Task 1: Examine the data

In this study a new fast-acting drug “AMP1050” is being developed. Individuals with P. falciparum-malaria were randomized into a parallel group study where all subjects were given combination therapy with piperakine (long-acting antimalarial drug) and, either artemisinin or AMP1050 at different dose levels. During this study, several blood samples were taken from the individuals. The blood samples were analyzed and the number of parasites in the blood was measured as parasites/µL (the pharmacodynamics or PD). Concentration of the administered drug in the blood samples was also analyzed and measured as ng/mL (the pharmacokinetics or PK).

rm(list = ls())
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
library(ggplot2)
library(PKNCA)
## 
## Attaching package: 'PKNCA'
## The following object is masked from 'package:stats':
## 
##     filter
data <- read.csv("data_13.csv")

data$DRUG <- as.factor(data$DRUG)
data$Model <- as.factor(data$Model)
data$AMT <- as.numeric(data$AMT)
## Warning: NAs introduced by coercion
data$DV <- as.numeric(data$DV)
## Warning: NAs introduced by coercion
data$DV[is.na(data$DV)] <- 0

data %<>%
  group_by(ID) %>%
  fill(AMT, .direction = "updown") %>%
  ungroup()

data %<>%
  group_by(ID) %>%
  mutate(DRUG_cat = ifelse(DRUG == 1 & AMT == 500, "ART-500",
                      ifelse(DRUG == 1 & AMT == 100, "ART-100",
                             ifelse(DRUG == 2 & AMT == 100, "AMP-100",
                                    ifelse(DRUG == 2 & AMT == 200, "AMP-200",
                                           ifelse(DRUG == 2 & AMT == 700, "AMP-700",
                                                  ifelse(DRUG == 2 & AMT == 800, "AMP-800", as.factor(DRUG)))))))) 
ggplot(data,aes(x = DRUG_cat, y = study_arm, color = DRUG_cat)) +
  geom_point(size = 3, shape = 16) +
  labs(title = "Scatter Plot of Study Arm by Drug ",
       x = "Drug ",
       y = "Study Arm / ID") +
  theme_classic()

data %>%
  group_by(study_arm) %>%
  select(study_arm, side_effect,DV,FEMALE,AGE,HT,WT) %>%
  ungroup() %>%
  GGally::ggpairs()
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

data %>%
  group_by(study_arm) %>%
  filter(study_arm < 2 | study_arm > 17, Model == 2) %>%
  ggplot(aes(x = study_arm, y = DV, fill = factor(study_arm))) +
  geom_col() +
  theme_light() +
  labs(title = "Comparison of Emax across Study Arms",
       x = "Study Arm",
       y = "Effect")

data_comb <- data %>%
  filter(Time %in% c(0, 6), Model == 2) %>% 
  pivot_wider(names_from = Time, values_from = DV, names_prefix = "PAR_") %>% 
  mutate(delta_log_PAR = log(PAR_0) - log(PAR_6)) %>% 
  select(-c("PAR_0", "PAR_6")) 
## Warning: There were 2 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `delta_log_PAR = log(PAR_0) - log(PAR_6)`.
## ℹ In group 63: `ID = 63`.
## Caused by warning in `log()`:
## ! NaNs produced
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
data_D2  <- data_comb %>%
  filter(DRUG == 2)

Task 2: Standard statistical analysis

Different study arms that are medicating with the various doses of AMP1050 were compared with Artemisinin (study arm 1), by t-test. Results indicates that the lowest dose of AMP1050 with significance effect is AMP1050 100 mg.

P-value Significance level
< 0.001 ***
< 0.01 **
< 0.05 * 
ggplot(data_comb,aes(x = DRUG_cat, y = delta_log_PAR)) +
  geom_boxplot() +
  ggtitle(" effect of different treatment doses ") +
  xlab("Study arms") +
  ylab("Effect") +
  theme_minimal()
## Warning: Removed 2 rows containing non-finite values (`stat_boxplot()`).

