Nothing
### R code from vignette source 'Ch_analysing_longitudinal_dataII.Rnw'
###################################################
### code chunk number 1: setup
###################################################
rm(list = ls())
s <- search()[-1]
s <- s[-match(c("package:base", "package:stats", "package:graphics", "package:grDevices",
"package:utils", "package:datasets", "package:methods", "Autoloads"), s)]
if (length(s) > 0) sapply(s, detach, character.only = TRUE)
if (!file.exists("tables")) dir.create("tables")
if (!file.exists("figures")) dir.create("figures")
set.seed(290875)
options(prompt = "R> ", continue = "+ ",
width = 63, # digits = 4,
SweaveHooks = list(leftpar = function()
par(mai = par("mai") * c(1, 1.05, 1, 1))))
HSAURpkg <- require("HSAUR")
if (!HSAURpkg) stop("cannot load package ", sQuote("HSAUR"))
rm(HSAURpkg)
a <- Sys.setlocale("LC_ALL", "C")
book <- TRUE
refs <- cbind(c("AItR", "SI", "CI", "ANOVA", "MLR", "GLM",
"DE", "RP", "SA", "ALDI", "ALDII", "MA", "PCA",
"MDS", "CA"), 1:15)
ch <- function(x, book = TRUE) {
ch <- refs[which(refs[,1] == x),]
if (book) {
return(paste("Chapter~\\\\ref{", ch[1], "}", sep = ""))
} else {
return(paste("Chapter~\\\\ref{", ch[2], "}", sep = ""))
}
}
###################################################
### code chunk number 2: ALDII-gee
###################################################
library("gee")
###################################################
### code chunk number 3: ALDII-BtheB-data
###################################################
data("BtheB", package = "HSAUR")
BtheB$subject <- factor(rownames(BtheB))
nobs <- nrow(BtheB)
BtheB_long <- reshape(BtheB, idvar = "subject",
varying = c("bdi.2m", "bdi.4m", "bdi.6m", "bdi.8m"),
direction = "long")
BtheB_long$time <- rep(c(2, 4, 6, 8), rep(nobs, 4))
###################################################
### code chunk number 4: ALDII-BtheB-geefit-indep
###################################################
osub <- order(as.integer(BtheB_long$subject))
BtheB_long <- BtheB_long[osub,]
btb_gee <- gee(bdi ~ bdi.pre + treatment + length + drug,
data = BtheB_long, id = subject, family = gaussian,
corstr = "independence")
###################################################
### code chunk number 5: ALDII-BtheB-geefit-ex
###################################################
btb_gee1 <- gee(bdi ~ bdi.pre + treatment + length + drug,
data = BtheB_long, id = subject, family = gaussian,
corstr = "exchangeable")
###################################################
### code chunk number 6: ALDII-BtheB-geesummary
###################################################
summary(btb_gee)
###################################################
### code chunk number 7: ALDII-BtheB-gee1summary
###################################################
summary(btb_gee1)
###################################################
### code chunk number 8: ALDII-respiratory-data
###################################################
data("respiratory", package = "HSAUR")
resp <- subset(respiratory, month > "0")
resp$baseline <- rep(subset(respiratory, month == "0")$status, rep(4, 111))
resp$nstat <- as.numeric(resp$status == "good")
###################################################
### code chunk number 9: ALDII-respiratory-fit
###################################################
resp_glm <- glm(status ~ centre + treatment + sex + baseline +
age, data = resp, family = "binomial")
resp_gee1 <- gee(nstat ~ centre + treatment + sex + baseline +
age, data = resp, family = "binomial", id = subject,
corstr = "independence", scale.fix = TRUE, scale.value = 1)
resp_gee2 <- gee(nstat ~ centre + treatment + sex + baseline +
age, data = resp, family = "binomial", id = subject,
corstr = "exchangeable", scale.fix = TRUE, scale.value = 1)
###################################################
### code chunk number 10: ALDII-resp-glm-summary
###################################################
summary(resp_glm)
###################################################
### code chunk number 11: ALDII-resp-gee1summary
###################################################
summary(resp_gee1)
###################################################
### code chunk number 12: ALDII-resp-gee2-summary
###################################################
summary(resp_gee2)
###################################################
### code chunk number 13: ALDII-resp-confint
###################################################
se <- summary(resp_gee2)$coefficients["treatmenttreatment",
"Robust S.E."]
coef(resp_gee2)["treatmenttreatment"] +
c(-1, 1) * se * qnorm(0.975)
###################################################
### code chunk number 14: ALDII-resp-confint-exp
###################################################
exp(coef(resp_gee2)["treatmenttreatment"] +
c(-1, 1) * se * qnorm(0.975))
###################################################
### code chunk number 15: ALDII-epilepsy
###################################################
data("epilepsy", package = "HSAUR")
itp <- interaction(epilepsy$treatment, epilepsy$period)
tapply(epilepsy$seizure.rate, itp, mean)
tapply(epilepsy$seizure.rate, itp, var)
###################################################
### code chunk number 16: ALDII-plot1
###################################################
layout(matrix(1:2, nrow = 1))
ylim <- range(epilepsy$seizure.rate)
placebo <- subset(epilepsy, treatment == "placebo")
progabide <- subset(epilepsy, treatment == "Progabide")
boxplot(seizure.rate ~ period, data = placebo,
ylab = "Number of seizures",
xlab = "Period", ylim = ylim, main = "Placebo")
boxplot(seizure.rate ~ period, data = progabide,
main = "Progabide", ylab = "Number of seizures",
xlab = "Period", ylim = ylim)
###################################################
### code chunk number 17: ALDII-plot2
###################################################
layout(matrix(1:2, nrow = 1))
ylim <- range(log(epilepsy$seizure.rate + 1))
boxplot(log(seizure.rate + 1) ~ period, data = placebo,
main = "Placebo", ylab = "Log number of seizures",
xlab = "Period", ylim = ylim)
boxplot(log(seizure.rate + 1) ~ period, data = progabide,
main = "Progabide", ylab = "Log number of seizures",
xlab = "Period", ylim = ylim)
###################################################
### code chunk number 18: ALDII-epilepsy-gee
###################################################
per <- rep(log(2),nrow(epilepsy))
epilepsy$period <- as.numeric(epilepsy$period)
fm <- seizure.rate ~ base + age + treatment + offset(per)
epilepsy_glm <- glm(fm, data = epilepsy, family = "poisson")
epilepsy_gee1 <- gee(fm, data = epilepsy, family = "poisson",
id = subject, corstr = "independence", scale.fix = TRUE,
scale.value = 1)
epilepsy_gee2 <- gee(fm, data = epilepsy, family = "poisson",
id = subject, corstr = "exchangeable", scale.fix = TRUE,
scale.value = 1)
epilepsy_gee3 <- gee(fm, data = epilepsy, family = "poisson",
id = subject, corstr = "exchangeable", scale.fix = FALSE,
scale.value = 1)
###################################################
### code chunk number 19: ALDII-espilepsy-glm-summary
###################################################
summary(epilepsy_glm)
###################################################
### code chunk number 20: ALDII-espilepsy-gee1-summary
###################################################
summary(epilepsy_gee1)
###################################################
### code chunk number 21: ALDII-espilepsy-gee2-summary
###################################################
summary(epilepsy_gee2)
###################################################
### code chunk number 22: ALDII-espilepsy-gee3-summary
###################################################
summary(epilepsy_gee3)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.