# updated 8/24/10
# updated 3/23/2010

# This S-PLUS script computes left-censored summary statistics for NAWQA Trace Element data
# It also computes summary statistics and quantiles for Barium and Strontium which do not have 
# censored data

# Text output will be written to the file specified in the sink command.
sink("report.trace.summary.out.txt")



# Read data file
# Nondetects are recensored to LT-MDLs
# Data have been revised to remove nondetects above 1 ug/L
attach.data.frame(trace.nawqa.censored)


#start of wilson function that we create as 
wilson<-function(target, target2, target3) {


cat("Processing target", target3, "\n")

#Parametric Summary Statistics for left censored data
print(mdlpar(lcens(target, target2 == "<")))

cat("Processing target", target3, "\n")

#Parametric Summary Statistics and Percentiles for left censored data
print(summary(mdlpar(lcens(target, target2 == "<"))))



# End of wilson function
}

# for every trace element generate a new line below
# following the example

wilson(aluminum, aluminum.rmk, "aluminum")
wilson(antimony, antimony.rmk, "antimony")
wilson(arsenic, arsenic.rmk, "arsenic")
wilson(beryllium, beryllium.rmk, "beryllium")
# wilson(boron, boron.rmk, "boron")
wilson(cadmium, cadmium.rmk, "cadmium")
wilson(chromium, chromium.rmk, "chromium")
wilson(cobalt, cobalt.rmk, "cobalt")
wilson(copper, copper.rmk, "copper")
wilson(iron, iron.rmk, "iron")
wilson(lead, lead.rmk, "lead")
# wilson(lithium, lithium.rmk, "lithium")
wilson(manganese, manganese.rmk, "manganese")
wilson(molybdenum, molybdenum.rmk, "molybdenum")
wilson(nickel, nickel.rmk, "nickel")
wilson(selenium, selenium.rmk, "selenium")
wilson(silver, silver.rmk, "silver")
wilson(thallium, thallium.rmk, "thallium")
wilson(uranium, uranium.rmk, "uranium")
wilson(vanadium, vanadium.rmk, "vanadium")
wilson(zinc, zinc.rmk, "zinc")


#start of wilson2 function that we create as 
wilson2<-function(target, target5) {


cat("Processing target", target5, "\n")

# do separate summary statistics for analytes with no censored data
# **Lithium only has one nondetect**
print(summary(target))

cat("Processing target", target5, "\n")

print(quantile(target,c(0.10,0.90), na.rm=T))



# End of wilson2 function
}



wilson2(barium, "barium")
wilson2(boron, "boron")
wilson2(lithium, "lithium")
wilson2(strontium, "strontium")



# do separate summary statistics for analytes with no censored data
# print(summary(barium))
# print(quantile(barium,c(0.10,0.90), na.rm=T))

# print(summary(strontium))
# print(quantile(strontium,c(0.10,0.90), na.rm=T))

#Detach the data frame
detach(trace.nawqa.censored)

# close output file
sink()