# required packages pacman::p_load(psych,dplyr,deSolve) # 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")) # 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")) # 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 and plotting out_IR_I<-subset(out_IR, select=c("I")) out_IU_I<-subset(out_IU, select=c("I")) out_NI_I<-subset(out_NI, select=c("I")) out_all <- cbind(out_IR_I,out_IU_I,out_NI_I) matplot(x = times, y = out_all, type = "l", xlab = "Time (days)", ylab = "Proportion Infected", main = "", lwd = 1, lty = c(1,2,3), bty = "l", col = c(2)) legend("topright", col = c(2), legend=c("Indigenous-Remote","Indigenous-Urban","Non-Indigenous"), lwd=1,lty =c(1,2,3),bty= "o",box.lty=3,box.lwd=1, box.col="white")