Perform a goodness-of-fit test to determine whether a data set appears to come from a normal distribution, lognormal distribution, or lognormal distribution (alternative parameterization) based on a sample of data that has been subjected to Type I or Type II censoring.

gofTestCensored(x, censored, censoring.side = "left", test = "sf", 
    distribution = "norm", est.arg.list = NULL, 
    prob.method = "hirsch-stedinger", plot.pos.con = 0.375, 
    keep.data = TRUE, data.name = NULL, censoring.name = NULL)

Arguments

x

numeric vector of observations. Missing (NA), undefined (NaN), and infinite (Inf, -Inf) values are allowed but will be removed.

censored

numeric or logical vector indicating which values of x are censored. This must be the same length as x. If the mode of censored is "logical", TRUE values correspond to elements of x that are censored, and FALSE values correspond to elements of x that are not censored. If the mode of censored is "numeric", it must contain only 1's and 0's; 1 corresponds to TRUE and 0 corresponds to FALSE. Missing (NA) values are allowed but will be removed.

censoring.side

character string indicating on which side the censoring occurs. The possible values are "left" (the default) and "right".

test

character string defining which goodness-of-fit test to perform. Possible values are: "sw" Shapiro-Wilk. "sf" Shapiro-Francia; the default. "ppcc" Probability Plot Correlation Coefficient.

The Shapiro-Wilk test is only available for singly censored data.

See the DETAILS section for more information.

distribution

