######################################
# Basic requirements                 # 
######################################

df <- subset(Data,complete=="1")

library(readxl)
library(ggplot2)
library(tidyverse)
library(bfsMaps)

######################################
# Table 1                            # 
######################################


table(df$lang)
round(100*table(df$lang)/259,1)

table(df$q01)
round(100*table(df$q01)/259,1)

table(df$q02)
round(100*table(df$q02)/259,1)

table(df$q03)
round(100*table(df$q03)/259,1)

table(df$q04)
round(100*table(df$q04)/259,1)


table(df$q05)
round(100*table(df$q05)/259,1)

table(df$q07b)
round(100*table(df$q07b)/259,1)

df$q07_clinpharm <- 0
df$q07_clinpharm[df$q07c==1 | df$q07d==1 | df$q07e==1] <- 1
table(df$q07_clinpharm)
round(100*table(df$q07_clinpharm)/259,1)

table(df$q07f)
round(100*table(df$q07f)/259,1)

table(df$q07g)
round(100*table(df$q07g)/259,1)

table(df$q07h)
round(100*table(df$q07h)/259,1)

table(df$q07i)
round(100*table(df$q07i)/259,1)

table(df$q07j)
round(100*table(df$q07j)/259,1)


table(df$q08e)
round(100*table(df$q08e)/259,1)

table(df$q09)
round(100*table(df$q09)/259,1)

table(df$q10)
round(100*table(df$q10)/259,1)

table(df$q12)
round(100*table(df$q12)/259,1)


df$q30[df$q30==0] <- NA
table(df$q30)


df$q21[df$q21==0] <- NA
table(df$q21)
sum(table(df$q21))



df$q32[df$q32==0] <- NA
table(df$q32)
sum(table(df$q32))


df$q34[df$q34==0 | df$q34 > 4] <- NA
table(df$q34)
sum(table(df$q34))


######################################
# Figure 2                           # 
######################################


x <- c(15,42,10,0,0,1,1,0,2,10,5,3,4,0,0,0,7,4,35,9,14,21,7,5,3,2)/c(239,178,35,1,16,2,3,3,16,83,26,75,49,14,6,1,54,46,127,25,205,254,124,57,180,20)

names(x) <- levels(d.bfsrg$kt_x)

round(100*x,1)

# define the a color ramp with 10 colors
cols <- colorRampPalette(colors = c("white","lightyellow2","darkolivegreen1","green3","springgreen4"))(500)
PlotKant(rownames(x), col=FindColor(x, cols = cols, min.x=0, max.x=.5),
         border="black",main="\nParticipation by cantons")
ColorLegend(x="left", width=15000, labels=paste0(seq(0, 50, 5),"%"),  
            cols=cols, cex=0.8, adj=c(1,0.5), frame="black", inset=c(-0.09, 0.35))
AddLakes(col="lightsteelblue1", border="black" )
AddRivers(col="steelblue1",lwd=1.25)
PlotCH(col=NA, add=TRUE, lwd=1.25, border="black")
points(sf::st_coordinates(GetMap("stkt.pnt")$geometry),
       pch=10, bg="black",cex=1.2, lwd=1.35)


######################################
# Figure 3                           # 
######################################

df$q13a[df$q13a==0] <- NA
df$q13b[df$q13b==0] <- NA
df$q13c[df$q13c==0] <- NA

table(df$q13a)
sum(table(df$q13a))
table(df$q13b)
sum(table(df$q13b))
table(df$q13c)
sum(table(df$q13c))

y <- c(rep("Brief counseling\nduring NRT sale\n(n = 253)", 5) , rep("Opportunistic tobacco\ncessation counseling\n(n = 250)" , 5) , rep("Dedicated smoking\ncessation counseling\n(n = 236)" , 5))
y <- factor(y, levels=c("Brief counseling\nduring NRT sale\n(n = 253)","Opportunistic tobacco\ncessation counseling\n(n = 250)","Dedicated smoking\ncessation counseling\n(n = 236)"))
x <- rep(c("Daily" , "Weekly","Monthly","Yearly", "Never"), 3)
x <- factor(x, c("Daily" ,"Weekly", "Monthly","Yearly", "Never"))
value <- c(table(df$q13a),table(df$q13b),c(0,4,45,72,115))
data <- data.frame(x,y,value)



