Предсказатели масштабирования в lme4 glmer не разрешают предупреждения о собственных значениях;как и альтернативная оптимизация - PullRequest
Новая
       37

Предсказатели масштабирования в lme4 glmer не разрешают предупреждения о собственных значениях;как и альтернативная оптимизация

0 голосов
/ 18 декабря 2018

Я анализирую данные (включены ниже) с использованием функции lme4 glmer в R. Модель, которую я строю, состоит из распределенной по Пуассону переменной отклика (obs), одного случайного фактора (area), одно непрерывное смещение (duration), пять непрерывных фиксированных эффектов (can_perc, can_n, time, temp, cloud_cover) и один биномиальный фиксированный коэффициент эффекта (burnt).Перед подгонкой модели я проверил коллинеарность и удалил все коллинеарные переменные.

Исходная модель: 

q1 = glmer(obs ~ can_perc + can_n  + time * temp + 
           cloud_cover + factor(burnt) + (1|area) + offset(dat$duration), 
           data=dat, family=poisson, na.action = na.fail)

(Примечание: мне нужно указать na.action как 'na.fail'как я хочу dredge() модель позже, и это требуется для этого.)

Запуск модели выдает следующее предупреждение:

"Гессиан численно сингулярен: параметрыопределяется не однозначно "

В аналогичных вариациях модели я также получил предупреждение:

" In checkConv (attr (opt, "производные"), opt$ par, ctrl = control $ checkConv,: Модель почти неопознаваема: большое отношение собственных значений - Масштабировать переменные? "

Из моего ограниченного понимания совета здесь https://rdrr.io/cran/lme4/man/troubleshooting.html и в других местах, обаиз этих предупреждений отражают аналогичную проблему: гессиан (матрица обратной кривизны) имеет большое собственное значение, что указывает на то, что (в пределах числовых допусков) поверхность в некотором направлении совершенно плоская.предупреждения и ссылку, я перемасштабировал все непрерывные переменные предиктора, используя scale().Я также масштабировал переменную смещения (я пробовал как с масштабированием, так и без него).Модель с масштабированными переменными предиктора здесь: 

q2 = glmer(obs ~ s.can_perc + s.can_n  + s.time * s.temp + 
           s.cloud_cover + factor(burnt) + (1|area) +
           offset(dat$s.duration), 
           data=dat, family=poisson, na.action = na.fail)

Однако я еще не избежал собственных значений!Масштабированная модель выдает два предупреждения:

"невозможно оценить масштабированный градиент"
"Модель не может сходиться: вырожденный гессиан с 1 отрицательным собственным значением"

У меня естьМного разыскивал в Интернете и не смог найти другой случай / решение того, как справиться с проблемами собственных значений после масштабирования предикторов, кроме проверки того, что модель не была ошибочно определена.

Попытки устранить предупреждения /улучшить оптимизацию:

На основании этих страниц / документов: https://cran.r -project.org / web / packages / lme4 / lme4.pdf

https://stats.stackexchange.com/questions/164457/r-glmer-warnings-model-fails-to-converge-model-is-nearly-unidentifiable

https://rdrr.io/cran/lme4/man/isSingular.html

https://stats.stackexchange.com/questions/242109/model-failed-to-converge-warning-in-lmer

и другие,

У меня есть:

  1. провереноспецификации модели и данные об ошибках (ничего, что я могу видеть - я что-то пропустил?)

  2. проверено на сингулярность с помощью is_singular(x, tol = 1e-05) (каким-то образом этот вызов функции эволюционировал с isSingular() дотекущая форма?): все модели дают ЛОЖЬ.

  3. проверенная сходимость, мeasure с converge_ok(q2, tolerance = 0.001): все модели дают FALSE, если я существенно не увеличил допуск;однако они существенно различаются по степени сходимости.

  4. пробовали различные оптимизаторы / методы оценки модели следующим образом:

    • a) glmerControl(optimizer = "bobyqa") and glmerControl(optimizer ="Nelder_Mead")
    • b) glmerControl(optimizer ='optimx', optCtrl=list(method='nlminb'))
    • c) bobyqa, Nelder_Mead, optimx.nlminb, optimx.L-BFGS-B, nloptwrap.NLOPT_LN_NELDERMEAD, nloptwrap.NLOPT_LN_BOBYQA и функция 1095 * с использованием функции N95пакет optimx.

