R random.effects GAM 中的随机效应

说明

为了估计的目的，GAM 的平滑分量可以被视为随机效应。这意味着可以通过两种方式将更传统的随机效应项合并到 GAM 中。第一种方法将所有平滑转换为适合通过标准混合建模软件估计的固定和随机分量。一旦 GAM 处于这种形式，那么传统的随机效应就很容易添加，并且整个模型被估计为一般的混合模型。 gamm4 包中的gamm 和gamm4 以这种方式运行。

第二种方法以与平滑表示相同的方式表示 GAM 中的传统随机效应 - 作为惩罚回归项。此方法可以与gam通过利用s(...,bs="re")模型中的项：参见smooth.construct.re.smooth.spec，了解完整详细信息。基本思想是，例如，s(x,z,g,bs="re")生成一个 i.i.d.模型矩阵的高斯随机效应由下式给出model.matrix(~x:z:g-1)——原则上，这些术语可以采用任意数量的参数。这种简单的方法足以实现各种常用的随机效应结构。例如如果g那么是一个因子s(g,bs="re")为每个级别生成一个随机系数g，随机系数全部建模为 i.i.d.普通的。如果g是一个因子并且x是数字，那么s(x,g,bs="re")产生一个 i.i.d.与响应相关的正态随机斜率x对于每个级别g.如果h那么是另一个因子s(h,g,bs="re")产生通常的 i.i.d.普通的g-h相互作用。请注意，基于 Wood (2013)，还对此类项实施了相当有用的零随机效应近似检验。如果精度矩阵在乘法常数内已知，则可以通过以下方式提供xt的论点s.看mgcv smooth.construct.re.smooth.spec了解详细信息和示例。一些模型需要相同随机效应的不同水平之间的差异：这些可以按照中所述来实现linear.functional.terms.

或者，但不太直接，可以使用 gam 的 paraPen 参数：请参阅 gam.models 。如果平滑参数估计是通过 ML 或 REML(例如 gam(...,method="REML") )进行的，那么这种方法是完全传统的基于似然的随机效应处理。

gam 对于拟合具有大量随机效应的模型可能会很慢，因为它没有利用稀疏性，而稀疏性通常是参数随机效应的一个特征。它不能用于系数多于数据的模型。然而，当随机效应的数量适中时，gam 通常比 gamm 或 gamm4 更快、更可靠。

为了方便使用 gam 的随机效应，gam.vcomp 是一个用于将平滑参数转换为方差分量的实用程序。如果平滑度估计是通过 ML 或 REML 进行的，它还提供置信区间。

请注意，将随机效应视为平滑并不能消除与测试方差分量是否等于零相关的常见问题：请参阅 summary.gam 和 anova.gam 。

例子

## see also examples for gam.models, gam.vcomp, gamm
## and smooth.construct.re.smooth.spec

## simple comparison of lme and gam
require(mgcv)
require(nlme)
b0 <- lme(travel~1,data=Rail,~1|Rail,method="REML") 

b <- gam(travel~s(Rail,bs="re"),data=Rail,method="REML")

intervals(b0)
gam.vcomp(b)
anova(b)
plot(b)

## simulate example...
dat <- gamSim(1,n=400,scale=2) ## simulate 4 term additive truth

fac <- sample(1:20,400,replace=TRUE)
b <- rnorm(20)*.5
dat$y <- dat$y + b[fac]
dat$fac <- as.factor(fac)

rm1 <- gam(y ~ s(fac,bs="re")+s(x0)+s(x1)+s(x2)+s(x3),data=dat,method="ML")
gam.vcomp(rm1)

fv0 <- predict(rm1,exclude="s(fac)") ## predictions setting r.e. to 0
fv1 <- predict(rm1) ## predictions setting r.e. to predicted values
## prediction setting r.e. to 0 and not having to provide 'fac'...
pd <- dat; pd$fac <- NULL
fv0 <- predict(rm1,pd,exclude="s(fac)",newdata.guaranteed=TRUE)

## Prediction with levels of fac not in fit data.
## The effect of the new factor levels (or any interaction involving them)
## is set to zero.
xx <- seq(0,1,length=10)
pd <- data.frame(x0=xx,x1=xx,x2=xx,x3=xx,fac=c(1:10,21:30))
fv <- predict(rm1,pd)
pd$fac <- NULL
fv0 <- predict(rm1,pd,exclude="s(fac)",newdata.guaranteed=TRUE)

作者

Simon Wood <simon.wood@r-project.org>

参考

Wood, S.N. (2013) A simple test for random effects in regression models. Biometrika 100:1005-1010

Wood, S.N. (2011) Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. Journal of the Royal Statistical Society (B) 73(1):3-36

Wood, S.N. (2008) Fast stable direct fitting and smoothness selection for generalized additive models. Journal of the Royal Statistical Society (B) 70(3):495-518

Wood, S.N. (2006) Low rank scale invariant tensor product smooths for generalized additive mixed models. Biometrics 62(4):1025-1036

也可以看看

gam.vcomp、gam.models、smooth.terms、smooth.construct.re.smooth.spec、gamm