Updated analysis script

Author: etpeterson

tests/IVIMmodels/unit_tests/analyze.r

@@ -1,6 +1,24 @@
+#!/usr/bin/env Rscript
+
+#Run like this:
+#Rscript --vanilla tests/IVIMmodels/unit_tests/analyze.r test_output_priors.csv test_duration_priors.csv
+args = commandArgs(trailingOnly=TRUE)
+output_name = "test_output.csv"
+duration_name = "test_duration.csv"
+runPrediction = FALSE
+if (length(args)>=1) {
+   output_name = args[1]
+}
+if (length(args)>=2) {
+    duration_name = args[2]
+}
+print(output_name)
+print(duration_name)
+
 library(plyr)  # order matters, we need dplyr to come later for grouping
 library(dplyr)
 library(tidyverse)
+library(data.table)
 
 
 plot_ivim <- function(data, fileExtension) {
@@ -15,33 +33,33 @@ plot_ivim <- function(data, fileExtension) {
     ggsave(paste("Dp", fileExtension, sep=""), plot=Dp_plot, width = 50, height = 50, units = "cm")
 }
 
-data <- read.csv("test_output.csv")
+data <- read.csv(output_name)
 data <- data %>% mutate_if(is.character, as.factor)
 plot_ivim(data, ".pdf")
 
 data_restricted <- data[data$Region %in% c("Liver", "spleen", "Right kydney cortex", "right kidney medulla"),]
 plot_ivim(data_restricted, "_limited.pdf")
 
-data_duration <- read.csv("test_duration.csv")
+data_duration <- read.csv(duration_name)
 data_duration <- data_duration %>% mutate_if(is.character, as.factor)
 data_duration$ms <- data_duration$Duration..us./data_duration$Count/1000
 ggplot(data_duration, aes(x=Algorithm, y=ms)) + geom_boxplot() + scale_x_discrete(guide = guide_axis(angle = 90)) + ggtitle("Fit Duration") + ylab("Time (ms)")
 ggsave("durations.pdf", width = 20, height = 20, units = "cm")
 
 
-# some kind of black voodoo magic to number unique groups
-# data_new <- data %>% group_by(Algorithm, Region, SNR) %>% mutate(index = row_number()) %>% ungroup()
-# Then this widens it so we can lm()
-data_wide <- data %>% pivot_wider(names_from=Algorithm, values_from=c(f_fitted, Dp_fitted, D_fitted), id_cols=c(Region, SNR, index, f, D, Dp))
-# linear fit for f
-f_model <- lm(f ~ f_fitted_IAR_LU_biexp + f_fitted_IAR_LU_modified_mix + f_fitted_IAR_LU_modified_topopro + f_fitted_IAR_LU_segmented_2step + f_fitted_IAR_LU_segmented_3step + f_fitted_IAR_LU_subtracted, data=data_wide)
-#f_model <- lm(f ~ f_fitted_ETP_SRI_LinearFitting + f_fitted_IAR_LU_biexp + f_fitted_IAR_LU_modified_mix + f_fitted_IAR_LU_modified_topopro + f_fitted_IAR_LU_segmented_2step + f_fitted_IAR_LU_segmented_3step + f_fitted_IAR_LU_subtracted, data=data_new_wide)
-D_model <- lm(D ~ D_fitted_ETP_SRI_LinearFitting + D_fitted_IAR_LU_biexp + D_fitted_IAR_LU_modified_topopro + D_fitted_IAR_LU_segmented_2step + D_fitted_IAR_LU_segmented_3step + D_fitted_IAR_LU_subtracted, data=data_wide)
-# D_model <- lm(D ~ D_fitted_ETP_SRI_LinearFitting + D_fitted_IAR_LU_biexp + D_fitted_IAR_LU_modified_mix + D_fitted_IAR_LU_modified_topopro + D_fitted_IAR_LU_segmented_2step + D_fitted_IAR_LU_segmented_3step + D_fitted_IAR_LU_subtracted, data=data_new_wide)
-Dp_model <- lm(Dp ~ Dp_fitted_IAR_LU_biexp + Dp_fitted_IAR_LU_modified_mix + Dp_fitted_IAR_LU_modified_topopro + Dp_fitted_IAR_LU_segmented_2step + Dp_fitted_IAR_LU_segmented_3step + Dp_fitted_IAR_LU_subtracted, data=data_wide)
-# Dp_model <- lm(Dp ~ Dp_fitted_ETP_SRI_LinearFitting + Dp_fitted_IAR_LU_biexp + Dp_fitted_IAR_LU_modified_mix + Dp_fitted_IAR_LU_modified_topopro + Dp_fitted_IAR_LU_segmented_2step + Dp_fitted_IAR_LU_segmented_3step + Dp_fitted_IAR_LU_subtracted, data=data_new_wide)
-# predict new data from existing model
-predict(object = f_model, newdata = data_wide)
+if (runPrediction) {
+    # Then this widens it so we can lm()
+    data_wide <- data %>% pivot_wider(names_from=Algorithm, values_from=c(f_fitted, Dp_fitted, D_fitted), id_cols=c(Region, SNR, index, f, D, Dp))
+    # linear fit for f
+    f_model <- lm(f ~ f_fitted_IAR_LU_biexp + f_fitted_IAR_LU_modified_mix + f_fitted_IAR_LU_modified_topopro + f_fitted_IAR_LU_segmented_2step + f_fitted_IAR_LU_segmented_3step + f_fitted_IAR_LU_subtracted, data=data_wide)
+    #f_model <- lm(f ~ f_fitted_ETP_SRI_LinearFitting + f_fitted_IAR_LU_biexp + f_fitted_IAR_LU_modified_mix + f_fitted_IAR_LU_modified_topopro + f_fitted_IAR_LU_segmented_2step + f_fitted_IAR_LU_segmented_3step + f_fitted_IAR_LU_subtracted, data=data_new_wide)
+    D_model <- lm(D ~ D_fitted_ETP_SRI_LinearFitting + D_fitted_IAR_LU_biexp + D_fitted_IAR_LU_modified_topopro + D_fitted_IAR_LU_segmented_2step + D_fitted_IAR_LU_segmented_3step + D_fitted_IAR_LU_subtracted, data=data_wide)
+    # D_model <- lm(D ~ D_fitted_ETP_SRI_LinearFitting + D_fitted_IAR_LU_biexp + D_fitted_IAR_LU_modified_mix + D_fitted_IAR_LU_modified_topopro + D_fitted_IAR_LU_segmented_2step + D_fitted_IAR_LU_segmented_3step + D_fitted_IAR_LU_subtracted, data=data_new_wide)
+    Dp_model <- lm(Dp ~ Dp_fitted_IAR_LU_biexp + Dp_fitted_IAR_LU_modified_mix + Dp_fitted_IAR_LU_modified_topopro + Dp_fitted_IAR_LU_segmented_2step + Dp_fitted_IAR_LU_segmented_3step + Dp_fitted_IAR_LU_subtracted, data=data_wide)
+    # Dp_model <- lm(Dp ~ Dp_fitted_ETP_SRI_LinearFitting + Dp_fitted_IAR_LU_biexp + Dp_fitted_IAR_LU_modified_mix + Dp_fitted_IAR_LU_modified_topopro + Dp_fitted_IAR_LU_segmented_2step + Dp_fitted_IAR_LU_segmented_3step + Dp_fitted_IAR_LU_subtracted, data=data_new_wide)
+    # predict new data from existing model
+    predict(object = f_model, newdata = data_wide)
+}
 
 
 ivim_decay <- function(f, D, Dp, bvalues) {
@@ -63,4 +81,4 @@ curves_restricted <- generate_curves(data_restricted, bvalues, "")
 data_curves_restricted <- cbind(curves_restricted, signal_fitted=curves_restricted_fitted$signal_fitted)
 curve_plot <- ggplot(data_curves_restricted, aes(x=bvalue))  + facet_grid(Region ~ SNR) + geom_line(alpha=0.2, aes(y=signal_fitted, group=interaction(Algorithm, index), color=Algorithm)) + geom_line(aes(y=signal)) + ylim(0, 1) + ylab("Signal (a.u.)")
 print(curve_plot)
-ggsave("curve_plot.pdf", plot=curve_plot, width = 30, height = 30, units = "cm")
+ggsave("curve_plot.pdf", plot=curve_plot, width = 30, height = 30, units = "cm")