R gamm 廣義加性混合模型

說明

通過在正常錯誤恒等鏈接情況下調用 lme 或調用 gammPQL(對 MASS 庫中的 glmmPQL 進行修改)，將指定的廣義加性混合模型 (GAMM) 擬合到數據否則。在後一種情況下，估計值僅約為 MLE。該例程通常比 gam 慢，並且在數值上不太穩健。

要使用 lme4 代替 nlme 作為底層擬合引擎，請參閱包 gamm4 中的 gamm4 。

平滑在對 gam 的調用中指定為固定效應模型公式的一部分，但平滑的擺動分量被視為隨機效應。 lme 可用的隨機效應結構和相關性結構用於指定其他隨機效應和相關性。

假設隨機效應和相關結構主要用於對數據中的殘差相關性進行建模，並且主要興趣在於對包括平滑在內的固定效應模型公式中的項進行推斷。因此，例程計算固定效應公式中所有項(包括平滑項)係數的後驗協方差矩陣。

為了有效地使用此函數，熟悉 gam 和 lme 的使用會有所幫助。

用法

gamm(formula,random=NULL,correlation=NULL,family=gaussian(),
data=list(),weights=NULL,subset=NULL,na.action,knots=NULL,
control=list(niterEM=0,optimMethod="L-BFGS-B",returnObject=TRUE),
niterPQL=20,verbosePQL=TRUE,method="ML",drop.unused.levels=TRUE,
mustart=NULL, etastart=NULL,...)

參數

formula	GAM 公式(另請參見 formula.gam 和 gam.models )。這類似於 glm 的公式，隻不過可以將平滑項( s 、 te 等)添加到公式的右側。請注意，不支持用於平滑和固定平滑參數的id。任何偏移量都應在公式中指定。
random	在調用 lme 時指定的(可選)隨機效應結構：僅允許使用 list 形式，以便於在 gamm 內操作隨機效應結構，以處理平滑項。請參閱下麵的示例。
correlation	可選的 corStruct 對象(請參閱 corClasses )，用於定義 lme 中的相關結構。假定該對象的公式中的任何分組因子都嵌套在任何隨機效應分組因子中，無需在公式中明確說明(這與 lme 的行為略有不同)。這是一般情況下的 GEE 關聯方法。請參閱下麵的示例。
family	family 用於調用 glm 或 gam 。如果沒有偏移項，則帶有恒等鏈接的默認 gaussian 會導致 gamm 通過直接調用 lme 進行擬合，否則使用 gammPQL。
data	包含模型響應變量和公式所需的協變量的 DataFrame 或列表。默認情況下，變量取自 environment(formula) ，通常是調用 gamm 的環境。
weights	在一般情況下，權重與glm權重含義相同。 lme 類型權重參數隻能在恒等鏈接高斯情況下使用，沒有偏移量(有關如何使用此類參數的詳細信息，請參閱 lme 的文檔)。
subset	一個可選向量，指定要在擬合過程中使用的觀測子集。
na.action	一個函數，指示當數據包含“NA”時應該發生什麽。默認值由 ‘options’ 的“na.action”設置設置，如果未設置，則為“na.fail”。 “factory-fresh” 默認為“na.omit”。
knots	這是一個可選列表，包含用戶指定的用於基礎構造的結值。不同的項可以使用不同數量的結，除非它們共享協變量。
control	lme 的擬合控製參數列表，用於替換 lmeControl 返回的默認值。請注意 lme 使用的 EM 迭代次數的設置：平滑是使用自定義 pdMat 類設置的，目前 EM 迭代代碼不支持這些類。如果您提供控製值列表，建議也包含 niterEM=0 ，並且如果您想擾亂模型擬合中使用的起始值(通常是更差的值！)，則僅從 0 開始增加。僅當您的 R 版本沒有 nlminb 優化器函數時，才使用 optimMethod 選項。
niterPQL	PQL 迭代的最大次數(如果有)。
verbosePQL	PQL 是否應該報告其進展情況？
method	當直接調用lme 時，在高斯加性混合模型情況下使用"ML" 或"REML" 中的哪一個。在一般情況下被忽略(或者如果模型有偏移)，在這種情況下使用gammPQL。
drop.unused.levels	默認情況下，未使用的級別會在擬合之前從因子中刪除。對於某些涉及因子變量的平滑，您可能需要將其關閉。僅當您知道自己在做什麽時才這樣做。
mustart	如果使用 PQL，則為平均值的起始值。
etastart	如果使用 PQL，則線性預測器的起始值(如果提供，則覆蓋 mustart)。
...	傳遞的進一步論點，例如至lme