pairwise.t.test(data_comb$delta_log_PAR, data_comb$DRUG_cat, p.adjust.method = "bonf", alternative = "less")
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  data_comb$delta_log_PAR and data_comb$DRUG_cat 
## 
##         AMP-100 AMP-200 AMP-700 AMP-800
## AMP-200 1       -       -       -      
## AMP-700 1       1       -       -      
## AMP-800 1       1       1       -      
## ART-500 1       1       5.3e-06 1.2e-07
## 
## P value adjustment method: bonferroni

AUC-calc

PK_data <- data %>%
  filter(Model == 1 ) %>%
  select(ID, DRUG, Time, DV, AMT, DRUG_cat)


data_conc <- PK_data %>%
  select(ID,Time, DV,DRUG) 

data_dose <- PK_data %>% filter(!is.na(AMT)) 
  

nca_conc <- PKNCAconc(data = data_conc,
                      formula = DV ~ Time | ID + DRUG)

nca_dose <- PKNCAdose(data = data_dose,
                      formula = AMT ~ Time | ID + DRUG)

nca_frame <- PKNCAdata(nca_conc,
                       nca_dose,
                       # NCA parameters & AUC interval specifications
                       intervals = data.frame(start = 0.,
                                             end = 6.,
                                             cmax = TRUE,
                                             auclast = TRUE))

nca_result <- pk.nca(nca_frame)

data_pk_params <- nca_result$result %>% 
  select("ID", "DRUG", Params = "PPTESTCD", Values = "PPORRES") %>% 
  pivot_wider(names_from = Params, values_from = Values)

data_comb %<>%
  group_by(ID) %>%
    merge(data_pk_params,by = c("ID" ,"DRUG"))
data_D2 %<>%
  group_by(ID) %>%
  merge(data_pk_params, by = c("ID" ,"DRUG"))

Task 3: Model based analysis

lm_auc seems to give best results. we proceed with it

library(rstanemax)
## Loading required package: Rcpp
#linear model:
lm_AMT <- lm(delta_log_PAR ~ AMT, data = data_D2)
summary(lm_AMT)
## 
## Call:
## lm(formula = delta_log_PAR ~ AMT, data = data_D2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.2333  -2.0342  -0.2412   1.8655  19.1917 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.863614   1.056533  -4.603 1.64e-05 ***
## AMT          0.028045   0.001923  14.585  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.355 on 76 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.7368, Adjusted R-squared:  0.7333 
## F-statistic: 212.7 on 1 and 76 DF,  p-value: < 2.2e-16
lm_cmax <- lm(delta_log_PAR ~ cmax, data = data_D2) 
summary(lm_cmax)  
## 
## Call:
## lm(formula = delta_log_PAR ~ cmax, data = data_D2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.6925  -1.8488  -0.6888   1.2744  18.2918 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.142592   0.791870  -3.969 0.000163 ***
## cmax         0.050062   0.002769  18.076  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.534 on 76 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.8113, Adjusted R-squared:  0.8088 
## F-statistic: 326.8 on 1 and 76 DF,  p-value: < 2.2e-16
lm_auc <- lm(delta_log_PAR ~ auclast, data = data_D2) 
summary(lm_auc)  
## 
## Call:
## lm(formula = delta_log_PAR ~ auclast, data = data_D2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.9209  -1.4356  -0.5528   0.7459  13.4071 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.9121046  0.6311807  -6.198 2.72e-08 ***
## auclast      0.0146749  0.0006108  24.024  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.56 on 76 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.8836, Adjusted R-squared:  0.8821 
## F-statistic: 577.2 on 1 and 76 DF,  p-value: < 2.2e-16
AIC(lm_AMT,lm_cmax,lm_auc)
##         df      AIC
## lm_AMT   3 487.0890
## lm_cmax  3 461.1277
## lm_auc   3 423.4138

Multilinear

