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

[Viechtbauer, Wolfgang (NP)]

The residuals are not the same as the sampling errors. Yes, the sampling errors have variance v_ij, but the *observed residuals* have a different variance which depends on the model being fitted, since Var(y_ij - fitted_ij) = Var(y_ij) + Var(fitted_ij) - 2*Cov(yij,fitted_ij).

Best,
Wolfgang

>-----Original Message-----
>From: Yuhang Hu [mailto:yh342 using nau.edu]
>Sent: Wednesday, 25 January, 2023 19:13
>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 there is a gap in my understanding, then, regarding
>the difference between sampling errors which I think for the current model (i.e.,
>yij = B0 + B1*yearij + uj + vij + eij) are defined as:
>
>eij = yij - theta_ij (true effect_ij), with Var(eij) = dat$vi for each eij ~
>Normal[0, dat$vi].
>
>and residuals i.e., Var(yij - fitted(yij)) .
>
>I thought that the variances of the sampling errors (i.e., dat$vi) which indicate
>that effect sizes collected from the literature are estimated with error
>eventually are the same as the Var(yij - fitted(yij)) thinking that fitted(yij)
>is taken as theta_ij.
>
>Best,
>Yuhang
>
>On Wed, Jan 25, 2023 at 9:48 AM Viechtbauer, Wolfgang (NP)
><wolfgang.viechtbauer using maastrichtuniversity.nl> wrote:
>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