Updates

Author: ngreifer

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: apm
 Title: Averaged Prediction Models
-Version: 0.0.0.9002
+Version: 0.1.0
 Authors@R: 
   c(person("Thomas", "Leavitt", email = "[email protected]", role = c("aut"),
            comment = c(ORCID = "0000-0002-3668-6409")),

R/apm_est.R

@@ -113,19 +113,15 @@ apm_est <- function(fits, post_time, M = 0, R = 1000L, all_models = FALSE, cl =
   BMA_weights <- fits$BMA_weights
 
   #Remove models that won't contribute
-  if (!all_models) {
+  if (all_models) {
+    models_to_keep <- seq_along(models)
+    fits_to_keep <- seq_along(fits$val_fits)
+  }
+  else {
     models_to_keep <- which(BMA_weights > 0)
     fits_to_keep <- which(grid$model %in% models_to_keep)
 
-    # models <- models[models_to_keep]
     BMA_weights <- BMA_weights[models_to_keep]
-    
-    # grid <- grid[fits_to_keep, , drop = FALSE]
-    # fits$val_fits <- fits$val_fits[fits_to_keep]
-  }
-  else {
-    models_to_keep <- seq_along(models)
-    fits_to_keep <- seq_along(fits$val_fits)
   }
 
   #Prep everything for bootstrap that doesn't involve weights
@@ -163,10 +159,8 @@ apm_est <- function(fits, post_time, M = 0, R = 1000L, all_models = FALSE, cl =
       ti <- grid[["time_ind"]][fi]
 
       d <- mods[[mi]]$data
-      
-      time <- val_times[ti]
-      
-      .subset_f_post_list[[fi]] <- which(d[[time_var]] == time)
+
+      .subset_f_post_list[[fi]] <- which(d[[time_var]] == val_times[ti])
 
       .val_data_f_val_list[[fi]] <- d[.subset_f_post_list[[fi]], , drop = FALSE]
 
@@ -194,10 +188,6 @@ apm_est <- function(fits, post_time, M = 0, R = 1000L, all_models = FALSE, cl =
     for (mi in models_to_keep) {
       model <- models[[mi]]
 
-      d <- mods[[mi]]$data
-      
-      time <- post_time
-      
       .subset_mi <- .subset_m_post_list[[mi]]
 
       .val_data_mi <- .val_data_m_post_list[[mi]]
@@ -260,14 +250,8 @@ apm_est <- function(fits, post_time, M = 0, R = 1000L, all_models = FALSE, cl =
       mi <- grid[["model"]][fi]
       ti <- grid[["time_ind"]][fi]
 
-      d <- mods[[mi]]$data
-      
-      time <- val_times[ti]
-      
       .subset_fi <- .subset_f_post_list[[fi]]
 
-      .val_data_fi <- .val_data_f_val_list[[fi]]
-      
       .val_weights_fi <- .weights[.subset_fi] * weights[.subset_fi]
 
       .val_groups_fi <- .val_groups_f_val_list[[fi]]

R/apm_pre.R

@@ -227,7 +227,7 @@ apm_pre <- function(models, data, weights = NULL, group_var, time_var,
   }
 
   #Difference in average prediction errors
-  pred_error_diffs_mat <- pred_errors_array[, , "1"] - pred_errors_array[, , "0"]
+  pred_error_diffs_mat <- pred_errors_array[, , 2L] - pred_errors_array[, , 1L]
 
   #Simulate to get BMA weights
 
@@ -277,15 +277,12 @@ apm_pre <- function(models, data, weights = NULL, group_var, time_var,
         p <- exp(p)
       }
 
-      predicted_val_means_s_i <- setNames(
-        vapply(group_levels, function(g) {
+      predicted_val_means_s_i <- vapply(group_levels, function(g) {
           .wtd_mean(p, val_weights[[f]], val_groups[[f]][[g]])
-        }, numeric(1L)),
-        group_levels
-      )
+        }, numeric(1L))
 
-      mat[ti, mi] <- (observed_val_means[[val_time_c]]["1"] - observed_val_means[[val_time_c]]["0"]) -
-        (predicted_val_means_s_i["1"] - predicted_val_means_s_i["0"])
+      mat[ti, mi] <- (observed_val_means[[val_time_c]][2L] - observed_val_means[[val_time_c]][1L]) -
+        (predicted_val_means_s_i[2L] - predicted_val_means_s_i[1L])
     }
 
     mat

R/plot.apm_pre_fits.R

@@ -14,15 +14,15 @@
 #' A `ggplot` object, which can be manipulated using `ggplot2` syntax (after loading `ggplot2`).
 #' 
 #' @details
-#' When `type = "weights"`, `plot()` displays a bar plot with a bar for each model with height equal to the BMA weight/posterior probability of selection for that model. (Note that the plot margins can sometimes cut off the models names; use `theme(plot.margins =)` after loading `ggplot2` to extend the left margin of the plot to ensure all text is visible. Alternatively, the axis text can be rotated using `theme(axis.text.x =)`.)
+#' When `type = "weights"`, `plot()` displays a bar plot with a bar for each model with height equal to the BMA weight/posterior probability of selection for that model. (Note that the plot margins can sometimes cut off the model names; use `theme(plot.margins =)` after loading `ggplot2` to extend the left margin of the plot to ensure all text is visible. Alternatively, the axis text can be rotated using `theme(axis.text.x =)`.)
 #' 
-#' When `type = "errors"`, `plot()` displays a lattice of bar plots, with a plot for each model displaying the difference in average prediction errors for each validation period. The period with the largest difference in average prediction errors will be shaded black. The model with the smallest maximum absolute difference in average prediction errors will have a gray label.
+#' When `type = "errors"`, `plot()` displays a lattice of bar plots with a plot for each model displaying the difference in average prediction errors for each validation period. The period with the largest difference in average prediction errors will be shaded black. The model with the smallest maximum absolute difference in average prediction errors will have a gray label.
 #' 
-#' When `type = "predict"`, `plot()` displays a lattice of line plots, with a plot for each model displaying the observed and predicted outcomes for each validation period under each model. The observed outcomes are displayed as points, while the predicted outcomes are displayed as lines.
+#' When `type = "predict"`, `plot()` displays a lattice of line plots with a plot for each model displaying the observed and predicted outcomes for each validation period under each model. The observed outcomes are displayed as points, while the predicted outcomes are displayed as lines.
 #' 
-#' When `type = "corrected"`, `plot()` displays a lattice of line plots, with a plot for each model displaying the observed and corrected predictions for the treated group for each validation period under each model. The observed outcomes are displayed as points, while the corrected predictions are displayed as lines. Corrected predictions are computed as the observed outcome in the treated group minus the prediction error in the treated group plus the prediction error in the control group.
+#' When `type = "corrected"`, `plot()` displays a lattice of line plots with a plot for each model displaying the observed and corrected predictions for the treated group for each validation period under each model. The observed outcomes are displayed as points, while the corrected predictions are displayed as lines. Corrected predictions are computed as the observed outcome in the treated group minus the prediction error in the treated group plus the prediction error in the control group.
 #' 
-#' @seealso [apm_pre()] to to compute the difference in average prediction errors and BMA weights; `ggplot2::geom_col()`, which is used to create the plots.
+#' @seealso [apm_pre()] to to compute the difference in average prediction errors and BMA weights; [ggplot2::geom_col()], which is used to create the plots.
 #' 
 #' @examples 
 #' data("ptpdata")

R/utils.R

@@ -2,10 +2,10 @@
 .wtd_mean <- function(x, w = NULL, subset = NULL) {
   if (!is.null(subset)) {
     if (is.null(w)) {
-      return(.wtd_mean(x[subset]))
+      return(Recall(x[subset]))
     }
 
-    return(.wtd_mean(x[subset], w = w[subset]))
+    return(Recall(x[subset], w = w[subset]))
   }
 
   if (is.null(w)) {
@@ -19,10 +19,10 @@
 .wtd_sd <- function(x, w = NULL, subset = NULL) {
   if (!is.null(subset)) {
     if (is.null(w)) {
-      return(.wtd_sd(x[subset]))
+      return(Recall(x[subset]))
     }
 
-    return(.wtd_sd(x[subset], w = w[subset]))
+    return(Recall(x[subset], w = w[subset]))
   }
 
   if (is.null(w)) {

README.Rmd

@@ -18,65 +18,26 @@ knitr::opts_chunk$set(
 <!-- badges: start -->
 <!-- badges: end -->
 
-```{r setup}
-library(apm)
-data("ptpdata")
-```
-
-### Supplying models
-
-We can specify the models to test using `apm_mod()`. This create a full cross of all supplied arguments, which include model formula, families, whether the outcome is logged or not, whether fixed effects are included or not, whether the outcome should be a difference, and whether outcome lags should appear as predictors. Below, we create a cross of 9 models.
-
-```{r}
-models <- apm_mod(deaths ~ 1,
-                  family = list("gaussian", "quasipoisson"),
-                  log = c(TRUE, FALSE),
-                  lag = 0, diff_k = 0, 
-                  time_trend = 0:2)
-
-models
-```
+## Introduction
 
-Normally, this cross would yield 12 = 3 (formulas) x 2 (families) x 2 (log T/F), but by default any models with non-linear links and `log = TRUE` are removed, leaving 9 models. If we want to manually add other models, we can so by creating a new models object and appending it to the current one.
+The `apm` package implements *Averaged Prediction Models (APM)*, a Bayesian model averaging approach for controlled pre-post designs. These designs compare differences over time between a group that becomes exposed (treated group) and one that remains unexposed (comparison group). With appropriate causal assumptions, they can identify the causal effect of the exposure/treatment.
 
-```{r}
-models2 <- apm_mod(list(deaths ~ 1),
-                   diff_k = 1)
+In APM, we specify a collection of models that predict untreated outcomes. Our causal identifying assumption is that the model's prediction errors would be equal (in expectation) in the treated and comparison groups in the absence of the exposure. This is a generalization of familiar methods like Difference-in-Differences (DiD) and Comparative Interrupted Time Series (CITS).
 
-models <- c(models, models2)
+Because many models may be plausible for this prediction task, we combine them using Bayesian model averaging. We weight each model by its robustness to violations of the causal assumption.
 
-models
-```
-
-This leaves us with 10 models.
+## Installation
 
-### Fitting the models
+To install the development version from GitHub, use:
 
-Next we fit all 10 models to the data. We do so once for each validation time to compute the average prediction error that will be used to select the optimal model. All models are fit simultaneously so the simulation can use the full joint distribution of model parameter estimates. For each validation time, each model is fit using a dataset that contains data points prior to that time.
+```{r eval = FALSE}
 
-We use `apm_fit()` to fit the models, and calculate the prediction errors and BMA weights.
+# Install devtools if not already installed
+install.packages("remotes")
 
-```{r}
-fits <- apm_pre(models,
-                data = ptpdata,
-                group_var = "group",
-                time_var = "year",
-                unit_var = "state",
-                val_times = 2004:2007)
+# Install apm package from GitHub if not already installed
+remotes::install_github("tl2624/apm")
 
-fits
 ```
 
-### Computing the ATT
-
-We compute the ATT using `apm_est()`, which uses bootstrapping to compute model uncertainty due to sampling along with uncertainty due to model selection.
-
-```{r}
-est <- apm_est(fits,
-               post_time = 2008,
-               M = 1)
-
-est
-
-summary(est)
-```
+See `vignette("apm")` for details on using the package.