#DOSE
mlm1_test = lm(delta_log_PAR ~ AMT + AGE + FEMALE + WT + HT, data = data_comb)
mlm_AMT = lm(delta_log_PAR ~ AMT + WT, data = data_comb)
summary(mlm1_test)
## 
## Call:
## lm(formula = delta_log_PAR ~ AMT + AGE + FEMALE + WT + HT, data = data_comb)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.0336  -2.9558  -0.1364   2.2479  15.2255 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 20.822434   9.635203   2.161 0.033375 *  
## AMT          0.027619   0.001748  15.802  < 2e-16 ***
## AGE          0.035272   0.086225   0.409 0.683472    
## FEMALE      -0.694994   1.106583  -0.628 0.531577    
## WT          -0.189334   0.054151  -3.496 0.000738 ***
## HT          -0.089463   0.066024  -1.355 0.178846    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.699 on 89 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.7757, Adjusted R-squared:  0.7631 
## F-statistic: 61.56 on 5 and 89 DF,  p-value: < 2.2e-16
summary(mlm_AMT)
## 
## Call:
## lm(formula = delta_log_PAR ~ AMT + WT, data = data_comb)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.8499 -2.6487  0.1112  2.2611 14.5479 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  9.813265   2.750187   3.568 0.000573 ***
## AMT          0.027219   0.001677  16.234  < 2e-16 ***
## WT          -0.235419   0.040616  -5.796 9.43e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.671 on 92 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.7709, Adjusted R-squared:  0.7659 
## F-statistic: 154.8 on 2 and 92 DF,  p-value: < 2.2e-16
#CMAX
mlm2_test = lm(delta_log_PAR ~ cmax + AGE + FEMALE + WT + HT, data = data_comb)
mlm_cmax = lm(delta_log_PAR ~ cmax + WT, data = data_comb)
summary(mlm2_test)
## 
## Call:
## lm(formula = delta_log_PAR ~ cmax + AGE + FEMALE + WT + HT, data = data_comb)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.8465  -1.9143  -0.7146   1.3682  17.1034 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.223534   8.713894  -0.026    0.980    
## cmax         0.050057   0.002733  18.316   <2e-16 ***
## AGE          0.011549   0.076752   0.150    0.881    
## FEMALE      -1.840488   0.985043  -1.868    0.065 .  
## WT          -0.030215   0.050501  -0.598    0.551    
## HT          -0.003129   0.058456  -0.054    0.957    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.197 on 89 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.821,  Adjusted R-squared:  0.811 
## F-statistic: 81.66 on 5 and 89 DF,  p-value: < 2.2e-16
summary(mlm_cmax)
## 
## Call:
## lm(formula = delta_log_PAR ~ cmax + WT, data = data_comb)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.8247  -1.7201  -0.6172   1.3194  18.1696 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.367666   2.765278  -0.856    0.394    
## cmax         0.050259   0.002716  18.505   <2e-16 ***
## WT          -0.013615   0.039301  -0.346    0.730    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.226 on 92 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.8125, Adjusted R-squared:  0.8084 
## F-statistic: 199.3 on 2 and 92 DF,  p-value: < 2.2e-16
#AUC
mlm3_test = lm(delta_log_PAR ~ auclast + AGE + FEMALE + WT + HT, data = data_comb)
mlm_auclast = lm(delta_log_PAR ~ auclast + WT, data = data_comb)
summary(mlm3_test)
## 
## Call:
## lm(formula = delta_log_PAR ~ auclast + AGE + FEMALE + WT + HT, 
##     data = data_comb)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.0332  -1.5049  -0.4785   1.0148  12.9415 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.2443336  6.9340319  -0.468   0.6410    
## auclast      0.0146980  0.0006083  24.161   <2e-16 ***
## AGE          0.0160783  0.0609360   0.264   0.7925    
## FEMALE      -1.3791925  0.7830259  -1.761   0.0816 .  
## WT           0.0008486  0.0403452   0.021   0.9833    
## HT          -0.0048706  0.0464298  -0.105   0.9167    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.334 on 89 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.8871, Adjusted R-squared:  0.8807 
## F-statistic: 139.8 on 5 and 89 DF,  p-value: < 2.2e-16
summary(mlm_auclast)
## 
## Call:
## lm(formula = delta_log_PAR ~ auclast + WT, data = data_comb)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.6071  -1.3818  -0.3812   0.8856  13.7214 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.9792915  2.2188323  -2.244   0.0272 *  
## auclast      0.0147574  0.0006023  24.502   <2e-16 ***
## WT           0.0120247  0.0313255   0.384   0.7020    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.347 on 92 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.8823, Adjusted R-squared:  0.8798 
## F-statistic: 344.9 on 2 and 92 DF,  p-value: < 2.2e-16
#comparison of all models 

