Suppose we have
• An lmer{lme4} model, MyModel, computed on dataframe SomeDATA;
• Dataframe NEWDATA, which contains the variables used in
computing MyModel, but different values for these variables.
In NEWDATA I would like to compute (in an automated, quick, easy,
non-error-prone manner) the fitted values of the fixed portion of the model.
Also I want to compute the standard errors for these fits.
One use of this would be to display graphically the fitted portion of an lmer model so as
to see the "effects" of selected variables. For instance, such a display would
be useful for judging the *practical* significance of variables that have tested as
*statistically* significant.
To this end, I would like to obtain (in an automated manner) the "design
matrix" of the model for the NEWDATA, say NEWDATA_designMatrix. Then the fits would
be
NEWDATA_designMatrix %*% cbind(fixef(MyModel))
and the variance of the fits
NEWDATA_designMatrix %*% vcov(MyModel) %*% t(NEWDATA_designMatrix)
Are not model.frame() and model.matrix() the way to obtain
NEWDATA_designMatrix? Below is an unsuccessful attempt to use these functions.
Can anyone explain how to obtain the desired result?
# Simulate unbalanced longitudinal data which include a factor
(categorical predictor).
set.seed(5)
obsPerID<-5
N_ID<-10
SomeDATA<-data.frame(
ID=rep( 1:N_ID, each=obsPerID)
, X=rnorm( obsPerID * N_ID )
, categ=factor(sample(LETTERS[1:3], size=obsPerID * N_ID, replace=T))
)
SomeDATA$Y <- SomeDATA$X + SomeDATA$ID + rnorm( obsPerID * N_ID ) / 2 +
ifelse(SomeDATA$categ=="C", 5, 0)
SomeDATA<-SomeDATA[order(SomeDATA$ID, SomeDATA$X),]
print(summary(SomeDATA))
library(lattice)
print(
xyplot( Y ~ X | categ, groups=ID, data=SomeDATA, type="l"
, layout=c(3,1)
)
)
# I would like to easily slap the fixed portion of the fits of a regression model on top of this plot.
library(lme4)
MyModel<-lmer( Y ~ X + categ + (1|ID), data=SomeDATA)
print(MyModel)
# Create an orderly NEWDATA.
N_NEWDATA<-7
NEWDATA<-data.frame(
X=seq(-2, 2, length=N_NEWDATA)
, ID=2
, categ=factor(
sample(c("A", "C"), replace=T, size=N_NEWDATA)
, levels=levels(SomeDATA$categ))
, Y=0
)
NEWDATA$nuisance<-sample(1:3, replace=T, size=nrow(NEWDATA))
print(NEWDATA)
NEWDATA_frame<-model.frame( formula=formula(MyModel), data=NEWDATA)
print(NEWDATA_frame)
# The following generates an error:
NEWDATA_matrix<-model.matrix( object=formula(MyModel), data=NEWDATA_frame)
# The following does not conform to NEWDATA but instead reverts to
SomeDATA:
NEWDATA_matrix<-model.matrix( object=MyModel, data=NEWDATA_frame)
print(NEWDATA_matrix)
Thanks for any insights
Jacob A. Wegelin
Assistant Professor
Department of Biostatistics
Virginia Commonwealth University
730 East Broad Street Room 3006
P. O. Box 980032
Richmond VA 23298-0032
APPENDIX
sessionInfo()
R version 2.10.1 (2009-12-14)
i386-apple-darwin9.8.0
locale:
[1] en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] lme4_0.999375-32 Matrix_0.999375-33 lattice_0.17-26
loaded via a namespace (and not attached):
[1] grid_2.10.1 tools_2.10.1
END APPENDIX
______________________________________________
R-help@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.
