-
Notifications
You must be signed in to change notification settings - Fork 0
/
chapter11.qmd
116 lines (82 loc) · 3.04 KB
/
chapter11.qmd
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
---
title: "Chapter 11 - Joint Analysis of Longitudinal and Survival Data"
---
## Slides
Lecture slides [here](chap11.html){target="_blank"}. (To convert html to pdf, press E $\to$ Print $\to$ Destination: Save to pdf)
## Base R Code
```{r}
#| code-fold: true
#| code-summary: "Show the code"
#| eval: false
##################################################################
# This code generates all numerical results in chapter 11. ##
##################################################################
############################################################
# Analysis of the anti-retroviral drug trial
############################################################
#load required packages
library(nlme) # for linear mixed effects model
library(survival)
# read in the study dataset
data <- read.table("Data//Anti-retroviral Trial//aids.txt")
head(data)
#load the JM package
install.packages("JM")
library(JM)
#####################################################
# Figure 11.2 Histograms of CD4 count and square-root
# transformation
#####################################################
par(mfrow=c(1,2))
hist(data$CD4,xlab="CD4 count", main = "", lwd=2)
hist(sqrt(data$CD4),xlab="Square root of CD4 count",main="",lwd=2)
# taking square root of CD4 count
data$y <- sqrt(data$CD4)
# create a de-duplicated data for survival sub-model
data.surv <- data[!duplicated(data$id),]
# Joint model fit for the HIV/AIDS dataset
# longitudinal sub-model
longit.sub <- lme(y ~ obsmo + obsmo:drug + sex + hist,
random = ~ obsmo|id, data = data)
# survival sub-model
surv.sub <- coxph(Surv(time, status) ~ drug + sex + hist,
data = data.surv, x = TRUE)
# combine the two models
# piecewise linear baseline with six internal knots
joint.model <- jointModel(longit.sub, surv.sub,
timeVar = "obsmo",
method = "piecewise-PH-aGH")
joint.model$coefficients
# print out the summary
summary(joint.model)
######################################
# Figure 11.3
######################################
# compute the mean trajectories based on the longitudinal
# sub-model results
t <- (0:180)/10
g0 <- 2.2123
g1 <- -0.0414
g3 <- 0.9153
g4 <- 0.0061
ddC.nonhist <- (g0+g3+g1*t)^2
ddC.hist <- (g0+g1*t)^2
ddI.nonhist <- (g0+g3+(g1+g4)*t)^2
ddI.hist <- (g0+(g1+g4)*t)^2
#######################################################
# Figure 11.3
# Model-based estimates of the mean trajectories
# of CD4 count for female patients by treatment
# group and previous infection status. Solid, ddC;
# dotted, ddI
######################################################
par(mfrow=c(1,2))
plot(t, ddC.nonhist,type='l',frame.plot=F,
main="No previous infection",xlim=c(0,20),ylim=c(0,10),
xlab="Time (months)", ylab="Mean CD4 cell count", lwd=2, cex.main = 1)
lines(t,ddI.nonhist,lty=3,lwd=2)
plot(t, ddC.hist,type='l',frame.plot=F,
main="Previous infection",xlim=c(0,20),ylim=c(0,10),
xlab="Time (months)", ylab="Mean CD4 cell count", lwd=2, cex.main = 1)
lines(t,ddI.hist,lty=3,lwd=2)
```