require(INLA)
require(spdep)
require(ggplot2)
require(BayesX)
m<-read.bnd("C:/p/malawi.bnd")
sp<-bnd2sp(m)
mg<-bnd2gra(m)
nb<-gra2nb(mg)
adjmat <- as(nb2mat(nb, style = "B"), "Matrix")
d<-read.csv("C:/p/covid19data.csv")
idareaint<-d$idarea
idtimeint<-d$idtime

prioriid<-c(1,0.001)
priorbesag<-c(1,0.001)
priorrw2<-c(1,0.00005)

formula <- y ~ f(idarea, model = "bym", param=priorbesag, graph = adjmat) + 
  f(idarea2, idtime, model = "iid",param=prioriid) + idtime

mod <- inla(formula,family="poisson",data=d,E=population, 
                 control.predictor=list(compute=TRUE), 
                 control.compute=list(dic=TRUE,cpo=TRUE))
mod$dic$dic
l<- -mean(log(mod$cpo))
round(l, digits = 3)

# spatial temporal model with type 1 interaction
formula1<- y ~ f(idarea,model="bym",param=priorbesag,graph=adjmat) + f(idtime,model="rw2",param=priorrw2) + f(idtime2,model="iid",param=prioriid) + f(idint,model="iid",param=prioriid)

mod1 <- inla(formula1,family="poisson",data=d,E=population, control.predictor=list(compute=TRUE),control.compute=list(dic=TRUE,cpo=TRUE))

mod1$dic$dic
l1<- -mean(log(mod1$cpo))
round(l1, digits = 3)

formula12<- y ~f(idarea,model="bym", param=priorbesag,graph=adjmat) +
    f(idtime,model="rw1") + f(idtime2,model="iid",param=prioriid) + f(idint,model="iid",param=prioriid)

mod12 <- inla(formula12,family="poisson",data=d,E=population, 
                 control.predictor=list(compute=TRUE), 
                 control.compute=list(dic=TRUE,cpo=TRUE))

mod12$dic$dic
l12<- -mean(log(mod12$cpo))
round(l12, digits = 3)

formula2<- y ~1+f(idarea,model="bym",param=priorbesag, graph=adjmat,adjust.for.con.comp = FALSE) +
    f(idtime,model="rw2",param=priorrw2) + f(idtime2,model="iid",param=prioriid) +f(idtimeint,
  model = "iid",param=prioriid,group = idareaint, control.group = list(model = "besag",param=priorbesag,graph = adjmat))

mod2 <- inla(formula2,family="poisson",data=d,E=population,
control.predictor=list(compute=TRUE),control.compute=list(dic=TRUE,cpo=TRUE))
mod2$dic$dic
l2<- -mean(log(mod2$cpo))
round(l2, digits = 3)

formula22<- y ~1+f(idarea,model="bym",param=priorbesag, graph=adjmat,adjust.for.con.comp = FALSE) +
    f(idtime,model="rw1") + f(idtime2,model="iid",param=prioriid) +f(idtimeint,
  model = "iid",param=prioriid,group = idareaint, control.group = list(model = "besag",param=priorbesag,graph = adjmat))

mod22 <- inla(formula22,family="poisson",data=d,E=population, 
                 control.predictor=list(compute=TRUE), 
                 control.compute=list(dic=TRUE,cpo=TRUE))
mod22$dic$dic
l22<- -mean(log(mod22$cpo))
round(l22, digits = 3)

formula3<- y ~ f(idarea,model="bym", param=priorbesag,graph=adjmat) +
    f(idtime,model="rw2",param=priorrw2) + f(idtime2,model="iid",param=prioriid)+f(idareaint,
  model = "iid",param=prioriid, group = idtimeint, control.group = list(model = "rw2",param=priorrw2))  

mod3 <- inla(formula3,family="poisson",data=d,E=population, 
                 control.predictor=list(compute=TRUE), 
                 control.compute=list(dic=TRUE,cpo=TRUE))
mod3$dic$dic
l3<- -mean(log(mod3$cpo))
round(l3, digits = 3)

formula33<- y ~ f(idarea,model="bym", param=priorbesag,graph=adjmat) +
    f(idtime,model="rw1") + f(idtime2,model="iid",param=prioriid)+f(idareaint,
  model = "iid",param=prioriid, group = idtimeint, control.group = list(model = "rw1"))  



mod33 <- inla(formula33,family="poisson",data=d,E=population, 
                 control.predictor=list(compute=TRUE), 
                 control.compute=list(dic=TRUE,cpo=TRUE))
mod33$dic$dic
l33<- -mean(log(mod33$cpo))
round(l33, digits = 3)


formula4<- y ~ f(idarea,model="bym", param=priorbesag,graph=adjmat) +
    f(idtime,model="rw2",param=priorrw2) + f(idtime2,model="iid",param=prioriid)+f(idareaint,
  model = "besag",param=priorbesag, graph = adjmat,group = idtimeint, control.group = list(model = "rw2",param=priorrw2)
)

mod4 <- inla(formula4,family="poisson",data=d,E=population,control.fixed= list(prec.intercept = 1),
                 control.predictor=list(compute=TRUE), 
                 control.compute=list(dic=TRUE,cpo=TRUE))
mod4$dic$dic
l4<- -mean(log(mod4$cpo))
round(l4, digits = 3)

formula44<- y ~ f(idarea,model="bym", param=priorbesag,graph=adjmat) +
    f(idtime,model="rw1") + f(idtime2,model="iid",param=prioriid)+f(idareaint,
  model = "besag",param=priorbesag, graph = adjmat,group = idtimeint, control.group = list(model = "rw1"))

