This page explains the details of estimating weights using a user-defined function. The function must take in arguments that are passed to it by weightit() or weightitMSM() and return a vector of weights or a list containing the weights.

To supply user-defined function, the function object should be entered directly to method; for example, for a function fun, method = fun.

Point Treatments

The following arguments are automatically passed to the user-defined function, which should have named parameters corresponding to them:

  • treat: a vector of treatment status for each unit. This comes directly from the left hand side of the formula passed to weightit and so will have it's type (e.g., numeric, factor, etc.), which may need to be converted.

  • covs: a data frame of covariate values for each unit. This comes directly from the right hand side of the formula passed to weightit. The covariates are processed so that all columns are numeric; all factor variables are split into dummies and all interactions are evaluated. All levels of factor variables are given dummies, so the matrix of the covariates is not full rank. Users can use make_full_rank, which accepts a numeric matrix or data frame and removes columns to make it full rank, if a full rank covariate matrix is desired.

  • s.weights: a numeric vector of sampling weights, one for each unit.

  • ps: a numeric vector of propensity scores.

  • subset: a logical vector the same length as treat that is TRUE for units to be included in the estimation and FALSE otherwise. This is used to subset the input objects when exact is used. treat, covs, s.weights, and ps, if supplied, will already have been subsetted by subset.

  • estimand: a character vector of length 1 containing the desired estimand. The characters will have been converted to uppercase. If "ATC" was supplied to estimand, weightit sets focal to the control level (usually 0 or the lowest level of treat) and sets estimand to "ATT".

  • focal: a character vector of length 1 containing the focal level of the treatment when the estimand is the ATT (or the ATC as detailed above). weightit() ensures the value of focal is a level of treat.

  • stabilize: a logical vector of length 1. It is not processed by weightit() before it reaches the fitting function.

  • moments: a numeric vector of length 1. It is not processed by weightit() before it reaches the fitting function except that as.integer is applied to it. This is used in other methods to determine whether polynomials of the entered covariates are to be used in the weight estimation.

  • int: a logical vector of length 1. It is not processed by weightit() before it reaches the fitting function. This is used in other methods to determine whether interactions of the entered covariates are to be used in the weight estimation.

None of these parameters are required to be in the fitting function. These are simply those that are automatically available.

In addition, any additional arguments supplied to weightit() will be passed on to the fitting function. weightit() ensures the arguments correspond to the parameters of the fitting function and throws an error if an incorrectly named argument is supplied and the fitting function doesn't include ... as a parameter.

The fitting function must output either a numeric vector of weights or a list (or list-like object) with an entry named wither "w" or "weights". If a list, the list can contain other named entries, but only entries named "w", "weights", "ps", and "fit.obj" will be processed. "ps" is a vector of propensity scores and "fit.obj" should be an object used in the fitting process that a user may want to examine and that is included in the weightit output object as "obj" when include.obj = TRUE. The "ps" and "fit.obj" components are optional, but "weights" or "w" is required.

Longitudinal Treatments

Longitudinal treatments can be handled either by running the fitting function for point treatments for each time point and multiplying the resulting weights together or by running a method that accommodates multiple time points and outputs a single set of weights. For the former, weightitMSM() can be used with the user-defined function just as it is with weightit(). The latter method is not yet accommodated by weightitMSM(), but will be someday, maybe.

Examples


library("cobalt")
data("lalonde", package = "cobalt")

#A user-defined version of method = "ps"
my.ps <- function(treat, covs, estimand, focal = NULL) {
  covs <- make_full_rank(covs)
  d <- data.frame(treat, covs)
  f <- formula(d)
  ps <- glm(f, data = d, family = "binomial")$fitted
  w <- get_w_from_ps(ps, treat = treat, estimand = estimand,
                     focal = focal)

  return(list(w = w, ps = ps))
}

#Balancing covariates between treatment groups (binary)
(W1 <- weightit(treat ~ age + educ + married +
                  nodegree + re74, data = lalonde,
                method = my.ps, estimand = "ATT"))
