I'm currently trying to calculate the effect size (Cohen's f^2) of a given effect but need to run a null and partial model with the random effects pre-specified in order to do so (according to Selya et al., 2012 for dealing with continuous predictors). Selya et al. outline the code necessary to do so in SAS, but I'm trying to figure out how one would do it in R.
When I went to use the last piece of code from a prior similar question, I kept getting the error that the thetas do not match ("3!=4"). I think the issue arises because my original model has a cross-level interaction term while the prior post has only random intercepts. I'm not just trying to hold constant the random effects on the intercept but also the random effects of a slope. How can I amend the code to make this run? I called the getME function with "theta" and saw indeed that my mod3b model has 4 theta values listed, so I imagine I will need to add another parameter term to the code. I just cannot figure out how to get the variance of the fourth theta term to show in my original output. I appreciate any and all help!
Here's my original cross-classified model with random slopes & random intercepts:
mod3b <- lmer(FitBelong~Condition*Gender +
(1+Condition|ResponseID) + (1|Stimuli),
data=LFS1Ensemble, REML=TRUE)
summary(mod3b)
I've adapted code from the answer on the link to be as follows:
#Effect size of interaction#
buildMM <- function(theta) {
dd <- as.function(mod3b)
ff <- dd(theta)
opt <- list(par=c(0,0,0),fval=ff,conv=0)
mm <- mkMerMod(environment(dd), opt, lmod$reTrms, fr = lmod$fr,
mc = quote(hacked_lmer()))
return(mm)
}
objfun <- function(x,target=c(3.92244,0.08805,0.09683)) {
mm <- buildMM(sqrt(x))
return(sum((unlist(VarCorr(mm))-target)^2))
}
s0 <- c(3.92244,0.08805,0.09683)/sigma(mod3b)^2
opt <- optim(fn=objfun,par=s0)
mm_final <- buildMM(sqrt(opt$par))
summary(mm_final)
The error it throws is: "Error: theta size mismatch"
When I do the traceback it gives me:
6 stop(structure(list(message = "theta size mismatch", call = NULL,
cppstack = NULL), .Names = c("message", "call", "cppstack"
), class = c("std::invalid_argument", "C++Error", "error", "condition"
)))
5 dd(theta)
4 buildMM(sqrt(x))
3 fn(par, ...)
2 (function (par)
fn(par, ...))(c(2.1571475240413, 0.0484231344499436, 0.0532516991344468
))
1 optim(fn = objfun, par = s0)
Any and all help much much appreciated!
