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)
m2 <- geese(y ~ offset(log(t)) + x + trt + x:trt, id = id,
data = seiz.l, subset = id!=49,
corstr = "exch", family=poisson)
summary(m2)
## Using fixed correlation matrix
cor.fixed <- matrix(c(1, 0.5, 0.25, 0.125, 0.125,
0.5, 1, 0.25, 0.125, 0.125,
0.25, 0.25, 1, 0.5, 0.125,
0.125, 0.125, 0.5, 1, 0.125,
0.125, 0.125, 0.125, 0.125, 1), 5, 5)
cor.fixed
zcor <- rep(cor.fixed[lower.tri(cor.fixed)], 59)
m3 <- geese(y ~ offset(log(t)) + x + trt + x:trt, id = id,
data = seiz.l, family = poisson,
corstr = "fixed", zcor = zcor)
summary(m3)
data(ohio)
fit <- geese(resp ~ age + smoke + age:smoke, id=id, data=ohio,
family=binomial, corstr="exch", scale.fix=TRUE)
summary(fit)
fit.ar1 <- geese(resp ~ age + smoke + age:smoke, id=id, data=ohio,
family=binomial, corstr="ar1", scale.fix=TRUE)
summary(fit.ar1)
###### simulated data
## a function to generate a dataset
gendat <- function() {
id <- gl(50, 4, 200)
visit <- rep(1:4, 50)
x1 <- rbinom(200, 1, 0.6) ## within cluster varying binary covariate
x2 <- runif(200, 0, 1) ## within cluster varying continuous covariate
phi <- 1 + 2 * x1 ## true scale model
## the true correlation coefficient rho for an ar(1)
## correlation structure is 0.667.
rhomat <- 0.667 ^ outer(1:4, 1:4, function(x, y) abs(x - y))
chol.u <- chol(rhomat)
noise <- as.vector(sapply(1:50, function(x) chol.u %*% rnorm(4)))
e <- sqrt(phi) * noise
y <- 1 + 3 * x1 - 2 * x2 + e
dat <- data.frame(y, id, visit, x1, x2)
dat
}
dat <- gendat()
fit <- geese(y ~ x1 + x2, id = id, data = dat, sformula = ~ x1,
corstr = "ar1", jack = TRUE, j1s = TRUE, fij = TRUE)
summary(fit)
#### create user-defined design matrix of unstrctured correlation.
#### in this case, zcor has 4*3/2 = 6 columns, and 50 * 6 = 300 rows
zcor <- genZcor(clusz = rep(4, 50), waves = dat$visit, "unstr")
zfit <- geese(y ~ x1 + x2, id = id, data = dat, sformula = ~ x1,
corstr = "userdefined", zcor = zcor,
jack = TRUE, j1s = TRUE, fij = TRUE)
summary(zfit)
#### Now, suppose that we want the correlation of 1-2, 2-3, and 3-4
#### to be the same. Then zcor should have 4 columns.
z2 <- matrix(NA, 300, 4)
z2[,1] <- zcor[,1] + zcor[,4] + zcor[,6]
z2[,2:4] <- zcor[, c(2, 3, 5)]
summary(geese(y ~ x1 + x2, id = id, data = dat, sformula = ~ x1,
corstr = "userdefined", zcor = z2,
jack = TRUE, j1s = TRUE, fij = TRUE))
#### Next, we introduce non-constant cluster sizes by
#### randomly selecting 60 percent of the data
good <- sort(sample(1:nrow(dat), .6 * nrow(dat)))
mdat <- dat[good,]
summary(geese(y ~ x1 + x2, id = id, data = mdat, waves = visit,
sformula = ~ x1, corstr="ar1",
jack = TRUE, j1s = TRUE, fij = TRUE))Run the code above in your browser using DataLab