R/BayesCTDesigncode.R

Defines functions plot.bayes_ctd_array print.bayes_ctd_array simple_sim historic_sim

Documented in historic_sim plot.bayes_ctd_array print.bayes_ctd_array simple_sim

#' Two Arm Bayesian Clinical Trial Simulation with Historical Data
#'
#' \code{historic_sim()} returns an S3 object of class \code{bayes_ctd_array}, which
#' will contain simulation results for power, statistic estimation, bias,
#' variance, and mse as requested by user.
#'
#' The object \code{bayes_ctd_array} has 6 elements: a list containing simulation
#' results (\code{data}), copies of the 4 function arguments \code{subj_per_arm},
#' \code{a0_vals}, \code{effect_vals}, and \code{rand_control_diff}, and finally
#' a \code{objtype} value indicating that \code{historic_sim()} was used. Each element of
#' \code{data} is a four-dimensional array, where each dimension is determined by the
#' length of parameters \code{subj_per_arm}, \code{a0_vals}, \code{effect_vals}, and
#' \code{rand_control_diff}. The size of \code{data} depends on which results are
#' requested by the user. At a minimum, at least one of \code{subj_per_arm},
#' \code{a0_vals}, \code{effect_vals}, or \code{rand_control_diff} must contain at
#' least 2 values, while the other three must contain at least 1 value. The \code{data}
#' list will always contain two elements: an array of power results (\code{power}) and
#' an array of estimation results (\code{est}). In addition to \code{power} and
#' \code{est}, data may also contain elements \code{var}, \code{bias}, or \code{mse},
#' depending on the values of \code{get_var}, \code{get_bias}, and \code{get_mse}. The
#' values returned in \code{est} are in the form of hazard ratios, mean ratios, odds
#' ratios, or mean differences depending on the value of \code{outcome_type}.  For a
#' Gaussian outcome, the estimation results are differences in group means (experimental
#' group minus control group). For a logistic outcome, the estimation results are odds
#' ratios (experimental group over control group). For lognormal and Poisson outcomes,
#' the estimation results are mean ratios (experimental group over control group). For a
#' piecewise exponential or a Weibull outcome, the estimation results are hazard
#' ratios (experimental group over control group).  The values returned in \code{bias},
#' \code{var}, and \code{mse} are on the scale of the values returned in
#' \code{est}.
#'
#' The object \code{bayes_ctd_array} has two primary methods, \code{print()} and
#' \code{plot()}, for printing and plotting slices of the arrays contained in
#' \code{bayes_ctd_array$data}.
#'
#' As dimensions of the four dimensional array increases, the time required to complete
#' the simulation will increase; however, it will be faster than a similar simulation
#' based on repeated calls to MCMC routines to analyze each simulated trial.
#'
#' The meaning of the estimation results, and the test used to generate power results,
#' depends on the outcome used. In all cases, power is based on a two-sided test
#' involving a (1-alpha)100\% credible interval, where the interval is used to determine
#' if the null hypothesis should be rejected (null value outside of the interval) or
#' not rejected (null value inside the interval). For a Gaussian outcome, the 95\%
#' credible interval is an interval for the difference in group means
#' (experimental group minus control group), and the test determines if 0 is in or
#' outside of the interval. For a Bernoulli outcome, the 95\% credible interval
#' is an interval for the odds ratio (experimental group over control group),
#' and the test determines if 1 is in or outside of the interval. For a lognormal or
#' a Poisson outcome, the 95\% credible interval is an interval for the mean ratio
#' (experimental group over control group), and the test determines if 1 is in or
#' outside of the interval. Finally, for a piecewise exponential or a Weibull outcome,
#' the 95\% credible interval is an interval for the hazard ratio (experimental group
#' over control group), and the test determines if 1 is in or outside of the interval.
#'
#' Please refer to the examples for illustration of package use.
#'
#' @param trial_reps Number of trials to replicate within each combination of
#'   \code{a0_vals}, \code{subj_per_arm}, \code{effect_vals}, and \code{rand_control_parms}.
#'   As the number of trials increases, the precision of the estimate will increase.
#'   Default is 100.
#' @param outcome_type Outcome distribution. Must be equal to \code{weibull},
#'   \code{lognormal}, \code{pwe} (Piecewise Exponential), \code{gaussian},
#'   \code{bernoulli}, or \code{poisson}.  Default is \code{weibull}.
#' @param subj_per_arm A vector of sample sizes, all of which must be positive
#'   integers.  Default is \code{c(50, 100, 150, 200, 250)}.
#' @param a0_vals A vector of power prior parameters ranging from 0 to 1, where 0
#'   implies no information from historical data should be used, and 1 implies all of
#'   the information from historical data should be used.  A value between 0 and 1
#'   implies that a proportion of the information from historical data will be used.
#'   Default is \code{c(0, 0.33, 0.67, 1)}.
#' @param effect_vals A vector of effects that should be reasonable for the
#'   outcome_type being studied, hazard ratios for Weibull, odds ratios for
#'   Bernoulli, mean ratios for Poisson, etc..  When \code{effect_vals} contain
#'   the null effect for a given \code{outcome_type}, the \code{power} component
#'   of \code{data} will contain an estimate of Type One Error.  In order to
#'   have a good set of Type One Error estimates, \code{trial_reps} need to be
#'   at least 10,000.  In such a case, if the total number of combinations
#'   made up from \code{subj_per_arm}, \code{a0_vals}, \code{effect_vals}, and
#'   \code{rand_control_diff} is very large, the time to complete the simulation
#'   can be substantial.  Default is \code{c(0.6, 1, 1.4)}.
#' @param rand_control_diff For piecewise exponential and Weibull outcomes, this is
#'   a vector of hazard ratios (randomized controls over historical controls)
#'   representing differences between historical and randomized controls.  For
#'   lognormal and Poisson outcomes, this is a vector of mean ratios (randomized
#'   controls over historical controls).  For a Bernoulli outcome, this is a vector
#'   of odds ratios (randomized controls over historical controls).  For a Gaussian
#'   outcome, this is a vector of mean differences (randomized minus historical
#'   controls). Default is \code{c(0.8, 1, 1.2)}.
#' @param hist_control_data A dataset of historical data.  Default is \code{NULL}.
#'   For survival outcomes, historical datasets must have 4 columns: id, treatment,
#'   event_time, and status.  The value of treatment should be 0.  For other
#'   outcomes, historical datasets must have columns: id, treatment, and y.
#' @param time_vec  A vector of time values which are used to create time periods
#'   within which the exponential hazard is constant.  Only used for piecewise
#'   exponential models.  Default is \code{NULL}.
#' @param censor_value A single value at which right censoring occurs when
#'   simulating randomized subject outcomes.  Used with survival outcomes.
#'   Default is \code{NULL}, where \code{NULL} implies no right censoring.
#' @param alpha A number ranging between 0 and 1 that defines the acceptable Type 1
#'   error rate. Default is 0.05.
#' @param get_var A TRUE/FALSE indicator of whether an array of variance
#'   estimates will be returned.  Default is \code{FALSE}.
#' @param get_bias A TRUE/FALSE indicator of whether an array of bias
#'   estimates will be returned.  Default is \code{FALSE}.
#' @param get_mse A TRUE/FALSE indicator of whether an array of MSE
#'   estimates will be returned.  Default is \code{FALSE}.
#' @param seedval A seed value for pseudo-random number generation.
#' @param quietly A TRUE/FALSE indicator of whether notes are printed
#'   to output about simulation progress as the simulation runs.  If
#'   running interactively in RStudio or running in the R console,
#'   \code{quietly} can be set to FALSE.  If running in a Notebook or
#'   knitr document, \code{quietly} needs to be set to TRUE.  Otherwise
#'   each note will be printed on a separate line and it will take up
#'   a lot of output space.  Default is \code{TRUE}.
#'
#' @return \code{historic_sim()} returns an S3 object of class \code{bayes_ctd_array}.
#'   As noted in details, an object of class \code{bayes_ctd_array }has 6 elements: a
#'   list of simulation results (\code{data}), copies of the 4 function arguments
#'   \code{subj_per_arm}, \code{a0_vals}, \code{effect_vals}, and
#'   \code{rand_control_diff}, and finally \code{objtype} indicating that \code{historic_sim()}
#'   was used. See details for a discussion about the contents of
#'   \code{data}. Results from the simulation contained in the \code{bayes_ctd_array}
#'   object can be printed or plotted using the \code{print()} and
#'   \code{plot()} methods. The results can also be accessed using basic list
#'   element identification and array slicing. For example, to get the 4-dimensional
#'   array of power results from a simulation, one could use the code
#'   \code{bayes_ctd_array$data$power}, where \code{bayes_ctd_array} is replaced
#'   with the name of the variable containing the \code{bayes_ctd_array} object. If
#'   one wanted a table of power for sample size by a0, while holding effect equal to
#'   the first considered value and control differences equal to the second considered
#'   value, then the code is \code{bayes_ctd_array$data$power[,,1,2]}, where
#'   \code{bayes_ctd_array} is replaced with the name of the variable containing the
#'   \code{bayes_ctd_array} object.
#'
#' @references Eggleston et al. BayesCTDesign: An R Package for Bayesian Trial Design Using Historical Control Data. \emph{Journal of Statistical Software}, November 2021, Volume 100, Issue 21,  <doi:10.18637/jss.v100.i21>.
#'
#' @examples
#' #Generate a sample of historical data for use in example.
#' set.seed(2250)
#' SampleHistData <- genweibulldata(sample_size=60, scale1=2.82487,
#'                                  hazard_ratio=0.6, common_shape=3,
#'                                  censor_value=3)
#' histdata <- subset(SampleHistData, subset=(treatment==0))
#' histdata$id <- histdata$id+10000
#'
#' #Run a Weibull simulation, using historic_sim().
#' #For meaningful results, trial_reps needs to be much larger than 2.
#' weibull_test <- historic_sim(trial_reps = 2, outcome_type = "weibull",
#'                              subj_per_arm = c(50, 100, 150),
#'                              a0_vals = c(0, 0.50, 1),
#'                              effect_vals = c(0.6, 1),
#'                              rand_control_diff = c(0.8, 1),
#'                              hist_control_data = histdata, time_vec = NULL,
#'                              censor_value = 3, alpha = 0.05, get_var = TRUE,
#'                              get_bias = TRUE, get_mse = TRUE, seedval=123,
#'                              quietly=TRUE)
#'
#' #Tabulate the simulation results for power.
#' test_table <- print(x=weibull_test, measure="power",
#'                     tab_type="WX|YZ", effect_val=0.6,
#'                     rand_control_diff_val=1.0)
#' print(test_table)
#'
#' \donttest{
#' #Create a plot of the power simulation results.
#' plot(x=weibull_test, measure="power", tab_type="WX|YZ",
#'      smooth=FALSE, plot_out=TRUE, effect_val=0.6,
#'      rand_control_diff_val=1.0)
#' #Create a plot of the estimated hazard ratio simulation results.
#' plot(x=weibull_test, measure="est", tab_type="WX|YZ",
#'      smooth=FALSE, plot_out=TRUE, effect_val=0.6,
#'      rand_control_diff_val=1.0)
#' #Create a plot of the hazard ratio variance simulation results.
#' plot(x=weibull_test, measure="var", tab_type="WX|YZ",
#'      smooth=FALSE, plot_out=TRUE, effect_val=0.6,
#'      rand_control_diff_val=1.0)
#' #Create a plot of the hazard ratio bias simulation results.
#' plot(x=weibull_test, measure="bias", tab_type="WX|YZ",
#'      smooth=FALSE, plot_out=TRUE, effect_val=0.6,
#'      rand_control_diff_val=1.0)
#' #Create a plot of the hazard ratio mse simulation results.
#' plot(x=weibull_test, measure="mse", tab_type="WX|YZ",
#'      smooth=FALSE, plot_out=TRUE, effect_val=0.6,
#'      rand_control_diff_val=1.0)
#'
#' #Create other power plots using different values for tab_type
#' plot(x=weibull_test, measure="power", tab_type="XY|WZ",
#'      smooth=FALSE, plot_out=TRUE, subj_per_arm_val=150,
#'      rand_control_diff_val=1.0)
#'
#' plot(x=weibull_test, measure="power", tab_type="XZ|WY",
#'      smooth=FALSE, plot_out=TRUE, subj_per_arm_val=150, effect_val=0.6)
#'
#' plot(x=weibull_test, measure="power", tab_type="YZ|WX",
#'      smooth=FALSE, plot_out=TRUE, subj_per_arm_val=150, a0_val=0.5)
#'
#' plot(x=weibull_test, measure="power", tab_type="WY|XZ",
#'      smooth=FALSE, plot_out=TRUE, rand_control_diff_val=1, a0_val=0.5)
#'
#' plot(x=weibull_test, measure="power", tab_type="WZ|XY",
#'      smooth=FALSE, plot_out=TRUE, effect_val=0.6, a0_val=0.5)
#' }
#'
#' \donttest{
#' #Run Poisson simulation, using historic_sim(), but set two design characteristic
#' # parameters to only 1 value.
#' #Note: historic_sim() can take a while to run.
#' #Generate a sample of historical poisson data for use in example.
#' set.seed(2250)
#' samplehistdata <- genpoissondata(sample_size=60, mu1=1, mean_ratio=1.0)
#' histdata <- subset(samplehistdata, subset=(treatment==0))
#' histdata$id <- histdata$id+10000
#'
#' #For meaningful results, trial_reps needs to be larger than 100.
#' poisson_test <- historic_sim(trial_reps = 100, outcome_type = "poisson",
#'                               subj_per_arm = c(50, 75, 100, 125, 150, 175, 200, 225, 250),
#'                               a0_vals = c(1),
#'                               effect_vals = c(0.6),
#'                               rand_control_diff = c(0.6, 1, 1.6),
#'                               hist_control_data = histdata, time_vec = NULL,
#'                               censor_value = 3, alpha = 0.05, get_var = TRUE,
#'                               get_bias = TRUE, get_mse = TRUE, seedval=123,
#'                               quietly=TRUE)
#'
#' #Tabulate the simulation results for power.
#' test_table <- print(x=poisson_test, measure="power",
#'                     tab_type=NULL)
#' print(test_table)
#'
#' #Create a plot of the power simulation results.
#' plot(x=poisson_test, measure="power", tab_type=NULL,
#'      smooth=FALSE, plot_out=TRUE)
#' }
#'
#' \donttest{
#' #At least one of subj_per_arm, a0_vals, effect_vals, or rand_control_diff
#' #must contain at least 2 values.
#' #Generate a sample of historical lognormal data for use in example.
#' set.seed(2250)
#' samplehistdata <- genlognormaldata(sample_size=60, mu1=1.06, mean_ratio=0.6, common_sd=1.25,
#'                                    censor_value=3)
#' histdata <- subset(samplehistdata, subset=(treatment==0))
#' histdata$id <- histdata$id+10000
#'
#' #Run a Lognormal simulation, using historic_sim().
#' #For meaningful results, trial_reps needs to be larger than 100.
#' lognormal_test <- historic_sim(trial_reps = 100, outcome_type = "lognormal",
#'                                subj_per_arm = c(25,50,75,100,125,150,175,200,225,250),
#'                                a0_vals = c(1.0),
#'                                effect_vals = c(0.6),
#'                                rand_control_diff = c(1.8),
#'                                hist_control_data = histdata, time_vec = NULL,
#'                                censor_value = 3, alpha = 0.05, get_var = TRUE,
#'                                get_bias = TRUE, get_mse = TRUE, seedval=123,
#'                                quietly=TRUE)
#'
#' test_table <- print(x=lognormal_test, measure="power",
#'                     tab_type=NULL)
#' print(test_table)
#' #Create a plot of the power simulation results.
#' plot(x=lognormal_test, measure="power", tab_type=NULL,
#'      smooth=TRUE, plot_out=TRUE)
#' }
#'
#' @export
historic_sim <- function(trial_reps = 100, outcome_type = "weibull", subj_per_arm = c(50, 100, 150, 200, 250), a0_vals = c(0,
    0.33, 0.67, 1), effect_vals = c(0.6, 1, 1.4), rand_control_diff = c(0.8, 1, 1.2), hist_control_data = NULL, time_vec = NULL,
    censor_value = NULL, alpha = 0.05, get_var = FALSE, get_bias = FALSE, get_mse = FALSE, seedval=NULL, quietly=TRUE){

	#set random seed
	set.seed(seedval)
    #------------- Go through all the high level checks for proper input. -------------#
    global_error_checks(outcome_type, a0_vals, subj_per_arm, hist_control_data, rand_control_diff, get_var, get_bias,
        get_mse)

    #-- Given outcome type go through low level checks for proper input and run simulation --#
    if (tolower(outcome_type) == "weibull") {
        weibull_error_checks(effect_vals, hist_control_data, rand_control_diff, censor_value, alpha)
        results <- weibull_sim(trial_reps = trial_reps, subj_per_arm = subj_per_arm, a0_vals = a0_vals, effect_vals = effect_vals,
            rand_control_diff = rand_control_diff, hist_control_data = hist_control_data, censor_value = censor_value,
            alpha = alpha, get_var = get_var, get_bias = get_bias, get_mse = get_mse, quietly=quietly)

    } else if (tolower(outcome_type) == "lognormal") {
        lognormal_error_checks(effect_vals, hist_control_data, rand_control_diff, censor_value, alpha)
        results <- lognormal_sim(trial_reps = trial_reps, subj_per_arm = subj_per_arm, a0_vals = a0_vals, effect_vals = effect_vals,
            rand_control_diff = rand_control_diff, hist_control_data = hist_control_data, censor_value = censor_value,
            alpha = alpha, get_var = get_var, get_bias = get_bias, get_mse = get_mse, quietly=quietly)

    } else if (tolower(outcome_type) == "pwe") {
        pwe_error_checks(effect_vals, hist_control_data, rand_control_diff, time_vec, censor_value, alpha)
        results <- pwe_sim(trial_reps = trial_reps, subj_per_arm = subj_per_arm, a0_vals = a0_vals, effect_vals = effect_vals,
            rand_control_diff = rand_control_diff, hist_control_data = hist_control_data, time_vec_val = time_vec, censor_value = censor_value,
            alpha = alpha, get_var = get_var, get_bias = get_bias, get_mse = get_mse, quietly=quietly)

    } else if (tolower(outcome_type) == "gaussian") {
        gaussian_error_checks(effect_vals, hist_control_data, rand_control_diff, alpha)
        results <- gaussian_sim(trial_reps = trial_reps, subj_per_arm = subj_per_arm, a0_vals = a0_vals, effect_vals = effect_vals,
            rand_control_diff = rand_control_diff, hist_control_data = hist_control_data, alpha = alpha, get_var = get_var,
            get_bias = get_bias, get_mse = get_mse, quietly=quietly)

    } else if (tolower(outcome_type) == "bernoulli") {
        bernoulli_error_checks(effect_vals, hist_control_data, rand_control_diff, alpha)
        results <- bernoulli_sim(trial_reps = trial_reps, subj_per_arm = subj_per_arm, a0_vals = a0_vals, effect_vals = effect_vals,
            rand_control_diff = rand_control_diff, hist_control_data = hist_control_data, alpha = alpha, get_var = get_var,
            get_bias = get_bias, get_mse = get_mse, quietly=quietly)

    } else if (tolower(outcome_type) == "poisson") {
        poisson_error_checks(effect_vals, hist_control_data, rand_control_diff, alpha)
        results <- poisson_sim(trial_reps = trial_reps, subj_per_arm = subj_per_arm, a0_vals = a0_vals, effect_vals = effect_vals,
            rand_control_diff = rand_control_diff, hist_control_data = hist_control_data, alpha = alpha, get_var = get_var,
            get_bias = get_bias, get_mse = get_mse, quietly=quietly)
    }

    if ((length(subj_per_arm) == 1 & length(a0_vals) == 1 & length(effect_vals) == 1) |
       (length(subj_per_arm) == 1 & length(a0_vals) == 1 & length(rand_control_diff) == 1) |
       (length(subj_per_arm) == 1 & length(effect_vals) == 1 & length(rand_control_diff) == 1) |
       (length(a0_vals) == 1 & length(effect_vals) == 1 & length(rand_control_diff) == 1)){
         results$objtype <- 'realsimple'
    }
    results
}

