+++ /dev/null
-##################################\r
-### Version 4.2 - October 2018 ###\r
-##################################\r
-### Any questions or problems please contact: Jens Nielsen, Phone: +45 32683965 E-mail: nls@ssi.dk\r
-\r
-# Install required packages\r
-#install.packages("ISOweeks")\r
-\r
-ptrend <- 0.05\r
-p26 <- 0.05\r
-p52 <- 0.10\r
-\r
-### START: deaths data ###\r
-if (A_MOMO == 1) {\r
- deaths <- try(read.table(paste0(indir,"/A-MOMO data.txt"), sep=",", dec=".", header=T)[,c("group","YoDi","WoDi","nb","nbc")])\r
- if (inherits(deaths,"try-error")) {\r
- stop(paste0("Error: ",indir,"/A-MOMO data.txt"))\r
- }\r
- deaths$deaths <- pmax(deaths$nb, deaths$nbc)\r
- deaths$agegrp <- match(deaths$group, c("0to4","5to14","15to64","65P","Total")) - 1 # Convert group to integer (0-5)\r
- deaths <- deaths[!is.na(deaths$agegrp), c("agegrp", "YoDi", "WoDi", "deaths")]\r
- names(deaths)[2:3] <- c("year","week")\r
-}\r
-if (A_MOMO != 1) {\r
- deaths <- try(read.table(paste0(indir,"/deaths.txt"), sep=";", dec=".", header=T)[,c("agegrp","year","week","deaths")])\r
- if (inherits(deaths,"try-error")) {\r
- stop(paste0("Error: ",indir,"/deaths.txt"))\r
- }\r
-}\r
-deaths <- subset(deaths, (start_year*100+start_week <= year*100+week) & (year*100+week <= end_year*100+end_week))\r
-deaths[is.na(deaths)] <- 0\r
-deaths <- deaths[order(deaths$agegrp,deaths$year,deaths$week),]\r
-### END: deaths data ###\r
-\r
-### Start: Temperature data ###\r
-\r
-#' Temperature data from WeatherData\r
-#' PARM ET data containing then variables: date, pop3, NUTS3, temp\r
-#' PARM start_year first year to be included. Must be >= 2000\r
-#' PARM start_week first week in start_year to be included\r
-#' PARM end_year last year to be included\r
-#' PARM end_week last week in end_year to be included\r
-TemperatureData <- function(ET, start_year, start_week, end_year, end_week) {\r
- ET$date <- as.Date(as.character(ET$date,"%Y%m%d"))\r
- ET$NUTS3 <- as.character(ET$NUTS3)\r
- ET[is.na(ET$pop3), "pop3"] <- 1\r
- ET <- ET[order(ET$NUTS3, ET$date),] # Ordering speeds summarizing a bit\r
- ET <- aggregate(cbind(temp, pop3) ~ NUTS3 + date, data=ET, function(x) mean(x, na.rm = TRUE))\r
- ET$pop3.sum <- with(ET, tapply(pop3, date, sum))[as.character(ET$date)]\r
- ET$temp <- ET$temp * ET$pop3 / ET$pop3.sum\r
- ET <- aggregate(temp ~ date, data = ET, function(x) sum(x, na.rm = TRUE))\r
- 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)\r
- ET$ISOweek <- ISOweek::ISOweek(ET$date)\r
- ET <- subset(ET, (paste0(start_year,"-W",formatC(start_week, width=2, flag="0")) <= ISOweek) &\r
- (ISOweek <= paste0(end_year,"-W",formatC(end_week, width=2, flag="0"))))\r
- ET <- aggregate(temp ~ ISOweek, data = ET, function(x) c(mean(x, na.rm = TRUE), range(x, na.rm = TRUE)))\r
- ET <- cbind(as.data.frame(ET$ISOweek), as.data.frame(ET$temp))\r
- colnames(ET) <- c("ISOweek","temp","tmin","tmax")\r
- ET$wk <- order(ET$ISOweek)\r
- ET$sin52 <- sin((2*pi/(365.25/7))*ET$wk)\r
- ET$cos52 <- cos((2*pi/(365.25/7))*ET$wk)\r
- ET$ptemp <- predict(glm(temp ~ sin52 + cos52, data=ET[!(is.na(ET$temp) | is.infinite(ET$temp)),]), ET)\r
- ET$ptmin <- predict(glm(tmin ~ sin52 + cos52, data=ET[!(is.na(ET$tmin) | is.infinite(ET$tmin)),]), ET)\r
- ET$ptmax <- predict(glm(tmax ~ sin52 + cos52, data=ET[!(is.na(ET$tmax) | is.infinite(ET$tmax)),]), ET)\r
- ET[,c("wk","cos52","sin52")] <- NULL\r
- ET$t = ifelse(is.na(ET$temp) | is.infinite(ET$temp), ET$ptemp, ET$temp)\r
- ET$ET <- with(ET, (t-ptmax)*(t>ptmax)+(t-ptmin)*(t<ptmin))\r
- ET$year <- as.numeric(substr(ET$ISOweek,1,4))\r
- ET$week <- as.numeric(substr(ET$ISOweek,7,8))\r
- ET[, c("ISOweek", "t")] <- NULL\r
- return(ET <- ET[order(ET$year,ET$week),])\r
-}\r
-\r
-# Read weather data\r
-if (WeatherData == 1) {\r
- ET <- try(read.table(paste0("http://euromomo.eu/methods/weather/wdata_",country.code,".txt"), header = TRUE, sep = ";", dec = ".")[, c("date","pop3","NUTS3","temp")])\r
- if (inherits(ET,"try-error")) {\r
- stop(paste0("Cannot download wdata_",country.code,".txt from http://euromomo.eu/methods/weather"))\r
- }\r
-}\r
-if (WeatherData != 1) {\r
- ET <- try(read.table(paste0(indir,"/wdata_",country.code,".txt"), header = TRUE, sep = ";", dec = ".")[, c("date","pop3","NUTS3","temp")])\r
- if (inherits(ET,"try-error")) {\r
- stop(paste0("Error: ",indir,"/wdata_",country.code,".txt"))\r
- }\r
-}\r
-\r
-ET <- TemperatureData(ET,start_year,start_week,end_year,end_week)\r
-\r
-write.table(ET, file=paste0(indir,"/temperature.txt"), row.names = FALSE, quote = FALSE, sep =";", dec=".")\r
-\r
-rm(TemperatureData)\r
-### END: Temperature data ###\r
-\r
-### START: Population data ###\r
-if (population) {\r
- pop_data <- try(read.table(paste0(indir,"/population.txt"), sep=";", dec=".", header=TRUE)[,c("agegrp","year","week","N")])\r
- if (inherits(pop_data,"try-error")) {\r
- stop(paste0("Error: ",indir,"/population.txt"))\r
- }\r
-} \r
-### END: Population data ###\r
-\r
-### START: Merge data ###\r
-results4 <- try(read.table(paste0(indir,"/IA.txt"), sep=";", dec=".", header=TRUE))\r
-if (inherits(results4,"try-error")) {\r
- stop(paste0("Error: ",indir,"/IA.txt"))\r
-}\r
-results4 <- merge(deaths, results4[,c("agegrp","year","week","IA")], by=c("agegrp","year","week"), all=FALSE)\r
-results4 <- merge(results4,ET[,c("year","week","ET")], by=c("year","week"), all=FALSE)\r
-if (population) {\r
- results4 <- merge(results4, pop_data, by=c("agegrp","year","week"), all=FALSE)\r
- rm(pop_data)\r
-} else {\r
- results4$N <- 1\r
-}\r
-rm(deaths, ET)\r
-### END: Merge data ###\r
-\r
-### START: Prepare data ###\r
-\r
-# Function to get lag, forward and backward\r
-vecshift <- function(x, shift=0) {\r
- if (shift==0) return(x)\r
- if (shift>0) return(c(x[-(1:shift)], rep(NA,shift)))\r
- if (shift<0) return(c(rep(NA, -shift), x[1:(length(x)+shift)]))\r
-}\r
-\r
-results4 <- subset(results4, (year*100+week >= start_year*100+start_week) & (year*100+week <= end_year*100+end_week))\r
-results4 <- results4[with(results4, order(agegrp, year, week)),]\r
-\r
-results4$season <- with(results4, year - (week<27))\r
-results4$summer <- with(results4, week>=21 & results4$week<=39)\r
-results4$winter <- !results4$summer\r
-\r
-# lags IA\r
-for (a in 0:4) {\r
- for (s in sort(unique(results4[results4$agegrp==a, "season"]))) {\r
- for (d in 0:IAlags) {\r
- 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)\r
- }\r
- }\r
-}\r
-# warm/cold summer/winter and lags\r
-for (s in c("summer","winter")) {\r
- results4[,paste0("cold_", s)] <- with(results4, -((ET<0) & get(s)) * ET)\r
- results4[,paste0("warm_", s)] <- with(results4, ((ET>0) & get(s)) * ET)\r
- for (t in c("cold","warm")) {\r
- for (d in 0:ETlags) {\r
- results4[,paste0("d", d, "_", t, "_", s)] <- NA\r
- for (a in 0:4) {\r
- results4[results4$agegrp==a,paste0("d", d, "_", t, "_", s)] <- vecshift(results4[results4$agegrp==a,paste0( t, "_", s)], -d)\r
- }\r
- }\r
- }\r
- results4[,paste0("cold_", s)] <- NULL\r
- results4[,paste0("warm_", s)] <- NULL\r
-}\r
-rm(a, s, t, d, vecshift)\r
-results4$summer <- NULL\r
-results4$winter <- NULL\r
-\r
-# Create wk\r
-results4 <- results4[with(results4, order(agegrp, year, week)),]\r
-results4$wk <- unlist(lapply(table(results4$agegrp), seq))\r
-\r
-results4$sin52 <- sin((2*pi/(365.25/7)) * results4$wk)\r
-results4$cos52 <- cos((2*pi/(365.25/7)) * results4$wk)\r
-results4$sin26 <- sin((4*pi/(365.25/7)) * results4$wk)\r
-results4$cos26 <- cos((4*pi/(365.25/7)) * results4$wk)\r
-\r
-results4[is.na(results4)] <- 0\r
-### END: Prepare data ###\r
-\r
-### START: Estimation ###\r
-results4$EB <- NA\r
-results4$VlogB <- NA\r
-results4$EIA <- NA\r
-results4$VlogIA <- NA\r
-results4$EET <- NA\r
-results4$VlogET <- NA\r
-results4$Vdeaths <- NA\r
-\r
-for (a in 0:4) {\r
- print(paste("### Age group ",a,"###"))\r
- wk = "wk"\r
- f <- paste(c("deaths ~ ",wk," sin52 + cos52 + sin26 + cos26", \r
- grep("^d[0-9]", names(results4), value=TRUE)), collapse=" + ")\r
- m <- try(glm(f, quasipoisson, offset = log(N), data=results4[results4$agegrp==a,]))\r
- if (!inherits(m,"try-error") & m$converge & (median(results4[results4$agegrp==a, "deaths"]) > 0)) {\r
- fa <- paste(c("deaths ~ sin52 + cos52 + sin26 + cos26",\r
- grep("^d[0-9]", names(results4), value=TRUE)), collapse=" + ")\r
- ma <- glm(fa, quasipoisson, offset = log(N), data=results4[results4$agegrp==a,])\r
- if (anova(m, ma, dispersion = max(1,sum(residuals(m, type = "deviance")^2)/df.residual(m)), test="LRT")$`Pr(>Chi)`[2]>ptrend) {\r
- wk = ""\r
- m <- ma\r
- }\r
- fa <- paste(c("deaths ~ ",wk, grep("^d[0-9]", names(results4), value=TRUE)), collapse=" + ")\r
- ma <- glm(fa, quasipoisson, offset = log(N), data=results4[results4$agegrp==a,])\r
- if (anova(m, ma, dispersion = max(1, sum(residuals(m, type = "deviance")^2)/df.residual(m)), test="LRT")$`Pr(>Chi)`[2]>max(p26,p52)) {\r
- m <- ma\r
- } else {\r
- fa <- paste(c("deaths ~ ",wk," + cos52 + sin52", grep("^d[0-9]", names(results4), value=TRUE)), collapse=" + ")\r
- ma <- glm(fa, quasipoisson, offset = log(N), data=results4[results4$agegrp==a,])\r
- if (anova(m, ma, dispersion = max(1, sum(residuals(m, type = "deviance")^2)/df.residual(m)), test="LRT")$`Pr(>Chi)`[2]>p26) {\r
- m <- ma\r
- } else {\r
- fa <- paste(c("deaths ~ ",wk, grep("^d[0-9]", names(results4), value=TRUE)), collapse=" + ")\r
- ma <- glm(fa, quasipoisson, offset = log(N), data=results4[results4$agegrp==a,])\r
- if (anova(m, ma, dispersion = max(1, sum(residuals(m, type = "deviance")^2)/df.residual(m)), test="LRT")$`Pr(>Chi)`[2]>p52) {\r
- m <- ma\r
- }\r
- }\r
- }\r
- } else {\r
- if (inherits(m,"try-error")) print("### Could not fit model ###")\r
- if (!m$converge) print("### Model did not converge ###")\r
- if (median(results4[results4$agegrp==a, "deaths"]) == 0) print("### Zero inflated ###")\r
- print("### Simple model with only trend used ###")\r
- f <- paste(c("deaths ~ wk"))\r
- m <- glm(f, quasipoisson, offset = log(N), data=results4[results4$agegrp==a,])\r
- }\r
- print(summary(m, dispersion = max(1, sum(residuals(m, type = "deviance")^2)/df.residual(m))))\r
- # Baseline\r
- results4.B <- results4[results4$agegrp==a,]\r
- for (d in 0:IAlags) {\r
- for (s in min(results4.B$season):max(results4.B$season)) {\r
- results4.B[,paste("d", d, "_IA", s, sep="")] <- 0\r
- }\r
- for (s in c("summer","winter")) {\r
- results4.B[,paste("d", d, "_warm_", s, sep="")] <- 0\r
- results4.B[,paste("d", d, "_cold_", s, sep="")] <- 0\r
- }\r
- }\r
- results4[results4$agegrp==a,]$EB <- exp(predict.glm(m, newdata=results4.B, se.fit=TRUE)$fit)\r
- 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\r
- # IA\r
- results4.IA <- results4[results4$agegrp==a,]\r
- for (s in c("summer","winter")) {\r
- for (d in 0:ETlags) {\r
- results4.IA[,paste("d", d, "_warm_", s, sep="")] <- 0\r
- results4.IA[,paste("d", d, "_cold_", s, sep="")] <- 0\r
- }\r
- }\r
- results4[results4$agegrp==a,]$EIA <- exp(predict.glm(m, newdata=results4.IA, se.fit=TRUE)$fit)\r
- 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\r
- # ET\r
- results4.ET <- results4[results4$agegrp==a,]\r
- for (d in 0:IAlags) {\r
- for (s in sort(unique(results4.ET$season))) {\r
- results4.ET[,paste("d", d, "_IA", s, sep="")] <- 0\r
- }\r
- }\r
- results4[results4$agegrp==a,]$EET <- exp(predict.glm(m, newdata=results4.ET, se.fit=T)$fit)\r
- 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\r
-\r
- results4[results4$agegrp==a,]$Vdeaths <- with(results4[results4$agegrp==a,], EB * max(1, sum(residuals(m, type = "deviance")^2)/df.residual(m))) \r
-}\r
-\r
-# Delete external objects\r
-rm(a, f, fa, s, d, m, ma, results4.B, results4.ET, results4.IA, ptrend, p26, p52, wk)\r
-### END: Estimation ###\r
-\r
-### START: Post estimation ###\r
-# Keep relevant variables #\r
-results4 <- results4[, c("agegrp", "year", "week", "deaths", "N", "IA", "ET", "Vdeaths", "EB", "VlogB", "EIA", "VlogIA", "EET", "VlogET")]\r
-\r
-# Baseline, EIA and EET estimation variances - not on log scale\r
-results4$VB <- with(results4, (exp(VlogB)-1)*exp(2*log(EB) + VlogB))\r
-results4$VIA <- with(results4, (exp(VlogIA)-1)*exp(2*log(EIA) + VlogIA)) \r
-results4$VET <- with(results4, (exp(VlogET)-1)*exp(2*log(EET) + VlogET))\r
-\r
-# Effects of IA and ET\r
-results4$EdIA <- with(results4, EIA - EB)\r
-results4[is.na(results4$EdIA),"EdIA"] <- 0\r
-results4$EdET <- with(results4, EET - EB)\r
-results4[is.na(results4$EdET),"EdET"] <- 0\r
-# Excess relative to baseline\r
-results4$excess <- with(results4, deaths - EB)\r
-# Unexplained excess\r
-results4$uexcess <- with(results4, deaths - (EB + EdIA + EdET))\r
-# Exclude negative IA effects\r
-if (IArest) results4$uexcess <- with(results4, deaths - (EB + pmax(0,EdIA) + EdET))\r
-\r
-# Baseline 2/3 residual confidence intervals\r
-results4$RVB <- with(results4, ((2/3)*(EB^(2/3-1))^2)*Vdeaths + ((2/3)*(EB^(2/3-1))^2)*VB)\r
-results4[with(results4, is.na(RVB) | is.infinite(RVB)), "RVB"] <- 0\r
-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)))\r
-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))\r
-\r
-# EdIA 2/3 residual confidence intervals\r
-results4$RVdIA <- with(results4, ((2/3)*(EB^(2/3-1))^2)*VB + ((2/3)*(EIA^(2/3-1))^2)*VIA)\r
-results4[with(results4, is.na(RVdIA) | is.infinite(RVdIA)), "RVdIA"] <- 0\r
-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))\r
-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)))\r
-\r
-# EdET 2/3 residual confidence intervals\r
-results4$RVdET <- with(results4, ((2/3)*(EB^(2/3-1))^2)*VB + ((2/3)*(EET^(2/3-1))^2)*VET)\r
-results4[with(results4, is.na(RVdET) | is.infinite(RVdET)), "RVdET"] <- 0\r
-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))\r
-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)))\r
-\r
-results4$summer <- with(results4, ifelse(week>=21 & week<=39, year, NA))\r
-results4$winter <- with(results4, ifelse(week<=20 | week>=40, year - (week<=20), NA))\r
-\r
-# Order\r
-results4 <- results4[with(results4, order(agegrp, year, week)),]\r
-\r
-cummulate <- function(x, pos=FALSE) {\r
- # x = variable, EdIA\r
- y <- x[,1]\r
- if (pos) y <- x[,1] * (x[,2] > 0)\r
- return(cumsum(y))\r
-}\r
-cCI <- function(x) {\r
- # x = cEB, cVB, cE*, cV*\r
- colnames(x) <- c("cEB", "cVB", "cE", "cV")\r
- x$cEd <- x$cE - x$cEB\r
- x$cRVd <- ((2/3)*(x$cEB^(2/3-1))^2)*x$cVB + ((2/3)*(x$cE^(2/3-1))^2)*x$cV\r
- x[with(x, is.na(cRVd) | is.infinite(cRVd)), "cRVd"] <- 0\r
- 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))\r
- 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))\r
- return(x[, c("cEd", "cCI_95L", "cCI_95U")])\r
-}\r
-\r
-# EdIA and EdET residual variances\r
-for (s in c("summer", "winter", "year")) {\r
- nn <- !is.na(results4[,s]) # logical vector indicating season\r
- for (a in 0:4) {\r
- ag <- (results4$agegrp==a)\r
- for (v in c("excess","uexcess")) {\r
- results4[nn & ag ,paste0("c", v, "_", s)] <- NA\r
- for (sy in sort(unique(results4[,s]))) {\r
- results4[(results4[,s]==sy) & nn & ag, paste0("c", v, "_", s)] <-\r
- cummulate(results4[(results4[,s]==sy) & nn & ag, c(v, "EdIA")], FALSE)\r
- }\r
- }\r
- for (v in c("B","IA")) {\r
- results4[nn & ag ,paste0("cE", v, "_", s)] <- NA\r
- results4[nn & ag ,paste0("cV", v, "_", s)] <- NA\r
- for (sy in sort(unique(results4[,s]))) {\r
- results4[(results4[,s]==sy) & nn & ag, paste0("cE", v, "_", s)] <-\r
- cummulate(results4[(results4[,s]==sy) & nn & ag, c(paste0("E",v),"EdIA")], IArest)\r
- results4[(results4[,s]==sy) & nn & ag, paste0("cV", v, "_", s)] <-\r
- cummulate(results4[(results4[,s]==sy) & nn & ag, c(paste0("V",v),"EdIA")], IArest)\r
- }\r
- }\r
- results4[nn & ag, paste0("cEd", v, "_", s)] <- NA\r
- results4[nn & ag, paste0("cEd", v, "_", s, "_95L")] <- NA\r
- results4[nn & ag, paste0("cEd", v, "_", s, "_95U")] <- NA\r
- results4[nn & ag, c(paste0("cEd", v, "_", s), paste0("cEd", v, "_", s, "_95L"), paste0("cEd", v,"_", s, "_95U"))] <-\r
- cCI(results4[nn & ag, c(paste0("cEB_", s), paste0("cVB_", s), paste0("cE",v,"_", s), paste0("cV",v,"_",s))])\r
- for (v in c("B","ET")) {\r
- results4[nn & ag ,paste0("cE", v, "_", s)] <- NA\r
- results4[nn & ag ,paste0("cV", v, "_", s)] <- NA\r
- for (sy in sort(unique(results4[,s]))) {\r
- results4[(results4[,s]==sy) & nn & ag, paste0("cE", v, "_", s)] <-\r
- cummulate(results4[(results4[,s]==sy) & nn & ag, c(paste0("E",v),"EdIA")], FALSE)\r
- results4[(results4[,s]==sy) & nn & ag, paste0("cV", v, "_", s)] <-\r
- cummulate(results4[(results4[,s]==sy) & nn & ag, c(paste0("V",v),"EdIA")], FALSE)\r
- }\r
- }\r
- results4[nn & ag ,paste0("cEd", v, "_", s)] <- NA\r
- results4[nn & ag ,paste0("cEd", v, "_", s, "_95L")] <- NA\r
- results4[nn & ag ,paste0("cEd", v, "_", s, "_95U")] <- NA\r
- results4[nn & ag ,c(paste0("cEd", v, "_", s), paste0("cEd", v, "_", s, "_95L"), paste0("cEd", v, "_", s, "_95U"))] <-\r
- cCI(results4[nn & ag , c(paste0("cEB_", s), paste0("cVB_", s), paste0("cE",v,"_", s), paste0("cV",v,"_",s))])\r
- }\r
-}\r
-\r
-results4$country <- country\r
-results4$IArestricted <- as.integer(IArest)\r
-results4 <- results4[, c("country", "IArestricted", "agegrp", "year", "week", "deaths", "Vdeaths", "N", "IA", "ET",\r
- "EB", "EB_95L", "EB_95U", "VB", \r
- "EIA", "VIA", "EET", "VET", \r
- "EdIA", "EdIA_95L", "EdIA_95U", \r
- "EdET", "EdET_95L", "EdET_95U", \r
- "cexcess_year", "cuexcess_year", \r
- "cEdIA_year", "cEdIA_year_95L", "cEdIA_year_95U", \r
- "cEdET_year", "cEdET_year_95L", "cEdET_year_95U", \r
- "summer", "cexcess_summer", "cuexcess_summer", \r
- "cEdIA_summer", "cEdIA_summer_95L", "cEdIA_summer_95U", \r
- "cEdET_summer", "cEdET_summer_95L", "cEdET_summer_95U", \r
- "winter", "cexcess_winter", "cuexcess_winter", \r
- "cEdIA_winter", "cEdIA_winter_95L", "cEdIA_winter_95U", \r
- "cEdET_winter", "cEdET_winter_95L", "cEdET_winter_95U")]\r
-\r
-if (IArest) { \r
- write.table(results4, file=paste0(outdir,"/",country,"_output_v4_IArestricted.txt"), \r
- row.names = FALSE, quote = FALSE, sep =";", dec=".", na="") \r
-} else { \r
- write.table(results4, file=paste0(outdir,"/",country,"_output_v4.txt"), \r
- row.names = FALSE, quote = FALSE, sep =";", dec=".", na="")\r
-}\r
-\r
-rm(a, s, sy, v, nn, ag, cummulate, cCI, results4)\r
-### END: Post estimation ###\r
projects / covid19.git / blobdiff