Search code examples
ranovascopingmixed-models

Scoping-related (?): anova() on list of created mixed-effects models


In a project where I'm performing mixed-effects modelling using lme, I'm trying to compare models with different correlation structures and equal fixed parts. As I'll be building a lot of these models (for different dependent variables), I tried to write a function to generate a list of models with different correlation structures, as in the example below (I really tried to keep it to a minimum working example).

If I run an anova() on the elements of this list, this works, but only if fixedPart is in my global environment. Why is this the case? Is there a way to circumvent this problem, so that I can just keep m and re-use/delete fixedPart?

I presume this problem is related to the (lexical) scoping in R, but I cannot find a way to actually fix it.

Thanks in advance!

#Dependencies
library(multilevel)
library(multcomp)

#Generate sample data
nVals = 100
sData = rnorm(nVals, mean = 1, sd = 1)
dF <- data.frame(nSubject = 1:nVals, 
                 v1data = sData + rnorm(nVals, mean = 0, sd = 0.1),
                 v2data = sData + rnorm(nVals, mean = 0, sd = 0.1),
                 v3data = sData + rnorm(nVals, mean = 0, sd = 0.4))
dLongF = reshape(data=dF, varying=c("v1data","v2data","v3data"), v.names='data', direction="long", idvar="nSubject", times=1:3)

#Define function to assess different covariance structures
doAllCorrModels <- function(dataF, subjVarName, visitVarName, fixedPart){
    mList <- vector("list",2)
    mList[[1]] <- lme(fixedPart, #Random intercept, homogeneous variance
             random=as.formula(paste("~1|", subjVarName)),
             data=dataF,
             weights=NULL)
    mList[[2]] <- lme(fixedPart, #Random intercept, heterogeneous variance
             random=as.formula(paste("~1|", subjVarName)),
             data=dataF,
             weights=varIdent(form = as.formula(paste("~1|", visitVarName)))
    )
    mList
}

#Get different covariance structures
dataF <- dLongF
subjVarName <- "nSubject"
visitVarName <- "time"
fixedPart <- data ~ time
m <- doAllCorrModels(dataF, subjVarName, visitVarName, fixedPart)

#This works:
a1 <- anova(m[[1]], m[[2]])

#But this does not:
rm(fixedPart)
a2 <- anova(m[[1]], m[[2]])

Solution

  • You can avoid this by using do.call:

    doAllCorrModels <- function(dataF, subjVarName, visitVarName, fixedPart){
      mList <- vector("list",2)
      mList[[1]] <- do.call(lme, list(fixed = fixedPart,
                                      random=as.formula(paste("~1|", subjVarName)),
                                      data=dataF,
                                      weights=NULL))
    
      mList[[2]] <- do.call(lme, list(fixed = fixedPart,
                                      random=as.formula(paste("~1|", subjVarName)),
                                      data=dataF,
                                      weights=varIdent(form = as.formula(paste("~1|", visitVarName)))))
      mList
    }