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

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