[R] gee and geepack: different results?
mendola@dssm.unipa.it
mendola at dssm.unipa.it
Fri Oct 24 18:22:05 CEST 2003
Hi, I downloaded both gee and geepack, and I am trying to understand the
differences between the two libraries.
I used the same data and estimated the same model, with a correlation
structure autoregressive of order 1. Surprisingly for me, I found very
different results. Coefficients are slightly different in value but
sometimes opposite in sign.
Moreover, the estimate of rho (correlation coefficient) given by gee is
0.5759 (see element 1.2 in the working correlation matrix) while the
estimate given by geese is 0.4519.
Could somebody explain me what happened?
Data are in dati22, which is only a part of the data for the study I
intend to perform. Here what I did ,using first gee and then geese:
> str(dati22)
`data.frame': 47 obs. of 7 variables:
$ TR : Factor w/ 4 levels "7","8","9","11": 1 1 1 1 1 1 1 1 1 1 ...
$ STAT : Factor w/ 3 levels "1","2","3": 1 1 1 1 2 2 2 2 3 3 ...
$ PIANTA: int 2 2 2 2 2 2 2 2 41 41 ...
$ ANLEP : int 1999 1998 1999 1998 1999 1998 1999 1998 1999 1998 ...
$ eta : int 12 11 10 9 11 10 11 10 14 13 ...
$ VCRE : num 5.3 6.9 11 9.9 7.9 9.2 14.2 11.9 10.5 10 ...
$ temp : num 19.7 20.0 19.7 20.0 19.7 ...
> mio1.2<- gee(VCRE ~ TR + STAT + temp + eta, id = PIANTA, +
+ data = dati22, family = Gamma, corstr = "AR-M", Mv = 1)
> summary(mio1.2)
GEE: GENERALIZED LINEAR MODELS FOR DEPENDENT DATA
gee S-function, version 4.13 modified 98/01/27 (1998)
Model:
Link: Reciprocal
Variance to Mean Relation: Gamma
Correlation Structure: AR-M , M = 1
Call:
gee(formula = VCRE ~ TR + STAT + temp + eta, id = PIANTA, data = dati22,
family = Gamma, corstr = "AR-M", Mv = 1)
Summary of Residuals:
Min 1Q Median 3Q Max
-4.5154238 -2.0766622 -0.2153521 0.9418182 6.3871421
Coefficients:
Estimate Naive S.E. Naive z Robust S.E. Robust z
(Intercept) 0.411995251 0.174052364 2.36707644 0.131632816 3.12988253
TR8 -0.001154422 0.021191593 -0.05447549 0.011807163 -0.09777305
TR9 0.019559907 0.024379471 0.80231056 0.008993803 2.17482050
TR11 -0.041092894 0.021609580 -1.90160537 0.015384050 -2.67113620
STAT2 0.023886745 0.014219390 1.67987130 0.013543550 1.76369899
STAT3 0.045749728 0.016844262 2.71604237 0.012862504 3.55682917
temp -0.020141682 0.008819798 -2.28368975 0.006851784 -2.93962600
eta 0.008021081 0.001465212 5.47434944 0.001129986 7.09838725
Estimated Scale Parameter: 0.1049601
Number of Iterations: 13
Working Correlation
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 1.0000 0.5758 0.3315 0.1909 0.1099 0.0633 0.0364 0.0210
[2,] 0.5758 1.0000 0.5758 0.3315 0.1909 0.1099 0.0633 0.0364
[3,] 0.3315 0.5758 1.0000 0.5758 0.3315 0.1909 0.1099 0.0633
[4,] 0.1909 0.3315 0.5758 1.0000 0.5758 0.3315 0.1909 0.1099
[5,] 0.1099 0.1909 0.3315 0.5758 1.0000 0.5758 0.3315 0.1909
[6,] 0.0633 0.1099 0.1909 0.3315 0.5758 1.0000 0.5758 0.3315
[7,] 0.0364 0.0633 0.1099 0.1909 0.3315 0.5758 1.0000 0.5758
[8,] 0.0210 0.0364 0.0633 0.1099 0.1909 0.3315 0.5758 1.0000
USING geese:
> mio1.2pack<-geese(VCRE ~ TR + STAT + temp + eta, id = PIANTA, data =
dati22, +
+ family = Gamma, corstr = "ar1")
Coefficients:
estimate san.se wald p
(Intercept) 0.44799 0.18279 6.007 1.43e-02
TR8 0.00444 0.00985 0.203 6.52e-01
TR9 0.02182 0.00632 11.923 5.54e-04
TR11 -0.03728 0.01407 7.024 8.04e-03
STAT2 0.02116 0.01476 2.056 1.52e-01
STAT3 0.04270 0.01245 11.758 6.06e-04
temp -0.02233 0.00973 5.273 2.17e-02
eta 0.00865 0.00136 40.683 1.79e-10
Scale Model:
Scale Link: identity
Estimated Scale Parameters:
estimate san.se wald p
(Intercept) 0.081 0.0287 7.98 0.00474
Correlation Model:
Correlation Structure: ar1
Correlation Link: identity
Estimated Correlation Parameters:
estimate san.se wald p
alpha 0.452 0.272 2.77 0.0963
Returned Error Value: 0
Number of clusters: 6 Maximum cluster size: 8
Thanks, Daria
****************************
Daria Mendola
Department of Statistics and Matematics
"Silvio Vianelli"
University di Palermo
Viale delle Scienze - Edificio 13
90128 Palermo, Italy
phone +39 091 6626210
fax +39 091 485726
email: mendola at dssm.unipa.it
More information about the R-help
mailing list