[R-meta] How to interpret sigma(model)^2 in metafor

Viechtbauer, Wolfgang (NP) wo||g@ng@v|echtb@uer @end|ng |rom m@@@tr|chtun|ver@|ty@n|
Tue Jan 24 09:21:27 CET 2023


Dear Yuhang,

sigma() is a generic function:

> sigma                                                                                                                 
function (object, ...)                                                                                                  
UseMethod("sigma")                                                                                                      
<bytecode: 0x2e1a3c8>                                                                                                   
<environment: namespace:stats>   

So, when calling sigma() on an object from the metafor package, the 'S3 dispatch mechanism' will first check if there is a method for the type (i.e., class) of object that you are passing to the sigma() function.

> library(metafor)
> methods(sigma)
[1] sigma.default* sigma.gls*     sigma.lmList*  sigma.lme*     sigma.mlm*
see '?methods' for accessing help and source code

Since there is none (there is no sigma.rma() function or anything like it), it will call sigma.default(). So let's look at what that does:

> sigma.default
Error: object 'sigma.default' not found

Hmmm, why can't we look at the code for this function? Note the * after sigma.default -- this indicates that the method definition is not exported. But we can still look at this with:

> getAnywhere(sigma.default)
A single object matching 'sigma.default' was found
It was found in the following places
  registered S3 method for sigma from namespace stats
  namespace:stats
with value

function (object, use.fallback = TRUE, ...)
sqrt(deviance(object, ...)/(nobs(object, use.fallback = use.fallback) -
    sum(!is.na(coef(object)))))
<bytecode: 0x9806a08>
<environment: namespace:stats>

(or, if you would happen to know that this function comes from the stats package, you could use stats:::sigma.default).

So, we can see what is happening. In essence:

dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)
res <- rma(yi, vi, data=dat)
sigma(res)
sqrt(deviance(res)/(nobs(res) - sum(!is.na(coef(res)))))

This has, as far as I am concerned, no logical meaning for rma objects.

Note that I try to be very explicit in the documentation what kind of methods are available (and meaningful) for a given object in the metafor package:

https://wviechtb.github.io/metafor/reference/rma.uni.html#methods

sigma() is not listed there. I cannot prevent default methods from being called unless I would actually put a sigma.rma() method into the metafor package. I have actually considered this, but I don't have a good idea what meaningful result this should return.

In any case, I hope this provides you with some idea how you can dig into the code (and the mechanisms of how it is being called) in general.

Best,
Wolfgang

>Hello Colleagues,
>
>By habit, I always check the variance of the residuals of my ordinary
>regression models using: sigma(model)^2, which is also printed in the
>output.
>
>I know that the variance of the residuals in meta-regression is not
>estimated but rather taken as being known and fixed by virtue of user's
>supplying the 'vi' or 'V' to functions such as rma.uni() and
>rma.mv().
>
>So, I was wondering what is the interpretation of sigma(rma.uni_model)^2
>and sigma(rma.mv_model)^2 and how they connect to the user-supplied
>'vi'.
>
>Thank you very much for your time.
>
>Best,
>Yuhang



More information about the R-sig-meta-analysis mailing list