Skip to content

mixed model fitting between R and SAS

2 messages · array chip, ONKELINX, Thierry

#
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
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: test.txt
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20110807/f9e1ff72/attachment.txt>
#
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