[R-meta] HKSJ adjusted error by hand

[Yefeng Yang]

Hi Wolfgang,

Thank you for clarifying the misconception I made about the relationship between the scale factor in the K&H method and H^2. Your illustration looks pretty good.

Interesting to see this adjustment does not improve the inference for multilevel/multivariate models. Regarding the Satterthwaite approximation to the dfs you mentioned in earlier email, actually I have done a little bit of investigation about it before. The main point I found is that using more than 500 meta-analyses across three disciplines, using the combination of Satterthwaite approximation with multilevel models (in my case, three levels) performs as well as the combination of cluster-robust variance estimation with multilevel models. I used James's clubSandwich package, which uses both robust error and Satterthwaite approximation of dfs for the test of coefficient. This makes me  think that maybe adjusting dfs is enough in terms of hypothesis testing? I am not aware of whether James has done any simulation studies to assess the separate contribution of robust errors and dfs adjustment.

Best,
Yefeng
________________________________
From: R-sig-meta-analysis <r-sig-meta-analysis-bounces using r-project.org> on behalf of Viechtbauer, Wolfgang (NP) via R-sig-meta-analysis <r-sig-meta-analysis using r-project.org>
Sent: Friday, September 27, 2024 19:58
To: R Special Interest Group for Meta-Analysis <r-sig-meta-analysis using r-project.org>
Cc: Viechtbauer, Wolfgang (NP) <wolfgang.viechtbauer using maastrichtuniversity.nl>
Subject: Re: [R-meta] HKSJ adjusted error by hand

Just to clarify one misconception: The scale factor in the K&H method is not H^2. Using the traditional definition of H^2 as Q/df, the computation of H^2 uses the weights from a fixed-effects model and the residuals are also computed accordingly:

dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)

rma_fit <- rma(yi, vi, mods = ~ ablat + year, data=dat, method="FE")
wi <- 1 / rma_fit$vi
ri <- residuals(rma_fit)
df <- rma_fit$k - rma_fit$p
# these are all identical
sum(wi * ri^2) / df
rma_fit$H2
rma_fit$QE / df

But in the K&H method, the scale factor can be thought of as a generalized H^2 statistic that uses the results from a random/mixed-effects model:

rma_fit <- rma(yi, vi, mods = ~ ablat + year, data=dat, method="DL")
wi <- 1 / (rma_fit$vi + rma_fit$tau2)
ri <- residuals(rma_fit)
df <- rma_fit$k - rma_fit$p
# this is not the same as H^2 or Q/df
sum(wi * ri^2) / df
rma_fit$H2
rma_fit$QE / df

(of course then sum(wi * ri^2) / df will change depending on which tau^2 estimator we use and also the equivalence between H2 and Q/df breaks down once we use a different estimator than DL, but this is a separate issue).

In any case, the computation of sum(wi * ri^2) / df can be easily generalized to multilevel/multivariate models, leading to a K&H method for such models. Now that I looked at this again, the df of k-p also follow in that case from the distribution of quadratic forms, so it is what it is. And based on simulation studies I have done, this approach did not work nearly as well as it does for standard random/mixed-effects models, so I never bothered to write this up. Might be useful though to at least put this out there, so that efforts to improve inference methods for such models can be focused elsewhere.

Best,
Wolfgang

