---
title: "Analýza cenzorovaných dat "
output: html_notebook
---

Organic plant protection products in the river Jagst (Germany) in 2013

```{r}

library(webchem)
jagst

```


# https://rpubs.com/xuyzhang/870450



Data pod DL

```{r}

library(NADA)

v = c('<1', 1, '<1', 1, 2)
splitQual(v)

v = c('<DL', 1, '<DL', 1, 2)
splitQual(v)

v = c('1', 1, '1', 1, 2)
splitQual(v)


obs = c(0.5, 0.5, 1.0, 1.5, 5.0, 10, 100)
censored = c(TRUE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE)

pctCen(obs, censored)

Cen(obs, censored, type = "left")

```


## Trend

Celkové využitelné koncentrace chromu v proudech potoka Gales Creek v Oregonu, USA.

```{r}

library(NADA2)

data(Gales_Creek) # NADA2 package
head(Gales_Creek)

NADA2::ATS(Gales_Creek$TCr,Gales_Creek$CrND,Gales_Creek$dectime,LOG=FALSE)

NADA2::ATSmini(Gales_Creek$TCr,Gales_Creek$CrND,Gales_Creek$dectime)

```


Koncentrace mědi a zinku v podzemních vodách ze dvou zón v San Joaquin Valley v Kalifornii. 

```{r}

# boxplot
NADA2::cboxplot(CuZn$Zn,CuZn$ZnCen,CuZn$Zone,minmax = TRUE,Xlab="Zone",Ylab="Zn")

NADA2::cboxplot(CuZn$Zn,CuZn$ZnCen,CuZn$Zone,LOG=TRUE,Xlab="Zone",Ylab="Zn")

NADA2::cboxplot(CuZn$Zn,CuZn$ZnCen,CuZn$Zone,LOG=TRUE,show = TRUE, minmax = TRUE, Xlab="Zone",Ylab="Zn")

```


Koncentrace pyrenu v miligramech na litr z 20 stanic monitorujících kvalitu vody v PugetSound ve státě Washington, USA.

```{r}

library(NADA) # nelze NA
library(EnvStats)
library(survival)

data(ShePyrene) # NADA package
head(ShePyrene)

# QQ plot
NADA2::cenQQ(ShePyrene$Pyrene,ShePyrene$PyreneCen)

NADA2::cenCompareQQ(ShePyrene$Pyrene,ShePyrene$PyreneCen, Yname = "", printrslt = TRUE)

EnvStats::qqPlotCensored(ShePyrene$Pyrene,ShePyrene$PyreneCen,add.line = TRUE)

```


Koncentrace pyrenu v miligramech na litr z 20 stanic monitorujících kvalitu vody v PugetSound ve státě Washington, USA.

