#Loading packages

require(nortest)
## Loading required package: nortest
require(tidyverse)
## Loading required package: tidyverse
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.3     v purrr   0.3.4
## v tibble  3.1.1     v dplyr   1.0.5
## v tidyr   1.1.3     v stringr 1.4.0
## v readr   1.4.0     v forcats 0.5.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
require(lubridate)
## Loading required package: lubridate
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
require(coin)
## Loading required package: coin
## Loading required package: survival
require(ggeasy)
## Loading required package: ggeasy
require(compareGroups)
## Loading required package: compareGroups
require(table1)
## Loading required package: table1
## 
## Attaching package: 'table1'
## The following objects are masked from 'package:base':
## 
##     units, units<-
require(pgirmess)
## Loading required package: pgirmess

#Select dababase

setwd("C:/Users/Dell/OneDrive - UMP/Research/Queensland team/WBE  Hanoi")
nithanoi=read.csv("nitinhanoi.csv",header=T,na.strings = ".")
nithanoi=nithanoi%>%
  mutate(Date=dmy(Date))%>%
  mutate(weekday= wday(Date,label=T,abbr = T))%>%
  mutate(Period=as.factor(Period))
attach(nithanoi)

#Compare amount of biomarkers detected

biomarkers=data.frame(NITavg,COTavg,HCOTavg)%>%
  gather("Biomarker", "Amount",c(NITavg,COTavg,HCOTavg))%>%
  mutate(Biomarker=recode(Biomarker,"NITavg"="Nicotine", "COTavg"="Cotinine","HCOTavg"="3-hydroxycotinine"))
ad.test(biomarkers$Amount)
## 
##  Anderson-Darling normality test
## 
## data:  biomarkers$Amount
## A = 10.607, p-value < 2.2e-16
kruskal.test(Amount~as.factor(Biomarker),data = biomarkers)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Amount by as.factor(Biomarker)
## Kruskal-Wallis chi-squared = 185.7, df = 2, p-value < 2.2e-16
kruskalmc(biomarkers$Amount,as.factor(biomarkers$Biomarker))
## Multiple comparison test after Kruskal-Wallis 
## p.value: 0.05 
## Comparisons
##                              obs.dif critical.dif difference
## 3-hydroxycotinine-Cotinine  47.93478      28.1738       TRUE
## 3-hydroxycotinine-Nicotine 108.57065      28.1738       TRUE
## Cotinine-Nicotine          156.50543      28.1738       TRUE
table1(~Amount|as.factor(Biomarker),data = biomarkers)
3-hydroxycotinine
(N=92)		 Cotinine
(N=92)		 Nicotine
(N=92)		 Overall
(N=276)
Amount
Mean (SD) 2.73 (0.542) 2.25 (0.434) 5.10 (1.01) 3.36 (1.43)
Median [Min, Max] 2.84 [1.25, 3.76] 2.30 [1.04, 3.06] 5.08 [1.42, 7.98] 2.86 [1.04, 7.98] 
createTable(compareGroups(Biomarker~Amount, data = biomarkers))
## 
## --------Summary descriptives table by 'Biomarker'---------
## 
## __________________________________________________________ 
##        3-hydroxycotinine  Cotinine    Nicotine   p.overall 
##              N=92           N=92        N=92               
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
## Amount    2.73 (0.54)    2.25 (0.43) 5.10 (1.01)  <0.001   
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
p=ggplot(biomarkers, aes(x=Biomarker, y=Amount, fill=Biomarker))
p+ geom_boxplot() + ggtitle("Biomarkers detected in wastewater")+ labs(y="Amount (µg/ml)",x=NULL)+
  #coord_flip()+ 
  theme_bw()+ 
  theme_classic()+
  easy_remove_legend_title()+
  easy_center_title()+
  theme(legend.position = "bottom")

#Weekly mass loads of three biomarkers

kruskal.test(NITavg~weekday,data=nithanoi)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  NITavg by weekday
## Kruskal-Wallis chi-squared = 3.5726, df = 6, p-value = 0.7343
kruskal.test(COTavg~weekday,data=nithanoi)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  COTavg by weekday
## Kruskal-Wallis chi-squared = 3.3779, df = 6, p-value = 0.7601
kruskal.test(HCOTavg~weekday,data = nithanoi)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  HCOTavg by weekday
## Kruskal-Wallis chi-squared = 2.61, df = 6, p-value = 0.856

#Biomarkers in three period

