Allowing correlation parameters in gls to depend on grouping factor

Question

Below there's a MWE of the issue I'm encountering. I'm using the Orthodont dataset from the nlme package which consists of 4 measurements for 27 children (16 girls, 11 boys). To model the correlation I'm using an unstructured covariance structure by specifying correlation = corSymm(form = ~1|Subject). I'm allowing for non-constant variance across different measurement occasions, but I would also like to allow the variance-covariance parameters to be different for boys and girls (e.g., because correlation between measurements is possibly higher/lower for boys than for girls). I can allow for this heterogeneity for the variance parameters by specifying weights = varIdent(form = ~1|age*Sex), but does anyone know how to allow for / specify this heterogeneity for the correlation parameters?

I know this is possible in SAS proc mixed by specifying the group option in the repeated statement (http://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_mixed_sect019.htm), but I haven't found a way to deal with this in R.

Many thanks in advance!

library(nlme)
head(Orthodont)
fit <- gls(distance ~ age * Sex, correlation = corSymm(form = ~1|Subject),
           weights = varIdent(form = ~1|age*Sex), data = Orthodont, na.action = na.exclude)

If you look at the description of form for ?corSymm, I'm not sure you can...as it says the following: When a grouping factor is present in form, the correlation structure is assumed to apply only to observations within the same grouping level; observations with different grouping levels are assumed to be uncorrelated. Maybe you can prespecify by using an initialized general correlation structure? The help function has some info on that, because I have no clue how to do that. — OFish

rcorty rcorty · Accepted Answer · 2015-03-31T17:12:07

I'm not completely sure this is what you're looking for, but give this a try:

fit <- gls(distance ~ age * Sex,
           correlation = corSymm(form = ~1|Subject/Sex),
           weights = varIdent(form = ~1|age*Sex), 
           data = Orthodont, na.action = na.exclude)

It fits a totally separate day-to-day correlation for males and females. If you run

summary(fit)

it shows this pretty clearly.

Allowing correlation parameters in gls to depend on grouping factor

1 Answers