## generate clustered data
gendat <- function(ncl, clsz) {
## ncl: number of clusters
## clsz: cluster size (all equal)
id <- rep(1:ncl, each = clsz)
visit <- rep(1:clsz, ncl)
n <- ncl * clsz
x1 <- rbinom(n, 1, 0.5) ## within cluster varying binary covariate
x2 <- runif(n, 0, 1) ## within cluster varying continuous covariate
## the true correlation coefficient rho for an ar(1)
## correlation structure is 2/3
rho <- 2/3
rhomat <- rho ^ outer(1:4, 1:4, function(x, y) abs(x - y))
chol.u <- chol(rhomat)
noise <- as.vector(sapply(1:ncl, function(x) chol.u %*% rnorm(clsz)))
y <- 1 + 3 * x1 - 2 * x2 + noise
dat <- data.frame(y, id, visit, x1, x2)
dat
}
simdat <- gendat(100, 4)
fit0 <- geese(y ~ x1, id = id, data = simdat, corstr = "un")
fit1 <- geese(y ~ x1 + x2, id = id, data = simdat, corstr = "un")
compCoef(fit0, fit1)
Run the code above in your browser using DataLab