CausalImpact source code (R) | NeuroΨ Webdesign

Click the link below to download the R script:
Download CausalImpact R Script

  o-o                  o     o-O-o                    o        O            o                    o--o                   
 /                     |       |                      |       / \           |          o         |   |                 
O     o-o o  o o-o  oo |       |   o-O-o o-o o-o o-o -o-     o---o o--o o-o | o  o o-o   o-o     O-Oo  
 \    |-| |  |  \  | | |       |   | | | | | |-| |    |      |   | |  | |-| | |  |  \  |  \      |  \ 
  o-o o-o-o--o o-o o-o-o     o-O-o o o o O-o o-o o-o  o      o   o o  o o-o o o--O o-o | o-o     o   o 
                                         |                                       |                       
                                         o                                    o--o                              

# Load necessary libraries with pacman
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  CausalImpact, bsts, zoo, magrittr, ggplot2, dplyr, tidyr, plotly, lmtest, tseries,
  heatmaply, RColorBrewer, moments
)

# Load data from URL
data <- read.csv("https://neuropsy-webdesign.de/arima/intervention.csv")
data$date <- as.Date(data$date)

# Check for duplicate dates and handle them
data <- data[!duplicated(data$date), ]

# Function to find index closest to date
find_index_by_date <- function(date) {
  idx <- which.min(abs(data$date - date))
  return(idx)
}

# Parameters
pre_start <- as.Date("2000-01-09")
pre_end <- as.Date("2020-12-26")
post_start <- as.Date("2020-12-27")
post_end <- as.Date("2021-01-03")
mcmc_iterations <- 1000
n_seasons <- 52
state_component <- "AddLocalLevel"
burn_in <- 500
ci_level <- 0.95
n_simulations <- 1000

# Create time series object
ts_data <- zoo(data[, c("mortality_rate", "pcr_rate", "intervention_rate")], order.by = data$date)

# Specify the state-space model
ss <- list()
if (state_component == "AddLocalLevel") {
  ss <- AddLocalLevel(ss, ts_data[, "mortality_rate"])
} else if (state_component == "AddLocalLinearTrend") {
  ss <- AddLocalLinearTrend(ss, ts_data[, "mortality_rate"])
} else if (state_component == "AddSeasonal") {
  ss <- AddSeasonal(ss, ts_data[, "mortality_rate"], nseasons = n_seasons)
}

# Fit the BSTS model with the additional seasonal component
bsts_model <- bsts(ts_data[, "mortality_rate"], state.specification = ss, niter = mcmc_iterations)

# Extract posterior samples for impact analysis
posterior_samples <- predict(bsts_model, burn = burn_in)

# Prepare the data for CausalImpact
impact_data <- zoo(
  cbind(
    ts_data[, "mortality_rate"],
    posterior_samples$mean,
    posterior_samples$interval[1, ],
    posterior_samples$interval[2, ]
  ),
  order.by = index(ts_data)
)

colnames(impact_data) <- c("response", "point.pred", "point.pred.lower", "point.pred.upper")

# Perform the CausalImpact analysis with the defined post-period
impact <- CausalImpact(
  impact_data,
  pre.period = c(pre_start, pre_end),
  post.period = c(post_start, post_end),
  alpha = 1 - ci_level
)

# Plot the results
plot(impact)
summary(impact)

# Create plots
impact_series <- as.data.frame(impact$series)
impact_series$index <- as.Date(row.names(impact_series))

p1 <- ggplot(impact_series, aes(x = index)) +
  geom_line(aes(y = response), color = "blue") +
  geom_line(aes(y = point.pred), color = "red", linetype = "dashed") +
  geom_ribbon(aes(ymin = point.pred.lower, ymax = point.pred.upper), fill = "red", alpha = 0.2) +
  labs(title = "Causal Impact Analysis", x = "Date", y = "Mortality Rate") +
  theme_minimal()

print(p1)

# Create interactive Plotly plots
p2 <- plot_ly(data = impact_series, x = ~index, y = ~response, type = 'scatter', mode = 'lines', name = 'Observed', line = list(color = 'blue')) %>%
  add_trace(y = ~point.pred, name = 'Predicted', line = list(color = 'red', dash = 'dash')) %>%
  add_ribbons(ymin = ~point.pred.lower, ymax = ~point.pred.upper, name = '95% CI', fillcolor = 'rgba(255, 0, 0, 0.2)', line = list(color = 'rgba(255, 0, 0, 0)')) %>%
  layout(title = 'Causal Impact Analysis', xaxis = list(title = 'Date'), yaxis = list(title = 'Mortality Rate'))

print(p2)

