[R] p values for a GEE model

[Renaud Lancelot]

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