#Script to obtain the figures

library(readr)
library(readxl)
library(xlsx)
library(tidyverse)
library(biscale)
library(cowplot)
library(sf)
library(ggpubr)
library(ggplot2)
library(deeptime)
library(poLCAParallel)
library(broom)
library(knitr)
library(flextable)
library(kableExtra)
library(useful)
library(tmap)
library(maptools)
library(regclass)
library(sf)
library(patchwork)
library(Cairo)
library(ggpubr)
library(moonBook)

#Figure 1

summary1 <- read.xlsx("Source Data File.xlsx",sheetIndex = 1)

p1 <- ggplot(data=summary1, aes(x=Age.Group, y=Frequency, fill=Conditions)) +
  geom_bar(stat="identity", position=position_dodge())+
  geom_text(aes(label=Frequency), vjust=1.6, color="white",
            position = position_dodge(0.9), size=3.5)+
  labs(x = "Age Group", y="", fill = "Conditions")+
  #scale_x_discrete(labels=c("<16", "16-44","45-64","65-79", ">79"))
  scale_fill_manual(values=c("#8da0cb", "#66c2a5","#FC8D62"))+
  #scale_color_brewer(palette = "Set2")+
  theme_light(base_size = 18)

p1

#Figure 2

labels <- c("Cardiometabolic + Cardiovascular + \nMusculoskeletal+ Respiratory \n+ Retinopathy + Renal (Complex)",
            "Dementia + Mental","Cardiometabolic + Musculoskeletal",
            "Cardiometabolic + Cardiovascular + \nMusculoskeletal + Respiratory \n(Complex)","Cardiovascular",
            "Dependence + Liver Disease","Cardiometabolic + Cardiovascular",
            "Cardiometabolic + Cardiovascular + \nMusculoskeletal + Mental \n(Complex)",
            "Cardiometabolic + Mental + \nMusculoskeletal + Respiratory \n(Complex)",
            "Cardiometabolic","Hypertension + Dyslipidemia \n+ Arthrosis","Asthma + Dependence",
            "Musculoskeletal","Hypthyroidism + Dyslipidemia","Mental + Developmental Problems",
            "Asthma + Obesity","Dermatitis + Food Intolerance","Digestive","Dermatitis + Digestive",
            "Developmental Problems \n+ Neurological","Mental","Respiratory","Asthma + Dermatitis",
            "Cardiometabolic","Asthma + Food Intolerance")
pcm = read.xlsx("Source Data File.xlsx",sheetIndex = 2)
colours = c( "#A54657",  "#582630", "#F7EE7F", "#4DAA57","#F1A66A", 
             "#A54657",  "#582630", "#F7EE7F", "#4DAA57","#F1A66A")
pcm$Group <- factor(pcm$Group,levels=unique(pcm$Group))

p2 <- ggplot(pcm, aes(x = Group, y = factor(Patterns,level = rev(unique(pcm$Patterns))))) + 
  geom_point(aes(size = Prevalence, fill = Group), alpha = 0.75, shape = 21) +
  scale_size_continuous(limits = c(0.000001, 52), range = c(1,30), breaks = c(1,25,50)) + 
  labs( x= "", y = "", size = "Relative Prevalence (%)", fill = "")  + 
  theme(legend.key=element_blank(), 
        axis.text.x = element_text(colour = "black", size = 16, face = "bold", angle = 60, hjust = 1), 
        axis.text.y = element_text(colour = "black", face = "bold", size = 17), 
        legend.text = element_text(size = 10, face ="bold", colour ="black"), 
        legend.title = element_text(size = 14, face = "bold"), 
        panel.background = element_blank(), panel.border = element_rect(colour = "black", fill = NA, size = 1.2), 
        legend.position = "right") +  
  scale_fill_manual(values = colours, guide = FALSE) + 
  scale_y_discrete(limits = rev(levels(unique(pcm$Patterns))),
                   labels=label_wrap(48)) + 
  geom_vline(xintercept = 5.5, linetype="dotted", 
             color = "black", size=2.5)+
  theme(panel.grid.major = element_line(size = 0.4,color="grey"))

p2

#Figure 5

df = read.xlsx("Source Data File.xlsx",sheetIndex = 5)
colnames(df) <- df[1,]
df <- df[-1, ]

df2 <- df[,1:4]
df2[,4] <- as.numeric(df2[,4])
df2$Age_Group <- factor(df2$Age_Group, levels=c("<16","16-44","45-64","65-79",">79"))
p1 <- ggplot(data=df2, aes(x=Age_Group, y=GP_visits, fill=Cluster)) +
  geom_bar(stat="identity",position="Dodge") +
  labs(title = "Mean number of GP visits by sex, age, and SES area", x = "Age Group", y="", fill = "Local Health Area")+
  facet_wrap(vars(Sex),labeller = labeller(Sex = c("Mujer"="Women","Hombre" ="Men")))+
  scale_fill_manual(values=c("#fc8d62", "#8da0cb", "#66c2a5")) +
  theme_light()

