[R-meta] Query model structure rma.mv

[Viechtbauer, Wolfgang (SP)]

Dear Sybryn,

Please see below for my responses.

Best,
Wolfgang

>-----Original Message-----
>From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-bounces using r-project.org] On
>Behalf Of Sybryn Maes
>Sent: Thursday, 08 July, 2021 16:28
>To: r-sig-meta-analysis using r-project.org
>Subject: [R-meta] Query model structure rma.mv
>
>Dear colleagues,
>
>I have three questions regarding a meta-analysis multivariate model
>structure I am using for my dataset that is complex in the following way:
>There are multiple sites, and within those sites, multiple years of data
>can occur. These years are not always the same, so one site can have data
>for eg 2003 and 2004, whereas another site could have data only for 2005.
>Another one can have data for 2001 to 2010 every year available.
>
>I propose the following structure for my data :
>
>*metamodel <- rma.mv(Hedges_SMD, Hedges_SMD_VARIANCE, random=list(~ 1 |
>Site_ID / Observation_ID , ~ Year |  Site_ID ),*
>
>*struct="CAR", data=SMD_SITELEVEL)*
>
>I am including the second part of the random structure list to account for
>potential autocorrelation and non-independence between the different years
>within a same site.
>
>However, I have the following questions when looking at the output and when
>trying to visualize the results (estimates + CIs) with a simple forest plot:
>
>1) Does the "*forest*" command also work with complex rma.mv structures -
>i.e. for making forest plots? I have an outcome and do not get errors, but
>I cannot find anywhere whether it is appropriate to use for these multivariate
>objects.

Yes, you can use forest() in combination with rma.mv() objects. The function doesn't do anything special for 'rma.mv' objects, so it just shows the individual estimates (plus it adds the pooled estimate at the bottom if the model doesn't include any moderators). If you want to indicate some kind of grouping visually, you have to take some extra steps. See, for example, the last example at:

https://wviechtb.github.io/metafor/reference/forest.rma.html

>2) When including my CAR structure, the *weights *change when I show them
>in the forest plot (using "showweights=TRUE" to visualize them). This I
>find strange as the sample sizes are not changing in any way. Why could this be?

Adding or changing the random effects structure impacts the weighting structure. This is to be expected.

>3) In order to *diagnostize *whether the non-independence and
>autocorrelation is appropriately taken into account, I thought of running a
>variogram with and without the CAR structure in the model, and an acf()
>command. However, none of these two tools actually shows me that the
>autocorrelation is reduced or removed. What could be the reason for this?
>Could I find out in another way whether it is a "better" model with the CAR
>structure included?

I don't know what commands you ran, but if you looked at the ACF of the residuals within studies, then those will still be autocorrelated, even if the autocorrelation is accounted for in the model. Here is an example (not really a meta-analysis, so 'yi' are raw data, but this illustrates the principle and simplifies things since we don't have to worry about a nested structure of observations within studies):

############################

library(metafor)

sigma <- 0.5
rho   <- 0.6
n     <- 1000
time  <- 1:n
const <- rep(1, n)

yi <- arima.sim(model=list(ar=rho), n=n) * sqrt(1-rho^2) * sigma

# treat all observations as independent
res1 <- rma.mv(yi, 0, random = ~ 1 | time)
res1

# assume an AR1 structure for the observations
res2 <- rma.mv(yi, 0, random = ~ time | const, struct="AR")
res2

acf(resid(res1))
acf(resid(res2))

# both ACFs show the same autocorrelation in the residuals

############################

To 'get rid' of the autocorrelation in the residuals, one has to compute what are sometimes called 'normalized residuals' (see Pinheiro & Bates, 2000) or 'Choleksy residuals'. Continuing with the example:

acf(resid(res1, type="cholesky"))
acf(resid(res2, type="cholesky"))

The ACF for res2 shows no noteworthy autocorrelation.

For meta-analytic data, you would have to look at the residuals within studies. Also, there is the issue of how one should deal with any random effects that are 'higher' up in the hierarchy of your model, that is, when computing the residuals, should one compute them only around the fixed effects or should one remove/partial out the contributions of any additional random effects? This gets complicated quickly, so instead, I would just suggest to run a LRT comparing the two models:

anova(res1, res2)

This will tell you whether the model with the autocorrelation structure fits significantly better than the one without.

>Thank you for any assistance/advice!
>
>Cheers,
>Sybryn