Hi all,
I have a dataset in which the output Y is observed on two groups of
patients (treatment factor T with 2 levels).
Every subject in each group is observed three times (not time points but
just technical replication).
I am interested in estimating the treatment effect and take into account
the fact that I have repeated measurements for every subject.
If I do this with repeated measures ANOVA (in which the patient is
considered a random effect) I got the following results:
library(nlme)
data<-read.table("http://146.9.88.18/uploads/dataGEE.txt",header=TRUE)
res<-lme(Y~T,random=~1|P,data=data)
summary(res)
So the p-value for significance of the treatment effect is 0.069.
I would like to use also as a variant analysis a Generalized Estimation
Equation Model, like
library(gee)
summary(gee(Y~T,id=P,data=data))
Questions:
A) Is the gee approach suitable in this case with the model formulae I
use?
B) Can I obtain a p-value for the fixed effect T ?
Thanks,
Laurentiu Tarca
[[alternative HTML version deleted]]
On Mon, 10 Apr 2006, Tarca, Adi wrote:> Hi all, > > I have a dataset in which the output Y is observed on two groups of > patients (treatment factor T with 2 levels). > > Every subject in each group is observed three times (not time points but > just technical replication). > > I am interested in estimating the treatment effect and take into account > the fact that I have repeated measurements for every subject. > > If I do this with repeated measures ANOVA (in which the patient is > considered a random effect) I got the following results: > > library(nlme) > > data<-read.table("http://146.9.88.18/uploads/dataGEE.txt",header=TRUE) > > res<-lme(Y~T,random=~1|P,data=data) > summary(res) > > So the p-value for significance of the treatment effect is 0.069. > > I would like to use also as a variant analysis a Generalized Estimation > Equation Model, like > > library(gee) > summary(gee(Y~T,id=P,data=data)) > > Questions: > > A) Is the gee approach suitable in this case with the model formulae I > use?Yes, but your choice of working correlation probably is not. Your results are basically least-squares here. I think you want corstr="exchangeable" or perhaps some form of AR if the repeated observations have a clear time base.> B) Can I obtain a p-value for the fixed effect T ?Well, the output is Coefficients: Estimate Naive S.E. Naive z Robust S.E. Robust z (Intercept) 14.781968 0.3686055 40.102410 0.3608597 40.963207 TPTLD 1.857166 0.5259206 3.531266 0.8570370 2.166961 and if you read up on the background you will know that the 'z' scores are approximately normally distributed. There are detailed comparisons of examples of the lme and gee approaches in MASS chapter 10. -- Brian D. Ripley, ripley at stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UK Fax: +44 1865 272595
2006/4/10, Tarca, Adi <atarca at med.wayne.edu>:> Hi all, > > I have a dataset in which the output Y is observed on two groups of > patients (treatment factor T with 2 levels). > > Every subject in each group is observed three times (not time points but > just technical replication). > > I am interested in estimating the treatment effect and take into account > the fact that I have repeated measurements for every subject. > > If I do this with repeated measures ANOVA (in which the patient is > considered a random effect) I got the following results: > > > > library(nlme) > > data<-read.table("http://146.9.88.18/uploads/dataGEE.txt",header=TRUE) > > res<-lme(Y~T,random=~1|P,data=data) > > summary(res) > > > > So the p-value for significance of the treatment effect is 0.069. > > I would like to use also as a variant analysis a Generalized Estimation > Equation Model, like > > library(gee) > > summary(gee(Y~T,id=P,data=data))Beware: the default within-group correlation structure is independence in the gee function (see argument corstr). I think you want an exchangeable correlation structure, i.e. the same within-group correlation for all the measurements:> Data <- read.table("http://146.9.88.18/uploads/dataGEE.txt", header = TRUE) > > library(gee) > fm1 <- gee(Y ~ T, id = P, data = Data, corstr = "exchangeable")[1] "Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27" [1] "running glm to get initial regression estimate" [1] 14.781968 1.857166> summary(fm1)GEE: GENERALIZED LINEAR MODELS FOR DEPENDENT DATA gee S-function, version 4.13 modified 98/01/27 (1998) Model: Link: Identity Variance to Mean Relation: Gaussian Correlation Structure: Exchangeable Call: gee(formula = Y ~ T, id = P, data = Data, corstr = "exchangeable") Summary of Residuals: Min 1Q Median 3Q Max -3.42512726 -0.99762726 -0.04887726 0.48282733 7.73737274 Coefficients: Estimate Naive S.E. Naive z Robust S.E. Robust z (Intercept) 14.853839 0.6069154 24.474317 0.3855683 38.524536 TPTLD 1.713788 0.8586943 1.995807 0.8423924 2.034429 Estimated Scale Parameter: 3.945557 Number of Iterations: 2 Working Correlation [,1] [,2] [,3] [1,] 1.000000 0.897846 0.897846 [2,] 0.897846 1.000000 0.897846 [3,] 0.897846 0.897846 1.000000> Questions: > > A) Is the gee approach suitable in this case with the model formulae I > use? > > B) Can I obtain a p-value for the fixed effect T ?You can extract the fixed-effect coefficients with coef:> coef(summary(fm1))Estimate Naive S.E. Naive z Robust S.E. Robust z (Intercept) 14.853839 0.6069154 24.474317 0.3855683 38.524536 TPTLD 1.713788 0.8586943 1.995807 0.8423924 2.034429 Then, get the P values using a normal approximation for the distribution of z:> 2 * pnorm(abs(coef(summary(fm1))[,5]), lower.tail = FALSE)(Intercept) TPTLD 0.00000000 0.04190831 Best, Renaud> > > > Thanks, > > > > Laurentiu Tarca > > > > > > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html >-- Renaud LANCELOT D?partement Elevage et M?decine V?t?rinaire (EMVT) du CIRAD Directeur adjoint charg? des affaires scientifiques CIRAD, Animal Production and Veterinary Medicine Department Deputy director for scientific affairs Campus international de Baillarguet TA 30 / B (B?t. B, Bur. 214) 34398 Montpellier Cedex 5 - France T?l +33 (0)4 67 59 37 17 Secr. +33 (0)4 67 59 39 04 Fax +33 (0)4 67 59 37 95
On Tue, 11 Apr 2006, Renaud Lancelot wrote:> 2006/4/10, Tarca, Adi <atarca at med.wayne.edu>: >> Hi all, >> >> I have a dataset in which the output Y is observed on two groups of >> patients (treatment factor T with 2 levels). >> >> Every subject in each group is observed three times (not time points but >> just technical replication). >> >> I am interested in estimating the treatment effect and take into account >> the fact that I have repeated measurements for every subject. >><<Snip>>>> I would like to use also as a variant analysis a Generalized Estimation >> Equation Model, like >> >> library(gee) >> >> summary(gee(Y~T,id=P,data=data)) > > Beware: the default within-group correlation structure is independence > in the gee function (see argument corstr). I think you want an > exchangeable correlation structure, i.e. the same within-group > correlation for all the measurements:He has a linear model with the same number of observations for each person and no covariates that vary within a person. The independence and exchangeable working correlations will give identical answers. -thomas Thomas Lumley Assoc. Professor, Biostatistics tlumley at u.washington.edu University of Washington, Seattle
On Tue, 11 Apr 2006, Thomas Lumley wrote:> On Tue, 11 Apr 2006, Renaud Lancelot wrote: > >> 2006/4/10, Tarca, Adi <atarca at med.wayne.edu>: >>> Hi all, >>> >>> I have a dataset in which the output Y is observed on two groups of >>> patients (treatment factor T with 2 levels). >>> >>> Every subject in each group is observed three times (not time points but >>> just technical replication). >>> >>> I am interested in estimating the treatment effect and take into account >>> the fact that I have repeated measurements for every subject. >>> > <<Snip>> >>> I would like to use also as a variant analysis a Generalized Estimation >>> Equation Model, like >>> >>> library(gee) >>> >>> summary(gee(Y~T,id=P,data=data)) >> >> Beware: the default within-group correlation structure is independence >> in the gee function (see argument corstr). I think you want an >> exchangeable correlation structure, i.e. the same within-group >> correlation for all the measurements: > > > He has a linear model with the same number of observations for each personNot so: some have 3 and some have 2, and the two levels of T are not quite balanced (29/28).> and no covariates that vary within a person. The independence and > exchangeable working correlations will give identical answers.They do not in this case (nor should they, I believe). -- Brian D. Ripley, ripley at stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UK Fax: +44 1865 272595
On Tue, 11 Apr 2006, Prof Brian Ripley wrote:> On Tue, 11 Apr 2006, Thomas Lumley wrote: >> >> He has a linear model with the same number of observations for each person > > Not so: some have 3 and some have 2, and the two levels of T are not quite > balanced (29/28). > >> and no covariates that vary within a person. The independence and >> exchangeable working correlations will give identical answers. > > They do not in this case (nor should they, I believe). >You are right, of course. I believed the message rather than looking at the data. -thomas Thomas Lumley Assoc. Professor, Biostatistics tlumley at u.washington.edu University of Washington, Seattle