Step 1: Develop the Propensity Model

The first step in developing the weights to make the population more “even” is to develop a propensity score for the treatment. Here I’ll run a logistic regression using glm() to predict where the Field Goal attempt will be iced based on the covariates that were unbalanced from before.

p_iced <- glm(is_iced ~ regulation_time_remaining + attempted_distance + 
             drive_is_home_offense + I(score_diff^2) + prior_miss + offense_win_prob, 
           data = fg_attempts, 
           family = 'binomial')

These are called propensity models because their output represent the propensity of a given attempt to get iced.

Step 2: Use the Propensity Scores to weight the non-Iced Field Goal Attempts

Then by using the augment() function from the {{broom}} package, I can add the predicted values from the model to the fg_attempts data set. The probabilities from this model can be used to re-weight the data in any number of ways. You can make adjust both the test and control to make them look like each other. You can also adjust test group to look like the control, and you can weight the control to look like the test.

In this case, since I want to understand the causal effect of icing the kicker on kicks that are likely to be iced, I’ll be re-weighting the control group to look like the test group. Thus, I will be looking for the average treatment effect among the treated (ATT) vs. the overall average treatment effect (ATE).

The formula for re-weighting the population for the ATT is:

Where p_i is the attempt’s propensity to be iced and Z_i is whether the attempt was actually iced. This winds up assigning each attempt in the test group to 1 and will upweight field goal attempts that had higher propensities for being iced from the non-iced group.

weighted_dt <- p_iced %>% 
  augment(type.predict = "response", data = fg_attempts) %>%
  mutate(
    w_att = ((.fitted * (is_iced=='Iced'))/.fitted) + 
      ((.fitted*(is_iced != 'Iced'))/(1-.fitted))
  )

Before showing the effects of the weighting let’s first look at the unweighted propensity scores:

ggplot(weighted_dt, aes(x = .fitted, fill = is_iced)) + 
  geom_density(alpha = .5) + 
  scale_x_continuous(labels = scales::percent) + 
  #scale_y_log10(labels = scales::comma) + 
  scale_fill_manual(values = c('Iced' = 'green', 'Not Iced' = 'blue')) + 
  labs(x = "P(Icing The Kicker)",
       y = "",
       title = "Probability of a FG Attempt Being Iced (Unweighted)",
       fill = "Kicker Iced?") + 
  cowplot::theme_cowplot() + 
  theme(
    legend.position = 'bottom',
    legend.justification = 'center',
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  )

This makes it very clear that the distributions of propensity to ice differs very heavily between the group that was actually iced and that which was not. Its also nice to see that the group that was actually iced generally have higher propensity scores that those that don’t.

Now let’s look at the distribution of propensity score when taking the weights into account. The distribution of the Iced group is shown in green and is unchanged from the pre and post weightings. On the bottom is the Non-Iced attempts. The overall distribution is shown in grey and the re-weighted distribution is shown in blue. Notice how it more closely reflects the distribution of the iced group.

weighted_dt %>%
  tidyr::spread(is_iced, .fitted, sep = "_p") %>%
  ggplot() +
  geom_histogram(bins = 50, aes(is_iced_pIced), alpha = 0.5) + 
  geom_histogram(bins = 50, aes(is_iced_pIced, weight = w_att), fill = "green", alpha = 0.5) + 
  geom_histogram(bins = 50, alpha = 0.5, aes(x = `is_iced_pNot Iced`, y = -..count..)) + 
  geom_histogram(bins = 50, aes(x = `is_iced_pNot Iced`, weight = w_att, y = -..count..), fill = "blue", alpha = 0.5) + 
  geom_hline(yintercept = 0, lwd = 0.5) +
  scale_y_continuous(label = abs) +
  scale_x_continuous(label = scales::percent) + 
  labs(title = "Post-Weighted Probability of FG Attempt Being Iced",
       subtitle = "grey is unweighted distribution",
       x = "P(Icing the Kicker)",
       y = "# of Attempts") + 
  theme_minimal() + 
  geom_rect(aes(xmin = 0.45, xmax = .47, ymin = 5, ymax = 100), fill = "#5DB854") + 
  geom_text(aes(x = 0.46, y = 50), label = "Iced", angle = 270, color = "white") + 
  geom_rect(aes(xmin = 0.45, xmax = .47, ymin = -100, ymax = -5), fill = "#5154B8") + 
  geom_text(aes(x = 0.46, y = -50), label = "Non-Iced", angle = 270, color = "white") + 
  coord_cartesian(ylim = c(-100, 100))

Step 3: Ensure the Post-Weighted Data is no longer Imbalanced

The next step is to run some diagnostics to make sure that the imbalance that we saw back at the beginning of this post in the standardized mean differences have gone away. I’m going to use the {{survey}} package and the tbl_svysummary() function from {{gtsummary}} to create a survey design object that incorporates the weights that were derived above. The ids = ~ 1 code tells the design object that there are no clusters in this.

svy_des <- svydesign(
  ids = ~ 1,
  data = weighted_dt,
  weights = ~ w_att
)

Then the syntax of tbl_svysummary() is identical to tbl_summary() just it uses the survey design object rather than a data frame. Like with the additional table, I’m adding in the standardized mean difference.

tbl_svysummary(
  svy_des,
  by = 'is_iced',
  include = c(regulation_time_remaining, attempted_distance, 
              drive_is_home_offense, score_diff, prior_miss, 
              offense_win_prob, is_iced),
  statistic = list(all_continuous() ~ "{mean} ({sd})",
                   all_dichotomous() ~ "{p}%"),
  digits = list(everything() ~ 2)
) %>% 
  add_difference(everything() ~ "smd")

