[R-sig-ME] [R] mixed model fitting between R and SAS
array chip
arrayprofile at yahoo.com
Tue Aug 9 18:46:40 CEST 2011
Dear Thierry,
Thanks very much again for your input! Pardon me that I still got 2 questions:
First regarding the argument "correlation=corAR1()", how does model correctly know the grouping structure to apply this particular correlation structure appropriately to observations?
With nested model specification (random=~1|pid/biopsy.site), I understand it's straightforward to figure out that it applies to the 2 time points within the same biopsy.site within the same subject.
But with random slope model specification (e.g. random=list(pid=pdIdent(~0 + biopsy.site))), how does the model figure out what observations it applies the correlation structure to? The random term only indicates that "pid' is the grouping variable, will it also add in the random slope variable as one of the grouping variable as well?
Also the "time" variable in the data is a factor, not a continuous variable, so I can't use corAR1(~time) syntax, just use corAR1() is ok?
The second question is when you say the 2 models are mathematically equal but different parameterisation, which particular pdClasses structure are you referring to, pdIdent()?
Actually, I tried all of pdClasses, but got different results by comparing fixed effect estimates:
> fixef(lme(y~age + time * trt, random=~1|pid/biopsy.site, data = test,correlation=corAR1()))
(Intercept) age timeDAY 8 trtB trtC
-0.0664413311 0.0001562844 0.1041869499 -0.0198187975 0.2225816296
trtD timeDAY 8:trtB timeDAY 8:trtC timeDAY 8:trtD
0.0600800217 -0.0691869499 -0.2417960052 -0.0906501595
> fixef(lme(y~age + time * trt, random=list(pid=pdIdent(~0+biopsy.site)), data = test,correlation=corAR1()))
(Intercept) age timeDAY 8 trtB trtC
-7.329312e-02 3.404107e-05 1.356288e-01 2.397692e-02 2.435217e-01
trtD timeDAY 8:trtB timeDAY 8:trtC timeDAY 8:trtD
6.953377e-02 -1.057646e-01 -2.983250e-01 -1.256355e-01
> fixef(lme(y~age + time * trt, random=list(pid=pdCompSymm(~0+biopsy.site)), data = test,correlation=corAR1()))
(Intercept) age timeDAY 8 trtB trtC
-7.370652e-02 2.550198e-05 1.371592e-01 2.763169e-02 2.446715e-01
trtD timeDAY 8:trtB timeDAY 8:trtC timeDAY 8:trtD
7.000318e-02 -1.075655e-01 -3.009948e-01 -1.271838e-01
> fixef(lme(y~age + time * trt, random=list(pid=pdSymm(~0+biopsy.site)), data = test,correlation=corAR1()))
(Intercept) age timeDAY 8 trtB trtC
-0.0855931906 0.0002516472 0.1364347017 0.0105537241 0.2457796991
trtD timeDAY 8:trtB timeDAY 8:trtC timeDAY 8:trtD
0.0709443633 -0.1120889667 -0.3029470038 -0.1273406900
> fixef(lme(y~age + time * trt, random=list(pid=pdDiag(~0+biopsy.site)), data = test,correlation=corAR1()))
(Intercept) age timeDAY 8 trtB trtC
-0.0854646844 0.0002610673 0.1351021118 0.0067663008 0.2451183499
trtD timeDAY 8:trtB timeDAY 8:trtC timeDAY 8:trtD
0.0707007033 -0.1102181249 -0.3006730214 -0.1259701144
Greatly appreciated!
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: Tuesday, August 9, 2011 3:20 AM
Subject: RE: [R] mixed model fitting between R and SAS
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
More information about the R-sig-mixed-models
mailing list