
###############################################################################################################################
###############################################################################################################################

rm(list=ls(all=TRUE))
graphics.off()
library(foreign)
library(Rmisc)
library(rstatix)
library(tidyverse)

###############################################################################################################################
###############################################################################################################################

#Emyria datasets
getwd()
data1=read.dta(file="Table2_dataset.dta")
data2=read.dta(file="Table3_dataset.dta")
data3=read.dta(file="Table4_dataset.dta")
data4=read.dta(file="Table5_and_Figure1_dataset.dta")
data5=read.dta(file="Figure2_dataset.dta")
data6=read.dta(file="Figure3_dataset.dta")
data7=read.dta(file="Figure4_dataset.dta")
data8=read.dta(file="Table6_dataset.dta")
data9=read.dta(file="Table7_dataset.dta")

###############################################################################################################################
###############################################################################################################################

#Create and re-label variables
#Educational Attainment
data1$EducationalAttainment = `levels<-`(addNA(data1$Education),c(levels(data1$Education),"Not reported"))
data1$EducationalAttainment = factor(data1$EducationalAttainment,levels=c("Postgraduate degree",
                                                                          "Bachelor or honors",
                                                                          "Advanced diploma or certificate",
                                                                          "Year 12",
                                                                          "High school not completed",
                                                                          "Other",
                                                                          "Not reported"))
#Labour force status
data1$LabourForceStatus = `levels<-`(addNA(data1$EmploymentStatus),c(levels(data1$EmploymentStatus),"Not reported"))
#Medication categories
data2$Antidepressants = factor(data2$Antidepressants,levels=c("Yes",
                                                              "No"))
data2$Antipsychotics = factor(data2$Antipsychotics,levels=c("Yes",
                                                            "No"))
data2$Benzodiazepines = factor(data2$Benzodiazepines,levels=c("Yes",
                                                              "No"))
data2$CompoundAnalgesics = factor(data2$CompoundAnalgesics,levels=c("Yes",
                                                                    "No"))
data2$GabaAnalogues = factor(data2$GabaAnalogues,levels=c("Yes",
                                                          "No"))
data2$Opioids = factor(data2$Opioids,levels=c("Yes",
                                              "No"))
data2$SimpleAnalgesics = factor(data2$SimpleAnalgesics,levels=c("Yes",
                                                                "No"))
data2$OtherPainMedications = factor(data2$OtherPainMedications,levels=c("Yes",
                                                                        "No"))
data2$Cannabinoids = factor(data2$Cannabinoids,levels=c("Yes",
                                                        "No"))
data2$Others = factor(data2$Others,levels=c("Yes",
                                            "No"))
#Primary diagnosis
data3$PrimaryDiagnosisCategory = factor(data3$PrimaryDiagnosisCategory,levels=c("Pain",
                                                                                "Psychiatric",
                                                                                "Neurological",
                                                                                "Other"))

data3$PrimaryDiagnosis = ifelse(data3$PrimaryDiagnosisCategory=="Other","Other",as.character(data3$PrimaryDiagnosis))
data3$PrimaryDiagnosis = factor(data3$PrimaryDiagnosis,levels=c("Chronic non-cancer pain",
                                                                "Cancer pain",
                                                                "Migraine/headache",
                                                                "Back pain",
                                                                "Complex regional pain syndrome",
                                                                "Insomnia",
                                                                "Anxiety",
                                                                "Post-traumatic stress disorder",
                                                                "Depression",
                                                                "Parkinson's disease",
                                                                "Epilepsy",
                                                                "Multiple sclerosis",
                                                                "Other"))
#Number of medical cannabis products
data4$NumberOfMedicalCannabisProducts = ifelse(data4$NumberOfMedicalCannabisProducts==1,"One",ifelse(data4$NumberOfMedicalCannabisProducts==2,"Two","Three or more"))
data4$NumberOfMedicalCannabisProducts = factor(data4$NumberOfMedicalCannabisProducts,levels=c("One",
                                                                                              "Two",
                                                                                              "Three or more"))
#Medical cannabis categories
data4$MedicalCannabisCategory = factor(data4$MedicalCannabisCategory,levels=c("Balanced",
                                                                              "CBD only",
                                                                              "THC-dominant",
                                                                              "CBD-dominant",
                                                                              "THC only"))


