Chapter 13, Section 13.4

 

Muscatine Coronary Risk Factor Study

Marginal Logistic Regression Model

 

 

 

library(foreign)

ds <- read.dta("c:/muscatine.dta")

musc <- na.omit(ds)

attach(musc)

 

cage <- age - 12

> cage2 <- cage*cage

 

library(geepack)

> model1 <- ordgee(ordered(y) ~ factor(gender) + cage + cage2 + 

+     factor(gender):cage + factor(gender):cage2, id = id, 

+     waves=occasion, mean.link="logit", corstr=("unstructured"))

 

Warning message:

In binomial(link) : use of binomial(link=link) is deprecated

 

summary(model1)

 

Call:

ordgee(formula = ordered(y) ~ factor(gender) + cage + cage2 + 

    factor(gender):cage + factor(gender):cage2, id = id, 

    waves = occasion, mean.link = "logit", corstr = ("unstructured"))

 

Mean Model:

 Mean Link:                 logit 

 Variance to Mean Relation: binomial 

 

 Coefficients:

                          estimate      san.se        wald            p

Inter:0               -1.214613103 0.050571150 576.8597850 0.000000e+00

factor(gender)1        0.115330450 0.071158497   2.6268450 1.050703e-01

cage                   0.037419375 0.013263832   7.9589357 4.785054e-03

cage2                 -0.017437692 0.003378786  26.6352422 2.457205e-07

factor(gender)1:cage   0.007510802 0.018268075   0.1690390 6.809673e-01

factor(gender)1:cage2  0.003860069 0.004632095   0.6944407 4.046580e-01

 

Scale is fixed.

 

Correlation Model:

 Correlation Structure:     unstructured 

 Correlation Link:          log 

 

 Estimated Correlation Parameters:

        estimate    san.se     wald p

alpha.1 3.130702 0.1535950 415.4599 0

alpha.2 2.408103 0.1455606 273.6921 0

alpha.3 2.793549 0.1351264 427.3978 0

 

Returned Error Value:    0 

Number of clusters:   4856   Maximum cluster size: 3

 

 

 

> model2 <- ordgee(ordered(y) ~ factor(gender) + cage + cage2, id = id, 

+    waves=occasion, mean.link="logit", corstr=("unstructured"))

 

Warning message:

In binomial(link) : use of binomial(link=link) is deprecated

 

summary(model2)

 

Call:

ordgee(formula = ordered(y) ~ factor(gender) + cage + cage2, 

    id = id, waves = occasion, mean.link = "logit", corstr = ("unstructured"))

 

Mean Model:

 Mean Link:                 logit 

 Variance to Mean Relation: binomial 

 

 Coefficients:

                   estimate      san.se      wald            p

Inter:0         -1.22930460 0.047708183 663.94676 0.000000e+00

factor(gender)1  0.14416724 0.062678388   5.29051 2.144194e-02

cage             0.04136873 0.009107800  20.63086 5.569105e-06

cage2           -0.01544435 0.002307421  44.80081 2.181333e-11

 

Scale is fixed.

 

Correlation Model:

 Correlation Structure:     unstructured 

 Correlation Link:          log 

 

 Estimated Correlation Parameters:

        estimate    san.se     wald p

alpha.1 3.129770 0.1535978 415.1976 0

alpha.2 2.406007 0.1455562 273.2326 0

alpha.3 2.793338 0.1350829 427.6084 0

 

Returned Error Value:    0 

Number of clusters:   4856   Maximum cluster size: 3

 

 

 

 

 

Clinical Trial of Antibiotics for Leprosy

Marginal Log-linear Regression Model

 

library(foreign)

ds <- read.dta("c:/leprosy.dta") 

dslong <- reshape(ds, idvar="id", varying=c("y1","y2"), 

+    v.names="y", timevar="time", time=0:1, direction="long")

dslong <- dslong[order(dslong$iddslong$time),]

