Re: [R] Using lme() inside a function

From: ONKELINX, Thierry <>
Date: Fri, 09 May 2008 15:45:55 +0200


This solutions works with R 2.7.0 under windows

PredRes <- function(cal, val){

    cal <<- cal
    lmemod <- lme(distance ~ age * Sex, random = ~1|Subject, data = cal, method="ML")

    themod <- stepAIC(lmemod, dir="both", trace = FALSE)     prs <- predict(themod, newdata = val)     obs <- val$distance
    print(mean(obs - prs))

PredRes(cal = subset(Orthodont, age!=14), val = subset(Orthodont, age==14))

HTH, Thierry

ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature and Forest
Cel biometrie, methodologie en kwaliteitszorg / Section biometrics, methodology and quality assurance
Gaverstraat 4
9500 Geraardsbergen
tel. + 32 54/436 185

To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of.
~ Sir Ronald Aylmer Fisher

The plural of anecdote is not data.
~ Roger Brinner

The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey

-----Oorspronkelijk bericht-----
Van: [] Namens Jorunn Slagstad
Verzonden: vrijdag 9 mei 2008 15:24
Aan: R-help
Onderwerp: [R] Using lme() inside a function

Dear R-help

I'm working on a large dataset which I have divided into 20 subsets based on similar features. Each subset consists of observations from different locations and I wish to use the location as a random effect. For each group I want to select regressors by a stepwise procedure and include a random effect on the intercept. I use stepAIC() and lme(). (The lmer()-function doesn't work with the stepAIC()-function.)

Since I have many groups, and I wish to do the same thing for each group, I have constructed a function which takes the dataset as input variable and gives a prediction result (here mean absolute error) as output.

This is an example using the Orthodont dataset:

 lmemod=lme(distance~age*Sex, random=~1|Subject, data=subset(D1,age!=14),method="ML")
 prs=predict(themod,newdata=subset(D1,age==14))  obs<-subset(D1,age==14)$distance

Using this function with D1=Orthodont gives:

> PredRes(Orthodont)
 Start: AIC=345.12
 distance ~ age * Sex
 Error in subset(D1, age != 14) : object "D1" not found

The code works when I take Orthodont in directly:

 lmemod=lme(distance~age*Sex, random=~1|Subject, data=subset(Orthodont,age!=14),method="ML")  themod=stepAIC(lmemod,dir="both")
Start: AIC=345.12
distance ~ age * Sex

          Df AIC
- age:Sex 1 344.49
<none> 345.12

Step: AIC=344.49
distance ~ age + Sex

          Df    AIC
<none>       344.49
+ age:Sex  1 345.12
- Sex      1 348.70
- age      1 371.77
 obs<-subset(Orthodont,age==14)$distance > print(mean(obs-prs))
[1] 0.2962963

How can I make this code work with dataset as input variable in a function?

I'm using R version:

 platform = x86_64-unknown-linux-gnu
 arch = x86_64
 os = linux-gnu
 system = x86_64, linux-gnu
 status =
 major = 2
 minor = 6.0
 year = 2007
 month = 10
 day = 03
 svn rev = 43063
 language = R
 version.string = R version 2.6.0 (2007-10-03)

 .GlobalEnv, package:MASS, package:nlme, package:stats, package:graphics, package:grDevices, package:utils, package:datasets, package:methods, Autoloads, package:base

By the way the R version 2.7.0 Patched (2008-05-08 r45647) gives the same error message.

Jorunn Slagstad

______________________________________________ mailing list
PLEASE do read the posting guide
and provide commented, minimal, self-contained, reproducible code.

______________________________________________ mailing list
PLEASE do read the posting guide
and provide commented, minimal, self-contained, reproducible code.
Received on Fri 09 May 2008 - 14:00:13 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Fri 09 May 2008 - 14:30:36 GMT.

Mailing list information is available at Please read the posting guide before posting to the list.

list of date sections of archive