library("gcmr") ################### Section 4.1 data("epilepsy", package = "gcmr") ## independence model mod.ind <- gcmr(counts ~ offset(log(time)) + visit + trt + visit:trt, data = epilepsy, subset = (id != 49), marginal = negbin.marg, cormat = cluster.cormat(id, type = "ind")) summary(mod.ind) ## sandwich standard errors library("sandwich") library("lmtest") coeftest(mod.ind, vcov. = sandwich(mod.ind)) ## AR(1) correlation mod.ar1 <- update(mod.ind, cormat = cluster.cormat(id, "ar1"), seed = 12345, nrep = 100) coeftest(mod.ar1, vcov. = sandwich(mod.ar1)) ## alternative analysis with GEE library("geepack") gee.ar1 <- geeglm(counts ~ offset(log(time)) + visit + trt + visit:trt, data = epilepsy, id = id, subset = (id != 49), family = poisson, corstr = "ar1") summary(gee.ar1) ## profile likelihood (time consuming) profile(mod.ar1, which = 4) ################### Section 4.2 data("HUR", package = "gcmr") plot(HUR, ylab = "rate", xlab = "time") ## beta regression with trend and ARMA(1,3) errors trend <- scale(time(HUR)) mod <- gcmr(HUR ~ trend | trend, marginal = beta.marg, cormat = arma.cormat(1, 3)) summary(mod) ## diagnostic plots par(mfrow = c(2, 2)) plot(mod) ################### Section 4.3 data("malaria", package = "gcmr") ## spatial distances library("sp") D <- spDists(cbind(malaria$x, malaria$y))/1000 ## model without area mod <- gcmr(cbind(cases, size - cases) ~ netuse + I(green/100) + phc, data = malaria, marginal = binomial.marg, cormat = matern.cormat(D), seed = 12345) summary(mod) ## model with area mod.area <- update(mod, . ~ . + area) summary(mod.area) AIC(mod, mod.area) ## diagnostic plots par(mfrow = c(2, 2)) plot(mod.area)