[R-sig-ME] [R] mixed model fitting between R and SAS

[ONKELINX, Thierry]

John,

Please keep sending the mails to the mixed models list which is more appropriate than R-help.

Both models have a different interpretation.
the nested model is indeed
y_ijk=fixed_effect + b_i + b_ij + e_ijk
the the one with the random slopes is rather
y_ijk=fixed_effect + b_1i + b_2i + e_ijk
with b_1i the effect of location 1 for pid i and b_2i the effect of location 2 for pid i.
The fits will be mathematically equal, but the parametrisation is different.

Best regards,

Thierry

> -----Oorspronkelijk bericht-----
> Van: array chip [mailto:arrayprofile at yahoo.com]
> Verzonden: dinsdag 9 augustus 2011 0:50
> Aan: ONKELINX, Thierry; r-help at stat.math.ethz.ch
> Onderwerp: Re: [R] mixed model fitting between R and SAS
> 
> 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