Fixed Interpolation/Extrapolation for Partial AUCs

- Fixed an 'subscript out of bounds' issue for CEST Interpolation/Extrapolation - Corrected the logic for Interpolation/Extrpolation - Identified Issue for OPTIMIZEKEL for SS dosing types and Partially Implemented it - Need to account for OPTIMIZEKEL for SS dosing types based on parameter list (Will be pushed in the next commit)
tensfeldt · Nov 13, 2019 · e3be44d · e3be44d
1 parent a06c1c7
commit e3be44d
Show file tree

Hide file tree

Showing 12 changed files with 88 additions and 42 deletions.
diff --git a/openNCA/R/auc_lin.R b/openNCA/R/auc_lin.R
@@ -114,6 +114,7 @@ auc_lin <- function(conc = NULL, time = NULL, exflag = NULL, interpolate = NULL,
 ##
       est_tmp <- estimate_missing_concentration(conc = conc, time = time, auc_method = "LIN", model = model, dosing_type = dosing_type, told = told, orig_conc = orig_conc, orig_time = orig_time)
       conc <- est_tmp[[1]]
+      tmp <- data.frame(time, conc)
     }
     for(i in 1:(nrow(tmp)-1)){
       auc_df[i] <- ((conc[i] + conc[i+1])/2)*(time[i+1]-time[i])

diff --git a/openNCA/R/auc_lin_log.R b/openNCA/R/auc_lin_log.R
@@ -121,6 +121,7 @@ auc_lin_log <- function(conc = NULL, time = NULL, exflag = NULL, t_max = NULL, i
 ##
         est_tmp <- estimate_missing_concentration(conc = conc, time = time, auc_method = "LIN", model = model, dosing_type = dosing_type, told = told, orig_conc = orig_conc, orig_time = orig_time)
         conc <- est_tmp[[1]]
+        tmp <- data.frame(time, conc)
       }
       for(i in 1:(nrow(tmp)-1)){
         if(tmp$time[i+1] <= t_max || tmp$conc[i] == 0 || tmp$conc[i+1] == 0 || tmp$conc[i] == tmp$conc[i+1]){

diff --git a/openNCA/R/auc_lin_up_log_down.R b/openNCA/R/auc_lin_up_log_down.R
@@ -117,6 +117,7 @@ auc_lin_up_log_down <- function(conc = NULL, time = NULL, exflag = NULL, interpo
 ##
       est_tmp <- estimate_missing_concentration(conc = conc, time = time, auc_method = "LIN", model = model, dosing_type = dosing_type, told = told, orig_conc = orig_conc, orig_time = orig_time)
       conc <- est_tmp[[1]]
+      tmp <- data.frame(time, conc)
     }
     for(i in 1:(nrow(tmp)-1)){
       curr_c <- as.numeric(conc[i])

diff --git a/openNCA/R/auc_log.R b/openNCA/R/auc_log.R
@@ -116,6 +116,7 @@ auc_log <- function(conc = NULL, time = NULL, exflag = NULL, interpolate = NULL,
 ##
       est_tmp <- estimate_missing_concentration(conc = conc, time = time, auc_method = "LIN", model = model, dosing_type = dosing_type, told = told, orig_conc = orig_conc, orig_time = orig_time)
       conc <- est_tmp[[1]]
+      tmp <- data.frame(time, conc)
     }
     if(!is.na(t_max)){
       for(i in 1:(nrow(tmp)-1)){

diff --git a/openNCA/R/auc_t1_t2.R b/openNCA/R/auc_t1_t2.R
@@ -150,7 +150,7 @@ auc_t1_t2 <- function(conc = NULL, time = NULL, t1 = NULL, t2 = NULL, method = 1
   } else if(is.na(t_max) && method==1) { # 2019-09-11/TGT/ - update 2019-09-25/TGT/ updated to also consider method
       return(NA)
   }
-
+  
   if(!(is.numeric(conc) && is.vector(conc)) ){
     stop("Error in auc_t1_t2: 'conc' is not a numeric vector")
   }
@@ -193,13 +193,27 @@ auc_t1_t2 <- function(conc = NULL, time = NULL, t1 = NULL, t2 = NULL, method = 1
   }
   conc <- conc[time >= t1 & time <= t2]
   time <- time[time >= t1 & time <= t2]
-
+
+##  2019-11-12/RD Added to Interpolate/Extrapolate data properly for missing time values
+## 
+  if(isTRUE(interpolate) || isTRUE(extrapolate)){
+    if(!(t1 %in% time)){
+      time <- c(t1, time)
+      conc <- c(NA, conc)
+    }
+    if(!(t2 %in% time)){
+      time <- c(time, t2)
+      conc <- c(conc, NA)
+    }
+  } else {
+    if(sum(conc, na.rm = T) == 0){
+      return(0)
+    }
+  }
+
   if(length(time) != length(conc)){
     stop("Error in auc_t1_t2: length of 'time' and 'conc' vectors are not equal")
   }
-  if(sum(conc, na.rm = T) == 0){
-    return(0)
-  }
 
   if(method == 1){
     return(auc_lin_log(conc = conc, time = time, exflag = exflag, t_max = t_max, interpolate = interpolate, extrapolate = extrapolate, model = model, dosing_type = dosing_type, told = told, orig_conc = orig_conc, orig_time = orig_time))

diff --git a/openNCA/R/estimate_missing_concentration.R b/openNCA/R/estimate_missing_concentration.R
@@ -52,11 +52,12 @@ estimate_missing_concentration <- function(conc = NULL, time = NULL, auc_method
           tmp$INT_EXT[1] <- "INT"
         }
       } else if((orig_time[1] < time[1]) && (time[1] < orig_time[length(orig_time)])){
-        idx <- which(time[1] == orig_time)
+        idx <- which(orig_time < time[1])
+        idx <- idx[length(idx)]
         if(auc_method == "LIN"){
-          conc[1] <- interpolate_lin(conc1 = orig_conc[idx-1], time1 = orig_time[idx-1], conc2 = orig_conc[idx+1], time2 = orig_time[i+1], est_time = time[1])
+          conc[1] <- interpolate_lin(conc1 = orig_conc[idx], time1 = orig_time[idx], conc2 = orig_conc[idx+1], time2 = orig_time[idx+1], est_time = time[1])
         } else if(auc_method == "LOG"){
-          conc[1] <- interpolate_log(conc1 = orig_conc[idx-1], time1 = orig_time[idx-1], conc2 = orig_conc[idx+1], time2 = orig_time[i+1], est_time = time[1])
+          conc[1] <- interpolate_log(conc1 = orig_conc[idx], time1 = orig_time[idx], conc2 = orig_conc[idx+1], time2 = orig_time[idx+1], est_time = time[1])
         }
         tmp$INT_EXT[1] <- "INT"
       } else if(time[1] >= orig_time[length(orig_time)]){
@@ -72,11 +73,12 @@ estimate_missing_concentration <- function(conc = NULL, time = NULL, auc_method
 ##  
     } else if(is.na(conc[nrow(tmp)]) && (i+1) == nrow(tmp)){
       if((orig_time[1] < time[nrow(tmp)]) && (time[nrow(tmp)] < orig_time[length(orig_time)])){
-        idx <- which(time[nrow(tmp)] == orig_time)
+        idx <- which(orig_time < time[nrow(tmp)])
+        idx <- idx[length(idx)]
         if(auc_method == "LIN"){
-          conc[nrow(tmp)] <- interpolate_lin(conc1 = orig_conc[idx-1], time1 = orig_time[idx-1], conc2 = orig_conc[idx+1], time2 = orig_time[i+1], est_time = time[nrow(tmp)])
+          conc[nrow(tmp)] <- interpolate_lin(conc1 = orig_conc[idx], time1 = orig_time[idx], conc2 = orig_conc[idx+1], time2 = orig_time[idx+1], est_time = time[nrow(tmp)])
         } else if(auc_method == "LOG"){
-          conc[nrow(tmp)] <- interpolate_log(conc1 = orig_conc[idx-1], time1 = orig_time[idx-1], conc2 = orig_conc[idx+1], time2 = orig_time[i+1], est_time = time[nrow(tmp)])
+          conc[nrow(tmp)] <- interpolate_log(conc1 = orig_conc[idx], time1 = orig_time[idx], conc2 = orig_conc[idx+1], time2 = orig_time[idx+1], est_time = time[nrow(tmp)])
         }
         tmp$INT_EXT[nrow(tmp)] <- "INT"
       } else {

diff --git a/openNCA/R/run_M1_SD_computation.R b/openNCA/R/run_M1_SD_computation.R
@@ -1321,10 +1321,6 @@ run_M1_SD_computation <- function(data = NULL, map = NULL, method = 1, model_reg
         } else {
           cest_kel <- rep(NA, length(conc))
         }
-##        2019-11-08/RD/ Added to add cest interpolation/extrapolation data as an output
-##
-
-        end_idx <- which(tmp_df[,map_data$TIME] == cest_tmp$TIME[length(cest_tmp$TIME)])
 
 ##################################################################################################################################
 ### Computations are completed at this point
@@ -1342,11 +1338,13 @@ run_M1_SD_computation <- function(data = NULL, map = NULL, method = 1, model_reg
                 while(cest_idx <= nrow(cest_tmp)){
                   if(cest_tmp[cest_idx, "TIME"] <= time[e]){
                     if(!(cest_tmp[cest_idx, "TIME"] %in% time)){
-                      pkdr_idx <- which(tmp_df[,map_data$TIME] == cest_tmp[cest_idx, "TIME"])
+##                      2019-11-12/RD/ Commented the way to retrieve the PKDATAROWID, need to find any alternative
+##                      
+                      #pkdr_idx <- which(tmp_df[,map_data$TIME] == cest_tmp[cest_idx, "TIME"])
                       if(cest_tmp[cest_idx, "INT_EXT"] == "INT"){
-                        tmp_est_row <- c(tmp_df[pkdr_idx,"PKDATAROWID"], unique(data_data[,map_data$SDEID])[i], cest_tmp[cest_idx, "TIME"], NA, cest_tmp[cest_idx, "CONC"], NA, NA, NA)
+                        tmp_est_row <- c(NA, unique(data_data[,map_data$SDEID])[i], cest_tmp[cest_idx, "TIME"], NA, cest_tmp[cest_idx, "CONC"], NA, NA, NA)
                       } else if(cest_tmp[cest_idx, "INT_EXT"] == "EXT"){
-                        tmp_est_row <- c(tmp_df[pkdr_idx,"PKDATAROWID"], unique(data_data[,map_data$SDEID])[i], cest_tmp[cest_idx, "TIME"], NA, NA, cest_tmp[cest_idx, "CONC"], NA, NA)
+                        tmp_est_row <- c(NA, unique(data_data[,map_data$SDEID])[i], cest_tmp[cest_idx, "TIME"], NA, NA, cest_tmp[cest_idx, "CONC"], NA, NA)
                       }
                       est_data[est_idx,] <- tmp_est_row
                       est_idx <- est_idx + 1
@@ -1360,6 +1358,8 @@ run_M1_SD_computation <- function(data = NULL, map = NULL, method = 1, model_reg
                         cest_idx <- cest_idx + 1
                       }
                     }
+                  } else {
+                    break
                   }
                 }
               }
@@ -1716,7 +1716,6 @@ run_M1_SD_computation <- function(data = NULL, map = NULL, method = 1, model_reg
     tmp_int_type <- computation_df[,names(computation_df) == as.character(regular_int_type[n])]
     if(!is.null(ncol(tmp_int_type))){
       for(r in 1:length(tmp_int_type)){
-        print(computation_df[,names(computation_df) == as.character(regular_int_type[n])][,r])
         suppressWarnings(computation_df[,names(computation_df) == as.character(regular_int_type[n])][,r] <- as.numeric(computation_df[,names(computation_df) == as.character(regular_int_type[n])][,r]))
       }
     } else {

diff --git a/openNCA/R/run_M1_SS_computation.R b/openNCA/R/run_M1_SS_computation.R
@@ -1011,7 +1011,7 @@ run_M1_SS_computation <- function(data = NULL, map = NULL, method = 1, model_reg
     warning("Kel optimization cannot be performed because 'TMAX', 'TLAST', 'CMAX', 'CLAST', 'AUCLAST' are not part of the calulcated parameters AND Flag 'FLGACCEPTKELCRIT' and Flag 'FLGXKEL' are not present in the dataset")
   }
 
-  if(optimize_kel && "TMAX" %in% parameter_list && "TLAST" %in% parameter_list && "CMAX" %in% parameter_list && "CLAST" %in% parameter_list && "AUCLAST" %in% parameter_list &&
+  if(optimize_kel && "TMAXi" %in% parameter_list && "TLAST" %in% parameter_list && "CMAXi" %in% parameter_list && "CLASTi" %in% parameter_list && "AUCLAST" %in% parameter_list &&
      "FLGACCEPTKELCRIT" %in% names(map_data) && "FLGEXKEL" %in% names(map_data) && map_data$FLGEXKEL %in% names(data_data)){
     kel_flag_optimized <- integer()
   }
@@ -1296,14 +1296,16 @@ run_M1_SS_computation <- function(data = NULL, map = NULL, method = 1, model_reg
         if(comp_required[["AUCLAST"]]){
           auclast <- auc_last(conc = tmp_df[,map_data$CONC], time = tmp_df[,map_data$TIME], method = method, exflag = auc_flag, t_last = t_last, t_max = t_max)
         }
-        if(optimize_kel && "TMAX" %in% parameter_list && "TLAST" %in% parameter_list && "CMAX" %in% parameter_list && "CLAST" %in% parameter_list && "AUCLAST" %in% parameter_list &&
+        if(optimize_kel && "TMAXi" %in% parameter_list && "TLAST" %in% parameter_list && "CMAXi" %in% parameter_list && "CLASTi" %in% parameter_list && "AUCLAST" %in% parameter_list &&
            "FLGACCEPTKELCRIT" %in% names(map_data) && "FLGEXKEL" %in% names(map_data) && map_data$FLGEXKEL %in% names(data_data)){
 ### 2019-09-04/TGT/
 ###          orig_time <- tmp_df[,map_data$TIME]
 ###          orig_conc <- tmp_df[,map_data$CONC]
           tmp_time <- orig_time
           tmp_conc <- orig_conc
 
+          #print(tmp_time)
+          #print(tmp_conc)
           if("FLGNOCMAX" %in% names(map_data) && (map_data$FLGNOCMAX == 1 || map_data$FLGNOCMAX == 0)){
             flg_no_cmax <- as.logical(map_data$FLGNOCMAX)
             if(flg_no_cmax){
@@ -1341,6 +1343,10 @@ run_M1_SS_computation <- function(data = NULL, map = NULL, method = 1, model_reg
               tmp_conc <- orig_conc[s_conc:e_conc]
             }
           }
+          #print("new data")
+          #print(tmp_time)
+          #print(tmp_conc)
+          #print("------")
 
           kel_n <- as.numeric(flag_df$CRIT[match("KELNOPT", flag_df$VAR)])
           kel_op <- flag_df$OPR[match("KELNOPT", flag_df$VAR)]
@@ -1357,6 +1363,7 @@ run_M1_SS_computation <- function(data = NULL, map = NULL, method = 1, model_reg
               ulist <- c(ulist,tlist)
             }
           }
+          #print(ulist)
 
           kelr_val <- kel_r(conc = tmp_df[,map_data$CONC], time = tmp_df[,map_data$TIME])[["KELRSQ"]]
           if("AUCXPCTO" %in% flag_df$VAR){
@@ -1367,11 +1374,14 @@ run_M1_SS_computation <- function(data = NULL, map = NULL, method = 1, model_reg
             stop("Error in optimize kel")
           }
 
+          #print("for loop")
           selected_idx <- NA
           saved_kel_opt <- 0
           for(k in 1:length(ulist)){
             sel_time <- ulist[[k]]
             sel_conc <- tmp_conc[match(sel_time, tmp_time)]
+            #print(sel_time)
+            #print(sel_conc)
 
             kelr_opt <- kel_r(conc = sel_conc, time = sel_time)[["KELRSQ"]]
             if("AUCXPCTO" %in% flag_df$VAR){
@@ -1398,6 +1408,8 @@ run_M1_SS_computation <- function(data = NULL, map = NULL, method = 1, model_reg
                 selected_idx <- match(sel_time, orig_time)
               }
             }
+            #print("selected KEL Flags")
+            #print(selected_idx)
           }
           tmp_kel_flag <- rep(1, length(kel_flag))
           tmp_kel_flag[selected_idx] <- 0
@@ -2059,11 +2071,13 @@ run_M1_SS_computation <- function(data = NULL, map = NULL, method = 1, model_reg
                 while(cest_idx <= nrow(cest_tmp)){
                   if(cest_tmp[cest_idx, "TIME"] <= time[e]){
                     if(!(cest_tmp[cest_idx, "TIME"] %in% time)){
-                      pkdr_idx <- which(tmp_df[,map_data$TIME] == cest_tmp[cest_idx, "TIME"])
+##                      2019-11-12/RD/ Commented the way to retrieve the PKDATAROWID, need to find any alternative
+##                      
+                      #pkdr_idx <- which(tmp_df[,map_data$TIME] == cest_tmp[cest_idx, "TIME"])
                       if(cest_tmp[cest_idx, "INT_EXT"] == "INT"){
-                        tmp_est_row <- c(tmp_df[pkdr_idx,"PKDATAROWID"], unique(data_data[,map_data$SDEID])[i], cest_tmp[cest_idx, "TIME"], NA, cest_tmp[cest_idx, "CONC"], NA, NA, NA)
+                        tmp_est_row <- c(NA, unique(data_data[,map_data$SDEID])[i], cest_tmp[cest_idx, "TIME"], NA, cest_tmp[cest_idx, "CONC"], NA, NA, NA)
                       } else if(cest_tmp[cest_idx, "INT_EXT"] == "EXT"){
-                        tmp_est_row <- c(tmp_df[pkdr_idx,"PKDATAROWID"], unique(data_data[,map_data$SDEID])[i], cest_tmp[cest_idx, "TIME"], NA, NA, cest_tmp[cest_idx, "CONC"], NA, NA)
+                        tmp_est_row <- c(NA, unique(data_data[,map_data$SDEID])[i], cest_tmp[cest_idx, "TIME"], NA, NA, cest_tmp[cest_idx, "CONC"], NA, NA)
                       }
                       est_data[est_idx,] <- tmp_est_row
                       est_idx <- est_idx + 1
@@ -2077,6 +2091,8 @@ run_M1_SS_computation <- function(data = NULL, map = NULL, method = 1, model_reg
                         cest_idx <- cest_idx + 1
                       }
                     }
+                  } else {
+                    break
                   }
                 }
               }
@@ -2619,7 +2635,6 @@ run_M1_SS_computation <- function(data = NULL, map = NULL, method = 1, model_reg
     tmp_int_type <- computation_df[,names(computation_df) == as.character(regular_int_type[n])]
     if(!is.null(ncol(tmp_int_type))){
       for(r in 1:length(tmp_int_type)){
-        print(computation_df[,names(computation_df) == as.character(regular_int_type[n])][,r])
         suppressWarnings(computation_df[,names(computation_df) == as.character(regular_int_type[n])][,r] <- as.numeric(computation_df[,names(computation_df) == as.character(regular_int_type[n])][,r]))
       }
     } else {
@@ -2671,7 +2686,7 @@ run_M1_SS_computation <- function(data = NULL, map = NULL, method = 1, model_reg
   results_list$data_out <- computation_df
   results_list$est_data <- est_data
 
-  if(optimize_kel && "TMAX" %in% parameter_list && "TLAST" %in% parameter_list && "CMAX" %in% parameter_list && "CLAST" %in% parameter_list && "AUCLAST" %in% parameter_list &&
+  if(optimize_kel && "TMAXi" %in% parameter_list && "TLAST" %in% parameter_list && "CMAXi" %in% parameter_list && "CLASTi" %in% parameter_list && "AUCLAST" %in% parameter_list &&
      "FLGACCEPTKELCRIT" %in% names(map_data) && "FLGEXKEL" %in% names(map_data) && map_data$FLGEXKEL %in% names(data_data)){
 ###    if("KEL" %in% parameter_list){
 ###    if("KEL" %in% parameter_list){

diff --git a/openNCA/R/run_M2_SD_computation.R b/openNCA/R/run_M2_SD_computation.R
@@ -1308,11 +1308,13 @@ run_M2_SD_computation <- function(data = NULL, map = NULL, method = 1, model_reg
                 while(cest_idx <= nrow(cest_tmp)){
                   if(cest_tmp[cest_idx, "TIME"] <= time[e]){
                     if(!(cest_tmp[cest_idx, "TIME"] %in% time)){
-                      pkdr_idx <- which(tmp_df[,map_data$TIME] == cest_tmp[cest_idx, "TIME"])
+##                      2019-11-12/RD/ Commented the way to retrieve the PKDATAROWID, need to find any alternative
+##                      
+                      #pkdr_idx <- which(tmp_df[,map_data$TIME] == cest_tmp[cest_idx, "TIME"])
                       if(cest_tmp[cest_idx, "INT_EXT"] == "INT"){
-                        tmp_est_row <- c(tmp_df[pkdr_idx,"PKDATAROWID"], unique(data_data[,map_data$SDEID])[i], cest_tmp[cest_idx, "TIME"], NA, cest_tmp[cest_idx, "CONC"], NA, NA, NA)
+                        tmp_est_row <- c(NA, unique(data_data[,map_data$SDEID])[i], cest_tmp[cest_idx, "TIME"], NA, cest_tmp[cest_idx, "CONC"], NA, NA, NA)
                       } else if(cest_tmp[cest_idx, "INT_EXT"] == "EXT"){
-                        tmp_est_row <- c(tmp_df[pkdr_idx,"PKDATAROWID"], unique(data_data[,map_data$SDEID])[i], cest_tmp[cest_idx, "TIME"], NA, NA, cest_tmp[cest_idx, "CONC"], NA, NA)
+                        tmp_est_row <- c(NA, unique(data_data[,map_data$SDEID])[i], cest_tmp[cest_idx, "TIME"], NA, NA, cest_tmp[cest_idx, "CONC"], NA, NA)
                       }
                       est_data[est_idx,] <- tmp_est_row
                       est_idx <- est_idx + 1
@@ -1326,6 +1328,8 @@ run_M2_SD_computation <- function(data = NULL, map = NULL, method = 1, model_reg
                         cest_idx <- cest_idx + 1
                       }
                     }
+                  } else {
+                    break
                   }
                 }
               }
@@ -1674,7 +1678,6 @@ run_M2_SD_computation <- function(data = NULL, map = NULL, method = 1, model_reg
     tmp_int_type <- computation_df[,names(computation_df) == as.character(regular_int_type[n])]
     if(!is.null(ncol(tmp_int_type))){
       for(r in 1:length(tmp_int_type)){
-        print(computation_df[,names(computation_df) == as.character(regular_int_type[n])][,r])
         suppressWarnings(computation_df[,names(computation_df) == as.character(regular_int_type[n])][,r] <- as.numeric(computation_df[,names(computation_df) == as.character(regular_int_type[n])][,r]))
       }
     } else {