Вот код: 

# singularity and convergence for first two models:
is_singular(s1, tol = 1e-05) # FALSE (a good thing?)
converge_ok(s1, tol = 1e-05) # FALSE (a bad thing?) 0.0259109730912352

is_singular(s2, tol = 1e-05) # FALSE (a good thing?)
converge_ok(s2, tol = 1e-05) # FALSE (a bad thing?) 0.0023434329028163
# I looked at singularity and converge measures for the others below, but omitted for brevity.

# Alternate optimisations for q1:
q1.bobyqa = glmer(obs ~ can_perc + can_n  + time * temp + cloud_cover + factor(burnt) + (1|area) + offset(dat$duration), data=dat, family=poisson, na.action = na.fail, glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e5)))
# Warning 1: unable to evaluate scaled gradient
# Warning 2: Model failed to converge: degenerate  Hessian with 1 negative eigenvalues

q1.neldermead = glmer(obs ~ can_perc + can_n  + time * temp + cloud_cover + factor(burnt) +  (1|area) + offset(dat$duration), data=dat, family=poisson, na.action = na.fail,  glmerControl(optimizer ="Nelder_Mead", optCtrl = list(maxfun = 2e5)))
# Warning: unable to evaluate scaled gradient Hessian is numerically singular: parameters are not uniquely determined

q1.nlminb = glmer(obs ~ can_perc + can_n  + time * temp + cloud_cover + factor(burnt) + (1|area) + offset(dat$duration), data=dat, family=poisson, na.action = na.fail, glmerControl(optimizer ='optimx', optCtrl=list(method='nlminb')))
# Warning: Parameters or bounds appear to have different scalings. This can cause poor performance in optimization. 
# It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.convergence code 9999 from optimxError in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev,  :   (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate

all_fit(q1)

# Alternate optimisations for q2:
q2.bobyqa = glmer(obs ~ s.can_perc + s.can_n  + s.time * s.temp + s.cloud_cover + factor(burnt) +  (1|area) + offset(dat$s.duration), data=dat, family=poisson, na.action = na.fail, glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e5)))
# Warning 1: unable to evaluate scaled gradient
# Warning 2: Model failed to converge: degenerate  Hessian with 1 negative eigenvalues

q2.neldermead = glmer(obs ~ s.can_perc + s.can_n  + s.time * s.temp + s.cloud_cover + factor(burnt) + (1|area) + offset(dat$s.duration), data=dat, family=poisson, na.action = na.fail, glmerControl(optimizer ="Nelder_Mead", optCtrl = list(maxfun = 2e5)))
# Warning: unable to evaluate scaled gradient Hessian is numerically singular: parameters are not uniquely determined

q2.nlminb = glmer(obs ~ s.can_perc + s.can_n  + s.time * s.temp + s.cloud_cover + factor(burnt) +  (1|area) + offset(dat$s.duration), data=dat, family=poisson, na.action = na.fail, control = glmerControl(optimizer ='optimx', optCtrl=list(method='nlminb')))
# Warning: Model is nearly unidentifiable: large eigenvalue ratio - Rescale variables?

all_fit(q2)

Выход из вышеприведенного кода для немасштабированной модели (q1): 

is_singular(s1, tol = 1e-05) # FALSE (a good thing?)
[1] FALSE
converge_ok(s1, tol = 1e-05) # FALSE (a bad thing?) 0.0259109730912352
0.0259109730912352 
             FALSE 
is_singular(s2, tol = 1e-05) # FALSE (a good thing?)
[1] FALSE
alternate optimisations for original model:
q1.bobyqa = glmer(obs ~ can_perc + can_n  + time * temp + cloud_cover + factor(burnt) + (1|area) + offset(dat$duration), data=dat, family=poisson, na.action = na.fail, glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e5)))
unable to evaluate scaled gradientModel failed to converge: degenerate  Hessian with 1 negative eigenvalues

