[R-sig-ME] different fits for geese and geeglm in geepack?

[Achaz von Hardenberg]

>>
Dear all,

I am trying to fit a GEE model on eagle productivity (number of  
hatched offspring per nest) using the geeglm function in the library  
geepack and I found an odd result.
My understanding is that the function geese and geeglm should give the  
same fits, as actually geeglm uses geese to fit the model, providing a  
glm style output.
However, if I fit the same model with geeglm and geese I get slightly  
different estimates of the parameters. Most striking  is the  
difference between the estimates of the correlation parameter where  
the differences are huge:
(geese: alpha= -0.0727, se= 0.0608; geeglm:  -0.219, se=  0.091).

Anybody knows why this is like that and which of the two I should  
rather trust? (I attach the outputs of the two models at the end of  
this mail)

thanks a lot for your hints!

Achaz von Hardenberg

PS:
One more question actually: anybody has got some code to calculate the  
QICu values to compare GEE models with and without specific fixed  
factors using geepack?

##################
GEESE OUTPUT:

Call:
geese(formula = PROD ~ as.factor(clas3) + Twinter + Tprecova,
     id = Nterr, waves = anno, data = aquile.dat2, family = poisson,
     corstr = "ar1")

Mean Model:
  Mean Link:                 log
  Variance to Mean Relation: poisson

  Coefficients:
                   estimate san.se  wald        p
(Intercept)        -1.7850 0.3661 23.77 1.08e-06
as.factor(clas3)2   0.7165 0.2475  8.38 3.79e-03
as.factor(clas3)3   0.5052 0.3368  2.25 1.34e-01
Twinter            -0.1066 0.0464  5.28 2.16e-02
Tprecova           -0.0549 0.0368  2.22 1.36e-01

Scale Model:
  Scale Link:                identity

  Estimated Scale Parameters:
             estimate san.se wald        p
(Intercept)    0.764  0.123 38.6 5.17e-10

Correlation Model:
  Correlation Structure:     ar1
  Correlation Link:          identity

  Estimated Correlation Parameters:
       estimate san.se wald     p
alpha  -0.0727 0.0608 1.43 0.232

Returned Error Value:    0
Number of clusters:   21   Maximum cluster size: 20

###############################################
GEEGLM OUTPUT:

Call:
geeglm(formula = PROD ~ as.factor(clas3) + Twinter + Tprecova,
     family = poisson, data = aquile.dat2, id = Nterr, waves = anno,
     corstr = "ar1")

  Coefficients:
                   Estimate Std.err  Wald Pr(>|W|)
(Intercept)        -1.7992  0.3751 23.01  1.6e-06 ***
as.factor(clas3)2   0.7412  0.2623  7.99   0.0047 **
as.factor(clas3)3   0.5253  0.3335  2.48   0.1152
Twinter            -0.1083  0.0473  5.24   0.0220 *
Tprecova           -0.0483  0.0354  1.86   0.1721
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Estimated Scale Parameters:
             Estimate Std.err
(Intercept)    0.773   0.117

Correlation: Structure = ar1  Link = identity

Estimated Correlation Parameters:
       Estimate Std.err
alpha   -0.219   0.091
Number of clusters:   21   Maximum cluster size: 20


Dr. Achaz von Hardenberg
--------------------------------------------------------------------------------------------------------
Centro Studi Fauna Alpina - Alpine Wildlife Research Centre
Servizio Sanitario e della Ricerca Scientifica
Parco Nazionale Gran Paradiso, Degioz, 11, 11010-Valsavarenche (Ao),  
Italy

Present address:
National Centre for Statistical Ecology
School of Mathematics, Statistics and Actuarial Science,
University of Kent,  Canterbury, UK

E-mail: achaz.hardenberg at pngp.it
	     fauna at pngp.it
Skype: achazhardenberg
Mobile: +44.(0)783.266.5995