Skip to content

Commit ddbfbbd

Browse files
danielhp777danielhp777
authored
Add files via upload
1 parent db95c4b commit ddbfbbd

4 files changed

Lines changed: 738 additions & 0 deletions

File tree

‎EDA_code.R‎

Lines changed: 141 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,141 @@
1+
library(ggplot2)
2+
library(dplyr)
3+
data=read.table("Data2.txt",header=TRUE)
4+
head(data)
5+
6+
# plot of logrna as a whole
7+
ggplot(data, aes(x = calwk, y = logrna, group = patid, color = factor(trtarm))) +
8+
geom_line(size = 1) +
9+
geom_point(size = 1) +
10+
facet_wrap(~ trtarm, ncol = 4, labeller = label_both) + # 3 columns: one for each treatment arm
11+
labs(title = "Longitudinal Log RNA Viral Load by Treatment Arm",
12+
x = "Calendar Week", y = "Log10 RNA Viral Load",
13+
) +
14+
theme_minimal() + theme(legend.position = "none")
15+
16+
# separating data (6 measurements per subject for most subjects, m=481)
17+
time <- c(0,2,4,8,16,24)
18+
calwk <- data$calwk
19+
logrna <- as.numeric(data$logrna)
20+
nnrti <- factor(data$nnrti)
21+
txday <- data$txday
22+
cens <- factor(data$cens)
23+
trtarm <- factor(data$trtarm)
24+
25+
# averages for each trtarm
26+
mean_logrna <- aggregate(logrna ~ trtarm + calwk, data = data, FUN = mean, na.rm = TRUE)
27+
ggplot(mean_logrna, aes(x = calwk, y = logrna, color = factor(trtarm), group = trtarm)) +
28+
geom_line(size = 1.2) +
29+
geom_point(size = 2) +
30+
labs(title = "Mean log10 Viral Load Over Time by Treatment", x = "Week",
31+
y = "Mean log10 Viral Load",
32+
color = "Treatment Arm") +
33+
scale_color_discrete(labels =
34+
c("Saquinavir", "Indinavir", "Nelfinavir", "Placebo")) +
35+
theme_minimal()
36+
37+
# creating new data frames that only contain one treatment arm
38+
trt1 <- NULL
39+
trt2 <- NULL
40+
trt3 <- NULL
41+
trt4 <- NULL
42+
for(i in 1:length(data$patid)){
43+
if(data$trtarm[i]==1){
44+
trt1 <- rbind(trt1, data[i,])
45+
}
46+
if(data$trtarm[i]==2){
47+
trt2 <- rbind(trt2, data[i,])
48+
}
49+
if(data$trtarm[i]==3){
50+
trt3 <- rbind(trt3, data[i,])
51+
}
52+
if(data$trtarm[i]==4){
53+
trt4 <- rbind(trt4, data[i,])
54+
}
55+
}
56+
57+
# calculate the mean and sd logrna values at each time point for each treatment arm
58+
trt1_means <- aggregate(logrna ~ calwk, data = trt1, FUN = mean)
59+
trt2_means <- aggregate(logrna ~ calwk, data = trt2, FUN = mean)
60+
trt3_means <- aggregate(logrna ~ calwk, data = trt3, FUN = mean)
61+
trt4_means <- aggregate(logrna ~ calwk, data = trt4, FUN = mean)
62+
trt1_sd <- aggregate(logrna ~ calwk, data = trt1, FUN= sd)
63+
trt2_sd <- aggregate(logrna ~ calwk, data = trt2, FUN= sd)
64+
trt3_sd <- aggregate(logrna ~ calwk, data = trt3, FUN= sd)
65+
trt4_sd <- aggregate(logrna ~ calwk, data = trt4, FUN= sd)
66+
trt1_mean_sd <- cbind(trt1_means, trt1_sd[,2])
67+
trt2_mean_sd <- cbind(trt2_means, trt2_sd[,2])
68+
trt3_mean_sd <- cbind(trt3_means, trt3_sd[,2])
69+
trt4_mean_sd <- cbind(trt4_means, trt4_sd[,2])
70+
71+
# plotting each treatment arm data (not included in dissertation)
72+
ggplot(trt1_mean_sd, aes(x=calwk, y=logrna)) + geom_line(size=1) + geom_point(size=1.5) +
73+
geom_ribbon(aes(ymin = trt1_means[,2] - trt1_sd[,2],
74+
ymax = trt1_means[,2] + trt1_sd[,2]),
75+
alpha = 0.4, fill = "blue") +
76+
labs(title = 'Mean log-10 HIV RNA values for patients receiving saquinavir (treatment 1)') + theme_minimal()
77+
78+
ggplot(trt2_mean_sd, aes(x=calwk, y=logrna)) + geom_line(size=1) + geom_point(size=1.5) +
79+
geom_ribbon(aes(ymin = trt2_means[,2] - trt2_sd[,2],
80+
ymax = trt2_means[,2] + trt2_sd[,2]),
81+
alpha = 0.4, fill = "red") +
82+
labs(title = 'Mean log-10 HIV RNA values for patients receiving indinavir (treatment 2)') + theme_minimal()
83+
84+
ggplot(trt3_mean_sd, aes(x=calwk, y=logrna)) + geom_line(size=1) + geom_point(size=1.5) +
85+
geom_ribbon(aes(ymin = trt3_means[,2] - trt3_sd[,2],
86+
ymax = trt3_means[,2] + trt3_sd[,2]),
87+
alpha = 0.4, fill = "green") +
88+
labs(title = 'Mean log-10 HIV RNA values for patients receiving nelfinavir (treatment 3)') + theme_minimal()
89+
90+
ggplot(trt4_mean_sd, aes(x=calwk, y=logrna)) + geom_line(size=1) + geom_point(size=1.5) +
91+
geom_ribbon(aes(ymin = trt4_means[,2] - trt4_sd[,2],
92+
ymax = trt4_means[,2] + trt4_sd[,2]),
93+
alpha = 0.4, fill = "purple") +
94+
labs(title = 'Mean log-10 HIV RNA values for patients receiving the placebo (treatment 4)') + theme_minimal()
95+
96+
# Proportions of groups where HIV cured after week 24
97+
# saquinavir:
98+
cured1 <- sum(trt1$logrna < 3 & trt1$calwk == 24)
99+
infected1 <- sum(trt1$logrna >= 3 & trt1$calwk == 24)
100+
cured_percentage1 <- 100 * cured1/(cured1 + infected1)
101+
102+
# indinavir:
103+
cured2 <- sum(trt2$logrna < 3 & trt2$calwk == 24)
104+
infected2 <- sum(trt2$logrna >= 3 & trt2$calwk == 24)
105+
cured_percentage2 <- 100 * cured2/(cured2 + infected2)
106+
107+
# nelfanavir:
108+
cured3 <- sum(trt3$logrna < 3 & trt3$calwk == 24)
109+
infected3 <- sum(trt3$logrna >= 3 & trt3$calwk == 24)
110+
cured_percentage3 <- 100 * cured3/(cured3 + infected3)
111+
112+
# placebo:
113+
cured4 <- sum(trt4$logrna < 3 & trt4$calwk == 24)
114+
infected4 <- sum(trt4$logrna >= 3 & trt4$calwk == 24)
115+
cured_percentage4 <- 100 * cured4/(cured4 + infected4)
116+
117+
# bar chart of patients cured after 24 weeks
118+
df_barchart <- data.frame(Treatment = c("Saquinavir", "Saquinavir",
119+
"Indinavir", "Indinavir",
120+
"Nelfinavir", "Nelfinavir",
121+
"Placebo", "Placebo"),
122+
State = c("Cured", "Infected",
123+
"Cured", "Infected",
124+
"Cured", "Infected",
125+
"Cured", "Infected"),
126+
count = c(cured1, infected1,
127+
cured2, infected2,
128+
cured3, infected3,
129+
cured4, infected4),
130+
percentage_cured = c(cured_percentage1, NA,
131+
cured_percentage2, NA,
132+
cured_percentage3, NA,
133+
cured_percentage4, NA))
134+
ggplot(df_barchart, aes(x= Treatment, y = count, fill = State)) +
135+
geom_bar(stat = "identity") +
136+
labs(title = "Proportion of participants cured from HIV after 24 weeks",
137+
x = "Treatment", y = "Number of participants",
138+
fill = "Cured or Infected?") +
139+
geom_text(aes(label = ifelse(!is.na(percentage_cured), paste0(round(percentage_cured,1), "%"), "")),
140+
position = position_stack(vjust = 0.5)) +
141+
theme_minimal()

