2307 lines
98 KiB
R
2307 lines
98 KiB
R
|
|
|
|
library(shiny)
|
|
library(shinyjs)
|
|
library(shinyAce)
|
|
library(purrr)
|
|
library(gslnls)
|
|
library(tidyverse)
|
|
library(ggplot2)
|
|
library(reshape2)
|
|
library(openxlsx)
|
|
library(DT)
|
|
library(ggpubr)
|
|
library(gridExtra)
|
|
library(drc)
|
|
library(twopartm)
|
|
library(car)
|
|
library(dplyr)
|
|
library(tinytex)
|
|
|
|
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)
|
|
}
|
|
|
|
plot_f <- function(dat, 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")
|
|
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(is.null(det_sig)) {
|
|
if (CORdat<0) {
|
|
startlist <- list(a=sigmoid[1], b=-sigmoid[5],cs=sigmoid[7],
|
|
d=sigmoid[3],r=sigmoid[8])
|
|
} else {
|
|
|
|
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*(log_dose-(cs-r*isSample)))),
|
|
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
|
|
})
|
|
|
|
a <- s_mr$coefficients[1,1]
|
|
b <- s_mr$coefficients[2,1]
|
|
cs <- s_mr$coefficients[3,1]
|
|
d <- s_mr$coefficients[4,1]
|
|
r <- s_mr$coefficients[5,1]
|
|
|
|
log_dose <- unique(all_l$log_dose)
|
|
seq_x <- seq(min(log_dose),max(log_dose),0.1)
|
|
SAMPLE <- a+(d-a)/(1+exp(b*(seq_x-(cs-r))))
|
|
REF <- a+(d-a)/(1+exp(b*(seq_x-(cs))))
|
|
|
|
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)))
|
|
}
|
|
#browser()
|
|
pl_df <- cbind(seq_x, SAMPLE, REF, SAMPLEtrue, REFtrue)
|
|
all_l2$readout[all_l2$readout < 0] <- 0.01
|
|
all_l2$readouttrans <- log(all_l2$readout)
|
|
slopeEC50 <- b*(a-d)/4
|
|
|
|
Xbendl3 <- cs-(1.31696/b)
|
|
Xbendu3 <- cs+(1.31696/b)
|
|
XbendlT <- cs-r-(1.31696/b)
|
|
XbenduT <- cs-r+(1.31696/b)
|
|
bendpoints <- c(bendREF_lower = round(Xbendl3,3), bendREF_upper=round(Xbendu3,3),
|
|
bendSAMPLE_lower = round(XbendlT,3), bendSAMPLE_upper=round(XbenduT,3))
|
|
Dat$bendpoints <- bendpoints
|
|
Dat$cfordils <- cs
|
|
|
|
p <- ggplot(all_l2, aes(x=log_dose, y=readout, color=factor(isRef))) +
|
|
geom_point(shape=factor(isRef), alpha=0.8) +
|
|
labs(title = paste("restricted 4pl; bendp:", round(Xbendl3,3),round(Xbendu3,3),round(XbendlT,3),round(XbenduT,3)),
|
|
color="product") +
|
|
scale_color_manual(labels=c("test","reference"), values=c("red","blue")) +
|
|
scale_shape_manual(labels=c("test","reference")) +
|
|
theme_bw() +
|
|
theme(axis.text = element_text(size=14))
|
|
|
|
p2 <- p + geom_line(data=as.data.frame(pl_df), aes(x=seq_x, y=SAMPLE), color="red",
|
|
inherit.aes = F) +
|
|
geom_line(data=as.data.frame(pl_df), aes(x=seq_x, y=REF), color="blue",
|
|
inherit.aes = F) +
|
|
geom_line(data=as.data.frame(pl_df), aes(x=seq_x, y=SAMPLEtrue), color="red", linetype=2, alpha=0.4,
|
|
inherit.aes = F) +
|
|
geom_line(data=as.data.frame(pl_df), aes(x=seq_x, y=REFtrue), color="blue", linetype=2, alpha=0.4,
|
|
inherit.aes = F) +
|
|
geom_vline(xintercept=c(Xbendl3, Xbendu3), col="blue",linetype=2) +
|
|
geom_vline(xintercept=c(XbendlT, XbenduT), col="red",linetype=2) +
|
|
annotate("text", x=cs, y=a+(d-a)/2, label="0", size=5) +
|
|
theme(legend.position="none")
|
|
Dat$p2 <- p2
|
|
|
|
# transformed plots
|
|
p_rt <- ggplot(all_l2, aes(x=log_dose, y=readouttrans, color=factor(isRef))) +
|
|
geom_point(shape=factor(isRef), alpha=0.8) +
|
|
labs(title = paste("restricted transformed 4pl"), color="product") +
|
|
scale_color_manual(labels=c("test","reference"), values=c("red","blue")) +
|
|
theme_bw()
|
|
|
|
mrt <- gsl_nls(fn = readouttrans ~ a+(d-a)/(1+exp(b*(log_dose-(cs-r*isSample)))),
|
|
data=all_l2,
|
|
start=startlist,
|
|
control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6))
|
|
s_mrt <- summary(mrt)
|
|
a_trans <- s_mrt$coefficients[1,1]
|
|
b_trans <- s_mrt$coefficients[2,1]
|
|
cs_trans <- s_mrt$coefficients[3,1]
|
|
d_trans <- s_mrt$coefficients[4,1]
|
|
r_trans <- s_mrt$coefficients[5,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)
|
|
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
|
|
SAMPLEtrans <- a_trans+(d_trans-a_trans)/(1+exp(b_trans*(seq_x-(cs_trans-r_trans))))
|
|
REFtrans <- a_trans+(d_trans-a_trans)/(1+exp(b_trans*(seq_x-(cs_trans))))
|
|
|
|
pl_df_trans <- cbind(seq_x, SAMPLEtrans, REFtrans)
|
|
p_rt2 <- p_rt + geom_line(data=as.data.frame(pl_df_trans), aes(x=seq_x, y=SAMPLEtrans), color="red",
|
|
inherit.aes = F) +
|
|
geom_line(data=as.data.frame(pl_df_trans), aes(x=seq_x, y=REFtrans), color="blue",
|
|
inherit.aes = F) +
|
|
geom_vline(xintercept=c(XbendlTrans, XbenduTrans), col="blue",linetype=2) +
|
|
geom_vline(xintercept=c(XbendlTransT, XbenduTransT), col="red",linetype=2) +
|
|
theme(legend.position = "none", axis.text=element_text(size=14))
|
|
|
|
if (is.null(det_sig)) {
|
|
unrestr <- drm(readout ~ exp(log_dose), isSample, data=all_l2, fct=LL.4(),
|
|
pmodels=data.frame(isSample, isSample,isSample,isSample))
|
|
Sum_u <- summary(unrestr)
|
|
ast <- Sum_u$coefficients[3,1]
|
|
ate <- Sum_u$coefficients[4,1]
|
|
bst <- Sum_u$coefficients[1,1]
|
|
bte <- Sum_u$coefficients[2,1]
|
|
cst <- log(Sum_u$coefficients[7,1])
|
|
cte <- log(Sum_u$coefficients[8,1])
|
|
dst <- Sum_u$coefficients[5,1]
|
|
dte <- Sum_u$coefficients[6,1]
|
|
} else {
|
|
ast <- det_sig[5]
|
|
ate <- det_sig[6]
|
|
bst <- det_sig[1]
|
|
bte <- det_sig[2]
|
|
cst <- det_sig[7]
|
|
cte <- det_sig[8]
|
|
dst <- det_sig[3]
|
|
dte <- det_sig[4]
|
|
}
|
|
REFu <- ast + (dst-ast)/(1+exp(bst*(seq_x-cst)))
|
|
SAMPLEu <- ate + (dte-ate)/(1+exp(bte*(seq_x-cte)))
|
|
pl_df2 <- cbind(seq_x, SAMPLEu, REFu)
|
|
#browser()
|
|
pu <- ggplot(all_l2, aes(x=log_dose, y=readout, color=factor(isRef))) +
|
|
geom_point() +
|
|
labs(title="unrestricted 4_pl-Model", color="product") +
|
|
scale_color_manual(labels = c("test","reference"), values=c("red","blue")) +
|
|
theme_bw()
|
|
pu2 <- pu + geom_line(data=as.data.frame(pl_df2), aes(x=seq_x, y=SAMPLEu),
|
|
color="red", inherit.aes = F) +
|
|
geom_line(data=as.data.frame(pl_df2), aes(x=seq_x, y=REFu),
|
|
color="blue", inherit.aes = 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() +
|
|
labs(title="unrestricted transformed 4_pl-Model", color="product") +
|
|
scale_color_manual(labels = c("test","reference"), values=c("red","blue")) +
|
|
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]
|
|
|
|
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)))
|
|
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),
|
|
color="red", inherit.aes = F) +
|
|
geom_line(data=as.data.frame(pl_df2u_t), aes(x=seq_x, y=REFu_trans),
|
|
color="blue", inherit.aes = F,
|
|
show.legend = F)
|
|
pu3_t <- pu2_t
|
|
grid.arrange(p2,p_rt2,pu2_,pu3_t, nrow=2)
|
|
}
|
|
|
|
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) {
|
|
yAxfac <- (ds-as)
|
|
log_dose <- log_conc
|
|
isRef <- rep(c(1,0),1,each=length(log_conc)*noDilSeries)
|
|
isSample <- rep(c(0,1),1,each=length(log_conc)*noDilSeries)
|
|
#browser()
|
|
av <- as*isRef + at*isSample + (ds*isRef + dt*isSample - as*isRef - at*isSample)/
|
|
(1+isRef*exp(bs*(cs - log_dose)) + isSample*exp(bt*(ct-log_dose)))
|
|
if (heteroNoise) {
|
|
# heterosc noise
|
|
ro_jit <- matrix(unlist(map(av, function(x) x+rnorm(1,0,x*sd_fac/100))), nrow=noDils, ncol=noDilSeries*2)
|
|
} else {
|
|
# homosc noise
|
|
ro_jit <- matrix(unlist(map(av, function(x) x+rnorm(1,0,sd_fac*yAxfac/100))), nrow=noDils, ncol=noDilSeries*2)
|
|
}
|
|
ro_jit <- abs(ro_jit)
|
|
|
|
ro_jit2 <- cbind(ro_jit, log_dose)
|
|
if (noDilSeries==3) {
|
|
colnames(ro_jit2) <- c("R_dil1","R_dil2","R_dil3","T_dil1","T_dil2","T_dil3", "log_dose")
|
|
} else {
|
|
colnames(ro_jit2) <- c("R_dil1","R_dil2","T_dil1","T_dil2", "log_dose")
|
|
}
|
|
return(ro_jit2)
|
|
}
|
|
|
|
LinPotTab <- function(circles, Lim, PureErrFlag) {
|
|
circ_ABl <- circles
|
|
circ_Al <- circ_ABl[circ_ABl$isSample ==1,]
|
|
circ_Al <- circ_ABl[circ_ABl$isSample ==0,]
|
|
# restr CSSI model
|
|
modAB <- lm(readout ~ log_dose + isSample, circ_ABl)
|
|
coeffs <- modAB$coefficients
|
|
SU_modAB <- tryCatch({
|
|
SU_modAB <- summary(modAB)
|
|
}, error = function(msg) {
|
|
return(NA)
|
|
})
|
|
# Intercept diff/slope modAB
|
|
linPot <- exp(modAB$coefficients[3]/modAB$coefficients[2])
|
|
|
|
if(PureErrFlag) {
|
|
FitAnova <- anova(lm(readout ~ factor(log_dose)*isSample, circ_ABl))
|
|
meanPureErr <- FitAnova[4,3]
|
|
DFsPure <- FitAnova[4,1]
|
|
VCOV <- vcov(modAB)
|
|
V_V <- VCOV/SU_modAB$sigma^2
|
|
VCOVpure <- V_V*meanPureErr
|
|
SEsPure <- sqrt(diag(V_V)*meanPureErr)
|
|
}
|
|
|
|
log_pot_delta <- deltaMethod(modAB, "isSample/log_dose")
|
|
if (PureErrFlag) {
|
|
V_ <- log_pot_delta$SE^2/SU_modAB$sigma^2
|
|
V_p <- V_*meanPureErr
|
|
potDeltaPureSE <- sqrt(V_p)
|
|
CI_log_low <- log_pot_delta$Estimate - qt(0.975, DFsPure)*potDeltaPureSE
|
|
CI_log_up <- log_pot_delta$Estimate + qt(0.975, DFsPure)*potDeltaPureSE
|
|
} else {
|
|
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
|
|
}
|
|
#browser()
|
|
ExpLinPot <- exp(c(log_pot_delta$Estimate, CI_log_low, CI_log_up))
|
|
if (ExpLinPot[2]*100>Lim[[9]] & ExpLinPot[3]*100<Lim[[10]]) test_potCI <- 0 else test_potCI <- 1
|
|
# Rel pot CI
|
|
relLinpotCI <- ExpLinPot/linPot*100
|
|
pottab <- cbind(round(linPot*100,3), round(ExpLinPot[2]*100,3), round(ExpLinPot[3]*100,3),
|
|
round(test_potCI,3), round(relLinpotCI[2],3),round(relLinpotCI[3],3))
|
|
colnames(pottab) <- c("Potency","lower 95%CI", "upper 95%CI", "test_result", "lowerRel95%CI","upperRel95%CI")
|
|
return(pottab)
|
|
}
|
|
|
|
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)
|
|
isSample <- rep(c(0,1),1,each=nrow(all_l)/2)
|
|
all_l$isRef <- isRef
|
|
all_l$isSample <- isSample
|
|
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
|
|
|
|
circ_ABl <- circles
|
|
circ_Al <- circ_ABl[circ_ABl$isSample ==1,]
|
|
circ_Bl <- circ_ABl[circ_ABl$isSample ==0,]
|
|
# restr CSSI model
|
|
modAB <- lm(readout ~ log_dose + isSample, circ_ABl)
|
|
# unrestr with interact term SSSI
|
|
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
|
|
ExpLinPot <- exp(c(log_pot_delta$Estimate, CI_log_low, CI_log_up))
|
|
if (ExpLinPot[2]*100>Lim[9] & ExpLinPot[3]*100>Lim[10]) test_potCI <- 0 else test_potCI <- 1
|
|
|
|
su_mod <- summary(modAB)$coefficients
|
|
su_mod2 <- cbind(data.frame(parameter = c("intercept REF","slope REF","intercepts diff.")), su_mod)
|
|
su_modU <- summary(modABu)$coefficients
|
|
su_modU2 <- cbind(data.frame(parameter = c("intercept REF","slope REF","intercepts diff.","slope difference")), su_modU)
|
|
|
|
uCI_SloDiff <- su_modU[4,1] + qt(0.975,8)*su_modU[4,2]
|
|
lCI_SloDiff <- su_modU[4,1] - qt(0.975,8)*su_modU[4,2]
|
|
SlopeDiffCI <- c(su_modU[4,1], lCI_SloDiff,uCI_SloDiff)
|
|
|
|
lenCirc <- nrow(circ_ABl)
|
|
dfTreat <- 3
|
|
dfPrep <- 1
|
|
dfReg <- 1
|
|
dfnonP <- 1
|
|
dfRMSE <- c(lenCirc-3-1)
|
|
dfTotal <- lenCirc-1
|
|
dfPureE <- lenCirc-6
|
|
dfNonLin <- dfRMSE-dfPureE
|
|
|
|
RSS <- sum(resid(lm(readout ~ log_dose*isSample, circ_ABl))^2)
|
|
MSE <- RSS/dfRMSE
|
|
SSE <- sum(resid(lm(readout ~ factor(log_dose)*isSample, circ_ABl))^2)
|
|
MSpure <- SSE/dfPureE
|
|
|
|
if (PureErrFlag) {
|
|
FitAnova <- anova(lm(readout ~ factor(log_dose)*isSample, circ_ABl))
|
|
meanPureErr <- FitAnova[4,3]
|
|
SU_modAB <- tryCatch({
|
|
SU_modAB <- summary(modAB)
|
|
}, error = function(msg) {
|
|
return(NA)
|
|
})
|
|
if (length(SU_modAB)>1) s_modABcoeffs <- summary(modAB)$coefficients
|
|
|
|
DFsPure <- FitAnova[4,1]
|
|
VCOV <- vcov(modAB)
|
|
V_V <- VCOV/SU_modAB$sigma^2
|
|
VCOVpure <- V_V*meanPureErr
|
|
SEsPure <- sqrt(diag(V_V)*meanPureErr)
|
|
su_mod2[,3] <- SEsPure
|
|
su_mod2[,4] <- su_mod2[,2]/su_mod2[,3]
|
|
su_mod2[,5] <- 2*(1-pt(abs(su_mod[,4]), FitAnova[4,1]))
|
|
|
|
s_mu <- summary(modABu)$coefficients
|
|
SU_modABu <- summary(modABu)
|
|
VCOVu <- vcov(modABu)
|
|
V_Vu <- VCOVu/SU_modABu$sigma^2
|
|
SEsPureU <- sqrt(diag(V_Vu)*meanPureErr)
|
|
|
|
su_modU2[,3] <- SEsPureU
|
|
su_modU2[,4] <- su_modU2[,2]/su_modU2[,3]
|
|
su_modU2[,5] <- 2*(1-pt(abs(su_modU2[,4]), FitAnova[4,1]))
|
|
|
|
uCI_SloDiffP <- su_modU[4,1] + qt(0.975,8)*SEsPureU[4]
|
|
lCI_SloDiffP <- su_modU[4,1] - qt(0.975,8)*SEsPureU[4]
|
|
SlopeDiffCI <- c(su_modU[4,1], lCI_SloDiffP,uCI_SloDiffP)
|
|
|
|
SSRes <- SSE
|
|
dfRes <- dfPureE
|
|
|
|
} else {
|
|
SSRes <- RSS
|
|
dfRes <- dfRMSE
|
|
}
|
|
|
|
# treatment
|
|
SStreat <- print(sum((predict(lm(readout ~ factor(log_dose)*isSample, circ_ABl))-mean(circ_ABl$readout))^2))
|
|
F_treat <- (SStreat/dfTreat)/(SSRes/dfRes)
|
|
# Preparation
|
|
SSprep <- print(sum((predict(lm(readout ~ isSample, circ_ABl))-mean(circ_ABl$readout))^2))
|
|
F_prep <- (SSprep/dfTreat)/(SSRes/dfRes)
|
|
# Regression
|
|
# ANOVA tape II SS of regression
|
|
SSreg <- Anova(lm(readout ~log_dose + isSample, circ_ABl))[1,1]
|
|
# Non-parallelism
|
|
# diff of RSS of restricted and unrestricted model
|
|
SSnonpar <- sum(resid(modAB)^2) - sum(resid(modABu)^2)
|
|
F_nonpar <- SSnonpar/(sum(resid(lm(readout ~ factor(log_dose)*isSample, circ_ABl))^2)/(lenCirc-4))
|
|
|
|
# non-linearity
|
|
SSnonlin <- sum((predict(modABu)-predict(lm(readout ~ as.factor(log_dose)*isSample, circ_ABl)))^2)
|
|
# = RSS-SSE
|
|
# Total SS
|
|
SStot <- sum((circ_ABl$readout-mean(circ_ABl$readout))^2)
|
|
# Significance of R^2 F-ratio
|
|
# MSR/MSE
|
|
# sample A
|
|
F_R2_A <- sum((predict(lm(readout ~ log_dose+ I(log_dose^2), circ_Al)) - mean(predict(modA)))^2 - (predict(modA) - mean(circ_Al$readout))^2)/
|
|
(sum((predict(lm(readout ~ log_dose+ I(log_dose^2), circ_Al)) - circ_Al$readout)^2)/(nrow(circ_Al)-3))
|
|
pFR2_A <- round(pf(F_R2_A,1,6),4)
|
|
# sample B
|
|
F_R2_B <- sum((predict(lm(readout ~ log_dose+ I(log_dose^2), circ_Bl)) - mean(predict(modB)))^2 - (predict(modB) - mean(circ_Bl$readout))^2)/
|
|
(sum((predict(lm(readout ~ log_dose+ I(log_dose^2), circ_Bl)) - circ_Bl$readout)^2)/(nrow(circ_Bl)-3))
|
|
pFR2_B <- round(pf(F_R2_B,1,6),4)
|
|
# sign of non-lin with pure error: MSSnonlin/MSSE
|
|
F_nonlin <- (SSnonlin/2)/(SSE/dfPureE)
|
|
|
|
# sign of slope
|
|
F_slope_B <- sum((predict(modB) - mean(circ_Bl$readout))^2)/(sum((circ_Bl$readout - predict(modB))^2)/(nrow(circ_Bl)-2))
|
|
F_slope_A <- sum((predict(modA) - mean(circ_Al$readout))^2)/(sum((circ_Al$readout - predict(modA))^2)/(nrow(circ_Al)-2))
|
|
# F-test on regression: MSSreg/MSSE
|
|
if (is.na(F_nonlin)) F_nonlin <- 0
|
|
if (F_nonlin > 0) {
|
|
p_F_nonlin <- round(pf(F_nonlin,2,dfPureE, lower.tail = F),5)
|
|
} else { p_F_nonlin <- "SSnonlin neg or 0"; }
|
|
|
|
# significances
|
|
F_regr <- (SSreg/1)/(SSRes/dfRes)
|
|
p_F_regr <- round(pf(F_regr,1,dfRes, lower.tail = F),3)
|
|
p_F_treat <- round(pf(F_treat,3,dfRes, lower.tail = F),3)
|
|
p_F_prep <- round(pf(F_prep,1,dfRes, lower.tail = F),3)
|
|
p_F_slope_A <- round(pf(F_slope_A,1,(nrow(circ_Al)-2), lower.tail = F),3)
|
|
p_F_slope_B <- round(pf(F_slope_B,1,(nrow(circ_Bl)-2), lower.tail = F),3)
|
|
p_F_nonp <- round(pf(F_nonpar,1,dfRes, lower.tail = F),3)
|
|
p_F_LoF <- p_F_nonlin
|
|
|
|
res_tab_lin <- data.frame(test = c("F-test on sign. of regression", "F_test on non-lin",
|
|
"F-test on R^2 A","F_test on R^2 B",
|
|
"F-test on slope A","F-test on slope B",
|
|
"F-test on non-parallelism","F-test on preparation"),
|
|
test_results = c(ifelse(p_F_regr<0.05,0,1),ifelse(p_F_nonlin<0.05,0,1),
|
|
ifelse(pFR2_A<0.05,0,1),ifelse(pFR2_B<0.05,0,1),
|
|
ifelse(p_F_slope_A<0.05,0,1),ifelse(p_F_slope_B<0.05,0,1),
|
|
ifelse(p_F_nonp<0.05,0,1),ifelse(p_F_prep<0.05,0,1)),
|
|
estimate = c(p_F_regr, p_F_nonlin,pFR2_A,pFR2_B,p_F_slope_A,
|
|
p_F_slope_B,p_F_nonp,p_F_prep),
|
|
Source = c("Treatment","Preparation","Regression","Non-parallelism",
|
|
"Resid Error","Non-linearity","Pure error", "Total"),
|
|
df = c(dfTreat,1,1,1,dfRMSE,2,dfPureE,lenCirc-1),
|
|
SumSquares = c(round(SStreat,5),round(SSprep,5),round(SSreg,5),
|
|
round(SSnonpar,5),round(RSS,5),round(SSnonlin,5),
|
|
round(SSE,5),round(SStot,5)),
|
|
MS = c(round(SStreat/dfTreat,5),round(SSprep,5),round(SSreg,5),
|
|
round(SSnonpar,5),round(RSS/dfRMSE,5),round(SSnonlin/2,5),
|
|
round(SSE/dfPureE,5),round(SStot/dfTotal,5)),
|
|
"F-value" = c(round(F_treat,5), round(F_prep,5),round(F_regr,5),
|
|
round(F_nonpar,5),"",round(F_nonlin,5),"",""),
|
|
"p-value" = c(p_F_treat, p_F_prep, p_F_regr, p_F_nonp, "", p_F_LoF, "",""))
|
|
RET <- list(res_tab_lin, su_modU2, SlopeDiffCI, su_mod2)
|
|
return(RET)
|
|
}
|
|
|
|
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)
|
|
isSample <- rep(c(0,1),1,each=nrow(all_l)/2)
|
|
all_l$isRef <- isRef
|
|
all_l$isSample <- isSample
|
|
all_l$Conc <- exp(all_l$log_dose)
|
|
all_l$readout[all_l$readout < 0] <- 0.01
|
|
all_l$readouttrans <- log(all_l$readout)
|
|
|
|
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
|
|
})
|
|
|
|
if (!PureErrFlag) {
|
|
pot_est <- exp(confintd(mr, "r", method="asymptotic"))
|
|
potU_est <- exp(confintd(mu, "r", method="asymptotic"))
|
|
} 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
|
|
SEsPure <- sqrt(diag(V_V)*meanPureErr)
|
|
pot_est <- c(exp(s_mr[5,1]),exp(s_mr[5,1]-qt(0.975,nrow(all_l)-5)*SEsPure[5]),
|
|
exp(s_mr[5,1]+qt(0.975,nrow(all_l)-5)*SEsPure[5]))
|
|
# unrestricted
|
|
s_mu <- summary(mu)$coefficients
|
|
SU_mu <- summary(mu)
|
|
VCOVu <- vcov(mu)
|
|
V_Vu <- VCOVu/SU_mu$sigma^2
|
|
SEsPureU <- sqrt(diag(V_Vu)*meanPureErr)
|
|
potU_est <- c(exp(s_mu[7,1]),exp(s_mu[7,1]-qt(0.975,nrow(all_l)-8)*SEsPureU[7]),
|
|
exp(s_mu[7,1]+qt(0.975,nrow(all_l)-8)*SEsPureU[7]))
|
|
} # 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
|
|
})
|
|
|
|
|
|
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"))
|
|
colnames(pot_est_log) <- c("estimate","lowerCI2","upperCI")
|
|
colnames(potU_est_log) <- c("estimate","lowerCI2","upperCI")
|
|
|
|
su_mr_log <- summary(mr_log)
|
|
Dat$RMSE_Rlog <- su_mr_log$sigma
|
|
su_mu_log <- 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)
|
|
|
|
potALL2 <- cbind(c("restricted","unrestricted","transformed restr","untransf restr"), potALL)
|
|
return(potALL2)
|
|
}
|
|
|
|
ParamCI_F <- function(xt,xs,se_xt, se_xs, CoVarlog,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_xt <- (se_xt/xt)^2
|
|
se_log_ratio <- sqrt(var_log_xs + var_log_xt) #-2*CoVarlog)
|
|
|
|
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
|
|
ci_ratio <- exp(c(lower_log_ratio, upper_log_ratio))
|
|
return(ci_ratio)
|
|
}
|
|
|
|
tests_FUNC <- function(ro_new, 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)
|
|
isSample <- rep(c(0,1),1,each=nrow(all_l)/2)
|
|
all_l$isRef <- isRef
|
|
all_l$isSample <- isSample
|
|
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))
|
|
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(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
|
|
})
|
|
|
|
smu <- tryCatch({ summary(mu) },
|
|
error=function(err) {
|
|
return(0)
|
|
})
|
|
|
|
POTr_CI <- potAll2[2:3]
|
|
|
|
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
|
|
SEsPure <- sqrt(diag(V_V)*meanPureErr)
|
|
VCOVpure <- V_V*meanPureErr
|
|
DFsPure <- FitAnova[4,1]
|
|
|
|
|
|
|
|
testPOTr <- logical()
|
|
if (POTr_CI[1]*100>Lim[[9]] & POTr_CI[2]*100<Lim[[10]] ) testPOTr <- 0 else testPOTr <- 1
|
|
|
|
potU <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")),
|
|
pmodels=data.frame(isSample, isSample,isSample,isSample))
|
|
potAllU <- EDcomp(potU, percVec=c(50,50), interval="delta", display=FALSE)
|
|
potAllU2 <- potAllU[1:3]
|
|
sum_potU <- summary(potU)
|
|
coeffs <- potU$coefficients
|
|
coeffs[1] <- ifelse(CORro<0, -coeffs[1], coeffs[1])
|
|
coeffs[2] <- ifelse(CORro<0, -coeffs[2], coeffs[2])
|
|
names(coeffs) <- c("bs","bt","ds","dt","as","at","cs","ct")
|
|
|
|
e_c_ref <- coeffs[7]
|
|
e_c_test <- coeffs[8]
|
|
coeffs[7:8] <- log(coeffs[7:8])
|
|
test_c <- logical()
|
|
if((potAllU2[2] >Lim[[9]]/100 & potAllU2[3] <Lim[[10]]/100)) test_c <- 0 else test_c <- 1
|
|
|
|
#### ANOVA ----
|
|
noConc <- length(unique(all_l$Conc))
|
|
nofitted <- noConc
|
|
AnovaDFs <- c(nofitted-1,1,3,nofitted-4-1,nrow(all_l)-nofitted, nofitted,nrow(all_l)-2*nofitted,nrow(all_l)-1)
|
|
SStreat <- round(sum((predict(potU)-mean(all_l$readout))^2),5)
|
|
SSregr <- round(sum((predict(pot)-mean(all_l$readout))^2),5)
|
|
# non-parallelism
|
|
SSnonparall <- round(sum(resid(pot)^2)-sum(resid(potU)^2),5)
|
|
SSprep <- round(sum((predict(lm(readout ~ isSample, all_l))-mean(all_l$readout))^2),5)
|
|
|
|
RSS <- round(sum(potU$predres[,2]^2),5)
|
|
RSS_df <- AnovaDFs[5]
|
|
MSEunr <- RSS/RSS_df
|
|
RMSEunr <- sqrt(RSS/RSS_df)
|
|
# Pure Err
|
|
FitAnova <- anova(lm(readout ~ factor(Conc)*isSample, all_l))
|
|
SSE <- sum(resid(lm(readout ~ factor(Conc)*isSample, all_l))^2) # =FitAnova[4,2]
|
|
SSE_df <- FitAnova[4,1]
|
|
PureMSE <- SSE/SSE_df
|
|
RMSE_pure <- sqrt(PureMSE)
|
|
## non-lin = LoF
|
|
if (PureErrFlag) { ERR <- PureMSE; ERR_df <- SSE_df } else { ERR <- MSEunr; ERR_df <- RSS_df }
|
|
SSnonlin <- sum((predict(lm(readout ~ factor(Conc)*isSample, all_l))-predict(potU))^2)
|
|
LoF_df <- FitAnova[1,1]+FitAnova[2,1]
|
|
F_regr <- (SSregr/AnovaDFs[3])/ERR
|
|
p_F_regr <- round(pf(F_regr, AnovaDFs[3], ERR_df, lower.tail = F),5)
|
|
if (ncol(ro_new)<4) F_nonlin <- 0 else F_nonlin <- (SSnonlin/AnovaDFs[6])/ERR
|
|
if (F_nonlin > 0) {
|
|
p_F_nonlin <- round(pf(F_nonlin, AnovaDFs[6], ERR_df, lower.tail = F),5)
|
|
} else { p_F_nonlin <- "SSnonlin neg or single dilutions" }
|
|
|
|
test_a <- test_b <- test_d <- test_ad <- logical()
|
|
|
|
RSS_r <- round(sum(pot$predres[,2]^2),5)
|
|
MSE_r <- RSS_r/(nrow(all_l)-5)
|
|
RMSE_r <- round(sqrt(MSE_r),6)
|
|
Dat$RMSE_r <- RMSE_r
|
|
Dat$RMSE_pure <- RMSE_pure
|
|
Dat$RMSE_unr <- round(RMSEunr,6)
|
|
#browser()
|
|
## EQ test on lower As diff
|
|
ds <- coeffs["ds"]
|
|
dt <- coeffs["dt"]
|
|
lAs_diff <- (dt-ds)
|
|
uCI_laDiff <- lAs_diff+qt(0.975,df.residual(mu))*sqrt(sum_potU$coefficients[3,2]^2+sum_potU$coefficients[4,2]^2)
|
|
lCI_laDiff <- lAs_diff-qt(0.975,df.residual(mu))*sqrt(sum_potU$coefficients[3,2]^2+sum_potU$coefficients[4,2]^2)
|
|
if (uCI_laDiff < Lim[[2]] & lCI_laDiff > 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])
|
|
#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 (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 <- smu$coefficients["bs",1]
|
|
bt <- smu$coefficients["bt",1]
|
|
if (PureErrFlag) se_bs <- sqrt(VCOVpure["bs","bs"]) else se_bs <- smu$coefficients["bs",2]
|
|
if (PureErrFlag) se_bt <- sqrt(VCOVpure["bt","bt"]) else se_bt <- smu$coefficients["bt",2]
|
|
if (PureErrFlag) CoVarlog_b <- VCOVpure["bt","bs"] else CoVarlog_b <- vcovMU["bt","bs"]
|
|
slopeCI2 <- ParamCI_F(bt,bs,se_bt, se_bs,CoVarlog_b, DFs, Conf=0.975)
|
|
if (slopeCI2[1] > Lim[[5]] & slopeCI2[2] < Lim[[6]]) test_b <- 0 else test_b <- 1
|
|
estUppA <- round(at/as,5)
|
|
|
|
Dat$slopeRatioCI <- slopeCI
|
|
|
|
#### 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])
|
|
|
|
as <- smu$coefficients["as",1]
|
|
at <- smu$coefficients["at",1]
|
|
if (PureErrFlag) se_as <- sqrt(VCOVpure["as","as"]) else se_as <- smu$coefficients["as",2]
|
|
if (PureErrFlag) se_at <- sqrt(VCOVpure["at","at"]) else se_at <- smu$coefficients["at",2]
|
|
if (PureErrFlag) CoVarlog_a <- VCOVpure["at","as"] else CoVarlog_a <- vcovMU["at","as"]
|
|
lAsCI2 <- ParamCI_F(at,as,se_at, se_as,CoVarlog_a, DFs, Conf=0.975)
|
|
if (lAsCI2[1] > Lim[[3]] & lAsCI2[2] < Lim[[4]]) test_d <- 0 else test_d <- 1
|
|
estLowA <- round(at/as,5)
|
|
|
|
Dat$lAsCI <- lAsCI2
|
|
|
|
#### EQtest on ratio of As difference ----
|
|
AsDiffRatio <- (at-dt)/(as-ds)
|
|
|
|
at_dt <- (at-dt)
|
|
as_ds <- (as-ds)
|
|
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"])
|
|
se_dt_atRMSE <- sqrt(vcovMU["at","at"]+vcovMU["dt","dt"]-2*vcovMU["at","dt"])
|
|
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)
|
|
if (AsDiffCI2[1] > Lim[[11]] & AsDiffCI2[2] < Lim[[12]]) test_ad <- 0 else test_ad <- 1
|
|
estLowA <- round(at/as,5)
|
|
|
|
Dat$up_lowAs <- (as-ds)
|
|
|
|
lowerCIlowerA <- lAsCI2[1]; lowerCIupperA <- uAsCI2[1]; upperCIlowerA <- lAsCI2[2]; upperCIupperA <- uAsCI2[2]
|
|
test_lowA <- test_d; test_uppA <- test_a
|
|
#browser()
|
|
res_tab <- data.frame(test= c("F-test on sign. of regression*",
|
|
"EQ test on significance of lower asymptotes",
|
|
"EQ test ratio of lower asymptotes",
|
|
"EQ test ratio of Hill slopes",
|
|
"EQ test ratio of upper asymptotes",
|
|
"F-test on non-linearity*",
|
|
"EQ test ratio of asymptote difference",
|
|
"geom. rel. CI restr. model",
|
|
"geom. rel. CI unrestr. model"),
|
|
test_results = c(ifelse(p_F_regr<0.05,0,1), test_la_diff, test_lowA, test_b, test_uppA,
|
|
ifelse(p_F_nonlin>1,1, ifelse(p_F_nonlin<0.05,1,0)), test_ad,
|
|
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)),
|
|
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),
|
|
round(uAsCI2[1],5), "-", round(AsDiffCI2[1],5), round(potAll2[2],2), round(potAllU2[2],2)),
|
|
upper_CI = c(RMSE_pure, round(uCI_laDiff,3), round(lAsCI2[2],5), round(slopeCI2[2],5),
|
|
round(uAsCI2[2],5), "-", round(AsDiffCI2[2],5), round(potAll2[3],2), round(potAllU2[3],2))
|
|
)
|
|
return(res_tab)
|
|
}
|
|
|
|
ANOVA4plUnresfunc <- function(ro_new, sigmoid) {
|
|
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)
|
|
isSample <- rep(c(0,1),1,each=all_len/2)
|
|
all_l$isRef <- isRef
|
|
all_l$isSample <- isSample
|
|
all_l$Conc <- exp(all_l$log_dose)
|
|
all_l$readout[all_l$readout < 0] <- 0.01
|
|
|
|
pot <- drm(readout ~ Conc, isSample, data=all_l2, 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)
|
|
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)
|
|
|
|
ANOVAtab <- data.frame(Source = c("Treatment","Preparation","Regression",
|
|
"Non-Parallelism","Residual Error","Non-linearity",
|
|
"Pure Error","Total"),
|
|
DF = AnovaDFs,
|
|
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,"","")
|
|
)
|
|
|
|
return(ANOVAtab)
|
|
}
|
|
|
|
perConcTab <- function(ro_new, noDilSeries) {
|
|
Reftab <- ro_new[,c(1:noDilSeries)]
|
|
Testtab <- ro_new[,c((noDilSeries+1):(2*noDilSeries))]
|
|
tReftab <- t(Reftab)
|
|
colnames(tReftab) <- round(ro_new[,ncol(ro_new)],5)
|
|
|
|
avs <- apply(tReftab,2,mean)
|
|
sds <- apply(tReftab,2,sd)
|
|
cv <- sds/avs*100
|
|
tReftab2 <- rbind(tReftab, avs,sds,cv)
|
|
|
|
tTesttab <- t(Testtab)
|
|
colnames(tTesttab) <- round(ro_new[,ncol(ro_new)],5)
|
|
|
|
avs_test <- apply(tTesttab,2,mean)
|
|
sds_test <- apply(tTesttab,2,sd)
|
|
cv_test <- sds_test/avs_test*100
|
|
tTesttab2 <- rbind(tTesttab, avs_test,sds_test,cv_test)
|
|
concTab <- rbind(tReftab2, tTesttab2)
|
|
return(concTab)
|
|
|
|
}
|
|
|
|
divFUN <- function(x,Div,N,res,noDil) {
|
|
N <- N+1
|
|
y <- x/Div
|
|
res <- c(res,y)
|
|
if (N==noDil) { return(res) }
|
|
divFUN(y,Div,N,res,noDil)
|
|
}
|
|
|
|
#_______________________________________________________________________________
|
|
#### Main function ----
|
|
#_______________________________________________________________________________
|
|
|
|
function(input,output, session) {
|
|
v <- reactiveValues(num_dose=0, next.dose.t=0)
|
|
|
|
sigmoid <- reactive({
|
|
sig <- c(input$lowAsymptREF, input$lowAsymptTEST,input$uppAsymptREF,input$uppAsymptREF,
|
|
input$slopeREF,input$slopeTEST,input$EC50,input$potDiff)
|
|
sig
|
|
})
|
|
|
|
CONC <- reactive({
|
|
Konz_ <- c(input$CONC1,input$CONC2,input$CONC3,input$CONC4,
|
|
input$CONC5,input$CONC6,input$CONC7,input$CONC8,
|
|
input$CONC9,input$CONC10,input$CONC11,input$CONC12)
|
|
if (any(na.omit(Konz_)==0)) Konz_[Konz_ ==0] <- 0.0000001
|
|
Konz <- na.omit(Konz_)
|
|
})
|
|
|
|
Dils <- reactive({
|
|
Dilutions <- c(input$ConcStart,input$dilutionFac,input$NoDil,input$NoDilSer)
|
|
})
|
|
|
|
#### input EXCEL file ----
|
|
|
|
observe({
|
|
if (!is.null(input$iFile)) {
|
|
inFile <- input$iFile
|
|
ext <- tools::file_ext(inFile$name)
|
|
file.rename(inFile$datapath, paste(inFile$datapath, ".xlsx",sep=""))
|
|
t.filelocation <- gsub('\\\\','/',paste(inFile$datapath, ext,sep="."))
|
|
sheets <- openxlsx::getSheetNames(t.filelocation)
|
|
dat <- lapply(sheets, openxlsx::read.xlsx, xlsxFile = t.filelocation)
|
|
names(dat) <- sheets
|
|
Dat$wb <- dat
|
|
names(Dat$wb) <- sheets
|
|
Dat$sheets <- sheets
|
|
Dat$FileName <- input$iFile[["name"]]
|
|
|
|
}
|
|
})
|
|
output$sheetName <- renderUI({
|
|
if (!is.null(Dat$wb)) {
|
|
#browser()
|
|
cnSheets <- Dat$sheets
|
|
cnSheets2 <- c("please choose", cnSheets)
|
|
selectInput(inputId = "sheet", label="Select a sheet:",choices = cnSheets2)
|
|
}
|
|
})
|
|
observeEvent(input$sign_out, {
|
|
unlink(input$iFile$datapath)
|
|
reset(id = "") # from shinyjs package
|
|
})
|
|
|
|
|
|
#### process XLSX file ----
|
|
observe({
|
|
if (!is.null(input$iFile)) {
|
|
if (!is.null(input$sheet)) {
|
|
if (input$sheet != "please choose") {
|
|
|
|
|
|
XLdat <- Dat$wb[input$sheet][[1]]
|
|
if (is.null(XLdat)) XLdat <- Dat$wb[Dat$sheets[1]][[1]]
|
|
cn <- colnames(XLdat)
|
|
logI <- grep("log", cn)
|
|
logDoseI <- grep("log_dose", cn)
|
|
if (length(logI)>0 & length(logDoseI)==0) {
|
|
XLdat$log_dose <- XLdat[,logI]
|
|
XLdat2 <- XLdat[,-logI]
|
|
CORro <- cor(XLdat$log_dose, XLdat[,3])
|
|
} else if (length(logI)==0 & length(logDoseI)==0) {
|
|
Ind <- grep("dilu|dose|Dose|Conc|conc",cn)
|
|
XLdat$log_dose <- log(XLdat[,Ind])
|
|
CORro <- cor(XLdat[,Ind], XLdat[,3])
|
|
XLdat2 <- XLdat[,-Ind]
|
|
} else if (length(logI)>0 & length(logDoseI)>0) {
|
|
XLdat2 <- XLdat
|
|
CORro <- cor(XLdat[,logI], XLdat[,3])
|
|
}
|
|
Dat$EXCEL <- XLdat2
|
|
PureErrFlag <- input$PureErr
|
|
warning_text2 <- reactive({
|
|
ifelse(PureErrFlag, 'Pure Error is selected', '')
|
|
})
|
|
output$PureErrW2 <- renderText(warning_text2())
|
|
|
|
noDilSeries <-(ncol(XLdat2)-1)/2
|
|
noDils <- nrow(XLdat2)
|
|
Dat$noDilSeriesXL <- noDilSeries
|
|
|
|
all_l <- melt(data.frame(XLdat2), id.vars="log_dose",variable.name = "replname", value.name = "readout")
|
|
isRef <- rep(c(1,0),1,each=nrow(XLdat2)*noDilSeries)
|
|
isSample <- rep(c(0,1),1,each=nrow(XLdat2)*noDilSeries)
|
|
all_l$isRef <- isRef
|
|
all_l$isSample <- isSample
|
|
all_l$Conc <- exp(all_l$log_dose)
|
|
# all_l$readout[all_l$readout < 0] <- 0.01
|
|
REP$all_l <- all_l
|
|
|
|
#### 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
|
|
})
|
|
|
|
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)
|
|
coeffsMR <- Smr$coefficients[,1]
|
|
coeffsMU <- Smu$coefficients[,1]
|
|
Dat$coeffsMRes <- coeffsMR
|
|
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"))
|
|
colnames(pot_est) <- c("estimate","lowerCI","upperCI")
|
|
colnames(potU_est) <- c("estimate","lowerCI","upperCI")
|
|
} else {
|
|
FitAnova <- anova(lm(readout ~ factor(log_dose)*isSample, circ_ABl))
|
|
meanPureErr <- FitAnova[4,3]
|
|
DFsPure <- FitAnova[4,1]
|
|
VCOV <- vcov(mr)
|
|
V_V <- 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
|
|
#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]),
|
|
upperCI = exp(coeffsMU[7]+qt(0.975,DFsPure)*SEsPureU[7]))
|
|
|
|
}
|
|
|
|
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]
|
|
ProzSD_diagn <- RMSE_unr_diagn*100/up_lowDiffDiagn
|
|
Dat$ProzSD_XL <- ProzSD_diagn
|
|
|
|
observe({
|
|
pot_est3 <- data.frame(pot_est*100)
|
|
MaxPl <- max(input$upperPot, pot_est3$upperCI)
|
|
MinPl <- max(input$lowerPot, pot_est3$lowerCI)
|
|
MaxPl_ <- MaxPl*1.1
|
|
MinPl_ <- MinPl*0.9
|
|
|
|
p_relCI <- ggplot(data=pot_est3, aes(xmin=lowerCI, xmax=upperCI, y=1)) +
|
|
geom_linerange(size=4, col="darkseagreen",alpha=0.5) +
|
|
geom_point(aes(x=estimate, y=1), col="grey15", shape=13, size=10) +
|
|
geom_vline(xintercept = c(input$lowerPot, input$upperPot), col="indianred") +
|
|
annotate("text", x=input$lowerPot-13, y=1.040, label=paste("lower EAC:", input$lowerPot), col="indianred") +
|
|
annotate("text", x=input$upperPot+13, y=1.040, label=paste("upper EAC:", input$lupperPot), col="indianred") +
|
|
annotate("text", x=pot_est3$lowerPot-10, y=1.020, label=paste("lower CL:", round(pot_est3$lowerCI,1)), col="darkgreen") +
|
|
annotate("text", x=pot_est3$upperPot+10, y=1.020, label=paste("upper CL:", round(pot_est3$upperCI,1)), col="darkgreen") +
|
|
annotate("text", x=pot_est3$estimate, y=0.98, label=paste("rel. potency:", round(pot_est3$estimate,1)), col="black") +
|
|
ylim(c(0.95, 1.05)) +
|
|
xlim(c(MinPl_,MaxPl_)) +
|
|
xlab("relative potency + confidence interval") +
|
|
theme_bw() +
|
|
theme(axis.title.y=element_blank(),
|
|
axis.text.y=element_blank(),
|
|
axis.ticks.y=element_blank())
|
|
output$relpotTestPlot <- renderPlot({
|
|
p_relCI
|
|
})
|
|
REP$relpotTestPlot <- p_relCI
|
|
})
|
|
|
|
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,"","")
|
|
)
|
|
|
|
output$ANOVAXLS <- renderTable({ ANOVAtab2 })
|
|
|
|
REP$ANOVAXLS <- ANOVAtab2
|
|
|
|
#browser()
|
|
|
|
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]
|
|
ProzSDlog_diagn <- RMSE_unrlog_diagn * 100 / up_lowDifflogDiagn
|
|
|
|
#### Diagnostic RMSE table ####
|
|
DiagnTable <- data.frame(parameter = c("RMSE unrestricted", "RMSE_restr.", "Diff_upper-lowerAsymp", "%SD (unrestricted)",
|
|
"RMSE log_unrestricted", "RMSE log_restr", "diff_up-lowAsymp_log", "%SD (log unrestricted)"),
|
|
result = c(round(RMSE_unr_diagn, 4), round(RMSE_res_diagn, 4),
|
|
round(up_lowDiffDiagn, 4), round(ProzSD_diagn, 4),
|
|
round(RMSE_unrlog_diagn, 4), round(RMSE_reslog_diagn, 4),
|
|
round(up_lowDifflogDiagn, 4), round(ProzSDlog_diagn, 4)))
|
|
|
|
DatSDiagnTable <- DiagnTable
|
|
REPSDiagnTable <- DiagnTable
|
|
|
|
logpotest <- exp(confintd(mrlog, "r", method = "asymptotic")) # compParm(logpot, "c")
|
|
logpotuest <- 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]))
|
|
colnames(logpotest) <- c("estimate", "lowerCI", "upperCI")
|
|
|
|
colnames(logpotuest) <- c("estimate", "lowerCI", "upperCI")
|
|
|
|
cnXL <- colnames(XLdat2)
|
|
Filesample <- data.frame(c("File name", "samples"), c(Dat$Filename, paste(cnXL[1], " vs ", cnXL[4], "\n")))
|
|
colnames(Filesample) <- c("", "")
|
|
outputIFilesampl <- renderTable({ Filesample })
|
|
|
|
UnRPLAausw <- data.frame(Information = c("model", "lower asymptote Ref", "Hill's slope Ref", "upper asymptote Ref","EC50 Ref",
|
|
"lower asymptote Test", "Hill's slope Test",
|
|
"upper asymptote Test","EC50 Difference",
|
|
"relative potency", "lower CI", "upper CI"),
|
|
Results = unlist(c("UNRESTRICTED", round(coeffsMU, 3), round(potU_est, 3)))) # von psl_nls
|
|
|
|
# "log relative potency", "log lower CI", "log upper CI", round(logpotest, 3), round(compParm(potu, "c", display = F), 3)
|
|
|
|
outputcoeffs_unr <- renderTable({
|
|
UnRPLAausw
|
|
})
|
|
|
|
# browser()
|
|
UnRPLAausw2 <- data.frame(Dat$bendpointsTRANS)
|
|
if (length(UnRPLAausw2) > 0) colnames(UnRPLAausw2) <- c("bendpoints log")
|
|
|
|
outputcoeffs_unr2 <- renderTable({
|
|
UnRRPLAausw2
|
|
})
|
|
|
|
REP$UnRPLAausw <- UnRPLAausw
|
|
REP$UnRPLAausw2 <- UnRPLAausw2
|
|
|
|
# browser()
|
|
coeffs_R <- coeffsMR # pot$coefficients
|
|
coeffs_R[5] <- coeffs_R[4] - coeffs_R[5]
|
|
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) ---
|
|
PLAAusw <- data.frame(
|
|
Information = c("model", "lower asymptote", "Hill's slope", "upper asymptote","EC50 Ref",
|
|
"EC50 Test", "relative potency",
|
|
"lower CI", "upper CI"),
|
|
Results = unlist(c("RESTRICTED", round(coeffs_R, 3),
|
|
round(pot_est[1, ] * 100, 3)))) # von gs1_nls
|
|
output$coeffs_r <- renderTable({ PLAAusw })
|
|
|
|
PLAAusw2 <- data.frame(Dat$bendpoints)
|
|
output$coeffs_r2 <- renderTable({ PLAAusw2 }, digits = 3, rownames = T)
|
|
|
|
REP$PLAausw <- PLAAusw
|
|
REP$PLBend <- PLAAusw2
|
|
|
|
# --- Koeffizienten-Extraktion ---
|
|
logcoeffs_R <- Smrlog[, 1] # logpot$coefficients
|
|
names(logcoeffs_R) <- c("lower A", "Hill's slope", "upper A", "EC50 REF","EC50 DIFF")
|
|
|
|
# --- Ergebnistabelle: LOG RESTRICTED ---
|
|
|
|
LogPLAAusw <- data.frame(
|
|
Information = c("model", "lower asymptote", "Hill's slope", "upper asymptote","EC50 Ref",
|
|
"EC50 difference", "log relative potency",
|
|
"log lower CI", "log upper CI"),
|
|
Results = unlist(c("LOG RESTRICTED", round(logcoeffs_R, 3),
|
|
round(logpotest * 100, 3)))) # von gs1_nls
|
|
|
|
output$logcoeffs_r <- renderTable({ LogPLAAusw })
|
|
REP$LogPLAausw <- LogPLAAusw
|
|
|
|
logcoeffs_UNR <- Smulog[,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"
|
|
)
|
|
|
|
# --- Ergebnistabelle: LOG UNRESTRICTED ---
|
|
|
|
LogUnrPLAAusw <- data.frame(
|
|
Information = c("model", "lower asymptote Ref", "Hill's slope Ref", "upper asymptote Ref","EC50 Ref",
|
|
"lower asymptote Test", "Hill's slope Test", "upper asymptote Test","EC50 Diff" ,
|
|
"relative potency", "lower CI", "upper CI"),
|
|
|
|
Results = unlist(c("LOG UNRESTRICTED", round(logcoeffs_UNR, 3),
|
|
round(logpotest * 100, 3)))) # von gs1_nls
|
|
|
|
output$logcoeffs_unr <- renderTable({
|
|
LogUnrPLAAusw
|
|
})
|
|
REP$LogUnrPLAausw <- LogUnrPLAAusw
|
|
#browser()
|
|
Dat$coeffs_UN <- coeffs_UN
|
|
|
|
|
|
if (exists("Ind")) {
|
|
Dat$dilution <- XLdat[,Ind]
|
|
} else Dat$dilution <- exp(XLdat[,logI])
|
|
# --- Plot-Ausgabe ---
|
|
output$XLplot <- renderPlot({
|
|
plot_f(XLdat2, sigmoid = NULL, det_sig = coeffs_UN)
|
|
})
|
|
|
|
REP$XLdat2 <- XLdat2
|
|
|
|
# --- Diagnose-Plots (Residualanalyse) ---
|
|
output$diagplot <- renderPlot({
|
|
op <- par(mfrow = c(1, 2), mar = c(3.2, 3.2, 2, .5), mgp = c(2, .7, 0))
|
|
|
|
# 1. Residuals vs Fitted
|
|
plot(residuals(pot) ~ fitted(pot), main = "Residuals restricted")
|
|
abline(h = 0)
|
|
qqnorm(residuals(pot))
|
|
qqline(residuals(pot))
|
|
plot(residuals(potU) ~ fitted(potU), main = "Residuals unrestricted")
|
|
abline(h = 0)
|
|
|
|
qqnorm(residuals(potU))
|
|
qqline(residuals(potU))
|
|
|
|
par(op) # Parameter zurücksetzen
|
|
})
|
|
|
|
output$AIC <- renderTable({
|
|
AIC <- AIC(pot, potU)
|
|
})
|
|
|
|
output$VarDiagn <- renderTable({
|
|
DiagnTable
|
|
}, digits=4)
|
|
|
|
output$relpotplot <- renderPlot({
|
|
relpot(potU, intervall="fieller", bty="l",
|
|
main="Quality of rel. potency over response")
|
|
})
|
|
|
|
|
|
|
|
} # !please choose
|
|
} # input$sheet
|
|
} # inout$iFile
|
|
})
|
|
|
|
#### make geomDils reactive ----
|
|
observe({
|
|
#browser()
|
|
if (!is.na(input$ConcStart) & !is.na(input$dilutionFac) &!is.na(input$NoDil) &!is.na(input$NoDilSer)) {
|
|
upR <- input$ConcStart
|
|
noDil <- input$NoDil
|
|
noDilSer <- input$NoDilSer
|
|
Div <- input$dilutionFac
|
|
res <- c()
|
|
N_ <- 1
|
|
Conc <- c(upR, divFUN(upR,Div,N=N_,res,noDil))
|
|
|
|
Dat$MetaConc <- Conc
|
|
}
|
|
})
|
|
|
|
#### updateSlider on XLSX ----
|
|
observe({
|
|
if (!is.null(Dat$potDiffXL)) {
|
|
updateSliderInput(session, "potencydiff",
|
|
value=round(as.numeric(Dat$potDiffXL[[1]]),5))
|
|
}
|
|
})
|
|
observeEvent(input$potencydiff, {
|
|
if (!is.null(Dat$potDiffXL)) {
|
|
updateSliderInput(session, "potencydiff",
|
|
value=round(as.numeric(input$potencydiff),5))
|
|
}
|
|
})
|
|
observe({
|
|
if (!is.null(Dat$ProzSD_XL)) {
|
|
updateSliderInput(session, "sdfacf",
|
|
value=round(as.numeric(Dat$ProzSD_XL[[1]]),5))
|
|
}
|
|
})
|
|
observeEvent(input$sdfac, {
|
|
if (!is.null(Dat$ProzSD_XL)) {
|
|
updateSliderInput(session, "sdfac",
|
|
value=round(as.numeric(Dat$ProzSD_XL[[1]]),5))
|
|
}
|
|
})
|
|
|
|
#### updaterNumeric Input ----
|
|
observe({
|
|
if(!is.null(Dat$coeffs_UN)) {
|
|
updateNumericInput(session, "lowAsymptREF",
|
|
value=round(as.numeric(Dat$coeffs_UN[3]),5), min=0)
|
|
updateNumericInput(session, "lowAsymptTEST",
|
|
value=round(as.numeric(Dat$coeffs_UN[4]),5), min=0)
|
|
updateNumericInput(session, "uppAsymptREF",
|
|
value=round(as.numeric(Dat$coeffs_UN[5]),5), min=0)
|
|
updateNumericInput(session, "uppAsymptTEST",
|
|
value=round(as.numeric(Dat$coeffs_UN[6]),5), min=0)
|
|
updateNumericInput(session, "slopeREF",
|
|
value=round(as.numeric(Dat$coeffs_UN[1]),5))
|
|
updateNumericInput(session, "slopeTEST",
|
|
value=round(as.numeric(Dat$coeffs_UN[2]),5))
|
|
updateNumericInput(session, "EC50",
|
|
value=round(as.numeric(Dat$coeffs_UN[7]),5))
|
|
updateNumericInput(session, "potDiff",
|
|
value=round(as.numeric(Dat$coeffs_UN[7])- as.numeric(Dat$coeffs_UN[8]),5))
|
|
}
|
|
})
|
|
|
|
observe({
|
|
if(!is.null(Dat$dilution)) {
|
|
updateNumericInput(session, "CONC1",
|
|
value=as.numeric(Dat$dilution[1]))
|
|
updateNumericInput(session, "CONC2",
|
|
value=as.numeric(Dat$dilution[2]))
|
|
updateNumericInput(session, "CONC3",
|
|
value=as.numeric(Dat$dilution[3]))
|
|
updateNumericInput(session, "CONC4",
|
|
value=as.numeric(Dat$dilution[4]))
|
|
updateNumericInput(session, "CONC5",
|
|
value=as.numeric(Dat$dilution[5]))
|
|
updateNumericInput(session, "CONC6",
|
|
value=as.numeric(Dat$dilution[6]))
|
|
updateNumericInput(session, "CONC7",
|
|
value=as.numeric(Dat$dilution[7]))
|
|
updateNumericInput(session, "CONC8",
|
|
value=as.numeric(Dat$dilution[8]))
|
|
updateNumericInput(session, "CONC9",
|
|
value=as.numeric(Dat$dilution[9]))
|
|
updateNumericInput(session, "CONC10",
|
|
value=as.numeric(Dat$dilution[10]))
|
|
updateNumericInput(session, "CONC11",
|
|
value=as.numeric(Dat$dilution[11]))
|
|
updateNumericInput(session, "CONC12",
|
|
value=as.numeric(Dat$dilution[12]))
|
|
|
|
}
|
|
})
|
|
|
|
observe({
|
|
if(!is.null(Dat$MetaConc)) {
|
|
updateNumericInput(session, "CONC1",
|
|
value=as.numeric(Dat$MetaConc[1]))
|
|
updateNumericInput(session, "CONC2",
|
|
value=as.numeric(Dat$MetaConc[2]))
|
|
updateNumericInput(session, "CONC3",
|
|
value=as.numeric(Dat$MetaConc[3]))
|
|
updateNumericInput(session, "CONC4",
|
|
value=as.numeric(Dat$MetaConc[4]))
|
|
updateNumericInput(session, "CONC5",
|
|
value=as.numeric(Dat$MetaConc[5]))
|
|
updateNumericInput(session, "CONC6",
|
|
value=as.numeric(Dat$MetaConc[6]))
|
|
updateNumericInput(session, "CONC7",
|
|
value=as.numeric(Dat$MetaConc[7]))
|
|
updateNumericInput(session, "CONC8",
|
|
value=as.numeric(Dat$MetaConc[8]))
|
|
updateNumericInput(session, "CONC9",
|
|
value=as.numeric(Dat$MetaConc[9]))
|
|
updateNumericInput(session, "CONC10",
|
|
value=as.numeric(Dat$MetaConc[10]))
|
|
updateNumericInput(session, "CONC11",
|
|
value=as.numeric(Dat$MetaConc[11]))
|
|
updateNumericInput(session, "CONC12",
|
|
value=as.numeric(Dat$MetaConc[12]))
|
|
|
|
}
|
|
})
|
|
|
|
#### render logDilsText ----
|
|
output$logdil <- renderText({
|
|
if (!is.null(Dat$MetaConc)) {
|
|
Conc <- Dat$MetaConc
|
|
} else Conc <- CONC()
|
|
logdilu <-log(Conc)
|
|
logdilu
|
|
})
|
|
|
|
|
|
#### reactive dataset sim ----
|
|
|
|
sim <- reactive({
|
|
|
|
sd_fac_ <- as.numeric(input$sdfac)
|
|
r_ <- log(as.numeric(input$potencydiff)/100)
|
|
as = sigmoid()[1]; bs = sigmoid()[5];cs = sigmoid()[7];ds = sigmoid()[3];at = sigmoid()[2];
|
|
bt = sigmoid()[6];r = sigmoid()[8]; ct = cs-r_; dt = sigmoid()[4];
|
|
if (!is.null(Dat$MetaConc)) Conc <- Dat$MetaConc else Conc <- CONC()
|
|
log_conc <- log(Conc)
|
|
av_test <- as + (ds-as)/(1+exp(bs*(cs - log_conc)))
|
|
av_ref <- at + (dt-at)/(1+exp(bt*(ct - log_conc)))
|
|
#browser()
|
|
if (!is.na(input$NoDilSer)) {
|
|
noDilSer <- input$NoDilSer
|
|
} else if (!is.null(Dat$noDilSeriesXL)) noDilSer <- Dat$noDilSeriesXL else noDilSer <- 3
|
|
if (!is.na(input$NoDil)) noDil <- input$NoDil else noDil <- length(log_conc)
|
|
isRef <- rep(c(1,0), 1,each=noDilSer*noDil)
|
|
isSample <- rep(c(0,1), 1,each=noDilSer*noDil)
|
|
|
|
if (is.null(Dat$EXCEL)) {
|
|
ro_new <- Calc_DilRes(as=as,at=at,ds=ds,dt=dt,cs=cs,ct=ct,r=r_,bt=bt,bs=bs, log_conc = log_conc,
|
|
sd_fac=sd_fac_,
|
|
# auslenkU=outlierU,
|
|
# auslenkM=outlierM,
|
|
# auslenkL=outlierL,
|
|
heteroNoise = input$heterosked, noDilSeries = noDilSer, noDils = noDil)
|
|
} else ro_new <- Dat$EXCEL
|
|
})
|
|
# })
|
|
|
|
####sim2 ----
|
|
sim2 <- reactive({
|
|
tab <- sim()
|
|
if (is.null(Dat$EXCEL)) return(tab) else return(Dat$EXCEL)
|
|
})
|
|
|
|
|
|
#### Plot 4pl ----
|
|
output$plot <- renderPlot({
|
|
#browser()
|
|
sigmoid <- sigmoid()
|
|
det_sig=NULL
|
|
plot_f(sim2(),sigmoid, det_sig)
|
|
})
|
|
|
|
#### Testergebnisse ----
|
|
#observe({
|
|
observeEvent(input$StartCalc,{
|
|
PureErrFlag <- input$PureErr
|
|
warning_text3 <- reactive({
|
|
ifelse(PureErrFlag, 'Pure error selected','')
|
|
})
|
|
|
|
|
|
output$PureErrW3 <- renderText(warning_text3())
|
|
|
|
Limite <- list(as.numeric(input$lEACdiffla), as.numeric(input$uEACdiffla),
|
|
as.numeric(input$lEACratiola), as.numeric(input$uEACratiola),
|
|
as.numeric(input$lEACratioSlope), as.numeric(input$uEACratioSlope),
|
|
as.numeric(input$lEACratioua), as.numeric(input$uEACratioua),
|
|
as.numeric(input$lowerPot), as.numeric(input$upperPot),
|
|
as.numeric(input$lEACratioAdiff), as.numeric(input$uEACratioAdiff))
|
|
Dat$limite <- Limite
|
|
|
|
|
|
tab <- tests_FUNC(sim2(), Limite, PureErrFlag = PureErrFlag)
|
|
|
|
tab[1,6:7] <- c("-","-")
|
|
Dat$tests_FUNC <- tab
|
|
REP$testsTab <- tab
|
|
tab2 <- tab[1:7,]
|
|
|
|
dat <- datatable(tab2,options = list(
|
|
paging=TRUE,
|
|
dom="t",
|
|
rownames=FALSE
|
|
)) %>% formatStyle("test_results",
|
|
target='row',
|
|
backgroundColor = styleEqual(c(-1,0,1),
|
|
c("pink",'lightgreen','lightgrey')))
|
|
browser()
|
|
output$secondary <- DT::renderDataTable({ dat})
|
|
|
|
output$EQtests <- DT::renderDataTable( dat, options = list(lengthChange=FALSE) )
|
|
|
|
}) # observe
|
|
|
|
####plot CIs ----
|
|
observe({
|
|
tab <- Dat$tests_FUNC
|
|
if (is.null(tab)) return(NULL)
|
|
|
|
tab2 <- tab[-c(1,2,3,6),]
|
|
tab2[,3:7] <- apply(tab2[,3:7],2,as.numeric)
|
|
tab2[4:5,3:7] <- tab2[4:5,3:7]/100
|
|
|
|
p_CIs <- ggplot(tab2,aes(x=test,y=estimate, color=test,group=test)) +
|
|
geom_point() +
|
|
geom_errorbar(aes(ymin=lower_CI, ymax=upper_CI), width=0.4) +
|
|
geom_crossbar(aes(ymin=lower_limit, ymax=upper_limit), size=0.8) +
|
|
coord_flip() +
|
|
theme_bw() +
|
|
theme(legend.position = "none",text = element_text(size=20))
|
|
|
|
output$CIplot <- renderPlot({ p_CIs}, height=200)
|
|
|
|
REP$CIplot <- p_CIs
|
|
})
|
|
|
|
output$simdat <- DT::renderDataTable({
|
|
tab <- sim2()
|
|
if (is.character(tab)) stop(tab)
|
|
|
|
tab2 <- round(tab, 5)
|
|
colnames(tab2) <- c(paste("T", seq(1,(ncol(tab2)-1)/2)),
|
|
paste("R", seq(1,(ncol(tab2)-1)/2)), "log_conc" )
|
|
dat <- datatable(tab2, options=list(
|
|
paging=T,
|
|
pageLength=20,
|
|
dom="t"
|
|
))
|
|
})
|
|
|
|
output$Conctab <- DT::renderDataTable({
|
|
if (!is.na(Dils()[1]) & is.na(Dils()[4])) return(NULL)
|
|
tab <- sim2()
|
|
if (is.character(tab)) stop(tab)
|
|
if (!is.na(Dils()[4])) {
|
|
noDilSer <- Dils()[4]
|
|
} else if (!is.null(Dat$noDilSeriesXL)) {
|
|
noDilSer <- Dat$noDilSeriesXL
|
|
} else { noDilSer <- 3 }
|
|
|
|
Conc <- CONC()
|
|
Conctab <- perConcTab(tab, noDilSeries = noDilSer)
|
|
Dat$Conctab <- Conctab
|
|
|
|
dat <- datatable(Conctab, options=list(
|
|
paging=T,
|
|
pageLength=12,
|
|
dom="t"
|
|
)) %>% formatStyle(0,
|
|
target='row',
|
|
backgroundColor = styleEqual(c("avs","sds","cv", "avs_test","sds_test","cv_test"),
|
|
c('lightgrey','lightgreen','pink','lightgrey','lightgreen','pink'))
|
|
) %>% formatRound(columns=colnames(Conctab), digits=3)
|
|
})
|
|
|
|
#### linear Plot output ----
|
|
|
|
output$plotLin <- renderPlot({
|
|
tab <- sim2()
|
|
if (is.character(tab)) stop(tab)
|
|
|
|
if (!is.na(Dils()[4])) noDilSer <- Dils()[4] else noDilSer = (ncol(tab)-1)/2
|
|
Conc <- CONC()
|
|
Conctab <- Dat$Conctab
|
|
if (!is.na(Dils()[3])) noDil <- Dils()[3] else noDil = length(Conc)
|
|
|
|
slopeSt <- slopeTe <- matrix(NA, nrow=noDil-2,ncol=2)
|
|
for (i in 1:(noDil-2)) {
|
|
avs <- Conctab[noDilSer+1,]
|
|
threes <- data.frame(lnC=log(Conc[i:(i+2)]), resp=avs[i:(i+2)])
|
|
lm3St <- lm(resp ~ lnC, data=threes)
|
|
slopeSt[i,] <- lm3St$coefficients
|
|
avt <- Conctab[noDilSer*2+4,]
|
|
threet <- data.frame(lnC=log(Conc[i:(i+2)]), resp=avt[i:(i+2)])
|
|
lm3Te <- lm(resp ~ lnC, data=threet)
|
|
slopeTe[i,] <- lm3Te$coefficients
|
|
}
|
|
|
|
indS <- which(abs(slopeSt[,2]) == max(abs(slopeSt[,2])))
|
|
indT <- which(abs(slopeTe[,2]) == max(abs(slopeTe[,2])))
|
|
|
|
pl_ <- slopeSt[indS,1]+slopeSt[indS,2]*log(Conc)
|
|
pl_T <- slopeTe[indT,1]+slopeTe[indT,2]*log(Conc)
|
|
pl_df <- data.frame(lnC=log(Conc), plotS=pl_, plotT=pl_T)
|
|
|
|
all_l <- melt(data.frame(tab), 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)
|
|
all_l2S <- all_l2[all_l2$isRef == 1,]
|
|
all_l2T <- all_l2[all_l2$isRef == 0,]
|
|
all_mS <- all_l2S[order(all_l2S$log_dose, decreasing=TRUE),]
|
|
all_mT <- all_l2T[order(all_l2T$log_dose, decreasing=TRUE),]
|
|
|
|
circleS <- all_mS[(indS*noDilSer-(noDilSer-1)):((indS+2)*noDilSer),]
|
|
circleT <- all_mT[(indT*noDilSer-(noDilSer-1)):((indT+2)*noDilSer),]
|
|
circle <- rbind(circleS,circleT)
|
|
|
|
Dat$circles <- circle
|
|
sigmoid <- sigmoid()
|
|
log_dose <- unique(all_l$log_dose)
|
|
seq_x <- seq(min(log_dose), max(log_dose),0.1)
|
|
SAMPLEtrue <- sigmoid[2] + (sigmoid[4]-sigmoid[2])/(1+exp(sigmoid[6]*((sigmoid[7]-log(input$potencydiff/100)-seq_x))))
|
|
REFtrue <- sigmoid[1] + (sigmoid[3]-sigmoid[1])/(1+exp(sigmoid[5]*((sigmoid[7]-seq_x))))
|
|
|
|
truePL_df <- cbind(seq_x,SAMPLEtrue, REFtrue)
|
|
|
|
p <- ggplot(all_l2,aes(x=log_dose,y=readout, color=factor(isRef))) +
|
|
geom_point() +
|
|
labs(title=paste("linear regression model", indS,indT), color="product") +
|
|
scale_colour_manual(labels = c("test","reference"), values=c("red","blue")) +
|
|
ylim(min(all_l2$readout),max(all_l2$readout)) +
|
|
theme_bw()
|
|
p2 <- p + geom_line(data=pl_df,aes(x=lnC,y=plotS),color="blue",
|
|
inherit.aes = F) +
|
|
geom_line(data=pl_df,aes(x=lnC,y=plotT),color="red",
|
|
inherit.aes = F) +
|
|
geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=SAMPLEtrue),color="red", linetype=2,alpha=0.4,
|
|
inherit.aes = F) +
|
|
geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=REFtrue),color="blue", linetype=2,alpha=0.4,
|
|
inherit.aes = F) +
|
|
labs(title = paste("unrestricted linear regression model",indS,indT), color="product") +
|
|
theme(legend.position="none", axis.text = element_text(size=14))
|
|
p3 <- p2 + geom_point(circle, mapping=aes(x=log_dose, y=readout, shape=factor(isRef),
|
|
size=5,alpha=0.2), inherit.aes = FALSE) +
|
|
scale_shape_manual(labels=c("test","reference"), values=c(21,21))
|
|
|
|
mLin <- gsl_nls(readout ~ (intS+r)*isSample + intS*isRef + k*log_dose,
|
|
data=circle,
|
|
start=list(intS = 0, k=1,r=0),
|
|
control = gsl_nls_control(xtol=1e-10,ftol=1e-10,gtol=1e-10))
|
|
# alternativ: modAB <- lm(readout ~ log_dose+isSample, circle)
|
|
sum_mLin <- summary(mLin)
|
|
|
|
pl_restS <- sum_mLin$coefficients[1,1] + sum_mLin$coefficients[2,1]*log(Conc)
|
|
pl_restT <- sum_mLin$coefficients[1,1] + sum_mLin$coefficients[3,1] + sum_mLin$coefficients[2,1]*log(Conc)
|
|
pl_rest <- data.frame(lnC=log(Conc), plotS=pl_restS, plotT=pl_restT)
|
|
|
|
pr2 <- p + geom_line(data=pl_rest,aes(x=lnC,y=plotS),color="blue",
|
|
inherit.aes = F) +
|
|
geom_line(data=pl_rest,aes(x=lnC,y=plotT),color="red",
|
|
inherit.aes = F) +
|
|
geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=SAMPLEtrue),color="red", linetype=2,alpha=0.4,
|
|
inherit.aes = F) +
|
|
geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=REFtrue),color="blue", linetype=2,alpha=0.4,
|
|
inherit.aes = F) +
|
|
labs(title = paste("restricted linear regression model",indS,indT), color="product") +
|
|
theme(legend.position="none", axis.text = element_text(size=14))
|
|
pr3 <- pr2 + geom_point(circle, mapping=aes(x=log_dose, y=readout, shape=factor(isRef),
|
|
size=5,alpha=0.2), inherit.aes = FALSE) +
|
|
scale_shape_manual(labels=c("test","reference"), values=c(21,21))
|
|
grid.arrange(p3,pr3,nrow=1)
|
|
})
|
|
|
|
#### linear PLA tests ----
|
|
output$TESTSlin <- DT::renderDataTable({
|
|
tab <- sim2()
|
|
if (is.character(tab)) stop(tab)
|
|
Conc <- CONC()
|
|
Limite <- list(as.numeric(input$lEACdiffla), as.numeric(input$uEACdiffla),
|
|
as.numeric(input$lEACratiola), as.numeric(input$uEACratiola),
|
|
as.numeric(input$lEACratioSlope), as.numeric(input$uEACratioSlope),
|
|
as.numeric(input$lEACratioua), as.numeric(input$uEACratioua),
|
|
as.numeric(input$lowerPot), as.numeric(input$upperPot),
|
|
as.numeric(input$lEACratioAdiff), as.numeric(input$uEACratioAdiff))
|
|
|
|
circles <- Dat$circles
|
|
PureErrFlag <- input$PureErr
|
|
warning_text <- reactive({
|
|
ifelse(PureErrFlag, 'Pure error is selected','')
|
|
})
|
|
output$PureErrW <- renderText(warning_text())
|
|
|
|
LIN <- ANOVAlintests(tab,circles,Limite,PureErrFlag=PureErrFlag)
|
|
df <- LIN[[1]]
|
|
su_modU <- LIN[[2]]
|
|
su_mod2 <- LIN[[4]]
|
|
|
|
output$SummaryModABu <- renderTable({ su_modU }, digits=5)
|
|
output$SummaryModAB <- renderTable({ su_mod2 }, digits=5)
|
|
|
|
slopeDiffCI <- t(data.frame(LIN[[3]]))
|
|
colnames(slopeDiffCI) <- c("slope difference","lower CI","upper CI")
|
|
output$SlopeDiffCI <- renderTable({ slopeDiffCI },digits=5)
|
|
#browser()
|
|
Dat$ANOVA <- df[,4:length(df)]
|
|
dat <- datatable(df[,1:3],
|
|
options=list(
|
|
paging=T, dom="t",rownames=F
|
|
)) %>% formatStyle("test_results", target="row",backgroundColor = styleEqual(c(-1,0,1),
|
|
c("pink","lightgreen","lightgrey")))
|
|
|
|
})
|
|
|
|
#### output 4PL ANOVA tests ---
|
|
output$ANOVA <- DT::renderDataTable({
|
|
sigmoid <- sigmoid()
|
|
tab <- ANOVA4plUnresfunc(sim2(),sigmid)
|
|
dat <- datatable(df[,1:3],
|
|
options=list(
|
|
dom="t",rownames=F
|
|
)) %>% formatStyle("p.value", target="row",
|
|
backgroundColor = styleEqual(c("p.value"),
|
|
c("lightgrey")))
|
|
})
|
|
|
|
#### output RMSEs ----
|
|
output$RMSE <- renderText({
|
|
paste("RMSE (unrestricted model):", Dat$RMSE_unr, "(~ entered % upper-lower asymptote)\n",
|
|
"RMSE restricted model:", Dat$RMSE_r, "\n",
|
|
"Pure RMSE unrestricted model:", Dat$RMSE_pure, "\n",
|
|
"%SD (unr model): ", Dat$RMSE_unr*100/Dat$up_lowAs, "(calculated as: RMSE/(upper-lower Asymptote)*100\n",
|
|
"RMSE (log restr. model): ", Dat$RMSE_Rlog, "\n",
|
|
"RMSE (log unrestr. model): ", Dat$RMSE_Ulog, "\n",
|
|
"%SDlog (unr model): ", Dat$RMSE_Ulog*100/Dat$up_lowAslog )
|
|
})
|
|
|
|
output$ANOVAlin <- DT::renderDataTable({
|
|
ANOVAlin <- Dat$ANOVA
|
|
dat <- datatable(ANOVAlin,
|
|
options=list(
|
|
dom="t",rownames=F
|
|
)) %>% formatStyle("p.value", target='cell',
|
|
backgroundColor = styleEqual(c("p.value"),
|
|
c("lightgrey")))
|
|
})
|
|
### output pot tab ----
|
|
output$pottab <- DT::renderDataTable({
|
|
|
|
Lim <- list(as.numeric(input$lEACdiffla), as.numeric(input$uEACdiffla),
|
|
as.numeric(input$lEACratiola), as.numeric(input$uEACratiola),
|
|
as.numeric(input$lEACratioSlope), as.numeric(input$uEACratioSlope),
|
|
as.numeric(input$lEACratioua), as.numeric(input$uEACratioua),
|
|
as.numeric(input$lowerPot), as.numeric(input$upperPot))
|
|
|
|
circles <- Dat$circles
|
|
PureErrFlag <- input$PureErr
|
|
pottab <- LinPotTab(circles,Lim,PureErrFlag = PureErrFlag)
|
|
#browser()
|
|
dat <- datatable(pottab,
|
|
options=list(
|
|
dom="t",rownames=F
|
|
)) %>% formatStyle("test_result", target='row',
|
|
backgroundColor = styleEqual(c(0,1), c("lightgrey")))
|
|
})
|
|
|
|
#### 4pl potency table ----
|
|
output$pottab4pl <- DT::renderDataTable({
|
|
ro_new <- sim2()
|
|
Dils_ <- Dils()
|
|
if (!is.na(Dils()[4])) noDilSer <- Dils()[4] else noDilSer <- 3
|
|
PureErrFl <- input$PureErr
|
|
pottab4 <- pot4plFUNC(ro_new = ro_new, PureErrFlag = PureErrFl)
|
|
Lim <- list(as.numeric(input$lEACdiffla), as.numeric(input$uEACdiffla),
|
|
as.numeric(input$lEACratiola), as.numeric(input$uEACratiola),
|
|
as.numeric(input$lEACratioSlope), as.numeric(input$uEACratioSlope),
|
|
as.numeric(input$lEACratioua), as.numeric(input$uEACratioua),
|
|
as.numeric(input$lowerPot), as.numeric(input$upperPot))
|
|
if (as.numeric(pottab4[1,3])*100>Lim[[9]] & as.numeric(pottab4[1,4])*100<Lim[[10]] ) {
|
|
test_potCI <- 0
|
|
} else {test_potCI <- 1 }
|
|
if (as.numeric(pottab4[2,3])*100>Lim[[9]] & as.numeric(pottab4[2,4])*100<Lim[[10]] ) {
|
|
test_potUCI <- 0
|
|
} else {test_potUCI <- 1 }
|
|
if (as.numeric(pottab4[3,3])*100>Lim[[9]] & as.numeric(pottab4[3,4])*100<Lim[[10]] ) {
|
|
test_potCI_t <- 0
|
|
} else {test_potCI_t <- 1 }
|
|
if (as.numeric(pottab4[4,3])*100>Lim[[9]] & as.numeric(pottab4[4,4])*100<Lim[[10]] ) {
|
|
test_potUCI_t <- 0
|
|
} else {test_potUCI_t <- 1 }
|
|
|
|
pottab4_ <- cbind(pottab4, data.frame(tests=c(test_potCI, test_potUCI,test_potCI_t,test_potUCI_t)))
|
|
colnames(pottab4_) <- c("model","potency","lower95%CI","upper95%CI","test_result")
|
|
pottab4_$potency <- round(as.numeric(pottab4_$potency)*100,3)
|
|
pottab4_$`lower95%CI` <- round(as.numeric(pottab4_$`lower95%CI`)*100,3)
|
|
pottab4_$`upper95%CI` <- round(as.numeric(pottab4_$`upper95%CI`)*100,3)
|
|
pottab4_$relative_lowerCL <- round(pottab4_[,3]/pottab4_[,2]*100,3)
|
|
pottab4_$relative_upperCL <- round(pottab4_[,4]/pottab4_[,2]*100,3)
|
|
dat <- datatable(pottab4_,
|
|
options=list(
|
|
paging=T, dom="t",rownames=F
|
|
)) %>% formatStyle("test_result", target="row",backgroundColor = styleEqual(c(0,1),
|
|
c("lightgreen","pink")))
|
|
|
|
})
|
|
|
|
#### Dilutions Simulator ----
|
|
output$plotfordilutions <- renderPlot({
|
|
tab <- sim()
|
|
|
|
tab <- as.data.frame(tab)
|
|
dils <- tab$log_dose
|
|
min_y <- min(tab[,1:3])
|
|
max_y <- max(tab[,1:3])
|
|
|
|
if (input$fixupper) {
|
|
dils_av <- dils-max(dils)
|
|
dils_av_ <- dils_av*(input$dilslider/100+1)
|
|
dils2 <- round(dils_av_ + max(dils),4)
|
|
dilfactors <- 1/exp(dils2-lag(dils2))
|
|
} else {
|
|
if (!is.null(Dat$cfordils)) {
|
|
av <- Dat$cfordils
|
|
} else { av <- (min(dils) + max(dils))/2 }
|
|
dils_av <- dils-av
|
|
dils_avsc <- dils_av*(input$dilslider/100+1)
|
|
dils2 <- dils_avsc+av
|
|
dilfactors <- 1/exp(dils2-lag(dils2))
|
|
}
|
|
|
|
Dat$newDils <- dils2
|
|
|
|
sigmoid <- sigmoid()
|
|
|
|
#browser()
|
|
BPs <- Dat$bendpoints
|
|
EC50REF <- (BPs[2]+BPs[1])/2
|
|
Einh <- abs((BPs[2]-BPs[1])/5)
|
|
asyml <- EC50REF-2*(EC50REF-BPs[1])
|
|
asymu <- EC50REF+2*(EC50REF-BPs[1])
|
|
|
|
det_sig <- Dat$coeffs_UN
|
|
|
|
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]
|
|
Xbend50 <- max(Xbend50l, Xbend50u)
|
|
Xbend200 <- min(Xbend200l, Xbend200u)
|
|
dummy <- plot_f(tab,sigmoid,det_sig=NULL)
|
|
} else {
|
|
|
|
#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]
|
|
Xbend50 <- max(Xbend50l, Xbend50u)
|
|
Xbend200 <- min(Xbend200l, Xbend200u)
|
|
dummy <- plot_f(tab,sigmoid=NULL,det_sig=det_sig)
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
pl_df <- cbind(dils2, SAMPLE50, SAMPLE200)
|
|
|
|
#browser()
|
|
# scenario2
|
|
eqSpac <- abs((BPs[1]-BPs[2])/5)
|
|
optdils <- c((asyml+BPs[1])/2, BPs[1], BPs[1]+1*eqSpac, BPs[1]+2*eqSpac,BPs[1]+3*eqSpac,BPs[1]+4*eqSpac,BPs[2], (asymu+BPs[2])/2)
|
|
# scenario 3
|
|
eqSpac_3 <- abs((BPs[1]-BPs[2])/3)
|
|
optdils_3 <- c(BPs[1]-2*eqSpac_3, BPs[1]-eqSpac_3, BPs[1], BPs[1]+1*eqSpac_3, BPs[1]+2*eqSpac_3,BPs[2], BPs[2]+eqSpac_3, BPs[2]+2*eqSpac_3)
|
|
# scenario 6
|
|
Einh2 <- abs(((BPs[2]-BPs[1])*0.7)/5)
|
|
eqSpac2 <- (2*0.7/Einh)/3
|
|
optdils2 <- c((asyml+BPs[1])/2, BPs[1], EC50REF-1.5*Einh2, EC50REF-0.5*Einh2,EC50REF+0.5*Einh2,EC50REF+1.5*Einh2, BPs[2], (asymu+BPs[2])/2)
|
|
# steep slope
|
|
eqSpac3 <- (abs(Xbend200-Xbend50))/5
|
|
optdils3 <- c(Xbend200-eqSpac3,Xbend200, Xbend200+1*eqSpac3, Xbend200+2*eqSpac3,Xbend200+3*eqSpac3,Xbend200+4*eqSpac3,Xbend50, Xbend50+eqSpac3)
|
|
|
|
output$extremebps <- renderTable({
|
|
ExtremeBPs <- c(Xbend50,Xbend200)
|
|
DF2 <- data.frame(sample=c("50% sample (right)", "200% sample (left)"), Extreme_BPs=ExtremeBPs)
|
|
DF2
|
|
})
|
|
|
|
optD <- data.frame(cbind(optdils, optdils_3,optdils2, optdils3))
|
|
colnames(optD) <- c("scenario2","scenario3","scenario6","steep slope")
|
|
|
|
output$optimalDils <- renderTable({ optD })
|
|
|
|
output$adjlogdil <- renderTable({
|
|
adjlogdilfactors <- round(dilfactors,3)
|
|
adjlogdils <- round(dils2,3)
|
|
adjdils <- round(exp(dils2),3)
|
|
DilsTable <- data.frame('adjusted ln(dilutions)' = adjlogdils,
|
|
'adjusted ln_dilution_factors' = adjlogdilfactors,
|
|
'adjusted dilutions' = adjdils)
|
|
DilsTable
|
|
})
|
|
|
|
if (!is.null(Dat$p2)) {
|
|
p2 <- Dat$p2
|
|
p_dil <- p2 +
|
|
annotate("pointrange",x=dils2,y=rep(min_y, length(dils2)), xmin=min(dils2), xmax=max(dils2)) +
|
|
annotate("text", x=dils2,y=rep(min_y+(max_y-min_y)*0.05, length(dils2)), label=as.character(round(dils2,3))) +
|
|
annotate("text", x=dils2[-1]+(max(dils2)-min(dils2))*0.05,
|
|
y=rep(min_y+(max_y-min_y)*0.1, length(dils2[-1])),
|
|
label=as.character(round(dilfactors[-1],3))) +
|
|
geom_line(data=as.data.frame(pl_df),aes(x=dils2,y=SAMPLE50), color="grey15", linetype=2,
|
|
inherit.aes = F) +
|
|
geom_line(data=as.data.frame(pl_df),aes(x=dils2,y=SAMPLE200), color="grey15", linetype=2,
|
|
inherit.aes = F) +
|
|
geom_vline(xintercept=c(Xbend50,Xbend200), col="grey15", linetype=2) +
|
|
{if (input$scenario =="scenario 6") annotate("pointrange",x=optdils2,y=rep(min_y+(max_y-min_y)*0.2, length(optdils2)),
|
|
xmin=min(optdils2), xmax=max(optdils2), color="seagreen")} +
|
|
{if (input$scenario =="scenario 6") annotate("text",x=optdils2,y=rep(min_y+(max_y-min_y)*0.25, length(optdils2)),
|
|
label=as.character(round(optdils2,3)), color="seagreen")} +
|
|
{if (input$scenario =="scenario 2") annotate("pointrange",x=optdils,y=rep(min_y+(max_y-min_y)*0.2, length(optdils)),
|
|
xmin=min(optdils), xmax=max(optdils), color="seagreen")} +
|
|
{if (input$scenario =="scenario 2") annotate("text",x=optdils,y=rep(min_y+(max_y-min_y)*0.25, length(optdils)),
|
|
label=as.character(round(optdils,3)), color="seagreen")} +
|
|
{if (input$scenario =="scenario 3") annotate("pointrange",x=optdils_3,y=rep(min_y+(max_y-min_y)*0.2, length(optdils_3)),
|
|
xmin=min(optdils_3), xmax=max(optdils_3), color="seagreen")} +
|
|
{if (input$scenario =="scenario 3") annotate("text",x=optdils_3,y=rep(min_y+(max_y-min_y)*0.25, length(optdils_3)),
|
|
label=as.character(round(optdils_3,3)), color="seagreen")} +
|
|
{if (input$scenario =="steep slope") annotate("pointrange",x=optdils3,y=rep(min_y+(max_y-min_y)*0.2, length(optdils3)),
|
|
xmin=min(optdils3), xmax=max(optdils3), color="seagreen")} +
|
|
{if (input$scenario =="steep slope") annotate("text",x=optdils3,y=rep(min_y+(max_y-min_y)*0.25, length(optdils3)),
|
|
label=as.character(round(optdils3,3)), color="seagreen")} +
|
|
annotate("text",x=optdils[1],y=(max_y+min_y)*0.5,
|
|
label=paste("in green: optimal \n dilutions acc. to Whitepaper\n", input$scenario), color="seagreen",
|
|
size=14/.pt,fontface="bold")
|
|
}
|
|
print(p_dil)
|
|
|
|
})
|
|
|
|
#### Dilutions CI table ----
|
|
observe({
|
|
|
|
output$CIs <- renderTable({
|
|
PureErrFlag <- input$PureErr
|
|
if (is.null(Dat$coeffs_UN)) {
|
|
# checks if an EXCEL was uploaded
|
|
sigmoid <- sigmoid()
|
|
det_sig=NULL
|
|
|
|
ast = sigmoid()[1];bst = sigmoid()[5];cst = sigmoid()[7];dst = sigmoid()[3];ate = sigmoid()[2];
|
|
bte = sigmoid()[6];r_ = sigmoid()[8];
|
|
cte = cst-r_;dte = sigmoid()[4];
|
|
} else {
|
|
sigmoid <- NULL
|
|
det_sig <- Dat$coeffs_UN
|
|
ast <- det_sig[3]
|
|
ate <- det_sig[4]
|
|
bst <- det_sig[1]
|
|
bte <- det_sig[2]
|
|
cst <- det_sig[7]
|
|
cte <- det_sig[7] -log(input$potencydiff/100)
|
|
dst <- det_sig[5]
|
|
dte <- det_sig[6]
|
|
r_ <- log(input$potencydiff/100)
|
|
|
|
}
|
|
if (!is.na(input$NoDilSer)) {
|
|
noDilSer <- input$NoDilSer
|
|
} else if (!is.null(Dat$NoDilSeriesXL)) noDilSer <- Dat$noDilSeriesXL else noDilSer <- 3
|
|
if (!is.na(input$NoDil)) noDil <- input$NoDil else noDil <- length(Dat$newDils)
|
|
#browser()
|
|
tab <- Calc_DilRes(as=ast,at=ate,ds=dst,dt=dte,cs=cst,ct=cte,r=r_,bt=bte,bs=bst,
|
|
sd_fac=input$sdfac,log_conc=Dat$newDils,
|
|
# auslenkU=outlierU,
|
|
# auslenkM=outlierM,
|
|
# auslenkL=outlierL,
|
|
heteroNoise = FALSE, noDilSeries = noDilSer, noDils = noDil)
|
|
Limite <- list(as.numeric(input$lEACdiffla), as.numeric(input$uEACdiffla),
|
|
as.numeric(input$lEACratiola), as.numeric(input$uEACratiola),
|
|
as.numeric(input$lEACratioSlope), as.numeric(input$uEACratioSlope),
|
|
as.numeric(input$lEACratioua), as.numeric(input$uEACratioua),
|
|
as.numeric(input$lowerPot), as.numeric(input$upperPot),
|
|
as.numeric(input$lEACratioAdiff), as.numeric(input$uEACratioAdiff))
|
|
|
|
CItable <- tests_FUNC(tab,Limite,PureErrFlag=PureErrFlag)
|
|
|
|
CItable_ <- CItable[-c(1,2,6,8,9),-c(2,4,5)]
|
|
potAll <- pot4plFUNC(tab, input$PureErr)
|
|
restrPot <- potAll[1,1:4]
|
|
restrPot[2:4] <- round(as.numeric(restrPot[2:4]),5)
|
|
potAll_ <- rbind(CItable_, restrPot)
|
|
potAll_$CIwidth <- as.numeric(potAll_[,4])-as.numeric(potAll_[,3])
|
|
potAll_[,1] <- c("ratio of lower asymptotes","ratio of slopes","ratio of upper asymptotes", "ratio of asympt. differences","restricted potency")
|
|
|
|
output$bps <- renderTable({
|
|
DF <- data.frame(sample=names(Dat$bendpoints),BPs=Dat$bendpoints)
|
|
DF
|
|
})
|
|
return(potAll_)
|
|
})
|
|
})
|
|
|
|
#### simulations ----
|
|
observe({
|
|
observeEvent(input$goSim,{
|
|
sd_fac_ <- as.numeric(input$sdfac)
|
|
|
|
r_ <- log(as.numeric(input$potencydiff)/100)
|
|
Conc <- Dat$MetaConc
|
|
as = sigmoid()[1]; bs = sigmoid()[5];cs = sigmoid()[7];ds = sigmoid()[3];at = sigmoid()[2];
|
|
bt = sigmoid()[6];r = sigmoid()[8]; ct = cs-r_; dt = sigmoid()[4]
|
|
|
|
if (!is.null(Dat$MetaConc)) {
|
|
Conc <- Dat$MetaConc
|
|
} else {
|
|
Conc <- CONC()
|
|
}
|
|
log_dose <- log(Conc)
|
|
yAxfac <- (ds-as)
|
|
|
|
if (!is.na(input$NoDilSer)) {
|
|
noDilSer <- input$NoDilSer
|
|
} else if (!is.null(Dat$NoDilSeriesXL)) noDilSer <- Dat$noDilSeriesXL else noDilSer <- 3
|
|
if (!is.na(input$NoDil)) noDil <- input$NoDil else noDil <- length(Conc)
|
|
isRef <- rep(c(1,0),1,each=noDilSer*noDil)
|
|
isSample <- rep(c(0,1),1,each=noDilSer*noDil)
|
|
N <- as.numeric(input$simN)
|
|
|
|
av <- as*isRef + at*isSample + (ds*isRef + dt*isSample - as*isRef - at*isSample)/
|
|
(1+isRef*exp(bs*(cs - log_dose)) + isSample*exp(bt*(ct-log_dose)))
|
|
|
|
resHist <- matrix(NA,nrow=N, ncol=13)
|
|
residualsList <- list()
|
|
start.time2 <- Sys.time()
|
|
withProgress(message = 'Making plot', value=0, {
|
|
for (i in 1:N) {
|
|
if (input$heterosked) {
|
|
# heterosc noise
|
|
ro_jit <- matrix(unlist(map(av, function(x) x+rnorm(1,0,x*sd_fac_/100))), nrow=noDil, ncol=noDilSer*2)
|
|
} else {
|
|
# homosc noise
|
|
ro_jit <- matrix(unlist(map(av, function(x) x+rnorm(1,0,sd_fac_*yAxfac/100))), nrow=noDil, ncol=noDilSer*2)
|
|
}
|
|
# browser()
|
|
ro_jit <- abs(ro_jit)
|
|
ro_new <- cbind(ro_jit, log_dose)
|
|
all_l <- melt(data.frame(ro_new), id.vars="log_dose", variable.name = "replname", value.name = "readout")
|
|
all_l$isRef <- isRef
|
|
all_l$isSample <- isSample
|
|
all_l$Conc <- exp(all_l$log_dose)
|
|
pot <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")),
|
|
pmodels=data.frame(1,1,1,isSample))
|
|
potAll <- EDcomp(pot, percVec=c(50,50), interval="delta", display=FALSE)
|
|
potAll2 <- potAll[1:3]
|
|
RSS <- sum(pot$predres[,2]^2)
|
|
dfreed <- nrow(all_l)-5
|
|
MSE <- RSS/dfreed
|
|
|
|
potU <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")),
|
|
pmodels=data.frame(isSample, isSample,isSample,isSample))
|
|
DF_U <- nrow(all_l)-8
|
|
|
|
uAsratio <- compParm(potU, "a",display=F)
|
|
uCIuAs <- uAsratio[1]+qt(0.975,DF_U)*uAsratio[2]
|
|
lCIuAs <- uAsratio[1]-qt(0.975,DF_U)*uAsratio[2]
|
|
lAsratio <- compParm(potU, "d",display=F)
|
|
uCIlAs <- lAsratio[1]+qt(0.975,DF_U)*lAsratio[2]
|
|
lCIlAs <- lAsratio[1]-qt(0.975,DF_U)*lAsratio[2]
|
|
Sloperatio <- compParm(potU, "b",display=F)
|
|
uCISlo <- Sloperatio[1]+qt(0.975,DF_U)*Sloperatio[2]
|
|
lCISlo <- Sloperatio[1]-qt(0.975,DF_U)*Sloperatio[2]
|
|
su <- summary(potU)
|
|
|
|
v <- vcov(potU)[c(5,6),c(5,6)]
|
|
Vd <- vcov(potU)[c(3,4),c(3,4)]
|
|
Va_d <- v+Vd
|
|
|
|
A_DTEST <- su$coefficients[6,1]-su$coefficients[4,1]
|
|
A_DREF <- su$coefficients[5,1]-su$coefficients[3,1]
|
|
if (abs(at/(sqrt(Va_d[2,2]/3))) > qt(0.95,2)) {
|
|
try(Fie_ad <- round(FiellerRatio(A_DREF,A_DTEST, Va_d),5))
|
|
}
|
|
if (!exists("Fie_ad")) Fie_ad <- NA
|
|
|
|
resHist[i,] <- c(potAll2, sqrt(MSE),Sloperatio[1],lCISlo, uCISlo,
|
|
uAsratio[1], lCIuAs, uCIuAs, Fie_ad[1],Fie_ad[2],Fie_ad[3])
|
|
colnames(resHist) <- c("pot4pl","lCI4pl","uCI4pl","RMSE","estSlope_ratio",
|
|
"lCISlope_ratio","uCISlope_ratio","estuAs_ratio",
|
|
"lCIuAs_ratio","uCIuAs_ratio","estAsyDiff_ratio",
|
|
"lCIAsyDiff_ratio", "uCIAsyDiff_ratio")
|
|
|
|
incProgress(1/N, detail=paste("Doing simulations",i))
|
|
} # withProgress
|
|
|
|
})
|
|
end.time2 <- Sys.time()
|
|
|
|
Dat$resHist <- resHist
|
|
})
|
|
})
|
|
|
|
|
|
#### simulation Histograms output ----
|
|
|
|
output$plotHistuAs <- renderPlot({
|
|
if (!is.null(Dat$resHist)) {
|
|
|
|
resHist <- Dat$resHist
|
|
#browser()
|
|
resHistuAs <- as.data.frame(resHist[,8:10])
|
|
resHistuAs_l <- melt(data.frame(resHistuAs), variable.name="ratio_CIs", value.name = "readout")
|
|
#browser()
|
|
lowquant_uAs <- quantile(resHistuAs[,2], probs=as.numeric(input$lowQuant)/100)
|
|
upquant_uAs <- quantile(resHistuAs[,3], probs=as.numeric(input$uppQuant)/100)
|
|
|
|
p_uAs <- ggplot(resHistuAs_l) +
|
|
geom_histogram(aes(readout, fill=ratio_CIs),alpha=0.5,position="identity") +
|
|
labs(title = paste("upper asymptote ratio EACs:", round(lowquant_uAs,3), " to ", round(upquant_uAs,3))) +
|
|
geom_vline(xintercept = c(lowquant_uAs, upquant_uAs), color="black", linetype="dashed", linewidth=1) +
|
|
geom_vline(xintercept = c(input$lEACratioua , input$uEACratioua), color="red", linetype="dashed", linewidth=1) +
|
|
theme_bw()
|
|
|
|
# asympt diff ratio
|
|
resHistAsDiff <- as.data.frame(resHist[,11:13])
|
|
resHistAsDiff_l <- melt(data.frame(resHistAsDiff), variable.name="ratio_CIs", value.name = "readout")
|
|
|
|
lowquant_AsDiff <- quantile(resHistAsDiff[,2], probs=as.numeric(input$lowQuant)/100)
|
|
upquant_AsDiff <- quantile(resHistAsDiff[,3], probs=as.numeric(input$uppQuant)/100)
|
|
|
|
p_AsDiff <- ggplot(resHistAsDiff_l, aes(readout, fill=ratio_CIs)) +
|
|
geom_histogram(alpha=0.5,position="identity") +
|
|
labs(title = paste("asymptote diff. ratio EACs:", round(lowquant_AsDiff,3), " to ", round(upquant_AsDiff,3))) +
|
|
geom_vline(xintercept = c(lowquant_AsDiff, upquant_AsDiff), color="black", linetype="dashed", linewidth=1) +
|
|
geom_vline(xintercept = c(input$lEACratioAdiff , input$uEACratioAdiff), color="red", linetype="dashed", linewidth=1) +
|
|
theme_bw()
|
|
|
|
# Slope ratio
|
|
resHistSlo <- as.data.frame(resHist[,5:7])
|
|
resHistSlo_l <- melt(data.frame(resHistSlo), variable.name="ratio_CIs", value.name = "readout")
|
|
|
|
lowquant_Slo <- quantile(resHistSlo[,2], probs=as.numeric(input$lowQuant)/100)
|
|
upquant_Slo <- quantile(resHistSlo[,3], probs=as.numeric(input$uppQuant)/100)
|
|
|
|
p_Slo <- ggplot(resHistSlo_l, aes(readout, fill=ratio_CIs)) +
|
|
geom_histogram(alpha=0.5,position="identity") +
|
|
labs(title = paste("Slope ratio EACs:", round(lowquant_Slo,3), " to ", round(upquant_Slo,3))) +
|
|
geom_vline(xintercept = c(lowquant_Slo, upquant_Slo), color="black", linetype="dashed", linewidth=1) +
|
|
geom_vline(xintercept = c(input$lEACratioSlope , input$uEACratioSlope), color="red", linetype="dashed", linewidth=1) +
|
|
theme_bw()
|
|
|
|
# poency ratio
|
|
resHistPot <- as.data.frame(resHist[,1:3])
|
|
resHistPot_l <- melt(data.frame(resHistPot), variable.name="ratio_CIs", value.name = "readout")
|
|
|
|
lowquant_Pot <- quantile(resHistPot[,2], probs=as.numeric(input$lowQuant)/100)
|
|
upquant_Pot <- quantile(resHistPot[,3], probs=as.numeric(input$uppQuant)/100)
|
|
#browser()
|
|
p_Pot <- ggplot(resHistPot_l, aes(readout, fill=ratio_CIs)) +
|
|
geom_histogram(alpha=0.5,position="identity") +
|
|
labs(title = paste("Poency ratio EACs:", round(lowquant_Pot,3), " to ", round(upquant_Pot,3))) +
|
|
geom_vline(xintercept = c(lowquant_Pot, upquant_Pot), color="black", linetype="dashed", linewidth=1) +
|
|
geom_vline(xintercept = c(input$lowerPot/100, input$upperPot/100), color="red", linetype="dashed", linewidth=1) +
|
|
theme_bw()
|
|
|
|
grid.arrange(p_Slo, p_AsDiff, p_uAs, p_Pot, nrow=1)
|
|
|
|
}
|
|
})
|
|
|
|
#### download XL report----
|
|
|
|
output$downloadXLReport <- downloadHandler(
|
|
filename= paste0("Report_4PLEvaluation", input$iFile,".pdf"),
|
|
|
|
content = function(file) {
|
|
tpdr <- tempdir()
|
|
tempReport <- file.path(tpdr,"Doc_BioassayReport.Rmd")
|
|
file.copy("Doc_BioassayReport.Rmd", tempReport, overwrite = T)
|
|
|
|
tempReportc <- file.path(tpdr,"logo.png")
|
|
file.copy("logo.png", tempReportc, overwrite = T)
|
|
|
|
rmarkdown::render(tempReport, output_file = file,
|
|
params = list(FileName = Dat$FileName,
|
|
author = Dat$author,
|
|
REP = REP,
|
|
coeffs = Dat$coeffs_UN),
|
|
envir = new.env(parent = globalenv()))
|
|
|
|
}
|
|
)
|
|
|
|
|
|
# observe({
|
|
# pfad <- getwd()
|
|
# Z <- read.csv('./Zaehlfile.csv')
|
|
# Z[,3] <- sapply(Z[,3], as.character)
|
|
# nrowsZ <- nrow(Z)
|
|
#
|
|
# Z[nrowsZ+1,1] <- Z[nrowsZ,1]+1
|
|
# Z[nrowsZ+1,1] <- Sys.Date()
|
|
# Z[nrowsZ+1,1] <- NA
|
|
# write.csv(Z, file= './Zaahlfile.csv', row.names = F)
|
|
# })
|
|
# observe({
|
|
# Z <- read.csv('./Zaehlfile.csv')
|
|
# Dat$Zaehl <- read.csv('./Zaehlfile.csv')
|
|
# nrowsZ <- nrow(Z)
|
|
# message2 <- paste("No. of calls:", Z[nrowsZ,1])
|
|
# output$stats <- renderPrint(cat(message2))
|
|
# })
|
|
# output$downloadZaehl <- downloadHandler(
|
|
# filename = function() {
|
|
# paste("Zaehlfile", Sys.Date(),".csv", sep="")
|
|
#
|
|
# },
|
|
# content = function(file) { write.csv(Dat$Zaehl, file = file, row.names = F)},
|
|
# contentType = "text/csv"
|
|
# )
|
|
} # end of session function
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|