Part-III-code-block-4.R

library(dplyr)
library(ggplot2)
library(splines)

days_grid <- 0:180
interior_knots <- c(2, 6, 25)
#interior_knots <- c(6, 25, 80)

my_spline <- ns(days_grid, knots = interior_knots)

perf_hat <- c()
eta_df <- data.frame()

for (j in 1:5) {
  subject_df <- subjects_df %>% filter(subject == j)
  training <- subject_df$training

  z_vars <- list()
  for (n in 1:nrow(subject_df)) {
    spline_pred  <- predict(my_spline, newx = (n - 1):1)
    spline_vars  <- colSums(spline_pred * training[1:(n - 1)]) # convolution
    spline_const <- sum(training[1:(n - 1)])
    z_vars[[n]]  <- c(spline_const, spline_vars)
  }

  z_vars_df <- Reduce(rbind.data.frame, z_vars)
  names(z_vars_df) <- paste0("z_", 1:ncol(z_vars_df))

  subject_aug_df <- cbind(subject_df, z_vars_df)

  spline_reg <- lm(perf_cp ~ z_1 + z_2 + z_3 + z_4 + z_5, data = subject_aug_df)
  summary(spline_reg)

  subject_aug_df$perf_hat <- predict(spline_reg, subject_aug_df)
  perf_hat <- c(perf_hat, subject_aug_df$perf_hat)

  spline_vars_grid <- cbind(1, predict(my_spline, newx = days_grid))
  eta <- as.numeric(spline_vars_grid %*% coef(spline_reg)[-1])
  eta_df <- rbind(eta_df, data.frame(subject = j, day = days_grid, eta = eta))

}
subjects_df$perf_hat_spline <- perf_hat

ggplot(eta_df, aes(x = day, y = eta)) +
  geom_line(color = "blue", size = 1.5) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "green", size = 1) +
  geom_point(data = data.frame(x = interior_knots), aes(x = x, y = 0),
             shape = 10, size = 3, color = "brown") +
  facet_grid(subject ~ .) +
  ggtitle("Spline-based estimation of impulse response \u03D5(t)") +
  xlab("Day (n)") + ylab("Lag distribution value") +
  theme(text = element_text(size = 16))

ggplot(subjects_df, aes(x = day, y = perf_cp)) +
  geom_col(data = subjects_df, aes(y = 10 * training), color = "grey",
           width = .2) +
  geom_point() +
  geom_line(aes(y = perf_hat), color = "blue", size = 1) +
  geom_line(aes(y = perf_hat_spline), color = "dark orange", size = 1) +
  facet_grid(subject ~ .) +
  annotate("text", label = "Relative training intensities (x10)",
           x = 70, y = 25, color = "grey32") +
  ggtitle("Adding convolution-based performance prediction") +
  xlab("Day (n)") + ylab("Performance Scale") +
  theme(text = element_text(size = 16))