AIC(lm_AMT,lm_cmax,lm_auc,mlm_AMT,mlm1_test,mlm_cmax,mlm2_test,mlm_auclast,mlm3_test)
## Warning in AIC.default(lm_AMT, lm_cmax, lm_auc, mlm_AMT, mlm1_test, mlm_cmax, :
## models are not all fitted to the same number of observations
##             df      AIC
## lm_AMT       3 487.0890
## lm_cmax      3 461.1277
## lm_auc       3 423.4138
## mlm_AMT      4 567.4158
## mlm1_test    7 571.3878
## mlm_cmax     4 548.3762
## mlm2_test    7 549.9427
## mlm_auclast  4 504.1082
## mlm3_test    7 506.1951
get_pred_ci <- function(fit_to_plot){
  if (class(fit_to_plot) == "nls") {
    if (!nzchar(system.file(package = "mvtnorm"))) install.packages("mvtnorm")
    for (i in 1:500) {
      if (i == 1) nls_pred <- c()
      params <- drop(mvtnorm::rmvnorm(1,coef(fit_to_plot),
                                      sigma = vcov(fit_to_plot)))
      nls_pred <- cbind(nls_pred,predict(fit_to_plot,params))
    }
    nls_pred <- t(apply(nls_pred,1,quantile,probs = c(0.025, 0.975)))
    nls_pred <- as.data.frame(cbind(fit = predict(fit_to_plot),nls_pred)) 
    mod_pred <- nls_pred
    names(mod_pred) <- c("fit","lwr","upr")
  } else if (class(fit_to_plot) == "lm") {
    mod_pred <- data.frame(predict(lm_fit, interval = "confidence"))
  }
  return(mod_pred)
}

None linear

data_comb %<>%
  group_by(ID) %>%
  na.omit(delta_log_PAR)

### Non-linear

##Dose

#nlm_AMT <- nls(delta_log_PAR ~ E0 + (EMAX * AMT) / (ED50 + AMT), data = data_comb, start = list(E0 = min(data_comb$delta_log_PAR, na.rm = TRUE), EMAX = max(data_comb$delta_log_PAR, na.rm = TRUE),ED50 = 500), trace = TRUE)
#summary(nlm_AMT)


nlm_AMT_GAMMA <- nls(delta_log_PAR ~ E0 + (EMAX * AMT^GAMMA) / (ED50^GAMMA + AMT^GAMMA),
               data = data_comb,
               start = list(E0 = min(data_comb$delta_log_PAR, na.rm = TRUE),
                            EMAX = max(data_comb$delta_log_PAR, na.rm = TRUE),
                            ED50 = 1000, 
                            GAMMA = 2),
               algorithm= "port" ,
               trace = TRUE)