data <- data %>% arrange(x)
data <- data %>% arrange(y)

data$x <- factor(x, c("Never", "Yearly","Monthly","Weekly","Daily","Always"))

data$group <- c(rep(1,5),rep(2,5),rep(3,5))
data$pos[data$group==1] <- cumsum(data$value[data$group==1]) - (data$value[data$group==1] / 2)
data$tot[data$group==1] <- sum(data$value[data$group==1])
data$prop[data$group==1] <- paste0(data$val[data$group==1]," (",round(100* (data$val[data$group==1] / data$tot[data$group==1]),1),"%)")
data$pos[data$group==2] <- cumsum(data$value[data$group==2]) - (data$value[data$group==2] / 2)
data$tot[data$group==2] <- sum(data$value[data$group==2])
data$prop[data$group==2] <- paste0(data$val[data$group==2]," (",round(100* (data$val[data$group==2] / data$tot[data$group==2]),1),"%)")
data$pos[data$group==3] <- cumsum(data$value[data$group==3]) - (data$value[data$group==3] / 2)
data$tot[data$group==3] <- sum(data$value[data$group==3])
data$prop[data$group==3] <- paste0(data$val[data$group==3]," (",round(100* (data$val[data$group==3] / data$tot[data$group==3]),1),"%)")

data$prop[9] <- paste0(data$val[9]," (","34.0%)" )

data$prop[data$value / data$tot < .025] <- NA


ggplot(data, aes(fill=x, y=value, x=y)) + 
  geom_bar(position="fill", stat="identity")+scale_y_continuous(labels = scales::percent, breaks = c(0,.1,.2,.3,.4,.5,.6,.7,.8,.9,1))+
  scale_fill_manual(values=c(
    "#FFDB7D",
    "#fffdba",
    "#d3e7cb",
    "#a8d7a5",
    "#AeD5F1"),name="") +
  labs(x = "",y="")+
  geom_text(aes(y= pos/tot, label=prop), vjust=1, color="black", size=3)+
  theme(legend.position="right")+ guides(fill = guide_legend(nrow = 5, byrow = TRUE))



table(df$q14a,df$q08e)
table(df$q14b,df$q08e)
table(df$q14c,df$q08e)




df$q24a[df$q24a==0] <- NA
df$q24b[df$q24b==0] <- NA
df$q24c[df$q24c==0] <- NA

table(df$q24a)
sum(table(df$q24a[df$q08e==0]))
table(df$q24b)
sum(table(df$q24b[df$q08e==0]))
table(df$q24c)
sum(table(df$q24c[df$q08e==0]))



plot2<-data.frame(Mean=c(mean(df$q24a[df$q08e==0],na.rm=T),mean(df$q24b[df$q08e==0],na.rm=T),mean(df$q24c[df$q08e==0],na.rm=T)), 
                  sd=c(sd(df$q24a[df$q08e==0],na.rm=T),sd(df$q24b[df$q08e==0],na.rm=T),sd(df$q24c[df$q08e==0],na.rm=T)), 
                  Category=c("Smoking cessation counseling\nas part of my activities","Motivation to conduct smoking\ncessation counseling","Perception of being\nsufficiently trained"))
plot2$Category <- factor(plot2$Category, levels=c("Perception of being\nsufficiently trained","Motivation to conduct smoking\ncessation counseling","Smoking cessation counseling\nas part of my activities"))


ggplot(plot2,  aes(x=Mean, y=Category, fill=Category)) + scale_x_continuous(breaks=(1:6*1),limits=c(0.995,6.005))+  
  geom_point(col="red",cex=.7)+ 
  geom_errorbar(aes(xmin=Mean-(1.96*(sd/sqrt(215))), xmax=Mean+(1.96*(sd/sqrt(215)))),col="red" , width=.2, 
                position=position_dodge(0.05)) + 
  #  facet_grid(~Category)+ 
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank(), panel.grid.minor = element_blank(),
        axis.text.x=element_text(angle = 0, hjust = 1),
        legend.position = "none")




df$q25a[df$q25a==0] <- NA
df$q25b[df$q25b==0] <- NA
df$q25c[df$q25c==0] <- NA
df$q25d[df$q25d==0] <- NA
df$q25e[df$q25e==0 | df$q25e>6] <- NA

