[R-meta] residuals in lmer versus rma.mv

[Viechtbauer, Wolfgang (NP)]

Hi Martin,

For rma.mv(), you should use:

m.rma <- rma.mv(Y  ~ X,
                V=0, random = list(~ 1|Paper/PubID, ~ 1|ID),
                R=list(ID = Rmat), Rscale=FALSE,
                data=dat, method = "REML", digits=5, sparse = TRUE)

But you still won't get the same residuals, because residuals(m.lmer) computes 'conditional residuals' (Y minus 'fitted values based on the fixed effects plus random effects'), while residuals(m.rma) computes 'marginal residuals' (Y minus 'fitted values based on the fixed effects').

For 'rma.uni' models, rstandard(<model>, type="conditional") also computes conditional residuals, but I have not yet implemented this for 'rma.mv' models. But you can get these manually as follows:

dat$`Paper/PubID` <- paste0(dat$Paper, "/", dat$PubID)
re <- ranef(m.rma)
Res.rma <- dat$Y - (fitted(m.rma) +
re$`Paper`$intrcpt[match(dat$`Paper`, row.names(re$Paper))] + 
re$`Paper/PubID`$intrcpt[match(dat$`Paper/PubID`, row.names(re$`Paper/PubID`))])
head(Res.rma)

The minor differences arise due to slightly different parameter estimates, but that can be ignored.

If you install the 'devel' version, you can do the slightly simpler:

re <- ranef(m.rma, expand=TRUE)
Res.rma <- dat$Y - (fitted(m.rma) + re$`Paper`$intrcpt + re$`Paper/PubID`$intrcpt)
head(Res.rma)

(note that the 'expand' argument is not currently documented).

Best,
Wolfgang

>-----Original Message-----
>From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-bounces using r-project.org] On
>Behalf Of Martineau, Roger (AAFC/AAC) via R-sig-meta-analysis
>Sent: Thursday, 02 March, 2023 12:12
>To: r-sig-meta-analysis using r-project.org
>Cc: Martineau, Roger (AAFC/AAC)
>Subject: [R-meta] residuals in lmer versus rma.mv
>
>Dear list-members,
>
>I extracted residuals from a lmer model and would like to get the same (or very
>similar) numbers from an equivalent rma.mv model. The example is reported below
>and the residuals from lmer and rma.mv are different:
>
>> head(Res.lmer)
>       1        2        3        4        5        6
> 0.20559  0.03306 -0.12163 -0.30483  0.10283  0.35307
>
>
>> head(Res.rma)
>
>        1         2         3         4         5         6
>
> 0.438639  0.238639  0.089245 -0.111968  0.338639  0.586820
>
>What should I do to get 0.205, 0.033, -0.122, etc as the residuals from the
>rma.mv model ?
>
>Thanks in advance,
>
>Roger  :)
>
>Dr Roger Martineau, DVM PhD
>Centre de recherche et de d�veloppement de Sherbrooke / Sherbrooke Research and
>Development Centre
>Agriculture et Agroalimentaire Canada/Agriculture and Agri-Food Canada
>2000 rue College / 2000 College Street
>Sherbrooke, Qu�bec, Canada,  J1M 0C8
>R�ceptionniste/Receptionist: 819-565-9171
>
># Example
>dat <- iris
>dat <- within(dat, {
>  Y       <- Sepal.Length
>  X       <- Petal.Length
>  TID     <- 1:nrow(dat)
>  PubID   <- rep(c(1:5), times = 30)
>  Paper   <- rep(c(1:10), each = 15)
>  TID     <- as.factor(TID)
>  PubID   <- as.factor(PubID)
>  Paper   <- as.factor(Paper)
>  Wgt     <- round(Sepal.Length)^0.5
>  ID      <- 1:nrow(dat)
>})
>dat <- subset (dat, select = -
>c(Sepal.Length,Sepal.Width,Petal.Length,Petal.Width,Species))
>str(dat)
>
># lmer model
>m.lmer <- lmer(Y ~ X + (1|Paper/PubID), weight = Wgt, data=dat)
>summary(m.lmer)
>Res.lmer <- residuals(m.lmer)
>head(Res.lmer)
>
># rma.mv model
>Rmat <- diag(1/dat$Wgt)
>rownames(Rmat) <- colnames(Rmat) <- dat$ID
>m.rma <- rma.mv(Y  ~ X,
>                V=0, random = list(~ 1|Paper, ~ 1|PubID, ~ 1|ID),
>                R=list(ID = Rmat), Rscale=FALSE,
>                data=dat, method = "REML", digits=5, sparse = TRUE)
>summary(m.rma)
>Res.rma <- residuals.rma(m.rma)
>head(Res.rma)