R gam.models 指定广义加性模型

说明

本页旨在提供有关如何指定 GAM 的更多信息。 GAM 是一种 GLM，其中线性预测器部分取决于预测器的平滑函数和(可能是虚拟)预测器的平滑函数的线性函数之和。

具体来说，让 y_i 表示具有均值 \mu_i 和指数族分布的独立随机变量，或者不适合使用 quasi-likelihood 方法的已知均值方差关系。那么 GAM 的线性预测器具有类似的结构

其中g是已知的平滑单调‘link’函数，{\bf X}_i\beta是线性预测器的参数部分，x_j是预测变量，f_j是平滑函数，L_i是f_3 。当然，可能有多个线性函数项，也可能没有。

这里的关键思想是，响应对预测变量的依赖性可以表示为参数 sub-model 加上一个或多个预测变量的一些平滑函数(函数)的总和。因此，该模型相对于严格参数线性或广义线性模型相当灵活，但仍然比完全通用模型具有更多的结构，后者表示响应只是所有协变量的一些平滑函数。

请注意重要的一点。为了使模型可识别，平滑函数通常必须被限制为具有零均值(通常接管协变量值的集合)。如果涉及平滑的项在其跨度中包含常数函数，则需要该约束。 gam 始终应用此类约束，除非存在 by 变量，在这种情况下，将评估是否需要约束(见下文)。

以下部分讨论为 gam 指定模型结构。分发和链接函数的规范是使用 gam 的 family 参数完成的，其工作方式与 glm 相同。因此，本页重点关注 gam 的模型公式。

具有简单平滑项的模型

考虑示例模型。

其中响应变量 y_i 具有期望 \mu_i 和 g 是链接函数。

gam其公式为
y ~ x1 + x2 + s(x3) + s(x4,x5).
这将使用平滑的默认基础(每个平滑的薄板回归样条基础)，并自动选择两个平滑的有效自由度。平滑基础的维度也被赋予默认值(基础的维度设置基础的最大可能自由度的上限 - 该限制通常比基础维度小一)。有关如何控制平滑的完整详细信息，请参见s和te，关于基本尺寸选择的进一步讨论可以在choose.k。目前假设我们希望将第一个平滑的基础更改为维度为 20 的三次回归样条基础，同时将第二项固定为 25 自由度。适当的公式是：
y ~ x1 + x2 + s(x3,bs="cr",k=20) + s(x4,x5,k=26,fx=TRUE).

上面假设x_{4}和x_5自然地具有相似的尺度(例如，它们可能是坐标)，因此各向同性平滑是合适的。如果这个假设是错误的，那么张量积平滑可能会更好(参见te)。
y ~ x1 + x2 + s(x3) + te(x4,x5)
将生成一个平滑的张量积x_{4}和x_5。默认情况下，该平滑的基础维度为 25，并使用三次回归样条边际。改变默认值很容易。例如
y ~ x1 + x2 + s(x3) + te(x4,x5,bs=c("cr","ps"),k=c(6,7))
指定张量积应使用 6 阶三次回归样条边际和 7 阶 P-spline 边际来创建基本维度为 42 的平滑。

嵌套项/函数方差分析

有时，指定具有主效应+交互结构的平滑模型很有趣，例如

或者

例如：此类模型应使用ti模型公式中的项。例如：
y ~ ti(x) + ti(z) + ti(x,z)， 或者
y ~ ti(x) + ti(z) + ti(v) + ti(x,z) + ti(x,v) + ti(z,v)+ti(x,z,v).
ti项与适当排除的成分主效应产生相互作用。 (其实没必要用ti这里的主要影响术语，s也可以使用术语。)

gam 允许 te 和 s 进行嵌套(或 ‘overlap’)，并自动生成辅助条件以使此类模型可识别，但生成的模型比使用 ti 项构建的模型稳定性和可解释性要差得多。

‘by’变量

by 变量是构建“varying-coefficient 模型”(地理回归模型)以及使用因子或参数项平滑 ‘interact’ 的方法。它们也是指定平滑的一般线性函数的关键。

