[R-meta] estimating tau2 under no effect

[Viechtbauer, Wolfgang (NP)]

Hi Filippo,

Pretty much all functions that fit 'normal-normal models' in metafor make use of the fact that the MLE of the fixed effects can be profiled out of the likelihood function and hence the optimization only pertains to the random effects. This has the advantage of simplifying the optimization problem, but precludes doing what you like to do, namely to fix a particular fixed effect (in this case, mu) to some value (in this case, 0). At one point, I considered adding this functionality, since this would also be useful for obtaining profile likelhood CIs for the fixed effects, but it would require a lot of internal changes that I didn't want to tackle. Plus, Wald-type CIs for the fixed effects (possibly with some adjustments, such as using a t-distribution for the critical values, the Knapp-Hartung method) work very well, so there would be little to no benefit for a whole lot of work.

However, when I added the possibility to fit location-scale models with rma(), I did also implement this possibility, although this is undocumented. This was mostly for my own research purposes, but hey, if this is useful for you, then here we go:

First, consider the standard random-effects model:

fit <- rma(yi, vi, data = dat)
fit

which we can reformulate as a location-scale model as follows:

fit <- rma(yi, vi, data = dat, scale = ~ 1)
fit
predict(fit, newscale=1, transf=exp)

You will find that the results are identical (except for minor numerical differences that can arise due to the different ways these models are fitted internally). By default, rma() also uses the profiling trick, but this can be switched off and then the optimization proceeds over the fixed and random effects simultaneously:

fit <- rma(yi, vi, data = dat, scale = ~ 1, optbeta=TRUE)
fit
predict(fit, newscale=1, transf=exp)

Again, results should be identical, give or take a few multiples of the Planck constant. Now that we are optimizing over beta (i.e., mu) as well, we can fix its value:

fit <- rma(yi, vi, data = dat, scale = ~ 1, optbeta=TRUE, beta=0)
fit
predict(fit, newscale=1, transf=exp)

There might be functions in metafor that will not behave nicely when such a model object is passed to it. But the above definitely works.

Best,
Wolfgang

> -----Original Message-----
> From: R-sig-meta-analysis <r-sig-meta-analysis-bounces using r-project.org> On Behalf
> Of Filippo Gambarota via R-sig-meta-analysis
> Sent: Friday, March 15, 2024 12:40
> To: R meta <r-sig-meta-analysis using r-project.org>
> Cc: Filippo Gambarota <filippo.gambarota using gmail.com>
> Subject: [R-meta] estimating tau2 under no effect
>
> Hi,
> Maybe a kind of unusual question but I would like to be able to estimate
> the variance of a random-effects model assuming that there is no real
> effect (mu_theta = 0).
> I did this manually maximising the log-likelihood of the random-effect
> model by fixing the average effect to 0 and the estimation is very close to
> the tau2 estimated with sign-flip permutations (that simulate the null
> hypothesis of no effect). I put the code below.
> My question is about, is there any way to do this directly in metafor? I
> was deeply reading the source code of rma.uni and rma.mv but I am not able
> to find a way to trick the model fixing beta to 0.
>
> ```
> ## SIMULATE A DATASET
>
> k <- 50
> n <- rpois(k, 40)
> d <- 0.5
> tau2 <- 0.1
> thetai <- rnorm(k, 0, sqrt(tau2))
>
> yi <- vi <- rep(NA, k)
>
> for(i in 1:k){
>   x0 <- rnorm(n[i], 0, 1)
>   x1 <- rnorm(n[i], d + thetai[i], 1)
>   yi[i] <- mean(x1) - mean(x0)
>   vi[i] <- var(x0)/n[i] + var(x1)/n[i]
> }
>
> dat <- data.frame(
>   id = 1:k,
>   yi,
>   vi
> )
>
> # random-effect model
> fit <- rma(yi, vi, data = dat)
>
> # sign-flip permutations
> nperm <- 1000
> tau2p <- rep(NA, nperm)
> tau2p[1] <- fit$tau2
>
> for(i in 2:nperm){
>   s <- sample(c(-1, 1), k, replace = TRUE)
>   yip <- dat$yi * s
>   fitp <- rma(yip, vi, data = dat)
>   tau2p[i] <- fitp$tau2
> }
>
> # restricted maximum likelihood equation
> meta_REML <- function(x, yi, vi){
>   theta <- x[1]
>   tau2 <- x[2]
>   -0.5 * sum(log(tau2 + vi)) - 0.5 * log(sum(1/(tau2 + vi))) - 0.5 *
> sum((yi - theta)^2 / (tau2 + vi))
> }
>
> # very crude approach, just to see
> tau2i <- seq(0, 1, 0.001)
> ll <- sapply(tau2i, function(t) meta_REML(c(0, t), dat$yi, dat$vi))
>
> # should be similar, tau2 under H0
> tau2i[which.max(ll)]
> median(tau2p)
> ```
> Thank you!
>
> --
> *Filippo Gambarota, PhD*
> Postdoctoral Researcher - University of Padova
> Department of Developmental and Social Psychology
> Website: filippogambarota.xyz
> Research Groups: Colab <http://colab.psy.unipd.it/>   Psicostat
> <https://psicostat.dpss.psy.unipd.it/>