```{r}

library(NADA) # nelze NA
library(EnvStats)

data(ShePyrene) # NADA package
head(ShePyrene)
NADA::censummary(ShePyrene$Pyrene,ShePyrene$PyreneCen)

# vyber distribuce
EnvStats::distChooseCensored(CuZn$Zn,CuZn$ZnCen, censoring.side = "left", alpha = 0.05,method = "sf", choices = c("norm", "gamma", "lnorm"),est.arg.list = NULL, prob.method = "hirsch-stedinger",plot.pos.con = 0.375, warn = TRUE, keep.data = TRUE,data.name = NULL, censoring.name = NULL)

#
sw.norm <- EnvStats::gofTestCensored(ShePyrene$Pyrene,ShePyrene$PyreneCen, censoring.side = "left", test = "sf",distribution = "lnorm", est.arg.list = NULL,prob.method = "hirsch-stedinger", plot.pos.con = 0.375,keep.data = TRUE, data.name = NULL, censoring.name = NULL)
#dist = c("norm","lnorm","gamma")
sw.norm
plot(sw.norm)

# MLE
NADA::cenmle(ShePyrene$Pyrene,ShePyrene$PyreneCen)
EnvStats::elnormAltCensored(ShePyrene$Pyrene,ShePyrene$PyreneCen, ci=TRUE, ci.method ="bootstrap",n.bootstraps = 5000)

# Kaplan-Maier (KM)
NADA::cenfit(ShePyrene$Pyrene,ShePyrene$PyreneCen)
EnvStats::enparCensored(ShePyrene$Pyrene,ShePyrene$PyreneCen, ci=TRUE, ci.method ="bootstrap",n.bootstraps = 5000)

# Regression on Order Statistics (ROS)
ros = NADA::cenros(ShePyrene$Pyrene,ShePyrene$PyreneCen)
mean(ros)
sd(ros)
quantile(ros)
plot(ros)
EnvStats::elnormAltCensored(ShePyrene$Pyrene,ShePyrene$PyreneCen, method="rROS", ci=TRUE, ci.method ="bootstrap", n.bootstraps = 5000)

# nonpar
EnvStats::enparCensored(ShePyrene$Pyrene,ShePyrene$PyreneCen, censoring.side = "left", correct.se = FALSE,left.censored.min = "DL", right.censored.max = "DL", ci = FALSE,ci.method = "normal.approx", ci.type = "two-sided", conf.level = 0.95,pivot.statistic = "z", ci.sample.size = NULL, n.bootstraps = 1000)


# vsechny metody
NADA::censtats(ShePyrene$Pyrene,ShePyrene$PyreneCen)


d1 <- data.frame(time = ShePyrene$Pyrene, ind = as.numeric(ShePyrene$PyreneCen))
censdata1 = fitdistrplus::Surv2fitdistcens(time = d1$time, event = d1$ind, type = "left")
normd = fitdistcens(censdata, distr="norm", start=NULL, fix.arg=NULL,keepdata = TRUE, keepdata.nb=100)
lnormd = fitdistcens(censdata, distr="lnorm", start=NULL, fix.arg=NULL,keepdata = TRUE, keepdata.nb=100)
gammd = fitdistcens(censdata, distr="gamma", start=NULL, fix.arg=NULL,keepdata = TRUE, keepdata.nb=100)
ft = list(normd,lnormd,gammd)

fitdistrplus::cdfcompcens(ft, xlogscale = FALSE, ylogscale = FALSE,Turnbull.confint = FALSE,NPMLE.method = "Wang",add = FALSE, plotstyle = "graphics")

fitdistrplus::qqcompcens(ft, xlogscale = FALSE, ylogscale = FALSE, addlegend = TRUE, ynoise = TRUE, NPMLE.method = "Wang", plotstyle = "graphics")

fitdistrplus::ppcompcens(ft, xlogscale = FALSE, ylogscale = FALSE, addlegend = TRUE, ynoise = TRUE, NPMLE.method = "Wang", plotstyle = "graphics")


d1 <- data.frame(time = ShePyrene$Pyrene, ind = as.numeric(ShePyrene$PyreneCen))
survival::Surv(time = d1$time, event = d1$ind, type = "left")
mn <- survival::survreg(formula=Surv(d1$time, d1$ind, type="left")~1, dist="gaussian")   
summary(mn)
mn[[1]]
mln <- survival::survreg(Surv(d1$time,d1$ind,type="left")~1,dist='loggaussian')
summary(mln)
exp(mln[[1]])


```


Koncentrace mědi a zinku v podzemních vodách ze dvou zón v San Joaquin Valley v Kalifornii. 

```{r}

library(NADA) # nelze NA
library(EnvStats)

data(CuZn) # NADA package
head(CuZn)
NADA::censummary(CuZn$Zn,CuZn$ZnCen)

# MLE
NADA::cenmle(CuZn$Zn,CuZn$ZnCen)
EnvStats::elnormAltCensored(CuZn$Zn,CuZn$ZnCen, ci=TRUE, ci.method ="bootstrap", n.bootstraps = 5000)

# Kaplan-Maier (KM)
NADA::cenfit(CuZn$Zn,CuZn$ZnCen) # chyba
EnvStats::enparCensored(CuZn$Zn,CuZn$ZnCen, ci=TRUE, ci.method ="bootstrap", n.bootstraps = 5000)

# Regression on Order Statistics (ROS)
EnvStats::elnormAltCensored(CuZn$Zn,CuZn$ZnCen, method="rROS", ci=TRUE, ci.method ="bootstrap",n.bootstraps = 5000)

# nonpar
EnvStats::enparCensored(CuZn$Zn,CuZn$ZnCen, censoring.side = "left", correct.se = FALSE,left.censored.min = "DL", right.censored.max = "DL", ci = FALSE,ci.method = "normal.approx", ci.type = "two-sided", conf.level = 0.95,pivot.statistic = "z", ci.sample.size = NULL, n.bootstraps = 1000)

```

### Párový test