#' Two Arm Bayesian Clinical Trial Simulation without Historical Data
#'
#' \code{simple_sim()} returns an S3 object of class \code{bayes_ctd_array}, which
#' will contain simulation results for power, statistic estimation, bias, variance,
#' and mse as requested by user.
#'
#' The object \code{bayes_ctd_array} has 6 elements: a list containing simulation
#' results (\code{data}), copies of the 4 function arguments \code{subj_per_arm},
#' \code{a0_vals}, \code{effect_vals}, and \code{rand_control_diff}, and finally
#' a \code{objtype} value indicating that \code{simple_sim()} was used. Each element of
#' \code{data} is a four-dimensional array, where each dimension is determined by the
#' length of parameters \code{subj_per_arm}, \code{a0_vals}, \code{effect_vals}, and
#' \code{rand_control_diff}. The size of \code{data} depends on which results are
#' requested by the user. At a minimum, at least one of \code{subj_per_arm},
#' \code{a0_vals}, \code{effect_vals}, or \code{rand_control_diff} must contain at
#' least 2 values, while the other three must contain at least 1 value.  The \code{data}
#' list will always contain two elements: an array of power results (\code{power}) and
#' an array of estimation results (\code{est}).  In addition to \code{power} and
#' \code{est}, \code{data} may also contain elements \code{var}, \code{bias}, or
#' \code{mse}, depending on the values of \code{get_var}, \code{get_bias}, and
#' \code{get_mse}. The values returned in \code{est} are in the form of hazard ratios,
#' mean ratios, odds ratios, or mean differences depending on the value of
#' \code{outcome_type}.   For a Gaussian outcome, the estimation results are
#' differences in group means (experimental group minus control group). For a
#' logistic outcome, the estimation results are odds ratios (experimental group over
#' control group). For lognormal and Poisson outcomes, the estimation results are mean
#' ratios (experimental group over control group). For a piecewise exponential or a
#' Weibull outcome, the estimation results are hazard ratios (experimental group over
#' control group).  The values returned in \code{bias}, \code{var}, and \code{mse} are
#' on the scale of the values returned in \code{est}.
#'
#' The object \code{bayes_ctd_array} has two primary methods, \code{print()} and
#' \code{plot()}, for printing and plotting slices of the arrays contained in
#' \code{bayes_ctd_array$data}.
#'
#' As dimensions of the four dimensional array increases, the time required to complete
#' the simulation will increase; however, it will be faster than a similar simulation
#' based on repeated calls to MCMC routines to analyze each simulated trial.
#'
#' The meaning of the estimation results, and the test used to generate power results,
#' depends on the outcome used. In all cases, power is based on a two-sided test
#' involving a (1-alpha)100\% credible interval, where the interval is used to determine
#' if the null hypothesis should be rejected (null value outside of the interval) or
#' not rejected (null value inside the interval). For a Gaussian outcome, the 95\%
#' credible interval is an interval for the difference in group means
#' (experimental group minus control group), and the test determines if 0 is in or
#' outside of the interval. For a Bernoulli outcome, the 95\% credible interval
#' is an interval for the odds ratio (experimental group over control group),
#' and the test determines if 1 is in or outside of the interval. For a lognormal or
#' a Poisson outcome, the 95\% credible interval is an interval for the mean ratio
#' (experimental group over control group), and the test determines if 1 is in or
#' outside of the interval. Finally, for a piecewise exponential or a Weibull outcome,
#' the 95\% credible interval is an interval for the hazard ratio (experimental group
#' over control group), and the test determines if 1 is in or outside of the interval.
#'
#' For a Gaussian outcome, the \code{control_parms} values should be \code{(mean, sd)},
#' where mean is the mean parameter for the control group used in a call to \code{rnorm()},
#' and sd is the common sd parameter for both groups used in a call to\code{rlnorm()}.
#'
#' For a Bernoulli outcome, the \code{control_parms} values should be \code{(prob)}, where
#' prob is the event probability for the control group used in a call to \code{rbinom()}.
#'
#' For a lognormal outcome, the \code{control_parms} values should be \code{(meanlog, sdlog)},
#' where meanlog is the meanlog parameter for the control group used in a call to
#' \code{rlnorm()}, and sdlog is the common sdlog parameter for both groups used in
#' a call to \code{rlnorm()}.
#'
#' For a Poisson outcome, the \code{control_parms} value should be \code{(lambda)}, where
#' lambda is the lambda parameter for the control group used in a call to \code{rpois()} and
#' is equal to the mean of a Poisson distribution.
#'
#' For a Weibull outcome, the \code{control_parms} values should be \code{(scale, shape)},
#' where scale is the scale parameter for the control group used in a call to
#' \code{rweibull()}, and shape is the common shape parameter for both groups used in
#' a call to \code{rweibull()}.
#'
#' For a piecewise exponential outcome, the \code{control_parms} values should be a vector
#' of lambdas used in a call to \code{eha::rpch()}.  Each element in \code{control_parms}
#' is a hazard for an interval defined by the \code{time_vec} parameter.
#'
#' Please refer to the examples for illustration of package use.
#'
#' @param trial_reps Number of trials to replicate within each combination of
#'   \code{a0_vals}, \code{subj_per_arm}, \code{effect_vals}, and \code{rand_control_parms}.
#'   As the number of trials increases, the precision of the estimate will increase.
#'   Default is 100.
#' @param outcome_type Outcome distribution. Must be equal to \code{weibull},
#'   \code{lognormal}, \code{pwe} (Piecewise Exponential), \code{gaussian},
#'   \code{bernoulli}, or \code{poisson}.  Default is \code{weibull}.
#' @param subj_per_arm A vector of sample sizes, all of which must be positive
#'   integers.  Default is \code{c(50, 100, 150, 200, 250)}.
#' @param effect_vals A vector of effects that should be reasonable for the
#'   outcome_type being studied, hazard ratios for Weibull, odds ratios for
#'   Bernoulli, mean ratios for Poisson, etc..  When \code{effect_vals} contain
#'   the null effect for a given \code{outcome_type}, the \code{power} component
#'   of \code{data} will contain an estimate of Type One Error.  In order to
#'   have a good set of Type One Error estimates, \code{trial_reps} need to be
#'   at least 10,000.  In such a case, if the total number of combinations
#'   made up from \code{subj_per_arm}, \code{a0_vals}, \code{effect_vals}, and
#'   \code{rand_control_diff} is very large, the time to complete the simulation
#'   can be substantial. Default is \code{c(0.6, 1, 1.4)}.
#' @param control_parms A vector of parameter values defining the outcome
#'   distribution for randomized controls. See Details for what is required for
#'   each \code{outcome_type}.
#' @param time_vec  A vector of time values that are used to create time periods
#'   within which the exponential hazard is constant.  Only used for piecewise
#'   exponential models.  Default is \code{NULL}.
#' @param censor_value A single value at which right censoring occurs when
#'   simulating randomized subject outcomes.  Used with survival outcomes.
#'   Default is \code{NULL}, where \code{NULL} implies no right censoring.
#' @param alpha A number ranging between 0 and 1 that defines the acceptable Type 1
#'   error rate. Default is 0.05.
#' @param get_var A TRUE/FALSE indicator of whether an array of variance
#'   estimates will be returned.  Default is \code{FALSE}.
#' @param get_bias A TRUE/FALSE indicator of whether an array of bias
#'   estimates will be returned.  Default is \code{FALSE}.
#' @param get_mse A TRUE/FALSE indicator of whether an array of MSE
#'   estimates will be returned.  Default is \code{FALSE}.
#' @param seedval A seed value for pseudo-random number generation.
#' @param quietly A TRUE/FALSE indicator of whether notes are printed
#'   to output about simulation progress as the simulation runs.  If
#'   running interactively in RStudio or running in the R console,
#'   \code{quietly} can be set to FALSE.  If running in a Notebook or
#'   knitr document, \code{quietly} needs to be set to TRUE.  Otherwise
#'   each note will be printed on a separate line and it will take up
#'   a lot of output space.  Default is \code{TRUE}.
#'
#' @return \code{simple_sim()} returns an S3 object of class \code{bayes_ctd_array}.
#'   As noted in Details, an object of class \code{bayes_ctd_array} has 6 elements: a
#'   list containing simulation results (\code{data}), copies of the 4 function
#'   arguments \code{subj_per_arm}, \code{a0_vals}, \code{effect_vals}, and
#'   \code{rand_control_diff}, and finally \code{objtype} indicating that \code{simple_sim()}
#'   was used. See Details for a discussion about the contents of
#'   \code{data}. Results from the simulation contained in the \code{bayes_ctd_array}
#'   object can be printed or plotted using the \code{print()} and
#'   \code{plot()} methods. The results can also be accessed using basic list
#'   element identification and array slicing. For example, to get the power results
#'   from a simulation, one could use the code \code{bayes_ctd_array$data$power}, where
#'   \code{bayes_ctd_array} is replaced with the name of the variable containing the
#'   \code{bayes_ctd_array} object. Even though this is a 4-dimensional array, the power
#'   results only occupy a single 2-dimensional table. To print this 2-dimensional table,
#'   one would use the code \code{bayes_ctd_array$data$power[,1,,1]}, where
#'   \code{bayes_ctd_array} is replaced with the name of the variable containing the
#'   \code{bayes_ctd_array} object.
#'
#' @references Eggleston et al. (2021). BayesCTDesign: An R Package for Bayesian Trial Design Using Historical Control Data. \emph{Journal of Statistical Software}, November 2021, Volume 100, Issue 21,  <doi:10.18637/jss.v100.i21>.
#'
#'
#' @examples
#' #Run a Weibull simulation, using simple_sim().
#' #For meaningful results, trial_reps needs to be much larger than 2.
#' weibull_test <- simple_sim(trial_reps = 2, outcome_type = "weibull",
#'                            subj_per_arm = c(50, 100, 150, 200),
#'                            effect_vals = c(0.6, 1, 1.4),
#'                            control_parms = c(2.82487,3), time_vec = NULL,
#'                            censor_value = NULL, alpha = 0.05,
#'                            get_var = TRUE, get_bias = TRUE, get_mse = TRUE,
#'                            seedval=123, quietly=TRUE)
#'
#' #Tabulate the simulation results for power.
#' test_table <- print(x=weibull_test, measure="power",
#'                     tab_type=NULL, subj_per_arm_val=NULL, a0_val=NULL,
#'                     effect_val=NULL, rand_control_diff_val=NULL)
#' print(test_table)
#'
#' #Create a plot of the power simulation results.
#' plot(x=weibull_test, measure="power", tab_type=NULL,
#'      smooth=FALSE, plot_out=TRUE)
#' #Create a plot of the estimated hazard ratio simulation results.
#' plot(x=weibull_test, measure="est", tab_type=NULL,
#'      smooth=FALSE, plot_out=TRUE)
#' #Create a plot of the hazard ratio variance simulation results.
#' plot(x=weibull_test, measure="var", tab_type=NULL,
#'      smooth=FALSE, plot_out=TRUE)
#' #Create a plot of the hazard ratio bias simulation results.
#' plot(x=weibull_test, measure="bias", tab_type=NULL,
#'      smooth=FALSE, plot_out=TRUE)
#' #Create a plot of the hazard ratio mse simulation results.
#' plot(x=weibull_test, measure="mse", tab_type=NULL,
#'      smooth=FALSE, plot_out=TRUE)
#' @export
simple_sim <- function(trial_reps = 100, outcome_type = "weibull", subj_per_arm = c(50, 100, 150, 200, 250), effect_vals = c(0.6,1,
    1.4), control_parms = NULL, time_vec = NULL, censor_value = NULL, alpha = 0.05, get_var = FALSE, get_bias = FALSE,
    get_mse = FALSE, seedval=NULL, quietly=TRUE) {

	#set random seed
	set.seed(seedval)
    #------------- Go through all the high level checks for proper input. -------------#
    global_error_checks_simple(outcome_type, subj_per_arm, get_var, get_bias, get_mse)

    #-- Given outcome type go through low level checks for proper input and run simulation --#
    if (tolower(outcome_type) == "weibull") {
        weibull_error_checks_simple(effect_vals, control_parms, censor_value, alpha)
        results <- simple_weibull_sim(trial_reps = trial_reps, subj_per_arm = subj_per_arm, effect_vals = effect_vals,
            scale1_value = control_parms[1], common_shape_value = control_parms[2], censor_value = censor_value, alpha = alpha,
            get_var = get_var, get_bias = get_bias, get_mse = get_mse, quietly=quietly)

    } else if (tolower(outcome_type) == "lognormal") {
        lognormal_error_checks_simple(effect_vals, control_parms, censor_value, alpha)
        results <- simple_lognormal_sim(trial_reps = trial_reps, subj_per_arm = subj_per_arm, effect_vals = effect_vals,
            mu1_val = control_parms[1], common_sd_val = control_parms[2], censor_value = censor_value, alpha = alpha,
            get_var = get_var, get_bias = get_bias, get_mse = get_mse, quietly=quietly)

    } else if (tolower(outcome_type) == "pwe") {
        pwe_error_checks_simple(effect_vals, time_vec, control_parms, censor_value, alpha)
        results <- simple_pwe_sim(trial_reps = trial_reps, subj_per_arm = subj_per_arm, effect_vals = effect_vals, time_vec_val = time_vec,
            rc_hazards = control_parms, censor_value = censor_value, alpha = alpha, get_var = get_var, get_bias = get_bias,
            get_mse = get_mse, quietly=quietly)

    } else if (tolower(outcome_type) == "gaussian") {
        gaussian_error_checks_simple(effect_vals, control_parms, alpha)
        results <- simple_gaussian_sim(trial_reps = trial_reps, subj_per_arm = subj_per_arm, effect_vals = effect_vals,
            mu1_val = control_parms[1], common_sd_val = control_parms[2], alpha = alpha, get_var = get_var, get_bias = get_bias,
            get_mse = get_mse, quietly=quietly)

    } else if (tolower(outcome_type) == "bernoulli") {
        bernoulli_error_checks_simple(effect_vals, control_parms, alpha)
        results <- simple_bernoulli_sim(trial_reps = trial_reps, subj_per_arm = subj_per_arm, effect_vals = effect_vals,
            prob1_val = control_parms[1], alpha = alpha, get_var = get_var, get_bias = get_bias, get_mse = get_mse, quietly=quietly)

    } else if (tolower(outcome_type) == "poisson") {
        poisson_error_checks_simple(effect_vals, control_parms, alpha)
        results <- simple_poisson_sim(trial_reps = trial_reps, subj_per_arm = subj_per_arm, effect_vals = effect_vals,
            mu1_val = control_parms[1], alpha = alpha, get_var = get_var, get_bias = get_bias, get_mse = get_mse, quietly=quietly)
    }
}



