R predict.gam 根據擬合的 GAM 模型進行預測

說明

采用 gam() 生成的擬合 gam 對象，並在給定模型協變量的一組新值或用於模型擬合的原始值的情況下生成預測。基於模型係數的後驗分布，預測可能伴隨有標準誤差。該例程可以選擇返回矩陣，模型係數必須預先乘以該矩陣，以便在提供的協變量值處產生線性預測變量的值：這對於獲得從模型導出的數量的可信區域(例如，平滑)，以及 R 之外的查找表預測(請參見下麵的示例代碼)。

用法

## S3 method for class 'gam'
predict(object,newdata,type="link",se.fit=FALSE,terms=NULL,
        exclude=NULL,block.size=NULL,newdata.guaranteed=FALSE,
        na.action=na.pass,unconditional=FALSE,iterms.type=NULL,...)

參數

object	由 gam() 生成的擬合 gam 對象。
newdata	包含需要預測的模型協變量值的 DataFrame 或列表。如果未提供，則返回與原始數據相對應的預測。如果提供了newdata，那麽它應該包含預測所需的所有變量：如果沒有，則會生成警告。請參閱與 link{linear.functional.terms} 一起使用的詳細信息。
type	當其值為"link"(默認)時，返回線性預測器(可能具有相關的標準誤差)。當type="terms" 線性預測器的每個分量單獨返回(可能帶有標準誤差)：這包括參數模型分量，後麵是每個平滑分量，但不包括任何偏移和任何截距。 type="iterms" 是相同的，隻是為平滑分量返回的任何標準誤差將包括截距/總體平均值的不確定性。當返回 type="response" 對響應規模的預測時(可能具有近似標準誤差)。當 type="lpmatrix" 時，返回一個矩陣，該矩陣在後乘於參數向量時產生線性預測器的值(減去任何偏移)(在這種情況下 se.fit 被忽略)。後一個選項對於獲取從模型導出的數量的方差估計最有用：例如積分數量或平滑的導數。線性預測矩陣也可用於在 R 之外實現近似預測(請參見下麵的示例代碼)。
se.fit	當此值為 TRUE(非默認值)時，將為每個預測返回標準誤差估計值。
terms	如果type=="terms" 或type="iterms" 則僅返回此數組中指定的項(平滑或參數)的結果。否則，該數組中未命名的任何術語都將設置為零。如果NULL則包含所有術語。
exclude	如果 type=="terms" 或 type="iterms" 則不會返回此數組中指定的項(平滑或參數)。否則，該數組中命名的任何術語都將設置為零。如果NULL，則不排除任何術語。請注意，這是模型摘要中顯示的術語名稱，請參閱示例。您可以通過設置 newdata.guaranteed=TRUE 來避免為排除的平滑項提供協變量，這將避免對 newdata 進行所有檢查(不能跳過參數項的協變量)。
block.size	每次調用底層代碼時要處理的最大預測數：越大速度越快，但內存消耗更大。設置為 < 1 以使用預測總數。如果NULL，則如果提供了新數據，則塊大小為 1000，否則為模型框架中的行數。
newdata.guaranteed	設置為 TRUE 以關閉對 newdata 的所有檢查，除了因子水平的健全性之外：這可以加快大型預測任務的速度，但 newdata 必須完整，並且不需要預測變量的 NA 值。模型。
na.action	如何處理 newdata 中的 NA 值。使用默認值 na.pass ，newdata 中包含所需預測變量的 NA 值的任何行都會產生 NA 預測(即使相關術語沒有 NA 預測變量)。 na.exclude 或 na.omit 會導致刪除 newdata 行(如果它們包含所需預測變量的任何 NA 值)。如果 newdata 缺失，則 NA 處理由 object$na.action 確定。
unconditional	如果TRUE，則使用平滑參數不確定性校正協方差矩陣(如果可用)，否則使用以估計平滑參數為條件的協方差矩陣。
iterms.type	如果 type="iterms" 則標準誤差可以包括總體平均值的不確定性(默認、包括固定效應和隨機效應)或僅包含非平滑固定效應平均值的不確定性(iterms.type=2)。
...	其他論點。

細節

predict.gam 產生的標準誤差基於擬合 gam 對象中參數 Vp 的貝葉斯後驗協方差矩陣。

當使用 linear.functional.terms 模型進行預測時，有兩種可能性。如果要在預測中使用求和約定，就像在擬合中一樣，那麽 newdata 應該是一個列表，其中命名的矩陣參數對應於擬合中作為矩陣的任何變量。或者，人們可以選擇簡單地評估特定值的構成平滑，在這種情況下，矩陣參數可以用向量替換(並且 newdata 可以是數據幀)。有關示例代碼，請參閱linear.functional.terms。