# Download animated *.gif
# gif_url <- "https://neuropsy-webdesign.de/arima/timeline-ani1.gif"
# download.file(gif_url, destfile = "animated_timeseries_plot.gif", mode = "wb")

100

101

102

103

104

105

106

107

108

o-o o o-O-o o O o o--o

/ | | | / \ | o | |

O o-o o o o-o oo | | o-O-o o-o o-o o-o -o- o---o o--o o-o | o o o-o o-o O-Oo

\ |-| | | \ | | | | | | | | | |-| | | | | | | |-| | | | \ | \ | \

o-o o-o-o--o o-o o-o-o o-O-o o o o O-o o-o o-o o o o o o o-o o o--O o-o | o-o o o

| |

o o--o

# Load necessary libraries with pacman

if (!require("pacman")) install.packages("pacman")

pacman::p_load(

CausalImpact, bsts, zoo, magrittr, ggplot2, dplyr, tidyr, plotly, lmtest, tseries,

heatmaply, RColorBrewer, moments

)

# Load data from URL

data <- read.csv("https://neuropsy-webdesign.de/arima/intervention.csv")

data$date <- as.Date(data$date)

# Check for duplicate dates and handle them

data <- data[!duplicated(data$date), ]

# Function to find index closest to date

find_index_by_date <- function(date) {

idx <- which.min(abs(data$date - date))

return(idx)

}

# Parameters

pre_start <- as.Date("2000-01-09")

pre_end <- as.Date("2020-12-26")

post_start <- as.Date("2020-12-27")

post_end <- as.Date("2021-01-03")

mcmc_iterations <- 1000

n_seasons <- 52

state_component <- "AddLocalLevel"

burn_in <- 500

ci_level <- 0.95

n_simulations <- 1000

# Create time series object

ts_data <- zoo(data[, c("mortality_rate", "pcr_rate", "intervention_rate")], order.by = data$date)

# Specify the state-space model

ss <- list()

if (state_component == "AddLocalLevel") {

ss <- AddLocalLevel(ss, ts_data[, "mortality_rate"])

} else if (state_component == "AddLocalLinearTrend") {

ss <- AddLocalLinearTrend(ss, ts_data[, "mortality_rate"])

} else if (state_component == "AddSeasonal") {

ss <- AddSeasonal(ss, ts_data[, "mortality_rate"], nseasons = n_seasons)

}

# Fit the BSTS model with the additional seasonal component

bsts_model <- bsts(ts_data[, "mortality_rate"], state.specification = ss, niter = mcmc_iterations)

# Extract posterior samples for impact analysis

posterior_samples <- predict(bsts_model, burn = burn_in)

# Prepare the data for CausalImpact

impact_data <- zoo(

cbind(

ts_data[, "mortality_rate"],

posterior_samples$mean,

posterior_samples$interval[1, ],

posterior_samples$interval[2, ]

order.by = index(ts_data)

)

colnames(impact_data) <- c("response", "point.pred", "point.pred.lower", "point.pred.upper")

# Perform the CausalImpact analysis with the defined post-period

impact <- CausalImpact(

impact_data,

pre.period = c(pre_start, pre_end),

post.period = c(post_start, post_end),

alpha = 1 - ci_level

)

# Plot the results

plot(impact)

summary(impact)

# Create plots

impact_series <- as.data.frame(impact$series)

impact_series$index <- as.Date(row.names(impact_series))

p1 <- ggplot(impact_series, aes(x = index)) +

geom_line(aes(y = response), color = "blue") +

geom_line(aes(y = point.pred), color = "red", linetype = "dashed") +

geom_ribbon(aes(ymin = point.pred.lower, ymax = point.pred.upper), fill = "red", alpha = 0.2) +

labs(title = "Causal Impact Analysis", x = "Date", y = "Mortality Rate") +

theme_minimal()

print(p1)

# Create interactive Plotly plots

p2 <- plot_ly(data = impact_series, x = ~index, y = ~response, type = 'scatter', mode = 'lines', name = 'Observed', line = list(color = 'blue')) %>%

add_trace(y = ~point.pred, name = 'Predicted', line = list(color = 'red', dash = 'dash')) %>%

add_ribbons(ymin = ~point.pred.lower, ymax = ~point.pred.upper, name = '95% CI', fillcolor = 'rgba(255, 0, 0, 0.2)', line = list(color = 'rgba(255, 0, 0, 0)')) %>%

layout(title = 'Causal Impact Analysis', xaxis = list(title = 'Date'), yaxis = list(title = 'Mortality Rate'))

print(p2)

# Download animated *.gif

# gif_url <- "https://neuropsy-webdesign.de/arima/timeline-ani1.gif"

# download.file(gif_url, destfile = "animated_timeseries_plot.gif", mode = "wb")