# Fitting model with unstructured covariance matrix with gls in R

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?

## 1 Answer

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.

• Many thanks for your help, Dimitris. That's unfortunate. Yes, for this specific toy example that is true, I could just split the data in two parts, but things are different in the actual data example I have. – user2960323 Jan 17 '19 at 11:40