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

[Viechtbauer, Wolfgang (NP)]

Depends on how one counts :)

But I would say this is a 3-level model.

Best,
Wolfgang

>-----Original Message-----
>From: Martineau, Roger (AAFC/AAC) [mailto:roger.martineau using AGR.GC.CA]
>Sent: Thursday, 02 March, 2023 16:25
>To: Viechtbauer, Wolfgang (NP)
>Cc: R Special Interest Group for Meta-Analysis
>Subject: RE: residuals in lmer versus rma.mv
>
>Dear Wolfgang,
>
>As written, is this a 3- or a 2-level meta-regression model ?
>
>Best regards,
>
>Roger ☺
>
>-----Original Message-----
>From: Viechtbauer, Wolfgang (NP) <wolfgang.viechtbauer using maastrichtuniversity.nl>
>Sent: Thursday, March 2, 2023 9:49 AM
>To: R Special Interest Group for Meta-Analysis <r-sig-meta-analysis using r-
>project.org>
>Cc: Martineau, Roger (AAFC/AAC) <roger.martineau using AGR.GC.CA>
>Subject: RE: residuals in lmer versus rma.mv
>
>As always, please respond to the mailing list, not just to me.
>
>I reran the code again (with random = list(~ 1|Paper/PubID, ~ 1|ID)) and I get
>the same results also using the devel version as described:
>
>> 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)
>          1           2           3           4           5           6
> 0.20558777  0.03305894 -0.12163497 -0.30482782  0.10283645  0.35307630
>>
>> re <- ranef(m.rma, expand=TRUE)
>> Res.rma <- dat$Y - (fitted(m.rma) + re$`Paper`$intrcpt +
>> re$`Paper/PubID`$intrcpt)
>> head(Res.rma)
>          1           2           3           4           5           6
> 0.20558777  0.03305894 -0.12163497 -0.30482782  0.10283645  0.35307630
>
>So not sure what the issue is.
>
>As for 2) That goes a bit beyond what I can offer here.
>
>Best,
>Wolfgang
>
>>-----Original Message-----
>>From: Martineau, Roger (AAFC/AAC) [mailto:roger.martineau using AGR.GC.CA]
>>Sent: Thursday, 02 March, 2023 14:35
>>To: Viechtbauer, Wolfgang (NP)
>>Subject: RE: residuals in lmer versus rma.mv
>>
>>ATTACHMENT(S) REMOVED: For Wolfgang.docx
>>
>>Dear Wolfgang,
>>
>>Thank you for your quick response.
>>
>>1. I ran the 2 scripts to get the conditional residuals; I get the same
>>results as the lmer model with the long script but not with the one
>>from the devel version, I didn’t find why (see below).
>>2. I saw you corrected the random argument, I will redo my stats
>>accordingly. My revised paper is due March 20th. A reviewer asked to
>>detail further the statistical methodology in the text. Would you be
>>kind enough to review the attached text that describes the statistical
>>methodology. Indeed, your contribution is reported in the
>>Acknowledgments section. I am not sure any more if it is a 3- or
>>2-level meta-regression and don’t want to publish something wrong.
>>
>>Thank you in advance,
>>
>>Roger ☺
>>
>>> 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)
>>        1         2         3         4         5         6
>> 0.205588  0.033059 -0.121635 -0.304828  0.102836  0.353076
>>> re <- ranef(m.rma, expand=TRUE)
>>> Res.rma <- dat$Y - (fitted(m.rma) + re$`Paper`$intrcpt +
>>re$`Paper/PubID`$intrcpt)
>>> head(Res.rma)
>>       1        2        3        4        5        6
>> 0.20559 -0.10791 -0.17307 -0.13438  0.45520  0.68216
>>
>>Thank you so much,
>>
>>Roger ☺
>>
>>-----Original Message-----
>>From: Viechtbauer, Wolfgang (NP)
>>wolfgang.viechtbauer using maastrichtuniversity.nl
>>Sent: Thursday, March 2, 2023 7:22 AM
>>To: R Special Interest Group for Meta-Analysis
>>r-sig-meta-analysis using r-project.org
>>Cc: Martineau, Roger (AAFC/AAC) roger.martineau using AGR.GC.CA
>>Subject: RE: residuals in lmer versus rma.mv
>>
>>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)