###############################################################################################################################
###############################################################################################################################

#Regression analysis
#Table 2
#Gender
df = data1
df %>%
  group_by(Gender) %>%
  summarise(n = n()) %>%
  mutate(freq = sprintf("%.1f",(n/sum(n))*100))
#Age
df %>%
  rstatix::get_summary_stats(Age,
                             type="mean_sd")
#Educational attainment
df %>%
  group_by(EducationalAttainment) %>%
  summarise(n = n()) %>%
  mutate(freq = sprintf("%.1f",(n/sum(n))*100))
#Labour force status
df %>%
  group_by(LabourForceStatus) %>%
  summarise(n = n()) %>%
  mutate(freq = sprintf("%.1f",(n/sum(n))*100))

#Table 3
#Medication categories
#Simple analgesics
df = data2
df %>%
  group_by(SimpleAnalgesics) %>%
  summarise(n = n()) %>%
  mutate(freq = sprintf("%.1f",(n/sum(n))*100))
#Opioids
df %>%
  group_by(Opioids) %>%
  summarise(n = n()) %>%
  mutate(freq = sprintf("%.1f",(n/sum(n))*100))
#Antidepressants
df %>%
  group_by(Antidepressants) %>%
  summarise(n = n()) %>%
  mutate(freq = sprintf("%.1f",(n/sum(n))*100))
#Benzodiazepines
df %>%
  group_by(Benzodiazepines) %>%
  summarise(n = n()) %>%
  mutate(freq = sprintf("%.1f",(n/sum(n))*100))
#GABA analogues
df %>%
  group_by(GabaAnalogues) %>%
  summarise(n = n()) %>%
  mutate(freq = sprintf("%.1f",(n/sum(n))*100))
#Other pain medications
df %>%
  group_by(OtherPainMedications) %>%
  summarise(n = n()) %>%
  mutate(freq = sprintf("%.1f",(n/sum(n))*100))
#Antipsychotics
df %>%
  group_by(Antipsychotics) %>%
  summarise(n = n()) %>%
  mutate(freq = sprintf("%.1f",(n/sum(n))*100))
#Compound analgesics
df %>%
  group_by(CompoundAnalgesics) %>%
  summarise(n = n()) %>%
  mutate(freq = sprintf("%.1f",(n/sum(n))*100))
#Cannabinoids
df %>%
  group_by(Cannabinoids) %>%
  summarise(n = n()) %>%
  mutate(freq = sprintf("%.1f",(n/sum(n))*100))
#Others
df %>%
  group_by(Others) %>%
  summarise(n = n()) %>%
  mutate(freq = sprintf("%.1f",(n/sum(n))*100))
#Number of medications taken daily
df %>%
  rstatix::get_summary_stats(NumberOfMedications,
                             type="common",
                             show=c("median","min","max"))

#Table 4
#Primary diagnosis
df = data3
df %>%
  group_by(PrimaryDiagnosisCategory,PrimaryDiagnosis) %>%
  summarise(n = n()) %>%
  mutate(freq = sprintf("%.1f",(n/3961)*100)) %>%
  print(n=Inf)
#Number of comorbidities
df %>%
  rstatix::get_summary_stats(NumberOfComorbidities,
                             type="mean_sd")

#Table 5
#Number of products
df = data4
df %>%
  group_by(NumberOfMedicalCannabisProducts) %>%
  summarise(n = n()) %>%
  mutate(freq = sprintf("%.1f",(n/sum(n))*100))
#Categories
df %>%
  group_by(MedicalCannabisCategory) %>%
  summarise(n = n()) %>%
  mutate(freq = sprintf("%.1f",(n/sum(n))*100))

