inst/doc/ITScontrol_demonstration_slope.R

In multipleITScontrol: Interrupted Time Series with a Control and Multiple Interventions

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE, 
  echo = FALSE,
  warning = FALSE,
  message = FALSE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(multipleITScontrol)
library(dplyr)
library(ggplot2)
library(lubridate)
library(stringi)
library(rlang)
library(purrr)

phei_calendar <- function(df,
                          date_column = NULL,
                          factor_column = NULL,
                          colours = NULL,
                          title = "Placeholder: Please supply title or 'element_blank()' to `title` argument",
                          subtitle = "Placeholder: Please supply subtitle or 'element_blank()' to `subtitle` argument",
                          caption = "PH.Intelligence@hertfordshire.gov.uk",
                          ncol,
                          ...) {


  date_column <- rlang::sym(date_column)
  factor_column <- rlang::sym(factor_column)

  df <- df |> dplyr::mutate(
    mon = lubridate::month(!!date_column, label = T, abbr = F),
    wkdy = weekdays(!!date_column,
                    abbreviate =
                      T
    ) |> forcats::fct_relevel("Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun"),
    day = lubridate::mday(!!date_column),
    week = stringi::stri_datetime_fields(!!date_column)$WeekOfMonth,
    year = lubridate::year(!!date_column),
    year_mon = zoo::as.yearmon(!!date_column, "%Y %m")
  ) |>
    dplyr::mutate(across(week, ~ dplyr::case_when(wkdy == "Sun" ~ week - 1,
                                           .default = as.numeric(week)
    )))
  
  df %>%
    ggplot2::ggplot(., ggplot2::aes(wkdy, week)) +
    # custom theme stuff below
    # geom_tile and facet_wrap will do all the heavy lifting
    ggplot2::geom_tile(
      alpha = 0.8,
      ggplot2::aes(fill = !!factor_column),
      color = "black", ...
    ) +
    ggplot2::facet_wrap(~year_mon, scales = "free_x", ncol = ncol) +
    ggplot2::geom_text(ggplot2::aes(label = day)) +
    # put your y-axis down, flip it, and reverse it
    ggplot2::scale_y_reverse(breaks = NULL) +
    # manually fill scale colors to something you like...
    ggplot2::scale_fill_manual(
      values = colours,
      na.value = "white",
      na.translate = FALSE
    ) +
    ggpubr::theme_pubclean() +
    ggplot2::theme(legend.position = "bottom") +
    ggplot2::labs(
      fill = "",
      x = "",
      y = "",
      title = element_blank(),
      caption = "PH.Intelligence@hertfordshire.gov.uk"
    )
}

## ----calendar, echo = FALSE, warning = FALSE, message = FALSE, fig.align="center", fig.height=10, fig.width=7, fig.retina=3----

tibble_data <- tibble::tibble(Date = seq(as.Date("2025-03-03"), as.Date("2026-08-30"), by = "day"),
  Period = dplyr::case_when(
    Date >= as.Date("2025-03-03") & Date <= as.Date("2025-08-31") ~ "Pre-intervention period",
    Date >= as.Date("2025-09-01") & Date <= as.Date("2026-03-01") ~ "Intervention 1) Reading Program",
    Date >= as.Date("2026-03-02") & Date <= as.Date("2026-08-30") ~ "Intervention 2) Peer Tutoring Sessions"
  ))

plot <- phei_calendar(
  tibble_data,
  date_column = "Date",
  "Period",
  colours = c("#3b5163", "#80bb77", "#afd0f0"),
  ncol = 3
) +
  theme(strip.text = element_text(size = rel(0.5)),
        axis.text = element_text(size = rel(0.5)),
        plot.caption = element_text(size = rel(0.5)),
        legend.text = element_text(size = rel(0.5)))
  

plot$layers[[2]]$aes_params$size <- 3

plot



## ----step_1_load_data---------------------------------------------------------
DT::datatable(its_data_school, options = list(dom = 'tip'), rownames = FALSE)