用于指定平滑的 s 和 te 术语接受参数 by ，该参数是与平滑协变量具有相同维度的数值或因子变量。如果 by 变量是数字，则其 i^{th} 元素乘以与相关平滑项对应的模型矩阵的 i^{th} 行。

影响平滑交互(另请参阅factor.smooth.interaction)。如果 by 变量是 factor，则它会为因子的每个级别生成指示向量，除非它是 ordered 因子。在无序情况下，然后为每个因子水平复制平滑项的模型矩阵，并且每个副本的行与其指示变量的相应行相乘。每个因子级别的平滑度惩罚也会重复。简而言之，为每个因子级别生成不同的平滑(s 和 te 的 id 参数可用于强制所有此类平滑具有相同的平滑参数)。 ordered by 变量以相同的方式处理，只是不会为有序因子的第一级生成平滑(请参见下面的 b3 示例)。当相同的平滑在模型中多次出现且具有不同的因子 by 变量时，这对于设置可识别模型非常有用。

作为一个例子，考虑模型

其中f是一个光滑函数，并且z_i是一个数值变量。适当的公式是：
y ~ s(x,by=z)
- 这by参数确保平滑函数乘以协变量z。请注意，当使用因子变量时，中心约束将应用于平滑，这通常意味着因子变量也应作为参数项包含在内。

下面的示例代码还说明了因子 by 变量的使用。

by 变量可以作为数值矩阵提供，作为指定一般线性函数项的一部分。

如果 by 变量存在并且是数字(而不是因子)，则相应的平滑仅在以下情况下受到可识别性约束：(i) by 变量是常量向量，或者 (ii) 对于矩阵 by 变量，L ，如果 L%*%rep(1,ncol(L)) 是常量，或者 (iii) 如果用户定义的平滑构造函数显式提供可识别性约束，并且该约束具有属性 "always.apply" 。

使用 ‘id’ 进行平滑链接

有时需要坚持不同的平滑项具有相同的平滑度。这可以通过使用 s 或 te 术语的 id 参数来完成。共享 id 的平滑将具有相同的平滑参数。实际上，只有当平滑使用相同的基函数时，这才有意义，并且默认行为是强制发生这种情况：共享 id 的所有平滑都具有与 id 发生的第一个平滑相同的基函数。请注意，如果您希望每个平滑函数具有完全相同的函数，那么最好通过使用“线性函数项”下涵盖的求和约定来实现。

作为示例，假设 E(y_i)\equiv\mu_i 和

但是那个f_1和f_3应具有相同的平滑参数(并且x_2和x_3是在不同的尺度上)。然后gam公式
y ~ s(x1,id=1) + te(x_2,x3) + s(x4,id=1)
就会达到想要的结果。id可以是数字或字符串。给予id到一个带有因子的项by变量导致因子每个级别的平滑具有相同的平滑参数。

gamm 不支持平滑项 id 。

线性泛函项

一般线性函数项在样条文献中有着悠久的历史，包括在惩罚 GLM 上下文中(例如参见 Wahba 1990)。这些术语涵盖不同的系数模型/地理回归、函数 GLM(即具有函数预测变量的 GLM)、GLASS 模型等，并允许例如针对聚合协变量值进行平滑。

此类条款的实施于mgcv对平滑项使用简单的“求和约定”：如果以矩阵形式提供平滑的协变量，则隐含对矩阵列上的评估平滑求和。每个协变量矩阵和任意by变量矩阵必须具有相同的维度。例如，考虑术语
s(X,Z,by=L)
其中X,Z和L是n \times p矩阵。让f表示指定的薄板回归样条。由此产生的贡献i^{\rm th}线性预测器的元素是

如果未提供 L，则其所有元素均取为 1。在 R 代码术语中，令 F 表示通过评估 X 和 Z 中的值的平滑度而获得的 n \times p 矩阵。那么该项对线性预测器的贡献是rowSums(L*F)(请注意，这里是逐个元素相乘！)。

求和约定适用于te 项和s 项。 linear.functional.terms 中提供了更多详细信息和示例。

随机效应