a character string denoting the abbreviation of the assumed distribution. Only continous distributions are allowed. See the help file for Distribution.df for a list of distributions and their abbreviations. Examples of possible values are:
distribution="norm" (Normal distribution; the default),
distribution="lnorm" (Lognormal distribution),
distribution="lnormAlt" (Lognormal distribution, alternative parameterization,
distribution="gamma" (Gamma distribution),
distribution="gammaAlt" (Gamma distribution, alternative parameterization).

The results for the goodness-of-fit test are identical for distribution="lnorm" and distribution="lnormAlt", the only difference in ouput is that when distribution="lnorm" the returned estimated parameters are the mean and standard deviation based on the log-scale of the data, whereas when distribution="lnormAlt" the returned estimated parameters are the mean and coefficient of variation based on the original scale of the data.

Also, the results for the goodness-of-fit test are identical for distribution="gamma" and distribution="gammaAlt", the only difference in ouput is that when distribution="gamma" the returned estimated parameters are the shape and scale, whereas when distribution="lnormAlt" the returned estimated parameters are the mean and coefficient of variation.

est.arg.list

a list of arguments to be passed to the function estimating the distribution parameters. For example, if distribution="lnormAlt" setting
est.arg.list=list(method="bcmle") indicates using the bias-corrected maximum likelihood estimators (see the help file for elnormAltCensored). The default value is est.arg.list=NULL so that all default values for the estimating function are used. The estimated parameters are provided in the output merely for information, and the choice of the method of estimation has no effect on the goodness-of-fit test statistic or p-value.

prob.method

character string indicating what method to use to compute the plotting positions (empirical probabilities) when test="sf" or test="ppcc". Possible values are:
"modified kaplan-meier" (modification of product-limit method of Kaplan and Meier (1958)),
"nelson" (hazard plotting method of Nelson (1972)),
"michael-schucany" (generalization of the product-limit method due to Michael and Schucany (1986)), and
"hirsch-stedinger" (generalization of the product-limit method due to Hirsch and Stedinger (1987)).
The default value is prob.method="hirsch-stedinger".

The "nelson" method is only available for censoring.side="right", and the "modified kaplan-meier" method is only available for censoring.side="left". See the DETAILS section and the help file for ppointsCensored for more information.

plot.pos.con

numeric scalar between 0 and 1 containing the value of the plotting position constant to use when test="sf" or test="ppcc". The default value is
plot.pos.con=0.375. See the DETAILS section and the help file for
ppointsCensored for more information.

keep.data

logical scalar indicating whether to return the original data. The default value is keep.data=TRUE.

data.name

optional character string indicating the name for the data used for argument x.

censoring.name

optional character string indicating the name for the data used for argument censored.

Details

Let \(\underline{x} = c(x_1, x_2, \ldots, x_N)\) denote a vector of \(N\) observations from from some distribution with cdf \(F\). Suppose we want to test the null hypothesis that \(F\) is the cdf of a normal (Gaussian) distribution with some arbitrary mean \(\mu\) and standard deviation \(\sigma\) against the alternative hypothesis that \(F\) is the cdf of some other distribution. The table below shows the random variable for which \(F\) is the assumed cdf, given the value of the argument distribution.

Value ofRandom Variable for
distributionDistribution Namewhich \(F\) is the cdf
"norm"Normal\(X\)
"lnorm"Lognormal (Log-space)\(log(X)\)
"lnormAlt"Lognormal (Untransformed)\(log(X)\)

Assume \(n\) (\(0 < n < N\)) of these observations are known and \(c\) (\(c=N-n\)) of these observations are all censored below (left-censored) or all censored above (right-censored) at \(k\) fixed censoring levels $$T_1, T_2, \ldots, T_k; \; k \ge 1 \;\;\;\;\;\; (1)$$ For the case when \(k \ge 2\), the data are said to be Type I multiply censored. For the case when \(k=1\), set \(T = T_1\). If the data are left-censored and all \(n\) known observations are greater than or equal to \(T\), or if the data are right-censored and all \(n\) known observations are less than or equal to \(T\), then the data are said to be Type I singly censored (Nelson, 1982, p.7), otherwise they are considered to be Type I multiply censored.

Let \(c_j\) denote the number of observations censored below or above censoring level \(T_j\) for \(j = 1, 2, \ldots, k\), so that $$\sum_{i=1}^k c_j = c \;\;\;\;\;\; (2)$$ Let \(x_{(1)}, x_{(2)}, \ldots, x_{(N)}\) denote the “ordered” observations, where now “observation” means either the actual observation (for uncensored observations) or the censoring level (for censored observations). For right-censored data, if a censored observation has the same value as an uncensored one, the uncensored observation should be placed first. For left-censored data, if a censored observation has the same value as an uncensored one, the censored observation should be placed first.

Note that in this case the quantity \(x_{(i)}\) does not necessarily represent the \(i\)'th “largest” observation from the (unknown) complete sample.

Note that for singly left-censored data: $$x_{(1)} = x_{(2)} = \cdots = x_{(c)} = T \;\;\;\;\;\; (3)$$ and for singly right-censored data: $$x_{(n+1)} = x_{(n+2)} = \cdots = x_{(N)} = T \;\;\;\;\;\; (4)$$

Finally, let \(\Omega\) (omega) denote the set of \(n\) subscripts in the “ordered” sample that correspond to uncensored observations.

Shapiro-Wilk Goodness-of-Fit Test for Singly Censored Data (test="sw")
Equation (8) in the help file for gofTest shows that for the case of complete ordered data \(\underline{x}\), the Shapiro-Wilk \(W\)-statistic is the same as the square of the sample product-moment correlation between the vectors \(\underline{a}\) and \(\underline{x}\): $$W = r(\underline{a}, \underline{x})^2 \;\;\;\;\;\; (5)$$ where $$r(\underline{x}, \underline{y}) = \frac{\sum_{i=1}^N (x_i - \bar{x})(y_i - \bar{y})}{[\sum_{i=1}^n (x_i - \bar{x})^2 \sum_{i=1}^n (y_i - \bar{y})^2]^{1/2}} \;\;\;\;\;\;\; (6)$$ and \(\underline{a}\) is defined by: $$\underline{a} = \frac{\underline{m}^T V^{-1}}{[\underline{m}^T V^{-1} V^{-1} \underline{m}]^{1/2}} \;\;\;\;\;\; (7)$$ where \(^T\) denotes the transpose operator, and \(\underline{m}\) is the vector of expected values and \(V\) is the variance-covariance matrix of the order statistics of a random sample of size \(N\) from a standard normal distribution. That is, the values of \(\underline{a}\) are the expected values of the standard normal order statistics weighted by their variance-covariance matrix, and normalized so that $$\underline{a}^T \underline{a} = 1 \;\;\;\;\;\; (8)$$

Computing Shapiro-Wilk W-Statistic for Singly Censored Data
For the case of singly censored data, following Smith and Bain (1976) and Verrill and Johnson (1988), Royston (1993) generalizes the Shapiro-Wilk \(W\)-statistic to: $$W = r(\underline{a}_{\Delta}, \underline{x}_{\Delta})^2 \;\;\;\;\;\; (9)$$ where for left singly-censored data: $$\underline{a}_{\Delta} = (a_{c+1}, a_{c+2}, \ldots, a_{N}) \;\;\;\;\;\; (10)$$ $$\underline{x}_{\Delta} = (x_{(c+1)}, x_{(c+2)}, \ldots, x_{(N)}) \;\;\;\;\;\; (11)$$ and for right singly-censored data: $$\underline{a}_{\Delta} = (a_1, a_2, \ldots, a_n) \;\;\;\;\;\; (12)$$ $$\underline{x}_{\Delta} = (x_{(1)}, x_{(2)}, \ldots, x_{(n)}) \;\;\;\;\;\; (13)$$ Just like the function gofTest, when test="sw", the function gofTestCensored uses Royston's (1992a) approximation for the coefficients \(\underline{a}\) (see the help file for gofTest).

Computing P-Values for the Shapiro-Wilk Test
Verrill and Johnson (1988) show that the asymptotic distribution of the statistic in Equation (9) above is normal, but the rate of convergence is “surprisingly slow” even for complete samples. They provide a table of empirical percentiles of the distribution for the \(W\)-statistic shown in Equation (9) above for several sample sizes and percentages of censoring.

Based on the tables given in Verrill and Johnson (1988), Royston (1993) approximated the 90'th, 95'th, and 99'th percentiles of the distribution of the z-statistic computed from the \(W\)-statistic. (The distribution of this z-statistic is assumed to be normal, but not necessarily a standard normal.) Denote these percentiles by \(Z_{0.90}\), \(Z_{0.95}\), and \(Z_{0.99}\). The true mean and standard deviation of the z-statistic are estimated by the intercept and slope, respectively, from the linear regression of \(Z_{\alpha}\) on \(\Phi^{-1}(\alpha)\) for \(\alpha\) = 0.9, 0.95, and 0.99, where \(\Phi\) denotes the cumulative distribution function of the standard normal distribution. The p-value associated with this test is then computed as: $$p = 1 - \Phi(\frac{z - \mu_z}{\sigma_z}) \;\;\;\;\;\; (14)$$ Note: Verrill and Johnson (1988) produced their tables based on Type II censoring. Royston's (1993) approximation to the p-value of these tests, however, should be fairly accurate for Type I censored data as well.

Testing Goodness-of-Fit for Any Continuous Distribution
The function gofTestCensored extends the Shapiro-Wilk test that accounts for censoring to test for goodness-of-fit for any continuous distribution by using the idea of Chen and Balakrishnan (1995), who proposed a general purpose approximate goodness-of-fit test based on the Cramer-von Mises or Anderson-Darling goodness-of-fit tests for normality. The function gofTestCensored modifies the approach of Chen and Balakrishnan (1995) by using the same first 2 steps, and then applying the Shapiro-Wilk test that accounts for censoring:

  1. Let \(\underline{x} = x_1, x_2, \ldots, x_n\) denote the vector of \(n\) ordered observations, ignoring censoring status. Compute cumulative probabilities for each \(x_i\) based on the cumulative distribution function for the hypothesized distribution. That is, compute \(p_i = F(x_i, \hat{\theta})\), where \(F(x, \theta)\) denotes the hypothesized cumulative distribution function with parameter(s) \(\theta\), and \(\hat{\theta}\) denotes the estimated parameter(s) using an estimation method that accounts for censoring (e.g., assuming a Gamma distribution with alternative parameterization, call the function link{egammaAltCensored}).

  2. Compute standard normal deviates based on the computed cumulative probabilities:
    \(y_i = \Phi^{-1}(p_i)\)

  3. Perform the Shapiro-Wilk goodness-of-fit test (that accounts for censoring) on the \(y_i\)'s.

Shapiro-Francia Goodness-of-Fit Test (test="sf")
Equation (15) in the help file for gofTest shows that for the complete ordered data \(\underline{x}\), the Shapiro-Francia \(W'\)-statistic is the same as the squared Pearson correlation coefficient associated with a normal probability plot.

Computing Shapiro-Francia W'-Statistic for Censored Data
For the case of singly censored data, following Smith and Bain (1976) and Verrill and Johnson (1988), Royston (1993) extends the computation of the Weisberg-Bingham Approximation to the \(W'\)-statistic to the case of singly censored data: $$\tilde{W}' = r(\underline{c}_{\Delta}, \underline{x}_{\Delta})^2 \;\;\;\;\;\; (14)$$ where for left singly-censored data: $$\underline{c}_{\Delta} = (c_{c+1}, c_{c+2}, \ldots, c_{N}) \;\;\;\;\;\; (15)$$ $$\underline{x}_{\Delta} = (x_{(c+1)}, x_{(c+2)}, \ldots, x_{(N)}) \;\;\;\;\;\; (16)$$ and for right singly-censored data: $$\underline{a}_{\Delta} = (a_1, a_2, \ldots, a_n) \;\;\;\;\;\; (17)$$ $$\underline{x}_{\Delta} = (x_{(1)}, x_{(2)}, \ldots, x_{(n)}) \;\;\;\;\;\; (18)$$ and \(\underline{c}\) is defined as: $$\underline{c} = \frac{\underline{\tilde{m}}}{[\underline{\tilde{m}}' \underline{\tilde{m}}]^{1/2}} \;\;\;\;\;\; (19)$$ where $$\tilde{m}_i = \Phi^{-1}(\frac{i - (3/8)}{n + (1/4)}) \;\;\;\;\;\; (20)$$ and \(\Phi\) denotes the standard normal cdf. Note: Do not confuse the elements of the vector \(\underline{c}\) with the scalar \(c\) which denotes the number of censored observations. We use \(\underline{c}\) here to be consistent with the notation in the help file for gofTest.

Just like the function gofTest, when test="sf", the function gofTestCensored uses Royston's (1992a) approximation for the coefficients \(\underline{c}\) (see the help file for gofTest).

In general, the Shapiro-Francia test statistic can be extended to multiply censored data using Equation (14) with \(\underline{c}_{\Delta}\) defined as the orderd values of \(c_i\) associated with uncensored observations, and \(\underline{x}_{\Delta}\) defined as the ordered values of \(x_i\) associated with uncensored observations: $$\underline{c}_{\Delta} = \cup_{i \in \Omega} \;\; c_{(i)} \;\;\;\;\;\; (21)$$ $$\underline{x}_{\Delta} = \cup_{i \in \Omega} \;\; x_{(i)} \;\;\;\;\;\; (22)$$ and where the plotting positions in Equation (20) are replaced with any of the plotting positions available in ppointsCensored (see the description for the argument prob.method).

Computing P-Values for the Shapiro-Francia Test
Verrill and Johnson (1988) show that the asymptotic distribution of the statistic in Equation (14) above is normal, but the rate of convergence is “surprisingly slow” even for complete samples. They provide a table of empirical percentiles of the distribution for the \(\tilde{W}'\)-statistic shown in Equation (14) above for several sample sizes and percentages of censoring.

As for the Shapiro-Wilk test, based on the tables given in Verrill and Johnson (1988), Royston (1993) approximated the 90'th, 95'th, and 99'th percentiles of the distribution of the z-statistic computed from the \(\tilde{W}'\)-statistic. (The distribution of this z-statistic is assumed to be normal, but not necessarily a standard normal.) Denote these percentiles by \(Z_{0.90}\), \(Z_{0.95}\), and \(Z_{0.99}\). The true mean and standard deviation of the z-statistic are estimated by the intercept and slope, respectively, from the linear regression of \(Z_{\alpha}\) on \(\Phi^{-1}(\alpha)\) for \(\alpha\) = 0.9, 0.95, and 0.99, where \(\Phi\) denotes the cumulative distribution function of the standard normal distribution. The p-value associated with this test is then computed as: $$p = 1 - \Phi(\frac{z - \mu_z}{\sigma_z}) \;\;\;\;\;\; (23)$$ Note: Verrill and Johnson (1988) produced their tables based on Type II censoring. Royston's (1993) approximation to the p-value of these tests, however, should be fairly accurate for Type I censored data as well, although this is an area that requires further investigation.

Testing Goodness-of-Fit for Any Continuous Distribution
The function gofTestCensored extends the Shapiro-Francia test that accounts for censoring to test for goodness-of-fit for any continuous distribution by using the idea of Chen and Balakrishnan (1995), who proposed a general purpose approximate goodness-of-fit test based on the Cramer-von Mises or Anderson-Darling goodness-of-fit tests for normality. The function gofTestCensored modifies the approach of Chen and Balakrishnan (1995) by using the same first 2 steps, and then applying the Shapiro-Francia test that accounts for censoring:

  1. Let \(\underline{x} = x_1, x_2, \ldots, x_n\) denote the vector of \(n\) ordered observations, ignoring censoring status. Compute cumulative probabilities for each \(x_i\) based on the cumulative distribution function for the hypothesized distribution. That is, compute \(p_i = F(x_i, \hat{\theta})\) where \(F(x, \theta)\) denotes the hypothesized cumulative distribution function with parameter(s) \(\theta\), and \(\hat{\theta}\) denotes the estimated parameter(s) using an estimation method that accounts for censoring (e.g., assuming a Gamma distribution with alternative parameterization, call the function link{egammaAltCensored}).

  2. Compute standard normal deviates based on the computed cumulative probabilities:
    \(y_i = \Phi^{-1}(p_i)\)

  3. Perform the Shapiro-Francia goodness-of-fit test (that accounts for censoring) on the \(y_i\)'s.

Probability Plot Correlation Coefficient (PPCC) Goodness-of-Fit Test (test="ppcc")
The function gofTestCensored computes the PPCC test statistic using Blom plotting positions. It can be shown that the square of this statistic is equivalent to the Weisberg-Bingham Approximation to the Shapiro-Francia \(W'\)-test (Weisberg and Bingham, 1975; Royston, 1993). Thus the PPCC goodness-of-fit test is equivalent to the Shapiro-Francia goodness-of-fit test.

Value

a list of class "gofCensored" containing the results of the goodness-of-fit test. See the help files for gofCensored.object for details.

References

Birnbaum, Z.W., and F.H. Tingey. (1951). One-Sided Confidence Contours for Probability Distribution Functions. Annals of Mathematical Statistics 22, 592-596.

Blom, G. (1958). Statistical Estimates and Transformed Beta Variables. John Wiley and Sons, New York.

Conover, W.J. (1980). Practical Nonparametric Statistics. Second Edition. John Wiley and Sons, New York.

Dallal, G.E., and L. Wilkinson. (1986). An Analytic Approximation to the Distribution of Lilliefor's Test for Normality. The American Statistician 40, 294-296.

D'Agostino, R.B. (1970). Transformation to Normality of the Null Distribution of \(g1\). Biometrika 57, 679-681.

D'Agostino, R.B. (1971). An Omnibus Test of Normality for Moderate and Large Size Samples. Biometrika 58, 341-348.

D'Agostino, R.B. (1986b). Tests for the Normal Distribution. In: D'Agostino, R.B., and M.A. Stephens, eds. Goodness-of Fit Techniques. Marcel Dekker, New York.

D'Agostino, R.B., and E.S. Pearson (1973). Tests for Departures from Normality. Empirical Results for the Distributions of \(b2\) and \(\sqrt{b1}\). Biometrika 60(3), 613-622.

D'Agostino, R.B., and G.L. Tietjen (1973). Approaches to the Null Distribution of \(\sqrt{b1}\). Biometrika 60(1), 169-173.

Fisher, R.A. (1950). Statistical Methods for Research Workers. 11'th Edition. Hafner Publishing Company, New York, pp.99-100.

Gibbons, R.D., D.K. Bhaumik, and S. Aryal. (2009). Statistical Methods for Groundwater Monitoring, Second Edition. John Wiley & Sons, Hoboken.

Kendall, M.G., and A. Stuart. (1991). The Advanced Theory of Statistics, Volume 2: Inference and Relationship. Fifth Edition. Oxford University Press, New York.

Royston, J.P. (1992a). Approximating the Shapiro-Wilk W-Test for Non-Normality. Statistics and Computing 2, 117-119.

Royston, J.P. (1992b). Estimation, Reference Ranges and Goodness of Fit for the Three-Parameter Log-Normal Distribution. Statistics in Medicine 11, 897-912.

Royston, J.P. (1992c). A Pocket-Calculator Algorithm for the Shapiro-Francia Test of Non-Normality: An Application to Medicine. Statistics in Medicine 12, 181-184.

Royston, P. (1993). A Toolkit for Testing for Non-Normality in Complete and Censored Samples. The Statistician 42, 37-43.

Ryan, T., and B. Joiner. (1973). Normal Probability Plots and Tests for Normality. Technical Report, Pennsylvannia State University, Department of Statistics.

Shapiro, S.S., and R.S. Francia. (1972). An Approximate Analysis of Variance Test for Normality. Journal of the American Statistical Association 67(337), 215-219.

Shapiro, S.S., and M.B. Wilk. (1965). An Analysis of Variance Test for Normality (Complete Samples). Biometrika 52, 591-611.

Verrill, S., and R.A. Johnson. (1987). The Asymptotic Equivalence of Some Modified Shapiro-Wilk Statistics – Complete and Censored Sample Cases. The Annals of Statistics 15(1), 413-419.

Verrill, S., and R.A. Johnson. (1988). Tables and Large-Sample Distribution Theory for Censored-Data Correlation Statistics for Testing Normality. Journal of the American Statistical Association 83, 1192-1197.

Weisberg, S., and C. Bingham. (1975). An Approximate Analysis of Variance Test for Non-Normality Suitable for Machine Calculation. Technometrics 17, 133-134.

Author

Steven P. Millard (EnvStats@ProbStatInfo.com)

Note

The Shapiro-Wilk test (Shapiro and Wilk, 1965) and the Shapiro-Francia test (Shapiro and Francia, 1972) are probably the two most commonly used hypothesis tests to test departures from normality. The Shapiro-Wilk test is most powerful at detecting short-tailed (platykurtic) and skewed distributions, and least powerful against symmetric, moderately long-tailed (leptokurtic) distributions. Conversely, the Shapiro-Francia test is more powerful against symmetric long-tailed distributions and less powerful against short-tailed distributions (Royston, 1992b; 1993).

In practice, almost any goodness-of-fit test will not reject the null hypothesis if the number of observations is relatively small. Conversely, almost any goodness-of-fit test will reject the null hypothesis if the number of observations is very large, since “real” data are never distributed according to any theoretical distribution (Conover, 1980, p.367). For most cases, however, the distribution of “real” data is close enough to some theoretical distribution that fairly accurate results may be provided by assuming that particular theoretical distribution. One way to asses the goodness of the fit is to use goodness-of-fit tests. Another way is to look at quantile-quantile (Q-Q) plots (see qqPlotCensored).

Examples

  # Generate 30 observations from a gamma distribution with 
  # parameters mean=10 and cv=1 and then censor observations less than 5.
  # Then test the hypothesis that these data came from a gamma 
  # distribution using the Shapiro-Wilk test.
  #
  # The p-value for the complete data is p = 0.86, while  
  # the p-value for the censored data is p = 0.52.
  # (Note:  the call to set.seed lets you reproduce this example.)

  set.seed(598)

  dat <- sort(rgammaAlt(30, mean = 10, cv = 1))
  dat
#>  [1]  0.5313509  1.4741833  1.9936208  2.7980636  3.4509840  3.7987348
#>  [7]  4.5542952  5.5207531  5.5253596  5.7177872  5.7513827  9.1086375
#> [13]  9.8444090 10.6247123 10.9304922 11.7925398 13.3432689 13.9562777
#> [19] 14.6029065 15.0563342 15.8730642 16.0039936 16.6910715 17.0288922
#> [25] 17.8507891 19.1105522 20.2657141 26.3815970 30.2912797 42.8726101
  # [1]  0.5313509  1.4741833  1.9936208  2.7980636  3.4509840
  # [6]  3.7987348  4.5542952  5.5207531  5.5253596  5.7177872
  #[11]  5.7513827  9.1086375  9.8444090 10.6247123 10.9304922
  #[16] 11.7925398 13.3432689 13.9562777 14.6029065 15.0563342
  #[21] 15.8730642 16.0039936 16.6910715 17.0288922 17.8507891
  #[26] 19.1105522 20.2657141 26.3815970 30.2912797 42.8726101

  dat.censored <- dat
  censored <- dat.censored < 5
  dat.censored[censored] <- 5

  # Results for complete data:
  #---------------------------
  gofTest(dat, test = "sw", dist = "gammaAlt")
#> 
#> Results of Goodness-of-Fit Test
#> -------------------------------
#> 
#> Test Method:                     Shapiro-Wilk GOF Based on
#>                                  Chen & Balakrisnan (1995)
#> 
#> Hypothesized Distribution:       Gamma
#> 
#> Estimated Parameter(s):          mean = 12.4248552
#>                                  cv   =  0.7901752
#> 
#> Estimation Method:               MLE
#> 
#> Data:                            dat
#> 
#> Sample Size:                     30
#> 
#> Test Statistic:                  W = 0.9814711
#> 
#> Test Statistic Parameter:        n = 30
#> 
#> P-value:                         0.8631802
#> 
#> Alternative Hypothesis:          True cdf does not equal the
#>                                  Gamma Distribution.

  #Results of Goodness-of-Fit Test
  #-------------------------------
  #
  #Test Method:                     Shapiro-Wilk GOF Based on 
  #                                 Chen & Balakrisnan (1995)
  #
  #Hypothesized Distribution:       Gamma
  #
  #Estimated Parameter(s):          mean = 12.4248552
  #                                 cv   =  0.7901752
  #
  #Estimation Method:               MLE
  #
  #Data:                            dat
  #
  #Sample Size:                     30
  #
  #Test Statistic:                  W = 0.981471
  #
  #Test Statistic Parameter:        n = 30
  #
  #P-value:                         0.8631802
  #
  #Alternative Hypothesis:          True cdf does not equal the
  #                                 Gamma Distribution.


  # Results for censored data:
  #---------------------------
  gof.list <- gofTestCensored(dat.censored, censored, test = "sw", 
    distribution = "gammaAlt")
  gof.list  
#> 
#> Results of Goodness-of-Fit Test
#> Based on Type I Censored Data
#> -------------------------------
#> 
#> Test Method:                     Shapiro-Wilk GOF
#>                                  (Singly Censored Data)
#>                                  Based on Chen & Balakrisnan (1995)
#> 
#> Hypothesized Distribution:       Gamma
#> 
#> Censoring Side:                  left
#> 
#> Censoring Level(s):              5 
#> 
#> Estimated Parameter(s):          mean = 12.4911448
#>                                  cv   =  0.7617343
#> 
#> Estimation Method:               MLE
#> 
#> Data:                            dat.censored
#> 
#> Censoring Variable:              censored
#> 
#> Sample Size:                     30
#> 
#> Percent Censored:                23.33333%
#> 
#> Test Statistic:                  W = 0.9613711
#> 
#> Test Statistic Parameters:       N     = 30.0000000
#>                                  DELTA =  0.2333333
#> 
#> P-value:                         0.522329
#> 
#> Alternative Hypothesis:          True cdf does not equal the
#>                                  Gamma Distribution.

  #Results of Goodness-of-Fit Test
  #Based on Type I Censored Data
  #-------------------------------
  #
  #Test Method:                     Shapiro-Wilk GOF
  #                                 (Singly Censored Data)
  #                                 Based on Chen & Balakrisnan (1995)
  #
  #Hypothesized Distribution:       Gamma
  #
  #Censoring Side:                  left
  #
  #Censoring Level(s):              5 
  #
  #Estimated Parameter(s):          mean = 12.4911448
  #                                 cv   =  0.7617343
  #
  #Estimation Method:               MLE
  #
  #Data:                            dat.censored
  #
  #Censoring Variable:              censored
  #
  #Sample Size:                     30
  #
  #Percent Censored:                23.3%
  #
  #Test Statistic:                  W = 0.9613711
  #
  #Test Statistic Parameters:       N     = 30.0000000
  #                                 DELTA =  0.2333333
  #
  #P-value:                         0.522329
  #
  #Alternative Hypothesis:          True cdf does not equal the
  #                                 Gamma Distribution.

  # Plot the results for the censored data
  #---------------------------------------
  dev.new()
  plot(gof.list)

  #==========

  # Continue the above example, but now test the hypothesis that 
  # these data came from a lognormal distribution 
  # (alternative parameterization) using the Shapiro-Wilk test.
  #
  # The p-value for the complete data is p = 0.056, while  
  # the p-value for the censored data is p = 0.11.

  # Results for complete data:
  #---------------------------
  gofTest(dat, test = "sw", dist = "lnormAlt")
#> 
#> Results of Goodness-of-Fit Test
#> -------------------------------
#> 
#> Test Method:                     Shapiro-Wilk GOF
#> 
#> Hypothesized Distribution:       Lognormal
#> 
#> Estimated Parameter(s):          mean = 13.757239
#>                                  cv   =  1.148872
#> 
#> Estimation Method:               mvue
#> 
#> Data:                            dat
#> 
#> Sample Size:                     30
#> 
#> Test Statistic:                  W = 0.9322226
#> 
#> Test Statistic Parameter:        n = 30
#> 
#> P-value:                         0.05626823
#> 
#> Alternative Hypothesis:          True cdf does not equal the
#>                                  Lognormal Distribution.

  #Results of Goodness-of-Fit Test
  #-------------------------------
  #
  #Test Method:                     Shapiro-Wilk GOF
  #
  #Hypothesized Distribution:       Lognormal
  #
  #Estimated Parameter(s):          mean = 13.757239
  #                                 cv   =  1.148872
  #
  #Estimation Method:               mvue
  #
  #Data:                            dat
  #
  #Sample Size:                     30
  #
  #Test Statistic:                  W = 0.9322226
  #
  #Test Statistic Parameter:        n = 30
  #
  #P-value:                         0.05626823
  #
  #Alternative Hypothesis:          True cdf does not equal the
  #                                 Lognormal Distribution.


  # Results for censored data:
  #---------------------------
  gof.list <- gofTestCensored(dat.censored, censored, test = "sw", 
    distribution = "lnormAlt")
  gof.list  
#> 
#> Results of Goodness-of-Fit Test
#> Based on Type I Censored Data
#> -------------------------------
#> 
#> Test Method:                     Shapiro-Wilk GOF
#>                                  (Singly Censored Data)
#> 
#> Hypothesized Distribution:       Lognormal
#> 
#> Censoring Side:                  left
#> 
#> Censoring Level(s):              5 
#> 
#> Estimated Parameter(s):          mean = 13.0382221
#>                                  cv   =  0.9129512
#> 
#> Estimation Method:               MLE
#> 
#> Data:                            dat.censored
#> 
#> Censoring Variable:              censored
#> 
#> Sample Size:                     30
#> 
#> Percent Censored:                23.33333%
#> 
#> Test Statistic:                  W = 0.9292406
#> 
#> Test Statistic Parameters:       N     = 30.0000000
#>                                  DELTA =  0.2333333
#> 
#> P-value:                         0.114511
#> 
#> Alternative Hypothesis:          True cdf does not equal the
#>                                  Lognormal Distribution.

  #Results of Goodness-of-Fit Test
  #Based on Type I Censored Data
  #-------------------------------
  #
  #Test Method:                     Shapiro-Wilk GOF
  #                                 (Singly Censored Data)
  #
  #Hypothesized Distribution:       Lognormal
  #
  #Censoring Side:                  left
  #
  #Censoring Level(s):              5 
  #
  #Estimated Parameter(s):          mean = 13.0382221
  #                                 cv   =  0.9129512
  #
  #Estimation Method:               MLE
  #
  #Data:                            dat.censored
  #
  #Censoring Variable:              censored
  #
  #Sample Size:                     30
  #
  #Percent Censored:                23.3%
  #
  #Test Statistic:                  W = 0.9292406
  #
  #Test Statistic Parameters:       N     = 30.0000000
  #                                 DELTA =  0.2333333
  #
  #P-value:                         0.114511
  #
  #Alternative Hypothesis:          True cdf does not equal the
  #                                 Lognormal Distribution.

  # Plot the results for the censored data
  #---------------------------------------
  dev.new()
  plot(gof.list)

  #----------

  # Redo the above example, but specify the quasi-minimum variance 
  # unbiased estimator of the mean.  Note that the method of 
  # estimating the parameters has no effect on the goodness-of-fit 
  # test (see the DETAILS section above).

  gofTestCensored(dat.censored, censored, test = "sw", 
    distribution = "lnormAlt", est.arg.list = list(method = "qmvue"))
#> 
#> Results of Goodness-of-Fit Test
#> Based on Type I Censored Data
#> -------------------------------
#> 
#> Test Method:                     Shapiro-Wilk GOF
#>                                  (Singly Censored Data)
#> 
#> Hypothesized Distribution:       Lognormal
#> 
#> Censoring Side:                  left
#> 
#> Censoring Level(s):              5 
#> 
#> Estimated Parameter(s):          mean = 12.8722749
#>                                  cv   =  0.8712549
#> 
#> Estimation Method:               Quasi-MVUE
#> 
#> Data:                            dat.censored
#> 
#> Censoring Variable:              censored
#> 
#> Sample Size:                     30
#> 
#> Percent Censored:                23.33333%
#> 
#> Test Statistic:                  W = 0.9292406
#> 
#> Test Statistic Parameters:       N     = 30.0000000
#>                                  DELTA =  0.2333333
#> 
#> P-value:                         0.114511
#> 
#> Alternative Hypothesis:          True cdf does not equal the
#>                                  Lognormal Distribution.

  #Results of Goodness-of-Fit Test
  #Based on Type I Censored Data
  #-------------------------------
  #
  #Test Method:                     Shapiro-Wilk GOF
  #                                 (Singly Censored Data)
  #
  #Hypothesized Distribution:       Lognormal
  #
  #Censoring Side:                  left
  #
  #Censoring Level(s):              5 
  #
  #Estimated Parameter(s):          mean = 12.8722749
  #                                 cv   =  0.8712549
  #
  #Estimation Method:               Quasi-MVUE
  #
  #Data:                            dat.censored
  #
  #Censoring Variable:              censored
  #
  #Sample Size:                     30
  #
  #Percent Censored:                23.3%
  #
  #Test Statistic:                  W = 0.9292406
  #
  #Test Statistic Parameters:       N     = 30.0000000
  #                                 DELTA =  0.2333333
  #
  #P-value:                         0.114511
  #
  #Alternative Hypothesis:          True cdf does not equal the
  #                                 Lognormal Distribution.

  #----------
  # Clean up

  rm(dat, dat.censored, censored, gof.list)
  graphics.off()

  #==========

  # Check the assumption that the silver data stored in Helsel.Cohn.88.silver.df 
  # follows a lognormal distribution and plot the goodness-of-fit test results.  
  # Note that the small p-value and the shape of the Q-Q plot 
  # (an inverted S-shape) suggests that the log transformation is not quite strong 
  # enough to "bring in" the tails (i.e., the log-transformed silver data has tails 
  # that are slightly too long relative to a normal distribution).  
  # Helsel and Cohn (1988, p.2002) note that the gross outlier of 560 mg/L tends to 
  # make the shape of the data resemble a gamma distribution.

  dum.list <- with(Helsel.Cohn.88.silver.df, 
    gofTestCensored(Ag, Censored, test = "sf", dist = "lnorm"))

  dum.list
#> 
#> Results of Goodness-of-Fit Test
#> Based on Type I Censored Data
#> -------------------------------
#> 
#> Test Method:                     Shapiro-Francia GOF
#>                                  (Multiply Censored Data)
#> 
#> Hypothesized Distribution:       Lognormal
#> 
#> Censoring Side:                  left
#> 
#> Censoring Level(s):               0.1  0.2  0.3  0.5  1.0  2.0  2.5  5.0  6.0 10.0 20.0 25.0 
#> 
#> Estimated Parameter(s):          meanlog = -1.040572
#>                                  sdlog   =  2.354847
#> 
#> Estimation Method:               MLE
#> 
#> Data:                            Ag
#> 
#> Censoring Variable:              Censored
#> 
#> Sample Size:                     56
#> 
#> Percent Censored:                60.71429%
#> 
#> Test Statistic:                  W = 0.8957198
#> 
#> Test Statistic Parameters:       N     = 56.0000000
#>                                  DELTA =  0.6071429
#> 
#> P-value:                         0.03490314
#> 
#> Alternative Hypothesis:          True cdf does not equal the
#>                                  Lognormal Distribution.
  #Results of Goodness-of-Fit Test
  #Based on Type I Censored Data
  #-------------------------------
  #
  #Test Method:                     Shapiro-Francia GOF
  #                                 (Multiply Censored Data)
  #
  #Hypothesized Distribution:       Lognormal
  #
  #Censoring Side:                  left
  #
  #Censoring Level(s):               0.1  0.2  0.3  0.5  1.0  2.0  2.5  5.0  
  #                                  6.0 10.0 20.0 25.0 
  #
  #Estimated Parameter(s):          meanlog = -1.040572
  #                                 sdlog   =  2.354847
  #
  #Estimation Method:               MLE
  #
  #Data:                            Ag
  #
  #Censoring Variable:              Censored
  #
  #Sample Size:                     56
  #
  #Percent Censored:                60.7%
  #
  #Test Statistic:                  W = 0.8957198
  #
  #Test Statistic Parameters:       N     = 56.0000000
  #                                 DELTA =  0.6071429
  #
  #P-value:                         0.03490314
  #
  #Alternative Hypothesis:          True cdf does not equal the
  #                                 Lognormal Distribution.

  dev.new()
  plot(dum.list)
#> Warning: Cannot construct histogram for multiply censored data when complete observations are between censoring levels

  #----------
  
  # Clean up
  #---------

  rm(dum.list)
  graphics.off()

  #==========

  # Chapter 15 of USEPA (2009) gives several examples of looking 
  # at normal Q-Q plots and estimating the mean and standard deviation 
  # for manganese concentrations (ppb) in groundwater at five background wells. 
  # In EnvStats these data are stored in the data frame 
  # EPA.09.Ex.15.1.manganese.df.

  # Here we will test whether the data appear to come from a normal 
  # distribution, then we will test to see whether they appear to come 
  # from a lognormal distribution.
  #--------------------------------------------------------------------


  # First look at the data:
  #-----------------------

  EPA.09.Ex.15.1.manganese.df
#>    Sample   Well Manganese.Orig.ppb Manganese.ppb Censored
#> 1       1 Well.1                 <5           5.0     TRUE
#> 2       2 Well.1               12.1          12.1    FALSE
#> 3       3 Well.1               16.9          16.9    FALSE
#> 4       4 Well.1               21.6          21.6    FALSE
#> 5       5 Well.1                 <2           2.0     TRUE
#> 6       1 Well.2                 <5           5.0     TRUE
#> 7       2 Well.2                7.7           7.7    FALSE
#> 8       3 Well.2               53.6          53.6    FALSE
#> 9       4 Well.2                9.5           9.5    FALSE
#> 10      5 Well.2               45.9          45.9    FALSE
#> 11      1 Well.3                 <5           5.0     TRUE
#> 12      2 Well.3                5.3           5.3    FALSE
#> 13      3 Well.3               12.6          12.6    FALSE
#> 14      4 Well.3              106.3         106.3    FALSE
#> 15      5 Well.3               34.5          34.5    FALSE
#> 16      1 Well.4                6.3           6.3    FALSE
#> 17      2 Well.4               11.9          11.9    FALSE
#> 18      3 Well.4                 10          10.0    FALSE
#> 19      4 Well.4                 <2           2.0     TRUE
#> 20      5 Well.4               77.2          77.2    FALSE
#> 21      1 Well.5               17.9          17.9    FALSE
#> 22      2 Well.5               22.7          22.7    FALSE
#> 23      3 Well.5                3.3           3.3    FALSE
#> 24      4 Well.5                8.4           8.4    FALSE
#> 25      5 Well.5                 <2           2.0     TRUE

  #   Sample   Well Manganese.Orig.ppb Manganese.ppb Censored
  #1       1 Well.1                 <5           5.0     TRUE
  #2       2 Well.1               12.1          12.1    FALSE
  #3       3 Well.1               16.9          16.9    FALSE
  #...
  #23      3 Well.5                3.3           3.3    FALSE
  #24      4 Well.5                8.4           8.4    FALSE
  #25      5 Well.5                 <2           2.0     TRUE

  longToWide(EPA.09.Ex.15.1.manganese.df, 
    "Manganese.Orig.ppb", "Sample", "Well",
    paste.row.name = TRUE)  
#>          Well.1 Well.2 Well.3 Well.4 Well.5
#> Sample.1     <5     <5     <5    6.3   17.9
#> Sample.2   12.1    7.7    5.3   11.9   22.7
#> Sample.3   16.9   53.6   12.6     10    3.3
#> Sample.4   21.6    9.5  106.3     <2    8.4
#> Sample.5     <2   45.9   34.5   77.2     <2

  #         Well.1 Well.2 Well.3 Well.4 Well.5
  #Sample.1     <5     <5     <5    6.3   17.9
  #Sample.2   12.1    7.7    5.3   11.9   22.7
  #Sample.3   16.9   53.6   12.6     10    3.3
  #Sample.4   21.6    9.5  106.3     <2    8.4
  #Sample.5     <2   45.9   34.5   77.2     <2


  # Now test whether the data appear to come from 
  # a normal distribution.  Note that these data 
  # are multiply censored, so we'll use the 
  # Shapiro-Francia test.  
  #----------------------------------------------

  gof.normal <- with(EPA.09.Ex.15.1.manganese.df,
    gofTestCensored(Manganese.ppb, Censored, test = "sf"))

  gof.normal
#> 
#> Results of Goodness-of-Fit Test
#> Based on Type I Censored Data
#> -------------------------------
#> 
#> Test Method:                     Shapiro-Francia GOF
#>                                  (Multiply Censored Data)
#> 
#> Hypothesized Distribution:       Normal
#> 
#> Censoring Side:                  left
#> 
#> Censoring Level(s):              2 5 
#> 
#> Estimated Parameter(s):          mean = 15.23508
#>                                  sd   = 30.62812
#> 
#> Estimation Method:               MLE
#> 
#> Data:                            Manganese.ppb
#> 
#> Censoring Variable:              Censored
#> 
#> Sample Size:                     25
#> 
#> Percent Censored:                24%
#> 
#> Test Statistic:                  W = 0.8368016
#> 
#> Test Statistic Parameters:       N     = 25.00
#>                                  DELTA =  0.24
#> 
#> P-value:                         0.004662658
#> 
#> Alternative Hypothesis:          True cdf does not equal the
#>                                  Normal Distribution.

  #Results of Goodness-of-Fit Test
  #Based on Type I Censored Data
  #-------------------------------
  #
  #Test Method:                     Shapiro-Francia GOF
  #                                 (Multiply Censored Data)
  #
  #Hypothesized Distribution:       Normal
  #
  #Censoring Side:                  left
  #
  #Censoring Level(s):              2 5 
  #
  #Estimated Parameter(s):          mean = 15.23508
  #                                 sd   = 30.62812
  #
  #Estimation Method:               MLE
  #
  #Data:                            Manganese.ppb
  #
  #Censoring Variable:              Censored
  #
  #Sample Size:                     25
  #
  #Percent Censored:                24%
  #
  #Test Statistic:                  W = 0.8368016
  #
  #Test Statistic Parameters:       N     = 25.00
  #                                 DELTA =  0.24
  #
  #P-value:                         0.004662658
  #
  #Alternative Hypothesis:          True cdf does not equal the
  #                                 Normal Distribution.

  # Plot the results:
  #------------------

  dev.new()
  plot(gof.normal)
#> Warning: Cannot construct histogram for multiply censored data when complete observations are between censoring levels

  #----------

  # Now test to see whether the data appear to come from 
  # a lognormal distribuiton.
  #-----------------------------------------------------

  gof.lognormal <- with(EPA.09.Ex.15.1.manganese.df,
    gofTestCensored(Manganese.ppb, Censored, test = "sf", 
    distribution = "lnorm"))

  gof.lognormal
#> 
#> Results of Goodness-of-Fit Test
#> Based on Type I Censored Data
#> -------------------------------
#> 
#> Test Method:                     Shapiro-Francia GOF
#>                                  (Multiply Censored Data)
#> 
#> Hypothesized Distribution:       Lognormal
#> 
#> Censoring Side:                  left
#> 
#> Censoring Level(s):              2 5 
#> 
#> Estimated Parameter(s):          meanlog = 2.215905
#>                                  sdlog   = 1.356291
#> 
#> Estimation Method:               MLE
#> 
#> Data:                            Manganese.ppb
#> 
#> Censoring Variable:              Censored
#> 
#> Sample Size:                     25
#> 
#> Percent Censored:                24%
#> 
#> Test Statistic:                  W = 0.9864426
#> 
#> Test Statistic Parameters:       N     = 25.00
#>                                  DELTA =  0.24
#> 
#> P-value:                         0.9767731
#> 
#> Alternative Hypothesis:          True cdf does not equal the
#>                                  Lognormal Distribution.

  #Results of Goodness-of-Fit Test
  #Based on Type I Censored Data
  #-------------------------------
  #
  #Test Method:                     Shapiro-Francia GOF
  #                                 (Multiply Censored Data)
  #
  #Hypothesized Distribution:       Lognormal
  #
  #Censoring Side:                  left
  #
  #Censoring Level(s):              2 5 
  #
  #Estimated Parameter(s):          meanlog = 2.215905
  #                                 sdlog   = 1.356291
  #
  #Estimation Method:               MLE
  #
  #Data:                            Manganese.ppb
  #
  #Censoring Variable:              Censored
  #
  #Sample Size:                     25
  #
  #Percent Censored:                24%
  #
  #Test Statistic:                  W = 0.9864426
  #
  #Test Statistic Parameters:       N     = 25.00
  #                                 DELTA =  0.24
  #
  #P-value:                         0.9767731
  #
  #Alternative Hypothesis:          True cdf does not equal the
  #                                 Lognormal Distribution.

  # Plot the results:
  #------------------

  dev.new()
  plot(gof.lognormal)
#> Warning: Cannot construct histogram for multiply censored data when complete observations are between censoring levels

  #----------

  # Clean up
  #---------

  rm(gof.normal, gof.lognormal)
  graphics.off()