R/validation_error_srsc.R

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 Error of estimates with respect to truth
#'@description
#'
#'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 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()
#'
#'
#'
# ####1#### ####2#### ####3#### ####4#### ####5#### ####6#### ####7#### ####8#### ####9####
#'#=============================  The second example ======================================
#'
#'#   If user do 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 estimates, we run the following:
#'
#'   fit <- fit_Bayesian_FROC(dataList.Chakra.1,ite = 1111,summary =FALSE,cha=3)
#'
#'
#'
#'
#'
#'
#'#  Secondly, extract the expected a posterior estimators (EAPs) from the object fit
#'
#'
#'   z <- rstan::get_posterior_mean(fit,par=c("z"))[,"mean-all chains"]
#'
#'
#'
#'
#'
#'#  Thirdly we use this z as a true values.
#'
#'
#'    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 do converge or not.
#'
#'     stan_trace( a$fit[[3]],pars = "A")
#'
#'
#'    # Check p value

#'     ppp( a$fit[[3]])
#'
#'
#'
#'     # In the above example, the posterior predictive p value is enough large,
#'     # but the model did not converge in R that criteria, which will cause
#'     # that the model does not fit to data. However p value is said
#'     # we can not reject the null hypothesis that the model does fit.
#'     # The author think this contradiction cause that the
#'     # number of MCMC iterations are too small which leads us to incorrect
#'     # Monte Carlo integral for p value. Thu p value is not correct.
#'     # Calculation of p value relies on the law of large number and thus
#'     # to obtain reliable posterior predictive p value, we need enough large
#'     # MCMC samples. 2019 August 29
#'
#'
#'
#'
#'
#'
#'                                           # Revised in 2019 August 29
#'
#'
#'
#'
#'
#'#========================================================================================
#'#            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 does not generate negative biases
#'#========================================================================================
#'
#'
#'  validation.dataset_srsc(absolute.errors = TRUE)
#'
#'
#'}# 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,
  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
    # Here fitting  ------

    fit[[seed]]   <- fit_Bayesian_FROC(
      dataList=list.of.dataList[[seed]],
      ite = ite,
      cha = cha,
      new.imaging.device = FALSE,
      summary = FALSE,
      ModifiedPoisson = ModifiedPoisson,
      DrawCurve = FALSE
    )


    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

  }#for seed

  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:", 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.
#'
#'
#'
#'
#'}# 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


){
 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 ,

      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. 13, 2021, 5:22 a.m.