可以使用 s(...,bs="re") 术语(请参阅 smooth.construct.re.smooth.spec )或下面介绍的 gam 的 paraPen 参数将随机效应添加到 gam 模型中。有关更多详细信息，请参阅gam.vcomp、random.effects 和smooth.construct.re.smooth.spec。另一种方法是使用 gamm 的方法。

惩罚参数项

如果添加平滑类、平滑恒等式、by 变量和求和约定的能力仍不足以准确实现您所需的惩罚 GLM，gam 还允许您惩罚模型公式中的参数项。这对于允许在公式中包含一个或多个矩阵项以及每个矩阵项的一系列二次罚矩阵非常有用。

假设你已经建立了一个模型矩阵\bf X，并想要惩罚相应的系数，\beta并有两项处罚\beta^T {\bf S}_1 \beta和\beta^T {\bf S}_2 \beta。那么像下面这样的东西是合适的：
gam(y ~ X - 1,paraPen=list(X=list(S1,S2)))
paraPen参数应该是一个列表，其中元素的名称与被惩罚的术语相对应。的每个元素paraPen本身是一个列表，带有可选元素L,rank和sp：所有其他元素必须是惩罚矩阵。如果存在的话，rank是一个向量，给出每个惩罚矩阵的排名(如果不存在，则以数字方式确定)。L是一个矩阵，它将基础对数平滑参数映射到实际乘以各个二次惩罚的对数平滑参数：如果未提供，则将其视为恒等式。sp是(基础)平滑参数值的向量：正值被视为固定值，负值表示应估计平滑参数。若不提供则视为全部阴性。

paraPen 的一个明显应用是合并随机效应，下面提供了一个示例。在这种情况下，提供的惩罚矩阵将是随机效应的(广义)逆协方差矩阵——即精度矩阵。对应于这些惩罚之一的协方差矩阵的最终估计由惩罚矩阵的(广义)逆乘以估计的尺度参数并除以惩罚的估计的平滑参数给出。例如，如果您使用单位矩阵来惩罚某些被视为 i.i.d 的系数。高斯随机效应，则对于此惩罚，它们的估计方差将是估计的尺度参数除以平滑参数的估计。请参阅下面的‘rail’ 示例。

应谨慎对待惩罚参数项的 P 值。如果您必须拥有它们，则在 anova.gam 和 summary.gam 中使用选项 freq=TRUE ，这往往会为以这种方式实现的随机效应提供合理的结果，但不适用于具有排名低效率惩罚的项(或具有宽eigen-spectrum)。

例子

require(mgcv)
set.seed(10)
## simulate date from y = f(x2)*x1 + error
dat <- gamSim(3,n=400)

b<-gam(y ~ s(x2,by=x1),data=dat)
plot(b,pages=1)
summary(b)

## Factor `by' variable example (with a spurious covariate x0)
## simulate data...

dat <- gamSim(4)

## fit model...
b <- gam(y ~ fac+s(x2,by=fac)+s(x0),data=dat)
plot(b,pages=1)
summary(b)

## note that the preceding fit is the same as....
b1<-gam(y ~ s(x2,by=as.numeric(fac==1))+s(x2,by=as.numeric(fac==2))+
            s(x2,by=as.numeric(fac==3))+s(x0)-1,data=dat)
## ... the `-1' is because the intercept is confounded with the 
## *uncentred* smooths here.
plot(b1,pages=1)
summary(b1)

## repeat forcing all s(x2) terms to have the same smoothing param
## (not a very good idea for these data!)
b2 <- gam(y ~ fac+s(x2,by=fac,id=1)+s(x0),data=dat)
plot(b2,pages=1)
summary(b2)

## now repeat with a single reference level smooth, and 
## two `difference' smooths...
dat$fac <- ordered(dat$fac)
b3 <- gam(y ~ fac+s(x2)+s(x2,by=fac)+s(x0),data=dat,method="REML")
plot(b3,pages=1)
summary(b3)


rm(dat)

## An example of a simple random effects term implemented via 
## penalization of the parametric part of the model...

dat <- gamSim(1,n=400,scale=2) ## simulate 4 term additive truth
## Now add some random effects to the simulation. Response is 
## grouped into one of 20 groups by `fac' and each groups has a
## random effect added....
fac <- as.factor(sample(1:20,400,replace=TRUE))
dat$X <- model.matrix(~fac-1)
b <- rnorm(20)*.5
dat$y <- dat$y + dat$X%*%b

