Tuesday, June 3, 2014

A few R functions to summarize lmer results

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