从R中的GAMM类随机效应中提取标准误差

Question

我建立了一个广义加性混合效应模型，包括一个固定效应和随机截距效应（这是一个分类变量）。在运行模型后，我可以使用ranef(m1$lme)$x[[1]]提取每个类别的随机截距。但是，当我尝试使用se.ranef(m1$lme)提取随机效应的标准误差时，该函数无法工作。尝试使用se.ranef(m1)和se.ranef(m1$gam)也不起作用。我不知道是否因为这些函数仅适用于lmer类的模型？

是否有人可以帮助我从“gamm”类中提取我的随机截距的标准误差？我想使用随机截距和标准误差来绘制我的gamm模型的最佳线性无偏预测器。

我的初始模型的形式为：gamm(y ~ s(z), random = list(x = ~1), data = dat)。

library(mgcv)
library(arm)

example <- gamm(mag ~ s(depth), random = list(stations = ~1), data = quakes)
summary(example$gam)

#Family: gaussian 
#Link function: identity 

#Formula:
 # mag ~ s(depth)

#Parametric coefficients:
#  Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  5.02300    0.04608     109   <2e-16 ***
#  ---
#  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Approximate significance of smooth terms:
#  edf Ref.df     F p-value    
#s(depth) 3.691  3.691 43.12  <2e-16 ***
#  ---
#  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
#R-sq.(adj) =  0.0725   
#Scale est. = 0.036163  n = 1000

ranef(example$lme)$stations[[1]] # extract random intercepts
#se.ranef(example$lme) # extract se of random intercepts - Problem line - doesn't work?

- bio123

Answer 1

我不是nlme::lme内部工作的专家，但我认为从该模型中获得所需内容并不容易——与lmer()等拟合的模型不同，ranef()方法不允许返回随机效应的后验或条件方差。

有两个选择：

1. 使用gamm4:gamm4()拟合模型，其中对象的混合模型面应该与se.ranef()一起使用；或者
2. 使用随机效应基础拟合gam()模型。

设置

library("mgcv")
library("gamm4")
library("arm")

library("ggplot2")
library("cowplot")
theme_set(theme_bw())

选项1：gamm4::gamm4()

这是将您的模型直接转换为gamm4::gamm4()所需语法的翻译。

quakes <- transform(quakes, fstations = factor(stations))

m1 <- gamm4::gamm4(mag ~ s(depth), random = ~ (1 | fstations),
                   data = quakes)
re1 <- ranef(m1$mer)[["fstations"]][,1]
se1 <- se.ranef(m1$mer)[["fstations"]][,1]

注意，我将stations转换为因子，因为mgcv::gam需要一个因子来拟合随机截距。

选项2：mgcv::gam() 对于这个选项，我们使用随机效应基础知识。惩罚样条模型理论表明，如果我们以特定方式书写数学公式，则该模型具有与混合模型相同的形式，基础部分的波动作为随机效应，无限平滑部分的基础用作固定效应。同样的理论也允许反向过程；我们可以制定一个完全受到惩罚的样条基础，这相当于一个随机效应。

m2 <- gam(mag ~ s(depth) + s(fstations, bs = "re"),
          data = quakes, method = "REML")

我们还需要做一些工作来获取“估计”的随机效应和标准误差。我们需要在fstations的水平上从模型中进行预测。我们还需要传入模型中其他术语的值，但由于模型是可加的，我们可以忽略它们的影响并提取随机效应。

newd <- with(quakes, data.frame(depth = mean(depth),
                                fstations = levels(fstations)))
p <- predict(m2, newd, type = "terms", se.fit = TRUE)
re2 <- p[["fit"]][ , "s(fstations)"]
se2 <- p[["se.fit"]][ , "s(fstations)"]

这些选项有何区别？

re <- data.frame(GAMM = re1, GAM = re2)
se <- data.frame(GAMM = se1, GAM = se2)

p1 <- ggplot(re, aes(x = GAMM, y = GAM)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1) +
  coord_equal() +
  labs(title = "Random effects")

p2 <- ggplot(se, aes(x = GAMM, y = GAM)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1) +
  coord_equal() +
  labs(title = "Standard errors")

plot_grid(p1, p2, nrow = 1, align = "hv")

“估计值”是等效的，但在GAM版本中，它们的标准误差略大。

- Gavin Simpson

---

据我所知，这里关于lme中随机效应值的标准误差的所有陈述都是正确的（在?ranef.lme或?intervals.lme中也没有提到可能性）。您可能希望评论一下为什么在“估计值”周围使用引号（即，随机效应的预测/BLUPs/条件模式不是以标准频率学意义上的估计值）。 - Ben Bolker