-
Notifications
You must be signed in to change notification settings - Fork 1
/
R-bios6643-L12.Rmd
executable file
·139 lines (92 loc) · 4.61 KB
/
R-bios6643-L12.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
---
title: "BIOS6643. L12 GEE and Marginal Models"
cheauthor: "EJC"
-- date: ""
header-includes:
- \usepackage{amsmath}
- \usepackage{float}
output: pdf_document
---
\newcommand{\bi}{\begin{itemize}}
\newcommand{\ei}{\end{itemize}}
\newcommand{\itt}{\item}
```{r setup, include=FALSE, echo=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
<!-- -->
```{r installp, echo = FALSE, eval=T, include=F}
##library(dplyr)
##library(tidyverse)
##library(ggplot2)
library(gee)
library(geepack)
```
Fitting population-averaged models using the gee() function and the geeglm() function in geepack. We will need to install.packages("gee")
and install.packages("geepack").
The syntax for each is similar. A general call looks like
fit.object <- gee(model formula, id, corstr, family, data)
- **corstr** is a specification for the working correlation matrix
- **family** is a specification for the scaled exponential family that
would be relevant under independence; the canonical link is the default
# Dental data
## Example of continuous outcome
```{r eval=T, echo = T, include=TRUE}
## Read in the data
dat.den <- read.csv("/Users/juarezce/Documents/OneDrive - The University of Colorado Denver/BIOS6643/BIOS6643_Notes/data/dental.csv")
head(dat.den,3)
##dat.den$id <- as.factor(dat.den$id)
##dat.den$gender <- as.factor(dat.den$gender)
ind.gee.den <- gee(distance ~ age*gender,
id = id,
data = dat.den,
family = gaussian,
corstr = "independence"
)
summary(ind.gee.den)
```
Recall the interpretation of the parameter estimates is at the population level. For instance for boys, the $\hat{\beta}_{age}=0.78$ may be interpreted as the average (expected) increase in distance per year of age increased.
# Seizure data using gee()
## Epileptic Seizure Study of a randomized trial reported in Thall and Vail (1990).
- 59 subjects with epilepsy suffering from simple or partial seizures were assigned at random to receive either the anti-epileptic drug progabide or a placebo in addition to a standard chemotherapy regimen all were taking.
- Because each individual might be prone to different rates of experiencing seizures, the investigators first tried to get a sense of this by recording the number of seizures suffered by each subject over the 8-week period prior to the start of administration of the assigned treatment.
- It is common in such studies to record such baseline measurements, so that the effect of treatment for each subject can be measured relative to how that subject behaved prior to treatment.
- Following initiation of treatment, the number of seizures for each subject was counted for each of 4 consecutive 2-week periods.
- The age of each subject at the start of the study was also recorded, as it was suspected that subject age might be associated with the effect of the treatment.
```{r eval=T, echo = T, include=T}
# Read in the data
dat.sz <- read.table("/Users/juarezce/Documents/OneDrive - The University of Colorado Denver/BIOS6643/BIOS6643_Notes/data/epilepsy.txt")
colnames(dat.sz) <- c("subj","seize","visit","trt","base","age")
## trt=0 corresponds to placebo
head(dat.sz,3)
# Create other covariates
dat.sz$o <- 8*(dat.sz$visit==0)+2*(dat.sz$visit>0)
dat.sz$logo <- log(dat.sz$o)
dat.sz$vm0 <- as.numeric(dat.sz$visit>0)
# Basic models -- the unstructured fit using gee() does not converge
# even with maxiter set to be much larger than the default
##un.gee <- gee(seize ~ trt + offset(logo), id=subj, family=poisson,
## corstr="unstructured", data=dat.sz, maxiter=100)
cs.gee <- gee(seize ~ trt + offset(logo),id=subj,family=poisson,
corstr="exchangeable", data=dat.sz)
summary(cs.gee)
## AR1
ar1.gee <- gee(seize ~ trt + offset(logo),id=subj,family=poisson,
corstr="AR-M",Mv=1,data=dat.sz)
summary(ar1.gee)
```
# Seizure data using geeglm()
```{r eval=T, echo = T, include=TRUE}
## UN using geeglm
un.geeglm <- geeglm(seize ~ trt, id=subj,family=poisson("log"),
offset=logo, corstr="unstructured",data=dat.sz)
summary(un.geeglm)
## CS using geeglm
cs.geeglm <- geeglm(seize ~ trt,id=subj,family=poisson("log"),
offset=logo, corstr="exchangeable",data=dat.sz)
summary(cs.geeglm)
ar1.geeglm <- geeglm(seize ~ trt,id=subj,family=poisson("log"),
offset=logo, corstr="ar1",data=dat.sz)
summary(ar1.geeglm)
```
# Investigate if there is a different different effect after baseline visit 0
This means include an interaction between the indicator variable for visit>0 and the treatment indicator. Use two different working correlation structures and compare the results.