euromomo/Estimation_v42.R

   1 ##################################
   2 ### Version 4.2 - October 2018 ###
   3 ##################################
   4 ### Any questions or problems please contact: Jens Nielsen, Phone: +45 32683965 E-mail: nls@ssi.dk
   5
   6 # Install required packages
   7 #install.packages("ISOweeks")
   8
   9 ptrend <- 0.05
  10 p26 <- 0.05
  11 p52 <- 0.10
  12
  13 ### START: deaths data ###
  14 if (A_MOMO == 1) {
  15   deaths <- try(read.table(paste0(indir,"/A-MOMO data.txt"), sep=",", dec=".", header=T)[,c("group","YoDi","WoDi","nb","nbc")])
  16   if (inherits(deaths,"try-error")) {
  17     stop(paste0("Error: ",indir,"/A-MOMO data.txt"))
  18   }
  19   deaths$deaths <- pmax(deaths$nb, deaths$nbc)
  20   deaths$agegrp <- match(deaths$group, c("0to4","5to14","15to64","65P","Total")) - 1  # Convert group to integer (0-5)
  21   deaths <- deaths[!is.na(deaths$agegrp), c("agegrp", "YoDi", "WoDi", "deaths")]
  22   names(deaths)[2:3] <- c("year","week")
  23 }
  24 if (A_MOMO != 1) {
  25   deaths <- try(read.table(paste0(indir,"/deaths.txt"), sep=";", dec=".", header=T)[,c("agegrp","year","week","deaths")])
  26   if (inherits(deaths,"try-error")) {
  27     stop(paste0("Error: ",indir,"/deaths.txt"))
  28   }
  29 }
  30 deaths <- subset(deaths, (start_year*100+start_week <= year*100+week) & (year*100+week <= end_year*100+end_week))
  31 deaths[is.na(deaths)] <- 0
  32 deaths <- deaths[order(deaths$agegrp,deaths$year,deaths$week),]
  33 ### END: deaths data ###
  34
  35 ### Start: Temperature data ###
  36
  37 #' Temperature data from WeatherData
  38 #' PARM ET data containing then variables: date, pop3, NUTS3, temp
  39 #' PARM start_year first year to be included. Must be >= 2000
  40 #' PARM start_week first week in start_year to be included
  41 #' PARM end_year last year to be included
  42 #' PARM end_week last week in end_year to be included
  43 TemperatureData <- function(ET, start_year, start_week, end_year, end_week) {
  44   ET$date <- as.Date(as.character(ET$date,"%Y%m%d"))
  45   ET$NUTS3 <- as.character(ET$NUTS3)
  46   ET[is.na(ET$pop3), "pop3"] <- 1
  47   ET <- ET[order(ET$NUTS3, ET$date),]  # Ordering speeds summarizing a bit
  48   ET <- aggregate(cbind(temp, pop3) ~ NUTS3 + date, data=ET, function(x) mean(x, na.rm = TRUE))
  49   ET$pop3.sum <- with(ET, tapply(pop3, date, sum))[as.character(ET$date)]
  50   ET$temp <- ET$temp * ET$pop3 / ET$pop3.sum
  51   ET <- aggregate(temp ~ date, data = ET, function(x) sum(x, na.rm = TRUE))
  52   ET <- merge(ET, data.frame(date = seq.Date(as.Date(paste0(start_year,"/01/01")), as.Date(paste0(end_year+1,"/01/01")), by = "day")), by = "date", all = TRUE)
  53   ET$ISOweek <- ISOweek::ISOweek(ET$date)
  54   ET <- subset(ET, (paste0(start_year,"-W",formatC(start_week, width=2, flag="0")) <= ISOweek) &
  55                    (ISOweek <= paste0(end_year,"-W",formatC(end_week, width=2, flag="0"))))
  56   ET <- aggregate(temp ~ ISOweek, data = ET, function(x) c(mean(x, na.rm = TRUE), range(x, na.rm = TRUE)))
  57   ET <- cbind(as.data.frame(ET$ISOweek), as.data.frame(ET$temp))
  58   colnames(ET) <- c("ISOweek","temp","tmin","tmax")
  59   ET$wk <- order(ET$ISOweek)
  60   ET$sin52 <- sin((2*pi/(365.25/7))*ET$wk)
  61   ET$cos52 <- cos((2*pi/(365.25/7))*ET$wk)
  62   ET$ptemp <- predict(glm(temp ~ sin52 + cos52, data=ET[!(is.na(ET$temp) | is.infinite(ET$temp)),]), ET)
  63   ET$ptmin <- predict(glm(tmin ~ sin52 + cos52, data=ET[!(is.na(ET$tmin) | is.infinite(ET$tmin)),]), ET)
  64   ET$ptmax <- predict(glm(tmax ~ sin52 + cos52, data=ET[!(is.na(ET$tmax) | is.infinite(ET$tmax)),]), ET)
  65   ET[,c("wk","cos52","sin52")] <- NULL
  66   ET$t = ifelse(is.na(ET$temp) | is.infinite(ET$temp), ET$ptemp, ET$temp)
  67   ET$ET <- with(ET, (t-ptmax)*(t>ptmax)+(t-ptmin)*(t<ptmin))
  68   ET$year <- as.numeric(substr(ET$ISOweek,1,4))
  69   ET$week <- as.numeric(substr(ET$ISOweek,7,8))
  70   ET[, c("ISOweek", "t")] <- NULL
  71   return(ET <- ET[order(ET$year,ET$week),])
  72 }
  73
  74 # Read weather data
  75 if (WeatherData == 1) {
  76   ET <- try(read.table(paste0("http://euromomo.eu/methods/weather/wdata_",country.code,".txt"), header = TRUE, sep = ";", dec = ".")[, c("date","pop3","NUTS3","temp")])
  77   if (inherits(ET,"try-error")) {
  78     stop(paste0("Cannot download wdata_",country.code,".txt from http://euromomo.eu/methods/weather"))
  79   }
  80 }
  81 if (WeatherData != 1) {
  82   ET <- try(read.table(paste0(indir,"/wdata_",country.code,".txt"), header = TRUE, sep = ";", dec = ".")[, c("date","pop3","NUTS3","temp")])
  83   if (inherits(ET,"try-error")) {
  84     stop(paste0("Error: ",indir,"/wdata_",country.code,".txt"))
  85   }
  86 }
  87
  88 ET <- TemperatureData(ET,start_year,start_week,end_year,end_week)
  89
  90 write.table(ET, file=paste0(indir,"/temperature.txt"), row.names = FALSE, quote = FALSE, sep =";", dec=".")
  91
  92 rm(TemperatureData)
  93 ### END: Temperature data ###
  94
  95 ### START: Population data ###
  96 if (population) {
  97   pop_data <- try(read.table(paste0(indir,"/population.txt"), sep=";", dec=".", header=TRUE)[,c("agegrp","year","week","N")])
  98   if (inherits(pop_data,"try-error")) {
  99     stop(paste0("Error: ",indir,"/population.txt"))
 100   }
 101 }
 102 ### END: Population data ###
 103
 104 ### START: Merge data ###
 105 results4 <- try(read.table(paste0(indir,"/IA.txt"), sep=";", dec=".", header=TRUE))
 106 if (inherits(results4,"try-error")) {
 107   stop(paste0("Error: ",indir,"/IA.txt"))
 108 }
 109 results4 <- merge(deaths, results4[,c("agegrp","year","week","IA")],  by=c("agegrp","year","week"), all=FALSE)
 110 results4 <- merge(results4,ET[,c("year","week","ET")], by=c("year","week"), all=FALSE)
 111 if (population) {
 112   results4 <- merge(results4, pop_data,  by=c("agegrp","year","week"), all=FALSE)
 113   rm(pop_data)
 114 } else {
 115   results4$N <- 1
 116 }
 117 rm(deaths, ET)
 118 ### END: Merge data ###
 119
 120 ### START: Prepare data ###
 121
 122 # Function to get lag, forward and backward
 123 vecshift <- function(x, shift=0) {
 124   if (shift==0) return(x)
 125   if (shift>0) return(c(x[-(1:shift)], rep(NA,shift)))
 126   if (shift<0) return(c(rep(NA, -shift), x[1:(length(x)+shift)]))
 127 }
 128
 129 results4 <- subset(results4, (year*100+week >= start_year*100+start_week) & (year*100+week <= end_year*100+end_week))
 130 results4 <- results4[with(results4, order(agegrp, year, week)),]
 131
 132 results4$season <- with(results4, year - (week<27))
 133 results4$summer <- with(results4, week>=21 & results4$week<=39)
 134 results4$winter <- !results4$summer
 135
 136 # lags IA
 137 for (a in 0:4) {
 138   for (s in sort(unique(results4[results4$agegrp==a, "season"]))) {
 139     for (d in 0:IAlags) {
 140       results4[results4$agegrp==a, paste("d", d, "_IA", s, sep="")] <- vecshift(results4[results4$agegrp==a, "IA"], -d) * as.numeric(results4[results4$agegrp==a, "season"]==s)
 141     }
 142   }
 143 }
 144 # warm/cold summer/winter and lags
 145 for (s in c("summer","winter")) {
 146   results4[,paste0("cold_", s)] <- with(results4, -((ET<0) & get(s)) * ET)
 147   results4[,paste0("warm_", s)] <- with(results4, ((ET>0) & get(s)) * ET)
 148   for (t in c("cold","warm")) {
 149     for (d in 0:ETlags) {
 150       results4[,paste0("d", d, "_", t, "_", s)] <- NA
 151       for (a in 0:4) {
 152         results4[results4$agegrp==a,paste0("d", d, "_", t, "_", s)] <- vecshift(results4[results4$agegrp==a,paste0( t, "_", s)], -d)
 153       }
 154     }
 155   }
 156   results4[,paste0("cold_", s)] <- NULL
 157   results4[,paste0("warm_", s)] <- NULL
 158 }
 159 rm(a, s, t, d, vecshift)
 160 results4$summer <- NULL
 161 results4$winter <- NULL
 162
 163 # Create wk
 164 results4 <- results4[with(results4, order(agegrp, year, week)),]
 165 results4$wk <- unlist(lapply(table(results4$agegrp), seq))
 166
 167 results4$sin52 <- sin((2*pi/(365.25/7)) * results4$wk)
 168 results4$cos52 <- cos((2*pi/(365.25/7)) * results4$wk)
 169 results4$sin26 <- sin((4*pi/(365.25/7)) * results4$wk)
 170 results4$cos26 <- cos((4*pi/(365.25/7)) * results4$wk)
 171
 172 results4[is.na(results4)] <- 0
 173 ### END: Prepare data ###
 174
 175 ### START: Estimation ###
 176 results4$EB <- NA
 177 results4$VlogB <- NA
 178 results4$EIA <- NA
 179 results4$VlogIA <- NA
 180 results4$EET <- NA
 181 results4$VlogET <- NA
 182 results4$Vdeaths <- NA
 183
 184 for (a in 0:4) {
 185   print(paste("### Age group ",a,"###"))
 186   wk = "wk"
 187   f <- paste(c("deaths ~ ",wk," sin52 + cos52 + sin26 + cos26",
 188     grep("^d[0-9]", names(results4), value=TRUE)), collapse=" + ")
 189   m <- try(glm(f, quasipoisson, offset = log(N), data=results4[results4$agegrp==a,]))
 190   if (!inherits(m,"try-error") & m$converge & (median(results4[results4$agegrp==a, "deaths"]) > 0)) {
 191     fa <- paste(c("deaths ~ sin52 + cos52 + sin26 + cos26",
 192                   grep("^d[0-9]", names(results4), value=TRUE)), collapse=" + ")
 193     ma <- glm(fa, quasipoisson, offset = log(N), data=results4[results4$agegrp==a,])
 194     if (anova(m, ma, dispersion = max(1,sum(residuals(m, type = "deviance")^2)/df.residual(m)), test="LRT")$`Pr(>Chi)`[2]>ptrend) {
 195       wk = ""
 196       m <- ma
 197     }
 198     fa <- paste(c("deaths ~ ",wk, grep("^d[0-9]", names(results4), value=TRUE)), collapse=" + ")
 199     ma <- glm(fa, quasipoisson, offset = log(N), data=results4[results4$agegrp==a,])
 200     if (anova(m, ma, dispersion = max(1, sum(residuals(m, type = "deviance")^2)/df.residual(m)), test="LRT")$`Pr(>Chi)`[2]>max(p26,p52)) {
 201       m <- ma
 202     } else {
 203       fa <- paste(c("deaths ~ ",wk," + cos52 + sin52", grep("^d[0-9]", names(results4), value=TRUE)), collapse=" + ")
 204       ma <- glm(fa, quasipoisson, offset = log(N), data=results4[results4$agegrp==a,])
 205       if (anova(m, ma, dispersion = max(1, sum(residuals(m, type = "deviance")^2)/df.residual(m)), test="LRT")$`Pr(>Chi)`[2]>p26) {
 206         m <- ma
 207       } else {
 208         fa <- paste(c("deaths ~ ",wk, grep("^d[0-9]", names(results4), value=TRUE)), collapse=" + ")
 209         ma <- glm(fa, quasipoisson, offset = log(N), data=results4[results4$agegrp==a,])
 210         if (anova(m, ma, dispersion = max(1, sum(residuals(m, type = "deviance")^2)/df.residual(m)), test="LRT")$`Pr(>Chi)`[2]>p52) {
 211           m <- ma
 212         }
 213       }
 214     }
 215   } else {
 216     if (inherits(m,"try-error")) print("### Could not fit model ###")
 217     if (!m$converge) print("### Model did not converge ###")
 218     if (median(results4[results4$agegrp==a, "deaths"]) == 0) print("### Zero inflated ###")
 219     print("### Simple model with only trend used ###")
 220     f <- paste(c("deaths ~ wk"))
 221     m <- glm(f, quasipoisson, offset = log(N), data=results4[results4$agegrp==a,])
 222   }
 223   print(summary(m, dispersion = max(1, sum(residuals(m, type = "deviance")^2)/df.residual(m))))
 224   # Baseline
 225   results4.B <- results4[results4$agegrp==a,]
 226   for (d in 0:IAlags) {
 227     for (s in min(results4.B$season):max(results4.B$season)) {
 228       results4.B[,paste("d", d, "_IA", s, sep="")] <- 0
 229     }
 230     for (s in c("summer","winter")) {
 231       results4.B[,paste("d", d, "_warm_", s, sep="")] <- 0
 232       results4.B[,paste("d", d, "_cold_", s, sep="")] <- 0
 233     }
 234   }
 235   results4[results4$agegrp==a,]$EB <- exp(predict.glm(m, newdata=results4.B, se.fit=TRUE)$fit)
 236   results4[results4$agegrp==a,]$VlogB <- predict.glm(m, newdata=results4.B, dispersion = max(1, sum(residuals(m, type = "deviance")^2)/df.residual(m)), se.fit=TRUE)$se.fit^2
 237   # IA
 238   results4.IA <- results4[results4$agegrp==a,]
 239   for (s in c("summer","winter")) {
 240     for (d in 0:ETlags) {
 241       results4.IA[,paste("d", d, "_warm_", s, sep="")] <- 0
 242       results4.IA[,paste("d", d, "_cold_", s, sep="")] <- 0
 243     }
 244   }
 245   results4[results4$agegrp==a,]$EIA <- exp(predict.glm(m, newdata=results4.IA, se.fit=TRUE)$fit)
 246   results4[results4$agegrp==a,]$VlogIA <- predict.glm(m, newdata=results4.IA, dispersion = max(1, sum(residuals(m, type = "deviance")^2)/df.residual(m)), se.fit=TRUE)$se.fit^2
 247   # ET
 248   results4.ET <- results4[results4$agegrp==a,]
 249   for (d in 0:IAlags) {
 250     for (s in sort(unique(results4.ET$season))) {
 251       results4.ET[,paste("d", d, "_IA", s, sep="")] <- 0
 252     }
 253   }
 254   results4[results4$agegrp==a,]$EET <-  exp(predict.glm(m, newdata=results4.ET, se.fit=T)$fit)
 255   results4[results4$agegrp==a,]$VlogET <- predict.glm(m, newdata=results4.ET, dispersion = max(1, sum(residuals(m, type = "deviance")^2)/df.residual(m)), se.fit=T)$se.fit^2
 256
 257   results4[results4$agegrp==a,]$Vdeaths <- with(results4[results4$agegrp==a,], EB * max(1, sum(residuals(m, type = "deviance")^2)/df.residual(m)))
 258 }
 259
 260 # Delete external objects
 261 rm(a, f, fa, s, d, m, ma, results4.B, results4.ET, results4.IA, ptrend, p26, p52, wk)
 262 ### END: Estimation ###
 263
 264 ### START: Post estimation ###
 265 # Keep relevant variables #
 266 results4 <- results4[, c("agegrp", "year", "week", "deaths", "N", "IA", "ET", "Vdeaths", "EB", "VlogB", "EIA", "VlogIA", "EET", "VlogET")]
 267
 268 # Baseline, EIA and EET estimation variances - not on log scale
 269 results4$VB <- with(results4, (exp(VlogB)-1)*exp(2*log(EB) + VlogB))
 270 results4$VIA <- with(results4, (exp(VlogIA)-1)*exp(2*log(EIA) + VlogIA))
 271 results4$VET <- with(results4, (exp(VlogET)-1)*exp(2*log(EET) + VlogET))
 272
 273 # Effects of IA and ET
 274 results4$EdIA <- with(results4, EIA - EB)
 275 results4[is.na(results4$EdIA),"EdIA"] <- 0
 276 results4$EdET <- with(results4, EET - EB)
 277 results4[is.na(results4$EdET),"EdET"] <- 0
 278 # Excess relative to baseline
 279 results4$excess <- with(results4, deaths - EB)
 280 # Unexplained excess
 281 results4$uexcess <- with(results4, deaths - (EB + EdIA + EdET))
 282 # Exclude negative IA effects
 283 if (IArest) results4$uexcess <- with(results4, deaths - (EB + pmax(0,EdIA) + EdET))
 284
 285 # Baseline 2/3 residual confidence intervals
 286 results4$RVB <- with(results4, ((2/3)*(EB^(2/3-1))^2)*Vdeaths + ((2/3)*(EB^(2/3-1))^2)*VB)
 287 results4[with(results4, is.na(RVB) | is.infinite(RVB)), "RVB"] <- 0
 288 results4$EB_95L <- with(results4, pmax(0,sign((sign(EB)*abs(EB)^(2/3))-1.96*sqrt(RVB))*abs((sign(EB)*abs(EB)^(2/3))-1.96*sqrt(RVB))^(3/2)))
 289 results4$EB_95U <- with(results4, sign((sign(EB)*abs(EB)^(2/3))+1.96*sqrt(RVB))*abs((sign(EB)*abs(EB)^(2/3))+1.96*sqrt(RVB))^(3/2))
 290
 291 # EdIA 2/3 residual confidence intervals
 292 results4$RVdIA <- with(results4, ((2/3)*(EB^(2/3-1))^2)*VB + ((2/3)*(EIA^(2/3-1))^2)*VIA)
 293 results4[with(results4, is.na(RVdIA) | is.infinite(RVdIA)), "RVdIA"] <- 0
 294 results4$EdIA_95L <- with(results4, sign((sign(EdIA)*abs(EdIA)^(2/3))-1.96*sqrt(RVdIA))*abs((sign(EdIA)*abs(EdIA)^(2/3))-1.96*sqrt(RVdIA))^(3/2))
 295 results4$EdIA_95U <- with(results4, pmax(EdIA_95L,sign((sign(EdIA)*abs(EdIA)^(2/3))+1.96*sqrt(RVdIA))*abs((sign(EdIA)*abs(EdIA)^(2/3))+1.96*sqrt(RVdIA))^(3/2)))
 296
 297 # EdET 2/3 residual confidence intervals
 298 results4$RVdET <- with(results4, ((2/3)*(EB^(2/3-1))^2)*VB + ((2/3)*(EET^(2/3-1))^2)*VET)
 299 results4[with(results4, is.na(RVdET) | is.infinite(RVdET)), "RVdET"] <- 0
 300 results4$EdET_95L <- with(results4, sign((sign(EdET)*abs(EdET)^(2/3))-1.96*sqrt(RVdET))*abs((sign(EdET)*abs(EdET)^(2/3))-1.96*sqrt(RVdET))^(3/2))
 301 results4$EdET_95U <- with(results4, pmax(EdET_95L,sign((sign(EdET)*abs(EdET)^(2/3))+1.96*sqrt(RVdET))*abs((sign(EdET)*abs(EdET)^(2/3))+1.96*sqrt(RVdET))^(3/2)))
 302
 303 results4$summer <- with(results4, ifelse(week>=21 & week<=39, year, NA))
 304 results4$winter <- with(results4, ifelse(week<=20 | week>=40, year - (week<=20), NA))
 305
 306 # Order
 307 results4 <- results4[with(results4, order(agegrp, year, week)),]
 308
 309 cummulate <- function(x, pos=FALSE) {
 310   # x = variable, EdIA
 311   y <- x[,1]
 312   if (pos) y <- x[,1] * (x[,2] > 0)
 313   return(cumsum(y))
 314 }
 315 cCI <- function(x) {
 316   # x = cEB, cVB, cE*, cV*
 317   colnames(x) <- c("cEB", "cVB", "cE", "cV")
 318   x$cEd <- x$cE - x$cEB
 319   x$cRVd <- ((2/3)*(x$cEB^(2/3-1))^2)*x$cVB + ((2/3)*(x$cE^(2/3-1))^2)*x$cV
 320   x[with(x, is.na(cRVd) | is.infinite(cRVd)), "cRVd"] <- 0
 321   x$cCI_95L <- with(x, sign((sign(cEd)*abs(cEd)^(2/3))-1.96*sqrt(cRVd))*abs((sign(cEd)*abs(cEd)^(2/3))-1.96*sqrt(cRVd))^(3/2))
 322   x$cCI_95U <- with(x, sign((sign(cEd)*abs(cEd)^(2/3))+1.96*sqrt(cRVd))*abs((sign(cEd)*abs(cEd)^(2/3))+1.96*sqrt(cRVd))^(3/2))
 323   return(x[, c("cEd", "cCI_95L", "cCI_95U")])
 324 }
 325
 326 # EdIA and EdET residual variances
 327 for (s in c("summer", "winter", "year")) {
 328   nn <- !is.na(results4[,s])  # logical vector indicating season
 329   for (a in 0:4) {
 330     ag <- (results4$agegrp==a)
 331     for (v in c("excess","uexcess")) {
 332       results4[nn & ag ,paste0("c", v, "_", s)] <- NA
 333       for (sy in sort(unique(results4[,s]))) {
 334         results4[(results4[,s]==sy) & nn & ag, paste0("c", v, "_", s)] <-
 335           cummulate(results4[(results4[,s]==sy) & nn & ag, c(v, "EdIA")], FALSE)
 336       }
 337     }
 338     for (v in c("B","IA")) {
 339       results4[nn & ag ,paste0("cE", v, "_", s)] <- NA
 340       results4[nn & ag ,paste0("cV", v, "_", s)] <- NA
 341       for (sy in sort(unique(results4[,s]))) {
 342         results4[(results4[,s]==sy) & nn & ag, paste0("cE", v, "_", s)] <-
 343           cummulate(results4[(results4[,s]==sy) & nn & ag, c(paste0("E",v),"EdIA")], IArest)
 344         results4[(results4[,s]==sy) & nn & ag, paste0("cV", v, "_", s)] <-
 345           cummulate(results4[(results4[,s]==sy) & nn & ag, c(paste0("V",v),"EdIA")],  IArest)
 346       }
 347     }
 348     results4[nn & ag, paste0("cEd", v, "_", s)] <- NA
 349     results4[nn & ag, paste0("cEd", v, "_", s, "_95L")] <- NA
 350     results4[nn & ag, paste0("cEd", v, "_", s, "_95U")] <- NA
 351     results4[nn & ag, c(paste0("cEd", v, "_", s), paste0("cEd", v, "_", s, "_95L"), paste0("cEd", v,"_", s, "_95U"))] <-
 352       cCI(results4[nn & ag, c(paste0("cEB_", s), paste0("cVB_", s), paste0("cE",v,"_", s), paste0("cV",v,"_",s))])
 353     for (v in c("B","ET")) {
 354       results4[nn & ag ,paste0("cE", v, "_", s)] <- NA
 355       results4[nn & ag ,paste0("cV", v, "_", s)] <- NA
 356       for (sy in sort(unique(results4[,s]))) {
 357         results4[(results4[,s]==sy) & nn & ag, paste0("cE", v, "_", s)] <-
 358           cummulate(results4[(results4[,s]==sy) & nn & ag, c(paste0("E",v),"EdIA")], FALSE)
 359         results4[(results4[,s]==sy) & nn & ag, paste0("cV", v, "_", s)] <-
 360           cummulate(results4[(results4[,s]==sy) & nn & ag, c(paste0("V",v),"EdIA")],  FALSE)
 361       }
 362     }
 363     results4[nn & ag ,paste0("cEd", v, "_", s)] <- NA
 364     results4[nn & ag ,paste0("cEd", v, "_", s, "_95L")] <- NA
 365     results4[nn & ag ,paste0("cEd", v, "_", s, "_95U")] <- NA
 366     results4[nn & ag ,c(paste0("cEd", v, "_", s), paste0("cEd", v, "_", s, "_95L"), paste0("cEd", v, "_", s, "_95U"))] <-
 367       cCI(results4[nn & ag , c(paste0("cEB_", s), paste0("cVB_", s), paste0("cE",v,"_", s), paste0("cV",v,"_",s))])
 368   }
 369 }
 370
 371 results4$country <- country
 372 results4$IArestricted <- as.integer(IArest)
 373 results4 <- results4[, c("country", "IArestricted", "agegrp", "year", "week", "deaths", "Vdeaths", "N", "IA", "ET",
 374                          "EB", "EB_95L", "EB_95U", "VB",
 375                          "EIA", "VIA", "EET", "VET",
 376                          "EdIA", "EdIA_95L", "EdIA_95U",
 377                          "EdET", "EdET_95L", "EdET_95U",
 378                          "cexcess_year", "cuexcess_year",
 379                          "cEdIA_year", "cEdIA_year_95L", "cEdIA_year_95U",
 380                          "cEdET_year", "cEdET_year_95L", "cEdET_year_95U",
 381                          "summer", "cexcess_summer", "cuexcess_summer",
 382                          "cEdIA_summer", "cEdIA_summer_95L", "cEdIA_summer_95U",
 383                          "cEdET_summer", "cEdET_summer_95L", "cEdET_summer_95U",
 384                          "winter", "cexcess_winter", "cuexcess_winter",
 385                          "cEdIA_winter", "cEdIA_winter_95L", "cEdIA_winter_95U",
 386                          "cEdET_winter", "cEdET_winter_95L", "cEdET_winter_95U")]
 387
 388 if (IArest) {
 389   write.table(results4, file=paste0(outdir,"/",country,"_output_v4_IArestricted.txt"),
 390       row.names = FALSE, quote = FALSE, sep =";", dec=".", na="")
 391 } else {
 392   write.table(results4, file=paste0(outdir,"/",country,"_output_v4.txt"),
 393       row.names = FALSE, quote = FALSE, sep =";", dec=".", na="")
 394 }
 395
 396 rm(a, s, sy, v, nn, ag, cummulate, cCI, results4)
 397 ### END: Post estimation ###