Introduction to the simITS package
In simITS: Analysis via Simulation of Interrupted Time Series (ITS) Data

knitr::opts_chunk$set(
  collapse = TRUE,
  fig.align='center',
  comment = "#>"
)
library( ggplot2 )
library( simITS )

Introduction

This vignette quickly outlines the primary method calls for conducting an analysis of an Interrupted Time Series using the simulation approach proposed in the companion paper.

We first cover a simple regression model, then show how to do smoothing, then seasonality. We also make a brief note about generating fake data for the purposes of conducting simulation studies.

Basic ITS analysis

We use the raw Mecklenberg data to illustrate the simITS package.

data(mecklenberg)
head( mecklenberg )
meck = mutate( mecklenberg, pbail = 100 * pbail )
ggplot( meck, aes( x=month, y=pbail)) +
  geom_rect(aes( ymin=-Inf, ymax=Inf, xmin=0.5, xmax=25, fill="lightgray"), col = "lightgray", alpha=0.25) +
  scale_fill_identity(name = "", guide = "none", labels = c('Post Policy era')) +
  geom_hline( yintercept = 0, col="black") +
  geom_line( col="black", lty=1, lwd=0.5) +
  geom_point() +
  scale_x_continuous( breaks = c(-29,-24,-18,-12,-6,1,6,12,18,24)) +    
  coord_cartesian(xlim=c(-29.5,24.5), ylim=c(0,100), expand=FALSE) +
  labs( title = " ", y = "Percent cases assigned bail", x = " " )

To have autoregressive errors we use lagged outcomes. We can add lagged outcomes (and covariates) as so:

meck = add_lagged_covariates( meck, outcomename = "pbail", covariates=NULL )
sample_n( meck, 5 ) %>% arrange( month )

This package passes functions for fitting the model, and then uses these functions for doing the extrapolation. For the default, we use the package's fit_model_default() which is a simple line (with lagged outcome as a covariate):

meck.pre = filter( meck, month <= 0 )    
mod = fit_model_default( meck.pre, "pbail" )
summary( mod )

To run the entire simulation and extrapolation as a call, we can directly do:

t0 = 0
envelope = process_outcome_model( "pbail", meck, 
                                  t0=t0, R = 100, 
                                  summarize = TRUE, smooth=FALSE )
sample_n( envelope, 5 ) %>% arrange( month )

And plotting our results:

ggplot( envelope, aes( month ) ) +
  geom_line( aes( y=Y ), alpha = 0.6 ) +  # original data
  geom_point( aes( y=Y ) ) +              # original data
  geom_ribbon( aes( ymin=Ymin, ymax=Ymax ), alpha=0.2 ) +
  geom_line( aes( y = Ystar ), col="darkgrey" ) +
  geom_vline( xintercept = t0+0.5)

We provide a nice utility function to generate these graphs:

make_envelope_graph(envelope, t0=t0)

Testing and Impact Intervals

We can aggregate impacts for several time points as follows. First call process_outcome_model() without summarizing:

predictions = process_outcome_model( "pbail", meck,
                                     t0=t0, R = 100,
                                     summarize = FALSE, smooth=FALSE )

Then use aggregate_simulation_results():

sstat = aggregate_simulation_results( orig.data = meck, outcomename = "pbail",
                                      predictions = predictions, months = 1:18 )

quantile( sstat$t, c( 0.025, 0.975 ))
sstat$t.obs
sstat$t.obs - quantile( sstat$t, c( 0.025, 0.975 ))

Generating fake data

For simulation we also offer a fake data generator. It works like this:

dat = generate_fake_data( t_min=-60, t_max=18, t0 = 0 )
qplot( month, Y, data=dat, geom = c( "point","line") )

Smoothing

Here we demonstrate summarizing and smoothing, using the fake data we just generated.

envelope = process_outcome_model( "Y", dat, t0=t0, R = 100, 
                                  summarize = TRUE, smooth=TRUE )
make_envelope_graph(envelope, t0 )

We can smooth to different degrees using the smooth_k parameter:

alphas = c( 6, 11, 20, 100 )
preds = purrr::map( alphas, function( alpha ) {
  pds = process_outcome_model( "Y", dat,
                               t0=t0, R = 20,
                               summarize = FALSE, smooth=TRUE,
                               smooth_k = alpha )
  pds
} )
names( preds ) = alphas
preds = bind_rows( preds, .id="alpha_k" )
ggplot( filter( preds, month >= t0 ), aes( month, Ysmooth ) ) +
  facet_wrap( ~ alpha_k ) +
  geom_line( aes( group=Run, col=alpha_k ), alpha=0.5, na.rm=TRUE) +
  geom_line( data=dat, aes( month, Y ), col="black", alpha=0.5 ) +
  geom_vline( xintercept=t0, col="red" ) +
  labs( x="month", y="proportion given bail")

Seasonality and covariates

A seasonality model on some fake data with a strong seasonality component is easy to fit. You just construct some code to fit the seasonality model via the make_fit_season_model() factory (you need to have the covariates pre-constructed in your data):

data( newjersey )
fit_season_model_qtemp =  make_fit_season_model( ~ temperature + Q2 + Q3 + Q4 )

envelope = process_outcome_model( "n.warrant", newjersey, t0=-7, R = 100, 
                                  summarize = TRUE, smooth=TRUE, 
                                  fit_model = fit_season_model_qtemp )
make_envelope_graph( envelope, t0=-7 )

Note how it will construct the lagged covariates automatically. The make_fit_season_model() method records what covariates are needed from the passed formula.

Smoothing and seasonality

We can smooth around a seasonality model either with a default smoother made from the specified seasonality model (as was done above) or, like the following, with a specified one of your choice:

smoother = make_model_smoother( fit_model = fit_season_model_sin, covariates = newjersey )
envelope_sin = process_outcome_model( "n.warrant", newjersey, t0=-7, R = 100,
                                  summarize = TRUE, smooth=TRUE, smoother = smoother, smooth_k = 11,
                                  fit_model = fit_season_model_qtemp )
envelope_sin$Ysmooth.base = envelope$Ysmooth
envelope_sin$Ysmooth1.base = envelope$Ysmooth1
make_envelope_graph( filter( envelope_sin, month > -30 ), t0=-7 ) +
  geom_line( aes( y=Ysmooth.base ), col="blue", na.rm=TRUE ) +
  geom_line( aes( y=Ysmooth1.base ), col="blue", lty=2, na.rm=TRUE )