if (FALSE) {
library(agridat)
data(usgs.herbicides)
dat <- usgs.herbicides
# create censored data for one trait
dat$y <- as.numeric(dat$atrazine)
dat$ycen <- is.na(dat$y)
dat$y[is.na(dat$y)] <- .05
# library(NADA)
# # percent censored
# with(dat, censummary(y, censored=ycen))
# # median/mean
# with(dat, cenmle(y, ycen, dist="lognormal"))
# # boxplot
# with(dat, cenboxplot(obs=y, cen=ycen, log=FALSE,
# main="usgs.herbicides"))
# # with(dat, boxplot(y))
# pp <- with(dat, ros(obs=y, censored=ycen, forwardT="log")) # default lognormal
# plot(pp)
# plotfun <- function(vv){
# dat$y <- as.numeric(dat[[vv]])
# dat$ycen <- is.na(dat$y)
# dat$y[is.na(dat$y)] <- .01
# # qqnorm(log(dat$y), main=vv) # ordinary qq plot shows censored values
# pp <- with(dat, ros(obs=y, censored=ycen, forwardT="log"))
# plot(pp, main=vv) # omits censored values
# }
# op <- par(mfrow=c(3,3))
# vnames <- c("acetochlor", "alachlor", "ametryn", "atrazine","CIAT", "CEAT", "cyanazine", #"CAM",
# "dimethenamid", "flufenacet")
# for(vv in vnames) plotfun(vv)
# par(op)
}
Run the code above in your browser using DataLab