gpt4 book ai didi

r - 在 gamm4 R 包中进行 GAM-GEE?

转载 作者:行者123 更新时间:2023-12-02 22:43:20 28 4
gpt4 key购买 nike

我正在尝试分析生物体的一些视觉横断面数据以生成栖息地分布模型。一旦发现生物,就会在给定的时间间隔收集点数据来跟踪它们。由于这些“跟随”之间的自相关,我希望利用类似于 Pirotta 等人的 GAM-GEE 方法。 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

这是我的 block 变量的摘要(指示每个 block 中存在自相关的组数

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

但是,我想使用'gamm4'包,因为我对Simon Wood教授的包和功能比较熟悉,看来gamm4可能是最合适的。重要的是要注意模型具有二元响应(生物体沿横断面存在或不存在),因此我认为 gamm4 比 gamm 更合适。在 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)

但是,通过检查输出 (summary(b$gam)),特别是 summary(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
>

我如何确保在“ block ”变量的每个唯一值中确实考虑了自相关?解释“summary(b$mer)”输出的最简单方法是什么?

结果确实不同于使用相同变量和参数但没有“correlation=...”项的普通 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 ”变量的每个值内的相关性。我如何制定模型来说明自相关可以存在于 block 的每个单个值中,但假设 block 之间是独立的?

另一方面,对于较大的模型(变量多于 2),我在完成模型后还收到以下警告消息:

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

最佳答案

  • gamm4 建立在 lme4 之上,它允许correlation 参数(相比之下到 nlme 包,它是 mgcv::gamm 的基础)。 mgcv::gamm 确实 处理二进制数据,尽管它使用 PQL,PQL 通常不如 gamm4/lme4 中的 Laplace/GHQ 近似准确.不幸的是 (!!) 你没有收到警告,告诉你 correlation 参数被忽略了(当我尝试一个使用 correlation 参数的简单示例时lme4,我确实收到警告,但额外的参数可能在 gamm4 中的某处被吞没了)。
  • 您想要的自相关结构(“自相关可以存在于 block 的每个单个值中,但假设 block 之间是独立的”)正是相关结构在 nlme 中编码的方式(因此在 mgcv::gamm).
  • 我会使用 mcgv::gamm,并建议如果可能的话,您可以在一些具有已知结构的模拟数据上尝试一下(或使用上面补充 Material 中提供的数据集和看看你是否可以用你的替代方法重现他们的定性结论)。
  • StackOverflow 很不错,但在 r-sig-mixed-models@r-project.org
  • 可能有更多混合模型专业知识

关于r - 在 gamm4 R 包中进行 GAM-GEE?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/10432671/

28 4 0
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com