#Figure 1
df = data4
df1 = select(df,MedicalCannabisTotalCBDDose,Followup)
names(df1) = c("X1","X2")
df1$X3 = "CBD"
df2 = df1
df1 = select(df,MedicalCannabisTotalTHCDose,Followup)
names(df1) = c("X1","X2")
df1$X3 = "THC"
df2 = bind_rows(df2,df1)
df1 = select(df,MedicalCannabisTotalCBDDose,MedicalCannabisTotalTHCDose,Followup)
names(df1) = c("X0","X1","X2")
df1$X1 = df1$X0+df1$X1
df1$X3 = "Total"
df2 = bind_rows(df2,select(df1,X1,X2,X3))
scaleFUN <- function(x) sprintf("%.0f", x)
df2$X2 = factor(df2$X2,levels=c("3 Month","6 Month","12 Month","18 Month","24 Month"))
df2_sum = summarySE(df2,measurevar="X1",groupvar=c("X3","X2"))
df2_sum$X4 = paste(df2_sum$X2," (N=",df2_sum$N,")",sep="")
colour_gradient = scales::seq_gradient_pal("white","black","Lab")(seq(0,1,length.out=6))[2:6]
figure1 = ggplot(subset(df2_sum,X3!="Total"),aes(x=X2,y=X1,group=X3))+
  annotate("rect",ymin=0,ymax=100,xmin=c(0.501,1.501,2.501,3.501,4.501),xmax=c(1.499,2.499,3.499,4.499,5.499),fill=colour_gradient,fill="white",alpha=0.1)+
  geom_bar(stat="identity",aes(fill=X3),position="stack",width=.5)+
  geom_point(data=subset(df2_sum,X3=="Total"),size=1,colour="black")+
  geom_errorbar(data=subset(df2_sum,X3=="Total"),aes(ymin=X1-ci,ymax=X1+ci),width=.3)+
  labs(x="Follow-up",y="Mean daily doses (in mg)")+
  theme(text=element_text(colour="black",size=10),panel.background=element_rect(fill="white"),axis.line=element_line(colour="black"),panel.grid.major.y=element_blank(),panel.grid.major.x=element_blank(),panel.grid.minor.x=element_blank(),legend.position="bottom",legend.box.background=element_rect(),legend.box.margin=margin(1,1,1,1),legend.key=element_rect(fill="white"),legend.title=element_blank(),legend.spacing.x=unit(0.1,'cm'))+
  scale_fill_manual(values=c("palegreen4","gray50"))+
  scale_y_continuous(breaks=seq(min(0),max(100),by=20),labels=scaleFUN)+
  scale_x_discrete(labels=df2_sum$X4)+
  coord_cartesian(ylim=c(0,100))
ggsave(filename="Figure1.tiff",plot=figure1,device="tiff",path=paste(dirname(getwd()),"Figures",sep="/"),width=2250,height=750,units="px",dpi=300)

#Figure 2
#Panel A
df = data5
df1 = select(df,SF36PhysicalFunctioningScore,Followup)
names(df1) = c("X1","X2")
df1$X3 = "Physical functioning"
df2 = df1
df1 = select(df,SF36RolePhysicalScore,Followup)
names(df1) = c("X1","X2")
df1$X3 = "Role-physical"
df2 = bind_rows(df2,df1)
df1 = select(df,SF36BodilyPainScore,Followup)
names(df1) = c("X1","X2")
df1$X3 = "Bodily pain"
df2 = bind_rows(df2,df1)
df1 = select(df,SF36GeneralHealthScore,Followup)
names(df1) = c("X1","X2")
df1$X3 = "General health"
df2 = bind_rows(df2,df1)
scaleFUN <- function(x) sprintf("%.0f", x)
df2$X2 = factor(df2$X2,levels=c("Baseline","3 Month","6 Month","12 Month","18 Month","24 Month"))
df2_sum = summarySE(df2,measurevar="X1",groupvar=c("X3","X2"))
df2_sum$X4 = paste(df2_sum$X2," (N=",df2_sum$N,")",sep="")
colour_gradient = scales::seq_gradient_pal("white","black","Lab")(seq(0,1,length.out=6))
figure2 = ggplot(df2_sum,aes(x=X2,y=X1,group=X3))+
  annotate("rect",ymin=0,ymax=60,xmin=c(0.501,1.501,2.501,3.501,4.501,5.501),xmax=c(1.499,2.499,3.499,4.499,5.499,6.499),fill=colour_gradient,fill="white",alpha=0.1)+
  geom_point(size=1,aes(colour=X3),position=position_dodge(width=.3))+
  geom_line(size=1,aes(colour=X3),position=position_dodge(width=.3))+
  geom_errorbar(aes(ymin=X1-ci,ymax=X1+ci),width=.3,position=position_dodge(width=.3))+
  labs(x="Follow-up",y="Mean survey scores")+
  theme(text=element_text(colour="black",size=10),panel.background=element_rect(fill="white"),axis.line=element_line(colour="black"),panel.grid.major.y=element_blank(),panel.grid.major.x=element_blank(),panel.grid.minor.x=element_blank(),legend.position="bottom",legend.box.background=element_rect(),legend.box.margin=margin(1,1,1,1),legend.key=element_rect(fill="white"),legend.title=element_blank(),legend.spacing.x=unit(0.1,'cm'))+
  scale_colour_manual(values=c("purple4","brown4","orange4","gray50"))+
  scale_y_continuous(breaks=seq(min(0),max(60),by=15),labels=scaleFUN)+
  scale_x_discrete(labels=df2_sum$X4)+
  coord_cartesian(ylim=c(0,60))
