Consider the following dataset:


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):
  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 

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?


migrated from stackoverflow.com Jan 14 at 13:33

This question came from our site for professional and enthusiast programmers.


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.

  • $\begingroup$ 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. $\endgroup$ – user2960323 Jan 17 at 11:40

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.