Fitting model with unstructured covariance matrix with gls in R

Question

Consider the following dataset:

library(nlme)
library(lme4)

data<-sleepstudy
data<-data[data$Days<5,]
data$b<-sample(0:1,nrow(data), rep=T)

I am trying to replicate the results of this Stata code in R:

use data, clear
xtmixed Reaction Days##b || Subject: , nocon reml residuals(un, t(Days) by(b) )

This is basically a model for longitudinal data with unstructured mean and covariance matrix, and heteroscedasticity by a grouping factor b.

Since, due to the nocon option, there is no random effect, I thought the right function to use was gls. Currently I was fitting the following model:

fit <- gls(Reaction ~ Days*b , data = data, corr=corSymm(form=~1|Subject),
weights=varIdent(form=~1|Days*b), method="REML", na.action=na.omit)

However, looking at the output:

Correlation Structure: General
Formula: ~1 | Subject 
Parameter estimate(s):
Correlation: 
  1     2     3     4    
2 0.859                  
3 0.703 0.876            
4 0.732 0.870 0.951      
5 0.625 0.743 0.815 0.931
Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | Days * b 
 Parameter estimates:
      0*1       1*0       2*1       3*1       4*0       0*0       1*1           2*0       4*1 
1.0000000 1.3521764 0.5289573 1.5353564 1.2546706 1.2076804 0.9908704 1.1038659 1.5339758 
      3*0 
1.0425169

this seems to correctly give a different variance for each level of b and time point, but a common correlation for both levels of grouping factor b. I tried with corr=corSymm(form=~1|Subject/b) but this did not make a difference.

Is it possible to allow for grouping level b-specific correlations?

Answer 1

As far as I know, the correlation structures that are available for gls(), via its correlation argument do not allow to estimate correlation parameters conditional on covariates.

However, for the specific model you want to fit, and because both the mean and variance-covariance structure are different for the two levels of b, you could split your dataset into two parts according to the levels of b and estimate the model in each part.