[R] nlme splines model.frame.default error

Jacob Wegelin jacobwegelin at fastmail.fm
Tue Jun 2 23:56:42 CEST 2015


On 2015-06-02 Tue 14:20, Jacob Wegelin wrote:
> I want to use a specialized function to compute knots on the fly to pass to 
> splines::bs within a call to nlme::lme.
>
> The specialized function will compute knots from the x-axis variable (the x 
> argument in a call to splines::bs).
>
> The syntax works with lm. But when I try it with lme, the following error is 
> returned:
>
> Error in model.frame.default(formula = ~why + eks + toydat + ID, data = list( 
> :
>  invalid type (list) for variable 'toydat'
>
> Below is a simplified example showing the "desired syntax," which returns the 
> error, and a workaround. I am mystified what inherently is wrong with the 
> "desired syntax", since it works fine with stats::lm.
>
> set.seed(5)
> library(splines)
> library(nlme)
> NID<-5
> toydat<-data.frame(eks=0:0:((2*NID)-1))
> toydat$ID<-factor(rep(LETTERS[1:NID], each=2))
> toydat
> toydat$why<-with(toydat, as.integer(100*sin((eks / 7 * pi)^2) + 
> rnorm(eks)/10))
> customKnotsFn<-function(some)3.5
> print(toydat)
> print(summary(toydat))
> print(customKnotsFn)
>
> # lm has no trouble:
>
> mylm<-lm(why~bs(eks,knots=customKnotsFn(some=toydat$eks)), data=toydat)
> print(mylm$call)
>
> #  The "desired syntax" returns an error:
>
> lme(fixed=why~bs(eks,knots=customKnotsFn(some=toydat$eks)),random=~1|ID, 
> data=toydat)
>
> #  Removing the argument from customKnotsFn eliminates the error, but then 
> customKnotsFn is useless for practical purposes:
>
> lme(fixed=why~bs(eks,knots=customKnotsFn()),random=~1|ID, data=toydat)
>
> #  Workaround returns no error:
>
> mylme<- with(toydat, 
> lme(fixed=why~bs(eks,knots=customKnotsFn(some=eks)),random=~1|ID))
> print(mylme$call)

# Here is a perhaps slicker workaround:

mybs <-function(x, df) {
      myknots<-customKnotsFn(x)
      bs(x, knots=myknots)
}

lme(fixed=why~mybs(eks, df=2), random=~1|ID, data=toydat)

# But why cannot lme handle the "desired syntax"?



More information about the R-help mailing list