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
different fits for geese and geeglm in geepack?
2 messages · Achaz von Hardenberg, Paul Johnson
2 days later
On Wed, Nov 25, 2009 at 5:35 AM, Achaz von Hardenberg
<achaz.hardenberg at gmail.com> wrote:
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.
If you gave us the exact commands that were run & the data, we could probably figure this out. I was interested to observe the problem you found, but can't reproduce it. I've installed geepack and run the examples for geese.fit, and then too exact same options and put them inside geeglm. I get the exact same results. As a result, I think you have either not run exactly the same model in the 2 cases, or geeglm is putting some settings in an unexpected way. I think you will find your answer if you study the defaults as they are set in geeglm. Type
geeglm
It will show the actual formula, you can see options that the function is setting on your behalf. It appears to me geeglm assumes your link function is identity, but the other might offer you the usual poisson with a log link. Then again, with a working example of geese.fit and geeglm that give different estimates, I could probably get further.
Paul E. Johnson Professor, Political Science 1541 Lilac Lane, Room 504 University of Kansas