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.