為了便於使用 termplot 進行繪圖，如果 object 擁有屬性 "para.only" 和 type=="terms"，則僅返回 1 階參數項(即 termplot 可以處理的參數項)。

請注意，與其他預測函數一樣，與 gam 模型公式中指定的偏移量不同，在預測時始終忽略作為參數提供給 gam 的任何偏移量。

請參閱示例，了解如何使用 lpmatrix 獲取從模型導出的數量的可信區域。

值

如果type=="lpmatrix"，則返回一個矩陣，當應用於模型係數向量時，該矩陣將在所提供的協變量值處給出線性預測變量值(減去任何偏移)的向量。否則，如果 se.fit 是 TRUE，則返回 2 項列表，其中項(兩個數組)fit 和 se.fit 包含預測和關聯的標準誤差估計，否則返回預測數組。返回數組的維度取決於 type 是否是 "terms"：如果是，則數組是二維的，線性預測器中的每個項都是單獨的，否則數組是一維的，包含線性預測器/預測值(或相應的 s.e.s)。按項返回的線性預測器將不包括偏移量或截距。

newdata 可以是 DataFrame 架、列表或 model.frame：如果它是模型框架，則必須提供所有變量。

警告

如果在平滑調用中使用協變量的數據相關變換，則預測可能不正確。請參閱示例。

請注意，此函數的行為與 Splus 中的predict.gam() 不同。

type=="terms" 與 predict.lm 對參數模型組件所做的並不完全匹配。

例子

library(mgcv)
n <- 200
sig <- 2
dat <- gamSim(1,n=n,scale=sig)

b <- gam(y~s(x0)+s(I(x1^2))+s(x2)+offset(x3),data=dat)

newd <- data.frame(x0=(0:30)/30,x1=(0:30)/30,x2=(0:30)/30,x3=(0:30)/30)
pred <- predict.gam(b,newd)
pred0 <- predict(b,newd,exclude="s(x0)") ## prediction excluding a term
## ...and the same, but without needing to provide x0 prediction data...
newd1 <- newd;newd1$x0 <- NULL ## remove x0 from `newd1'
pred1 <- predict(b,newd1,exclude="s(x0)",newdata.guaranteed=TRUE)

## custom perspective plot...

m1 <- 20;m2 <- 30; n <- m1*m2
x1 <- seq(.2,.8,length=m1);x2 <- seq(.2,.8,length=m2) ## marginal grid points
df <- data.frame(x0=rep(.5,n),x1=rep(x1,m2),x2=rep(x2,each=m1),x3=rep(0,n))
pf <- predict(b,newdata=df,type="terms")
persp(x1,x2,matrix(pf[,2]+pf[,3],m1,m2),theta=-130,col="blue",zlab="")

#############################################
## difference between "terms" and "iterms"
#############################################
nd2 <- data.frame(x0=c(.25,.5),x1=c(.25,.5),x2=c(.25,.5),x3=c(.25,.5))
predict(b,nd2,type="terms",se=TRUE)
predict(b,nd2,type="iterms",se=TRUE)

#########################################################
## now get variance of sum of predictions using lpmatrix
#########################################################

Xp <- predict(b,newd,type="lpmatrix") 

## Xp %*% coef(b) yields vector of predictions

a <- rep(1,31)
Xs <- t(a) %*% Xp ## Xs %*% coef(b) gives sum of predictions
var.sum <- Xs %*% b$Vp %*% t(Xs)


#############################################################
## Now get the variance of non-linear function of predictions
## by simulation from posterior distribution of the params
#############################################################

rmvn <- function(n,mu,sig) { ## MVN random deviates
  L <- mroot(sig);m <- ncol(L);
  t(mu + L%*%matrix(rnorm(m*n),m,n)) 
}

br <- rmvn(1000,coef(b),b$Vp) ## 1000 replicate param. vectors
res <- rep(0,1000)
for (i in 1:1000)
{ pr <- Xp %*% br[i,] ## replicate predictions
  res[i] <- sum(log(abs(pr))) ## example non-linear function
}
mean(res);var(res)

## loop is replace-able by following .... 

res <- colSums(log(abs(Xp %*% t(br))))



