在新mgcv包中使用update（）函數與gam（）聯合使用

Question

爲了測試一系列gam模型的捕獲計數爲研究物種的七個性/生命階段，我使用的函數僱用update()迭代地將一系列預測變量代入基礎模型中，併爲每個預測變量產生一個AIC得分。在新的mgcv包中使用gam模型實現相同的代碼似乎存在問題。這是一個功能簡化的數據子集。

repex=structure(list(Day = c(183L, 190L, 197L, 204L, 211L, 218L, 225L, 
232L, 239L, 246L, 175L, 182L), M = c(18L, 43L, 22L, 20L, 
7L, 1L, 1L, 0L, 0L, 0L, 0L, 17L), Solar = c(77L, 59L, 20L, 
55L, -3L, -44L, 13L, 58L, 8L, 6L, -28L, 12L)), .Names = c("Day", "M", "Solar"), 
class = "data.frame", row.names = c(NA, -12L)) 

我更新mgcv到1.8-9之前，該功能以這種形式成功地跑了（注意，我編輯從今以後直接引用變量，而不是附加repex）：

MAIC<-function(x){ 
    m<-gam(repex$M~s(repex$Day),data=repex,family=poisson) 
    m<-update(m,.~.+x) 
    return(AIC(m)) 
} 

我然後將產生AIC評分列表像這樣的東西：

lapply(c('Solar'),function(x) MAIC(repex[ , x])) 

我更新R鍵3.2.3和mgcv到1.8-9後，我跑了上述腳本以及簡單地測試功能：

MAIC(repex$Solar) 

，並收到此消息：

Error in eval(expr, envir, enclos) : object 'x' not found 

我一直在把玩這個四周，已經確定的是，違背了一些建議，沒有代碼行m<-update(m,.~.+x)基本上有問題。我簡化了上述MAIC功能，以試圖找到故障源，並與一個glm()和lm()呼叫上述repex數據子集成功地運行它：

MAIC1 <- function(x){ 
    m <- glm(repex$M~repex$Day,data=repex,family=poisson) 
    m <- update(m,.~.+x) 
    return(AIC(m)) 
} 
MAIC1(repex$Solar) 

MAIC2 <- function(x){ 
    m <- lm(repex$M~repex$Day,data=repex) 
    m <- update(m,.~.+x) 
    return(AIC(m)) 
} 
MAIC2(repex$Solar) 

但是，當我的模式更改爲GAM，我收到上述錯誤：

MAIC3 <- function(x){ 
    m <- gam(repex$M~s(repex$Day),data=repex) 
    m <- update(m,.~.+x) 
    return(AIC(m)) 
} 
MAIC(repex$Solar) 

這發生不管基部GAM是如何構造的，除非省略update()如下：

MAIC4<-function(x){ 
    m<-gam(repex$M~s(repex$Day)+x,data=repex) 
    return(AIC(m)) 
} 
MAIC4(repex$Solar) 

我sessionInfo()召喚着：

R version 3.2.3 (2015-12-10) 
Platform: x86_64-w64-mingw32/x64 (64-bit) 
Running under: Windows 7 x64 (build 7601) Service Pack 1 

我ls()召喚着：

[1] "MAIC" "MAIC1" "MAIC2" "MAIC3" "MAIC4" "repex" 

我find("x")召喚着：

character(0) 

最後，我反覆覈對的MAIC1(repex$Solar)的AIC分數，MAIC2(repex$Solar)和MAIC4(repex$Solar)，分別爲：

AIC(glm(repex$M~repex$Day+repex$Solar,data=repex,family=poisson)) 
AIC(lm(repex$M~repex$Day+repex$Solar,data=repex)) 
AIC(gam(repex$M~s(repex$Day)+repex$Solar,data=repex)) 

希望這有助於清理一些事情。

Answer 1

reformulate()就是做這個明智的方法：

repex <- data.frame(Day = c(183, 190, 197, 204, 211, 
218, 225, 232, 239, 246, 175, 182), 
    M = c(18, 43, 22, 20, 7, 1, 1, 0, 0, 0, 0, 17), 
    Solar = c(77, 59, 20, 55, -3, -44, 13, 58, 
    8, 6, -28, 12)) 

library(mgcv) 
MAIC <- function(x) { 
    form <- reformulate(c("s(Day)",x),response="M") 
    m <- gam(form, data=repex,family=poisson) 
    return(AIC(m)) 
} 
MAIC("Solar") 

這使得它更容易，例如，對列名的載體上操作。

如果你真的希望能夠用「原始」變量名（例如Solar而不是"Solar"，你可以使用deparse(substitute())

MAIC2 <- function(x) { 
    xx <- deparse(substitute(x)) 
    form <- reformulate(c("s(Day)",xx),response="M") 
    m <- gam(form, data=repex,family=poisson) 
    return(AIC(m)) 
} 
MAIC2(Solar) 

（或者，你可以只寫MAIC2 <- function(x) MAIC(deparse(substitute(x))) ...）

如果你真的想使用update用生變量，它需要多一點的魔術......

MAIC3 <- function(x) { 
    m <- gam(M~s(Day),data=repex,family=poisson) 
    m2 <- update(m,bquote(.~.+.(substitute(x)))) 
    return(AIC(m2)) 
} 
MAIC3(Solar)