从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?

Answer 1

我不是对内部工作的期望,nlme::lme但我认为从模型中得到你想要的东西并不容易 - 该ranef()方法不允许返回随机效应的后验或条件方差,不同于lmer()和安装的模型的方法.

我想到了两个选择

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版本中它们的标准误差稍大.