[R-meta] residuals in lmer versus rma.mv
Martineau, Roger (AAFC/AAC)
roger@m@rt|ne@u @end|ng |rom AGR@GC@CA
Thu Mar 2 16:25:26 CET 2023
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
CAUTION: This email originated from outside of the organization. Do not click links or open attachments unless you recognize the sender and know the content is safe.
ATTENTION: Ce courriel provient de l’extérieur de l’organisation. Ne cliquez pas sur les liens et n’ouvrez pas les pièces jointes à moins que vous ne reconnaissiez l’expéditeur et que vous sachiez que le contenu est sûr.
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)
More information about the R-sig-meta-analysis
mailing list