[R-meta] Deviance and it d.f. in pairwise meta analysis

[Viechtbauer, Wolfgang (NP)]

Ah, now I understand. Interesting question, since this is not something that is discussed in any textbooks or articles on meta-analysis as far as I know.

First, it may help to clarify how the deviance is computed. Often, deviance is just computed as -2 * log likelihood, but metafor takes the slightly more accurate definition of the deviance as described here:

https://en.wikipedia.org/wiki/Deviance_(statistics)

that is, deviance is 2 times the difference of the log likelihood of the 'saturated model' compared to the fitted model of interest. The saturated model is the model that has a parameter for every observation and hence fits perfectly. So let's manually compute the deviance for an equal-effects model:

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

# fit the equal-effects model
res.ee <- rma(yi, vi, data=dat, method="EE")
res.ee

# fit the saturated model
res.s <- rma(yi, vi, mods = ~ 0 + factor(trial), data=dat, method="FE")
res.s

# computation of the deviance
2 * c(logLik(res.s) - logLik(res.ee))

# compare with fitstats()
fitstats(res.ee)

# deviance = Q statistic
res.ee$QE

Note that in this case, the deviance is identical to the Q statistic. Hence, under the null hypothesis of homogeneity, the deviance follows a chi-square distribution with k-1 degrees of freedom (just like the Q statistic). To be precise, this assumes unbiased estimates, normal sampling distributions, and known sampling variances, so the above is only really true when the sample sizes of the individual studies are sufficiently large. On the other hand, k is not relevant here, so the above is true even when k=2.

Under a random-effects model, the computation is the same:

# fit the random-effects model
res.re <- rma(yi, vi, data=dat, method="ML")
res.re

# computation of the deviance
2 * c(logLik(res.s) - logLik(res.re))

# compare with fitstats()
fitstats(res.re)

The saturated model is also res.s (one cannot use method="ML" for the saturated model, since then the model is overparameterized, as tau^2 must be 0 in this model by definition).

The distribution of the deviance is complicated now, since even when the null hypothesis of homogeneity is true, the estimate of tau^2 will be larger than 0 in about half of the cases, leading to a deviance that is closer to res.s than that of res.ee. However, as k gets larger, the estimate of tau^2 also gets increasingly closer to 0, so asymptotically (i.e., when k is large) the deviance will also follow a chi^2 distribution with df = k-1.

As for your last question: df.residual() always returns k-p, irrespective of the model you have fit (with or without tau^2).

Best,
Wolfgang

> -----Original Message-----
> From: Marimuthu S <sm using mcmaster.ca>
> Sent: Wednesday, March 5, 2025 16:15
> To: Viechtbauer, Wolfgang (NP) <wolfgang.viechtbauer using maastrichtuniversity.nl>; R
> Special Interest Group for Meta-Analysis <r-sig-meta-analysis using r-project.org>
> Subject: Re: Deviance and it d.f. in pairwise meta analysis
>
> Hello Dr. Wolfgang,
>
> Thank you for your response. I apologize for any confusion. I meant to ask about
> the theoretical distribution of deviance statistics.
> Additionally, I am wondering, do the residual degrees of freedom account for the
> heterogeneity parameter (tau)?
>
> Marimuthu
> ________________________________________
> From: Viechtbauer, Wolfgang (NP)
> <mailto:wolfgang.viechtbauer using maastrichtuniversity.nl>
> Sent: Wednesday, March 5, 2025 8:36 AM
> To: R Special Interest Group for Meta-Analysis <mailto:r-sig-meta-analysis using r-
> project.org>
> Cc: Marimuthu S <mailto:sm using mcmaster.ca>
> Subject: RE: Deviance and it d.f. in pairwise meta analysis
>
> Dear Marimuthu,
>
> You can use df.residual() to obtain these:
>
> > df.residual(res1)
> [1] 12
>
> I am not sure what you mean by 'its distribution'. The df are not a statistic,
> so they do not have a distribution.
>
> By the way, no need for REML=F in fitstats() if the model was fitted with ML
> (fitstats() automatically then provides the ML values).
>
> Best,
> Wolfgang
>
> > -----Original Message-----
> > From: R-sig-meta-analysis <mailto:r-sig-meta-analysis-bounces using r-project.org>
> On Behalf
> > Of Marimuthu S via R-sig-meta-analysis
> > Sent: Monday, March 3, 2025 22:12
> > To: Marimuthu S via R-sig-meta-analysis <mailto:r-sig-meta-analysis using r-
> project.org>
> > Cc: Marimuthu S <mailto:sm using mcmaster.ca>
> > Subject: [R-meta] Deviance and it d.f. in pairwise meta analysis
> >
> > Hello Everyone,
> >
> > I fitted meta-analysis model using rma() function. I can extract all the fit
> > statistics using fitstats() function. But it doesn't give degrees of freedom
> > (d.f) for deviance.
> >
> > Could you please let me know if it's possible to calculate the deviance
> degrees
> > of freedom (d.f.) and derive its distribution in pairwise meta-analysis? If
> it's
> > not feasible, I'd appreciate understanding the reasons behind it. Any insights
> > or leads would be greatly appreciated. Thank you for your time!
> >
> > Here is my code:
> >
> > library(metafor)
> >
> > dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)
> > ### fit random-effects model
> > res1 <- rma(yi, vi, data=dat, method="ML")
> >
> > > fitstats(res1, REML=F)
> >                  ML
> > logLik:   -12.66508
> > deviance:  37.11602
> > AIC:       29.33015
> > BIC:       30.46005
> > AICc:      30.53015
> >
> > Regards,
> >
> > Marimuthu,
> > Ph.D. candidate
> > McMaster University, Canada