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

Tue, Jan 24, 2023 12:21 AM #

Dear Yuhang,

sigma() is a generic function:

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.

[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:

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:

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

Yuhang Hu

Tue, Jan 24, 2023 6:23 PM #

Thank you so much Wolfgang, for your response.

On the same note, is there a way to extract the variance of the residuals
for the model, say from a fitted model like below?

dat <- dat.konstantopoulos2011
res <- rma.mv(yi ~ I(year-mean(year)), vi, random = ~ 1 |district/study,
data= dat)

Thank you.
Yuhang


On Tue, Jan 24, 2023 at 1:21 AM Viechtbauer, Wolfgang (NP) <

wolfgang.viechtbauer at maastrichtuniversity.nl> wrote:

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