[R-meta] estimating tau2 under no effect

[Filippo Gambarota]

Thank you Wolfgang! That is extremely helpful.

On Fri, 15 Mar 2024 at 15:38, Viechtbauer, Wolfgang (NP) <
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 <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/>
>


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

	[[alternative HTML version deleted]]