# required packages pacman::p_load(dplyr,RColorBrewer,psych,ggplot2,reshape2,deSolve,nleqslv,cowplot) ABC_C<-read.csv("/Users/tvino/Documents/PhD/Analysis/ABC_Household_Analysis/PeerJ/FINAL/Codes/ABC_household_raw_de-identified.csv",header=T) # Figure 3 - Household size distributions #To create the Bar chart for Household size with Age Category classification with A,C,B proportion #Creating the proportion ABC_C["Prop_A"]<-ABC_C$all_adult/ABC_C$overall_no ABC_C["Prop_C"]<-ABC_C$all_child/ABC_C$overall_no ABC_C["Prop_B"]<-ABC_C$all_baby/ABC_C$overall_no #View(ABC_C) ABC_C_Selected<-subset(ABC_C,select = c("overall_no","rem_urb","Prop_A","Prop_C","Prop_B")) #View(ABC_C_Selected) Remote_Sel<-ABC_C_Selected[which(ABC_C$rem_urb=="ABC remote"),] Urban_Sel<-ABC_C_Selected[which(ABC_C$rem_urb=="ABC urban"),] #To melt, removing the urban/remote variables Remote_Sel2<-subset(Remote_Sel,select = c("overall_no","Prop_A","Prop_C","Prop_B")) Urban_Sel2<-subset(Urban_Sel,select = c("overall_no","Prop_A","Prop_C","Prop_B")) #View(Remote_Sel2) Remote_Sel2.m <- melt(Remote_Sel2,id.vars = "overall_no") Urban_Sel2.m <- melt(Urban_Sel2,id.vars = "overall_no") #View(Remote_Sel2.m) #To add the urban/remote in same axis Remote_Sel2.m["rem_urb"]<-"Remote" Urban_Sel2.m["rem_urb"]<-"Urban" Merged_Rem_Urb <- rbind(Remote_Sel2.m, Urban_Sel2.m) #View(Merged_Rem_Urb) # Renaming the colum names and variable names Merged_Rem_Urb<-rename(Merged_Rem_Urb, Household_Size=overall_no,Proportion_of_Age_Category=variable, Count=value) Merged_Rem_Urb$Proportion_of_Age_Category <- as.character(Merged_Rem_Urb$Proportion_of_Age_Category) Merged_Rem_Urb$Proportion_of_Age_Category[Merged_Rem_Urb$Proportion_of_Age_Category == "Prop_A"] <- "Adult" Merged_Rem_Urb$Proportion_of_Age_Category[Merged_Rem_Urb$Proportion_of_Age_Category == "Prop_C"] <- "Child" Merged_Rem_Urb$Proportion_of_Age_Category[Merged_Rem_Urb$Proportion_of_Age_Category == "Prop_B"] <- "Baby" Merged_Rem_Urb$Proportion_of_Age_Category<-as.factor(Merged_Rem_Urb$Proportion_of_Age_Category) #View(Merged_Rem_Urb) Merged_Rem_Urb$Proportion_of_Age_Category<- factor(Merged_Rem_Urb$Proportion_of_Age_Category, levels = c("Baby", "Child", "Adult")) #renaming column values to match variables written in paper Merged_Rem_Urb$Proportion_of_Age_Category<-as.character(Merged_Rem_Urb$Proportion_of_Age_Category) Merged_Rem_Urb[Merged_Rem_Urb=="Child"]<-"School aged Child" Merged_Rem_Urb[Merged_Rem_Urb=="Baby"]<-"Pre-school aged Child" Merged_Rem_Urb$Proportion_of_Age_Category<- factor(Merged_Rem_Urb$Proportion_of_Age_Category, levels = c("Pre-school aged Child", "School aged Child", "Adult")) #Seperate graphs to label as A, B and then combine Remote<-Merged_Rem_Urb[Merged_Rem_Urb$rem_urb == 'Remote',] Urban<-Merged_Rem_Urb[Merged_Rem_Urb$rem_urb == 'Urban',] Remote_p<- ggplot(Remote, aes(x = Household_Size, y = Count, fill=Proportion_of_Age_Category)) + geom_bar(stat='identity') + #facet_grid(. ~rem_urb ) + scale_fill_brewer(palette = "Set1") + theme_bw(base_size = 17) + theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank()) + labs(list( x = "Household Size", y = "Count",fill= "Age Category")) +xlim(0,25)+ylim(0,27) Urban_p<- ggplot(Urban, aes(x = Household_Size, y = Count, fill=Proportion_of_Age_Category)) + geom_bar(stat='identity') + #facet_grid(. ~rem_urb ) + scale_fill_brewer(palette = "Set1") + theme_bw(base_size = 17) + theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank()) + labs(list( x = "Household Size", y = "Count",fill= "Age Category")) + xlim(0,25)+ylim(0,27) Urban_p #To combine (using package cowplot) plot_grid(Remote_p + theme(legend.position = 'none'), Urban_p + theme(legend.position = 'none'), get_legend(Remote_p), labels = c('A', 'B'), label_size=25, nrow=1, rel_widths = c(1, 1, 0.45)) ggsave(filename="/Users/tvino/Documents/PhD/Analysis/ABC_Household_Analysis/PeerJ/Final/Revised_Main/Figure_03.png", width = 38, height = 20, units = "cm") # Figure 4 - Room occupancy rates #renaming the values in Rem_Urb and column names ABC_C$rem_urb<-as.character(ABC_C$rem_urb) ABC_C[ABC_C=="ABC remote"]<-"Remote" ABC_C[ABC_C=="ABC urban"]<-"Urban" #Reducing some shires as Other levels(ABC_C$shire)[table(ABC_C$shire) < 10] <- 'Other' #Apply sorting ABC_C$shire <- factor(ABC_C$shire, levels(ABC_C$shire)[order(summary(ABC_C$shire),decreasing = TRUE)]) #Seperate graphs to label as A, B and then combine Remote_Occupancy<-ABC_C[ABC_C$rem_urb == 'Remote',] Urban_Occupancy<-ABC_C[ABC_C$rem_urb == 'Urban',] # Converting No.of.Rooms.in.a.house.hold to factor ABC_C$No.of.Rooms.in.a.house.hold<-as.factor(ABC_C$No.of.Rooms.in.a.house.hold) Remote_Occupancy$No.of.Rooms.in.a.house.hold<-as.factor(Remote_Occupancy$No.of.Rooms.in.a.house.hold) Urban_Occupancy$No.of.Rooms.in.a.house.hold<-as.factor(Urban_Occupancy$No.of.Rooms.in.a.house.hold) Remote_O_p<-ggplot(Remote_Occupancy, aes(No.of.Rooms.in.a.house.hold,Average.no.of.ppl.in.a.room)) + geom_violin(fill="grey87",colour="white",scale="count") + theme_bw(base_size = 17) + theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank())+ geom_jitter(aes(colour=shire),width = 0.25) + #facet_grid(. ~rem_urb ) + scale_colour_manual(values = c("#e41a1c","#377eb8", "#4daf4a","#984ea3")) + labs(list( x = "No of Rooms", y = "Occupancy Rate")) Urban_O_p<- ggplot(Urban_Occupancy, aes(No.of.Rooms.in.a.house.hold,Average.no.of.ppl.in.a.room)) + geom_violin(fill="grey87",colour="white",scale="count") + theme_bw(base_size = 17) + theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank())+ geom_jitter(aes(colour=shire),width = 0.25) + #facet_grid(. ~rem_urb ) + scale_colour_manual(values = c("#e41a1c","#984ea3","#ff7f00","#e78ac3")) + labs(list( x = "No of Rooms", y = "Occupancy Rate")) + scale_x_discrete(name="No of Rooms", limits = c(1,2,3,4,5)) Combined<-ggplot(ABC_C, aes(No.of.Rooms.in.a.house.hold,Average.no.of.ppl.in.a.room)) + geom_violin(fill="grey87",colour="white",scale="count") + theme_bw(base_size = 17) + theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank())+ geom_jitter(aes(colour=shire),width = 0.25) + #facet_grid(. ~rem_urb ) + scale_colour_manual(values = c("#e41a1c","#377eb8", "#4daf4a","#984ea3","#ff7f00","#e78ac3")) + labs(list( x = "No of Rooms", y = "Occupancy Rate"))+ylim(0,6) #To combine (using package cowplot) plot_grid(Remote_O_p + theme(legend.position = 'none'), Urban_O_p + theme(legend.position = 'none'), get_legend(Combined), labels = c('A', 'B'), label_size=25, nrow=1, rel_widths = c(1, 1, 0.28)) ggsave(filename="/Users/tvino/Documents/PhD/Analysis/ABC_Household_Analysis/PeerJ/Final/Revised_Main/Figure_04.png", width = 35, height = 20, units = "cm") # Figure 5 - Household contact matrices # Figure 6 - Household contact matrices – Indigenous and non-Indigenous # Figure 7 - Proportion of infected population - Indigenous and non-Indigenous # Calculating Community mixing Avg_contact<-13.35 Avg_contact_Community<-Avg_contact*0.78 ACC<-round(Avg_contact_Community, digits = 0) #Other parameters phi<-0.67 #(Rate of change from Exposure to infected (1.5 days)) gamma<-0.67 # (Rate of change from Infected to Recovered (1.5 days)) q1 <-0.14 #probability of transmission at a contact at HH level q2 <-0.5 *q1 #probability of transmission at a contact at community level Ea<-1e-6 # Number of Exposed adults, In this model, I'm assuming only one adult gets exposed Ec<-1e-6 # Number of Exposed children Eb<-1e-6 # Number of Exposed babies # SEIR Model- Indigenous-Remote #As per ABS Census data from Table builder, Under persons and relationships Age in single years, Melbourne, IREA locations of Northern Territory and section of State/urban centres and localities geography -NT-SOS level classification #Age proportion of Indigenous people IRA<-0.77 IRC<-0.14 IRB<-0.09 #From HH contact matrix C1bb_ir<-0.5 C1bc_ir<-1.7 C1ba_ir<-4.2 C1cb_ir<-1.7 C1cc_ir<-3.4 C1ca_ir<-11.4 C1ab_ir<-4.2 C1ac_ir<-11.4 C1aa_ir<-11.9 S = matrix(c("Sa","Sc","Sb"), nrow = 3, ncol=1) E = matrix(c("Ea","Ec","Eb"), nrow = 3, ncol=1) I = matrix(c("Ia","Ic","Ib"), nrow = 3, ncol=1) R = matrix(c("Ra","Rc","Rb"), nrow = 3, ncol=1) C1h = matrix (c("C1bb","C1bc","C1ba","C1cb","C1cc","C1ca","C1ab","C1ac","C1aa"), nrow = 3, ncol=3) C2c = matrix (c("C2bb","C2bc","C2ba","C2cb","C2cc","C2ca","C2ab","C2ac","C2aa"), nrow = 3, ncol=3) SEIR_MODEL <- function(time, state, parameters) { with(as.list(c(state, parameters)), { S = matrix(state[1:3], nrow = 3, ncol=1) E = matrix(state[4:6], nrow = 3, ncol=1) I = matrix(state[7:9], nrow = 3, ncol=1) R = matrix(state[10:12], nrow = 3, ncol=1) C1h = matrix (c(C1bb,C1bc,C1ba,C1cb,C1cc,C1ca,C1ab,C1ac,C1aa), nrow=3, ncol=3) C2c = matrix (c(C2bb,C2bc,C2ba,C2cb,C2cc,C2ca,C2ab,C2ac,C2aa), nrow=3, ncol=3) dS <- (-((q1*C1h) %+% (q2*C2c) %*% I)) * as.vector(S) dE <- ((q1*C1h) %+% (q2*C2c) %*% I) * as.vector(S) - phi * E dI <- phi * E - gamma * I dR <- gamma * I return(list(c(dS, dE, dI, dR))) }) } init <- c(Sa=IRA-Ea, Sc=IRC, Sb=IRB, Ea=Ea, Ec=0, Eb=0, Ia=0, Ic=0, Ib=0, Ra=0.0, Rc=0.0, Rb=0.0) parameters <- c(C1bb=C1bb_ir, C1bc=C1bc_ir, C1ba=C1ba_ir, C1cb=C1cb_ir, C1cc=C1cc_ir, C1ca=C1ca_ir, C1ab=C1ab_ir, C1ac=C1ac_ir, C1aa=C1aa_ir, C2bb=IRB*ACC, C2bc=IRC*ACC, C2ba=IRA*ACC, C2cb=IRB*ACC, C2cc=IRC*ACC, C2ca=IRA*ACC, C2ab=IRB*ACC, C2ac=IRC*ACC, C2aa=IRA*ACC, phi=phi,gamma=gamma, q1=q1, q2=q2) times <- seq(0, 250, by = 1) out6 <- ode(y=init, times=times, func=SEIR_MODEL, parms=parameters) out7 <- as.data.frame(out6) out7$time <- NULL #Combining the age proportion in the output out8<-out7%>%mutate(S=Sa+Sb+Sc)%>%mutate(E=Ea+Eb+Ec)%>%mutate(I=Ia+Ib+Ic)%>%mutate(R=Ra+Rb+Rc) out_IR<-subset(out8, select=c("S","E","I","R")) # FINDING R0 C1h = matrix (c(C1bb_ir, C1bc_ir, C1ba_ir, C1cb_ir, C1cc_ir, C1ca_ir,C1ab_ir, C1ac_ir, C1aa_ir), nrow=3, ncol=3) C2c = matrix (c(IRB*ACC, IRC*ACC, IRA*ACC, IRB*ACC,IRC*ACC,IRA*ACC, IRB*ACC, IRC*ACC, IRA*ACC), nrow=3, ncol=3) C_COM<-(q1*C1h) + (q2*C2c) C_COM #Using next generation matrix dslnex_IR <- function(x) { y <- numeric(3) y[1] <- R0_IR[1]*x[1]+ R0_IR[4]*x[2]+R0_IR[7]*(1-x[1]-x[2])-x[3]*x[1] y[2] <- R0_IR[2]*x[1]+ R0_IR[5]*x[2]+R0_IR[8]*(1-x[1]-x[2])-x[3]*x[2] y[3] <- R0_IR[3]*x[1]+ R0_IR[6]*x[2]+R0_IR[9]*(1-x[1]-x[2])-x[3]*(1-x[1]-x[2]) y } xstart <- c(0.1,0.2,1) fstart <- dslnex(xstart) xstart fstart nleqslv(xstart, dslnex, control=list(btol=.01)) # SEIR Model- Indigenous-Urban #As per ABS Census data from Table builder #Age proportion of Indigenous people IUA<-0.72 IUC<-0.17 IUB<-0.11 #From HH contact matrix C1bb_iu<-0.3 C1bc_iu<-0.8 C1ba_iu<-2 C1cb_iu<-0.8 C1cc_iu<-0.5 C1ca_iu<-3 C1ab_iu<-2 C1ac_iu<-3 C1aa_iu<-4.8 S = matrix(c("Sa","Sc","Sb"), nrow = 3, ncol=1) E = matrix(c("Ea","Ec","Eb"), nrow = 3, ncol=1) I = matrix(c("Ia","Ic","Ib"), nrow = 3, ncol=1) R = matrix(c("Ra","Rc","Rb"), nrow = 3, ncol=1) C1h = matrix (c("C1bb","C1bc","C1ba","C1cb","C1cc","C1ca","C1ab","C1ac","C1aa"), nrow = 3, ncol=3) C2c = matrix (c("C2bb","C2bc","C2ba","C2cb","C2cc","C2ca","C2ab","C2ac","C2aa"), nrow = 3, ncol=3) SEIR_MODEL <- function(time, state, parameters) { with(as.list(c(state, parameters)), { S = matrix(state[1:3], nrow = 3, ncol=1) E = matrix(state[4:6], nrow = 3, ncol=1) I = matrix(state[7:9], nrow = 3, ncol=1) R = matrix(state[10:12], nrow = 3, ncol=1) C1h = matrix (c(C1bb,C1bc,C1ba,C1cb,C1cc,C1ca,C1ab,C1ac,C1aa), nrow=3, ncol=3) C2c = matrix (c(C2bb,C2bc,C2ba,C2cb,C2cc,C2ca,C2ab,C2ac,C2aa), nrow=3, ncol=3) dS <- (-((q1*C1h) %+% (q2*C2c) %*% I)) *as.vector(S) dE <- ((q1*C1h) %+% (q2*C2c) %*% I) *as.vector(S) - phi * E dI <- phi * E - gamma * I dR <- gamma * I return(list(c(dS, dE, dI, dR))) }) } init <- c(Sa=IUA-Ea, Sc=IUC, Sb=IUB, Ea=Ea, Ec=0, Eb=0, Ia=0, Ic=0, Ib=0, Ra=0.0, Rc=0.0, Rb=0.0) parameters <- c(C1bb=C1bb_iu, C1bc=C1bc_iu, C1ba=C1ba_iu, C1cb=C1cb_iu, C1cc=C1cc_iu, C1ca=C1ca_iu, C1ab=C1ab_iu, C1ac=C1ac_iu, C1aa=C1aa_iu, C2bb=IUB*ACC, C2bc=IUC*ACC, C2ba=IUA*ACC, C2cb=IUB*ACC, C2cc=IUC*ACC, C2ca=IUA*ACC, C2ab=IUB*ACC, C2ac=IUC*ACC, C2aa=IUA*ACC, phi=phi,gamma=gamma, q1=q1, q2=q2) times <- seq(0, 250, by = 1) out6 <- ode(y=init, times=times, func=SEIR_MODEL, parms=parameters) out7 <- as.data.frame(out6) out7$time <- NULL #Combining the age proportion in the output out8<-out7%>%mutate(S=Sa+Sb+Sc)%>%mutate(E=Ea+Eb+Ec)%>%mutate(I=Ia+Ib+Ic)%>%mutate(R=Ra+Rb+Rc) out_IU<-subset(out8, select=c("S","E","I","R")) # FINDING R0 C1h = matrix (c(C1bb_iu, C1bc_iu, C1ba_iu, C1cb_iu, C1cc_iu, C1ca_iu,C1ab_iu, C1ac_iu, C1aa_iu), nrow=3, ncol=3) C2c = matrix (c(IUB*ACC, IUC*ACC, IUA*ACC, IUB*ACC,IUC*ACC,IUA*ACC, IUB*ACC, IUC*ACC, IUA*ACC), nrow=3, ncol=3) C_COM<-(q1*C1h) + (q2*C2c) R0_IU<-C_COM/gamma R0_IU<-round(R0_IU,2) R0_IU dslnex_IU <- function(x) { y <- numeric(3) y[1] <- R0_IU[1]*x[1]+ R0_IU[4]*x[2]+R0_IU[7]*(1-x[1]-x[2])-x[3]*x[1] y[2] <- R0_IU[2]*x[1]+ R0_IU[5]*x[2]+R0_IU[8]*(1-x[1]-x[2])-x[3]*x[2] y[3] <- R0_IU[3]*x[1]+ R0_IU[6]*x[2]+R0_IU[9]*(1-x[1]-x[2])-x[3]*(1-x[1]-x[2]) y } xstart <- c(0.1,0.2,1) fstart <- dslnex(xstart) xstart fstart nleqslv(xstart, dslnex_IU, control=list(btol=.01)) # SEIR Model- Indigenous-Whole IA<-0.75 IC<-0.15 IB<-0.10 #From HH contact matrix C1bb_i<-0.5 C1bc_i<-1.5 C1ba_i<-3.8 C1cb_i<-1.5 C1cc_i<-2.8 C1ca_i<-10 C1ab_i<-3.8 C1ac_i<-10 C1aa_i<-10.6 S = matrix(c("Sa","Sc","Sb"), nrow = 3, ncol=1) E = matrix(c("Ea","Ec","Eb"), nrow = 3, ncol=1) I = matrix(c("Ia","Ic","Ib"), nrow = 3, ncol=1) R = matrix(c("Ra","Rc","Rb"), nrow = 3, ncol=1) C1h = matrix (c("C1bb","C1bc","C1ba","C1cb","C1cc","C1ca","C1ab","C1ac","C1aa"), nrow = 3, ncol=3) # Contact matrix of household (Nic created) C2c = matrix (c("C2bb","C2bc","C2ba","C2cb","C2cc","C2ca","C2ab","C2ac","C2aa"), nrow = 3, ncol=3) # Contact matrix of community SEIR_MODEL <- function(time, state, parameters) { with(as.list(c(state, parameters)), { S = matrix(state[1:3], nrow = 3, ncol=1) E = matrix(state[4:6], nrow = 3, ncol=1) I = matrix(state[7:9], nrow = 3, ncol=1) R = matrix(state[10:12], nrow = 3, ncol=1) C1h = matrix (c(C1bb,C1bc,C1ba,C1cb,C1cc,C1ca,C1ab,C1ac,C1aa), nrow=3, ncol=3) C2c = matrix (c(C2bb,C2bc,C2ba,C2cb,C2cc,C2ca,C2ab,C2ac,C2aa), nrow=3, ncol=3) dS <- (-((q1*C1h) %+% (q2*C2c) %*% I)) *as.vector(S) dE <- ((q1*C1h) %+% (q2*C2c) %*% I) *as.vector(S) - phi * E dI <- phi * E - gamma * I dR <- gamma * I return(list(c(dS, dE, dI, dR))) }) } init <- c(Sa=IA-Ea, Sc=IC, Sb=IB, Ea=Ea, Ec=0, Eb=0, Ia=0, Ic=0, Ib=0, Ra=0.0, Rc=0.0, Rb=0.0) parameters <- c(C1bb=C1bb_i, C1bc=C1bc_i, C1ba=C1ba_i, C1cb=C1cb_i, C1cc=C1cc_i, C1ca=C1ca_i, C1ab=C1ab_i, C1ac=C1ac_i, C1aa=C1aa_i, C2bb=IB*ACC, C2bc=IC*ACC, C2ba=IA*ACC, C2cb=IB*ACC, C2cc=IC*ACC, C2ca=IA*ACC, C2ab=IB*ACC, C2ac=IC*ACC, C2aa=IA*ACC, phi=phi,gamma=gamma, q1=q1, q2=q2) times <- seq(0, 250, by = 1) out6 <- ode(y=init, times=times, func=SEIR_MODEL, parms=parameters) out7 <- as.data.frame(out6) out7$time <- NULL #Combining the age proportion in the output out8<-out7%>%mutate(S=Sa+Sb+Sc)%>%mutate(E=Ea+Eb+Ec)%>%mutate(I=Ia+Ib+Ic)%>%mutate(R=Ra+Rb+Rc) out_I<-subset(out8, select=c("S","E","I","R")) # SEIR Model- Non-Indigenous Urban NIA<-0.80 NIC<-0.12 NIB<-0.08 #From HH contact matrix C1bb_ni<-0.02 C1bc_ni<-0.06 C1ba_ni<-0.16 C1cb_ni<-0.06 C1cc_ni<-0.23 C1ca_ni<-0.97 C1ab_ni<-0.16 C1ac_ni<-0.97 C1aa_ni<-1.62 S = matrix(c("Sa","Sc","Sb"), nrow = 3, ncol=1) E = matrix(c("Ea","Ec","Eb"), nrow = 3, ncol=1) I = matrix(c("Ia","Ic","Ib"), nrow = 3, ncol=1) R = matrix(c("Ra","Rc","Rb"), nrow = 3, ncol=1) C1h = matrix (c("C1bb","C1bc","C1ba","C1cb","C1cc","C1ca","C1ab","C1ac","C1aa"), nrow = 3, ncol=3) C2c = matrix (c("C2bb","C2bc","C2ba","C2cb","C2cc","C2ca","C2ab","C2ac","C2aa"), nrow = 3, ncol=3) SEIR_MODEL <- function(time, state, parameters) { with(as.list(c(state, parameters)), { S = matrix(state[1:3], nrow = 3, ncol=1) E = matrix(state[4:6], nrow = 3, ncol=1) I = matrix(state[7:9], nrow = 3, ncol=1) R = matrix(state[10:12], nrow = 3, ncol=1) C1h = matrix (c(C1bb,C1bc,C1ba,C1cb,C1cc,C1ca,C1ab,C1ac,C1aa), nrow=3, ncol=3) C2c = matrix (c(C2bb,C2bc,C2ba,C2cb,C2cc,C2ca,C2ab,C2ac,C2aa), nrow=3, ncol=3) dS <- (-((q1*C1h) %+% (q2*C2c) %*% I)) *as.vector(S) dE <- ((q1*C1h) %+% (q2*C2c) %*% I) *as.vector(S) - phi * E dI <- phi * E - gamma * I dR <- gamma * I return(list(c(dS, dE, dI, dR))) }) } init <- c(Sa=NIA-Ea, Sc=NIC, Sb=NIB, Ea=Ea, Ec=0, Eb=0, Ia=0, Ic=0, Ib=0, Ra=0.0, Rc=0.0, Rb=0.0) parameters <- c(C1bb=C1bb_ni, C1bc=C1bc_ni, C1ba=C1ba_ni, C1cb=C1cb_ni, C1cc=C1cc_ni, C1ca=C1ca_ni, C1ab=C1ab_ni, C1ac=C1ac_ni, C1aa=C1aa_ni, C2bb=NIB*ACC, C2bc=NIC*ACC, C2ba=NIA*ACC, C2cb=NIB*ACC, C2cc=NIC*ACC, C2ca=NIA*ACC, C2ab=NIB*ACC, C2ac=NIC*ACC, C2aa=NIA*ACC, phi=phi,gamma=gamma, q1=q1, q2=q2) times <- seq(0, 250, by = 1) out6 <- ode(y=init, times=times, func=SEIR_MODEL, parms=parameters) out7 <- as.data.frame(out6) out7$time <- NULL #Combining the age proportion in the output out8<-out7%>%mutate(S=Sa+Sb+Sc)%>%mutate(E=Ea+Eb+Ec)%>%mutate(I=Ia+Ib+Ic)%>%mutate(R=Ra+Rb+Rc) out_NI<-subset(out8, select=c("S","E","I","R")) # Combining the "I" of IR,IU and NI, renaming and plotting out_IR_I<-subset(out_IR, select=c("I")) out_IR_I<-rename(out_IR_I, "Indigenous-Remote"=I) out_IU_I<-subset(out_IU, select=c("I")) out_IU_I<-rename(out_IU_I, "Indigenous-Urban"=I) out_NI_I<-subset(out_NI, select=c("I")) out_NI_I<-rename(out_NI_I, "Non-Indigenous"=I) out_all <- cbind(out_IR_I,out_IU_I,out_NI_I) matplot(x = times, y = out_all, type = "l", cex.lab = 1.5, cex.axis = 1.5, xlab = "Time (days)", ylab = "Proportion Infected", main = "", lwd = 2, lty = c(1,2,3), bty = "l", col = c(2)) legend("topright", col = c(2), legend=c("Indigenous-Remote","Indigenous-Urban","Non-Indigenous"), lwd=2,lty =c(1,2,3),bty= "o",box.lty=3,box.lwd=1, box.col="white",cex = 1.5) # plotting with ggplot times_f_07 <- as.data.frame(times) out_all<-as.data.frame(out_all) out_all_time<-cbind(times_f_07,out_all) # Melt the data to long format out_all_time_melt <- melt(out_all_time, id="times") #rename the titles out_all_time_melt<-rename(out_all_time_melt, "Population"=variable, "Infected_proportion"=value) ggplot(out_all_time_melt, aes(x=times, y=Infected_proportion, lty=Population)) + geom_line(col="red") + scale_linetype_manual(values=c("solid","dashed", "dotdash"))+ labs(list( x = "Time (days)", y = "Proportion Infected"))+ylim(0,0.15) + theme_bw(base_size = 17)+ theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank()) ggsave(filename="/Users/tvino/Documents/PhD/Analysis/ABC_Household_Analysis/PeerJ/Final/Revised_Main/Figure_07.png", width = 35, height = 20, units = "cm") # FINDING R0 C1h = matrix (c(C1bb_ni, C1bc_ni, C1ba_ni, C1cb_ni, C1cc_ni, C1ca_ni,C1ab_ni, C1ac_ni, C1aa_ni), nrow=3, ncol=3) C2c = matrix (c(NIB*ACC, NIC*ACC, NIA*ACC, NIB*ACC, NIC*ACC, NIA*ACC, NIB*ACC, NIC*ACC, NIA*ACC), nrow=3, ncol=3) C_COM<-(q1*C1h)+(q2*C2c) C_COM R0_NI<-C_COM/gamma R0_NI<-round(R0_NI,2) R0_NI dslnex_NI <- function(x) { y <- numeric(3) y[1] <- R0_NI[1]*x[1]+ R0_NI[4]*x[2]+R0_NI[7]*(1-x[1]-x[2])-x[3]*x[1] y[2] <- R0_NI[2]*x[1]+ R0_NI[5]*x[2]+R0_NI[8]*(1-x[1]-x[2])-x[3]*x[2] y[3] <- R0_NI[3]*x[1]+ R0_NI[6]*x[2]+R0_NI[9]*(1-x[1]-x[2])-x[3]*(1-x[1]-x[2]) y } xstart <- c(0.1,0.2,1) fstart <- dslnex(xstart) xstart fstart nleqslv(xstart, dslnex_NI, control=list(btol=.01)) # Supp Figure 1 - Household Size distribution for Census data Census<-read.csv("/Users/tvino/Documents/PhD/Analysis/ABC_Household_Analysis/PeerJ/Census_2011_Comparison.csv",header=T) View(Census) ggplot(Census, aes(x = HH_Size, y = Percentage)) + geom_bar(stat='identity',fill="#377eb8")+ facet_grid(. ~Rem_urb )+ scale_colour_manual(values = c("red"))+ theme_bw(base_size = 17)+ theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank())+ labs(list( x = "Household Size", y = "Percentage")) # Supp Figure 2 - Contact matrices for each shire