Lymph node position variation for rectal cancer patients treated with long course radiotherapy: Supplementary data analysis
Dennis Tideman Arp
Ane L. Appelt
Rasmus Froberg Brøndum
Rasa Mikalone
Martin Skovmos Nielsen
Laurids Østergaard Poulsen
29/05/2024
Read and transform data
data <- read_tsv("../ExternalData/Data_lymphnodes_141223.txt") %>%
mutate(Patient = factor(sub("Pt", "#", Patient), levels = paste0("#", 1:14)),
patient_node = factor(sub("Pt", "#", sub(" MR[0-9]", "", Unq_ID)),
levels = paste0(rep(paste0("#", 1:14), each = 8),"_", 1:8)),
`Pretreat/during` = factor(`Pretreat/during`, levels = c("Pretreat", "During")),
Rel_bladder = Rel_bladder / 100,
Baseline_bladder = Baseline_bladder / 100,
Shift_LR = Shift_X * 10,
Shift_AP = Shift_Y * 10,
Shift_CC = Shift_Z * 10,
Shift_LR_GTV = Shift_X_GTV * 10,
Shift_AP_GTV = Shift_Y_GTV * 10,
Shift_CC_GTV = Shift_Z_GTV * 10,
Location_plot = case_when(Location == "Pre" ~ "Presacral",
Location == "LN" ~ "Lateral nodes",
Location == "Meso" ~ "Mesorectum")) Compare the shifs for bony and soft tissue match
bind_rows(data.frame(Location = data$Location_plot, Shift = data$Shift_LR, Shift_GTV = data$Shift_LR_GTV, Direction = "LR"),
data.frame(Location = data$Location_plot, Shift = data$Shift_AP, Shift_GTV = data$Shift_AP_GTV, Direction = "AP"),
data.frame(Location = data$Location_plot, Shift = data$Shift_CC, Shift_GTV = data$Shift_CC_GTV, Direction = "CC")) %>%
ggplot(aes(x = Shift, y = Shift_GTV)) +
geom_point() +
xlab("Bony match") +
ylab("Soft tissue match") +
facet_grid(cols = vars(Location), rows = vars(Direction))
Position variation relative to bony structures
Visualize data
First we try to plot the data for each patient and Shift direction
fig1 <- data %>%
dplyr::select(Patient, Location_plot, Shift_LR:Shift_CC) %>%
pivot_longer(cols = Shift_LR:Shift_CC) %>%
mutate(name = factor(sub("Shift_","", name), levels = c("LR", "AP", "CC"))) %>%
ggplot(aes(x = as.factor(Patient), y = value, fill = as.factor(Patient))) +
geom_boxplot() +
xlab("Patient") +
ylab("Shift (mm)") +
scale_fill_discrete(name = "Patient") +
#theme_bw() +
theme_classic() +
facet_grid(rows = vars(name), cols = vars(Location_plot), space = "free_x", scales = "free_x") +
theme(legend.position = "none", axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
strip.text.y.right = element_text(angle = 0)) +
scale_fill_manual(values = BloodCol)
fig1
ggsave(plot = fig1, filename = "../Output/fig1.pdf")Then we plot the data for individual nodes
fig2 <- data %>% dplyr::select(Patient, patient_node, Location_plot, Shift_LR:Shift_CC) %>%
pivot_longer(cols = Shift_LR:Shift_CC) %>%
mutate(name = factor(sub("Shift_","", name), levels = c("LR", "AP", "CC"))) %>%
ggplot(aes(x = as.factor(patient_node), y = value, fill = as.factor(Patient))) +
geom_boxplot() +
xlab("Patient / Node") +
ylab("Shift (mm)") +
scale_fill_discrete(name = "Patient") +
theme_classic() +
#theme_bw() +
facet_grid(rows = vars(name), cols = vars(Location_plot), space = "free_x", scales = "free_x") +
theme(legend.position = "none", axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
strip.text.y.right = element_text(angle = 0)) +
scale_fill_manual(values = BloodCol)
fig2
ggsave(plot = fig2, filename = "../Output/fig2.pdf")Systematic and random var
A lot of the variance seems to be between nodes within patients, i.e. modelling the variance with only an overall and patient offsets, might miss some systematic variance, that would be taken into account by including a node effect. We do this by utilizing a nested mixed model, with individual nodes nested within patients. This is done for all locations and directions in order to estimate the overall mean, systematic, and random variance. The systematic variance is calculated as the sum of variances for the patient and node effects.
sys_rand_results <- list()
for(location in c("Combined", "LN", "Meso", "Pre")){
for(direction in c("LR","AP","CC")){
if(location == "Combined"){
sys_rand_results[[paste0(location,"_",direction)]] <-
lmer(as.formula(paste0("Shift_", direction, " ~ (1|Patient/patient_node)")), data)
} else{
sys_rand_results[[paste0(location,"_",direction)]] <-
lmer(as.formula(paste0("Shift_", direction, " ~ (1|Patient/patient_node)")), data %>% filter(Location == location))
}
}
}summary_results <- data.frame("location" = NULL, "direction" = NULL, "M"=NULL, Patient_Node_SD = NULL, Patient_SD = NULL, "Sys_SD"=NULL, "Rand_SD" = NULL)
for(name in names(sys_rand_results)){
location <- str_split(name, "_")[[1]][1]
direction <- str_split(name, "_")[[1]][2]
#cat(location, direction, "\n")
results <- broom.mixed::tidy(sys_rand_results[[name]])
M <- results$estimate[results$term == "(Intercept)"]
Patient_Node_SD <- results$estimate[results$group == "patient_node:Patient" & results$term == "sd__(Intercept)"]
Patient_SD <- results$estimate[results$group == "Patient" & results$term == "sd__(Intercept)"]
Sys_SD <- sqrt(sum(results$estimate[results$term == "sd__(Intercept)"])^2)
Rand_SD <- results$estimate[results$term == "sd__Observation"]
summary_results <- rbind(summary_results, data.frame(location, direction, M, Patient_Node_SD, Patient_SD, Sys_SD, Rand_SD))
}summary_results %>%
pivot_longer(cols = M:Rand_SD) %>%
pivot_wider(names_from = direction) %>%
kable()| location | name | LR | AP | CC |
|---|---|---|---|---|
| Combined | M | 0.2010718 | -0.6743243 | -0.3872820 |
| Combined | Patient_Node_SD | 1.1650586 | 1.4291795 | 1.0566977 |
| Combined | Patient_SD | 0.7981935 | 0.9786845 | 1.8069732 |
| Combined | Sys_SD | 1.9632521 | 2.4078641 | 2.8636709 |
| Combined | Rand_SD | 1.2373878 | 1.7499741 | 1.8351158 |
| LN | M | 0.4556550 | -0.8218578 | -0.0777707 |
| LN | Patient_Node_SD | 0.3765537 | 0.0000000 | 0.1482093 |
| LN | Patient_SD | 1.8799117 | 0.5911154 | 0.4147415 |
| LN | Sys_SD | 2.2564654 | 0.5911154 | 0.5629508 |
| LN | Rand_SD | 0.7944124 | 0.7028148 | 0.6418723 |
| Meso | M | 0.2675469 | -1.0468786 | -0.0084708 |
| Meso | Patient_Node_SD | 1.2938815 | 1.3142692 | 0.9713403 |
| Meso | Patient_SD | 0.8592124 | 1.4895386 | 1.2306529 |
| Meso | Sys_SD | 2.1530939 | 2.8038078 | 2.2019931 |
| Meso | Rand_SD | 1.3001927 | 1.8160538 | 2.0644052 |
| Pre | M | -0.1300072 | 0.1150739 | -1.5381713 |
| Pre | Patient_Node_SD | 0.4472965 | 0.0000000 | 0.0000000 |
| Pre | Patient_SD | 0.5663485 | 1.9273261 | 3.6762997 |
| Pre | Sys_SD | 1.0136450 | 1.9273261 | 3.6762997 |
| Pre | Rand_SD | 1.3660713 | 2.1588453 | 1.7674469 |
Model control
sys_rand_results %>%
purrr::map(~ {.x %>% residuals}) %>%
bind_rows(.id = "id") %>%
column_to_rownames("id") %>%
t() %>%
data.frame() %>%
pivot_longer(everything()) %>%
mutate(location = sapply(strsplit(name,"_"), function(x) x[[1]]),
direction = sapply(strsplit(name,"_"), function(x) x[[2]])) %>%
ggplot(aes(sample = value)) +
stat_qq() +
stat_qq_line() +
facet_grid(cols = vars(direction), rows = vars(location))
Effect of location
We test if the location of the nodes have an effect on the magnitude of the shift in either direction
loc_results <- list()
for(direction in c("LR","AP","CC")){
loc_results[[direction]] <-
lmer(as.formula(paste0("Shift_", direction, " ~ Location + (1|Patient/patient_node)")), data)
}We extract the estimates and p.values of the location effects. The model design uses contrasts, i.e. “LN” is included in the intercept, and effects for “Meso” and “Pre” is given as the difference in mean compared to “LN”. Results indicate no significant difference in the average magnitude of the shift across the different locations for all directions.
lapply(loc_results, broom.mixed::tidy) %>%
bind_rows(.id = "Direction") %>%
filter(effect == "fixed") %>%
dplyr::select(Direction, estimate, term, p.value) %>%
pivot_wider(names_from = "Direction", values_from = c("estimate", "p.value")) %>%
dplyr::select(term, estimate_LR, p.value_LR, estimate_AP, p.value_AP, estimate_CC, p.value_CC) %>%
kable()| term | estimate_LR | p.value_LR | estimate_AP | p.value_AP | estimate_CC | p.value_CC |
|---|---|---|---|---|---|---|
| (Intercept) | 0.2645534 | 0.5943285 | -0.1741210 | 0.7773523 | -0.8446492 | 0.2327226 |
| LocationMeso | -0.0085068 | 0.9864788 | -0.8031499 | 0.1914590 | 0.5896335 | 0.2789507 |
| LocationPre | -0.3256012 | 0.5974848 | 0.2386066 | 0.7498443 | 0.3815595 | 0.5661269 |
Model control
loc_results %>%
purrr::map(~ {.x %>% residuals}) %>%
bind_rows(.id = "id") %>%
column_to_rownames("id") %>%
t() %>%
data.frame() %>%
pivot_longer(everything()) %>%
ggplot(aes(sample = value)) +
stat_qq() +
stat_qq_line() +
facet_wrap(~name, ncol = 3)
Effect of Pretreat/During
We test if the radiation time (Pretreat / During) of the nodes have an effect on the magnitude of the shift in either direction
predur_results <- list()
for(direction in c("LR","AP","CC")){
predur_results[[direction]] <-
lmer(as.formula(paste0("Shift_", direction, " ~ `Pretreat/during` + (1|Patient/patient_node)")), data)
}We extract the estimates and p.values of the effects. The model design uses contrasts, i.e. “Pretreat” is included in the intercept, and the effect of “During” is given as the difference in mean compared to “Pretreat”. Results indicate a significant difference in the magnitude of the shift in the AP direction.
lapply(predur_results, broom.mixed::tidy) %>%
bind_rows(.id = "Direction") %>%
filter(effect == "fixed") %>%
dplyr::select(Direction, estimate, term, p.value) %>%
pivot_wider(names_from = "Direction", values_from = c("estimate", "p.value")) %>%
dplyr::select(term, estimate_LR, p.value_LR, estimate_AP, p.value_AP, estimate_CC, p.value_CC) %>%
kable()| term | estimate_LR | p.value_LR | estimate_AP | p.value_AP | estimate_CC | p.value_CC |
|---|---|---|---|---|---|---|
| (Intercept) | 0.1266453 | 0.6737760 | -0.1479231 | 0.6938437 | -0.4067930 | 0.4653648 |
Pretreat/duringDuring |
0.1240437 | 0.3921726 | -0.8773224 | 0.0000130 | 0.0325137 | 0.8798785 |
Model control
predur_results %>%
purrr::map(~ {.x %>% residuals}) %>%
bind_rows(.id = "id") %>%
column_to_rownames("id") %>%
t() %>%
data.frame() %>%
pivot_longer(everything()) %>%
ggplot(aes(sample = value)) +
stat_qq() +
stat_qq_line() +
facet_wrap(~name, ncol = 3)
Effect of initial bladder volume
We test if the baseline bladder volume has an effect on the magnitude of the shift in either direction
baseline_bladder_results <- list()
for(direction in c("LR","AP","CC")){
baseline_bladder_results[[direction]] <-
lmer(as.formula(paste0("Shift_", direction, " ~ Baseline_bladder + (1|Patient/patient_node)")), data)
}We extract the estimates and p.values of the effects. We see no significant effect in either direction.
lapply(baseline_bladder_results, broom.mixed::tidy) %>%
bind_rows(.id = "Direction") %>%
filter(effect == "fixed") %>%
dplyr::select(Direction, estimate, term, p.value) %>%
pivot_wider(names_from = "Direction", values_from = c("estimate", "p.value")) %>%
dplyr::select(term, estimate_LR, p.value_LR, estimate_AP, p.value_AP, estimate_CC, p.value_CC) %>%
kable()| term | estimate_LR | p.value_LR | estimate_AP | p.value_AP | estimate_CC | p.value_CC |
|---|---|---|---|---|---|---|
| (Intercept) | -0.3992842 | 0.5117812 | -1.3139761 | 0.1071471 | -0.3653089 | 0.7683312 |
| Baseline_bladder | 0.3496578 | 0.2697830 | 0.3712344 | 0.3525229 | -0.0142343 | 0.9817177 |
Model control
baseline_bladder_results %>%
purrr::map(~ {.x %>% residuals}) %>%
bind_rows(.id = "id") %>%
column_to_rownames("id") %>%
t() %>%
data.frame() %>%
pivot_longer(everything()) %>%
ggplot(aes(sample = value)) +
stat_qq() +
stat_qq_line() +
facet_wrap(~name, ncol = 3)
Effect of relative bladder
Here we investigate the effect of the relative volume of the bladder on the magnitude of the shift in either direction and location. First we do a plot of Rel_bladder vs shift, and add simple regression lines for each patient/node combination.
data %>% dplyr::select(Patient, patient_node, Location, Rel_bladder, Shift_LR:Shift_CC) %>%
pivot_longer(cols = Shift_LR:Shift_CC) %>%
mutate(name = factor(sub("Shift_","", name), levels = c("LR", "AP", "CC"))) %>%
ggplot(aes(x = Rel_bladder, y = value, col = patient_node)) +
geom_point() +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
facet_grid(rows = vars(Location), cols = vars(name)) +
theme(legend.position = "none")
We test the effect of the relative fill of the bladder combined and stratified by location and direction using the same nested mixed model as before, i.e. each patient/node gets an indvidual intercept.
#bladder_strat_results <- list()
#for(location in c("LN", "Meso", "Pre")){
# for(direction in c("LR","AP","CC")){
# bladder_strat_results[[paste0(location,"_",direction)]] <-
# lmer(as.formula(paste0("Shift_", direction, " ~ Rel_bladder + (1|Patient/patient_node)")), data %>% filter(Location == location))
# }
#}
bladder_strat_results <- list()
for(location in c("Combined", "LN", "Meso", "Pre")){
for(direction in c("LR","AP","CC")){
if(location == "Combined"){
bladder_strat_results[[paste0(location,"_",direction)]] <-
lmer(as.formula(paste0("Shift_", direction, " ~ Rel_bladder + (1|Patient/patient_node)")), data)
} else{
bladder_strat_results[[paste0(location,"_",direction)]] <-
lmer(as.formula(paste0("Shift_", direction, " ~ Rel_bladder + (1|Patient/patient_node)")), data %>% filter(Location == location))
}
}
}Results show a tendency for the relative bladder volume to increase the magnitude of the shift in the LR direction for LN, and decrease the magnitude of the shift in the AP direction for Pre. For the combined model we see no significance. Effect sizes are given in mm/(100 ml).
lapply(bladder_strat_results, broom.mixed::tidy) %>%
bind_rows(.id = "test") %>%
filter(term == "Rel_bladder") %>%
mutate(location = sapply(strsplit(test,"_"), function(x) x[[1]]),
direction = sapply(strsplit(test,"_"), function(x) x[[2]])) %>%
dplyr::select(location, direction, estimate, p.value) %>%
pivot_wider(names_from = "direction", values_from = c("estimate", "p.value")) %>%
dplyr::select(location, estimate_LR, p.value_LR, estimate_AP, p.value_AP, estimate_CC, p.value_CC) %>%
kable()| location | estimate_LR | p.value_LR | estimate_AP | p.value_AP | estimate_CC | p.value_CC |
|---|---|---|---|---|---|---|
| Combined | 0.1210409 | 0.4384574 | -0.2126829 | 0.3163010 | -0.0600586 | 0.8035887 |
| LN | 1.0306305 | 0.0000757 | -0.1629788 | 0.4385251 | -0.0433228 | 0.8072436 |
| Meso | -0.0145070 | 0.9386773 | -0.0864561 | 0.7448413 | -0.1855976 | 0.5018719 |
| Pre | 0.2362193 | 0.6491276 | -2.0525819 | 0.0160118 | 1.2142631 | 0.1067566 |
Model control
bladder_strat_results %>%
purrr::map(~ {.x %>% residuals}) %>%
bind_rows(.id = "id") %>%
column_to_rownames("id") %>%
t() %>%
data.frame() %>%
pivot_longer(everything()) %>%
mutate(location = sapply(strsplit(name,"_"), function(x) x[[1]]),
direction = sapply(strsplit(name,"_"), function(x) x[[2]])) %>%
ggplot(aes(sample = value)) +
stat_qq() +
stat_qq_line() +
facet_grid(cols = vars(direction), rows = vars(location))
Multivariable model
## model for bladder vs shift with interaction on location. Seems to disrupt the effect on shift_x in pre, seen in the stratified model.
bladder_comb_results <- list()
for(direction in c("LR","AP","CC")){
bladder_comb_results[[paste0(direction)]] <-
lmer(as.formula(paste0("Shift_", direction, " ~ Rel_bladder*Location + Baseline_bladder + `Pretreat/during` + (1|Patient/patient_node)")), data)
}
multires <- lapply(bladder_comb_results, broom.mixed::tidy)
kable(multires$LR, caption = "Results of multivariate model for LR")| effect | group | term | estimate | std.error | statistic | df | p.value |
|---|---|---|---|---|---|---|---|
| fixed | NA | (Intercept) | -0.4483684 | 0.8267699 | -0.5423134 | 25.21826 | 0.5923632 |
| fixed | NA | Rel_bladder | 0.5403734 | 0.3721851 | 1.4518942 | 272.99282 | 0.1476793 |
| fixed | NA | LocationMeso | -0.2665242 | 0.6190836 | -0.4305142 | 80.93659 | 0.6679659 |
| fixed | NA | LocationPre | -0.4169470 | 0.7475431 | -0.5577565 | 86.44987 | 0.5784522 |
| fixed | NA | Baseline_bladder | 0.5902214 | 0.3386970 | 1.7426239 | 16.55700 | 0.0999399 |
| fixed | NA | Pretreat/duringDuring |
0.1602233 | 0.1465819 | 1.0930633 | 239.20257 | 0.2754656 |
| fixed | NA | Rel_bladder:LocationMeso | -0.3383002 | 0.3902254 | -0.8669354 | 234.24621 | 0.3868645 |
| fixed | NA | Rel_bladder:LocationPre | -0.0975009 | 0.6150846 | -0.1585162 | 296.27146 | 0.8741580 |
| ran_pars | patient_node:Patient | sd__(Intercept) | 1.2238461 | NA | NA | NA | NA |
| ran_pars | Patient | sd__(Intercept) | 0.7225389 | NA | NA | NA | NA |
| ran_pars | Residual | sd__Observation | 1.2332279 | NA | NA | NA | NA |
kable(multires$AP, caption = "Results of multivariate model for AP")| effect | group | term | estimate | std.error | statistic | df | p.value |
|---|---|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 0.3782706 | 1.0819588 | 0.3496165 | 26.50545 | 0.7293892 |
| fixed | NA | Rel_bladder | 0.2279167 | 0.4817259 | 0.4731253 | 248.20310 | 0.6365396 |
| fixed | NA | LocationMeso | -1.2175390 | 0.7473219 | -1.6292029 | 85.03680 | 0.1069692 |
| fixed | NA | LocationPre | -0.9596202 | 0.8937768 | -1.0736688 | 87.64266 | 0.2859184 |
| fixed | NA | Baseline_bladder | 0.1445068 | 0.4643286 | 0.3112166 | 18.21203 | 0.7591702 |
| fixed | NA | Pretreat/duringDuring |
-0.9193817 | 0.2001493 | -4.5934798 | 240.11055 | 0.0000070 |
| fixed | NA | Rel_bladder:LocationMeso | -0.4942605 | 0.4966695 | -0.9951497 | 191.45429 | 0.3209195 |
| fixed | NA | Rel_bladder:LocationPre | -1.9714140 | 0.8088849 | -2.4371998 | 287.48609 | 0.0154080 |
| ran_pars | patient_node:Patient | sd__(Intercept) | 1.3021503 | NA | NA | NA | NA |
| ran_pars | Patient | sd__(Intercept) | 1.1231176 | NA | NA | NA | NA |
| ran_pars | Residual | sd__Observation | 1.6839764 | NA | NA | NA | NA |
kable(multires$CC, caption = "Results of multivariate model for CC")| effect | group | term | estimate | std.error | statistic | df | p.value |
|---|---|---|---|---|---|---|---|
| fixed | NA | (Intercept) | -0.3909911 | 1.3644571 | -0.2865544 | 17.12426 | 0.7778916 |
| fixed | NA | Rel_bladder | -0.0805933 | 0.4862616 | -0.1657406 | 205.42844 | 0.8685241 |
| fixed | NA | LocationMeso | 0.2747550 | 0.6956535 | 0.3949596 | 74.62727 | 0.6939993 |
| fixed | NA | LocationPre | 1.3354288 | 0.8212550 | 1.6260831 | 79.68728 | 0.1078805 |
| fixed | NA | Baseline_bladder | -0.2008768 | 0.6373896 | -0.3151554 | 13.54160 | 0.7574477 |
| fixed | NA | Pretreat/duringDuring |
-0.0020683 | 0.2192924 | -0.0094315 | 240.01235 | 0.9924827 |
| fixed | NA | Rel_bladder:LocationMeso | -0.2737144 | 0.4851239 | -0.5642154 | 128.79906 | 0.5735887 |
| fixed | NA | Rel_bladder:LocationPre | 1.9314657 | 0.8320990 | 2.3211970 | 248.98927 | 0.0210843 |
| ran_pars | patient_node:Patient | sd__(Intercept) | 0.9528923 | NA | NA | NA | NA |
| ran_pars | Patient | sd__(Intercept) | 1.8512367 | NA | NA | NA | NA |
| ran_pars | Residual | sd__Observation | 1.8441888 | NA | NA | NA | NA |
Position variation relative to primary tumor
Visualize data
First we try to plot the data for each patient and Shift direction
data %>% dplyr::select(Patient, Location_plot, Shift_LR_GTV:Shift_CC_GTV) %>%
pivot_longer(cols = Shift_LR_GTV:Shift_CC_GTV) %>%
mutate(name = factor(sapply(str_split(name, "_"), function(x) x[2]), levels = c("LR", "AP", "CC"))) %>%
ggplot(aes(x = as.factor(Patient), y = value, fill = as.factor(Patient))) +
geom_boxplot() +
xlab("Patient") +
ylab("Shift (mm)") +
scale_fill_discrete(name = "Patient") +
#theme_bw() +
theme_classic() +
facet_grid(rows = vars(name), cols = vars(Location_plot), space = "free_x", scales = "free_x") +
theme(legend.position = "none", axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
strip.text.y.right = element_text(angle = 0)) +
scale_fill_manual(values = BloodCol) 
Then we plot the data for individual nodes
data %>% dplyr::select(Patient, patient_node, Location_plot, Shift_LR_GTV:Shift_CC_GTV) %>%
pivot_longer(cols = Shift_LR_GTV:Shift_CC_GTV) %>%
mutate(name = factor(sapply(str_split(name, "_"), function(x) x[2]), levels = c("LR", "AP", "CC"))) %>%
ggplot(aes(x = as.factor(patient_node), y = value, fill = as.factor(Patient))) +
geom_boxplot() +
xlab("Patient / Node") +
ylab("Shift (mm)") +
scale_fill_discrete(name = "Patient") +
#theme_bw() +
theme_classic() +
facet_grid(rows = vars(name), cols = vars(Location_plot), space = "free_x", scales = "free_x") +
theme(legend.position = "none", axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
strip.text.y.right = element_text(angle = 0)) +
scale_fill_manual(values = BloodCol)
Systematic and random var
We use a similar model for estimating the systematic and random variance
sys_rand_results <- list()
for(location in c("Combined", "LN", "Meso", "Pre")){
for(direction in c("LR","AP","CC")){
if(location == "Combined"){
sys_rand_results[[paste0(location,"_",direction)]] <-
lmer(as.formula(paste0("Shift_", direction, "_GTV ~ (1|Patient/patient_node)")), data)
} else{
sys_rand_results[[paste0(location,"_",direction)]] <-
lmer(as.formula(paste0("Shift_", direction, "_GTV ~ (1|Patient/patient_node)")), data %>% filter(Location == location))
}
}
}summary_results <- data.frame("location" = NULL, "direction" = NULL, "M"=NULL, Patient_Node_SD = NULL, Patient_SD = NULL, "Sys_SD"=NULL, "Rand_SD" = NULL)
for(name in names(sys_rand_results)){
location <- str_split(name, "_")[[1]][1]
direction <- str_split(name, "_")[[1]][2]
#cat(location, direction, "\n")
results <- broom.mixed::tidy(sys_rand_results[[name]])
M <- results$estimate[results$term == "(Intercept)"]
Patient_Node_SD <- results$estimate[results$group == "patient_node:Patient" & results$term == "sd__(Intercept)"]
Patient_SD <- results$estimate[results$group == "Patient" & results$term == "sd__(Intercept)"]
Sys_SD <- sqrt(sum(results$estimate[results$term == "sd__(Intercept)"])^2)
Rand_SD <- results$estimate[results$term == "sd__Observation"]
summary_results <- rbind(summary_results, data.frame(location, direction, M, Patient_Node_SD, Patient_SD, Sys_SD, Rand_SD))
}summary_results %>%
pivot_longer(cols = M:Rand_SD) %>%
pivot_wider(names_from = direction) %>%
kable()| location | name | LR | AP | CC |
|---|---|---|---|---|
| Combined | M | -0.2717650 | -0.1162171 | -0.0361357 |
| Combined | Patient_Node_SD | 1.0775641 | 1.3152756 | 0.8151997 |
| Combined | Patient_SD | 1.4395725 | 1.9131843 | 2.2620974 |
| Combined | Sys_SD | 2.5171366 | 3.2284599 | 3.0772971 |
| Combined | Rand_SD | 1.6623952 | 2.1772837 | 2.3860259 |
| LN | M | -0.4199141 | -0.0084884 | 0.6535926 |
| LN | Patient_Node_SD | 0.0000000 | 0.0000000 | 0.0000000 |
| LN | Patient_SD | 1.8853376 | 1.7554829 | 2.0785069 |
| LN | Sys_SD | 1.8853376 | 1.7554829 | 2.0785069 |
| LN | Rand_SD | 1.3275719 | 1.8435798 | 2.3356596 |
| Meso | M | -0.1853583 | -0.5154813 | 0.4862730 |
| Meso | Patient_Node_SD | 1.2058598 | 1.2535963 | 0.7943324 |
| Meso | Patient_SD | 1.5743252 | 2.2235336 | 1.7825555 |
| Meso | Sys_SD | 2.7801851 | 3.4771299 | 2.5768879 |
| Meso | Rand_SD | 1.8013233 | 2.0756091 | 2.4121591 |
| Pre | M | -0.4360391 | -0.3053377 | -1.1587052 |
| Pre | Patient_Node_SD | 0.2707124 | 0.0000000 | 0.0000000 |
| Pre | Patient_SD | 0.6023547 | 2.6858341 | 3.1918241 |
| Pre | Sys_SD | 0.8730670 | 2.6858341 | 3.1918241 |
| Pre | Rand_SD | 1.4189305 | 2.6386131 | 2.1843062 |
Model control
sys_rand_results %>%
purrr::map(~ {.x %>% residuals}) %>%
bind_rows(.id = "id") %>%
column_to_rownames("id") %>%
t() %>%
data.frame() %>%
pivot_longer(everything()) %>%
mutate(location = sapply(strsplit(name,"_"), function(x) x[[1]]),
direction = sapply(strsplit(name,"_"), function(x) x[[2]])) %>%
ggplot(aes(sample = value)) +
stat_qq() +
stat_qq_line() +
facet_grid(cols = vars(direction), rows = vars(location))
Effect of location
We test if the location of the nodes have an effect on the magnitude of the shift in either direction
loc_results <- list()
for(direction in c("LR","AP","CC")){
loc_results[[direction]] <-
lmer(as.formula(paste0("Shift_", direction, "_GTV ~ Location + (1|Patient/patient_node)")), data)
}We extract the estimates and p.values of the location effects. The model design uses contrasts, i.e. “LN” is included in the intercept, and effects for “Meso” and “Pre” is given as the difference in mean compared to “LN”. Results indicate no significant difference in the magnitude of the shift across the different locations for all directions.
lapply(loc_results, broom.mixed::tidy) %>%
bind_rows(.id = "Direction") %>%
filter(effect == "fixed") %>%
dplyr::select(Direction, estimate, term, p.value) %>%
pivot_wider(names_from = "Direction", values_from = c("estimate", "p.value")) %>%
dplyr::select(term, estimate_LR, p.value_LR, estimate_AP, p.value_AP, estimate_CC, p.value_CC) %>%
kable()| term | estimate_LR | p.value_LR | estimate_AP | p.value_AP | estimate_CC | p.value_CC |
|---|---|---|---|---|---|---|
| (Intercept) | -0.3422084 | 0.5804914 | 0.5798518 | 0.4613936 | -0.5184555 | 0.5160629 |
| LocationMeso | 0.1407248 | 0.7910635 | -1.0232053 | 0.1136012 | 0.6225594 | 0.2584632 |
| LocationPre | -0.1141043 | 0.8608102 | -0.0675561 | 0.9312341 | 0.3983822 | 0.5539322 |
Model control
loc_results %>%
purrr::map(~ {.x %>% residuals}) %>%
bind_rows(.id = "id") %>%
column_to_rownames("id") %>%
t() %>%
data.frame() %>%
pivot_longer(everything()) %>%
ggplot(aes(sample = value)) +
stat_qq() +
stat_qq_line() +
facet_wrap(~name, ncol = 3)
Effect of Pretreat/During
We test if the radiation time (Pretreat / During) of the nodes have an effect on the magnitude of the shift in either direction
predur_results <- list()
for(direction in c("LR","AP","CC")){
predur_results[[direction]] <-
lmer(as.formula(paste0("Shift_", direction, "_GTV ~ `Pretreat/during` + (1|Patient/patient_node)")), data)
}We extract the estimates and p.values of the effects. The model design uses contrasts, i.e. “Pretreat” is included in the intercept, and the effect of “During” is given as the difference in mean compared to “Pretreat”. Results indicate a significant difference in the magnitude of the shift in the AP and CC directions.
lapply(predur_results, broom.mixed::tidy) %>%
bind_rows(.id = "Direction") %>%
filter(effect == "fixed") %>%
dplyr::select(Direction, estimate, term, p.value) %>%
pivot_wider(names_from = "Direction", values_from = c("estimate", "p.value")) %>%
dplyr::select(term, estimate_LR, p.value_LR, estimate_AP, p.value_AP, estimate_CC, p.value_CC) %>%
kable()| term | estimate_LR | p.value_LR | estimate_AP | p.value_AP | estimate_CC | p.value_CC |
|---|---|---|---|---|---|---|
| (Intercept) | -0.2125844 | 0.6431913 | 0.5296932 | 0.3816160 | -0.5323658 | 0.4353069 |
Pretreat/duringDuring |
-0.0986339 | 0.6127151 | -1.0765027 | 0.0000173 | 0.8270492 | 0.0028483 |
Model control
predur_results %>%
purrr::map(~ {.x %>% residuals}) %>%
bind_rows(.id = "id") %>%
column_to_rownames("id") %>%
t() %>%
data.frame() %>%
pivot_longer(everything()) %>%
ggplot(aes(sample = value)) +
stat_qq() +
stat_qq_line() +
facet_wrap(~name, ncol = 3)
Effect of initial bladder volume
We test if the baseline bladder volume has an effect on the magnitude of the shift in either direction
baseline_bladder_results <- list()
for(direction in c("LR","AP","CC")){
baseline_bladder_results[[direction]] <-
lmer(as.formula(paste0("Shift_", direction, "_GTV ~ Baseline_bladder + (1|Patient/patient_node)")), data)
}We extract the estimates and p.values of the effects. We see no significant effect in either direction.
lapply(baseline_bladder_results, broom.mixed::tidy) %>%
bind_rows(.id = "Direction") %>%
filter(effect == "fixed") %>%
dplyr::select(Direction, estimate, term, p.value) %>%
pivot_wider(names_from = "Direction", values_from = c("estimate", "p.value")) %>%
dplyr::select(term, estimate_LR, p.value_LR, estimate_AP, p.value_AP, estimate_CC, p.value_CC) %>%
kable()| term | estimate_LR | p.value_LR | estimate_AP | p.value_AP | estimate_CC | p.value_CC |
|---|---|---|---|---|---|---|
| (Intercept) | 0.1177918 | 0.9065926 | -2.148502 | 0.0798561 | -0.4448488 | 0.7684851 |
| Baseline_bladder | -0.2225138 | 0.6611792 | 1.135978 | 0.0671120 | 0.2291203 | 0.7622322 |
Model control
baseline_bladder_results %>%
purrr::map(~ {.x %>% residuals}) %>%
bind_rows(.id = "id") %>%
column_to_rownames("id") %>%
t() %>%
data.frame() %>%
pivot_longer(everything()) %>%
ggplot(aes(sample = value)) +
stat_qq() +
stat_qq_line() +
facet_wrap(~name, ncol = 3)
Effect of relative bladder
Here we investigate the effect of the relative volume of the bladder on the magnitude of the shift in either direction and location. First we do a plot of Rel_bladder vs shift, and add simple regression lines for each patient/node combination.
data %>% dplyr::select(Patient, patient_node, Location, Rel_bladder, Shift_LR_GTV:Shift_CC_GTV) %>%
pivot_longer(cols = Shift_LR_GTV:Shift_CC_GTV) %>%
mutate(name = factor(sapply(str_split(name, "_"), function(x) x[2]), levels = c("LR", "AP", "CC"))) %>%
ggplot(aes(x = Rel_bladder, y = value, col = patient_node)) +
geom_point() +
stat_smooth(method = "lm",
formula = y ~ x,
geom = "smooth") +
facet_grid(rows = vars(Location), cols = vars(name)) +
theme(legend.position = "none")
We test the effect of the relative fill of the bladder combined and stratified by location and direction using the same nested mixed model as before, i.e. each patient/node gets an indvidual intercept.
#bladder_strat_results <- list()
#for(location in c("LN", "Meso", "Pre")){
# for(direction in c("LR","AP","CC")){
# bladder_strat_results[[paste0(location,"_",direction)]] <-
# lmer(as.formula(paste0("Shift_", direction, "_GTV ~ Rel_bladder + (1|Patient/patient_node)")), data %>% filter(Location == location))
# }
#}
bladder_strat_results <- list()
for(location in c("Combined", "LN", "Meso", "Pre")){
for(direction in c("LR","AP","CC")){
if(location == "Combined"){
bladder_strat_results[[paste0(location,"_",direction)]] <-
lmer(as.formula(paste0("Shift_", direction, "_GTV ~ Rel_bladder + (1|Patient/patient_node)")), data)
} else{
bladder_strat_results[[paste0(location,"_",direction)]] <-
lmer(as.formula(paste0("Shift_", direction, "_GTV ~ Rel_bladder + (1|Patient/patient_node)")), data %>% filter(Location == location))
}
}
}Results show a tendency for the relative bladder volume to increase the magnitude of the shift in the LR and AP directions for LN, while combining across locations only gives significance for AP. Effect sizes are given as mm/(100 ml)
lapply(bladder_strat_results, broom.mixed::tidy) %>%
bind_rows(.id = "test") %>%
filter(term == "Rel_bladder") %>%
mutate(location = sapply(strsplit(test,"_"), function(x) x[[1]]),
direction = sapply(strsplit(test,"_"), function(x) x[[2]])) %>%
dplyr::select(location, direction, estimate, p.value) %>%
pivot_wider(names_from = "direction", values_from = c("estimate", "p.value")) %>%
dplyr::select(location, estimate_LR, p.value_LR, estimate_AP, p.value_AP, estimate_CC, p.value_CC) %>%
kable()| location | estimate_LR | p.value_LR | estimate_AP | p.value_AP | estimate_CC | p.value_CC |
|---|---|---|---|---|---|---|
| Combined | 0.3769833 | 0.0784180 | -0.6809947 | 0.0138516 | 0.2991526 | 0.3361438 |
| LN | 1.4067589 | 0.0012831 | -1.3309466 | 0.0138977 | 0.7930715 | 0.2468077 |
| Meso | 0.2005602 | 0.4480756 | -0.5287316 | 0.0908914 | -0.0232290 | 0.9448980 |
| Pre | 0.4055876 | 0.4439214 | -1.1283503 | 0.2961191 | 1.4343397 | 0.1161893 |
Model control
bladder_strat_results %>%
purrr::map(~ {.x %>% residuals}) %>%
bind_rows(.id = "id") %>%
column_to_rownames("id") %>%
t() %>%
data.frame() %>%
pivot_longer(everything()) %>%
mutate(location = sapply(strsplit(name,"_"), function(x) x[[1]]),
direction = sapply(strsplit(name,"_"), function(x) x[[2]])) %>%
ggplot(aes(sample = value)) +
stat_qq() +
stat_qq_line() +
facet_grid(cols = vars(direction), rows = vars(location))
Multivariable model
## model for bladder vs shift with interaction on location. Seems to disrupt the effect on shift_x in pre, seen in the stratified model.
bladder_comb_results <- list()
for(direction in c("LR","AP","CC")){
bladder_comb_results[[paste0(direction)]] <-
lmer(as.formula(paste0("Shift_", direction, "_GTV ~ Rel_bladder*Location + Baseline_bladder + `Pretreat/during` + (1|Patient/patient_node)")), data)
}
multires <- lapply(bladder_comb_results, broom.mixed::tidy)
kable(multires$LR, caption = "Results of multivariate model for LR")| effect | group | term | estimate | std.error | statistic | df | p.value |
|---|---|---|---|---|---|---|---|
| fixed | NA | (Intercept) | -0.0812829 | 1.2175389 | -0.0667600 | 18.82398 | 0.9474768 |
| fixed | NA | Rel_bladder | 0.5571822 | 0.4653913 | 1.1972339 | 230.69708 | 0.2324442 |
| fixed | NA | LocationMeso | -0.1028153 | 0.7092340 | -0.1449667 | 79.30152 | 0.8851050 |
| fixed | NA | LocationPre | -0.2108367 | 0.8372974 | -0.2518062 | 80.24619 | 0.8018346 |
| fixed | NA | Baseline_bladder | 0.1259385 | 0.5493045 | 0.2292691 | 13.55801 | 0.8220822 |
| fixed | NA | Pretreat/duringDuring |
-0.0445372 | 0.1972378 | -0.2258045 | 238.00351 | 0.8215473 |
| fixed | NA | Rel_bladder:LocationMeso | -0.2226693 | 0.4760226 | -0.4677703 | 165.27437 | 0.6405651 |
| fixed | NA | Rel_bladder:LocationPre | 0.0022740 | 0.7843740 | 0.0028991 | 274.76641 | 0.9976890 |
| ran_pars | patient_node:Patient | sd__(Intercept) | 1.1322225 | NA | NA | NA | NA |
| ran_pars | Patient | sd__(Intercept) | 1.5251324 | NA | NA | NA | NA |
| ran_pars | Residual | sd__Observation | 1.6587175 | NA | NA | NA | NA |
kable(multires$AP, caption = "Results of multivariate model for AP")| effect | group | term | estimate | std.error | statistic | df | p.value |
|---|---|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 0.2569344 | 1.4455307 | 0.1777440 | 22.02002 | 0.8605490 |
| fixed | NA | Rel_bladder | -0.4483474 | 0.5654427 | -0.7929139 | 215.87984 | 0.4286985 |
| fixed | NA | LocationMeso | -1.2396934 | 0.8306429 | -1.4924505 | 76.21745 | 0.1397090 |
| fixed | NA | LocationPre | -0.6231960 | 0.9828147 | -0.6340930 | 79.42560 | 0.5278422 |
| fixed | NA | Baseline_bladder | 0.4187503 | 0.6585909 | 0.6358276 | 16.47529 | 0.5336223 |
| fixed | NA | Pretreat/duringDuring |
-1.1921816 | 0.2479840 | -4.8074949 | 238.31113 | 0.0000027 |
| fixed | NA | Rel_bladder:LocationMeso | -0.3361974 | 0.5710881 | -0.5886962 | 143.58663 | 0.5569901 |
| fixed | NA | Rel_bladder:LocationPre | -0.9764154 | 0.9619141 | -1.0150754 | 263.27626 | 0.3110018 |
| ran_pars | patient_node:Patient | sd__(Intercept) | 1.2378443 | NA | NA | NA | NA |
| ran_pars | Patient | sd__(Intercept) | 1.8256947 | NA | NA | NA | NA |
| ran_pars | Residual | sd__Observation | 2.0857877 | NA | NA | NA | NA |
kable(multires$CC, caption = "Results of multivariate model for CC")| effect | group | term | estimate | std.error | statistic | df | p.value |
|---|---|---|---|---|---|---|---|
| fixed | NA | (Intercept) | -1.3962200 | 1.6291271 | -0.8570357 | 15.38721 | 0.4045759 |
| fixed | NA | Rel_bladder | 0.8306066 | 0.5620517 | 1.4778119 | 182.60204 | 0.1411809 |
| fixed | NA | LocationMeso | 0.0035166 | 0.7368717 | 0.0047723 | 65.06397 | 0.9962069 |
| fixed | NA | LocationPre | 1.1271695 | 0.8792061 | 1.2820310 | 74.73726 | 0.2037944 |
| fixed | NA | Baseline_bladder | 0.6455385 | 0.7835178 | 0.8238977 | 13.32497 | 0.4244997 |
| fixed | NA | Pretreat/duringDuring |
0.8873404 | 0.2789002 | 3.1815699 | 240.20608 | 0.0016579 |
| fixed | NA | Rel_bladder:LocationMeso | -0.5764313 | 0.5348672 | -1.0777092 | 92.21149 | 0.2839748 |
| fixed | NA | Rel_bladder:LocationPre | 1.7084852 | 0.9771164 | 1.7484971 | 202.04164 | 0.0818968 |
| ran_pars | patient_node:Patient | sd__(Intercept) | 0.7089587 | NA | NA | NA | NA |
| ran_pars | Patient | sd__(Intercept) | 2.3093311 | NA | NA | NA | NA |
| ran_pars | Residual | sd__Observation | 2.3459189 | NA | NA | NA | NA |