‎LMM _code.R‎

Lines changed: 79 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,79 @@
1+
library(nlme)
2+
library(dplyr)
3+
library(ggplot2)
4+
5+
# load and process data
6+
data <- read.table("Data2.txt",header=TRUE)
7+
8+
# separating data (6 measurements per subject for most subjects, m=481)
9+
time <- c(0,2,4,8,16,24)
10+
calwk <- data$calwk
11+
logrna <- as.numeric(data$logrna)
12+
nnrti <- factor(data$nnrti)
13+
txday <- data$txday
14+
cens <- factor(data$cens)
15+
trtarm <- factor(data$trtarm)
16+
patid <- data$patid
17+
18+
hiv <- data.frame(patid, nnrti, calwk, txday, logrna, cens, trtarm)
19+
20+
set.seed(158)
21+
patients <- seq(1:481)
22+
test_ids <- sample(patients, size = 96, replace = FALSE) # choosing the patient ID's for test data
23+
test_ids <- sort(test_ids)
24+
train_ids <- setdiff(patients, test_ids)
25+
26+
hiv_train <- hiv %>% filter(patid %in% train_ids)
27+
hiv_test <- hiv %>% filter(patid %in% test_ids)
28+
29+
# Linear Mixed-Effect Models: choosing the best LMM to use (random slope models)
30+
# independent within subject correlation structure
31+
lme1 <- lme(logrna ~ nnrti + calwk + txday + cens + trtarm,
32+
data = hiv_train, random = ~ calwk | patid, method = 'ML')
33+
summary(lme1)
34+
35+
# compound symmetry structure for within subject correlation
36+
lme2 <- update(lme1, correlation = corCompSymm(form = ~ 1 | patid))
37+
summary(lme2)
38+
39+
# fewer fixed effects, random slope with independent within subject correlation
40+
lme3 <- lme(logrna ~ trtarm * calwk, data = hiv_train,
41+
random = ~calwk | patid, method = 'ML')
42+
summary(lme3)
43+
44+
lme4 <- update(lme3, correlation = corCompSymm(form = ~ 1 | patid))
45+
summary(lme4)
46+
47+
anova(lme1, lme2, lme3, lme4) # we choose lme1 based on AIC/BIC
48+
49+
# improve parameter estimates by using REML instead of ML
50+
lme5 <- update(lme1, method = 'REML')
51+
summary(lme5) # shows the fitted parameters for the chosen LMM
52+
53+
# predictions
54+
lme5_pred <- predict(lme5, newdata = hiv_test, level = 0) # level asserts that we only use fixed effects in prediction, since expectation of random components is 0
55+
56+
# plot predictions vs actual responses
57+
D5 <- data.frame(real_values = hiv_test$logrna, Predictions = lme5_pred)
58+
59+
ggplot(D5, aes(x = Predictions, y = real_values)) +
60+
geom_point(alpha = 0.6, aes(color = 'Data points')) +
61+
geom_smooth(method = "lm", se = FALSE, aes(color = 'Regression line')) + # best fit line
62+
geom_line(linewidth = 1, data = data.frame(x = range(D5$Predictions), y = range(D5$Predictions)), aes(x=x, y=y, color = 'Identity line')) +
63+
scale_color_manual(name = 'Legend', values = c('Data points'= 'deeppink3', 'Regression line' = 'black', 'Identity line' = 'green')) +
64+
labs(
65+
x = "Predicted viral load", y = "Actual viral load",
66+
title =
67+
"Predicted vs Actual viral loads for the testing data (LMM)") +
68+
theme_minimal()
69+
70+
# mse
71+
mse_lmm <- mean((lme5_pred - hiv_test$logrna)^2)
72+
mse_lmm
73+
74+
# mae
75+
mae_lmm <- mean(abs(lme5_pred - hiv_test$logrna))
76+
mae_lmm
77+
78+
# correlation between predicted and actual responses
79+
cor(hiv_test$logrna, lme5_pred)

0 commit comments

Comments
 (0)