+++ /dev/null
-### Graphs - calendar time ###\r
-\r
-if (IArest) {\r
- results4 <- read.table(file=paste0(outdir,"/",country,"_output_v4_IArestricted.txt"), header = TRUE, sep =";", dec=".")\r
- results4$EIA <- results4$EB + pmax(0,results4$EdIA)\r
- results4$EIAET <- results4$EB + pmax(0,results4$EdIA) + results4$EdET\r
- note <- "Note: Negative effects of IA are not shown"\r
-} else {\r
- results4 <- read.table(file=paste0(outdir,"/",country,"_output_v4.txt"), header = TRUE, sep =";", dec=".")\r
- results4$EIA <- results4$EB + results4$EdIA\r
- results4$EIAET <- results4$EB + results4$EdIA + results4$EdET\r
- note <- ""\r
-}\r
-results4 <- results4[with(results4, order(agegrp, year, week)),]\r
-results4$wk <- unlist(lapply(table(results4$agegrp), seq))\r
-results4$agegrp.labels <- ordered(results4$agegrp, \r
- labels = c("0-4 years","5-14 years","15-64 years","Aged 65","Total"))\r
-results4$yw <- sprintf("%04dw%02d", results4$year, results4$week)\r
-\r
-\r
-plotDeaths <- function(a, mortality=FALSE, annot=TRUE) {\r
- if (mortality) {\r
- results4[,c("deaths","EB","EB_95L","EB_95U","EIAET","EIA")] <- \r
- 100000 * results4[,c("deaths","EB","EB_95L","EB_95U","EIAET","EIA")] / results4[,"N"]\r
- }\r
- with(results4[(results4$agegrp==a),], {\r
- plot(wk, EB, type="n", xaxt="n", bty="l", xlab=NA, ylab=NA, \r
- ylim=range(pretty(range(c(EIAET, EIA, EB_95U, EB_95L)))))\r
- if (annot) {\r
- axis(1, at=wk[as.integer(substr(yw,6,7)) %in% c(1,27)], \r
- labels=yw[as.integer(substr(yw,6,7)) %in% c(1,27)], las=2, cex.axis=0.9)\r
- }\r
- grid(nx=NA, ny=NULL)\r
- abline(v=wk[as.integer(substr(yw,6,7)) %in% c(1,27)], col="lightgray", lty="dotted")\r
- polygon(x=c(wk,rev(wk)), y=c(deaths,rev(EB)), \r
- col=do.call(rgb, c(as.list(col2rgb("gray")/255), alpha=0.7)),\r
- border=do.call(rgb, c(as.list(col2rgb("gray")/255), alpha=0.7)))\r
- points(wk, EIAET, col="darkgreen", lwd=2, type="l")\r
- points(wk, EIA, col="red", lwd=2, type="l")\r
- points(wk, EB, lwd=2, type="l")\r
- points(wk, EB_95L, lty="dashed", lwd=2, type="l")\r
- points(wk, EB_95U, lty="dashed", lwd=2, type="l")\r
- })\r
- if (annot) {\r
- mtext(sprintf("FluMOMO v4 - week %s, %s", end_week, end_year), side=3, line=1, adj=1, cex=0.8)\r
- mtext("year-week", side=1, line=4.5)\r
- if (mortality) {\r
- mtext(paste("Mortality by age group:",unique(results4$agegrp.labels)[a+1]), side=3, line=1, adj=0, cex=1.2)\r
- mtext("Deaths / 100.000 per week", side=2, line=2.5)\r
- } else {\r
- mtext(paste("Number of deaths, age group:",unique(results4$agegrp.labels)[a+1]), side=3, line=1, adj=0, cex=1.2)\r
- mtext("Number of deaths per week", side=2, line=2.5)\r
- }\r
- legend("bottom", c("Effect of Influenza Activity (IA)",\r
- "Effect of excess temperature (on top of IA)", "Baseline", "95% confidence interval"), \r
- lwd=2, col=c("red","darkgreen","black","black"), lty=c(rep("solid",3),"dashed"), \r
- bty="n", ncol=2, xpd=NA, inset=c(0,-0.5), seg.len=4)\r
- }\r
-}\r
-\r
-\r
-plotDeathsMultiple <- function(mortality=FALSE) {\r
- par(mar=c(3,2,3,2), mfrow=c(5,1), oma=c(8,3.5,3,0))\r
- for(a in 0:4) {\r
- plotDeaths(a, mortality=mortality, annot=FALSE)\r
- mtext(unique(results4$agegrp.labels)[a+1], side=3, line=1)\r
- if (a==4) {\r
- with(results4[(results4$agegrp==a),], \r
- axis(1, at=wk[as.integer(substr(yw,6,7)) %in% c(1,27)], \r
- labels=yw[as.integer(substr(yw,6,7)) %in% c(1,27)], las=2))\r
- }\r
- }\r
- mtext(sprintf("FluMOMO v4 - week %s, %s", end_week, end_year), side=3, line=0, adj=0.9, cex=0.8, outer=TRUE)\r
- mtext("year-week", side=1, line=5)\r
- if (mortality) {\r
- mtext("Mortality by age group", side=3, line=0, adj=0, cex=1.2, outer=TRUE)\r
- mtext("Deaths / 100.000 per week", side=2.5, line=1, outer=TRUE)\r
- } else {\r
- mtext("Number of deaths by age group", side=3, line=0, adj=0, cex=1.2, outer=TRUE)\r
- mtext("Number of deaths per week", side=2.5, line=1, outer=TRUE)\r
- }\r
- legend("bottom", c("Effect of Influenza Activity (IA)",\r
- "Effect of excess temperature (on top of IA)", "Baseline", "95% confidence interval"), \r
- lwd=2, col=c("red","darkgreen","black","black"), lty=c(rep("solid",3),"dashed"), \r
- bty="n", ncol=2, xpd=NA, inset=c(0,-0.9), seg.len=4)\r
-}\r
-\r
-\r
-png(paste0(outdir,"/deaths_agegroups_v4", ifelse(IArest, "_IArestricted", ""), ".png"), width=1200, height=1600, res=130)\r
-plotDeathsMultiple(mortality=FALSE)\r
-dev.off()\r
-for (a in 0:4) {\r
- png(paste0(outdir, "/deaths_agegroup_", a, "_v4", ifelse(IArest, "_IArestricted", ""), ".png"), width=1200, height=800, res=130)\r
- par(mar=c(10,4,3,2))\r
- plotDeaths(a, mortality=FALSE)\r
- dev.off()\r
-}\r
-\r
-if (population) {\r
- png(paste0(outdir,"/mr_agegroups_v4", ifelse(IArest, "_IArestricted", ""), ".png"), width=1200, height=1600, res=130)\r
- plotDeathsMultiple(mortality=TRUE)\r
- dev.off() \r
- for (a in 0:4) {\r
- png(paste0(outdir,"/mr_agegroup", a, "_v4", ifelse(IArest, "_IArestricted", ""), ".png"), width=1200, height=800, res=130)\r
- par(mar=c(10,4,3,2))\r
- plotDeaths(a, mortality=TRUE)\r
- dev.off()\r
- }\r
-}\r
-rm(plotDeaths, plotDeathsMultiple, results4, a, note)\r
projects / covid19.git / blobdiff