table(df$q25a)
sum(table(df$q25a))
table(df$q25b)
sum(table(df$q25b))
table(df$q25c)
sum(table(df$q25c))
table(df$q25d)
sum(table(df$q25d))
table(df$q25e)
sum(table(df$q25e))



plot2<-data.frame(Mean=c(mean(df$q25a,na.rm=T),mean(df$q25b,na.rm=T),mean(df$q25c,na.rm=T),mean(df$q25d,na.rm=T),mean(df$q25e,na.rm=T)), 
                  sd=c(sd(df$q25a,na.rm=T),sd(df$q25b,na.rm=T),sd(df$q25c,na.rm=T),sd(df$q25d,na.rm=T),sd(df$q25e,na.rm=T)), 
                  Category=c("Financial compensation","More training","A decision aid","Higher demand","Better collaboration"))
plot2$Category <- factor(plot2$Category, levels=c("Better collaboration","Higher demand","A decision aid","More training","Financial compensation"))

n <- c(250,254,252,250,229)
ggplot(plot2,  aes(x=Mean, y=Category, fill=Category)) + scale_x_continuous(breaks=(1:6*1),limits=c(0.995,6.005))+  
  geom_point(col="red",cex=.7)+ 
  geom_errorbar(aes(xmin=Mean-(1.96*(sd/sqrt(n))), xmax=Mean+(1.96*(sd/sqrt(n)))),col="red" , width=.2, 
                position=position_dodge(0.05)) + 
  #  facet_grid(~Category)+ 
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank(), panel.grid.minor = element_blank(),
        axis.text.x=element_text(angle = 0, hjust = 1),
        legend.position = "none")




######################################
# Figure 4                           # 
######################################


df$q36a[df$q36a==0] <- NA
df$q36b[df$q36b==0] <- NA
df$q36c[df$q36c==0] <- NA
df$q36d[df$q36d==0] <- NA
df$q36e[df$q36e==0] <- NA
df$q36f[df$q36f==0] <- NA

df$q37a[df$q37a==0] <- NA
df$q37b[df$q37b==0] <- NA
df$q37c[df$q37c==0] <- NA
df$q37d[df$q37d==0] <- NA
df$q37e[df$q37e==0] <- NA
df$q37f[df$q37f==0] <- NA


df$q38a[df$q38a==0] <- NA
df$q38b[df$q38b==0] <- NA
df$q38c[df$q38c==0] <- NA
df$q38d[df$q38d==0] <- NA
df$q38e[df$q38e==0] <- NA
df$q38e[df$q38f==0] <- NA
df$q38e[df$q38g==0] <- NA
df$q38e[df$q38h==0] <- NA
df$q38e[df$q38i==0] <- NA
df$q38e[df$q38j==0] <- NA
df$q38e[df$q38k==0] <- NA


plot1<-data.frame(Mean=c(mean(df$q36a,na.rm=T),mean(df$q36b,na.rm=T),mean(df$q36c,na.rm=T),mean(df$q36d,na.rm=T),mean(df$q36e,na.rm=T),mean(df$q36f,na.rm=T),
                         mean(df$q37a,na.rm=T),mean(df$q37b,na.rm=T),mean(df$q37c,na.rm=T),mean(df$q37d,na.rm=T),mean(df$q37e,na.rm=T),mean(df$q37f,na.rm=T),
                         mean(df$q38a,na.rm=T),mean(df$q38b,na.rm=T),mean(df$q38c,na.rm=T),mean(df$q38d,na.rm=T),mean(df$q38e,na.rm=T),mean(mean(df$q38f,na.rm=T),mean(df$q38g,na.rm=T),mean(df$q38h,na.rm=T),mean(df$q38i,na.rm=T),mean(df$q38k,na.rm=T))), 
                  sd=c(sd(df$q36a,na.rm=T),sd(df$q36b,na.rm=T),sd(df$q36c,na.rm=T),sd(df$q36d,na.rm=T),sd(df$q36e,na.rm=T),sd(df$q36f,na.rm=T),
                       sd(df$q37a,na.rm=T),sd(df$q37b,na.rm=T),sd(df$q37c,na.rm=T),sd(df$q37d,na.rm=T),sd(df$q37e,na.rm=T),sd(df$q37f,na.rm=T),
                       sd(df$q38a,na.rm=T),sd(df$q38b,na.rm=T),sd(df$q38c,na.rm=T),sd(df$q38d,na.rm=T),sd(df$q38e,na.rm=T),mean(sd(df$q38f,na.rm=T),sd(df$q38g,na.rm=T),sd(df$q38h,na.rm=T),sd(df$q38i,na.rm=T),sd(df$q38k,na.rm=T))), 
                  n = c(sum(table(df$q36a)),sum(table(df$q36b)),sum(table(df$q36c)),sum(table(df$q36d)),sum(table(df$q36e)),sum(table(df$q36f)),
                        sum(table(df$q37a)),sum(table(df$q37b)),sum(table(df$q37c)),sum(table(df$q37d)),sum(table(df$q37e)),sum(table(df$q37f)),
          
                                      sum(table(df$q38a)),sum(table(df$q38b)),sum(table(df$q38c)),sum(table(df$q38d)),sum(table(df$q38e)),sum(table(df$q38f))),
                  Category=c("Cigarettes","E-cigarettes","Tobacco heaters","Snus","Nicotine pouches","Pharmaceutical NRTs"),
                  Perceptions=c(rep("Health hazard",6),rep("Cancerogenicity",6),rep("Addictiveness",6))) 