ggsave(filename="Figure2PanelA.tiff",plot=figure2,device="tiff",path=paste(dirname(getwd()),"Figures",sep="/"),width=2250,height=750,units="px",dpi=300)
#Panel B
df = data5
df1 = select(df,SF36VitalityScore,Followup)
names(df1) = c("X1","X2")
df1$X3 = "Vitality"
df2 = df1
df1 = select(df,SF36SocialFunctioningScore,Followup)
names(df1) = c("X1","X2")
df1$X3 = "Social functioning"
df2 = bind_rows(df2,df1)
df1 = select(df,SF36RoleEmotionalScore,Followup)
names(df1) = c("X1","X2")
df1$X3 = "Role-emotional"
df2 = bind_rows(df2,df1)
df1 = select(df,SF36MentalHealthScore,Followup)
names(df1) = c("X1","X2")
df1$X3 = "Mental health"
df2 = bind_rows(df2,df1)
scaleFUN <- function(x) sprintf("%.0f", x)
df2$X2 = factor(df2$X2,levels=c("Baseline","3 Month","6 Month","12 Month","18 Month","24 Month"))
df2_sum = summarySE(df2,measurevar="X1",groupvar=c("X3","X2"))
df2_sum$X4 = paste(df2_sum$X2," (N=",df2_sum$N,")",sep="")
colour_gradient = scales::seq_gradient_pal("white","black","Lab")(seq(0,1,length.out=6))
figure2 = ggplot(df2_sum,aes(x=X2,y=X1,group=X3))+
  annotate("rect",ymin=20,ymax=80,xmin=c(0.501,1.501,2.501,3.501,4.501,5.501),xmax=c(1.499,2.499,3.499,4.499,5.499,6.499),fill=colour_gradient,fill="white",alpha=0.1)+
  geom_point(size=1,aes(colour=X3),position=position_dodge(width=.3))+
  geom_line(size=1,aes(colour=X3),position=position_dodge(width=.3))+
  geom_errorbar(aes(ymin=X1-ci,ymax=X1+ci),width=.3,position=position_dodge(width=.3))+
  labs(x="Follow-up",y="Mean survey scores")+
  theme(text=element_text(colour="black",size=10),panel.background=element_rect(fill="white"),axis.line=element_line(colour="black"),panel.grid.major.y=element_blank(),panel.grid.major.x=element_blank(),panel.grid.minor.x=element_blank(),legend.position="bottom",legend.box.background=element_rect(),legend.box.margin=margin(1,1,1,1),legend.key=element_rect(fill="white"),legend.title=element_blank(),legend.spacing.x=unit(0.1,'cm'))+
  scale_colour_manual(values=c("red4","blue4","palegreen4","yellow4"))+
  scale_y_continuous(breaks=seq(min(20),max(80),by=15),labels=scaleFUN)+
  scale_x_discrete(labels=df2_sum$X4)+
  coord_cartesian(ylim=c(20,80))
ggsave(filename="Figure2PanelB.tiff",plot=figure2,device="tiff",path=paste(dirname(getwd()),"Figures",sep="/"),width=2250,height=750,units="px",dpi=300)

#Figure 3
#Panel A
df = subset(data6,!is.na(BPIPainSeverityScore)
            &!is.na(BPIPainInterferenceScore))
