[R-meta] estimating tau2 under no effect

Tue, Mar 19, 2024 8:05 AM

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/>

[R-meta] estimating tau2 under no effect

Thread (5 messages)