##   0:     2663.1734: -3.67594  33.9593  1000.00  2.00000
##   1:     1359.9953: -2.15533  40.5630  846.604  1.89827
##   2:     1336.8143: -1.46793  38.8960  806.545  2.44619
##   3:     1330.3561: -1.71462  32.9168  708.475  2.64572
##   4:     1327.5171: -1.64536  31.3328  683.649  2.86339
##   5:     1324.1098: -1.46470  27.7847  633.435  3.24159
##   6:     1319.3641: -1.34237  24.9678  594.834  3.89511
##   7:     1314.8378: -1.26195  21.9780  557.175  4.50493
##   8:     1309.3286: -1.19746  20.5236  540.389  5.32937
##   9:     1305.4864: -1.13725  18.8616  523.141  6.74302
##  10:     1302.0435: -1.10881  17.8756  514.762  9.12361
##  11:     1297.5749: -1.11562  18.0377  515.740  10.6775
##  12:     1297.4863: -1.11604  18.0634  515.274  10.9762
##  13:     1297.4861: -1.11602  18.0631  515.259  10.9963
##  14:     1297.4861: -1.11602  18.0630  515.258  10.9969
summary(nlm_AMT_GAMMA)
## 
## Formula: delta_log_PAR ~ E0 + (EMAX * AMT^GAMMA)/(ED50^GAMMA + AMT^GAMMA)
## 
## Parameters:
##       Estimate Std. Error t value Pr(>|t|)    
## E0     -1.1160     0.8552  -1.305    0.195    
## EMAX   18.0630     1.9702   9.168 1.41e-14 ***
## ED50  515.2579    28.4275  18.125  < 2e-16 ***
## GAMMA  10.9969    12.8410   0.856    0.394    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.34 on 91 degrees of freedom
## 
## Algorithm "port", convergence message: relative convergence (4)
##Cmax

nlm_cmax <- nls(delta_log_PAR ~ E0 + (EMAX * cmax) / (ED50 + cmax), data = data_comb,
                start = list(E0 = min(data_comb$delta_log_PAR, na.rm = TRUE), 
                             EMAX = max(data_comb$delta_log_PAR, na.rm = TRUE),
                             ED50 = 1000), 
                              trace = TRUE)
## 7042.434    (2.19e+00): par = (-3.675938 33.95934 1000)
## 4314.954    (1.64e+00): par = (-5.593921 42.83437 209.0485)
## 2254.314    (9.23e-01): par = (-3.60822 53.41605 490.1278)
## 1222.140    (3.73e-02): par = (-7.131415 65.67036 646.2288)
## 1220.526    (5.53e-03): par = (-7.423884 65.03921 616.3817)
## 1220.495    (1.03e-03): par = (-7.390536 65.43856 624.5661)
## 1220.494    (2.50e-04): par = (-7.400444 65.3414 622.6502)
## 1220.494    (5.99e-05): par = (-7.398178 65.36499 623.1142)
## 1220.494    (1.44e-05): par = (-7.39873 65.35933 623.0028)
## 1220.494    (3.46e-06): par = (-7.398598 65.36069 623.0295)
summary(nlm_cmax)
## 
## Formula: delta_log_PAR ~ E0 + (EMAX * cmax)/(ED50 + cmax)
## 
## Parameters:
##      Estimate Std. Error t value Pr(>|t|)    
## E0     -7.399      1.203  -6.148 1.99e-08 ***
## EMAX   65.361     10.237   6.385 6.88e-09 ***
## ED50  623.030    194.104   3.210  0.00183 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.642 on 92 degrees of freedom
## 
## Number of iterations to convergence: 9 
## Achieved convergence tolerance: 3.462e-06
nlm_cmax_GAMMA <- nls(delta_log_PAR ~ E0 + (EMAX * cmax^GAMMA) / (ED50^GAMMA + cmax^GAMMA),
               data = data_comb,
               start = list(E0 = min(data_comb$delta_log_PAR, na.rm = TRUE),
                            EMAX = max(data_comb$delta_log_PAR, na.rm = TRUE),
                            ED50 = 1000, 
                            GAMMA = 2),
               algorithm = "port" ,
               trace = TRUE)
