euromomo/Output_cumulated_v42.R

   1 ### Graphs - cumulated ###
   2
   3 results4 <- read.table(
   4       paste0(outdir, "/", country, "_output_v4", ifelse(IArest, "_IArestricted", ""), ".txt"),
   5       header = TRUE, sep =";", dec=".")
   6 note <- ifelse(IArest, "Note: Negative effects of IA are not shown", "")
   7 results4$agegrp.labels <- factor(results4$agegrp,
   8     labels = c("0-4 years", "5-14 years", "15-64 years", "Aged 65", "Total"))
   9
  10
  11 pal <- c("#E41A1C", "#377EB8", "#4DAF4A", "#984EA3", "#FF7F00", "#A65628", "#F781BF")
  12
  13 plotCumulSingle <- function(a, mortality=FALSE, lwd=2, ylab=NULL) {
  14   if (mortality) {
  15     yvar <- "cm"
  16     ylb <- "Deaths / 100,000"
  17     title <- paste0("Cumulated ",s," mortality attributable to ",t,", age group: ",unique(results4$agegrp.labels)[a+1])
  18   } else {
  19     yvar <- "cd"
  20     ylb <- "Number of deaths"
  21     title <- paste0("Cumulated ",s," number of deaths attributable to ",t,", age group: ",unique(results4$agegrp.labels)[a+1])
  22   }
  23   if (!is.null(ylab)) ylb <- ylab
  24   xlb <- "Week number"
  25   if (is.na(ylb)) xlb <- NA
  26
  27   with(out[(out$agegrp==a),], {
  28     plot(x=wk, y=get(yvar), type="n", xaxt="n", bty="l",
  29       xlab=xlb, ylab=ylb,
  30       col=pal[1:length(levels(season))][as.integer(season)]
  31     )
  32     if (!is.na(ylb)) {
  33       mtext(title, side=3, line=2, adj=0)
  34       mtext(sprintf("FluMOMO v4 - week %s, %s", end_week, end_year), side=3, line=2, adj=1, cex=0.8)
  35     }
  36     axis(1, at=seq(1, max(wk), 4), labels = unique(floor(week))[seq(1, max(wk), 4)])
  37     grid(nx=NA, ny=NULL)
  38     abline(v=seq(1, max(wk), 4), col="lightgray", lty="dotted")
  39     mapply(function(s, col, l) {
  40       points(x=wk[season==s], y=get(yvar)[season==s], col=col, type="l", lty=l, lwd=lwd)
  41     }, s=levels(season), col=rev(pal[1:length(levels(season))]),
  42     l=c(rep("dotted", length(levels(season))-1), "solid"))
  43   })
  44 }
  45
  46
  47 plotCumulMultiple <- function(mortality=FALSE, lwd=2, ylab=NULL) {
  48   par(mfrow=c(5,1), oma=c(6,3,2,0))
  49   for (a in 0:4) {
  50     par(mar=c(3,2,3,2))
  51     plotCumulSingle(a, mortality, ylab=NA)
  52     legend("topleft", levels(results4$agegrp.labels)[a+1], bty="n", xpd=NA, cex=1.3, inset=c(0,-0.2))
  53   }
  54   mtext(paste0("Cumulated ",s, ifelse(mortality, " mortality", " number of deaths"), " attributable to ",t,", age group: ",unique(results4$agegrp.labels)[a+1]), side=3, line=0, adj=0, outer=TRUE)
  55   mtext(sprintf("FluMOMO v4 - week %s, %s", end_week, end_year), side=3, line=0, adj=0.95, cex=0.8, outer=TRUE)
  56   mtext(ifelse(mortality, "Deaths / 100,000", "Number of deaths"), side=2, line=1, outer=TRUE)
  57   mtext("Week number", side=1, line=3, cex=0.8)
  58   legend("bottom", legend=rev(levels(out$season)),
  59     lty = c("solid", rep("dotted", length(levels(out$season))-1)),
  60     col = pal[1:length(levels(out$season))], lwd=2, cex=1.2,
  61     bty="n", horiz=TRUE, xpd=NA, inset=c(0,-0.59))
  62 }
  63
  64
  65 for (s in c("summer","winter","year")) {
  66   for (t in c("IA")) {
  67     out <- results4[!is.na(results4[,s]),
  68         c("agegrp","agegrp.labels","year","week",paste0("cEd",t,"_",s),"N")]
  69     colnames(out) <- c("agegrp","agegrp.labels","year","week","cd","N")
  70     out$cm <- 100000 * out$cd / out$N
  71     out[out$week==53,"week"] <- 52.5
  72     if (s == "winter") {
  73       out$wk <- with(out, 1 + (week-40)*(week>=40) + (week+12)*(week<=20))
  74       out$season <- factor(paste0(as.character(with(out, year*(week>=40) + (year-1)*(week<=20))),"/",substr(as.character(with(out, year*(week>=40) + (year-1)*(week<=20))+1),3,4)))
  75     }
  76     if (s == "summer") {
  77       out$wk <- out$week-20
  78       out$season <- factor(out$year)
  79     }
  80     if (s == "year") {
  81       out$wk <- out$week
  82       out$season <- factor(out$year)
  83     }
  84     out <- out[with(out, order(agegrp, wk)),]
  85     png(paste0(outdir,"/cumulated_",t,"_deaths_agegroups_",s,"_v4", ifelse(IArest, "_IArestricted", ""), ".png"), width=1200, height=1600, res=130)
  86     plotCumulMultiple()
  87     dev.off()
  88     if (population) {
  89       png(paste0(outdir,"/cumulated_",t,"_mr_agegroups_",s,"_v4", ifelse(IArest, "_IArestricted", ""), ".png"), width=1200, height=1600, res=130)
  90       plotCumulMultiple(mortality=TRUE)
  91       dev.off()
  92     }
  93     for (a in 0:4) {
  94       png(paste0(outdir,"/cumulated_",t,"_deaths_agegroup_",a,"_",s,"_v4", ifelse(IArest, "_IArestricted", ""), ".png"), width=1200, height=800, res=130)
  95       par(oma=c(3,0,0,0))
  96       plotCumulSingle(a)
  97       legend("bottom", legend=rev(levels(out$season)),
  98         lty = c("solid", rep("dotted", length(levels(out$season))-1)),
  99         col = pal[1:length(levels(out$season))], lwd=2,
 100         bty="n", horiz=TRUE, xpd=NA, inset=c(0,-0.35))
 101       dev.off()
 102       if (population) {
 103         png(paste0(outdir,"/cumulated_",t,"_mr_agegroup_",a,"_",s,"_v4", ifelse(IArest, "_IArestricted", ""), ".png"), width=1200, height=800, res=130)
 104         par(oma=c(3,0,0,0))
 105         plotCumulSingle(a, mortality=TRUE)
 106         legend("bottom", legend=rev(levels(out$season)),
 107           lty = c("solid", rep("dotted", length(levels(out$season))-1)),
 108           col = pal[1:length(levels(out$season))], lwd=2,
 109           bty="n", horiz=TRUE, xpd=NA, inset=c(0,-0.35))
 110         dev.off()
 111       }
 112     }
 113   }
 114 }
 115 rm(plotCumulSingle, plotCumulMultiple, results4, a, s, t, out, note, pal)