[R-sig-ME] lme: count the number extra parameters estimated for variance or covariances

[James Pustejovsky]

Hi Simon,

Yes, the estimated residual variance among males is as you've calculated.
Using lmeInfo again:

VC <- extract_varcomp(fit)
with(VC, var_params * sigma_sq)

Female is the reference level here, so the residual variance there is
simply VC$sigma_sq.

James

On Sun, Nov 1, 2020 at 12:19 PM Simon Harmel <sim.harmel using gmail.com> wrote:

> By the way, Wolfgang, for my model above the variance function returns
> 1.047655 for males. I know this is a multiplicative constant, but should
> I multiply this constant by the estimated residual variance for the entire
> model to obtain the estimated residual variance for males? That is:
>
> sigma(fit)^2 * 1.047655 ## > [1] 38.91265 ## what would be the residual
> variance for females?
>
> ===============
> library(nlme)
>
> hsb <- read.csv('
> https://raw.githubusercontent.com/rnorouzian/e/master/hsb.csv')
> hsb$female <- as.factor(hsb$female)
>
> fit <- lme(math ~ female, random = ~ 1|sch.id, data = hsb, weights =
> varIdent(form = ~1 |female))
>
> Variance function:
> Formula: ~1 | female
>  Parameter estimates:
>        1         0
> 1.000000 1.047655
>
>
>
> On Sun, Nov 1, 2020 at 11:15 AM Simon Harmel <sim.harmel using gmail.com> wrote:
>
>> Great, many thanks to all!
>>
>> On Sun, Nov 1, 2020 at 11:14 AM Viechtbauer, Wolfgang (SP) <
>> wolfgang.viechtbauer using maastrichtuniversity.nl> wrote:
>>
>>> Ah, very cool. Did not know about the lmeInfo package.
>>>
>>> And yes, one can think of the parameter/coefficient for females as a
>>> binary predictor that allows the error variance to differ for females from
>>> that of the males.
>>>
>>> Best,
>>> Wolfgang
>>>
>>> >-----Original Message-----
>>> >From: Simon Harmel [mailto:sim.harmel using gmail.com]
>>> >Sent: Sunday, 01 November, 2020 18:06
>>> >To: James Pustejovsky
>>> >Cc: Viechtbauer, Wolfgang (SP); Harold Doran; r-sig-mixed-models
>>> >Subject: Re: [R-sig-ME] lme: count the number extra parameters
>>> estimated for
>>> >variance or covariances
>>> >
>>> >Thanks, James!
>>> >
>>> >On Sun, Nov 1, 2020 at 11:03 AM James Pustejovsky <jepusto using gmail.com>
>>> wrote:
>>> >Simon,
>>> >
>>> >Here is a quick way to accomplish the same thing that Wolfgang
>>> demonstrated,
>>> >using the lmeInfo package:
>>> >VC <- lmeInfo::extract_varcomp(fit)   # get all the variance components
>>> >lengths(VC)                                        # count the number of
>>> >estimated parameters in each component
>>> >sum(lengths(VC))                               # total number of
>>> variance
>>> >component parameters
>>> >
>>> >Kind Regards,
>>> >James
>>> >
>>> >On Sun, Nov 1, 2020 at 10:45 AM Simon Harmel <sim.harmel using gmail.com>
>>> wrote:
>>> >Thank you, Wolfgang (sorry for misspelling). So, by: " The coefficient
>>> for
>>> >the females is of course [one additional parameter]" you mean for the
>>> >variance coefficient of `female == 1` as a binary predictor, right?
>>> >
>>> >On Sun, Nov 1, 2020 at 10:31 AM Viechtbauer, Wolfgang (SP) <
>>> >wolfgang.viechtbauer using maastrichtuniversity.nl> wrote:
>>> >
>>> >> I meant that the 1 for males is not an estimated parameter. The
>>> >> coefficient for the females is of course (and hence one additional
>>> >> parameter). Apologies for the confusion.
>>> >>
>>> >> For correlation structures, there will indeed be a 'corStruct' element
>>> >> under 'modelStruct'.
>>> >>
>>> >> Best,
>>> >> Wolfgang
>>> >>
>>> >> >-----Original Message-----
>>> >> >From: Simon Harmel [mailto:sim.harmel using gmail.com]
>>> >> >Sent: Sunday, 01 November, 2020 17:15
>>> >> >To: Viechtbauer, Wolfgang (SP)
>>> >> >Cc: Harold Doran; r-sig-mixed-models
>>> >> >Subject: Re: [R-sig-ME] lme: count the number extra parameters
>>> estimated
>>> >> for
>>> >> >variance or covariances
>>> >> >
>>> >> >Thank you Wolfang. That was exactly what I was looking for. If an
>>> lme()
>>> >> >model uses  `correlation = corAR1()`, then I'm assuming something new
>>> >> >will appear for the question mark in the following:
>>> `m2$modelStruct$???`,
>>> >> >right?
>>> >> >
>>> >> >Wolfang, on the one the hand you mentioned: "you will also get the
>>> >> >coefficient (= 1) for the males. But that is not actually an
>>> estimated
>>> >> >parameter",
>>> >> >
>>> >> >On the other hand you mentioned: "And these *are* parameters
>>> (besides the
>>> >> >fixed effects and the vars/covs of the random effects)."
>>> >> >
>>> >> >Multiple software I used show that my model with `varIdent(form = ~1
>>> >> >|female)` estimates one additional parameter compared to a
>>> corresponding
>>> >> >model without `varIdent(form = ~1 |female)`.
>>> >> >
>>> >> >Books (e.g., Mixed Effects Models and Extensions in Ecology with R by
>>> >Zuur
>>> >> >et al, 2009; Pinheiro & Bates, 2000, p. 209) also clearly mention
>>> >> >`varIdent(form = ~1 |female)` estimates one more parameter.
>>> >> >
>>> >> >Would you please clarify?
>>> >> >
>>> >> >On Sun, Nov 1, 2020 at 9:44 AM Viechtbauer, Wolfgang (SP)
>>> >> ><wolfgang.viechtbauer using maastrichtuniversity.nl> wrote:
>>> >> >Dear Simon,
>>> >> >
>>> >> >For variance structures, you can use:
>>> >> >
>>> >> >coef(fit$modelStruct$varStruct)
>>> >> >
>>> >> >That will give you the parameter estimates involved in the variance
>>> >> >structure (in their constrained form as used during the
>>> optimization).
>>> >> With:
>>> >> >
>>> >> >coef(fit$modelStruct$varStruct, unconstrained=FALSE)
>>> >> >
>>> >> >you can get the unconstrained estimates. Only coefficients that are
>>> >> actually
>>> >> >estimated are returned by default. With:
>>> >> >
>>> >> >coef(fit$modelStruct$varStruct, unconstrained=FALSE, allCoef=TRUE)
>>> >> >
>>> >> >you will also get the coefficient (= 1) for the males. But that is
>>> not
>>> >> >actually an estimated parameter. For more details, see:
>>> >> >
>>> >> >help(coef.varFunc)
>>> >> >
>>> >> >And these *are* parameters (besides the fixed effects and the
>>> vars/covs
>>> >of
>>> >> >the random effects).
>>> >> >
>>> >> >Best,
>>> >> >Wolfgang
>>> >> >
>>> >> >>-----Original Message-----
>>> >> >>From: R-sig-mixed-models [mailto:
>>> >> r-sig-mixed-models-bounces using r-project.org]
>>> >> >>On Behalf Of Harold Doran
>>> >> >>Sent: Sunday, 01 November, 2020 16:26
>>> >> >>To: Simon Harmel
>>> >> >>Cc: r-sig-mixed-models
>>> >> >>Subject: Re: [R-sig-ME] lme: count the number extra parameters
>>> estimated
>>> >> >for
>>> >> >>variance or covariances
>>> >> >>
>>> >> >>I think you need to understand what a reproducible example is
>>> intended
>>> >to
>>> >> >>do. Your data estimates a model and yields a fiitted model object.
>>> What
>>> >> >>parameter from that object using an extractor are you intending to
>>> find?
>>> >> >>
>>> >> >>For example, a well posed question would be something like. I want
>>> to
>>> >> >>extract the fixed effects from a fitted model object. How do I get
>>> them?
>>> >> >>
>>> >> >>To say I want the “parameters estimated for modeling residual
>>> variances”
>>> >> >etc
>>> >> >>makes no sense. The parameters of a mixed model are the fixed
>>> effects
>>> >and
>>> >> >>the marginal variances (and covariances) of the random effects.
>>> >> >>
>>> >> >>So, specifically what parameters do you think exist in a model that
>>> you
>>> >> >>want?
>>> >> >>
>>> >> >>From: Simon Harmel <sim.harmel using gmail.com>
>>> >> >>Sent: Sunday, November 1, 2020 9:58 AM
>>> >> >>To: Harold Doran <harold.doran using cambiumassessment.com>
>>> >> >>Cc: r-sig-mixed-models <r-sig-mixed-models using r-project.org>
>>> >> >>Subject: Re: [R-sig-ME] lme: count the number extra parameters
>>> estimated
>>> >> >for
>>> >> >>variance or covariances
>>> >> >>
>>> >> >>Dear Harold,
>>> >> >>
>>> >> >>My question "specifically" is:  is there a way (e.g., via an
>>> extractor
>>> >> >>function) to obtain parameters estimated for modeling residual
>>> variances
>>> >> or
>>> >> >>covariances from an "lme" model?
>>> >> >>
>>> >> >>For concreteness, please consider the reproducible model I provided
>>> in
>>> >my
>>> >> >>original post in which a variance function has been used.
>>> >> >>
>>> >> >>Thanks,
>>> >> >>
>>> >> >>On Sun, Nov 1, 2020, 4:43 AM Harold Doran
>>> >> >><harold.doran using cambiumassessment.com<mailto:
>>> >> harold.doran using cambiumassessment.c
>>> >> >o
>>> >> >>m>> wrote:
>>> >> >>In order to answer that you need to specify what "thing" you want.
>>> The
>>> >> >>object itself has many things and there are extractor functions to
>>> grab
>>> >> >many
>>> >> >>of them. I say "thing" because the *parameters* of a mixed model
>>> are the
>>> >> >>fixed effects and the variance components. Random effects etc are
>>> not
>>> >> >>parameters of a mixed model.
>>> >> >>
>>> >> >>You can always look at the structure of a fitted model object in R
>>> to
>>> >see
>>> >> >>what things are generally in it.
>>> >> >>
>>> >> >>-----Original Message-----
>>> >> >>From: R-sig-mixed-models <r-sig-mixed-models-bounces using r-
>>> >> >project.org<mailto:r-
>>> >> >>sig-mixed-models-bounces using r-project.org>> On Behalf Of Simon Harmel
>>> >> >>Sent: Sunday, November 1, 2020 4:02 AM
>>> >> >>To: r-sig-mixed-models <r-sig-mixed-models using r-project.org<mailto:
>>> r-sig-
>>> >> >mixed-
>>> >> >>models using r-project.org>>
>>> >> >>Subject: [R-sig-ME] lme: count the number extra parameters
>>> estimated for
>>> >> >>variance or covariances
>>> >> >>
>>> >> >>External email alert: Be wary of links & attachments.
>>> >> >>
>>> >> >>Hello All,
>>> >> >>
>>> >> >>In addition to fixed and random effects, is there a way to extract
>>> how
>>> >> many
>>> >> >>other parameters (for modeling residual variances or covariances) an
>>> >> >"lme()"
>>> >> >>object has estimated?
>>> >> >>
>>> >> >>Here is a reproducible example:
>>> >> >>
>>> >> >>library(nlme)
>>> >> >>
>>> >> >>hsb <- read.csv('
>>> >> >>https://raw.githubusercontent.com/rnorouzian/e/master/hsb.csv')
>>> >> >>hsb$female <- as.factor(hsb$female)
>>> >> >>
>>> >> >>fit <- lme(math ~ female, random = ~ 1|sch.id<http://sch.id>, data
>>> =
>>> >> hsb,
>>> >> >>weights = varIdent(form = ~1 |female))
>>>
>>

	[[alternative HTML version deleted]]