table1(~NIT1+NIT2+COT1+COT2+HCOT1+HCOT2|Period,data = biomarkers)
1
(N=30)		 2
(N=31)		 3
(N=31)		 Overall
(N=92)
NIT1
Mean (SD) 5.23 (1.17) 5.49 (1.16) 5.79 (1.58) 5.51 (1.33)
Median [Min, Max] 5.25 [1.86, 7.94] 5.37 [3.23, 8.77] 5.56 [2.92, 9.90] 5.36 [1.86, 9.90]
NIT2
Mean (SD) 4.07 (1.01) 5.07 (1.03) 4.88 (1.15) 4.68 (1.14)
Median [Min, Max] 4.31 [0.970, 5.45] 4.93 [3.51, 7.42] 4.99 [1.71, 6.98] 4.70 [0.970, 7.42]
Missing 0 (0%) 1 (3.2%) 0 (0%) 1 (1.1%)
COT1
Mean (SD) 2.14 (0.261) 2.34 (0.525) 2.40 (0.465) 2.29 (0.444)
Median [Min, Max] 2.13 [1.34, 2.54] 2.41 [0.970, 3.27] 2.40 [1.44, 3.83] 2.29 [0.970, 3.83]
COT2
Mean (SD) 1.85 (0.539) 2.47 (0.470) 2.32 (0.455) 2.21 (0.550)
Median [Min, Max] 1.99 [0.730, 2.45] 2.58 [1.10, 3.40] 2.41 [0.910, 2.87] 2.32 [0.730, 3.40]
Missing 0 (0%) 1 (3.2%) 0 (0%) 1 (1.1%)
HCOT1
Mean (SD) 2.88 (0.352) 2.30 (0.505) 3.14 (0.598) 2.77 (0.607)
Median [Min, Max] 2.93 [1.89, 3.52] 2.31 [1.09, 3.04] 3.16 [1.94, 5.07] 2.80 [1.09, 5.07]
HCOT2
Mean (SD) 2.54 (0.722) 2.46 (0.465) 3.06 (0.622) 2.69 (0.663)
Median [Min, Max] 2.80 [1.17, 3.57] 2.47 [1.13, 3.26] 3.19 [1.23, 4.19] 2.80 [1.13, 4.19]
Missing 0 (0%) 1 (3.2%) 0 (0%) 1 (1.1%)

#Biomarkers detected in 2 sites

#Nicotin
nit=data.frame(Period,NIT1,NIT2)%>%
  gather("Location", "AmountNIT",c(-Period))
createTable(compareGroups(Location~AmountNIT,data=nit))
## 
## --------Summary descriptives table by 'Location'---------
## 
## ___________________________________________ 
##              NIT1        NIT2     p.overall 
##              N=92        N=91               
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
## AmountNIT 5.51 (1.33) 4.68 (1.14)  <0.001   
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
oneway_test(nit$AmountNIT~as.factor(nit$Location),distribution=approximate(nresample = 10000))
## 
##  Approximative Two-Sample Fisher-Pitman Permutation Test
## 
## data:  nit$AmountNIT by as.factor(nit$Location) (NIT1, NIT2)
## Z = 4.3138, p-value = 1e-04
## alternative hypothesis: true mu is not equal to 0
#Cotinine
cot=data.frame(Period,COT1,COT2)%>%
  gather("Location", "AmountCOT",c(-Period))
createTable(compareGroups(Location~AmountCOT,data=cot))
## 
## --------Summary descriptives table by 'Location'---------
## 
## ___________________________________________ 
##              COT1        COT2     p.overall 
##              N=92        N=91               
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
## AmountCOT 2.29 (0.44) 2.21 (0.55)   0.267   
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
oneway_test(cot$AmountCOT~as.factor(cot$Location),distribution=approximate(nresample = 10000))
## 
##  Approximative Two-Sample Fisher-Pitman Permutation Test
## 
## data:  cot$AmountCOT by as.factor(cot$Location) (COT1, COT2)
## Z = 1.1132, p-value = 0.2687
## alternative hypothesis: true mu is not equal to 0
#Hydroxycotinine
hcot=data.frame(Period,HCOT1,HCOT2)%>%
  gather("Location", "AmountHCOT",c(-Period))
