Hi,

I am trying to fit a ordinal GEE model using ordgee {geepack}. In order to check the validity of the function, I specified the correlation structure as independence (i.e. constr = "independence") and compared the result with that using polr {MASS}. 

Because a GEE model with an independent working correlation structure is equivalent to an ordinary GLM model, we would expect the same results from both functions (there might be small differences due to dissimilarity in approximation methods). However, I got a very different results from the two approaches.

Here is my R codes and ouputs for the two approaches:
1.) 
library(geepack)summary(ordgee(UACR~HR+SP+bmi+diabetes+age+HDL+ht.med+male+SCHV+SCRV+smoke+statin+AP,id=x18,mean.link="logit",corstr="independence"))

Call:
ordgee(formula = UACR ~ HR + SP + bmi + diabetes + age + HDL + ht.med + male + SCHV + SCRV + smoke + statin + AP, id = x18, mean.link = "logit", corstr = "independence")

Mean Model:
 Mean Link:                 logit 
 Variance to Mean Relation: binomial 

 Coefficients:
              estimate      san.se         wald           p
Inter:0  -1.5843906608 2.408716706 4.326663e-01 0.510683514
Inter:1  -2.6094525748 2.401539294 1.180645e+00 0.277224675
Inter:2  -3.6006918946 2.410172287 2.231905e+00 0.135187061
HR        0.0491831800 0.018913783 6.762020e+00 0.009311831
SP       -0.0287228998 0.011465183 6.276168e+00 0.012237244
bmi       0.0363585224 0.029345236 1.535102e+00 0.215348334
diabetes  0.2557417872 0.464391207 3.032742e-01 0.581837057
age       0.0188629344 0.022226239 7.202555e-01 0.396060127
HDL      -0.0210835675 0.014822770 2.023156e+00 0.154916842
ht.med    0.0089636059 0.383824776 5.453804e-04 0.981368392
male      1.2806017705 0.473921757 7.301547e+00 0.006889526
SCHV     -0.0008670472 0.005862042 2.187699e-02 0.882414803
SCRV      0.5667065765 0.879731084 4.149699e-01 0.519457682
smoke    -0.6064323366 0.352556881 2.958737e+00 0.085414655
statin   -0.0008867941 0.408462933 4.713465e-06 0.998267753
AP        0.1159595043 0.036258436 1.022809e+01 0.001383178

Scale is fixed.

Correlation Model:
 Correlation Structure:     independence 

Returned Error Value:    0 
Number of clusters:   285   Maximum cluster size: 6 

2.)
library(MASS)
summary(polr(UACR~HR+SP+bmi+diabetes+age+HDL+ht.med+male+SCHV+SCRV+smoke+statin+AP,method="logistic"))

Re-fitting to get Hessian

Call:
polr(formula = UACR ~ HR + SP + bmi + diabetes + age + HDL + ht.med + male + SCHV + SCRV + smoke + statin + AP, method = "logistic")

Coefficients:
                 Value  Std. Error      t value
HR        0.0004674288 0.009690505  0.048235757
SP       -0.0106449514 0.005305501 -2.006398816
bmi       0.0298158386 0.015501658  1.923396794
diabetes  0.6563205176 0.251591908  2.608671009
age       0.0034051085 0.010467630  0.325298895
HDL      -0.0153783914 0.006227563 -2.469407699
ht.med   -0.1146770998 0.194130303 -0.590722303
male      0.1798842955 0.227476463  0.790782014
SCHV      0.0022770371 0.002675591  0.851040793
SCRV      0.0010954808 0.440261143  0.002488252
smoke     0.0207174113 0.165287870  0.125341390
statin    0.4320833418 0.206676105  2.090630370
AP        0.0419551861 0.015799215  2.655523401

Intercepts:
    Value   Std. Error t value
0|1 -0.9883  1.1258    -0.8778
1|2  0.1712  1.1246     0.1522
2|3  1.3306  1.1258     1.1819

Residual Deviance: 1376.931 
AIC: 1408.931 

Is there anything wrong with my code? I guess my real question should be, could I safely use the ordgee program to fit a ordinal GEE model with exchangeable correlation structure?

Thanks very much for your time and professional opinion,
High Seng


	[[alternative HTML version deleted]]


