проводя GAM-GEE в пакете gamm4 R?

Question

Я пытаюсь проанализировать некоторые визуальные данные обрезания организмов для создания модели распределения среды обитания. Как только организмы видны, они следуют за тем, как собираются данные о точке за определенный промежуток времени. Из-за автокорреляции среди этих «следует», я хочу использовать подход GAM-GEE, аналогичный подходу Pirotta et al. 2011, используя пакеты «yags» и «сплайны» (http://www.int-res.com/abstracts/meps/v436/p257-272/). Их R-скрипты показаны здесь (http://www.int-res.com/articles/suppl/m436p257_supp/m436p257_supp1-code.r). Я использовал этот код с ограниченным успехом и несколько проблем, с которыми модели не сходятся.

Ниже представлена ​​структура моих данных:

> str(dat2) 

'data.frame': 10792 obs. of 4 variables: 

$ dist_slag  : num 26475 26340 25886 25400 24934 ... 
$ Depth   : num -10.1 -10.5 -16.6 -22.2 -29.7 ... 
$ dolphin_presence: int 0 0 0 0 0 0 0 0 0 0 ... 


$ block   : int 1 1 1 1 1 1 1 1 1 1 ... 


> head(dat2) 

    dist_slag Depth dolphin_presence block 
1 26475.47 -10.0934    0  1 
2 26340.47 -10.4870    0  1 
3 25886.33 -16.5752    0  1 
4 25399.88 -22.2474    0  1 



5 24934.29 -29.6797    0  1 
6 24519.90 -26.2370    0  1 

Здесь приводится краткое изложение моей переменной блока (с указанием количества групп, для которых существует автокорреляции в пределах каждого блока

> summary(dat2$block) 
    Min. 1st Qu. Median Mean 3rd Qu. Max. 
    1.00 39.00 76.00 73.52 111.00 148.00 

Однако я хотел бы использовать пакет «gamm4», так как я больше знаком с пакетами и функциями профессора Саймона Вуда, и кажется, что gamm4 может быть наиболее подходящим. Важно отметить, что модели имеют двоичный ответ (отсутствие присутствия организма вдоль разреза), и поэтому я думаю, что gamm4 более уместен, чем gamm. В Гамм помощи, она обеспечивает следующий пример для автокорреляции в пределах факторов:

## more complicated autocorrelation example - AR errors 
## only within groups defined by `fac' 
e <- rnorm(n,0,sig) 
for (i in 2:n) e[i] <- 0.6*e[i-1]*(fac[i-1]==fac[i]) + e[i] 
y <- f + e 
b <- gamm(y~s(x,k=20),correlation=corAR1(form=~1|fac)) 

Следуя этому примеру, следующий код, который я использовал для моего набора данных

b <- gamm4(dolphin_presence~s(dist_slag)+s(Depth),random=(form=~1|block), family=binomial(),data=dat) 

Однако, рассматривая выход ((b $ gam)) и, в частности, резюме (b $ mer)), я либо не уверен, как интерпретировать результаты, либо не считаю, что автокорреляция внутри группы принимается во внимание.

> summary(b$gam) 

Family: binomial 
Link function: logit 

Formula: 
dolphin_presence ~ s(dist_slag) + s(Depth) 

Parametric coefficients: 


      Estimate Std. Error z value Pr(>|z|) 
    (Intercept) -13.968  5.145 -2.715 0.00663 ** 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Approximate significance of smooth terms: 


       edf Ref.df Chi.sq p-value  
s(dist_slag) 4.943 4.943 70.67 6.85e-14 *** 
s(Depth)  6.869 6.869 115.59 < 2e-16 *** 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 



R-sq.(adj) = 0.317glmer.ML score = 10504 Scale est. = 1   n = 10792 
> 

> summary(b$mer) 
Generalized linear mixed model fit by the Laplace approximation 


    AIC BIC logLik deviance 
10514 10551 -5252 10504 
Random effects: 
Groups Name   Variance Std.Dev. 
Xr  s(dist_slag) 1611344 1269.39 
Xr.0 s(Depth)  98622 314.04 
Number of obs: 10792, groups: Xr, 8; Xr.0, 8 



Fixed effects: 
       Estimate Std. Error z value Pr(>|z|) 
X(Intercept)  -13.968  5.145 -2.715 0.00663 ** 
Xs(dist_slag)Fx1 -35.871  33.944 -1.057 0.29063 
Xs(Depth)Fx1  3.971  3.740 1.062 0.28823 


--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Correlation of Fixed Effects: 
      X(Int) X(_)F1 
Xs(dst_s)F1 0.654  
Xs(Dpth)Fx1 -0.030 0.000 
> 

Как убедиться, что автокорреляция действительно учитываются в пределах каждого уникального значения «блок» переменный? Каков самый простой способ интерпретировать вывод для «summary (b $ mer)»?

Результаты отличаются от обычных gam (пакет mgcv) с использованием тех же переменных и параметров без условия «корреляция = ...», что указывает на то, что происходит что-то другое.

Однако, когда я использую другую переменную для термина корреляции (сезона), я получаю тот же результат:

> dat2 <- data.frame(dist_slag = dat$dist_slag, Depth = dat$Depth, dolphin_presence = dat$dolphin_presence, 

+ block = dat$block, season=dat$season) 
> head(dat2) 
     dist_slag Depth dolphin_presence block season 
1 26475.47 -10.0934    0  1  F 
2 26340.47 -10.4870    0  1  F 

3 25886.33 -16.5752    0  1  F 
4 25399.88 -22.2474    0  1  F 
5 24934.29 -29.6797    0  1  F 
6 24519.90 -26.2370    0  1  F 

> summary(dat2$season) 

    F S 
3224 7568 


> b <- gamm4(dolphin_presence~s(dist_slag)+s(Depth),correlation=corAR1(1, form=~1 | season), family=binomial(),data=dat2) 
> summary(b$gam) 

Family: binomial 
Link function: logit 


Formula: 
dolphin_presence ~ s(dist_slag) + s(Depth) 

Parametric coefficients: 
      Estimate Std. Error z value Pr(>|z|) 
    (Intercept) -13.968  5.145 -2.715 0.00663 ** 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 


Approximate significance of smooth terms: 
       edf Ref.df Chi.sq p-value  
s(dist_slag) 4.943 4.943 70.67 6.85e-14 *** 
s(Depth)  6.869 6.869 115.59 < 2e-16 *** 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 


R-sq.(adj) = 0.317glmer.ML score = 10504 Scale est. = 1   n = 10792 
> summary(b$mer) 
Generalized linear mixed model fit by the Laplace approximation 
    AIC BIC logLik deviance 

10514 10551 -5252 10504 
Random effects: 
Groups Name   Variance Std.Dev. 
Xr  s(dist_slag) 1611344 1269.39 
Xr.0 s(Depth)  98622 314.04 
Number of obs: 10792, groups: Xr, 8; Xr.0, 8 


Fixed effects: 
       Estimate Std. Error z value Pr(>|z|) 
X(Intercept)  -13.968  5.145 -2.715 0.00663 ** 
Xs(dist_slag)Fx1 -35.871  33.944 -1.057 0.29063 
Xs(Depth)Fx1  3.971  3.740 1.062 0.28823 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Correlation of Fixed Effects: 
      X(Int) X(_)F1 
Xs(dst_s)F1 0.654  
Xs(Dpth)Fx1 -0.030 0.000 
> 

Я просто хочу, чтобы убедиться, что он правильно позволяет корреляции в пределах каждого значения для переменная «block». Как сформулировать модель, чтобы сказать, что автокорреляция может существовать в пределах каждого отдельного значения для блока, но взять независимость от блоков?

На другой ноте, я также получаю следующее предупреждающее сообщение после завершения модели для больших моделей (с многими переменными, чем 2):

Warning message: 
In mer_finalize(ans) : false convergence (8) 

Answer 1

• gamm4 построен на вершине lme4, который делает не Разрешить параметр correlation (в отличие от nlme, пакет, который лежит в основе mgcv::gamm).mgcv::gammделает обрабатывает двоичные данные, хотя использует PQL, который, как правило, менее точен, чем приближения Лапласа/GHQ, как в gamm4/lme4. Несчастливо (!!), что вы не получаете предупреждение о том, что аргумент correlation игнорируется (когда я пробую простой пример, используя аргумент correlation с lme4, я получаю предупреждение, но возможно, что дополнительный аргумент проглатывается где-то внутри gamm4).
• Ваша желаемая структура автокорреляции («автокорреляция может существовать в пределах каждого отдельного значения для блока, но допускать независимость от блоков») в точности соответствует тому, как корреляционные структуры закодированы в nlme (и, следовательно, в mgcv::gamm).
• Я бы использовал mcgv::gamm и предположил бы, что, если это вообще возможно, вы попробуете его на некоторых смоделированных данных с известной структурой (или используйте набор данных, указанный в дополнительном материале выше, и посмотрите, сможете ли вы воспроизвести свои качественные выводы с помощью своих альтернативный подход).
• StackOverflow это хорошо, но, вероятно, более смешанная модель экспертиза на r-sig-mixed-models@r-project.org