####################################### #TRANSFER FUNCTIONS: REGRESSION TABLES# ####################################### #loading libraries library(TSA) library(doBy) library(stargazer) library(lmtest) library(plyr) library(forecast) library(reshape) #setting wd (NOTE: CHANGE THIS TO THE DIRECTORY OF THE REPLICATION DATA ON YOUR MACHINE) setwd("C:/Users/mlc872/Dropbox/Terrorism/Natural_Experiments/Daylight saving/Epidemiology/Revision/Data") ############## #OPENING DATA# ############## full <- read.csv("full_data_v2.csv", as.is=T) ############### #ARIMAX TABLES# ############### #TABLE 2 #dependent variables dep <- ts(full$depress_count, frequency = 365.25/7) mani.bip <- ts(full$manic_bipolar_count, frequency = 365.25/7) #transfer dataframes from <- data.frame(full$from_clock_change_week) to <- data.frame(full$to_clock_change_week) both <- data.frame(full$from_clock_change_week , full$to_clock_change_week) #running models fit1 <- arimax(log(dep),order=c(1,1,1), xtransf=both, transfer=list(c(1,0),c(1,0))) fit2 <- arimax(log(mani.bip),order=c(1,1,1), xtransf=both, transfer=list(c(1,0),c(1,0))) #building table and writing to html stargazer(coeftest(fit1),coeftest(fit2), type="html", title = "The effect of daylight saving time on affective disorders", dep.var.caption='Dependent variables in logs', dep.var.labels=c("Depression", "Bipolar disorder"), column.labels=c("ARIMA(1,1,1)", "ARIMA(1,1,1)"), covariate.labels = c("Summertime (AR1)","Summertime (MA0)", "Wintertime (AR1)", "Wintertime(MA0)"), add.lines=list(c("Obs.", sum(!is.na(fit1$residuals))-fit1["arma"][[1]][6], sum(!is.na(fit2$residuals))-fit2["arma"][[1]][6]), c("Log-likelihood",round(fit1$loglik,d=2), round(fit2$loglik,d=2)), c("AIC", round(fit1$aic,d=2), round(fit2$aic,d=2)) ), ci = T, report="vcs", multicolumn = F, notes.append = F, omit=seq(from=1, to=2, by=1), notes = ("ARMA-coefficients omitted"), out="daylight_saving_tab2.htm") #TABLE A1 #running models fit0 <- arimax(log(dep),order=c(1,1,1), xtransf=both, transfer=list(c(1,0),c(1,0))) fit1 <- arimax(log(dep),order=c(1,1,1), xreg = fourier(log(dep), K = 1), xtransf=both, transfer=list(c(1,0),c(1,0))) fit2 <- arimax(log(dep),order=c(1,1,1), xreg = fourier(log(dep), K = 2), xtransf=both, transfer=list(c(1,0),c(1,0))) fit3 <- arimax(log(dep),order=c(1,1,1), xreg = fourier(log(dep), K = 3), xtransf=both, transfer=list(c(1,0),c(1,0))) fit4 <- arimax(log(dep),order=c(1,1,1), xreg = fourier(log(dep), K = 4), xtransf=both, transfer=list(c(1,0),c(1,0))) fit5 <- arimax(log(dep),order=c(1,1,1), xreg = fourier(log(dep), K = 26), xtransf=both, transfer=list(c(1,0),c(1,0))) #building table and writing to html stargazer(coeftest(fit0),coeftest(fit1), coeftest(fit2),coeftest(fit3), coeftest(fit4),coeftest(fit5), type="html", title = "The effect of daylight saving time transitions, including k Fourier terms", dep.var.caption='Dependent variable: log(major depressive episodes)', dep.var.labels=c("k = 0 (original)", "k = 1", "k = 2", "k = 3", "k = 4", "k = 26"), column.labels=c("ARIMA(1,1,1)", "ARIMA(1,1,1)", "ARIMA(1,1,1)", "ARIMA(1,1,1)", "ARIMA(1,1,1)"), covariate.labels = c("Summertime (AR1)","Summertime (MA0)", "Wintertime (AR1)", "Wintertime(MA0)"), add.lines=list(c("Obs.", sum(!is.na(fit0$residuals))-fit0["arma"][[1]][6], sum(!is.na(fit1$residuals))-fit1["arma"][[1]][6], sum(!is.na(fit2$residuals))-fit2["arma"][[1]][6], sum(!is.na(fit3$residuals))-fit3["arma"][[1]][6], sum(!is.na(fit4$residuals))-fit4["arma"][[1]][6], sum(!is.na(fit5$residuals))-fit5["arma"][[1]][6]), c("Log-likelihood", round(fit0$loglik,d=2), round(fit1$loglik,d=2), round(fit2$loglik,d=2), round(fit3$loglik,d=2), round(fit4$loglik,d=2), round(fit5$loglik,d=2)), c("AIC", round(fit0$aic,d=2), round(fit1$aic,d=2), round(fit2$aic,d=2), round(fit3$aic,d=2), round(fit4$aic,d=2), round(fit5$aic,d=2))), ci = T, report="vcs", multicolumn = F, notes.append = F, omit=c("ar[0-9]","ma[0-9]","S[0-9]+\\.","C[0-9]+\\."), notes = ("ARMA-coefficients and Fourier terms omitted"), out="daylight_saving_taba1.htm") #TABLE A2 #conditions variables precip <- ts(full$precipitation, frequency = 365.25/7) temp <- ts(full$temperature, frequency = 365.25/7) light <- ts(full$daylight, frequency = 365.25/7) #running models fit1 <- arimax(log(dep),order=c(1,1,1), xreg = precip, xtransf=both, transfer=list(c(1,0),c(1,0))) fit2 <- arimax(log(dep),order=c(1,1,1), xreg = temp, xtransf=both, transfer=list(c(1,0),c(1,0))) fit3 <- arimax(log(dep),order=c(1,1,1), xreg = light, xtransf=both, transfer=list(c(1,0),c(1,0))) stargazer(coeftest(fit1), coeftest(fit2), coeftest(fit3), type="html", title = "Controlling for weather and minutes of daylight", dep.var.caption='Dependent variable: log(major depressive episodes)', dep.var.labels=c("Precipitation", "Temperature", "Min. of daylight"), column.labels=c("ARIMA(1,1,1)", "ARIMA(1,1,1)","ARIMA(1,1,1)"), covariate.labels = c("Conditions","Summertime (AR1)","Summertime (MA0)", "Wintertime (AR1)", "Wintertime(MA0)"), add.lines=list(c("Obs.", sum(!is.na(fit1$residuals))-fit1["arma"][[1]][6], sum(!is.na(fit2$residuals))-fit2["arma"][[1]][6], sum(!is.na(fit3$residuals))-fit3["arma"][[1]][6]), c("Log-likelihood", round(fit1$loglik,d=2), round(fit2$loglik,d=2), round(fit3$loglik,d=2)), c("AIC", round(fit1$aic,d=2), round(fit2$aic,d=2), round(fit3$aic,d=2))), ci = T, report="vcs", multicolumn = F, notes.append = F, omit=seq(from=1, to=2, by=1), notes = ("ARMA-coefficients omitted"), out="daylight_saving_taba2.htm") ############# ###FIGURES### ############# #INTERVENTION EFFECT #running models fit1 <- arimax(log(dep),order=c(1,1,1), xtransf=both, transfer=list(c(1,0),c(1,0))) intv.effect <- 1*(seq_along(dep)== 1) intv.effect <-ts( filter(intv.effect, filter = coef(fit1)["full.to_clock_change_week-AR1"], method ='rec',sides = '1' ) *coef(fit1)["full.to_clock_change_week-MA0"]) intv.effect = exp(intv.effect) tsp(intv.effect) <- tsp(dep) #win.metafile("figure3.emf", width = 12) plot(100*(intv.effect-1)[1:40],type='l', main = "Figure 3. The depressogenic effect of the transition from summer time to standard time", ylab="% change in the number of cases",xlab="Week from transition to standard time") abline(h = 0, lty = 2) #dev.off() #RAW TIME SERIES #time resetting dep <- ts(full$depress_count) mani.bip <- ts(full$manic_bipolar_count) #depression par(family="serif") plot(dep, xaxt = "n", xlab = "Year", ylab = "Average no. of daily incidents in a week", main = "Figure A3: Incidence rate for major depressive episodes, 1995-2012") axis(1, at= append(seq(from = 1, to = 940, by = 52.18),940), labels=seq(from = 1995, to = 2013, by = 1)) #bipolar disorder par(family = "serif") plot(mani.bip, xaxt = "n", xlab = "Year", ylab = "Average no. of daily incidents in a week", main = "Figure A4: Incidence rate for bipolar disorder, 1995-2012") axis(1, at= append(seq(from = 1, to = 940, by = 52.18),940), labels=seq(from = 1995, to = 2013, by = 1)) #time resetting dep <- ts(full$depress_count, frequency = 365.25/7) mani.bip <- ts(full$manic_bipolar_count, frequency = 365.25/7) #PLACEBO TESTS #Placebo test with constant p, d, and q #No seasonal corrections #setting up vectors coef <- NULL lower <- NULL upper <- NULL int_time <- NULL y <- 1 for (z in seq(-45,45, by=5)) { cat("Week of intervention: ", z, "\n") int_time[y]<-z to_temp <- numeric(nrow(to)) for (x in 1:length(which(to==1))) { i <- which(to==1)[x]+z to_temp[i]<-1 } to_temp <- to_temp[1:940] fit1 <- arimax(log(dep),order=c(1,1,1), xtransf=to_temp, transfer=list(c(1,0))) cat("Coeff. on MA0: ", coef(fit1)["T1-MA0"],"\n") coef[y]<-coef(fit1)["T1-MA0"] cat("Lower: ",confint(fit1)["T1-MA0",1],"\n") lower[y]<-confint(fit1)["T1-MA0",1] cat("Upper: ",confint(fit1)["T1-MA0",2],"\n") upper[y]<-confint(fit1)["T1-MA0",2] y <- y+1 } #plotting results par(family = "serif") plot(int_time, coef,type='n', ylim=c(min(lower[is.finite(lower)]),max(upper[is.finite(upper)])), ylab="Coefficient on intervention pulse", xlab="Intervention time (weeks from actual intervention)", main = "Figure A1: Placebo Interventions", cex.lab=1.5, cex.axis=1.5, cex.main=1.5) polygon(c(int_time,rev(int_time)),c(lower,rev(upper)), col="grey",border=NA) lines(int_time, coef,type='l') abline(a=0,b=0,lty=3,col='grey50') abline(v=0, lty=3,col='grey50') legend(5,-0.7,"95 % confidence intervals", density= 10000 ,bty="n", col="grey", fill= "grey", border ="NA", cex=1.3, pt.cex=1.3) #Placebo test with constant p, d, and q #26 Fourier paris correcting for seasonality #setting up vectors coef <- NULL lower <- NULL upper <- NULL int_time <- NULL se <- NULL y <- 1 for (z in seq(-45,45, by=5)) { cat("Week of intervention: ", z, "\n") int_time[y]<-z to_temp <- numeric(nrow(to)) for (x in 1:length(which(to==1))) { i <- which(to==1)[x]+z to_temp[i]<-1 } to_temp <- to_temp[1:940] fit1 <- arimax(log(dep),order=c(1,1,1), xreg = fourier(log(dep), K = 26), xtransf=to_temp, transfer=list(c(1,0))) cat("Coeff. on MA0: ", coef(fit1)["T1-MA0"],"\n") coef[y]<-coef(fit1)["T1-MA0"] cat("Lower: ",confint(fit1)["T1-MA0",1],"\n") lower[y]<-confint(fit1)["T1-MA0",1] cat("Upper: ",confint(fit1)["T1-MA0",2],"\n") upper[y]<-confint(fit1)["T1-MA0",2] se[y] <- sqrt(fit1$var.coef[ncol(fit1$var.coef),nrow(fit1$var.coef)]) y <- y+1 } #interpolating NAN's nan_time <- int_time[which(is.nan(upper))] upper_ipo <- na.approx(upper) lower_ipo <- na.approx(lower) #plotting results jumps <- 5 par(family="serif") plot(int_time, coef,type='n', ylim=c(min(lower[is.finite(lower)]),max(upper[is.finite(upper)])), ylab="Coefficient on intervention pulse", xlab="Intervention time (weeks from actual intervention)", main = "Placebo Interventions, 26 Fourier terms", cex.lab=1.5, cex.axis=1.5, cex.main=1.5) polygon(c(int_time,rev(int_time)),c(lower_ipo,rev(upper_ipo)), col="grey",border=NA) for (i in 1:length(nan_time)) { #making x-points x.points <- c(int_time[int_time <= nan_time[i]+jumps & int_time >= nan_time[i]-jumps]) x.points[1] <- mean(x.points[1:2]) x.points[length(x.points)] <- mean(x.points[(length(x.points)-1):length(x.points)]) x.points.rev <- rev(x.points) #making y-points y.points <- c(lower_ipo[int_time <= nan_time[i]+jumps & int_time >= nan_time[i]-jumps]) y.points[1] <- mean(y.points[1:2]) y.points[length(y.points)] <- mean(y.points[(length(y.points)-1):length(y.points)]) y.points.rev <- rev(c(upper_ipo[int_time <= nan_time[i]+jumps & int_time >= nan_time[i]-jumps])) y.points.rev[1] <- mean(y.points.rev[1:2]) y.points.rev[length(y.points.rev)] <- mean(y.points.rev[(length(y.points.rev)-1):length(y.points.rev)]) #drawing white polygon polygon(c(x.points,x.points.rev), c(y.points, y.points.rev), col="white" ,border=NA) #crossing white polygon out polygon(c(x.points,x.points.rev), c(y.points, y.points.rev), col="grey50", density = 20 ,border=NA) polygon(c(x.points,x.points.rev), c(y.points, y.points.rev), col="grey50", density = 20 ,border=NA,angle=135) } lines(int_time, coef,type='l') abline(a=0,b=0,lty=3,col='grey50') abline(v=0, lty=3,col='grey50') legend(5,-0.15,c("95 % confidence intervals","Standard errors not calculated, \nconfidence intervals invalid"), angle=c(90,135),density=c(10000,40),bty="n", col="grey", fill=c("grey","grey"), border =c("NA","NA"), cex=1-3, pt.cex=1.3) legend(5,-0.15,c("95 % confidence intervals","Standard errors not calculated, \nconfidence intervals invalid"), angle=c(90,45),density=c(10000,40),bty="n", col="grey", fill=c("grey","grey"), border =c("NA","NA"), cex=1.3, pt.cex=1.3) ##################################### #COMPARING ALL THE ESTIMATES BY RANK# ##################################### #discarding nan's se2 <- se[!is.nan(se)] coef2 <- coef[!is.nan(se)] #sampling from each distribution draws <- 10000 for (i in 1:length(se2)) { s <- rnorm(draws, coef2[i], se2[i]) if (i > 1) { sim <- cbind(sim, s) } else { sim <- data.frame(s) } names(sim)[names(sim) == "s"] <- paste0("s",i) } #s9 is the true intervention here orders <- apply(sim, 1, function(x) rank(-x)[9]) #testing whether this is different from a uniform distribution chi <- chisq.test(table(orders)) chi.text <- ifelse(chi$p.value==0, ) #adding wilcoxon test sim2 <- melt(sim) sim2$true <- as.numeric(sim2$variable=="s9") w1 <- wilcox.test(sim2$value[sim2$true==1],sim2$value[sim2$true==0]) w2 <- wilcox.test(sim2$value[sim2$true==1],sim2$value[sim2$true==0], alternative = "g") #creating bar plot #counts counts <- table(orders) while (length(counts) < 17) { counts <- append(counts,0) } names(counts) <- 1:17 #plot itself png(file="bar_rank.png",width=5000, height=4500, res = 500) par(family="serif") bars <- barplot(counts/sum(counts), xlab = "Rank", ylab = "Proportion", axes = T, axis.lty = 1) title("efigure 3: Distribution of ranks of the true intervention", line = 2) text(x = bars, y = counts/sum(counts)+0.02, labels=as.character(round(counts/sum(counts),digits = 2)), xpd=TRUE) mtext(paste0("(",draws, " draws, Pearson test for uniform distribution, \nMann-Whitney test for identical distributions)"), cex = 0.9) if (chi$p.value==0) { text(x = 15, y = 0.3, bquote(atop(chi^2 == .(chi$statistic), italic(p) < 0.001))) } else { text(x = 15, y = 0.3, bquote(atop(chi^2 == .(chi$statistic), italic(p) == .(chi$p.value)))) } if (w2$p.value==0) { text(x = 15, y = 0.2, bquote(atop(italic(U) == .(w2$statistic), italic(p) < 0.001))) } else { text(x = 15, y = 0.2, bquote(atop(italic(U) == .(w2$statistic), italic(p) == .(w2$p.value)))) } dev.off() #drawing plot for overlapping distributions means <- summaryBy(value ~ variable, data = sim2) dens <- apply(sim, 2, function(x) max(density(x)[[2]])) md <- cbind(means, dens) cols <- append(rep("grey", 8), "green") cols <- append(cols, rep("grey", 8)) windowsFonts(Times=windowsFont("TT Times New Roman")) p1 <- ggplot(sim2, aes(value, fill=variable)) + geom_density(alpha=.4, color=NA) + ggtitle("efigure4: Simulated estimate distributions") + xlab("Estimate") + ylab("Density") + geom_segment(data = md, aes(x = value.mean, y = 0, xend = value.mean, yend = dens), linetype = 2, color = "white") + scale_fill_manual(values = cols, breaks = c("s1","s9"), labels = c("Placebo interventions","True intervention")) + theme_classic() + theme(legend.title=element_blank(), legend.position=c(0.9,0.75), plot.title = element_text(face="bold", family = "Times"), axis.text.x = element_text(family = "Times"), axis.text.y = element_text(family = "Times"), axis.title.x = element_text(family = "Times"), axis.title.y = element_text(family = "Times"), legend.text = element_text(family = "Times")) p1 png(file="densities.png",width=5000, height=4500, res = 500) par(mar=c(5,3,2,2)+0.1) par(family="serif") p1 dev.off()