# 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)