library(geepack)
## Load in Seizure Data
seizures <- read.csv("data/Epilepsy/seizures_full.csv")
# Should re-shape this to be long
seizures_long <- reshape(
data = seizures,
varying = c("y0","y1","y2","y3","y4","offset0","offset1","offset2","offset3","offset4"),
sep = "",
direction = "long",
idvar = "ID"
)
# Plot the data by treatment group
lattice::xyplot(y ~ time|as.factor(Trt),
groups = ID,
data = seizures_long,
panel = function(x, y){
lattice::panel.xyplot(x, y, type='p')
lattice::panel.linejoin(x, y, fun=mean, horizontal=F, lwd=2, col=1)
})
# Maybe an outlier?
plot(y ~ time, data = seizures_long, subset = (ID == 49))
# Remove ID: 49 from the data.
seizures_long <- seizures_long[which(seizures_long$ID != 49), ]
### GEEGLM is very picky
### Data *must* be ordered by ID, then by time
seizures_long <- seizures_long[order(seizures_long$ID, seizures_long$time), ]
# We are interested in the following questions:
## What is the relative rate, in the placebo group, post-vs-pre treatment.
## Is there a significant difference in pre-vs-post RRs between the two treatment arms?
## Is there a treatment effect on the rates of seizures, overall?
## Is there evidence of overdisperson?
## What is the estimated correlation?
## Does baseline age change any of these conclusions?
# First Let's try to Decide on Correlation
uns <- geeglm(
y ~ offset(log(offset)) + Age*Trt*as.factor(time),
data = seizures_long,
family = poisson,
id = ID,
corstr = "uns"
)
exch <- geeglm(
y ~ offset(log(offset)) + Age*Trt*as.factor(time),
data = seizures_long,
family = poisson,
id = ID,
corstr = "exch"
)
ar1 <- geeglm(
y ~ offset(log(offset)) + Age*Trt*as.factor(time),
data = seizures_long,
family = poisson,
id = ID,
corstr = "ar1"
)
cbind(QIC(uns),QIC(exch),QIC(ar1))
## Exchangeable Seems Preferable!
# We want pre-vs-post treatment indicator, rather than time
seizures_long$post <- as.numeric(seizures_long$time != 0)
# We want questions regarding RR in placebo and Trt group, pre and post
model1 <- geeglm(
y ~ offset(log(offset)) + Trt*post,
family = poisson,
data = seizures_long,
id = ID,
corstr = "exch"
)
summary(model1)
# What is the relative rate, in the placebo group, post-vs-pre treatment.
## placebo group, pre-treatment: b0
## placebo group, post-treatment: b0 + b2
## => log[Relative Rate]: b2
RR1 <- exp(coef(model1)['post']) # Is this significant? (Summary says no. But)
Z <- coef(model1)['post']/sqrt(diag(vcov(model1))['post'])
2*(1 - pnorm(Z))
# Is there a significant difference in pre-vs-post RRs between the two treatment arms?
## treatment group, pre-treatment: b0 + b1
## treatment group, post-treatment: b0 + b1 + b2 + b3
## => log[Relative Rate]: b2 + b3
RR2 <- exp(coef(model1)['post'] + coef(model1)['Trt:post'])
RR2
# If they are the same between Trt = 0 and Trt = 1, then b3 = 0
Z <- abs(coef(model1)['Trt:post']/sqrt(diag(vcov(model1))['Trt:post']))
2*(1 - pnorm(Z))
## Is there a treatment effect on the rates of seizures, overall?
# H0: b1 = b3 = 0
# Can compute this directly as we did last time, LB = 0
# Can also use anova
model1.reduced <- update(model1, formula = y ~ offset(log(offset)) + post)
anova(model1.reduced, model1) # p = 0.16, do not reject H0.
## Is there evidence of overdisperson?
# variance = phi*mu_{ij}
summary(model1)$dispersion # 10.4 => Std.err 2.28
## What is the estimated correlation?
summary(model1)$corr # 0.598 => Std.err 0.0811
## Does baseline age change any of these conclusions?
# IDEA: try to work through the factors again.
# However, we can see whether adding age at all helps, easily
model2 <- geeglm(
y ~ offset(log(offset)) + Trt*post + Trt*Age + Age*post,
family = poisson,
data = seizures_long,
id = ID,
corstr = "exch"
)
# This adds 3 additional factors compared to model1
# Can use a Wald based test of H0: Age = Trt:Age = post:Age = 0
# Can use ANOVA in R
anova(model2, model1) # Reject H0! Something here matters!
summary(model2) # Trt:Age does not seem to matter; drop that
model2.reduced <- geeglm(
y ~ offset(log(offset)) + Trt*post + Age*post,
family = poisson,
data = seizures_long,
id = ID,
corstr = "exch"
)
# Now we could re-answer the questions above, using this new model!
# Note, we will have to indicate that all conclusions are "at fixed age!"
# Placebo baseline vs. post: b0 + b3*Age vs. b0 + b2 + b3*Age
# still just H0: b2 = 0
RR1.new <- exp(coef(model2.reduced)['post']) # Is this significant? (Summary says no. But)
Z.new <- abs(coef(model2.reduced)['post']/sqrt(diag(vcov(model2.reduced))['post']))
2*(1 - pnorm(Z.new)) # Reject H0!
# etc.
QIC(model2.reduced)
QIC(model1.reduced)