ICT case study: Phenobarbital, groep 6
Pharmacokinetic data analysis
Ilse Jansen, Natalia Khan, Parsa Shameli, Laurens Bakker
Introduction
Phenobarbital is a benzodiazepine that is used to prevent convulsions in newborns. Pediatricians wonder whether the pharmacokinetics of phenobarbital in newborns that are treated with therapeutic hypothermia (TH) are the same for normal (normo-thermic) newborns. They have concentration (mg/L) – time (h) data available from patients that were treated with TH (=1) and patients of similar age that are not treated with TH (=0). All patients received a dose of 40 mg in this study.
Notes:
- Make sure that for all questions below, the R-code required to obtain the answers are shown (you can use code folding however)
- Knit this file, so that you get a more readable case study
- You can use this file to fill in your answers, don’t forget to change the authors at the top.
- Write a short (1 paragraph) conclusion at the end, where you summarize your findings
- Knit this file to HTML to upload to Brightspace. Other formats will not be graded! (Tip: Knit early and often to prevent knitting problems just before the deadline)
Exercises
Exercise 1
Read in the dataset Phenobarbital_TH.csv and describe
the design of the study in terms of number of individuals, treatments in
the study arms and number and timing of PK samples. Show the table using
DT (https://rstudio.github.io/DT/).
NB: Don’t change the file; all changes should be done
with R code
Solution
This study was done with 30 individuals, 15 in each treatment arm. The first treatment arm (TH = 1) consist of newborns with therapeutic hypothermia, while the second treatment arm (TH = 0) consist of newborns without this condition. Both groups of patients received a dose of 40 mg of phenobarbital. The concentration (mg/mL) was measured at several timepoints: 0h, 0.5h, 2h, 4h, 8h, 12h, 24h and 48h. This gives 7 observations per patient, which adds up to a total of 210 observations.
# write your code here
data_phen <- read.csv("Phenobarbital_TH.csv")
data_phen## ID TIME DV TH
## 1 2017005 0.0 0.0000 1
## 2 2017005 0.5 32.9500 1
## 3 2017005 2.0 28.3500 1
## 4 2017005 4.0 27.4400 1
## 5 2017005 8.0 13.7400 1
## 6 2017005 24.0 5.1400 1
## 7 2017005 48.0 0.4900 1
## 8 2017043 0.0 0.0000 1
## 9 2017043 0.5 49.5200 1
## 10 2017043 2.0 38.1100 1
## 11 2017043 4.0 36.5100 1
## 12 2017043 8.0 29.8400 1
## 13 2017043 24.0 7.3800 1
## 14 2017043 48.0 2.8900 1
## 15 2017057 0.0 0.0000 1
## 16 2017057 0.5 32.0700 1
## 17 2017057 2.0 35.0900 1
## 18 2017057 4.0 25.6500 1
## 19 2017057 8.0 17.3300 1
## 20 2017057 24.0 7.5800 1
## 21 2017057 48.0 1.3400 1
## 22 2017107 0.0 0.0000 1
## 23 2017107 0.5 38.4600 1
## 24 2017107 2.0 34.6300 1
## 25 2017107 4.0 30.1600 1
## 26 2017107 8.0 29.6200 1
## 27 2017107 24.0 27.8100 1
## 28 2017107 48.0 17.1900 1
## 29 2017154 0.0 0.0000 1
## 30 2017154 0.5 39.1600 1
## 31 2017154 2.0 31.4400 1
## 32 2017154 4.0 35.2600 1
## 33 2017154 8.0 50.3400 1
## 34 2017154 24.0 10.9800 1
## 35 2017154 48.0 4.9800 1
## 36 2017159 0.0 0.0000 1
## 37 2017159 0.5 36.3200 1
## 38 2017159 2.0 34.3200 1
## 39 2017159 4.0 25.9200 1
## 40 2017159 8.0 22.9400 1
## 41 2017159 24.0 5.6200 1
## 42 2017159 48.0 11.7600 1
## 43 2017261 0.0 0.0000 1
## 44 2017261 0.5 42.3100 1
## 45 2017261 2.0 42.0800 1
## 46 2017261 4.0 35.2100 1
## 47 2017261 8.0 32.2600 1
## 48 2017261 24.0 18.7400 1
## 49 2017261 48.0 4.5700 1
## 50 2018004 0.0 0.0000 1
## 51 2018004 0.5 46.2500 1
## 52 2018004 2.0 42.2100 1
## 53 2018004 4.0 39.6500 1
## 54 2018004 8.0 28.7200 1
## 55 2018004 24.0 12.3600 1
## 56 2018004 48.0 4.9500 1
## 57 2018048 0.0 0.0000 1
## 58 2018048 0.5 49.6700 1
## 59 2018048 2.0 36.8800 1
## 60 2018048 4.0 35.4200 1
## 61 2018048 8.0 20.1900 1
## 62 2018048 24.0 9.9500 1
## 63 2018048 48.0 1.4600 1
## 64 2018058 0.0 0.0000 1
## 65 2018058 0.5 41.3800 1
## 66 2018058 2.0 36.0700 1
## 67 2018058 4.0 40.0300 1
## 68 2018058 8.0 32.4900 1
## 69 2018058 24.0 23.3200 1
## 70 2018058 48.0 9.2100 1
## 71 2018114 0.0 0.0000 1
## 72 2018114 0.5 45.8000 1
## 73 2018114 2.0 40.2100 1
## 74 2018114 4.0 39.9500 1
## 75 2018114 8.0 33.1800 1
## 76 2018114 24.0 19.8100 1
## 77 2018114 48.0 12.7000 1
## 78 2018143 0.0 0.0000 1
## 79 2018143 0.5 37.8700 1
## 80 2018143 2.0 41.0300 1
## 81 2018143 4.0 32.4200 1
## 82 2018143 8.0 21.8300 1
## 83 2018143 24.0 9.4100 1
## 84 2018143 48.0 1.3300 1
## 85 2018181 0.0 0.0000 1
## 86 2018181 0.5 44.1000 1
## 87 2018181 2.0 34.1700 1
## 88 2018181 4.0 27.4800 1
## 89 2018181 8.0 25.8400 1
## 90 2018181 24.0 10.6200 1
## 91 2018181 48.0 3.3200 1
## 92 2018203 0.0 0.0000 1
## 93 2018203 0.5 31.7300 1
## 94 2018203 2.0 37.0000 1
## 95 2018203 4.0 31.5300 1
## 96 2018203 8.0 30.7400 1
## 97 2018203 24.0 14.4400 1
## 98 2018203 48.0 3.8000 1
## 99 2018221 0.0 0.0000 1
## 100 2018221 0.5 50.2200 1
## 101 2018221 2.0 49.1200 1
## 102 2018221 4.0 39.7300 1
## 103 2018221 8.0 27.9900 1
## 104 2018221 24.0 16.5500 1
## 105 2018221 48.0 5.0300 1
## 106 2017003 0.0 0.0000 0
## 107 2017003 1.0 32.9500 0
## 108 2017003 2.0 28.3500 0
## 109 2017003 6.0 17.4400 0
## 110 2017003 12.0 10.7400 0
## 111 2017003 24.0 4.1400 0
## 112 2017003 48.0 0.3240 0
## 113 2017012 0.0 0.0000 0
## 114 2017012 1.0 49.5400 0
## 115 2017012 2.0 38.1100 0
## 116 2017012 6.0 33.4100 0
## 117 2017012 12.0 19.8400 0
## 118 2017012 24.0 5.3800 0
## 119 2017012 48.0 2.0900 0
## 120 2017062 0.0 0.0000 0
## 121 2017062 1.0 35.0700 0
## 122 2017062 2.0 32.0900 0
## 123 2017062 6.0 22.6500 0
## 124 2017062 12.0 15.3300 0
## 125 2017062 24.0 2.5800 0
## 126 2017062 48.0 1.1400 0
## 127 2017064 0.0 0.0000 0
## 128 2017064 1.0 38.4600 0
## 129 2017064 2.0 34.6300 0
## 130 2017064 6.0 30.1600 0
## 131 2017064 12.0 19.6200 0
## 132 2017064 24.0 13.8100 0
## 133 2017064 48.0 3.1900 0
## 134 2017096 0.0 0.0000 0
## 135 2017096 1.0 39.1600 0
## 136 2017096 2.0 31.4400 0
## 137 2017096 6.0 28.2600 0
## 138 2017096 12.0 10.3400 0
## 139 2017096 24.0 11.9800 0
## 140 2017096 48.0 1.7800 0
## 141 2017133 0.0 0.0000 0
## 142 2017133 1.0 36.3200 0
## 143 2017133 2.0 30.3200 0
## 144 2017133 6.0 24.9200 0
## 145 2017133 12.0 18.9400 0
## 146 2017133 24.0 5.6200 0
## 147 2017133 48.0 0.7600 0
## 148 2017201 0.0 0.0000 0
## 149 2017201 1.0 42.0100 0
## 150 2017201 2.0 32.5100 0
## 151 2017201 6.0 25.2100 0
## 152 2017201 12.0 15.2600 0
## 153 2017201 24.0 3.8700 0
## 154 2017201 48.0 0.4570 0
## 155 2017211 0.0 0.0000 0
## 156 2017211 1.0 43.5400 0
## 157 2017211 2.0 36.1400 0
## 158 2017211 6.0 26.6500 0
## 159 2017211 12.0 12.7500 0
## 160 2017211 24.0 2.3600 0
## 161 2017211 48.0 0.1900 0
## 162 2018018 0.0 0.0000 0
## 163 2018018 1.0 50.6700 0
## 164 2018018 2.0 38.6800 0
## 165 2018018 6.0 34.2500 0
## 166 2018018 12.0 19.1900 0
## 167 2018018 24.0 9.8500 0
## 168 2018018 48.0 1.3700 0
## 169 2018032 0.0 0.0000 0
## 170 2018032 1.0 40.8900 0
## 171 2018032 2.0 32.7100 0
## 172 2018032 6.0 24.3500 0
## 173 2018032 12.0 12.4900 0
## 174 2018032 24.0 8.3200 0
## 175 2018032 48.0 0.5180 0
## 176 2018097 0.0 0.0000 0
## 177 2018097 1.0 44.8000 0
## 178 2018097 2.0 39.2400 0
## 179 2018097 6.0 36.5400 0
## 180 2018097 12.0 20.8500 0
## 181 2018097 24.0 9.8300 0
## 182 2018097 48.0 0.7210 0
## 183 2018106 0.0 0.0000 0
## 184 2018106 1.0 42.3100 0
## 185 2018106 2.0 37.8200 0
## 186 2018106 6.0 25.2100 0
## 187 2018106 12.0 12.2600 0
## 188 2018106 24.0 6.7400 0
## 189 2018106 48.0 0.4570 0
## 190 2018135 0.0 0.0000 0
## 191 2018135 1.0 45.3400 0
## 192 2018135 2.0 66.1400 0
## 193 2018135 6.0 22.6500 0
## 194 2018135 12.0 13.7200 0
## 195 2018135 24.0 1.3600 0
## 196 2018135 48.0 0.0937 0
## 197 2018147 0.0 0.0000 0
## 198 2018147 1.0 40.2400 0
## 199 2018147 2.0 34.7300 0
## 200 2018147 6.0 27.5300 0
## 201 2018147 12.0 17.0740 0
## 202 2018147 24.0 7.4400 0
## 203 2018147 48.0 0.2780 0
## 204 2018178 0.0 0.0000 0
## 205 2018178 1.0 47.0200 0
## 206 2018178 2.0 43.6590 0
## 207 2018178 6.0 33.7300 0
## 208 2018178 12.0 17.9900 0
## 209 2018178 24.0 8.6500 0
## 210 2018178 48.0 1.0300 0library(DT)
datatable(data_phen)Exercise 2
How many unique individuals are included in the data set for both treatment arms?
Solution
In total there are 30 unique patients for both treatment arms thus 15 patients per study arm.
# write your code here
ID_phen <- data_phen %>% group_by(TH) %>% count(ID)
ID_phen## # A tibble: 30 × 3
## # Groups: TH [2]
## TH ID n
## <int> <int> <int>
## 1 0 2017003 7
## 2 0 2017012 7
## 3 0 2017062 7
## 4 0 2017064 7
## 5 0 2017096 7
## 6 0 2017133 7
## 7 0 2017201 7
## 8 0 2017211 7
## 9 0 2018018 7
## 10 0 2018032 7
## # ℹ 20 more rowsnrow(ID_phen)## [1] 30Exercise 3
Are there any differences between the pharmacokinetic profiles of phenobarbital patients treated with TH and patients that are not treated with TH? Base your assessment on a plot of the concentration-time data that includes all patients and uses different colors for patients with and without TH. Do this on a linear scale and semi-logarithmic scale.
Solution
The concentration-time profiles suggest that TH-treated patients have slightly slower elimination of the drug than normothermic patients, while peak concentrations appear similar. This indicates potential differences in clearance but similar distribution volume
TH_1_data <- data_phen %>%
group_by(TH == 1)
TH_0_data <- data_phen %>%
group_by(TH == 0)# write your code here
graph_norm <- ggplot(data_phen, aes(x = TIME, y = DV, group = ID, color = factor(TH), linetype = factor(TH))) + geom_line() + geom_point() + labs(
title = "Pharmacokinetic profiles of phenobarbital patients",
x = "Time (hours)",
y = "Concentration (mg/L)",
color = "Treatment (TH)",
linetype = "TH"
)
graph_log <- ggplot(data_phen, aes(x = TIME, y = DV, group = ID, color = factor(TH), linetype = factor(TH))) +
geom_line() + geom_point() + scale_y_log10() + labs(
title = "Pharmacokinetic profiles of phenobarbital patients on semi log scale",
x = "Time (hours)",
y = "Concentration (mg/L)",
color = "Treatment (TH)",
linetype = "TH"
)
graph_norm
graph_log## Warning: Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
# Linear scale interactive plot
graph_norm_plotly <- ggplotly(graph_norm)
graph_norm_plotly# Semi-logarithmic scale interactive plot
graph_log_plotly <- ggplotly(graph_log)## Warning: Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axisgraph_log_plotlyExercise 4
Can you detect any outlying data points? Make sure that you can differentiate ID numbers and treatment arms in the visualization. Do this on a linear scale and semi-logarithmic scale.
Hint: Fancy statistical outlier analysis not needed, with the correct visualization you should be able to identify the outliers by eye (also don’t worry too much about getting the correct ones, this question and the next are more about the removal and traceability than the actual points, however, you should explain why you think the points you removed are outliers)
outliers: 2017159, 48h, TH=1 → increases from 24h to 48h
2017154, 8h, TH=1 → very high value
2018135, 2h, TH=0 → very high value
2017096, 24h, TH=0 → increases from 12h to 48h
By visual inspection, the following points were identified as outliers: 2017159 (48h, TH=1), 2017154 (8h, TH=1), 2018135 (2h, TH=0), and 2017096 (24h, TH=0). These points either increase unexpectedly or have extreme values compared to the rest of the dataset.
Solution
Exercise 5
Make a subset of the dataset that does not include the outlying data points and use that dataset for the following questions. Make sure that it is clear in your script what subset of the data you have used (or what data points you have taken out), because traceability is very important!
Solution
data_phen_clean <- data_phen %>%
filter(!( (ID == 2017159 & TIME == 48 & TH == 1) |
(ID == 2017154 & TIME == 8 & TH == 1) |
(ID == 2018135 & TIME == 2 & TH == 0) |
(ID == 2017096 & TIME == 24 & TH == 0) )) %>%
filter(!(DV == 0))
data_phen_clean## ID TIME DV TH
## 1 2017005 0.5 32.9500 1
## 2 2017005 2.0 28.3500 1
## 3 2017005 4.0 27.4400 1
## 4 2017005 8.0 13.7400 1
## 5 2017005 24.0 5.1400 1
## 6 2017005 48.0 0.4900 1
## 7 2017043 0.5 49.5200 1
## 8 2017043 2.0 38.1100 1
## 9 2017043 4.0 36.5100 1
## 10 2017043 8.0 29.8400 1
## 11 2017043 24.0 7.3800 1
## 12 2017043 48.0 2.8900 1
## 13 2017057 0.5 32.0700 1
## 14 2017057 2.0 35.0900 1
## 15 2017057 4.0 25.6500 1
## 16 2017057 8.0 17.3300 1
## 17 2017057 24.0 7.5800 1
## 18 2017057 48.0 1.3400 1
## 19 2017107 0.5 38.4600 1
## 20 2017107 2.0 34.6300 1
## 21 2017107 4.0 30.1600 1
## 22 2017107 8.0 29.6200 1
## 23 2017107 24.0 27.8100 1
## 24 2017107 48.0 17.1900 1
## 25 2017154 0.5 39.1600 1
## 26 2017154 2.0 31.4400 1
## 27 2017154 4.0 35.2600 1
## 28 2017154 24.0 10.9800 1
## 29 2017154 48.0 4.9800 1
## 30 2017159 0.5 36.3200 1
## 31 2017159 2.0 34.3200 1
## 32 2017159 4.0 25.9200 1
## 33 2017159 8.0 22.9400 1
## 34 2017159 24.0 5.6200 1
## 35 2017261 0.5 42.3100 1
## 36 2017261 2.0 42.0800 1
## 37 2017261 4.0 35.2100 1
## 38 2017261 8.0 32.2600 1
## 39 2017261 24.0 18.7400 1
## 40 2017261 48.0 4.5700 1
## 41 2018004 0.5 46.2500 1
## 42 2018004 2.0 42.2100 1
## 43 2018004 4.0 39.6500 1
## 44 2018004 8.0 28.7200 1
## 45 2018004 24.0 12.3600 1
## 46 2018004 48.0 4.9500 1
## 47 2018048 0.5 49.6700 1
## 48 2018048 2.0 36.8800 1
## 49 2018048 4.0 35.4200 1
## 50 2018048 8.0 20.1900 1
## 51 2018048 24.0 9.9500 1
## 52 2018048 48.0 1.4600 1
## 53 2018058 0.5 41.3800 1
## 54 2018058 2.0 36.0700 1
## 55 2018058 4.0 40.0300 1
## 56 2018058 8.0 32.4900 1
## 57 2018058 24.0 23.3200 1
## 58 2018058 48.0 9.2100 1
## 59 2018114 0.5 45.8000 1
## 60 2018114 2.0 40.2100 1
## 61 2018114 4.0 39.9500 1
## 62 2018114 8.0 33.1800 1
## 63 2018114 24.0 19.8100 1
## 64 2018114 48.0 12.7000 1
## 65 2018143 0.5 37.8700 1
## 66 2018143 2.0 41.0300 1
## 67 2018143 4.0 32.4200 1
## 68 2018143 8.0 21.8300 1
## 69 2018143 24.0 9.4100 1
## 70 2018143 48.0 1.3300 1
## 71 2018181 0.5 44.1000 1
## 72 2018181 2.0 34.1700 1
## 73 2018181 4.0 27.4800 1
## 74 2018181 8.0 25.8400 1
## 75 2018181 24.0 10.6200 1
## 76 2018181 48.0 3.3200 1
## 77 2018203 0.5 31.7300 1
## 78 2018203 2.0 37.0000 1
## 79 2018203 4.0 31.5300 1
## 80 2018203 8.0 30.7400 1
## 81 2018203 24.0 14.4400 1
## 82 2018203 48.0 3.8000 1
## 83 2018221 0.5 50.2200 1
## 84 2018221 2.0 49.1200 1
## 85 2018221 4.0 39.7300 1
## 86 2018221 8.0 27.9900 1
## 87 2018221 24.0 16.5500 1
## 88 2018221 48.0 5.0300 1
## 89 2017003 1.0 32.9500 0
## 90 2017003 2.0 28.3500 0
## 91 2017003 6.0 17.4400 0
## 92 2017003 12.0 10.7400 0
## 93 2017003 24.0 4.1400 0
## 94 2017003 48.0 0.3240 0
## 95 2017012 1.0 49.5400 0
## 96 2017012 2.0 38.1100 0
## 97 2017012 6.0 33.4100 0
## 98 2017012 12.0 19.8400 0
## 99 2017012 24.0 5.3800 0
## 100 2017012 48.0 2.0900 0
## 101 2017062 1.0 35.0700 0
## 102 2017062 2.0 32.0900 0
## 103 2017062 6.0 22.6500 0
## 104 2017062 12.0 15.3300 0
## 105 2017062 24.0 2.5800 0
## 106 2017062 48.0 1.1400 0
## 107 2017064 1.0 38.4600 0
## 108 2017064 2.0 34.6300 0
## 109 2017064 6.0 30.1600 0
## 110 2017064 12.0 19.6200 0
## 111 2017064 24.0 13.8100 0
## 112 2017064 48.0 3.1900 0
## 113 2017096 1.0 39.1600 0
## 114 2017096 2.0 31.4400 0
## 115 2017096 6.0 28.2600 0
## 116 2017096 12.0 10.3400 0
## 117 2017096 48.0 1.7800 0
## 118 2017133 1.0 36.3200 0
## 119 2017133 2.0 30.3200 0
## 120 2017133 6.0 24.9200 0
## 121 2017133 12.0 18.9400 0
## 122 2017133 24.0 5.6200 0
## 123 2017133 48.0 0.7600 0
## 124 2017201 1.0 42.0100 0
## 125 2017201 2.0 32.5100 0
## 126 2017201 6.0 25.2100 0
## 127 2017201 12.0 15.2600 0
## 128 2017201 24.0 3.8700 0
## 129 2017201 48.0 0.4570 0
## 130 2017211 1.0 43.5400 0
## 131 2017211 2.0 36.1400 0
## 132 2017211 6.0 26.6500 0
## 133 2017211 12.0 12.7500 0
## 134 2017211 24.0 2.3600 0
## 135 2017211 48.0 0.1900 0
## 136 2018018 1.0 50.6700 0
## 137 2018018 2.0 38.6800 0
## 138 2018018 6.0 34.2500 0
## 139 2018018 12.0 19.1900 0
## 140 2018018 24.0 9.8500 0
## 141 2018018 48.0 1.3700 0
## 142 2018032 1.0 40.8900 0
## 143 2018032 2.0 32.7100 0
## 144 2018032 6.0 24.3500 0
## 145 2018032 12.0 12.4900 0
## 146 2018032 24.0 8.3200 0
## 147 2018032 48.0 0.5180 0
## 148 2018097 1.0 44.8000 0
## 149 2018097 2.0 39.2400 0
## 150 2018097 6.0 36.5400 0
## 151 2018097 12.0 20.8500 0
## 152 2018097 24.0 9.8300 0
## 153 2018097 48.0 0.7210 0
## 154 2018106 1.0 42.3100 0
## 155 2018106 2.0 37.8200 0
## 156 2018106 6.0 25.2100 0
## 157 2018106 12.0 12.2600 0
## 158 2018106 24.0 6.7400 0
## 159 2018106 48.0 0.4570 0
## 160 2018135 1.0 45.3400 0
## 161 2018135 6.0 22.6500 0
## 162 2018135 12.0 13.7200 0
## 163 2018135 24.0 1.3600 0
## 164 2018135 48.0 0.0937 0
## 165 2018147 1.0 40.2400 0
## 166 2018147 2.0 34.7300 0
## 167 2018147 6.0 27.5300 0
## 168 2018147 12.0 17.0740 0
## 169 2018147 24.0 7.4400 0
## 170 2018147 48.0 0.2780 0
## 171 2018178 1.0 47.0200 0
## 172 2018178 2.0 43.6590 0
## 173 2018178 6.0 33.7300 0
## 174 2018178 12.0 17.9900 0
## 175 2018178 24.0 8.6500 0
## 176 2018178 48.0 1.0300 0ggplot(data_phen_clean, aes(x = TIME, y = DV, group = ID, color = factor(TH))) +
geom_line() + geom_point() +
labs(
title = "Pharmacokinetic profiles of phenobarbital patients without outliers",
x = "Time (hours)",
y = "Concentration (mg/L)",
color = "Treatment (TH)",
linetype = "TH"
)
ggplot(data_phen_clean, aes(x = TIME, y = DV, group = ID, color = factor(TH))) +
geom_line() + geom_point() + scale_y_log10() +
labs(
title = "Pharmacokinetic profiles of phenobarbital patients without outliers on semi log scale",
x = "Time (hours)",
y = "ln(Concentration (mg/L)",
color = "Treatment (TH)",
linetype = "TH"
)
Exercise 6
What are the average Cmax (peak concentration), average Cmin (trough concentration), and average Cav (average concentration) for both treatment arms? Does it make sense to compare these values between both groups given the study design?
Hint: You need to summarize the output from summarize. (First summarize is per patient, second is summarizing over the patients)
Solution
These values provide a first impression of concentration differences between the TH and non-TH groups, but due to limited sampling and study design, direct comparison may not fully reflect differences in drug distribution or clearance.
# write your code here
library(dplyr)
pk_per_patient <- data_phen_clean %>%
group_by(ID, TH) %>%
summarise(
Cmax = max(DV),
Cmin = min(DV),
Cav = mean(DV),
)## `summarise()` has grouped output by 'ID'. You can override using the `.groups`
## argument.pk_per_patient## # A tibble: 30 × 5
## # Groups: ID [30]
## ID TH Cmax Cmin Cav
## <int> <int> <dbl> <dbl> <dbl>
## 1 2017003 0 33.0 0.324 15.7
## 2 2017005 1 33.0 0.49 18.0
## 3 2017012 0 49.5 2.09 24.7
## 4 2017043 1 49.5 2.89 27.4
## 5 2017057 1 35.1 1.34 19.8
## 6 2017062 0 35.1 1.14 18.1
## 7 2017064 0 38.5 3.19 23.3
## 8 2017096 0 39.2 1.78 22.2
## 9 2017107 1 38.5 17.2 29.6
## 10 2017133 0 36.3 0.76 19.5
## # ℹ 20 more rows pk_per_arm <- pk_per_patient %>%
group_by(TH) %>%
summarise(
avg_Cmax = mean(Cmax),
avg_Cmin = mean(Cmin),
avg_Cav = mean(Cav),
)
pk_per_arm## # A tibble: 2 × 4
## TH avg_Cmax avg_Cmin avg_Cav
## <int> <dbl> <dbl> <dbl>
## 1 0 41.9 0.960 21.2
## 2 1 42.0 5.26 26.3Exercise 7
The graph of concentration-time data reveals a linear profile on a semi-logarithmic scale. In order to perform a linear regression on the data, add a new column of (natural) log-transformed values of the concentrations to the data set.
Solution
# write your code here
data_lnDV <- data_phen_clean %>%
mutate(DV_ln = log(DV))
DV_ln_graph <- ggplot(data_lnDV, aes(x = TIME, y= DV, group = ID, color=factor(TH))) + geom_line() +
labs(
title = "Pharmacokinetic profiles of phenobarbital patients on ln scale",
x = "Time (hours)",
y = "ln(Concentration (mg/L)",
color = "Treatment (TH)",
linetype = "TH"
)
DV_ln_graph
Exercise 8
Fit a linear model to the ln-transformed concentration data, using both time and treatment as covariates. Make models with and without interaction and explain which model is better.
The peak concentration is mainly determined by the distribution volume of phenobarbital in this population, while the slope of the profiles is mainly determined by the clearance. Do you anticipate differences in distribution volume and clearance between patients on TH treatment and patients that are not on TH treatment, given the outcome of the model fit?
Hint: Which model should you use for your anticipation?
Bonus: Using the data_grid() and
gather_predictions() functions from the modelr package, can
you show the difference between the (predictions of the) two models in a
plot?
Solution
The model which takes interaction into account should be used since there is a significance shown in comparsion to the model which does not take interaction into account.
# write your code here
library(dplyr)
library(modelr)
# Zonder interaction
lm_add <- lm(DV_ln ~ TIME + TH, data = data_lnDV)
# Met interaction
lm_int <- lm(DV_ln ~ TIME * TH, data = data_lnDV)
summary(lm_add)##
## Call:
## lm(formula = DV_ln ~ TIME + TH, data = data_lnDV)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.55269 -0.26527 0.00459 0.27912 2.20259
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.499811 0.071485 48.959 < 2e-16 ***
## TIME -0.069058 0.002535 -27.241 < 2e-16 ***
## TH 0.456700 0.084399 5.411 2.06e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5593 on 173 degrees of freedom
## Multiple R-squared: 0.8195, Adjusted R-squared: 0.8175
## F-statistic: 392.8 on 2 and 173 DF, p-value: < 2.2e-16summary(lm_int)##
## Call:
## lm(formula = DV_ln ~ TIME * TH, data = data_lnDV)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.9713 -0.1687 0.0199 0.1887 1.5864
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.788327 0.068503 55.302 < 2e-16 ***
## TIME -0.087604 0.003016 -29.045 < 2e-16 ***
## TH -0.088320 0.094650 -0.933 0.352
## TIME:TH 0.036728 0.004244 8.653 3.46e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4682 on 172 degrees of freedom
## Multiple R-squared: 0.8743, Adjusted R-squared: 0.8721
## F-statistic: 398.7 on 3 and 172 DF, p-value: < 2.2e-16anova(lm_add, lm_int)## Analysis of Variance Table
##
## Model 1: DV_ln ~ TIME + TH
## Model 2: DV_ln ~ TIME * TH
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 173 54.119
## 2 172 37.705 1 16.414 74.877 3.46e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1grid <- data_lnDV %>%
data_grid(TIME = seq(min(TIME), max(TIME), length.out = 100),
TH = unique(TH)) %>%
gather_predictions(lm_add, lm_int)
ggplot(data_lnDV, aes(x = TIME, y = DV_ln, color = factor(TH))) +
geom_point(alpha = 0.5) +
geom_line(data = grid, aes(y = pred, linetype = model), size = 1) +
labs(y = "ln(Concentration)", color = "TH", linetype = "Model") +
labs(
title = "Linear Regression Analysis of ln-Transformed Phenobarbital Concentrations: Assessing the Effect of Therapeutic Hypothermia",
x = "Time (hours)",
y = "ln(Concentration (mg/L)",
color = "Treatment (TH)",
linetype = "TH"
)## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Exercise 9
Fit a linear model to the ln-transformed concentration data of all individual patients in both treatment arms and make a table that includes ID number, treatment arm and the concentration at time = 0 (time0).
Hints:
- Have a look at rowwise (https://dplyr.tidyverse.org/articles/rowwise.html#modelling-1) or at https://purrr.tidyverse.org/ and the slides
- Depending on your approach you might need to make a tibble out of the results, and join this with the original tibble
Solution
# write your code here
per_patient_models <- data_lnDV %>%
group_by(ID, TH) %>%
nest()
per_patient_models## # A tibble: 30 × 3
## # Groups: ID, TH [30]
## ID TH data
## <int> <int> <list>
## 1 2017005 1 <tibble [6 × 3]>
## 2 2017043 1 <tibble [6 × 3]>
## 3 2017057 1 <tibble [6 × 3]>
## 4 2017107 1 <tibble [6 × 3]>
## 5 2017154 1 <tibble [5 × 3]>
## 6 2017159 1 <tibble [5 × 3]>
## 7 2017261 1 <tibble [6 × 3]>
## 8 2018004 1 <tibble [6 × 3]>
## 9 2018048 1 <tibble [6 × 3]>
## 10 2018058 1 <tibble [6 × 3]>
## # ℹ 20 more rowsper_patient_coef <- per_patient_models %>%
mutate(model = map(data, ~ lm(DV_ln ~ TIME, data = .)),
intercept = map_dbl(model, ~ coef(.x)[1]),
time0 = exp(intercept)) %>%
select(ID, TH, time0)
datatable(per_patient_coef)Exercise 10
The concentration at time zero, \(time0 \approx \frac{dose}{distribution volume}\). Assume that the distribution volume in the population is normally distributed. Can you determine whether there is a statistically significant difference in the distribution volume of phenobarbital for patients with TH treatment and without TH treatment?
Solution
The t-test shows that the distribution volumes are not significantly different between the TH and non-TH patients (p > 0.05). In conclusion, TH does not appear to affect distribution of phenobarbital.
# write your code here
dose <- 40
per_patient_vd <- per_patient_coef %>%
mutate(Vd = dose / time0)
t_test <- t.test(Vd ~ TH, data = per_patient_vd)
t_test##
## Welch Two Sample t-test
##
## data: Vd by TH
## t = -1.4623, df = 25.683, p-value = 0.1558
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
## -0.16983603 0.02868954
## sample estimates:
## mean in group 0 mean in group 1
## 0.9136266 0.9841998Conclusion
This study compared the pharmacokinetics of phenobarbital in newborns treated with therapeutic hypothermia (TH) versus normothermic controls. Peak concentrations (Cmax) were similar between groups, but TH-treated patients showed a slightly slower elimination, suggesting reduced clearance. Time-zero concentrations and derived distribution volumes were comparable, and statistical testing did not show significant differences in distribution volume between groups. Overall, therapeutic hypothermia appears to modestly affect phenobarbital elimination without substantially altering distribution volume, highlighting the need for careful monitoring in these patients.