r - 从 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() 拟合的模型的方法不同。和公司。

两个选项浮现在脑海

使用 gamm4:gamm4() 拟合模型，其中对象的混合模型面应与 se.ranef() 一起使用, 或

用 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 版本中稍大一些。