Re: [R] Repeated measures using lme

From: Douglas Bates <bates_at_stat.wisc.edu>
Date: Fri, 07 Mar 2008 08:17:02 -0600

On Thu, Mar 6, 2008 at 3:51 AM, Tom Cohen <tom.cohen78_at_yahoo.se> wrote:
> Dear list,

> I am trying to do a repeated analysis using lme in R and a little bit unsure if I have set up the right statement. The problem is the IL6 (interleukin 6) was measured 5 times on each individual in each of 6 companies. The hypotheses are to see whether there is a relationship between IL6 and the total dust in each of the companies and if there is any change in IL6 across time points. So the fixed effects are total dust, company and random effects is individual. The model would be like this

> lme(IL6~ dust + time*company, random=~time|individual,
> correlation=corAR1(form=~time|individual), data=dat)

> with time as a repeated measure.

> the analysis in SAS would be

> proc mixed data=dat;
> class time individual company;
> model IL6=dust company time company*time;
> repeated time/Sub=individual(company) type=AR(1) r rcorr;
> random individual;
> lsmeans company time company*time/slice=time;
> run;

I'm not sure that those are equivalent specifications. If I read the SAS code correctly (and I don't have a lot of experience with SAS) the equivalent call to lme would have random = ~ 1 | individual

I would be more inclined to start with a model that does have a random effect for time but does not have the additional correlation structure. In lme this would be as in your specification but omitting the correlation argument. In lmer it would be

lmer(IL6 ~ dust + time * company + (time | individual), data)

In general it is difficult to estimate both a random effect for the slope w.r.t. time and an AR(1) correlation pattern over time within individual. The motivation for these two types of models is not identical but the end results in terms of the marginal covariance structure of the responses are similar. Generally it is better to stick with just one such correlation-inducing term, as you have done in the SAS specification.

Before fitting models to such data I recommend extensive plotting to see what types of models are indicated. With longitudinal data collected on individuals I would always begin with

library(lattice)
xyplot(IL6 ~ time | individual, data)

and probably enhance that by taking account of the grouping of individuals within company by grouping the panels according to company and ordering according to some characteristic of the individual within company.

Let me take the opportunity to plug Deepayan Sarkar's book "Lattice: Multivariate Data Visualization with R" due for release "real soon" and the web site http://lmdvr.r-forge.r-project.org/ for the book. It is a great experience just to browse the figures for each chapter at the web site.

I enclose a plot of the sort that I describe. In this case there is a single group of subjects.

> Am I writing the right code in R that would give me the same results if doing the analysis in SAS. Also is there any command in R that does the same thing as SLICE in SAS does, to test when in time there is difference between companies?
>
> Thanks for any help,
> Tom
>
>
> ---------------------------------
> Sök efter kärleken!
>
> [[alternative HTML version deleted]]
>
>
> ______________________________________________
> R-help_at_r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
>



R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Received on Fri 07 Mar 2008 - 14:20:39 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Fri 07 Mar 2008 - 15:30:20 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive