# [R] Using lme() inside a function

Date: Fri, 09 May 2008 15:23:43 +0200

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

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

 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:

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