#> A weightit object
#>  - method: "my.ps" (a user-defined method)
#>  - number of obs.: 614
#>  - sampling weights: none
#>  - treatment: 2-category
#>  - estimand: ATT (focal: 1)
#>  - covariates: age, educ, married, nodegree, re74
summary(W1)
#>                  Summary of weights
#> 
#> - Weight ranges:
#> 
#>            Min                                  Max
#> treated 1.0000               ||              1.0000
#> control 0.0222 |---------------------------| 2.0438
#> 
#> - Units with 5 most extreme weights by group:
#>                                            
#>               7      5      4      3      1
#>  treated      1      1      1      1      1
#>             411    595    269    409    296
#>  control 1.3303 1.4365 1.5005 1.6369 2.0438
#> 
#> - Weight statistics:
#> 
#>         Coef of Var   MAD Entropy # Zeros
#> treated       0.000 0.000   -0.00       0
#> control       0.823 0.701    0.33       0
#> 
#> - Effective Sample Sizes:
#> 
#>            Control Treated
#> Unweighted  429.       185
#> Weighted    255.99     185
bal.tab(W1)
#> Call
#>  weightit(formula = treat ~ age + educ + married + nodegree + 
#>     re74, data = lalonde, method = my.ps, estimand = "ATT")
#> 
#> Balance Measures
#>                Type Diff.Adj
#> prop.score Distance   0.0199
#> age         Contin.   0.0459
#> educ        Contin.  -0.0360
#> married      Binary   0.0044
#> nodegree     Binary   0.0080
#> re74        Contin.  -0.0275
#> 
#> Effective sample sizes
#>            Control Treated
#> Unadjusted  429.       185
#> Adjusted    255.99     185
#Balancing covariates for longitudinal treatments
# my.ps is used at each time point.
library("twang")
#> To reproduce results from prior versions of the twang package, please see the version="legacy" option described in the documentation.
data("iptwExWide", package = "twang")
(W2 <- weightitMSM(list(tx1 ~ age + gender + use0,
                        tx2 ~ tx1 + use1 + age + gender + use0,
                        tx3 ~ tx2 + use2 + tx1 + use1 + age + gender + use0),
                   data = iptwExWide,
                   method = my.ps))
