euromomo/Output_calendar_v42.R

   1 ### Graphs - calendar time ###
   2
   3 if (IArest) {
   4   results4 <- read.table(file=paste0(outdir,"/",country,"_output_v4_IArestricted.txt"), header = TRUE, sep =";", dec=".")
   5   results4$EIA <- results4$EB + pmax(0,results4$EdIA)
   6   results4$EIAET <- results4$EB + pmax(0,results4$EdIA) + results4$EdET
   7   note <- "Note: Negative effects of IA are not shown"
   8 } else {
   9   results4 <- read.table(file=paste0(outdir,"/",country,"_output_v4.txt"), header = TRUE, sep =";", dec=".")
  10   results4$EIA <- results4$EB + results4$EdIA
  11   results4$EIAET <- results4$EB + results4$EdIA + results4$EdET
  12   note <- ""
  13 }
  14 results4 <- results4[with(results4, order(agegrp, year, week)),]
  15 results4$wk <- unlist(lapply(table(results4$agegrp), seq))
  16 results4$agegrp.labels <- ordered(results4$agegrp,
  17     labels = c("0-4 years","5-14 years","15-64 years","Aged 65","Total"))
  18 results4$yw <- sprintf("%04dw%02d", results4$year, results4$week)
  19
  20
  21 plotDeaths <- function(a, mortality=FALSE, annot=TRUE) {
  22   if (mortality) {
  23     results4[,c("deaths","EB","EB_95L","EB_95U","EIAET","EIA")] <-
  24       100000 * results4[,c("deaths","EB","EB_95L","EB_95U","EIAET","EIA")] / results4[,"N"]
  25   }
  26   with(results4[(results4$agegrp==a),], {
  27     plot(wk, EB, type="n", xaxt="n", bty="l", xlab=NA, ylab=NA,
  28         ylim=range(pretty(range(c(EIAET, EIA, EB_95U, EB_95L)))))
  29     if (annot) {
  30       axis(1, at=wk[as.integer(substr(yw,6,7)) %in% c(1,27)],
  31           labels=yw[as.integer(substr(yw,6,7)) %in% c(1,27)], las=2, cex.axis=0.9)
  32     }
  33     grid(nx=NA, ny=NULL)
  34     abline(v=wk[as.integer(substr(yw,6,7)) %in% c(1,27)], col="lightgray", lty="dotted")
  35     polygon(x=c(wk,rev(wk)), y=c(deaths,rev(EB)),
  36       col=do.call(rgb, c(as.list(col2rgb("gray")/255), alpha=0.7)),
  37       border=do.call(rgb, c(as.list(col2rgb("gray")/255), alpha=0.7)))
  38     points(wk, EIAET, col="darkgreen", lwd=2, type="l")
  39     points(wk, EIA, col="red", lwd=2, type="l")
  40     points(wk, EB, lwd=2, type="l")
  41     points(wk, EB_95L, lty="dashed", lwd=2, type="l")
  42     points(wk, EB_95U, lty="dashed", lwd=2, type="l")
  43   })
  44   if (annot) {
  45     mtext(sprintf("FluMOMO v4 - week %s, %s", end_week, end_year), side=3, line=1, adj=1, cex=0.8)
  46     mtext("year-week", side=1, line=4.5)
  47     if (mortality) {
  48       mtext(paste("Mortality by age group:",unique(results4$agegrp.labels)[a+1]), side=3, line=1, adj=0, cex=1.2)
  49       mtext("Deaths / 100.000 per week", side=2, line=2.5)
  50     } else {
  51       mtext(paste("Number of deaths, age group:",unique(results4$agegrp.labels)[a+1]), side=3, line=1, adj=0, cex=1.2)
  52       mtext("Number of deaths per week", side=2, line=2.5)
  53     }
  54     legend("bottom", c("Effect of Influenza Activity (IA)",
  55         "Effect of excess temperature (on top of IA)", "Baseline", "95% confidence interval"),
  56         lwd=2, col=c("red","darkgreen","black","black"), lty=c(rep("solid",3),"dashed"),
  57         bty="n", ncol=2, xpd=NA, inset=c(0,-0.5), seg.len=4)
  58   }
  59 }
  60
  61
  62 plotDeathsMultiple <- function(mortality=FALSE) {
  63   par(mar=c(3,2,3,2), mfrow=c(5,1), oma=c(8,3.5,3,0))
  64   for(a in 0:4) {
  65     plotDeaths(a, mortality=mortality, annot=FALSE)
  66     mtext(unique(results4$agegrp.labels)[a+1], side=3, line=1)
  67     if (a==4) {
  68       with(results4[(results4$agegrp==a),],
  69         axis(1, at=wk[as.integer(substr(yw,6,7)) %in% c(1,27)],
  70           labels=yw[as.integer(substr(yw,6,7)) %in% c(1,27)], las=2))
  71     }
  72   }
  73   mtext(sprintf("FluMOMO v4 - week %s, %s", end_week, end_year), side=3, line=0, adj=0.9, cex=0.8, outer=TRUE)
  74   mtext("year-week", side=1, line=5)
  75   if (mortality) {
  76     mtext("Mortality by age group", side=3, line=0, adj=0, cex=1.2, outer=TRUE)
  77     mtext("Deaths / 100.000 per week", side=2.5, line=1, outer=TRUE)
  78   } else {
  79     mtext("Number of deaths by age group", side=3, line=0, adj=0, cex=1.2, outer=TRUE)
  80     mtext("Number of deaths per week", side=2.5, line=1, outer=TRUE)
  81   }
  82   legend("bottom", c("Effect of Influenza Activity (IA)",
  83       "Effect of excess temperature (on top of IA)", "Baseline", "95% confidence interval"),
  84       lwd=2, col=c("red","darkgreen","black","black"), lty=c(rep("solid",3),"dashed"),
  85       bty="n", ncol=2, xpd=NA, inset=c(0,-0.9), seg.len=4)
  86 }
  87
  88
  89 png(paste0(outdir,"/deaths_agegroups_v4", ifelse(IArest, "_IArestricted", ""), ".png"), width=1200, height=1600, res=130)
  90 plotDeathsMultiple(mortality=FALSE)
  91 dev.off()
  92 for (a in 0:4) {
  93   png(paste0(outdir, "/deaths_agegroup_", a, "_v4", ifelse(IArest, "_IArestricted", ""), ".png"), width=1200, height=800, res=130)
  94   par(mar=c(10,4,3,2))
  95   plotDeaths(a, mortality=FALSE)
  96   dev.off()
  97 }
  98
  99 if (population) {
 100   png(paste0(outdir,"/mr_agegroups_v4", ifelse(IArest, "_IArestricted", ""), ".png"), width=1200, height=1600, res=130)
 101   plotDeathsMultiple(mortality=TRUE)
 102   dev.off()
 103   for (a in 0:4) {
 104     png(paste0(outdir,"/mr_agegroup", a, "_v4", ifelse(IArest, "_IArestricted", ""), ".png"), width=1200, height=800, res=130)
 105     par(mar=c(10,4,3,2))
 106     plotDeaths(a, mortality=TRUE)
 107     dev.off()
 108   }
 109 }
 110 rm(plotDeaths, plotDeathsMultiple, results4, a, note)