[R-sig-ME] reference category in varIdent function (once again)

[Sebastian Meyer]

I think it is unfortunate that varIdent() selects the reference category 
based on how the data is sorted (using the first level that occurs). 
There is also a corresponding bug report 
<https://bugs.r-project.org/show_bug.cgi?id=18505>.

Your example illustrates an additional issue with this behaviour, 
because lme() internally reorders the data by the grouping factor(s) 
before initialization, so the order of that factor ('Subject' in your 
example, with first level "M16", so male) will indirectly determine the 
reference level of varIdent(), regardless of how your data or strata 
levels are ordered originally.

It seems that only if there are more than two strata, the reference 
level for varIdent() can be chosen via a named initial parameter 
'value', leaving out the desired reference level, but I haven't tested 
if this works as intended with both gls() and lme(). It would then make 
sense to support such a choice also in the case of only two strata, so 
for a single parameter.

Best wishes,

	Sebastian Meyer


Am 20.12.24 um 20:39 schrieb James Pustejovsky:
> Greetings,
> 
> I'm trying to fit a linear mixed effect model using nlme::lme() that
> includes heteroskedastic variances by group, where the variance for one
> group is fixed to a specific value. I *think* I can do this using
> varIdent() and setting lmeControl(sigma = x) for a specified value x.
> However, I've run into what is evidently an old problem: that there doesn't
> appear to be any way to control the reference level of varIdent().
> 
> As was noted in a very old post to R-SIG-mixed-models (
> https://stat.ethz.ch/pipermail/r-sig-mixed-models/2013q3/020582.html), the
> reference level of varIdent() within an lme() call is set to the most
> common level of the grouping factor. For reasons that are too involved to
> get into, I need to set the reference level to something else (in fact, the
> least frequently occurring level), so it seems I'm stuck.
> 
> Curiously, the behavior of varIdent() within a gls() call is a bit
> different. As documented in this StackOverflow (
> https://stackoverflow.com/questions/75929847/reference-category-in-varident-function-nlme-package-depends-on-data-order),
> the reference level used within gls() depends on the sort order of the
> data. I've included below a minimal working example that demonstrates the
> difference in behavior.
> 
> Does anyone know of any work-around for setting the reference level of
> varIdent() within an lme() call? I'd be grateful for any leads/suggestions
> on what I could do (up to and including to just use a different package
> already).
> 
> Thanks,
> James
> 
> 
> ----------------------------------------------------------------------------
> library(nlme)
> 
> # Make reversed factor
> table(Orthodont$Sex)
> Orthodont$Sex_rev <- Orthodont$Sex
> levels(Orthodont$Sex_rev) <- rev(levels(Orthodont$Sex))
> table(Orthodont$Sex_rev)
> 
> Ortho_male_first <- Orthodont[order(Orthodont$Sex),]
> Ortho_female_first <- Orthodont[order(Orthodont$Sex, decreasing = TRUE),]
> 
> # Fit gls models with heteroskedastic variance
> gls1 <- gls(
>    distance ~ age + Sex,
>    weights = varIdent(form = ~ 1 | Sex),
>    data = Ortho_male_first
> )
> gls2 <- update(gls1, data = Ortho_female_first)
> gls3 <- update(gls1, weights = varIdent(form = ~ 1 |Sex_rev))
> gls4 <- update(gls3, data = Ortho_female_first)
> 
> # Parameterization depends on sort order
> coef(gls1$modelStruct$varStruct)
> coef(gls2$modelStruct$varStruct)
> # But not on factor level order
> coef(gls3$modelStruct$varStruct)
> coef(gls4$modelStruct$varStruct)
> 
> # Fit lme models with heteroskedastic level-1 variance
> 
> lme1 <- lme(
>    distance ~ age + Sex,
>    random = ~ 1,
>    weights = varIdent(form = ~ 1 | Sex),
>    data = Ortho_male_first
> )
> lme2 <- update(lme1, data = Ortho_female_first)
> lme3 <- update(lme1, weights = varIdent(form = ~ 1 |Sex_rev))
> lme4 <- update(lme3, data = Ortho_female_first)
> 
> # Parameterization depends neither on sort order nor on factor level order
> coef(lme1$modelStruct$varStruct)
> coef(lme2$modelStruct$varStruct)
> coef(lme3$modelStruct$varStruct)
> coef(lme4$modelStruct$varStruct)
> 
> 	[[alternative HTML version deleted]]
> 
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models