Проводим GAM-GEE в пакете gamm4 R?
Я пытаюсь проанализировать некоторые визуальные данные трансектов организмов, чтобы создать модель распределения среды обитания. Как только организмы обнаружены, за ними следуют, когда точечные данные собираются за определенный промежуток времени. Из-за автокорреляции между этими "следует", я хочу использовать подход GAM-GEE, подобный подходу Pirotta et al. 2011, с использованием пакетов 'yags' и 'splines' (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 может быть наиболее подходящим. Важно отметить, что модели имеют бинарный ответ (присутствие организма в отсутствии поперечного разреза), и поэтому я считаю, что гамма-4 более уместна, чем гамма. В справке 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
>
Как я могу гарантировать, что автокорреляция действительно учитывается в каждом уникальном значении переменной "блок"? Какой самый простой способ интерпретировать вывод для "резюме (b$mer)"?
Результаты отличаются от обычной игры (пакет mgcv) с использованием тех же переменных и параметров без термина "correlation=...", указывающего на то, что происходит что-то другое.
Однако, когда я использую другую переменную для термина корреляции (сезона), я получаю ОДНО ЖЕ вывод:
> 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
>
Я просто хочу убедиться, что он правильно учитывает корреляцию внутри каждого значения для переменной "блок". Как мне сформулировать модель, чтобы сказать, что автокорреляция может существовать в каждом отдельном значении для блока, но предполагать независимость между блоками?
С другой стороны, я также получаю следующее предупреждающее сообщение после завершения модели для более крупных моделей (с большим количеством переменных, чем 2):
Warning message:
In mer_finalize(ans) : false convergence (8)
1 ответ
gamm4построен на вершинеlme4, что не позволяетcorrelationпараметр (в отличие отnlme, пакет, который лежит в основеmgcv::gamm).mgcv::gammобрабатывает двоичные данные, хотя использует PQL, который, как правило, менее точен, чем приближения Лапласа /GHQ, как вgamm4/lme4, К сожалению (!!) вы не получаете предупреждение о том, чтоcorrelationаргумент игнорируется (когда я пытаюсь простой пример, используяcorrelationспорить сlme4Я получаю предупреждение, но возможно, что дополнительный аргумент проглатывается где-то внутриgamm4).- Желаемая автокорреляционная структура ("автокорреляция может существовать в каждом отдельном значении для блока, но предполагать независимость между блоками") - это именно тот способ, которым кодируются корреляционные структуры.
nlme(и, следовательно, вmgcv::gamm). - я хотел бы использовать
mcgv::gammи предложил бы, если это вообще возможно, попробовать его на некоторых смоделированных данных с известной структурой (или использовать набор данных, предоставленный в дополнительном материале выше, и посмотреть, сможете ли вы воспроизвести их качественные выводы с помощью альтернативного подхода). - Stackru - это хорошо, но, возможно, в
r-sig-mixed-models@r-project.org