#> A weightitMSM object
#>  - method: "my.ps" (a user-defined method)
#>  - number of obs.: 1000
#>  - sampling weights: none
#>  - number of time points: 3 (tx1, tx2, tx3)
#>  - treatment: 
#>     + time 1: 2-category
#>     + time 2: 2-category
#>     + time 3: 2-category
#>  - covariates: 
#>     + baseline: age, gender, use0
#>     + after time 1: tx1, use1, age, gender, use0
#>     + after time 2: tx2, use2, tx1, use1, age, gender, use0
summary(W2)
#>                  Summary of weights
#> 
#>                        Time 1                       
#>                  Summary of weights
#> 
#> - Weight ranges:
#> 
#>            Min                                   Max
#> treated 1.6667 |---|                         15.5292
#> control 2.6123 |---------------------------| 82.2349
#> 
#> - Units with 5 most extreme weights by group:
#>                                                 
#>              477     442     307     409     906
#>  treated 13.4464 13.6233 13.7059  14.436 15.5292
#>              206     282     641     547     980
#>  control 44.8027 46.1355 58.4164 79.8919 82.2349
#> 
#> - Weight statistics:
#> 
#>         Coef of Var   MAD Entropy # Zeros
#> treated       0.481 0.376   0.117       0
#> control       0.768 0.494   0.222       0
#> 
#> - Effective Sample Sizes:
#> 
#>            Control Treated
#> Unweighted  294.     706. 
#> Weighted    185.18   573.6
#> 
#>                        Time 2                       
#>                  Summary of weights
#> 
#> - Weight ranges:
#> 
#>            Min                                   Max
#> treated 1.6667 |---------------------------| 82.2349
#> control 2.6123 |--------------------------|  79.8919
#> 
#> - Units with 5 most extreme weights by group:
#>                                                 
#>               34     561     553     641     980
#>  treated  33.371 34.7405 42.9442 58.4164 82.2349
#>              109      95     206     282     547
#>  control 39.7862 42.2256 44.8027 46.1355 79.8919
#> 
#> - Weight statistics:
#> 
#>         Coef of Var   MAD Entropy # Zeros
#> treated       0.960 0.653   0.336       0
#> control       0.737 0.384   0.159       0
#> 
#> - Effective Sample Sizes:
#> 
#>            Control Treated
#> Unweighted   492.   508.  
#> Weighted     318.9  264.49
#> 
#>                        Time 3                       
#>                  Summary of weights
#> 
#> - Weight ranges:
#> 
#>            Min                                   Max
#> treated 1.6667 |-------------|               44.8027
#> control 2.6123 |---------------------------| 82.2349
#> 
#> - Units with 5 most extreme weights by group:
#>                                                 
#>              936     650     763     109     206
#>  treated  26.003  26.471 28.0582 39.7862 44.8027
#>              553     282     641     547     980
#>  control 42.9442 46.1355 58.4164 79.8919 82.2349
#> 
#> - Weight statistics:
#> 
#>         Coef of Var   MAD Entropy # Zeros
#> treated       0.773 0.565   0.245       0
#> control       0.873 0.469   0.227       0
#> 
#> - Effective Sample Sizes:
#> 
#>            Control Treated
#> Unweighted  415.     585. 
#> Weighted    235.67   366.4
#> 
bal.tab(W2)
#> Call
#>  weightitMSM(formula.list = list(tx1 ~ age + gender + use0, tx2 ~ 
#>     tx1 + use1 + age + gender + use0, tx3 ~ tx2 + use2 + tx1 + 
#>     use1 + age + gender + use0), data = iptwExWide, method = my.ps)
#> 
#> Balance summary across all time points
#>              Times     Type Max.Diff.Adj
#> prop.score 1, 2, 3 Distance       0.0251
#> age        1, 2, 3  Contin.       0.0703
#> gender     1, 2, 3   Binary       0.0263
#> use0       1, 2, 3  Contin.       0.0558
#> tx1           2, 3   Binary       0.0171
#> use1          2, 3  Contin.       0.0316
#> tx2              3   Binary       0.0085
#> use2             3  Contin.       0.0315
#> 
#> Effective sample sizes
#>  - Time 1
#>            Control Treated
#> Unadjusted  294.     706. 
#> Adjusted    185.18   573.6
#>  - Time 2
#>            Control Treated
#> Unadjusted   492.   508.  
#> Adjusted     318.9  264.49
#>  - Time 3
#>            Control Treated
#> Unadjusted  415.     585. 
#> Adjusted    235.67   366.4
# Kernel balancing using the KBAL package, available
# using devtools::install_github("chadhazlett/KBAL").
# Only the ATT and ATC are available. Use 'kbal.method'
# instead of 'method' in weightit() to choose between
# "ebal" and "el".

if (FALSE) {
kbal.fun <- function(treat, covs, estimand, focal, ...) {
    args <- list(...)
    if (is_not_null(focal))
        treat <- as.numeric(treat == focal)
    else if (estimand != "ATT")
        stop("estimand must be 'ATT' or 'ATC'.", call. = FALSE)
    if ("kbal.method" %in% names(args)) {
        names(args)[names(args) == "kbal.method"] <- "method"
    }
    args[!names(args) %in% setdiff(names(formals(KBAL::kbal)),
        c("X", "D"))] <- NULL
    k.out <- do.call(KBAL::kbal, c(list(X = covs, D = treat),
        args))
    w <- k.out$w
    return(list(w = w))
}

(Wk <- weightit(treat ~ age + educ + married +
                nodegree + re74, data = lalonde,
                method = kbal.fun, estimand = "ATT",
                kbal.method = "ebal"))
summary(Wk)
bal.tab(Wk, disp.ks = TRUE)
}