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

[Simon Harmel]

Thank you James. Both SPSS and the HLM software estimate the variance for
males to be ~40.

But that aside, given the larger estimated variance for males (coded 0), I
thought I should see a "wider" distribution for males compared to females.
I don't see that to be the case, though (see ggplot below)?

hsb <- read.csv('
https://raw.githubusercontent.com/rnorouzian/e/master/hsb.csv')
hsb$female <- as.factor(hsb$female)

ggplot(hsb) + aes(x = math, color = female) + facet_wrap(~female, scales =
"free")+
  geom_histogram()

Of course, when I plot the "female" variable (a factor of 1s & 0s) against
the residuals of the model I see that the model has correctly worked
showing males are more scattered than females in their residual.

hsb <- within(hsb, resid <- resid(fit))
ggplot(hsb) + aes(x = female, y = resid) + facet_wrap(~female)+geom_point()

On Sun, Nov 1, 2020 at 3:41 PM James Pustejovsky <jepusto using gmail.com> wrote:

> 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]]