0
$\begingroup$

For the longitudinal data provided below, we have the following variables: the response variable y, the time variable 'week', 'grp' (with two levels: grp1 and grp2), and 'subject'. My objective is to construct a model to evaluate the mean response at each week for each group and to incorporate an appropriate correlation structure for the repeated measures, such as the unstructured covariance. How can this be implemented in R? An example of the model specification in SAS is provided below.


proc mixed data = df; class grp week; model y = grp week grp*week / ddfm = km; repeated week / subject type = un r; run;

An example data 'df' has the following format (actual data can contain many more observations):

df

A tibble: 12 × 4

subject grp week y 1 1 grp1 2 0.692 2 1 grp1 4 1.63 3 1 grp1 6 2.46 4 3 grp1 2 0.753 5 3 grp1 4 1.48 6 3 grp1 6 1.94 7 2 grp2 2 0.356 8 2 grp2 4 0.637 9 2 grp2 6 1.03 10 4 grp2 2 0.596 11 4 grp2 4 0.907 12 4 grp2 6 1.48

Thank you very much!

$\endgroup$
8
  • $\begingroup$ Thanks so much for the answer. I usually use lme4 package for fitting mixed models in R. But I struggled to find a way to specify unstructured correlation for repeated measures using lmer. Is there a recommended reading that give a clear and concise summary on how to specify difference mixed effect models using nlme? $\endgroup$
    – user13154
    Mar 24 at 13:59
  • $\begingroup$ One more question, before specifying the gls model as you did below, besides setting week as a factor, it is necessary to set id as a factor as well, right? $\endgroup$
    – user13154
    Mar 24 at 14:03
  • $\begingroup$ No, that is not needed afaik, you could try. $\endgroup$
    – BenP
    Mar 24 at 14:09
  • $\begingroup$ Btw I feel your question would be more interesting for readers if you would ask about how to choose a suitable structure, and which statistical measures exist, if any, that help one choosing. I.e. there are different patterns than only the unstructired. $\endgroup$
    – BenP
    Mar 24 at 14:16
  • $\begingroup$ ok. Agree that is more useful. Thanks a lot!! $\endgroup$
    – user13154
    Mar 24 at 14:17

1 Answer 1

1
$\begingroup$

Although the question is formulated as a programming/syntax problem I will answer it, because I'm not sure if you are aware of the many other options to model (co)variance matrices over "weeks". And also because I've struggled quite a while myself to figure out the (many) possibilities of package "nlme" and those of "gls", which is often a great alternative for random effect models, in case of longitudinal data.

Before giving the syntax, a "warning". Are you sure you need un unstructed covariance pattern, which is quite parameter-consuming? There are many other structures in the "gls" procedure from "nlme" which you could consider. Also you "should" study your correlations and variances over "weeks", how they differ, if these differences make sense to you. And for any "structure" you choose, compare the variances and correlations over the weeks as produced by the model with the ones you observed for your data. Do the model-produced variances and correlations make sense to you and are they in line with the observed ones? But you may know all this, why am I telling you this (probably because it's a rainy sunday morning).

It's been a while since I used SAS, but I'll give it a shot. You could use "gls" from package "nlme" to do this. Let's say your data are in data frame "da". Variable "week" should be an integer variable running from 1,2,3, ... etc. If not, make a new one, with name "weeknr"e.g. It's important that this new variable "weeknr" is integer valued and starting at value 1, not zero! For the interaction you must have a factor for week.

I would start with a model without "grp", to first study the variances and correlations over the weeks.

da$weekf <- factor(da$week)
mymodel <- gls(y ~ weekf, 
               correlation = corSymm(form = ~ weeknr|id),
               weights = varIdent(form = ~1|weeknr), data=da)  
summary(mymodel)

The summary shows the effects of the week-dummies, how they differ. To see the variances and correlations produced by the gls model, you can use a script, which I will put at the end of this answer. Before running your "full" model, I would first choose the right pattern for the variances / correlation over the weeks, not only look at the "unstructured" pattern. Use AIC and BIC to compare the different structures in terms of "fit", and use theory/common sense to judge if the model predicted vars and corrs are in line with your (theoretical) expectations.

Then to your "full" model:

mymodel <- gls(y ~ grp*weekf, 
               correlation = corSymm(form = ~ weeknr|id),
               weights = varIdent(form = ~1|weeknr), data=da)  
summary(mymodel)

To see the model-produced variances and correlations, here is a R script that you must run first, and next you can use:

corandcov.gls(mymodel)

That's it!

R script for shown variances and correlations over time for gls models.

# corandcov_gls
#
# Function found on the internet for showing (co)variances and 
# correlations for a generalised least squares (GLS) model 
# estimated with gls form nlme package.
# To call the function use:
#
# corandcov_gls(modelname)
#
# where "modelname" is the name of the GLS model you estimated.
#--------------------------------------------------------------

corandcov.gls <- function(glsob, cov = TRUE, ...) {
  
  corm <- nlme::corMatrix(glsob$modelStruct$corStruct)[[5]]
  
  varstruct <- glsob$modelStruct$varStruct
  
  if (!is.null(varstruct)) {
  varests <- stats::coef(varstruct, uncons = FALSE, allCoef = TRUE)
  covm <- corm * glsob$sigma ^ 2 * t(t(varests)) %*% t(varests)
  
  return(list(correlation_matrix = corm,
         variance_structure = varstruct,
         covariance_matrix = covm))
  }
  
  if (is.null(varstruct)) {
    # The homogeneous variance models CS, AR1 en Toep only have 1 variance
    # and no varStruct which equal NULL then.
    var <- getVarCov(glsob)[1,1]
    list(correlation_matrix = corm, variance = var)
  }
}
$\endgroup$

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