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

[Viechtbauer, Wolfgang (NP)]

They are not quite the same. dat$vi are not the sampling variances of the residuals, but of the sampling errors. Those are indeed assumed to be known and fixed. However, the sampling variances of the residuals depend on the model you are fitting.

Best,
Wolfgang

>-----Original Message-----
>From: Yuhang Hu [mailto:yh342 using nau.edu]
>Sent: Wednesday, 25 January, 2023 17:35
>To: Viechtbauer, Wolfgang (NP)
>Cc: r-sig-meta-analysis using r-project.org
>Subject: Re: Re: [R-meta] How to interpret sigma(model)^2 in metafor
>
>Thank you, Wolfgang. I guess my expectation was that `rstandard(res)$se^2` should
>be equal to `dat$vi` since the model takes `dat$vi` as given (known and fixed)
>for the sampling distribution of each residual in each row (e_ij ~ Normal[0,
>dat$vi])?
>
>Thank you very much,
>Yuhang
>
>On Wed, Jan 25, 2023 at 1:40 AM Viechtbauer, Wolfgang (NP)
><wolfgang.viechtbauer using maastrichtuniversity.nl> wrote:
>Depends on what you mean by 'variance of the residuals'. With:
>
>rstandard(res)$se^2
>
>you can obtain the *sampling variance* of the residuals. Note that each residual
>has its own sampling variance.
>
>If you just want the *sample variance* of the residuals, then
>
>var(resid(res))
>
>would give you that.
>
>Best,
>Wolfgang
>
>>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 using 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