attach(dslong

 

drugn <- as.numeric(drug)

timeA <- time*(drugn == 2)

timeB <- time*(drugn == 3)

timeAB <- time*I(drugn != 1)

  

library(geepack)

summary(geeglm(y ~ time + timeA + timeB, data=dslong, id = id, 

+    waves=time, family=poisson("log"), corstr="exch", std.err="san.se"))

 

Call:

geeglm(formula = y ~ time + timeA + timeB, family = poisson("log"), 

    data = dslong, id = id, waves = time, corstr = "exch", std.err = "san.se")

 

 Coefficients:

             Estimate   Std.err    Wald Pr(>|W|)    

(Intercept)  2.373354  0.082801 821.582   <2e-16 ***

time        -0.002877  0.151947   0.000   0.9849    

timeA       -0.562571  0.218491   6.630   0.0100 *  

timeB       -0.495284  0.230981   4.598   0.0320 *  

---

Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

 

Estimated Scale Parameters:

            Estimate Std.err

(Intercept)    3.214  0.5199

 

Correlation: Structure = exchangeable  Link = identity 

 

Estimated Correlation Parameters:

      Estimate Std.err

alpha   0.7384 0.07194

Number of clusters:   30   Maximum cluster size: 2 

 

 

summary(geeglm(y ~ time + timeAB, data=dslong, id = id, 

+    waves=time, family=poisson("log"), corstr="exch", std.err="san.se"))

 

Call:

geeglm(formula = y ~ time + timeAB, family = poisson("log"), 

    data = dslong, id = id, waves = time, corstr = "exch", std.err = "san.se")

 

 Coefficients:

            Estimate  Std.err   Wald Pr(>|W|)    

(Intercept)  2.37335  0.08280 821.58   <2e-16 ***

time        -0.00286  0.15195   0.00   0.9850    

timeAB      -0.52783  0.19498   7.33   0.0068 ** 

---

Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

 

Estimated Scale Parameters:

            Estimate Std.err

(Intercept)     3.23   0.541

 

Correlation: Structure = exchangeable  Link = identity 

 

Estimated Correlation Parameters:

      Estimate Std.err

alpha    0.738  0.0715

Number of clusters:   30   Maximum cluster size: 2 

 

 

 

 

 

Arthritis Clinical Trial

Marginal Proportional Odds (Ordinal) Regression Model

 

 

library(foreign)

ds <- read.dta("c:/arthritis.dta")

dslong <- reshape(ds, idvar="id", varying=c("y1","y2","y3","y4"), 

+    v.names="y", timevar="time", time=1:4, direction="long")

dslong <- dslong[order(dslong$iddslong$time),]

dslong <- dslong[complete.cases(dslong$y),] 

attach(dslong

 

month <- 2*(time-1)

sqrtmonth <- month^0.5

 

library(geepack)

 

## Note: Results from ordgee below are discernibly different 

## from results obtained using alternative statistical 

## software (e.g., SAS, Stata). The estimated coefficients do not 

## correspond to those under a "working independence" assumption.

## Commands and output presented for illustrative purposes only.

 

summary(ordgee(ordered(y) ~ trt + sqrtmonth + trt:sqrtmonth

+    id = id, waves=time, corstr=("independence")))

 

Call:

ordgee(formula = ordered(y) ~ trt + sqrtmonth + trt:sqrtmonth

    id = id, waves = time, corstr = ("independence"))

 

Mean Model:

 Mean Link:                 logit 

 Variance to Mean Relation: binomial 

 

 Coefficients:

                 estimate     san.se        wald            p

Inter:1        3.18050940 0.22094412 207.2184734 0.000000e+00

Inter:2        1.14013965 0.17883503  40.6453531 1.825208e-10

Inter:3       -0.63240562 0.16823727  14.1301589 1.705868e-04

Inter:4       -2.56786564 0.22013455 136.0718133 0.000000e+00

trt           -0.14795869 0.24712397   0.3584687 5.493579e-01

sqrtmonth     -0.24775817 0.08521025   8.4541958 3.642023e-03

trt:sqrtmonth -0.08530304 0.12045483   0.5015111 4.788369e-01

 

Scale is fixed.

 

Correlation Model:

 Correlation Structure:     independence 

 

Returned Error Value:    0 

Number of clusters:   303   Maximum cluster size: 4 

 

Warning message:

In binomial(link) : use of binomial(link=link) is deprecated