## now fit appropriate random effect model...
PP <- list(X=list(rank=20,diag(20)))
rm <- gam(y~ X+s(x0)+s(x1)+s(x2)+s(x3),data=dat,paraPen=PP)
plot(rm,pages=1)
## Get estimated random effects standard deviation...
sig.b <- sqrt(rm$sig2/rm$sp[1]);sig.b 

## a much simpler approach uses "re" terms...

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

## Simple comparison with lme, using Rail data. 
## See ?random.effects for a simpler method
require(nlme)
b0 <- lme(travel~1,data=Rail,~1|Rail,method="ML") 
Z <- model.matrix(~Rail-1,data=Rail,
     contrasts.arg=list(Rail="contr.treatment"))
b <- gam(travel~Z,data=Rail,paraPen=list(Z=list(diag(6))),method="ML")

b0 
(b$reml.scale/b$sp)^.5 ## `gam' ML estimate of Rail sd
b$reml.scale^.5         ## `gam' ML estimate of residual sd

b0 <- lme(travel~1,data=Rail,~1|Rail,method="REML") 
Z <- model.matrix(~Rail-1,data=Rail,
     contrasts.arg=list(Rail="contr.treatment"))
b <- gam(travel~Z,data=Rail,paraPen=list(Z=list(diag(6))),method="REML")

b0 
(b$reml.scale/b$sp)^.5 ## `gam' REML estimate of Rail sd
b$reml.scale^.5         ## `gam' REML estimate of residual sd

################################################################
## Approximate large dataset logistic regression for rare events
## based on subsampling the zeroes, and adding an offset to
## approximately allow for this.
## Doing the same thing, but upweighting the sampled zeroes
## leads to problems with smoothness selection, and CIs.
################################################################
n <- 50000  ## simulate n data 
dat <- gamSim(1,n=n,dist="binary",scale=.33)
p <- binomial()$linkinv(dat$f-6) ## make 1's rare
dat$y <- rbinom(p,1,p)      ## re-simulate rare response

## Now sample all the 1's but only proportion S of the 0's
S <- 0.02                   ## sampling fraction of zeroes
dat <- dat[dat$y==1 | runif(n) < S,] ## sampling

## Create offset based on total sampling fraction
dat$s <- rep(log(nrow(dat)/n),nrow(dat))

lr.fit <- gam(y~s(x0,bs="cr")+s(x1,bs="cr")+s(x2,bs="cr")+s(x3,bs="cr")+
              offset(s),family=binomial,data=dat,method="REML")

## plot model components with truth overlaid in red
op <- par(mfrow=c(2,2))
fn <- c("f0","f1","f2","f3");xn <- c("x0","x1","x2","x3")
for (k in 1:4) {
       plot(lr.fit,select=k,scale=0)
       ff <- dat[[fn[k]]];xx <- dat[[xn[k]]]
       ind <- sort.int(xx,index.return=TRUE)$ix
       lines(xx[ind],(ff-mean(ff))[ind]*.33,col=2)
}
par(op)
rm(dat)

## A Gamma example, by modify `gamSim' output...
 
dat <- gamSim(1,n=400,dist="normal",scale=1)
dat$f <- dat$f/4 ## true linear predictor 
Ey <- exp(dat$f);scale <- .5 ## mean and GLM scale parameter
## Note that `shape' and `scale' in `rgamma' are almost
## opposite terminology to that used with GLM/GAM...
dat$y <- rgamma(Ey*0,shape=1/scale,scale=Ey*scale)
bg <- gam(y~ s(x0)+ s(x1)+s(x2)+s(x3),family=Gamma(link=log),
          data=dat,method="REML")
plot(bg,pages=1,scheme=1)

作者

Simon N. Wood simon.wood@r-project.org

参考

Wahba (1990) Spline Models of Observational Data SIAM.

Wood S.N. (2017) Generalized Additive Models: An Introduction with R (2nd edition). Chapman and Hall/CRC Press.