[R] How can one make stepAIC and lme

[Prof Brian Ripley]

You are using subset() as a function -- lme has a subset argument, and 
subset() has non-standard semantics.

The dataset you are using is subset(D1,age!=14): it is that which you need 
to make sure is visible.  You didn't make D1 visible, just its columns.

On Tue, 13 May 2008, Jorunn Slagstad wrote:

> Hi again,
>
> I've tried using with() inside my function in the following manner:
>
> library(nlme)
> library(MASS)
> PredRes<-function(D1)
> {
> with(D1, {lmemod<-lme(distance~age*Sex, random=~1|Subject,
> data=subset(D1,age!=14), method="ML")
> themod<-stepAIC(lmemod,dir="both")
> summary(themod)
> prs=predict(themod,newdata=subset(D1,age==14))
> obs<-subset(D1,age==14)$distance
> print(mean(obs-prs))
> } )
> }
>
> but the error message still reads:
>
>> PredRes(Orthodont)
> Start:  AIC=345.12
> distance ~ age * Sex
>
> Error in subset(D1, age != 14) : object "D1" not found
>
> Does this mean that the only option for making a function of this type
> work is to use the "<<-" operator?
> --
> Regards
> Jorunn
>
>
>
> 2008/5/9 Prof Brian Ripley <ripley at stats.ox.ac.uk>:
>> It's a known scoping issue in lme -- you are doing this from a function.
>>  Make sure your dataset is visible -- e.g. use with().
>>
>>
>>
>>  On Fri, 9 May 2008, Jorunn Slagstad wrote:
>>
>>
>>> 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 av 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
>>>
>>> ______________________________________________
>>> R-help at r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>>> and provide commented, minimal, self-contained, reproducible code.
>>>
>>>
>>
>>  --
>>  Brian D. Ripley,                  ripley at stats.ox.ac.uk
>>  Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
>>  University of Oxford,             Tel:  +44 1865 272861 (self)
>>  1 South Parks Road,                     +44 1865 272866 (PA)
>>  Oxford OX1 3TG, UK                Fax:  +44 1865 272595
>>
>>
>>
>>  ______________________________________________
>>  R-help at r-project.org mailing list
>>  https://stat.ethz.ch/mailman/listinfo/r-help
>>  PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>>  and provide commented, minimal, self-contained, reproducible code.
>>
>

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595