df1 = select(df,BPIPainSeverityScore,Followup)
names(df1) = c("X1","X2")
df1$X3 = "Pain severity"
df2 = df1
df1 = select(df,BPIPainInterferenceScore,Followup)
names(df1) = c("X1","X2")
df1$X3 = "Pain interference"
df2 = bind_rows(df2,df1)
scaleFUN <- function(x) sprintf("%.0f", x)
df2$X2 = factor(df2$X2,levels=c("Baseline","3 Month","6 Month","12 Month","18 Month","24 Month"))
df2_sum = summarySE(df2,measurevar="X1",groupvar=c("X3","X2"))
df2_sum$X4 = paste(df2_sum$X2," (N=",df2_sum$N,")",sep="")
colour_gradient = scales::seq_gradient_pal("white","black","Lab")(seq(0,1,length.out=6))
figure3 = ggplot(df2_sum,aes(x=X2,y=X1,group=X3))+
  annotate("rect",ymin=3,ymax=7,xmin=c(0.501,1.501,2.501,3.501,4.501,5.501),xmax=c(1.499,2.499,3.499,4.499,5.499,6.499),fill=colour_gradient,fill="white",alpha=0.1)+
  geom_point(size=1,aes(colour=X3),position=position_dodge(width=.3))+
  geom_line(size=1,aes(colour=X3),position=position_dodge(width=.3))+
  geom_errorbar(aes(ymin=X1-ci,ymax=X1+ci),width=.3,position=position_dodge(width=.3))+
  labs(x="Follow-up",y="Mean survey scores")+
  theme(text=element_text(colour="black",size=10),panel.background=element_rect(fill="white"),axis.line=element_line(colour="black"),panel.grid.major.y=element_blank(),panel.grid.major.x=element_blank(),panel.grid.minor.x=element_blank(),legend.position="bottom",legend.box.background=element_rect(),legend.box.margin=margin(1,1,1,1),legend.key=element_rect(fill="white"),legend.title=element_blank(),legend.spacing.x=unit(0.1,'cm'))+
  scale_colour_manual(values=c("purple4","brown4"))+
  scale_y_continuous(breaks=seq(min(3),max(7),by=1),labels=scaleFUN)+
  scale_x_discrete(labels=df2_sum$X4)+
  coord_cartesian(ylim=c(3,7))
ggsave(filename="Figure3PanelA.tiff",plot=figure3,device="tiff",path=paste(dirname(getwd()),"Figures",sep="/"),width=2250,height=750,units="px",dpi=300)
#Panel B
df = subset(data6,!is.na(DASS21AnxietyScore)
            &!is.na(DASS21DepressionScore)
            &!is.na(DASS21StressScore))
df1 = select(df,DASS21AnxietyScore,Followup)
names(df1) = c("X1","X2")
df1$X3 = "Anxiety"
df2 = df1
df1 = select(df,DASS21DepressionScore,Followup)
names(df1) = c("X1","X2")
df1$X3 = "Depression"
df2 = bind_rows(df2,df1)
df1 = select(df,DASS21StressScore,Followup)
names(df1) = c("X1","X2")
df1$X3 = "Stress"
df2 = bind_rows(df2,df1)
scaleFUN <- function(x) sprintf("%.0f", x)
df2$X2 = factor(df2$X2,levels=c("Baseline","3 Month","6 Month","12 Month","18 Month","24 Month"))
df2_sum = summarySE(df2,measurevar="X1",groupvar=c("X3","X2"))
df2_sum$X4 = paste(df2_sum$X2," (N=",df2_sum$N,")",sep="")
colour_gradient = scales::seq_gradient_pal("white","black","Lab")(seq(0,1,length.out=6))
figure3 = ggplot(df2_sum,aes(x=X2,y=X1,group=X3))+
  annotate("rect",ymin=4,ymax=20,xmin=c(0.501,1.501,2.501,3.501,4.501,5.501),xmax=c(1.499,2.499,3.499,4.499,5.499,6.499),fill=colour_gradient,fill="white",alpha=0.1)+
  geom_point(size=1,aes(colour=X3),position=position_dodge(width=.3))+
  geom_line(size=1,aes(colour=X3),position=position_dodge(width=.3))+
  geom_errorbar(aes(ymin=X1-ci,ymax=X1+ci),width=.3,position=position_dodge(width=.3))+
  labs(x="Follow-up",y="Mean survey scores")+
  theme(text=element_text(colour="black",size=10),panel.background=element_rect(fill="white"),axis.line=element_line(colour="black"),panel.grid.major.y=element_blank(),panel.grid.major.x=element_blank(),panel.grid.minor.x=element_blank(),legend.position="bottom",legend.box.background=element_rect(),legend.box.margin=margin(1,1,1,1),legend.key=element_rect(fill="white"),legend.title=element_blank(),legend.spacing.x=unit(0.1,'cm'))+
  scale_colour_manual(values=c("red4","blue4","palegreen4"))+
  scale_y_continuous(breaks=seq(min(4),max(20),by=4),labels=scaleFUN)+
  scale_x_discrete(labels=df2_sum$X4)+
  coord_cartesian(ylim=c(4,20))
