[R] Using lme() inside a function

[Jorunn Slagstad]

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:

 library(MASS)
 library(nlme)
 PredRes<-function(D1)
 {
 lmemod=lme(distance~age*Sex, random=~1|Subject,
data=subset(D1,age!=14),method="ML")
 themod=stepAIC(lmemod,dir="both")
 prs=predict(themod,newdata=subset(D1,age==14))
 obs<-subset(D1,age==14)$distance
 print(mean(obs-prs))
 }

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
 prs=predict(themod,newdata=subset(Orthodont,age==14))
 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:

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)

Locale:
LC_CTYPE=en_US.UTF-8;LC_NUMERIC=C;LC_TIME=en_US.UTF-8;LC_COLLATE=en_US.UTF-8;LC_MONETARY=en_US.UTF-8;LC_MESSAGES=en_US.UTF-8;LC_PAPER=en_US.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=en_US.UTF-8;LC_IDENTIFICATION=C

Search Path:
 .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.

--
Regards
Jorunn Slagstad