[R-meta] estimating tau2 under no effect

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 at 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 at r-project.org>
Cc: Filippo Gambarota <filippo.gambarota at 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/>