ggsave(filename="Figure3PanelB.tiff",plot=figure3,device="tiff",path=paste(dirname(getwd()),"Figures",sep="/"),width=2250,height=750,units="px",dpi=300)
#Panel C
df = subset(data6,!is.na(ISIInsomniaScore))
df1 = select(df,ISIInsomniaScore,Followup)
names(df1) = c("X1","X2")
df1$X3 = "Insomnia"
df2 = df1
scaleFUN <- function(x) sprintf("%.0f", x)
df2$X2 = factor(df2$X2,levels=c("Baseline","3 Month","6 Month","12 Month","18 Month","24 Month"))
df2_sum = summarySE(df2,measurevar="X1",groupvar=c("X3","X2"))
df2_sum$X4 = paste(df2_sum$X2," (N=",df2_sum$N,")",sep="")
colour_gradient = scales::seq_gradient_pal("white","black","Lab")(seq(0,1,length.out=6))
figure3 = ggplot(df2_sum,aes(x=X2,y=X1,group=X3))+
  annotate("rect",ymin=6,ymax=18,xmin=c(0.501,1.501,2.501,3.501,4.501,5.501),xmax=c(1.499,2.499,3.499,4.499,5.499,6.499),fill=colour_gradient,fill="white",alpha=0.1)+
  geom_point(size=1,aes(colour=X3),position=position_dodge(width=.3))+
  geom_line(size=1,aes(colour=X3),position=position_dodge(width=.3))+
  geom_errorbar(aes(ymin=X1-ci,ymax=X1+ci),width=.3,position=position_dodge(width=.3))+
  labs(x="Follow-up",y="Mean survey scores")+
  theme(text=element_text(colour="black",size=10),panel.background=element_rect(fill="white"),axis.line=element_line(colour="black"),panel.grid.major.y=element_blank(),panel.grid.major.x=element_blank(),panel.grid.minor.x=element_blank(),legend.position="bottom",legend.box.background=element_rect(),legend.box.margin=margin(1,1,1,1),legend.key=element_rect(fill="white"),legend.title=element_blank(),legend.spacing.x=unit(0.1,'cm'))+
  scale_colour_manual(values=c("yellow4"))+
  scale_y_continuous(breaks=seq(min(6),max(18),by=3),labels=scaleFUN)+
  scale_x_discrete(labels=df2_sum$X4)+
  coord_cartesian(ylim=c(6,18))
ggsave(filename="Figure3PanelC.tiff",plot=figure3,device="tiff",path=paste(dirname(getwd()),"Figures",sep="/"),width=2250,height=750,units="px",dpi=300)

