Nothing
R/validation_error_srsc.R
In BayesianFROC: FROC Analysis by Bayesian Approaches
Defines functions
validation.draw_srsc
error_srsc_variance_visualization
error_srsc_error_visualization
error_srsc
validation.dataset_srsc
Documented in
error_srsc
error_srsc_error_visualization
error_srsc_variance_visualization
validation.dataset_srsc
validation.draw_srsc
#' @title Errors of Estimator for any Given true parameter
#@title ----
#'@description
#'
#' This function replicates many datasets from a model with a given model parameter as truth.
#' Then fit a model to these datasets and get estimates.
#' Then calculates mean or variance of the difference of estimates and truth.
#'
#'
#' This function gives a validation under the assumption that
#' we obtain a fixed true model parameter drawn from the prior.
#'
#'But, you know, the author noticed that it is not sufficient because
#' , in Bayesian, such a true parameter is merely one time sample from
#' prior. So, what we should do is to draw randomly many samples from priors.
#'
#' Well, I know, but, I am being such a couch potato. In the future, I would do,
#' if I was alive.
#' Even if I try this, this effort never take account by someone.
#' But, to live, it is required. No money, no life. 2020 NOv 29
#'
#'
#'
#' In the future, what I should do is that the following calculations.
#'
#'
#' 1. Draw a sample of parameter \eqn{\theta} from prior, namely, \eqn{\theta \sim \pi(\theta)}
#'
#' 2. Draw a sample of data \eqn{x=x(\theta)} from a model with the sample of prior \eqn{ x \sim \pi(x|\theta)}
#'
#' 3. Draw a sample of parameter \eqn{\theta'=\theta'(x(_\theta))} from a posterior with the sample of data \eqn{ \theta' \sim \pi(\theta|x)}
#'
#' 4. Calculates the integral \eqn{ \int | \theta' - \theta|^2 \pi(\theta)d\theta} as the error of estimates among all priors..
#'
#'
#'
#'
#'
#'
#'
#' @details
#'
#'Let us denote a model parameter by \eqn{\theta_0}
#' \eqn{N_I} by a number of
#' images and number
#' of lesions by \eqn{N_L} which are specified
#' by user as the variables of the function.
#'
#' \describe{
#' \item{ \strong{(I)} Replicates models for datasets \eqn{D_1,D_2,...,D_k,...,D_K}. }{ }
#' \item{ Draw a dataset \eqn{D_k} from a likelihood (model), namely }{ \eqn{D_k \sim likelihood(|\theta_0)}. }
#' \item{ Draw a MCMC samples \eqn{\{ \theta_i (D_k)\}} from a posterior, namely }{ \eqn{ \theta _i \sim \pi(|D_k)}. }
#' \item{ Calculate a posterior mean, namely }{ \eqn{ \bar{\theta}(D_k) := \sum_i \theta_i(D_k) }. }
#' \item{ Calculates error for \eqn{D_k} }{ \eqn{\epsilon_k}:=Trueth - posterior mean estimates of \eqn{D_k} = \eqn{|\theta_0 - \bar{\theta}(D_k)|} (or = \eqn{\theta_0 - \bar{\theta}(D_k)}, accordinly by the user specified \code{absolute.errors} ). }
#' \item{ \strong{(II)} Calculates mean of errors }{ mean of errors \eqn{ \bar{\epsilon}(\theta_0,N_I,N_L)}= \eqn{ \frac{1}{K} \sum \epsilon_k } }
#' }
#'
#' Running this function, we can see that the error \eqn{ \bar{\epsilon}(\theta_0,N_I,N_L)} decreases
#' monotonically as a given number of images \eqn{N_I} or a given number of lesions \eqn{N_L} increases.
#'
#' Also, the scale of error also will be found. Thus this function can show how our estimates are correct.
#' Scale of error differs for each componenet of model parameters.
#'
#'
#'
#' Revised 2019 August 28
#'
#'
#'
#'
#'
#'
#'
#'@return Return values is,
#'
#' \describe{
#' \item{ Stanfit objects }{ for each Replicated datasets }
#' \item{ Errors }{ EAPs minus true values, in the above notations, it is \eqn{ \bar{\epsilon}(\theta_0,N_I,N_L)} }
#' \item{ Variances of estimators. }{ This calculates the vaiance of posterior means over all replicated datasets }
#' }
#'
#'@param mean.truth This is a parameter
#'of the latent Gaussian assumption
#'for the noise distribution.
#'@param sd.truth This is a parameter of the latent
#' Gaussian assumption for the noise distribution.
#'@param z.truth This is a parameter of the
#' latent Gaussian assumption for the noise distribution.
# @param seed_sampling to be passed to rstan::sampling()
#'@param NL Number of Lesions.
#'@param NI Number of Images.
#'@param replicate.datset A Number indicate
#'that how many you replicate dataset
#' from user's specified dataset.
#'@param serial.number A positive integer
#'or Character. This is for programming perspective.
#' The author use this to print the serial
#' numbre of validation. This will be used
#' in the validation function.
#'@param base_size An numeric for size of object, this is for the package developer.
#'@param absolute.errors A logical specifying whether mean of errors is defined by
#'
#' \describe{
#' \item{ \code{TRUE} }{ \eqn{ \bar{\epsilon}(\theta_0,N_I,N_L)}= \eqn{ \frac{1}{K} \sum | \epsilon_k | } }
#' \item{ \code{FALSE} }{ \eqn{ \bar{\epsilon}(\theta_0,N_I,N_L)}= \eqn{ \frac{1}{K} \sum \epsilon_k } }
#' }
#'
#'
#'
#'@inheritParams fit_Bayesian_FROC
#'@inheritParams DrawCurves
#'@examples
#' \dontrun{
# ####1#### ####2#### ####3#### ####4#### ####5#### ####6#### ####7#### ####8#### ####9####
#'#=========================== The first example ======================================
#'
#'
#'# It is sufficient to run the function with default variable
#'
#' datasets <- validation.dataset_srsc()
#'
#'
#'# Today 2020 Nov 29 I have completely forgotten this function, oh it helps me. Thank me.
# ####1#### ####2#### ####3#### ####4#### ####5#### ####6#### ####7#### ####8#### ####9####
#'#============================= The second example ======================================
#'
#'# If user does not familiar with the values of thresholds, then
#'# it would be better to use the actual estimated values
#'# as an example of true parameters. In the following,
#'# I explain this.
#'# First, to get posterior mean estimates, we run the following:
#'
#'
#'
#' fit <- fit_Bayesian_FROC(dataList.Chakra.1,ite = 1111,summary =FALSE,cha=3)
#'
#'
#'
#'
#'
#'
#'# Secondly, extract the expected a posterior estimates (EAPs) from the object fit
#'# Note that EAP is also called "posterior mean"
#'
#' z <- rstan::get_posterior_mean(fit,par=c("z"))[,"mean-all chains"]
#'
#'
#'
#'
#'
#'# Thirdly we use this z as a true value.
#'
#'
#' datasets <- validation.dataset_srsc(z.truth = z)
#'
#'
#'
#'#========================================================================================
#'# 1) extract replicated fitted model object
#'#========================================================================================
#'
#'
#' # Replicates models
#'
#' a <- validation.dataset_srsc(replicate.datset = 3,ite = 111)
#'
#'
#'
#' # Check convergence, in the above MCMC iterations = 111 which is too small to get
#' # a convergence MCMC chain, and thus the following example will the example
#' # of a non-convergent model in the r hat criteria.
#'
#' ConfirmConvergence( a$fit[[3]])
#'
#'
#' # Check trace plot to confirm whether MCMC chain converge or not.
#'
#' stan_trace( a$fit[[3]],pars = "A")
#'
#'
#' # Check p value, for chi square goodness of fit whose null hypothesis is that
#' # the model is fitted well.
#'
#' fit@posterior_predictive_pvalue_for_chi_square_goodness_of_fit
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#' # Revised in 2019 August 29
#'
#' # Revised in 2020 Nov 28
#' # It is weird, awesome,
#' # What a fucking English,...I fix it.
#'
#'
#'
#'
#'#========================================================================================
#'# 1) Histogram of error of postrior means for replicated datasets
#'#========================================================================================
#'#'
#'
#' a<- validation.dataset_srsc(replicate.datset = 100)
#' hist(a$error.of.AUC,breaks = 111)
#' hist(a$error.of.AUC,breaks = 30)
#'
#'
#'
#'
#'
#'
# ####1#### ####2#### ####3#### ####4#### ####5#### ####6#### ####7#### ####8#### ####9####
#'#========================================================================================
#'# absolute.errors = FALSE generates negative biases
#'#========================================================================================
#'
#'
#' validation.dataset_srsc(absolute.errors = FALSE)
#'
#'
#'
# ####1#### ####2#### ####3#### ####4#### ####5#### ####6#### ####7#### ####8#### ####9####
#'#========================================================================================
#'# absolute.errors = TRUE coerce negative biases to positives, i.e., L^2 norm
#'#========================================================================================
#'
#'
#' validation.dataset_srsc(absolute.errors = TRUE)
#'
#'
#'
#'
#'
#'
#'
#'
# ####1#### ####2#### ####3#### ####4#### ####5#### ####6#### ####7#### ####8#### ####9####
#'#========================================================================================
#'# Check each fitted model object
#'#========================================================================================
#'
#'
#'
#'
#' a <- validation.dataset_srsc(verbose = TRUE)
#'
#' a$fit[[2]]
#'
#'
#'
#'class(a$fit[[2]])
#'
#' rstan::traceplot(a$fit[[2]], pars = c("A"))
#'
#'
#'
#'
#'
#'
# ####1#### ####2#### ####3#### ####4#### ####5#### ####6#### ####7#### ####8#### ####9####
#'#========================================================================================
#'# NaN ... why? 2021 Dec
#'#========================================================================================
#'
#' fits <- validation.dataset_srsc()
#'
#' f <-fits$fit[[1]]
#' rstan::extract(f)$dl
#' sum(rstan::extract(f)$dl)
#' Is.nan.in.MCMCsamples <- as.logical(!prod(!is.nan(rstan::extract(f)$dl)))
#' rstan::extract(f)$A[525]
#' a<-rstan::extract(f)$a[525]
#' b<-rstan::extract(f)$b[525]
#'
#' Phi( a/sqrt(b^2+1) )
#' x<-rstan::extract(f)$dl[2]
#'
#' a<-rstan::extract(f)$a
#' b<-rstan::extract(f)$b
#'
#' a/(b^2+1)
#' Phi(a/(b^2+1))
#' mean( Phi(a/(b^2+1)) )
#'
#' #'
#'
#'
#'
#'
#'}# dontrun
#' @export
validation.dataset_srsc <-function(
replicate.datset =3,
ModifiedPoisson = FALSE,
mean.truth=0.6,
sd.truth=5.3,
z.truth =c(-0.8,0.7,2.38),
NL=259,
NI=57,
ite =1111,
cha =1,
summary=TRUE,
see = 1234,
verbose = FALSE,
serial.number=1,
base_size=0,
absolute.errors = TRUE
# , convergence.model.only = FALSE
){
C <- length(z.truth)
c<-C:1
dz.truth <-vector()
for (cd in 1:C-1) {
dz.truth[cd] <- z.truth[cd+1]-z.truth[cd]
}
if(ModifiedPoisson==F){ NX <- NI}
if(ModifiedPoisson==T){ NX <- NL}
p.truth <- vector()
list.of.dataList <- list()
vector.of.dataList <- list()
name.vector.of.dataList <- vector()
for (cd in 1:C) {
name.vector.of.dataList[cd] <- paste("f[",C-cd+1,"]",sep = "")
}
for (cd in 1:C) {
name.vector.of.dataList[C+cd] <- paste("h[",C-cd+1,"]",sep = "")
}
name.vector.of.dataList[2*C+1] <- "NL"
name.vector.of.dataList[2*C+2] <- "NI"
l.truth <- -log( stats::pnorm(z.truth))
for (seed in 1:replicate.datset ) {
f.inv <- vector()#these should in for sentence
h.inv <- vector()
f <- vector()
h <- vector()
for (cd in 1:C) {
if (cd==C) {p.truth[cd]= 1 - stats::pnorm((z.truth[cd]-mean.truth)/sd.truth)
}else{
p.truth[cd]<- stats::pnorm((z.truth[cd+1]-mean.truth)/sd.truth) - stats::pnorm((z.truth[cd]-mean.truth)/sd.truth)
}}
# here is hitrate modified -----
deno <-vector()
hit_rate.truth <- vector()
deno[C-1] = 1-p.truth[C];
for(cd in 3:C){ deno[c[cd]] = deno[c[cd-1]]-p.truth[c[cd-1]]; }
hit_rate.truth[C] = p.truth[C];
for(cd in 1:C-1) hit_rate.truth[cd] = p.truth[cd]/deno[cd];
s <-0# 2019 Sept 30. This is a key idea
c <-C:1
for (cd in 1:C) {
if(ModifiedPoisson==F){
if (cd==C) {
set.seed(seed = seed);
f.inv[cd] <- stats::rpois(n= 1, lambda = (l.truth[cd]-0)*NI )
}else{
set.seed(seed = seed);
f.inv[cd] <- stats::rpois(n= 1, lambda = (l.truth[cd]-l.truth[cd+1])*NI )
}#else
}# if ModifiedPoisson==F
if(ModifiedPoisson==T){
if (cd==C) {
set.seed(seed = seed);
f.inv[cd] <- stats::rpois(n= 1, lambda = (l.truth[cd]-0)*NL )
}else{
set.seed(seed = seed);
f.inv[cd] <- stats::rpois(n= 1, lambda = (l.truth[cd]-l.truth[cd+1])*NL)
}#else set.seed(seed = seed); hits <- stats::rbinom(n=1,size = NL,prob = p.truth[cd])
}# if ModifiedPoisson==T
# set.seed(seed = seed); h.inv[cd] <- stats::rbinom(n=1,
# size = NL-s, # 2019 Sept 30. This is a key idea # In order to avoid the sum of hits may be greater than NL
# prob = p.truth[cd])
# s <- s + h.inv[cd] # 2019 Sept 30. This is a key idea # In order to avoid the sum of hits may be greater than NL
#
#
set.seed(seed = seed); h[cd] <- stats::rbinom(n=1,
size = NL-s, # 2019 Sept 30. This is a key idea # In order to avoid the sum of hits may be greater than NL
prob = hit_rate.truth[c[cd]])
s <- s + h[cd] # 2019 Sept 30. This is a key idea # In order to avoid the sum of hits may be greater than NL
f[C-cd+1] <- f.inv[cd]
# h[C-cd+1] <- h.inv[cd]
}# for cd in 1:C
list.of.dataList[[seed]] <- list(
NL=NL,
NI=NI,
f=f,
h=h,
C=C
)
a<-append( list.of.dataList[[seed]]$f, list.of.dataList[[seed]]$h);
b <-append( list.of.dataList[[seed]]$NL, list.of.dataList[[seed]]$NI);
vector.of.dataList[[seed]] <-append(a,b)
# browser()
names(vector.of.dataList[[seed]] ) <- name.vector.of.dataList
}#for seed
#Here
#
# * Datasets were created !!
#
# * Next Fitting !!
#
#Here
fit <- list()
validation.of.AUC.EAP <-vector()
validation.of.AUC.CI.lower <-vector()
validation.of.AUC.CI.upper <-vector()
validation.of.z.Threshold.EAP <- list()
validation.of.z.Threshold.CI.lower <- list()
validation.of.z.Threshold.CI.upper <- list()
validation.of.dz.Threshold.EAP <- list()
validation.of.dz.Threshold.CI.lower <- list()
validation.of.dz.Threshold.CI.upper <- list()
validation.of.mean.of.latent.EAP <- vector()
validation.of.mean.of.latent.CI.lower <- vector()
validation.of.mean.of.latent.CI.upper <- vector()
validation.of.Standard.Deviation.latent.EAP <- vector()
validation.of.Standard.Deviation.latent.CI.lower <- vector()
validation.of.Standard.Deviation.CI.upper <- vector()
validation.of.p.Hitrate.EAP <- list()
validation.of.p.Hitrate.CI.lower <- list()
validation.of.p.Hitrate.CI.upper <- list()
validation.of.l.FalseRate.EAP <- list()
validation.of.l.FalseRate.CI.lower <- list()
validation.of.l.FalseRate.CI.upper <- list()
validation.of.loglikelihood.of.model <- vector()
convergent.dataList <- list()
convergent.dataList.as.dataframe <- data.frame(row.names = name.vector.of.dataList)
# rownames(convergent.dataList.as.dataframe) <- name.vector.of.dataList
s<-0
# base_size <-0
for (seed in 1:replicate.datset ) {
#It seems better to make a indicator and
#do not print the result of rstan
# fit ------
color_message(seed,"-th fitting STARTS-----------------------------------------------------------------------")
fit[[seed]] <- fit_Bayesian_FROC(
dataList=list.of.dataList[[seed]],
ite = ite,
cha = cha,
see = see,#seed ---------
new.imaging.device = FALSE,
summary = FALSE,
verbose = verbose,
ModifiedPoisson = ModifiedPoisson,
DrawCurve = FALSE
)
# };browser(); for (seed in 1:replicate.datset ) {
print(stanfit_from_its_inherited_class( fit[[seed]] ) )
base_size <- base_size + size_of_return_value(fit[[seed]] ,is_return_value = FALSE,summary = FALSE)
# here col of size ----
print(size_of_return_value(fit[[seed]] ,is_return_value = FALSE, base_size =base_size,col=TRUE ))
# view data ----
print(viewdata(list.of.dataList[[seed]]))
message("\n----------------------------------------\n",sep = "")
message("\n* ",serial.number, " -th validation")
message("\n* The (",seed," / ", replicate.datset, " ) -th fitting finished.\n",sep = "")
message("\n----------------------------------------\n",sep = "")
# Here ------
if ( fit[[seed]]@convergence ==TRUE) {
s <-s+1
convergent.dataList[[s]] <-fit[[seed]]@dataList
a<-append( convergent.dataList[[s]]$f, convergent.dataList[[s]]$h);
b <-append( convergent.dataList[[s]]$NL, convergent.dataList[[s]]$NI);
convergent.dataList.as.dataframe[,s] <-append(a,b)
# browser()
#ssss -----
ssss <- summary_EAP_CI_srsc( fit[[seed]],summary = FALSE,dig =11)
validation.of.AUC.EAP[s] <- ssss$AUC.EAP
validation.of.AUC.CI.lower[s] <-ssss$AUC.CI.lower
validation.of.AUC.CI.upper[s] <- ssss$AUC.CI.upper
validation.of.z.Threshold.EAP[[s]] <- ssss$z.Threshold.EAP
validation.of.z.Threshold.CI.lower[[s]] <- ssss$z.Threshold.CI.lower
validation.of.z.Threshold.CI.upper[[s]] <- ssss$z.Threshold.CI.upper
validation.of.dz.Threshold.EAP[[s]] <- ssss$dz.Threshold.EAP
validation.of.dz.Threshold.CI.lower[[s]] <- ssss$dz.Threshold.CI.lower
validation.of.dz.Threshold.CI.upper[[s]] <- ssss$dz.Threshold.CI.upper
validation.of.mean.of.latent.EAP[s] <- ssss$mean.of.latent.EAP
validation.of.mean.of.latent.CI.lower[s] <- ssss$mean.of.latent.CI.lower
validation.of.mean.of.latent.CI.upper[s] <- ssss$mean.of.latent.CI.upper
validation.of.Standard.Deviation.latent.EAP[s] <- ssss$Standard.Deviation.latent.EAP
validation.of.Standard.Deviation.latent.CI.lower[s] <- ssss$Standard.Deviation.latent.CI.lower
validation.of.Standard.Deviation.CI.upper[s] <- ssss$Standard.Deviation.CI.upper
validation.of.p.Hitrate.EAP[[s]] <- ssss$p.Hitrate.EAP
validation.of.p.Hitrate.CI.lower[[s]] <- ssss$p.Hitrate.CI.lower
validation.of.p.Hitrate.CI.upper[[s]] <- ssss$p.Hitrate.CI.upper
validation.of.l.FalseRate.EAP[[s]] <- ssss$l.FalseRate.EAP
validation.of.l.FalseRate.CI.lower[[s]] <- ssss$l.FalseRate.CI.lower
validation.of.l.FalseRate.CI.upper[[s]] <- ssss$l.FalseRate.CI.upper
validation.of.loglikelihood.of.model[s] <- ssss$loglikelihood.of.model
}# if convergence is true
color_message(seed,"-th fitting FINISHED-----------------------------------------------------------------------")
}#for seed
# browser()
col.names <- vector()
for (sd in 1:s) {
col.names[sd] <- paste(sd,"-th model",sep = "")
}
colnames(convergent.dataList.as.dataframe) <- col.names
number.of.convergent.model <-s#Here
message("\n ----- Comments for Validation -----\n")
message("\n* Number of all replicated models:",replicate.datset,"\n")
message("\n* Number of convergent models:",s,"\n")
message("\n* Convergence rate:= convergent models / convergent and non convergent models =", round( (s/replicate.datset)*100 , digits = 2),"% \n")
a.truth <- mean.truth/sd.truth
b.truth <- 1/sd.truth
AUC.truth <- stats::pnorm(a.truth/ sqrt( 1+b.truth^2))
l.truth <-l.truth[1:C]
# here error.of.AUC------
error.of.AUC <- validation.of.AUC.EAP - AUC.truth
if (absolute.errors) error.of.AUC <- abs(validation.of.AUC.EAP - AUC.truth)
l1 <-validation.of.z.Threshold.EAP
# l2 <-rep( list(z.truth) , length(unlist( validation.of.z.Threshold.EAP))/length(z.truth) )
l2 <-rep( list(z.truth) , number.of.convergent.model)
# here error.of.z.Threshold.EAP------
# browser()
error.of.z.Threshold.EAP <- Map( `-`,l1 , l2 )
if (absolute.errors) error.of.z.Threshold.EAP <- lapply( Map( `-`,l1 , l2 ),abs)
mean.error.of.z.Threshold.EAP<- rep(NaN, length(error.of.z.Threshold.EAP[[1]]))
sd.error.of.z.Threshold.EAP<- rep(NaN, length(error.of.z.Threshold.EAP[[1]]))
name.z.Threshold.EAP<- rep(NaN, length(error.of.z.Threshold.EAP[[1]]))
for(i.th.column in 1: length(error.of.z.Threshold.EAP[[1]]) ){
mean.error.of.z.Threshold.EAP[i.th.column]<-mean(as.numeric(lapply(error.of.z.Threshold.EAP,"[",n=i.th.column)))
# sd.error.of.z.Threshold.EAP[i.th.column]<-stats::sd(as.numeric(lapply(error.of.z.Threshold.EAP,"[",n=i.th.column)))
sd.error.of.z.Threshold.EAP[i.th.column]<-stats::sd(as.numeric(lapply(validation.of.z.Threshold.EAP,"[",n=i.th.column)))
name.z.Threshold.EAP[i.th.column]<- paste("z[", i.th.column, "]")
}
# browser()
l1 <-validation.of.dz.Threshold.EAP
# l2 <-rep( list(z.truth) , length(unlist( validation.of.z.Threshold.EAP))/length(z.truth) )
l2 <-rep( list(dz.truth) , number.of.convergent.model)
# error ----
# browser()
error.of.dz.Threshold.EAP <- Map( `-`,l1 , l2 )
if (absolute.errors) error.of.dz.Threshold.EAP <- lapply( Map( `-`,l1 , l2 ),abs)
mean.error.of.dz.Threshold.EAP<- rep(NaN, length(error.of.dz.Threshold.EAP[[1]]))
sd.error.of.dz.Threshold.EAP<- rep(NaN, length(error.of.dz.Threshold.EAP[[1]]))
name.dz.Threshold.EAP<- rep(NaN, length(error.of.dz.Threshold.EAP[[1]]))
for(i.th.column in 1: length(error.of.dz.Threshold.EAP[[1]]) ){
mean.error.of.dz.Threshold.EAP[i.th.column]<-mean(as.numeric(lapply(error.of.dz.Threshold.EAP,"[",n=i.th.column)))
# sd.error.of.z.Threshold.EAP[i.th.column]<-stats::sd(as.numeric(lapply(error.of.z.Threshold.EAP,"[",n=i.th.column)))
sd.error.of.dz.Threshold.EAP[i.th.column]<-stats::sd(as.numeric(lapply(validation.of.dz.Threshold.EAP,"[",n=i.th.column)))
name.dz.Threshold.EAP[i.th.column]<- paste("dz[", i.th.column, "]")
}
# error ----
# browser()
if (absolute.errors) error.of.mean.of.latent.EAP <- abs( validation.of.mean.of.latent.EAP - mean.truth)
error.of.mean.of.latent.EAP <- validation.of.mean.of.latent.EAP - mean.truth
if (absolute.errors) error.of.Standard.Deviation.latent.EAP <-abs( validation.of.Standard.Deviation.latent.EAP -sd.truth)
error.of.Standard.Deviation.latent.EAP <- validation.of.Standard.Deviation.latent.EAP -sd.truth
#S.E.
sd.error.of.mean.of.latent.EAP <- stats::sd(validation.of.mean.of.latent.EAP)
sd.error.of.Standard.Deviation.latent.EAP <-stats::sd(validation.of.mean.of.latent.EAP)
sd.error.of.AUC <-stats::sd(validation.of.AUC.EAP)
# sd.error.of.AUC <-stats::sd(error.of.AUC)
name.Gaussian <- append(name.z.Threshold.EAP, c("mean.of.signal", "sd.of.signal","AUC") )
Bias.Gaussian <- append(mean.error.of.z.Threshold.EAP, c(mean(error.of.mean.of.latent.EAP),mean(error.of.Standard.Deviation.latent.EAP),mean(error.of.AUC ) ) )
Std.Err.Gaussian <- append(sd.error.of.z.Threshold.EAP, c(sd.error.of.mean.of.latent.EAP, sd.error.of.Standard.Deviation.latent.EAP, sd.error.of.AUC) )
name.Gaussian <- append(name.Gaussian, name.dz.Threshold.EAP )
Bias.Gaussian <- append(Bias.Gaussian, mean.error.of.dz.Threshold.EAP )
Std.Err.Gaussian <- append(Std.Err.Gaussian,sd.error.of.dz.Threshold.EAP )
if (summary==TRUE) {
message("\n* This is a model parameter for the latent Gaussian distribution.\n")
print(knitr::kable( data.frame(Param.Name=name.Gaussian, ############ Here !!
Bias = Bias.Gaussian, ########### Here !!
S.E. = Std.Err.Gaussian )########### Here !!
))
}
l1 <-validation.of.p.Hitrate.EAP
l2 <-rep( list(p.truth) , length(unlist( validation.of.p.Hitrate.EAP))/length(p.truth) )
error.of.p.Hitrate.EAP <- Map( `-`,l1 , l2 )
if (absolute.errors) error.of.p.Hitrate.EAP <- lapply( Map( `-`,l1 , l2 ),abs)
mean.error.p.Hitrate.EAP<- rep(NaN, length(error.of.p.Hitrate.EAP[[1]]))
sd.error.p.Hitrate.EAP<- rep(NaN, length(error.of.p.Hitrate.EAP[[1]]))
name.p.Hitrate.EAP<- rep(NaN, length(error.of.p.Hitrate.EAP[[1]]))
for(i.th.column in 1: length(error.of.p.Hitrate.EAP[[1]]) ){
# mean ----
mean.error.p.Hitrate.EAP[i.th.column] <- mean(as.numeric(lapply(error.of.p.Hitrate.EAP,"[",n=i.th.column)))
sd.error.p.Hitrate.EAP[i.th.column] <- stats::sd(as.numeric(lapply(error.of.p.Hitrate.EAP,"[",n=i.th.column)))
name.p.Hitrate.EAP[i.th.column] <- paste("p[", i.th.column, "]")
}
if (summary==TRUE) {
print(knitr::kable( data.frame(Param.Name=name.p.Hitrate.EAP,########### Here !!
Bias = mean.error.p.Hitrate.EAP,########### Here !!
S.E. = sd.error.p.Hitrate.EAP )########### Here !!
)
)
}
# browser()
l1 <-validation.of.l.FalseRate.EAP
l2 <-rep( list(l.truth) , length(unlist( validation.of.l.FalseRate.EAP))/length(l.truth) )
# error ----
error.of.l.FalseRate.EAP <- Map( `-`,l1 , l2 )
if (absolute.errors) error.of.l.FalseRate.EAP <- lapply( Map( `-`,l1 , l2 ),abs)
# mean ----
mean.error.l.FalseRate.EAP<- rep(NaN, length(error.of.l.FalseRate.EAP[[1]]))
sd.error.l.FalseRate.EAP<- rep(NaN, length(error.of.l.FalseRate.EAP[[1]]))
name.l.FalseRate.EAP<- rep(NaN, length(error.of.l.FalseRate.EAP[[1]]))
for(i.th.column in 1: length(error.of.l.FalseRate.EAP[[1]]) ){
# mean ----
mean.error.l.FalseRate.EAP[i.th.column]<-mean(as.numeric(lapply(error.of.p.Hitrate.EAP,"[",n=i.th.column)))
sd.error.l.FalseRate.EAP[i.th.column]<-stats::sd(as.numeric(lapply(error.of.p.Hitrate.EAP,"[",n=i.th.column)))
name.l.FalseRate.EAP[i.th.column]<- paste("lambda[", i.th.column, "]")
}
if (summary==TRUE) {
print(knitr::kable( data.frame(Param.Name=name.l.FalseRate.EAP,########### Here !!
Bias = mean.error.l.FalseRate.EAP,########### Here !!
S.E. = sd.error.l.FalseRate.EAP )) )########### Here !!
# browser()
}
Name.of.all.param <- append(name.Gaussian,
append(name.p.Hitrate.EAP,
name.l.FalseRate.EAP)
)
# error ----
Bias.of.all.param <- append(
Bias.Gaussian,
append(
mean.error.p.Hitrate.EAP,
mean.error.l.FalseRate.EAP)
)
S.E.of.all.param <- append(Std.Err.Gaussian,
append(
sd.error.p.Hitrate.EAP,
sd.error.l.FalseRate.EAP)
)
results <- data.frame( Name.of.all.param=Name.of.all.param,
Bias.of.all.param=Bias.of.all.param,
S.E.of.all.param=S.E.of.all.param
)
# error ----
error.of.mean.of.latent.EAP <- validation.of.mean.of.latent.EAP - mean.truth
error.of.Standard.Deviation.latent.EAP <- validation.of.Standard.Deviation.latent.EAP -sd.truth
return_value<- list(
results=results,
Name.of.all.param = Name.of.all.param,
Bias.of.all.param = Bias.of.all.param,
S.E.of.all.param = S.E.of.all.param,
vector.of.dataList=vector.of.dataList,
convergent.dataList = convergent.dataList,
convergent.dataList.as.dataframe = convergent.dataList.as.dataframe,
############################################################# Convergence number and all (non conv and conv) number
number.of.convergent.model =number.of.convergent.model,
replicate.datset=replicate.datset,
#Here
fit=fit,
ModifiedPoisson=ModifiedPoisson,
mean.truth=mean.truth,
sd.truth=sd.truth,
p.truth=p.truth,
l.truth=l.truth,
z.truth=z.truth,
AUC.truth=AUC.truth,
validation.of.AUC.EAP =validation.of.AUC.EAP,
validation.of.AUC.CI.lower =validation.of.AUC.CI.lower,
validation.of.AUC.CI.upper =validation.of.AUC.CI.upper,
validation.of.z.Threshold.EAP = validation.of.z.Threshold.EAP,
validation.of.z.Threshold.CI.lower = validation.of.z.Threshold.CI.lower,
validation.of.z.Threshold.CI.upper = validation.of.z.Threshold.CI.upper,
validation.of.mean.of.latent.EAP = validation.of.mean.of.latent.EAP,
validation.of.mean.of.latent.CI.lower = validation.of.mean.of.latent.CI.lower,
validation.of.mean.of.latent.CI.upper = validation.of.mean.of.latent.CI.upper,
validation.of.Standard.Deviation.latent.EAP = validation.of.Standard.Deviation.latent.EAP,
validation.of.Standard.Deviation.latent.CI.lower = validation.of.Standard.Deviation.latent.CI.lower,
validation.of.Standard.Deviation.CI.upper = validation.of.Standard.Deviation.CI.upper,
validation.of.p.Hitrate.EAP = validation.of.p.Hitrate.EAP,
validation.of.p.Hitrate.CI.lower = validation.of.p.Hitrate.CI.lower,
validation.of.p.Hitrate.CI.upper = validation.of.p.Hitrate.CI.upper,
validation.of.l.FalseRate.EAP = validation.of.l.FalseRate.EAP,
validation.of.l.FalseRate.CI.lower = validation.of.l.FalseRate.CI.lower,
validation.of.l.FalseRate.CI.upper = validation.of.l.FalseRate.CI.upper,
validation.of.loglikelihood.of.model = validation.of.loglikelihood.of.model,
# here errors ------
error.of.AUC=error.of.AUC,
error.of.z.Threshold.EAP = error.of.z.Threshold.EAP,
error.of.mean.of.latent.EAP =error.of.mean.of.latent.EAP,
error.of.Standard.Deviation.latent.EAP =error.of.Standard.Deviation.latent.EAP,
error.of.p.Hitrate.EAP =error.of.p.Hitrate.EAP,
error.of.l.FalseRate.EAP =error.of.l.FalseRate.EAP,
dataList=list.of.dataList
)
invisible( return_value )
}#function
#'@title Validation via replicated datasets from a model at a given model parameter
#'@description
#' Print for a given true parameter,
#' a errors of estimates from replicated dataset.
#'
#' Also print a standard error which is the variance of estimates.
#'
#'
#'
#' Suppose that \eqn{\theta_0} is a given true model parameter with a given number of images \eqn{N_I} and a given number of lesions \eqn{N_L}, specified by user.
#' \describe{
#' \item{ \strong{(I)} }{}
#' \item{ \strong{(I.1)} Synthesize a collection of dataset \eqn{D_k} (\eqn{k=1,2,...,K}) from a likelihood (model) at a given parameter \eqn{\theta_0}, namely }{ \eqn{D_k \sim} likelihood( \eqn{\theta_0}). }
#' \item{ \strong{(I.2)} Replicates \eqn{K} models fitted to each dataset \eqn{D_k} (\eqn{k=1,2,...,K}), namely, draw MCMC samples \eqn{\{ \theta_i (D_k);i=1,...,I\}} from each posterior of the dataset \eqn{D_k}, namely }{ \eqn{ \theta _i(D_k) \sim \pi(|D_k)}. }
#' \item{ \strong{(I.3)} Calculate posterior means for the set of data \eqn{D_k} (\eqn{k=1,2,...,K}), namely }{ \eqn{ \bar{\theta}(D_k) := \frac{1}{I} \sum_i \theta_i(D_k) }. }
#' \item{ \strong{(I.4)} Calculates error for each dataset \eqn{D_k} }{ \eqn{\epsilon_k}:= Truth - estimates = \eqn{\theta_0 - \bar{\theta}(D_k)}. }
#' \item{ \strong{(II)} Calculates mean of errors over all datasets \eqn{D_k} (\eqn{k=1,2,...,K}) }{ mean of errors \eqn{ \bar{\epsilon}(\theta_0,N_I,N_L)}= \eqn{ \frac{1}{K} \sum \epsilon_k }. }
#' \item{ NOTE }{ We note that if a fitted model does not converge,( namely R hat is far from one), then it is omiited from this calculation.}
#' \item{ \strong{(III) } Calculates mean of errors for various number of lesions and images }{ mean of errors \eqn{ \bar{\epsilon}(\theta_0,N_I,N_L)} }
#' }
#'
#' For example, if \eqn{(N_I^1,N_L^1),(N_I^2,N_L^2),(N_I^3,N_L^3),...,(N_I^m,N_L^m)}, then
#' \eqn{ \bar{\epsilon}(\theta_0,N_I^1,N_L^1)},
#' \eqn{ \bar{\epsilon}(\theta_0,N_I^2,N_L^2)},
#' \eqn{ \bar{\epsilon}(\theta_0,N_I^3,N_L^3)},...,
#' \eqn{ \bar{\epsilon}(\theta_0,N_I^m,N_L^m)} are calculated.
#'
#' To obtain precise error,
#' The number of replicated fitted models (denoted by \eqn{K}) should be large enough.
#' If \eqn{K} is small, then it causes a bias.
#' \eqn{K} = \code{replicate.datset}: a variable of the function \code{error_srsc}.
#'
#'
#'
#'
#'
#'
#'
#'
#describe
#'
#'
#' Running this function, we can see that the error
#' \eqn{ \bar{\epsilon}(\theta_0,N_I,N_L)} decreases
#' monotonically as a given number of images \eqn{N_I}
#' or a given number of lesions \eqn{N_L} increases.
#'
#' Also, the scale of error also will be found. Thus this function can show how our estimates are correct.
#' Scale of error differs for each componenet of model parameters.
#'
#'
#'
#' Revised 2019 August 28
#'
#'
#'
#' @details
#'
#'
#' In Bayesian inference,
#' if sample size is large,
#' then posterior tends to the Dirac measure.
#' So, the error and variance of estimates
#' should be tends to zero as sample size tends to infinity.
#'
#' This function check this phenomenen.
#'
#' If model has problem, then it contains some non-decreasing vias
#' with respect to sample size.
#'
#' Revised 2019 Nov 1
#'
#'
#'
#' Provides a reliability of our posterior mean estimates.
#' Using this function, we can
#' find what digit makes sence.
#'
#' In the real world, the data for modality
#' comparison or observer performan evaluation is
#' 100 images or 200 images. In such scale data,
#' any estimate of AUC will contain error at most 0.0113....
#' So, the value of AUC should round in 0.XXX
#' and not 0.XXXX or 0.XXXXX or more. Since
#' error is 0.00113... and hence 4 digit
#' or more digit is meaningless. In such manner,
#' we can analyize the errors.
#'
#'We note that if we increase the
#'number of images or lesinons,
#' the errors decrease.
#'
#' For example, if we use 20000 images in FROC trial,
#' then the error of AUC will be 0.0005... and thus,
#' and so on. Thus large number of images gives
#' us more reliable AUC. However the radiologist
#' cannot read such large (20000) images.
#'
#' Thus, the error will be 0.00113...
#'
#'
#' If the number of images are given before hand and moreover if we obtains the estimates,
#' then we can run this function using these two, we can find the estimated errors by simulation.
#' Of course, the esimates is not the truth, but roughly speaking, if we assume that
#' the estimates is not so far from truth, and the error analysis is rigid with respect to
#' changing the truth, then we can say using estimates as truth, the result of this error analysis can be regarded as an
#' actual error.
#'
#'
#' I want to go home. Unfortunatly, my house is ...
#'
#'
#'
#'@return Replicated datasets, estimates,
#' errors,...etc I made this program 1 years ago?
#' and now I forget ... the precise return values.
#'When I see today, 2019 August. It retains too many return
#'values to explain all of them.
#@param ----
#'@param NLvector A vector of positive integers,
#' indicating a collection of numbers of Lesions.
#@param NIvector Number of Images.
#'@param ratio A positive \strong{\emph{rational}} number,
#' with which Number of Images is determined by the formula:
#' (number of images) = \code{ratio} times (numbser of lesions).
#' Note that in calculation, it rounds \code{ratio * NLvector } to an integer.
#'@inheritParams fit_Bayesian_FROC
#'@inheritParams DrawCurves
#'@inheritParams validation.dataset_srsc
#'@examples
#' \dontrun{
#'#========================================================================================
#'# 0) 0-th example
#'#========================================================================================
#'
#'
#' datasets <-error_srsc(
#' NLvector = c(100,10000,1000000),
#' ite = 2222
#' )
#'
#'
#' # By the following, we can extract only datasets whose
#' # model has converged.
#' datasets$convergent.dataList.as.dataframe
#'
#'
#'
#'
#'#========================================================================================
#'# 1) 1-st example
#'#========================================================================================
#'# Long width is required in R console.
#'
#'
#'
#' datasets <-error_srsc(NLvector = c(
#' 50L,
#' 111L,
#' 11111L
#' ),
#' # NIvector,
#' ratio=2,
#' replicate.datset =3,
#' ModifiedPoisson = FALSE,
#' mean.truth=0.6,
#' sd.truth=5.3,
#' z.truth =c(-0.8,0.7,2.38),
#' ite =2222
#' )
#'
#'
#'
#'#========================================================================================
#'# 2) Plot the error of AUC with respect to NI
#'#========================================================================================
#'
#'
#' a <-error_srsc(NLvector = c(
#' 33L,
#' 50L,
#' 111L,
#' 11111L
#' ),
#' # NIvector,
#' ratio=2,
#' replicate.datset =3,
#' ModifiedPoisson = FALSE,
#' mean.truth=0.6,
#' sd.truth=5.3,
#' z.truth =c(-0.8,0.7,2.38),
#' ite =2222
#' )
#'
#'
#'
#'
#'
#'
#' aa <- a$Bias.for.various.NL
#'
#'
#' error.of.AUC <- aa[8,]
#' y <- subset(aa[8,], select = 2:length(aa[8,]))
#' y <- as.numeric(y)
#' y <- abs(y)
#' upper_y <- max(y)
#' lower_y <- min(y)
#'
#' x <- 1:length(y)
#'
#'
#'
#' plot(x,y, ylim=c(lower_y, upper_y))
#'
#'# From this plot, we cannot see whether the error has decreased or not.
#'# Thus, we replot with the log y-axis, the we will see that the error
#'# has decreased with respect to number of images and lesions.
#'
#' ggplot(data.frame(x=x,y=y), aes(x = x, y = y)) +
#' geom_line() +
#' geom_point() +
#' scale_y_log10()
#'
#'# Revised 2019 Sept 25
#'
#'
#' # General print of log scale
#' df<-data.frame(x=c(10,100,1000,10,100,1000),
#' y=c(1100,220000,33000000,1300,240000,36000000),
#' group=c("1","1","1","2","2","2")
#' )
#'
#' ggplot2::ggplot(df, aes(x = x, y = y, shape = group)) +
#' ggplot2::geom_line(position = position_dodge(0.2)) + # Dodge lines by 0.2
#' ggplot2::geom_point(position = position_dodge(0.2), size = 4)+ # Dodge points by 0.2
#' ggplot2::scale_y_log10()+
#' ggplot2::scale_x_log10()
#'
#'
#'#========================================================================================
#'# 2) Add other param into plot plain of the error of AUC with respect to NI
#'#========================================================================================
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#' a <-error_srsc(NLvector = c(
#' 111L,
#' 11111L
#' ),
#' # NIvector,
#' ratio=2,
#' replicate.datset =3,
#' ModifiedPoisson = FALSE,
#' mean.truth=0.6,
#' sd.truth=5.3,
#' z.truth =c(-0.8,0.7,2.38),
#' ite =2222
#' )
#' aa <- a$Bias.for.various.NL
#'
#'
#' error.of.AUC <- aa[8,]
#' y1 <- subset(aa[8,], select = 2:length(aa[8,]))
#' y1 <- as.numeric(y1)
#' y1 <- abs(y1)
#'
#' LLL <-length(y1)
#'
#' y2 <- subset(aa[7,], select = 2:length(aa[7,]))
#' y2 <- as.numeric(y2)
#' y2 <- abs(y2)
#'
#' y <- c(y1,y2)
#'
#'
#' upper_y <- max(y)
#' lower_y <- min(y)
#'
#' group <- rep(seq(1,2,1),1 , each=LLL)
#' x <- rep(seq(1,LLL,1),2 , each=1)
#' group <- as.character(group)
#' df <- data.frame(x=x,y=y,group=group)
#'
#'
#' ggplot2::ggplot(df, aes(x = x, y = y, shape = group)) +
#' ggplot2::geom_line(position = position_dodge(0.2)) + # Dodge lines by 0.2
#' ggplot2::geom_point(position = position_dodge(0.2), size = 4)+ # Dodge points by 0.2
#' ggplot2::scale_y_log10()
#' # ggplot2::scale_x_log10()
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'#========================================================================================
#'# Confidence level = 4
#'#========================================================================================
#'
#'
#'
#'
#'
#' datasets <-error_srsc(NLvector = c(
#' 111L,
#' 11111L
#' ),
#' # NIvector,
#' ratio=2,
#' replicate.datset =3,
#' ModifiedPoisson = FALSE,
#' mean.truth=-0.22,
#' sd.truth=5.72,
#' z.truth =c(-0.46,-0.20,0.30,1.16),
#' ite =2222
#' )
#'
#'
#'
#'
#'
#' error_srsc_variance_visualization(datasets)
#'
#'# The parameter of model is 7 in which the ggplot2 fails with the following warning:
#'
#'
#'# The shape palette can deal with a maximum of 6 discrete values because more than 6
#'# becomes difficult to
#'# discriminate; you have 7. Consider specifying shapes manually if you must have them.
#'
#'
# ####1#### ####2#### ####3#### ####4#### ####5#### ####6#### ####7#### ####8#### ####9####
#'#========================================================================================
#'# NaN ... why? 2021 Dec
#'#========================================================================================
#'
#' fits <- validation.dataset_srsc()
#'
#' f <-fits$fit[[2]]
#' rstan::extract(f)$dl
#' sum(rstan::extract(f)$dl)
#' Is.nan.in.MCMCsamples <- as.logical(!prod(!is.nan(rstan::extract(f)$dl)))
#' rstan::extract(f)$A[525]
#' a<-rstan::extract(f)$a
#' b<-rstan::extract(f)$b
#'
#'
#'
#' Phi( a[525]/sqrt(b[525]^2+1) )
#' a[525]/sqrt(b[525]^2+1)
#' Phi( a/sqrt(b^2+1) )
#' x<-rstan::extract(f)$dl[2]
#'
#' a<-rstan::extract(f)$a
#' b<-rstan::extract(f)$b
#'
#' a/(b^2+1)
#' Phi(a/(b^2+1))
#' mean( Phi(a/(b^2+1)) )
#'
#' #'
#'
#'
#'}# dontrun
#' @export
# OUt put can be made by TeX source.
error_srsc <-function(NLvector = c(
100L,
10000L,
1000000L
),
# NIvector,
ratio=2,# NUmber of image is not so important, since when we consider ModifiedPoisson = TRUE
replicate.datset =3,
ModifiedPoisson = FALSE,
mean.truth=0.6,
sd.truth=5.3,
z.truth =c(-0.8,0.7,2.38),
ite =2222,
cha =1,
verbose = FALSE
){
C <- length(z.truth)
d <- list()
Bias.of.all.param <- list()
S.E.of.all.param <- list()
Bias.of.all.param.with.NL.NI <- list()
S.E.of.all.param.with.NL.NI <- list()
Name.of.all.param.with.NL.NI <- vector()
col.names_for_SE <- vector()
col.names_for_SE[1] <-paste("Variance of estimates")
col.names <- vector()
col.names[1] <-paste("Error of estimates")
convergent.dataList.as.dataframe <- list()
number.of.convergent.model<- vector()
number.of.replicate.datset<- vector()
# col.names[1] <-paste("")
s <- 0
base_size <- 0
for (nl in NLvector) {
s<- s+1
ni<- floor(ratio*nl)
#Create the replicating datasets for different number of lesions
# validation.dataset_srsc -----
d[[s]]<- validation.dataset_srsc(
replicate.datset =replicate.datset,
ModifiedPoisson = ModifiedPoisson,
mean.truth=mean.truth,
sd.truth=sd.truth,
z.truth =z.truth,
NL=ni,
NI=nl,
ite =ite,
cha =cha,
base_size = base_size ,
verbose = verbose,
summary=FALSE,
serial.number= paste( "The number of Lesion = ",nl, " and The number of images",ni, ". \n\n* The (" , s, "/", length(NLvector), ")")
)
base_size <- base_size + size_of_return_value(d[[s]], is_return_value = FALSE,base_size = base_size,summary = FALSE)
# Here ------
number.of.convergent.model[s] <- d[[s]]$number.of.convergent.model#Here
number.of.replicate.datset[s] <- d[[s]]$replicate.datset#Here
Bias.of.all.param[[s]] <-d[[s]]$Bias.of.all.param
S.E.of.all.param[[s]] <-d[[s]]$S.E.of.all.param
Bias.of.all.param.with.NL.NI[[s]] <- append(append(append(append(ni,nl),
Bias.of.all.param[[s]]),
number.of.convergent.model[s]),
number.of.replicate.datset[s])
S.E.of.all.param.with.NL.NI[[s]] <- append(append(append(append(ni,nl),# 2019 Sept 6
S.E.of.all.param[[s]]),# 2019 Sept 6
number.of.convergent.model[s]),# 2019 Sept 6
number.of.replicate.datset[s])# 2019 Sept 6
col.names[s+1] <-paste(s,"-th validation",sep = "")
if(s==1) col.names[s+1] <-paste(s,"-st validation",sep = "")
if(s==2) col.names[s+1] <-paste(s,"-nd validation",sep = "")
if(s==3) col.names[s+1] <-paste(s,"-rd validation",sep = "")
col.names_for_SE[s+1] <-paste(s,"-th validation",sep = "")# 2019 Sept 6
if(s==1) col.names_for_SE[s+1] <-paste(s,"-st validation",sep = "")# 2019 Sept 6
if(s==2) col.names_for_SE[s+1] <-paste(s,"-nd validation",sep = "")# 2019 Sept 6
if(s==3) col.names_for_SE[s+1] <-paste(s,"-rd validation",sep = "")# 2019 Sept 6
convergent.dataList.as.dataframe[[s]] <- d[[s]]$convergent.dataList.as.dataframe
}# for nl in NLvector
list.names <- vector()
for (sd in 1:s) {
list.names [sd] <- paste(sd,"-th dataset of convergent model ",sep = "")
}
names(convergent.dataList.as.dataframe) <- list.names
# Name.of.all.param.with.NL.NI<- append(append("Number of Images",
# "Number of Lesions"),
# d[[1]]$Name.of.all.param )
Name.of.all.param.with.NL.NI<- append( append( append(append("Number of Images",
"Number of Lesions"),
d[[1]]$Name.of.all.param ),
"number of convergent models"),#Here
"number of replicated datsets")#Here
# browser()
l <-list(
Name.of.all.param.with.NL.NI = Name.of.all.param.with.NL.NI,
Bias.of.all.param.with.NL.NI = Bias.of.all.param.with.NL.NI
)#list
ll <-list(# 2019 Sept 6
Name.of.all.param.with.NL.NI = Name.of.all.param.with.NL.NI,# 2019 Sept 6
S.E.of.all.param.with.NL.NI = S.E.of.all.param.with.NL.NI# 2019 Sept 6
)#list# 2019 Sept 6
# browser()
options(scipen=100)
#Here, list l is converted, and each component of list become column vector!! This is very good!
Bias.for.various.NL<-as.data.frame(l)
S.E.for.various.NL<-as.data.frame(ll)# 2019 Sept 6
names(Bias.for.various.NL) <- col.names
names(S.E.for.various.NL)<- col.names_for_SE # 2019 Sept 6
message(
"
==================================================================
* Each column is independent.
* Each row means the mean of the difference of [ truth - estimate ]
* If each number is small, then it means that each error is small.
* In my model, the most important quantity is a parameter of AUC,
whose error is the most important.
")
message("\n* Variance of estimates, (Standard Error) \n")# 2019 Sept 6
print(knitr::kable(S.E.for.various.NL) )# 2019 Sept 6
message("\n* Error between estimates and specified true parameter \n")
print(knitr::kable(Bias.for.various.NL) )
return_value <- list(
Bias.for.various.NL =Bias.for.various.NL,
S.E.for.various.NL=S.E.for.various.NL,
Name.of.all.param.with.NL.NI = Name.of.all.param.with.NL.NI,
Bias.of.all.param.with.NL.NI = Bias.of.all.param.with.NL.NI,
S.E.of.all.param.with.NL.NI = S.E.of.all.param.with.NL.NI,
C=C,
length.of.NLvector = length( NLvector),
NLvector=NLvector,
ratio=ratio,# NUmber of image is not so important, since when we consider ModifiedPoisson = TRUE
replicate.datset =replicate.datset,
ModifiedPoisson = ModifiedPoisson,
mean.truth=mean.truth,
sd.truth=sd.truth,
z.truth =z.truth,
Number.of.MCMC.iteration =ite,
number.of.convergent.model =number.of.convergent.model,#Here
number.of.replicate.datset =number.of.replicate.datset, #Here
datasets =d,
convergent.dataList.as.dataframe=convergent.dataList.as.dataframe
)
size_of_return_value(
return_value
)
invisible(
return_value
)#invisible
}
# datasets$datasets[[2]]$dataList
#' @title Visualization for Error of Estimator
#' @description
#' The function plot the graph of errors with respect to sample sizes.
#'
#'
#'\strong{\emph{Error plot}}
#'\describe{
#'\item{ \strong{\emph{x-axis}} }{ Sample sizes }
#'\item{ \strong{\emph{y-axis}} }{ Error for each parameter }
#'}
#'
#' @param return.value.of_error_srsc A return value
#' of the function \code{\link{error_srsc}()}.
#' @param log_scale_x.axis A logical,
#' whether x axis is log scale or not.
#' @return A long format dataframe of
#' error and its parameter name
#' @export
#' @seealso \link{error_srsc_variance_visualization}
#'
#' @examples
#'
#' # General plot
#'
#' df <- data.frame(x=runif(100),y=runif(100),g= as.factor(rep(1:5,10)))
#'
#' ggplot(df, aes(x = x, y = y, shape = g)) +
#' geom_point(size = 3) +
#' scale_shape_manual(values = c(1,2,3,4,5,6,7,8,9))
#'
#'
#'
#'
#'
#'
#'
#'
#'
#' df <- data.frame(x=runif(100),y=runif(100),g= as.factor(rep(1:25,4)))
#'
#' # Use slightly larger points and use custom values for the shape scale
#'
#'
#' ggplot(df, aes(x = x, y = y, shape = g)) +
#' geom_point(size = 3) +
#' scale_shape_manual(values = c(1,2,3,4,5,6,7,8,9,10,
#' 11,12,13,14,15,16,17,18,19,20,21,22,23,24,25))
#'
#' \dontrun{
#' a <- error_srsc()
#'
#' error_srsc_error_visualization(a)
#'
#'
# ####1#### ####2#### ####3#### ####4#### ####5#### ####6#### ####7#### ####8#### ####9####
#'#========================================================================================
#'# In case of C = 4, arbitrary C is available.
#'#========================================================================================
#'
#' a <-error_srsc(NLvector = c(
#' 100,
#' 10000,
#' 1000000
#' ),
#' ratio=2,
#' replicate.datset =2,
#' ModifiedPoisson = FALSE,
#' mean.truth=0.6,
#' sd.truth=5.3,
#' z.truth =c(-0.8,0.7,2.38,3), # Here we use the C=4
#' ite =500
#' )
#'
#' error_srsc_error_visualization(a)
#' error_srsc_variance_visualization(a)
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
# ####1#### ####2#### ####3#### ####4#### ####5#### ####6#### ####7#### ####8#### ####9####
#'#========================================================================================
#'# In case of C = 7, arbitrary C is available.
#'#========================================================================================
#'
#'
#'
#'
#'
#'
#'
#' #'
#' a <-error_srsc(NLvector = c(
#' 100,
#' 10000,
#' 100000
#' ),
#' ratio=2,
#' replicate.datset =2,
#' ModifiedPoisson = FALSE,
#' mean.truth=0.6,
#' sd.truth=5.3,
#' z.truth =c(-0.8,0.7,2.38,3,3.4,3.6,3.8), # Here we use the C=7
#' ite =500
#' )
#'
#' error_srsc_error_visualization(a)
#' error_srsc_variance_visualization(a)
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#' }
#'
error_srsc_error_visualization <- function(return.value.of_error_srsc,
log_scale_x.axis = TRUE
){
a<-return.value.of_error_srsc
aa <- a$Bias.for.various.NL
NLvector <- a$NLvector
length.of.NLvector <- a$length.of.NLvector
C <- a$C
LLL <-length.of.NLvector
y.list <- list()
group.list <-list()
for (cd in 1:(C+3)) {
y.list[[cd]] <- vector()
y.list[[cd]] <- abs(
as.numeric(
subset(aa[cd+2,], select = 2:length(aa[cd+2,]))
)
)
# group.list[[cd]] <- subset(aa[,1], subset = 2:(C+2),select =1 )
}
for (cd in 1:C) {
group.list[[cd]] <- vector()
for (nnn in 1:LLL) {
group.list[[cd]][nnn] <- paste("z[",cd,"]")
}
}
group.list[[C+1]] <- vector()
group.list[[C+2]] <- vector()
group.list[[C+3]] <- vector()
for (nnn in 1:LLL) {
group.list[[C+1]][nnn] <- "mean.of.signal"
group.list[[C+2]][nnn] <- "sd.of.signal"
group.list[[C+3]][nnn] <- "AUC"
}
y <- unlist(y.list)
group <- unlist(group.list)
# upper_y <- max(y)
# lower_y <- min(y)
# group <- rep(seq(1,LLL,1),1 , each=C+2)
# group <- as.character(group)
if(log_scale_x.axis==FALSE) x <- rep(seq(1,LLL,1),C+3 , each=1)
if(log_scale_x.axis==TRUE) x <- rep(NLvector ,C+3 , each=1)
# browser()
# browser()
# browser()
df <- data.frame(x=x,y=y,group=group)
# browser()
if (log_scale_x.axis==TRUE) {
print(
ggplot2::ggplot(df, ggplot2::aes(x = x, y = y, shape = group)) +
ggplot2::geom_line(position = ggplot2::position_dodge(0.2)) + # Dodge lines by 0.2
ggplot2::geom_point(position = ggplot2::position_dodge(0.2), size = 4)+ # Dodge points by 0.2
ggplot2::scale_shape_manual(values = c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25))+ # This is a valiation of dots in the plot, and important in case of the number of confidence level is greater than 4
ggplot2::scale_y_log10()+
ggplot2::scale_x_log10()
)
}
if (log_scale_x.axis==FALSE) {
print( ggplot2::ggplot(df, ggplot2::aes(x = x, y = y, shape = group)) +
ggplot2::geom_line(position = ggplot2::position_dodge(0.2)) + # Dodge lines by 0.2
ggplot2::geom_point(position = ggplot2::position_dodge(0.2), size = 4)+ # Dodge points by 0.2
ggplot2::scale_shape_manual(values = c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25))+ # This is a valiation of dots in the plot, and important in case of the number of confidence level is greater than 4
ggplot2::scale_y_log10()
#+ ggplot2::scale_x_log10()
)
}
# browser()
}
#' @title Visualization Of variance Analysis
#'
#' @param return.value.of_error_srsc A return value
#' of the function \code{\link{error_srsc}()}.
#' @param log_scale_x.axis A logical, whether x axis is log scale.
#' @return A long format dataframe of error and its parameter name
#' @export
#'
#' @examples
#' \dontrun{
#' a <- error_srsc()
#'
#' error_srsc_variance_visualization(a)
#'
#'
#'
#' a <- error_srsc(replicate.datset = 10)
#' error_srsc_variance_visualization(a)
#'
#' }
#'
error_srsc_variance_visualization <- function(return.value.of_error_srsc,
log_scale_x.axis = TRUE
){
a<-return.value.of_error_srsc
aa <- a$S.E.for.various.NL
NLvector <- a$NLvector
length.of.NLvector <- a$length.of.NLvector
C <- a$C
LLL <-length.of.NLvector
y.list <- list()
group.list <-list()
for (cd in 1:(C+3)) {
y.list[[cd]] <- vector()
y.list[[cd]] <- abs(
as.numeric(
subset(aa[cd+2,], select = 2:length(aa[cd+2,]))
)
)
# group.list[[cd]] <- subset(aa[,1], subset = 2:(C+2),select =1 )
}
for (cd in 1:C) {
group.list[[cd]] <- vector()
for (nnn in 1:LLL) {
group.list[[cd]][nnn] <- paste("z[",cd,"]")
}
}
group.list[[C+1]] <- vector()
group.list[[C+2]] <- vector()
group.list[[C+3]] <- vector()
for (nnn in 1:LLL) {
group.list[[C+1]][nnn] <- "mean.of.signal"
group.list[[C+2]][nnn] <- "sd.of.signal"
group.list[[C+3]][nnn] <- "AUC"
}
y <- unlist(y.list)
group <- unlist(group.list)
# upper_y <- max(y)
# lower_y <- min(y)
# group <- rep(seq(1,LLL,1),1 , each=C+2)
# group <- as.character(group)
if(log_scale_x.axis==FALSE) x <- rep(seq(1,LLL,1),C+3 , each=1)
if(log_scale_x.axis==TRUE) x <- rep(NLvector ,C+3 , each=1)
# browser()
# browser()
# browser()
df <- data.frame(x=x,y=y,group=group)
# browser()
if (log_scale_x.axis==TRUE) {
print(
ggplot2::ggplot(df, ggplot2::aes(x = x, y = y, shape = group)) +
ggplot2::geom_line(position = ggplot2::position_dodge(0.2)) + # Dodge lines by 0.2
ggplot2::geom_point(position = ggplot2::position_dodge(0.2), size = 4)+ # Dodge points by 0.2
ggplot2::scale_shape_manual(values = c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25))+ # This is a valiation of dots in the plot, and important in case of the number of confidence level is greater than 4
ggplot2::scale_y_log10()+
ggplot2::scale_x_log10()
)
}
if (log_scale_x.axis==FALSE) {
print( ggplot2::ggplot(df, ggplot2::aes(x = x, y = y, shape = group)) +
ggplot2::geom_line(position = ggplot2::position_dodge(0.2)) + # Dodge lines by 0.2
ggplot2::geom_point(position = ggplot2::position_dodge(0.2), size = 4)+ # Dodge points by 0.2
ggplot2::scale_shape_manual(values = c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25))+ # This is a valiation of dots in the plot, and important in case of the number of confidence level is greater than 4
ggplot2::scale_y_log10()
#+ ggplot2::scale_x_log10()
)
}
# browser()
}
#'@title Draw Curves for validation dataset
#'@description drawing curves.
#'
#'
#' \strong{Red curve} indicates an FROC curve of truth parameter.
#'
#' \strong{Other curves} are drawn using replicated estimates.
#'
#'
#'@inheritParams DrawCurves
#'@inheritParams DrawCurves_MRMC_pairwise
#'@inheritParams fit_Bayesian_FROC
#'@param validation.data This is a return value of
#' the function \code{validation.dataset_srsc}.
#'
#'@return \code{NULL}
#'
#'
#'@examples
#'
#' \dontrun{
#'#--------------------------------------------------------------------------------------
#'# 1) Draw the curve for each replicated dataset
#'#--------------------------------------------------------------------------------------
#'
#' datasets <- validation.dataset_srsc()
#'
#' validation.draw_srsc(datasets)
#'
#'
#'
#'
#'#--------------------------------------------------------------------------------------
#'# 1) Draw the curve for each replicated dataset
#'#--------------------------------------------------------------------------------------
#'
#'
#'
#' datasets <- validation.dataset_srsc(replicate.datset = 5)
#'
#'
#' validation.draw_srsc(datasets)
#'
#'
#'}# dottest
#' @export
validation.draw_srsc <-function(validation.data,
mesh.for.drawing.curve=11111,
upper_y=1,
DrawFROCcurve=TRUE
){
replicate.datset <- as.integer(validation.data$replicate.datset)
s<-vector()
for (seed in 1:replicate.datset ) {if (validation.data$fit[[seed]]@convergence == TRUE) {
s <- append(s, validation.data$fit[[seed]]@metadata$ff )
}}
upper_x <-max(s) ### Nice! I was great 2019 September 1
s<-0
for (seed in 1:replicate.datset ) {
if (validation.data$fit[[seed]]@convergence == TRUE) {
s<-s+1
validation.data$fit[[seed]]@index <-s
DrawCurves(validation.data$fit[[seed]],
new.imaging.device = FALSE,
title = FALSE,
DrawFROCcurve=DrawFROCcurve,
indexCFPCTP = s,
upper_x=upper_x,
upper_y=upper_y,
# cex=0.04
)
}
if (validation.data$fit[[seed]]@convergence == FALSE) {
message("\n* The ",seed,"-th curve are not drawn, since the model did not converge.\n",sep = "")
}
}
mean.truth <- validation.data$mean.truth
sd.truth <- validation.data$sd.truth
a.truth <- mean.truth/sd.truth
b.truth <- 1/sd.truth
if(validation.data$fit[[1]]@studyDesign == "srsc.per.image"){ xlabel = 'mean of cumulative false positives per image' }
if(validation.data$fit[[1]]@studyDesign == "srsc.per.lesion"){ xlabel = 'mean of cumulative false positives per lesion' }
set.seed(1);ll<- stats::rchisq(mesh.for.drawing.curve, 1)
lll<- 0.99+ll
x.FROC<-append(ll,lll)
y.FROC.truth <- 1 - stats::pnorm( b.truth*stats::qnorm(exp(-x.FROC)) - a.truth )
if( missing(upper_y)){ upper_y <- 1.0 }
suppressWarnings(graphics::par(new=TRUE));
plot(x.FROC,y.FROC.truth,
col = 'red',
cex= 0.1 ,
xlim = c(0,upper_x ),ylim = c(0,upper_y),
xlab = xlabel,
ylab = 'cumulative hit per lesion'
,main =""
);
}#Def of function
Try the BayesianFROC package in your browser
Any scripts or data that you put into this service are public.
BayesianFROC documentation built on Jan. 23, 2022, 9:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions validation.draw_srsc error_srsc_variance_visualization error_srsc_error_visualization error_srsc validation.dataset_srsc
Documented in error_srsc error_srsc_error_visualization error_srsc_variance_visualization validation.dataset_srsc validation.draw_srsc
#' @title Errors of Estimator for any Given true parameter
#@title ----
#'@description
#'
#' This function replicates many datasets from a model with a given model parameter as truth.
#' Then fit a model to these datasets and get estimates.
#' Then calculates mean or variance of the difference of estimates and truth.
#'
#'
#' This function gives a validation under the assumption that
#' we obtain a fixed true model parameter drawn from the prior.
#'
#'But, you know, the author noticed that it is not sufficient because
#' , in Bayesian, such a true parameter is merely one time sample from
#' prior. So, what we should do is to draw randomly many samples from priors.
#'
#' Well, I know, but, I am being such a couch potato. In the future, I would do,
#' if I was alive.
#' Even if I try this, this effort never take account by someone.
#' But, to live, it is required. No money, no life. 2020 NOv 29
#'
#'
#'
#' In the future, what I should do is that the following calculations.
#'
#'
#' 1. Draw a sample of parameter \eqn{\theta} from prior, namely, \eqn{\theta \sim \pi(\theta)}
#'
#' 2. Draw a sample of data \eqn{x=x(\theta)} from a model with the sample of prior \eqn{ x \sim \pi(x|\theta)}
#'
#' 3. Draw a sample of parameter \eqn{\theta'=\theta'(x(_\theta))} from a posterior with the sample of data \eqn{ \theta' \sim \pi(\theta|x)}
#'
#' 4. Calculates the integral \eqn{ \int | \theta' - \theta|^2 \pi(\theta)d\theta} as the error of estimates among all priors..
#'
#'
#'
#'
#'
#'
#'
#' @details
#'
#'Let us denote a model parameter by \eqn{\theta_0}
#' \eqn{N_I} by a number of
#' images and number
#' of lesions by \eqn{N_L} which are specified
#' by user as the variables of the function.
#'
#' \describe{
#' \item{ \strong{(I)} Replicates models for datasets \eqn{D_1,D_2,...,D_k,...,D_K}. }{ }
#' \item{ Draw a dataset \eqn{D_k} from a likelihood (model), namely }{ \eqn{D_k \sim likelihood(|\theta_0)}. }
#' \item{ Draw a MCMC samples \eqn{\{ \theta_i (D_k)\}} from a posterior, namely }{ \eqn{ \theta _i \sim \pi(|D_k)}. }
#' \item{ Calculate a posterior mean, namely }{ \eqn{ \bar{\theta}(D_k) := \sum_i \theta_i(D_k) }. }
#' \item{ Calculates error for \eqn{D_k} }{ \eqn{\epsilon_k}:=Trueth - posterior mean estimates of \eqn{D_k} = \eqn{|\theta_0 - \bar{\theta}(D_k)|} (or = \eqn{\theta_0 - \bar{\theta}(D_k)}, accordinly by the user specified \code{absolute.errors} ). }
#' \item{ \strong{(II)} Calculates mean of errors }{ mean of errors \eqn{ \bar{\epsilon}(\theta_0,N_I,N_L)}= \eqn{ \frac{1}{K} \sum \epsilon_k } }
#' }
#'
#' Running this function, we can see that the error \eqn{ \bar{\epsilon}(\theta_0,N_I,N_L)} decreases
#' monotonically as a given number of images \eqn{N_I} or a given number of lesions \eqn{N_L} increases.
#'
#' Also, the scale of error also will be found. Thus this function can show how our estimates are correct.
#' Scale of error differs for each componenet of model parameters.
#'
#'
#'
#' Revised 2019 August 28
#'
#'
#'
#'
#'
#'
#'
#'@return Return values is,
#'
#' \describe{
#' \item{ Stanfit objects }{ for each Replicated datasets }
#' \item{ Errors }{ EAPs minus true values, in the above notations, it is \eqn{ \bar{\epsilon}(\theta_0,N_I,N_L)} }
#' \item{ Variances of estimators. }{ This calculates the vaiance of posterior means over all replicated datasets }
#' }
#'
#'@param mean.truth This is a parameter
#'of the latent Gaussian assumption
#'for the noise distribution.
#'@param sd.truth This is a parameter of the latent
#' Gaussian assumption for the noise distribution.
#'@param z.truth This is a parameter of the
#' latent Gaussian assumption for the noise distribution.
# @param seed_sampling to be passed to rstan::sampling()
#'@param NL Number of Lesions.
#'@param NI Number of Images.
#'@param replicate.datset A Number indicate
#'that how many you replicate dataset
#' from user's specified dataset.
#'@param serial.number A positive integer
#'or Character. This is for programming perspective.
#' The author use this to print the serial
#' numbre of validation. This will be used
#' in the validation function.
#'@param base_size An numeric for size of object, this is for the package developer.
#'@param absolute.errors A logical specifying whether mean of errors is defined by
#'
#' \describe{
#' \item{ \code{TRUE} }{ \eqn{ \bar{\epsilon}(\theta_0,N_I,N_L)}= \eqn{ \frac{1}{K} \sum | \epsilon_k | } }
#' \item{ \code{FALSE} }{ \eqn{ \bar{\epsilon}(\theta_0,N_I,N_L)}= \eqn{ \frac{1}{K} \sum \epsilon_k } }
#' }
#'
#'
#'
#'@inheritParams fit_Bayesian_FROC
#'@inheritParams DrawCurves
#'@examples
#' \dontrun{
# ####1#### ####2#### ####3#### ####4#### ####5#### ####6#### ####7#### ####8#### ####9####
#'#=========================== The first example ======================================
#'
#'
#'# It is sufficient to run the function with default variable
#'
#' datasets <- validation.dataset_srsc()
#'
#'
#'# Today 2020 Nov 29 I have completely forgotten this function, oh it helps me. Thank me.
# ####1#### ####2#### ####3#### ####4#### ####5#### ####6#### ####7#### ####8#### ####9####
#'#============================= The second example ======================================
#'
#'# If user does not familiar with the values of thresholds, then
#'# it would be better to use the actual estimated values
#'# as an example of true parameters. In the following,
#'# I explain this.
#'# First, to get posterior mean estimates, we run the following:
#'
#'
#'
#' fit <- fit_Bayesian_FROC(dataList.Chakra.1,ite = 1111,summary =FALSE,cha=3)
#'
#'
#'
#'
#'
#'
#'# Secondly, extract the expected a posterior estimates (EAPs) from the object fit
#'# Note that EAP is also called "posterior mean"
#'
#' z <- rstan::get_posterior_mean(fit,par=c("z"))[,"mean-all chains"]
#'
#'
#'
#'
#'
#'# Thirdly we use this z as a true value.
#'
#'
#' datasets <- validation.dataset_srsc(z.truth = z)
#'
#'
#'
#'#========================================================================================
#'# 1) extract replicated fitted model object
#'#========================================================================================
#'
#'
#' # Replicates models
#'
#' a <- validation.dataset_srsc(replicate.datset = 3,ite = 111)
#'
#'
#'
#' # Check convergence, in the above MCMC iterations = 111 which is too small to get
#' # a convergence MCMC chain, and thus the following example will the example
#' # of a non-convergent model in the r hat criteria.
#'
#' ConfirmConvergence( a$fit[[3]])
#'
#'
#' # Check trace plot to confirm whether MCMC chain converge or not.
#'
#' stan_trace( a$fit[[3]],pars = "A")
#'
#'
#' # Check p value, for chi square goodness of fit whose null hypothesis is that
#' # the model is fitted well.
#'
#' fit@posterior_predictive_pvalue_for_chi_square_goodness_of_fit
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#' # Revised in 2019 August 29
#'
#' # Revised in 2020 Nov 28
#' # It is weird, awesome,
#' # What a fucking English,...I fix it.
#'
#'
#'
#'
#'#========================================================================================
#'# 1) Histogram of error of postrior means for replicated datasets
#'#========================================================================================
#'#'
#'
#' a<- validation.dataset_srsc(replicate.datset = 100)
#' hist(a$error.of.AUC,breaks = 111)
#' hist(a$error.of.AUC,breaks = 30)
#'
#'
#'
#'
#'
#'
# ####1#### ####2#### ####3#### ####4#### ####5#### ####6#### ####7#### ####8#### ####9####
#'#========================================================================================
#'# absolute.errors = FALSE generates negative biases
#'#========================================================================================
#'
#'
#' validation.dataset_srsc(absolute.errors = FALSE)
#'
#'
#'
# ####1#### ####2#### ####3#### ####4#### ####5#### ####6#### ####7#### ####8#### ####9####
#'#========================================================================================
#'# absolute.errors = TRUE coerce negative biases to positives, i.e., L^2 norm
#'#========================================================================================
#'
#'
#' validation.dataset_srsc(absolute.errors = TRUE)
#'
#'
#'
#'
#'
#'
#'
#'
# ####1#### ####2#### ####3#### ####4#### ####5#### ####6#### ####7#### ####8#### ####9####
#'#========================================================================================
#'# Check each fitted model object
#'#========================================================================================
#'
#'
#'
#'
#' a <- validation.dataset_srsc(verbose = TRUE)
#'
#' a$fit[[2]]
#'
#'
#'
#'class(a$fit[[2]])
#'
#' rstan::traceplot(a$fit[[2]], pars = c("A"))
#'
#'
#'
#'
#'
#'
# ####1#### ####2#### ####3#### ####4#### ####5#### ####6#### ####7#### ####8#### ####9####
#'#========================================================================================
#'# NaN ... why? 2021 Dec
#'#========================================================================================
#'
#' fits <- validation.dataset_srsc()
#'
#' f <-fits$fit[[1]]
#' rstan::extract(f)$dl
#' sum(rstan::extract(f)$dl)
#' Is.nan.in.MCMCsamples <- as.logical(!prod(!is.nan(rstan::extract(f)$dl)))
#' rstan::extract(f)$A[525]
#' a<-rstan::extract(f)$a[525]
#' b<-rstan::extract(f)$b[525]
#'
#' Phi( a/sqrt(b^2+1) )
#' x<-rstan::extract(f)$dl[2]
#'
#' a<-rstan::extract(f)$a
#' b<-rstan::extract(f)$b
#'
#' a/(b^2+1)
#' Phi(a/(b^2+1))
#' mean( Phi(a/(b^2+1)) )
#'
#' #'
#'
#'
#'
#'
#'}# dontrun
#' @export
validation.dataset_srsc <-function(
replicate.datset =3,
ModifiedPoisson = FALSE,
mean.truth=0.6,
sd.truth=5.3,
z.truth =c(-0.8,0.7,2.38),
NL=259,
NI=57,
ite =1111,
cha =1,
summary=TRUE,
see = 1234,
verbose = FALSE,
serial.number=1,
base_size=0,
absolute.errors = TRUE
# , convergence.model.only = FALSE
){
C <- length(z.truth)
c<-C:1
dz.truth <-vector()
for (cd in 1:C-1) {
dz.truth[cd] <- z.truth[cd+1]-z.truth[cd]
}
if(ModifiedPoisson==F){ NX <- NI}
if(ModifiedPoisson==T){ NX <- NL}
p.truth <- vector()
list.of.dataList <- list()
vector.of.dataList <- list()
name.vector.of.dataList <- vector()
for (cd in 1:C) {
name.vector.of.dataList[cd] <- paste("f[",C-cd+1,"]",sep = "")
}
for (cd in 1:C) {
name.vector.of.dataList[C+cd] <- paste("h[",C-cd+1,"]",sep = "")
}
name.vector.of.dataList[2*C+1] <- "NL"
name.vector.of.dataList[2*C+2] <- "NI"
l.truth <- -log( stats::pnorm(z.truth))
for (seed in 1:replicate.datset ) {
f.inv <- vector()#these should in for sentence
h.inv <- vector()
f <- vector()
h <- vector()
for (cd in 1:C) {
if (cd==C) {p.truth[cd]= 1 - stats::pnorm((z.truth[cd]-mean.truth)/sd.truth)
}else{
p.truth[cd]<- stats::pnorm((z.truth[cd+1]-mean.truth)/sd.truth) - stats::pnorm((z.truth[cd]-mean.truth)/sd.truth)
}}
# here is hitrate modified -----
deno <-vector()
hit_rate.truth <- vector()
deno[C-1] = 1-p.truth[C];
for(cd in 3:C){ deno[c[cd]] = deno[c[cd-1]]-p.truth[c[cd-1]]; }
hit_rate.truth[C] = p.truth[C];
for(cd in 1:C-1) hit_rate.truth[cd] = p.truth[cd]/deno[cd];
s <-0# 2019 Sept 30. This is a key idea
c <-C:1
for (cd in 1:C) {
if(ModifiedPoisson==F){
if (cd==C) {
set.seed(seed = seed);
f.inv[cd] <- stats::rpois(n= 1, lambda = (l.truth[cd]-0)*NI )
}else{
set.seed(seed = seed);
f.inv[cd] <- stats::rpois(n= 1, lambda = (l.truth[cd]-l.truth[cd+1])*NI )
}#else
}# if ModifiedPoisson==F
if(ModifiedPoisson==T){
if (cd==C) {
set.seed(seed = seed);
f.inv[cd] <- stats::rpois(n= 1, lambda = (l.truth[cd]-0)*NL )
}else{
set.seed(seed = seed);
f.inv[cd] <- stats::rpois(n= 1, lambda = (l.truth[cd]-l.truth[cd+1])*NL)
}#else set.seed(seed = seed); hits <- stats::rbinom(n=1,size = NL,prob = p.truth[cd])
}# if ModifiedPoisson==T
# set.seed(seed = seed); h.inv[cd] <- stats::rbinom(n=1,
# size = NL-s, # 2019 Sept 30. This is a key idea # In order to avoid the sum of hits may be greater than NL
# prob = p.truth[cd])
# s <- s + h.inv[cd] # 2019 Sept 30. This is a key idea # In order to avoid the sum of hits may be greater than NL
#
#
set.seed(seed = seed); h[cd] <- stats::rbinom(n=1,
size = NL-s, # 2019 Sept 30. This is a key idea # In order to avoid the sum of hits may be greater than NL
prob = hit_rate.truth[c[cd]])
s <- s + h[cd] # 2019 Sept 30. This is a key idea # In order to avoid the sum of hits may be greater than NL
f[C-cd+1] <- f.inv[cd]
# h[C-cd+1] <- h.inv[cd]
}# for cd in 1:C
list.of.dataList[[seed]] <- list(
NL=NL,
NI=NI,
f=f,
h=h,
C=C
)
a<-append( list.of.dataList[[seed]]$f, list.of.dataList[[seed]]$h);
b <-append( list.of.dataList[[seed]]$NL, list.of.dataList[[seed]]$NI);
vector.of.dataList[[seed]] <-append(a,b)
# browser()
names(vector.of.dataList[[seed]] ) <- name.vector.of.dataList
}#for seed
#Here
#
# * Datasets were created !!
#
# * Next Fitting !!
#
#Here
fit <- list()
validation.of.AUC.EAP <-vector()
validation.of.AUC.CI.lower <-vector()
validation.of.AUC.CI.upper <-vector()
validation.of.z.Threshold.EAP <- list()
validation.of.z.Threshold.CI.lower <- list()
validation.of.z.Threshold.CI.upper <- list()
validation.of.dz.Threshold.EAP <- list()
validation.of.dz.Threshold.CI.lower <- list()
validation.of.dz.Threshold.CI.upper <- list()
validation.of.mean.of.latent.EAP <- vector()
validation.of.mean.of.latent.CI.lower <- vector()
validation.of.mean.of.latent.CI.upper <- vector()
validation.of.Standard.Deviation.latent.EAP <- vector()
validation.of.Standard.Deviation.latent.CI.lower <- vector()
validation.of.Standard.Deviation.CI.upper <- vector()
validation.of.p.Hitrate.EAP <- list()
validation.of.p.Hitrate.CI.lower <- list()
validation.of.p.Hitrate.CI.upper <- list()
validation.of.l.FalseRate.EAP <- list()
validation.of.l.FalseRate.CI.lower <- list()
validation.of.l.FalseRate.CI.upper <- list()
validation.of.loglikelihood.of.model <- vector()
convergent.dataList <- list()
convergent.dataList.as.dataframe <- data.frame(row.names = name.vector.of.dataList)
# rownames(convergent.dataList.as.dataframe) <- name.vector.of.dataList
s<-0
# base_size <-0
for (seed in 1:replicate.datset ) {
#It seems better to make a indicator and
#do not print the result of rstan
# fit ------
color_message(seed,"-th fitting STARTS-----------------------------------------------------------------------")
fit[[seed]] <- fit_Bayesian_FROC(
dataList=list.of.dataList[[seed]],
ite = ite,
cha = cha,
see = see,#seed ---------
new.imaging.device = FALSE,
summary = FALSE,
verbose = verbose,
ModifiedPoisson = ModifiedPoisson,
DrawCurve = FALSE
)
# };browser(); for (seed in 1:replicate.datset ) {
print(stanfit_from_its_inherited_class( fit[[seed]] ) )
base_size <- base_size + size_of_return_value(fit[[seed]] ,is_return_value = FALSE,summary = FALSE)
# here col of size ----
print(size_of_return_value(fit[[seed]] ,is_return_value = FALSE, base_size =base_size,col=TRUE ))
# view data ----
print(viewdata(list.of.dataList[[seed]]))
message("\n----------------------------------------\n",sep = "")
message("\n* ",serial.number, " -th validation")
message("\n* The (",seed," / ", replicate.datset, " ) -th fitting finished.\n",sep = "")
message("\n----------------------------------------\n",sep = "")
# Here ------
if ( fit[[seed]]@convergence ==TRUE) {
s <-s+1
convergent.dataList[[s]] <-fit[[seed]]@dataList
a<-append( convergent.dataList[[s]]$f, convergent.dataList[[s]]$h);
b <-append( convergent.dataList[[s]]$NL, convergent.dataList[[s]]$NI);
convergent.dataList.as.dataframe[,s] <-append(a,b)
# browser()
#ssss -----
ssss <- summary_EAP_CI_srsc( fit[[seed]],summary = FALSE,dig =11)
validation.of.AUC.EAP[s] <- ssss$AUC.EAP
validation.of.AUC.CI.lower[s] <-ssss$AUC.CI.lower
validation.of.AUC.CI.upper[s] <- ssss$AUC.CI.upper
validation.of.z.Threshold.EAP[[s]] <- ssss$z.Threshold.EAP
validation.of.z.Threshold.CI.lower[[s]] <- ssss$z.Threshold.CI.lower
validation.of.z.Threshold.CI.upper[[s]] <- ssss$z.Threshold.CI.upper
validation.of.dz.Threshold.EAP[[s]] <- ssss$dz.Threshold.EAP
validation.of.dz.Threshold.CI.lower[[s]] <- ssss$dz.Threshold.CI.lower
validation.of.dz.Threshold.CI.upper[[s]] <- ssss$dz.Threshold.CI.upper
validation.of.mean.of.latent.EAP[s] <- ssss$mean.of.latent.EAP
validation.of.mean.of.latent.CI.lower[s] <- ssss$mean.of.latent.CI.lower
validation.of.mean.of.latent.CI.upper[s] <- ssss$mean.of.latent.CI.upper
validation.of.Standard.Deviation.latent.EAP[s] <- ssss$Standard.Deviation.latent.EAP
validation.of.Standard.Deviation.latent.CI.lower[s] <- ssss$Standard.Deviation.latent.CI.lower
validation.of.Standard.Deviation.CI.upper[s] <- ssss$Standard.Deviation.CI.upper
validation.of.p.Hitrate.EAP[[s]] <- ssss$p.Hitrate.EAP
validation.of.p.Hitrate.CI.lower[[s]] <- ssss$p.Hitrate.CI.lower
validation.of.p.Hitrate.CI.upper[[s]] <- ssss$p.Hitrate.CI.upper
validation.of.l.FalseRate.EAP[[s]] <- ssss$l.FalseRate.EAP
validation.of.l.FalseRate.CI.lower[[s]] <- ssss$l.FalseRate.CI.lower
validation.of.l.FalseRate.CI.upper[[s]] <- ssss$l.FalseRate.CI.upper
validation.of.loglikelihood.of.model[s] <- ssss$loglikelihood.of.model
}# if convergence is true
color_message(seed,"-th fitting FINISHED-----------------------------------------------------------------------")
}#for seed
# browser()
col.names <- vector()
for (sd in 1:s) {
col.names[sd] <- paste(sd,"-th model",sep = "")
}
colnames(convergent.dataList.as.dataframe) <- col.names
number.of.convergent.model <-s#Here
message("\n ----- Comments for Validation -----\n")
message("\n* Number of all replicated models:",replicate.datset,"\n")
message("\n* Number of convergent models:",s,"\n")
message("\n* Convergence rate:= convergent models / convergent and non convergent models =", round( (s/replicate.datset)*100 , digits = 2),"% \n")
a.truth <- mean.truth/sd.truth
b.truth <- 1/sd.truth
AUC.truth <- stats::pnorm(a.truth/ sqrt( 1+b.truth^2))
l.truth <-l.truth[1:C]
# here error.of.AUC------
error.of.AUC <- validation.of.AUC.EAP - AUC.truth
if (absolute.errors) error.of.AUC <- abs(validation.of.AUC.EAP - AUC.truth)
l1 <-validation.of.z.Threshold.EAP
# l2 <-rep( list(z.truth) , length(unlist( validation.of.z.Threshold.EAP))/length(z.truth) )
l2 <-rep( list(z.truth) , number.of.convergent.model)
# here error.of.z.Threshold.EAP------
# browser()
error.of.z.Threshold.EAP <- Map( `-`,l1 , l2 )
if (absolute.errors) error.of.z.Threshold.EAP <- lapply( Map( `-`,l1 , l2 ),abs)
mean.error.of.z.Threshold.EAP<- rep(NaN, length(error.of.z.Threshold.EAP[[1]]))
sd.error.of.z.Threshold.EAP<- rep(NaN, length(error.of.z.Threshold.EAP[[1]]))
name.z.Threshold.EAP<- rep(NaN, length(error.of.z.Threshold.EAP[[1]]))
for(i.th.column in 1: length(error.of.z.Threshold.EAP[[1]]) ){
mean.error.of.z.Threshold.EAP[i.th.column]<-mean(as.numeric(lapply(error.of.z.Threshold.EAP,"[",n=i.th.column)))
# sd.error.of.z.Threshold.EAP[i.th.column]<-stats::sd(as.numeric(lapply(error.of.z.Threshold.EAP,"[",n=i.th.column)))
sd.error.of.z.Threshold.EAP[i.th.column]<-stats::sd(as.numeric(lapply(validation.of.z.Threshold.EAP,"[",n=i.th.column)))
name.z.Threshold.EAP[i.th.column]<- paste("z[", i.th.column, "]")
}
# browser()
l1 <-validation.of.dz.Threshold.EAP
# l2 <-rep( list(z.truth) , length(unlist( validation.of.z.Threshold.EAP))/length(z.truth) )
l2 <-rep( list(dz.truth) , number.of.convergent.model)
# error ----
# browser()
error.of.dz.Threshold.EAP <- Map( `-`,l1 , l2 )
if (absolute.errors) error.of.dz.Threshold.EAP <- lapply( Map( `-`,l1 , l2 ),abs)
mean.error.of.dz.Threshold.EAP<- rep(NaN, length(error.of.dz.Threshold.EAP[[1]]))
sd.error.of.dz.Threshold.EAP<- rep(NaN, length(error.of.dz.Threshold.EAP[[1]]))
name.dz.Threshold.EAP<- rep(NaN, length(error.of.dz.Threshold.EAP[[1]]))
for(i.th.column in 1: length(error.of.dz.Threshold.EAP[[1]]) ){
mean.error.of.dz.Threshold.EAP[i.th.column]<-mean(as.numeric(lapply(error.of.dz.Threshold.EAP,"[",n=i.th.column)))
# sd.error.of.z.Threshold.EAP[i.th.column]<-stats::sd(as.numeric(lapply(error.of.z.Threshold.EAP,"[",n=i.th.column)))
sd.error.of.dz.Threshold.EAP[i.th.column]<-stats::sd(as.numeric(lapply(validation.of.dz.Threshold.EAP,"[",n=i.th.column)))
name.dz.Threshold.EAP[i.th.column]<- paste("dz[", i.th.column, "]")
}
# error ----
# browser()
if (absolute.errors) error.of.mean.of.latent.EAP <- abs( validation.of.mean.of.latent.EAP - mean.truth)
error.of.mean.of.latent.EAP <- validation.of.mean.of.latent.EAP - mean.truth
if (absolute.errors) error.of.Standard.Deviation.latent.EAP <-abs( validation.of.Standard.Deviation.latent.EAP -sd.truth)
error.of.Standard.Deviation.latent.EAP <- validation.of.Standard.Deviation.latent.EAP -sd.truth
#S.E.
sd.error.of.mean.of.latent.EAP <- stats::sd(validation.of.mean.of.latent.EAP)
sd.error.of.Standard.Deviation.latent.EAP <-stats::sd(validation.of.mean.of.latent.EAP)
sd.error.of.AUC <-stats::sd(validation.of.AUC.EAP)
# sd.error.of.AUC <-stats::sd(error.of.AUC)
name.Gaussian <- append(name.z.Threshold.EAP, c("mean.of.signal", "sd.of.signal","AUC") )
Bias.Gaussian <- append(mean.error.of.z.Threshold.EAP, c(mean(error.of.mean.of.latent.EAP),mean(error.of.Standard.Deviation.latent.EAP),mean(error.of.AUC ) ) )
Std.Err.Gaussian <- append(sd.error.of.z.Threshold.EAP, c(sd.error.of.mean.of.latent.EAP, sd.error.of.Standard.Deviation.latent.EAP, sd.error.of.AUC) )
name.Gaussian <- append(name.Gaussian, name.dz.Threshold.EAP )
Bias.Gaussian <- append(Bias.Gaussian, mean.error.of.dz.Threshold.EAP )
Std.Err.Gaussian <- append(Std.Err.Gaussian,sd.error.of.dz.Threshold.EAP )
if (summary==TRUE) {
message("\n* This is a model parameter for the latent Gaussian distribution.\n")
print(knitr::kable( data.frame(Param.Name=name.Gaussian, ############ Here !!
Bias = Bias.Gaussian, ########### Here !!
S.E. = Std.Err.Gaussian )########### Here !!
))
}
l1 <-validation.of.p.Hitrate.EAP
l2 <-rep( list(p.truth) , length(unlist( validation.of.p.Hitrate.EAP))/length(p.truth) )
error.of.p.Hitrate.EAP <- Map( `-`,l1 , l2 )
if (absolute.errors) error.of.p.Hitrate.EAP <- lapply( Map( `-`,l1 , l2 ),abs)
mean.error.p.Hitrate.EAP<- rep(NaN, length(error.of.p.Hitrate.EAP[[1]]))
sd.error.p.Hitrate.EAP<- rep(NaN, length(error.of.p.Hitrate.EAP[[1]]))
name.p.Hitrate.EAP<- rep(NaN, length(error.of.p.Hitrate.EAP[[1]]))
for(i.th.column in 1: length(error.of.p.Hitrate.EAP[[1]]) ){
# mean ----
mean.error.p.Hitrate.EAP[i.th.column] <- mean(as.numeric(lapply(error.of.p.Hitrate.EAP,"[",n=i.th.column)))
sd.error.p.Hitrate.EAP[i.th.column] <- stats::sd(as.numeric(lapply(error.of.p.Hitrate.EAP,"[",n=i.th.column)))
name.p.Hitrate.EAP[i.th.column] <- paste("p[", i.th.column, "]")
}
if (summary==TRUE) {
print(knitr::kable( data.frame(Param.Name=name.p.Hitrate.EAP,########### Here !!
Bias = mean.error.p.Hitrate.EAP,########### Here !!
S.E. = sd.error.p.Hitrate.EAP )########### Here !!
)
)
}
# browser()
l1 <-validation.of.l.FalseRate.EAP
l2 <-rep( list(l.truth) , length(unlist( validation.of.l.FalseRate.EAP))/length(l.truth) )
# error ----
error.of.l.FalseRate.EAP <- Map( `-`,l1 , l2 )
if (absolute.errors) error.of.l.FalseRate.EAP <- lapply( Map( `-`,l1 , l2 ),abs)
# mean ----
mean.error.l.FalseRate.EAP<- rep(NaN, length(error.of.l.FalseRate.EAP[[1]]))
sd.error.l.FalseRate.EAP<- rep(NaN, length(error.of.l.FalseRate.EAP[[1]]))
name.l.FalseRate.EAP<- rep(NaN, length(error.of.l.FalseRate.EAP[[1]]))
for(i.th.column in 1: length(error.of.l.FalseRate.EAP[[1]]) ){
# mean ----
mean.error.l.FalseRate.EAP[i.th.column]<-mean(as.numeric(lapply(error.of.p.Hitrate.EAP,"[",n=i.th.column)))
sd.error.l.FalseRate.EAP[i.th.column]<-stats::sd(as.numeric(lapply(error.of.p.Hitrate.EAP,"[",n=i.th.column)))
name.l.FalseRate.EAP[i.th.column]<- paste("lambda[", i.th.column, "]")
}
if (summary==TRUE) {
print(knitr::kable( data.frame(Param.Name=name.l.FalseRate.EAP,########### Here !!
Bias = mean.error.l.FalseRate.EAP,########### Here !!
S.E. = sd.error.l.FalseRate.EAP )) )########### Here !!
# browser()
}
Name.of.all.param <- append(name.Gaussian,
append(name.p.Hitrate.EAP,
name.l.FalseRate.EAP)
)
# error ----
Bias.of.all.param <- append(
Bias.Gaussian,
append(
mean.error.p.Hitrate.EAP,
mean.error.l.FalseRate.EAP)
)
S.E.of.all.param <- append(Std.Err.Gaussian,
append(
sd.error.p.Hitrate.EAP,
sd.error.l.FalseRate.EAP)
)
results <- data.frame( Name.of.all.param=Name.of.all.param,
Bias.of.all.param=Bias.of.all.param,
S.E.of.all.param=S.E.of.all.param
)
# error ----
error.of.mean.of.latent.EAP <- validation.of.mean.of.latent.EAP - mean.truth
error.of.Standard.Deviation.latent.EAP <- validation.of.Standard.Deviation.latent.EAP -sd.truth
return_value<- list(
results=results,
Name.of.all.param = Name.of.all.param,
Bias.of.all.param = Bias.of.all.param,
S.E.of.all.param = S.E.of.all.param,
vector.of.dataList=vector.of.dataList,
convergent.dataList = convergent.dataList,
convergent.dataList.as.dataframe = convergent.dataList.as.dataframe,
############################################################# Convergence number and all (non conv and conv) number
number.of.convergent.model =number.of.convergent.model,
replicate.datset=replicate.datset,
#Here
fit=fit,
ModifiedPoisson=ModifiedPoisson,
mean.truth=mean.truth,
sd.truth=sd.truth,
p.truth=p.truth,
l.truth=l.truth,
z.truth=z.truth,
AUC.truth=AUC.truth,
validation.of.AUC.EAP =validation.of.AUC.EAP,
validation.of.AUC.CI.lower =validation.of.AUC.CI.lower,
validation.of.AUC.CI.upper =validation.of.AUC.CI.upper,
validation.of.z.Threshold.EAP = validation.of.z.Threshold.EAP,
validation.of.z.Threshold.CI.lower = validation.of.z.Threshold.CI.lower,
validation.of.z.Threshold.CI.upper = validation.of.z.Threshold.CI.upper,
validation.of.mean.of.latent.EAP = validation.of.mean.of.latent.EAP,
validation.of.mean.of.latent.CI.lower = validation.of.mean.of.latent.CI.lower,
validation.of.mean.of.latent.CI.upper = validation.of.mean.of.latent.CI.upper,
validation.of.Standard.Deviation.latent.EAP = validation.of.Standard.Deviation.latent.EAP,
validation.of.Standard.Deviation.latent.CI.lower = validation.of.Standard.Deviation.latent.CI.lower,
validation.of.Standard.Deviation.CI.upper = validation.of.Standard.Deviation.CI.upper,
validation.of.p.Hitrate.EAP = validation.of.p.Hitrate.EAP,
validation.of.p.Hitrate.CI.lower = validation.of.p.Hitrate.CI.lower,
validation.of.p.Hitrate.CI.upper = validation.of.p.Hitrate.CI.upper,
validation.of.l.FalseRate.EAP = validation.of.l.FalseRate.EAP,
validation.of.l.FalseRate.CI.lower = validation.of.l.FalseRate.CI.lower,
validation.of.l.FalseRate.CI.upper = validation.of.l.FalseRate.CI.upper,
validation.of.loglikelihood.of.model = validation.of.loglikelihood.of.model,
# here errors ------
error.of.AUC=error.of.AUC,
error.of.z.Threshold.EAP = error.of.z.Threshold.EAP,
error.of.mean.of.latent.EAP =error.of.mean.of.latent.EAP,
error.of.Standard.Deviation.latent.EAP =error.of.Standard.Deviation.latent.EAP,
error.of.p.Hitrate.EAP =error.of.p.Hitrate.EAP,
error.of.l.FalseRate.EAP =error.of.l.FalseRate.EAP,
dataList=list.of.dataList
)
invisible( return_value )
}#function
#'@title Validation via replicated datasets from a model at a given model parameter
#'@description
#' Print for a given true parameter,
#' a errors of estimates from replicated dataset.
#'
#' Also print a standard error which is the variance of estimates.
#'
#'
#'
#' Suppose that \eqn{\theta_0} is a given true model parameter with a given number of images \eqn{N_I} and a given number of lesions \eqn{N_L}, specified by user.
#' \describe{
#' \item{ \strong{(I)} }{}
#' \item{ \strong{(I.1)} Synthesize a collection of dataset \eqn{D_k} (\eqn{k=1,2,...,K}) from a likelihood (model) at a given parameter \eqn{\theta_0}, namely }{ \eqn{D_k \sim} likelihood( \eqn{\theta_0}). }
#' \item{ \strong{(I.2)} Replicates \eqn{K} models fitted to each dataset \eqn{D_k} (\eqn{k=1,2,...,K}), namely, draw MCMC samples \eqn{\{ \theta_i (D_k);i=1,...,I\}} from each posterior of the dataset \eqn{D_k}, namely }{ \eqn{ \theta _i(D_k) \sim \pi(|D_k)}. }
#' \item{ \strong{(I.3)} Calculate posterior means for the set of data \eqn{D_k} (\eqn{k=1,2,...,K}), namely }{ \eqn{ \bar{\theta}(D_k) := \frac{1}{I} \sum_i \theta_i(D_k) }. }
#' \item{ \strong{(I.4)} Calculates error for each dataset \eqn{D_k} }{ \eqn{\epsilon_k}:= Truth - estimates = \eqn{\theta_0 - \bar{\theta}(D_k)}. }
#' \item{ \strong{(II)} Calculates mean of errors over all datasets \eqn{D_k} (\eqn{k=1,2,...,K}) }{ mean of errors \eqn{ \bar{\epsilon}(\theta_0,N_I,N_L)}= \eqn{ \frac{1}{K} \sum \epsilon_k }. }
#' \item{ NOTE }{ We note that if a fitted model does not converge,( namely R hat is far from one), then it is omiited from this calculation.}
#' \item{ \strong{(III) } Calculates mean of errors for various number of lesions and images }{ mean of errors \eqn{ \bar{\epsilon}(\theta_0,N_I,N_L)} }
#' }
#'
#' For example, if \eqn{(N_I^1,N_L^1),(N_I^2,N_L^2),(N_I^3,N_L^3),...,(N_I^m,N_L^m)}, then
#' \eqn{ \bar{\epsilon}(\theta_0,N_I^1,N_L^1)},
#' \eqn{ \bar{\epsilon}(\theta_0,N_I^2,N_L^2)},
#' \eqn{ \bar{\epsilon}(\theta_0,N_I^3,N_L^3)},...,
#' \eqn{ \bar{\epsilon}(\theta_0,N_I^m,N_L^m)} are calculated.
#'
#' To obtain precise error,
#' The number of replicated fitted models (denoted by \eqn{K}) should be large enough.
#' If \eqn{K} is small, then it causes a bias.
#' \eqn{K} = \code{replicate.datset}: a variable of the function \code{error_srsc}.
#'
#'
#'
#'
#'
#'
#'
#'
#describe
#'
#'
#' Running this function, we can see that the error
#' \eqn{ \bar{\epsilon}(\theta_0,N_I,N_L)} decreases
#' monotonically as a given number of images \eqn{N_I}
#' or a given number of lesions \eqn{N_L} increases.
#'
#' Also, the scale of error also will be found. Thus this function can show how our estimates are correct.
#' Scale of error differs for each componenet of model parameters.
#'
#'
#'
#' Revised 2019 August 28
#'
#'
#'
#' @details
#'
#'
#' In Bayesian inference,
#' if sample size is large,
#' then posterior tends to the Dirac measure.
#' So, the error and variance of estimates
#' should be tends to zero as sample size tends to infinity.
#'
#' This function check this phenomenen.
#'
#' If model has problem, then it contains some non-decreasing vias
#' with respect to sample size.
#'
#' Revised 2019 Nov 1
#'
#'
#'
#' Provides a reliability of our posterior mean estimates.
#' Using this function, we can
#' find what digit makes sence.
#'
#' In the real world, the data for modality
#' comparison or observer performan evaluation is
#' 100 images or 200 images. In such scale data,
#' any estimate of AUC will contain error at most 0.0113....
#' So, the value of AUC should round in 0.XXX
#' and not 0.XXXX or 0.XXXXX or more. Since
#' error is 0.00113... and hence 4 digit
#' or more digit is meaningless. In such manner,
#' we can analyize the errors.
#'
#'We note that if we increase the
#'number of images or lesinons,
#' the errors decrease.
#'
#' For example, if we use 20000 images in FROC trial,
#' then the error of AUC will be 0.0005... and thus,
#' and so on. Thus large number of images gives
#' us more reliable AUC. However the radiologist
#' cannot read such large (20000) images.
#'
#' Thus, the error will be 0.00113...
#'
#'
#' If the number of images are given before hand and moreover if we obtains the estimates,
#' then we can run this function using these two, we can find the estimated errors by simulation.
#' Of course, the esimates is not the truth, but roughly speaking, if we assume that
#' the estimates is not so far from truth, and the error analysis is rigid with respect to
#' changing the truth, then we can say using estimates as truth, the result of this error analysis can be regarded as an
#' actual error.
#'
#'
#' I want to go home. Unfortunatly, my house is ...
#'
#'
#'
#'@return Replicated datasets, estimates,
#' errors,...etc I made this program 1 years ago?
#' and now I forget ... the precise return values.
#'When I see today, 2019 August. It retains too many return
#'values to explain all of them.
#@param ----
#'@param NLvector A vector of positive integers,
#' indicating a collection of numbers of Lesions.
#@param NIvector Number of Images.
#'@param ratio A positive \strong{\emph{rational}} number,
#' with which Number of Images is determined by the formula:
#' (number of images) = \code{ratio} times (numbser of lesions).
#' Note that in calculation, it rounds \code{ratio * NLvector } to an integer.
#'@inheritParams fit_Bayesian_FROC
#'@inheritParams DrawCurves
#'@inheritParams validation.dataset_srsc
#'@examples
#' \dontrun{
#'#========================================================================================
#'# 0) 0-th example
#'#========================================================================================
#'
#'
#' datasets <-error_srsc(
#' NLvector = c(100,10000,1000000),
#' ite = 2222
#' )
#'
#'
#' # By the following, we can extract only datasets whose
#' # model has converged.
#' datasets$convergent.dataList.as.dataframe
#'
#'
#'
#'
#'#========================================================================================
#'# 1) 1-st example
#'#========================================================================================
#'# Long width is required in R console.
#'
#'
#'
#' datasets <-error_srsc(NLvector = c(
#' 50L,
#' 111L,
#' 11111L
#' ),
#' # NIvector,
#' ratio=2,
#' replicate.datset =3,
#' ModifiedPoisson = FALSE,
#' mean.truth=0.6,
#' sd.truth=5.3,
#' z.truth =c(-0.8,0.7,2.38),
#' ite =2222
#' )
#'
#'
#'
#'#========================================================================================
#'# 2) Plot the error of AUC with respect to NI
#'#========================================================================================
#'
#'
#' a <-error_srsc(NLvector = c(
#' 33L,
#' 50L,
#' 111L,
#' 11111L
#' ),
#' # NIvector,
#' ratio=2,
#' replicate.datset =3,
#' ModifiedPoisson = FALSE,
#' mean.truth=0.6,
#' sd.truth=5.3,
#' z.truth =c(-0.8,0.7,2.38),
#' ite =2222
#' )
#'
#'
#'
#'
#'
#'
#' aa <- a$Bias.for.various.NL
#'
#'
#' error.of.AUC <- aa[8,]
#' y <- subset(aa[8,], select = 2:length(aa[8,]))
#' y <- as.numeric(y)
#' y <- abs(y)
#' upper_y <- max(y)
#' lower_y <- min(y)
#'
#' x <- 1:length(y)
#'
#'
#'
#' plot(x,y, ylim=c(lower_y, upper_y))
#'
#'# From this plot, we cannot see whether the error has decreased or not.
#'# Thus, we replot with the log y-axis, the we will see that the error
#'# has decreased with respect to number of images and lesions.
#'
#' ggplot(data.frame(x=x,y=y), aes(x = x, y = y)) +
#' geom_line() +
#' geom_point() +
#' scale_y_log10()
#'
#'# Revised 2019 Sept 25
#'
#'
#' # General print of log scale
#' df<-data.frame(x=c(10,100,1000,10,100,1000),
#' y=c(1100,220000,33000000,1300,240000,36000000),
#' group=c("1","1","1","2","2","2")
#' )
#'
#' ggplot2::ggplot(df, aes(x = x, y = y, shape = group)) +
#' ggplot2::geom_line(position = position_dodge(0.2)) + # Dodge lines by 0.2
#' ggplot2::geom_point(position = position_dodge(0.2), size = 4)+ # Dodge points by 0.2
#' ggplot2::scale_y_log10()+
#' ggplot2::scale_x_log10()
#'
#'
#'#========================================================================================
#'# 2) Add other param into plot plain of the error of AUC with respect to NI
#'#========================================================================================
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#' a <-error_srsc(NLvector = c(
#' 111L,
#' 11111L
#' ),
#' # NIvector,
#' ratio=2,
#' replicate.datset =3,
#' ModifiedPoisson = FALSE,
#' mean.truth=0.6,
#' sd.truth=5.3,
#' z.truth =c(-0.8,0.7,2.38),
#' ite =2222
#' )
#' aa <- a$Bias.for.various.NL
#'
#'
#' error.of.AUC <- aa[8,]
#' y1 <- subset(aa[8,], select = 2:length(aa[8,]))
#' y1 <- as.numeric(y1)
#' y1 <- abs(y1)
#'
#' LLL <-length(y1)
#'
#' y2 <- subset(aa[7,], select = 2:length(aa[7,]))
#' y2 <- as.numeric(y2)
#' y2 <- abs(y2)
#'
#' y <- c(y1,y2)
#'
#'
#' upper_y <- max(y)
#' lower_y <- min(y)
#'
#' group <- rep(seq(1,2,1),1 , each=LLL)
#' x <- rep(seq(1,LLL,1),2 , each=1)
#' group <- as.character(group)
#' df <- data.frame(x=x,y=y,group=group)
#'
#'
#' ggplot2::ggplot(df, aes(x = x, y = y, shape = group)) +
#' ggplot2::geom_line(position = position_dodge(0.2)) + # Dodge lines by 0.2
#' ggplot2::geom_point(position = position_dodge(0.2), size = 4)+ # Dodge points by 0.2
#' ggplot2::scale_y_log10()
#' # ggplot2::scale_x_log10()
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'#========================================================================================
#'# Confidence level = 4
#'#========================================================================================
#'
#'
#'
#'
#'
#' datasets <-error_srsc(NLvector = c(
#' 111L,
#' 11111L
#' ),
#' # NIvector,
#' ratio=2,
#' replicate.datset =3,
#' ModifiedPoisson = FALSE,
#' mean.truth=-0.22,
#' sd.truth=5.72,
#' z.truth =c(-0.46,-0.20,0.30,1.16),
#' ite =2222
#' )
#'
#'
#'
#'
#'
#' error_srsc_variance_visualization(datasets)
#'
#'# The parameter of model is 7 in which the ggplot2 fails with the following warning:
#'
#'
#'# The shape palette can deal with a maximum of 6 discrete values because more than 6
#'# becomes difficult to
#'# discriminate; you have 7. Consider specifying shapes manually if you must have them.
#'
#'
# ####1#### ####2#### ####3#### ####4#### ####5#### ####6#### ####7#### ####8#### ####9####
#'#========================================================================================
#'# NaN ... why? 2021 Dec
#'#========================================================================================
#'
#' fits <- validation.dataset_srsc()
#'
#' f <-fits$fit[[2]]
#' rstan::extract(f)$dl
#' sum(rstan::extract(f)$dl)
#' Is.nan.in.MCMCsamples <- as.logical(!prod(!is.nan(rstan::extract(f)$dl)))
#' rstan::extract(f)$A[525]
#' a<-rstan::extract(f)$a
#' b<-rstan::extract(f)$b
#'
#'
#'
#' Phi( a[525]/sqrt(b[525]^2+1) )
#' a[525]/sqrt(b[525]^2+1)
#' Phi( a/sqrt(b^2+1) )
#' x<-rstan::extract(f)$dl[2]
#'
#' a<-rstan::extract(f)$a
#' b<-rstan::extract(f)$b
#'
#' a/(b^2+1)
#' Phi(a/(b^2+1))
#' mean( Phi(a/(b^2+1)) )
#'
#' #'
#'
#'
#'}# dontrun
#' @export
# OUt put can be made by TeX source.
error_srsc <-function(NLvector = c(
100L,
10000L,
1000000L
),
# NIvector,
ratio=2,# NUmber of image is not so important, since when we consider ModifiedPoisson = TRUE
replicate.datset =3,
ModifiedPoisson = FALSE,
mean.truth=0.6,
sd.truth=5.3,
z.truth =c(-0.8,0.7,2.38),
ite =2222,
cha =1,
verbose = FALSE
){
C <- length(z.truth)
d <- list()
Bias.of.all.param <- list()
S.E.of.all.param <- list()
Bias.of.all.param.with.NL.NI <- list()
S.E.of.all.param.with.NL.NI <- list()
Name.of.all.param.with.NL.NI <- vector()
col.names_for_SE <- vector()
col.names_for_SE[1] <-paste("Variance of estimates")
col.names <- vector()
col.names[1] <-paste("Error of estimates")
convergent.dataList.as.dataframe <- list()
number.of.convergent.model<- vector()
number.of.replicate.datset<- vector()
# col.names[1] <-paste("")
s <- 0
base_size <- 0
for (nl in NLvector) {
s<- s+1
ni<- floor(ratio*nl)
#Create the replicating datasets for different number of lesions
# validation.dataset_srsc -----
d[[s]]<- validation.dataset_srsc(
replicate.datset =replicate.datset,
ModifiedPoisson = ModifiedPoisson,
mean.truth=mean.truth,
sd.truth=sd.truth,
z.truth =z.truth,
NL=ni,
NI=nl,
ite =ite,
cha =cha,
base_size = base_size ,
verbose = verbose,
summary=FALSE,
serial.number= paste( "The number of Lesion = ",nl, " and The number of images",ni, ". \n\n* The (" , s, "/", length(NLvector), ")")
)
base_size <- base_size + size_of_return_value(d[[s]], is_return_value = FALSE,base_size = base_size,summary = FALSE)
# Here ------
number.of.convergent.model[s] <- d[[s]]$number.of.convergent.model#Here
number.of.replicate.datset[s] <- d[[s]]$replicate.datset#Here
Bias.of.all.param[[s]] <-d[[s]]$Bias.of.all.param
S.E.of.all.param[[s]] <-d[[s]]$S.E.of.all.param
Bias.of.all.param.with.NL.NI[[s]] <- append(append(append(append(ni,nl),
Bias.of.all.param[[s]]),
number.of.convergent.model[s]),
number.of.replicate.datset[s])
S.E.of.all.param.with.NL.NI[[s]] <- append(append(append(append(ni,nl),# 2019 Sept 6
S.E.of.all.param[[s]]),# 2019 Sept 6
number.of.convergent.model[s]),# 2019 Sept 6
number.of.replicate.datset[s])# 2019 Sept 6
col.names[s+1] <-paste(s,"-th validation",sep = "")
if(s==1) col.names[s+1] <-paste(s,"-st validation",sep = "")
if(s==2) col.names[s+1] <-paste(s,"-nd validation",sep = "")
if(s==3) col.names[s+1] <-paste(s,"-rd validation",sep = "")
col.names_for_SE[s+1] <-paste(s,"-th validation",sep = "")# 2019 Sept 6
if(s==1) col.names_for_SE[s+1] <-paste(s,"-st validation",sep = "")# 2019 Sept 6
if(s==2) col.names_for_SE[s+1] <-paste(s,"-nd validation",sep = "")# 2019 Sept 6
if(s==3) col.names_for_SE[s+1] <-paste(s,"-rd validation",sep = "")# 2019 Sept 6
convergent.dataList.as.dataframe[[s]] <- d[[s]]$convergent.dataList.as.dataframe
}# for nl in NLvector
list.names <- vector()
for (sd in 1:s) {
list.names [sd] <- paste(sd,"-th dataset of convergent model ",sep = "")
}
names(convergent.dataList.as.dataframe) <- list.names
# Name.of.all.param.with.NL.NI<- append(append("Number of Images",
# "Number of Lesions"),
# d[[1]]$Name.of.all.param )
Name.of.all.param.with.NL.NI<- append( append( append(append("Number of Images",
"Number of Lesions"),
d[[1]]$Name.of.all.param ),
"number of convergent models"),#Here
"number of replicated datsets")#Here
# browser()
l <-list(
Name.of.all.param.with.NL.NI = Name.of.all.param.with.NL.NI,
Bias.of.all.param.with.NL.NI = Bias.of.all.param.with.NL.NI
)#list
ll <-list(# 2019 Sept 6
Name.of.all.param.with.NL.NI = Name.of.all.param.with.NL.NI,# 2019 Sept 6
S.E.of.all.param.with.NL.NI = S.E.of.all.param.with.NL.NI# 2019 Sept 6
)#list# 2019 Sept 6
# browser()
options(scipen=100)
#Here, list l is converted, and each component of list become column vector!! This is very good!
Bias.for.various.NL<-as.data.frame(l)
S.E.for.various.NL<-as.data.frame(ll)# 2019 Sept 6
names(Bias.for.various.NL) <- col.names
names(S.E.for.various.NL)<- col.names_for_SE # 2019 Sept 6
message(
"
==================================================================
* Each column is independent.
* Each row means the mean of the difference of [ truth - estimate ]
* If each number is small, then it means that each error is small.
* In my model, the most important quantity is a parameter of AUC,
whose error is the most important.
")
message("\n* Variance of estimates, (Standard Error) \n")# 2019 Sept 6
print(knitr::kable(S.E.for.various.NL) )# 2019 Sept 6
message("\n* Error between estimates and specified true parameter \n")
print(knitr::kable(Bias.for.various.NL) )
return_value <- list(
Bias.for.various.NL =Bias.for.various.NL,
S.E.for.various.NL=S.E.for.various.NL,
Name.of.all.param.with.NL.NI = Name.of.all.param.with.NL.NI,
Bias.of.all.param.with.NL.NI = Bias.of.all.param.with.NL.NI,
S.E.of.all.param.with.NL.NI = S.E.of.all.param.with.NL.NI,
C=C,
length.of.NLvector = length( NLvector),
NLvector=NLvector,
ratio=ratio,# NUmber of image is not so important, since when we consider ModifiedPoisson = TRUE
replicate.datset =replicate.datset,
ModifiedPoisson = ModifiedPoisson,
mean.truth=mean.truth,
sd.truth=sd.truth,
z.truth =z.truth,
Number.of.MCMC.iteration =ite,
number.of.convergent.model =number.of.convergent.model,#Here
number.of.replicate.datset =number.of.replicate.datset, #Here
datasets =d,
convergent.dataList.as.dataframe=convergent.dataList.as.dataframe
)
size_of_return_value(
return_value
)
invisible(
return_value
)#invisible
}
# datasets$datasets[[2]]$dataList
#' @title Visualization for Error of Estimator
#' @description
#' The function plot the graph of errors with respect to sample sizes.
#'
#'
#'\strong{\emph{Error plot}}
#'\describe{
#'\item{ \strong{\emph{x-axis}} }{ Sample sizes }
#'\item{ \strong{\emph{y-axis}} }{ Error for each parameter }
#'}
#'
#' @param return.value.of_error_srsc A return value
#' of the function \code{\link{error_srsc}()}.
#' @param log_scale_x.axis A logical,
#' whether x axis is log scale or not.
#' @return A long format dataframe of
#' error and its parameter name
#' @export
#' @seealso \link{error_srsc_variance_visualization}
#'
#' @examples
#'
#' # General plot
#'
#' df <- data.frame(x=runif(100),y=runif(100),g= as.factor(rep(1:5,10)))
#'
#' ggplot(df, aes(x = x, y = y, shape = g)) +
#' geom_point(size = 3) +
#' scale_shape_manual(values = c(1,2,3,4,5,6,7,8,9))
#'
#'
#'
#'
#'
#'
#'
#'
#'
#' df <- data.frame(x=runif(100),y=runif(100),g= as.factor(rep(1:25,4)))
#'
#' # Use slightly larger points and use custom values for the shape scale
#'
#'
#' ggplot(df, aes(x = x, y = y, shape = g)) +
#' geom_point(size = 3) +
#' scale_shape_manual(values = c(1,2,3,4,5,6,7,8,9,10,
#' 11,12,13,14,15,16,17,18,19,20,21,22,23,24,25))
#'
#' \dontrun{
#' a <- error_srsc()
#'
#' error_srsc_error_visualization(a)
#'
#'
# ####1#### ####2#### ####3#### ####4#### ####5#### ####6#### ####7#### ####8#### ####9####
#'#========================================================================================
#'# In case of C = 4, arbitrary C is available.
#'#========================================================================================
#'
#' a <-error_srsc(NLvector = c(
#' 100,
#' 10000,
#' 1000000
#' ),
#' ratio=2,
#' replicate.datset =2,
#' ModifiedPoisson = FALSE,
#' mean.truth=0.6,
#' sd.truth=5.3,
#' z.truth =c(-0.8,0.7,2.38,3), # Here we use the C=4
#' ite =500
#' )
#'
#' error_srsc_error_visualization(a)
#' error_srsc_variance_visualization(a)
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
# ####1#### ####2#### ####3#### ####4#### ####5#### ####6#### ####7#### ####8#### ####9####
#'#========================================================================================
#'# In case of C = 7, arbitrary C is available.
#'#========================================================================================
#'
#'
#'
#'
#'
#'
#'
#' #'
#' a <-error_srsc(NLvector = c(
#' 100,
#' 10000,
#' 100000
#' ),
#' ratio=2,
#' replicate.datset =2,
#' ModifiedPoisson = FALSE,
#' mean.truth=0.6,
#' sd.truth=5.3,
#' z.truth =c(-0.8,0.7,2.38,3,3.4,3.6,3.8), # Here we use the C=7
#' ite =500
#' )
#'
#' error_srsc_error_visualization(a)
#' error_srsc_variance_visualization(a)
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#' }
#'
error_srsc_error_visualization <- function(return.value.of_error_srsc,
log_scale_x.axis = TRUE
){
a<-return.value.of_error_srsc
aa <- a$Bias.for.various.NL
NLvector <- a$NLvector
length.of.NLvector <- a$length.of.NLvector
C <- a$C
LLL <-length.of.NLvector
y.list <- list()
group.list <-list()
for (cd in 1:(C+3)) {
y.list[[cd]] <- vector()
y.list[[cd]] <- abs(
as.numeric(
subset(aa[cd+2,], select = 2:length(aa[cd+2,]))
)
)
# group.list[[cd]] <- subset(aa[,1], subset = 2:(C+2),select =1 )
}
for (cd in 1:C) {
group.list[[cd]] <- vector()
for (nnn in 1:LLL) {
group.list[[cd]][nnn] <- paste("z[",cd,"]")
}
}
group.list[[C+1]] <- vector()
group.list[[C+2]] <- vector()
group.list[[C+3]] <- vector()
for (nnn in 1:LLL) {
group.list[[C+1]][nnn] <- "mean.of.signal"
group.list[[C+2]][nnn] <- "sd.of.signal"
group.list[[C+3]][nnn] <- "AUC"
}
y <- unlist(y.list)
group <- unlist(group.list)
# upper_y <- max(y)
# lower_y <- min(y)
# group <- rep(seq(1,LLL,1),1 , each=C+2)
# group <- as.character(group)
if(log_scale_x.axis==FALSE) x <- rep(seq(1,LLL,1),C+3 , each=1)
if(log_scale_x.axis==TRUE) x <- rep(NLvector ,C+3 , each=1)
# browser()
# browser()
# browser()
df <- data.frame(x=x,y=y,group=group)
# browser()
if (log_scale_x.axis==TRUE) {
print(
ggplot2::ggplot(df, ggplot2::aes(x = x, y = y, shape = group)) +
ggplot2::geom_line(position = ggplot2::position_dodge(0.2)) + # Dodge lines by 0.2
ggplot2::geom_point(position = ggplot2::position_dodge(0.2), size = 4)+ # Dodge points by 0.2
ggplot2::scale_shape_manual(values = c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25))+ # This is a valiation of dots in the plot, and important in case of the number of confidence level is greater than 4
ggplot2::scale_y_log10()+
ggplot2::scale_x_log10()
)
}
if (log_scale_x.axis==FALSE) {
print( ggplot2::ggplot(df, ggplot2::aes(x = x, y = y, shape = group)) +
ggplot2::geom_line(position = ggplot2::position_dodge(0.2)) + # Dodge lines by 0.2
ggplot2::geom_point(position = ggplot2::position_dodge(0.2), size = 4)+ # Dodge points by 0.2
ggplot2::scale_shape_manual(values = c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25))+ # This is a valiation of dots in the plot, and important in case of the number of confidence level is greater than 4
ggplot2::scale_y_log10()
#+ ggplot2::scale_x_log10()
)
}
# browser()
}
#' @title Visualization Of variance Analysis
#'
#' @param return.value.of_error_srsc A return value
#' of the function \code{\link{error_srsc}()}.
#' @param log_scale_x.axis A logical, whether x axis is log scale.
#' @return A long format dataframe of error and its parameter name
#' @export
#'
#' @examples
#' \dontrun{
#' a <- error_srsc()
#'
#' error_srsc_variance_visualization(a)
#'
#'
#'
#' a <- error_srsc(replicate.datset = 10)
#' error_srsc_variance_visualization(a)
#'
#' }
#'
error_srsc_variance_visualization <- function(return.value.of_error_srsc,
log_scale_x.axis = TRUE
){
a<-return.value.of_error_srsc
aa <- a$S.E.for.various.NL
NLvector <- a$NLvector
length.of.NLvector <- a$length.of.NLvector
C <- a$C
LLL <-length.of.NLvector
y.list <- list()
group.list <-list()
for (cd in 1:(C+3)) {
y.list[[cd]] <- vector()
y.list[[cd]] <- abs(
as.numeric(
subset(aa[cd+2,], select = 2:length(aa[cd+2,]))
)
)
# group.list[[cd]] <- subset(aa[,1], subset = 2:(C+2),select =1 )
}
for (cd in 1:C) {
group.list[[cd]] <- vector()
for (nnn in 1:LLL) {
group.list[[cd]][nnn] <- paste("z[",cd,"]")
}
}
group.list[[C+1]] <- vector()
group.list[[C+2]] <- vector()
group.list[[C+3]] <- vector()
for (nnn in 1:LLL) {
group.list[[C+1]][nnn] <- "mean.of.signal"
group.list[[C+2]][nnn] <- "sd.of.signal"
group.list[[C+3]][nnn] <- "AUC"
}
y <- unlist(y.list)
group <- unlist(group.list)
# upper_y <- max(y)
# lower_y <- min(y)
# group <- rep(seq(1,LLL,1),1 , each=C+2)
# group <- as.character(group)
if(log_scale_x.axis==FALSE) x <- rep(seq(1,LLL,1),C+3 , each=1)
if(log_scale_x.axis==TRUE) x <- rep(NLvector ,C+3 , each=1)
# browser()
# browser()
# browser()
df <- data.frame(x=x,y=y,group=group)
# browser()
if (log_scale_x.axis==TRUE) {
print(
ggplot2::ggplot(df, ggplot2::aes(x = x, y = y, shape = group)) +
ggplot2::geom_line(position = ggplot2::position_dodge(0.2)) + # Dodge lines by 0.2
ggplot2::geom_point(position = ggplot2::position_dodge(0.2), size = 4)+ # Dodge points by 0.2
ggplot2::scale_shape_manual(values = c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25))+ # This is a valiation of dots in the plot, and important in case of the number of confidence level is greater than 4
ggplot2::scale_y_log10()+
ggplot2::scale_x_log10()
)
}
if (log_scale_x.axis==FALSE) {
print( ggplot2::ggplot(df, ggplot2::aes(x = x, y = y, shape = group)) +
ggplot2::geom_line(position = ggplot2::position_dodge(0.2)) + # Dodge lines by 0.2
ggplot2::geom_point(position = ggplot2::position_dodge(0.2), size = 4)+ # Dodge points by 0.2
ggplot2::scale_shape_manual(values = c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25))+ # This is a valiation of dots in the plot, and important in case of the number of confidence level is greater than 4
ggplot2::scale_y_log10()
#+ ggplot2::scale_x_log10()
)
}
# browser()
}
#'@title Draw Curves for validation dataset
#'@description drawing curves.
#'
#'
#' \strong{Red curve} indicates an FROC curve of truth parameter.
#'
#' \strong{Other curves} are drawn using replicated estimates.
#'
#'
#'@inheritParams DrawCurves
#'@inheritParams DrawCurves_MRMC_pairwise
#'@inheritParams fit_Bayesian_FROC
#'@param validation.data This is a return value of
#' the function \code{validation.dataset_srsc}.
#'
#'@return \code{NULL}
#'
#'
#'@examples
#'
#' \dontrun{
#'#--------------------------------------------------------------------------------------
#'# 1) Draw the curve for each replicated dataset
#'#--------------------------------------------------------------------------------------
#'
#' datasets <- validation.dataset_srsc()
#'
#' validation.draw_srsc(datasets)
#'
#'
#'
#'
#'#--------------------------------------------------------------------------------------
#'# 1) Draw the curve for each replicated dataset
#'#--------------------------------------------------------------------------------------
#'
#'
#'
#' datasets <- validation.dataset_srsc(replicate.datset = 5)
#'
#'
#' validation.draw_srsc(datasets)
#'
#'
#'}# dottest
#' @export
validation.draw_srsc <-function(validation.data,
mesh.for.drawing.curve=11111,
upper_y=1,
DrawFROCcurve=TRUE
){
replicate.datset <- as.integer(validation.data$replicate.datset)
s<-vector()
for (seed in 1:replicate.datset ) {if (validation.data$fit[[seed]]@convergence == TRUE) {
s <- append(s, validation.data$fit[[seed]]@metadata$ff )
}}
upper_x <-max(s) ### Nice! I was great 2019 September 1
s<-0
for (seed in 1:replicate.datset ) {
if (validation.data$fit[[seed]]@convergence == TRUE) {
s<-s+1
validation.data$fit[[seed]]@index <-s
DrawCurves(validation.data$fit[[seed]],
new.imaging.device = FALSE,
title = FALSE,
DrawFROCcurve=DrawFROCcurve,
indexCFPCTP = s,
upper_x=upper_x,
upper_y=upper_y,
# cex=0.04
)
}
if (validation.data$fit[[seed]]@convergence == FALSE) {
message("\n* The ",seed,"-th curve are not drawn, since the model did not converge.\n",sep = "")
}
}
mean.truth <- validation.data$mean.truth
sd.truth <- validation.data$sd.truth
a.truth <- mean.truth/sd.truth
b.truth <- 1/sd.truth
if(validation.data$fit[[1]]@studyDesign == "srsc.per.image"){ xlabel = 'mean of cumulative false positives per image' }
if(validation.data$fit[[1]]@studyDesign == "srsc.per.lesion"){ xlabel = 'mean of cumulative false positives per lesion' }
set.seed(1);ll<- stats::rchisq(mesh.for.drawing.curve, 1)
lll<- 0.99+ll
x.FROC<-append(ll,lll)
y.FROC.truth <- 1 - stats::pnorm( b.truth*stats::qnorm(exp(-x.FROC)) - a.truth )
if( missing(upper_y)){ upper_y <- 1.0 }
suppressWarnings(graphics::par(new=TRUE));
plot(x.FROC,y.FROC.truth,
col = 'red',
cex= 0.1 ,
xlim = c(0,upper_x ),ylim = c(0,upper_y),
xlab = xlabel,
ylab = 'cumulative hit per lesion'
,main =""
);
}#Def of function
Try the BayesianFROC package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.