knitr::opts_chunk$set(echo = TRUE)
library(readr)
library(magrittr)
library(ggplot2)
library(dplyr)
library(knitr)
hiv_data <- read_csv(
"https://raw.githubusercontent.com/ComputationalModeling/IPML-Data/master/01HIVseries/HIVseries.csv",
col_names = c("time_in_days", "viral_load")
)
hiv_data %<>%
mutate(log_viral_load = log(viral_load)) %>%
mutate(index = rownames(hiv_data))
hiv_data %>% knitr::kable()
| time_in_days | viral_load | log_viral_load | index |
|---|---|---|---|
| 0.0000 | 106100.0 | 11.572137 | 1 |
| 0.0831 | 93240.0 | 11.442932 | 2 |
| 0.1465 | 166720.0 | 12.024071 | 3 |
| 0.2587 | 153780.0 | 11.943278 | 4 |
| 0.4828 | 118800.0 | 11.685197 | 5 |
| 0.7448 | 116900.0 | 11.669074 | 6 |
| 0.9817 | 109570.0 | 11.604319 | 7 |
| 1.2563 | 111350.0 | 11.620434 | 8 |
| 1.4926 | 74388.0 | 11.217050 | 9 |
| 1.7299 | 83291.0 | 11.330096 | 10 |
| 1.9915 | 66435.0 | 11.103979 | 11 |
| 3.0011 | 35408.0 | 10.474693 | 12 |
| 4.0109 | 21125.0 | 9.958213 | 13 |
| 5.0090 | 20450.0 | 9.925738 | 14 |
| 5.9943 | 15798.0 | 9.667639 | 15 |
| 7.0028 | 4785.2 | 8.473283 | 16 |
time <- seq(from = 0, to = 10, length.out = 101)
get_model <- function(A, alpha, B, beta, time) {
model <- (A * exp(-alpha * time)) + (B * exp(-beta * time))
return (model)
}
get_slope_and_intercept <- function(
x_values,
y_values,
index1,
index2)
{
rise <- y_values[index2] - y_values[index1]
run <- x_values[index2] - x_values[1]
slope <- rise / run
intercept <- y_values[index1] - (slope * x_values[index1])
output <- data.frame(
rise,
run,
slope,
intercept
)
return (output)
}
early_stage_clearance <- get_slope_and_intercept(
x_values = hiv_data$time_in_days,
y_values = hiv_data$log_viral_load,
index1 = 1,
index2 = 10
) %>%
data.frame(clearance_stage="early stage")
late_stage_clearance <- get_slope_and_intercept(
x_values = hiv_data$time_in_days,
y_values = hiv_data$log_viral_load,
index1 = 11,
index2 = 16
) %>%
data.frame(clearance_stage="later stage")
clearance_rates <- merge(
early_stage_clearance,
late_stage_clearance,
all = TRUE
)
clearance_rates %>% kable()
| rise | run | slope | intercept | clearance_stage |
|---|---|---|---|---|
| -2.6306962 | 7.0028 | -0.3756635 | 11.85211 | later stage |
| -0.2420415 | 1.7299 | -0.1399165 | 11.57214 | early stage |
We can try plotting the two distinct linear models
raw_data_plot <- ggplot()
raw_data_plot <- raw_data_plot + geom_text(
aes(
x = time_in_days,
y = log_viral_load,
label = index
),
data = hiv_data
)
full_line_plot <- geom_abline(
aes(
slope = slope,
intercept = intercept,
group = clearance_stage,
color = clearance_stage
),
data = clearance_rates
)
print(
raw_data_plot +
full_line_plot +
scale_color_brewer(type = "qual")
)
clearance_1 <- with(
hiv_data,
data.frame(
x = time_in_days[1],
xend = time_in_days[10],
y = log_viral_load[1],
yend = log_viral_load[10],
clearance_stage = "early"
)
)
clearance_2 <- with(
hiv_data,
data.frame(
x = time_in_days[11],
xend = time_in_days[16],
y = log_viral_load[11],
yend = log_viral_load[16],
clearance_stage = "late"
)
)
clearance_rates <- merge(
clearance_1,
clearance_2,
all = TRUE
)
line_segment_plot <- raw_data_plot
line_segment_plot <- line_segment_plot + geom_segment(
aes(
x = x,
xend = xend,
y = y,
yend = yend,
group = clearance_stage,
color = clearance_stage
),
data = clearance_rates
)
line_segment_plot <- line_segment_plot +
scale_color_brewer(type = "qual") +
theme_bw()
print(line_segment_plot)
)