mod44 <- inla(formula44,family="poisson",data=d,E=population, 
                 control.predictor=list(compute=TRUE), 
                 control.compute=list(dic=TRUE,cpo=TRUE))
mod44$dic$dic
l44<- -mean(log(mod44$cpo))
round(l44, digits = 3)

# plotting cumlative cases over time
ggplot(d) +   geom_line(aes(x=idtime,y=cumulative, group=1, colour = 'Confirmed cases')) +    geom_line(aes(x=idtime,y=Dead, group = 2, colour = 'Dead')) +   ylab('Cumulative values')+xlab('Time (Week)')

12
# overall region risk by time
rate.est <- matrix(mod1$summary.fitted.values[,1], nrow = 31, byrow = FALSE)
rate.est <- as.data.frame(rate.est)
ST1<-rate.est[,1]
ST2<-rate.est[,2]
ST3<-rate.est[,3]
ST4<-rate.est[,4]
ST5<-rate.est[,5]
ST6<-rate.est[,6]
spatial<- data.frame(ID=seq(1,31),ST1,ST2,ST3,ST4,ST5,ST6)

par(mar=c(1,1,1,1))
par(mfrow=c(2,3))


plot1<-drawmap(data=spatial,map=m,regionvar="ID",
plotvar="ST1",swapcolors=T,cols="grey",density=36,drawnames=F,main="Week 1", cex.legend=1,cex.names=1)

plot2<-drawmap(data=spatial,map=m,regionvar="ID",
plotvar="ST2",swapcolors=T,cols="grey",density=36,drawnames=F,main="Week 2",cex.legend=1,cex.names=1)

plot3<-drawmap(data=spatial,map=m,regionvar="ID",
plotvar="ST3",swapcolors=T,cols="grey",density=36,drawnames=F,main="Week 3",cex.legend=1,cex.names=1)

plot4<-drawmap(data=spatial,map=m,regionvar="ID",
plotvar="ST4",swapcolors=T,cols="grey",density=36,drawnames=F,main="Week 4",cex.legend=1,cex.names=1)

plot5<-drawmap(data=spatial,map=m,regionvar="ID",
plotvar="ST5",swapcolors=T,cols="grey",density=36,drawnames=F,main="Week 5",cex.legend=1,cex.names=1)

plot6<-drawmap(data=spatial,map=m,regionvar="ID",
plotvar="ST6",swapcolors=T,cols="grey",density=36,drawnames=F,main="Week 6",cex.legend=1,cex.names=1)

dev.off()

#Overall Risk trend by region
theme_set(theme_gray(base_size = 12)) # font size
rate.est <- matrix(mod1$summary.fitted.values[,1], nrow = 186, byrow = TRUE)
rate.est <- as.data.frame(rate.est)
d<- data.frame(ID=seq(1,186),d,rate.est)
d$Time<-d$idtime
d$District<-d$idarea
d$Time <- factor(d$Time,labels = c("Week 1", "Week 2", "Week 3", "Week 4", "Week 5","Week 6"))
d$District <- factor(d$District,labels = c("DZ", "DA", "KU", "LL", "MC","KK",
"NU","NS","SA","LL CITY","CP","KA","MZ","NB","RU","LA","MZ CITY","BT","CK","CZ","MHG","MH","MJ","MN","NE","TO","ZA","PE","BK","BT CITY","ZA CITY"))
p2<-ggplot(d, aes(x = Time, y = V1, group = District)) + geom_line(aes(colour = District))+ylab("Estimated Relative Risk")
p2 


# Interaction effect
STint <- matrix(exp(mod1$summary.random$idint[,2]), nrow = 31, byrow = FALSE)
STint <- as.data.frame(STint)
ST1<-STint[,1]
ST2<-STint[,2]
ST3<-STint[,3]
ST4<-STint[,4]
ST5<-STint[,5]
ST6<-STint[,6]
spatial<- data.frame(ID=seq(1,31),ST1,ST2,ST3,ST4,ST5,ST6)

par(mar=c(1,1,1,1))
par(mfrow=c(2,3))
plot1<-drawmap(data=spatial,map=m,regionvar="ID",
plotvar="ST1",swapcolors=T,cols="grey",density=36,drawnames=F,main="Week 1", cex.legend=1,cex.names=1)

plot2<-drawmap(data=spatial,map=m,regionvar="ID",
plotvar="ST2",swapcolors=T,cols="grey",density=36,drawnames=F,main="Week 2",cex.legend=1,cex.names=1)

plot3<-drawmap(data=spatial,map=m,regionvar="ID",
plotvar="ST3",swapcolors=T,cols="grey",density=36,drawnames=F,main="Week 3",cex.legend=1,cex.names=1)

plot4<-drawmap(data=spatial,map=m,regionvar="ID",
plotvar="ST4",swapcolors=T,cols="grey",density=36,drawnames=F,main="Week 4",cex.legend=1,cex.names=1)

plot5<-drawmap(data=spatial,map=m,regionvar="ID",
plotvar="ST5",swapcolors=T,cols="grey",density=36,drawnames=F,main="Week 5",cex.legend=1,cex.names=1)

plot6<-drawmap(data=spatial,map=m,regionvar="ID",
plotvar="ST6",swapcolors=T,cols="grey",density=36,drawnames=F,main="Week 6",cex.legend=1,cex.names=1)