Notice that all of the SMDs are now below the 0.1 threshold.

Another way to visualize the changes in SMDs between is using a Love Plot. Named after Dr. Thomas E. Love, the Love Plot is a way of summarizing covariate balance before and after weighting. In the first code block, I calculate both the weighted and unweighted standardized mean differences using the {{smd}} package. In the smd() code blocks, I pass in each variable, the group variable, and in the case of the weighted version, the weights.

smds <- weighted_dt %>%
  # Calculate the SMD for Each Variable
  summarise(
    # List Variables to run functions
    across(c(regulation_time_remaining, attempted_distance, 
             drive_is_home_offense, score_diff, prior_miss,
             offense_win_prob),
           # List functions
           list(
             unweighted = ~smd(.x, is_iced, na.rm = T)$estimate, 
             weighted = ~smd(.x, is_iced, w_att, na.rm = T)$estimate 
           ),
           # Assign how the naming will show up in the output
           # Assign placeholder _zzz_ to split on in the next step
           .names = "{.col}_zzz_{.fn}")
  )

smds %>% 
  pivot_longer( 
    everything(),
    values_to = "SMD", 
    names_to = c("variable", "Method"), 
    names_sep = "_zzz_"
  ) %>%
  ggplot(
    aes(x = abs(SMD), y = variable, group = Method, color = Method)
  ) +  
  geom_line(orientation = "y") +
  geom_point() + 
  geom_vline(xintercept = 0.1, color = "black", size = 0.1) + 
  labs(title = "Love Plot Pre/Post Weighting",
       subtitle = "Post-Weighted Variables are All Balanced",
       y= "") + 
  cowplot::theme_cowplot()

The Love Plot clearly shows that the weighted version of the data has corrected the imbalances that we’ve seen in the unweighted version since all variables are now below 0.1. So it looks like the propensity score weighting

Step 4: Calculate the ATT

The final step is to calculate the average treatment effect on the treatment by regressing our outcome variable (Fields Goal Made) by the “treatment” (Whether the kick was iced or not) weighted by the weighting scheme we came up with above. I’m using a linear probability model for convenience so that the coefficient is more human interpretable. But there is a case to be made for using a logistic regression since fg_made is binary.

final.model <- lm(fg_made ~ is_iced, data = weighted_dt, weights = w_att)

tidy(final.model, conf.int = T) %>%
  select(term, estimate, conf.low, conf.high) %>% 
  mutate(across(where(is.numeric), ~scales::percent(.x, accuracy = .01))) %>% 
  knitr::kable(align = 'c')

From this model the results look like that icing the kicker results in a decreased success rate of 2.85%. The confidence intervals from the linear model would suggest that its statistically significant. However, the confidence intervals generated above are overly optimistic as the weights are treated as separate individuals rather than actual weights. In order to get more robust confidence intervals, I’ll use bootstrapping to redo the entire process 1000 times. The following function does the entire process from above (propensity score -> weights -> output model).

#### Bootstrapping Estimates
fit_ipw <- function(split, ...) { 
  .df <- analysis(split) 
  
  # fit propensity score model
  propensity_model <- glm(
    is_iced ~ regulation_time_remaining + attempted_distance + 
             drive_is_home_offense + I(score_diff^2)  + prior_miss + offense_win_prob, 
    family = binomial(), 
    data = .df
  )
  
  # calculate inverse probability weights
  .df <- propensity_model %>% 
    augment(type.predict = "response", data = .df) %>%
    mutate(
      w_att = ((.fitted * (is_iced=='Iced'))/.fitted) + 
      ((.fitted*(is_iced != 'Iced'))/(1-.fitted))
    )
  
  # fit correctly bootstrapped ipw model
  lm(fg_made ~ is_iced, data = .df, weights = w_att) %>%
    tidy()
}

The bootstrapping will be done using the {{rsample}} package and the bootstraps() function. In the function I ask for 1,000 bootstrapped samples (the apparent option includes a 1001st sample that’s the entire data set). Then I apply the above function to every bootstrapped sample through {{purrr}}’s map() function.

# fit ipw model to bootstrapped samples
set.seed(20220130)
ipw_results <- bootstraps(fg_attempts, 1000, apparent = TRUE) %>% 
  mutate(results = map(splits, fit_ipw))

Finally, the int_t() function generates confidence intervals from the t-distribution based on the results of the 1,000 bootstrapped samples.

# get t-statistic-based CIs
int_t(ipw_results, results) %>%
  filter(term == "is_icedIced") %>% 
  select(term, .lower, .estimate, .upper) %>% 
  mutate(across(where(is.numeric), ~scales::percent(.x, accuracy = .01))) %>%
  knitr::kable(align = 'c')

From the bootstrapped results, we have the a similar point estimate of -2.82% which is much smaller than the 8.3% that was seen in the naive analysis. but the confidence intervals now spans from -5.88% to 0.50% making the results not significantly different from zero.

So in conclusion, we can’t definitively say that icing the kicker is actually harmful to the kicker’s success which seems consistent with the other studies that say that either its ineffective or only mildly effective at best.

In the next post in this series, I’ll be looking at alternative causal inference methodologies like G-computation and targeted maximum likelihood estimation (TLME) to see if the results are similar or different to the results from this post.