plot1$Category <- factor(plot1$Category, levels=c("Pharmaceutical NRTs","Nicotine pouches","Snus","E-cigarettes","Tobacco heaters","Cigarettes"))

plot1$Perceptions <- factor(plot1$Perceptions, levels=c("Addictiveness","Health hazard","Cancerogenicity"))

ggplot(plot1,  aes(x=Mean, y=Category, fill=Perceptions)) + scale_x_continuous(breaks=(1:6*1),limits=c(0.995,6.005))+ 
  geom_point(col="red",cex=.7)+ 
  geom_errorbar(aes(xmin=Mean-(1.96*(sd/sqrt(n))), xmax=Mean+(1.96*(sd/sqrt(n)))), width=.2, 
                position=position_dodge(0.05),color="red") + 
  facet_wrap(~Perceptions,nrow=3)+ 
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank(), panel.grid.minor = element_blank(),
        legend.position = "none")



######################################
# Figure 5                           # 
######################################


df$q29a[df$q29a==0] <- NA
df$q29b[df$q29b==0] <- NA
df$q29c[df$q29c==0] <- NA
df$q29d[df$q29d==0] <- NA
df$q29e[df$q29e==0] <- NA

sum(table(df$q29a))
sum(table(df$q29b))
sum(table(df$q29c))
sum(table(df$q29d))
sum(table(df$q29e))

y <- c(rep("E-cigarettes\n(n = 249)" , 6) , rep("Tobacco heaters\n(n = 249)" , 6) , rep("Snus\n(n = 247)" , 6) , rep("Nicotine pouches\n(n = 238)" , 6), rep("NRTs\n(n = 248)" , 6) )
y <- factor(y, levels=c("NRTs\n(n = 248)","Nicotine pouches\n(n = 238)","Snus\n(n = 247)","Tobacco heaters\n(n = 249)","E-cigarettes\n(n = 249)"))
x <- rep(c("Always\n(100%)" , "Very often\n(>80%)","Often\n(50%-80%)","Sometimes\n(20%-50%)","Rarely\n(<20%)", "Never\n(0%)"), 5)
x <- factor(x, c("Always\n(100%)" , "Very often\n(>80%)","Often\n(50%-80%)","Sometimes\n(20%-50%)","Rarely\n(<20%)", "Never\n(0%)"))
value <- c(table(df$q29a),table(df$q29b),c(1,0,1,2,7,236),table(df$q29d),table(df$q29e))
data <- data.frame(x,y,value)



y <- factor(y, levels=c("Tobacco heaters\n(n = 249)","E-cigarettes\n(n = 249)","Snus\n(n = 247)","Nicotine pouches\n(n = 238)","NRTs\n(n = 248)"))
data <- data.frame(x,y,value)
data <- data %>% arrange(x)
data <- data %>% arrange(y)