##   0:     6615.0972: -3.67594  33.9593  1000.00  2.00000
##   1:     2969.0238: -1.06033  35.7427  932.706  1.42406
##   2:     888.35234: -6.44603  27.1111  157.890  1.56444
##   3:     659.85163: -6.98374  37.8367  249.412  1.27171
##   4:     640.30574: -6.93451  41.2176  270.195  1.26588
##   5:     609.23233: -6.67062  41.0647  285.730  1.33110
##   6:     573.68407: -3.33017  36.7639  315.590  1.66749
##   7:     555.30936: -3.13558  33.9460  285.753  1.94549
##   8:     550.86367: -2.32040  28.9239  254.383  2.37670
##   9:     545.48804: -2.33508  29.4246  261.490  2.49025
##  10:     545.40435: -2.29358  29.1440  259.561  2.53746
##  11:     545.40243: -2.28750  29.1276  259.503  2.54121
##  12:     545.39998: -2.27901  29.0946  259.441  2.54830
##  13:     545.39986: -2.27657  29.0847  259.402  2.55012
##  14:     545.39985: -2.27598  29.0824  259.394  2.55055
##  15:     545.39985: -2.27584  29.0818  259.392  2.55065
##  16:     545.39985: -2.27580  29.0817  259.392  2.55067
summary(nlm_cmax_GAMMA)
## 
## Formula: delta_log_PAR ~ E0 + (EMAX * cmax^GAMMA)/(ED50^GAMMA + cmax^GAMMA)
## 
## Parameters:
##       Estimate Std. Error t value Pr(>|t|)    
## E0     -2.2758     0.8399  -2.709  0.00805 ** 
## EMAX   29.0817     3.1757   9.158 1.48e-14 ***
## ED50  259.3919    22.7239  11.415  < 2e-16 ***
## GAMMA   2.5507     0.4705   5.421 4.83e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.462 on 91 degrees of freedom
## 
## Algorithm "port", convergence message: relative convergence (4)
## AUC

nlm_auclast <- nls(delta_log_PAR ~ E0 + (EMAX * auclast) / (ED50 + auclast), data = data_comb,
                start = list(E0 = min(data_comb$delta_log_PAR, na.rm = TRUE), 
                             EMAX = max(data_comb$delta_log_PAR, na.rm = TRUE),
                             ED50 = 1000), 
                              trace = TRUE)
## 2198.824    (1.46e+00): par = (-3.675938 33.95934 1000)
## 1405.527    (9.12e-01): par = (-2.620688 66.79992 3248.3)
## 798.1182    (1.93e-01): par = (-6.608581 96.84047 4441.392)
## 770.0040    (2.25e-02): par = (-6.825371 94.58061 4045.365)
## 769.6409    (2.77e-03): par = (-6.779813 96.33516 4175.566)
## 769.6363    (6.94e-04): par = (-6.796859 95.86418 4144.375)
## 769.6361    (1.78e-04): par = (-6.792978 95.98655 4152.534)
## 769.6360    (4.59e-05): par = (-6.794006 95.95512 4150.437)
## 769.6360    (1.18e-05): par = (-6.793743 95.96323 4150.978)
## 769.6360    (3.06e-06): par = (-6.79381 95.96114 4150.839)
summary(nlm_auclast)
## 
## Formula: delta_log_PAR ~ E0 + (EMAX * auclast)/(ED50 + auclast)
## 
## Parameters:
##       Estimate Std. Error t value Pr(>|t|)    
## E0     -6.7938     0.8129  -8.357 6.48e-13 ***
## EMAX   95.9611    18.6927   5.134 1.58e-06 ***
## ED50 4150.8390  1229.6321   3.376  0.00108 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.892 on 92 degrees of freedom
## 
## Number of iterations to convergence: 9 
## Achieved convergence tolerance: 3.059e-06
nlm_auclast_GAMMA <- nls(delta_log_PAR ~ E0 + (EMAX * auclast^GAMMA) / (ED50^GAMMA + auclast^GAMMA),
               data = data_comb,
               start = list(E0 = min(data_comb$delta_log_PAR, na.rm = TRUE),
                            EMAX = max(data_comb$delta_log_PAR, na.rm = TRUE),
                            ED50 = 1000, 
                            GAMMA = 2),
               algorithm = "port" ,
               trace = TRUE)