alternate optimisations for original model:
q1.bobyqa = glmer(obs ~ can_perc + can_n  + time * temp + cloud_cover + factor(burnt) + (1|area) + offset(dat$duration), data=dat, family=poisson, na.action = na.fail, glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e5)))
    unable to evaluate scaled gradientModel failed to converge: degenerate      Hessian with 1 negative eigenvalues
    q1.neldermead = glmer(obs ~ can_perc + can_n  + time * temp + cloud_cover + factor(burnt) + (1|area) + offset(dat$duration), data=dat, family=poisson, na.action = na.fail, glmerControl(optimizer ="Nelder_Mead", optCtrl = list(maxfun = 2e5)))
    unable to evaluate scaled gradient Hessian is numerically singular: parameters are not uniquely determined

all_fit(q1)
bobyqa. : unable to evaluate scaled gradientModel failed to converge:     degenerate  Hessian with 1 negative eigenvalues[OK]
Nelder_Mead. : unable to evaluate scaled gradient Hessian is numerically singular: parameters are not uniquely determined[OK]
optimx.nlminb : Parameters or bounds appear to have different scalings.
This can cause poor performance in optimization. 
It is important for derivative free methods like BOBYQA, UOBYQA,     NEWUOA.convergence code 9999 from optimxParameters or bounds appear to have different scalings.
This can cause poor performance in optimization. 
It is important for derivative free methods like BOBYQA, UOBYQA,     NEWUOA.convergence code 9999 from optimx[ERROR]
optimx.L-BFGS-B : Parameters or bounds appear to have different scalings.
This can cause poor performance in optimization. 
It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.convergence code 9999 from optimxParameters or bounds appear to have different scalings.
This can cause poor performance in optimization. 
It is important for derivative free methods like BOBYQA, UOBYQA,     NEWUOA.convergence code 9999 from optimx[ERROR]
nloptwrap.NLOPT_LN_NELDERMEAD : [ERROR]
nloptwrap.NLOPT_LN_BOBYQA : [ERROR]
nmkbw. : [ERROR]

$`bobyqa.`
    Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
     Family: poisson  ( log )
    Formula: obs ~ can_perc + can_n + time * temp + cloud_cover + factor(burnt) +  (1 | area) + offset(dat$duration)
       Data: dat
      AIC       BIC    logLik  deviance  df.resid 
     311.0473  330.3356 -146.5237  293.0473        54 
    Random effects:
     Groups Name        Std.Dev.
     area   (Intercept) 1.992   
    Number of obs: 63, groups:  area, 8
    Fixed Effects:
(Intercept)  can_perc  can_n   time     temp  
-67.4998    -1.3180    0.0239    4.8025    1.7793  
cloud_cover  factor(burnt)unburnt             time:temp  
             -0.3813               18.5676               -0.1748  
convergence code 0; 2 optimizer warnings; 0 lme4 warnings 

$Nelder_Mead.
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
     Family: poisson  ( log )
    Formula: obs ~ can_perc + can_n + time * temp + cloud_cover + factor(burnt) +      (1 | area) + offset(dat$duration)
       Data: dat
      AIC       BIC    logLik  deviance  df.resid 
     311.0473  330.3356 -146.5237  293.0473        54 
    Random effects:
     Groups Name        Std.Dev.
     area   (Intercept) 1.992   
    Number of obs: 63, groups:  area, 8
    Fixed Effects:
             (Intercept)              
can_perc      can_n       time      temp  
-67.48057    -1.31791    0.02389    4.80463    1.78012  
cloud_cover  factor(burnt)unburnt    time:temp  
-0.38118    18.52637      -0.17483  
convergence code 0; 2 optimizer warnings; 0 lme4 warnings 

$optimx.nlminb
<std::runtime_error in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat,    compDev = compDev,     grpFac = fac, maxit = maxit, verbose = verbose): (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate>

$`optimx.L-BFGS-B`
<std::runtime_error in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev,     grpFac = fac, maxit = maxit, verbose = verbose): (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate>

$nloptwrap.NLOPT_LN_NELDERMEAD
<simpleError in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev,     grpFac = fac, maxit = maxit, verbose = verbose): Downdated VtV is not positive definite>

$nloptwrap.NLOPT_LN_BOBYQA
<simpleError in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev,     grpFac = fac, maxit = maxit, verbose = verbose): Downdated VtV is not positive definite>

