| Title: | Propensity Score Sampling for quasi-experimental studies with an external comparison group |
|---|---|
| Description: | IN DEVELOPMENT! Propensity score sampling is a method for finding a matched comparison group in a quasi-experimental study where an 'external' comparison cohort is being used. By 'external', we mean a cohort for which the intervention is not available to anyone (for example, a different time or region from the treatment group). |
| Authors: | Rachel Oughton [aut, cre] |
| Maintainer: | Rachel Oughton <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.0.0.9000 |
| Built: | 2026-10-19 13:54:18 UTC |
| Source: | https://github.com/RachelOughton/pss |
Create a summary table for a vector of seeds
candidates( df, propscores, seeds, smd_cols, smd_order, replace = TRUE, downsample = FALSE, drop_int = FALSE, use_attr = TRUE, cov_cols = NULL, arm_col = NULL, intervention_level = NULL )candidates( df, propscores, seeds, smd_cols, smd_order, replace = TRUE, downsample = FALSE, drop_int = FALSE, use_attr = TRUE, cov_cols = NULL, arm_col = NULL, intervention_level = NULL )
df |
A data frame from which to generate matched samples |
propscores |
A data frame (created by |
seeds |
A vector of random seeds to be used |
smd_cols |
A character vector of columns for which to assess balance. This can be the same as |
smd_order |
An integer denoting the highest order of interaction for which to calculate the standardised mean difference (SMD) |
replace |
logical. Passed to |
downsample |
logical. Passed to |
drop_int |
logical. Passed to |
use_attr |
logical. Should the attributes of |
cov_cols |
(only specify if |
arm_col |
(only specify if |
intervention_level |
(only specify if |
candidates() returns a data frame with a row for each random seed and columns
seed
n_int: the number of treated / intervention cases included in the matched sample for that seed
n_comp: the number of comparison cases included in the matched sample for that seed
SMD_...: column for each covariate and covariate interaction, giving the calculated SMD, found using smd::smd()
matched_sample(), expand_matched_df()
## First use `propscores` to create a propensity score data frame ## This example uses the joint approach and specifies p(Z=1) such ## that the expected size of the comparison group is 300 eg_j_exp300 <- propscores( df = eg_data, cov_cols = c("risk_ass", "category", "sus_age_bin"), arm_col = "Arm", intervention_level = "Intervention", cov_dist = "joint", pz1_fn = pz1_expectedN(300) ) ## Now we can use `eg_j_exp300` to find a matched comparison group ## for each seed in the vector `seeds`. [candidates()] calls ## [matched_sample()] each time it does this, and we specify that ## it should downsample the intervention cases (`downsample = TRUE`) ## rather than sampling from the comparison cases with replacement. ## In the output of [candidates()] we include SMD values for ## four variables and include interactions up to order 3. Because ## `suspect_age` is numeric it will not be included in the interactions. candidates( df = eg_data, propscores = eg_j_exp300, seeds = 1:20, smd_cols = c("suspect_age", "category", "risk_ass", "sus_age_bin"), smd_order = 3, downsample = TRUE, replace = FALSE )## First use `propscores` to create a propensity score data frame ## This example uses the joint approach and specifies p(Z=1) such ## that the expected size of the comparison group is 300 eg_j_exp300 <- propscores( df = eg_data, cov_cols = c("risk_ass", "category", "sus_age_bin"), arm_col = "Arm", intervention_level = "Intervention", cov_dist = "joint", pz1_fn = pz1_expectedN(300) ) ## Now we can use `eg_j_exp300` to find a matched comparison group ## for each seed in the vector `seeds`. [candidates()] calls ## [matched_sample()] each time it does this, and we specify that ## it should downsample the intervention cases (`downsample = TRUE`) ## rather than sampling from the comparison cases with replacement. ## In the output of [candidates()] we include SMD values for ## four variables and include interactions up to order 3. Because ## `suspect_age` is numeric it will not be included in the interactions. candidates( df = eg_data, propscores = eg_j_exp300, seeds = 1:20, smd_cols = c("suspect_age", "category", "risk_ass", "sus_age_bin"), smd_order = 3, downsample = TRUE, replace = FALSE )
Display the treatment cases that will be dropped under a particular matched sample
dropped_df(matched_df, use_attr = TRUE, cov_cols = NULL)dropped_df(matched_df, use_attr = TRUE, cov_cols = NULL)
matched_df |
A data frame output by |
use_attr |
logical. Should the attributes of |
cov_cols |
(only specify if |
dropped_df() returns a data frame with the following columns
A column for each of cov_cols
n_dropped: how many treated / intervention cases with this combination of cov_cols have been dropped / downsampled
n_kept: how many treated / intervention cases with this combination of cov_cols have been kept
The rows will be ordered in descending order of n_dropped
## First use [propscores()] to create a propensity score data frame ## This example uses the joint approach and specifies p(Z=1) such ## that the expected size of the comparison group is 300 eg_j_exp300 <- propscores( df = eg_data, cov_cols = c("risk_ass", "category", "sus_age_bin"), arm_col = "Arm", intervention_level = "Intervention", cov_dist = "joint", pz1_fn = pz1_expectedN(300) ) ## Now we can use `eg_j_exp300` to find a matched comparison group ## In this version we sample from the comparison group with replacement ## and do not drop any intervention cases for which there are no ## comparison cases to sample match_j_exp300_replace <- matched_sample( df = eg_data, propscores = eg_j_exp300, seed = 20, replace = TRUE, downsample = FALSE, drop_int = FALSE ) ## Finally we use [dropped_df()] to summarise the treated / intervention ## cases that have been dropped or downsampled ## included `include` times match_j_exp300_replace |> dropped_df()## First use [propscores()] to create a propensity score data frame ## This example uses the joint approach and specifies p(Z=1) such ## that the expected size of the comparison group is 300 eg_j_exp300 <- propscores( df = eg_data, cov_cols = c("risk_ass", "category", "sus_age_bin"), arm_col = "Arm", intervention_level = "Intervention", cov_dist = "joint", pz1_fn = pz1_expectedN(300) ) ## Now we can use `eg_j_exp300` to find a matched comparison group ## In this version we sample from the comparison group with replacement ## and do not drop any intervention cases for which there are no ## comparison cases to sample match_j_exp300_replace <- matched_sample( df = eg_data, propscores = eg_j_exp300, seed = 20, replace = TRUE, downsample = FALSE, drop_int = FALSE ) ## Finally we use [dropped_df()] to summarise the treated / intervention ## cases that have been dropped or downsampled ## included `include` times match_j_exp300_replace |> dropped_df()
An anonymised and slightly altered version of data from a study investigating forensic marking to protect victims of domestic abuse in the UK
eg_dataeg_data
eg_dataA data frame with 999 rows and 20 columns:
Whether the victim was given forensic marking equipment
Whether the perpetrator was warned about the use of forensic marking equipment (can only be Y if FM_equip is TRUE)
The principal crime committed at the reference incident, according to the home office lookup table
...
Expand a matched comparison data frame
expand_matched_df(matched_df)expand_matched_df(matched_df)
matched_df |
A data frame output by |
expand_matched_df() returns a data frame with all the columns of matched_df except include. Each row (row i) is now included include[i] times. That is, rows with matched_df$include[i]==0 are omitted from the returned data frame, and rows with matched_df$include[i] > 1 are included multiple times.
## First use `propscores` to create a propensity score data frame ## This example uses the joint approach and specifies p(Z=1) such ## that the expected size of the comparison group is 300 eg_j_exp300 <- propscores( df = eg_data, cov_cols = c("risk_ass", "category", "sus_age_bin"), arm_col = "Arm", intervention_level = "Intervention", cov_dist = "joint", pz1_fn = pz1_expectedN(300) ) ## Now we can use `eg_j_exp300` to find a matched comparison group ## In this version we sample from the comparison group with replacement ## and do not drop any intervention cases for which there are no ## comparison cases to sample match_j_exp300_replace <- matched_sample( df = eg_data, propscores = eg_j_exp300, seed = 20, replace = TRUE, downsample = FALSE, drop_int = FALSE ) ## Finally we use `expand_matched_df` to expand this so that each row is ## included `include` times match_j_exp300_replace |> expand_matched_df()## First use `propscores` to create a propensity score data frame ## This example uses the joint approach and specifies p(Z=1) such ## that the expected size of the comparison group is 300 eg_j_exp300 <- propscores( df = eg_data, cov_cols = c("risk_ass", "category", "sus_age_bin"), arm_col = "Arm", intervention_level = "Intervention", cov_dist = "joint", pz1_fn = pz1_expectedN(300) ) ## Now we can use `eg_j_exp300` to find a matched comparison group ## In this version we sample from the comparison group with replacement ## and do not drop any intervention cases for which there are no ## comparison cases to sample match_j_exp300_replace <- matched_sample( df = eg_data, propscores = eg_j_exp300, seed = 20, replace = TRUE, downsample = FALSE, drop_int = FALSE ) ## Finally we use `expand_matched_df` to expand this so that each row is ## included `include` times match_j_exp300_replace |> expand_matched_df()
matched_sample() uses propensity score sampling to generate a matched sample.
matched_sample( df, propscores, seed, replace = TRUE, downsample = FALSE, drop_int = FALSE, use_attr = TRUE, cov_cols = NULL, arm_col = NULL, intervention_level = NULL )matched_sample( df, propscores, seed, replace = TRUE, downsample = FALSE, drop_int = FALSE, use_attr = TRUE, cov_cols = NULL, arm_col = NULL, intervention_level = NULL )
df |
A data frame containing treatment group and comparison cohort data. These should |
propscores |
A data frame with a column for each variable in |
seed |
The random seed to be used, for reproducibility |
replace |
logical. Should the comparison cases be sampled with replacement? |
downsample |
logical. Should the treatment cases be downsampled when there aren't enough comparison cases? Only one of |
drop_int |
logical. Should cases in the intervention/treated group be dropped if there are no equivalent cases in the comparison cohort to sample from? |
use_attr |
logical. Should the attributes of |
cov_cols |
(only specify if |
arm_col |
(only specify if |
intervention_level |
(only specify if |
matched_sample() returns a data frame that is an expanded form of df, with the same number of rows, in the same order, but with some extra columns added:
PropScore: the propensity score for that row, from propscores
include: how many copies of this row are included in the matched dataset. Zero means the row has been dropped. Values greater than one mean the row will be duplicated.
seed: the random seed that was used. This will be the same for all rows and is included in the output for reproducibility.
The data frame can be expanded into a usable data frame (that is, with row i included include[i] times) using expand_matched_df().
matched_sample() will also issue messages indicating:
How many comparison cases have been repeated (and how many times) if replace == TRUE
How many intervention cases have been downsampled, if downsample == TRUE
How many (if any) intervention cases have been dropped because there are no similar comparison cases to sample from, if drop_int == TRUE
How many (if any) intervention cases have been kept even though there are no similar comparison cases to sample from, if drop_int == TRUE
## First use `propscores` to create a propensity score data frame ## This example uses the joint approach and specifies p(Z=1) such ## that the expected size of the comparison group is 300 eg_j_exp300 <- propscores( df = eg_data, cov_cols = c("risk_ass", "category", "sus_age_bin"), arm_col = "Arm", intervention_level = "Intervention", cov_dist = "joint", pz1_fn = pz1_expectedN(300) ) ## Now we can use `eg_j_exp300` to find a matched comparison group ## In this example we use the attribute information in `propscores`, ## we sample from the comparison group with replacement ## and do not drop any intervention cases for which there are no ## comparison cases to sample match_j_exp300_replace <- matched_sample( df = eg_data, propscores = eg_j_exp300, seed = 20, replace = TRUE, downsample = FALSE, drop_int = FALSE )## First use `propscores` to create a propensity score data frame ## This example uses the joint approach and specifies p(Z=1) such ## that the expected size of the comparison group is 300 eg_j_exp300 <- propscores( df = eg_data, cov_cols = c("risk_ass", "category", "sus_age_bin"), arm_col = "Arm", intervention_level = "Intervention", cov_dist = "joint", pz1_fn = pz1_expectedN(300) ) ## Now we can use `eg_j_exp300` to find a matched comparison group ## In this example we use the attribute information in `propscores`, ## we sample from the comparison group with replacement ## and do not drop any intervention cases for which there are no ## comparison cases to sample match_j_exp300_replace <- matched_sample( df = eg_data, propscores = eg_j_exp300, seed = 20, replace = TRUE, downsample = FALSE, drop_int = FALSE )
candidates() data frame. If a threshold value for the SMD is given (via smd_thresh) then the output for those seeds with any SMD value exceeding the threshold are greyed, and a message is generated listing the successful seeds.Plots standardised mean difference (SMD) for each covariate and interaction in a candidates() data frame. If a threshold value for the SMD is given (via smd_thresh) then the output for those seeds with any SMD value exceeding the threshold are greyed, and a message is generated listing the successful seeds.
plot_smd(summary_df, smd_thresh = NULL)plot_smd(summary_df, smd_thresh = NULL)
summary_df |
A data frame generated by |
smd_thresh |
Numerical. The threshold above which the SMD is considered 'high'. If |
A plot with covariates and covariate interactions on the x axis and SMD on the y axis. Lines connect SMD values for the same seed.
## First use `propscores` to create a propensity score data frame ## This example uses the joint approach and specifies p(Z=1) such ## that the expected size of the comparison group is 300 eg_j_exp300 <- propscores( df = eg_data, cov_cols = c("risk_ass", "category", "sus_age_bin"), arm_col = "Arm", intervention_level = "Intervention", cov_dist = "joint", pz1_fn = pz1_expectedN(300) ) ## Now we can use [candidates()] to tabulate the SMD values for ## four variables and include interactions up to order 3. Because ## `suspect_age` is numeric it will not be included in the interactions. sdf20 = candidates( df = eg_data, propscores = eg_j_exp300, seeds = 1:20, smd_cols = c("suspect_age", "category", "risk_ass", "sus_age_bin"), smd_order = 3, downsample = TRUE, replace = FALSE ) ## Finally we use [plot_smd()] to display the SMD values for each ## covariate and interaction. If a numeric value is given for `smd_thresh` ## then points and lines for seeds for which at lease one SMD values ## is above the threshold are plotted with transparency `alpha = 0.5`. plot_smd(sdf20, smd_thresh = 0.15)## First use `propscores` to create a propensity score data frame ## This example uses the joint approach and specifies p(Z=1) such ## that the expected size of the comparison group is 300 eg_j_exp300 <- propscores( df = eg_data, cov_cols = c("risk_ass", "category", "sus_age_bin"), arm_col = "Arm", intervention_level = "Intervention", cov_dist = "joint", pz1_fn = pz1_expectedN(300) ) ## Now we can use [candidates()] to tabulate the SMD values for ## four variables and include interactions up to order 3. Because ## `suspect_age` is numeric it will not be included in the interactions. sdf20 = candidates( df = eg_data, propscores = eg_j_exp300, seeds = 1:20, smd_cols = c("suspect_age", "category", "risk_ass", "sus_age_bin"), smd_order = 3, downsample = TRUE, replace = FALSE ) ## Finally we use [plot_smd()] to display the SMD values for each ## covariate and interaction. If a numeric value is given for `smd_thresh` ## then points and lines for seeds for which at lease one SMD values ## is above the threshold are plotted with transparency `alpha = 0.5`. plot_smd(sdf20, smd_thresh = 0.15)
propscores() generates a data frame of propensity score estimates from study / trial data, to be used in the propensity score sampling method, using the Bayes Theorem factorisation of the propensity score
Of the functions listed, only propscores() should be used directly.
.props_fun() is called by propscores()
to find the proportions of the data within each category for each arm. One of pz1_estimate(), pz1_rescale() or pz1_expectedN() should be given (with its own argument n_pz1) as the pz1_fn argument, to determine which method is used to specify .
propscores(df, cov_cols, arm_col, intervention_level, cov_dist, pz1_fn) .props_fun(df, cov_cols, cov_list) pz1_estimate(n_pz1) pz1_rescale(n_pz1) pz1_expectedN(n_pz1)propscores(df, cov_cols, arm_col, intervention_level, cov_dist, pz1_fn) .props_fun(df, cov_cols, cov_list) pz1_estimate(n_pz1) pz1_rescale(n_pz1) pz1_expectedN(n_pz1)
df |
A data frame containing treatment group and comparison cohort data |
cov_cols |
A character vector containing the column names of the covariates to be matched on. These should all be factor / categorical data. |
arm_col |
The name of the column indicating which rows are treated cases and which are comparison cases. These should be factor or character, with only two levels / options. |
intervention_level |
The value in the |
cov_dist |
Either "joint" for fully joint distribution, "marginal" for fully marginal distribution, or a list of of character vectors indicating groups of covariates that should be considered jointly. |
pz1_fn |
One of the functions used to specify p(Z=1): |
cov_list |
An argument to |
n_pz1 |
An argument to |
propscores() returns a data frame with a row for each combination of values of cov_cols and the following columns
a column for each element of cov_cols
px_int, the estimate of
px_comp, the estimate of
count_comp, the number of cases in the comparison cohort
count_int, the number of cases in the intervention cohort
PropScore, which contains the estimated propensity score
Some further columns, depending on which method is used
The data frame is assigned attributes cov_cols, arm_col and intervention_level, storing the information for use in subsequent functions such as matched_sample() (if use_attr is set to TRUE).
.props_fun returns a function that creates a data frame containing px_int ( ), px_comp and ratio (px_int / px_comp) for every combination of the levels of the covariates in cov_cols, using the approach specified via cov_dist.
pz1_estimate returns the estimate of using the 'estimate' method (eg. a value provided by an expert). The argument is used as the estimate.
pz1_rescale finds the value of such that the largest propensity score will be n_pz1 (the argument given to pz1_fn)
pz1_expectedN finds the value of such that the expected size of the matched comparison group will be n_pz1 (the argument to pz1_fn).
.props_fun(): Find table of estimates of , and their ratio using fully joint approach
pz1_estimate(): Specify using the estimate method. The argument given is the estimate used
pz1_rescale(): Specify using the rescale method. The argument given is the value the largest estimated propensity score should take
pz1_expectedN(): Specify using the expectedN method. The argument given should be the expected size of the matched comparison group
# Creates a propensity score data frame from `eg_data` # using the joint approach, and rescaling so that the # highest propensity score is 1. propscores( df = eg_data, cov_cols = c("risk_ass", "category", "sus_age_bin"), arm_col = "Arm", intervention_level = "Intervention", cov_dist = "joint", pz1_fn = pz1_rescale(1) ) # Creates a propensity score data frame from `eg_data` # using the fully marginal approach, and setting p(Z=1) # so that the expected size of the matched comparison # group will be 300 propscores( df = eg_data, cov_cols = c("risk_ass", "category", "sus_age_bin"), arm_col = "Arm", intervention_level = "Intervention", cov_dist = "marginal", pz1_fn = pz1_expectedN(300) ) # Creates a propensity score data frame from `eg_data` # for which "risk_ass" and "category" are considered # jointly and "sus_age_bin" marginally, and setting p(Z=1) # so that the expected size of the matched comparison # group will be 300 propscores( df = eg_data, cov_cols = c("risk_ass", "category", "sus_age_bin"), arm_col = "Arm", intervention_level = "Intervention", cov_dist = "marginal", pz1_fn = pz1_expectedN(300) )# Creates a propensity score data frame from `eg_data` # using the joint approach, and rescaling so that the # highest propensity score is 1. propscores( df = eg_data, cov_cols = c("risk_ass", "category", "sus_age_bin"), arm_col = "Arm", intervention_level = "Intervention", cov_dist = "joint", pz1_fn = pz1_rescale(1) ) # Creates a propensity score data frame from `eg_data` # using the fully marginal approach, and setting p(Z=1) # so that the expected size of the matched comparison # group will be 300 propscores( df = eg_data, cov_cols = c("risk_ass", "category", "sus_age_bin"), arm_col = "Arm", intervention_level = "Intervention", cov_dist = "marginal", pz1_fn = pz1_expectedN(300) ) # Creates a propensity score data frame from `eg_data` # for which "risk_ass" and "category" are considered # jointly and "sus_age_bin" marginally, and setting p(Z=1) # so that the expected size of the matched comparison # group will be 300 propscores( df = eg_data, cov_cols = c("risk_ass", "category", "sus_age_bin"), arm_col = "Arm", intervention_level = "Intervention", cov_dist = "marginal", pz1_fn = pz1_expectedN(300) )
pss_check_data checks that the data frame provided is suitable. Will return either a very slightly cleaned up version, or an error message.
The function does the following:
checks the arm_col column has two unique values, one of which is intervention_level
adds an Arm column, with two levels ("Intervention" and "Comparison")
checks the columns named in cov_cols are present and in a suitable form (ideally factor, but if character it will coerce to factor and give a warning)
Removes rows for which any of cov_cols or arm_col is missing.
pss_check_data(df, cov_cols, arm_col, intervention_level)pss_check_data(df, cov_cols, arm_col, intervention_level)
df |
A data frame containing treatment group and comparison cohort data |
cov_cols |
A character vector containing the column names of the covariates to be matched on. These should all be factor / categorical data. |
arm_col |
The name of the column indicating which rows are treated cases and which are comparison cases. These should be factor or character, with only two levels / options. |
intervention_level |
The value in the |
pss_check_data returns a (possibly slightly tidied up) data frame or an error message. The data frame will have an Arm column with levels Intervention and Comparison.
pss_check_data( df = eg_data, cov_cols = c("risk_ass", "category", "sus_age_bin"), arm_col = "Arm", intervention_level = "Intervention" )pss_check_data( df = eg_data, cov_cols = c("risk_ass", "category", "sus_age_bin"), arm_col = "Arm", intervention_level = "Intervention" )