ROUT report cleaned up, outlier table and Q-valus added
This commit is contained in:
Guxa committed
1 parent
0eba5831b6
commit
e263ba66e4
3 files changed
+105
-78
No files matched your search
+2
-2
@@ -1382,8 +1382,8 @@ tests_FUNC <- function(ro_new, Lim, PureErrFlag) {
|
||||
"EQ test ratio of upper asymptotes",
|
||||
"F-test on Lack-of-Fit*",
|
||||
"EQ test ratio of asymptote difference",
|
||||
"geom. rel. CI restr. model",
|
||||
"geom. rel. CI unrestr. model"
|
||||
"Potency CL restr. model",
|
||||
"Potency CL unrestr. model (fio)"
|
||||
),
|
||||
test_results = c(
|
||||
ifelse(p_F_regr < 0.05, 0, 1), test_la_diff, test_lowA, test_b, test_uppA,
|
||||
|
||||
+49
-31
@@ -23,7 +23,7 @@ params:
|
||||
author: "Author: `r params$author`"
|
||||
title: |
|
||||
| {width=1in}
|
||||
| 4PL bioassay evaluation
|
||||
| Robust outlier detection report
|
||||
subtitle: |
|
||||
`r params$FileName`
|
||||
|
||||
@@ -69,6 +69,8 @@ coeffsMU_ROUT <- REP$UnRPLAauswROUT
|
||||
CIplot <- REP$CIplot
|
||||
testsTabROUT <- REP$testsTabROUT
|
||||
relpotTestPlot <- REP$relpotTestPlot
|
||||
Qslider <- REP$Qslider
|
||||
OUTs_ <- REP$OUTs_
|
||||
|
||||
#browser()
|
||||
|
||||
@@ -87,17 +89,12 @@ All data used for the 4PL evaluation is shown in table 1:
|
||||
```{r alll, echo=FALSE, warning=FALSE, results='asis'}
|
||||
|
||||
kable(XLdat2, format = "markdown", caption= "Uploaded data (test and reference) ", digits=3)
|
||||
|
||||
if (!is.null(ro_newROUT)) {
|
||||
kable(ro_newROUT, format = "markdown", caption= "Data after exclusion of suspected outliers (see [6]) ", digits=3)
|
||||
}
|
||||
|
||||
```
|
||||
|
||||
|
||||
# Results
|
||||
|
||||
## Overall result
|
||||
## Results before exclusion
|
||||
|
||||
```{r Over_all, echo=FALSE, comment=NA, warning=NA, message=NA}
|
||||
|
||||
@@ -124,9 +121,9 @@ colFmt <- function() {
|
||||
`r colFmt()`
|
||||
|
||||
|
||||
## 4pl-regression
|
||||
### 4pl-regression
|
||||
|
||||
Relative potency (absolute and relative confidence limits) BEFORE EXCLUSION of potential outliers are shown in Table 3 for information. `r if(PureErr) {"Pure Error is used for calculations."}`
|
||||
Relative potency (absolute and relative confidence limits) BEFORE EXCLUSION of potential outliers are shown in Table 4 for information. `r if(PureErr) {"Pure Error is used for calculations."}`
|
||||
`r if (!PureErr) {"RMSE of restricted model is used for confidence limit calculation."}`
|
||||
|
||||
```{r Pot_tab4pl, echo=FALSE, comment=NA, warning=NA, message=NA}
|
||||
@@ -139,12 +136,12 @@ kable(pottab4plXL, format = "markdown", caption= "Relative potency with absolute
|
||||
|
||||
|
||||
```
|
||||
NOTE: results of unrestricted model for Information only.
|
||||
NOTE: results of unrestricted model for information only.
|
||||
|
||||
|
||||
## Plot of the data and models
|
||||
## Robust non-linear fitting and exclusion plot
|
||||
|
||||
Plots in Figure 1 shows the restricted model.
|
||||
Figure 1 shows the data with suspected outliers in red (reference in circles, test sample in triangles). The blue and red lines mark the 4PL fit AFTER exclusion of suspected outliers.
|
||||
|
||||
|
||||
```{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'}
|
||||
@@ -155,11 +152,33 @@ plot(ROUTplot)
|
||||
```
|
||||
|
||||
|
||||
# Results after exclusion of potential outliers
|
||||
Q-value set for exclusion of results (see [6]): `r Qslider`
|
||||
|
||||
## ANOVA table
|
||||
```{r ronewROUT, echo=FALSE, warning=FALSE, results='asis'}
|
||||
|
||||
The ANOVA of the unconstrained model is listed in table 4. Bates and Watts [4] proposed a test on parallelism which compares the residual sum of squares of the restricted model (ResRSSE) with the residual sum of squares of the unrestricted model (UnresRSSE). If the UnresRSSE is significantly smaller than the ResRSSE, the p-value of "Non-parallelism" is smaller than 0.05 (line 4 in table 4). This test is for information only as it may be overly sensitive in case of small overall variability of the data.
|
||||
if (!is.null(ro_newROUT)) {
|
||||
kable(ro_newROUT, format = "markdown", caption= "Data after exclusion of suspected outliers ", digits=3)
|
||||
}
|
||||
|
||||
```
|
||||
|
||||
|
||||
```{r label='OUTs_', echo=FALSE, warning=FALSE, results='asis'}
|
||||
|
||||
if (!is.null(OUTs_)) {
|
||||
kable(OUTs_, format = "markdown", caption= "Excluded results", digits=3)
|
||||
}
|
||||
|
||||
```
|
||||
|
||||
|
||||
|
||||
|
||||
## Results after exclusion of potential outliers
|
||||
|
||||
### ANOVA table
|
||||
|
||||
The ANOVA of the unconstrained model is listed in table 5. Bates and Watts [4] proposed a test on parallelism which compares the residual sum of squares of the restricted model (ResRSSE) with the residual sum of squares of the unrestricted model (UnresRSSE). If the UnresRSSE is significantly smaller than the ResRSSE, the p-value of "Non-parallelism" is smaller than 0.05 (line 4 in table 4). This test is for information only as it may be overly sensitive in case of small overall variability of the data.
|
||||
|
||||
```{r anovaxls, echo=FALSE, warning=FALSE, results='asis'}
|
||||
|
||||
@@ -169,9 +188,9 @@ kable(ANOVA_ROUT, format = "markdown", caption= "Analysis of variance", digits=3
|
||||
|
||||
```
|
||||
|
||||
## Assay suitability tests
|
||||
### Assay suitability tests
|
||||
|
||||
Table 5 lists the chosen suitability test results with confidence limits, where applicable. F-tests should be read with caution, if the overall variability is small, as the test gets overly sensitive.
|
||||
Table 6 lists the chosen suitability test results with confidence limits, where applicable. F-tests should be read with caution. If the overall variability is small, the tests get overly sensitive.
|
||||
|
||||
|
||||
```{r SST_ergebn, echo=FALSE, cache=FALSE, warning=FALSE, message=FALSE, tidy=TRUE}
|
||||
@@ -192,7 +211,7 @@ 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);")
|
||||
|
||||
cat(" Potency CL for information only")
|
||||
|
||||
|
||||
```
|
||||
@@ -224,9 +243,9 @@ colFmt2 <- function() {
|
||||
|
||||
|
||||
|
||||
## Fitting results with curve points
|
||||
### 4PL results with curve points
|
||||
|
||||
The results of the non-linear fitting procedure for the restricted model (5 parameters) is listed in Table 6:
|
||||
The results of the non-linear fitting procedure for the restricted model (5 parameters) is listed in Table 7:
|
||||
|
||||
```{r PLAausw, echo=FALSE, warning=FALSE, results='asis'}
|
||||
|
||||
@@ -237,7 +256,7 @@ kable(coeffsMR_ROUT , format = "markdown", caption= "Restricted 4PL model", digi
|
||||
|
||||
|
||||
|
||||
The results of the non-linear fitting procedure for the unrestricted model (8 parameters) is listed in table 7:
|
||||
The results of the non-linear fitting procedure for the unrestricted model (8 parameters) is listed in table 8:
|
||||
|
||||
```{r UnRPLAausw, echo=FALSE, warning=FALSE, results='asis'}
|
||||
|
||||
@@ -246,9 +265,8 @@ kable(coeffsMU_ROUT , format = "markdown", caption= "Unrestricted 4PL model", di
|
||||
```
|
||||
|
||||
|
||||
## Bend points and asymptote points
|
||||
|
||||
Sebaugh et al. proposed bend points for test and reference samples, that define the points with highest turning behavior. Table 6 lists these bendpoints as well as asymptote points ~ twice as far from the center as the bendpoints.
|
||||
Sebaugh et al. proposed bend points for test and reference samples, that define the points with highest turning behavior. Table 9 lists these bend-points as well as asymptote-points (i.e. twice as far from the center as the bendpoints). Asymptote points may serve as indicators for dilution setting: only one dilution is needed outside the asymptote point, if the Hill's slope is between 0.5 and 2.
|
||||
|
||||
```{r PLBend, echo=FALSE, warning=FALSE, results='asis'}
|
||||
|
||||
@@ -267,7 +285,7 @@ kable(BendsROUT2, format = "markdown", caption= "Bendpoints and asymptote points
|
||||
\noindent
|
||||
\begin{tabular}{p{6cm}p{1cm}p{6cm}}
|
||||
\cline{1-1} \cline{3-3}
|
||||
Date & & Signature
|
||||
Date & & Signature (`r author`)
|
||||
\end{tabular}
|
||||
|
||||
|
||||
@@ -312,20 +330,20 @@ $$
|
||||
B_{true}=B*\frac{D-A}{4}
|
||||
$$
|
||||
|
||||
## Confidence intervals
|
||||
## Confidence limits
|
||||
|
||||
In general, the confidence intervals are calculated as follows:
|
||||
In general, the confidence limits are calculated as follows:
|
||||
$$
|
||||
CI = \hat\theta\pm se(\hat\theta)*q^{t_{n-p}}_{1-\frac{\alpha}{2}}
|
||||
CL = \hat\theta\pm se(\hat\theta)*q^{t_{n-p}}_{1-\frac{\alpha}{2}}
|
||||
$$
|
||||
…where $\hat\theta$ is a fitted parameter or a linear combination thereof, q is the 1-alpha/2 quantile of the Student’s t-distribution with n-p degrees of freedom and se is the standard error derived from any covariance matrix.
|
||||
|
||||
Let $\theta$ be the 4+1 parameters of the fit (a, b, d, EC50 of reference and EC50 difference). It can be shown that the least squares estimator $\hat\theta$ is normally distributed with asymptotic covariance matrix. The gradient method provides one of several ways to calculate the covariance matrix:
|
||||
Let $\theta$ be the 4+1 parameters of the fit (A,B,D, EC50 of reference and EC50 difference). It can be shown that the least squares estimator $\hat\theta$ is normally distributed with asymptotic covariance matrix. The gradient method provides one of several ways to calculate the covariance matrix:
|
||||
|
||||
$$
|
||||
\hat{V(\theta)}= \sigma^2(A(\hat\theta)^T*A(\hat\theta))^{-1}
|
||||
\widehat{V(\theta)}= \sigma^2(A(\hat\theta)^T*A(\hat\theta))^{-1}
|
||||
$$
|
||||
where A($\theta$) is the n x p matrix of the first partial derivatives for each parameter (i.e. gradient) realized at the fitted parameter estimates. The RMSE of the model or the pure error is used as estimate of $\sigma$. The square root of the diagonals of $\hat{V(\theta)}$ gives the standard errors and with that confidence intervals (CI) can be computed.
|
||||
where A($\theta$) is the n x p matrix of the first partial derivatives for each parameter (i.e. gradient) realized at the fitted parameter estimates. The RMSE of the model or the pure error is used as estimate of $\sigma$. The square root of the diagonals of $\widehat{V(\theta)}$ gives the standard errors and with that confidence intervals (CI) can be computed.
|
||||
|
||||
# Literature
|
||||
|
||||
@@ -339,4 +357,4 @@ where A($\theta$) is the n x p matrix of the first partial derivatives for each
|
||||
|
||||
[5] Bates, D.M., Watts, D.G. (1988) 2. In: Nonlinear Regression Analysis and Its Applications. New York: Wiley, pp 52-58
|
||||
|
||||
[6] Motulsky, Brown Outlier testing
|
||||
[6] Motulsky HJ, Brown RE. Detecting outliers when fitting data with nonlinear regression - a new method based on robust nonlinear regression and the false discovery rate. BMC Bioinformatics. 2006 Mar 9;7:123. doi: 10.1186/1471-2105-7-123. PMID: 16526949; PMCID: PMC1472692.
|
||||
@@ -338,7 +338,8 @@ server <- function(input, output, session) {
|
||||
"Adjust the slider to mark the suspected outliers. Mostly, a Q-value of 2% is sufficient.",
|
||||
"If the general variability of the data is high, many datapoints will be flagged, also ones that are not deemed to be outliers",
|
||||
"An indicator for high variability is, when the Q-value needs to be increased above 5%, to flag suspected outliers.",
|
||||
"Then, please re-consider, if the suspected outlier is not 'just' normal variability."
|
||||
"Then, please re-consider, if the suspected outlier is not 'just' normal variability.",
|
||||
tableOutput("ANOVA_ROUT")
|
||||
),
|
||||
tabPanel(
|
||||
"parameter estimates",
|
||||
@@ -576,13 +577,13 @@ server <- function(input, output, session) {
|
||||
" CI for difference of slopes:",
|
||||
tableOutput("SlopeDiffCIMeta"),
|
||||
)
|
||||
),
|
||||
tabPanel(
|
||||
"Report",
|
||||
h4("Settings for report"),
|
||||
downloadButton("downloadXLReportMeta", label = "Download PDF report", class = "butt"),
|
||||
tags$style(type = "text/css", "#downloadXLReportMeta {background-color: orange; color: black;font-family: COurier New}"),
|
||||
)
|
||||
) #,
|
||||
# tabPanel(
|
||||
# "Report",
|
||||
# h4("Settings for report"),
|
||||
# downloadButton("downloadXLReportMeta", label = "Download PDF report", class = "butt"),
|
||||
# tags$style(type = "text/css", "#downloadXLReportMeta {background-color: orange; color: black;font-family: COurier New}"),
|
||||
# )
|
||||
)
|
||||
)
|
||||
# )
|
||||
@@ -822,6 +823,9 @@ server <- function(input, output, session) {
|
||||
##### ROUT outlier testing ----
|
||||
#browser()
|
||||
if(!is.null(input$Qslider)) {
|
||||
|
||||
|
||||
|
||||
all_lROUT <- all_l[complete.cases(all_l),]
|
||||
colnames(all_lROUT) <- c("log_dose","sample","y","isRef","isSample","conc")
|
||||
res <- rout_4pl_potency(all_lROUT, Q=input$Qslider/100)
|
||||
@@ -870,6 +874,9 @@ server <- function(input, output, session) {
|
||||
REP$ROUTplot <- ROUTplot
|
||||
|
||||
ANOVA_ROUT <- ANOVA4plUnresfunc(ro_new = ro_newROUT)
|
||||
output$ANOVA_ROUT <- renderTable({
|
||||
ANOVA_ROUT
|
||||
})
|
||||
REP$ANOVA_ROUT <- ANOVA_ROUT
|
||||
|
||||
Limite <- list(
|
||||
@@ -915,6 +922,8 @@ server <- function(input, output, session) {
|
||||
|
||||
|
||||
REP$BendsROUT <- BendsROUT
|
||||
REP$Qslider <- input$Qslider
|
||||
REP$OUTs_ <- OUTs_
|
||||
}
|
||||
|
||||
|
||||
@@ -2959,43 +2968,43 @@ server <- function(input, output, session) {
|
||||
|
||||
#### download Meta 4PL report----
|
||||
|
||||
observeEvent(input$btn2, {
|
||||
if(!Dat$FITsFlag) {
|
||||
runjs("$('#downloadXLReportMeta')[0].click();")
|
||||
} else {
|
||||
showModal(modalDialog(
|
||||
title = "No 4PL model to Download",
|
||||
"Please select other data before download.",
|
||||
easyClose = TRUE,
|
||||
footer = NULL
|
||||
))
|
||||
}
|
||||
})
|
||||
|
||||
output$downloadXLReportMeta <- downloadHandler(
|
||||
filename = paste0("Report_4PLEvaluation", Dat$RepIdentifier, ".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, "logov2.png")
|
||||
file.copy("logov2.png", tempReportc, overwrite = T)
|
||||
|
||||
rmarkdown::render(tempReport,
|
||||
output_file = file,
|
||||
params = list(
|
||||
FileName = Dat$FileName,
|
||||
author = Dat$Author,
|
||||
NoP = Dat$NoP,
|
||||
Assay = Dat$Assay,
|
||||
REP = REP,
|
||||
coeffs = Dat$coeffs_UN
|
||||
),
|
||||
envir = new.env(parent = globalenv())
|
||||
)
|
||||
}
|
||||
)
|
||||
# observeEvent(input$btn2, {
|
||||
# if(!Dat$FITsFlag) {
|
||||
# runjs("$('#downloadXLReportMeta')[0].click();")
|
||||
# } else {
|
||||
# showModal(modalDialog(
|
||||
# title = "No 4PL model to Download",
|
||||
# "Please select other data before download.",
|
||||
# easyClose = TRUE,
|
||||
# footer = NULL
|
||||
# ))
|
||||
# }
|
||||
# })
|
||||
#
|
||||
# output$downloadXLReportMeta <- downloadHandler(
|
||||
# filename = paste0("Report_4PLEvaluation", Dat$RepIdentifier, ".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, "logov2.png")
|
||||
# file.copy("logov2.png", tempReportc, overwrite = T)
|
||||
#
|
||||
# rmarkdown::render(tempReport,
|
||||
# output_file = file,
|
||||
# params = list(
|
||||
# FileName = Dat$FileName,
|
||||
# author = Dat$Author,
|
||||
# NoP = Dat$NoP,
|
||||
# Assay = Dat$Assay,
|
||||
# REP = REP,
|
||||
# coeffs = Dat$coeffs_UN
|
||||
# ),
|
||||
# envir = new.env(parent = globalenv())
|
||||
# )
|
||||
# }
|
||||
# )
|
||||
|
||||
|
||||
#### download Wizard report ----
|
||||
|
||||
Reference in new issue
Block a user