df2 <- df[,5:8]
df2[,4] <- as.numeric(df2[,4])
df2$Age_Group <- factor(df2$Age_Group, levels=c("<16","16-44","45-64","65-79",">79"))
p2 <- ggplot(data=df2, aes(x=Age_Group, y=Hospital_visits, fill=Cluster)) +
  geom_bar(stat="identity",position="Dodge") +
  labs(title = "Mean number of hospital visits by sex, age, and SES area", x = "Age Group", y="", fill = "Local Health Area")+
  facet_wrap(vars(Sex),labeller = labeller(Sex = c("Mujer"="Women","Hombre" ="Men")))+
  scale_fill_manual(values=c("#fc8d62", "#8da0cb", "#66c2a5")) +
  theme_light()

df2 <- df[,9:12]
df2[,4] <- as.numeric(df2[,4])
df2$Age_Group <- factor(df2$Age_Group, levels=c("<16","16-44","45-64","65-79",">79"))
p3 <- ggplot(data=df2, aes(x=Age_Group, y=AE_visits, fill=Cluster)) +
  geom_bar(stat="identity",position="Dodge") +
  labs(title = "Mean number of A&E visits by sex, age, and SES area", x = "Age Group", y="", fill = "Local Health Area")+
  facet_wrap(vars(Sex),labeller = labeller(Sex = c("Mujer"="Women","Hombre" ="Men")))+
  scale_fill_manual(values=c("#fc8d62", "#8da0cb", "#66c2a5")) +
  theme_light()

ggarrange(p1,p2,p3,nrow=3,ncol=1)

#Figure 6

df = read.xlsx("Source Data File.xlsx",sheetIndex = 6)
colnames(df) <- df[1,]
df <- df[-1, ]

df2 <- df[,1:3]
df2[,3] <- as.numeric(df2[,3])
df2 <- df2[order(df2$GP_visits,decreasing=TRUE),]
df2$Pattern <- factor(df2$Pattern, levels = rev(unique(df2$Pattern)))
p1 <- df2 %>%
  ggplot( aes(x=Pattern, y=GP_visits, group=Cluster, color=Cluster)) +
  labs(title = "Mean number of GP visits by pattern and SES area", x= "Pattern", y="",color="Local Health Area") +
  theme(axis.text.x = element_text(face="bold", size=7)) +
  geom_line() +
  scale_fill_manual(values=c("#fc8d62", "#8da0cb", "#66c2a5")) +
  theme_light()

df2 <- df[,4:6]
df2[,3] <- as.numeric(df2[,3])
df2 <- df2[order(df2$Hospital_visits,decreasing=TRUE),]
df2$Pattern <- factor(df2$Pattern, levels = rev(unique(df2$Pattern)))
p2 <- df2 %>%
  ggplot( aes(x=Pattern, y=Hospital_visits, group=Cluster, color=Cluster)) +
  labs(title = "Mean number of hospital visits by pattern and SES area", x= "Pattern", y="",color="Local Health Area") +
  theme(axis.text.x = element_text(face="bold", size=7)) +
  geom_line() +
  scale_fill_manual(values=c("#fc8d62", "#8da0cb", "#66c2a5")) +
  theme_light()

df2 <- df[,7:9]
df2[,3] <- as.numeric(df2[,3])
df2 <- df2[order(df2$AE_visits,decreasing=TRUE),]
df2$Pattern <- factor(df2$Pattern, levels = rev(unique(df2$Pattern)))
p3 <- df2 %>%
  ggplot( aes(x=Pattern, y=AE_visits, group=Cluster, color=Cluster)) +
  labs(title = "Mean number of A&E visits by pattern and SES area", x= "Pattern", y="",color="Local Health Area") +
  theme(axis.text.x = element_text(face="bold", size=7)) +
  geom_line() + 
  scale_fill_manual(values=c("#fc8d62", "#8da0cb", "#66c2a5")) +
  theme_light()

ggarrange(p1,p2,p3,nrow=3,ncol=1)

#Figure 7