細節

Wahba (1983) 和 Silverman (1985) 引入的樣條平滑貝葉斯模型開啟了估計廣義加性模型中項的平滑程度的可能性，作為被視為隨機效應的平滑項的擺動分量的方差。一些作者已經認識到這一點(參見 Wang 1998；Ruppert、Wand 和 Carroll，2003)，並且在正常誤差中，可以使用通用線性混合效應建模軟件(例如 lme)執行恒等鏈接情況估計。在廣義情況下，迄今為止僅可進行近似推斷，例如使用 Breslow 和 Clayton (1993) 的懲罰 Quasi-Likelihood 方法，如 Venables 和 Ripley (2002) 在 glmmPQL 中實現的方法。這種方法的一個優點是，它允許通過隨機效應或 nlme 庫中可用的相關結構來處理相關錯誤(使用超出嚴格相加情況的相關結構相當於使用 GEE 方法進行擬合)。

有關 GAM 如何表示為混合模型以及如何使用 gamm 中的 lme 或 gammPQL 進行估計的一些細節可以在 Wood (2004 ,2006a,b) 中找到。此外gamm 獲得所有固定效應和平滑項參數的後驗協方差矩陣。該方法類似於 Lin & Zhang (1999) 中說明的方法 - 基於參數的估計，獲得由權重、相關性和隨機效應結構隱含的數據(或廣義情況下的偽數據)的協方差矩陣這些項用於獲得固定和平滑效果的後驗協方差矩陣。

用於表示平滑項的基數與 gam 中使用的基數相同，但自適應平滑基數不可用。使用 predict.gam 從返回的 gam 對象進行預測很簡單，但這會將隨機效應設置為零。如果您想使用設置為預測值的隨機效應進行預測，那麽您可以調整下麵示例中給出的預測代碼。

如果 lme 收斂失敗，請考慮修改 options(mgcv.vc.logrange) ：減少它有助於消除可能性的不確定性(如果這是問題所在)，但減少太大可能會強製過度平滑或欠平滑。有關此選項的更多信息，請參閱notExp2。如果失敗，您可以嘗試增加 control 中的 niterEM 選項：這會擾亂擬合中使用的起始值，但通常是可能性較低的值！請注意，此版本的 gamm 最適合與 R 2.2.0 或更高版本以及 nlme 、3.1-62 及更高版本一起使用，因為它們使用改進的優化器。

值

返回包含兩項的列表：

gam	類 gam 的對象，與 GCV/UBRE 模型選擇相關的信息較少。目前，這包含足夠的信息來使用 predict 、 summary 和 print 方法以及 vis.gam ，但不能使用例如anova 方法函數用於比較模型。這是基於使用 gammPQL 時的工作模型。
lme	lme 或 gammPQL 返回的擬合模型對象。請注意，由於 GAMM 的拆分方式以及對 lme 和 gammPQL 的調用的構造方式，模型公式和分組結構可能看起來相當奇怪。

警告

gamm 的參數列表與 gam 略有不同，gam 參數(例如提供給 gamm 的 gamma 將被忽略)。

gamm 在處理二進製數據時表現不佳，因為它使用 PQL。最好將 gam 與 s(...,bs="re") 術語或 gamm4 一起使用。

gamm 假設您知道自己在做什麽！例如，與 MASS 中的 glmmPQL 不同，它將在 PQL 迭代收斂時從工作模型返回完整的 lme 對象，包括“對數似然”，即使這不是擬合 GAMM 的似然。

如果相關結構用於非常大的數據組，則該例程將非常慢並且占用大量內存。例如絕對不建議嘗試使用數千個數據運行示例部分中的空間示例：通常相關性應該僅適用於可以由分組因子定義的集群內，並且如果這些集群不會變得太大，則擬合通常是合適的可能的。

模型必須包含至少一種隨機效應：具有非零平滑參數的平滑，或參數 random 中指定的隨機效應。