##   0:     398.24279: -3.67594  33.9593  1000.00  2.00000
##   1:     369.31703: -3.45006  33.9331  1051.49  2.10049
##   2:     345.51968: -2.58271  34.6090  1151.63  2.22651
##   3:     338.65860: -2.52060  37.3512  1253.01  2.13641
##   4:     337.25954: -2.72211  40.0857  1333.16  2.00672
##   5:     337.19394: -2.75527  40.4317  1343.93  2.00072
##   6:     337.19393: -2.75577  40.4283  1343.75  2.00081
##   7:     337.19393: -2.75577  40.4284  1343.75  2.00081
summary(nlm_auclast_GAMMA)
## 
## Formula: delta_log_PAR ~ E0 + (EMAX * auclast^GAMMA)/(ED50^GAMMA + auclast^GAMMA)
## 
## Parameters:
##        Estimate Std. Error t value Pr(>|t|)    
## E0      -2.7558     0.7563  -3.644 0.000447 ***
## EMAX    40.4284     5.4168   7.464 4.87e-11 ***
## ED50  1343.7548   170.9371   7.861 7.42e-12 ***
## GAMMA    2.0008     0.2829   7.073 3.03e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.722 on 91 degrees of freedom
## 
## Algorithm "port", convergence message: relative convergence (4)
AIC(nlm_cmax,nlm_auclast,nlm_AMT_GAMMA,nlm_cmax_GAMMA,nlm_auclast_GAMMA)
##                   df      AIC
## nlm_cmax           4 520.1460
## nlm_auclast        4 476.3422
## nlm_AMT_GAMMA      5 593.8065
## nlm_cmax_GAMMA     5 511.4733
## nlm_auclast_GAMMA  5 465.7915
ART_500 <- data_comb %>%
  filter(DRUG_cat == "ART-500")
t.test(ART_500$delta_log_PAR, conf.level = 0.95)
## 
##  One Sample t-test
## 
## data:  ART_500$delta_log_PAR
## t = 4.9577, df = 16, p-value = 0.0001425
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  3.683843 9.187715
## sample estimates:
## mean of x 
##  6.435779
plot1 <- ggplot(data = data_comb, aes(x = auclast, y = delta_log_PAR)) +
  geom_point() + geom_ribbon(aes(ymin = 3.68 , ymax = 9.18), 
              fill = "lightgreen", alpha = 0.5) + 
  geom_hline(yintercept = 6.435779, color = "red", alpha = 0.5)


plot1

plot1 + geom_smooth(method = "lm") + theme_classic()
## `geom_smooth()` using formula = 'y ~ x'

fit_to_plot = nlm_auclast
mod_pred_ci = get_pred_ci(fit_to_plot)

plot1 +   geom_ribbon(aes(ymin = mod_pred_ci$lwr, ymax = mod_pred_ci$upr), 
              fill = "lightblue", alpha = 0.5) +
  geom_line(aes(y = mod_pred_ci$fit), col = "blue", lwd = 1)

plot2 <- ggplot(data = data_comb, aes(x = AMT, y = auclast)) +  geom_point()
plot2 + geom_smooth(method = "lm") + theme_classic() + 
  geom_hline(yintercept = 1100, color = "red", alpha = 0.5)
## `geom_smooth()` using formula = 'y ~ x'

