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