diff --git a/Plot4PL.png b/Plot4PL.png new file mode 100644 index 0000000..83285dd Binary files /dev/null and b/Plot4PL.png differ diff --git a/TestF4PLGHZ1.numbers b/TestF4PLGHZ1.numbers new file mode 100644 index 0000000..05fdcdd Binary files /dev/null and b/TestF4PLGHZ1.numbers differ diff --git a/TestF4PLGHZ1.xlsx b/TestF4PLGHZ1.xlsx new file mode 100644 index 0000000..60dcc1d Binary files /dev/null and b/TestF4PLGHZ1.xlsx differ diff --git a/TestF4PLGHZ2.xlsx b/TestF4PLGHZ2.xlsx new file mode 100644 index 0000000..7980463 Binary files /dev/null and b/TestF4PLGHZ2.xlsx differ diff --git a/TestFileRoxygen.R b/TestFileRoxygen.R new file mode 100644 index 0000000..e2537cd --- /dev/null +++ b/TestFileRoxygen.R @@ -0,0 +1,219 @@ +dat <- data.frame(REF1=c(1.1,1.2,2.1,3,5,6,8.1,9), REF1=c(1.2,1.5,2.1,3.1,4.9,6.1,8.3,9.1), + TEST1=c(1,1.3,2.5,3.5,5.9,6.9,8.1,9.1), TEST2=c(1.4,1.2,2.6,3.4,5.8,6.7,8.6,9.3), log_dose=c(1,0,-1,-2,-3,-4,-5,-6)) + +all_l <- melt(data.frame(dat), id.vars="log_dose", variable.name="replname", value.name = "readout") +isRef <- rep(c(1,0),1,each=nrow(all_l)/2) +isSample <- rep(c(0,1),1,each=nrow(all_l)/2) +all_l2 <- cbind(all_l, isRef, isSample) + +startlistmu <- list(as=1, bs=-1,cs=-3, + ds=10,at=1, bt=-1, + dt=10,r=0) + + +mu <- gsl_nls(fn = readout ~ as*isRef + at*isSample + (ds*isRef + dt*isSample - as*isRef - at*isSample)/ + (1+isRef*exp(bs*(cs - log_dose)) + isSample*exp(bt*(cs-r*isSample-log_dose))), + data=all_l, + start=startlistmu, + control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) +s_mu <- summary(mu)$coefficients[,1] +s_mu + + +CIRC <- data.frame(log_dose = c(-2.5,-2.5,-2.5, -3.2,-3.2,-3.2,-3.9,-3.9,-3.9, +-3.2,-3.2,-3.2,-3.9,-3.9,-3.9,-4.7,-4.7,-4.7), +replname= c("R_dil1","R_dil1","R_dil1", "R_dil2","R_dil2","R_dil2", "R_dil3","R_dil3","R_dil3", + "T_dil1","T_dil1","T_dil1", "T_dil2","T_dil2","T_dil2", "T_dil3","T_dil3","T_dil3"), + readout = c(72.1,75.8,76.04,59.8,61,62.7,43.6,45,41.5,53.5,62.2,65.9,48.3,43.8,43.14,28.17,29.2,31.2), + isRef=c(rep(1,9), rep(0,9)), + isSample = c(rep(0,9), rep(1,9))) +Lim <- c(rep(0,8), 70,130) +PureErrF <- TRUE + + +LinPotTab(circles=CIRC, Lim, PureErrF) + +# Potency lower 95%CI upper 95%CI test_result lowerRel95%CI upperRel95%CI +# isSample 104.959 87.994 125.196 0 83.836 119.281 + + +as=3; bs=1; cs=-4; ds=10; at=3; bt=1; dt=10;r=0.0001;ct=cs-r;sd_fac=0.1; gt=1; gs=1; +lnConc=c(1,0,-1,-2,-3,-4,-5,-6) +heteroNoise <- FALSE +noDilS <- 3 +noD <- 8 +Calc_DilRes(as,bs,cs,ds,at,bt,dt,r,ct,sd_fac,gt,gs,log_conc=lnConc, heteroNoise, noDilS, noD) + +# R_dil1 R_dil2 R_dil3 T_dil1 T_dil2 T_dil3 log_dose +# [1,] 9.953324 9.958311 9.952765 9.964621 9.951204 9.962239 1 +# [2,] 9.878524 9.877231 9.884573 9.866374 9.874134 9.878356 0 +# [3,] 9.654425 9.671862 9.680914 9.665887 9.670981 9.680247 -1 +# [4,] 9.174863 9.167155 9.143036 9.177929 9.165979 9.163067 -2 +# [5,] 8.120441 8.122103 8.110990 8.107848 8.124945 8.109104 -3 +# [6,] 6.495730 6.497168 6.501670 6.503780 6.497561 6.500025 -4 +# [7,] 4.886209 4.877895 4.883961 4.890778 4.874703 4.889223 -5 +# [8,] 3.829636 3.819950 3.830449 3.836558 3.831859 3.828217 -6 + +dat <- data.frame(REF1=c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2=c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497), +REF3=c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665), +TEST1=c(1405, 1523, 1502, 1474, 2383, 3221, 3589, 3445), TEST2=c(1420, 1516, 1544, 1512, 2226, 3219, 3327, 3591), +TEST3=c(1399, 1376, 1588, 1475, 2148, 3083, 2942, 3466), log_dose=c(5.01,3.401,2.708,2.015,1.32176,0.62861,-0.0645385,-1.6739764)) +TransF<- FALSE +Dat<- list() +p <- plot_f(dat,TransF) +print(p) +# Plot4PL.png +# Error in chol2inv(object$m$Rmat()) : +# element (2, 2) is zero, so the inverse cannot be computed +# In addition: +# Warning messages: +# 1: In nlsModel(formula, mf, cFit, .swts, jac) : +# singular gradient matrix at parameter estimates +# 2: In nlsModel(formula, mf, cFit, .swts, jac) : +# singular gradient matrix at parameter estimates + + +dat <- data.frame(R_dil1 =c(10221, 18258, 31993, 49336, 68332, 83527, 95584, 102229), +R_dil2=c( 10136, 19078, 31925, 49003, 68034, 83776, 95495, 101608), +T_dil1=c( 10830, 19891, 33915, 52131, 70617, 85784, 95937, 102791), +T_dil2=c( 11169, 20153, 34007, 52179, 69962, 85543, 96439, 102655), +log_dose=c( -1.2029, -1.89712, -2.590267, -3.2834, -3.97656, -4.66917, -5.362323, -6.05334)) +CIRC <- data.frame(log_dose=c( -2.590267, -2.590267, -3.2834 , -3.2834,-3.97656, -3.97656, -2.590267, -2.590267,-3.2834, -3.2834, -3.97656, -3.97656), +replname= c("R_dil1","R_dil2","R_dil1", "R_dil2","R_dil1","R_dil2", "T_dil1","T_dil2","T_dil1", "T_dil2","T_dil1","T_dil2"), + readout = c(31993, 31925, 49336, 49003 , 68332, 68034 , 33915, 34007, 52131, 52179 , 70617, 69962), + isRef=c(rep(1,6), rep(0,6)), + isSample = c(rep(0,6), rep(1,6))) +Lim <- c(rep(0,8),70,130) # only Lim 9 and 10 relevant +PureErrF <- TRUE + + +ANOVAlintests(ro_new=dat, circles=CIRC, Lim=Lim, PureErrF) + +# [[1]] +# test test_results estimate Source df SumSquares MS F.value p.value +# 1 F-test on sign. of regression 0 0.00000 Treatment 3 3889.33189 1296.44396 127.75497 0 +# 2 F_test on non-lin 0 0.58832 Preparation 1 969.90761 969.90761 31.85908 0 +# 3 F-test on R^2 A 0 0.25260 Regression 1 2903.63083 2903.63083 286.13136 0 +# 4 F_test on R^2 B 0 0.85290 Non-parallelism 1 4.53646 4.53646 0.52154 0.484 +# 5 F-test on slope A 0 0.00000 Resid Error 14 133.03173 9.50227 +# 6 F-test on slope B 0 0.00000 Non-linearity 2 11.25699 5.62850 0.55465 0.58832 +# 7 F-test on non-parallelism 0 0.48400 Pure error 12 121.77473 10.14789 +# 8 F-test on preparation 0 0.00000 Total 17 4011.10663 235.94745 +# +# [[2]] +# parameter Estimate Std. Error t value Pr(>|t|) +# (Intercept) intercept REF 131.223810 6.039254 21.7284793 5.288214e-11 +# log_dose slope REF 22.342857 1.857866 12.0260887 4.720297e-08 +# isSample intercepts diff. -4.977419 9.167857 -0.5429207 5.971246e-01 +# log_dose:isSample slope difference -1.698577 2.540472 -0.6686068 5.164021e-01 +# +# [[3]] +# log_dose:isSample log_dose:isSample +# -1.698577 -7.556917 4.159763 +# +# [[4]] +# parameter Estimate Std. Error t value Pr(>|t|) +# (Intercept) intercept REF 128.316878 4.191623 30.6126966 1.0000000 +# log_dose slope REF 21.434441 1.267154 16.9154178 1.0000000 +# isSample intercepts diff. 1.037479 1.765952 0.5874899 0.5951663 + + + + +ro_new <- data.frame(R_dil1 =c(10221, 18258, 31993, 49336, 68332, 83527, 95584, 102229), +R_dil2=c( 10136, 19078, 31925, 49003, 68034, 83776, 95495, 101608), +T_dil1=c( 10830, 19891, 33915, 52131, 70617, 85784, 95937, 102791), +T_dil2=c( 11169, 20153, 34007, 52179, 69962, 85543, 96439, 102655), +log_dose=c( -1.2029, -1.89712, -2.590267, -3.2834, -3.97656, -4.66917, -5.362323, -6.05334)) +PureErrF <- TRUE + + +pot4plFUNC(ro_new, PureErrF) +# pot_est "restricted" "0.911577062498655" "0.902232125904588" "0.921018789971077" +#potU_est "unrestricted" "0.900218901786783" "0.874169901564974" "0.92704412458425" +#r "ln-transformed restr" "0.914213577701921" "0.902850729973349" "0.925719432800612" +#r "ln-transformed unrestr" "0.990574387419529" "0.75948098671897" "1.29198443959817" + +xs=2; xt=3.2; se_xt=0.34;se_xs=0.23; DFs=32-16 + + +ParamCI_F(xt,xs,se_xt, se_xs, CoVar,DFs, Conf=0.975) +# [1] 1.14808 2.22981 + + +dat <- data.frame(REF1=c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2=c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497), +REF3=c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665), +TEST1=c(1405, 1523, 1502, 1474, 2383, 3221, 3589, 3445), TEST2=c(1420, 1516, 1544, 1512, 2226, 3219, 3327, 3591), +TEST3=c(1399, 1376, 1588, 1475, 2148, 3083, 2942, 3466), log_dose=c(5.01,3.401,2.708,2.015,1.32176,0.62861,-0.0645385,-1.6739764)) +Lim <- c(-1, 1, 0.005 , 2, 0.5, 2, 0.5,2,75,133, 75, 133) +PureErrF <- FALSE +#' +tests_FUNC(dat, Lim,PureErrF) +# test test_results estimate lower_limit upper_limit lower_CI upper_CI +# 1 F-test on sign. of regression* 0 0.00000 - - 166.34405 153.405644181258 +# 2 EQ test on lower asymptotes difference 1 -76.27934 -1 1 -270.281 117.723 +# 3 EQ test ratio of lower asymptotes 0 0.95010 0.005 2 0.84509 1.06817 +# 4 EQ test ratio of Hill slopes 1 0.92581 0.5 2 0.58023 2.01077 +# 5 EQ test ratio of upper asymptotes 0 0.95010 0.5 2 0.92865 1.03996 +# 6 F-test on non-linearity* 0 0.12180 - - - - +# 7 EQ test ratio of asymptote difference 1 1.00851 75 133 0.8744 1.16319 +# 8 geom. rel. CI restr. model 1 211.16000 75 133 1.81 2.46 +# 9 geom. rel. CI unrestr. model 1 197.72000 75 133 1.63 2.4 + +tests_FUNC(dat, Lim,PureErrFlag=T) +# test test_results estimate lower_limit upper_limit lower_CI upper_CI +# 1 F-test on sign. of regression* 0 0.00000 - - 166.34405 153.405644181258 +# 2 EQ test on lower asymptotes difference 1 -76.27935 -1 1 -270.281 117.723 +# 3 EQ test ratio of lower asymptotes 0 0.95010 0.005 2 0.85388 1.05716 +# 4 EQ test ratio of Hill slopes 0 0.92581 0.5 2 0.61299 1.90329 +# 5 EQ test ratio of upper asymptotes 0 0.95010 0.5 2 0.93331 1.03477 +# 6 F-test on non-linearity* 0 0.06397 - - - - +# 7 EQ test ratio of asymptote difference 1 1.00851 75 133 0.8855 1.14861 +# 8 geom. rel. CI restr. model 1 211.16000 75 133 1.81 2.46 +# 9 geom. rel. CI unrestr. model 1 197.72000 75 133 1.63 2.4 + + + +ro_new <- data.frame(REF1=c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2=c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497), +REF3=c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665), +TEST1=c(1405, 1523, 1502, 1474, 2383, 3221, 3589, 3445), TEST2=c(1420, 1516, 1544, 1512, 2226, 3219, 3327, 3591), +TEST3=c(1399, 1376, 1588, 1475, 2148, 3083, 2942, 3466), log_dose=c(5.01,3.401,2.708,2.015,1.32176,0.62861,-0.0645385,-1.6739764)) + + +ANOVA4plUnresfunc(ro_new) +# Source DF SumSquares MeanSquares F.value p_value +# 1 Treatment 7 36070704.40 5152957.772 179.14139 0 +# 2 Preparation 1 1131295.02 1131295.02083 39.32921 0 +# 3 Regression 3 36031469.73 12023568.134 417.54192 0 +# 4 Non-Parallelism 3 39234.65 13078.216 0.45466 0.715 +# 5 Residual Error 40 1150590.09 28764.75235 +# 6 Non-linearity 8 397524.76 49690.59506 67.56807 0 +# 7 Pure Error 32 23533.29 735.41536 +# 8 Total 47 37221294.48 + + +ro_new <- data.frame(REF1=c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2=c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497), +REF3=c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665), +TEST1=c(1405, 1523, 1502, 1474, 2383, 3221, 3589, 3445), TEST2=c(1420, 1516, 1544, 1512, 2226, 3219, 3327, 3591), +TEST3=c(1399, 1376, 1588, 1475, 2148, 3083, 2942, 3466), log_dose=c(5.01,3.401,2.708,2.015,1.32176,0.62861,-0.0645385,-1.6739764)) + + +perConcTab(ro_new, noDilSeries=3) +# 5.01 3.401 2.708 2.015 1.32176 0.62861 -0.06454 -1.67398 +# REF1 1547.0000000 1620.000000 1644.000000 2504.000000 3426.000000 3512.000000 3401.000000 3787.000000 +# REF2 1492.0000000 1536.000000 1384.000000 2286.000000 3046.000000 3479.000000 3516.000000 3497.000000 +# REF3 1468.0000000 1827.000000 1558.000000 2252.000000 3002.000000 3349.000000 2945.000000 3665.000000 +# avs 1502.3333333 1661.000000 1528.666667 2347.333333 3158.000000 3446.666667 3287.333333 3649.666667 +# sds 40.5010288 149.769823 132.458799 136.738193 233.135154 86.176176 301.993929 145.606776 +# cv 2.6958750 9.016847 8.664989 5.825257 7.382367 2.500276 9.186593 3.989591 +# TEST1 1405.0000000 1523.000000 1502.000000 1474.000000 2383.000000 3221.000000 3589.000000 3445.000000 +# TEST2 1420.0000000 1516.000000 1544.000000 1512.000000 2226.000000 3219.000000 3327.000000 3591.000000 +# TEST3 1399.0000000 1376.000000 1588.000000 1475.000000 2148.000000 3083.000000 2942.000000 3466.000000 +# avs_test 1408.0000000 1471.666667 1544.666667 1487.000000 2252.333333 3174.333333 3286.000000 3500.666667 +# sds_test 10.8166538 82.923660 43.003876 21.656408 119.692662 79.103308 325.442775 78.932461 +# cv_test 0.7682283 5.634677 2.784023 1.456383 5.314163 2.491966 9.903919 2.254784 + +x <- 1; Div <- 3;N <- 0; res <- c(); noDil <- 7 + +divFUN(x,Div,N,res,noDil) +# [1] 0.3333333333 0.1111111111 0.0370370370 0.0123456790 0.0041152263 0.0013717421 0.0004572474 diff --git a/app.R b/app.R index 45f14b1..1f72925 100644 --- a/app.R +++ b/app.R @@ -25,16 +25,143 @@ Dat <- reactiveValues() REP <- reactiveValues() #### functions ---- -dilFUN2 <- function(cs_,dils,Faktor) { - av <- cs_ - dils_av <- dils_av - dils_avsc <- dils_av*Faktor - dils2 <- dils_avsc+av - dilfactors <- 1/exp(dils2-lag(dils2)) - return(dilfactors) + + + +# dilFUN2 <- function(cs_,dils,Faktor) { +# av <- cs_ +# dils_av <- dils +# dils_avsc <- dils_av*Faktor +# dils2 <- dils_avsc+av +# dilfactors <- 1/exp(dils2-lag(dils2)) +# return(dilfactors) +# } + +#' Levenberg Marquard fit of 4 pl +#' +#' Returns list of following: summary of restricted and unrestricted model fits of the 4pl function, +#' potency estimates and respective CIs of restricted and unrestricted models, and the predictions thereof. +#' +#' @param ro_new The dilution series readouts for REF and TEST + 1 column log(loncentration) or concentration. +#' @param TransFlag A boolean if the plot should be with natural log as readout (TRUE) or not (FALSE). +#' @returns Returns list of following: summary of restricted and unrestricted model fits of the 4pl function, +#' potency estimates and respective CIs of restricted and unrestricted models, and the predictions thereof. +#' @export +#' @examples +#' dat <- data.frame(REF1=c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2=c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497), +#' REF3=c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665), +#' TEST1=c(1405, 1523, 1502, 1474, 2383, 3221, 3589, 3445), TEST2=c(1420, 1516, 1544, 1512, 2226, 3219, 3327, 3591), +#' TEST3=c(1399, 1376, 1588, 1475, 2148, 3083, 2942, 3466), log_dose=c(5.01,3.401,2.708,2.015,1.32176,0.62861,-0.0645385,-1.6739764)) +#' TransF<- FALSE +#' Dat<- list() +#' te <- Fitting_FUNC(dat,TransF) +#' print(te) + + +Fitting_FUNC <- function(ro_new, TransFlag=F) { + CORro <- cor(ro_new[,1],ro_new[,ncol(ro_new)]) + #browser() + all_l <- melt(data.frame(ro_new), id.vars="log_dose", variable.name="replname", value.name = "readout") + isRef <- rep(c(1,0),1,each=nrow(all_l)/2) + isSample <- rep(c(0,1),1,each=nrow(all_l)/2) + all_l2 <- cbind(all_l, isRef, isSample) +#browser() + if (CORro<0) SLOPE <- -1 else SLOPE <- 1 + if (!TransFlag) { + startlist <- list(a=min(ro_new[,2]), b=SLOPE, d=max(ro_new[,2]), cs=mean(all_l$log_dose),r=0) + + mr <- gsl_nls(fn = readout ~ a+(d-a)/(1+exp(b*((cs-r*isSample)-log_dose))), + data=all_l2, + start=startlist,#trace=T, + control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) + s_mr <- tryCatch({ + s_mr <- summary(mr) + }, + error = function(err) { + s_mr <- NULL + }) + } else { + startlist <- list(a=log(min(ro_new[,2])), b=SLOPE, d=log(max(ro_new[,2])), cs=mean(all_l$log_dose),r=0) + mrT <- gsl_nls(fn = log(readout) ~ a+(d-a)/(1+exp(b*((cs-r*isSample)-log_dose))), + data=all_l2, + start=startlist, + control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) + s_mr <- summary(mrT) + } + + if (!TransFlag) { + startlistmu <- list(as=min(ro_new[,2]), bs=SLOPE, ds=max(ro_new[,2]), cs=mean(all_l$log_dose), + at=min(ro_new[,2]), bt=SLOPE, dt=max(ro_new[,2]), r=0) + tryCatch({ + mu <- gsl_nls(fn = readout ~ as*isRef + at*isSample + (ds*isRef + dt*isSample - as*isRef - at*isSample)/ + (1+isRef*exp(bs*(cs - log_dose)) + isSample*exp(bt*(cs-r*isSample-log_dose))), + data=all_l, + start=startlistmu, + control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) + }, + error = function(msg){ + return(0) }) + + Sum_u <- tryCatch({ summary(mu) }, + error=function(msg){ + return(0) }) + } else { + startlistmu <- list(as=log(min(ro_new[,2])), bs=SLOPE, ds=log(max(ro_new[,2])), cs=mean(all_l$log_dose), + at=log(min(ro_new[,2])), bt=SLOPE, dt=log(max(ro_new[,2])), r=0) + tryCatch({ + muT <- gsl_nls(fn = log(readout) ~ as*isRef + at*isSample + (ds*isRef + dt*isSample - as*isRef - at*isSample)/ + (1+isRef*exp(bs*(cs - log_dose)) + isSample*exp(bt*(cs-r*isSample-log_dose))), + data=all_l, + start=startlistmu, + control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) + }, + error = function(msg){ + return(0) }) + + Sum_u <- tryCatch({ summary(muT) }, + error=function(msg){ + return(0) }) + } + if (!TransFlag) { + pot_est <- exp(confintd(mr, "r", method="asymptotic")) + potU_est <- exp(confintd(mu, "r", method="asymptotic")) + PRED <- predict(mr) + PREDu <- predict(mu) + } else { + pot_est <- exp(confintd(mrT, "r", method="asymptotic")) + potU_est <- exp(confintd(muT, "r", method="asymptotic")) + PRED <- predict(mrT) + PREDu <- predict(muT) + } + return(list(s_mr, Sum_u, pot_est, potU_est, PRED, PREDu)) } -plot_f <- function(dat, sigmoid,det_sig, TransFlag=F) { + +#' Plot sigmoidal curve +#' +#' Returns the final plots of the 4pl function as sigmoidal lines, and the single readouts as scatter, with REF in blue and TEST in red. +#' Two plots are returned in a grid object: the unrestricted model and the restricted model. +#' +#' @param dat The dilution series readouts for REF and TEST + 1 column log(loncentration) or concentration. +#' @param sigmoid The parameters of the 4pl fit of unrestricted model from EXCEL input. +#' @param det_sig The parameters of the 4pl fit of unrestricted model from the metadata input. +#' @param TransFlag A boolean if the plot should be with natural log as readout (TRUE) or not (FALSE). +#' @returns A grid object either of the original scale or the natural log of the readouts. +#' @export +#' @examples +#' dat <- data.frame(REF1=c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2=c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497), +#' REF3=c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665), +#' TEST1=c(1405, 1523, 1502, 1474, 2383, 3221, 3589, 3445), TEST2=c(1420, 1516, 1544, 1512, 2226, 3219, 3327, 3591), +#' TEST3=c(1399, 1376, 1588, 1475, 2148, 3083, 2942, 3466), log_dose=c(5.01,3.401,2.708,2.015,1.32176,0.62861,-0.0645385,-1.6739764)) +#' sigmoid <- c(0.7163324, 0.5636804, 10.6156340 , 9.9784160, -0.7504673, -0.7108692, -3.5788141, -0.6662962) +#' det_sig <- FALSE +#' TransF<- FALSE +#' Dat<- list() +#' p <- plot_f(dat,sigmoid,det_sig,TransF) +#' print(p) + + +plot_f <- function(dat, TransFlag=F) { #sigmoid,det_sig, CORdat <- cor(dat[,1],dat[,ncol(dat)]) #browser() all_l <- melt(data.frame(dat), id.vars="log_dose", variable.name="replname", value.name = "readout") @@ -42,27 +169,8 @@ plot_f <- function(dat, sigmoid,det_sig, TransFlag=F) { isSample <- rep(c(0,1),1,each=nrow(all_l)/2) all_l2 <- cbind(all_l, isRef, isSample) #browser() - if(is.null(det_sig)) { - - startlist <- list(a=sigmoid[1], b=sigmoid[5],cs=sigmoid[7], - d=sigmoid[3],r=sigmoid[8]) - - } else { - startlist <- list(a=det_sig[5], b=det_sig[1],cs=det_sig[7], - d=det_sig[3],r=det_sig[7] - det_sig[8]) - } - - mr <- gsl_nls(fn = readout ~ a+(d-a)/(1+exp(b*((cs-r*isSample)-log_dose))), - data=all_l, - start=startlist, - control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) - s_mr <- tryCatch({ - s_mr <- summary(mr) - }, - error = function(err) { - s_mr <- NULL - }) - + MODLS <- Fitting_FUNC(dat, TransFlag=F) + s_mr <- MODLS[[1]] a <- s_mr$coefficients["a",1] b <- s_mr$coefficients["b",1] cs <- s_mr$coefficients["cs",1] @@ -74,13 +182,22 @@ plot_f <- function(dat, sigmoid,det_sig, TransFlag=F) { SAMPLE <- a+(d-a)/(1+exp(b*((cs-r)-seq_x))) REF <- a+(d-a)/(1+exp(b*((cs)-seq_x))) - if (is.null(det_sig)) { - SAMPLEtrue <- sigmoid[2] + (sigmoid[4] -sigmoid[2])/(1+exp(sigmoid[6]*((sigmoid[7]-sigmoid[8]-seq_x)))) - REFtrue <- sigmoid[1] + (sigmoid[3] -sigmoid[1])/(1+exp(sigmoid[5]*((sigmoid[7]-seq_x)))) - } else { - SAMPLEtrue <- det_sig[4] + (det_sig[6] -det_sig[4])/(1+exp(det_sig[2]*(det_sig[8]-seq_x))) - REFtrue <- det_sig[3] + (det_sig[5] -det_sig[3])/(1+exp(det_sig[1]*(det_sig[7]-seq_x))) - } + s_mu <- MODLS[[2]] + + as <- s_mu$coefficients["as",1] + bs <- s_mu$coefficients["bs",1] + cs <- s_mu$coefficients["cs",1] + ds <- s_mu$coefficients["ds",1] + at <- s_mu$coefficients["at",1] + bt <- s_mu$coefficients["bt",1] + ct <- s_mu$coefficients["cs",1]-s_mu$coefficients["r",1] + dt <- s_mu$coefficients["dt",1] + r <- s_mu$coefficients["r",1] + + + SAMPLEtrue <- at + (dt -at)/(1+exp(bt*((cs-r-seq_x)))) + REFtrue <- as + (ds -as)/(1+exp(bs*((cs-seq_x)))) + #browser() pl_df <- cbind(seq_x, SAMPLE, REF, SAMPLEtrue, REFtrue) all_l2$readout[all_l2$readout < 0] <- 0.01 @@ -88,10 +205,10 @@ plot_f <- function(dat, sigmoid,det_sig, TransFlag=F) { slopeEC50 <- b*(d-a)/4 Intercept <- a+(d-a)/2-b*(d-a)/4*cs #browser() - Xbendl3 <- cs-(1.31696/b) - Xbendu3 <- cs+(1.31696/b) - XbendlT <- cs-r-(1.31696/b) - XbenduT <- cs-r+(1.31696/b) + Xbendl3 <- cs-(1.5434/b) + Xbendu3 <- cs+(1.5434/b) + XbendlT <- cs-r-(1.5434/b) + XbenduT <- cs-r+(1.5434/b) XasymplS <- cs-(3/b) XasympuS <- cs+(3/b) XasymplT <- cs-r-(3/b) @@ -104,7 +221,7 @@ plot_f <- function(dat, sigmoid,det_sig, TransFlag=F) { Dat$cfordils <- cs p <- ggplot(all_l2, aes(x=log_dose, y=readout, color=factor(isRef))) + - geom_point(shape=factor(isRef), alpha=0.8) + + geom_point(shape=factor(isRef), size=3,alpha=0.8) + labs(title = paste("restricted 4pl"), color="product") + scale_color_manual(labels=c("test","reference"), values=c("#C2173F","#4545BA")) + @@ -131,26 +248,23 @@ plot_f <- function(dat, sigmoid,det_sig, TransFlag=F) { # transformed plots p_rt <- ggplot(all_l2, aes(x=log_dose, y=readouttrans, color=factor(isRef))) + - geom_point(shape=factor(isRef), alpha=0.8) + + geom_point(shape=factor(isRef), size=3,alpha=0.8) + labs(title = paste("restricted transformed 4pl"), color="product") + scale_color_manual(labels=c("test","reference"), values=c("#C2173F","#4545BA")) + theme_bw() - mrt <- gsl_nls(fn = readouttrans ~ a+(d-a)/(1+exp(b*((cs-r*isSample)-log_dose))), - data=all_l2, - start=startlist, - control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) - s_mrt <- summary(mrt) + MODLStrans <- Fitting_FUNC(dat,TransFlag = TRUE) + s_mrt <- MODLStrans[[1]] a_trans <- s_mrt$coefficients["a",1] b_trans <- s_mrt$coefficients["b",1] cs_trans <- s_mrt$coefficients["cs",1] d_trans <- s_mrt$coefficients["d",1] r_trans <- s_mrt$coefficients["r",1] - XbendlTrans <- cs_trans-(1.31696/b_trans) - XbenduTrans <- cs_trans+(1.31696/b_trans) - XbendlTransT <- cs_trans-r_trans-(1.31696/b_trans) - XbenduTransT <- cs_trans-r_trans+(1.31696/b_trans) + XbendlTrans <- cs_trans-(1.5434/b_trans) + XbenduTrans <- cs_trans+(1.5434/b_trans) + XbendlTransT <- cs_trans-r_trans-(1.5434/b_trans) + XbenduTransT <- cs_trans-r_trans+(1.5434/b_trans) bendpointsTRANS <- c(bendREF_lower = round(XbendlTrans,3), bendREF_upper=round(XbenduTrans,3), bendSAMPLE_lower = round(XbendlTransT,3), bendSAMPLE_upper=round(XbenduTransT,3)) Dat$bendpointsTRANS <- bendpointsTRANS @@ -167,30 +281,7 @@ plot_f <- function(dat, sigmoid,det_sig, TransFlag=F) { theme(legend.position = "none", axis.text=element_text(size=14)) #browser() - if(is.null(det_sig)) { - - startlistmu <- list(as=sigmoid[1], bs=sigmoid[5],cs=sigmoid[7], - ds=sigmoid[3], at=sigmoid[2], bt=sigmoid[6], - dt=sigmoid[4],r=sigmoid[8]) - - } else { - startlistmu <- list(as=det_sig[5], bs=det_sig[1],cs=det_sig[7], - ds=det_sig[3],at=det_sig[6], bt=det_sig[2], - dt=det_sig[4],r=det_sig[7] - det_sig[8]) - } - tryCatch({ - mu <- gsl_nls(fn = readout ~ as*isRef + at*isSample + (ds*isRef + dt*isSample - as*isRef - at*isSample)/ - (1+isRef*exp(bs*(cs - log_dose)) + isSample*exp(bt*(cs-r*isSample-log_dose))), - data=all_l, - start=startlistmu, - control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) - }, - error = function(msg){ - return(0) }) - - Sum_u <- tryCatch({ summary(mu) }, - error=function(msg){ - return(0) }) + Sum_u <- MODLS[[2]] #browser() #if (is.null(det_sig)) { ast <- Sum_u$coefficients["as",1] @@ -216,9 +307,9 @@ plot_f <- function(dat, sigmoid,det_sig, TransFlag=F) { pl_df2 <- cbind(seq_x, SAMPLEu, REFu) #browser() pu <- ggplot(all_l2, aes(x=log_dose, y=readout, color=factor(isRef))) + - geom_point() + + geom_point(size=3) + labs(title="unrestricted 4_pl-Model", color="product") + - scale_color_manual(labels = c("test","reference"), values=c("#C2173F","#4545BA")) + + scale_color_manual(labels = c("test","reference"), values=c("#C2173F88","#4545BA88")) + theme_bw() pu2 <- pu + geom_line(data=as.data.frame(pl_df2), aes(x=seq_x, y=SAMPLEu), color="#C2173F", inherit.aes = F) + @@ -227,26 +318,26 @@ plot_f <- function(dat, sigmoid,det_sig, TransFlag=F) { show.legend = F) pu2_ <- pu2 + theme(legend.position = "none", axis.text = element_text(size=14)) + putrans <- ggplot(all_l2, aes(x=log_dose, y=readouttrans, color=factor(isRef))) + - geom_point() + + geom_point( size=3) + labs(title="unrestricted transformed 4_pl-Model", color="product") + - scale_color_manual(labels = c("test","reference"), values=c("#C2173F","#4545BA")) + + scale_color_manual(labels = c("test","reference"), values=c("#C2173FAA","#4545BAAA")) + theme_bw() - unrestr_trans <- drm(readouttrans ~ exp(log_dose), isSample, data=all_l2, fct=LL.4(), - pmodels=data.frame(isSample, isSample,isSample,isSample)) - Sum_ut <- summary(unrestr_trans) - ast_t <- Sum_ut$coefficients[3,1] - ate_t <- Sum_ut$coefficients[4,1] - bst_t <- Sum_ut$coefficients[1,1] - bte_t <- Sum_ut$coefficients[2,1] - cst_t <- log(Sum_ut$coefficients[7,1]) - cte_t <- log(Sum_ut$coefficients[8,1]) - dst_t <- Sum_ut$coefficients[5,1] - dte_t <- Sum_ut$coefficients[6,1] + Sum_ut <- MODLStrans[[2]] - REFu_trans <- ast_t + (dst_t-ast_t)/(1+exp(bst_t*(seq_x-cst_t))) - SAMPLEu_trans <- ate_t + (dte_t-ate_t)/(1+exp(bte_t*(seq_x-cte_t))) + ast_t <- Sum_ut$coefficients["as",1] + ate_t <- Sum_ut$coefficients["at",1] + bst_t <- Sum_ut$coefficients["bs",1] + bte_t <- Sum_ut$coefficients["bt",1] + cst_t <- Sum_ut$coefficients["cs",1] + cte_t <- Sum_ut$coefficients["cs",1]-Sum_ut$coefficients["r",1] + dst_t <- Sum_ut$coefficients["ds",1] + dte_t <- Sum_ut$coefficients["dt",1] + + REFu_trans <- ast_t + (dst_t-ast_t)/(1+exp(bst_t*(cst_t-seq_x))) + SAMPLEu_trans <- ate_t + (dte_t-ate_t)/(1+exp(bte_t*(cte_t-seq_x))) pl_df2u_t <- cbind(seq_x, SAMPLEu_trans, REFu_trans) pu2_t <- putrans + geom_line(data=as.data.frame(pl_df2u_t), aes(x=seq_x, y=SAMPLEu_trans), @@ -258,6 +349,27 @@ plot_f <- function(dat, sigmoid,det_sig, TransFlag=F) { if (TransFlag) grid.arrange(p_rt2,pu3_t, nrow=1) else grid.arrange(p2,pu2_, nrow=1) } +#' Simulate a well-plate readout +#' +#' Returns a possible well-plate result based on the metadata provided and a variability. +#' Homo- and heteroscedasticity can be simulated. +#' +#' @param as, bs, cs, ds, at, bt, dt,r,ct, gt, gs are the parameter of the 5pl logistic regression. +#' @param sd_fac The standard deviation with units in % of upper-lower asymptote. +#' @param log_conc The natural log of the concentrations. +#' @param heteroNoise Boolean if heteroscedstic noise should be added. +#' @param noDilSeries A number >1 indicating how many dilution series per product should be simulated. +#' @param noDils A number >7 indicating how many dilutions steps per product should be simulated. +#' @returns A data-frame with readouts and natural log of concentrations. +#' @export +#' @examples +#' as=3; bs=1; cs=-4; ds=10; at=3; bt=1; dt=10;r=0.0001;ct=cs-r;sd_fac=0.1; gt=1; gs=1; +#' lnConc=c(1,0,-1,-2,-3,-4,-5,-6) +#' heteroNoise <- FALSE +#' noDilS <- 3 +#' noD <- 8 +#' Calc_DilRes(as,bs,cs,ds,at,bt,dt,r,ct,sd_fac,gt,gs,log_conc=lnConc, heteroNoise, noDilS, noD) + Calc_DilRes <- function(as=3, bs=1, cs=-4, ds=10, at=3, bt=1, dt=10,r=0.0001,ct=cs-r, sd_fac=0.1, gt=1, gs=1, log_conc, heteroNoise=FALSE, noDilSeries, noDils) { @@ -286,6 +398,31 @@ Calc_DilRes <- function(as=3, bs=1, cs=-4, ds=10, at=3, bt=1, dt=10,r=0.0001,ct return(ro_jit2) } +#' Calculate the potency of linear PLA +#' +#' The delta Method is used for calculating the potency using a restricted linear model. +#' Absolute and relative CIs are calculated. The model RMSE or the Pure Error can be used. +#' +#' @param circles Data frame with the 3 consecutive concentrations and the respective readouts. +#' @param Lim The limits to check if the relative CI is within the limits. +#' @param PureErrFlag Boolean if Pure Error should be used. +#' @returns A data-frame with potency estimate, absolute CIs, test result, relative CIs. +#' @export +#' @examples +#' CIRC <- data.frame(log_dose = c(-2.5,-2.5,-2.5, -3.2,-3.2,-3.2,-3.9,-3.9,-3.9, +#' -3.2,-3.2,-3.2,-3.9,-3.9,-3.9,-4.7,-4.7,-4.7), +#' replname= c("R_dil1","R_dil1","R_dil1", "R_dil2","R_dil2","R_dil2", "R_dil3","R_dil3","R_dil3", +#' "T_dil1","T_dil1","T_dil1", "T_dil2","T_dil2","T_dil2", "T_dil3","T_dil3","T_dil3"), +#' readout = c(72.1,75.8,76.04,59.8,61,62.7,43.6,45,41.5,53.5,62.2,65.9,48.3,43.8,43.14,28.17,29.2,31.2), +#' isRef=c(rep(1,9), rep(0,9)), +#' isSample = c(rep(0,9), rep(1,9))) +#' Lim <- c(rep(0,8),70,130) # only Lim 9 and 10 relevant +#' PureErrF <- TRUE +#' +#' +#' LinPotTab(circles=CIRC, Lim, PureErrF) + + LinPotTab <- function(circles, Lim, PureErrFlag) { circ_ABl <- circles circ_Al <- circ_ABl[circ_ABl$isSample ==1,] @@ -333,6 +470,37 @@ LinPotTab <- function(circles, Lim, PureErrFlag) { return(pottab) } + +#' Calculate the ANOVA of linear PLA and restricted and unrestr. model summaries +#' +#' The delta Method is used for calculating the potency using a restricted linear model. +#' Absolute and relative CIs are calculated. The model RMSE or the Pure Error can be used. +#' +#' @param circles Data frame with the 3 consecutive concentrations and the respective readouts. +#' @param Lim The limits to check if the relative CI is within the limits. +#' @param PureErrFlag Boolean if Pure Error should be used. +#' @returns A data-frame with 1) data.frame with F-tests ANOVA and test results +#' 2) summary of unrestricted linear model 3) Slope difference +CIs +#' 4) summary of restricted linear model. +#' @export +#' @examples +#' dat <- data.frame(R_dil1 =c(10221, 18258, 31993, 49336, 68332, 83527, 95584, 102229), +#' R_dil2=c( 10136, 19078, 31925, 49003, 68034, 83776, 95495, 101608), +#' T_dil1=c( 10830, 19891, 33915, 52131, 70617, 85784, 95937, 102791), +#' T_dil2=c( 11169, 20153, 34007, 52179, 69962, 85543, 96439, 102655), +#' log_dose=c( -1.2029, -1.89712, -2.590267, -3.2834, -3.97656, -4.66917, -5.362323, -6.05334)) +#' CIRC <- data.frame(log_dose=c( -2.590267, -2.590267, -3.2834 , -3.2834, -3.97656, -3.97656, -2.590267, -2.590267,-3.2834, -3.2834, -3.97656, -3.97656), +#' replname= c("R_dil1","R_dil2","R_dil1", "R_dil2","R_dil1","R_dil2", "T_dil1","T_dil2","T_dil1", "T_dil2","T_dil1","T_dil2"), +#' readout = c(31993, 31925, 49336, 49003 , 68332, 68034 , 33915, 34007, 52131, 52179 , 70617, 69962), +#' isRef=c(rep(1,6), rep(0,6)), +#' isSample = c(rep(0,6), rep(1,6))) +#' Lim <- c(rep(0,8),70,130) # only Lim 9 and 10 relevant +#' PureErrF <- TRUE +#' +#' +#' ANOVAlintests(ro_new, circles, Lim, PureErrF) + + ANOVAlintests <- function(ro_new, circles, Lim, PureErrFlag) { all_l <- melt(data.frame(ro_new), id.vars="log_dose", variable.name = "replname", value.name = "readout") isRef <- rep(c(1,0),1,each=nrow(all_l)/2) @@ -342,7 +510,7 @@ ANOVAlintests <- function(ro_new, circles, Lim, PureErrFlag) { all_l$Conc <- exp(all_l$log_dose) all_lA <- all_l[all_l$isSample == 1,] # TEST all_lB <- all_l[all_l$isSample == 0,] # REF - + #browser() circ_ABl <- circles circ_Al <- circ_ABl[circ_ABl$isSample ==1,] circ_Bl <- circ_ABl[circ_ABl$isSample ==0,] @@ -352,7 +520,7 @@ ANOVAlintests <- function(ro_new, circles, Lim, PureErrFlag) { modABu <- lm(readout ~ log_dose + isSample + log_dose*isSample, circ_ABl) modA <- lm(readout ~ log_dose + isSample, circ_Al) modB <- lm(readout ~ log_dose + isSample, circ_Bl) - + log_pot_delta <- deltaMethod(modAB, "isSample/log_dose") CI_log_low <- log_pot_delta$Estimate - qt(0.975, df.residual(modAB))*log_pot_delta$SE CI_log_up <- log_pot_delta$Estimate + qt(0.975, df.residual(modAB))*log_pot_delta$SE @@ -501,6 +669,30 @@ ANOVAlintests <- function(ro_new, circles, Lim, PureErrFlag) { return(RET) } +#' Calculates the potency of 4PL PLA of all models model +#' +#' The gradient method is used for calculating the potency for a restricted model, an unrestricteed model, +#' and the log-transformations thereof. +#' Absolute and relative CIs are calculated. The model RMSE or the Pure Error can be used. +#' +#' @param ro_new The dilution series readouts for REF and TEST + 1 column log(concentration). +#' @param PureErrFlag Boolean if Pure Error should be used. +#' @returns A data-frame with 1) data.frame with F-tests ANOVA and test results +#' 2) summary of unrestricted linear model 3) Slope difference +CIs +#' 4) summary of restricted linear model. +#' @export +#' @examples +#' ro_new <- data.frame(R_dil1 =c(10221, 18258, 31993, 49336, 68332, 83527, 95584, 102229), +#' R_dil2=c( 10136, 19078, 31925, 49003, 68034, 83776, 95495, 101608), +#' T_dil1=c( 10830, 19891, 33915, 52131, 70617, 85784, 95937, 102791), +#' T_dil2=c( 11169, 20153, 34007, 52179, 69962, 85543, 96439, 102655), +#' log_dose=c( -1.2029, -1.89712, -2.590267, -3.2834, -3.97656, -4.66917, -5.362323, -6.05334)) +#' PureErrF <- TRUE +#' +#' +#' pot4plFUNC(ro_new, PureErrF) + + pot4plFUNC <- function(ro_new, PureErrFlag) { all_l <- melt(data.frame(ro_new), id.vars="log_dose", variable.name="replname", value.name = "readout") isRef <- rep(c(1,0),1,each=nrow(all_l)/2) @@ -510,116 +702,94 @@ pot4plFUNC <- function(ro_new, PureErrFlag) { all_l$Conc <- exp(all_l$log_dose) all_l$readout[all_l$readout < 0] <- 0.01 all_l$readouttrans <- log(all_l$readout) - +#browser() CORdat <- cor(ro_new[,1],ro_new[,ncol(ro_new)]) if (CORdat<0) SLOPE <- -1 else SLOPE <- 1 - - startlist <- list(a=min(ro_new[,2]), b=SLOPE, d=max(ro_new[,2]), cs=mean(all_l$log_dose),r=0) - tryCatch({ - mr <- gsl_nls(fn = readout ~ a+(d-a)/(1+exp(b*(cs-r*isSample-log_dose))), - data=all_l, - start=startlist, - control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) - }, - error = function(err) { - err$message - }) - - startlistmu <- list(as=max(ro_new[,2]), bs=SLOPE, ds=min(ro_new[,2]), cs=mean(all_l$log_dose), - at=max(ro_new[,2]), bt=SLOPE, dt=min(ro_new[,2]), r=0) - - tryCatch({ - mu <- gsl_nls(fn = readout ~ as*isRef + at*isSample + (ds*isRef + dt*isSample - as*isRef - at*isSample)/ - (1+isRef*exp(bs*(cs - log_dose)) + isSample*exp(bt*(cs-r*isSample-log_dose))), - data=all_l, - start=startlistmu, - control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) - }, - error = function(err) { - err$message - }) - + # + FITs <- Fitting_FUNC(ro_new, TransFlag = F) if (!PureErrFlag) { - pot_est <- exp(confintd(mr, "r", method="asymptotic")) - potU_est <- exp(confintd(mu, "r", method="asymptotic")) + pot_est <- FITs[[3]] + potU_est <- FITs[[4]] } else { FitAnova <- anova(lm(readout ~ factor(Conc)*isSample, all_l)) meanPureErr <- FitAnova[4,3] - SU_mr <- tryCatch({ - SU_mr <- summary(mr) - }, error = function(msg) { - return() - }) - - browser() - if (length(SU_mr)>1) { - s_mr <- SU_mr$coefficients - } else { SU_mr <- rep(NA,5) } - - VCOV <- vcov(mr) - V_V <- VCOV/SU_mr$sigma^2 + # SU_mr <- tryCatch({ + # SU_mr <- summary(mr) + # }, error = function(msg) { + # return() + # }) + SU_mr <- FITs[[1]] + SU_mrCoeff <- SU_mr$coefficients + # if (length(SU_mr)>1) { + # SU_mrCoeff <- SU_mr$coefficients + # } else { SU_mrCoeff <- rep(NA,5) } + #te2$cov.unscaled[1,1]*te2$sigma^2 + #VCOV <- vcov(mr) + #V_V <- VCOV/SU_mr$sigma^2 + V_V <- SU_mr$cov.unscaled SEsPure <- sqrt(diag(V_V)*meanPureErr) - pot_est <- c(exp(s_mr['r',1]),exp(s_mr['r',1]-qt(0.975,nrow(all_l)-5)*SEsPure['r']), - exp(s_mr['r',1]+qt(0.975,nrow(all_l)-5)*SEsPure['r'])) + pot_est <- c(exp(SU_mrCoeff['r',1]),exp(SU_mrCoeff['r',1]-qt(0.975,nrow(all_l)-5)*SEsPure['r']), + exp(SU_mrCoeff['r',1]+qt(0.975,nrow(all_l)-5)*SEsPure['r'])) # unrestricted - s_mu <- summary(mu)$coefficients - SU_mu <- summary(mu) - VCOVu <- vcov(mu) - V_Vu <- VCOVu/SU_mu$sigma^2 + SU_mu <- FITs[[2]] + s_muCoeff <- SU_mu$coefficients + #SU_mu <- summary(mu) + #VCOVu <- vcov(mu) + V_Vu <- SU_mu$cov.unscaled SEsPureU <- sqrt(diag(V_Vu)*meanPureErr) - potU_est <- c(exp(s_mu['r',1]),exp(s_mu['r',1]-qt(0.975,nrow(all_l)-8)*SEsPureU['r']), - + exp(s_mu['r',1]+qt(0.975,nrow(all_l)-8)*SEsPureU['r'])) + potU_est <- c(exp(s_muCoeff['r',1]),exp(s_muCoeff['r',1]-qt(0.975,nrow(all_l)-8)*SEsPureU['r']), + + exp(s_muCoeff['r',1]+qt(0.975,nrow(all_l)-8)*SEsPureU['r'])) } # PureErrFlag - startlistmr_log <- list(a=max(all_l$readouttrans), b=SLOPE, d=min(all_l$readouttrans), cs=mean(all_l$log_dose),r=0) - - tryCatch({ - mr_log <- gsl_nls(fn = readouttrans ~ a+(d-a)/(1+exp(b*(cs-r*isSample-log_dose))), - data=all_l, - start=startlistmr_log, - control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) - }, - error = function(err) { - err$message - }) + FITstrans <- Fitting_FUNC(ro_new, TransFlag = TRUE) - startlistmu_log <- list(as=max(ro_new[,2]), bs=SLOPE, ds=min(ro_new[,2]), cs=mean(all_l$log_dose), - at=max(ro_new[,2]), bt=SLOPE, dt=min(ro_new[,2]), r=0) - - tryCatch({ - mu_log <- gsl_nls(fn = readouttrans ~ as*isRef + at*isSample + (ds*isRef + dt*isSample - as*isRef - at*isSample)/ - (1+isRef*exp(bs*(cs - log_dose)) + isSample*exp(bt*(cs-r*isSample-log_dose))), - data=all_l, - start=startlistmu_log, - control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) - }, - error = function(err) { - err$message - }) - - pot_est_log <- exp(confintd(mr_log, "r", method="asymptotic")) - potU_est_log <- exp(confintd(mu_log, "r", method="asymptotic")) + # pot_est_log <- exp(confintd(mr_log, "r", method="asymptotic")) + # potU_est_log <- exp(confintd(mu_log, "r", method="asymptotic")) + pot_est_log <- FITstrans[[3]] + potU_est_log <- FITstrans[[4]] colnames(pot_est_log) <- c("estimate","lowerCI2","upperCI") colnames(potU_est_log) <- c("estimate","lowerCI2","upperCI") #browser() - su_mr_log <- summary(mr_log) + su_mr_log <- FITstrans[[1]] #summary(mr_log) Dat$RMSE_Rlog <- su_mr_log$sigma - su_mu_log <- summary(mu_log) + su_mu_log <- FITstrans[[2]] # summary(mu_log) Dat$RMSE_Ulog <- su_mu_log$sigma Dat$up_lowAslog <- su_mu_log$coefficients[1,1] - su_mu_log$coefficients[4,1] - potALL <- rbind(pot_est, potU_est, pot_est_log, potU_est_log) + potALL <- rbind(round(pot_est,6), round(potU_est,6), round(pot_est_log,6), round(potU_est_log,6)) potALL2 <- cbind(c("restricted","unrestricted","ln-transformed restr","ln-transformed unrestr"), potALL) return(potALL2) } -ParamCI_F <- function(xt,xs,se_xt, se_xs, CoVarlog,DFs, Conf=0.975) { +#' Calculates logarithmized CIs for ratios +#' +#' Ratio CIs can be unbound, and not calulatable. Therefore the ratios are logarithmized. +#' The Covariance is not used, as it is 0 for comparisons between products. +#' +#' @param xs The estimate of the reference sample. +#' @param xt The estimate of the test sample. +#' @param se_xs The standard error of the estimate of the reference sample. +#' @param se_xt The standard error of the estimate of the test sample. +#' @param CoVar The covariance between test and reference estimate (set to 0). +#' @param DFs Degrees of freedom, either from pure error term (no of readouts - no of concentrations*2), or RMSE term (no of readouts - 8 parameters). +#' @param Conf The confidence of the Confidence Interval (default 95% which requires 0.975). +#' @returns A data-frame with the lower and upper CI in anti-log form. +#' @export +#' @examples +#' xs=2; xt=3.2; se_xt=0.34;se_xs=0.23; DFs=32-16 +#' +#' +#' ParamCI_F(xt,xs,se_xt, se_xs, CoVar,DFs, Conf=0.975) + + + +ParamCI_F <- function(xt,xs,se_xt, se_xs, CoVar,DFs, Conf=0.975) { log_xs <- log(abs(xs)) log_xt <- log(abs(xt)) - var_log_xs <- (se_xs/xs)^2 # approximate variance of log(bs) + var_log_xs <- (se_xs/xs)^2 # approximate variance of log(xs) var_log_xt <- (se_xt/xt)^2 - se_log_ratio <- sqrt(var_log_xs + var_log_xt) #-2*CoVarlog) + se_log_ratio <- sqrt(var_log_xs + var_log_xt) #-2*CoVar/(xs*xt) lower_log_ratio <- log_xt-log_xs - qt(Conf,DFs)*se_log_ratio upper_log_ratio <- log_xt-log_xs + qt(Conf,DFs)*se_log_ratio @@ -627,8 +797,31 @@ ParamCI_F <- function(xt,xs,se_xt, se_xs, CoVarlog,DFs, Conf=0.975) { return(ci_ratio) } + +#' Calculates EQ tests and F-Tests for 4P +#' +#' CIs of asmptote ratios,slope ratios, asympt-difference ratio, and lower asympt-diff. calculated. +#' In addition the F-test on significance of regression and non-linearity are calculated. +#' +#' @param ro_new Data frame with the concentrations and the respective readouts. +#' @param Lim The limits to check if the relative CI is within the limits. +#' @param PureErrFlag Boolean if Pure Error should be used. +#' @returns A data-frame with the lower and upper CI in anti-log form. +#' @export +#' @examples +#' +#' dat <- data.frame(REF1=c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2=c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497), +#' REF3=c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665), +#' TEST1=c(1405, 1523, 1502, 1474, 2383, 3221, 3589, 3445), TEST2=c(1420, 1516, 1544, 1512, 2226, 3219, 3327, 3591), +#' TEST3=c(1399, 1376, 1588, 1475, 2148, 3083, 2942, 3466), log_dose=c(5.01,3.401,2.708,2.015,1.32176,0.62861,-0.0645385,-1.6739764)) +#' Lim <- c(-1, 1, 0.005 , 2, 0.5, 2, 0.5,2,75,133, 75, 133) +#' PureErrF <- FALSE +#' +#' tests_FUNC(ro_new, Lim,PureErrF) + + tests_FUNC <- function(ro_new, Lim, PureErrFlag) { - #browser() + all_l <- melt(data.frame(ro_new), id.vars="log_dose", variable.name="replname", value.name = "readout") isRef <- rep(c(1,0),1,each=nrow(all_l)/2) isSample <- rep(c(0,1),1,each=nrow(all_l)/2) @@ -636,56 +829,20 @@ tests_FUNC <- function(ro_new, Lim, PureErrFlag) { all_l$isSample <- isSample all_l$Conc <- exp(all_l$log_dose) all_l$readout[all_l$readout < 0] <- 0.01 - tryCatch({ - pot <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")), - pmodels=data.frame(1,1,1,isSample)) - }, - error = function(msg){ - return(0) }) - compParm(pot, "c",display=T) - ED50 <- ED(pot,c(50), interval="delta") - PotEst <- ED50[1,1]/ED50[2,1] - potAll <- EDcomp(pot, percVec=c(50,50), interval="delta", display=FALSE) - potAll2 <- potAll[1:3] - - CORro <- cor(ro_new[,1], ro_new[,ncol(ro_new)]) - - if (CORro<0) SLOPE <- -1 else SLOPE <- 1 - startlist <- list(a=max(ro_new[,2]), b=SLOPE, d=min(ro_new[,2]), cs=mean(all_l$log_dose),r=0) - tryCatch({ - mr <- gsl_nls(fn = readout ~ a+(d-a)/(1+exp(b*(cs-r*isSample-log_dose))), - data=all_l, - start=startlist, - control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) - }, - error = function(msg){ - return(0) }) - - startlistmu <- list(as=max(ro_new[,2]), bs=SLOPE, ds=min(ro_new[,2]), cs=mean(all_l$log_dose), - at=max(ro_new[,2]), bt=SLOPE, dt=min(ro_new[,2]), r=0) - tryCatch({ - mu <- gsl_nls(fn = readout ~ as*isRef + at*isSample + (ds*isRef + dt*isSample - as*isRef - at*isSample)/ - (1+isRef*exp(bs*(cs - log_dose)) + isSample*exp(bt*(cs-r*isSample-log_dose))), - data=all_l, - start=startlistmu, - control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) - }, - error = function(msg){ - return(0) }) - - smu <- tryCatch({ summary(mu) }, - error=function(msg){ - return(0) }) - - POTr_CI <- potAll2[2:3] + #browser() + FITs <- Fitting_FUNC(ro_new = ro_new, TransFlag = FALSE) + POTr_CI <- FITs[[3]][2:3] + potAll2 <- FITs[[3]] + smr <- FITs[[1]] + smu <- FITs[[2]] FitAnova <- anova(lm(readout ~ factor(Conc)*isSample, all_l)) # pure error pureSS <- FitAnova[4,2] pureSS_df <- FitAnova[4,1] meanPureErr <- FitAnova[4,3] - vcovMU <- vcov(mu) - V_V <- vcovMU/smu$sigma^2 + #vcovMU <- vcov(mu) + V_V <- smu$cov.unscaled SEsPure <- sqrt(diag(V_V)*meanPureErr) VCOVpure <- V_V*meanPureErr DFsPure <- FitAnova[4,1] @@ -695,33 +852,25 @@ tests_FUNC <- function(ro_new, Lim, PureErrFlag) { testPOTr <- logical() if (POTr_CI[1]*100>Lim[[9]] & POTr_CI[2]*100Lim[[9]]/100 & potAllU2[3] Lim[[1]]) test_la_diff <- 0 else test_la_diff <- 1 #### EQ test on upper asymptote ratio ---- - as <- coeffs["as"] - at <- coeffs["at"] - uAsRatio <- compParm(potU, "a","/",display=F) - uAsCI <- c(uAsRatio[1]-qt(0.975,RSS_df)*uAsRatio[2], uAsRatio[1]+qt(0.975,RSS_df)*uAsRatio[2]) + # as <- coeffs["as"] + # at <- coeffs["at"] + # uAsRatio <- compParm(potU, "a","/",display=F) + # uAsCI <- c(uAsRatio[1]-qt(0.975,RSS_df)*uAsRatio[2], uAsRatio[1]+qt(0.975,RSS_df)*uAsRatio[2]) #browser() ds <- smu$coefficients["ds",1] dt <- smu$coefficients["dt",1] if (PureErrFlag) se_ds <- sqrt(VCOVpure["ds","ds"]) else se_ds <- smu$coefficients["ds",2] if (PureErrFlag) se_dt <- sqrt(VCOVpure["dt","dt"]) else se_dt <- smu$coefficients["dt",2] if (PureErrFlag) CoVarlog_d <- VCOVpure["dt","ds"] else CoVarlog_d <- vcovMU["dt","ds"] - if (PureErrFlag) DFs <- DFsPure else DFs <- nrow(all_l)-noConc - uAsCI2 <- ParamCI_F(dt,ds,se_dt, se_ds,CoVarlog_d, DFs, Conf=0.9975) + if (PureErrFlag) DFs <- DFsPure else DFs <- nrow(all_l)-8 + uAsCI2 <- ParamCI_F(dt,ds,se_dt, se_ds,CoVarlog_d, DFs, Conf=0.975) if (uAsCI2[1] > Lim[[7]] & uAsCI2[2] < Lim[[8]]) test_a <- 0 else test_a <- 1 estUppA <- round(at/as,5) Dat$uAsCI <- uAsCI2 #### EQ test on slope ratio ---- - bs <- coeffs["bs"] - bt <- coeffs["bt"] - slopeRatio <- compParm(potU, "b","/",display=F) - slopeCI <- c(slopeRatio[1,1]-qt(0.975,RSS_df)*slopeRatio[1,2], slopeRatio[1,1]+qt(0.975,RSS_df)*slopeRatio[1,2]) + # bs <- coeffs["bs"] + # bt <- coeffs["bt"] + # slopeRatio <- compParm(potU, "b","/",display=F) + # slopeCI <- c(slopeRatio[1,1]-qt(0.975,RSS_df)*slopeRatio[1,2], slopeRatio[1,1]+qt(0.975,RSS_df)*slopeRatio[1,2]) bs <- smu$coefficients["bs",1] bt <- smu$coefficients["bt",1] @@ -792,12 +943,12 @@ tests_FUNC <- function(ro_new, Lim, PureErrFlag) { if (slopeCI2[1] > Lim[[5]] & slopeCI2[2] < Lim[[6]]) test_b <- 0 else test_b <- 1 estUppA <- round(at/as,5) - Dat$slopeRatioCI <- slopeCI + Dat$slopeRatioCI <- slopeCI2 #### EQ test on lower As ratio ---- - lAsRatio <- compParm(potU, "d","/",display=F) - slopeCI <- c(lAsRatio[1,1]-qt(0.975,RSS_df)*lAsRatio[1,2], lAsRatio[1,1]+qt(0.975,RSS_df)*lAsRatio[1,2]) + # lAsRatio <- compParm(potU, "d","/",display=F) + # slopeCI <- c(lAsRatio[1,1]-qt(0.975,RSS_df)*lAsRatio[1,2], lAsRatio[1,1]+qt(0.975,RSS_df)*lAsRatio[1,2]) as <- smu$coefficients["as",1] at <- smu$coefficients["at",1] @@ -811,10 +962,10 @@ tests_FUNC <- function(ro_new, Lim, PureErrFlag) { Dat$lAsCI <- lAsCI2 #### EQtest on ratio of As difference ---- - AsDiffRatio <- (at-dt)/(as-ds) + AsDiffRatio <- (dt-at)/(ds-as) - at_dt <- (at-dt) - as_ds <- (as-ds) + dt_at <- (dt-at) + ds_as <- (ds-as) se_ds_asPure <- sqrt(VCOVpure["as","as"]+VCOVpure["ds","ds"]-2*VCOVpure["as","ds"]) se_dt_atPure <- sqrt(VCOVpure["at","at"]+VCOVpure["dt","dt"]-2*VCOVpure["at","dt"]) se_ds_asRMSE <- sqrt(vcovMU["as","as"]+vcovMU["ds","ds"]-2*vcovMU["as","ds"]) @@ -822,7 +973,7 @@ tests_FUNC <- function(ro_new, Lim, PureErrFlag) { if (PureErrFlag) se_ds_as <- se_ds_asPure else se_ds_as <- se_ds_asRMSE if (PureErrFlag) se_dt_at <- se_dt_atPure else se_dt_at <- se_dt_atRMSE - AsDiffCI2 <- ParamCI_F(at_dt,as_ds,se_dt_at, se_ds_as,CoVarlog=0, DFs, Conf=0.975) + AsDiffCI2 <- ParamCI_F(dt_at,ds_as,se_dt_at, se_ds_as,CoVar=0, DFs, Conf=0.975) if (AsDiffCI2[1] > Lim[[11]] & AsDiffCI2[2] < Lim[[12]]) test_ad <- 0 else test_ad <- 1 estLowA <- round(at/as,5) @@ -845,7 +996,7 @@ tests_FUNC <- function(ro_new, Lim, PureErrFlag) { testPOTr, test_c), estimate = c(round(p_F_regr, 3), round(lAs_diff, 5), estLowA, round(bs/bt,5), estUppA, p_F_nonlin, - round(at_dt/as_ds, 5), round(potAll2[1]*100,2),round(potAllU2[1]*100,2)), + round(dt_at/ds_as, 5), round(potAll2[1]*100,2),round(potAllU2[1]*100,2)), lower_limit = c("-",Lim[[1]],Lim[[3]],Lim[[5]],Lim[[7]],"-",Lim[[11]],Lim[[9]],Lim[[9]]), upper_limit = c("-",Lim[[2]],Lim[[4]],Lim[[6]],Lim[[8]],"-",Lim[[12]],Lim[[10]],Lim[[10]]), lower_CI = c(RMSE_r, round(lCI_laDiff,3), round(lAsCI2[1],5), round(slopeCI2[1],5), @@ -856,7 +1007,26 @@ tests_FUNC <- function(ro_new, Lim, PureErrFlag) { return(res_tab) } -ANOVA4plUnresfunc <- function(ro_new, sigmoid) { +#' ANOVA unrestricted model +#' +#' ANOVA with unrestricted model. +#' +#' +#' @param ro_new Data frame with the concentrations and the respective readouts. +#' @returns A data-frame with the analysis of variance. +#' @export +#' @examples +#' +#' ro_new <- data.frame(REF1=c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2=c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497), +#' REF3=c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665), +#' TEST1=c(1405, 1523, 1502, 1474, 2383, 3221, 3589, 3445), TEST2=c(1420, 1516, 1544, 1512, 2226, 3219, 3327, 3591), +#' TEST3=c(1399, 1376, 1588, 1475, 2148, 3083, 2942, 3466), log_dose=c(5.01,3.401,2.708,2.015,1.32176,0.62861,-0.0645385,-1.6739764)) +#' +#' +#' ANOVA4plUnresfunc(ro_new) + + +ANOVA4plUnresfunc <- function(ro_new) { all_l <- melt(data.frame(ro_new), id.vars="log_dose", variable.name="replname", value.name = "readout") all_len <- nrow(all_l) isRef <- rep(c(1,0),1,each=all_len/2) @@ -866,26 +1036,31 @@ ANOVA4plUnresfunc <- function(ro_new, sigmoid) { all_l$Conc <- exp(all_l$log_dose) all_l$readout[all_l$readout < 0] <- 0.01 - pot <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")), - pmodels=data.frame(1,1,1,isSample)) - potU <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")), - pmodels=data.frame(isSample, isSample,isSample,isSample)) - SStreat <- round(sum((potU$predres[,1] - mean(all_l$readout))^2),5) + FITs <- Fitting_FUNC(ro_new = ro_new, TransFlag = FALSE) + smr <- FITs[[1]] + smu <- FITs[[2]] + smrPREDs <- FITs[[5]] + smuPREDs <- FITs[[6]] + # pot <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")), + # pmodels=data.frame(1,1,1,isSample)) + # potU <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")), + # pmodels=data.frame(isSample, isSample,isSample,isSample)) + SStreat <- round(sum((smuPREDs - mean(all_l$readout))^2),5) SStreat_df <- length(unique(all_l$log_dose))-1 - SSregr <- round(sum((predict(pot)-mean(all_l$readout))^2),5) + SSregr <- round(sum((smrPREDs-mean(all_l$readout))^2),5) ## Non-parallel - SSnonparallel <- round(sum(resid(pot)^2) - sum(resid(potU)^2),5) + SSnonparallel <- round(sum(smr$residuals^2) - sum(smu$residuals^2),5) ## Preparation SSprep <- round(sum((predict(lm(readout ~ isSample, all_l)) - mean(all_l$readout))^2),5) ## Resid Err - RSS <- round(sum(potU$predres[,2]^2),5) + RSS <- round(sum(smu$residuals^2),5) RSS_df <- nrow(all_l)-SStreat_df-1 FitAnova <- anova(lm(readout ~ factor(Conc)*isSample, all_l)) # PureErr - SSE <- FitAnova[4,3] + SSE <- round(FitAnova[4,3],5) SSE_df <- FitAnova[4,1] # Non-Linearity - SSnonlin <- round(sum((predict(lm(readout ~ factor(Conc)*isSample, all_l)) - predict(potU))^2),4) + SSnonlin <- round(sum((predict(lm(readout ~ factor(Conc)*isSample, all_l)) - smuPREDs)^2),4) LoF_df <- FitAnova[1,1]+FitAnova[2,1] ## Total SStot <- round(sum((all_l$readout -mean(all_l$readout))^2),5) @@ -915,6 +1090,26 @@ ANOVA4plUnresfunc <- function(ro_new, sigmoid) { return(ANOVAtab) } +#' Calculates descriptive statistics per concentration +#' +#' Mean, SD, and CV% per concentration and product.. +#' +#' +#' @param ro_new Data frame with the concentrations and the respective readouts. +#' @returns A data-frame with the concentrations and metadata. +#' @export +#' @examples +#' +#' ro_new <- data.frame(REF1=c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2=c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497), +#' REF3=c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665), +#' TEST1=c(1405, 1523, 1502, 1474, 2383, 3221, 3589, 3445), TEST2=c(1420, 1516, 1544, 1512, 2226, 3219, 3327, 3591), +#' TEST3=c(1399, 1376, 1588, 1475, 2148, 3083, 2942, 3466), log_dose=c(5.01,3.401,2.708,2.015,1.32176,0.62861,-0.0645385,-1.6739764)) +#' +#' +#' perConcTab(ro_new, noDilSeries=3) + + + perConcTab <- function(ro_new, noDilSeries) { Reftab <- ro_new[,c(1:noDilSeries)] Testtab <- ro_new[,c((noDilSeries+1):(2*noDilSeries))] @@ -938,6 +1133,25 @@ perConcTab <- function(ro_new, noDilSeries) { } +#' Calculates dilution series. +#' +#' Recursive function to calcuate a geometric dilution series. +#' +#' @param x Starting concentration. +#' @param Div Divisor to get next concentration. +#' @param N Number of dilution step, increases. +#' @param res Vector of results. +#' @param noDil Final number of concentrations -1 (the initial one). +#' @returns A data-frame with the concentrations and metadata. +#' @export +#' @examples +#' +#' x <- 1; Div <- 3;N <- 0; res <- c(); noDil <- 7 +#' +#' divFUN(x,Div,N,res,noDil) + + + divFUN <- function(x,Div,N,res,noDil) { N <- N+1 y <- x/Div @@ -963,9 +1177,9 @@ ui <- dashboardPage( # menuSubItem(icon = NULL, tags$li(a("Document", target="self",href="UserManual.pdf"))) # ), menuItem("EXCEL upload", tabName="Dataupload", icon=icon("magnet", lib="glyphicon")), - menuItem("4PL + report", tabName="fourPL", icon=icon("chart-line", lib="font-awesome")), + menuItem("4PL metadata + report", tabName="fourPL", icon=icon("chart-line", lib="font-awesome")), #menuItem("XLSX diagnostics", tabName="XLdiagn", icon=icon("chart-bar", lib="font-awesome")), - menuItem("Linear regression + report", tabName="pla", icon=icon("pencil", lib="glyphicon")), + # menuItem("Linear regression + report", tabName="pla", icon=icon("pencil", lib="glyphicon")), menuItem("Wizard", tabName="wizard", icon=icon("chart-column", lib="font-awesome")), menuItem("Documentation", tabName="documentation", icon=icon("chart-area", lib="font-awesome")) ), @@ -1053,14 +1267,26 @@ server <- function(input, output, session) { div(checkboxInput("PureErr", "Should pure error be used for calculation of CIs?", FALSE), style = "font-size: 24px !important;color: #C2173F"), - checkboxGroupInput("selectedSSTs", "Which suitability tests to be used?", choices= c("F-test on Regr."="1", - "EQ-test on lower asymptote difference"= "2", - "EQ-test on ratio of lower asymptote"= "3","EQ-test on ratio of Hill slopes"= "4", - "EQ-test on ratio of upper asymptote"= "5", "F-test on non-linearity"="6", - "EQ-test on ratio of asymptote differences"= "7"), - selected= c("1","2","3","4","5","6","7")), + #actionLink("selectall","SelectAll"), h5("\n\n\n Author: Franz Innerbichler, InnerAnalytics")), + column(6, + h4("Suitability tests for 4-parametric logistic regression"), + checkboxGroupInput("selectedSSTs", "Which suitability tests to be used?", choices= c("F-test on Regr."="1", + "EQ-test on lower asymptote difference"= "2", + "EQ-test on ratio of lower asymptote"= "3","EQ-test on ratio of Hill slopes"= "4", + "EQ-test on ratio of upper asymptote"= "5", "F-test on non-linearity"="6", + "EQ-test on ratio of asymptote differences"= "7"), + selected= c("1","2","3","4","5","6","7")), + h4("Suitability tests for Parallel Line Assay"), + checkboxGroupInput("selectedSSTsLinear", "Which suitability tests to be used?", + choices= c("F-test on Regr."="1", + "F-test on non-linearity"= "2", + "F-test on R^2 A"= "3","F-test on R^2 B"= "4", + "F-test on slope A"= "5", "F-test on slope B"="6", + "F-test on non-parallelism"= "7", "F-test on preparation"="8"), + selected= c("1","2","3","4","5","6","7","8")) + ) ), tabPanel("4pl-Analysis", @@ -1088,9 +1314,11 @@ server <- function(input, output, session) { plotOutput("XLplot"), plotOutput("diagnplot"), - - tableOutput("ANOVAXLS"), - DTOutput("EQtests")) + DTOutput("EQtests"), + DTOutput("pottab4pl"), + DTOutput("pottab4plTrans"), + tableOutput("ANOVAXLS") + ) ))), tabPanel("linear Analysis", @@ -1101,13 +1329,7 @@ server <- function(input, output, session) { column(12, numericInput("Limits",p("limit to be >", bsButton("q4",label="", icon=icon("info"), style="primary", size="extra-small")), bsPopover(id="q4", title="", content="The calculated limits ...")), - checkboxGroupInput("selectedSSTsLinear", "Which suitability tests to be used?", - choices= c("F-test on Regr."="1", - "F-test on non-linearity"= "2", - "F-test on R^2 A"= "3","F-test on R^2 B"= "4", - "F-test on slope A"= "5", "F-test on slope B"="6", - "F-test on non-parallelism"= "7", "F-test on preparation"="8"), - selected= c("1","2","3","4","5","6","7","8")), + ) )), mainPanel( @@ -1266,7 +1488,7 @@ server <- function(input, output, session) { htmlOutput("PureErrW3"), plotOutput("plot", width = "80%"), - DT::dataTableOutput("pottab4pl"), + DTOutput("pottab4pl"), "Footnote: test performed on relative CIs.", DTOutput("EQtests4pl"), # SSTs @@ -1467,56 +1689,35 @@ server <- function(input, output, session) { #### XLSX eval ---- if (CORro<0) SLOPE <- -1 else SLOPE <- 1 - ec50est <- (max(all_l$log_dose)+min(all_l$log_dose))/2 - startlist <- list(a=min(all_l$readout), b=SLOPE, d=max(all_l$readout), cs=ec50est,r=0) - tryCatch({ - mr <- gsl_nls(fn = readout ~ a+(d-a)/(1+exp(b*(cs-r*isSample-log_dose))), - data=all_l, - start=startlist, - control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) - }, - error = function(err) { - err$message - }) + FITs <- Fitting_FUNC(XLdat2, TransFlag = FALSE) + - startlistmu <- list(as=min(all_l$readout), bs=SLOPE, ds=max(all_l$readout), cs=ec50est, - at=min(all_l$readout), bt=SLOPE, dt=max(all_l$readout), r=0) - tryCatch({ - mu <- gsl_nls(fn = readout ~ as*isRef + at*isSample + (ds*isRef + dt*isSample - as*isRef - at*isSample)/ - (1+isRef*exp(bs*(cs - log_dose)) + isSample*exp(bt*(cs-r*isSample-log_dose))), - data=all_l, - start=startlistmu, - control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) - }, - error = function(err) { - err$message - }) - - Smr <- summary(mr) - Smu <- summary(mu) + Smr <- FITs[[1]] # summary(mr) + Smu <- FITs[[2]] # summary(mu) coeffsMR <- Smr$coefficients[,1] coeffsMU <- Smu$coefficients[,1] Dat$coeffsMRes <- coeffsMR Dat$coeffsMUnr <- coeffsMU + Dat$coeffs_UN <- coeffsMU names(coeffsMU) <- c("lowAsym REF", "slope REF","upperAsym REF","EC50 REF","lowAsym TEST","slope TEST","upperAsym TEST","r") if (!PureErrFlag) { - pot_est <- exp(confintd(mr, "r", method="asymptotic")) - potU_est <- exp(confintd(mu, "r", method="asymptotic")) + pot_est <- FITs[[3]] + potU_est <- FITs[[4]] colnames(pot_est) <- c("estimate","lowerCI","upperCI") colnames(potU_est) <- c("estimate","lowerCI","upperCI") } else { FitAnova <- anova(lm(readout ~ factor(log_dose)*isSample, all_l)) meanPureErr <- FitAnova[4,3] DFsPure <- FitAnova[4,1] - VCOV <- vcov(mr) - V_V <- VCOV/Smr$sigma^2 + #VCOV <- vcov(mr) + V_V <- Smr$cov.unscaled #VCOV/Smr$sigma^2 VCOVpure <- V_V*meanPureErr SEsPure <- sqrt(diag(V_V)*meanPureErr) pot_est <- data.frame(estimate=exp(coeffsMR[5]), lowerCI = exp(coeffsMR[5]-qt(0.975,DFsPure)*SEsPure[5]), - upperCI = exp(coeffsMR[5]+qt(0.975,DFsPure)*SEsPure[5])) - VCOVu <- vcov(mu) - V_Vu <- VCOVu/Smu$sigma^2 + upperCI = exp(coeffsMR[5]+qt(0.975,DFsPure)*SEsPure[5])) + #VCOVu <- vcov(mu) + V_Vu <- Smu$cov.unscaled #VCOVpure <- V_Vu*meanPureErr SEsPureU <- sqrt(diag(V_Vu)*meanPureErr) potU_est <- data.frame(estimate=exp(coeffsMU[7]), lowerCI = exp(coeffsMU[7]-qt(0.975,DFsPure)*SEsPureU[7]), @@ -1524,21 +1725,12 @@ server <- function(input, output, session) { } - pot <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")), - pmodels=data.frame(1,1,1,isSample)) - potU <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")), - pmodels=data.frame(isSample, isSample,isSample,isSample)) - SR <- summary(pot) - SU <- summary(potU) - coeffs_UN <- potU$coefficients - coeffs_UN[1] <- ifelse(xor(CORro>0, coeffs_UN[1]>0), -coeffs_UN[1],coeffs_UN[1]) - coeffs_UN[2] <- ifelse(xor(CORro>0, coeffs_UN[2]>0), -coeffs_UN[2],coeffs_UN[2]) - coeffs_UN[7:8] <- log(coeffs_UN[7:8]) - POTU <- EDcomp(potU, percVec = c(50,50), interval="delta",display=F) - Dat$potDiffXL <- POTU[1]*100 - RMSE_unr_diagn <- sqrt(SU$resVar) - RMSE_res_diagn <- sqrt(SR$resVar) - up_lowDiffDiagn <- SU$coefficients[5,1] - SU$coefficients[3,1] + + #browser() + Dat$potDiffXL <- potU_est[1]*100 + RMSE_unr_diagn <- Smu$sigma # sqrt(SU$resVar) + RMSE_res_diagn <- Smr$sigma #sqrt(SR$resVar) + up_lowDiffDiagn <- Smu$coefficients["ds",1]-Smu$coefficients["as",1] ProzSD_diagn <- RMSE_unr_diagn*100/up_lowDiffDiagn Dat$ProzSD_XL <- ProzSD_diagn @@ -1573,89 +1765,22 @@ server <- function(input, output, session) { output$relpotTestTab <- renderTable({ pot_est3 }) }) - SStreat <- round(sum((potU$predres[,1] - mean(all_l$readout))^2),5) - SStreat_df <- length(unique(all_l$log_dose))-1 - SSregr <- round(sum((predict(pot)-mean(all_l$readout))^2),5) - ## Non-parallel - SSnonparallel <- round(sum(resid(pot)^2) - sum(resid(potU)^2),5) - ## Preparation - SSprep <- round(sum((predict(lm(readout ~ isSample, all_l)) - mean(all_l$readout))^2),5) - ## Resid Err - RSS <- round(sum(potU$predres[,2]^2),5) - RSS_df <- nrow(all_l)-SStreat_df-1 - FitAnova <- anova(lm(readout ~ factor(Conc)*isSample, all_l)) - # PureErr - SSE <- FitAnova[4,3] - SSE_df <- FitAnova[4,1] - # Non-Linearity - SSnonlin <- round(sum((predict(lm(readout ~ factor(Conc)*isSample, all_l)) - predict(potU))^2),4) - LoF_df <- FitAnova[1,1]+FitAnova[2,1] - ## Total - SStot <- round(sum((all_l$readout -mean(all_l$readout))^2),5) - MSE <- RSS/RSS_df - noConc <- length(unique(all_l$Conc)) - AnovaDFs <- c(noConc-1, 1,3,noConc-4-1, nrow(all_l)-noConc, noConc, nrow(all_l)-noConc-noConc, nrow(all_l)-1) - p_SStreat <- round(pf((SStreat/AnovaDFs[1])/MSE, AnovaDFs[1],RSS_df, lower.tail = F),3) - p_SSprep <- round(pf((SSprep/AnovaDFs[2])/MSE, AnovaDFs[2],RSS_df, lower.tail = F),3) - p_SSregr <- round(pf((SSregr/AnovaDFs[3])/MSE, AnovaDFs[3],RSS_df, lower.tail = F),3) - p_SSnonp <- round(pf((SSnonparallel/AnovaDFs[4])/MSE, AnovaDFs[3],RSS_df, lower.tail = F),3) - p_SSLoF <- round(pf((SSnonlin/LoF_df)/(SSE/SSE_df), LoF_df,SSE_df, lower.tail = F),5) - - ANOVAtab2 <- data.frame(Source = c("Treatment","Preparation","Regression", - "Non-Parallelism","Residual Error","Non-linearity", - "Pure Error","Total"), - DF = round(AnovaDFs,0), - SumSquares = c(SStreat, SSprep,SSregr, SSnonparallel, - RSS, SSnonlin,SSE, SStot), - MeanSquares = c(round(SStreat/AnovaDFs[1],3), SSprep, round(SStreat/AnovaDFs[3],3),round(SSnonparallel/AnovaDFs[4],3), - round(MSE,5), round(SSnonlin/LoF_df,5), round(SSE/SSE_df,5),""), - "F-value" = c(round((SStreat/AnovaDFs[1])/MSE,5), round((SSprep/AnovaDFs[2])/MSE,5), - round((SSregr/AnovaDFs[3])/MSE,5),round((SSnonparallel/AnovaDFs[4])/MSE,5), - "",round((SSnonlin/LoF_df)/(SSE/SSE_df),5),"",""), - "p_value" = c(round(p_SStreat,3), p_SSprep, round(p_SSregr,3), p_SSnonp,"",p_SSLoF,"","") - ) - + + ANOVAtab2 <- ANOVA4plUnresfunc(ro_new = XLdat2) output$ANOVAXLS <- renderTable({ ANOVAtab2 }) REP$ANOVAXLS <- ANOVAtab2 - #browser() + #browser() + FITsTrans <- Fitting_FUNC(XLdat2, TransFlag = TRUE) + # + + SUlog <- FITsTrans[[2]] + SRlog <- FITsTrans[[1]] + RMSE_unrlog_diagn <- SUlog$sigma + RMSE_reslog_diagn <- SRlog$sigma - startlistlog <- list(a=min(log(all_l$readout)), b=SLOPE, d=max(log(all_l$readout)), cs=ec50est,r=0) - tryCatch({ - mrlog <- gsl_nls(fn = log(readout) ~ a+(d-a)/(1+exp(b*(cs-r*isSample-log_dose))), - data=all_l, - start=startlistlog, - control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) - }, - error = function(err) { - print("Error in mrlog gsl_nls") - }) - - startlistmulog <- list(as=min(log(all_l$readout)), bs=SLOPE, ds=max(log(all_l$readout)), cs=ec50est, - at=min(log(all_l$readout)), bt=SLOPE, dt=max(log(all_l$readout)), r=0) - tryCatch({ - mulog <- gsl_nls(fn = log(readout) ~ as*isRef + at*isSample + (ds*isRef + dt*isSample - as*isRef - at*isSample)/ - (1+isRef*exp(bs*(cs - log_dose)) + isSample*exp(bt*(cs-r*isSample-log_dose))), - data=all_l, - start=startlistmulog, - control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) - }, - error = function(err) { - print("Error in murlog gsl_nls") - }) - logpot <- drm(log(readout) ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")), - pmodels=data.frame(1,1,1,isSample)) - logpotU <- drm(log(readout) ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")), - pmodels=data.frame(isSample, isSample,isSample,isSample)) - Smrlog <- summary(mrlog)$coefficients - Smulog <- summary(mulog)$coefficients - SUlog <- summary(logpotU) - SRlog <- summary(logpot) - RMSE_unrlog_diagn <- sqrt(SUlog$resVar) - RMSE_reslog_diagn <- sqrt(SRlog$resVar) - - up_lowDifflogDiagn <- SUlog$coefficients[5, 1] - SUlog$coefficients[3, 1] + up_lowDifflogDiagn <- SUlog$coefficients["ds",1]-SUlog$coefficients["as",1] ProzSDlog_diagn <- RMSE_unrlog_diagn * 100 / up_lowDifflogDiagn #### Diagnostic RMSE table #### @@ -1669,8 +1794,8 @@ server <- function(input, output, session) { Dat$DiagnTable <- DiagnTable REP$DiagnTable <- DiagnTable - logpotest <- exp(confintd(mrlog, "r", method = "asymptotic")) # compParm(logpot, "c") - logpotuest <- exp(confintd(mulog, "r", method = "asymptotic")) # compParm(logpotu, "c") + logpotest <- FITsTrans[[3]] #exp(confintd(mrlog, "r", method = "asymptotic")) # compParm(logpot, "c") + logpotuest <- FITsTrans[[4]] # exp(confintd(mulog, "r", method = "asymptotic")) # compParm(logpotu, "c") # Berechnung der Konfidenzintervalle (CI) # logpotCI <- c(exp(Smrlog[5,1] - qt(0.975, nrow(all_1)-5) * Smrlog[5,2]), exp(Smrlog[5,1]), exp(Smrlog[5,1] + qt(0.975, nrow(all_1)-5) * Smrlog[5,2])) @@ -1716,7 +1841,7 @@ server <- function(input, output, session) { names(coeffs_R) <- c("lower A", "slope", "upper A", "EC50 REF", "EC50 TEST") # coeffs_R[4] <- log(coeffs_R[4]) # coeffs_R[5] <- log(coeffs_R[5]) - # --- Ergebnistabelle: RESTRICTED (Eingeschränktes Modell) --- + # --- Ergebnistabelle: RESTRICTED --- PLAAusw <- data.frame( Information = c("model", "lower asymptote", "Hill's slope", "upper asymptote","EC50 Ref", "EC50 Test", "relative potency", @@ -1732,7 +1857,7 @@ server <- function(input, output, session) { REP$PLBend <- PLAAusw2 # --- Koeffizienten-Extraktion --- - logcoeffs_R <- Smrlog[, 1] # logpot$coefficients + logcoeffs_R <- SRlog$coefficients[, 1] # logpot$coefficients names(logcoeffs_R) <- c("lower A", "Hill's slope", "upper A", "EC50 REF","EC50 DIFF") # --- Ergebnistabelle: LOG RESTRICTED --- @@ -1747,7 +1872,7 @@ server <- function(input, output, session) { output$logcoeffs_r <- renderTable({ LogPLAAusw }) REP$LogPLAausw <- LogPLAAusw - logcoeffs_UNR <- Smulog[,1] + logcoeffs_UNR <- SUlog$coefficients[,1] names(logcoeffs_UNR) <- c("lower asymptote Ref", "Hill's slope Ref", "upper asymptote Ref","EC50 Ref", "lower asymptote Test", "Hill's slope Test", "upper asymptote Test","EC50 Diff" ) @@ -1767,7 +1892,7 @@ server <- function(input, output, session) { }) REP$LogUnrPLAausw <- LogUnrPLAAusw #browser() - Dat$coeffs_UN <- coeffs_UN + if (exists("Ind")) { @@ -1775,7 +1900,7 @@ server <- function(input, output, session) { } else Dat$dilution <- exp(XLdat[,logI]) # --- Plot-Ausgabe --- output$XLplot <- renderPlot({ - plot_f(XLdat2, sigmoid = NULL, det_sig = coeffs_UN, TransFlag=F) + plot_f(XLdat2, sigmoid = NULL, det_sig = coeffsMU, TransFlag=F) }) REP$XLdat2 <- XLdat2 @@ -1783,21 +1908,27 @@ server <- function(input, output, session) { # --- Diagnose-Plots (Residualanalyse) --- output$diagnplot <- renderPlot({ op <- par(mfrow = c(2, 2), mar = c(3.2, 3.2, 2, .5), mgp = c(2, .7, 0)) - + PREDs <- FITs[[5]] + PREDsU <- FITs[[6]] # 1. Residuals vs Fitted - plot(residuals(pot) ~ fitted(pot), main = "Residuals restricted") + plot(Smr$residuals ~ PREDs, main = "Residuals restricted") abline(h = 0) - qqnorm(residuals(pot)) - qqline(residuals(pot)) - plot(residuals(potU) ~ fitted(potU), main = "Residuals unrestricted") + qqnorm(Smr$residuals) + qqline(Smr$residuals) + plot(Smu$residuals ~ PREDsU, main = "Residuals unrestricted") abline(h = 0) - qqnorm(residuals(potU)) - qqline(residuals(potU)) + qqnorm(Smu$residuals) + qqline(Smu$residuals) par(op) # Parameter zurücksetzen }) + + pot <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")), + pmodels=data.frame(1,1,1,isSample)) + potU <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")), + pmodels=data.frame(isSample, isSample,isSample,isSample)) output$AIC <- renderTable({ AIC <- AIC(pot, potU) }) @@ -2389,9 +2520,12 @@ server <- function(input, output, session) { slopeDiffCI <- t(data.frame(LIN[[3]])) colnames(slopeDiffCI) <- c("slope difference","lower CI","upper CI") output$SlopeDiffCI <- renderTable({ slopeDiffCI },digits=5) - #browser() + + SelTestsL <- as.numeric(input$selectedSSTsLinear) + df2 <- df[SelTestsL,] + Dat$ANOVA <- df[,4:length(df)] - dat <- datatable(df[,1:3], + dat <- datatable(df2[,1:3], options=list( paging=T, dom="t",rownames=F )) %>% formatStyle("test_results", target="row",backgroundColor = styleEqual(c(-1,0,1), @@ -2402,7 +2536,7 @@ server <- function(input, output, session) { #### output 4PL ANOVA tests --- output$ANOVA <- DT::renderDataTable({ sigmoid <- sigmoid() - tab <- ANOVA4plUnresfunc(sim2(),sigmoid) + tab <- ANOVA4plUnresfunc(sim2()) # ,sigmoid dat <- datatable(tab, options=list( dom="t",rownames=F @@ -2469,11 +2603,11 @@ server <- function(input, output, session) { pottab4_ <- data.frame(pottab4) - pottab4_$potency <- as.numeric(pottab4[,2])*100 - pottab4_$`lower95%CI` <- as.numeric(pottab4[,3])*100 - pottab4_$`upper95%CI` <- as.numeric(pottab4[,4])*100 - pottab4_$relative_lowerCL <- round(pottab4_[,6]/pottab4_[,5]*100,3) - pottab4_$relative_upperCL <- round(pottab4_[,7]/pottab4_[,5]*100,3) + pottab4_$potency <- round(as.numeric(pottab4[,2])*100,1) + pottab4_$`lower95%CI` <- round(as.numeric(pottab4[,3])*100,2) + pottab4_$`upper95%CI` <- round(as.numeric(pottab4[,4])*100,2) + pottab4_$relative_lowerCL <- round(pottab4_[,6]/pottab4_[,5]*100,2) + pottab4_$relative_upperCL <- round(pottab4_[,7]/pottab4_[,5]*100,2) if (as.numeric(pottab4_$relative_lowerCL[1]) > Lim[[9]] & as.numeric(pottab4_$relative_upperCL[1]) < Lim[[10]] ) { test_potCI <- 0 @@ -2492,17 +2626,17 @@ server <- function(input, output, session) { output$pottab4pl <- DT::renderDataTable({ dat <- datatable(pottab4_[1:2,], - options=list( + options=list(digits=3, paging=T, dom="t",rownames=F - )) %>% formatStyle("test_result", target="row",backgroundColor = styleEqual(c(0,1), - c("lightgreen","pink"))) + )) %>% formatStyle("test_result", target="row",backgroundColor = styleEqual(c(0,1), + c("lightgreen","pink"))) }) output$pottab4plTrans <- DT::renderDataTable({ dat <- datatable(pottab4_[3:4,], - options=list( + options=list(digits=3, paging=T, dom="t",rownames=F )) %>% formatStyle("test_result", target="row",backgroundColor = styleEqual(c(0,1), - c("lightgreen","pink"))) + c("lightgreen","pink"))) }) }) @@ -2546,10 +2680,10 @@ server <- function(input, output, session) { if (is.null(Dat$coeffs_UN)) { SAMPLE50 <- sigmoid[1] + (sigmoid[3] - sigmoid[1])/(1+exp(sigmoid[5]*( (sigmoid[7]+0.693147)- dils2))) SAMPLE200 <- sigmoid[1] + (sigmoid[3] - sigmoid[1])/(1+exp(sigmoid[5]*( (sigmoid[7]-0.693147)-dils2))) - Xbend50l <- sigmoid[7] + 0.693147-1.31696/sigmoid[5] - Xbend200l <- sigmoid[7] - 0.693147-1.31696/sigmoid[5] - Xbend50u <- sigmoid[7] + 0.693147+1.31696/sigmoid[5] - Xbend200u <- sigmoid[7] - 0.693147+1.31696/sigmoid[5] + Xbend50l <- sigmoid[7] + 0.693147-1.5434/sigmoid[5] + Xbend200l <- sigmoid[7] - 0.693147-1.5434/sigmoid[5] + Xbend50u <- sigmoid[7] + 0.693147+1.5434/sigmoid[5] + Xbend200u <- sigmoid[7] - 0.693147+1.5434/sigmoid[5] Xbend50 <- max(Xbend50l, Xbend50u) Xbend200 <- min(Xbend200l, Xbend200u) dummy <- plot_f(tab,sigmoid,det_sig=NULL) @@ -2558,10 +2692,10 @@ server <- function(input, output, session) { #browser() SAMPLE50 <- det_sig[3] + (det_sig[5] - det_sig[3])/(1+exp(det_sig[1]*(det_sig[7]+0.693147-dils2))) SAMPLE200 <- det_sig[3] + (det_sig[5] - det_sig[3])/(1+exp(det_sig[1]*(det_sig[7]-0.693147-dils2))) - Xbend50l <- det_sig[7] + 0.693147-1.31696/det_sig[1] - Xbend200l <- det_sig[7] - 0.693147-1.31696/det_sig[1] - Xbend50u <- det_sig[7] + 0.693147+1.31696/det_sig[1] - Xbend200u <- det_sig[7] - 0.693147+1.31696/det_sig[1] + Xbend50l <- det_sig[7] + 0.693147-1.5434/det_sig[1] + Xbend200l <- det_sig[7] - 0.693147-1.5434/det_sig[1] + Xbend50u <- det_sig[7] + 0.693147+1.5434/det_sig[1] + Xbend200u <- det_sig[7] - 0.693147+1.5434/det_sig[1] Xbend50 <- max(Xbend50l, Xbend50u) Xbend200 <- min(Xbend200l, Xbend200u) dummy <- plot_f(tab,sigmoid=NULL,det_sig=det_sig)