$nmkbw.
<std::runtime_error in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev,     grpFac = fac, maxit = maxit, verbose = verbose): (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate>

Выход из вышеприведенного кода для масштабированной модели (q2): 

alternate optimisations for q2:
q2.bobyqa = glmer(obs ~ s.can_perc + s.can_n  + s.time * s.temp + s.cloud_cover + factor(burnt) + (1|area) + offset(dat$s.duration), data=dat, family=poisson, na.action = na.fail, glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e5)))
Model is nearly unidentifiable: large eigenvalue ratio - Rescale variables?
q2.neldermead = glmer(obs ~ s.can_perc + s.can_n  + s.time * s.temp + s.cloud_cover + factor(burnt) + (1|area) + offset(dat$s.duration), data=dat, family=poisson, na.action = na.fail, glmerControl(optimizer ="Nelder_Mead", optCtrl = list(maxfun = 2e5)))
unable to evaluate scaled gradientModel failed to converge: degenerate  Hessian with 1 negative eigenvalues

all_fit(q2)
bobyqa. : Model is nearly unidentifiable: large eigenvalue ratio
 - Rescale variables?[OK]
Nelder_Mead. : unable to evaluate scaled gradientModel failed to converge: degenerate  Hessian with 1 negative eigenvalues[OK]
optimx.nlminb : Model is nearly unidentifiable: large eigenvalue ratio
 - Rescale variables?[OK]
optimx.L-BFGS-B : unable to evaluate scaled gradientModel failed to converge: degenerate  Hessian with 1 negative eigenvalues[OK]
nloptwrap.NLOPT_LN_NELDERMEAD : [ERROR]
nloptwrap.NLOPT_LN_BOBYQA : [ERROR]
nmkbw. : [ERROR]
$`bobyqa.`
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: poisson  ( log )
Formula: n_shreiberi ~ s.can_perc + s.can_n + s.time * s.temp + s.cloud_cover +  
    factor(burnt) + (1 | area) + offset(dat$s.duration)
   Data: dat
      AIC       BIC    logLik  deviance  df.resid 
 316.8412  336.1294 -149.4206  298.8412        54 
Random effects:
 Groups Name        Std.Dev.
 area   (Intercept) 2.523   
Number of obs: 63, groups:  area, 8
Fixed Effects:
(Intercept)    s.can_perc    s.can_n    s.time    s.temp  
-18.19816    -0.22152    0.45839    0.05239    -0.24983  
       s.cloud_cover  factor(burnt)unburnt         s.time:s.temp  
            -0.19691              17.92390              -0.13948  
convergence code 0; 1 optimizer warnings; 0 lme4 warnings 

$Nelder_Mead.
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: poisson  ( log )
Formula: n_shreiberi ~ s.can_perc + s.can_n + s.time * s.temp + s.cloud_cover +  
    factor(burnt) + (1 | area) + offset(dat$s.duration)
   Data: dat
      AIC       BIC    logLik  deviance  df.resid 
 316.8408  336.1290 -149.4204  298.8408        54 
Random effects:
 Groups Name        Std.Dev.
 area   (Intercept) 2.524   
Number of obs: 63, groups:  area, 8
Fixed Effects:
         (Intercept)            s.can_perc               s.can_n                s.time                s.temp  
           -19.29632              -0.22153               0.45840               0.05241              -0.24990  
       s.cloud_cover  factor(burnt)unburnt         s.time:s.temp  
            -0.19692              19.02091              -0.13949  
convergence code 0; 2 optimizer warnings; 0 lme4 warnings 

$optimx.nlminb
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: poisson  ( log )
Formula: n_shreiberi ~ s.can_perc + s.can_n + s.time * s.temp + s.cloud_cover +  
    factor(burnt) + (1 | area) + offset(dat$s.duration)
   Data: dat
      AIC       BIC    logLik  deviance  df.resid 
 316.8412  336.1294 -149.4206  298.8412        54 
Random effects:
 Groups Name        Std.Dev.
 area   (Intercept) 2.523   
Number of obs: 63, groups:  area, 8
Fixed Effects:
         (Intercept)            s.can_perc               s.can_n                s.time                s.temp  
           -18.23626              -0.22152               0.45839               0.05239              -0.24983  
       s.cloud_cover  factor(burnt)unburnt         s.time:s.temp  
            -0.19691              17.96199              -0.13948  