#' Print Data from Two Arm Bayesian Clinical Trial Simulation.
#'
#' \code{print.bayes_ctd_array()} takes an S3 object of class \code{bayes_ctd_array}, and
#' prints a two dimensional slice from the data generated by a clinical trial simulation
#' using \code{historic_sim()} or \code{simple_sim()}.
#'
#' If the object of class \code{bayes_ctd_array} is created by \code{historic_sim()},
#' then the function \code{print()} allows the user to print user-specified 1- and 2-
#' dimensional slices of the simulation results based on slicing code described
#' below.  If the object of class \code{bayes_ctd_array} is created by
#' \code{simple_sim()}, a basic table of characteristic by sample size and effect is created.
#'
#' If the object of class \code{bayes_ctd_array} is created by \code{simple_sim()}, then
#' all four trial characteristics (\code{subj_per_arm_val}, \code{a0_vals},
#' \code{effect_val}, and \code{rand_control_diff_val}) can be ignored, as can the
#' parameter defining what type of table to print, \code{tab_type}.  A call to
#' \code{print()} will require the user to specify a measure (power, est, var, bias,
#' or mse).
#'
#' If the object of class \code{bayes_ctd_array} is created by \code{historic_sim()},
#' a call to \code{print()} will require the user to specify a measure
#' (power, est, var, bias, or mse) and may require the user to specify a table type.
#' A table type, \code{tab_type}, will be required if 3 of the 4 trial characteristics
#' are equal to a vector of 2 or more values.  The table type specification
#' uses the letters W, X, Y, and Z.  The letter W represents the subject per arm
#' dimension.  The letter X represents the a0 dimension.  The letter Y represents
#' the effect dimension.  The letter Z represents the control difference dimension.
#' To define a slice of the 4-dimensional array, these letters are put into an AB|CD
#' pattern.  The two letters to the right of the vertical bar define which variables
#' are held constant.  The two letters to the left of the vertical bar define which
#' variables are going to show up in the rows (first letter) and in the columns (second
#' letter).  For example if tab_type equals \code{WX|YZ}, then effect and control
#' differences will be held constant, while sample size will be represented by the rows
#' in the generated table and a0 values will be represented by the columns.  The actual
#' values that are printed in the tables depend on what measure is requested in the
#' parameter \code{measure}.
#'
#' \itemize{
#'   \item \code{tab_type='WX|YZ'}, Sample Size by a0
#'   \item \code{tab_type='WY|XZ'}, Sample Size by Effect
#'   \item \code{tab_type='WZ|XY'}, Sample Size by Control Differences
#'   \item \code{tab_type='XY|WZ'}, a0 by Effect
#'   \item \code{tab_type='XZ|WY'}, a0 by Control Differences
#'   \item \code{tab_type='YZ|WX'}, Effect by Control Differences
#'   \item \code{tab_type='ZX|WY'}, Control Differences by a0
#'   \item \code{tab_type='XW|YZ'}, a0 by Sample Size
#'   \item \code{tab_type='YW|XZ'}, Effect by Sample Size
#'   \item \code{tab_type='YX|WZ'}, Effect by a0
#'   \item \code{tab_type='ZW|XY'}, Control Differences by Sample Size
#'   \item \code{tab_type='ZY|WX'}, Control Differences by Effect
#' }
#'
#' It is very important to populate the values of \code{subj_per_arm_val},
#' \code{a0_vals}, \code{effect_val}, and \code{rand_control_diff_val} correctly given
#' the value of tab_type, when the object of class \code{bayes_ctd_array} is created by
#' \code{historic_sim()} and at least 3 of the four parameters have more than one
#' value.  On the other hand, if 2 or more of the four parameters have only one value,
#' then \code{subj_per_arm_val}, \code{a0_vals}, \code{effect_val},
#' \code{rand_control_diff_val}, as well as \code{tab_type} can be ignored.  If the last
#' two letters are \code{YZ}, then \code{effect_val} and \code{rand_control_diff_val}
#' must be populated.  If the last two letters are \code{XZ}, then \code{a0_vals} and
#' \code{rand_control_diff_val} must be populated.  If the last two letters are
#' \code{XY}, then \code{a0_vals} and \code{effect_val} must be populated.  If the last
#' two letters are \code{WZ}, then \code{sample_val} and \code{rand_control_diff_val}
#' must be populated.  If the last two letters are \code{WY}, then \code{sample_size_val}
#' and \code{effect_val} must be populated.  If the last two letters are \code{WX}, then
#' \code{sample_size_val} and \code{a0_vals} must be populated.
#'
#' If the object of class \code{bayes_ctd_array} is created by \code{simple_sim()}, the
#' parameters \code{tab_type}, \code{subj_per_arm_val}, \code{a0_vals}, \code{effect_val},
#' and \code{rand_control_diff_val} are ignored.
#'
#' @param x Name of object of class \code{bayes_ctd_array} containing
#'   data from clinical trial simulation.
#' @param measure Must be equal to \code{power}, \code{est}, \code{var}, \code{bias},
#'   or \code{mse}.  Default is \code{power}.  Case does not matter.
#' @param tab_type A character string that must equal \code{WX|YZ}, \code{WY|XZ},
#'   \code{WZ|XY}, \code{XY|WZ}, \code{XZ|WY}, \code{YZ|WX}, \code{ZX|WY}, \code{XW|YZ},
#'   \code{YW|XZ}, \code{YX|WZ}, \code{ZW|XY}, \code{ZX|WY}, \code{ZY|WX} when
#'   \code{x} is generated by \code{historic_sim()}.  Default is
#'   \code{WX|YZ}.  When \code{x} is generated by \code{simple_sim()},
#'   \code{tab_type} is ignored.
#' @param subj_per_arm_val Must be non-missing, if \code{x} is generated
#'   by \code{historic_sim()} and sample size is being held constant.
#'   If \code{x} is generated by \code{historic_sim()} and sample size
#'   is being held constant, \code{subj_per_arm_val} must equal a value submitted
#'   to \code{historic_sim()} within the \code{subj_per_arm} parameter.  When
#'   \code{x} is generated by \code{simple_sim()}, \code{subj_per_arm_val}
#'   is ignored.
#' @param a0_val Must be non-missing, if \code{x} is generated
#'   by \code{historic_sim()} and a0, the power prior parameter, is being held
#'   constant.  If \code{x} is generated by \code{historic_sim()} and
#'   a0 is being held constant, \code{a0_val} must equal a value submitted
#'   to \code{historic_sim()} within the \code{a0_vals} parameter.  When
#'   \code{x} is generated by \code{simple_sim()}, \code{a0_val} is
#'   ignored.
#' @param effect_val  Must be non-missing, if \code{x} is generated
#'   by \code{historic_sim()} and effect is being held constant.  If
#'   \code{x} is generated by \code{historic_sim()} and effect is being
#'   held constant, \code{effect_val} must equal a value submitted to
#'   \code{historic_sim()} within the \code{effect_vals} parameter.  When
#'   \code{x} is generated by \code{simple_sim()}, \code{effect_val} is
#'   ignored.
#' @param rand_control_diff_val Must be non-missing, if \code{x} is
#'   generated by \code{historic_sim()} and differences between randomized
#'   and historical controls are being held constant.  If \code{x}
#'   is generated by \code{historic_sim()} and control differences are being
#'   held constant, \code{rand_control_diff_val} must equal a value submitted to
#'   \code{historic_sim()} within the \code{rand_control_diff} parameter.  When
#'   \code{x} is generated by \code{simple_sim()},
#'   \code{rand_control_diff_val} is ignored.
#' @param print_chg_warn A parameter not used by the user, but is used by
#'   \code{plot()} to ensure warnings are not printed twice.
#' @param ...	further arguments passed to or from other methods.
#'
#' @return \code{print()} returns a two dimensional array of simulation results.
#'
#' @references Eggleston et al. BayesCTDesign: An R Package for Bayesian Trial Design Using Historical Control Data.  \emph{Journal of Statistical Software}, November 2021, Volume 100, Issue 21,  <doi:10.18637/jss.v100.i21>.
#'
#' @examples
#' #Run a Weibull simulation, using simple_sim().
#' #For meaningful results, trial_reps needs to be much larger than 2.
#' weibull_test <- simple_sim(trial_reps = 2, outcome_type = "weibull",
#'                            subj_per_arm = c(50, 100, 150, 200),
#'                            effect_vals = c(0.6, 1, 1.4),
#'                            control_parms = c(2.82487,3),
#'                            time_vec = NULL, censor_value = NULL,
#'                            alpha = 0.05, get_var = TRUE,
#'                            get_bias = TRUE, get_mse = TRUE,
#'                            seedval=123, quietly=TRUE)
#'
#' #Tabulate the simulation results for power.
#' test_table <- print(x=weibull_test, measure="power",
#'                     tab_type=NULL, subj_per_arm_val=NULL, a0_val=NULL,
#'                     effect_val=NULL, rand_control_diff_val=NULL)
#' print(test_table)
#'
#' #Tabulate the simulation results for estimates.
#' print(x=weibull_test, measure="est")
#'
#' #Tabulate the simulation results for variance.
#' print(x=weibull_test, measure="var")
#'
#' #Tabulate the simulation results for bias.
#' print(x=weibull_test, measure="bias")
#'
#' #Tabulate the simulation results for mse.
#' print(x=weibull_test, measure="mse")
#'
#' \donttest{
#' #Run another weibull simulation, using historic_sim().
#' #Note: historic_sim() can take a while to run.
#' #Generate a sample of historical data for use in example.
#' set.seed(2250)
#' SampleHistData <- genweibulldata(sample_size=60, scale1=2.82487,
#'                                  hazard_ratio=0.6, common_shape=3,
#'                                  censor_value=3)
#' histdata <- subset(SampleHistData, subset=(treatment==0))
#' histdata$id <- histdata$id+10000
#'
#' #For meaningful results, trial_reps needs to be larger than 100.
#' weibull_test2 <- historic_sim(trial_reps = 100, outcome_type = "weibull",
#'                               subj_per_arm = c(50, 100, 150, 200, 250),
#'                               a0_vals = c(0, 0.33, 0.67, 1),
#'                               effect_vals = c(0.6, 1, 1.4),
#'                               rand_control_diff = c(0.8, 1, 1.2),
#'                               hist_control_data = histdata, time_vec = NULL,
#'                               censor_value = 3, alpha = 0.05, get_var = TRUE,
#'                               get_bias = TRUE, get_mse = TRUE, seedval=123,
#'                               quietly=TRUE)
#'
#' #Tabulate the simulation results for power.
#' test_table <- print(x=weibull_test2, measure="power",
#'                     tab_type="WX|YZ", effect_val=0.6,
#'                     rand_control_diff_val=1.0)
#' print(test_table)
#'
#' #Tabulate the simulation results for estimates.
#' print(x=weibull_test2, measure="est", tab_type="WX|YZ",
#'       effect_val=0.6, rand_control_diff_val=1.0)
#'
#' #Tabulate the simulation results for variance.
#' print(x=weibull_test2, measure="var", tab_type="WX|YZ",
#'       effect_val=0.6, rand_control_diff_val=1.0)
#'
#' #Tabulate the simulation results for bias.
#' print(x=weibull_test2, measure="bias", tab_type="WX|YZ",
#'       effect_val=0.6, rand_control_diff_val=1.0)
#'
#' #Tabulate the simulation results for mse.
#' print(x=weibull_test2, measure="mse", tab_type="WX|YZ",
#'       effect_val=0.6, rand_control_diff_val=1.0)
#' }
#'
#' \donttest{
#' #Run a Bernoulli simulation, using historic_sim().
#' #Generate a sample of historical Bernoulli data for use in example.
#' set.seed(2250)
#' samplehistdata <- genbernoullidata(sample_size=60, prob1=0.6, odds_ratio=0.6)
#' histdata <- subset(samplehistdata, subset=(treatment==0))
#' histdata$id <- histdata$id+10000
#'
#' #For meaningful results, trial_reps needs to be larger than 100.
#' bernoulli_test <- historic_sim(trial_reps = 100, outcome_type = "bernoulli",
#'                               subj_per_arm = c(150),
#'                               a0_vals = c(1.0),
#'                               effect_vals = c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0),
#'                               rand_control_diff = c(1.8),
#'                               hist_control_data = histdata, time_vec = NULL,
#'                               censor_value = 3, alpha = 0.05, get_var = TRUE,
#'                               get_bias = TRUE, get_mse = TRUE, seedval=123,
#'                               quietly=TRUE)
#' test_table <- print(x=bernoulli_test, measure="power",
#'                     tab_type=NULL, effect_val=NULL,
#'                     subj_per_arm_val=NULL)
#' print(test_table)
#'
#' #If only one or two of the subj_per_arm, a0_vals, effect_vals, or
#' #rand_control_diff parameters have length greater than 1, then
#' #only bayes_ctd_array and measure parameters are needed.
#' #Tabulate the simulation results for estimates.
#' print(x=bernoulli_test, measure="est")
#'
#' #Tabulate the simulation results for variance.
#' print(x=bernoulli_test, measure="var")
#'
#' #Tabulate the simulation results for bias.
#' print(x=bernoulli_test, measure="bias")
#'
#' #Tabulate the simulation results for mse.
#' print(x=bernoulli_test, measure="mse")
#' }
#'
#' @export
print.bayes_ctd_array <- function(x = NULL, measure = "power", tab_type = "WX|YZ", subj_per_arm_val = NULL,
                                        a0_val = NULL, effect_val = NULL, rand_control_diff_val = NULL, print_chg_warn = 1, ...) {

  # (W): subj_per_arm,
  # (X): a0_val,
  # (Y): effect_val,
  # (Z): rand_control_diff
  # tab_type='WX|YZ', Table of Design Characteristics: Sample Size by a0
  # tab_type='WY|XZ', Table of Design Characteristics: Sample Size by Effect
  # tab_type='WZ|XY', Table of Design Characteristics: Sample Size by Control Differences
  # tab_type='XY|WZ', Table of Design Characteristics: a0 by Effect
  # tab_type='XZ|WY', Table of Design Characteristics: a0 by Control Diffferences
  # tab_type='YZ|WX', Table of Design Characteristics: Effect by Control Diffferences
  # tab_type='ZX|WY', Table of Design Characteristics: Control Diffferences by a0
  # tab_type='XW|YZ', Table of Design Characteristics: a0 by Sample Size
  # tab_type='YW|XZ', Table of Design Characteristics: Effect by Sample Size
  # tab_type='YX|WZ', Table of Design Characteristics: Effect by a0
  # tab_type='ZW|XY', Table of Design Characteristics: Control Differences by Sample Size
  # tab_type='ZY|WX', Table of Design Characteristics: Control Differences by Effect

  #Add simple code to help out the user with tab_type on matters where order is
  #not important and can be easily fixed.
  tab_type = toupper(tab_type)
  if (!((length(tab_type)== 0) && (typeof(tab_type) == "character"))){
    if (x$objtype == 'historic') {
      if (tab_type=='WX|ZY') tab_type <- 'WX|YZ'
      if (tab_type=='WY|ZX') tab_type <- 'WY|XZ'
      if (tab_type=='WZ|YX') tab_type <- 'WZ|XY'
      if (tab_type=='XY|ZW') tab_type <- 'XY|WZ'
      if (tab_type=='XZ|YW') tab_type <- 'XZ|WY'
      if (tab_type=='YZ|XW') tab_type <- 'YZ|WX'
      if (tab_type=='ZX|YW') tab_type <- 'ZX|WY'
      if (tab_type=='XW|ZY') tab_type <- 'XW|YZ'
      if (tab_type=='YW|ZX') tab_type <- 'YW|XZ'
      if (tab_type=='YX|ZW') tab_type <- 'YX|WZ'
      if (tab_type=='ZW|YX') tab_type <- 'ZW|XY'
      if (tab_type=='ZY|XW') tab_type <- 'ZY|WX'
    }
  }
  #Do the same for measure
  measure = tolower(measure)
  #The above code was added very late in initial production, so time was not
  # available to go through all code and remove code no longer necessary.


  #------------- Go through all the high level checks for proper input. -------------#

  if (x$objtype == 'simple') {
    tab_type <- "WY|XZ"
    a0_val <- 0
    rand_control_diff_val <- 1
    subj_per_arm_val <- x$subj_per_arm[1]
    effect_val <- x$effect_vals[1]
    if (print_chg_warn == 1) {
      print("Since simple_sim was used, tab_type was set to WY|XZ")
      print("Values for tab_type, subj_per_arm_val, a0_val, effect_val, and rand_control_diff_val were ignored")
      print("This works towards putting all results in a single table, effect by sample size")
    }
  }

  if (x$objtype == 'realsimple') {
    if (length(x$subj_per_arm) > 1) {
      if (print_chg_warn == 1) {
        print("Since only subj_per_arm vector has more than 1 element, tab_type was set to WX|YZ")
        print("This works towards putting all results in a single table")
      }
      tab_type <- "WX|YZ"
      a0_val <- x$a0_vals[1]
      effect_val <- x$effect_vals[1]
      rand_control_diff_val <- x$rand_control_diff[1]
    }
    if (length(x$a0_vals) > 1) {
      if (print_chg_warn == 1) {
        print("Since only a0_vals vector has more than 1 element, tab_type was set to XY|WZ")
        print("This works towards putting all results in a single table")
      }
      tab_type <- "XY|WZ"
      subj_per_arm_val <- x$subj_per_arm[1]
      effect_val <- x$effect_vals[1]
      rand_control_diff_val <- x$rand_control_diff[1]
    }
    if (length(x$effect_vals) > 1) {
      if (print_chg_warn == 1) {
        print("Since only effect_vals vector has more than 1 element, tab_type was set to YZ|WX")
        print("This works towards putting all results in a single table")
      }
      tab_type <- "YZ|WX"
      subj_per_arm_val <- x$subj_per_arm[1]
      a0_val <- x$a0_vals[1]
      rand_control_diff_val <- x$rand_control_diff[1]
    }
    if (length(x$rand_control_diff) > 1) {
      if (print_chg_warn == 1) {
        print("Since only rand_control_diff vector has more than 1 element, tab_type was set to ZX|WY")
        print("This works towards putting all results in a single table")
      }
      tab_type <- "ZX|WY"
      subj_per_arm_val <- x$subj_per_arm[1]
      a0_val <- x$a0_vals[1]
      effect_val <- x$effect_vals[1]
    }
  }

  # Check to see if two dimensions of bayes_ctd_array have only 1 element each. Is so set tab_type to an appropriate
  # value.
  if (x$objtype == 'historic') {
    if (length(x$effect_vals) == 1 & length(x$rand_control_diff) == 1) {
      chkflag1 <- 0
      if(!((length(tab_type)== 0) && (typeof(tab_type) == "character"))){
        if ((tab_type != "WX|YZ") & (tab_type != "XW|YZ")) {
          tab_type <- "WX|YZ"
          chkflag1 <- 1
        }
      }
      if( (length(tab_type)== 0) && (typeof(tab_type) == "character") ){
        tab_type <- "WX|YZ"
        chkflag1 <- 1
      }
      effect_val <- x$effect_vals[1]
      rand_control_diff_val <- x$rand_control_diff[1]
      if (print_chg_warn == 1 & chkflag1 == 1) {
        print("Since effect_vals and rand_control_diff vectors only have 1 element each and tab_type not equal to 'WX|YZ' or 'XW|YZ', tab_type was set to WX|YZ")
        print("This works towards putting all results in a single table")
      }
    }
    if (length(x$a0_vals) == 1 & length(x$rand_control_diff) == 1) {
      chkflag1 <- 0
      if(!((length(tab_type)== 0) && (typeof(tab_type) == "character"))){
        if ((tab_type != "WY|XZ") & (tab_type != "YW|XZ")) {
          tab_type <- "WY|XZ"
          chkflag1 <- 1
        }
      }
      if((length(tab_type)== 0) && (typeof(tab_type) == "character")){
        tab_type <- "WY|XZ"
        chkflag1 <- 1
      }
      a0_val <- x$a0_vals[1]
      rand_control_diff_val <- x$rand_control_diff[1]
      if (print_chg_warn == 1 & chkflag1 == 1) {
        print("Since a0 and rand_control_diff vectors only have 1 element each and tab_type not equal to 'WY|XZ' or 'YW|XZ', tab_type was set to WY|XZ")
        print("This works towards putting all results in a single table")
      }
    }
    if (length(x$subj_per_arm) == 1 & length(x$rand_control_diff) == 1) {
      chkflag1 <- 0
      if(!((length(tab_type)== 0) && (typeof(tab_type) == "character"))){
        if ((tab_type != "XY|WZ") & (tab_type != "YX|WZ")) {
          tab_type <- "XY|WZ"
          chkflag1 <- 1
        }
      }
      if((length(tab_type)== 0) && (typeof(tab_type) == "character")){
        tab_type <- "XY|WZ"
        chkflag1 <- 1
      }
      subj_per_arm_val <- x$subj_per_arm
      rand_control_diff_val <- x$rand_control_diff
      if (print_chg_warn == 1 & chkflag1 == 1) {
        print("Since sample size and rand_control_diff vectors only have 1 element each and tab_type not equal to 'XY|WZ' or 'YX|WZ', tab_type was set to XY|WZ")
        print("This works towards putting all results in a single table")
      }
    }
    if (length(x$subj_per_arm) == 1 & length(x$effect_vals) == 1) {
      chkflag1 <- 0
      if(!((length(tab_type)== 0) && (typeof(tab_type) == "character"))){
        if ((tab_type != "XZ|WY") & (tab_type != "ZX|WY")) {
          tab_type <- "XZ|WY"
          chkflag1 <- 1
        }
      }
      if((length(tab_type)== 0) && (typeof(tab_type) == "character")){
        tab_type <- "XZ|WY"
        chkflag1 <- 1
      }
      subj_per_arm_val <- x$subj_per_arm
      effect_val <- x$effect_vals
      if (print_chg_warn == 1 & chkflag1 == 1) {
        print("Since sample size and effect_vals vectors only have 1 element each and tab_type not equal to 'XZ|WY' or 'ZX|WY', tab_type was set to XZ|WY")
        print("This works towards putting all results in a single table")
      }
    }
    if (length(x$subj_per_arm) == 1 & length(x$a0_vals) == 1) {
      chkflag1 <- 0
      if(!((length(tab_type)== 0) && (typeof(tab_type) == "character"))){
        if ((tab_type != "YZ|WX") & (tab_type != "ZY|WX")) {
          tab_type <- "YZ|WX"
          chkflag1 <- 1
        }
      }
      if((length(tab_type)== 0) && (typeof(tab_type) == "character")){
        tab_type <- "YZ|WX"
        chkflag1 <- 1
      }
      subj_per_arm_val <- x$subj_per_arm
      a0_val <- x$a0_vals
      if (print_chg_warn == 1 & chkflag1 == 1) {
        print("Since sample size and a0 vectors only have 1 element each and tab_type not equal to 'YZ|WX' or 'ZY|WX', tab_type was set to YZ|WX")
        print("This works towards putting all results in a single table")
      }
    }
    if (length(x$a0_vals) == 1 & length(x$effect_vals) == 1) {
      chkflag1 <- 0
      if(!((length(tab_type)== 0) && (typeof(tab_type) == "character"))){
        if ((tab_type != "WZ|XY") & (tab_type != "ZW|XY")) {
          tab_type <- "WZ|XY"
          chkflag1 <- 1
        }
      }
      if((length(tab_type)== 0) && (typeof(tab_type) == "character")){
        tab_type <- "WZ|XY"
        chkflag1 <- 1
      }
      a0_val <- x$a0_vals
      effect_val <- x$effect_vals
      if (print_chg_warn == 1 & chkflag1 == 1) {
        print("Since a0 and effect_vals vectors only have 1 element each and tab_type not equal to 'WZ|XY' or 'ZW|XY', tab_type was set to WZ|XY")
        print("This works towards putting all results in a single table")
      }
    }
  }

  print_error_checks(x, measure, tab_type, subj_per_arm_val, a0_val, effect_val, rand_control_diff_val)

  #Select the data array depending on user defined value of measure.
  if (tolower(measure) == "power")
    dataarray <- x$data$power
  if (tolower(measure) == "est")
    dataarray <- x$data$est
  if (tolower(measure) == "var")
    dataarray <- x$data$var
  if (tolower(measure) == "bias")
    dataarray <- x$data$bias
  if (tolower(measure) == "mse")
    dataarray <- x$data$mse
  if (toupper(tab_type) == "WX|YZ") {
    # Table of Design Characteristics: Sample Size by a0
    effect_val_levels <- x$effect_vals
    rand_control_diff_levels <- x$rand_control_diff
    # check to make sure effect_val is in effect_vals and
    # rand_control_diff_val is in rand_control_diff
    effect_val_loc <- match(effect_val, effect_val_levels)
    rand_control_diff_loc <- match(rand_control_diff_val, rand_control_diff_levels)
    return_table <- dataarray[,, effect_val_loc, rand_control_diff_loc]
  }
  if (toupper(tab_type) == "WY|XZ") {
    # Table of Design Characteristics: Sample Size by Effect
    a0_val_levels <- x$a0_vals
    rand_control_diff_levels <- x$rand_control_diff
    # check to make sure a0_val is in a0_vals and
    # rand_control_diff_val is in rand_control_diff
    a0_val_loc <- match(a0_val, a0_val_levels)
    rand_cntrl_diff_loc <- match(rand_control_diff_val, rand_control_diff_levels)
    return_table <- dataarray[, a0_val_loc,, rand_cntrl_diff_loc]
  }
  if (toupper(tab_type) == "WZ|XY") {
    # Table of Design Characteristics: Sample Size by Control Differences
    a0_val_levels <- x$a0_vals
    effect_val_levels <- x$effect_vals
    # check to make sure effect_val is in effect_vals and
    # a0_val is in a0_vals
    a0_val_loc <- match(a0_val, a0_val_levels)
    effect_val_loc <- match(effect_val, effect_val_levels)
    return_table <- dataarray[, a0_val_loc, effect_val_loc, ]
  }
  if (toupper(tab_type) == "XY|WZ") {
    # Table of Design Characteristics: a0 by Effect
    subj_per_arm_levels <- x$subj_per_arm
    rand_control_diff_levels <- x$rand_control_diff
    # check to make sure subj_per_arm_val is in subj_per_arm and
    # rand_control_diff_val is in rand_control_diff
    subj_per_arm_loc <- match(subj_per_arm_val, subj_per_arm_levels)
    rand_control_diff_loc <- match(rand_control_diff_val, rand_control_diff_levels)
    return_table <- dataarray[subj_per_arm_loc,,, rand_control_diff_loc]
  }
  if (toupper(tab_type) == "XZ|WY") {
    # Table of Design Characteristics: a0 by Control Diffferences
    subj_per_arm_levels <- x$subj_per_arm
    effect_val_levels <- x$effect_vals
    # check to make sure effect_val is in effect_vals and
    # subj_per_arm_val is in subj_per_arm
    subj_per_arm_loc <- match(subj_per_arm_val, subj_per_arm_levels)
    effect_val_loc <- match(effect_val, effect_val_levels)
    return_table <- dataarray[subj_per_arm_loc,, effect_val_loc, ]
  }
  if (toupper(tab_type) == "YZ|WX") {
    # Table of Design Characteristics: Effect by Control Diffferences
    subj_per_arm_levels <- x$subj_per_arm
    a0_val_levels <- x$a0_vals
    # check to make sure subj_per_arm_val is in subj_per_arm and
    # a0_val is in a0_vals
    subj_per_arm_loc <- match(subj_per_arm_val, subj_per_arm_levels)
    a0_val_loc <- match(a0_val, a0_val_levels)
    return_table <- dataarray[subj_per_arm_loc, a0_val_loc,, ]
  }
  if (toupper(tab_type) == "XW|YZ") {
    # Table of Design Characteristics: a0 by sample size
    effect_val_levels <- x$effect_vals
    rand_control_diff_levels <- x$rand_control_diff
    # check to make sure rand_control_diff_val is in rand_control_diff and
    # effect_val is in effect_vals
    effect_val_loc <- match(effect_val, effect_val_levels)
    rand_control_diff_loc <- match(rand_control_diff_val, rand_control_diff_levels)
    return_table <- t(dataarray[,,effect_val_loc,rand_control_diff_loc])
  }
  if (toupper(tab_type) == "YW|XZ") {
    # Table of Design Characteristics: Effect by sample size
    a0_val_levels <- x$a0_vals
    rand_control_diff_levels <- x$rand_control_diff
    # check to make sure rand_control_diff_val is in rand_control_diff and
    # a0_val is in a0_vals
    a0_val_loc <- match(a0_val, a0_val_levels)
    rand_control_parm_loc <- match(rand_control_diff_val, rand_control_diff_levels)
    return_table <- t(dataarray[,a0_val_loc,,rand_control_parm_loc])
  }
  if (toupper(tab_type) == "YX|WZ") {
    # Table of Design Characteristics: Effect by a0
    subj_per_arm_levels <- x$subj_per_arm
    rand_control_diff_levels <- x$rand_control_diff
    # check to make sure subj_per_arm_val is in subj_per_arm and
    # rand_control_diff_val is in rand_control_diff
    subj_per_arm_loc <- match(subj_per_arm_val, subj_per_arm_levels)
    rand_control_parm_loc <- match(rand_control_diff_val, rand_control_diff_levels)
    return_table <- dataarray[subj_per_arm_loc,,,rand_control_parm_loc]
  }
  if (toupper(tab_type) == "ZW|XY") {
    # Table of Design Characteristics: Control Differences by sample size
    a0_val_levels <- x$a0_vals
    effect_val_levels <- x$effect_vals
    # check to make sure a0_val is in a0_vals and
    # effect_val is in effect_vals
    a0_val_loc <- match(a0_val, a0_val_levels)
    effect_val_loc <- match(effect_val, effect_val_levels)
    return_table <- t(dataarray[,a0_val_loc,effect_val_loc, ])
  }
  if (toupper(tab_type) == "ZX|WY") {
    # Table of Design Characteristics: Control Differences by a0
    subj_per_arm_levels <- x$subj_per_arm
    effect_val_levels <- x$effect_vals
    # check to make sure subj_per_arm_val is in subj_per_arm and
    # effect_val is in effect_vals
    subj_per_arm_loc <- match(subj_per_arm_val, subj_per_arm_levels)
    effect_val_loc <- match(effect_val, effect_val_levels)
    return_table <- t(dataarray[subj_per_arm_loc,,effect_val_loc, ])
  }
  if (toupper(tab_type) == "ZY|WX") {
    # Table of Design Characteristics: Control Differences by Effect
    subj_per_arm_levels <- x$subj_per_arm
    a0_val_levels <- x$a0_vals
    # check to make sure subj_per_arm_val is in subj_per_arm and
    # a0_val is in a0_vals
    subj_per_arm_loc <- match(subj_per_arm_val, subj_per_arm_levels)
    a0_val_loc <- match(a0_val, a0_val_levels)
    return_table <- t(dataarray[subj_per_arm_loc,a0_val_loc,, ])
  }

  if (toupper(tab_type) == "YX|WZ") {
    # Table of Design Characteristics: Effect by a0
    subj_per_arm_levels <- x$subj_per_arm
    rand_control_diff_levels <- x$rand_control_diff
    # check to make sure subj_per_arm_val is in subj_per_arm and
    # rand_control_diff_val is in rand_control_diff
    subj_per_arm_loc <- match(subj_per_arm_val, subj_per_arm_levels)
    rand_control_diff_val_loc <- match(rand_control_diff_val, rand_control_diff_levels)
    return_table <- t(dataarray[subj_per_arm_loc,,,rand_control_diff_val_loc ])
  }

  return_table
}



