################################################################################ # Code for Gehringer, Rode, Schomaker # # The effect of load shedding on pediatric hospital admissions in South Africa # ################################################################################ ################################################################################ # To obtain the data to run the code, and reproduce the results, please contact# # Christian Gehringer # # Michael Schomaker # # # # Note that analyses of this data can only be published after approval from # # Red Cross Hospital # ################################################################################ ################################################################################ # Main variables: # # Y = number admissions = "admissions" # # A = load shedding on respective day = "loads3" # # Abar2 = load shedding on same or previous 2 days = "loads3_last2" # # L= weather summaries = sunshine_ap, wind_tot, rain_tot, temp_tot, # # humidity_tot, pressure_tot # # and t (time), payweek, month, weekend # # "_lag1" refers to any Y,A,L measured the day before, similarly "_lag2" etc. # # "_last2" refers to any binary Y,A,L occurring on the same or prev. 2 days # # "ICD_" variables are number admissions per ICD 10 group # # "area_" variables indicate whether load shedding occurred in the respective # # area, i.e. 1-16 # ################################################################################ ########################################################### # Part 1: Prelimiaries - Load libraries, define functions # ########################################################### # Specify (YOUR) Working Directory setwd("C://Users//01429265//Dropbox//Documents//Projects_Other//Loadshedding//Results_reproduce") # packages library(ggplot2) # graphs library(gridExtra) # graphs library(MuMIn) # Frequentist Model Averaging [FMA] require(parallel) # parallel computing [used for FMA] library(MASS) # many, for example negative binomial model library(survival) library(tscount) # INGARCH model library(tmle) # TMLE: causal inference library(SuperLearner) # learners for TMLE library(MAMI) # from http://mami.r-forge.r-project.org/ # - additional learners for ltmle/SuperLearner library(xlsx) # save output as Excel file library(xtable) # save output as Latex table library(vcd) # for Cramer's V [association matrix for areas] library(corrplot) # plot association matrix # SESSION INFO (i.e. version of packages used) - sessionInfo() # We used R 3.4.4 # #attached base packages: # [1] grid parallel stats graphics grDevices utils datasets methods base # #other attached packages: # [1] leaflet_1.1.0 rgdal_1.2-18 sp_1.2-7 corrplot_0.84 vcd_1.4-4 xtable_1.8-2 xlsx_0.5.7 xlsxjars_0.6.1 rJava_0.9-9 MAMI_0.9.10 tmle_1.2.0-5 # [12] SuperLearner_2.0-22 nnls_1.4 tscount_1.4.1 survival_2.41-3 MASS_7.3-49 MuMIn_1.40.1 gridExtra_2.3 ggplot2_2.2.1 #loaded via a namespace (and not attached): # [1] Rcpp_0.12.14 mvtnorm_1.0-6 lattice_0.20-35 zoo_1.8-1 glmnet_2.0-13 digest_0.6.13 lmtest_0.9-35 foreach_1.4.5 mime_0.5 R6_2.2.2 plyr_1.8.4 stats4_3.4.4 # [13] pcaPP_1.9-72.1 Amelia_1.7.4 pillar_1.1.0 rlang_0.1.6 lazyeval_0.2.1 minqa_1.2.4 nloptr_1.0.4 rpart_4.1-13 Matrix_1.2-12 splines_3.4.4 lme4_1.1-15 foreign_0.8-69 # [25] htmlwidgets_1.0 munsell_0.4.3 shiny_1.0.5 httpuv_1.3.6.2 compiler_3.4.4 htmltools_0.3.6 nnet_7.3-12 tibble_1.4.1 quadprog_1.5-5 codetools_0.2-15 BMA_3.18.7 rrcov_1.4-3 # [37] leaps_3.0 nlme_3.1-131.1 ltsa_1.4.6 gtable_0.2.0 magrittr_1.5 scales_0.5.0 mice_2.46.0 robustbase_0.92-8 boot_1.3-20 iterators_1.0.9 DEoptimR_1.0-8 crosstalk_1.0.0 # [49] inline_0.3.14 colorspace_1.3-2 cluster_2.0.6 # seed set.seed(24121980) # Define some refined super learners for Sections 4ff. below here SL.gbm2 <- function(Y, X, newX, family, obsWeights, gbm.trees = 10000, interaction.depth = 2, shrinkage = 0.001, ...) { gbm.model <- as.formula(paste("Y~", paste(colnames(X), collapse = "+"))) if (family$family == "gaussian") { fit.gbm <- gbm::gbm(formula = gbm.model, data = X, distribution = "gaussian", n.trees = gbm.trees, interaction.depth = interaction.depth, shrinkage = shrinkage, cv.folds = 5, keep.data = TRUE, bag.fraction=0.75, weights = obsWeights, verbose = FALSE) } if (family$family == "binomial" & all(Y%in%c(0,1))) { fit.gbm <- gbm::gbm(formula = gbm.model, data = X, distribution = "bernoulli", n.trees = gbm.trees, interaction.depth = interaction.depth, shrinkage = shrinkage, cv.folds = 5, keep.data = TRUE, bag.fraction=0.75, weights = obsWeights, verbose = FALSE) } if (family$family == "binomial" & all(Y%in%c(0,1))==FALSE) { fit.gbm <- gbm::gbm(formula = gbm.model, data = X, distribution = "gaussian", n.trees = gbm.trees, interaction.depth = interaction.depth, shrinkage = shrinkage, cv.folds = 5, keep.data = TRUE, bag.fraction=0.75, weights = obsWeights, verbose = FALSE) } best.iter <- gbm::gbm.perf(fit.gbm, method = "cv", plot.it = FALSE) pred <- predict(fit.gbm, newdata = newX, best.iter, type = "response") fit <- list(object = fit.gbm, n.trees = best.iter) out <- list(pred = pred, fit = fit) class(out$fit) <- c("SL.gbm") return(out) } SL.earth2 <- function(Y, X, newX, family, obsWeights, id, degree = 2, penalty = 3, nk = max(21, 2 * ncol(X) + 1), pmethod = "backward", nfold = 0, ncross = 1, minspan = 0, endspan = 0, ...) { if (family$family == "gaussian") { fit.earth <- earth::earth(x = X, y = Y, degree = degree, nk = nk, penalty = penalty, pmethod = pmethod, nfold = nfold, ncross = ncross, minspan = minspan, endspan = endspan) } if (family$family == "binomial" & all(Y%in%c(0,1))){ fit.earth <- earth::earth(x = X, y = Y, degree = degree, nk = nk, penalty = penalty, pmethod = pmethod, nfold = nfold, ncross = ncross, minspan = minspan, endspan = endspan, glm = list(family = binomial)) } if (family$family == "binomial" & all(Y%in%c(0,1))==FALSE) { fit.earth <- earth::earth(x = X, y = Y, degree = degree, nk = nk, penalty = penalty, pmethod = pmethod, nfold = nfold, ncross = ncross, minspan = minspan, endspan = endspan) } pred <- predict(fit.earth, newdata = newX, type = "response") fit <- list(object = fit.earth) out <- list(pred = pred, fit = fit) class(out$fit) <- c("SL.earth") return(out) } screen.glmnet3 <- function(Y, X, family, alpha = 1, minscreen = 2, pw=F, maxtries=4, nfolds = 10, nlambda = 200, ...){ if(family$family == "binomial" & all(Y%in%c(0,1))==FALSE){myfamily<- "gaussian"}else{myfamily<-family$family} if (!is.matrix(X)) {X <- try(model.matrix(~-1 + ., data=X),silent=T)} successfulfit <- FALSE fitCV <- try(glmnet::cv.glmnet(x = X, y = Y, lambda = NULL, type.measure = "deviance", nfolds = nfolds, family = myfamily, alpha = alpha, nlambda = nlambda, keep=T),silent=T) if(class(fitCV)=="try-error"){ i <- 2 while(successfulfit==FALSE & i<=maxtries){ if(pw==T){cat(paste("glmnet failed, new try #",i,"\n"))} fitCV <- try(glmnet::cv.glmnet(x = X, y = Y, lambda = NULL, type.measure = "deviance", nfolds = nfolds, family = myfamily, alpha = alpha, nfolds=4, nlambda = nlambda*(i+3), keep=T, foldid=sample(fitCV$foldid)),silent=T) i <- i+1 if(class(fitCV)=="try-error"){successfulfit <- FALSE }else{successfulfit <- TRUE} } }else{successfulfit <- TRUE} whichVariable <- NULL if(successfulfit==TRUE){ si <- try(abs((max(fitCV$nzero)- (dim(X)[2]/2))-fitCV$nzero),silent=T) whichVariable2 <- try((as.numeric(coef(fitCV$glmnet.fit, s = fitCV$lambda.min))[-1] != 0),silent=T) whichVariable3 <- try(as.numeric(glmnet::coef.glmnet(fitCV, s = fitCV$lambda[si==min(si)][1]))[-1] != 0,silent=T) if(sum(whichVariable2)>1){whichVariable<-whichVariable2}else{if(sum(whichVariable3)>1 & sum(whichVariable3)<(dim(X)[2]/2)){whichVariable<-whichVariable3}} } if(is.null(whichVariable)){ whichVariable<-screen.cramersv(Y,X) if(pw==T){cat("Lasso failed and screening was based on Cramer's V\n")}} if(pw==T){cat(paste("Number of included variables:",sum(whichVariable),"\n"))} return(whichVariable) } screen.cramersv <- function(Y, X, nscreen=6, cts.num=10, ...){ var_cont <- apply(X, 2, function(x) (length(unique(x)) > cts.num)) cvdata <- X cutf <- function(x){cut(x, unique(quantile(x, prob = c(0, 0.2, 0.4, 0.6, 0.8, 1))),include.lowest=T)} if(sum(var_cont>0)){cvdata[,var_cont] <- apply(cvdata[,var_cont],2,cutf)} cvf <- function(xx){assocstats(table(Y,xx))$cramer} whichVariable <- order(apply(cvdata,2,cvf),decreasing=T)<=nscreen return(whichVariable) } ######################## ptm2 <- proc.time() ######################## # Part 2: Descriptives # ######################## adm <- read.csv('final_data_all.csv') adm$date <- as.Date(adm$date, "%Y-%m-%d") attach(adm) mycolors <- c("red","black") # Figure 2 pdf(file=paste(getwd(),"/Figure2.pdf",sep=""), width=9) ggplot(adm, aes(admissions, colour=Status, fill=Status)) + geom_density(alpha=0.15) +theme_bw() + facet_grid(. ~ weekend2) + scale_colour_manual(values = mycolors)+ scale_fill_manual(values = mycolors) + scale_x_continuous("Number of admissions", limits=c(0,115)) + scale_y_continuous(" ") + theme(axis.title.x = element_text(size=14), axis.text.x = element_text(size=14),axis.title.y = element_text(size=14, angle = 90), axis.text.y = element_text(size=14), legend.text = element_text(size=13), legend.title = element_text(size=13, face = "bold", hjust = 0),legend.position = "bottom", legend.key.width = unit(1.75, "cm")) dev.off() # eFigure 1 ggtemp <- as.data.frame(cbind(c(temp_ap,temp_kb,temp_obz,temp_mo,temp_yc,temp_tot))) colnames(ggtemp) <- c("Temperature") ggtemp$Date <- rep(adm$date,6) ggtemp$Station <- c(rep("Airport",365),rep("Kirstenbosch",365),rep("Observatory",365),rep("Reservoir",365),rep("Yacht Club",365),rep("Total",365)) gghum <- as.data.frame(cbind(c(humidity_ap,humidity_kb,humidity_obz,humidity_mo,humidity_yc,humidity_tot))) colnames(gghum) <- c("Humidity") gghum$Date <- rep(adm$date,6) gghum$Station <- c(rep("Airport",365),rep("Kirstenbosch",365),rep("Observatory",365),rep("Reservoir",365),rep("Yacht Club",365),rep("Total",365)) ggpressure <- as.data.frame(cbind(c(pressure_ap,presssure_kb,pressure_obz,pressure_mo,pressure_yc,pressure_tot))) colnames(ggpressure) <- c("Pressure") ggpressure$Date <- rep(adm$date,6) ggpressure$Station <- c(rep("Airport",365),rep("Kirstenbosch",365),rep("Observatory",365),rep("Reservoir",365),rep("Yacht Club",365),rep("Total",365)) ggrain <- as.data.frame(cbind(c(rain_ap,rain_kb,rain_obz,rain_mo,rain_yc,rain_tot))) colnames(ggrain) <- c("Rainfall") ggrain$Date <- rep(adm$date,6) ggrain$Station <- c(rep("Airport",365),rep("Kirstenbosch",365),rep("Observatory",365),rep("Reservoir",365),rep("Yacht Club",365),rep("Total",365)) ggwind <- as.data.frame(cbind(c(wind_ap,wind_kb,wind_obz,wind_mo,wind_yc,wind_tot))) colnames(ggwind) <- c("Wind") ggwind$Date <- rep(adm$date,6) ggwind$Station <- c(rep("Airport",365),rep("Kirstenbosch",365),rep("Observatory",365),rep("Reservoir",365),rep("Yacht Club",365),rep("Total",365)) g1 <- ggplot(ggtemp, aes(Date,Temperature, colour=Station)) + geom_smooth(se=F) +theme_bw() + scale_colour_brewer(palette = "Set1") + theme(axis.title.x = element_text(size=14), axis.text.x = element_text(size=14),axis.title.y = element_text(size=14, angle = 90), axis.text.y = element_text(size=14), legend.text = element_text(size=13), legend.title = element_text(size=13, face = "bold", hjust = 0),legend.position = "bottom", legend.key.width = unit(1.75, "cm")) g2 <- ggplot(ggpressure, aes(Date,Pressure, colour=Station)) + geom_smooth(se=F) +theme_bw() + scale_colour_brewer(palette = "Set1") + theme(axis.title.x = element_text(size=14), axis.text.x = element_text(size=14),axis.title.y = element_text(size=14, angle = 90), axis.text.y = element_text(size=14), legend.text = element_text(size=13), legend.title = element_text(size=13, face = "bold", hjust = 0),legend.position = "bottom", legend.key.width = unit(1.75, "cm")) g3 <- ggplot(gghum, aes(Date,Humidity, colour=Station)) + geom_smooth(se=F) +theme_bw() + scale_colour_brewer(palette = "Set1") +theme(axis.title.x = element_text(size=14), axis.text.x = element_text(size=14),axis.title.y = element_text(size=14, angle = 90), axis.text.y = element_text(size=14), legend.text = element_text(size=13), legend.title = element_text(size=13, face = "bold", hjust = 0),legend.position = "bottom", legend.key.width = unit(1.75, "cm")) g4 <- ggplot(ggrain, aes(Date,Rainfall, colour=Station)) + geom_smooth(se=F) +theme_bw() + scale_colour_brewer(palette = "Set1") + theme(axis.title.x = element_text(size=14), axis.text.x = element_text(size=14),axis.title.y = element_text(size=14, angle = 90), axis.text.y = element_text(size=14), legend.text = element_text(size=13), legend.title = element_text(size=13, face = "bold", hjust = 0),legend.position = "bottom", legend.key.width = unit(1.75, "cm")) g5 <- ggplot(ggwind, aes(Date,Wind, colour=Station)) + geom_smooth(se=F) +theme_bw() + scale_colour_brewer(palette = "Set1") + theme(axis.title.x = element_text(size=14), axis.text.x = element_text(size=14),axis.title.y = element_text(size=14, angle = 90), axis.text.y = element_text(size=14), legend.text = element_text(size=13), legend.title = element_text(size=13, face = "bold", hjust = 0),legend.position = "bottom", legend.key.width = unit(1.75, "cm")) g_legend<-function(a.gplot){ tmp <- ggplot_gtable(ggplot_build(a.gplot)) leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box") legend <- tmp$grobs[[leg]] return(legend)} mylegend <- g_legend(g1) eF1 <- grid.arrange(arrangeGrob(g1 + theme(legend.position="none"), g2 + theme(legend.position="none"), g3 + theme(legend.position="none"), g4 + theme(legend.position="none"), g5 + theme(legend.position="none"), mylegend, nrow=3)) # save manually for nicer look pdf(file=paste(getwd(),"/eFigure1.pdf",sep=""), width=12, height=12) plot(eF1) dev.off() # eFigure5 areas <- adm[,c("area1","area2","area3","area4","area5","area6","area7","area8","area9","area10","area11","area12","area13","area14","area15","area16")] cmat <- matrix(NA,ncol=16,nrow=16) # association matrix measured by Cramer's V for(i in 1:16){for(j in 1:16){ if(i!=j){ cmat[i,j] <- assocstats(table(areas[,i],areas[,j]))$cramer cmat[j,i] <- assocstats(table(areas[,i],areas[,j]))$cramer } }} diag(cmat) <- 1 areanames <- c("Belville","Maitland/Milnerton","Somerset West/Strand","Mitchell's Plain/Philippi Farmlands","Hanover Park/Lansdown/Newlands","Durbanville/Tyger Valley", "City Bowl/Camps Bay","Southern Peninsula","Pinelands/Langa/Industrial","Brackenfell/Kuils River", "Plumstead/Ottery/Hout Bay/Constantia","Athlone/Manenberg/Gugulethu","Goodwood/Parow/Plattekloof","Atlantis/Melkbosstrand","Obz/Rondebosch/Brooklyn","Philippi/Retreat") colnames(cmat) <- areanames rownames(cmat) <- areanames col3 <- colorRampPalette(c("green", "brown", "white", "blue", "orange", "red")) pdf(file=paste(getwd(),"/eFigure5.pdf", sep=""), width=10) corrplot::corrplot(cmat, diag=F, col = col3(200),tl.col = "black", cl.lim = c(0,1),type="upper") dev.off() detach(adm) ptm3 <- proc.time() ################## # Part 3: Models # ################## adm <- read.csv('final_data_all.csv') adm$date <- as.Date(adm$date, "%Y-%m-%d") ########################### # i) pre considerations # ########################### # Model family m_adm0 <- lm(admissions ~ loads3, data=adm) m_adm1 <- glm(admissions ~ loads3, family="poisson", data=adm) m_adm2 <- glm(admissions ~ loads3, family="quasipoisson", data=adm) summary(m_adm0) #plot(m_adm0) # assumptions violated (normality of residuals) summary(m_adm1) # assumptions violated (overdispersion: deviance/df >1) summary(m_adm2) # o.k. library(mgcv) # Weekend hypothesis m_we2 <- gam(admissions ~ weekend,family=quasipoisson,data=adm) summary(m_we2) m_we2.1 <- gam(loads3 ~ weekend,family=binomial,data=adm) summary(m_we2.1) ####################################################### # ii) selection of a priori considered variables # # with model averaging # ####################################################### attach(adm) # "wrapper", so that AIC can be used by MuMIn for quasi-likelihood methods, i.e. family=quasipoisson dfun <- function(object){ with(object,sum((weights * residuals^2)[weights > 0])/df.residual) } x.quasipoisson <- function(...) { res <- quasipoisson(...) res$aic <- poisson(...)$aic res } # parallelization to speed up calculation # note that this can take long when the number of cores in your computer is small # on a computer with 8 cores (of which we use 7) model averaging took about 1 hour clust <- makeCluster(getOption("cl.cores", max(detectCores()-1,1)), type = "PSOCK") clusterExport(clust, list("adm", "dfun","x.quasipoisson"), envir=environment()) # a) how many lags of y? global_model.i <- glm(admissions ~ admissions_lag1+admissions_lag2+admissions_lag3+admissions_lag4+admissions_lag5+admissions_lag6+admissions_lag7+admissions_lag8+admissions_lag9+admissions_lag10 +admissions_lag11+admissions_lag12+admissions_lag13+admissions_lag14+admissions_lag21+admissions_lag28 +I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)) +loads3_last2 + month+ sunshine_ap+wind_tot+rain_tot+rain_tot_lag1+temp_tot_lag1+humidity_tot_lag1+wind_tot_lag2+payweek, family="quasipoisson",na.action=na.fail, data=adm) global_model2.i <- glm(admissions ~ admissions_lag1+admissions_lag2+admissions_lag3+admissions_lag4+admissions_lag5+admissions_lag6+admissions_lag7+admissions_lag8+admissions_lag9+admissions_lag10+admissions_lag28 +admissions_lag11+admissions_lag12+admissions_lag13+admissions_lag14+admissions_lag21+admissions_lag28 +I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)) + loads3_last2 + month+ sunshine_ap+wind_tot+rain_tot+rain_tot_lag1+temp_tot_lag1+humidity_tot_lag1+wind_tot_lag2+payweek, family="poisson",na.action=na.fail, data=adm) global_model3.i <- update(global_model2.i,family="x.quasipoisson",na.action=na.fail) m_adm_QAIC.i <- pdredge(global_model3.i,rank="QAIC",chat=dfun(global_model2.i),cluster=clust,fixed=~loads3_last2 + month+ sunshine_ap+wind_tot+rain_tot+rain_tot_lag1+temp_tot_lag1+humidity_tot_lag1+wind_tot_lag2+payweek+I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)) ) i1 <- importance(model.avg(m_adm_QAIC.i)) rm(list=c("m_adm_QAIC.i")) i1 # variable importance (VI) for admission lags # b) seasonal component of t? global_model.ii <- glm(admissions ~ admissions_lag1+admissions_lag3+admissions_lag7+admissions_lag9 +I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)) +I(sin(pi*2*3*t/7)) +I(cos(pi*2*3*t/7)) +I(sin(pi*2*4*t/7)) +I(cos(pi*2*4*t/7)) +loads3_last2 + month+ sunshine_ap+wind_tot+rain_tot+rain_tot_lag1+temp_tot_lag1+humidity_tot_lag1+wind_tot_lag2+payweek, family="quasipoisson",na.action=na.fail, data=adm) global_model2.ii <- glm(admissions ~ admissions_lag1+admissions_lag3+admissions_lag7+admissions_lag9 +I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)) +I(sin(pi*2*3*t/7)) +I(cos(pi*2*3*t/7)) +I(sin(pi*2*4*t/7)) +I(cos(pi*2*4*t/7)) + loads3_last2 + month+ sunshine_ap+wind_tot+rain_tot+rain_tot_lag1+temp_tot_lag1+humidity_tot_lag1+wind_tot_lag2+payweek, family="poisson",na.action=na.fail, data=adm) global_model3.ii <- update(global_model2.ii,family="x.quasipoisson",na.action=na.fail) m_adm_QAIC.ii <- pdredge(global_model3.ii,rank="QAIC",chat=dfun(global_model2.ii),cluster=clust,fixed=~admissions_lag1+admissions_lag3+admissions_lag7+admissions_lag9 + loads3_last2 + month+ sunshine_ap+wind_tot+rain_tot+rain_tot_lag1+temp_tot_lag1+humidity_tot_lag1+wind_tot_lag2+payweek) i2 <- importance(model.avg(m_adm_QAIC.ii)) rm(list=c("m_adm_QAIC.ii")) i2 # variable importance (VI) for seasonal effect; only for amplitudes 1 and 2 both sine and cosine are selected # c) lags of load shedding global_model.iii <- glm(admissions ~ loads3 + loads3_lag1 + loads3_lag2 + loads3_lag3 + loads3_lag4 + loads3_lag5 + loads3_lag6 + loads3_lag7 + month+ sunshine_ap+wind_tot+rain_tot+rain_tot_lag1+temp_tot_lag1+humidity_tot_lag1+wind_tot_lag2+payweek ,family="quasipoisson",na.action=na.fail, data=adm) global_model2.iii <- glm(admissions ~ loads3 + loads3_lag1 + loads3_lag2 + loads3_lag3 + loads3_lag4 + loads3_lag5 + loads3_lag6 + loads3_lag7 + month+ sunshine_ap+wind_tot+rain_tot+rain_tot_lag1+temp_tot_lag1+humidity_tot_lag1+wind_tot_lag2+payweek ,family="poisson",na.action=na.fail, data=adm) global_model3.iii <- update(global_model2.iii,family="x.quasipoisson",na.action=na.fail) m_adm_QAIC.iii <- pdredge(global_model3.iii,rank="QAIC",chat=dfun(global_model2.iii),cluster=clust ,fixed=~month+ sunshine_ap+wind_tot+rain_tot+rain_tot_lag1+temp_tot_lag1+humidity_tot_lag1+wind_tot_lag2+payweek ) i3 <- importance(model.avg(m_adm_QAIC.iii)) rm(list=c("m_adm_QAIC.iii")) i3 # variable importance (VI) for lags of load shedding # d) weather indicators in relation to admissions [pressure lag2 and humidity lag2 omitted for computational feasability (based on initial screening)] global_model.iv <- glm(admissions ~ admissions_lag1+admissions_lag3+admissions_lag7+admissions_lag9 +I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)) +loads3_last2 + month+ payweek + wind_tot+pressure_tot+rain_tot+temp_tot+humidity_tot+sunshine_ap +wind_tot_lag1+pressure_tot_lag1+rain_tot_lag1+temp_tot_lag1+humidity_tot_lag1+sunshine_ap_lag1 +wind_tot_lag2+rain_tot_lag2+temp_tot_lag2+sunshine_ap_lag2 , family="quasipoisson",na.action=na.fail, data=adm) global_model2.iv <- glm(admissions ~ admissions_lag1+admissions_lag3+admissions_lag7+admissions_lag9 +I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)) +loads3_last2 + month+ payweek + wind_tot+pressure_tot+rain_tot+temp_tot+humidity_tot+sunshine_ap +wind_tot_lag1+pressure_tot_lag1+rain_tot_lag1+temp_tot_lag1+humidity_tot_lag1+sunshine_ap_lag1 +wind_tot_lag2+rain_tot_lag2+temp_tot_lag2+sunshine_ap_lag2 , family="poisson",na.action=na.fail, data=adm) global_model3.iv <- update(global_model2.iv,family="x.quasipoisson",na.action=na.fail) m_adm_QAIC.iv <- pdredge(global_model3.iv,rank="QAIC",chat=dfun(global_model2.iv),cluster=clust,fixed=~admissions_lag1+admissions_lag3+admissions_lag7+admissions_lag9 + loads3_last2 + month + payweek + I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)) ) i4 <- importance(model.avg(m_adm_QAIC.iv)) rm(list=c("m_adm_QAIC.iv")) i4 # variable importance (VI) for lags of weather indicators (w.r.t outcome "admissions") # e) weather indicators in relation to load shedding global_model.v <- glm(loads ~ month +loads3_lag1 + loads3_lag2 + payweek + wind_tot+pressure_tot+rain_tot+temp_tot+humidity_tot+sunshine_ap + wind_tot_lag1+pressure_tot_lag1+rain_tot_lag1+temp_tot_lag1+humidity_tot_lag1+sunshine_ap_lag1 + wind_tot_lag2+pressure_tot_lag2+rain_tot_lag2+temp_tot_lag2+humidity_tot_lag2+sunshine_ap_lag2 , family="binomial",na.action=na.fail, data=adm) m_adm_QAIC.v <- pdredge(global_model.v,rank="AIC",cluster=clust,fixed=~month +loads3_lag1 + loads3_lag2 + payweek) i5 <- importance(model.avg(m_adm_QAIC.v)) rm(list=c("m_adm_QAIC.v")) i5 # variable importance (VI) for lags of weather indicators (w.r.t outcome "load shedding") threshold <- 0.5 i4[i4>threshold & i4<0.99] # VI >0.5 weather indicators i5[i5>threshold & i5<0.99] # VI >0.5 weather indicators stopCluster(clust) # eFigure2 (looking visually at the selected weather indicators in terms of confounding) m_we <- gam(loads3 ~ +s(sunshine_ap)+s(wind_tot)+s(rain_tot)+s(rain_tot_lag1)+s(temp_tot_lag1)+s(humidity_tot_lag1)+s(wind_tot_lag2)+month,family=binomial,data=adm) summary(m_we) m_we3 <- gam(admissions ~ +s(sunshine_ap)+s(wind_tot)+s(rain_tot)+s(rain_tot_lag1)+s(temp_tot_lag1)+s(humidity_tot_lag1)+s(wind_tot_lag2) +s(admissions_lag1)+s(admissions_lag3)+s(admissions_lag7)+s(admissions_lag9) +I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)) ,family=quasipoisson,data=adm) summary(m_we3) pdf(file=paste(getwd(),"/eFigure2.pdf",sep=""), width=10,height=12) par(mfrow=c(4,3)) plot(m_we3, select=1, scale=0, trans=exp, xlab="sunshine (hours)", ylab="IRR for number of admissions", cex.axis=1.5,cex.lab=1.5,cex=1.5,cex.main=1.5) plot(m_we3, select=2, scale=0, trans=exp, xlab="wind speed", ylab="IRR for number of admissions", cex.axis=1.5,cex.lab=1.5,cex=1.5,cex.main=1.5) plot(m_we3, select=3, scale=0, trans=exp, xlab="rainfall", ylab="IRR for number of admissions", cex.axis=1.5,cex.lab=1.5,cex=1.5,cex.main=1.5) plot(m_we3, select=4, scale=0, trans=exp, xlab="rainfall (lag 1)", ylab="IRR for number of admissions", cex.axis=1.5,cex.lab=1.5,cex=1.5,cex.main=1.5) plot(m_we3, select=5, scale=0, trans=exp, xlab="temperature (lag 1)", ylab="IRR for number of admissions", cex.axis=1.5,cex.lab=1.5,cex=1.5,cex.main=1.5) plot(m_we3, select=6, scale=0, trans=exp, xlab="humidity (lag 1)", ylab="IRR for number of admissions", cex.axis=1.5,cex.lab=1.5,cex=1.5,cex.main=1.5) plot(m_we, select=1, scale=0, trans=plogis, xlab="sunshine (hours)", ylab="Probability of load shedding", cex.axis=1.5,cex.lab=1.5,cex=1.5,cex.main=1.5) plot(m_we, select=2, scale=0, trans=plogis, xlab="wind speed", ylab="Probability of load shedding", cex.axis=1.5,cex.lab=1.5,cex=1.5,cex.main=1.5) plot(m_we, select=3, scale=0, trans=plogis, xlab="rainfall", ylab="Probability of load shedding", cex.axis=1.5,cex.lab=1.5,cex=1.5,cex.main=1.5) plot(m_we, select=4, scale=0, trans=plogis, xlab="rainfall (lag 1)", ylab="Probability of load shedding", cex.axis=1.5,cex.lab=1.5,cex=1.5,cex.main=1.5) plot(m_we, select=5, scale=0, trans=plogis, xlab="temperature (lag 1)", ylab="Probability of load shedding", cex.axis=1.5,cex.lab=1.5,cex=1.5,cex.main=1.5) plot(m_we, select=6, scale=0, trans=plogis, xlab="humidity (lag 1)", ylab="Probability of load shedding", cex.axis=1.5,cex.lab=1.5,cex=1.5,cex.main=1.5) dev.off() ######################################## # iii) Main models (parts of Table 1) # ######################################## mainformula <- as.formula(admissions ~ loads3_last2 +s(sunshine_ap)+s(wind_tot)+s(rain_tot)+s(rain_tot_lag1)+s(temp_tot_lag1)+s(humidity_tot_lag1)+s(wind_tot_lag2) +month+payweek +s(admissions_lag1)+s(admissions_lag3)+s(admissions_lag7)+s(admissions_lag9) +I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7))) covariates <- "loads3_last2 +s(sunshine_ap)+s(wind_tot)+s(rain_tot)+s(rain_tot_lag1)+s(temp_tot_lag1)+s(humidity_tot_lag1)+s(wind_tot_lag2)+month+payweek +s(admissions_lag1)+s(admissions_lag3)+s(admissions_lag7)+s(admissions_lag9)+I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7))" covariates2 <- "loads3_last2 + sunshine_ap+wind_tot+rain_tot+rain_tot_lag1+temp_tot_lag1+humidity_tot_lag1+wind_tot_lag2+month+payweek +admissions_lag1+admissions_lag3+admissions_lag7+admissions_lag9+I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7))" # a) loadshedding on same day or preceding 2 days m_adm6 <- gam(mainformula, family="quasipoisson") summary(m_adm6) # Odds Ratio and cofidence Interval exp(m_adm6$coefficients["loads3_last2"]) exp(m_adm6$coefficients["loads3_last2"] - qnorm(0.975)*summary(m_adm6)[[22]]["loads3_last2",2]) exp(m_adm6$coefficients["loads3_last2"] + qnorm(0.975)*summary(m_adm6)[[22]]["loads3_last2",2]) # b) loadshedding on same day or preceding 2 days + Interaction m_adm6.1 <- gam(admissions ~ loads3_last2*weekend +s(sunshine_ap)+s(wind_tot)+s(rain_tot)+s(rain_tot_lag1)+s(temp_tot_lag1)+s(humidity_tot_lag1)+s(wind_tot_lag2)+month+payweek +s(admissions_lag1)+s(admissions_lag3)+s(admissions_lag7)+s(admissions_lag9) +I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)), family="quasipoisson") summary(m_adm6.1) #Weekday exp(m_adm6.1$coefficients["loads3_last2"]) exp(m_adm6.1$coefficients["loads3_last2"] - qnorm(0.975)*summary(m_adm6.1)[[22]]["loads3_last2",2]) exp(m_adm6.1$coefficients["loads3_last2"] + qnorm(0.975)*summary(m_adm6.1)[[22]]["loads3_last2",2]) #Weekend sigma <- sqrt(sqrt(m_adm6.1$Vp[2,2])^2 + sqrt(m_adm6.1$Vp[20,20])^2 + 2*(m_adm6.1$Vp[2,20])^2) # see Heumann, Schomaker, Shalabh "Introduction so Statisics", Springer 2016, p.279 for formula exp(m_adm6.1$coefficients["loads3_last2"]+m_adm6.1$coefficients["loads3_last2:weekendTRUE"]) exp(m_adm6.1$coefficients["loads3_last2"]+m_adm6.1$coefficients["loads3_last2:weekendTRUE"]-qnorm(0.975)*sigma) exp(m_adm6.1$coefficients["loads3_last2"]+m_adm6.1$coefficients["loads3_last2:weekendTRUE"]+qnorm(0.975)*sigma) # c) estimates for different lags m_adm8 <- gam(admissions ~ loads3 +s(sunshine_ap)+s(wind_tot)+s(rain_tot)+s(rain_tot_lag1)+s(temp_tot_lag1)+s(humidity_tot_lag1)+s(wind_tot_lag2)+month+payweek +s(admissions_lag1)+s(admissions_lag3)+s(admissions_lag7)+s(admissions_lag9) +I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)), family="quasipoisson") m_adm9 <- gam(admissions ~ loads3_lag1 +s(sunshine_ap)+s(wind_tot)+s(rain_tot)+s(rain_tot_lag1)+s(temp_tot_lag1)+s(humidity_tot_lag1)+s(wind_tot_lag2)+month+payweek +s(admissions_lag1)+s(admissions_lag3)+s(admissions_lag7)+s(admissions_lag9) +I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)), family="quasipoisson") m_adm10 <- gam(admissions ~ loads3_lag2 +s(sunshine_ap)+s(wind_tot)+s(rain_tot)+s(rain_tot_lag1)+s(temp_tot_lag1)+s(humidity_tot_lag1)+s(wind_tot_lag2)+month+payweek +s(admissions_lag1)+s(admissions_lag3)+s(admissions_lag7)+s(admissions_lag9) +I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)), family="quasipoisson") m_adm11 <- gam(admissions ~ loads3_lag3 +s(sunshine_ap)+s(wind_tot)+s(rain_tot)+s(rain_tot_lag1)+s(temp_tot_lag1)+s(humidity_tot_lag1)+s(wind_tot_lag2)+month+payweek +s(admissions_lag1)+s(admissions_lag3)+s(admissions_lag7)+s(admissions_lag9) +I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)), family="quasipoisson") m_adm12 <- gam(admissions ~ loads3_lag4 +s(sunshine_ap)+s(wind_tot)+s(rain_tot)+s(rain_tot_lag1)+s(temp_tot_lag1)+s(humidity_tot_lag1)+s(wind_tot_lag2)+month+payweek +s(admissions_lag1)+s(admissions_lag3)+s(admissions_lag7)+s(admissions_lag9) +I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)), family="quasipoisson") summary(m_adm8) summary(m_adm9) summary(m_adm10) summary(m_adm11) summary(m_adm12) # Odds Ratio and cofidence interval exp(m_adm8$coefficients["loads3"]) exp(m_adm8$coefficients["loads3"] - qnorm(0.975)*summary(m_adm8)[[22]]["loads3",2]) exp(m_adm8$coefficients["loads3"] + qnorm(0.975)*summary(m_adm8)[[22]]["loads3",2]) exp(m_adm9$coefficients["loads3_lag1"]) exp(m_adm9$coefficients["loads3_lag1"] - qnorm(0.975)*summary(m_adm9)[[22]]["loads3_lag1",2]) exp(m_adm9$coefficients["loads3_lag1"] + qnorm(0.975)*summary(m_adm9)[[22]]["loads3_lag1",2]) exp(m_adm10$coefficients["loads3_lag2"]) exp(m_adm10$coefficients["loads3_lag2"] - qnorm(0.975)*summary(m_adm10)[[22]]["loads3_lag2",2]) exp(m_adm10$coefficients["loads3_lag2"] + qnorm(0.975)*summary(m_adm10)[[22]]["loads3_lag2",2]) exp(m_adm11$coefficients["loads3_lag3"]) exp(m_adm11$coefficients["loads3_lag3"] - qnorm(0.975)*summary(m_adm11)[[22]]["loads3_lag3",2]) exp(m_adm11$coefficients["loads3_lag3"] + qnorm(0.975)*summary(m_adm11)[[22]]["loads3_lag3",2]) exp(m_adm12$coefficients["loads3_lag4"]) exp(m_adm12$coefficients["loads3_lag4"] - qnorm(0.975)*summary(m_adm12)[[22]]["loads3_lag4",2]) exp(m_adm12$coefficients["loads3_lag4"] + qnorm(0.975)*summary(m_adm12)[[22]]["loads3_lag4",2]) # Models for subwards and ICD10 codes m_adm13 <- gam(formula(paste("SpecgroupMedicine ~", covariates)), family="quasipoisson") exp(m_adm13$coefficients["loads3_last2"]) exp(m_adm13$coefficients["loads3_last2"] - qnorm(0.975)*summary(m_adm13)[[22]]["loads3_last2",2]) exp(m_adm13$coefficients["loads3_last2"] + qnorm(0.975)*summary(m_adm13)[[22]]["loads3_last2",2]) m_adm14 <- gam(formula(paste("SpecgroupSurgery ~", covariates)), family="quasipoisson") exp(m_adm14$coefficients["loads3_last2"]) exp(m_adm14$coefficients["loads3_last2"] - qnorm(0.975)*summary(m_adm14)[[22]]["loads3_last2",2]) exp(m_adm14$coefficients["loads3_last2"] + qnorm(0.975)*summary(m_adm14)[[22]]["loads3_last2",2]) m_adm15 <- gam(formula(paste("ICD_Infections ~", covariates)), family="quasipoisson") exp(m_adm15$coefficients["loads3_last2"]) exp(m_adm15$coefficients["loads3_last2"] - qnorm(0.975)*summary(m_adm15)[[22]]["loads3_last2",2]) exp(m_adm15$coefficients["loads3_last2"] + qnorm(0.975)*summary(m_adm15)[[22]]["loads3_last2",2]) m_adm16 <- gam(formula(paste("ICD_Respiration ~", covariates)), family="quasipoisson") exp(m_adm16$coefficients["loads3_last2"]) exp(m_adm16$coefficients["loads3_last2"] - qnorm(0.975)*summary(m_adm16)[[22]]["loads3_last2",2]) exp(m_adm16$coefficients["loads3_last2"] + qnorm(0.975)*summary(m_adm16)[[22]]["loads3_last2",2]) m_adm17 <- gam(formula(paste("ICD_InjuryPoisoning ~", covariates)), family="quasipoisson") exp(m_adm17$coefficients["loads3_last2"]) exp(m_adm17$coefficients["loads3_last2"] - qnorm(0.975)*summary(m_adm17)[[22]]["loads3_last2",2]) exp(m_adm17$coefficients["loads3_last2"] + qnorm(0.975)*summary(m_adm17)[[22]]["loads3_last2",2]) m_adm18 <- gam(formula(paste("ICD_EyeEar ~", covariates)), family="quasipoisson") exp(m_adm18$coefficients["loads3_last2"]) exp(m_adm18$coefficients["loads3_last2"] - qnorm(0.975)*summary(m_adm18)[[22]]["loads3_last2",2]) exp(m_adm18$coefficients["loads3_last2"] + qnorm(0.975)*summary(m_adm18)[[22]]["loads3_last2",2]) m_adm19 <- gam(formula(paste("ICD_Skin ~", covariates)), family="quasipoisson") exp(m_adm19$coefficients["loads3_last2"]) exp(m_adm19$coefficients["loads3_last2"] - qnorm(0.975)*summary(m_adm19)[[22]]["loads3_last2",2]) exp(m_adm19$coefficients["loads3_last2"] + qnorm(0.975)*summary(m_adm19)[[22]]["loads3_last2",2]) m_adm20 <- gam(formula(paste("ICD_Digestion ~", covariates)), family="quasipoisson") exp(m_adm20$coefficients["loads3_last2"]) exp(m_adm20$coefficients["loads3_last2"] - qnorm(0.975)*summary(m_adm20)[[22]]["loads3_last2",2]) exp(m_adm20$coefficients["loads3_last2"] + qnorm(0.975)*summary(m_adm20)[[22]]["loads3_last2",2]) adm$Other <- ICD_Endocrine + ICD_BehavNeurodev + ICD_Neurology + ICD_Circulation + ICD_Musculoscel + ICD_Urogenital + ICD_Perinatal + ICD_Congenital + ICD_NoClass + ICD_BehavNeurodev.1 + ICD_NeoOnco + ICD_External m_adm21 <- gam(formula(paste("adm$Other ~", covariates)), family="quasipoisson") exp(m_adm21$coefficients["loads3_last2"]) exp(m_adm21$coefficients["loads3_last2"] - qnorm(0.975)*summary(m_adm21)[[22]]["loads3_last2",2]) exp(m_adm21$coefficients["loads3_last2"] + qnorm(0.975)*summary(m_adm21)[[22]]["loads3_last2",2]) ########################### # iv) Alternative models # ########################### # negative binomial m_adm39 <- glm.nb(admissions ~ loads3_last2 +sunshine_ap+wind_tot+rain_tot+rain_tot_lag1+temp_tot_lag1+humidity_tot_lag1+wind_tot_lag2+month+payweek +admissions_lag1+admissions_lag3+admissions_lag7+admissions_lag9 +I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7))) exp(m_adm39$coefficients["loads3_last2"]) exp(m_adm39$coefficients["loads3_last2"] - qnorm(0.975)*summary(m_adm39)[[12]]["loads3_last2",2]) exp(m_adm39$coefficients["loads3_last2"] + qnorm(0.975)*summary(m_adm39)[[12]]["loads3_last2",2]) # INGARCH adm$z1 <- I(sin(pi*2*1*adm$t/7)) adm$z2 <- I(cos(pi*2*1*adm$t/7)) adm$z3 <- I(sin(pi*2*2*adm$t/7)) adm$z4 <- I(cos(pi*2*2*adm$t/7)) adm <- cbind(adm,model.matrix(~-1+adm$month)[,-1]) colnames(adm)[dim(adm)[2]:(dim(adm)[2]-10)]<- c("Aug","Dec","Feb","Jan","Jul","Jun","Mar","May","Nov","Oct","Sep") m_t1 <- tsglm(adm$admissions, model=list(past_obs=c(1,3,7,9), past_mean=c(1:3,7,28)), xreg=adm[,c("loads3_last2","Aug","Dec","Feb","Jan","Jul","Jun","Mar","May","Nov","Oct","Sep","sunshine_ap","wind_tot","rain_tot","rain_tot_lag1","temp_tot_lag1","humidity_tot_lag1","wind_tot_lag2","payweek","z1","z2","z3","z4","admissions_lag1","admissions_lag3","admissions_lag7","admissions_lag9")], link="log", distr="nbinom") m_t2 <- tsglm(adm$admissions, model=list(past_obs=c(1,3,7,9), past_mean=c(1:3)), xreg=adm[,c("loads3_last2","Aug","Dec","Feb","Jan","Jul","Jun","Mar","May","Nov","Oct","Sep","sunshine_ap","wind_tot","rain_tot","rain_tot_lag1","temp_tot_lag1","humidity_tot_lag1","wind_tot_lag2","payweek","z1","z2","z3","z4","admissions_lag1","admissions_lag3","admissions_lag7","admissions_lag9")], link="log", distr="nbinom") m_t3 <- tsglm(adm$admissions, model=list(past_obs=c(1,3,7,9), past_mean=c(1:3,7)), xreg=adm[,c("loads3_last2","Aug","Dec","Feb","Jan","Jul","Jun","Mar","May","Nov","Oct","Sep","sunshine_ap","wind_tot","rain_tot","rain_tot_lag1","temp_tot_lag1","humidity_tot_lag1","wind_tot_lag2","payweek","z1","z2","z3","z4","admissions_lag1","admissions_lag3","admissions_lag7","admissions_lag9")], link="log", distr="nbinom") m_t4 <- tsglm(adm$admissions, model=list(past_obs=c(1,3,7,9), past_mean=c(1,7,14,21,28)), xreg=adm[,c("loads3_last2","Aug","Dec","Feb","Jan","Jul","Jun","Mar","May","Nov","Oct","Sep","sunshine_ap","wind_tot","rain_tot","rain_tot_lag1","temp_tot_lag1","humidity_tot_lag1","wind_tot_lag2","payweek","z1","z2","z3","z4","admissions_lag1","admissions_lag3","admissions_lag7","admissions_lag9")], link="log", distr="nbinom") m_t5 <- tsglm(adm$admissions, model=list(past_obs=c(1,3,7,9), past_mean=c(1,3,7,10,14)), xreg=adm[,c("loads3_last2","Aug","Dec","Feb","Jan","Jul","Jun","Mar","May","Nov","Oct","Sep","sunshine_ap","wind_tot","rain_tot","rain_tot_lag1","temp_tot_lag1","humidity_tot_lag1","wind_tot_lag2","payweek","z1","z2","z3","z4","admissions_lag1","admissions_lag3","admissions_lag7","admissions_lag9")], link="log", distr="nbinom") m_t6 <- tsglm(adm$admissions, model=list(past_obs=c(1,3,7,9), past_mean=c(1:5)), xreg=adm[,c("loads3_last2","Aug","Dec","Feb","Jan","Jul","Jun","Mar","May","Nov","Oct","Sep","sunshine_ap","wind_tot","rain_tot","rain_tot_lag1","temp_tot_lag1","humidity_tot_lag1","wind_tot_lag2","payweek","z1","z2","z3","z4","admissions_lag1","admissions_lag3","admissions_lag7","admissions_lag9")], link="log", distr="nbinom") m_t7 <- tsglm(adm$admissions, model=list(past_obs=c(1,3,7,9), past_mean=c(1)), xreg=adm[,c("loads3_last2","Aug","Dec","Feb","Jan","Jul","Jun","Mar","May","Nov","Oct","Sep","sunshine_ap","wind_tot","rain_tot","rain_tot_lag1","temp_tot_lag1","humidity_tot_lag1","wind_tot_lag2","payweek","z1","z2","z3","z4","admissions_lag1","admissions_lag3","admissions_lag7","admissions_lag9")], link="log", distr="nbinom") m_t8 <- tsglm(adm$admissions, model=list(past_obs=c(1,3,7,9), past_mean=NULL), xreg=adm[,c("loads3_last2","Aug","Dec","Feb","Jan","Jul","Jun","Mar","May","Nov","Oct","Sep","sunshine_ap","wind_tot","rain_tot","rain_tot_lag1","temp_tot_lag1","humidity_tot_lag1","wind_tot_lag2","payweek","z1","z2","z3","z4","admissions_lag1","admissions_lag3","admissions_lag7","admissions_lag9")], link="log", distr="nbinom") m_t9 <- tsglm(adm$admissions, model=list(past_obs=c(1,3,7,9), past_mean=c(1:3,7,14,28)), xreg=adm[,c("loads3_last2","Aug","Dec","Feb","Jan","Jul","Jun","Mar","May","Nov","Oct","Sep","sunshine_ap","wind_tot","rain_tot","rain_tot_lag1","temp_tot_lag1","humidity_tot_lag1","wind_tot_lag2","payweek","z1","z2","z3","z4","admissions_lag1","admissions_lag3","admissions_lag7","admissions_lag9")], link="log", distr="nbinom") m_t10 <- tsglm(adm$admissions, model=list(past_obs=c(1,3,7,9), past_mean=c(1:3,7,14,21,28)), xreg=adm[,c("loads3_last2","Aug","Dec","Feb","Jan","Jul","Jun","Mar","May","Nov","Oct","Sep","sunshine_ap","wind_tot","rain_tot","rain_tot_lag1","temp_tot_lag1","humidity_tot_lag1","wind_tot_lag2","payweek","z1","z2","z3","z4","admissions_lag1","admissions_lag3","admissions_lag7","admissions_lag9")], link="log", distr="nbinom") AIC(m_t1) AIC(m_t2) # lowest AIC AIC(m_t3) AIC(m_t4) AIC(m_t5) AIC(m_t6) AIC(m_t7) AIC(m_t8) AIC(m_t9) AIC(m_t10) summary(m_t2) exp(m_t2$coefficients["loads3_last2"]) exp(m_t2$coefficients["loads3_last2"] - qnorm(0.975)*summary(m_t2)[[5]]["loads3_last2",2]) exp(m_t2$coefficients["loads3_last2"] + qnorm(0.975)*summary(m_t2)[[5]]["loads3_last2",2]) pdf(file=paste(getwd(),"/eFigure6.pdf", sep=""), width=10) par(mfrow=c(2,3)) plot(m_t4, ask=F) dev.off() detach(adm) detach("package:mgcv", unload=TRUE) ptm4 <- proc.time() ########################### # PART 4: Causal Analysis # ########################### adm <- read.csv('final_data_all.csv') adm$date <- as.Date(adm$date, "%Y-%m-%d") attach(adm) adm$z1 <- I(sin(pi*2*1*adm$t/7)) adm$z2 <- I(cos(pi*2*1*adm$t/7)) adm$z3 <- I(sin(pi*2*2*adm$t/7)) adm$z4 <- I(cos(pi*2*2*adm$t/7)) W <- as.matrix(cbind(adm$month, adm$sunshine_ap, adm$wind_tot, adm$rain_tot, adm$rain_tot_lag1, adm$temp_tot_lag1, adm$humidity_tot_lag1, adm$wind_tot_lag2, adm$payweek, adm$z1,adm$z2,adm$z3,adm$z4, adm$admissions_lag1, adm$admissions_lag3, adm$admissions_lag7, adm$admissions_lag9)) # TMLE (with super learning) mylibrary <- list("SL.glm", "SL.stepAIC","SL.step.forward", "SL.step", "SL.bayesglm","SL.mean", "SL.step.interaction", "SL.glmnet","SL.lae","SL.lae2","SL.lae.int","SL.gbm","SL.earth","SL.mma","SL.jma", c("SL.glm","screen.glmnet3"),c("SL.earth2","screen.glmnet3"),c("SL.gbm2","screen.cramersv")) mylibrary2 <- list("SL.glm", "SL.stepAIC" , "SL.step.forward", "SL.step", "SL.bayesglm","SL.mean","SL.lae","SL.lae2","SL.lae.int", "SL.step.interaction","SL.gbm","SL.earth","SL.knn", c("SL.glm","screen.glmnet3"),c("SL.earth2","screen.glmnet3"),c("SL.gbm2","screen.cramersv")) M_tmle_final <- tmle(admissions,loads3_last2, W, family="gaussian", gbound = c(0.01, 1), V=10, verbose=T, Q.SL.library = mylibrary, g.SL.library = mylibrary2) M_tmle_final M_tmle_final$Qinit$coef M_tmle_final$g$coef # Naive estimate (linear model) naive <- glm(admissions ~ loads3_last2 + sunshine_ap+wind_tot+rain_tot+rain_tot_lag1+temp_tot_lag1+humidity_tot_lag1+wind_tot_lag2+month+payweek +admissions_lag1+admissions_lag3+admissions_lag7+admissions_lag9 +I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)), data=adm, family="gaussian") summary(naive) confint(naive) # Table 1 (summary of Part 3 and Part 4) mr1 <- c(exp(m_adm6$coefficients["loads3_last2"]), exp(m_adm6$coefficients["loads3_last2"] - qnorm(0.975)*summary(m_adm6)[[22]]["loads3_last2",2]), exp(m_adm6$coefficients["loads3_last2"] + qnorm(0.975)*summary(m_adm6)[[22]]["loads3_last2",2])) mr2 <- c(exp(m_adm6.1$coefficients["loads3_last2"]), exp(m_adm6.1$coefficients["loads3_last2"] - qnorm(0.975)*summary(m_adm6.1)[[22]]["loads3_last2",2]), exp(m_adm6.1$coefficients["loads3_last2"] + qnorm(0.975)*summary(m_adm6.1)[[22]]["loads3_last2",2])) mr3 <- c(exp(m_adm6.1$coefficients["loads3_last2"]+m_adm6.1$coefficients["loads3_last2:weekendTRUE"]), exp(m_adm6.1$coefficients["loads3_last2"]+m_adm6.1$coefficients["loads3_last2:weekendTRUE"]-qnorm(0.975)*sigma), exp(m_adm6.1$coefficients["loads3_last2"]+m_adm6.1$coefficients["loads3_last2:weekendTRUE"]+qnorm(0.975)*sigma)) mr4 <- c(exp(m_adm8$coefficients["loads3"]), exp(m_adm8$coefficients["loads3"] - qnorm(0.975)*summary(m_adm8)[[22]]["loads3",2]), exp(m_adm8$coefficients["loads3"] + qnorm(0.975)*summary(m_adm8)[[22]]["loads3",2])) mr5 <- c(exp(m_adm9$coefficients["loads3_lag1"]), exp(m_adm9$coefficients["loads3_lag1"] - qnorm(0.975)*summary(m_adm9)[[22]]["loads3_lag1",2]), exp(m_adm9$coefficients["loads3_lag1"] + qnorm(0.975)*summary(m_adm9)[[22]]["loads3_lag1",2])) mr6 <- c(exp(m_adm10$coefficients["loads3_lag2"]), exp(m_adm10$coefficients["loads3_lag2"] - qnorm(0.975)*summary(m_adm10)[[22]]["loads3_lag2",2]), exp(m_adm10$coefficients["loads3_lag2"] + qnorm(0.975)*summary(m_adm10)[[22]]["loads3_lag2",2])) mr7 <- c(exp(m_adm11$coefficients["loads3_lag3"]), exp(m_adm11$coefficients["loads3_lag3"] - qnorm(0.975)*summary(m_adm11)[[22]]["loads3_lag3",2]), exp(m_adm11$coefficients["loads3_lag3"] + qnorm(0.975)*summary(m_adm11)[[22]]["loads3_lag3",2])) mr8 <- c(exp(m_adm12$coefficients["loads3_lag4"]), exp(m_adm12$coefficients["loads3_lag4"] - qnorm(0.975)*summary(m_adm12)[[22]]["loads3_lag4",2]), exp(m_adm12$coefficients["loads3_lag4"] + qnorm(0.975)*summary(m_adm12)[[22]]["loads3_lag4",2])) mr9 <- c(exp(m_adm14$coefficients["loads3_last2"]), exp(m_adm14$coefficients["loads3_last2"] - qnorm(0.975)*summary(m_adm14)[[22]]["loads3_last2",2]), exp(m_adm14$coefficients["loads3_last2"] + qnorm(0.975)*summary(m_adm14)[[22]]["loads3_last2",2])) mr10 <- c(exp(m_adm13$coefficients["loads3_last2"]) , exp(m_adm13$coefficients["loads3_last2"] - qnorm(0.975)*summary(m_adm13)[[22]]["loads3_last2",2]), exp(m_adm13$coefficients["loads3_last2"] + qnorm(0.975)*summary(m_adm13)[[22]]["loads3_last2",2])) mr11 <- c(exp(m_adm15$coefficients["loads3_last2"]), exp(m_adm15$coefficients["loads3_last2"] - qnorm(0.975)*summary(m_adm15)[[22]]["loads3_last2",2]), exp(m_adm15$coefficients["loads3_last2"] + qnorm(0.975)*summary(m_adm15)[[22]]["loads3_last2",2])) mr12 <- c(exp(m_adm18$coefficients["loads3_last2"]), exp(m_adm18$coefficients["loads3_last2"] - qnorm(0.975)*summary(m_adm18)[[22]]["loads3_last2",2]), exp(m_adm18$coefficients["loads3_last2"] + qnorm(0.975)*summary(m_adm18)[[22]]["loads3_last2",2])) mr13 <- c(exp(m_adm16$coefficients["loads3_last2"]), exp(m_adm16$coefficients["loads3_last2"] - qnorm(0.975)*summary(m_adm16)[[22]]["loads3_last2",2]), exp(m_adm16$coefficients["loads3_last2"] + qnorm(0.975)*summary(m_adm16)[[22]]["loads3_last2",2])) mr14 <- c(exp(m_adm20$coefficients["loads3_last2"]), exp(m_adm20$coefficients["loads3_last2"] - qnorm(0.975)*summary(m_adm20)[[22]]["loads3_last2",2]), exp(m_adm20$coefficients["loads3_last2"] + qnorm(0.975)*summary(m_adm20)[[22]]["loads3_last2",2])) mr15 <- c(exp(m_adm19$coefficients["loads3_last2"]), exp(m_adm19$coefficients["loads3_last2"] - qnorm(0.975)*summary(m_adm19)[[22]]["loads3_last2",2]), exp(m_adm19$coefficients["loads3_last2"] + qnorm(0.975)*summary(m_adm19)[[22]]["loads3_last2",2])) mr16 <- c(exp(m_adm17$coefficients["loads3_last2"]), exp(m_adm17$coefficients["loads3_last2"] - qnorm(0.975)*summary(m_adm17)[[22]]["loads3_last2",2]), exp(m_adm17$coefficients["loads3_last2"] + qnorm(0.975)*summary(m_adm17)[[22]]["loads3_last2",2])) mr17 <- c(exp(m_adm21$coefficients["loads3_last2"]), exp(m_adm21$coefficients["loads3_last2"] - qnorm(0.975)*summary(m_adm21)[[22]]["loads3_last2",2]) , exp(m_adm21$coefficients["loads3_last2"] + qnorm(0.975)*summary(m_adm21)[[22]]["loads3_last2",2])) mr18 <- c(M_tmle_final$estimates$ATE$psi,M_tmle_final$estimates$ATE$CI) mr19 <- c(coef(naive)[2],confint(naive)[2,]) mi <- i3 mr20 <- c(mi[names(mi)=="loads3"],mi[names(mi)=="loads3_lag1"],mi[names(mi)=="loads3_lag2"],mi[names(mi)=="loads3_lag3"],mi[names(mi)=="loads3_lag4"]) table1 <- matrix(NA,nrow=30,ncol=5) table1[1,2:4] <- mr1 table1[3,2:4] <- mr2 table1[4,2:4] <- mr3 table1[8,2:4] <- mr4 table1[9,2:4] <- mr5 table1[10,2:4] <- mr6 table1[11,2:4] <- mr7 table1[12,2:4] <- mr8 table1[15,2:4] <- mr9 table1[16,2:4] <- mr10 table1[19,2:4] <- mr11 table1[20,2:4] <- mr12 table1[21,2:4] <- mr13 table1[22,2:4] <- mr14 table1[23,2:4] <- mr15 table1[24,2:4] <- mr16 table1[25,2:4] <- mr17 table1[28,2:4] <- mr18 table1[29,2:4] <- mr19 table1[8:12,5] <- mr20 table1 <- round(table1, digits=2) t1names <- c("LS: same day or up to 2 days prior","Interaction:","Weekday","Weekend","Interaction: length of load shedding", "see eFigure3","Other models:","LS: same day","LS: 1 day prior","LS: 2 days prior","LS: 3 days prior","LS: 4 days prior","","By specialty:", "Surgical cases","Medical cases","","By ICD-10 codes (code range):","Certain Infections (A00-B99)","Eye & Ear (H00-H95)","Respiratory system (J00-J99)", "Digestive system (K00-K93)","Skin (L00-L99)","Injuries, Poisoning (S00-Y98)","Other","","","ATE","naive regression","") tab1_save <- data.frame(t1names,table1[,2],paste(table1[,3],table1[,4], sep=";"),table1[,5]) colnames(tab1_save) <- c("","est","95%CI","VI") tab1_save$est <- as.character(tab1_save$est) tab1_save$VI <- as.character(tab1_save$VI) write.xlsx(tab1_save,file=paste(getwd(),"/table1.xlsx",sep="")) detach(adm) detach("package:tmle", unload=TRUE) ptm5 <- proc.time() ############################################ # PART 5: results by area and length of LS # ############################################ adm <- read.csv('final_data_all.csv') adm$date <- as.Date(adm$date, "%Y-%m-%d") # association by length of loadshedding, i.e. interaction hypothesis (non-linear and linear) library(mgcv) attach(adm) # non-linear m_adm_a1 <- gam(admissions ~ s(area_all_time, by=loads3_last2) +s(sunshine_ap)+s(wind_tot)+s(rain_tot)+s(rain_tot_lag1)+s(temp_tot_lag1)+s(humidity_tot_lag1)+s(wind_tot_lag2) + month + payweek + s(admissions_lag1) + s(admissions_lag3) + s(admissions_lag7) + s(admissions_lag9) + I(sin(pi * 2 * 1 * t/7)) + I(cos(pi * 2 * 1 * t/7)) + I(sin(pi * 2 * 2 * t/7)) + I(cos(pi * 2 * 2 * t/7)), family="quasipoisson") summary(m_adm_a1) #eFigure3 pdf(file=paste(getwd(),"/eFigure3.pdf",sep=""), width=9) plot(m_adm_a1, shade=T,select=1,lwd=2, scale=0, xlab="Length of loadshedding (average over all areas)", ylab="IRR", trans=exp, ylim=c(0.8,1.5)) dev.off() # linear m_adm_a2 <- gam(admissions ~ area_all_time*loads3_last2 +s(sunshine_ap)+s(wind_tot)+s(rain_tot)+s(rain_tot_lag1)+s(temp_tot_lag1)+s(humidity_tot_lag1)+s(wind_tot_lag2) + month + payweek + s(admissions_lag1) + s(admissions_lag3) + s(admissions_lag7) + s(admissions_lag9) + I(sin(pi * 2 * 1 * t/7)) + I(cos(pi * 2 * 1 * t/7)) + I(sin(pi * 2 * 2 * t/7)) + I(cos(pi * 2 * 2 * t/7)), family="quasipoisson") summary(m_adm_a2) detach(adm) # by area m_adm_area1 <- gam(admissions ~ area1_last2 +s(sunshine_ap)+s(wind_tot)+s(rain_tot)+s(rain_tot_lag1)+s(temp_tot_lag1)+s(humidity_tot_lag1)+s(wind_tot_lag2)+month+payweek +s(admissions_lag1)+s(admissions_lag3)+s(admissions_lag7)+s(admissions_lag9)+I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)), family="quasipoisson", data=adm) m_adm_area2 <- gam(admissions ~ area2_last2 +s(sunshine_ap)+s(wind_tot)+s(rain_tot)+s(rain_tot_lag1)+s(temp_tot_lag1)+s(humidity_tot_lag1)+s(wind_tot_lag2)+month+payweek +s(admissions_lag1)+s(admissions_lag3)+s(admissions_lag7)+s(admissions_lag9)+I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)), family="quasipoisson", data=adm) m_adm_area3 <- gam(admissions ~ area3_last2 +s(sunshine_ap)+s(wind_tot)+s(rain_tot)+s(rain_tot_lag1)+s(temp_tot_lag1)+s(humidity_tot_lag1)+s(wind_tot_lag2)+month+payweek +s(admissions_lag1)+s(admissions_lag3)+s(admissions_lag7)+s(admissions_lag9)+I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)), family="quasipoisson", data=adm) m_adm_area4 <- gam(admissions ~ area4_last2 +s(sunshine_ap)+s(wind_tot)+s(rain_tot)+s(rain_tot_lag1)+s(temp_tot_lag1)+s(humidity_tot_lag1)+s(wind_tot_lag2)+month+payweek +s(admissions_lag1)+s(admissions_lag3)+s(admissions_lag7)+s(admissions_lag9)+I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)), family="quasipoisson", data=adm) m_adm_area5 <- gam(admissions ~ area5_last2 +s(sunshine_ap)+s(wind_tot)+s(rain_tot)+s(rain_tot_lag1)+s(temp_tot_lag1)+s(humidity_tot_lag1)+s(wind_tot_lag2)+month+payweek +s(admissions_lag1)+s(admissions_lag3)+s(admissions_lag7)+s(admissions_lag9)+I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)), family="quasipoisson", data=adm) m_adm_area6 <- gam(admissions ~ area6_last2 +s(sunshine_ap)+s(wind_tot)+s(rain_tot)+s(rain_tot_lag1)+s(temp_tot_lag1)+s(humidity_tot_lag1)+s(wind_tot_lag2)+month+payweek +s(admissions_lag1)+s(admissions_lag3)+s(admissions_lag7)+s(admissions_lag9)+I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)), family="quasipoisson", data=adm) m_adm_area7 <- gam(admissions ~ area7_last2 +s(sunshine_ap)+s(wind_tot)+s(rain_tot)+s(rain_tot_lag1)+s(temp_tot_lag1)+s(humidity_tot_lag1)+s(wind_tot_lag2)+month+payweek +s(admissions_lag1)+s(admissions_lag3)+s(admissions_lag7)+s(admissions_lag9)+I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)), family="quasipoisson", data=adm) m_adm_area8 <- gam(admissions ~ area8_last2 +s(sunshine_ap)+s(wind_tot)+s(rain_tot)+s(rain_tot_lag1)+s(temp_tot_lag1)+s(humidity_tot_lag1)+s(wind_tot_lag2)+month+payweek +s(admissions_lag1)+s(admissions_lag3)+s(admissions_lag7)+s(admissions_lag9)+I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)), family="quasipoisson", data=adm) m_adm_area9 <- gam(admissions ~ area9_last2 +s(sunshine_ap)+s(wind_tot)+s(rain_tot)+s(rain_tot_lag1)+s(temp_tot_lag1)+s(humidity_tot_lag1)+s(wind_tot_lag2)+month+payweek +s(admissions_lag1)+s(admissions_lag3)+s(admissions_lag7)+s(admissions_lag9)+I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)), family="quasipoisson", data=adm) m_adm_area10 <- gam(admissions ~ area10_last2 +s(sunshine_ap)+s(wind_tot)+s(rain_tot)+s(rain_tot_lag1)+s(temp_tot_lag1)+s(humidity_tot_lag1)+s(wind_tot_lag2)+month+payweek +s(admissions_lag1)+s(admissions_lag3)+s(admissions_lag7)+s(admissions_lag9)+I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)), family="quasipoisson", data=adm) m_adm_area11 <- gam(admissions ~ area11_last2 +s(sunshine_ap)+s(wind_tot)+s(rain_tot)+s(rain_tot_lag1)+s(temp_tot_lag1)+s(humidity_tot_lag1)+s(wind_tot_lag2)+month+payweek +s(admissions_lag1)+s(admissions_lag3)+s(admissions_lag7)+s(admissions_lag9)+I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)), family="quasipoisson", data=adm) m_adm_area12 <- gam(admissions ~ area12_last2 +s(sunshine_ap)+s(wind_tot)+s(rain_tot)+s(rain_tot_lag1)+s(temp_tot_lag1)+s(humidity_tot_lag1)+s(wind_tot_lag2)+month+payweek +s(admissions_lag1)+s(admissions_lag3)+s(admissions_lag7)+s(admissions_lag9)+I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)), family="quasipoisson", data=adm) m_adm_area13 <- gam(admissions ~ area13_last2 +s(sunshine_ap)+s(wind_tot)+s(rain_tot)+s(rain_tot_lag1)+s(temp_tot_lag1)+s(humidity_tot_lag1)+s(wind_tot_lag2)+month+payweek +s(admissions_lag1)+s(admissions_lag3)+s(admissions_lag7)+s(admissions_lag9)+I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)), family="quasipoisson", data=adm) m_adm_area14 <- gam(admissions ~ area14_last2 +s(sunshine_ap)+s(wind_tot)+s(rain_tot)+s(rain_tot_lag1)+s(temp_tot_lag1)+s(humidity_tot_lag1)+s(wind_tot_lag2)+month+payweek +s(admissions_lag1)+s(admissions_lag3)+s(admissions_lag7)+s(admissions_lag9)+I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)), family="quasipoisson", data=adm) m_adm_area15 <- gam(admissions ~ area15_last2 +s(sunshine_ap)+s(wind_tot)+s(rain_tot)+s(rain_tot_lag1)+s(temp_tot_lag1)+s(humidity_tot_lag1)+s(wind_tot_lag2)+month+payweek +s(admissions_lag1)+s(admissions_lag3)+s(admissions_lag7)+s(admissions_lag9)+I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)), family="quasipoisson", data=adm) m_adm_area16 <- gam(admissions ~ area16_last2 +s(sunshine_ap)+s(wind_tot)+s(rain_tot)+s(rain_tot_lag1)+s(temp_tot_lag1)+s(humidity_tot_lag1)+s(wind_tot_lag2)+month+payweek +s(admissions_lag1)+s(admissions_lag3)+s(admissions_lag7)+s(admissions_lag9)+I(sin(pi*2*1*t/7)) +I(cos(pi*2*1*t/7)) +I(sin(pi*2*2*t/7)) +I(cos(pi*2*2*t/7)), family="quasipoisson", data=adm) # table of results used for Figure 3 fig3 <- data.frame(matrix(NA,ncol=2,nrow=16)) colnames(fig3) <- c("area","IRR") IRRS <- c(exp(m_adm_area1$coefficients["area1_last2"]), exp(m_adm_area2$coefficients["area2_last2"]), exp(m_adm_area3$coefficients["area3_last2"]), exp(m_adm_area4$coefficients["area4_last2"]), exp(m_adm_area5$coefficients["area5_last2"]), exp(m_adm_area6$coefficients["area6_last2"]), exp(m_adm_area7$coefficients["area7_last2"]), exp(m_adm_area8$coefficients["area8_last2"]), exp(m_adm_area9$coefficients["area9_last2"]), exp(m_adm_area10$coefficients["area10_last2"]), exp(m_adm_area11$coefficients["area11_last2"]), exp(m_adm_area12$coefficients["area12_last2"]), exp(m_adm_area13$coefficients["area13_last2"]), exp(m_adm_area14$coefficients["area14_last2"]), exp(m_adm_area15$coefficients["area15_last2"]), exp(m_adm_area16$coefficients["area16_last2"])) fig3$IRR <- IRRS areanames <- c("Belville","Maitland/Milnerton","Somerset West/Strand","Mitchell's Plain/Philippi Farmlands","Hanover Park/Lansdown/Newlands","Durbanville/Tyger Valley", "City Bowl/Camps Bay","Southern Peninsula","Pinelands/Langa/Industrial","Brackenfell/Kuils River", "Plumstead/Ottery/Hout Bay/Constantia","Athlone/Manenberg/Gugulethu","Goodwood/Parow/Plattekloof","Atlantis/Melkbosstrand","Obz/Rondebosch/Brooklyn","Philippi/Retreat") fig3$area <- areanames fig3[order(fig3$IRR),] write.xlsx(fig3,file=paste(getwd(),"/figure3_table_of_results.xlsx",sep="")) ptm6 <- proc.time() ################################### # PART 6: results by diagnosis # ################################### adm <- read.csv('final_data_all.csv') attach(adm) # all diagnoses with number cases >25 index3 <- grep("A05",colnames(adm)) ind25 <- rep(NA,dim(adm)[2]) for(i in index3:dim(adm)[2])try({ if(sum(adm[,i])>25){ind25[i]<-TRUE}else{ind25[i]<-FALSE} }) ind25[is.na(ind25)]<-FALSE ind25 <- c(seq(1,dim(adm)[2]))[ind25] # ATE (TAKES LONG) - skip or reduce library if needed library(tmle) tresults1 <- NULL for(i in c(seq(1,dim(adm)[2]))[colnames(adm)%in%colnames(adm)[ind25]])try({ M_tmle <- suppressWarnings(tmle(as.matrix(adm[,i],ncol=1),loads3_last2, W, family="gaussian", gbound = c(0.01, 1), V=10, verbose=T, Q.SL.library = mylibrary, g.SL.library = mylibrary2)) ress <- data.frame(t(c(colnames(adm)[i],M_tmle$estimates$ATE$psi,M_tmle$estimates$ATE$CI))) tresults1 <- rbind(tresults1,ress) }) colnames(tresults1) <- c("ICD10","ATE","l95","u95") tresults1 <- data.frame(tresults1) tresults1$width <- as.numeric(paste(tresults1$u95))-as.numeric(paste(tresults1$l95)) tresults1[,2:4] <- round(apply(tresults1[,2:4],2,as.numeric),digits=4) tresults1[,5] <- round(tresults1[,5],digits=4) tresults1 <- tresults1[order(tresults1$width),] detach("package:tmle", unload=TRUE) library(mgcv) # IRR tresults2 <- NULL for(i in c(seq(1,dim(adm)[2]))[colnames(adm)%in%colnames(adm)[ind25]])try({ M_IRR <- try(mgcv::gam(formula(paste("adm[,i] ~", covariates2)), family="quasipoisson"),silent=T) p1 <- exp(M_IRR$coefficients["loads3_last2"]) p2 <- exp(M_IRR$coefficients["loads3_last2"] - qnorm(0.975)*summary(M_IRR)[[22]]["loads3_last2",2]) p3 <- exp(M_IRR$coefficients["loads3_last2"] + qnorm(0.975)*summary(M_IRR)[[22]]["loads3_last2",2]) ress <- data.frame(t(c(colnames(adm)[i],p1,p2,p3))) tresults2<- rbind(tresults2,ress) }) colnames(tresults2) <- c("ICD10","IRR","l95","u95") tresults2 <- data.frame(tresults2) tresults2$width <- as.numeric(paste(tresults2$u95))-as.numeric(paste(tresults2$l95)) tresults2[,2:4] <- round(apply(tresults2[,2:4],2,as.numeric),digits=4) tresults2[,5] <- round(tresults2[,5],digits=4) tresults2 <- tresults2[order(tresults2$width),] for(i in 1:dim(tresults2)[1]){for(j in 2:dim(tresults2)[2]){ if(as.numeric(tresults2[i,j])>1000){tresults2[i,j]<-NA} }} #eTable1 cl <- data.frame(rbind( c("S06","Intracranial injury"), c("J05","Acute obstructive laryngitis"), c("R11","Nausea and vomiting"), c("G03","Meningitis"), c("T23","Burn & corrosion of wrist and hand"), c("T22","Burn & corrosion of shoulder & upper limb"), c("T24","Burn & corrosion of lower limb"), c("A87","Viral meningitis"), c("S82","Fracture of lower leg"), c("J22","Acute lower respiratory infection"), c("Z00","General examination"), c("S42","Fracture of shoulder and upper arm"), c("T74","Child abuse, neglect & other maltreatment"), c("S52","Fracture of forearm"), c("T20","Burn & corrosion of head, face, and neck"), c("A41","Sepsis"), c("T29","Burns of multiple regions"), c("S72","Fracture of femur"), c("L02","Cutaneous abscess"), c("J44","Chronic obstructive pulmonary disease"), c("J45","Asthma"), c("R56","Convulsions"), c("Q21","Congenital malformations of cardiac septa"), c("G40","Epilepsy"), c("J18","Pneumonia"), c("S09","Unspecified injuries of head"), c("J21","Acute bronchiolitis"), c("A09","Infectious gastroenteritis"), c("T21","Burn & corrosion of trunk") )) colnames(cl) <- c("ICD10","Diagnosis") ef1 <- merge(cl,tresults1) ef2 <- merge(cl,tresults2) ef3 <- merge(ef2,ef1,by=1) ef4 <- ef3[order(ef3$IRR, decreasing=T),] ef4 <- ef4[,c("ICD10","Diagnosis.x","IRR","l95.x","u95.x","ATE","l95.y","u95.y")] colnames(ef4)<-c("ICD-10","Diagnosis","IRR","95%","CI","ATE","95%","CI") mytable1 <- xtable(ef4, align = "lr|l|ccc|ccc") print(mytable1,file=paste(getwd(),"/eTable1.tex",sep=""),table.placement="h",include.rownames=F,floating=F,comment=F) write.xlsx(ef4,file=paste(getwd(),"/eTable1.xlsx",sep="")) # FINISH and save R workspace ptm7 <- proc.time() sd2 <- round(((ptm3-ptm2)/60)[1],digits=2) sd3 <- round(((ptm4-ptm3)/60)[1],digits=2) sd4 <- round(((ptm5-ptm4)/60)[1],digits=2) sd5 <- round(((ptm6-ptm5)/60)[1],digits=2) sd6 <- round(((ptm7-ptm6)/60)[1],digits=2) sd7 <- round(((ptm7-ptm2)/60/60)[1],digits=2) cat(paste("The time for producing graphics was", sd2, "minutes \n")) cat(paste("The time for the main associative analyses was", sd3, "minutes \n")) cat(paste("The time for the main causal analyses was", sd4, "minutes \n")) cat(paste("The time for analyses by area was", sd5, "minutes \n")) cat(paste("The time for analyses by diagnosis was", sd6, "minutes \n")) cat(paste("The total analysis time was", sd7, "hours \n")) save.image(file=paste(getwd(),"/results.Rdata",sep=""))