As I am wrapping up with the growth curve models at hand, here are a couple of R functions to share with whoever is still using lmer from the pre-1.0 version lme4 (not ready to upgrade as yet--too many codes would require updating). These functions were written to summarize results from different lmer models.
The first three functions separately extract the model summary statistics (lmer.stats), the fixed effect parameters (lmer.fixef), and the random effect parameters (lmer.ranef) into data frames. The last function lmer.append can combine these results into an aggregated data frame, which can then be saved as a spreadsheet using the xlsx package.
Note: the if condition in lmer.ranef needs revision to make the columns consistent if you have more than one covariance term in any of the models. Otherwise R won't be able to aggregate the data frames.
lmer.stats<-function(lmer.name) {
A<-AIC(lmer.name)
B<-BIC(lmer.name)
ll<-logLik(lmer.name)
dg<-attr(ll,"df")
dv<-deviance(lmer.name)
obs.TIME<-length(lmer.name@y)
obs.CHILD<-sapply(ranef(lmer.name),nrow)[1]
names(obs.CHILD)<-NULL
obs.SCHOOL<-sapply(ranef(lmer.name),nrow)[2]
names(obs.SCHOOL)<-NULL
label<-deparse(substitute(lmer.name)) # identifier
df<-data.frame(label, "AIC"=A, "BIC"=B, "LL"=ll, "DEV"=dv, "df"=dg, "N"=obs.TIME, "CHILD"=obs.CHILD, "SCHOOL"=obs.SCHOOL)
}
lmer.ranef<-function(lmer.name){
re<-data.frame(summary(lmer.name)@REmat)
re<-subset(re,select=-Name)
label<-deparse(substitute(lmer.name)) # identifier
nr<-nrow(summary(lmer.name)@REmat)
md<-data.frame(rep(label,nr))
colnames(md)<-"Model"
dfr<-data.frame(cbind(md,re))
if (ncol(dfr)==4) { # random slope models have more columns
corr.col<-data.frame(rep(NA,nr))
colnames(corr.col)<-"Corr"
V6.col<-data.frame(rep(NA,nr))
colnames(V6.col)<-"V6"
dfr<-data.frame(cbind(dfr,corr.col,V6.col))
} else {
dfr<-dfr
}
}
lmer.fixef<-function(lmer.name){
beta<-data.frame("Beta"=fixef(lmer.name))
se<-data.frame("S.E."=sqrt(diag(vcov(lmer.name))))
vars<-data.frame(row.names(beta))
colnames(vars)<-"Variable"
vars$Variable<-gsub("\\)", "", vars$Variable) # deal with (Intercept)
vars$Variable<-gsub("\\(", "", vars$Variable)
label<-deparse(substitute(lmer.name)) # identifier
md<-data.frame(rep(label,length(lmer.name@fixef)))
colnames(md)<-"Model"
row.names(beta)<-NULL
dff<-data.frame(cbind(md,vars,beta,se))
}
lmer.append<-function(...,append=TRUE) {
label<<-deparse(substitute(...))
if (!append){
L.stats<<-lmer.stats(...)
L.ranef<<-lmer.ranef(...)
L.fixef<<-lmer.fixef(...)
} else {
L.stats<<-rbind(L.stats, lmer.stats(...))
L.ranef<<-rbind(L.ranef, lmer.ranef(...))
L.fixef<<-rbind(L.fixef, lmer.fixef(...))
}
}
Added 06/09/2014:
Someone just reminded me that the functionality of lmer.stat is similar to the base routine anova(). Say you have lmer model estimates A1, A2, and A3, anova(A1,A2,A3) returns a data frame that summarizes the degrees of freedom, AIC, BIC, log-likelihood, and results of deviance tests in relation to the first model, in this case, A1. Still, that does not include the other statistics that lmer.stats provides.
No comments:
Post a Comment