#Figure 4
df = data7
df1 = select(df,CGIEfficacyIndexScore,Followup)
names(df1) = c("X1","X2")
df1$X3 = "Efficacy index"
df2 = df1
df1 = select(df,CGIGlobalImprovementScore,Followup)
names(df1) = c("X1","X2")
df1$X3 = "Global improvement"
df2 = bind_rows(df2,df1)
df1 = select(df,CGISeverityOfIllnessScore,Followup)
names(df1) = c("X1","X2")
df1$X3 = "Severity of illness"
df2 = bind_rows(df2,df1)
scaleFUN <- function(x) sprintf("%.0f", x)
df2$X2 = factor(df2$X2,levels=c("Baseline","3 Month","6 Month","12 Month","18 Month","24 Month"))
df2_sum = summarySE(df2,measurevar="X1",groupvar=c("X3","X2"))
df2_sum$X4 = paste(df2_sum$X2," (N=",df2_sum$N,")",sep="")
colour_gradient = scales::seq_gradient_pal("white","black","Lab")(seq(0,1,length.out=6))
figure4 = ggplot(df2_sum,aes(x=X2,y=X1,group=X3))+
  annotate("rect",ymin=1,ymax=9,xmin=c(0.501,1.501,2.501,3.501,4.501,5.501),xmax=c(1.499,2.499,3.499,4.499,5.499,6.499),fill=colour_gradient,fill="white",alpha=0.1)+
  geom_point(size=1,aes(colour=X3),position=position_dodge(width=.3))+
  geom_line(size=1,aes(colour=X3),position=position_dodge(width=.3))+
  geom_errorbar(aes(ymin=X1-ci,ymax=X1+ci),width=.3,position=position_dodge(width=.3))+
  labs(x="Follow-up",y="Mean survey scores")+
  theme(text=element_text(colour="black",size=10),panel.background=element_rect(fill="white"),axis.line=element_line(colour="black"),panel.grid.major.y=element_blank(),panel.grid.major.x=element_blank(),panel.grid.minor.x=element_blank(),legend.position="bottom",legend.box.background=element_rect(),legend.box.margin=margin(1,1,1,1),legend.key=element_rect(fill="white"),legend.title=element_blank(),legend.spacing.x=unit(0.1,'cm'))+
  scale_colour_manual(values=c("red4","blue4","palegreen4"))+
  scale_y_continuous(breaks=seq(min(1),max(9),by=2),labels=scaleFUN)+
  scale_x_discrete(labels=df2_sum$X4)+
  coord_cartesian(ylim=c(1,9))
ggsave(filename="Figure4.tiff",plot=figure4,device="tiff",path=paste(dirname(getwd()),"Figures",sep="/"),width=2250,height=750,units="px",dpi=300)

#Table 6
df = data8
df1 = select(df,SF36PhysicalFunctioningScore,Followup)
names(df1)[1] = "Score"
df1$Variable = "Physical functioning"
df2 = df1
df1 = select(df,SF36RolePhysicalScore,Followup)
names(df1)[1] = "Score"
df1$Variable = "Role-physical"
df2 = bind_rows(df2,df1)
df1 = select(df,SF36BodilyPainScore,Followup)
names(df1)[1] = "Score"
df1$Variable = "Bodily pain"
df2 = bind_rows(df2,df1)
df1 = select(df,SF36GeneralHealthScore,Followup)
names(df1)[1] = "Score"
df1$Variable = "General health"
df2 = bind_rows(df2,df1)
df1 = select(df,SF36VitalityScore,Followup)
names(df1)[1] = "Score"
df1$Variable = "Vitality"
df2 = bind_rows(df2,df1)
df1 = select(df,SF36SocialFunctioningScore,Followup)
names(df1)[1] = "Score"
df1$Variable = "Social functioning"
df2 = bind_rows(df2,df1)
df1 = select(df,SF36RoleEmotionalScore,Followup)
names(df1)[1] = "Score"
df1$Variable = "Role-emotional"
df2 = bind_rows(df2,df1)
df1 = select(df,SF36MentalHealthScore,Followup)
names(df1)[1] = "Score"
df1$Variable = "Mental health"
df2 = bind_rows(df2,df1)
df1 = select(df,BPIPainSeverityScore,Followup)
names(df1)[1] = "Score"
df1$Variable = "Pain severity"
df2 = bind_rows(df2,df1)
df1 = select(df,BPIPainInterferenceScore,Followup)
names(df1)[1] = "Score"
df1$Variable = "Pain interference"
df2 = bind_rows(df2,df1)
df1 = select(df,DASS21AnxietyScore,Followup)
names(df1)[1] = "Score"
df1$Variable = "Anxiety"
df2 = bind_rows(df2,df1)
df1 = select(df,DASS21DepressionScore,Followup)
names(df1)[1] = "Score"
df1$Variable = "Depression"
df2 = bind_rows(df2,df1)
df1 = select(df,DASS21StressScore,Followup)
names(df1)[1] = "Score"
df1$Variable = "Stress"
df2 = bind_rows(df2,df1)
df1 = select(df,ISIInsomniaScore,Followup)
names(df1)[1] = "Score"
df1$Variable = "Insomnia"
df2 = bind_rows(df2,df1)
df1 = select(df,CGIEfficacyIndexScore,Followup)
names(df1)[1] = "Score"
df1$Variable = "Efficacy index"
df2 = bind_rows(df2,df1)
df1 = select(df,CGIGlobalImprovementScore,Followup)
names(df1)[1] = "Score"
df1$Variable = "Global improvement"
df2 = bind_rows(df2,df1)
df1 = select(df,CGISeverityOfIllnessScore,Followup)
names(df1)[1] = "Score"
df1$Variable = "Severity of illness"
df2 = bind_rows(df2,df1)
df1 = select(df,NumberOfMedications,Followup)
names(df1)[1] = "Score"
df1$Variable = "Number of medications"
df2 = bind_rows(df2,df1)
df2$Variable = factor(df2$Variable,levels=c("Bodily pain",
                                            "General health",
                                            "Physical functioning",
                                            "Role-physical",
                                            "Mental health",
                                            "Role-emotional",
                                            "Social functioning",
                                            "Vitality",
                                            "Anxiety",
                                            "Depression",
                                            "Stress",
                                            "Pain severity",
                                            "Pain interference",
                                            "Insomnia",
                                            "Efficacy index",
                                            "Global improvement",
                                            "Severity of illness",
                                            "Number of medications"))
