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

3 messages · Martineau, Roger, Wolfgang Viechtbauer

Thu, Mar 2, 2023 6:49 AM #

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:

re$`Paper`$intrcpt[match(dat$`Paper`, row.names(re$Paper))] +
  re$`Paper/PubID`$intrcpt[match(dat$`Paper/PubID`, row.names(re$`Paper/PubID`))])

1           2           3           4           5           6
 0.20558777  0.03305894 -0.12163497 -0.30482782  0.10283645  0.35307630

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 at 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 at maastrichtuniversity.nl
Sent: Thursday, March 2, 2023 7:22 AM
To: R Special Interest Group for Meta-Analysis r-sig-meta-analysis at r-project.org
Cc: Martineau, Roger (AAFC/AAC) roger.martineau at 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 at 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 at 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)

Martineau, Roger

Thu, Mar 2, 2023 7:25 AM #

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 at maastrichtuniversity.nl> 
Sent: Thursday, March 2, 2023 9:49 AM
To: R Special Interest Group for Meta-Analysis <r-sig-meta-analysis at r-project.org>
Cc: Martineau, Roger (AAFC/AAC) <roger.martineau at 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:

re$`Paper`$intrcpt[match(dat$`Paper`, row.names(re$Paper))] +
  re$`Paper/PubID`$intrcpt[match(dat$`Paper/PubID`, row.names(re$`Paper/PubID`))])

1           2           3           4           5           6
 0.20558777  0.03305894 -0.12163497 -0.30482782  0.10283645  0.35307630

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 at 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 at maastrichtuniversity.nl
Sent: Thursday, March 2, 2023 7:22 AM
To: R Special Interest Group for Meta-Analysis 
r-sig-meta-analysis at r-project.org
Cc: Martineau, Roger (AAFC/AAC) roger.martineau at 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 at 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 at 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)

Wolfgang Viechtbauer

Thu, Mar 2, 2023 7:34 AM #

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 at 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 at maastrichtuniversity.nl>
Sent: Thursday, March 2, 2023 9:49 AM
To: R Special Interest Group for Meta-Analysis <r-sig-meta-analysis at r-
project.org>
Cc: Martineau, Roger (AAFC/AAC) <roger.martineau at 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 at 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 at maastrichtuniversity.nl
Sent: Thursday, March 2, 2023 7:22 AM
To: R Special Interest Group for Meta-Analysis
r-sig-meta-analysis at r-project.org
Cc: Martineau, Roger (AAFC/AAC) roger.martineau at 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 at 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 at 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)