[R-meta] estimating tau2 under no effect

[Viechtbauer, Wolfgang (NP)]

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 using gmail.com>
> Sent: Sunday, March 17, 2024 08:00
> To: Viechtbauer, Wolfgang (NP) <wolfgang.viechtbauer using maastrichtuniversity.nl>
> Cc: R Special Interest Group for Meta-Analysis <r-sig-meta-analysis using 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 using 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 using 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 using r-project.org>
> > Cc: Filippo Gambarota <mailto: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 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/>