Test koncentrace arsenu v podzemní vodě oproti normě pro pitnou vodu 10 ug/l.

```{r}

library(NADA2)

data(Example1) # NADA2 package
head(Example1)

NADA2::cen_paired(Example1$Arsenic, Example1$NDisTRUE, 10, alt = "greater")

```

### Párový test

Testování, zda byly koncentrace atrazinu v podzemních vodách v červnu a v září ve stejných termínech v různých vrtech stejné

```{r}

library(NADA)
library(NADA2)

data(Atra) # NADA package
head(Atra)

NADA2::cen_paired(Atra$June, Atra$JuneCen, Atra$Sept, Atra$SeptCen)

NADA2::cen_signtest(Atra$June, Atra$JuneCen, Atra$Sept, Atra$SeptCen)

NADA2::cen_signedranktest(Atra$June, Atra$JuneCen, Atra$Sept, Atra$SeptCen) # silnejsi test

# Paired Prentice-Wilcoxon test 
NADA2::ppw.test(Atra$June, Atra$JuneCen, Atra$Sept, Atra$SeptCen, alternative = "two.sided", printstat = TRUE)

```


Koncentrace arsenu v podzemní vodě.

```{r}

library(NADA2)
library(EnvStats)
library(fitdistrplus)

data(Example1) # NADA2 package
head(Example1)

d1 = data.frame(time = Example1$Arsenic, ind = Example1$NDis1)
censdata = fitdistrplus::Surv2fitdistcens(time = d1$time, event = d1$ind, type = "left")
# type = c('right','left','interval','interval2')


NADA2::cboxplot(Example1$Arsenic, Example1$NDisTRUE, Ylab="Arsenic Conc", show = TRUE)

NADA2::cenCompareQQ(Example1$Arsenic, Example1$NDisTRUE, Yname = "", printrslt = TRUE)

fitdistrplus::plotdistcens(censdata, NPMLE = TRUE, NPMLE.method = "Wang")
fitdistrplus::plotdistcens(censdata, NPMLE = TRUE, NPMLE.method = "Turnbull.middlepoints")
fitdistrplus::plotdistcens(censdata, NPMLE = TRUE, NPMLE.method = "Turnbull.middlepoints", Turnbull.confint = TRUE)
fitdistrplus::plotdistcens(censdata,NPMLE = FALSE)


# gamma distribuce
EnvStats::egammaAltCensored(Example1$Arsenic, Example1$NDisTRUE, method = "mle", censoring.side = "left",ci = TRUE, ci.method = "profile.likelihood", ci.type = "two-sided",conf.level = 0.95, n.bootstraps = 1000, pivot.statistic = "z")
EnvStats::egammaCensored(Example1$Arsenic, Example1$NDisTRUE, method = "mle", censoring.side = "left",ci = TRUE, ci.method = "profile.likelihood", ci.type = "two-sided",conf.level = 0.95, pivot.statistic = "z", n.bootstraps = 1000)

gammd = fitdistrplus::fitdistcens(censdata, distr="gamma", start=NULL, fix.arg=NULL,keepdata = TRUE, keepdata.nb=100)
summary(gammd)
gammd$estimate
gammd$sd
fitdistrplus::quantile(gammd, probs = c(0.25, 0.5, 0.75))
fitdistrplus::mledist(censdata, distr="gamma", start = NULL, fix.arg = NULL, optim.method = "default",lower = -Inf, upper = Inf, custom.optim = NULL, weights = NULL, silent = TRUE,gradient = NULL, checkstartfix=FALSE)

fitdistrplus::plotdistcens(censdata, distr="gamma", para=list(gammd$estimate), leftNA = -Inf, rightNA = Inf,NPMLE = TRUE, Turnbull.confint = FALSE,NPMLE.method = "Wang")


# normalni distribuce
EnvStats::enormCensored(Example1$Arsenic, Example1$NDisTRUE, method = "mle", censoring.side = "left",ci = TRUE, ci.method = "profile.likelihood", ci.type = "two-sided",conf.level = 0.95, n.bootstraps = 1000, pivot.statistic = "z", seed = NULL)

normd = fitdistrplus::fitdistcens(censdata, distr="norm", start=NULL, fix.arg=NULL,keepdata = TRUE, keepdata.nb=100)
summary(normd)
normd$estimate
normd$sd
fitdistrplus::quantile(normd, probs = c(0.25, 0.5, 0.75))
fitdistrplus::mledist(censdata, distr="norm", start = NULL, fix.arg = NULL, optim.method = "default",lower = -Inf, upper = Inf, custom.optim = NULL, weights = NULL, silent = TRUE,gradient = NULL, checkstartfix=FALSE)

fitdistrplus::plotdistcens(censdata, distr="norm", para=list(normd$estimate), leftNA = -Inf, rightNA = Inf,NPMLE = TRUE, Turnbull.confint = FALSE,NPMLE.method = "Wang")


# log normalni distribuce
EnvStats::elnormAltCensored(Example1$Arsenic, Example1$NDisTRUE, method = "mle", censoring.side = "left",ci = TRUE, ci.method = "profile.likelihood", ci.type = "two-sided",conf.level = 0.95, n.bootstraps = 1000, pivot.statistic = "z")
EnvStats::elnormCensored(Example1$Arsenic, Example1$NDisTRUE, method = "mle", censoring.side = "left",ci = FALSE, ci.method = "profile.likelihood", ci.type = "two-sided",conf.level = 0.95, n.bootstraps = 1000, pivot.statistic = "z",nmc = 1000, seed = NULL)

lnormd = fitdistrplus::fitdistcens(censdata, distr="lnorm", start=NULL, fix.arg=NULL,keepdata = TRUE, keepdata.nb=100)
summary(lnormd)
lnormd$estimate
exp(lnormd$estimate)
lnormd$sd
fitdistrplus::quantile(lnormd, probs = c(0.25, 0.5, 0.75))
fitdistrplus::mledist(censdata, distr="lnorm", start = NULL, fix.arg = NULL, optim.method = "default",lower = -Inf, upper = Inf, custom.optim = NULL, weights = NULL, silent = TRUE,gradient = NULL, checkstartfix=FALSE)

fitdistrplus::plotdistcens(censdata, distr="lnorm", para=list(lnormd$estimate), leftNA = -Inf, rightNA = Inf,NPMLE = TRUE, Turnbull.confint = FALSE,NPMLE.method = "Wang")


# neparam. odhady
EnvStats::enparCensored(Example1$Arsenic, Example1$NDisTRUE, ci=TRUE, ci.method="bootstrap", ci.type="upper", n.bootstraps=5000)

# Box Cox
bcc = EnvStats::boxcoxCensored(Example1$Arsenic, Example1$NDisTRUE)
bcc
plot(bcc)
     
```