data$x <- factor(x, c("Never\n(0%)", "Rarely\n(<20%)","Sometimes\n(20%-50%)","Often\n(50%-80%)","Very often\n(>80%)","Always\n(100%)"))
data$y <- factor(y, levels=c("Tobacco heaters\n(n = 249)","E-cigarettes\n(n = 249)","Snus\n(n = 247)","Nicotine pouches\n(n = 238)","NRTs\n(n = 248)"))
data$group <- c(rep(1,6),rep(2,6),rep(3,6),rep(4,6),rep(5,6))
data$pos[data$group==1] <- cumsum(data$value[data$group==1]) - (data$value[data$group==1] / 2)
data$tot[data$group==1] <- sum(data$value[data$group==1])
data$prop[data$group==1] <- paste0(data$val[data$group==1]," (",round(100* (data$val[data$group==1] / data$tot[data$group==1]),1),"%)")
data$pos[data$group==2] <- cumsum(data$value[data$group==2]) - (data$value[data$group==2] / 2)
data$tot[data$group==2] <- sum(data$value[data$group==2])
data$prop[data$group==2] <- paste0(data$val[data$group==2]," (",round(100* (data$val[data$group==2] / data$tot[data$group==2]),1),"%)")
data$pos[data$group==3] <- cumsum(data$value[data$group==3]) - (data$value[data$group==3] / 2)
data$tot[data$group==3] <- sum(data$value[data$group==3])
data$prop[data$group==3] <- paste0(data$val[data$group==3]," (",round(100* (data$val[data$group==3] / data$tot[data$group==3]),1),"%)")
data$pos[data$group==4] <- cumsum(data$value[data$group==4]) - (data$value[data$group==4] / 2)
data$tot[data$group==4] <- sum(data$value[data$group==4])
data$prop[data$group==4] <- paste0(data$val[data$group==4]," (",round(100* (data$val[data$group==4] / data$tot[data$group==4]),1),"%)")
data$pos[data$group==5] <- cumsum(data$value[data$group==5]) - (data$value[data$group==5] / 2)
data$tot[data$group==5] <- sum(data$value[data$group==5])
data$prop[data$group==5] <- paste0(data$val[data$group==5]," (",round(100* (data$val[data$group==5] / data$tot[data$group==5]),1),"%)")

data$prop[4] <- paste0(data$val[4]," (","8.0%)" )
data$prop[11] <- paste0(data$val[11]," (","6.0%)" )
data$prop[25] <- paste0(data$val[25]," (","27.0%)" )
data$pos<- data$pos *(249/data$tot)

data$prop[data$value / data$tot < .025] <- NA


ggplot(data, aes(fill=x, y=value, x=y)) + 
  geom_bar(position="fill", stat="identity") + scale_y_continuous(labels = scales::percent, breaks = c(0,.1,.2,.3,.4,.5,.6,.7,.8,.9,1))+
  scale_fill_manual(values=c("#efbd66",
                             "#FFDB7D",
                             "#fffdba",
                             "#d3e7cb",
                             "#a8d7a5",
                             "#AeD5F1"),name="") +
  labs(x = "",y="")+
  geom_text(aes(y= pos/tot, label=prop), vjust=1, color="black", size=3)+
  theme(legend.position="right")+ guides(fill = guide_legend(nrow = 6, byrow = TRUE))




######################################
# Suppl. Figure 1                    # 
######################################


x <- c(239,178,35,1,16,2,3,3,16,83,26,75,49,14,6,1,54,46,127,25,205,254,124,57,180,20)

names(x) <- levels(d.bfsrg$kt_x)


# define the a color ramp with 10 colors
cols <- colorRampPalette(colors = c("lightyellow","darkolivegreen1","green3","springgreen4"))(10)
PlotKant(rownames(x), col=FindColor(x, cols = cols, min.x=0, max.x=250),
         border="black",main="\nPharmacies in each canton")
ColorLegend(x="left", width=15000, labels=paste0(seq(0, 250, 25)),  
            cols=cols, cex=0.8, adj=c(1,0.5), frame="black", inset=c(-0.09, 0.35))
AddLakes(col="lightsteelblue1", border="black" )
AddRivers(col="steelblue1",lwd=1.25)
PlotCH(col=NA, add=TRUE, lwd=1.25, border="black")
points(sf::st_coordinates(GetMap("stkt.pnt")$geometry),
       pch=10, bg="black",cex=1.2, lwd=1.35)