seizure.RdThe seizure data frame has 59 rows and 7 columns. The dataset has the
number of epiliptic seizures in each of four two-week intervals, and in a
baseline eight-week inverval, for treatment and control groups with a total
of 59 individuals.
seizureThis data frame contains the following columns:
the number of epiliptic seizures in the 1st 2-week interval
the number of epiliptic seizures in the 2nd 2-week interval
the number of epiliptic seizures in the 3rd 2-week interval
the number of epiliptic seizures in the 4th 2-week interval
an indicator of treatment
the number of epilitic seizures in a baseline 8-week interval
a numeric vector of subject age
Thall, P.F. and Vail S.C. (1990) Some covariance models for longitudinal count data with overdispersion. Biometrics 46: 657–671.
Diggle, P.J., Liang, K.Y., and Zeger, S.L. (1994) Analysis of Longitudinal Data. Clarendon Press.
data(seizure)
## Diggle, Liang, and Zeger (1994) pp166-168, compare Table 8.10
seiz.l <- reshape(seizure,
varying=list(c("base","y1", "y2", "y3", "y4")),
v.names="y", times=0:4, direction="long")
seiz.l <- seiz.l[order(seiz.l$id, seiz.l$time),]
seiz.l$t <- ifelse(seiz.l$time == 0, 8, 2)
seiz.l$x <- ifelse(seiz.l$time == 0, 0, 1)
m1 <- geese(y ~ offset(log(t)) + x + trt + x:trt, id = id,
data=seiz.l, corstr="exch", family=poisson)
summary(m1)
#>
#> Call:
#> geese(formula = y ~ offset(log(t)) + x + trt + x:trt, id = id,
#> data = seiz.l, family = poisson, corstr = "exch")
#>
#> Mean Model:
#> Mean Link: log
#> Variance to Mean Relation: poisson
#>
#> Coefficients:
#> estimate san.se wald p
#> (Intercept) 1.34761 0.1574 73.34238 0.0000
#> x 0.11184 0.1159 0.93061 0.3347
#> trt 0.02753 0.2218 0.01541 0.9012
#> x:trt -0.10473 0.2134 0.24073 0.6237
#>
#> Scale Model:
#> Scale Link: identity
#>
#> Estimated Scale Parameters:
#> estimate san.se wald p
#> (Intercept) 19.41 8.689 4.992 0.02546
#>
#> Correlation Model:
#> Correlation Structure: exch
#> Correlation Link: identity
#>
#> Estimated Correlation Parameters:
#> estimate san.se wald p
#> alpha 0.7767 0.07522 106.6 0
#>
#> Returned Error Value: 0
#> Number of clusters: 59 Maximum cluster size: 5
#>
m2 <- geese(y ~ offset(log(t)) + x + trt + x:trt, id = id,
data = seiz.l, subset = id!=49,
corstr = "exch", family=poisson)
summary(m2)
#>
#> Call:
#> geese(formula = y ~ offset(log(t)) + x + trt + x:trt, id = id,
#> data = seiz.l, subset = id != 49, family = poisson, corstr = "exch")
#>
#> Mean Model:
#> Mean Link: log
#> Variance to Mean Relation: poisson
#>
#> Coefficients:
#> estimate san.se wald p
#> (Intercept) 1.3476 0.1574 73.3424 0.00000
#> x 0.1118 0.1159 0.9306 0.33470
#> trt -0.1068 0.1937 0.3041 0.58130
#> x:trt -0.3024 0.1711 3.1248 0.07711
#>
#> Scale Model:
#> Scale Link: identity
#>
#> Estimated Scale Parameters:
#> estimate san.se wald p
#> (Intercept) 10.38 2.282 20.71 5.347e-06
#>
#> Correlation Model:
#> Correlation Structure: exch
#> Correlation Link: identity
#>
#> Estimated Correlation Parameters:
#> estimate san.se wald p
#> alpha 0.5978 0.08107 54.38 1.653e-13
#>
#> Returned Error Value: 0
#> Number of clusters: 58 Maximum cluster size: 5
#>
## Thall and Vail (1990)
seiz.l <- reshape(seizure, varying=list(c("y1","y2","y3","y4")),
v.names="y", direction="long")
seiz.l <- seiz.l[order(seiz.l$id, seiz.l$time),]
seiz.l$lbase <- log(seiz.l$base / 4)
seiz.l$lage <- log(seiz.l$age)
seiz.l$v4 <- ifelse(seiz.l$time == 4, 1, 0)
m3 <- geese(y ~ lbase + trt + lbase:trt + lage + v4,
sformula = ~ as.factor(time) - 1, id = id,
data = seiz.l, corstr = "exchangeable", family=poisson)
## compare to Model 13 in Table 4, noticeable difference
summary(m3)
#>
#> Call:
#> geese(formula = y ~ lbase + trt + lbase:trt + lage + v4, sformula = ~as.factor(time) -
#> 1, id = id, data = seiz.l, family = poisson, corstr = "exchangeable")
#>
#> Mean Model:
#> Mean Link: log
#> Variance to Mean Relation: poisson
#>
#> Coefficients:
#> estimate san.se wald p
#> (Intercept) -2.3979 0.81644 8.626 0.0033144
#> lbase 0.9222 0.08367 121.463 0.0000000
#> trt -1.5271 0.43832 12.138 0.0004940
#> lage 0.8181 0.23979 11.640 0.0006454
#> v4 -0.1647 0.06547 6.331 0.0118642
#> lbase:trt 0.6123 0.18043 11.517 0.0006894
#>
#> Scale Model:
#> Scale Link: identity
#>
#> Estimated Scale Parameters:
#> estimate san.se wald p
#> as.factor(time)1 3.378 0.8214 16.913 3.914e-05
#> as.factor(time)2 4.465 1.5305 8.510 3.531e-03
#> as.factor(time)3 7.420 3.6417 4.152 4.160e-02
#> as.factor(time)4 2.318 0.4554 25.923 3.552e-07
#>
#> Correlation Model:
#> Correlation Structure: exchangeable
#> Correlation Link: identity
#>
#> Estimated Correlation Parameters:
#> estimate san.se wald p
#> alpha 0.3866 0.06383 36.68 1.391e-09
#>
#> Returned Error Value: 0
#> Number of clusters: 59 Maximum cluster size: 4
#>
## set up a design matrix for the correlation
z <- model.matrix(~ age, data = seizure) # data is not seiz.l
## just to illustrate the scale link and correlation link
m4 <- geese(y ~ lbase + trt + lbase:trt + lage + v4,
sformula = ~ as.factor(time)-1, id = id,
data = seiz.l, corstr = "ar1", family = poisson,
zcor = z, cor.link = "fisherz", sca.link = "log")
summary(m4)
#>
#> Call:
#> geese(formula = y ~ lbase + trt + lbase:trt + lage + v4, sformula = ~as.factor(time) -
#> 1, id = id, data = seiz.l, zcor = z, family = poisson, cor.link = "fisherz",
#> sca.link = "log", corstr = "ar1")
#>
#> Mean Model:
#> Mean Link: log
#> Variance to Mean Relation: poisson
#>
#> Coefficients:
#> estimate san.se wald p
#> (Intercept) -2.4473 0.81219 9.080 0.0025847
#> lbase 0.9256 0.08600 115.852 0.0000000
#> trt -1.6537 0.45010 13.499 0.0002387
#> lage 0.8334 0.23607 12.463 0.0004151
#> v4 -0.1580 0.07104 4.945 0.0261617
#> lbase:trt 0.6545 0.18544 12.458 0.0004161
#>
#> Scale Model:
#> Scale Link: log
#>
#> Estimated Scale Parameters:
#> estimate san.se wald p
#> as.factor(time)1 1.2126 0.2423 25.04 5.617e-07
#> as.factor(time)2 1.5462 0.3656 17.88 2.351e-05
#> as.factor(time)3 2.0142 0.4795 17.64 2.664e-05
#> as.factor(time)4 0.8731 0.2074 17.72 2.553e-05
#>
#> Correlation Model:
#> Correlation Structure: ar1
#> Correlation Link: fisherz
#>
#> Estimated Correlation Parameters:
#> estimate san.se wald p
#> (Intercept) 2.54347 0.70686 12.947 0.0003204
#> age -0.04604 0.02285 4.061 0.0438922
#>
#> Returned Error Value: 0
#> Number of clusters: 59 Maximum cluster size: 4
#>