[R] mixed model fitting between R and SAS

[array chip]

Dear Thierry,

Thanks a lot for pointing me to the right direction! I still have some questions, really appreciate if you could provide any help:


Is there a relationship between the nested mixed model that I used vs. the model that you gave (using biop location as random slope of pid), i.e. 


lme(y~age + time * trt, random=~1|pid/biopsy.site, data = test, correlation=corAR1())

vs.
lme(y~age + time * trt, random=~0  + biopsy.site|pid, data = test,  correlation=corAR1())

In my nested mixed model, a variance of pid and a variance of biopsy.site within pid will be estimated. In your mixed model, there is no variance estimated at the pid level , instead variance for each biopsy.site is given. I thought the statistical model was always the same for both mixed models:

y_ijk=fixed_effect + b_i + b_ij + e_ijk

where b_i is the random effect for pid, and v_ij is the random effect for biopsy.site within pid. I thought that the difference is in my nested mixed model, b_ij is independent of each other and has the same variance, whereas in your mixed model, b_ij is modeled by the pdClasses() chosen. If my thought was correct, then my model should be the same as 


lme(y~age + time * trt, random=list(pid=pdIdent(~0  + biopsy.site)), data = test,  correlation=corAR1())

But they are not the same. What did I misunderstand here?

Many thanks

John




----- Original Message -----
From: "ONKELINX, Thierry" <Thierry.ONKELINX at inbo.be>
To: array chip <arrayprofile at yahoo.com>; "r-sig-mixed-models at r-project.org" <r-sig-mixed-models at r-project.org>
Cc: 
Sent: Monday, August 8, 2011 1:19 AM
Subject: RE: [R] mixed model fitting between R and SAS

Please don't cross-post.

- corAR1() models the correlation between the residuals of the two time points. 
- if you want a specific correlation structure for biopsy locations the you must use on of the pdClasses() and use biopsy location as random slope of pid rather than random effect nested in pid.

#basic structure = positive definitive symmetrical variance/covariance matrix
lme(y~age + time * trt, random=~0  + biopsy.site|pid, data = test,  correlation=corAR1(~time))
#no correlation between biopsy location and different variance
lme(y~age + time * trt, random=list(pid  = pdDiag(~0  + biopsy.site), data = test,  correlation=corAR1(~time))
#no correlation between biopsy location and equal variance
lme(y~age + time * trt, random=list(pid  = pdIdent(~0  + biopsy.site), data = test,  correlation=corAR1(~time))

Note that since you have only two biopsy locations there will be no difference between pdSymm (the default) and pdCompSymm

Best regards,

Thierry
----------------------------------------------------------------------------
ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek
team Biometrie & Kwaliteitszorg
Gaverstraat 4
9500 Geraardsbergen
Belgium

Research Institute for Nature and Forest
team Biometrics & Quality Assurance
Gaverstraat 4
9500 Geraardsbergen
Belgium

tel. + 32 54/436 185
Thierry.Onkelinx at inbo.be
www.inbo.be

To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of.
~ Sir Ronald Aylmer Fisher

The plural of anecdote is not data.
~ Roger Brinner

The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey

> -----Oorspronkelijk bericht-----
> Van: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
> Namens array chip
> Verzonden: maandag 8 augustus 2011 8:48
> Aan: r-sig-mixed-models at r-project.org
> CC: r-help
> Onderwerp: [R] mixed model fitting between R and SAS
> 
> Hi al,
> 
> I have a dataset (see attached), which basically involves 4 treatments for a
> chemotherapy drug. Samples were taken from 2 biopsy locations, and biopsy
> were taken at 2 time points. So each subject has 4 data points (from 2 biopsy
> locations and 2 time points). The objective is to study treatment difference.
> 
> I used lme to fit a mixed model that uses "biopsy.site nested within pid" as a
> random term, and used corAR1() as the correlation structure for between the 2
> time points:
> 
> 
> library(nlme)
> 
> test<-read.table("test.txt",sep='\t',header=T,row.names=1)
> fit<-lme(y~age + time * trt, random=~1|pid/biopsy.site, data = test,
> correlation=corAR1())
> 
> First, by above model specification, corAR1() is used for the correlation between
> the 2 time points; what is the correlation structure implicitly used for between
> biopsy locations? How do I specify a particular correlation structure for between
> biopsy locations in this situation?
> 
> Second, does anyone know how to write the above mixed model in SAS? One of
> my colleagues wrote the following, but it gave me different results:
> 
> proc mixed data=test;
> 
> class time trt pid biopsysite;
> model y=age time trt time*trt;
> random biopsysite
> repeated pid / type=ar(1)
> run;
> 
> Is there anyone familiar with SAS and know if the above SAS code does what the
> R code does?
> 
> Many thanks
> 
> John