mixed model fitting between R and SAS

Tue, Aug 9, 2011 3:20 AM

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

mixed model fitting between R and SAS

Thread (2 messages)