#' Plot Data from Two Arm Bayesian Clinical Trial Simulation.
#'
#' \code{plot.bayes_ctd_array()} takes an S3 object of class \code{bayes_ctd_array}, and
#' creates a line plot from a one or two dimensional slice of the data generated by a
#' clinical trial simulation using \code{historic_sim()} or \code{simple_sim()}.  The
#' plotted results can be smoothed or unsmoothed.
#'
#' If the object of class \code{bayes_ctd_array} is created by \code{historic_sim()},
#' the function \code{plot()} allows the user to create line plots of user-specified
#' 1- or 2- dimensional slices of the simulation results based on slicing code
#' described below. If the object of class \code{bayes_ctd_array} is created by
#' \code{simple_sim()}, a basic plot of characteristic by sample size and effect is created.
#'
#' If the object of class \code{bayes_ctd_array} is created by \code{simple_sim()}, then
#' all four trial characteristics (\code{subj_per_arm_val}, \code{a0_vals},
#' \code{effect_val}, and \code{rand_control_diff_val}) can be ignored as can the
#' parameter defining what type of plot to create through the parameter \code{tab_type}.
#' A call to \code{plot()} will require the user to specify a measure (power, est,
#' var, bias, or mse).
#'
#' If the object of class \code{bayes_ctd_array} is created by \code{historic_sim()},
#' when calling \code{plot()} the user must specify a measure to plot
#' (power, est, var, bias, or mse) and may be required to specify a plot type through
#' the \code{tab_type} parameter.  A plot type, \code{tab_type}, will be required if
#' 3 of the 4 trial characteristics are equal to a vector of 2 or more values.  This
#' plot type specification uses the letters W, X, Y, and Z.  The letter W represents
#' the subject per arm dimension.  The letter X represents the a0 dimension.  The
#' letter Y represents the effect dimension.  The letter Z represents the control
#' difference dimension.  To plot a slice of the 4-dimensional array, these letters
#' are put into an AB|CD pattern just like in \code{print()}.  The two letters
#' to the right of the vertical bar define which variables are held constant.  The two
#' letters to the left of the vertical bar define which variables are going to show up
#' in the plot.  The first letter defines the x-axis variable and the second letter
#' defines the stratification variable.  The result is a plot of power, estimate,
#' variance, bias, or mse by the trial characteristic represented by the first letter.
#' On this plot, one line will be created for each value of the trial characteristic
#' represented by the second letter.  For example if tab_type equals \code{WX|YZ},
#' then effect and control differences will be held constant, while sample size will be
#' represented along the horizontal axis and a0 values will be represented by separate
#' lines.  The actual values that are plotted on the y-axis depend on what measure is
#' requested in the parameter \code{measure}.
#'
#' \itemize{
#'   \item \code{tab_type='WX|YZ'}, Sample Size by a0
#'   \item \code{tab_type='WY|XZ'}, Sample Size by Effect
#'   \item \code{tab_type='WZ|XY'}, Sample Size by Control Differences
#'   \item \code{tab_type='XY|WZ'}, a0 by Effect
#'   \item \code{tab_type='XZ|WY'}, a0 by Control Differences
#'   \item \code{tab_type='YZ|WX'}, Effect by Control Differences
#'   \item \code{tab_type='ZX|WY'}, Control Differences by a0
#'   \item \code{tab_type='XW|YZ'}, a0 by Sample Size
#'   \item \code{tab_type='YW|XZ'}, Effect by Sample Size
#'   \item \code{tab_type='YX|WZ'}, Effect by a0
#'   \item \code{tab_type='ZW|XY'}, Control Differences by Sample Size
#'   \item \code{tab_type='ZY|WX'}, Control Differences by Effect
#' }
#' It is very important to populate the values of \code{subj_per_arm_val},
#' \code{a0_val}, \code{effect_val}, and \code{rand_control_diff_val} correctly given
#' the value of tab_type, when the object of class \code{bayes_ctd_array} is created by
#' \code{historic_sim()} and at least 3 of the four parameters have more than one
#' value.  On, the other hand, if 2 or more of the four parameters have only one value,
#' then \code{subj_per_arm_val}, \code{a0_vals}, \code{effect_val},
#' \code{rand_control_diff_val}, as well as \code{tab_type} can be ignored.  If the last
#' two letters are \code{YZ}, then \code{effect_val} and \code{rand_control_diff_val}
#' must be populated.  If the last two letters are \code{XZ}, then \code{a0_val} and
#' \code{rand_control_diff_val} must be populated.  If the last two letters are \code{XY},
#' then \code{a0_val} and \code{effect_val} must be populated.  If the last two letters
#' are \code{WZ}, then \code{sample_val} and \code{rand_control_diff_val} must be
#' populated.  If the last two letters are \code{WY}, then \code{sample_size_val} and
#' \code{effect_val} must be populated.  If the last two letters are \code{WX}, then
#' \code{sample_size_val} and \code{a0_val} must be populated.
#'
#' If the object of class \code{bayes_ctd_array} is created by \code{simple_sim()}, the
#' parameters \code{tab_type}, \code{subj_per_arm_val}, \code{a0_val}, \code{effect_val},
#' and \code{rand_control_diff_val} are ignored.
#'
#' @param x Name of object of class \code{bayes_ctd_array} containing
#'   data from clinical trial simulation.
#' @param measure Must be equal to \code{power}, \code{est}, \code{var}, \code{bias},
#'   or \code{mse}.  Default is \code{power}.  Case does not matter.
#' @param tab_type A character string that must equal \code{WX|YZ}, \code{WY|XZ},
#'   \code{WZ|XY}, \code{XY|WZ}, \code{XZ|WY}, \code{YZ|WX}, \code{ZX|WY}, \code{XW|YZ},
#'   \code{YW|XZ}, \code{YX|WZ}, \code{ZW|XY}, \code{ZX|WY}, \code{ZY|WX} when
#'   \code{x} is generated by \code{historic_sim()}.  Default is
#'   \code{WX|YZ}.  When \code{x} is generated by \code{simple_sim()},
#'   \code{tab_type} is ignored.
#' @param smooth A true/false parameter indicating whether smoothed results
#'   should be plotted. Note, smoothing of simulation results requires the length of
#'   \code{subj_per_arm_val} or \code{a0_val} or \code{effect_val} or
#'   \code{rand_control_diff_val}, whichever populates the x-axis on the graph to
#'   contain enough elements to justify the smoothing.  No checking occurs to
#'   determine if enough elements are present to justify smoothing.  Default is
#'   \code{FALSE}.
#' @param plot_out A true/false parameter indicating whether the plot should be
#'   produced.  This parameter is useful if the user only wants a table of smoothed
#'   values.  Default is \code{TRUE}.
#' @param subj_per_arm_val Must be non-missing, if \code{x} is generated
#'   by \code{historic_sim()} and sample size is being held constant.
#'   If \code{x} is generated by \code{historic_sim()} and sample size
#'   is being held constant, \code{subj_per_arm_val} must equal a value submitted
#'   to \code{historic_sim()} within the \code{subj_per_arm} parameter.  When
#'   \code{x} is generated by \code{simple_sim()}, \code{subj_per_arm_val}
#'   is ignored.
#' @param a0_val Must be non-missing, if \code{x} is generated
#'   by \code{historic_sim()} and a0, the power prior parameter, is being held
#'   constant.  If \code{x} is generated by \code{historic_sim()} and
#'   a0 is being held constant, \code{a0_val} must equal a value submitted
#'   to \code{historic_sim()} within the \code{a0_val} parameter.  When
#'   \code{x} is generated by \code{simple_sim()}, \code{a0_val} is
#'   ignored.
#' @param effect_val  Must be non-missing, if \code{x} is generated
#'   by \code{historic_sim()} and effect is being held constant.  If
#'   \code{x} is generated by \code{historic_sim()} and effect is being
#'   held constant, \code{effect_val} must equal a value submitted to
#'   \code{historic_sim()} within the \code{effect_vals} parameter.  When
#'   \code{x} is generated by \code{simple_sim()}, \code{effect_val} is
#'   ignored.
#' @param rand_control_diff_val Must be non-missing, if \code{x} is
#'   generated by \code{historic_sim()} and differences between randomized
#'   and historical controls are being held constant.  If \code{x}
#'   is generated by \code{historic_sim()} and control differences are being
#'   held constant, \code{rand_control_diff_val} must equal a value submitted to
#'   \code{historic_sim()} within the \code{rand_control_diff} parameter.  When
#'   \code{x} is generated by \code{simple_sim()},
#'   \code{rand_control_diff_val} is ignored.
#' @param span The \code{span} parameter value for a call \code{loess()}.  Default is 0.75.  If
#'   \code{span} is a single number, then that value will be used to smooth the data in all
#'   columns of the table being plotted.  If \code{span} is a vector, then it must have length
#'   equal to the number of columns being plotted.
#' @param degree The \code{degree} parameter value for a call \code{loess()}.  Default is 2.
#'   The value of \code{degree} will be used for all columns being plotted.
#' @param family The \code{family} parameter value for a call \code{loess()}.  Default is
#'   "\code{gaussian}".  The value of \code{family} will be used for all columns being plotted.
#' @param title Title for the plot.
#' @param ylim Lower and upper limits for y-axis of plot.
#' @param ...	further arguments passed to or from other methods.
#'
#' @return \code{plot()} returns a plot for a two dimensional array of simulation
#'   results.  If \code{smooth} is \code{TRUE}, then the plot is based on a smoothed
#'   version of the simulation results. If \code{smooth} is \code{FALSE}, then the plot
#'   is based on the raw data from the simulation results.  What actually is plotted
#'   depends on the value of \code{measure}.  If \code{plot_out} is \code{FALSE}, the
#'   plot is not created.  This option is useful when the user wants a table of smoothed
#'   simulation results but does not want the plot. Smoothing of simulation results
#'   requires the length of \code{subj_per_arm_val} or \code{a0_val} or \code{effect_val}
#'   or \code{rand_control_diff_val}, whichever populates the x-axis on the graph to
#'   contain enough elements to justify the smoothing.  No checking occurs to
#'   determine if enough elements are present to justify smoothing.
#'
#' @references Eggleston et al. (2021). BayesCTDesign: An R Package for Bayesian Trial Design Using Historical Control Data. \emph{Journal of Statistical Software}, November 2021, Volume 100, Issue 21,  <doi:10.18637/jss.v100.i21>.
#'
#' @examples
#' #Run a Weibull simulation, using simple_sim().
#' #For meaningful results, trial_reps needs to be much larger than 2.
#' weibull_test <- simple_sim(trial_reps = 2, outcome_type = "weibull",
#'                            subj_per_arm = c(50, 100, 150, 200),
#'                            effect_vals = c(0.6, 1),
#'                            control_parms = c(2.82487,3), time_vec = NULL,
#'                            censor_value = NULL, alpha = 0.05,
#'                            get_var = TRUE, get_bias = TRUE, get_mse = TRUE,
#'                            seedval=123, quietly=TRUE)
#'
#' #Create a plot of the power simulation results.
#' plot(x=weibull_test, measure="power", tab_type=NULL,
#'      smooth=FALSE, plot_out=TRUE, subj_per_arm_val=NULL, a0_val=NULL,
#'      effect_val=NULL, rand_control_diff_val=NULL)
#' #Create a plot of the hazard ratio simulation results.
#' plot(x=weibull_test, measure="est", tab_type=NULL,
#'      smooth=FALSE, plot_out=TRUE, subj_per_arm_val=NULL, a0_val=NULL,
#'      effect_val=NULL, rand_control_diff_val=NULL)
#' #Create a plot of the hazard ratio variance simulation results.
#' plot(x=weibull_test, measure="var", tab_type=NULL,
#'      smooth=FALSE, plot_out=TRUE, subj_per_arm_val=NULL, a0_val=NULL,
#'      effect_val=NULL, rand_control_diff_val=NULL)
#' #Create a plot of the hazard ratio bias simulation results.
#' plot(x=weibull_test, measure="bias", tab_type=NULL,
#'      smooth=FALSE, plot_out=TRUE, subj_per_arm_val=NULL, a0_val=NULL,
#'      effect_val=NULL, rand_control_diff_val=NULL)
#' #Create a plot of the hazard ratio mse simulation results.
#' plot(x=weibull_test, measure="mse", tab_type=NULL,
#'      smooth=FALSE, plot_out=TRUE, subj_per_arm_val=NULL, a0_val=NULL,
#'      effect_val=NULL, rand_control_diff_val=NULL)
#'
#' \donttest{
#' #Run a second Weibull simulation, using simple_sim() and smooth the plot.
#' #For meaningful results, trial_reps needs to be larger than 100.
#' weibull_test2 <- simple_sim(trial_reps = 100, outcome_type = "weibull",
#'                             subj_per_arm = c(50, 75, 100, 125, 150, 175, 200, 225, 250),
#'                             effect_vals = c(0.6, 1, 1.4),
#'                             control_parms = c(2.82487,3), time_vec = NULL,
#'                             censor_value = NULL, alpha = 0.05, get_var = TRUE,
#'                             get_bias = TRUE, get_mse = TRUE, seedval=123,
#'                             quietly=TRUE)
#'
#' #Create a plot of the power simulation results.
#' plot(x=weibull_test2, measure="power", tab_type=NULL,
#'      smooth=TRUE, plot_out=TRUE, subj_per_arm_val=NULL, a0_val=NULL,
#'      effect_val=NULL, rand_control_diff_val=NULL, span=c(1,1,1))
#' }
#'
#' \donttest{
#' #Run a third weibull simulation, using historic_sim().
#' #Note: historic_sim() can take a while to run.
#' #Generate a sample of historical data for use in example.
#' set.seed(2250)
#' SampleHistData <- genweibulldata(sample_size=60, scale1=2.82487,
#'                                  hazard_ratio=0.6, common_shape=3,
#'                                  censor_value=3)
#' histdata <- subset(SampleHistData, subset=(treatment==0))
#' histdata$id <- histdata$id+10000
#'
#' #For meaningful results, trial_reps needs to be larger than 100.
#' weibull_test3 <- historic_sim(trial_reps = 100, outcome_type = "weibull",
#'                               subj_per_arm = c(50, 100, 150, 200, 250),
#'                               a0_vals = c(0, 0.33, 0.67, 1),
#'                               effect_vals = c(0.6, 1, 1.4),
#'                               rand_control_diff = c(0.8, 1, 1.2),
#'                               hist_control_data = histdata, time_vec = NULL,
#'                               censor_value = 3, alpha = 0.05, get_var = TRUE,
#'                               get_bias = TRUE, get_mse = TRUE, seedval=123,
#'                               quietly=TRUE)
#'
#' #Create a plot of the power simulation results.
#' plot(x=weibull_test3, measure="power", tab_type="WX|YZ",
#'      smooth=FALSE, plot_out=TRUE, effect_val=0.6,
#'      rand_control_diff_val=1.0)
#' }
#'
#' \donttest{
#' #Run a Gaussian simulation, using historic_sim()
#' #Generate a sample of historical Gaussian data for use in example.
#' set.seed(2250)
#' samplehistdata <- gengaussiandata(sample_size=60, mu1=25, mean_diff=0, common_sd=3)
#' histdata <- subset(samplehistdata, subset=(treatment==0))
#' histdata$id <- histdata$id+10000
#'
#' #For meaningful results, trial_reps needs to be larger than 100.
#' gaussian_test <- historic_sim(trial_reps = 100, outcome_type = "gaussian",
#'                              subj_per_arm = c(150),
#'                              a0_vals = c(1.0),
#'                              effect_vals = c(0.15),
#'                              rand_control_diff = c(-4.0,-3.5,-3.0,-2.5,-2.0,
#'                                                    -1.5,-1.0,-0.5,0,0.5,1.0),
#'                              hist_control_data = histdata, time_vec = NULL,
#'                              censor_value = 3, alpha = 0.05, get_var = TRUE,
#'                              get_bias = TRUE, get_mse = TRUE, seedval=123,
#'                              quietly=TRUE)
#' test_table <- print(x=gaussian_test, measure="power",
#'                          tab_type=NULL, effect_val=NULL,
#'                          subj_per_arm_val=NULL)
#' print(test_table)
#' #Create a plot of the power simulation results.
#' plot(x=gaussian_test, measure="power", tab_type=NULL,
#'      smooth=TRUE, plot_out=TRUE, effect_val=NULL,
#'      rand_control_diff_val=NULL)
#' }
#'
#' \donttest{
#' #Generate a sample of historical pwe data for use in example.
#' set.seed(2250)
#' nvalHC <- 60
#' time.vec <- c(0.3,0.9,1.5,2.1,2.4)
#' lambdaHC.vec <- c(0.19,0.35,0.56,0.47,0.38,0.34)
#' censor.value <- 3
#'
#' SampleHistData <- genpwedata(nvalHC, lambdaHC.vec, 1.0, time.vec, censor.value)
#' histdata <- subset(SampleHistData, subset=(treatment==0))
#' histdata$indicator <- 2 #If set to 2, then historical controls will be collapsed with
#' #randomized controls, when time_vec is re-considered and
#' #potentially restructured.  If set to 1, then historical
#' #controls will be treated as a separated cohort when
#' #time_vec is being assessed for restructuring.
#' histdata$id <- histdata$id+10000
#'
#' #Run a pwe simulation, using historic_sim().
#' #For meaningful results, trial_reps needs to be larger than 100.
#' pwe_test <- historic_sim(trial_reps = 100, outcome_type = "pwe",
#'                         subj_per_arm = c(25,50,75,100,125,150,175,200,225,250),
#'                         a0_vals = c(1.0),
#'                         effect_vals = c(0.6),
#'                         rand_control_diff = c(1.8),
#'                         hist_control_data = histdata, time_vec = time.vec,
#'                         censor_value = 3, alpha = 0.05, get_var = TRUE,
#'                         get_bias = TRUE, get_mse = TRUE, seedval=123,
#'                         quietly=TRUE)
#'
#' #Create a plot of the power simulation results.
#' plot(x=pwe_test, measure="power", tab_type=NULL,
#'      smooth=TRUE, plot_out=TRUE, effect_val=NULL,
#'      rand_control_diff_val=NULL)
#' }
#'
#' @export
plot.bayes_ctd_array <- function(x = NULL, measure = "power", tab_type = "WX|YZ", smooth = FALSE,
                                       plot_out = TRUE, subj_per_arm_val = NULL, a0_val = NULL, effect_val = NULL, rand_control_diff_val = NULL, span = 0.75,
                                       degree = 2, family = "gaussian", title=NULL, ylim=NULL, ...) {
  # (W): subj_per_arm,
  # (X): a0_val,
  # (Y): effect_val,
  # (Z): rand_control_diff
  # tab_type='WX|YZ', Line Plot of Design Characteristics: Sample Size by a0
  # tab_type='WY|XZ', Line Plot of Design Characteristics: Sample Size by Effect
  # tab_type='WZ|XY', Line Plot of Design Characteristics: Sample Size by Control Differences
  # tab_type='XY|WZ', Line Plot of Design Characteristics: a0 by Effect
  # tab_type='XZ|WY', Line Plot of Design Characteristics: a0 by Control Diffferences
  # tab_type='YZ|WX', Line Plot of Design Characteristics: Effect by Control Diffferences
  # tab_type='ZX|WY', Line Plot of Design Characteristics: Control Diffferences by a0
  # tab_type='XW|YZ', Line Plot of Design Characteristics: a0 by Sample Size
  # tab_type='YW|XZ', Line Plot of Design Characteristics: Effect by Sample Size
  # tab_type='YX|WZ', Line Plot of Design Characteristics: Effect by a0
  # tab_type='ZW|XY', Line Plot of Design Characteristics: Control Differences by Sample Size
  # tab_type='ZY|WX', Line Plot of Design Characteristics: Control Differences by Effect

  #Add simple code to help out the user with tab_type on matters where order is
  #not important and can be easily fixed.
  tab_type = toupper(tab_type)
  if (!((length(tab_type)== 0) && (typeof(tab_type) == "character"))){
    if (x$objtype == 'historic') {
      if (tab_type=='WX|ZY') tab_type <- 'WX|YZ'
      if (tab_type=='WY|ZX') tab_type <- 'WY|XZ'
      if (tab_type=='WZ|YX') tab_type <- 'WZ|XY'
      if (tab_type=='XY|ZW') tab_type <- 'XY|WZ'
      if (tab_type=='XZ|YW') tab_type <- 'XZ|WY'
      if (tab_type=='YZ|XW') tab_type <- 'YZ|WX'
      if (tab_type=='ZX|YW') tab_type <- 'ZX|WY'
      if (tab_type=='XW|ZY') tab_type <- 'XW|YZ'
      if (tab_type=='YW|ZX') tab_type <- 'YW|XZ'
      if (tab_type=='YX|ZW') tab_type <- 'YX|WZ'
      if (tab_type=='ZW|YX') tab_type <- 'ZW|XY'
      if (tab_type=='ZY|XW') tab_type <- 'ZY|WX'
    }
  }
  #Do the same for measure
  measure = tolower(measure)
  #The above code was added very late in initial production, so time was not
  # available to go through all code and remove code no longer necessary.

  #------------- Go through all the high level checks for proper input. -------------#
  if (x$objtype == 'simple') {
    tab_type <- "WY|XZ"
    a0_val <- 0
    rand_control_diff_val <- 1
    subj_per_arm_val <- x$subj_per_arm[1]
    effect_val <- x$effect_vals[1]
    print("Since simple_sim was used, tab_type was set to WY|XZ")
    print("Values for tab_type, subj_per_arm_val, a0_val, effect_val, and rand_control_diff_val were ignored")
    print("This works towards putting all results on a single plot, effect by sample size")
  }

  if (x$objtype == 'realsimple') {
    if (length(x$subj_per_arm) > 1) {
      print("Since only subj_per_arm vector has more than 1 element, tab_type was set to WX|YZ")
      print("This works towards putting all results in a single table")
      tab_type <- "WX|YZ"
      a0_val <- x$a0_vals[1]
      effect_val <- x$effect_vals[1]
      rand_control_diff_val <- x$rand_control_diff[1]
    }
    if (length(x$a0_vals) > 1) {
      print("Since only a0_vals vector has more than 1 element, tab_type was set to XY|WZ")
      print("This works towards putting all results in a single table")
      tab_type <- "XY|WZ"
      subj_per_arm_val <- x$subj_per_arm[1]
      effect_val <- x$effect_vals[1]
      rand_control_diff_val <- x$rand_control_diff[1]
    }
    if (length(x$effect_vals) > 1) {
      print("Since only effect_vals vector has more than 1 element, tab_type was set to YZ|WX")
      print("This works towards putting all results in a single table")
      tab_type <- "YZ|WX"
      subj_per_arm_val <- x$subj_per_arm[1]
      a0_val <- x$a0_vals[1]
      rand_control_diff_val <- x$rand_control_diff[1]
    }
    if (length(x$rand_control_diff) > 1) {
      print("Since only rand_control_diff vector has more than 1 element, tab_type was set to ZX|WY")
      print("This works towards putting all results in a single table")
      tab_type <- "ZX|WY"
      subj_per_arm_val <- x$subj_per_arm[1]
      a0_val <- x$a0_vals[1]
      effect_val <- x$effect_vals[1]
    }
  }

  if (x$objtype == 'historic') {
    # Check to see if two dimensions of x have only 1 element each. Is so set tab_type to an appropriate
    # value.
    if (length(x$effect_vals) == 1 & length(x$rand_control_diff) == 1) {
      chkflag1 <- 0
      if(!((length(tab_type)== 0) && (typeof(tab_type) == "character"))){
        if ((tab_type != "WX|YZ") & (tab_type != "XW|YZ")) {
          tab_type <- "WX|YZ"
          chkflag1 <- 1
        }
      }
      if( (length(tab_type)== 0) && (typeof(tab_type) == "character")){
        tab_type <- "WX|YZ"
        chkflag1 <- 1
      }
      effect_val <- x$effect_vals[1]
      rand_control_diff_val <- x$rand_control_diff[1]
      if (chkflag1 == 1) {
        print("Since effect_vals and rand_control_diff vectors only have 1 element each and tab_type not equal to 'WX|YZ' or 'XW|YZ', tab_type was set to WX|YZ")
        print("This works towards putting all results on a single plot")
      }
    }
    if (length(x$a0_vals) == 1 & length(x$rand_control_diff) == 1) {
      chkflag1 <- 0
      if(!((length(tab_type)== 0) && (typeof(tab_type) == "character"))){
        if ((tab_type != "WY|XZ") & (tab_type != "YW|XZ")) {
          tab_type <- "WY|XZ"
          chkflag1 <- 1
        }
      }
      if( (length(tab_type)== 0) && (typeof(tab_type) == "character") ){
        tab_type <- "WY|XZ"
        chkflag1 <- 1
      }
      a0_val <- x$a0_vals[1]
      rand_control_diff_val <- x$rand_control_diff[1]
      if (chkflag1 == 1) {
        print("Since a0 and rand_control_diff vectors only have 1 element each and tab_type not equal to 'WY|XZ' or 'YW|XZ', tab_type was set to WY|XZ")
        print("This works towards putting all results on a single plot")
      }
    }
    if (length(x$subj_per_arm) == 1 & length(x$rand_control_diff) == 1) {
      chkflag1 <- 0
      if(!((length(tab_type)== 0) && (typeof(tab_type) == "character"))){
        if ((tab_type != "XY|WZ") & (tab_type != "YX|WZ")) {
          tab_type <- "XY|WZ"
          chkflag1 <- 1
        }
      }
      if( (length(tab_type)== 0) && (typeof(tab_type) == "character") ){
        tab_type <- "XY|WZ"
        chkflag1 <- 1
      }
      subj_per_arm_val <- x$subj_per_arm
      rand_control_diff_val <- x$rand_control_diff
      if (chkflag1 == 1) {
        print("Since sample size and rand_control_diff vectors only have 1 element each and tab_type not equal to 'XY|WZ' or 'YX|WZ', tab_type was set to XY|WZ")
        print("This works towards putting all results on a single plot")
      }
    }
    if (length(x$subj_per_arm) == 1 & length(x$effect_vals) == 1) {
      chkflag1 <- 0
      if(!((length(tab_type)== 0) && (typeof(tab_type) == "character"))){
        if ((tab_type != "XZ|WY") & (tab_type != "ZX|WY")) {
          tab_type <- "XZ|WY"
          chkflag1 <- 1
        }
      }
      if( (length(tab_type)== 0) && (typeof(tab_type) == "character") ){
        tab_type <- "XZ|WY"
        chkflag1 <- 1
      }
      subj_per_arm_val <- x$subj_per_arm
      effect_val <- x$effect_vals
      if (chkflag1 == 1) {
        print("Since sample size and effect_vals vectors only have 1 element each and tab_type not equal to 'XZ|WY' or 'ZX|WY', tab_type was set to XZ|WY")
        print("This works towards putting all results on a single plot")
      }
    }
    if (length(x$subj_per_arm) == 1 & length(x$a0_vals) == 1) {
      chkflag1 <- 0
      if(!((length(tab_type)== 0) && (typeof(tab_type) == "character"))){
        if ((tab_type != "YZ|WX") & (tab_type != "ZY|WX")) {
          tab_type <- "YZ|WX"
          chkflag1 <- 1
        }
      }
      if( (length(tab_type)== 0) && (typeof(tab_type) == "character") ){
        tab_type <- "YZ|WX"
        chkflag1 <- 1
      }
      subj_per_arm_val <- x$subj_per_arm
      a0_val <- x$a0_vals
      if (chkflag1 == 1) {
        print("Since sample size and a0 vectors only have 1 element each and tab_type not equal to 'YZ|WX' or 'ZY|WX', tab_type was set to YZ|WX")
        print("This works towards putting all results on a single plot")
      }
    }
    if (length(x$a0_vals) == 1 & length(x$effect_vals) == 1) {
      chkflag1 <- 0
      if(!((length(tab_type)== 0) && (typeof(tab_type) == "character"))){
        if ((tab_type != "WZ|XY") & (tab_type != "ZW|XY")) {
          tab_type <- "WZ|XY"
          chkflag1 <- 1
        }
      }
      if( (length(tab_type)== 0) && (typeof(tab_type) == "character") ){
        tab_type <- "WZ|XY"
        chkflag1 <- 1
      }
      a0_val <- x$a0_vals
      effect_val <- x$effect_vals
      if (chkflag1 == 1) {
        print("Since a0 and effect_vals vectors only have 1 element each and tab_type not equal to 'WZ|XY' or 'ZW|XY', tab_type was set to WZ|XY")
        print("This works towards putting all results on a single plot")
      }
    }
  }

  plot_error_checks(x, measure, tab_type, smooth, plot_out, subj_per_arm_val, a0_val, effect_val,
                          rand_control_diff_val, span)

  #If title is NULL, create a title from the two parameters that are fixed.
  if (is.null(title)){
    # (W): subj_per_arm_val,
    # (X): a0_val,
    # (Y): effect_val,
    # (Z): rand_control_diff_val
    #WX|YZ, WY|XZ, WZ|XY, XY|WZ, XZ|WY, YZ|WX, ZX|WY, XW|YZ, YW|XZ, YX|WZ, ZW|XY, ZX|WY, or ZY|WX

    if (tab_type == 'WX|YZ'){
      title = paste0('effect_val =', effect_val,' and rand_control_diff_val =', rand_control_diff_val)
    }
    if (tab_type == 'WY|XZ'){
      title = paste0('a0_val =', a0_val,' and rand_control_diff_val =', rand_control_diff_val)
    }
    if (tab_type == 'WZ|XY'){
      title = paste0('a0_val =', a0_val,' and effect_val =', effect_val)
    }
    if (tab_type == 'XY|WZ'){
      title = paste0('subj_per_arm_val =', subj_per_arm_val,' and rand_control_diff_val =', rand_control_diff_val)
    }
    if (tab_type == 'XZ|WY'){
      title = paste0('subj_per_arm_val =', subj_per_arm_val,' and effect_val =', effect_val)
    }
    if (tab_type == 'YZ|WX'){
      title = paste0('subj_per_arm_val =', subj_per_arm_val,' and a0_val =', a0_val)
    }
    if (tab_type == 'ZX|WY'){
      title = paste0('subj_per_arm_val =', subj_per_arm_val,' and effect_val =', effect_val)
    }
    if (tab_type == 'XW|YZ'){
      title = paste0('effect_val =', effect_val,' and rand_control_diff_val =', rand_control_diff_val)
    }
    if (tab_type == 'YW|XZ'){
      title = paste0('a0_val =', a0_val,' and rand_control_diff_val =', rand_control_diff_val)
    }
    if (tab_type == 'YX|WZ'){
      title = paste0('subj_per_arm_val =', subj_per_arm_val,' and rand_control_diff_val =', rand_control_diff_val)
    }
    if (tab_type == 'ZW|XY'){
      title = paste0('a0_val =', a0_val,' and effect_val =', effect_val)
    }
    if (tab_type == 'ZX|WY'){
      title = paste0('subj_per_arm_val =', subj_per_arm_val,' and effect_val =', effect_val)
    }
    if (tab_type == 'ZY|WX'){
      title = paste0('subj_per_arm_val =', subj_per_arm_val,' and a0_val =', a0_val)
    }
  }

  if (tab_type == "WX|YZ") {
    # Line Plot of Design Characteristics: Sample Size by a0
    # First generate the table using a call to print.
    table_mat <- print(x = x, measure = measure, tab_type = tab_type, effect_val = effect_val,
                             rand_control_diff_val = rand_control_diff_val, print_chg_warn = 0)
    table_mat <- data.frame(table_mat)
    if (dim(table_mat)[2]==1){
      tempcolname <- paste('X', as.character(x$a0_vals), sep='')
      colnames(table_mat) <- tempcolname
    }
    # Use a call to stats::loess to smooth simulation results is requested.
    if (smooth == TRUE) {
      colcnt <- dim(table_mat)[2]
      if (length(span) == 1)  #If span is single number, then use it to smooth all columns.
        span <- rep(span, colcnt)
      y <- as.numeric(rownames(table_mat))
      for (i in 1:colcnt) {
        loessfit <- stats::loess(table_mat[, i] ~ y, span = span[i], degree = degree, family = family)
        table_mat[, i] <- stats::predict(loessfit, data.frame(y = y))
      }
    }
    # If no plot needed, then return the table.  A user should only use this
    # option when they want to use the plot method to get a table of
    # smoothed results.
    if (plot_out == FALSE) {
      return(table_mat)
    }
    if (plot_out == TRUE) {
      #Shape the data for plotting.
      table_mat$row_names <- as.numeric(rownames(table_mat))
      table_dat <- reshape2::melt(table_mat, id = "row_names")
      colnames(table_dat) <- c("sample_size", "a0_val", measure)
      table_dat$sample_size <- as.numeric(table_dat$sample_size)

      #Construct proper values for the line variable.
      avec <- NULL
      if (table_dat$a0_val[1] == "table_mat") {
        table_dat$a0_val <- "0"
      } else if (table_dat$a0_val[1] != "table_mat") {
        for (i in 1:dim(table_dat)[1]) {
          avec[i] <- as.numeric(unlist(strsplit(as.character(table_dat[i, 2]), split = "X", fixed = TRUE)))[2]
        }
        table_dat$a0_val <- avec
        table_dat$a0_val <- as.character(table_dat$a0_val)
      }

      #Create the appropriate graph depending on the value of measure.
      if (measure == "power") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'sample_size', y = 'power', group = 'a0_val', color = 'a0_val')) + ggplot2::geom_line() +
          ggplot2::geom_point() + ggplot2::labs(color='a0 Value', x='Sample Size', y='Power') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "est") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'sample_size', y = 'est', group = 'a0_val', color = 'a0_val')) + ggplot2::geom_line() +
          ggplot2::geom_point() + ggplot2::labs(color='a0 Value', x='Sample Size', y='Estimate') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "var") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'sample_size', y = 'var', group = 'a0_val', color = 'a0_val')) + ggplot2::geom_line() +
          ggplot2::geom_point() + ggplot2::labs(color='a0 Value', x='Sample Size', y='Variance') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "bias") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'sample_size', y = 'bias', group = 'a0_val', color = 'a0_val')) + ggplot2::geom_line() +
          ggplot2::geom_point() + ggplot2::labs(color='a0 Value', x='Sample Size', y='Bias')  +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "mse") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'sample_size', y = 'mse', group = 'a0_val', color = 'a0_val')) + ggplot2::geom_line() +
          ggplot2::geom_point() + ggplot2::labs(color='a0 Value', x='Sample Size', y='MSE')  +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
    }
  }
  if (tab_type == "WY|XZ") {
    # Line Plot of Design Characteristics: Sample Size by Effect
    # First generate the table using a call to print.
    table_mat <- print(x = x, measure = measure, tab_type = "WY|XZ", a0_val = a0_val,
                             rand_control_diff_val = rand_control_diff_val, print_chg_warn = 0)
    table_mat <- data.frame(table_mat)
    if (dim(table_mat)[2]==1){
      tempcolname <- paste('X', as.character(x$effect_vals), sep='')
      colnames(table_mat) <- tempcolname
    }
    # Use a call to stats::loess to smooth simulation results is requested.
    if (smooth == TRUE) {
      colcnt <- dim(table_mat)[2]
      if (length(span) == 1) #If span is single number, then use it to smooth all columns.
        span <- rep(span, colcnt)
      y <- as.numeric(rownames(table_mat))
      for (i in 1:colcnt) {
        loessfit <- stats::loess(table_mat[, i] ~ y, span = span[i], degree = degree, family = family)
        table_mat[, i] <- stats::predict(loessfit, data.frame(y = y))
      }
    }
    # If no plot needed, then return the table.  A user should only use this
    # option when they want to use the plot method to get a table of
    # smoothed results.
    if (plot_out == FALSE) {
      return(table_mat)
    }
    if (plot_out == TRUE) {
      #Shape the data for plotting.
      table_mat$row_names <- as.numeric(rownames(table_mat))
      table_dat <- reshape2::melt(table_mat, id = "row_names")
      colnames(table_dat) <- c("sample_size", "effect_val", measure)
      table_dat$sample_size <- as.numeric(table_dat$sample_size)

      #Construct proper values for the line variable.
      avec <- NULL
      for (i in 1:dim(table_dat)[1]) {
        avec[i] <- as.numeric(unlist(strsplit(as.character(table_dat[i, 2]), split = "X", fixed = TRUE)))[2]
      }
      table_dat$effect_val <- avec
      table_dat$effect_val <- as.character(table_dat$effect_val)

      #Create the appropriate graph depending on the value of measure.
      if (measure == "power") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'sample_size', y = 'power', group = 'effect_val', color = 'effect_val')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Effect Value', x='Sample Size', y='Power') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "est") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'sample_size', y = 'est', group = 'effect_val', color = 'effect_val')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Effect Value', x='Sample Size', y='Estimate') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "var") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'sample_size', y = 'var', group = 'effect_val', color = 'effect_val')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Effect Value', x='Sample Size', y='Variance') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "bias") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'sample_size', y = 'bias', group = 'effect_val', color = 'effect_val')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Effect Value', x='Sample Size', y='Bias') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "mse") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'sample_size', y = 'mse', group = 'effect_val', color = 'effect_val')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Effect Value', x='Sample Size', y='MSE') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
    }
  }
  if (tab_type == "WZ|XY") {
    # Line Plot of Design Characteristics: Sample Size by Control Differences
    # First generate the table using a call to print.
    table_mat <- print(x = x, measure = measure, tab_type = "WZ|XY", a0_val = a0_val,
                             effect_val = effect_val, print_chg_warn = 0)
    table_mat <- data.frame(table_mat)
    if (dim(table_mat)[2]==1){
      tempcolname <- paste('X', as.character(x$rand_control_diff), sep='')
      colnames(table_mat) <- tempcolname
    }
    # Use a call to stats::loess to smooth simulation results is requested.
    if (smooth == TRUE) {
      colcnt <- dim(table_mat)[2]
      if (length(span) == 1) #If span is single number, then use it to smooth all columns.
        span <- rep(span, colcnt)
      y <- as.numeric(rownames(table_mat))
      for (i in 1:colcnt) {
        loessfit <- stats::loess(table_mat[, i] ~ y, span = span[i], degree = degree, family = family)
        table_mat[, i] <- stats::predict(loessfit, data.frame(y = y))
      }
    }
    # If no plot needed, then return the table.  A user should only use this
    # option when they want to use the plot method to get a table of
    # smoothed results.
    if (plot_out == FALSE) {
      return(table_mat)
    }
    if (plot_out == TRUE) {
      #Shape the data for plotting.
      table_mat$row_names <- as.numeric(rownames(table_mat))
      table_dat <- reshape2::melt(table_mat, id = "row_names")
      colnames(table_dat) <- c("sample_size", "control_differences", measure)
      table_dat$sample_size <- as.numeric(table_dat$sample_size)

      #Construct proper values for the line variable.
      avec <- NULL
      if (table_dat$control_differences[1] == "table_mat") {
        table_dat$control_differences <- "1"
      } else if (table_dat$control_differences[1] != "table_mat") {
        for (i in 1:dim(table_dat)[1]) {
          avec[i] <- as.numeric(unlist(strsplit(as.character(table_dat[i, 2]), split = "X", fixed = TRUE)))[2]
        }
        table_dat$control_differences <- avec
        table_dat$control_differences <- as.character(table_dat$control_differences)
      }

      #Create the appropriate graph depending on the value of measure.
      if (measure == "power") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'sample_size', y = 'power', group = 'control_differences', color = 'control_differences')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Control\nDifferences', x='Sample Size', y='Power') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "est") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'sample_size', y = 'est', group = 'control_differences', color = 'control_differences')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Control\nDifferences', x='Sample Size', y='Estimate') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "var") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'sample_size', y = 'var', group = 'control_differences', color = 'control_differences')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Control\nDifferences', x='Sample Size', y='Variance')+
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "bias") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'sample_size', y = 'bias', group = 'control_differences', color = 'control_differences')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Control\nDifferences', x='Sample Size', y='Bias') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "mse") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'sample_size', y = 'mse', group = 'control_differences', color = 'control_differences')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Control\nDifferences', x='Sample Size', y='MSE') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }

    }
  }

  if (tab_type == "XY|WZ") {
    # Line Plot of Design Characteristics: a0 by Effect
    # First generate the table using a call to print.
    table_mat <- print(x = x, measure = measure, tab_type = "XY|WZ", subj_per_arm_val = subj_per_arm_val,
                             rand_control_diff_val = rand_control_diff_val, print_chg_warn = 0)
    table_mat <- data.frame(table_mat)
    if (dim(table_mat)[2]==1){
      tempcolname <- paste('X', as.character(x$effect_vals), sep='')
      colnames(table_mat) <- tempcolname
    }
    # Use a call to stats::loess to smooth simulation results is requested.
    if (smooth == TRUE) {
      colcnt <- dim(table_mat)[2]
      if (length(span) == 1) #If span is single number, then use it to smooth all columns.
        span <- rep(span, colcnt)
      y <- as.numeric(rownames(table_mat))
      for (i in 1:colcnt) {
        loessfit <- stats::loess(table_mat[, i] ~ y, span = span[i], degree = degree, family = family)
        table_mat[, i] <- stats::predict(loessfit, data.frame(y = y))
      }
    }
    # If no plot needed, then return the table.  A user should only use this
    # option when they want to use the plot method to get a table of
    # smoothed results.
    if (plot_out == FALSE) {
      return(table_mat)
    }
    if (plot_out == TRUE) {
      #Shape the data for plotting.
      table_mat$row_names <- as.numeric(rownames(table_mat))
      table_dat <- reshape2::melt(table_mat, id = "row_names")
      colnames(table_dat) <- c("a0_val", "effect_val", measure)
      table_dat$a0_val <- as.numeric(table_dat$a0_val)

      #Construct proper values for the line variable.
      avec <- NULL
      for (i in 1:dim(table_dat)[1]) {
        avec[i] <- as.numeric(unlist(strsplit(as.character(table_dat[i, 2]), split = "X", fixed = TRUE)))[2]
      }
      table_dat$effect_val <- avec
      table_dat$effect_val <- as.character(table_dat$effect_val)

      #Create the appropriate graph depending on the value of measure.
      if (measure == "power") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'a0_val', y = 'power', group = 'effect_val', color = 'effect_val')) + ggplot2::geom_line() +
          ggplot2::geom_point() + ggplot2::labs(color='Effect Value', x='a0 Value', y='Power') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "est") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'a0_val', y = 'est', group = 'effect_val', color = 'effect_val')) + ggplot2::geom_line() +
          ggplot2::geom_point() + ggplot2::labs(color='Effect Value', x='a0 Value', y='Estimate') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "var") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'a0_val', y = 'var', group = 'effect_val', color = 'effect_val')) + ggplot2::geom_line() +
          ggplot2::geom_point() + ggplot2::labs(color='Effect Value', x='a0 Value', y='Variance') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "bias") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'a0_val', y = 'bias', group = 'effect_val', color = 'effect_val')) + ggplot2::geom_line() +
          ggplot2::geom_point() + ggplot2::labs(color='Effect Value', x='a0 Value', y='Bias') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "mse") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'a0_val', y = 'mse', group = 'effect_val', color = 'effect_val')) + ggplot2::geom_line() +
          ggplot2::geom_point() + ggplot2::labs(color='Effect Value', x='a0 Value', y='MSE') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
    }
  }

  if (tab_type == "XZ|WY") {
    # Line Plot of Design Characteristics: a0 by Control Diffferences
    # First generate the table using a call to print.
    table_mat <- print(x = x, measure = measure, tab_type = "XZ|WY", subj_per_arm_val = subj_per_arm_val,
                             effect_val = effect_val, print_chg_warn = 0)
    table_mat <- data.frame(table_mat)
    if (dim(table_mat)[2]==1){
      tempcolname <- paste('X', as.character(x$rand_control_diff), sep='')
      colnames(table_mat) <- tempcolname
    }
    # Use a call to stats::loess to smooth simulation results is requested.
    if (smooth == TRUE) {
      colcnt <- dim(table_mat)[2]
      if (length(span) == 1) #If span is single number, then use it to smooth all columns.
        span <- rep(span, colcnt)
      y <- as.numeric(rownames(table_mat))
      for (i in 1:colcnt) {
        loessfit <- stats::loess(table_mat[, i] ~ y, span = span[i], degree = degree, family = family)
        table_mat[, i] <- stats::predict(loessfit, data.frame(y = y))
      }
    }
    # If no plot needed, then return the table.  A user should only use this
    # option when they want to use the plot method to get a table of
    # smoothed results.
    if (plot_out == FALSE) {
      return(table_mat)
    }
    if (plot_out == TRUE) {
      #Shape the data for plotting.
      table_mat$row_names <- as.numeric(rownames(table_mat))
      table_dat <- reshape2::melt(table_mat, id = "row_names")
      colnames(table_dat) <- c("a0_val", "control_differences", measure)
      table_dat$a0_val <- as.numeric(table_dat$a0_val)

      #Construct proper values for the line variable.
      avec <- NULL
      for (i in 1:dim(table_dat)[1]) {
        avec[i] <- as.numeric(unlist(strsplit(as.character(table_dat[i, 2]), split = "X", fixed = TRUE)))[2]
      }
      table_dat$control_differences <- avec
      table_dat$control_differences <- as.character(table_dat$control_differences)

      #Create the appropriate graph depending on the value of measure.
      if (measure == "power") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'a0_val', y = 'power', group = 'control_differences', color = 'control_differences')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Control\nDifferences', x='a0 Value', y='Power') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "est") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'a0_val', y = 'est', group = 'control_differences', color = 'control_differences')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Control\nDifferences', x='a0 Value', y='Estimate') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "var") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'a0_val', y = 'var', group = 'control_differences', color = 'control_differences')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Control\nDifferences', x='a0 Value', y='Variance') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "bias") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'a0_val', y = 'bias', group = 'control_differences', color = 'control_differences')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Control\nDifferences', x='a0 Value', y='Bias') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "mse") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'a0_val', y = 'mse', group = 'control_differences', color = 'control_differences')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Control\nDifferences', x='a0 Value', y='MSE') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }

    }
  }

  if (tab_type == "YZ|WX") {
    # Line Plot of Design Characteristics: Effect by Control Diffferences
    # First generate the table using a call to print.
    table_mat <- print(x = x, measure = measure, tab_type = "YZ|WX", subj_per_arm_val = subj_per_arm_val,
                             a0_val = a0_val, print_chg_warn = 0)
    table_mat <- data.frame(table_mat)
    if (dim(table_mat)[2]==1){
      tempcolname <- paste('X', as.character(x$rand_control_diff), sep='')
      colnames(table_mat) <- tempcolname
    }
    # Use a call to stats::loess to smooth simulation results is requested.
    if (smooth == TRUE) {
      colcnt <- dim(table_mat)[2]
      if (length(span) == 1) #If span is single number, then use it to smooth all columns.
        span <- rep(span, colcnt)
      y <- as.numeric(rownames(table_mat))
      for (i in 1:colcnt) {
        loessfit <- stats::loess(table_mat[, i] ~ y, span = span[i], degree = degree, family = family)
        table_mat[, i] <- stats::predict(loessfit, data.frame(y = y))
      }
    }
    # If no plot needed, then return the table.  A user should only use this
    # option when they want to use the plot method to get a table of
    # smoothed results.
    if (plot_out == FALSE) {
      return(table_mat)
    }
    if (plot_out == TRUE) {
      #Shape the data for plotting.
      table_mat$row_names <- as.numeric(rownames(table_mat))
      table_dat <- reshape2::melt(table_mat, id = "row_names")
      colnames(table_dat) <- c("effect_val", "control_differences", measure)
      table_dat$effect_val <- as.numeric(table_dat$effect_val)

      #Construct proper values for the line variable.
      avec <- NULL
      if (table_dat$control_differences[1] == "table_mat") {
        table_dat$control_differences <- "1"
      } else if (table_dat$control_differences[1] != "table_mat") {
        for (i in 1:dim(table_dat)[1]) {
          avec[i] <- as.numeric(unlist(strsplit(as.character(table_dat[i, 2]), split = "X", fixed = TRUE)))[2]
        }
        table_dat$control_differences <- avec
        table_dat$control_differences <- as.character(table_dat$control_differences)
      }

      #Create the appropriate graph depending on the value of measure.
      if (measure == "power") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'effect_val', y = 'power', group = 'control_differences', color = 'control_differences')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Control\nDifferences', x='Effect Value', y='Power') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "est") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'effect_val', y = 'est', group = 'control_differences', color = 'control_differences')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Control\nDifferences', x='Effect Value', y='Estimate') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "var") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'effect_val', y = 'var', group = 'control_differences', color = 'control_differences')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Control\nDifferences', x='Effect Value', y='Variance') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "bias") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'effect_val', y = 'bias', group = 'control_differences', color = 'control_differences')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Control\nDifferences', x='Effect Value', y='Bias') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "mse") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'effect_val', y = 'mse', group = 'control_differences', color = 'control_differences')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Control\nDifferences', x='Effect Value', y='MSE') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
    }
  }

  if (tab_type == "ZX|WY") {
    # Line Plot of Design Characteristics: Control Diffferences by a0
    # First generate the table using a call to print.
    table_mat <- print(x = x, measure = measure, tab_type = "ZX|WY", subj_per_arm_val = subj_per_arm_val,
                             effect_val = effect_val, print_chg_warn = 0)
    #The following if block makes the 1D case work, if it occurs.
    #The return of a 1D vector from print(tab_type="ZX|WY")
    #  which is a 1byK row vector with class=matrix causes problems
    #when I make it into a dataframe.  The following code re-classifies
    #the object as a vector if necessary.
    if (is.matrix(table_mat) & dim(table_mat)[1]==1){
      colnamevec <- colnames(table_mat)
      table_mat <- as.vector(table_mat)
      names(table_mat) <- colnamevec
    }
    table_mat <- data.frame(table_mat)
    if (dim(table_mat)[2]==1){
      tempcolname <- paste('X', as.character(x$a0_vals), sep='')
      colnames(table_mat) <- tempcolname
    }
    # Use a call to stats::loess to smooth simulation results is requested.
    if (smooth == TRUE) {
      colcnt <- dim(table_mat)[2]
      if (length(span) == 1) #If span is single number, then use it to smooth all columns.
        span <- rep(span, colcnt)
      y <- as.numeric(rownames(table_mat))
      for (i in 1:colcnt) {
        loessfit <- stats::loess(table_mat[, i] ~ y, span = span[i], degree = degree, family = family)
        table_mat[, i] <- stats::predict(loessfit, data.frame(y = y))
      }
    }
    # If no plot needed, then return the table.  A user should only use this
    # option when they want to use the plot method to get a table of
    # smoothed results.
    if (plot_out == FALSE) {
      return(table_mat)
    }
    if (plot_out == TRUE) {
      #Shape the data for plotting.
      table_mat$row_names <- as.numeric(rownames(table_mat))
      table_dat <- reshape2::melt(table_mat, id = "row_names")
      colnames(table_dat) <- c("control_differences", "a0_val", measure)
      table_dat$control_differences <- as.numeric(table_dat$control_differences)

      #Construct proper values for the line variable.
      avec <- NULL
      if (table_dat$a0_val[1] == "table_mat") {
        table_dat$a0_val <- "1"
      } else if (table_dat$a0_val[1] != "table_mat") {
        for (i in 1:dim(table_dat)[1]) {
          avec[i] <- as.numeric(unlist(strsplit(as.character(table_dat[i, 2]), split = "X", fixed = TRUE)))[2]
        }
        table_dat$a0_val <- avec
        table_dat$a0_val <- as.character(table_dat$a0_val)
      }

      #Create the appropriate graph depending on the value of measure.
      if (measure == "power") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'control_differences', y = 'power', group = 'a0_val', color = 'a0_val')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='a0 Value', x='Control Differences', y='Power') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "est") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'control_differences', y = 'est', group = 'a0_val', color = 'a0_val')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='a0 Value', x='Control Differences', y='Estimate') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "var") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'control_differences', y = 'var', group = 'a0_val', color = 'a0_val')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='a0 Value', x='Control Differences', y='Variance') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "bias") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'control_differences', y = 'bias', group = 'a0_val', color = 'a0_val')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='a0 Value', x='Control Differences', y='Bias') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "mse") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'control_differences', y = 'mse', group = 'a0_val', color = 'a0_val')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='a0 Value', x='Control Differences', y='MSE') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
    }
  }

  if (tab_type == "XW|YZ") {
    # Line Plot of Design Characteristics: a0 by sample size
    # First generate the table using a call to print.
    table_mat <- print(x = x, measure = measure, tab_type = "XW|YZ", effect_val = effect_val,
                             rand_control_diff_val = rand_control_diff_val, print_chg_warn = 0)
    #The following if block makes the 1D case work, if it occurs.
    #The return of a 1D vector from print(tab_type="ZX|WY")
    #  which is a 1byK row vector with class=matrix causes problems
    #when I make it into a dataframe.  The following code re-classifies
    #the object as a vector if necessary.
    if (is.matrix(table_mat) & dim(table_mat)[1]==1){
      colnamevec <- colnames(table_mat)
      table_mat <- as.vector(table_mat)
      names(table_mat) <- colnamevec
    }
    table_mat <- data.frame(table_mat)
    if (dim(table_mat)[2]==1){
      tempcolname <- paste('X', as.character(x$subj_per_arm), sep='')
      colnames(table_mat) <- tempcolname
    }
    # Use a call to stats::loess to smooth simulation results is requested.
    if (smooth == TRUE) {
      colcnt <- dim(table_mat)[2]
      if (length(span) == 1) #If span is single number, then use it to smooth all columns.
        span <- rep(span, colcnt)
      y <- as.numeric(rownames(table_mat))
      for (i in 1:colcnt) {
        loessfit <- stats::loess(table_mat[, i] ~ y, span = span[i], degree = degree, family = family)
        table_mat[, i] <- stats::predict(loessfit, data.frame(y = y))
      }
    }
    # If no plot needed, then return the table.  A user should only use this
    # option when they want to use the plot method to get a table of
    # smoothed results.
    if (plot_out == FALSE) {
      return(table_mat)
    }
    if (plot_out == TRUE) {
      #Shape the data for plotting.
      table_mat$row_names <- as.numeric(rownames(table_mat))
      table_dat <- reshape2::melt(table_mat, id = "row_names")
      colnames(table_dat) <- c("a0_val", "subj_per_arm", measure)
      table_dat$a0_val <- as.numeric(table_dat$a0_val)

      #Construct proper values for the line variable.
      avec <- NULL
      if (table_dat$subj_per_arm[1] == "table_mat") {
        table_dat$subj_per_arm <- "1"
      } else if (table_dat$subj_per_arm[1] != "table_mat") {
        for (i in 1:dim(table_dat)[1]) {
          avec[i] <- as.numeric(unlist(strsplit(as.character(table_dat[i, 2]), split = "X", fixed = TRUE)))[2]
        }
        table_dat$subj_per_arm <- avec
        table_dat$subj_per_arm <- as.character(table_dat$subj_per_arm)
      }

      #Create the appropriate graph depending on the value of measure.
      if (measure == "power") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'a0_val', y = 'power', group = 'subj_per_arm', color = 'subj_per_arm')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Subjects per Arm', x='a0_val', y='Power') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "est") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'a0_val', y = 'est', group = 'subj_per_arm', color = 'subj_per_arm')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Subjects per Arm', x='a0_val', y='Estimate') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "var") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'a0_val', y = 'var', group = 'subj_per_arm', color = 'subj_per_arm')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Subjects per Arm', x='a0_val', y='Variance') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "bias") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'a0_val', y = 'bias', group = 'subj_per_arm', color = 'subj_per_arm')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Subjects per Arm', x='a0_val', y='Bias') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "mse") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'a0_val', y = 'mse', group = 'subj_per_arm', color = 'subj_per_arm')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Subjects per Arm', x='a0_val', y='MSE') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
    }
  }

  if (tab_type == "YW|XZ") {
    # Line Plot of Design Characteristics: Effect by sample size
    # First generate the table using a call to print.
    table_mat <- print(x = x, measure = measure, tab_type = "YW|XZ",
                             a0_val = a0_val, rand_control_diff_val=rand_control_diff_val, print_chg_warn = 0)
    #The following if block makes the 1D case work, if it occurs.
    #The return of a 1D vector from print(tab_type="ZX|WY")
    #  which is a 1byK row vector with class=matrix causes problems
    #when I make it into a dataframe.  The following code re-classifies
    #the object as a vector if necessary.
    if (is.matrix(table_mat) & dim(table_mat)[1]==1){
      colnamevec <- colnames(table_mat)
      table_mat <- as.vector(table_mat)
      names(table_mat) <- colnamevec
    }
    table_mat <- data.frame(table_mat)
    if (dim(table_mat)[2]==1){
      tempcolname <- paste('X', as.character(x$subj_per_arm), sep='')
      colnames(table_mat) <- tempcolname
    }
    # Use a call to stats::loess to smooth simulation results is requested.
    if (smooth == TRUE) {
      colcnt <- dim(table_mat)[2]
      if (length(span) == 1) #If span is single number, then use it to smooth all columns.
        span <- rep(span, colcnt)
      y <- as.numeric(rownames(table_mat))
      for (i in 1:colcnt) {
        loessfit <- stats::loess(table_mat[, i] ~ y, span = span[i], degree = degree, family = family)
        table_mat[, i] <- stats::predict(loessfit, data.frame(y = y))
      }
    }
    # If no plot needed, then return the table.  A user should only use this
    # option when they want to use the plot method to get a table of
    # smoothed results.
    if (plot_out == FALSE) {
      return(table_mat)
    }
    if (plot_out == TRUE) {
      #Shape the data for plotting.
      table_mat$row_names <- as.numeric(rownames(table_mat))
      table_dat <- reshape2::melt(table_mat, id = "row_names")
      colnames(table_dat) <- c("effect_val", "subj_per_arm", measure)
      table_dat$effect_val <- as.numeric(table_dat$effect_val)

      #Construct proper values for the line variable.
      avec <- NULL
      if (table_dat$subj_per_arm[1] == "table_mat") {
        table_dat$subj_per_arm <- "1"
      } else if (table_dat$subj_per_arm[1] != "table_mat") {
        for (i in 1:dim(table_dat)[1]) {
          avec[i] <- as.numeric(unlist(strsplit(as.character(table_dat[i, 2]), split = "X", fixed = TRUE)))[2]
        }
        table_dat$subj_per_arm <- avec
        table_dat$subj_per_arm <- as.character(table_dat$subj_per_arm)
      }

      #Create the appropriate graph depending on the value of measure.
      if (measure == "power") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'effect_val', y = 'power', group = 'subj_per_arm', color = 'subj_per_arm')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Subjects per Arm', x='Effect Value', y='Power') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "est") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'effect_val', y = 'est', group = 'subj_per_arm', color = 'subj_per_arm')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Subjects per Arm', x='Effect Value', y='Estimate') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "var") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'effect_val', y = 'var', group = 'subj_per_arm', color = 'subj_per_arm')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Subjects per Arm', x='Effect Value', y='Variance') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "bias") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'effect_val', y = 'bias', group = 'subj_per_arm', color = 'subj_per_arm')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Subjects per Arm', x='Effect Value', y='Bias') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "mse") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'effect_val', y = 'mse', group = 'control_differences', color = 'control_differences')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Control\nDifferences', x='Effect Value', y='MSE') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
    }
  }

  if (tab_type == "ZW|XY") {
    # Line Plot of Design Characteristics: Control differences by sample size
    # First generate the table using a call to print.
    table_mat <- print(x = x, measure = measure, tab_type = "ZW|XY",
                             a0_val = a0_val, effect_val=effect_val, print_chg_warn = 0)
    #The following if block makes the 1D case work, if it occurs.
    #The return of a 1D vector from print(tab_type="ZX|WY")
    #  which is a 1byK row vector with class=matrix causes problems
    #when I make it into a dataframe.  The following code re-classifies
    #the object as a vector if necessary.
    if (is.matrix(table_mat) & dim(table_mat)[1]==1){
      colnamevec <- colnames(table_mat)
      table_mat <- as.vector(table_mat)
      names(table_mat) <- colnamevec
    }
    table_mat <- data.frame(table_mat)
    if (dim(table_mat)[2]==1){
      tempcolname <- paste('X', as.character(x$subj_per_arm), sep='')
      colnames(table_mat) <- tempcolname
    }
    # Use a call to stats::loess to smooth simulation results is requested.
    if (smooth == TRUE) {
      colcnt <- dim(table_mat)[2]
      if (length(span) == 1) #If span is single number, then use it to smooth all columns.
        span <- rep(span, colcnt)
      y <- as.numeric(rownames(table_mat))
      for (i in 1:colcnt) {
        loessfit <- stats::loess(table_mat[, i] ~ y, span = span[i], degree = degree, family = family)
        table_mat[, i] <- stats::predict(loessfit, data.frame(y = y))
      }
    }
    # If no plot needed, then return the table.  A user should only use this
    # option when they want to use the plot method to get a table of
    # smoothed results.
    if (plot_out == FALSE) {
      return(table_mat)
    }
    if (plot_out == TRUE) {
      #Shape the data for plotting.
      table_mat$row_names <- as.numeric(rownames(table_mat))
      table_dat <- reshape2::melt(table_mat, id = "row_names")
      colnames(table_dat) <- c("control_differences", "subj_per_arm", measure)
      table_dat$control_differences <- as.numeric(table_dat$control_differences)

      #Construct proper values for the line variable.
      avec <- NULL
      if (table_dat$subj_per_arm[1] == "table_mat") {
        table_dat$subj_per_arm <- "1"
      } else if (table_dat$subj_per_arm[1] != "table_mat") {
        for (i in 1:dim(table_dat)[1]) {
          avec[i] <- as.numeric(unlist(strsplit(as.character(table_dat[i, 2]), split = "X", fixed = TRUE)))[2]
        }
        table_dat$subj_per_arm <- avec
        table_dat$subj_per_arm <- as.character(table_dat$subj_per_arm)
      }

      #Create the appropriate graph depending on the value of measure.
      if (measure == "power") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'control_differences', y = 'power', group = 'subj_per_arm', color = 'subj_per_arm')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Subjects per Arm', x='Control Differences', y='Power') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "est") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'control_differences', y = 'est', group = 'subj_per_arm', color = 'subj_per_arm')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Subjects per Arm', x='Control Differences', y='Estimate') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "var") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'control_differences', y = 'var', group = 'subj_per_arm', color = 'subj_per_arm')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Subjects per Arm', x='Control Differences', y='Variance') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "bias") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'control_differences', y = 'bias', group = 'subj_per_arm', color = 'subj_per_arm')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Subjects per Arm', x='Control Differences', y='Bias') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "mse") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'control_differences', y = 'mse', group = 'control_differences', color = 'control_differences')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='Control\nDifferences', x='Control Differences', y='MSE') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
    }
  }

  if (tab_type == "ZY|WX") {
    # Line Plot of Design Characteristics: a0 by Effect
    # First generate the table using a call to print.
    table_mat <- print(x = x, measure = measure, tab_type = "ZY|WX", subj_per_arm_val = subj_per_arm_val,
                             a0_val = a0_val, print_chg_warn = 0)
    #The following if block makes the 1D case work, if it occurs.
    #The return of a 1D vector from print(tab_type="ZX|WY")
    #  which is a 1byK row vector with class=matrix causes problems
    #when I make it into a dataframe.  The following code re-classifies
    #the object as a vector if necessary.
    if (is.matrix(table_mat) & dim(table_mat)[1]==1){
      colnamevec <- colnames(table_mat)
      table_mat <- as.vector(table_mat)
      names(table_mat) <- colnamevec
    }
    table_mat <- data.frame(table_mat)
    if (dim(table_mat)[2]==1){
      tempcolname <- paste('X', as.character(x$effect_vals), sep='')
      colnames(table_mat) <- tempcolname
    }
    # Use a call to stats::loess to smooth simulation results is requested.
    if (smooth == TRUE) {
      colcnt <- dim(table_mat)[2]
      if (length(span) == 1) #If span is single number, then use it to smooth all columns.
        span <- rep(span, colcnt)
      y <- as.numeric(rownames(table_mat))
      for (i in 1:colcnt) {
        loessfit <- stats::loess(table_mat[, i] ~ y, span = span[i], degree = degree, family = family)
        table_mat[, i] <- stats::predict(loessfit, data.frame(y = y))
      }
    }
    # If no plot needed, then return the table.  A user should only use this
    # option when they want to use the plot method to get a table of
    # smoothed results.
    if (plot_out == FALSE) {
      return(table_mat)
    }
    if (plot_out == TRUE) {
      #Shape the data for plotting.
      table_mat$row_names <- as.numeric(rownames(table_mat))
      table_dat <- reshape2::melt(table_mat, id = "row_names")
      colnames(table_dat) <- c("control_differences", "effect_val", measure)
      table_dat$control_differences <- as.numeric(table_dat$control_differences)

      #Construct proper values for the line variable.
      avec <- NULL
      for (i in 1:dim(table_dat)[1]) {
        avec[i] <- as.numeric(unlist(strsplit(as.character(table_dat[i, 2]), split = "X", fixed = TRUE)))[2]
      }
      table_dat$effect_val <- avec
      table_dat$effect_val <- as.character(table_dat$effect_val)

      #Create the appropriate graph depending on the value of measure.
      if (measure == "power") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'control_differences', y = 'power', group = 'effect_val', color = 'effect_val')) + ggplot2::geom_line() +
          ggplot2::geom_point() + ggplot2::labs(color='Effect Value', x='Control Differences', y='Power') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "est") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'control_differences', y = 'est', group = 'effect_val', color = 'effect_val')) + ggplot2::geom_line() +
          ggplot2::geom_point() + ggplot2::labs(color='Effect Value', x='Control Differences', y='Estimate') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "var") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'control_differences', y = 'var', group = 'effect_val', color = 'effect_val')) + ggplot2::geom_line() +
          ggplot2::geom_point() + ggplot2::labs(color='Effect Value', x='Control Differences', y='Variance') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "bias") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'control_differences', y = 'bias', group = 'effect_val', color = 'effect_val')) + ggplot2::geom_line() +
          ggplot2::geom_point() + ggplot2::labs(color='Effect Value', x='Control Differences', y='Bias') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "mse") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'control_differences', y = 'mse', group = 'effect_val', color = 'effect_val')) + ggplot2::geom_line() +
          ggplot2::geom_point() + ggplot2::labs(color='Effect Value', x='Control Differences', y='MSE') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
    }
  }

  if (tab_type == "YX|WZ") {
    # Line Plot of Design Characteristics: Effect by a0
    # First generate the table using a call to print.
    table_mat <- print(x = x, measure = measure, tab_type = "YX|WZ", subj_per_arm_val = subj_per_arm_val,
                             rand_control_diff_val = rand_control_diff_val, print_chg_warn = 0)
    #The following if block makes the 1D case work, if it occurs.
    #The return of a 1D vector from print(tab_type="ZX|WY")
    #  which is a 1byK row vector with class=matrix causes problems
    #when I make it into a dataframe.  The following code re-classifies
    #the object as a vector if necessary.
    if (is.matrix(table_mat) & dim(table_mat)[1]==1){
      colnamevec <- colnames(table_mat)
      table_mat <- as.vector(table_mat)
      names(table_mat) <- colnamevec
    }
    table_mat <- data.frame(table_mat)
    if (dim(table_mat)[2]==1){
      tempcolname <- paste('X', as.character(x$a0_vals), sep='')
      colnames(table_mat) <- tempcolname
    }
    # Use a call to stats::loess to smooth simulation results is requested.
    if (smooth == TRUE) {
      colcnt <- dim(table_mat)[2]
      if (length(span) == 1) #If span is single number, then use it to smooth all columns.
        span <- rep(span, colcnt)
      y <- as.numeric(rownames(table_mat))
      for (i in 1:colcnt) {
        loessfit <- stats::loess(table_mat[, i] ~ y, span = span[i], degree = degree, family = family)
        table_mat[, i] <- stats::predict(loessfit, data.frame(y = y))
      }
    }
    # If no plot needed, then return the table.  A user should only use this
    # option when they want to use the plot method to get a table of
    # smoothed results.
    if (plot_out == FALSE) {
      return(table_mat)
    }
    if (plot_out == TRUE) {
      #Shape the data for plotting.
      table_mat$row_names <- as.numeric(rownames(table_mat))
      table_dat <- reshape2::melt(table_mat, id = "row_names")
      colnames(table_dat) <- c("effect_val", "a0_val", measure)
      table_dat$effect_val <- as.numeric(table_dat$effect_val)

      #Construct proper values for the line variable.
      avec <- NULL
      if (table_dat$a0_val[1] == "table_mat") {
        table_dat$a0_val <- "1"
      } else if (table_dat$a0_val[1] != "table_mat") {
        for (i in 1:dim(table_dat)[1]) {
          avec[i] <- as.numeric(unlist(strsplit(as.character(table_dat[i, 2]), split = "X", fixed = TRUE)))[2]
        }
        table_dat$a0_val <- avec
        table_dat$a0_val <- as.character(table_dat$a0_val)
      }

      #Create the appropriate graph depending on the value of measure.
      if (measure == "power") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'effect_val', y = 'power', group = 'a0_val', color = 'a0_val')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='a0 Value', x='Effect Value', y='Power') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "est") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'effect_val', y = 'est', group = 'a0_val', color = 'a0_val')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='a0 Value', x='Effect Value', y='Estimate') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "var") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'effect_val', y = 'var', group = 'a0_val', color = 'a0_val')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='a0 Value', x='Effect Value', y='Variance') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "bias") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'effect_val', y = 'bias', group = 'a0_val', color = 'a0_val')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='a0 Value', x='Effect Value', y='Bias') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
      if (measure == "mse") {
        p <- ggplot2::ggplot(data = table_dat, ggplot2::aes_string(x = 'effect_val', y = 'mse', group = 'a0_val', color = 'a0_val')) +
          ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::labs(color='a0 Value', x='Effect Value', y='MSE') +
          ggplot2::ggtitle(title)
        if (!is.null(ylim)){
          p <- p + ggplot2::ylim(ylim)
        }
        return(p)
      }
    }
  }
}
begglest/BayesCTDesign documentation built on Nov. 29, 2021, 10:34 p.m.