Marc in the box
Marc in the box

Reputation: 11995

gls object produced with wrapper function fails for prediction with new data

In the following example, I'm using a wrapper function to fit a gls object. I can successfully return the model prediction for the original data, but not when using newdata, which returns the error "Error in eval(expr, envir, enclos) : 'nthcdr' needs a list to CDR down":

# library -----------------------------------------------------------------
library(nlme)


# wrapper function --------------------------------------------------------
test.gls <- function(data, ...){
  fit <- gls(data=data, ...)
  return(fit)
}


# make data ---------------------------------------------------------------
set.seed(1)
n <- 20
x <- sort(rlnorm(n, meanlog = 0.25, sdlog = 1.5))
b <- 100
cv <- 0.5
y <- x*b * rlnorm(n, 0, cv)
dat <- data.frame(x,y)


# fit model ---------------------------------------------------------------
fit <- test.gls(data=dat, model=y~x-1, weights=varExp())
class(fit)
# [1] "gls"
plot(y~x, dat)
lines(dat$x, predict(fit)) # works


# prediction --------------------------------------------------------------
newdat <- data.frame(x=sort(rlnorm(n, meanlog = 0.25, sdlog = 1.5)))
pred <- predict(fit)
newpred <- predict(fit, newdata = newdat)
# Error in eval(expr, envir, enclos) : 'nthcdr' needs a list to CDR down

Upvotes: 1

Views: 366

Answers (1)

Marc in the box
Marc in the box

Reputation: 11995

The following answer seems to provide one solution using do.call: https://stackoverflow.com/a/7668846/1199289

Example:

# library -----------------------------------------------------------------
library(nlme)


# wrapper function --------------------------------------------------------
test.gls <- function(data, ...){
  fit <- gls(data=data, ...)
  return(fit)
}
test.gls2 <- function(argList=NULL){
  fit <- do.call("gls", args = argList)
  return(fit)
}


# make data ---------------------------------------------------------------
set.seed(1)
n <- 20
x <- sort(rlnorm(n, meanlog = 0.25, sdlog = 1.5))
b <- 100
cv <- 0.5
y <- x*b * rlnorm(n, 0, cv)
dat <- data.frame(x,y)


# fit model ---------------------------------------------------------------
fit <- test.gls(data=dat, model=y~x-1, weights=varExp())
fit2 <- test.gls2(argList=list(data=dat, model=y~x-1, weights=varExp()))
class(fit)
# [1] "gls"
plot(y~x, dat)
lines(dat$x, predict(fit)) # works


# prediction --------------------------------------------------------------
newdat <- data.frame(x=sort(rlnorm(n, meanlog = 0.25, sdlog = 1.5)))
pred <- predict(fit)
newpred <- predict(fit, newdata = newdat)
# Error in eval(expr, envir, enclos) : 'nthcdr' needs a list to CDR down
newpred2 <- predict(fit2, newdata = newdat) # works

Upvotes: 3

Related Questions