Jorunn,
This solutions works with R 2.7.0 under windows
library(MASS)
library(nlme)
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
Belgium
tel. + 32 54/436 185
Thierry.Onkelinx op inbo.be
www.inbo.be
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: r-help-bounces op r-project.org [mailto:r-help-bounces op r-project.org]
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:
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.U
TF-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_ID
ENTIFICATION=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 op 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.