Skip to content

different fits for geese and geeglm in geepack?

2 messages · Achaz von Hardenberg, Paul Johnson

#
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
2 days later
#
On Wed, Nov 25, 2009 at 5:35 AM, Achaz von Hardenberg
<achaz.hardenberg at gmail.com> wrote:
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
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.