[R-meta] Query model structure rma.mv

[Sybryn Maes]

Dear Wolfgang,

Thanks a lot for your reply and clarificaitons. About part 3): anova shows
there is a significant improvement in the model fit when including the CAR
structure and acf() now indeed shows reductions in AC. I now feel confident
of the final model.

All the best,
Sybryn




Op di 13 jul. 2021 om 14:17 schreef Viechtbauer, Wolfgang (SP) <
wolfgang.viechtbauer using maastrichtuniversity.nl>:

> 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
>

	[[alternative HTML version deleted]]