> -----Original Message-----
> From: R-sig-meta-analysis <r-sig-meta-analysis-bounces using r-project.org> On Behalf
> Of Viechtbauer, Wolfgang (NP) via R-sig-meta-analysis
> Sent: Thursday, September 26, 2024 17:47
> To: R Special Interest Group for Meta-Analysis <r-sig-meta-analysis using r-
> project.org>; James Pustejovsky <jepusto using gmail.com>
> Cc: Viechtbauer, Wolfgang (NP) <wolfgang.viechtbauer using maastrichtuniversity.nl>
> Subject: Re: [R-meta] HKSJ adjusted error by hand
>
> I have played around with this idea and this is even implemented in rma.mv():
>
> dat <- dat.konstantopoulos2011
> rma.mv(yi, vi, random = ~ 1 | district/school, data=dat, test="knha")
>
> As you will see in the warning message, this is experimental. And based on some
> simulations I have done, this doesn't actually work that great. Maybe this could
> be improved by using a Satterthwaite approximation to the dfs. But even then, I
> am not sure if this will do all that much better than:
>
> rma.mv(yi, vi, random = ~ 1 | district/school, data=dat, test="t",
> dfs="contain")
>
> Best,
> Wolfgang
>
> > -----Original Message-----
> > From: R-sig-meta-analysis <r-sig-meta-analysis-bounces using r-project.org> On
> Behalf
> > Of Yefeng Yang via R-sig-meta-analysis
> > Sent: Thursday, September 26, 2024 17:35
> > To: James Pustejovsky <jepusto using gmail.com>; R Special Interest Group for Meta-
> > Analysis <r-sig-meta-analysis using r-project.org>
> > Cc: Yefeng Yang <yefeng.yang1 using unsw.edu.au>
> > Subject: Re: [R-meta] HKSJ adjusted error by hand
> >
> > Hi James,
> >
> > Thank you for you quick and helpful code.
> >
> > Interesting to see the scale factor (S_sq  in you notation) is just the
> > heterogeneity measure H^2, which is equal to Q / df. If this is the case, I am
> > wondering why this adjustment has not been extended to multilevel/multivariate
> > meta-analysis, given that one can also derive the scale factor from these
> > relative complex meta-analyses?
> >
> > Best,
> > Yefeng
> >
> > ________________________________
> > From: James Pustejovsky <jepusto using gmail.com>
> > Sent: Friday, September 27, 2024 0:41
> > To: R Special Interest Group for Meta-Analysis <r-sig-meta-analysis using r-
> > project.org>
> > Cc: Yefeng Yang <yefeng.yang1 using unsw.edu.au>
> > Subject: Re: [R-meta] HKSJ adjusted error by hand
> >
> > The adjustment is just a scalar that can be calculated "by hand" from the
> > weights and residuals of the model. Here's an example, which could be turned
> > into a function:
> >
> > library(metafor)
> >
> > dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)
> >
> > rma_fit <- rma(yi, vi, mods = ~ ablat + year, data=dat)
> >
> > # Calculate variance factor
> > wi <- 1 / (rma_fit$vi + rma_fit$tau2)
> > ri <- residuals(rma_fit)
> > df <- rma_fit$k - rma_fit$p
> > S_sq <- sum(wi * ri^2) / df
> >
> > # Scale the variance-covariance matrix
> > V_hksj <- S_sq * vcov(rma_fit)
> >
> > # Check against rma.uni()
> > rma_hksj <- update(rma_fit, test = "hksj")
> > all.equal(V_hksj, vcov(rma_hksj))
> >
> > Depending on your preferred settings, you might need to adjust by taking
> > max(S_sq, 1).
> >
> > James
> >
> > On Thu, Sep 26, 2024 at 9:03 AM Yefeng Yang via R-sig-meta-analysis <r-sig-
> meta-
> > analysis using r-project.org<mailto:r-sig-meta-analysis using r-project.org>> wrote:
> > Dear community,
> >
> > I am wondering whether there is a post-hoc way to calculate sampling variances
> > of the estimated regression coefficients from meta-analysis models based on
> the
> > Hartung-Knapp-Sidik-Jonkman method.
> >
> > To be more precise,
> >
> >   1.
> > I first fit a normal random-effect meta-analysis via `rma()` in metafor
> package:
> >
> > library(metafor)
> > dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)
> > mod <- rma(yi, vi, data=dat),
> >
> >   1.
> > then, how can I get the adjusted standard error of the estimated coefficients
> > (in this case, it is an intercept) based on the model object `mod`?
> >
> > Of course, we can get it using an on-the-fly way:
> > rma(yi, vi, data=dat, test = "hksj")
> >
> > But I have a couple of big datasets and want to report both the original
> > standard and adjusted errors. I do not have a high-performance PC and I would
> > like to avoid 're-fitting' the meta-analysis.
> >
> > Best,
> > Yefeng
_______________________________________________
R-sig-meta-analysis mailing list @ R-sig-meta-analysis using r-project.org
To manage your subscription to this mailing list, go to:
https://stat.ethz.ch/mailman/listinfo/r-sig-meta-analysis

	[[alternative HTML version deleted]]