df2$Followup = factor(df2$Followup,levels=c("Baseline",
                                            "3 Month",
                                            "12 Month",
                                            "24 Month"))
df3 = df2 %>%
  group_by(Variable,Followup) %>%
  rstatix::get_summary_stats(Score,type ="mean_ci") %>%
  select(Variable,Followup,n,mean,ci)
df2 = subset(df2,!Variable%in%c("Efficacy index","Global improvement"))
df2 = df2 %>% 
  group_by(Variable) %>% 
  t_test(Score~Followup,ref.group="Baseline") %>%
  mutate(Followup = group2) %>%
  select(Variable,Followup,p)
left_join(df3,df2) %>% 
  pivot_wider(names_from=Followup,
              names_glue = "{Followup} {.value}",
              values_from=c(n,mean,ci,p),
              names_vary = "slowest")

#Table 7
df = cbind(data9 %>%
             select(!matches("^(AdverseEventSeverity)")) %>% 
             gather(na.rm = TRUE) %>%
             mutate(AdverseEvent = value) %>%
             select(AdverseEvent),
           data9 %>% 
             select(matches("^(AdverseEventSeverity)")) %>% 
             gather(na.rm = TRUE) %>%
             mutate(Severity = as.factor(value)) %>%
             select(Severity))
df1 = df %>%
  group_by(AdverseEvent) %>%
  summarise(n = n()) %>%
  subset(n>=40) %>%
  select(AdverseEvent)
df$AdverseEvent = as.factor(ifelse(df$AdverseEvent%in%as.vector(df1$AdverseEvent),df$AdverseEvent,"Other"))
df$AdverseEvent = factor(df$AdverseEvent,levels=c("Sedation/sleepiness",
                                                  "Dry mouth",
                                                  "Lethargy/tiredness",
                                                  "Dizziness",
                                                  "Nausea",
                                                  "Concentration difficulty",
                                                  "Feeling high",
                                                  "Diarrhoea/loose stools",
                                                  "Increased appetite",
                                                  "Headache",
                                                  "Anxiety/panic attack",
                                                  "Vivid dreams",
                                                  "Hallucination",
                                                  "Impaired coordination",
                                                  "Other"))
df1 = df %>%
  group_by(AdverseEvent,Severity) %>%
  summarise(n = n()) %>%
  pivot_wider(names_from=Severity,
              names_glue = "{Severity} {.value}",
              values_from=c(n),
              names_vary = "slowest")
df2 = df %>%
  group_by(AdverseEvent) %>%
  summarise(n = n()) %>%
  mutate(freq = sprintf("%.1f",(n/3961)*100),
         Severity = "Total") %>%
  pivot_wider(names_from=Severity,
              names_glue = "{Severity} {.value}",
              values_from=c(n,freq),
              names_vary = "slowest")
left_join(df1,df2)

###############################################################################################################################
###############################################################################################################################
