[R-meta] estimating tau2 under no effect

Just a quick follow-up: It just occurred to me that I also hacked the
possibility to fix 'beta' values into rma.mv() at one point (had totally
forgotten about this). So this should also work:

dat$id <- 1:k
rma.mv(yi, vi, random = ~ 1 | id, data=dat, beta=0)

This is done in a more crude way, but should work at least for such a case
where there really is only an intercept term and we want to fix the
intercept to 0.

Best,
Wolfgang

-----Original Message-----
From: Filippo Gambarota <filippo.gambarota at gmail.com>
Sent: Sunday, March 17, 2024 08:00
To: Viechtbauer, Wolfgang (NP) <

wolfgang.viechtbauer at maastrichtuniversity.nl>

Cc: R Special Interest Group for Meta-Analysis <r-sig-meta-analysis at r-
project.org>
Subject: Re: [R-meta] estimating tau2 under no effect

Thank you Wolfgang! That is extremely helpful.

On Fri, 15 Mar 2024 at 15:38, Viechtbauer, Wolfgang (NP)
<mailto:wolfgang.viechtbauer at maastrichtuniversity.nl> wrote:
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 <mailto:

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 <mailto:r-sig-meta-analysis at r-project.org>
Cc: Filippo Gambarota <mailto: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 http://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: http://filippogambarota.xyz
Research Groups: Colab <http://colab.psy.unipd.it/>   Psicostat
<https://psicostat.dpss.psy.unipd.it/>