# inst/scripts/Manual/Chapter08.R In EnvStats: Package for Environmental Statistics, Including US EPA Guidance

```#####
# File:     Chapter08.R
#
# Purpose:  Reproduce Examples in Chapter 8 of the book:
#
#           Millard, SP. (2013).  EnvStats: An R Package for
#             Environmental Statistics.  Springer-Verlag.
#
# Author:   Steven P. Millard
#
# Last
# Updated:  2013/08/27
#####

#######################################################################################

library(EnvStats)

############################################################

# 8.4.1 Quantile (Empirical CDF) Plots for Censored Data
#-----------------------------------------------------

# Figures 8.1 & 8.2
#------------------

windows()
with(Helsel.Cohn.88.silver.df,
ecdfPlotCensored(Ag, Censored,
xlab = expression(paste("Ag (", mu, "g/L)", sep = "")),
main = "Empirical CDF Plot for Silver\n(Method of Michael-Schucany)",
include.cen = TRUE))

windows()
with(Helsel.Cohn.88.silver.df,
ecdfPlotCensored(log.Ag, Censored,
xlab = expression(paste("log [ Ag (", mu, "g/L) ]", sep = "")),
main = "Empirical CDF Plot for Silver\n(Method of Michael-Schucany)",
include.cen = TRUE))

#######################################################################################

# 8.4.2 Comparing an Empirical and Hypothesized CDF
#--------------------------------------------------

# Figures 8.3 & 8.4
#------------------

windows()
with(Helsel.Cohn.88.silver.df,
cdfCompareCensored(Ag, Censored, distribution = "lnorm",
xlab = expression(paste("Ag (", mu, "g/L)", sep = ""))))

windows()
with(Helsel.Cohn.88.silver.df,
cdfCompareCensored(log(Ag), Censored, distribution="norm",
xlab = expression(paste("log [ Ag (", mu, "g/L) ]", sep = ""))))

#########################################################################################

# 8.4.3 Comparing Two Empirical CDFs
#-----------------------------------

# Figure 8.5
#-----------

attach(Millard.Deverel.88.df)
Cu.Alluvial <- Cu[Zone == "Alluvial.Fan"]
Cu.Alluvial.cen <- Cu.censored[Zone == "Alluvial.Fan"]
Cu.Basin <- Cu[Zone == "Basin.Trough"]
Cu.Basin.cen <- Cu.censored[Zone == "Basin.Trough"]

windows()
cdfCompareCensored(x = Cu.Alluvial, censored = Cu.Alluvial.cen,
y = Cu.Basin, y.censored = Cu.Basin.cen)

detach("Millard.Deverel.88.df")
rm(Cu.Alluvial, Cu.Alluvial.cen, Cu.Basin, Cu.Basin.cen)

#####################################################################################

# 8.4.4 Q-Q Plots for Censored Data
#----------------------------------

# Figure 8.6
#-----------

windows()
with(Modified.TcCB.df,
qqPlotCensored(TcCB, Censored, distribution = "lnorm",
add.line = TRUE, points.col = "blue"))

###################################################################################

# 8.4.5 Box-Cox Transformations for Censored Data
#------------------------------------------------

# Figure 8.7
#-----------

boxcox.list <- with(Modified.TcCB.df,
boxcoxCensored(TcCB, Censored))

windows()
plot(boxcox.list)

rm(boxcox.list)

#################################################################################

# 8.5 Estimating Distribution Parameters
#---------------------------------------

# 8.5.1 The Normal Distribution
#------------------------------

# Estimating Parameters of a Normal Distribution
#-----------------------------------------------

with(EPA.09.Ex.15.1.manganese.df,
elnormCensored(Manganese.ppb, Censored))

with(EPA.09.Ex.15.1.manganese.df,
elnormCensored(Manganese.ppb, Censored, method = "qq.reg"))

with(EPA.09.Ex.15.1.manganese.df,
elnormCensored(Manganese.ppb, Censored,
method = "impute.w.qq.reg", prob.method = "hirsch-stedinger"))

with(EPA.09.Ex.15.1.manganese.df,
elnormCensored(Manganese.ppb, Censored,
method = "impute.w.qq.reg", prob.method = "hirsch-stedinger",
plot.pos.con = 0))

# Confidence Intervals for the Mean of the Normal Distribution
#-------------------------------------------------------------

with(EPA.09.Ex.15.1.manganese.df,
elnormCensored(Manganese.ppb, Censored, ci = TRUE))

#################################################################################

# 8.5.2 The Lognormal Distribution
#---------------------------------

with(Modified.TcCB.df,
elnormAltCensored(TcCB, Censored, method = "mle", ci = TRUE))

#################################################################################

# 8.5.3 The Gamma Distribution
#-----------------------------

with(Modified.TcCB.df,
egammaAltCensored(TcCB, Censored, ci = TRUE, ci.method = "profile.likelihood"))

#################################################################################

# 8.5.4 Estimating the Mean Nonparametrically
#--------------------------------------------

with(EPA.09.Ex.15.1.manganese.df,
enparCensored(log(Manganese.ppb), Censored, ci = TRUE))

with(EPA.09.Ex.15.1.manganese.df,
enparCensored(Manganese.ppb, Censored, ci = TRUE))

#################################################################################

# 8.6 Estimating Distribution Quantiles
#--------------------------------------

# 8.6.1 Parametric Estimates of Quantiles
#----------------------------------------

#--------------------------------------
# The Normal and Lognormal Distribution
#--------------------------------------

# First compute the estimated 90th percentile and
# associated 95% confidence interval using the complete data.
#------------------------------------------------------------
with(EPA.94b.tccb.df, eqlnorm(TcCB[Area == "Reference"], p = 0.9, ci = TRUE))

# Now compute the estimated 90th percentile and confidence interval
# using the censored data
#-------------------------------------------------------------------

# Estimate the mean and stanrdard deviation using the MLE method,
# and compute the confidence interval using the formula as if
#-----------------------------------------------------------------
with(Modified.TcCB.df, eqlnormCensored(TcCB, Censored, p = 0.9, ci = TRUE))

# Estimate the mean and stanrdard deviation using imputation
# with Q-Q regression (robust ROS), and compute the confidence
# interval using the gpq method.
#-----------------------------------------------------------------

with(Modified.TcCB.df, eqlnormCensored(TcCB, Censored, p = 0.9,
method = "impute.w.qq.reg", ci = TRUE,
ci.method = "gpq", seed = 47))

#--------------------
# Other Distributions
#--------------------

egamma.list <- with(Modified.TcCB.df, egammaCensored(TcCB, Censored))

eqgamma(egamma.list, p = 0.9)

# Clean up
#---------
rm(egamma.list)

##########################################################################

# 8.7 Prediction Intervals
#-------------------------

# 8.7.1 Parametric Prediction Intervals
#--------------------------------------

enorm.list <- with(EPA.09.Ex.18.3.TCE.df,
enormCensored(TCE.ppb[Well.type == "Background"],
Censored[Well.type == "Background"]))

predIntNorm(enorm.list, k = 4, pi.type = "upper", conf.level = 0.95, method = "exact")

# Clean up
#---------
rm(enorm.list)

##############################################################################

# 8.8 Tolerance Intervals
#------------------------

# 8.8.1 Parametric Tolerance Intervals
#-------------------------------------

# Estimate the mean and stanrdard deviation using the MLE method,
# and compute the tolerance interval using the formula as if
#-----------------------------------------------------------------
with(Modified.TcCB.df, tolIntLnormCensored(TcCB, Censored, ti.type = "upper"))

# Estimate the mean and stanrdard deviation using imputation
# with Q-Q regression (robust ROS), and compute the tolerence
# interval using the gpq method.
#-----------------------------------------------------------------
with(Modified.TcCB.df, tolIntLnormCensored(TcCB, Censored, ti.type = "upper",
method = "impute.w.qq.reg", ti.method = "gpq", seed = 47))

#############################################################################

# 8.9 Hypothesis Tests
#---------------------

# 8.9.1 Goodness-of-Fit Tests
#----------------------------

sw.list.norm <- with(Modified.TcCB.df, gofTestCensored(TcCB, Censored))
sw.list.norm

sw.list.lnorm <- with(Modified.TcCB.df, gofTestCensored(TcCB, Censored, dist = "lnorm"))
sw.list.lnorm

windows()
plot(sw.list.norm)

windows()
plot(sw.list.lnorm)

# Clean up
#---------
rm(sw.list.norm, sw.list.lnorm)

##############################################################################

# 8.9.2 Nonparametric Tests to Compare Two Groups
#------------------------------------------------

attach(Millard.Deverel.88.df)

Cu.AF <- Cu[Zone=="Alluvial.Fan"]
Cu.AF.cen <- Cu.censored[Zone=="Alluvial.Fan"]
Cu.BT <- Cu[Zone=="Basin.Trough"]
Cu.BT.cen <- Cu.censored[Zone=="Basin.Trough"]
twoSampleLinearRankTestCensored(Cu.AF, Cu.AF.cen, Cu.BT, Cu.BT.cen,
test = "normal.scores.2", var = "hypergeometric")

# Clean up
#---------
rm(Cu.AF, Cu.AF.cen, Cu.BT, Cu.BT.cen)

detach("Millard.Deverel.88.df")
```

## Try the EnvStats package in your browser

Any scripts or data that you put into this service are public.

EnvStats documentation built on Oct. 10, 2017, 1:05 a.m.