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.hardenb...@pngp.it
fa...@pngp.it
Skype: achazhardenberg
Mobile: +44.(0)783.266.5995
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.hardenb...@pngp.it
fa...@pngp.it
Skype: achazhardenberg
Mobile: +44.(0)783.266.5995
--------------------------------------------------------------------------------------------------------
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
