library(mrgsolve)
##
## Attaching package: 'mrgsolve'
## The following object is masked from 'package:stats':
##
## filter
library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --
## v ggplot2 3.3.3 v purrr 0.3.4
## v tibble 3.0.6 v dplyr 1.0.4
## v tidyr 1.1.2 v stringr 1.4.0
## v readr 1.4.0 v forcats 0.5.0
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks mrgsolve::filter(), stats::filter()
## x dplyr::lag() masks stats::lag()
library(ggplot2)
library(readr)
dat <- read_csv("data/011921.csv")
##
## -- Column specification --------------------------------------------------------
## cols(
## time = col_double(),
## conc = col_double(),
## evid = col_double(),
## cmt = col_double(),
## ID = col_double(),
## amt = col_double()
## )
dat
## # A tibble: 32 x 6
## time conc evid cmt ID amt
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 4.96 1 2 1 500
## 2 3 11.2 0 0 1 NA
## 3 10 11.5 0 0 1 NA
## 4 15 9.67 0 0 1 NA
## 5 30 8.51 0 0 1 NA
## 6 60 10.4 0 0 1 NA
## 7 1440 2.71 0 0 1 NA
## 8 2880 0.642 0 0 1 NA
## 9 0 0.436 1 2 2 500
## 10 3 9.32 0 0 2 NA
## # ... with 22 more rows
theme_set(
theme_bw() +
theme(legend.position = "top")
)
dat <- read_csv("data/011921.csv")
##
## -- Column specification --------------------------------------------------------
## cols(
## time = col_double(),
## conc = col_double(),
## evid = col_double(),
## cmt = col_double(),
## ID = col_double(),
## amt = col_double()
## )
# dat_1 <- dat[which(dat$ID<4),]
p <- ggplot(dat, aes(x=time, y=conc, colour=factor(ID))) +
geom_point(size=2) +
geom_smooth(method = lm, formula = y ~ splines::bs(x, 9), se = FALSE) +
labs(title = "Plot all data points")
p
## Warning in predict.lm(model, newdata = new_data_frame(list(x = xseq)), se.fit =
## se, : prediction from a rank-deficient fit may be misleading
## Warning in predict.lm(model, newdata = new_data_frame(list(x = xseq)), se.fit =
## se, : prediction from a rank-deficient fit may be misleading
## Warning in predict.lm(model, newdata = new_data_frame(list(x = xseq)), se.fit =
## se, : prediction from a rank-deficient fit may be misleading
## Warning in predict.lm(model, newdata = new_data_frame(list(x = xseq)), se.fit =
## se, : prediction from a rank-deficient fit may be misleading
dat_2 <- dat[which(dat$time<61),]
p <- ggplot(dat_2, aes(x=time, y=conc, colour=factor(ID))) +
geom_point(size=2) +
geom_smooth(method="auto", se=FALSE) +
labs(title = "Only keep data within 60 minutes") +
theme_bw()
p
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
source("previous_implementation.R")
fitting_func <- function(dat_to_fit, num_of_compartments=1, approach='bolus') {
formula_df <- formula_func()
experiment <- formula_df[formula_df[["num_of_compartments"]]==num_of_compartments &
formula_df[["delivery_approach"]]==approach,]
nls_formula <- as.formula(experiment[['formula_str']])
paras <- init_para(dat_to_fit, num_of_compartments)
nls_out <- NULL
predicted_t <- seq(from = dat_to_fit$t[1], to = dat_to_fit$t[nrow(dat_to_fit)], length.out = 100)
predicted_c <- NULL
if (num_of_compartments==1 & approach=='bolus'){
start_point <- c(paras[['init_para']][['init_val']])
start_point <- list(a = start_point[1], b = start_point[2])
}
else if (num_of_compartments==2 & approach=='bolus') {
start_point <- c(paras[['init_para']][['init_val']])
start_point <- list(a = start_point[1],
b = start_point[2],
c = start_point[3],
d = start_point[4])
}
else if (num_of_compartments==3 & approach=='bolus') {
start_point <- c(paras[['init_para']][['init_val']])
start_point <- list(a = start_point[1],
b = start_point[2],
c = start_point[3],
d = start_point[4],
e = start_point[5],
g = start_point[6])
}
nls_out <- try(nls(nls_formula, data = dat_to_fit, lower = paras[['init_para']][['min_val']], upper = paras[['init_para']][['max_val']],
start = start_point, algorithm = "default"), silent = TRUE)
if (class(nls_out)=="try-error") {
fn <- nls_bolus_error_solution(dat_to_fit, num_of_compartments, approach)
# print(fn)
nls_out <- try(nloptr::bobyqa(as.vector(unlist(start_point)), fn, lower = paras[['init_para']][['min_val']], upper = paras[['init_para']][['max_val']]),
silent = TRUE)
}
return(nls_out)
}
Confidenti Interval
One compartment
# colnames(dat_2)[1:9] <- c('t(h)', "con(mg/L)", "dose(mg)", 'parameters', "min", "max", "initial",
# 'weighting_approach', "specified_weight")
data_fit<- dat[c('time', 'conc', 'amt', 'ID')]
colnames(data_fit)[1:4]<-c('t',"con",'dose',"ID")
data_fit
## # A tibble: 32 x 4
## t con dose ID
## <dbl> <dbl> <dbl> <dbl>
## 1 0 4.96 500 1
## 2 3 11.2 NA 1
## 3 10 11.5 NA 1
## 4 15 9.67 NA 1
## 5 30 8.51 NA 1
## 6 60 10.4 NA 1
## 7 1440 2.71 NA 1
## 8 2880 0.642 NA 1
## 9 0 0.436 500 2
## 10 3 9.32 NA 2
## # ... with 22 more rows
myList <- unlist(unique(data_fit['ID']))
fit_output <- lapply(myList, function(id){
print(id)
temp_output <- fitting_func(data_fit[which(data_fit['ID']==id),])
temp_output <- coef(temp_output)
print(temp_output)
return(temp_output)
})
## [1] 1
## Warning in nls(nls_formula, data = dat_to_fit, lower = paras[["init_para"]]
## [["min_val"]], : upper and lower bounds ignored unless algorithm = "port"
## a b
## 9.5132539162 0.0008483735
## [1] 2
## Warning in nls(nls_formula, data = dat_to_fit, lower = paras[["init_para"]]
## [["min_val"]], : upper and lower bounds ignored unless algorithm = "port"
## a b
## 8.1330106540 0.0003043876
## [1] 3
## Warning in nls(nls_formula, data = dat_to_fit, lower = paras[["init_para"]]
## [["min_val"]], : upper and lower bounds ignored unless algorithm = "port"
## a b
## 1.201149e+01 6.970941e-04
## [1] 4
## Warning in nls(nls_formula, data = dat_to_fit, lower = paras[["init_para"]]
## [["min_val"]], : upper and lower bounds ignored unless algorithm = "port"
## a b
## 1.715843e+01 9.937707e-04
fit_output
## $ID1
## a b
## 9.5132539162 0.0008483735
##
## $ID2
## a b
## 8.1330106540 0.0003043876
##
## $ID3
## a b
## 1.201149e+01 6.970941e-04
##
## $ID4
## a b
## 1.715843e+01 9.937707e-04
fit_output <- data.frame(fit_output)
one_cp <- data.frame(t(apply(as.matrix(fit_output), 1, function(x){
temp_result <- c(mean(x), sd(x))
temp_result <- c(temp_result, mean(x)+c(-1.96, 1.96)*sd(x)/sqrt(length(x)))
})))
colnames(one_cp) <- c('Mean', "STD", 'LCI 95%', 'UCI 95')
rownames(one_cp) <- c('A', 'k10')
one_cp
## Mean STD LCI 95% UCI 95
## A 1.170405e+01 3.9747828691 7.8087597707 15.599334194
## k10 7.109065e-04 0.0002968488 0.0004199946 0.001001818
Two compartment
# colnames(dat_2)[1:9] <- c('t(h)', "con(mg/L)", "dose(mg)", 'parameters', "min", "max", "initial",
# 'weighting_approach', "specified_weight")
myList <- unlist(unique(data_fit['ID']))
fit_output <- lapply(myList, function(id){
print(id)
temp_output <- fitting_func(data_fit[which(data_fit['ID']==id),], num_of_compartments=2)
print(temp_output$par)
return(temp_output$par)
})
## [1] 1
## Warning in if (is.na(f0)) {: the condition has length > 1 and only the first
## element will be used
## [1] 3.251383 0.000000 1.703862 0.000000
## [1] 2
## Warning in if (is.na(f0)) {: the condition has length > 1 and only the first
## element will be used
## [1] 0.2880040 0.3162706 0.1482972 0.0000000
## [1] 3
## Warning in if (is.na(f0)) {: the condition has length > 1 and only the first
## element will be used
## [1] 0.2722648781 0.1102392367 0.0005073736 0.0138861430
## [1] 4
## Warning in if (is.na(f0)) {: the condition has length > 1 and only the first
## element will be used
## [1] 0.277072952 0.187542965 0.007364205 0.000000000
fit_output
## $ID1
## [1] 3.251383 0.000000 1.703862 0.000000
##
## $ID2
## [1] 0.2880040 0.3162706 0.1482972 0.0000000
##
## $ID3
## [1] 0.2722648781 0.1102392367 0.0005073736 0.0138861430
##
## $ID4
## [1] 0.277072952 0.187542965 0.007364205 0.000000000
fit_output <- data.frame(fit_output)
two_cp <- data.frame(t(apply(as.matrix(fit_output), 1, function(x){
temp_result <- c(mean(x), sd(x))
temp_result <- c(temp_result, mean(x)+c(-1.96, 1.96)*sd(x)/sqrt(length(x)))
})))
colnames(two_cp) <- c('Mean', "STD", 'LCI 95%', 'UCI 95')
rownames(two_cp) <- c('A', 'Lambda1', 'B', 'lambda2')
two_cp
## Mean STD LCI 95% UCI 95
## A 1.022181297 1.486149291 -0.434245008 2.47860760
## Lambda1 0.153513203 0.133024983 0.023148720 0.28387769
## B 0.465007759 0.828706683 -0.347124790 1.27714031
## lambda2 0.003471536 0.006943071 -0.003332674 0.01027575
Three compartment
# colnames(dat_2)[1:9] <- c('t(h)', "con(mg/L)", "dose(mg)", 'parameters', "min", "max", "initial",
# 'weighting_approach', "specified_weight")
data_fit<- dat[c('time', 'conc', 'amt', 'ID')]
colnames(data_fit)[1:4]<-c('t',"con",'dose',"ID")
data_fit
## # A tibble: 32 x 4
## t con dose ID
## <dbl> <dbl> <dbl> <dbl>
## 1 0 4.96 500 1
## 2 3 11.2 NA 1
## 3 10 11.5 NA 1
## 4 15 9.67 NA 1
## 5 30 8.51 NA 1
## 6 60 10.4 NA 1
## 7 1440 2.71 NA 1
## 8 2880 0.642 NA 1
## 9 0 0.436 500 2
## 10 3 9.32 NA 2
## # ... with 22 more rows
myList <- unlist(unique(data_fit['ID']))
fit_output <- lapply(myList, function(id){
print(id)
temp_output <- fitting_func(data_fit[which(data_fit['ID']==id),], num_of_compartments=3)
print(temp_output$par)
return(temp_output$par)
})
## [1] 1
## Warning in nls(nls_formula, data = dat_to_fit, lower = paras[["init_para"]]
## [["min_val"]], : upper and lower bounds ignored unless algorithm = "port"
## Warning in if (is.na(f0)) {: the condition has length > 1 and only the first
## element will be used
## [1] 2.2415476248 0.0000000000 2.6495788425 0.0008751636 0.0641192393
## [6] 0.0540710000
## [1] 2
## Warning in nls(nls_formula, data = dat_to_fit, lower = paras[["init_para"]]
## [["min_val"]], : upper and lower bounds ignored unless algorithm = "port"
## Warning in nls(nls_formula, data = dat_to_fit, lower = paras[["init_para"]]
## [["min_val"]], : the condition has length > 1 and only the first element will be
## used
## [1] 0.43630121 0.00000000 0.00000000 0.00240476 0.00000000 0.05407100
## [1] 3
## Warning in nls(nls_formula, data = dat_to_fit, lower = paras[["init_para"]]
## [["min_val"]], : upper and lower bounds ignored unless algorithm = "port"
## Warning in nls(nls_formula, data = dat_to_fit, lower = paras[["init_para"]]
## [["min_val"]], : the condition has length > 1 and only the first element will be
## used
## [1] 0.2723774174 0.2024437342 0.0003949654 0.0004422807 0.0000000000
## [6] 0.0540710000
## [1] 4
## Warning in if (is.na(f0)) {: the condition has length > 1 and only the first
## element will be used
## [1] 2.724895e-01 2.746801e+00 1.102961e-04 3.507245e-05 1.183740e-02
## [6] 5.407100e-02
fit_output
## $ID1
## [1] 2.2415476248 0.0000000000 2.6495788425 0.0008751636 0.0641192393
## [6] 0.0540710000
##
## $ID2
## [1] 0.43630121 0.00000000 0.00000000 0.00240476 0.00000000 0.05407100
##
## $ID3
## [1] 0.2723774174 0.2024437342 0.0003949654 0.0004422807 0.0000000000
## [6] 0.0540710000
##
## $ID4
## [1] 2.724895e-01 2.746801e+00 1.102961e-04 3.507245e-05 1.183740e-02
## [6] 5.407100e-02
fit_output <- data.frame(fit_output)
three_cp <- data.frame(t(apply(as.matrix(fit_output), 1, function(x){
temp_result <- c(mean(x), sd(x))
temp_result <- c(temp_result, mean(x)+c(-1.96, 1.96)*sd(x)/sqrt(length(x)))
})))
colnames(three_cp) <- c('Mean', "STD", 'LCI 95%', 'UCI 95')
rownames(three_cp) <- c('A', 'Lambda1', 'B', 'lambda2', 'C', 'lambda3')
three_cp
## Mean STD LCI 95% UCI 95
## A 0.8056789270 0.960357629 -1.354715e-01 1.74682940
## Lambda1 0.7373110903 1.343054550 -5.788824e-01 2.05350455
## B 0.6625210260 1.324705221 -6.356901e-01 1.96073214
## lambda2 0.0009393191 0.001035429 -7.540152e-05 0.00195404
## C 0.0189891591 0.030599827 -1.099867e-02 0.04897699
## lambda3 0.0540710000 0.000000000 5.407100e-02 0.05407100