gamm 在數值上不如 gam 穩定：lme 調用偶爾會失敗。請參閱詳細信息部分以獲取建議，或嘗試 ‘gamm4’ 包。

gamm 通常比 gam 慢得多，並且在某些平台上，您可能需要增加 R 可用的內存才能將其用於大型數據集(請參閱 memory.limit )。

請注意，擬合的 GAM 對象中返回的權重是虛擬的，而不是 PQL 迭代使用的權重：這使得部分殘差圖看起來很奇怪。

請注意，返回對象的 gam 對象部分在 gamObject 中定義的所有元素層麵並不完整，並且不繼承自 glm ：因此例如multi-model anova 調用將不起作用。它也是基於使用 PQL 時的工作模型。

gamm 中用於平滑參數的參數化通過有效無窮大和有效零將它們限製在上方和下方。有關如何更改此設置的詳細信息，請參閱notExp2。

不支持鏈接平滑參數和自適應平滑。

例子

library(mgcv)
## simple examples using gamm as alternative to gam
set.seed(0) 
dat <- gamSim(1,n=200,scale=2)
b <- gamm(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat)
plot(b$gam,pages=1)
summary(b$lme) # details of underlying lme fit
summary(b$gam) # gam style summary of fitted model
anova(b$gam) 
gam.check(b$gam) # simple checking plots

b <- gamm(y~te(x0,x1)+s(x2)+s(x3),data=dat) 
op <- par(mfrow=c(2,2))
plot(b$gam)
par(op)
rm(dat)


## Add a factor to the linear predictor, to be modelled as random
dat <- gamSim(6,n=200,scale=.2,dist="poisson")
b2 <- gamm(y~s(x0)+s(x1)+s(x2),family=poisson,
           data=dat,random=list(fac=~1))
plot(b2$gam,pages=1)
fac <- dat$fac
rm(dat)
vis.gam(b2$gam)

## In the generalized case the 'gam' object is based on the working
## model used in the PQL fitting. Residuals for this are not
## that useful on their own as the following illustrates...

gam.check(b2$gam) 

## But more useful residuals are easy to produce on a model
## by model basis. For example...

fv <- exp(fitted(b2$lme)) ## predicted values (including re)
rsd <- (b2$gam$y - fv)/sqrt(fv) ## Pearson residuals (Poisson case)
op <- par(mfrow=c(1,2))
qqnorm(rsd);plot(fv^.5,rsd)
par(op)

## now an example with autocorrelated errors....
n <- 200;sig <- 2
x <- 0:(n-1)/(n-1)
f <- 0.2*x^11*(10*(1-x))^6+10*(10*x)^3*(1-x)^10
e <- rnorm(n,0,sig)
for (i in 2:n) e[i] <- 0.6*e[i-1] + e[i]
y <- f + e
op <- par(mfrow=c(2,2))
## Fit model with AR1 residuals
b <- gamm(y~s(x,k=20),correlation=corAR1())
plot(b$gam);lines(x,f-mean(f),col=2)
## Raw residuals still show correlation, of course...
acf(residuals(b$gam),main="raw residual ACF")
## But standardized are now fine...
acf(residuals(b$lme,type="normalized"),main="standardized residual ACF")
## compare with model without AR component...
b <- gam(y~s(x,k=20))
plot(b);lines(x,f-mean(f),col=2)

## more complicated autocorrelation example - AR errors
## only within groups defined by `fac'
e <- rnorm(n,0,sig)
for (i in 2:n) e[i] <- 0.6*e[i-1]*(fac[i-1]==fac[i]) + e[i]
y <- f + e
b <- gamm(y~s(x,k=20),correlation=corAR1(form=~1|fac))
plot(b$gam);lines(x,f-mean(f),col=2)
par(op) 

## more complex situation with nested random effects and within
## group correlation 