createTable(compareGroups(Location~AmountHCOT,data=hcot))
## 
## --------Summary descriptives table by 'Location'---------
## 
## ____________________________________________ 
##               HCOT1       HCOT2    p.overall 
##               N=92        N=91               
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
## AmountHCOT 2.77 (0.61) 2.69 (0.66)   0.384   
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
oneway_test(hcot$AmountHCOT~as.factor(hcot$Location),distribution=approximate(nresample = 10000))
## 
##  Approximative Two-Sample Fisher-Pitman Permutation Test
## 
## data:  hcot$AmountHCOT by as.factor(hcot$Location) (HCOT1, HCOT2)
## Z = 0.87349, p-value = 0.3769
## alternative hypothesis: true mu is not equal to 0

#Pattern of nicotine consumption

shapiro.test(nithanoi$NICabs)
## 
##  Shapiro-Wilk normality test
## 
## data:  nithanoi$NICabs
## W = 0.95754, p-value = 0.004442
ad.test(nithanoi$NICabs)
## 
##  Anderson-Darling normality test
## 
## data:  nithanoi$NICabs
## A = 1.0924, p-value = 0.006925
qqnorm(nithanoi$NICabs)
qqline(nithanoi$NICabs,col="red")

#Overall use of nicotine
table1(~NICabs+ Tobacco,data=nithanoi)
Overall
(N=92)
NICabs
Mean (SD) 1.24 (0.239)
Median [Min, Max] 1.26 [0.569, 1.68]
Tobacco
Mean (SD) 9.62 (1.85)
Median [Min, Max] 9.79 [4.42, 13.1] 
#Weekly patterns of nicotine consumption
kruskal.test(NICabs~weekday, data=nithanoi)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  NICabs by weekday
## Kruskal-Wallis chi-squared = 3.3779, df = 6, p-value = 0.7601
av=aov(nithanoi$NICabs~nithanoi$weekday)
summary(av)
##                  Df Sum Sq Mean Sq F value Pr(>F)
## nithanoi$weekday  6  0.189 0.03158   0.538  0.778
## Residuals        85  4.992 0.05873
nithanoi$weekday=factor(nithanoi$weekday, levels=c("Mon","Tue","Wed","Thu","Fri","Sat","Sun"))
p=ggplot(nithanoi, aes(x=weekday, y=NICabs, fill=weekday))
p+ geom_boxplot() + ggtitle("Weekly profile of nicotine consumption in Hanoi")+ labs(y="Nicotine consumption (mg/person/day)",x=NULL)+
  coord_flip()+ 
  theme_bw()+ 
  theme_classic()+
  easy_remove_legend()+
  easy_center_title()

#Compare nicotine consumption in 3 periods

nithanoi%>%
  mutate(Period=recode(Period, "1"="Sep 2018","2"="Dec 2018-Jan 2019", "3"= "Dec 2019- Jan 2020"))%>%
  ggplot(aes(x= Period,y=NICabs, fill=Period))+
  geom_boxplot()+ theme_bw()+ theme_classic()+ labs(y="Nicotine consumption (mg/person/day)",x=NULL) + 
  ggtitle("Estimation of nicotine consumption over 3 periods")+ 
  easy_remove_legend()+
  easy_center_title()

createTable(compareGroups(Period~NICabs,data = nithanoi))
## 
## --------Summary descriptives table by 'Period'---------
## 
## ____________________________________________________ 
##             1           2           3      p.overall 
##           N=30        N=31        N=31               
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
## NICabs 1.09 (0.20) 1.32 (0.25) 1.30 (0.20)  <0.001   
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
table1(~NICabs|Period, data=nithanoi)
1
(N=30)		 2
(N=31)		 3
(N=31)		 Overall
(N=92)
NICabs
Mean (SD) 1.09 (0.201) 1.32 (0.253) 1.30 (0.196) 1.24 (0.239)
Median [Min, Max] 1.18 [0.569, 1.31] 1.33 [0.637, 1.68] 1.33 [0.646, 1.57] 1.26 [0.569, 1.68] 
kruskal.test(NICabs~as.factor(Period),data=nithanoi)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  NICabs by as.factor(Period)
## Kruskal-Wallis chi-squared = 19.834, df = 2, p-value = 4.933e-05
kruskalmc(NICabs,as.factor(Period))
## Multiple comparison test after Kruskal-Wallis 
## p.value: 0.05 
## Comparisons
##       obs.dif critical.dif difference
## 1-2 27.701613     16.37151       TRUE
## 1-3 24.975806     16.37151       TRUE
## 2-3  2.725806     16.23676      FALSE