convergence code 0; 1 optimizer warnings; 0 lme4 warnings 

$`optimx.L-BFGS-B`
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: poisson  ( log )
Formula: n_shreiberi ~ s.can_perc + s.can_n + s.time * s.temp + s.cloud_cover +  
    factor(burnt) + (1 | area) + offset(dat$s.duration)
   Data: dat
      AIC       BIC    logLik  deviance  df.resid 
 316.8412  336.1294 -149.4206  298.8412        54 
Random effects:
 Groups Name        Std.Dev.
 area   (Intercept) 2.524   
Number of obs: 63, groups:  area, 8
Fixed Effects:
         (Intercept)            s.can_perc               s.can_n                s.time                s.temp  
           -18.23581              -0.22155               0.45841               0.05242              -0.24997  
       s.cloud_cover  factor(burnt)unburnt         s.time:s.temp  
            -0.19694              17.96246              -0.13943  
convergence code 0; 2 optimizer warnings; 0 lme4 warnings 

$nloptwrap.NLOPT_LN_NELDERMEAD
<simpleError in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev,     grpFac = fac, maxit = maxit, verbose = verbose): Downdated VtV is not positive definite>

$nloptwrap.NLOPT_LN_BOBYQA
<simpleError in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev,     grpFac = fac, maxit = maxit, verbose = verbose): Downdated VtV is not positive definite>

$nmkbw.
<simpleError in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev,     grpFac = fac, maxit = maxit, verbose = verbose): Downdated VtV is not positive definite>

Данные:

Набор данных доступен по этой ссылке: https://www.dropbox.com/s/ud50uatztjq4bh9/20181217%20Surveys%20simplified%20data%20for%20stackX.xlsx?dl=0

Заключение и запрос:

Мне кажется, что ни один из этих альтернативных методов оптимизации также не удался;на самом деле, некоторые из них, похоже, вызвали другие предупреждения / ошибки, которые унесли бы меня в другую кроличью нору.

Может кто-нибудь посоветовать, как я могу продвинуться в подборе этих моделей?У меня нет намерения, чтобы они были конечными моделями, а скорее дноуглубить их, а затем выбрать оптимальные / лучшие модели из различных альтернативных моделей подмножеств.

1 Ответ

0 голосов
/ 21 декабря 2018

tl; dr Это похоже на случай полного разделения ;у вас нет никаких положительных результатов в вашем «обожженном» состоянии.Вам не нужно обязательно беспокоиться об этом - сравнения AIC все еще должны быть достаточно надежными - но вы, возможно, захотите понять, что происходит, прежде чем продолжить.Эта проблема (и способы ее устранения) обсуждаются в соответствующем разделе FAQ по GLMM (и есть множество соответствующих вопросов / ответов по CrossValidated ).

Откуда я знаю?Вот коэффициенты: 

  (Intercept)       s.can_perc               s.can_n                s.time                s.temp  
   -19.29632          -0.22153               0.45840               0.05241              -0.24990  
       s.cloud_cover  factor(burnt)unburnt         s.time:s.temp  
            -0.19692              19.02091              -0.13949

Каждый раз, когда у вас есть коэффициенты в (биномиальном или пуассоновском) GLM, которые больше (в абсолютном значении), чем 8-10, вам нужно беспокоиться (если вы не смотрите накоэффициент числового ковариата, который измеряется в очень больших единицах, например, если вы смотрите на количество углерода на заднем дворе в единицах гигатонн).Это означает, что изменение на одну единицу в предикторе вызывает (скажем) изменение на 10 единиц в log-odds (для модели биномиального / логит-линкового соединения), например, с вероятностью от 0,006 (plogis(-5)) до 0,994 (plogis(5)).В вашем случае, перехват составляет -19,29, поэтому при нулевых значениях всех предикторов в сожженном состоянии вы получите вероятность 4,2e-9.Другой огромный коэффициент для unburnt (19.02), поэтому при нулевых значениях всех предикторов в несгоревшем (несгоревшем?) Состоянии вы получите plogis(-19.29+19.02) = 0,43.

...