Better alternatives than normalizing to control using R

Better alternatives than normalizing to control using R Version 1 Christian Ritz ([email protected]) Last compiled: 2025-10-30 Introduction This document provides all Rcode and output used for the paper entitled “Better alternatives than normalizing to control: Case studies with algae toxicity and dose-response analysis”. There are three parts. The first part is a re-analysis of the data by Andersen et al. (1998). Specifically, data are imported from the accompanying file “andersen_etal1998-appendixB.csv”; data are from 4 experiments, but only data from experiment 1 and experiment 2 are used although similar dose-response analyses could be carried out for the two other experiments. Then dose-response models are fitted to the data using the three approaches described in detail in the paper: the normalization, correlation, and growth rate-based (raw data) approaches (using the Rpackage “drc”). Model fits are summarized and estimated ED10, ED20, and ED50 are reported with delta method-based confidence intervals. The fitted dose-response curves are also shown (plots obtained using the Rpackage “ggplot2”). In the second part simulations are carried out, based on datasets generated using the parameter estimates obtained from fitting dose-response models to the data from experiment 1 and experiment 2. For a number of scenarios with varying numbers of measurements for the control group 1000 datasets are obtained and bias, coefficient of variation (CV), coverage and length of confidence intervals are calculated for ED10, ED20, and ED50. Results are provided as LaTeX tables. The final part produces a figure showing the bias in the fitted dose-response curve (and also in the ED50) for the normalization approach compared to the growth rate-based approach, using simulated data to capture an average discrepancy between the two fitted curves. Re-analysis of data from Andersen et al. R packages used Loading relevant extension packages: library(drc) library(ggplot2) library(plyr) library(xtable) 1 Information on the Rsession: sessionInfo() ## R version 4.5.1 (2025-06-13) ## Platform: x86_64-pc-linux-gnu ## Running under: Ubuntu 24.04.3 LTS ## ## Matrix products: default ## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.12.0 ## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0 ## ## locale: ## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C ## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 ## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 ## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C ## [9] LC_ADDRESS=C LC_TELEPHONE=C ## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C ## ## time zone: Europe/Copenhagen ## tzcode source: system (glibc) ## ## attached base packages: ## [1] stats graphics grDevices utils datasets methods base ## ## other attached packages: ## [1] xtable_1.8-4 plyr_1.8.9 ggplot2_3.5.2 drc_3.2-0 drcData_1.1-3 ## [6] MASS_7.3-65 ## ## loaded via a namespace (and not attached): ## [1] Matrix_1.7-3 gtable_0.3.6 dplyr_1.1.4 compiler_4.5.1 ## [5] gtools_3.9.5 Rcpp_1.0.14 tidyselect_1.2.1 plotrix_3.8-4 ## [9] splines_4.5.1 scales_1.4.0 yaml_2.3.10 fastmap_1.2.0 ## [13] lattice_0.22-5 TH.data_1.1-3 R6_2.6.1 generics_0.1.4 ## [17] Formula_1.2-5 knitr_1.50 tibble_3.2.1 car_3.1-3 ## [21] pillar_1.10.2 RColorBrewer_1.1-3 rlang_1.1.6 multcomp_1.4-28 ## [25] xfun_0.52 cli_3.6.5 withr_3.0.2 magrittr_2.0.3 ## [29] digest_0.6.37 grid_4.5.1 mvtnorm_1.3-3 sandwich_3.1-1 ## [33] lifecycle_1.0.4 vctrs_0.6.5 evaluate_1.0.3 glue_1.8.0 ## [37] farver_2.1.2 codetools_0.2-20 zoo_1.8-14 survival_3.8-3 ## [41] abind_1.4-8 carData_3.0-5 rmarkdown_2.29 pkgconfig_2.0.3 ## [45] tools_4.5.1 htmltools_0.5.8.1 Data import and formatting Data from Andersen et al. (1998) consist of results from 4 separate dose-response experiments. In all experiments there are 6 replicates for the control (concentration 0) and for 2 of the experiments there are 3 replicates for each non-zero concentration (experiment 1 and 3). For the 2 other experiments concentrations were measured and there are no replicates. For each growth rate, which is an estimated slope, the corresponding standard error is also provided and there is a trend: the smaller the growth rate the smaller is the corresponding standard error. Note that the reciprocal standard errors will be weights in the nonlinear least squares estimation. ## Conc1 GR1 SE1 Conc2 GR2 SE2 Conc3 2 ## 1 0.000 2.3842697 0.022789 0.0000 2.014286855 0.002339547 0.0000 ## 2 0.000 2.3816736 0.045754 0.0000 1.777592604 0.021759548 0.0000 ## 3 0.000 2.3295015 0.020538 0.0000 2.007431568 0.003522097 0.0000 ## 4 0.000 2.3238405 0.019038 0.0000 1.974795254 0.008190391 0.0000 ## 5 0.000 2.4361314 0.023957 0.0000 2.181544605 0.007758268 0.0000 ## 6 0.000 2.2403203 0.059913 0.0000 2.021203652 0.006898431 0.0000 ## 7 5.625 2.5607525 0.030831 2.8937 2.003873121 0.005528977 6.1295 ## 8 5.625 2.6528099 0.039230 5.2087 1.900418733 0.002117418 6.1295 ## 9 5.625 2.2725164 0.057254 9.3757 1.940659369 0.005787856 6.1295 ## 10 11.250 2.1843847 0.061625 16.8760 2.132061593 0.001597769 12.2590 ## 11 11.250 2.2830658 0.024091 30.3770 2.034806770 0.005771680 12.2590 ## 12 11.250 2.3271210 0.013814 54.6790 1.845920248 0.002535440 12.2590 ## 13 22.500 1.9531553 0.058811 98.4220 1.506089761 0.006155441 24.5180 ## 14 22.500 1.9671423 0.025465 177.1600 0.748736268 0.003714306 24.5180 ## 15 22.500 2.0739440 0.015759 318.8800 0.002244535 0.003714306 24.5180 ## 16 45.000 1.7810941 0.033720 574.0000 0.001423340 0.003714306 49.0360 ## 17 45.000 1.6997967 0.037853 1033.2000 0.003312330 0.003714306 49.0360 ## 18 45.000 1.6374060 0.008258 NA NA NA 49.0360 ## 19 90.000 0.9344162 0.107094 NA NA NA 98.0720 ## 20 90.000 0.8480270 0.200411 NA NA NA 98.0720 ## 21 90.000 1.0391678 0.181461 NA NA NA 98.0720 ## 22 NA NA NA NA NA NA 196.1400 ## 23 NA NA NA NA NA NA 196.1400 ## 24 NA NA NA NA NA NA 196.1400 ## GR3 SE3 Conc4 GR4 SE4 ## 1 0.094201605 0.005671437 0.0 1.62862155 1.93982e-04 ## 2 0.095821757 0.003082070 0.0 1.70842215 3.84363e-04 ## 3 0.101719661 0.003377968 0.0 1.67516097 1.27243e-04 ## 4 0.097677733 0.002735039 0.0 1.62390305 1.37185e-04 ## 5 0.092875246 0.003953754 0.0 1.57240369 1.59674e-04 ## 6 0.091549900 0.003794118 0.0 1.61266088 1.60233e-04 ## 7 0.092682595 0.003190673 1.0 1.68462686 1.21073e-04 ## 8 0.089956861 0.003725896 1.0 1.63151726 7.90277e-05 ## 9 0.088366498 0.003953586 1.0 1.62652671 1.09256e-04 ## 10 0.100161797 0.001930250 13.3 1.46795912 2.55040e-04 ## 11 0.101425667 0.002236926 13.3 1.50856517 1.15686e-04 ## 12 0.098504437 0.002148053 13.3 1.33799425 5.69584e-05 ## 13 0.087521470 0.005298999 66.5 1.10908289 1.52591e-04 ## 14 0.087939234 0.004450987 66.5 1.04919619 9.48687e-05 ## 15 0.084547475 0.003122399 66.5 1.02152777 1.43305e-04 ## 16 0.065590789 0.001887717 332.5 0.10254159 3.28521e-04 ## 17 0.072841333 0.003088261 332.5 0.05138762 1.37442e-04 ## 18 0.071835786 0.002707159 332.5 -0.00804090 1.53620e-04 ## 19 0.031687968 0.006872033 665.0 0.04640898 5.49214e-04 ## 20 0.035072356 0.005021663 665.0 -0.05319582 8.41774e-05 ## 21 0.035348704 0.006886412 665.0 -0.05339642 9.02568e-05 ## 22 0.001044420 0.005337032 NA NA NA ## 23 0.001370118 0.005372703 NA NA NA ## 24 0.001040508 0.006152045 NA NA NA andersen.data <- read.csv2("andersen_etal1998-appendixB.csv") names(andersen.data) ## [1] "Conc1" "GR1" "SE1" "Conc2" "GR2" "SE2" "Conc3" "GR3" "SE3" 3 ## [10] "Conc4" "GR4" "SE4" summary(andersen.data) ## Conc1 GR1 SE1 Conc2 ## Min. : 0.00 Min. :0.848 Min. :0.008258 Min. : 0.000 ## 1st Qu.: 0.00 1st Qu.:1.781 1st Qu.:0.022789 1st Qu.: 0.000 ## Median :11.25 Median :2.240 Median :0.033720 Median : 9.376 ## Mean :24.91 Mean :2.015 Mean :0.051794 Mean : 136.534 ## 3rd Qu.:45.00 3rd Qu.:2.330 3rd Qu.:0.058811 3rd Qu.: 98.422 ## Max. :90.00 Max. :2.653 Max. :0.200411 Max. :1033.200 ## NA's :3 NA's :3 NA's :3 NA's :7 ## GR2 SE2 Conc3 GR3 ## Min. :0.001423 Min. :0.001598 Min. : 0.000 Min. :0.001041 ## 1st Qu.:1.506090 1st Qu.:0.003522 1st Qu.: 4.597 1st Qu.:0.058030 ## Median :1.940659 Median :0.003714 Median : 18.389 Median :0.088153 ## Mean :1.535082 Mean :0.005578 Mean : 48.269 Mean :0.071699 ## 3rd Qu.:2.014287 3rd Qu.:0.006155 3rd Qu.: 61.295 3rd Qu.:0.094607 ## Max. :2.181545 Max. :0.021760 Max. :196.140 Max. :0.101720 ## NA's :7 NA's :7 ## SE3 Conc4 GR4 SE4 ## Min. :0.001888 Min. : 0.0 Min. :-0.0534 Min. :5.696e-05 ## 1st Qu.:0.002995 1st Qu.: 0.0 1st Qu.: 0.1025 1st Qu.:1.093e-04 ## Median :0.003760 Median : 13.3 Median : 1.4680 Median :1.374e-04 ## Mean :0.004000 Mean :154.0 Mean : 1.0640 Mean :1.730e-04 ## 3rd Qu.:0.005309 3rd Qu.:332.5 3rd Qu.: 1.6265 3rd Qu.:1.602e-04 ## Max. :0.006886 Max. :665.0 Max. : 1.7084 Max. :5.492e-04 ## NA's :3 NA's :3 NA's :3 andersen.data ## Conc1 GR1 SE1 Conc2 GR2 SE2 Conc3 ## 1 0.000 2.3842697 0.022789 0.0000 2.014286855 0.002339547 0.0000 ## 2 0.000 2.3816736 0.045754 0.0000 1.777592604 0.021759548 0.0000 ## 3 0.000 2.3295015 0.020538 0.0000 2.007431568 0.003522097 0.0000 ## 4 0.000 2.3238405 0.019038 0.0000 1.974795254 0.008190391 0.0000 ## 5 0.000 2.4361314 0.023957 0.0000 2.181544605 0.007758268 0.0000 ## 6 0.000 2.2403203 0.059913 0.0000 2.021203652 0.006898431 0.0000 ## 7 5.625 2.5607525 0.030831 2.8937 2.003873121 0.005528977 6.1295 ## 8 5.625 2.6528099 0.039230 5.2087 1.900418733 0.002117418 6.1295 ## 9 5.625 2.2725164 0.057254 9.3757 1.940659369 0.005787856 6.1295 ## 10 11.250 2.1843847 0.061625 16.8760 2.132061593 0.001597769 12.2590 ## 11 11.250 2.2830658 0.024091 30.3770 2.034806770 0.005771680 12.2590 ## 12 11.250 2.3271210 0.013814 54.6790 1.845920248 0.002535440 12.2590 ## 13 22.500 1.9531553 0.058811 98.4220 1.506089761 0.006155441 24.5180 ## 14 22.500 1.9671423 0.025465 177.1600 0.748736268 0.003714306 24.5180 ## 15 22.500 2.0739440 0.015759 318.8800 0.002244535 0.003714306 24.5180 ## 16 45.000 1.7810941 0.033720 574.0000 0.001423340 0.003714306 49.0360 ## 17 45.000 1.6997967 0.037853 1033.2000 0.003312330 0.003714306 49.0360 ## 18 45.000 1.6374060 0.008258 NA NA NA 49.0360 ## 19 90.000 0.9344162 0.107094 NA NA NA 98.0720 ## 20 90.000 0.8480270 0.200411 NA NA NA 98.0720 ## 21 90.000 1.0391678 0.181461 NA NA NA 98.0720 ## 22 NA NA NA NA NA NA 196.1400 ## 23 NA NA NA NA NA NA 196.1400 ## 24 NA NA NA NA NA NA 196.1400 4 ## GR3 SE3 Conc4 GR4 SE4 ## 1 0.094201605 0.005671437 0.0 1.62862155 1.93982e-04 ## 2 0.095821757 0.003082070 0.0 1.70842215 3.84363e-04 ## 3 0.101719661 0.003377968 0.0 1.67516097 1.27243e-04 ## 4 0.097677733 0.002735039 0.0 1.62390305 1.37185e-04 ## 5 0.092875246 0.003953754 0.0 1.57240369 1.59674e-04 ## 6 0.091549900 0.003794118 0.0 1.61266088 1.60233e-04 ## 7 0.092682595 0.003190673 1.0 1.68462686 1.21073e-04 ## 8 0.089956861 0.003725896 1.0 1.63151726 7.90277e-05 ## 9 0.088366498 0.003953586 1.0 1.62652671 1.09256e-04 ## 10 0.100161797 0.001930250 13.3 1.46795912 2.55040e-04 ## 11 0.101425667 0.002236926 13.3 1.50856517 1.15686e-04 ## 12 0.098504437 0.002148053 13.3 1.33799425 5.69584e-05 ## 13 0.087521470 0.005298999 66.5 1.10908289 1.52591e-04 ## 14 0.087939234 0.004450987 66.5 1.04919619 9.48687e-05 ## 15 0.084547475 0.003122399 66.5 1.02152777 1.43305e-04 ## 16 0.065590789 0.001887717 332.5 0.10254159 3.28521e-04 ## 17 0.072841333 0.003088261 332.5 0.05138762 1.37442e-04 ## 18 0.071835786 0.002707159 332.5 -0.00804090 1.53620e-04 ## 19 0.031687968 0.006872033 665.0 0.04640898 5.49214e-04 ## 20 0.035072356 0.005021663 665.0 -0.05319582 8.41774e-05 ## 21 0.035348704 0.006886412 665.0 -0.05339642 9.02568e-05 ## 22 0.001044420 0.005337032 NA NA NA ## 23 0.001370118 0.005372703 NA NA NA ## 24 0.001040508 0.006152045 NA NA NA To be able to reproduce results from Andersen et al. (1998) growth inhibition is calculated and included in the dataset using the averages of the growth rates in the control groups, i.e., growth inhibition data are obtained using normalization of the original growth rates provided in the data. For each experiment a new variable is added to the dataset. Initially averages of the control groups are calculated (they will also be used later). control1 <- with(andersen.data, mean(GR1[1:6])) control2 <- with(andersen.data, mean(GR2[1:6])) control3 <- with(andersen.data, mean(GR3[1:6])) control4 <- with(andersen.data, mean(GR4[1:6])) andersen.data[["gi1"]] <- with(andersen.data, 1-GR1 /control1) andersen.data[["gi2"]] <- with(andersen.data, 1-GR2 /control2) andersen.data[["gi3"]] <- with(andersen.data, 1-GR3 /control3) andersen.data[["gi4"]] <- with(andersen.data, 1-GR4 /control4) By normalizing with the averages in the control groups, the resulting inhibition growth rates become correlated with each experiment as they were all normalized using the same average. Fitting naive models - the normalization approach Naive models ignore the correlation introduced through normalization, i.e., the growth inhibition data are treated as independent data. By construction growth inhibition data mostly be in the range 0 to 1 and therefore two-parameter doseresponse models assuming lower and upper limits equal to 0 and 1, respectively, are commonly used. Following Andersen et al. two-parameter Weibull type 2 models are fitted (possibly reflecting some asymmetry in the dose-response trend). The standard errors are incorporated as weights; for experiment 1 standard errors are increased with concentration whereas for the other experiments no clear trends are visible. 5 Initially, models are fitted without weights, using the function drm() from the package “drc”: gi1.W2.2nw <- drm(gi1 ~Conc1, data = andersen.data, fct = W2.2()) gi2.W2.2nw <- drm(gi2 ~Conc2, data = andersen.data, fct = W2.2()) Parameter estimated and estimated ED10, ED20, and ED50 are obtained using the summary() and ED() methods (also from the package “drc”): summary(gi1.W2.2nw) ## ## Model fitted: Weibull (type 2) with lower limit at 0 and upper limit at 1 (2 parms) ## ## Parameter estimates: ## ## Estimate Std. Error t-value p-value ## b:(Intercept) 1.46636 0.13826 10.606 2.023e-09 *** ## e:(Intercept) 95.44528 5.29151 18.037 2.064e-13 *** ## --- ## Signif. codes: 0 '***'0.001 '**'0.01 '*'0.05 '.'0.1 ' ' 1 ## ## Residual standard error: ## ## 0.05107554 (19 degrees of freedom) ED(gi1.W2.2nw, c(10,20,50), interval = "delta") ## ## Estimated effective doses ## ## Estimate Std. Error Lower Upper ## e:1:10 20.5712 2.4732 15.3946 25.7477 ## e:1:20 34.3174 2.6595 28.7509 39.8839 ## e:1:50 74.3368 3.3843 67.2533 81.4203 summary(gi2.W2.2nw) ## ## Model fitted: Weibull (type 2) with lower limit at 0 and upper limit at 1 (2 parms) ## ## Parameter estimates: ## ## Estimate Std. Error t-value p-value ## b:(Intercept) 2.30405 0.30192 7.6313 1.532e-06 *** ## e:(Intercept) 176.22304 8.49450 20.7455 1.855e-12 *** ## --- ## Signif. codes: 0 '***'0.001 '**'0.01 '*'0.05 '.'0.1 ' ' 1 ## ## Residual standard error: ## ## 0.04578018 (15 degrees of freedom) ED(gi2.W2.2nw, c(10,20,50), interval = "delta") ## ## Estimated effective doses ## ## Estimate Std. Error Lower Upper 6 ## e:1:10 66.3570 8.0854 49.1233 83.5907 ## e:1:20 91.9044 7.7057 75.4801 108.3287 ## e:1:50 150.3066 6.9367 135.5213 165.0919 The resulting model fits may be used to obtain starting values if needed for the weighted model fits below. Now models including the weights are fitted: gi1.W2.2 <- drm(gi1 ~Conc1, data = andersen.data, fct = W2.2(), weights = SE1 /control1, start = coef(gi1.W2.2nw)) # provision of starting values needed to avoid sub-optimal fit gi2.W2.2 <- drm(gi2 ~Conc2, data = andersen.data, fct = W2.2(), weights = SE2 /control2) The normalized standard errors are provided as weights, which in the function drm() are understood in this way: the smaller the value the more weight is put of the corresponding growth inhibition value in the least squares estimation. The resulting parameter estimates and estimated ED values are: summary(gi1.W2.2) ## ## Model fitted: Weibull (type 2) with lower limit at 0 and upper limit at 1 (2 parms) ## ## Parameter estimates: ## ## Estimate Std. Error t-value p-value ## b:(Intercept) 1.62894 0.17929 9.0857 2.409e-08 *** ## e:(Intercept) 84.89321 6.74838 12.5798 1.163e-10 *** ## --- ## Signif. codes: 0 '***'0.001 '**'0.01 '*'0.05 '.'0.1 ' ' 1 ## ## Residual standard error: ## ## 3.506035 (19 degrees of freedom) ED(gi1.W2.2,c(10,20,50), interval = "delta") ## ## Estimated effective doses ## ## Estimate Std. Error Lower Upper ## e:1:10 21.3254 1.7557 17.6507 25.0002 ## e:1:20 33.8042 1.3032 31.0765 36.5319 ## e:1:50 67.7885 3.8596 59.7103 75.8667 summary(gi2.W2.2) ## ## Model fitted: Weibull (type 2) with lower limit at 0 and upper limit at 1 (2 parms) ## ## Parameter estimates: ## ## Estimate Std. Error t-value p-value ## b:(Intercept) 2.42276 0.36571 6.6247 8.085e-06 *** ## e:(Intercept) 176.58826 9.52198 18.5453 9.396e-12 *** ## --- 7 ## Signif. codes: 0 '***'0.001 '**'0.01 '*'0.05 '.'0.1 ' ' 1 ## ## Residual standard error: ## ## 27.5425 (15 degrees of freedom) ED(gi2.W2.2,c(10,20,50), interval = "delta") ## ## Estimated effective doses ## ## Estimate Std. Error Lower Upper ## e:1:10 69.7542 9.7769 48.9153 90.5930 ## e:1:20 95.0799 9.3616 75.1262 115.0336 ## e:1:50 151.7967 8.2493 134.2136 169.3797 Compiling estimated EC values in to a tabel Based on the model fits EC10, EC20, and EC50 are estimated with two types of confidence intervals, based on the delta method and based on inverse regression; the latter results in asymmetric intervals. Only delta method-based intervals are reported in the paper. lapplyFct1 <- function(modelF) { cbind(ED(modelF, c(10,20,50), interval = "delta",display = FALSE), ED(modelF, c(10,20,50), interval = "inv",display = FALSE)[, 2:3]) } naiveEDmat1 <- matrix(NA,6,6) naiveEDmat1 <- ldply(lapply(list(gi1.W2.2, gi2.W2.2), lapplyFct1), function(lstElt) { colnames(lstElt) <- c("Est","SE","LCIdel","UCIdel","LCIinv","UCIinv") lstElt }) rownames(naiveEDmat1) <- c("Exp1.EC10","Exp1.EC20","Exp1.EC50","Exp2.EC10","Exp2.EC20","Exp2.EC50") naiveEDmat1 ## Est SE LCIdel UCIdel LCIinv UCIinv ## Exp1.EC10 21.32542 1.755714 17.65067 25.00017 17.72415 24.90412 ## Exp1.EC20 33.80420 1.303243 31.07648 36.53192 30.86826 36.37229 ## Exp1.EC50 67.78852 3.859576 59.71033 75.86670 61.50735 79.13302 ## Exp2.EC10 69.75416 9.776851 48.91530 90.59302 52.94491 92.71534 ## Exp2.EC20 95.07989 9.361550 75.12622 115.03356 76.62345 114.86085 ## Exp2.EC50 151.79665 8.249336 134.21361 169.37970 134.30111 170.44666 The confidence intervals agree quite well between the two methods for determining confidence intervals. Fitting models that accommodate the correlation The naive models may be extended to capture the correlation due to normalization. Andersen et al. (1998) derived formulas for estimating the variance-covariance matrix between all growth inhibition values, resulting in a matrix that may be input for an extended dose-response model. Of note these formulas also include the weights. 8 First functions for estimation of the variance-covariance matrix are defined and subsequently the matrices are calculated: sigma1 <- with(andersen.data, mean(SE1[1:6])) sigma2 <- with(andersen.data, mean(SE2[1:6])) vcovFct <- function(GR, SE, conVal, sigmaVal) { vcovTemp <- matrix(GR, ncol = 1)%*% matrix(GR, nrow = 1)*(sigmaValˆ2)/conValˆ4 diag(vcovTemp) <- SEˆ2/conValˆ2+GRˆ2*sigmaValˆ2/conValˆ4 vcovTemp } vcov1 <- with(andersen.data, vcovFct(GR1, SE1, control1, sigma1)) vcov2 <- with(andersen.data, vcovFct(GR2, SE2, control2, sigma2)) vcov1[1:10,1:10] ## [,1] [,2] [,3] [,4] [,5] ## [1,] 0.0002851772 0.0001908718 0.0001866906 0.0001862369 0.0001952361 ## [2,] 0.0001908718 0.0005699659 0.0001864873 0.0001860341 0.0001950235 ## [3,] 0.0001866906 0.0001864873 0.0002588285 0.0001819589 0.0001907514 ## [4,] 0.0001862369 0.0001860341 0.0001819589 0.0002471871 0.0001902879 ## [5,] 0.0001952361 0.0001950235 0.0001907514 0.0001902879 0.0003034729 ## [6,] 0.0001795434 0.0001793479 0.0001754192 0.0001749929 0.0001834488 ## [7,] 0.0002052235 0.0002050000 0.0002005094 0.0002000221 0.0002096874 ## [8,] 0.0002126011 0.0002123696 0.0002077175 0.0002072128 0.0002172255 ## [9,] 0.0001821237 0.0001819254 0.0001779402 0.0001775078 0.0001860852 ## [10,] 0.0001750607 0.0001748700 0.0001710394 0.0001706237 0.0001788685 ## [,6] [,7] [,8] [,9] [,10] ## [1,] 0.0001795434 0.0002052235 0.0002126011 0.0001821237 0.0001750607 ## [2,] 0.0001793479 0.0002050000 0.0002123696 0.0001819254 0.0001748700 ## [3,] 0.0001754192 0.0002005094 0.0002077175 0.0001779402 0.0001710394 ## [4,] 0.0001749929 0.0002000221 0.0002072128 0.0001775078 0.0001706237 ## [5,] 0.0001834488 0.0002096874 0.0002172255 0.0001860852 0.0001788685 ## [6,] 0.0008190864 0.0001928332 0.0001997654 0.0001711280 0.0001644914 ## [7,] 0.0001928332 0.0003926414 0.0002283378 0.0001956044 0.0001880186 ## [8,] 0.0001997654 0.0002283378 0.0005153918 0.0002026363 0.0001947777 ## [9,] 0.0001711280 0.0001956044 0.0002026363 0.0007675219 0.0001668554 ## [10,] 0.0001644914 0.0001880186 0.0001947777 0.0001668554 0.0008484674 Now a dose-response model, a two-parameter Weibull type 2 model, can be fitted using drm() and the additional argument varcov where the estimated variance-covariance matrix is supplied. gi1.W2.2vcov <- drm(gi1 ~Conc1, data = andersen.data, fct = W2.2(), varcov = vcov1[1:21,1:21], start = coef(gi1.W2.2)) # there are only 21 non-NA values for experiment 1 gi2.W2.2vcov <- drm(gi2 ~Conc2, data = andersen.data, fct = W2.2(), varcov = vcov2[1:17,1:17]) Checking the model fits: predMat1vcov <- predict(gi1.W2.2vcov, data.frame(Conc1 = c(c(1:10)/10,2:1000)), interval = "confidence") plot(gi1.W2.2vcov, broken = TRUE,bp=0.1,xlim = c(0,1000), ylim = c(-0.15,1.15), type = "all") lines(c(c(1:10)/10,2:1000), predMat1vcov[, 2], lty = 2)# lower limit 9 ## Exp1.EC50 1.0209225 1.020487 0.9749809 ## Exp2.EC10 0.9155723 1.043306 0.9821547 ## Exp2.EC20 0.9440913 1.039377 1.0146517 ## Exp2.EC50 0.9885376 1.032498 1.0269967 Overall the confidence intervals become wider compared to the naive approach; this observation is more pronounced for the delta method compared to inverse regression. For the delta method standard errors are slightly larger for model fits using raw data compared to models accommodating correlation and somewhat larger compared to the naive models. Normalization leads to too small standard errors and too narrow confidence intervals; the uncertainty is under-estimated. Generating table for the LaTeX manuscript: print(xtable(cbind(naiveEDmat1[, c(1,3,4)], naiveEDmat2[, c(1,3,4,7)], naiveEDmat3[, c(1,3,4,7)]), digits = c(1,1,1,1,1,1,1,2,1,1,1,2))) % latex table generated in R 4.5.1 by xtable 1.8-4 package % Thu Oct 30 13:24:58 2025 \begin{table}[ht] \centering \begin{tabular}{rrrrrrrrrrrr} \hline & Est & LCIdel & UCIdel & Est & LCIdel & UCIdel & ratiodel & Est & LCIdel & UCIdel & CIratiodel \\ \hline Exp1.EC10 & 21.3 & 17.7 & 25.0 & 18.7 & 14.2 & 23.2 & 1.22 & 18.9 & 14.2 & 23.6 & 1.28 \\ Exp1.EC20 & 33.8 & 31.1 & 36.5 & 31.3 & 27.4 & 35.3 & 1.46 & 31.5 & 27.3 & 35.6 & 1.52 \\ Exp1.EC50 & 67.8 & 59.7 & 75.9 & 68.2 & 60.5 & 76.0 & 0.96 & 68.2 & 60.3 & 76.1 & 0.98 \\ Exp2.EC10 & 69.8 & 48.9 & 90.6 & 63.5 & 43.9 & 83.1 & 0.94 & 64.2 & 43.6 & 84.8 & 0.99 \\ Exp2.EC20 & 95.1 & 75.1 & 115.0 & 88.9 & 69.4 & 108.4 & 0.98 & 89.6 & 69.2 & 110.0 & 1.02 \\ Exp2.EC50 & 151.8 & 134.2 & 169.4 & 147.9 & 130.4 & 165.4 & 0.99 & 148.4 & 130.2 & 166.5 & 1.03 \\ \hline \end{tabular} \end{table} Plotting fitted curves for growth inhibition It is possible to transform the data and the fitted dose-response curve for the growth rates to obtain a visualization in terms of growth inhibition. As a first step a helper function is defined, constructing growth inhibition from the growth rate-base model fit: gr2giFct <- function(modelFit, concVec, level = 0.95) { conLevel <- predict(modelFit, data.frame(0)) infLevel <- predict(modelFit, data.frame(Inf)) obsgi <- 1-(modelFit[["data"]][, 2]-infLevel) /(conLevel -infLevel) predgi <- 1-(predict(modelFit, data.frame(concVec)) -infLevel) /(conLevel -infLevel) parest <- coef(modelFit) dfVec <- modelFit$fct$deriv1(concVec, matrix(parest, length(concVec), length(parest), byrow = TRUE)) if (modelFit$fct$name %in% c("LL.3","W1.3","W2.3")) { parest[2]<- 1# fixing upper limit 16 beInd <- c(1,3) } if (modelFit$fct$name %in% c("LL.4","W1.4","W2.4")) { parest[c(2,3)] <- c(0,1)# fixing lower and upper limits beInd <- c(1,4) } dfVec <- modelFit$fct$deriv1(concVec, matrix(parest, length(concVec), length(parest), byrow = TRUE)) dfVec <- dfVec[, beInd] sepred <- sqrt(diag(dfVec %*% vcov(modelFit)[beInd, beInd] %*% t(dfVec))) # note that joint confidence bands could also be worked out tquan <- qt(1-(1-level)/2,df.residual(modelFit)) #confidence interval return(list(observed = 100 *obsgi, predicted = data.frame(conc = concVec, pred = 100 *predgi, se=100 *sepred, lli = 100 *(predgi -tquan *sepred), uli = 100 *(predgi +tquan *sepred)))) } Generating plot data for experiment 1: concGrid1 <- seq(0,500,by=1) plotData1 <- gr2giFct(gr1.W1.3, concGrid1) origData1 <- data.frame(conc = andersen.data$Conc1[1:21], gi=plotData1[["observed"]]) predData1raw <- as.data.frame(plotData1[["predicted"]]) predData1.0 <- 100 *predict(gi1.W2.2,data.frame(Conc1 = concGrid1), interval = "confidence") predData1nor <- as.data.frame(cbind(Conc = concGrid1, predData1.0)) predData1.1 <- 100 *predict(gi1.W2.2vcov, data.frame(Conc1 = concGrid1), interval = "confidence") predData1cor <- as.data.frame(cbind(Conc = concGrid1, predData1.1)) Making and displaying the plot: p1 <- ggplot(origData1, aes(x=conc, y=gi)) +theme_bw() + xlab(expression(paste("Dose (cl/l)"))) + ylab("Growth inhibition (%)")+ scale_x_log10(limits = c(1,500), breaks = c(1,10,100), labels = c("1","10","100")) + scale_y_continuous(limits = c(-15,100)) + coord_cartesian(clip = "off")+ geom_line(data = predData1cor, aes(x=Conc, y=Prediction), colour = "black", linetype = "dotted",size = 0.5)+ geom_line(data = predData1cor, aes(x=Conc, y=Lower), colour = "black", linetype = "dotted",size = 0.5)+ geom_line(data = predData1cor, aes(x=Conc, y=Upper), colour = "black", 17 linetype = "dotted",size = 0.5)+ geom_ribbon(data = predData1nor, aes(x=Conc, y=Prediction, ymin = Lower, ymax = pmin(Upper, 100)), alpha = 0.85,fill = "darkgrey")+ geom_line(data = predData1nor, aes(x=Conc, y=Prediction), colour = "black", linetype = "dashed")+ geom_ribbon(data = predData1raw, aes(x=conc, y=pred, ymin = lli, ymax = pmin(uli, 100)), alpha = 0.5,fill = "lightgrey")+ geom_line(data = predData1raw, aes(x=conc, y=pred), colour = "black")+ geom_point(alpha = 0.5,show.legend = FALSE)+ annotate("text",x=0,y = -15,size = 2,label = " 23%")+ annotate("text",x=5.625,y = -15,size = 2,label = "5%")+ annotate("text",x=11.25,y = -15,size = 2,label = "18%")+ annotate("text",x=22.5,y = -15,size = 2,label = "14%")+ annotate("text",x=45,y = -15,size = 2,label = "40%")+ annotate("text",x=90,y = -15,size = 2,label = "0.3%") p1 23% 5% 18% 14% 40% 0.3% 0 25 50 75 100 1 10 100 Dose (cl/l) Growth inhibition (%) Saving it as a pdf: 18 pdf("p1-v2.pdf") p1 dev.off() ## pdf ## 2 Generating plot data for experiment 2: concGrid2 <- seq(0,1500,by=10) plotData2 <- gr2giFct(gr2.W1.3, concGrid2) origData2 <- data.frame(conc = andersen.data$Conc2[1:17], gi=plotData2[["observed"]]) predData2raw <- as.data.frame(plotData2[["predicted"]]) predData2.0 <- 100 *predict(gi2.W2.2,data.frame(Conc2 = concGrid2), interval = "confidence") predData2nor <- as.data.frame(cbind(Conc = concGrid2, predData2.0)) predData2.1 <- 100 *predict(gi2.W2.2vcov, data.frame(Conc2 = concGrid2), interval = "confidence") predData2cor <- as.data.frame(cbind(Conc = concGrid2, predData2.1)) Making and displaying the plot: p2 <- ggplot(origData2, aes(x=conc, y=gi)) +theme_bw() + xlab(expression(paste("Dose (cl/l)"))) + ylab("Growth inhibition (%)")+ scale_x_log10(limits = c(1,1500), breaks = c(1,10,100,1000), labels = c("1","10","100","1000")) + scale_y_continuous(limits = c(-15,100)) + coord_cartesian(clip = "off")+ geom_line(data = predData2cor, aes(x=Conc, y=Prediction), colour = "black", linetype = "dotted",size = 0.5)+ geom_line(data = predData2cor, aes(x=Conc, y=Lower), colour = "black", linetype = "dotted",size = 0.5)+ geom_line(data = predData2cor, aes(x=Conc, y=Upper), colour = "black", linetype = "dotted",size = 0.5)+ geom_ribbon(data = predData2nor, aes(x=Conc, y=Prediction, ymin = Lower, ymax = pmin(Upper, 100)), alpha = 0.85,fill = "darkgrey")+ geom_line(data = predData2nor, aes(x=Conc, y=Prediction), colour = "black", linetype = "dashed")+ geom_ribbon(data = predData2raw, aes(x=conc, y=pred, ymin = lli, ymax = pmin(uli, 100)), alpha = 0.5,fill = "lightgrey")+ geom_line(data = predData2raw, aes(x=conc, y=pred), colour = "black")+ geom_point(alpha = 0.5,show.legend = FALSE)+ annotate("text",x=0,y = -15,size = 2,label = " 21%")+ annotate("text",x=2.89,y = -15,size = 2,label = "2%")+ annotate("text",x=5.21,y = -15,size = 2,label = "15%")+ annotate("text",x=9.38,y = -15,size = 2,label = "2%")+ 19 annotate("text",x=16.877,y = -15,size = 2,label = "26%")+ annotate("text",x=30.377,y = -15,size = 2,label = "2%")+ annotate("text",x=54.679,y = -15,size = 2,label = "10%")+ annotate("text",x=98.422,y = -15,size = 2,label = "2%")+ annotate("text",x=177.16,y = -15,size = 2,label = "5%")+ annotate("text",x=318.88,y = -15,size = 2,label = "5%")+ annotate("text",x=574,y = -15,size = 2,label = "5%")+ annotate("text",x=1033,y = -15,size = 2,label = "5%") p2 21% 2% 15% 2% 26% 2% 10% 2% 5% 5% 5% 5% 0 25 50 75 100 1 10 100 1000 Dose (cl/l) Growth inhibition (%) Saving the plot as a pdf: pdf("p2-v2.pdf") p2 dev.off() ## pdf ## 2 Simulation study By means of simulated dose-response data, which are generated from the estimates obtained for experiment 1, it’s investigated how the delta method and inverse regression perform in terms of coverage for ED10, ED20, and ED50 when the model based on raw data is fitted. 20 Preparations Fitting the three-parameter log-logistic model to data from experiments 1 and 2 to get parameter estimates and estimated ED10, ED20, and ED50 to be used for generation of simulated dose-response data in the simulations. gr1.LL.3 <- drm(GR1 ~Conc1, data = andersen.data, fct = LL.3()) plot(gr1.LL.3,broken = TRUE,bp=0.1,xlim = c(0,1000), type = "all") Conc1 GR1 0110 100 1000 1.0 1.5 2.0 2.5 summary(gr1.LL.3) ## ## Model fitted: Log-logistic (ED50 as parameter) with lower limit at 0 (3 parms) ## ## Parameter estimates: ## ## Estimate Std. Error t-value p-value ## b:(Intercept) 1.683675 0.192596 8.742 6.776e-08 *** ## d:(Intercept) 2.387384 0.040625 58.766 < 2.2e-16 *** ## e:(Intercept) 71.509086 4.156773 17.203 1.275e-12 *** ## --- ## Signif. codes: 0 '***'0.001 '**'0.01 '*'0.05 '.'0.1 ' ' 1 ## ## Residual standard error: ## ## 0.1201848 (18 degrees of freedom) edval1 <- ED(gr1.LL.3,c(10,20,50))[, 1] ## ## Estimated effective doses ## ## Estimate Std. Error ## e:1:10 19.3910 2.9954 ## e:1:20 31.3888 3.3089 21 ## e:1:50 71.5091 4.1568 edval1 ## e:1:10 e:1:20 e:1:50 ## 19.39098 31.38879 71.50909 gr2.LL.3 <- drm(GR2 ~Conc2, data = andersen.data, fct = LL.3()) plot(gr2.LL.3,broken = TRUE,bp=0.1,xlim = c(0,1000), type = "all") Conc2 GR2 0110 100 1000 0.0 0.5 1.0 1.5 2.0 summary(gr2.LL.3) ## ## Model fitted: Log-logistic (ED50 as parameter) with lower limit at 0 (3 parms) ## ## Parameter estimates: ## ## Estimate Std. Error t-value p-value ## b:(Intercept) 3.232339 0.469569 6.8836 7.518e-06 *** ## d:(Intercept) 1.993873 0.030944 64.4345 < 2.2e-16 *** ## e:(Intercept) 143.610098 8.290154 17.3230 7.470e-11 *** ## --- ## Signif. codes: 0 '***'0.001 '**'0.01 '*'0.05 '.'0.1 ' ' 1 ## ## Residual standard error: ## ## 0.1030675 (14 degrees of freedom) edval2 <- ED(gr2.LL.3,c(10,20,50))[, 1] ## ## Estimated effective doses ## ## Estimate Std. Error ## e:1:10 72.7726 9.0315 22 ## e:1:20 93.5241 8.7117 ## e:1:50 143.6101 8.2902 edval2 ## e:1:10 e:1:20 e:1:50 ## 72.77256 93.52411 143.61010 Helper function definitions Several helper functions are defined to facilitate organizing the simulation results. Defining a helper function for obtaining inverse regression-based confidence intervals for ED levels: EDintinv <- function(modelFit, p, invFct = NULL) # invFct = function(y, parameter vector) # only works for a single dose-response curve { if (!is.null(invFct)) {invFct <- inv} else {invFct <- modelFit[["fct"]][["inversion"]]} decrea <- coef(lm(resp~dose, data=as.data.frame(modelFit$dataList)))[2]<0 EDp <- as.vector(ED(modelFit, p, display = FALSE)[, 1]) fEDp <- predict(modelFit, data.frame(EDp), se.fit = FALSE,interval = "confidence") parmVec <- coef(modelFit) intMat <- t(apply(fEDp, 1, function(rowElt) { lowL <- invFct(rowElt[2], parmVec) uppL <- invFct(rowElt[3], parmVec) ciVec <- c(lowL, uppL) if (decrea) {ciVec <- rev(ciVec)} ciVec[is.nan(ciVec)] <- NA ciVec })) colnames(intMat) <- c("Lower CI","Upper CI") rownames(intMat) <- p return(intMat) } Now helper functions for running the simulations and subsequently calculating bias, CV, coverage probabilities, and length of confidence intervals are defined: simFct <- function(noSim, modelFit, concVec, seedVal, selInd, coInd, rseVal = NULL,fitOut = FALSE) { sigmaVal <- ifelse(is.null(rseVal), sqrt(summary(modelFit)$resVar), rseVal) noConc <- length(selInd) EDmat10 <- matrix(NA, noSim, 3*4) EDmat20 <- matrix(NA, noSim, 3*4) EDmat50 <- matrix(NA, noSim, 3*4) EDmat10inv <- matrix(NA, noSim, 3*3) EDmat20inv <- matrix(NA, noSim, 3*3) EDmat50inv <- matrix(NA, noSim, 3*3) EDmat10invX <- matrix(NA, noSim, 3*2) 23 EDmat20invX <- matrix(NA, noSim, 3*2) EDmat50invX <- matrix(NA, noSim, 3*2) if (fitOut) {fitList <- list()} set.seed(seedVal) for (i in 1:noSim) { respVec <- fitted(modelFit)[selInd] +rnorm(noConc, 0, sigmaVal) controlMean <- mean(respVec[coInd]) normVec <- 1-respVec /controlMean concVecTemp <- concVec[selInd] tempFit1 <- drm(normVec ~concVecTemp, fct = LL.2()) tempFit2 <- drm(normVec ~concVecTemp, fct = LL.2(), varcov = vcovFct(respVec, 1, controlMean, 1), start = coef(tempFit1)) tempFit3 <- drm(respVec ~concVecTemp, fct = LL.3()) if (fitOut) {fitList[[i]] <- list(tempFit1, tempFit2, tempFit3)} else{ EDtemp1 <- ED(tempFit1, c(10,20,50), interval = "delta",display = FALSE) EDtemp2 <- ED(tempFit2, c(10,20,50), interval = "delta",display = FALSE) EDtemp3 <- ED(tempFit3, c(10,20,50), interval = "delta",display = FALSE) EDmat10[i, 1:4]<- EDtemp1[1,1:4] EDmat20[i, 1:4]<- EDtemp1[2,1:4] EDmat50[i, 1:4]<- EDtemp1[3,1:4] EDmat10[i, 5:8]<- EDtemp2[1,1:4] EDmat20[i, 5:8]<- EDtemp2[2,1:4] EDmat50[i, 5:8]<- EDtemp2[3,1:4] EDmat10[i, 9:12]<- EDtemp3[1,1:4] EDmat20[i, 9:12]<- EDtemp3[2,1:4] EDmat50[i, 9:12]<- EDtemp3[3,1:4] EDtemp1inv <- ED(tempFit1, c(10,20,50), interval = "inv",display = FALSE) EDtemp2inv <- ED(tempFit2, c(10,20,50), interval = "inv",display = FALSE) EDtemp3inv <- ED(tempFit3, c(10,20,50), interval = "inv",display = FALSE) EDmat10inv[i, 1:3]<- EDtemp1inv[1,1:3] EDmat20inv[i, 1:3]<- EDtemp1inv[2,1:3] EDmat50inv[i, 1:3]<- EDtemp1inv[3,1:3] EDmat10inv[i, 4:6]<- EDtemp2inv[1,1:3] EDmat20inv[i, 4:6]<- EDtemp2inv[2,1:3] EDmat50inv[i, 4:6]<- EDtemp2inv[3,1:3] EDmat10inv[i, 7:9]<- EDtemp3inv[1,1:3] EDmat20inv[i, 7:9]<- EDtemp3inv[2,1:3] EDmat50inv[i, 7:9]<- EDtemp3inv[3,1:3] EDtemp1invX <- EDintinv(tempFit1, c(10,20,50)) 24 EDtemp2invX <- EDintinv(tempFit2, c(10,20,50)) EDtemp3invX <- EDintinv(tempFit3, c(10,20,50)) EDmat10invX[i, 1:2]<- EDtemp1invX[1, ] EDmat20invX[i, 1:2]<- EDtemp1invX[2, ] EDmat50invX[i, 1:2]<- EDtemp1invX[3, ] EDmat10invX[i, 3:4]<- EDtemp2invX[1, ] EDmat20invX[i, 3:4]<- EDtemp2invX[2, ] EDmat50invX[i, 3:4]<- EDtemp2invX[3, ] EDmat10invX[i, 5:6]<- EDtemp3invX[1, ] EDmat20invX[i, 5:6]<- EDtemp3invX[2, ] EDmat50invX[i, 5:6]<- EDtemp3invX[3, ] } } if (fitOut) {return(fitList)} else { return(list(EDmat10, EDmat20, EDmat50, EDmat10inv, EDmat20inv, EDmat50inv, EDmat10invX, EDmat20invX, EDmat50invX))} } biasFct0 <- function(simResElt, edval, meanSum, sdSum, logT) { biasMat0 <- matrix(NA,1,1+4*3+2) if (meanSum) {cenFct <- mean} else {cenFct <- median} if (sdSum) {scaFct <- sd} else {scaFct <- mad} biasMat0[1,1]<- edval if (logT) { estMat <- log(simResElt[, c(1,5,9)]) -log(edval) }else { estMat <- simResElt[, c(1,5,9)] -edval } biasMat0[1,c(2,6,10)] <- apply(estMat, 2, cenFct) biasMat0[1,c(4,8,12)] <- apply(estMat, 2, scaFct) if (logT) { biasMat0[1,c(3,7,11)] <- 100 *(exp(biasMat0[1,c(2,6,10)]) -1) biasMat0[1,c(5,9,13)] <- 100 *sqrt(exp(biasMat0[1,c(4,8,12)]ˆ2)-1) }else { biasMat0[1,c(3,7,11)] <- 100 *biasMat0[1,c(2,6,10)] /edval biasMat0[1,c(5,9,13)] <- 100 *biasMat0[1,c(4,8,12)] /edval } biasMat0[1,14]<- mean(simResElt[, 9]<simResElt[, 1]) biasMat0[1,15]<- mean(simResElt[, 9]<simResElt[, 5]) 25 ## [13,] 82.03751 ## [14,] 75.42575 ## [15,] 77.62648 # MAD and log scale biasFct(simRes1.test, edval1) ## True NORMd NORMdp NORMs NORMsp CORRd CORRdp ## ED10 19.39098 0.2632082 1.35737479 3.281954 16.925159 -0.04670665 -0.2408680 ## ED20 31.38879 0.1841568 0.58669630 3.544579 11.292501 -0.16847418 -0.5367337 ## ED50 71.50909 0.0420065 0.05874289 3.761001 5.259473 -0.24876805 -0.3478831 ## CORRs CORRsp RAWd RAWdp RAWs RAWsp NORMBig CORRBig ## ED10 2.940078 15.162091 0.5488003 2.8301840 2.868089 14.790841 0.44 0 ## ED20 3.134908 9.987350 0.5153057 1.6416873 3.035996 9.672230 0.44 0 ## ED50 3.543679 4.955564 0.1970671 0.2755832 3.539000 4.949021 0.44 0 biasFct(simRes1.test, edval1, TRUE,TRUE,TRUE) ## True NORMd NORMdp NORMs NORMsp CORRd ## ED10 19.39098 -0.0001121834 -0.01121771 0.16694455 16.811454 -0.013588638 ## ED20 31.38879 -0.0003576621 -0.03575981 0.11286944 11.322987 -0.010309820 ## ED50 71.50909 -0.0007773106 -0.07770085 0.05288678 5.292378 -0.004704643 ## CORRdp CORRs CORRsp RAWd RAWdp RAWs RAWsp ## ED10 -1.3496730 0.15092996 15.179359 0.017880750 1.8041568 0.14295561 14.368911 ## ED20 -1.0256856 0.10031414 10.056703 0.011854493 1.1925036 0.09509947 9.531489 ## ED50 -0.4693594 0.04998268 5.001392 0.001552538 0.1553744 0.04954543 4.957585 ## NORMBig CORRBig ## ED10 0.44 0 ## ED20 0.44 0 ## ED50 0.44 0 # Coverage and length covFct(simRes1.test, edval1) ## NORMcov NORMle NORMNA CORRcov CORRle CORRNA RAWcov RAWle RAWNA ## ED10 0.90 9.645798 0 0.98 12.16583 0 0.98 12.90117 0 ## ED20 0.92 10.268611 0 0.98 13.38473 0 0.98 14.10302 0 ## ED50 0.96 16.328786 0 1.00 17.19586 0 1.00 17.74233 0 ## NORMInv NORMil NORMiNA CORRInv CORRil CORRiNA RAWInv RAWil RAWiNA ## ED10 0.90 9.418613 0 0.98 11.78679 0 0.94 9.918046 0 ## ED20 0.92 10.132070 0 0.98 13.01963 0 0.92 10.706253 0 ## ED50 1.00 16.897137 0 1.00 17.78332 0 0.98 17.321477 0 ## NORMInv2 NORMi2l NORMi2NA CORRInv2 CORRi2l CORRiNA2 RAWInv2 RAWi2l ## ED10 0.90 9.775731 0 0.98 12.44079 0 0.96 10.13587 ## ED20 0.92 10.315847 0 0.98 13.49049 0 0.92 10.95223 ## ED50 1.00 16.423496 0 1.00 17.30607 0 1.00 16.72571 ## RAWi2NA ## ED10 0 ## ED20 0 ## ED50 0 # another test scenario based on experiment 1 (only 1 control) simRes2.test <- simFct(50, gr1.LL.3, andersen.data[["Conc1"]], 202409211, c(1,7:21), 1) biasFct(simRes2.test, edval1) 32 ## True NORMd NORMdp NORMs NORMsp CORRd CORRdp ## ED10 19.39098 2.1117308 10.8902754 6.280580 32.38919 0.1985295 1.0238242 ## ED20 31.38879 1.9048456 6.0685545 6.919063 22.04311 -0.1663657 -0.5300163 ## ED50 71.50909 0.2234496 0.3124772 5.285837 7.39184 -1.3567773 -1.8973495 ## CORRs CORRsp RAWd RAWdp RAWs RAWsp NORMBig CORRBig ## ED10 4.199597 21.657481 0.8914559 4.597272 3.849589 19.852476 0.58 0.26 ## ED20 4.644788 14.797603 0.6490386 2.067741 4.211304 13.416587 0.58 0.26 ## ED50 4.179660 5.844936 -0.8177515 -1.143563 4.169209 5.830321 0.58 0.26 biasFct(simRes2.test, edval1, TRUE,TRUE,TRUE) ## True NORMd NORMdp NORMs NORMsp CORRd ## ED10 19.39098 0.0577862649 5.94885216 0.31357778 32.144665 -0.01422774 ## ED20 31.38879 0.0366343959 3.73137053 0.21665588 21.922336 -0.01667800 ## ED50 71.50909 0.0004750369 0.04751497 0.07339615 7.349511 -0.02086675 ## CORRdp CORRs CORRsp RAWd RAWdp RAWs RAWsp ## ED10 -1.412701 0.22875814 23.178377 0.02630165 2.665059 0.1985332 20.050566 ## ED20 -1.653969 0.15438158 15.530603 0.01173129 1.180037 0.1348164 13.543132 ## ED50 -2.065055 0.05887137 5.892242 -0.01317691 -1.309048 0.0582440 5.829343 ## NORMBig CORRBig ## ED10 0.58 0.26 ## ED20 0.58 0.26 ## ED50 0.58 0.26 covFct(simRes2.test, edval1) ## NORMcov NORMle NORMNA CORRcov CORRle CORRNA RAWcov RAWle RAWNA ## ED10 0.54 11.04742 0 0.90 14.42050 0 0.92 15.76290 0 ## ED20 0.54 11.51847 0 0.88 16.09902 0 0.96 17.45509 0 ## ED50 0.92 17.55977 0 0.94 18.34671 0 0.96 19.05616 0 ## NORMInv NORMil NORMiNA CORRInv CORRil CORRiNA RAWInv RAWil RAWiNA ## ED10 0.52 10.69970 0 0.92 13.84821 0 0.80 10.59163 0 ## ED20 0.54 11.31623 0 0.88 15.39224 0 0.78 11.06356 0 ## ED50 0.92 18.28078 0 0.94 19.02148 0 0.94 17.30330 0 ## NORMInv2 NORMi2l NORMi2NA CORRInv2 CORRi2l CORRiNA2 RAWInv2 RAWi2l ## ED10 0.60 11.38338 0 0.86 14.77920 0 0.82 10.34955 ## ED20 0.54 11.62285 0 0.88 16.31081 0 0.78 11.38008 ## ED50 0.92 17.70407 0 0.94 18.49353 0 0.94 16.57424 ## RAWi2NA ## ED10 0 ## ED20 0 ## ED50 0 Merging results into a single table: bicoMerge(list(simRes1.test), edval1) # list of one ## NORMdp CORRdp RAWdp NORMsp CORRsp RAWsp NORMcov ## ED10-1 1.35737479 -0.2408680 2.8301840 16.925159 15.162091 14.790841 0.90 ## ED20-1 0.58669630 -0.5367337 1.6416873 11.292501 9.987350 9.672230 0.92 ## ED50-1 0.05874289 -0.3478831 0.2755832 5.259473 4.955564 4.949021 0.96 ## CORRcov RAWcov NORMle CORRle RAWle ## ED10-1 0.98 0.98 9.645798 12.16583 12.90117 ## ED20-1 0.98 0.98 10.268611 13.38473 14.10302 ## ED50-1 1.00 1.00 16.328786 17.19586 17.74233 33 bicoMat <- bicoMerge(list(simRes1.test, simRes2.test), edval1) bicoMat ## NORMdp CORRdp RAWdp NORMsp CORRsp RAWsp NORMcov ## ED10-1 1.35737479 -0.2408680 2.8301840 16.925159 15.162091 14.790841 0.90 ## ED10-2 10.89027544 1.0238242 4.5972718 32.389186 21.657481 19.852476 0.54 ## ED20-1 0.58669630 -0.5367337 1.6416873 11.292501 9.987350 9.672230 0.92 ## ED20-2 6.06855455 -0.5300163 2.0677405 22.043107 14.797603 13.416587 0.54 ## ED50-1 0.05874289 -0.3478831 0.2755832 5.259473 4.955564 4.949021 0.96 ## ED50-2 0.31247720 -1.8973495 -1.1435631 7.391840 5.844936 5.830321 0.92 ## CORRcov RAWcov NORMle CORRle RAWle ## ED10-1 0.98 0.98 9.645798 12.16583 12.90117 ## ED10-2 0.90 0.92 11.047416 14.42050 15.76290 ## ED20-1 0.98 0.98 10.268611 13.38473 14.10302 ## ED20-2 0.88 0.96 11.518465 16.09902 17.45509 ## ED50-1 1.00 1.00 16.328786 17.19586 17.74233 ## ED50-2 0.94 0.96 17.559773 18.34671 19.05616 Generating LaTeX code of the table: print(xtable(bicoMat)) % latex table generated in R 4.5.1 by xtable 1.8-4 package % Thu Oct 30 13:25:16 2025 \begin{table}[ht] \centering \begin{tabular}{rrrrrrrrrrrrr} \hline & NORMdp & CORRdp & RAWdp & NORMsp & CORRsp & RAWsp & NORMcov & CORRcov & RAWcov & NORMle & CORRle & RAWle \\ \hline ED10-1 & 1.36 & -0.24 & 2.83 & 16.93 & 15.16 & 14.79 & 0.90 & 0.98 & 0.98 & 9.65 & 12.17 & 12.90 \\ ED10-2 & 10.89 & 1.02 & 4.60 & 32.39 & 21.66 & 19.85 & 0.54 & 0.90 & 0.92 & 11.05 & 14.42 & 15.76 \\ ED20-1 & 0.59 & -0.54 & 1.64 & 11.29 & 9.99 & 9.67 & 0.92 & 0.98 & 0.98 & 10.27 & 13.38 & 14.10 \\ ED20-2 & 6.07 & -0.53 & 2.07 & 22.04 & 14.80 & 13.42 & 0.54 & 0.88 & 0.96 & 11.52 & 16.10 & 17.46 \\ ED50-1 & 0.06 & -0.35 & 0.28 & 5.26 & 4.96 & 4.95 & 0.96 & 1.00 & 1.00 & 16.33 & 17.20 & 17.74 \\ ED50-2 & 0.31 & -1.90 & -1.14 & 7.39 & 5.84 & 5.83 & 0.92 & 0.94 & 0.96 & 17.56 & 18.35 & 19.06 \\ \hline \end{tabular} \end{table} print(xtable(bicoMat), file = "testRes.txt") Experiment 1: Simulation results Simulations based on experiment 1 where there were 3 replicates per non-zero concentration. # 6 control measurements simRes1.6 <- simFct(1000, gr1.LL.3, andersen.data[["Conc1"]], 202409231,1:21,1:6) bicoMerge(list(simRes1.6), edval1) NORMdp CORRdp RAWdp NORMsp CORRsp RAWsp NORMcov ED10-1 0.39908002 -2.1659948 0.5854608 17.158873 15.813674 15.82218 0.828 ED20-1 0.06416267 -1.7779236 0.1824267 11.716948 10.757256 10.70525 0.817 ED50-1 0.34597428 -0.3135628 0.2562561 6.125883 5.946344 5.94285 0.934 CORRcov RAWcov NORMle CORRle RAWle 34 ED10-1 0.921 0.938 9.380123 11.67202 12.36773 ED20-1 0.926 0.944 10.007651 12.91356 13.60797 ED50-1 0.943 0.953 16.180481 16.98157 17.52696 print(xtable(bicoMerge(list(simRes1.6), edval1)), file = "simRes1-6.txt") # 1 control measurement simRes1.1 <- simFct(1000, gr1.LL.3, andersen.data[["Conc1"]], 202409232,c(1,7:21), 1) bicoMerge(list(simRes1.1), edval1) NORMdp CORRdp RAWdp NORMsp CORRsp RAWsp NORMcov ED10-1 2.7340337 -3.997766 0.3893737 31.808699 21.062389 19.30030 0.572 ED20-1 1.0346333 -3.432319 -0.2149269 22.540967 14.934043 13.35782 0.541 ED50-1 0.2997428 -1.212291 -0.2452038 8.718556 6.682077 6.31025 0.856 CORRcov RAWcov NORMle CORRle RAWle ED10-1 0.875 0.928 10.33237 14.12560 15.50051 ED20-1 0.886 0.937 10.98018 16.03298 17.41955 ED50-1 0.935 0.960 17.69472 18.64535 19.35072 print(xtable(bicoMerge(list(simRes1.1), edval1)), file = "simRes1-1.txt") # 3 control measurements simRes1.3 <- simFct(1000, gr1.LL.3, andersen.data[["Conc1"]], 202409233,c(1:3,7:21), 1:3) bicoMerge(list(simRes1.3), edval1) NORMdp CORRdp RAWdp NORMsp CORRsp RAWsp NORMcov ED10-1 0.5479027 -2.6513363 0.7124780 21.033329 18.229694 18.088154 0.750 ED20-1 0.0107314 -2.2624996 0.1444838 14.398477 12.411719 12.186125 0.729 ED50-1 0.2656506 -0.5781913 0.1263378 6.357336 5.898423 5.867177 0.920 CORRcov RAWcov NORMle CORRle RAWle ED10-1 0.891 0.931 9.548562 12.65577 13.66049 ED20-1 0.902 0.939 10.192509 14.14441 15.17225 ED50-1 0.929 0.945 16.511070 17.54822 18.27012 print(xtable(bicoMerge(list(simRes1.3), edval1)), file = "simRes1-3.txt") # 9 control measurements simRes1.9 <- simFct(1000, gr1.LL.3, andersen.data[["Conc1"]], 202409234, c(1,1,2,2,3,3,4,5,6,7:21), 1:9) bicoMerge(list(simRes1.9), edval1) NORMdp CORRdp RAWdp NORMsp CORRsp RAWsp NORMcov ED10-1 1.05910612 -1.3344055 1.0391455 14.971842 14.092655 14.079365 0.883 ED20-1 0.39157829 -1.3261244 0.3622014 10.109135 9.527545 9.478498 0.869 ED50-1 -0.06738513 -0.6693585 -0.1700105 5.608293 5.530988 5.534504 0.933 CORRcov RAWcov NORMle CORRle RAWle ED10-1 0.942 0.953 9.268770 11.02243 11.54526 ED20-1 0.932 0.950 9.870478 12.10165 12.61563 ED50-1 0.942 0.953 15.781412 16.42359 16.86089 print(xtable(bicoMerge(list(simRes1.9), edval1)), file = "simRes1-9.txt") # 12 control measurements simRes1.12 <- simFct(1000, gr1.LL.3, andersen.data[["Conc1"]], 202409235, c(1,1,2,2,3,3,4,4,5,5,6,6,7:21), 1:12) 35 bicoMerge(list(simRes1.12), edval1) NORMdp CORRdp RAWdp NORMsp CORRsp RAWsp NORMcov ED10-1 -0.03647479 -2.262355 -0.04742255 14.532883 14.116132 14.166109 0.886 ED20-1 -0.12369137 -1.733781 -0.15268420 9.885059 9.607884 9.602101 0.877 ED50-1 0.39145403 -0.158222 0.31050521 6.073844 6.025582 6.025693 0.931 CORRcov RAWcov NORMle CORRle RAWle ED10-1 0.927 0.945 9.237832 10.68949 11.12399 ED20-1 0.926 0.949 9.868287 11.73391 12.15423 ED50-1 0.937 0.947 16.028498 16.57520 16.94624 print(xtable(bicoMerge(list(simRes1.12), edval1)), file = "simRes1-12.txt") # 18 control measurements simRes1.18 <- simFct(1000, gr1.LL.3, andersen.data[["Conc1"]], 202409735,# 202409335, c(1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,6,6,6,7:21), 1:18) bicoMerge(list(simRes1.18), edval1) NORMdp CORRdp RAWdp NORMsp CORRsp RAWsp NORMcov ED10-1 0.6702582 -1.2729723 0.7429826 13.524698 13.303317 13.367104 0.906 ED20-1 0.4358704 -0.9653266 0.4688307 9.097648 8.949316 8.953595 0.906 ED50-1 0.6228156 0.1575785 0.5803782 5.479626 5.434090 5.433840 0.957 CORRcov RAWcov NORMle CORRle RAWle ED10-1 0.942 0.953 9.241050 10.31418 10.63964 ED20-1 0.947 0.956 9.850238 11.23334 11.53793 ED50-1 0.955 0.963 15.902291 16.31402 16.59000 print(xtable(bicoMerge(list(simRes1.18), edval1)), file = "simRes1-18.txt") bicoMat1 <- bicoMerge(list(simRes1.1, simRes1.3, simRes1.6, simRes1.9, simRes1.12, simRes1.18), edval1) bicoMat1 NORMdp CORRdp RAWdp NORMsp CORRsp RAWsp NORMcov ED10-1 2.73403371 -3.9977665 0.38937373 31.808699 21.062389 19.300298 0.572 ED10-2 0.54790270 -2.6513363 0.71247801 21.033329 18.229694 18.088154 0.750 ED10-3 0.39908002 -2.1659948 0.58546080 17.158873 15.813674 15.822182 0.828 ED10-4 1.05910612 -1.3344055 1.03914551 14.971842 14.092655 14.079365 0.883 ED10-5 -0.03647479 -2.2623548 -0.04742255 14.532883 14.116132 14.166109 0.886 ED10-6 0.67025820 -1.2729723 0.74298256 13.524698 13.303317 13.367104 0.906 ED20-1 1.03463327 -3.4323190 -0.21492686 22.540967 14.934043 13.357818 0.541 ED20-2 0.01073140 -2.2624996 0.14448384 14.398477 12.411719 12.186125 0.729 ED20-3 0.06416267 -1.7779236 0.18242665 11.716948 10.757256 10.705247 0.817 ED20-4 0.39157829 -1.3261244 0.36220138 10.109135 9.527545 9.478498 0.869 ED20-5 -0.12369137 -1.7337809 -0.15268420 9.885059 9.607884 9.602101 0.877 ED20-6 0.43587039 -0.9653266 0.46883066 9.097648 8.949316 8.953595 0.906 ED50-1 0.29974276 -1.2122912 -0.24520382 8.718556 6.682077 6.310250 0.856 ED50-2 0.26565056 -0.5781913 0.12633784 6.357336 5.898423 5.867177 0.920 ED50-3 0.34597428 -0.3135628 0.25625609 6.125883 5.946344 5.942850 0.934 ED50-4 -0.06738513 -0.6693585 -0.17001053 5.608293 5.530988 5.534504 0.933 ED50-5 0.39145403 -0.1582220 0.31050521 6.073844 6.025582 6.025693 0.931 ED50-6 0.62281558 0.1575785 0.58037817 5.479626 5.434090 5.433840 0.957 CORRcov RAWcov NORMle CORRle RAWle ED10-1 0.875 0.928 10.332366 14.12560 15.50051 ED10-2 0.891 0.931 9.548562 12.65577 13.66049 36 ED10-3 0.921 0.938 9.380123 11.67202 12.36773 ED10-4 0.942 0.953 9.268770 11.02243 11.54526 ED10-5 0.927 0.945 9.237832 10.68949 11.12399 ED10-6 0.942 0.953 9.241050 10.31418 10.63964 ED20-1 0.886 0.937 10.980178 16.03298 17.41955 ED20-2 0.902 0.939 10.192509 14.14441 15.17225 ED20-3 0.926 0.944 10.007651 12.91356 13.60797 ED20-4 0.932 0.950 9.870478 12.10165 12.61563 ED20-5 0.926 0.949 9.868287 11.73391 12.15423 ED20-6 0.947 0.956 9.850238 11.23334 11.53793 ED50-1 0.935 0.960 17.694717 18.64535 19.35072 ED50-2 0.929 0.945 16.511070 17.54822 18.27012 ED50-3 0.943 0.953 16.180481 16.98157 17.52696 ED50-4 0.942 0.953 15.781412 16.42359 16.86089 ED50-5 0.937 0.947 16.028498 16.57520 16.94624 ED50-6 0.955 0.963 15.902291 16.31402 16.59000 print(xtable(bicoMat1)) % latex table generated in R 4.5.1 by xtable 1.8-4 package % Thu Oct 30 13:38:13 2025 \begin{table}[ht] \centering \begin{tabular}{rrrrrrrrrrrrr} \hline & NORMdp & CORRdp & RAWdp & NORMsp & CORRsp & RAWsp & NORMcov & CORRcov & RAWcov & NORMle & CORRle & RAWle \\ \hline ED10-1 & 2.73 & -4.00 & 0.39 & 31.81 & 21.06 & 19.30 & 0.57 & 0.88 & 0.93 & 10.33 & 14.13 & 15.50 \\ ED10-2 & 0.55 & -2.65 & 0.71 & 21.03 & 18.23 & 18.09 & 0.75 & 0.89 & 0.93 & 9.55 & 12.66 & 13.66 \\ ED10-3 & 0.40 & -2.17 & 0.59 & 17.16 & 15.81 & 15.82 & 0.83 & 0.92 & 0.94 & 9.38 & 11.67 & 12.37 \\ ED10-4 & 1.06 & -1.33 & 1.04 & 14.97 & 14.09 & 14.08 & 0.88 & 0.94 & 0.95 & 9.27 & 11.02 & 11.55 \\ ED10-5 & -0.04 & -2.26 & -0.05 & 14.53 & 14.12 & 14.17 & 0.89 & 0.93 & 0.94 & 9.24 & 10.69 & 11.12 \\ ED10-6 & 0.67 & -1.27 & 0.74 & 13.52 & 13.30 & 13.37 & 0.91 & 0.94 & 0.95 & 9.24 & 10.31 & 10.64 \\ ED20-1 & 1.03 & -3.43 & -0.21 & 22.54 & 14.93 & 13.36 & 0.54 & 0.89 & 0.94 & 10.98 & 16.03 & 17.42 \\ ED20-2 & 0.01 & -2.26 & 0.14 & 14.40 & 12.41 & 12.19 & 0.73 & 0.90 & 0.94 & 10.19 & 14.14 & 15.17 \\ ED20-3 & 0.06 & -1.78 & 0.18 & 11.72 & 10.76 & 10.71 & 0.82 & 0.93 & 0.94 & 10.01 & 12.91 & 13.61 \\ ED20-4 & 0.39 & -1.33 & 0.36 & 10.11 & 9.53 & 9.48 & 0.87 & 0.93 & 0.95 & 9.87 & 12.10 & 12.62 \\ ED20-5 & -0.12 & -1.73 & -0.15 & 9.89 & 9.61 & 9.60 & 0.88 & 0.93 & 0.95 & 9.87 & 11.73 & 12.15 \\ ED20-6 & 0.44 & -0.97 & 0.47 & 9.10 & 8.95 & 8.95 & 0.91 & 0.95 & 0.96 & 9.85 & 11.23 & 11.54 \\ ED50-1 & 0.30 & -1.21 & -0.25 & 8.72 & 6.68 & 6.31 & 0.86 & 0.94 & 0.96 & 17.69 & 18.65 & 19.35 \\ ED50-2 & 0.27 & -0.58 & 0.13 & 6.36 & 5.90 & 5.87 & 0.92 & 0.93 & 0.94 & 16.51 & 17.55 & 18.27 \\ ED50-3 & 0.35 & -0.31 & 0.26 & 6.13 & 5.95 & 5.94 & 0.93 & 0.94 & 0.95 & 16.18 & 16.98 & 17.53 \\ ED50-4 & -0.07 & -0.67 & -0.17 & 5.61 & 5.53 & 5.53 & 0.93 & 0.94 & 0.95 & 15.78 & 16.42 & 16.86 \\ ED50-5 & 0.39 & -0.16 & 0.31 & 6.07 & 6.03 & 6.03 & 0.93 & 0.94 & 0.95 & 16.03 & 16.58 & 16.95 \\ ED50-6 & 0.62 & 0.16 & 0.58 & 5.48 & 5.43 & 5.43 & 0.96 & 0.95 & 0.96 & 15.90 & 16.31 & 16.59 \\ \hline \end{tabular} \end{table} print(xtable(bicoMat1), file = "simRes1.txt") Showing histograms of estimated ED10, ED20, and ED50 values for 3 scenarios (1, 3, and 6 control measurements). For the normalization approach: 37 # ED10 par(mfrow = c(2,2)) hist(simRes1.1[[1]][, 1]) hist(simRes1.3[[1]][, 1]) hist(simRes1.6[[1]][, 1]) par(mfrow = c(1,1)) Histogram of simRes1.1[[1]][, 1] simRes1.1[[1]][, 1] Frequency 0 10 20 30 40 0 150 Histogram of simRes1.3[[1]][, 1] simRes1.3[[1]][, 1] Frequency 10 15 20 25 30 35 0 100 Histogram of simRes1.6[[1]][, 1] simRes1.6[[1]][, 1] Frequency 10 15 20 25 30 0 100 # ED20 par(mfrow = c(2,2)) hist(simRes1.1[[2]][, 1]) hist(simRes1.3[[2]][, 1]) hist(simRes1.6[[2]][, 1]) par(mfrow = c(1,1)) 38 Histogram of simRes1.1[[2]][, 1] simRes1.1[[2]][, 1] Frequency 10 20 30 40 50 0 150 Histogram of simRes1.3[[2]][, 1] simRes1.3[[2]][, 1] Frequency 15 20 25 30 35 40 45 0 100 Histogram of simRes1.6[[2]][, 1] simRes1.6[[2]][, 1] Frequency 20 25 30 35 40 45 0 100 # ED50 par(mfrow = c(2,2)) hist(simRes1.1[[3]][, 1]) hist(simRes1.3[[3]][, 1]) hist(simRes1.6[[3]][, 1]) par(mfrow = c(1,1)) 39 Histogram of simRes1.1[[3]][, 1] simRes1.1[[3]][, 1] Frequency 50 60 70 80 90 0 150 Histogram of simRes1.3[[3]][, 1] simRes1.3[[3]][, 1] Frequency 60 70 80 90 0 200 Histogram of simRes1.6[[3]][, 1] simRes1.6[[3]][, 1] Frequency 60 70 80 90 0 200 For the growth rate-based approach: # ED10 par(mfrow = c(2,2)) hist(simRes1.1[[1]][, 9]) hist(simRes1.3[[1]][, 9]) hist(simRes1.6[[1]][, 9]) par(mfrow = c(1,1)) 40 Histogram of simRes1.1[[1]][, 9] simRes1.1[[1]][, 9] Frequency 5 10 15 20 25 30 0 100 Histogram of simRes1.3[[1]][, 9] simRes1.3[[1]][, 9] Frequency 10 15 20 25 30 0 100 Histogram of simRes1.6[[1]][, 9] simRes1.6[[1]][, 9] Frequency 10 15 20 25 30 0 100 250 # ED20 par(mfrow = c(2,2)) hist(simRes1.1[[2]][, 9]) hist(simRes1.3[[2]][, 9]) hist(simRes1.6[[2]][, 9]) par(mfrow = c(1,1)) 41 ED50-4 0.922 0.942 46.25279 35.99998 33.46842 ED50-5 0.925 0.937 41.99686 522.02100 32.94533 ED50-6 0.938 0.956 39.25654 39.08761 32.38579 print(xtable(bicoMat2)) % latex table generated in R 4.5.1 by xtable 1.8-4 package % Thu Oct 30 13:53:26 2025 \begin{table}[ht] \centering \begin{tabular}{rrrrrrrrrrrrr} \hline & NORMdp & CORRdp & RAWdp & NORMsp & CORRsp & RAWsp & NORMcov & CORRcov & RAWcov & NORMle & CORRle & RAWle \\ \hline ED10-1 & 2.60 & 1.65 & 0.91 & 25.01 & 20.27 & 13.58 & 0.89 & 0.93 & 0.95 & 61.77 & 50.39 & 41.70 \\ ED10-2 & 5.10 & 3.67 & 1.03 & 26.99 & 24.51 & 13.33 & 0.88 & 0.92 & 0.94 & 108.22 & 58.74 & 39.45 \\ ED10-3 & 6.77 & 4.69 & 2.07 & 26.38 & 23.55 & 12.65 & 0.90 & 0.92 & 0.95 & 230.71 & 53.37 & 38.77 \\ ED10-4 & 5.42 & 3.46 & 0.89 & 26.23 & 23.18 & 12.62 & 0.91 & 0.93 & 0.95 & 136.28 & 49.19 & 37.51 \\ ED10-5 & 5.31 & 2.75 & 0.70 & 26.01 & 21.98 & 12.40 & 0.91 & 0.94 & 0.96 & 60.69 & 477.14 & 36.73 \\ ED10-6 & 6.84 & 5.09 & 1.38 & 27.64 & 25.91 & 12.43 & 0.91 & 0.92 & 0.95 & 96.48 & 62.79 & 36.17 \\ ED20-1 & 1.15 & 0.66 & 0.45 & 17.53 & 13.75 & 10.02 & 0.88 & 0.93 & 0.95 & 54.98 & 44.83 & 39.39 \\ ED20-2 & 2.81 & 1.95 & 0.54 & 17.75 & 16.05 & 9.75 & 0.89 & 0.92 & 0.94 & 86.22 & 50.51 & 37.17 \\ ED20-3 & 3.97 & 2.68 & 1.32 & 17.07 & 15.29 & 9.17 & 0.90 & 0.92 & 0.94 & 184.81 & 46.46 & 36.39 \\ ED20-4 & 3.08 & 1.86 & 0.51 & 17.03 & 15.14 & 9.24 & 0.90 & 0.92 & 0.95 & 109.51 & 43.22 & 35.28 \\ ED20-5 & 2.99 & 1.38 & 0.36 & 16.88 & 14.40 & 9.12 & 0.90 & 0.94 & 0.95 & 51.37 & 492.36 & 34.54 \\ ED20-6 & 4.01 & 2.88 & 0.81 & 17.75 & 16.70 & 9.06 & 0.90 & 0.92 & 0.94 & 78.28 & 53.01 & 33.87 \\ ED50-1 & -0.21 & -0.26 & 0.05 & 8.69 & 6.78 & 5.99 & 0.92 & 0.93 & 0.94 & 48.38 & 38.48 & 36.61 \\ ED50-2 & 0.24 & 0.11 & 0.09 & 7.82 & 7.01 & 5.84 & 0.90 & 0.92 & 0.94 & 52.71 & 38.65 & 34.88 \\ ED50-3 & 0.56 & 0.31 & 0.41 & 7.82 & 7.02 & 5.60 & 0.92 & 0.93 & 0.95 & 45.57 & 37.24 & 34.34 \\ ED50-4 & 0.38 & 0.14 & 0.20 & 7.74 & 7.02 & 5.72 & 0.91 & 0.92 & 0.94 & 46.25 & 36.00 & 33.47 \\ ED50-5 & 0.30 & -0.05 & 0.12 & 7.81 & 7.00 & 5.75 & 0.90 & 0.93 & 0.94 & 42.00 & 522.02 & 32.95 \\ ED50-6 & 0.60 & 0.31 & 0.17 & 7.48 & 7.03 & 5.48 & 0.92 & 0.94 & 0.96 & 39.26 & 39.09 & 32.39 \\ \hline \end{tabular} \end{table} print(xtable(bicoMat2), file = "simRes2.txt") Again, showing histograms of estimated ED10, ED20, and ED50 values for 3 scenarios (1, 3, and 6 control measurements). For the normalization approach: # ED10 par(mfrow = c(2,2)) hist(simRes2.1[[1]][, 1]) hist(simRes2.3[[1]][, 1]) hist(simRes2.6[[1]][, 1]) par(mfrow = c(1,1)) 48 Histogram of simRes2.1[[1]][, 1] simRes2.1[[1]][, 1] Frequency 50 100 150 0 150 Histogram of simRes2.3[[1]][, 1] simRes2.3[[1]][, 1] Frequency 40 60 80 120 160 0 150 350 Histogram of simRes2.6[[1]][, 1] simRes2.6[[1]][, 1] Frequency 40 60 80 120 160 0 200 # ED20 par(mfrow = c(2,2)) hist(simRes2.1[[2]][, 1]) hist(simRes2.3[[2]][, 1]) hist(simRes2.6[[2]][, 1]) par(mfrow = c(1,1)) 49 Histogram of simRes2.1[[2]][, 1] simRes2.1[[2]][, 1] Frequency 60 80 100 140 180 0 150 Histogram of simRes2.3[[2]][, 1] simRes2.3[[2]][, 1] Frequency 60 80 100 140 180 0 150 350 Histogram of simRes2.6[[2]][, 1] simRes2.6[[2]][, 1] Frequency 60 80 100 140 180 0 200 # ED50 par(mfrow = c(2,2)) hist(simRes2.1[[3]][, 1]) hist(simRes2.3[[3]][, 1]) hist(simRes2.6[[3]][, 1]) par(mfrow = c(1,1)) 50 Histogram of simRes2.1[[3]][, 1] simRes2.1[[3]][, 1] Frequency 100 120 140 160 180 0 150 350 Histogram of simRes2.3[[3]][, 1] simRes2.3[[3]][, 1] Frequency 100 120 140 160 180 0 100 Histogram of simRes2.6[[3]][, 1] simRes2.6[[3]][, 1] Frequency 100 120 140 160 180 0 200 For the growth rate-based approach: # ED10 par(mfrow = c(2,2)) hist(simRes2.1[[1]][, 9]) hist(simRes2.3[[1]][, 9]) hist(simRes2.6[[1]][, 9]) par(mfrow = c(1,1)) 51 Histogram of simRes2.1[[1]][, 9] simRes2.1[[1]][, 9] Frequency 40 60 80 100 0 100 Histogram of simRes2.3[[1]][, 9] simRes2.3[[1]][, 9] Frequency 50 60 70 80 90 110 0 100 Histogram of simRes2.6[[1]][, 9] simRes2.6[[1]][, 9] Frequency 50 60 70 80 90 110 0 100 # ED20 par(mfrow = c(2,2)) hist(simRes2.1[[2]][, 9]) hist(simRes2.3[[2]][, 9]) hist(simRes2.6[[2]][, 9]) par(mfrow = c(1,1)) 52 Histogram of simRes2.1[[2]][, 9] simRes2.1[[2]][, 9] Frequency 60 80 100 120 0 100 250 Histogram of simRes2.3[[2]][, 9] simRes2.3[[2]][, 9] Frequency 70 80 90 100 120 0 100 Histogram of simRes2.6[[2]][, 9] simRes2.6[[2]][, 9] Frequency 70 80 90 100 120 0 100 250 # ED50 par(mfrow = c(2,2)) hist(simRes2.1[[3]][, 9]) hist(simRes2.3[[3]][, 9]) hist(simRes2.6[[3]][, 9]) par(mfrow = c(1,1)) 53 Histogram of simRes2.1[[3]][, 9] simRes2.1[[3]][, 9] Frequency 120 140 160 180 0 100 Histogram of simRes2.3[[3]][, 9] simRes2.3[[3]][, 9] Frequency 120 140 160 0 100 Histogram of simRes2.6[[3]][, 9] simRes2.6[[3]][, 9] Frequency 120 130 140 150 160 170 0 100 Comparison of normalization and growth rate-based approaches By means of simulated data a graphical comparison of the two approaches is obtained, which also shows how the ED50 values differ. 100 simulated datasets are generated based on parameter estimates from a 3-parameter log-logistic model fitted to the data from experiment 1: noSim <- 1000 simRes1.test <- simFct(noSim, gr1.LL.3, andersen.data[["Conc1"]], 202510283,1:21,1:6, rseVal = 0.3,fitOut = TRUE) Differences between normalization and growth rate-based approaches over a grid for each of the 100 simulated datasets are calculated: concGridX <- seq(10,500,by=2) lapres1 <- lapply(simRes1.test, function(listElt){100*predict(listElt[[1]], data.frame(concGridX))}) lapres2 <- as.data.frame(lapres1) colnames(lapres2) <- 1:noSim lapFct <- function(listElt) { plotDataX <- gr2giFct(listElt[[3]], concGridX) plotDataX[["predicted"]][, "pred"] } lapres3 <- as.data.frame(lapply(simRes1.test, lapFct)) colnames(lapres3) <- 1:noSim lapres23 <- lapres3 -lapres2 # differences between the two approaches 54 f0raw <- gr2giFct(gr1.LL.3, concGridX)[["predicted"]][, "pred"] # the true curve based on growth rate data p0Data <- data.frame(conc = concGridX, y0raw = f0raw, diff = apply(lapres23, 1, mean)) The fitted dose-response curve based on growth rates as well as the differences to the normalization approach are plotted: p0 <- ggplot(origData1, aes(x=conc, y=gi)) +theme_bw() +theme(legend.position = "none")+ xlab(expression(paste("Dose (cl/l)"))) + ylab("Growth inhibition (%)")+ scale_x_log10(limits = c(45,210), breaks = c(1,10,50,100,200), labels = c("1","10","50","100","200")) + scale_y_continuous(limits = c(0,100)) + coord_cartesian(clip = "off")+ geom_line(data = p0Data, aes(x=conc, y=y0raw), colour = "black", linetype = "solid")+ geom_line(data = p0Data, aes(x=conc, y=y0raw -diff), colour = "darkgrey", linetype = "dashed")+ annotate("text",x=85,y=75,size = 4,label = "Growth rate-based approach")+ annotate("text",x=125,y=55,size = 4,label = "Normalization approach", colour = "darkgrey")+ geom_segment(aes(x=0,y=50,xend = 71.5091,yend = 50), linetype = "dotted",size = 0.5)+ geom_segment(aes(x=71.5091,y=50,xend = 71.5091,yend = 0), linetype = "dashed",size = 0.5)+ annotate("text",x=59,y=0,size = 4,label = "ED50 = 71.5 cl/l")+ geom_segment(aes(x=74,y=50,xend = 74,yend = 0), linetype = "dashed",size = 0.5, colour = "darkgrey")+ annotate("text",x=89,y=0,size = 4,label = "ED50 = 74.0 cl/l", colour = "darkgrey")+ geom_point(data = data.frame(x = c(71.5,74), y = c(0,0)), aes(x=x, y=y), colour = "black",shape = 19) p0 55 Growth rate−based approach Normalization approach ED50 = 71.5 cl/l ED50 = 74.0 cl/l 0 25 50 75 100 50 100 200 Dose (cl/l) Growth inhibition (%) Saving the plot as a pdf: pdf("p0-v1.pdf") p0 dev.off() ## pdf ## 2 References • Andersen, J. S., Holst, H., Spliid, H., Andersen, H., Baun, A., Nyholm, N. (1998). Continuous Ecotoxicological Data Evaluated Relative to a Control Response. Journal of Agricultural, Biological, and Environmental Statistics,3, 405-420. http://www.jstor.org/stable/1400573 • R Core Team (2025). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.r-project.org/. 56