# Version 15 March 2022

library(covid19swiss)
dates<-unique(covid19swiss$date)[-1]

# prepare case data 
sel<-covid19swiss$data_type=="cases_total"
totcases<-data.frame(t(tapply(covid19swiss$value[sel],INDEX=list(covid19swiss$location[sel],covid19swiss$date[sel]), FUN = function(x){x})))
totcases[1,]<-0
totcases <- apply(totcases,2,function(x){for (i in 1:length(x)){if (is.na(x[i])) {x[i]<-x[i-1]}};return(t(x))})
cases<-as.data.frame(apply(totcases,2,diff))

# prepare death data 
sel<-covid19swiss$data_type=="deaths_total"
totdeaths<-data.frame(t(tapply(covid19swiss$value[sel],INDEX=list(covid19swiss$location[sel],covid19swiss$date[sel]), FUN = function(x){x})))
totdeaths[1,]<-0
totdeaths <- apply(totdeaths,2,function(x){for (i in 1:length(x)){if (is.na(x[i])) {x[i]<-x[i-1]}};return(t(x))})
deaths<-as.data.frame(apply(totdeaths,2,diff))

###########################################################
# 1. How did pandemic develop in Ticino?

# plot raw data
plot(dates,cases$TI, type="l")
title("Covid cases TI")
lines(smooth.spline(dates,y=cases$TI),col="blue")

plot(dates,deaths$TI, type="l")
title("Covid deaths TI")
lines(smooth.spline(dates,y=deaths$TI),col="blue")

# plot SQRT transformed data
plot(dates,sqrt(cases$TI),type="l",xlab="",ylab="",main = "Sqrt Cases TI")
plot(dates,sqrt(deaths$TI),type="l",xlab="",ylab="",main = "Sqrt Deaths TI")

# put transformed data in matrix
lcases<-as.data.frame(sqrt(as.matrix(cases)))
colnames(lcases)<-colnames(cases)


###########################################################
# 2. How did cantons infect each other?

# plot relationship between cantons
plot(lcases$VS,lcases$TI,xlab="number of cases in VS",
      ylab="number of cases in TI")

plot(lcases$LU,lcases$TI,xlab="number of cases in LU",
     ylab="number of cases in TI")

# Fit linear regressions
fit1<-lm(lcases$TI~lcases$VS)
summary(fit1)
plot(lcases$VS,lcases$TI,xlab="number of cases in VS",
     ylab="number of cases in TI")
abline(coef(fit1),col="blue")

fit2<-lm(lcases$TI~lcases$LU)
summary(fit2)
plot(lcases$LU,lcases$TI,xlab="number of cases in LU",
     ylab="number of cases in TI")
abline(coef(fit2),col="blue")


# Multiple regression
fit<-lm(lcases$TI~lcases$VS+lcases$LU)
summary(fit)


# network of Swiss relationships of infections
library(huge)
cases.hg<-huge(as.matrix(lcases),method = "mb")
cases.shg<-huge.select(cases.hg,criterion = "stars")
plot(cases.shg)

# plot infection network
library(igraph)
adj<-cases.hg$path[[cases.shg$opt.index]]
g<-graph_from_adjacency_matrix(adj,mode="undirected")
g<-set_vertex_attr(g,name="label", value=colnames(lcases))
plot(g,vertex.size=25,vertex.color="white", edge.width=2, curved=TRUE)

# Swiss map with most infectious cantons
library(RSwissMaps)

# colour each canton by number of links
dt <- can.template(2016)
dt$values <- rowSums(as.matrix(adj))[match(dt$name,colnames(lcases))]
can.plot(dt$bfs_nr, dt$values, 2016)

# Causal network
library(pcalg)
suffStat = list(C=cor(lcases),n=ncol(lcases))
cases.pc<-pc(suffStat = suffStat, indepTest = gaussCItest,labels=colnames(lcases),alpha=0.3)
plot(cases.pc,main="directed causal graph")



###########################################################
# 3. Did the mortality rate improve during the pandemic?

#deal with recording issues
wkcases<-as.data.frame(apply(cases,2,function(x){tapply(x,rep(1:33,each=7),sum)}))
wkdeaths<-as.data.frame(apply(deaths,2,function(x){tapply(x,rep(1:33,each=7),sum)}))


plot(wkcases$TI,xlab="Week",ylab="number Covid19 cases",type="l")
plot(wkdeaths$TI,xlab="Week",ylab="number Covid19 cases",type="l")

# raw mortality probability 
plot(wkdeaths$TI/wkcases$TI,type="b",ylab="mortality rate",xlab="week")
title("No delay")
plot(wkdeaths$TI[-(1:2)]/wkcases$TI[-(32:33)],type="b",ylab="mortality rate",xlab="week")
title("Two week delay")

#logistic regression of mortality ratio
wks<-rep(1:31,27)
y<-cbind(unlist(wkdeaths[-c(1,2),]),unlist(wkcases[-c(32,33),]-wkdeaths[-c(1,2),]))
sel<-(apply(y,1,sum)>0)&(y[,2]>=0)
mortality.glm<-glm(y~wks,family=binomial(link = "logit"),subset = sel)
summary(mortality.glm)

par(mfrow=c(1,1))
plot(sort(wks[sel]),predict(mortality.glm,type = "response")[order(wks[sel])],type="l",xlab="Week",main = "Prb of dying if you get Covid",ylab="")
abline(0,0)