[R-meta] estimating tau2 under no effect

[Filippo Gambarota]

Thanks! This approach is even better given that I can use rma.mv with the
standard syntax!

On Tue, 19 Mar 2024 at 16:05, Viechtbauer, Wolfgang (NP) <
wolfgang.viechtbauer using maastrichtuniversity.nl> wrote:

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


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