Generalized Pivotal Quantity for Confidence Interval for the Mean of a Normal Distribution Based on Censored Data

Generate a generalized pivotal quantity (GPQ) for a confidence interval for the mean of a Normal distribution based on singly or multiply censored data.

gpqCiNormSinglyCensored(n, n.cen, probs, nmc, method = "mle", 
    censoring.side = "left", seed = NULL, names = TRUE)

  gpqCiNormMultiplyCensored(n, cen.index, probs, nmc, method = "mle", 
    censoring.side = "left", seed = NULL, names = TRUE)

Arguments

n: positive integer $\ge 3$ indicating the sample size.
n.cen: for the case of singly censored data, a positive integer indicating the number of censored observations. The value of n.cen must be between 1 and n-2, inclusive.
cen.index: for the case of multiply censored data, a sorted vector of unique integers indicating the indices of the censored observations when the observations are “ordered”. The length of cen.index must be between 1 and n-2, inclusive, and the values of cen.index must be between 1 and n.
probs: numeric vector of values between 0 and 1 indicating the confidence level(s) associated with the GPQ(s).
nmc: positive integer $\ge 10$ indicating the number of Monte Carlo trials to run in order to compute the GPQ(s).
method: character string indicating the method to use for parameter estimation.

For singly censored data, possible values are "mle" (the default), "bcmle", "qq.reg", "qq.reg.w.cen.level", "impute.w.qq.reg",
"impute.w.qq.reg.w.cen.level", "impute.w.mle",
"iterative.impute.w.qq.reg", "m.est", and "half.cen.level". See the help file for enormCensored for details.

For multiply censored data, possible values are "mle" (the default), "qq.reg", "impute.w.qq.reg", and "half.cen.level". See the help file for enormCensored for details.
censoring.side: character string indicating on which side the censoring occurs. The possible values are "left" (the default) and "right".
seed: positive integer to pass to the function set.seed. This argument is ignored if seed=NULL (the default). Using the seed argument lets you reproduce the exact same result if all other arguments stay the same.
names: a logical scalar passed to quantile indicating whether to add a names attribute to the resulting GPQ(s). The default value is names=TRUE.

Details

The functions gpqCiNormSinglyCensored and gpqCiNormMultiplyCensored are called by
enormCensored when ci.method="gpq". They are used to construct generalized pivotal quantities to create confidence intervals for the mean $\mu$ of an assumed normal distribution.

This idea was introduced by Schmee et al. (1985) in the context of Type II singly censored data. The function gpqCiNormSinglyCensored generates GPQs using a modification of Algorithm 12.1 of Krishnamoorthy and Mathew (2009, p. 329). Algorithm 12.1 is used to generate GPQs for a tolerance interval. The modified algorithm for generating GPQs for confidence intervals for the mean $\mu$ is as follows:

Generate a random sample of $n$ observations from a standard normal (i.e., N(0,1)) distribution and let $z_{(1)}, z_{(2)}, \ldots, z_{(n)}$ denote the ordered (sorted) observations.
Set the smallest n.cen observations as censored.
Compute the estimates of $\mu$ and $\sigma$ by calling enormCensored using the method specified by the method argument, and denote these estimates as $\hat{\mu}^*, \; \hat{\sigma}^*$.
Compute the t-like pivotal quantity $\hat{t} = \hat{\mu}^*/\hat{\sigma}^*$.
Repeat steps 1-4 nmc times to produce an empirical distribution of the t-like pivotal quantity.

A two-sided $(1-\alpha)100\%$ confidence interval for $\mu$ is then computed as: $$[\hat{\mu} - \hat{t}_{1-(\alpha/2)} \hat{\sigma}, \; \hat{\mu} - \hat{t}_{\alpha/2} \hat{\sigma}]$$ where $\hat{t}_p$ denotes the $p$'th empirical quantile of the nmc generated $\hat{t}$ values.

Schmee at al. (1985) derived this method in the context of Type II singly censored data (for which these limits are exact within Monte Carlo error), but state that according to Regal (1982) this method produces confidence intervals that are close apporximations to the correct limits for Type I censored data.

The function gpqCiNormMultiplyCensored is an extension of this idea to multiply censored data. The algorithm is the same as for singly censored data, except Step 2 changes to:

2. Set observations as censored for elements of the argument cen.index that have the value TRUE.

