Covariate distributions

If cov_dist = "joint" then the joint density is used to estimate \(p(\mathbf{x}|Z=1)\) and \(p(\mathbf{x})\). That is, the number of cases is found for every combination of values of the covariates given in cov_cols for both trial arms. This is the most accurate approach, but may perform poorly if the number of cases is small relative to the number of covariates.

If cov_dist = "marginal" then \(p(\mathbf{x}|Z=1)\) and \(p(\mathbf{x})\) are estimated using the marginal approximation

\[\begin{align*} p\left(\mathbf{x}\mid Z=1\right) & = \prod\limits_{j=1}^k p\left(x^{(j)}\mid Z=1\right) \\ p\left(\mathbf{x}\right) & = \prod\limits_{j=1}^k p\left(x^{(j)}\right). \end{align*}\]

This is a strong assumption, since there will likely be some dependence between covariates, but it is more robust to empty combinations.

Finally, one can specify sub-groups of covariates to be treated jointly by entering a list to cov_dist. For example, if we set

cov_cols = c("category", "risk_ass", "sus_age_bin", "nprev_sus_bin")

and

cov_dist = list(c("category", "risk_ass"))

then proportions will be found for all combinations of category and risk_ass (treating these two covariates jointly) for each arm. Proportions will also be found for each value of sus_age_bin and nprev_sus_bin. These proportions will be multiplied together to estimate \(p(\mathbf{x}|Z=1)\) and \(p(\mathbf{x})\). This approach gives a middle ground between the fully marginal and fully joint approaches.

Specifying \(p(Z=1)\)

Unlike \(p(\mathbf{x}|Z=1)\) and \(p(\mathbf{x})\), \(p(Z=1)\) cannot be estimated from the data in df, but it does need to be specified to estimate the propensity scores. The function propscores() provides three approaches, specified via pz1_fn and its argument n_pz1. These are described in more detail in [the paper] but in summary

• pz1_fn = pz1_rescale(n_pz1) will set \(p(Z=1)\) such that the largest propensity score (unless any are infinite) is n_pz1
• pz1_fn = pz1_estimate(n_pz1) allows the user to specify an estimate of \(p(Z=1)\) directly, ie. \(p(Z=1) =\)n_pz1
• pz1_fn = pz1_expectedN(n_pz1) will set \(p(Z=1)\) such that the expected size of the matched comparison group will be n_pz1

Examples

We will first estimate propensity scores using the fully joint approach and an estimate of \(p(Z=1) = 0.1\):

ps_joint_est025 = propscores(
  df = eg_data,
  cov_cols = c("sus_age_bin", "category", "risk_ass"),
  arm_col = "Arm",
  intervention_level = "Intervention",
  cov_dist = "joint",
  pz1_fn = pz1_estimate(0.25)
)
#> There are 1 NA values in the sus_age_bin column. These rows will be lost.
#> There are 1 NA values in the category column. These rows will be lost.

Notice that this warns us that there are NA values in two of the columns, leading to rows being lost.

Each of these has columns for every element of cov_cols, and columns px_int, px_comp and PropScore. They will also have other columns, depending on which methods were used.

Below we display the data frame ps_joint_est025. Notice that in one row (HIGH, Theft_CD, 40+) the propensity score is Inf. This is because there are some such cases in the intervention group, but none in the comparison group, and so although the match sampling process requires some comparison cases, there are none to be sampled. There are also some rows (MEDIUM, Theft / CD, Under 29 and MEDIUM, Theft / CD, 33-39) for which there are no cases in the intervention cohort, and so no comparison cases will be sampled (and the propensity score is zero).

ps_joint_est025 |>
  dplyr::select(c("risk_ass", "category", "sus_age_bin", "px_int", "px_comp", "count_comp", "PropScore")) |>
#  dplyr::select(c("risk_ass", "category", "sus_age_bin", "count_int", "count_comp", "PropScore")) |>
  knitr::kable(digits = 3)

To avoid the arbitrary results of small numbers, but preserve a potential link between category and risk_ass, we might choose the partially marginal approach instead (for demonstration purposes we also switch to specifying \(p(z=1)\) so that the expected size of the comparison group is 250),

ps_mix_exp250 = propscores(
  df = eg_data,
  cov_cols = c("sus_age_bin", "category", "risk_ass"),
  arm_col = "Arm",
  intervention_level = "Intervention",
  cov_dist = list(c("risk_ass", "category")),
  pz1_fn = pz1_expectedN(250)
)
#> There are 1 NA values in the sus_age_bin column. These rows will be lost.
#> There are 1 NA values in the category column. These rows will be lost.

This has no infinite propensity scores

ps_mix_exp250 |>
  dplyr::select(c("risk_ass", "category", "sus_age_bin", "px_int", "px_comp", "PropScore")) |>
  #  dplyr::select(c("risk_ass", "category", "sus_age_bin", "count_int", "count_comp", "PropScore")) |>
  knitr::kable(digits = 3)

An alternative would be to adopt the fully marginal approach (changing the method for \(p(Z=1)\) again for demonstration purposes, now fixing \(p(Z=1)\) so that the largest propensity score is 2)

ps_marg_resc2 = propscores(
  df = eg_data,
  cov_cols = c("risk_ass", "category", "sus_age_bin", "nprev_sus_bin"),
  arm_col = "Arm",
  intervention_level = "Intervention",
  cov_dist = "marginal",
  pz1_fn = pz1_rescale(2)
)
#> There are 1 NA values in the category column. These rows will be lost.
#> There are 1 NA values in the sus_age_bin column. These rows will be lost.

Examples

candidates(
  df = eg_data,
  propscores = ps_joint_est025,
  seeds = 1:10,
  smd_cols = c("risk_ass", "category", "sus_age_bin"),
  smd_order = 2
)
#>    seed n_int n_comp SMD_risk_ass SMD_category SMD_sus_age_bin
#> 1     1   314    168   0.23308149   0.02636891      0.11576866
#> 2     2   314    171   0.04256640   0.05493603      0.07528036
#> 3     3   314    170   0.03389044   0.06374867      0.03572941
#> 4     4   314    185   0.02598238   0.04409182      0.09550639
#> 5     5   314    178   0.08440418   0.11249425      0.06965982
#> 6     6   314    172   0.03861228   0.12868680      0.07227845
#> 7     7   314    162   0.09666017   0.09146735      0.07851967
#> 8     8   314    172   0.02605045   0.02221172      0.09039220
#> 9     9   314    166   0.01098257   0.09910340      0.11214305
#> 10   10   314    162   0.05372981   0.15983827      0.06603256
#>    SMD_category.risk_ass SMD_risk_ass.sus_age_bin SMD_category.sus_age_bin
#> 1             0.28132402                0.2832541                0.1828508
#> 2             0.10018936                0.1447606                0.2225351
#> 3             0.15439603                0.1408354                0.1618457
#> 4             0.11952593                0.1372929                0.1725952
#> 5             0.19985826                0.1698178                0.2036004
#> 6             0.16773534                0.1146762                0.1860480
#> 7             0.19176077                0.1342730                0.2682980
#> 8             0.05807104                0.1154419                0.1484787
#> 9             0.14409534                0.1461870                0.2080970
#> 10            0.16829690                0.2032046                0.2425481