## Two-Group Tests

Koncentrace mědi a zinku v podzemních vodách ze dvou zón v San Joaquin Valley v Kalifornii.

```{r}

library(NADA)
library(EnvStats)
library(NADA2)

data(CuZn) # NADA package
head(CuZn)

NADA2::cboxplot(CuZn$Zn,CuZn$ZnCen,CuZn$Zone,Ylab="Zinc Concentration, in ug/L")

# MLE "t-test"
NADA2::cen2means(CuZn$Zn,CuZn$ZnCen,CuZn$Zone,LOG=FALSE, printstat = TRUE)

# perm test
NADA2::cenperm2(CuZn$Zn,CuZn$ZnCen,CuZn$Zone, R = 9999, alternative = "two.sided", printstat = TRUE)

# Wilcoxon rank sum test
EnvStats::twoSampleLinearRankTestCensored(CuZn$Zn[which(CuZn$Zone=="AlluvialFan")],CuZn$ZnCen[which(CuZn$Zone=="AlluvialFan")],CuZn$Zn[which(CuZn$Zone=="BasinTrough")],CuZn$ZnCen[which(CuZn$Zone=="BasinTrough")], test = "normal.scores.2", var = "hypergeometric")

```

## ANOVA

Koncentrace olova v krvi a některých orgánech volavek ve Virginii.

```{r}

library(NADA)
library(NADA2)

data(Golden) # NADA package
head(Golden)

NADA2::cboxplot(Golden$Liver, Golden$LiverCen, Golden$Dosage,Ylab = "Lead concentrations in liver, in ppb")

# MLE ANOVA
NADA2::cenanova(Golden$Liver, Golden$LiverCen, Golden$Dosage)

# perm ANOVA
NADA2::cenpermanova(Golden$Liver, Golden$LiverCen, Golden$Dosage)

# neparam. Peto - Peto test
NADA2::cen1way(Golden$Liver, Golden$LiverCen, Golden$Dosage)

```


