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]])
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
}