set.seed(0)
n.g <- 10
n<-n.g*10*4
## simulate smooth part...
dat <- gamSim(1,n=n,scale=2)
f <- dat$f
## simulate nested random effects....
fa <- as.factor(rep(1:10,rep(4*n.g,10)))
ra <- rep(rnorm(10),rep(4*n.g,10))
fb <- as.factor(rep(rep(1:4,rep(n.g,4)),10))
rb <- rep(rnorm(4),rep(n.g,4))
for (i in 1:9) rb <- c(rb,rep(rnorm(4),rep(n.g,4)))
## simulate auto-correlated errors within groups
e<-array(0,0)
for (i in 1:40) {
  eg <- rnorm(n.g, 0, sig)
  for (j in 2:n.g) eg[j] <- eg[j-1]*0.6+ eg[j]
  e<-c(e,eg)
}
dat$y <- f + ra + rb + e
dat$fa <- fa;dat$fb <- fb
## fit model .... 
b <- gamm(y~s(x0,bs="cr")+s(x1,bs="cr")+s(x2,bs="cr")+
  s(x3,bs="cr"),data=dat,random=list(fa=~1,fb=~1),
  correlation=corAR1())
plot(b$gam,pages=1)
summary(b$gam)
vis.gam(b$gam)

## Prediction from gam object, optionally adding 
## in random effects. 

## Extract random effects and make names more convenient...
refa <- ranef(b$lme,level=5)
rownames(refa) <- substr(rownames(refa),start=9,stop=20)
refb <- ranef(b$lme,level=6)
rownames(refb) <- substr(rownames(refb),start=9,stop=20)

## make a prediction, with random effects zero...
p0 <- predict(b$gam,data.frame(x0=.3,x1=.6,x2=.98,x3=.77))

## add in effect for fa = "2" and fb="2/4"...
p <- p0 + refa["2",1] + refb["2/4",1] 

## and a "spatial" example...
library(nlme);set.seed(1);n <- 100
dat <- gamSim(2,n=n,scale=0) ## standard example
attach(dat)
old.par<-par(mfrow=c(2,2))
contour(truth$x,truth$z,truth$f)  ## true function
f <- data$f                       ## true expected response
## Now simulate correlated errors...
cstr <- corGaus(.1,form = ~x+z)  
cstr <- Initialize(cstr,data.frame(x=data$x,z=data$z))
V <- corMatrix(cstr) ## correlation matrix for data
Cv <- chol(V)
e <- t(Cv) %*% rnorm(n)*0.05 # correlated errors
## next add correlated simulated errors to expected values
data$y <- f + e ## ... to produce response
b<- gamm(y~s(x,z,k=50),correlation=corGaus(.1,form=~x+z),
         data=data)
plot(b$gam) # gamm fit accounting for correlation
# overfits when correlation ignored.....  
b1 <- gamm(y~s(x,z,k=50),data=data);plot(b1$gam) 
b2 <- gam(y~s(x,z,k=50),data=data);plot(b2)
par(old.par)

作者

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

參考

Breslow, N. E. and Clayton, D. G. (1993) Approximate inference in generalized linear mixed models. Journal of the American Statistical Association 88, 9-25.

Lin, X and Zhang, D. (1999) Inference in generalized additive mixed models by using smoothing splines. JRSSB. 55(2):381-400

Pinheiro J.C. and Bates, D.M. (2000) Mixed effects Models in S and S-PLUS. Springer

Ruppert, D., Wand, M.P. and Carroll, R.J. (2003) Semiparametric Regression. Cambridge

Silverman, B.W. (1985) Some aspects of the spline smoothing approach to nonparametric regression. JRSSB 47:1-52

Venables, W. N. and Ripley, B. D. (2002) Modern Applied Statistics with S. Fourth edition. Springer.

Wahba, G. (1983) Bayesian confidence intervals for the cross validated smoothing spline. JRSSB 45:133-150

Wood, S.N. (2004) Stable and efficient multiple smoothing parameter estimation for generalized additive models. Journal of the American Statistical Association. 99:673-686

Wood, S.N. (2003) Thin plate regression splines. J.R.Statist.Soc.B 65(1):95-114

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

Wood S.N. (2006b) Generalized Additive Models: An Introduction with R. Chapman and Hall/CRC Press.

Wang, Y. (1998) Mixed effects smoothing spline analysis of variance. J.R. Statist. Soc. B 60, 159-174

https://www.maths.ed.ac.uk/~swood34/

也可以看看

magic 用於相關數據的替代方案，te , s , predict.gam , plot.gam , summary.gam , negbin , vis.gam , pdTens , gamm4 ( https://cran.r-project.org/package=gamm4 )