I would love to use emmeans with a complicated model from rma.mv() in metafor, and for this there is the function emmprep(). However, the model I have is unbalanced and nested, leading to redundant predictors being dropped from the model due to rank deficiencies. In this case, emmprep() does not work. As per the helpfile: "In this case, one should remove any redundancies in the original model fitted before using this function."
How do I do that? Do I have to use dummy variables to avoid rank deficiencies, or is there a simpler solution?
Here's a reproducible example:
library(metafor)
dat <- dat.bcg #take some example data to modify
#modify to make it work for this example (so now we have two factors)
dat$ablat<-c("a","b","c","a","b","b","a","b","c","a","b","c","a")
#get the yis & vis
dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat,
slab=paste0(author, ", ", year))
#run a nested model
res4 <- rma.mv(yi, vi, mods = ~ alloc/ablat, data=dat)
res4
library(emmeans) #load emmeans
sav <- emmprep(res4) #this returns the redundant predictors error
#"Error in emmprep(res4) :
# Cannot use function when some redundant predictors were dropped from the model."
(PS. I asked this question badly before and have deleted that version to avoid confusion)
It seems you don't really need emmprep() all that much. The qdrg() function in emmeans almost works out-of-the-box. The only hitch is that rma models name the intercept intrcpt instead of (Intercept), and we have to fix that. Here's an example with a dataset I created for the feature request I posted on metafor's GitHub site
library(metafor)
## Loading required package: Matrix
## Loading required package: metadat
## Loading required package: numDeriv
##
## Loading the 'metafor' package (version 4.2-0). For an
## introduction to the package please type: help(metafor)
library(emmeans)
dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)
dat <- transform(dat,
acat = factor(ablat < 40, labels = c("high", "low")),
era = factor(1+as.integer((year-1940)/15),
labels = c("early","mid","late")))
with(dat, table(acat, era)) # there's an empty cell
## era
## acat early mid late
## high 3 2 1
## low 0 2 5
mod <- rma(yi, vi, mods = ~ acat * era, data = dat)
## Warning: Redundant predictors dropped from the model.
emmprep(mod) # fails
## Error in emmprep(mod): Cannot use function when some redundant predictors were dropped from the model.
# However ...
rownames(mod$beta)[1] = "(Intercept)"
RG = qdrg(object = mod) # named argument 'object' is essential
RG
## 'emmGrid' object with variables:
## acat = high, low
## era = early, mid, late
## Some things are non-estimable (null space dim = 1)
emmeans(RG, ~era | acat)
## acat = high:
## era emmean SE df lower.CL upper.CL
## early -0.971 0.244 8 -1.533 -0.409
## mid -1.366 0.347 8 -2.166 -0.566
## late -1.442 0.325 8 -2.190 -0.693
##
## acat = low:
## era emmean SE df lower.CL upper.CL
## early nonEst NA NA NA NA
## mid -0.299 0.340 8 -1.082 0.485
## late -0.268 0.161 8 -0.641 0.104
##
## Confidence level used: 0.95
Created on 2023-07-22 with reprex v2.0.2