## ----echo = TRUE, results='hide'----------------------------------------------
transformed_data <- 
  multipleITScontrol::transform_data(df = its_data_school,
                                     time_var = "Date",
                                     group_var = "group_var",
                                     outcome_var = "score",
                                     intervention_dates = as.Date(c("2025-09-05", "2026-03-06")))

## -----------------------------------------------------------------------------
transformed_data

## ----echo = TRUE, results='hide'----------------------------------------------
fitted_ITS_model <- 
  multipleITScontrol::fit_its_model(transformed_data = transformed_data,
                                    impact_model = "slope",
                                    num_interventions = 2)

fitted_ITS_model

## -----------------------------------------------------------------------------
fitted_ITS_model

## ----echo = TRUE, results='hide'----------------------------------------------
my_summary_its_model <- multipleITScontrol::summary_its(fitted_ITS_model)

my_summary_its_model

## -----------------------------------------------------------------------------
my_summary_its_model

## ----echo = TRUE, results='hide'----------------------------------------------
sjPlot::tab_model(
  my_summary_its_model,
  dv.labels = "Average School Result",
  show.se = TRUE,
  collapse.se = TRUE,
  linebreak = FALSE,
  string.est = "Estimate (std. error)",
  string.ci = "95% CI",
  p.style = "numeric_stars"
)

## -----------------------------------------------------------------------------
sjPlot::tab_model(
  my_summary_its_model,
  dv.labels = "Average School Result",
  show.se = TRUE,
  collapse.se = TRUE,
  linebreak = FALSE,
  string.est = "Estimate (std. error)",
  string.ci = "95% CI",
  p.style = "numeric_stars"
)

a <- coef(my_summary_its_model)[[which(names(coef(my_summary_its_model)) == "A) Control y-axis intercept")]] |> round(2)
c <- coef(my_summary_its_model)[[which(names(coef(my_summary_its_model)) == "C) Control pre-intervention slope")]] |> round(2)
d <- coef(my_summary_its_model)[[which(names(coef(my_summary_its_model)) == "D) Pilot pre-intervention slope difference to control")]] |> round(2)
e <- coef(my_summary_its_model)[[which(names(coef(my_summary_its_model)) == "E) Control intervention 1 slope")]] |> round(2)
f <- coef(my_summary_its_model)[[which(names(coef(my_summary_its_model)) == "F) Pilot intervention 1 slope")]] |> round(2)
i <- coef(my_summary_its_model)[[which(names(coef(my_summary_its_model)) == "I) Control intervention 2 slope")]] |> round(2)
j <- coef(my_summary_its_model)[[which(names(coef(my_summary_its_model)) == "J) Pilot intervention 2 slope")]] |> round(2)

## ----echo = TRUE, results='hide'----------------------------------------------
slope_difference(model = my_summary_its_model, intervention = 1)

## ----echo = FALSE-------------------------------------------------------------
slope_difference(model = my_summary_its_model, intervention = 1)


## ----echo = TRUE, results= 'hide'---------------------------------------------
slope_difference(model = my_summary_its_model, intervention = 2)

## ----echo = FALSE-------------------------------------------------------------
slope_difference(model = my_summary_its_model, intervention = 2)


## ----echo = TRUE, results='hide'----------------------------------------------
transformed_data_with_predictions <- generate_predictions(transformed_data, fitted_ITS_model)

transformed_data_with_predictions

## -----------------------------------------------------------------------------
DT::datatable(transformed_data_with_predictions, options = list(dom = 'tip', scrollX = TRUE), rownames = FALSE)

## ----echo = TRUE, results='hide', fig.align="center", fig.width=7, fig.height=6, fig.retina=3----
its_plot(model = my_summary_its_model,
         data_with_predictions = transformed_data_with_predictions,
         time_var = "time",
         intervention_dates = as.Date(c("2025-09-05", "2026-03-06")), 
         y_axis = "Reading Comprehension Score")