The functions gpqCiNormSinglyCensored and gpqCiNormMultiplyCensored are computationally intensive and provided to the user to allow you to create your own tables.

Value

a numeric vector containing the GPQ(s).

References

Krishnamoorthy K., and T. Mathew. (2009). Statistical Tolerance Regions: Theory, Applications, and Computation. John Wiley and Sons, Hoboken.

Regal, R. (1982). Applying Order Statistic Censored Normal Confidence Intervals to Time Censored Data. Unpublished manuscript, University of Minnesota, Duluth, Department of Mathematical Sciences.

Schmee, J., D.Gladstein, and W. Nelson. (1985). Confidence Limits for Parameters of a Normal Distribution from Singly Censored Samples, Using Maximum Likelihood. Technometrics 27(2) 119–128.

Author

Steven P. Millard (EnvStats@ProbStatInfo.com)

Examples

  # Reproduce the entries for n=10 observations with n.cen=6 in Table 4 
  # of Schmee et al. (1985, p.122).
  #
  # Notes: 
  # 1. This table applies to right-censored data, and the 
  #    quantity "r" in this table refers to the number of 
  #    uncensored observations.
  #
  # 2. Passing a value for the argument "seed" simply allows 
  #    you to reproduce this example.  

  # NOTE:  Here to save computing time for the sake of example, we will specify 
  #        just 100 Monte Carlos, whereas Krishnamoorthy and Mathew (2009) 
  #        suggest *10,000* Monte Carlos.

  # Here are the values given in Schmee et al. (1985):
  Schmee.values <- c(-3.59, -2.60, -1.73, -0.24, 0.43, 0.58, 0.73)
  probs <- c(0.025, 0.05, 0.1, 0.5, 0.9, 0.95, 0.975)
  names(Schmee.values) <- paste(probs * 100, "%", sep = "")

  Schmee.values
#>  2.5%    5%   10%   50%   90%   95% 97.5% 
#> -3.59 -2.60 -1.73 -0.24  0.43  0.58  0.73 
  # 2.5%    5%   10%   50%   90%   95% 97.5% 
  #-3.59 -2.60 -1.73 -0.24  0.43  0.58  0.73

  gpqs <- gpqCiNormSinglyCensored(n = 10, n.cen = 6, probs = probs, 
    nmc = 100, censoring.side = "right", seed = 529)

  round(gpqs, 2)
#>  2.5%    5%   10%   50%   90%   95% 97.5% 
#> -2.46 -2.03 -1.38 -0.14  0.54  0.65  0.84 
  # 2.5%    5%   10%   50%   90%   95% 97.5% 
  #-2.46 -2.03 -1.38 -0.14  0.54  0.65  0.84 

  # This is what you get if you specify nmc = 1000 with the 
  # same value for seed:
  #-----------------------------------------------
  # 2.5%    5%   10%   50%   90%   95% 97.5% 
  #-3.50 -2.49 -1.67 -0.25  0.41  0.57  0.71


  # Clean up
  #---------
  rm(Schmee.values, probs, gpqs)

  #==========

  # Example of using gpqCiNormMultiplyCensored
  #-------------------------------------------

  # Consider the following set of multiply left-censored data:
  dat <- 12:16
  censored <- c(TRUE, FALSE, TRUE, FALSE, FALSE)

  # Since the data are "ordered" we can identify the indices of the 
  # censored observations in the ordered data as follow:

  cen.index <- (1:length(dat))[censored]
  cen.index
#> [1] 1 3
  #[1] 1 3

  # Now we can generate a GPQ using gpqCiNormMultiplyCensored.
  # Here we'll generate a GPQs to use to create a  
  # 95% confidence interval for left-censored data.

  # NOTE:  Here to save computing time for the sake of example, we will specify 
  #        just 100 Monte Carlos, whereas Krishnamoorthy and Mathew (2009) 
  #        suggest *10,000* Monte Carlos.

  gpqCiNormMultiplyCensored(n = 5, cen.index = cen.index,  
    probs = c(0.025, 0.975), nmc = 100, seed = 237)
#> Warning: NA/NaN function evaluation
#> Warning: NA/NaN function evaluation
#> Warning: NA/NaN function evaluation
#>      2.5%     97.5% 
#> -1.315592  1.848513 
  #     2.5%     97.5% 
  #-1.315592  1.848513 

  #----------

  # Clean up
  #---------
  rm(dat, censored, cen.index)