nls_fit_to_plot <- nlm_auclast_GAMMA
if (!nzchar(system.file(package = "mvtnorm"))) install.packages("mvtnorm")
for (i in 1:500) {
  if (i == 1) nls_pred <- c()
  params <- drop(mvtnorm::rmvnorm(1,coef(nls_fit_to_plot),
    sigma = vcov(nls_fit_to_plot)))
  nls_pred <- cbind(nls_pred,predict(nls_fit_to_plot,params))
}
nls_pred <- t(apply(nls_pred,1,quantile,probs = c(0.025, 0.975)))
nls_pred <- as.data.frame(cbind(fit = predict(nls_fit_to_plot),nls_pred))
plot1 +
geom_ribbon(aes(ymin = nls_pred$`2.5%`,ymax = nls_pred$`97.5%`), fill = "lightblue",alpha = 0.5) +
geom_line(aes(y = fitted(nls_fit_to_plot)),col = "blue",lwd = 1)

fit.emax <- stan_emax(delta_log_PAR ~ auclast, data = data_comb,
                      # the next line is only to make the output short
                      chains = 1, iter = 500, seed = 12345)
## 
## SAMPLING FOR MODEL 'emax' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 6.2e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.62 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:   1 / 500 [  0%]  (Warmup)
## Chain 1: Iteration:  50 / 500 [ 10%]  (Warmup)
## Chain 1: Iteration: 100 / 500 [ 20%]  (Warmup)
## Chain 1: Iteration: 150 / 500 [ 30%]  (Warmup)
## Chain 1: Iteration: 200 / 500 [ 40%]  (Warmup)
## Chain 1: Iteration: 250 / 500 [ 50%]  (Warmup)
## Chain 1: Iteration: 251 / 500 [ 50%]  (Sampling)
## Chain 1: Iteration: 300 / 500 [ 60%]  (Sampling)
## Chain 1: Iteration: 350 / 500 [ 70%]  (Sampling)
## Chain 1: Iteration: 400 / 500 [ 80%]  (Sampling)
## Chain 1: Iteration: 450 / 500 [ 90%]  (Sampling)
## Chain 1: Iteration: 500 / 500 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.347 seconds (Warm-up)
## Chain 1:                0.19 seconds (Sampling)
## Chain 1:                0.537 seconds (Total)
## Chain 1:
## Warning: The largest R-hat is NA, indicating chains have not mixed.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#r-hat
## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess
fit.emax
## ---- Emax model fit with rstanemax ----
## 
##          mean se_mean     sd    2.5%     25%     50%     75%   97.5%  n_eff
## emax    80.87    0.91   8.41   66.01   74.70   80.44   86.80   96.29  86.16
## e0      -7.22    0.09   0.73   -8.69   -7.65   -7.23   -6.74   -5.86  65.96
## ec50  3210.76   62.96 525.94 2203.35 2843.97 3191.13 3632.56 4084.57  69.78
## gamma    1.00     NaN   0.00    1.00    1.00    1.00    1.00    1.00    NaN
## sigma    2.96    0.02   0.21    2.58    2.82    2.95    3.08    3.44 116.87
##       Rhat
## emax  1.01
## e0    1.00
## ec50  1.01
## gamma  NaN
## sigma 1.00
## 
## * Use `extract_stanfit()` function to extract raw stanfit object
## * Use `extract_param()` function to extract posterior draws of key parameters
## * Use `plot()` function to visualize model fit
## * Use `posterior_predict()` or `posterior_predict_quantile()` function to get
##   raw predictions or make predictions on new data
## * Use `extract_obs_mod_frame()` function to extract raw data 
##   in a processed format (useful for plotting)
plot(fit.emax)

Are there any indications that the dose you chose will lead to side effects at a greater rate than those seen for artemisinin?

library(binom)
data_ae = data_comb %>% group_by(DRUG_cat) %>% 
  summarise(AE_tot = sum(side_effect),
            binom.confint(x = AE_tot, n = n(), methods = "wilson"),
            .groups = "drop_last")

ggplot(data = data_ae, aes(x = DRUG_cat, y = mean)) +
  geom_point(shape = 18, size = 7) +
  geom_errorbar(aes(x = DRUG_cat, y = mean, ymin = lower, ymax = upper), width = .3) +
  ylab("Proportion with side effect") +
  xlab("DRUG-Dose") +theme_test()