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

#Select dababase

t=file.choose()
nithanoi=read.csv(t,header=T,na.strings = ".")
nithanoi=nithanoi%>%
  mutate(Date=dmy(Date))%>%
  mutate(weekday= wday(Date,label=T,abbr = T))%>%
  mutate(Period=as.factor(Period))

#Biomarkers in two period

createTable(compareGroups(Period~NIT1+COT1+HCOT1+NIT2+COT2+HCOT2+NITavg+COTavg+HCOTavg, data = nithanoi))
## 
## --------Summary descriptives table by 'Period'---------
## 
## _________________________________________ 
##              1           2      p.overall 
##            N=30        N=31               
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
## NIT1    5.23 (1.17) 5.49 (1.16)   0.382   
## COT1    2.13 (0.26) 2.34 (0.53)   0.059   
## HCOT1   2.88 (0.35) 2.30 (0.51)  <0.001   
## NIT2    4.07 (1.01) 5.07 (1.03)  <0.001   
## COT2    1.85 (0.54) 2.47 (0.47)  <0.001   
## HCOT2   2.54 (0.72) 2.46 (0.46)   0.629   
## NITavg  4.63 (0.98) 5.28 (0.92)   0.010   
## COTavg  1.99 (0.37) 2.41 (0.46)  <0.001   
## HCOTavg 2.70 (0.49) 2.39 (0.44)   0.010   
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

#Biomarkers detected in 2 sites

attach(nithanoi)
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=61        N=60               
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
## AmountNIT 5.36 (1.17) 4.57 (1.13)  <0.001   
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
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=61        N=60               
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
## AmountCOT 2.24 (0.43) 2.16 (0.59)   0.384   
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
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=61        N=60               
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
## AmountHCOT 2.59 (0.52) 2.50 (0.60)   0.399   
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

#Weekly pattern of nicotine consumption

shapiro.test(nithanoi$NITuse)
## 
##  Shapiro-Wilk normality test
## 
## data:  nithanoi$NITuse
## W = 0.96428, p-value = 0.07216
ad.test(nithanoi$NITuse)
## 
##  Anderson-Darling normality test
## 
## data:  nithanoi$NITuse
## A = 0.78477, p-value = 0.03946
qqnorm(nithanoi$NITuse)
qqline(nithanoi$NITuse,col="red")

kruskal.test(NITuse~weekday, data=nithanoi)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  NITuse by weekday
## Kruskal-Wallis chi-squared = 4.6762, df = 6, p-value = 0.586
av=aov(nithanoi$NITuse~nithanoi$weekday)
summary(av)
##                  Df    Sum Sq   Mean Sq F value Pr(>F)
## nithanoi$weekday  6 2.360e-07 3.932e-08   0.573   0.75
## Residuals        54 3.702e-06 6.856e-08
p=ggplot(nithanoi, aes(x=weekday, y=NITuse, 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_title()+ 
  easy_center_title()+
  theme(legend.position = "bottom")

#Compare nicotine consumption in 2 periods

nithanoi%>%ggplot(aes(x= Period,y=NITuse, fill=Period))+
  geom_boxplot()+ theme_bw()+ theme_classic()+ labs(y="Nicotine consumption (mg/person/day)",x=NULL) + 
  ggtitle("Estimation of nicotine consumption over 2 periods")+ 
  easy_center_title()+ easy_remove_legend_title()+
  theme(legend.position = "bottom")

createTable(compareGroups(Period~NITuse,data = nithanoi))
## 
## --------Summary descriptives table by 'Period'---------
## 
## ________________________________________ 
##             1           2      p.overall 
##           N=30        N=31               
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
## NITuse 0.00 (0.00) 0.00 (0.00)  <0.001   
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
table1(~NITuse|Period, data=nithanoi)
1
(N=30)		 2
(N=31)		 Overall
(N=61)
NITuse
Mean (SD) 0.00109 (0.000204) 0.00132 (0.000253) 0.00121 (0.000256)
Median [Min, Max] 0.00118 [0.000563, 0.00131] 0.00134 [0.000636, 0.00168] 0.00123 [0.000563, 0.00168]