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

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

```{r}

library(webchem)
jagst

```


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)

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(CuZn$Zn,CuZn$ZnCen, 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)

```


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)

data(Example1) # NADA2 package
head(Example1)

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

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

# 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",ci.sample.size = sum(!censored))
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, n.bootstraps = 1000, pivot.statistic = "z",ci.sample.size = sum(!censored))

# 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",nmc = 1000, seed = NULL)

# 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)

# 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
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

```

Regrese s LOD

```{r}

library(lodr)

data(lod_data_ex)
head(lod_data_ex)

fitl <- lod_lm(data=lod_data_ex, frmla=y~x1+x2+x3, lod=c(0,0),var_LOD=c("x2", "x3"), nSamples=100, boots=0)
summary(fitl)
coef(fit)
residuals(fitl)

```


Hodnoty 17-alfa-hydroxyprogesteronu v krvi novorozencù

```{r}

library(openxlsx)
ohpx <- read.xlsx("c:\\Users\\lubop\\Documents\\kurs\\OHP.xlsx",
                  sheet = "OHP2", startRow = 1, colNames = TRUE, rowNames = FALSE,detectDates = FALSE, skipEmptyRows = TRUE,rows = NULL, cols = NULL,check.names = FALSE)
ohp2 <- rep(0, length(ohpx[,1]))
ohp2[which(ohpx[,1]=="<DL")]<- 1
ohp1 <- ohpx[,1]
ohp1[which(ohpx[,1]=="<DL")]<- 30
ohp1  <- as.numeric(ohp1)
ohp  <- cbind(ohp1,ohp2); ohp
colnames(ohp) = c("OHP","CENZ")
head(ohp)

as.logical(ohp[,"CENZ"])
pctCen(obs=ohp[,"OHP"], censored=as.logical(ohp[,"CENZ"]))
cenfit(Cen(obs=ohp[,"OHP"], censored=as.logical(ohp[,"CENZ"]))~1)


# podil hodnot pod DL (DL = 30)
t<-table(ohp[,"CENZ"])
names(t)<-c("necenzorováno","cenzorováno")
100*t/sum(t) # procenta cenzorovaných

# prumer vypocteny pouze z necenzorovanych hodnot
mean(ohp[ohp[,"CENZ"]==0,"OHP"]) 
median(ohp[ohp[,"CENZ"]==0,"OHP"]) 

# prumer vypocteny ze vsech hodnot pokud misto cenzorovani nechame detekcni limit
mean(ohp[,"OHP"])
median(ohp[,"OHP"]) 

# prumer vypocteny ze vsech hodnot pokud cenzorovane hodnoty nahradime nulou
xp<-ohp
xp[xp[,"CENZ"]==1,"OHP"]<-0
mean(xp[,"OHP"])
median(xp[,"OHP"]) 

# prumer vypocteny ze vsech hodnot pomoci survival analysis
library(survival)
# za predpokladu normalniho rozlozeni.
m<-survreg(formula=Surv(ohp[,"OHP"], ohp[,"CENZ"], type="left")~1, dist="gaussian")   
summary(m)
m[[1]]
# za predpokladu logaritmicko normalniho rozdeleni
m<- survreg(Surv(ohp[,"OHP"],ohp[,"CENZ"],type="left")~1,dist='loggaussian')
summary(m)
exp(m[[1]])

boxplot(ohp[,"OHP"])
hist(ohp[,"OHP"],freq = F)

```