#Function for figure 7
OR.sub=function(result,type=1,xlab="",ylab="",show.OR=TRUE,show.CI=FALSE,
                sig.level=1,cex=1.2,lwd=2,pch=18,col=NULL,...){
  result=result[result[[4]]<=sig.level,]
  count=length(result[,1])
  if(count<1) {
    cat("No variable to be plotted found")
    return(invisible())
  }
  if(is.null(col) | (length(col)!=2) ) {
    if(type==3) {col1="salmon";col2="darkturquoise"}
    else {col1="firebrick2";col2="dodgerblue3"}
  } else {
    col1=col[1];col2=col[2]
  }
  result=result[order(result[[1]],decreasing=TRUE),]
  
  max=max(result[,-4],na.rm=T)+0.1
  min=min(result[,-4],na.rm=T)
  
  x=log10(result[,1])*5/log10(max)
  x1=log10(result[,2])*5/log10(max)
  x2=log10(result[,3])*5/log10(max)
  opar<-par(no.readonly=TRUE)
  if(show.CI) par(mar=c(5,8,4,8))
  else par(mar=c(5,8,4,2))
  if(xlab=="") xlab=ifelse(colnames(result)[1]=="OR","Odds Ratios","Harzard Ratios")
  plot(result[,1],count:1,type="n",axes=FALSE,ylim=c(0.5,count+0.5),
       xlim=c(log10(min)*5/log10(max),5),xlab="",ylab=ylab,cex.axis=0.5,...)
  rect(par("usr")[1],par("usr")[3],par("usr")[2],par("usr")[4], border=FALSE,
       col = "white")
  
  abline(h=count:1,col="darkgray",lwd=1)
  abline(v=0,lty=2,col="darkgray",lwd=1)
  range=c(min,seq(1,max,1))
  range2=c(min+0.25,seq(1+0.25,max,1))
  #abline(v=log10(range2)*5/log10(max),col="white",lwd=1)
  abline(v=log10(range)*5/log10(max),col="darkgray",lwd=1)
  yscale=par("usr")[3:4]
  xscale=par("usr")[1:2]
  text(x=log10(range)*5/log10(max),y=par("usr")[3],range,pos=1,cex=0.7,xpd=TRUE)
  if(length(range2)<4) text(x=log10(range2)*5/log10(max),y=par("usr")[3],range2,pos=1,cex=0.6,xpd=TRUE)
  text(y=count:1,par("usr")[1],labels=rownames(result),pos=2,cex=0.9,xpd=TRUE)
  #text(x=5,y=par("usr")[3],xlab,pos=1,cex=1,xpd=TRUE)
  text(x=mean(range(xscale)),y=par("usr")[3],xlab,pos=1,cex=1,offset=2,xpd=TRUE)
  p=c()
  for(i in 1:count){
    if(is.nan(result[i,4])) p[i]=result[i,1]
    else if(result[i,4]<0.001) p[i]=paste(result[i,1],"***",sep="")
    else if(result[i,4]<0.01) p[i]=paste(result[i,1],"**",sep="")
    else if(result[i,4]<0.05) p[i]=paste(result[i,1],"*",sep="")
    else p[i]=result[i,1]
  }
  if(type<3){
    segments(x1,count:1,x2,count:1,lty=1,
             col=ifelse(x>0,col1,col2),lwd=lwd)
    if(show.OR) text(x,count:1,p,pos=3,cex=0.8)
    if(type==2) {
      segments(x1,(count:1)-0.1,x1,(count:1)+0.1,lty=1,
               col=ifelse(x>0,col1,col2),lwd=lwd)
      segments(x2,(count:1)-0.1,x2,(count:1)+0.1,lty=1,
               col=ifelse(x>0,col1,col2),lwd=lwd)
      points(x,count:1,col=ifelse(x>0,col1,col2),pch=pch,cex=cex)
    } else{
      points(x,count:1,col=ifelse(x>0,col1,col2),pch=pch,cex=cex)
      
    }
    
  }
  else{
    left=ifelse(x>0,0,x)
    right=ifelse(x>0,x,0)
    height=0.25
    rect(left,(count:1)+height,right,(count:1)-height,
         col=ifelse(x>0,col1,col2))
    segments(x1,count:1,x2,count:1,lty=1,lwd=lwd)
    if(show.OR) text(0,count:1,p,pos=ifelse(x>0,2,4),cex=0.8)
  }
  if(show.CI){
    text(y=count:1,par("usr")[2],
         paste(p," (",result[,2],"-",result[,3],")",sep=""),
         pos=4,xpd=TRUE)
    CI.title=paste(colnames(result)[1],"(95% C.I.)")
    text(y=max(yscale),par("usr")[2],CI.title,pos=4,xpd=TRUE)
  }
  #par(opar)
}

df = read.xlsx("Source Data File.xlsx",sheetIndex = 7)
colnames(df) <- df[1,]
df <- df[-1, ]

par(mfrow=c(3,1))
mod1 <- df[,2:5]
mod1 <- mod1 %>% mutate_if(is.character, as.numeric)
rownames(mod1) <- df[,1]
OR.sub(mod1,type=2,show.CI=TRUE,main="Local Health Area with Low SES",cex=1.1,col=c("coral1","cornflowerblue"))
mod2 <- df[,7:10]
mod2 <- mod2 %>% mutate_if(is.character, as.numeric)
rownames(mod1) <- df[,1]
OR.sub(mod2,type=2,show.CI=TRUE,main="Local Health Area with Medium SES",cex=1.1,col=c("coral1","cornflowerblue"))
mod3 <- df[,12:15]
mod3 <- mod3 %>% mutate_if(is.character, as.numeric)
rownames(mod3) <- df[,1]
OR.sub(mod3,type=2,show.CI=TRUE,main="Local Health Area with High SES",cex=1.1,col=c("coral1","cornflowerblue"))