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 your canton?

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

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

# plot SQRT transformed data
plot(dates,sqrt(cases$AG),type="l")
plot(dates,sqrt(deaths$AG),type="l")

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


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

# plot relationship between cantons
plot(lcases[,27],lcases[,1],xlab="number of cases in ZH",
      ylab="number of cases in AG")

plot(lcases[,4],lcases[,1],xlab="number of cases in BE",
     ylab="number of cases in AG")

# Fit linear regressions
fit1<-lm(lcases[,1]~lcases[,27])
summary(fit1)
plot(lcases[,27],lcases[,1],xlab="number of cases in ZH",
     ylab="number of cases in AG")
abline(coef(fit1),col="blue")

fit2<-lm(lcases[,1]~lcases[,4])
summary(fit2)
plot(lcases[,4],lcases[,1],xlab="number of cases in BE",
     ylab="number of cases in AG")
abline(coef(fit2),col="blue")


# Multiple regression
fit<-lm(lcases[,1]~lcases[,27]+lcases[,4])
summary(fit)


# network of Swiss relationships of infections
library(huge)
cases.hg<-huge(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)



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

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

plot(wkcases[,1],xlab="Week",ylab="number Covid19 cases",type="l")
plot(wkdeaths[,1],xlab="Week",ylab="number Covid19 cases",type="l")

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

#logistic regression of mortality ratio
wks<-rep(1:31,27)
y<-cbind(as.vector(wkdeaths[-c(1,2),]),as.vector(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)