##################################################################
## The following shows how to use use an "lpmatrix" as a lookup 
## table for approximate prediction. The idea is to create 
## approximate prediction matrix rows by appropriate linear 
## interpolation of an existing prediction matrix. The additivity 
## of a GAM makes this possible. 
## There is no reason to ever do this in R, but the following 
## code provides a useful template for predicting from a fitted 
## gam *outside* R: all that is needed is the coefficient vector 
## and the prediction matrix. Use larger `Xp'/ smaller `dx' and/or 
## higher order interpolation for higher accuracy.  
###################################################################

xn <- c(.341,.122,.476,.981) ## want prediction at these values
x0 <- 1         ## intercept column
dx <- 1/30      ## covariate spacing in `newd'
for (j in 0:2) { ## loop through smooth terms
  cols <- 1+j*9 +1:9      ## relevant cols of Xp
  i <- floor(xn[j+1]*30)  ## find relevant rows of Xp
  w1 <- (xn[j+1]-i*dx)/dx ## interpolation weights
  ## find approx. predict matrix row portion, by interpolation
  x0 <- c(x0,Xp[i+2,cols]*w1 + Xp[i+1,cols]*(1-w1))
}
dim(x0)<-c(1,28) 
fv <- x0%*%coef(b) + xn[4];fv    ## evaluate and add offset
se <- sqrt(x0%*%b$Vp%*%t(x0));se ## get standard error
## compare to normal prediction
predict(b,newdata=data.frame(x0=xn[1],x1=xn[2],
        x2=xn[3],x3=xn[4]),se=TRUE)

##################################################################
# illustration of unsafe scale dependent transforms in smooths....
##################################################################

b0 <- gam(y~s(x0)+s(x1)+s(x2)+x3,data=dat) ## safe
b1 <- gam(y~s(x0)+s(I(x1/2))+s(x2)+scale(x3),data=dat) ## safe
b2 <- gam(y~s(x0)+s(scale(x1))+s(x2)+scale(x3),data=dat) ## unsafe
pd <- dat; pd$x1 <- pd$x1/2; pd$x3 <- pd$x3/2
par(mfrow=c(1,2))
plot(predict(b0,pd),predict(b1,pd),main="b0 and b1 predictions match")
abline(0,1,col=2)
plot(predict(b0,pd),predict(b2,pd),main="b2 unsafe, doesn't match")
abline(0,1,col=2)


####################################################################
## Differentiating the smooths in a model (with CIs for derivatives)
####################################################################

## simulate data and fit model...
dat <- gamSim(1,n=300,scale=sig)
b<-gam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat)
plot(b,pages=1)

## now evaluate derivatives of smooths with associated standard 
## errors, by finite differencing...
x.mesh <- seq(0,1,length=200) ## where to evaluate derivatives
newd <- data.frame(x0 = x.mesh,x1 = x.mesh, x2=x.mesh,x3=x.mesh)
X0 <- predict(b,newd,type="lpmatrix") 

eps <- 1e-7 ## finite difference interval
x.mesh <- x.mesh + eps ## shift the evaluation mesh
newd <- data.frame(x0 = x.mesh,x1 = x.mesh, x2=x.mesh,x3=x.mesh)
X1 <- predict(b,newd,type="lpmatrix")

Xp <- (X1-X0)/eps ## maps coefficients to (fd approx.) derivatives
colnames(Xp)      ## can check which cols relate to which smooth

par(mfrow=c(2,2))
for (i in 1:4) {  ## plot derivatives and corresponding CIs
  Xi <- Xp*0 
  Xi[,(i-1)*9+1:9+1] <- Xp[,(i-1)*9+1:9+1] ## Xi%*%coef(b) = smooth deriv i
  df <- Xi%*%coef(b)              ## ith smooth derivative 
  df.sd <- rowSums(Xi%*%b$Vp*Xi)^.5 ## cheap diag(Xi%*%b$Vp%*%t(Xi))^.5
  plot(x.mesh,df,type="l",ylim=range(c(df+2*df.sd,df-2*df.sd)))
  lines(x.mesh,df+2*df.sd,lty=2);lines(x.mesh,df-2*df.sd,lty=2)
}

作者

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

The design is inspired by the S function of the same name described in Chambers and Hastie (1993) (but is not a clone).

參考

Chambers and Hastie (1993) Statistical Models in S. Chapman & Hall.

Marra, G and S.N. Wood (2012) Coverage Properties of Confidence Intervals for Generalized Additive Model Components. Scandinavian Journal of Statistics, 39(1), 53-74. doi:10.1111/j.1467-9469.2011.00760.x

Wood S.N. (2017, 2nd ed) Generalized Additive Models: An Introduction with R. Chapman and Hall/CRC Press. doi:10.1201/9781315370279

也可以看看

gam , gamm , plot.gam