## Korelace

Koncentrace olova v krvi a některých orgánech volavek ve Virginii.

```{r}

library(NADA)
library(NADA2)

data(PbHeron)
head(PbHeron)

# scatter plot

NADA2::kenplot(PbHeron$Blood, PbHeron$BloodCen, PbHeron$Kidney, PbHeron$KidneyCen,atsline = FALSE,xnam = NULL,ynam = NULL,Title = "")

NADA::cenxyplot(PbHeron$Blood, PbHeron$BloodCen, PbHeron$Kidney, PbHeron$KidneyCen, log="xy", lty="dashed")


# Kendall corel
NADA::cenken(PbHeron$Blood, PbHeron$BloodCen, PbHeron$Kidney, PbHeron$KidneyCen)


NADA2::kenplot(PbHeron$Blood, PbHeron$BloodCen, PbHeron$Kidney, PbHeron$KidneyCen, atsline = TRUE, xnam = "Cu",ynam = "Zn",Title = "Akritas - Theil - Sen line")

NADA2::ATS(PbHeron$Blood, PbHeron$BloodCen, PbHeron$Kidney, PbHeron$KidneyCen)

```


Koncentrace rozpuštěného železa v průběhu několika let v řece Brazos v Texasu. Byly použity letní koncentrace.

```{r}

library(NADA)
library(NADA2)

data(DFe) # NADA
head(DFe)

NADA::cenxyplot(DFe$Year, xcen=FALSE, DFe$Summer, DFe$SummerCen, log="", lty="dashed")

  
DFeReg = NADA::cenreg(Cen(DFe$Summer, DFe$SummerCen)~DFe$Year, dist="gaussian")
# dist = "lognormal", "gaussian"
DFeReg
summary(DFeReg)


DFeReg <-NADA2::cencorreg(DFe$Summer, DFe$SummerCen, DFe$Year)
DFeReg
summary(DFeReg)


DFeReg <-NADA2::ATS(DFe$Summer, DFe$SummerCen, DFe$Year,retrans=FALSE)
DFeReg
DFeReg

DFeReg <-NADA2::ATS(DFe$Summer, DFe$SummerCen, DFe$Year,retrans=TRUE)
DFeReg

DFeReg <-NADA2::ATSmini(DFe$Summer, DFe$SummerCen, DFe$Year)
DFeReg

```



Hodnoty 17-alfa-hydroxyprogesteronu v krvi novorozencù

```{r}

library(openxlsx)
ohpx <- read.xlsx("c:\\Users\\lubop\\Dropbox\\kurs\\OHP.xlsx",
                  sheet = "OHP2", startRow = 1, colNames = TRUE, rowNames = FALSE,detectDates = FALSE, skipEmptyRows = TRUE,rows = NULL, cols = NULL,check.names = FALSE)
head(ohpx)

# podil hodnot pod DL (DL = 30)
100*length(which(ohpx[,"OHP"]=="<30"))/length(ohpx[,"OHP"]) # procenta cenzorovaných

# prumer vypocteny pouze z necenzorovanych hodnot
nohp = as.numeric(ohpx[which(ohpx[,"OHP"]!="<30"),1])
mean(nohp) 
median(nohp) 

# prumer vypocteny ze vsech hodnot pokud misto cenzorovani nechame detekcni limit
lohp = ohpx[,"OHP"]
lohp[which(lohp=="<30")] = 30
lohp = as.numeric(lohp)
mean(lohp)
median(lohp) 

# prumer vypocteny ze vsech hodnot pokud cenzorovane hodnoty nahradime nulou
zohp = ohpx[,"OHP"]
zohp[which(zohp=="<30")] = 0
zohp = as.numeric(zohp)
mean(zohp)
median(zohp) 

# prumery vypoctene ze vsech hodnot pomoci survival analysis
library(NADA)
sohpx = splitQual(ohpx[,"OHP"])
library(survival)
# za predpokladu normalniho rozlozeni.
mg<-survreg(formula=Surv(sohpx$obs, sohpx$cen, type="left")~1, dist="gaussian")   
summary(mg)
mg[[1]]
# za predpokladu logaritmicko normalniho rozdeleni
mlg<- survreg(Surv(sohpx$obs, sohpx$cen,type="left")~1, dist='loggaussian')
summary(mlg)
exp(mlg[[1]])

```

