CODEWITHSUNDEEP
rfunction
user3277394
user3277394

Reputation: 41

Using update() function in conjunction with gam() in new mgcv package

In order to test a series of gam models for capture counts for each of seven sex/life-stages of a study species, I used a function that employed update() to iteratively substitute into the base model a sequence of predictors and produce an AIC score for each. There appears to be a problem implementing the same code with a gam model in the new mgcv package. Here is a functional, simplified subset of data to proceed.

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))

Before I updated mgcv to 1.8-9, the function ran successfully in this form (note, I edited henceforth to directly reference the variables, as opposed to attaching repex):

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

I would then produce a list of AIC scores with something like this:

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

After I updated R to 3.2.3 and mgcv to 1.8-9, I ran the above scripts as well as simply testing the function with:

MAIC(repex$Solar)

and receive this message:

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

I have been futzing around with this and have determined that, contrary to some suggestions, there isn't anything essentially problematic with the line of code m<-update(m,.~.+x). I simplified the above MAIC function to try to locate the source of trouble, and ran it successfully with the above repex data subset with a glm() and lm() call:

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)

But when I change the model to a gam, I receive the above error:

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

This happens no matter how the base gam is constructed, unless the update() is omitted as follows:

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

My sessionInfo() call brings:

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

My ls() call brings:

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

My find("x") call brings:

character(0)

Finally, I cross-checked the AIC scores of MAIC1(repex$Solar), MAIC2(repex$Solar), and MAIC4(repex$Solar), respectively, with:

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))

Hopefully this helps clear things up.

Upvotes: 2

Views: 1204

Answers (1)

Ben Bolker
Ben Bolker

Reputation: 226027

reformulate() is a sensible way to do this:

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")

This makes it easier, e.g., to operate over a vector of column names.

If you really want to be able to use "raw" variable names (e.g. Solar rather than "Solar", you can use 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)

(alternatively, you could just write MAIC2 <- function(x) MAIC(deparse(substitute(x))) ...)

If you really want to use update with raw variables it takes a little more magic ...

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

Upvotes: 2

Related Questions