4PL report update

This commit is contained in:
2026-05-14 18:38:27 +02:00
parent 9861af5fba
commit 9422490f25
4 changed files with 338 additions and 455 deletions
+192 -191
View File
@@ -90,198 +90,199 @@ You can also embed plots, for example:
```{r XLplot, echo=FALSE, warning=FALSE, fig.height=4, fig.width=6, fig.cap="Plot of models", fig.align='left'} ```{r XLplot, echo=FALSE, warning=FALSE, fig.height=4, fig.width=6, fig.cap="Plot of models", fig.align='left'}
plot_f <- function(dat, sigmoid,det_sig) { # plot_f <- function(dat, sigmoid,det_sig) {
CORdat <- cor(dat[,1],dat[,ncol(dat)]) # CORdat <- cor(dat[,1],dat[,ncol(dat)])
#
# 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)
#
# if(is.null(det_sig)) {
# if (CORdat<0) {
# startlist <- list(a=sigmoid[3], b=-sigmoid[5],cs=sigmoid[7],
# d=sigmoid[1],r=sigmoid[8])
# } else {
# startlist <- list(a=sigmoid[3],b=sigmoid[5],cs=sigmoid[7],
# d=sigmoid[1],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])
# }
# #browser()
# tryCatch({
# mr <- gsl_nls(fn = readout ~ 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))
# },
# error = function(err) {
# err$message
# })
# s_mr <- summary(mr)
# 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[4] + (sigmoid[2] -sigmoid[4])/(1+exp(sigmoid[6]*(seq_x-(sigmoid[7]-sigmoid[8]))))
# REFtrue <- sigmoid[3] + (sigmoid[1] -sigmoid[3])/(1+exp(sigmoid[5]*(seq_x-(sigmoid[7]))))
# } else {
# SAMPLEtrue <- det_sig[4] + (det_sig[6] -det_sig[4])/(1+exp(-det_sig[2]*(seq_x-(det_sig[8]))))
# REFtrue <- det_sig[3] + (det_sig[5] -det_sig[3])/(1+exp(-det_sig[1]*(seq_x-(det_sig[7]))))
# }
#
# 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))
#
# 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")
#
#
# # 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))
#
# 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)
# }
#
# plot_f(XLdat2, sigmoid=NULL, det_sig=coeffs)
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)
if(is.null(det_sig)) {
if (CORdat<0) {
startlist <- list(a=sigmoid[3], b=-sigmoid[5],cs=sigmoid[7],
d=sigmoid[1],r=sigmoid[8])
} else {
startlist <- list(a=sigmoid[3],b=sigmoid[5],cs=sigmoid[7],
d=sigmoid[1],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])
}
#browser()
tryCatch({
mr <- gsl_nls(fn = readout ~ 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))
},
error = function(err) {
err$message
})
s_mr <- summary(mr)
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[4] + (sigmoid[2] -sigmoid[4])/(1+exp(sigmoid[6]*(seq_x-(sigmoid[7]-sigmoid[8]))))
REFtrue <- sigmoid[3] + (sigmoid[1] -sigmoid[3])/(1+exp(sigmoid[5]*(seq_x-(sigmoid[7]))))
} else {
SAMPLEtrue <- det_sig[4] + (det_sig[6] -det_sig[4])/(1+exp(-det_sig[2]*(seq_x-(det_sig[8]))))
REFtrue <- det_sig[3] + (det_sig[5] -det_sig[3])/(1+exp(-det_sig[1]*(seq_x-(det_sig[7]))))
}
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))
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")
# 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))
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)
}
plot_f(XLdat2, sigmoid=NULL, det_sig=coeffs)
``` ```
+104 -237
View File
@@ -42,6 +42,7 @@ knitr::opts_chunk$set(echo = TRUE)
library(knitr) library(knitr)
library(DT) library(DT)
library(kableExtra)
REP <- params$REP REP <- params$REP
author <- params$author author <- params$author
@@ -49,13 +50,14 @@ coeffs <- params$coeffs
all_l <- REP$all_l all_l <- REP$all_l
ANOVAXLS <- REP$ANOVAXLS ANOVAXLS <- REP$ANOVAXLS
XLplot4pl <- REP$XLplot4pl
DiagnTable <- REP$DiagnTable DiagnTable <- REP$DiagnTable
UnRPLAausw <- REP$UnRPLAausw UnRPLAausw <- REP$UnRPLAausw
UnRPLBend <- REP$UnRPLBend UnRPLBend <- REP$UnRPLBend
PLAausw <- REP$PLAausw PLAausw <- REP$PLAausw
PLBend <- REP$PLBend PLBend <- REP$PLBend
LogPLAausw <- REP$LogPLAausw pottab4plXL <- REP$pottab4plXL
LogUnrPLAausw <- REP$LogUnrPLAausw Lim <- REP$Lim
XLdat2 <- REP$XLdat2 XLdat2 <- REP$XLdat2
@@ -70,262 +72,144 @@ relpotTestPlot <- REP$relpotTestPlot
# Introduction # Introduction
Bioassay potency estimation uses statistical methods to quantify the strength of a biological product or drug by comparing its response to that of a reference standard. Because biological responses are inherently variable, affected by assay conditions, cell systems or organisms, and measurement noise, the 4-parametric logistic regression is used to obtain reliable potency values. The variance for confidence interval calculation is coming from the regression procedure itself and is an excellent predictor for the variability of any future potency determinations. Bioassay potency estimation uses statistical methods to quantify the strength of a biological product or drug by comparing its response to that of a reference standard. Because biological responses are inherently variable, affected by assay conditions, cell systems or organisms, and measurement noise, the 4-parametric logistic regression is used to obtain reliable potency values.
USP<1034> recommends calculation of standard errors of ratios of the parameters using Fieller's theorem [Finney D.J. 1978] or using the "delta" method (for a discussion about the "delta" method see [Ver Hoef 2012]). However, the presented gradient approach using the differences on the log-scale is methematically more stable und thus preferable compared to any ratio approach ([Franz, V.H. 2007]). USP<1034> recommends calculation of standard errors of ratios of the parameters using Fieller's theorem [Finney DJ 1978] or using the "delta" method (for a discussion about the "delta" method see [Ver Hoef 2012]). However, the presented gradient approach using the differences on the log-scale is mathematically more stable und thus preferable compared to a ratio approach ([Franz VH 2007]).
# Results # Raw data
All data used for the 4PL evaluation is shown in table 1: All data used for the 4PL evaluation is shown in table 1:
```{r alll, echo=FALSE, warning=FALSE, results='asis'} ```{r alll, echo=FALSE, warning=FALSE, results='asis'}
kable(all_l, format = "markdown", caption= "Uploaded data (test and reference) in long format", digits=3) kable(XLdat2, format = "markdown", caption= "Uploaded data (test and reference) ", digits=3)
``` ```
The following 4 plots show all 4 models: restricted and unrestricted, and log transformed, respectively.
You can also embed plots, for example: # Results
```{r XLplot, echo=FALSE, warning=FALSE, fig.height=4, fig.width=6, fig.cap="Plot of models", fig.align='left'} ## Overall result
plot_f <- function(dat, sigmoid,det_sig) { ```{r Over_all, echo=FALSE, comment=NA, warning=NA, message=NA}
CORdat <- cor(dat[,1],dat[,ncol(dat)])
all_l <- melt(data.frame(dat), id.vars="log_dose", variable.name="replname", value.name = "readout") potFlag <- 0
isRef <- rep(c(1,0),1,each=nrow(all_l)/2) if (pottab4plXL["test_result"][[1]][1]==1) potFlag <- 1
isSample <- rep(c(0,1),1,each=nrow(all_l)/2) AnalysisFlag <- FALSE
all_l2 <- cbind(all_l, isRef, isSample) if (potFlag==1 | sum(testsTab$test_results)>0) AnalysisFlag <- TRUE
if(is.null(det_sig)) { colFmt <- function() {
if (CORdat<0) {
startlist <- list(a=sigmoid[3], b=-sigmoid[5],cs=sigmoid[7], outputFormat <- knitr::opts_knit$get("rmarkdown.pandoc.to")
d=sigmoid[1],r=sigmoid[8]) if(AnalysisFlag) {
} else { text <- paste("\\textcolor{red}{Analysis failed}",sep="")
startlist <- list(a=sigmoid[3],b=sigmoid[5],cs=sigmoid[7], } else {
d=sigmoid[1],r=sigmoid[8]) text <- paste("\\textcolor{black}{Analysis succeeded}>",sep="")
}
} 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])
} }
#browser() return(text)
tryCatch({
mr <- gsl_nls(fn = readout ~ 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))
},
error = function(err) {
err$message
})
s_mr <- summary(mr)
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[4] + (sigmoid[2] -sigmoid[4])/(1+exp(sigmoid[6]*(seq_x-(sigmoid[7]-sigmoid[8]))))
REFtrue <- sigmoid[3] + (sigmoid[1] -sigmoid[3])/(1+exp(sigmoid[5]*(seq_x-(sigmoid[7]))))
} else {
SAMPLEtrue <- det_sig[4] + (det_sig[6] -det_sig[4])/(1+exp(-det_sig[2]*(seq_x-(det_sig[8]))))
REFtrue <- det_sig[3] + (det_sig[5] -det_sig[3])/(1+exp(-det_sig[1]*(seq_x-(det_sig[7]))))
}
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))
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")
# 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))
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)
} }
plot_f(XLdat2, sigmoid=NULL, det_sig=coeffs)
```
`r colFmt()`
## 4pl-regression
Relative potency (absolute and relative confidence limits) are shown in Table 2:
```{r Pot_tab4pl, echo=FALSE, comment=NA, warning=NA, message=NA}
#browser()
if (pottab4plXL["test_result"][[1]][1]==1) { cat(paste("FAILED: relative potency CL result of restricted model outside limits: ", Lim[[9]], "and" ,Lim[[10]] ))}
kable(pottab4plXL, format = "markdown", caption= "Relative potency with absolute and relative CLs ", digits=3, row.names = F) %>%
kable_styling(latex_options = "hold_position")
``` ```
## Plot of the data and models
The following plots show the restricted and unrestricted model, respectively.
```{r XLplot, echo=FALSE, warning=FALSE, fig.height=4, fig.width=6, fig.cap="Plot of models", fig.align='left', comment=F, message=F, results='asis', fig.pos='H'}
library(cowplot)
plot_grid(XLplot4pl)
```
## ANOVA table
The ANOVA of the unconstrained model is listed in table 2: The ANOVA of the unconstrained model is listed in table 2:
```{r anovaxls, echo=FALSE, warning=FALSE, results='asis'} ```{r anovaxls, echo=FALSE, warning=FALSE, results='asis'}
kable(ANOVAXLS, format = "markdown", caption= "ANOVA table unrestricted", digits=3) kable(ANOVAXLS, format = "markdown", caption= "ANOVA table unrestricted", digits=3) %>%
kable_styling(latex_options = "hold_position")
```
## Analysis suitability tests
The following table lists the chosen suitabilit test results with confidence limits, where applicable:
```{r SST_ergebn, echo=FALSE, cache=FALSE, warning=FALSE, message=FALSE, tidy=TRUE}
kable(testsTab, row.names = F, format = "markdown", caption="SST results")
``` ```
```{r SST_ergebn, fig.align='center', fig.pos='htb!', echo=FALSE, cache=FALSE, warning=FALSE, message=FALSE, tidy=TRUE} \footnotesize
kable(testsTab[1:7,], row.names = F, format = "markdown", caption="SST results") ```{r Fussnote, echo=F, comment=NA}
cat("*...The estimate for F-test on regression and on non-linearity is the p-value")
cat( "F-test on regression passes if F-value > F-crit and thus p < 0.05")
cat( "F-test on non-linearity passes if F-value < F-crit and thus p > 0.05")
cat( "Test results outcome:")
cat(" 0 ... test passed (for EQ tests: CL within limits);")
cat(" 1 ... test failed (for EQ tests: CL not within limits);")
```
\normalsize
```{r AST_Ergebn, echo=FALSE, cache=FALSE, warning=FALSE, message=FALSE, tidy=TRUE}
TestsTabFlag <- FALSE
if (sum(testsTab$test_results)>0) TestsTabFlag <- TRUE
colFmt2 <- function() {
outputFormat <- knitr::opts_knit$get("rmarkdown.pandoc.to")
if(TestsTabFlag) {
text <- paste("\\textcolor{red}{Analysis suitability tests failed}",sep="")
} else {
text <- paste("\\textcolor{black}{Analysis suitability tests succeeded}",sep="")
}
return(text)
}
``` ```
*...The estimate for F-test on regression and on non-linearity is the p-value `r colFmt2()`
F-test on regression passes if F-value > F-crit and thus p < 0.05
F-test on non-linearity passes if F-value < F-crit and thus p > 0.05
Test results outcome:
0 ... test passed (for EQ tests: CI within limits);
1 ... test failed (for EQ tests CI not within limits);
-1 ... calculations unbound/denominator too close to 0
<!-- ```{r, label= 'CIplot', echo=FALSE, warning=FALSE, fig.width=100, fig.cap='Selected SSt confidence intervals with entered limits', fig.align='center'} -->
<!-- png("CIplot.png") -->
<!-- print(CIplot) -->
<!-- dev.off() -->
<!-- ``` -->
<!-- ![](CIplot.png){width=60%} -->
## Fitting results of the 4 models with bend points ## Fitting results of the 4 models with bend points
@@ -377,23 +261,6 @@ kable(UnRPLBend, format = "markdown", caption= "Bend points of 4PL unrestricted"
```{r LogPLAausw, echo=FALSE, warning=FALSE, results='asis'}
kable(LogPLAausw, format = "markdown", caption= "Restricted 4PL evaluation with log-transformed response", digits=3)
```
```{r LogUnRPLAausw, echo=FALSE, warning=FALSE, results='asis'}
kable(LogUnrPLAausw, format = "markdown", caption= "Unrestricted 4PL evaluation with log-transformed response", digits=3)
```
# Appendix: Formulas # Appendix: Formulas
## 4PL regression ## 4PL regression
+8
View File
@@ -49,3 +49,11 @@ expect_equal(TestTab, SolTab, check.attributes = FALSE)
``` ```
Note that the `echo = FALSE` parameter was added to the code chunk to prevent printing of the R code that generated the plot. Note that the `echo = FALSE` parameter was added to the code chunk to prevent printing of the R code that generated the plot.
``` {r}
y<-1:4;mean(y)
#> [1] 2.5
```
+33 -26
View File
@@ -776,7 +776,7 @@ server <- function(input, output, session) {
REP$PLAausw <- PLAAusw REP$PLAausw <- PLAAusw
REP$PLBend <- PLAAusw2 REP$PLBend <- PLAAusw2
#### Koeffizienten-Extraktion ---- #### Parameter extraktion ----
logcoeffs_R <- SRlog$coefficients[, 1] # logpot$coefficients logcoeffs_R <- SRlog$coefficients[, 1] # logpot$coefficients
names(logcoeffs_R) <- c("lower A", "Hill's slope", "upper A", "EC50 REF","EC50 DIFF") names(logcoeffs_R) <- c("lower A", "Hill's slope", "upper A", "EC50 REF","EC50 DIFF")
@@ -818,9 +818,13 @@ server <- function(input, output, session) {
if (exists("Ind")) { if (exists("Ind")) {
Dat$dilution <- XLdat[,Ind] Dat$dilution <- XLdat[,Ind]
} else Dat$dilution <- exp(XLdat[,logI]) } else Dat$dilution <- exp(XLdat[,logI])
# --- Plot-Ausgabe ---
##### Plot XL 4PL ----
output$XLplot <- renderPlot({ output$XLplot <- renderPlot({
plot_f(XLdat2, TransFlag=F) XLplot4pl <- plot_f(XLdat2, TransFlag=F)
REP$XLplot4pl <- XLplot4pl
XLplot4pl
}) })
REP$XLdat2 <- XLdat2 REP$XLdat2 <- XLdat2
@@ -1140,9 +1144,10 @@ server <- function(input, output, session) {
tab <- tests_FUNC(Dat$EXCEL, Limite, PureErrFlag = PureErrFlag) tab <- tests_FUNC(Dat$EXCEL, Limite, PureErrFlag = PureErrFlag)
tab[1,6:7] <- c("-","-") tab[1,6:7] <- c("-","-")
Dat$tests_FUNC <- tab
REP$testsTab <- tab
tab2 <- tab[SelTests,] tab2 <- tab[SelTests,]
Dat$tests_FUNC <- tab2
REP$testsTab <- tab2
dat <- datatable(tab2,options = list( dat <- datatable(tab2,options = list(
paging=TRUE, paging=TRUE,
@@ -1158,26 +1163,26 @@ server <- function(input, output, session) {
}) # observe }) # observe
#### plot CIs XL---- #### plot CIs XL----
observe({ # observe({
tab <- Dat$tests_FUNC # tab <- Dat$tests_FUNC
if (is.null(tab)) return(NULL) # if (is.null(tab)) return(NULL)
#
tab2 <- tab[-c(1,2,3,6),] # tab2 <- tab[-c(1,2,3,6),]
tab2[,3:7] <- apply(tab2[,3:7],2,as.numeric) # tab2[,3:ncol(tab2)] <- apply(tab2[,3:ncol(tab2)],2,as.numeric)
tab2[4:5,3:7] <- tab2[4:5,3:7]/100 # tab2[4:5,3:7] <- tab2[4:5,3:7]/100
#
p_CIs <- ggplot(tab2,aes(x=test,y=estimate, color=test,group=test)) + # p_CIs <- ggplot(tab2,aes(x=test,y=estimate, color=test,group=test)) +
geom_point() + # geom_point() +
geom_errorbar(aes(ymin=lower_CI, ymax=upper_CI), width=0.4) + # geom_errorbar(aes(ymin=lower_CI, ymax=upper_CI), width=0.4) +
geom_crossbar(aes(ymin=lower_limit, ymax=upper_limit), size=0.8) + # geom_crossbar(aes(ymin=lower_limit, ymax=upper_limit), size=0.8) +
coord_flip() + # coord_flip() +
theme_bw() + # theme_bw() +
theme(legend.position = "none",text = element_text(size=20)) # theme(legend.position = "none",text = element_text(size=20))
#
output$CIplot <- renderPlot({ p_CIs}, height=200) # output$CIplot <- renderPlot({ p_CIs}, height=200)
#
REP$CIplot <- p_CIs # REP$CIplot <- p_CIs
}) # })
#### simulated data tab Meta ---- #### simulated data tab Meta ----
output$simdat <- DT::renderDataTable({ output$simdat <- DT::renderDataTable({
@@ -1643,7 +1648,7 @@ server <- function(input, output, session) {
as.numeric(input$lEACratioSlope), as.numeric(input$uEACratioSlope), as.numeric(input$lEACratioSlope), as.numeric(input$uEACratioSlope),
as.numeric(input$lEACratioua), as.numeric(input$uEACratioua), as.numeric(input$lEACratioua), as.numeric(input$uEACratioua),
as.numeric(input$lowerPot), as.numeric(input$upperPot)) as.numeric(input$lowerPot), as.numeric(input$upperPot))
REP$Lim <- Lim
pottab4_ <- data.frame(pottab4) pottab4_ <- data.frame(pottab4)
pottab4_$potency <- round(as.numeric(pottab4[,2])*100,1) pottab4_$potency <- round(as.numeric(pottab4[,2])*100,1)
@@ -1667,6 +1672,8 @@ server <- function(input, output, session) {
pottab4_ <- cbind(pottab4_[,-(2:4)], data.frame(tests=c(test_potCI, test_potUCI,test_potCI_t,test_potUCI_t))) pottab4_ <- cbind(pottab4_[,-(2:4)], data.frame(tests=c(test_potCI, test_potUCI,test_potCI_t,test_potUCI_t)))
colnames(pottab4_) <- c("model","potency","lower95%CI","upper95%CI","relative_lower95%CI","relative_upper95%CI","test_result") colnames(pottab4_) <- c("model","potency","lower95%CI","upper95%CI","relative_lower95%CI","relative_upper95%CI","test_result")
REP$pottab4plXL <- pottab4_[1:2,]
output$pottab4plXL <- DT::renderDataTable({ output$pottab4plXL <- DT::renderDataTable({
dat <- datatable(pottab4_[1:2,], dat <- datatable(pottab4_[1:2,],
options=list(digits=3, options=list(digits=3,