Trimming and truncation improvements

NAMESPACE

@@ -11,13 +11,22 @@ S3method("[",ps_trim)
 S3method("[",ps_trim_matrix)
 S3method("[",ps_trunc)
 S3method("[",ps_trunc_matrix)
+S3method("[",psw)
+S3method("is.na<-",ps_trim)
+S3method("is.na<-",ps_trunc)
+S3method("is.na<-",psw)
 S3method(Summary,ps_calib)
 S3method(Summary,ps_trim)
 S3method(Summary,ps_trunc)
 S3method(anyDuplicated,ps_trim)
 S3method(anyDuplicated,ps_trunc)
 S3method(as.data.frame,ipw_diagnostic_se)
+S3method(as.data.frame,ps_trim_matrix)
+S3method(as.data.frame,ps_trunc_matrix)
 S3method(as.double,ps_calib)
+S3method(c,ps_trim)
+S3method(c,ps_trunc)
+S3method(c,psw)
 S3method(causalgenerics::ipw,default)
 S3method(causalgenerics::ipw,glm)
 S3method(causalgenerics::ipw,joint_wt_models)
@@ -74,6 +83,10 @@ S3method(is_unit_truncated,default)
 S3method(is_unit_truncated,ps_trunc)
 S3method(is_unit_truncated,ps_trunc_matrix)
 S3method(is_unit_truncated,psw)
+S3method(is_unit_wt_truncated,default)
+S3method(is_unit_wt_truncated,psw)
+S3method(is_wt_truncated,default)
+S3method(is_wt_truncated,psw)
 S3method(max,ps_calib)
 S3method(max,ps_trim)
 S3method(max,ps_trunc)
@@ -88,6 +101,8 @@ S3method(model_fits_levels,multinom)
 S3method(model_levels,multinom)
 S3method(na.omit,ps_trim)
 S3method(obj_print_footer,psw)
+S3method(predict_binary_ps,default)
+S3method(predict_binary_ps,multinom)
 S3method(print,ipw_diagnostic_se)
 S3method(print,propensity_density)
 S3method(print,propensity_density_meta)
@@ -107,12 +122,14 @@ S3method(ps_tilt,numeric)
 S3method(ps_trim,data.frame)
 S3method(ps_trim,default)
 S3method(ps_trim,glm)
+S3method(ps_trim,lm)
 S3method(ps_trim,matrix)
 S3method(ps_trim,multinom)
 S3method(ps_trim,ps_trim)
 S3method(ps_trunc,data.frame)
 S3method(ps_trunc,default)
 S3method(ps_trunc,glm)
+S3method(ps_trunc,lm)
 S3method(ps_trunc,matrix)
 S3method(ps_trunc,multinom)
 S3method(ps_trunc,ps_trunc)
@@ -127,6 +144,8 @@ S3method(sort,ps_trim)
 S3method(sort,ps_trunc)
 S3method(summary,ps_trim)
 S3method(summary,ps_trunc)
+S3method(tibble::as_tibble,ps_trim_matrix)
+S3method(tibble::as_tibble,ps_trunc_matrix)
 S3method(unique,ps_trim)
 S3method(unique,ps_trunc)
 S3method(vec_arith,ps_calib)
@@ -291,6 +310,12 @@ S3method(wt_entropy,numeric)
 S3method(wt_entropy,ps_calib)
 S3method(wt_entropy,ps_trim)
 S3method(wt_entropy,ps_trunc)
+S3method(wt_trunc,default)
+S3method(wt_trunc,numeric)
+S3method(wt_trunc,ps_calib)
+S3method(wt_trunc,ps_trim)
+S3method(wt_trunc,ps_trunc)
+S3method(wt_trunc,psw)
 export("estimand<-")
 export(as_conditional)
 export(as_marginal)
@@ -317,6 +342,8 @@ export(is_refit)
 export(is_stabilized)
 export(is_unit_trimmed)
 export(is_unit_truncated)
+export(is_unit_wt_truncated)
+export(is_wt_truncated)
 export(joint_wt_meta)
 export(joint_wt_models)
 export(new_psw)
@@ -341,6 +368,8 @@ export(wt_atu)
 export(wt_cens)
 export(wt_entropy)
 export(wt_joint)
+export(wt_trunc)
+export(wt_trunc_sensitivity)
 import(vctrs)
 importFrom(causalgenerics,
   "estimand<-",

R/categorical-models.R

@@ -166,7 +166,7 @@ model_levels.multinom <- function(model) {
 # that leave nothing to read as the probability of the exposure, which is what a
 # binary exposure needs.
 #
-# `remedy` goes unread here. What is wrong with a `multinom` of more than two
+# `remedy` and `problem` go unread here. What is wrong with a `multinom` of more than two
 # levels is the number of levels rather than the family, so the remedy is the
 # one written below whichever argument the model arrived in. A `multinom`
 # supplied as a numerator model never reaches this method at all:
@@ -178,6 +178,7 @@ check_binary_model_family.multinom <- function(
   model,
   arg = ".propensity",
   remedy = NULL,
+  problem = NULL,
   call = rlang::caller_env()
 ) {
   n_levels <- length(model$lev)
@@ -211,6 +212,14 @@ extract_binary_ps.multinom <- function(model, call = rlang::caller_env()) {
   as.numeric(stats::fitted(model))
 }
 
+# `predict()` on a `multinom` offers no response type. Its probabilities for a
+# two-level fit are the one column `fitted()` reports, the probability of the
+# second level.
+#' @export
+predict_binary_ps.multinom <- function(model, newdata) {
+  as.numeric(stats::predict(model, newdata = newdata, type = "probs"))
+}
+
 # A `multinom` keeps no model frame, so `model.frame()` rebuilds one by
 # re-evaluating the fitting call in the environment the formula came from. That
 # environment no longer holds the arguments of a fit made inside a function, and

R/continuous-models.R

@@ -36,8 +36,15 @@ extract_continuous_ps.lm <- function(model, call = rlang::caller_env()) {
 # route says which, so the fit is refused here rather than downstream, where a
 # matrix with a row for each unit is only ever described by a length that the
 # caller never wrote.
+#
+# `problem` and `remedy` say what the caller needed the single mean for and what
+# to do instead, since trimming a dose model cannot take the means on their own.
 check_continuous_model_response <- function(
   model,
+  problem = "Weights for a continuous exposure need a model of one conditional
+             mean for each unit.",
+  remedy = "Fit the exposure on its own, or pass the conditional means of this
+            exposure to {.arg .propensity} as a numeric vector.",
   call = rlang::caller_env()
 ) {
   fitted_values <- stats::fitted(model)
@@ -50,13 +57,11 @@ check_continuous_model_response <- function(
 
   abort(
     c(
-      "Weights for a continuous exposure need a model of one conditional mean
-       for each unit.",
+      problem,
       x = "{.arg .propensity} is {.cls {class(model)[[1]]}}, a fit of
            {dims[[2]]} response{?s}, whose fitted values are {dims[[1]]} by
            {dims[[2]]}.",
-      i = "Fit the exposure on its own, or pass the conditional means of this
-           exposure to {.arg .propensity} as a numeric vector."
+      i = remedy
     ),
     error_class = "propensity_ps_shape_error",
     call = call
@@ -101,12 +106,16 @@ extract_continuous_ps.gam <- function(model, call = rlang::caller_env()) {
 # `arg` is the argument the model arrived in, and `remedy` the advice that goes
 # with it: a propensity score model's fitted means can be passed directly
 # instead, and a numerator model a caller stabilizes on is a model or nothing.
+# `problem` states what the caller needed the model for, since a refit of a
+# trimmed dose model needs the same model for a different purpose.
 check_continuous_model_family <- function(
   model,
   arg = ".propensity",
   remedy = "Fit the propensity score model with {.fun gaussian}, {.fun lm},
             {.fun mgcv::gam}, or {.fun MASS::rlm}, or pass fitted conditional
             means to {.arg .propensity} directly.",
+  problem = "Weights for a continuous exposure need a model of its conditional
+             mean with a single spread.",
   call = rlang::caller_env()
 ) {
   family <- model[["family"]]
@@ -137,8 +146,7 @@ check_continuous_model_family <- function(
 
   abort(
     c(
-      "Weights for a continuous exposure need a model of its conditional mean
-       with a single spread.",
+      problem,
       x = fit_with,
       i = remedy
     ),
@@ -154,14 +162,19 @@ check_continuous_model_family <- function(
 # own.
 #
 # `arg` names the argument the model arrived in and `remedy` says what to do
-# about it, because the same question is asked of two models: the propensity
-# score model of a binary exposure, and the numerator model that estimates the
-# probability the weights are stabilized on.
+# about it, because the same question is asked of more than one model: the
+# propensity score model of a binary exposure, the numerator model that
+# estimates the probability the weights are stabilized on, and the model a
+# modified score is read from or refit with. `problem` states what the caller
+# needed a probability for, which differs between building weights, bounding
+# scores, and refitting them.
 check_binary_model_family <- function(
   model,
   arg = ".propensity",
   remedy = "Fit the propensity score model with {.fun binomial}, or pass
             fitted probabilities to {.arg .propensity} directly.",
+  problem = "Weights for a binary exposure need a model of the probability of
+             the exposure.",
   call = rlang::caller_env()
 ) {
   UseMethod("check_binary_model_family")
@@ -178,23 +191,20 @@ check_binary_model_family.default <- function(
   arg = ".propensity",
   remedy = "Fit the propensity score model with {.fun binomial}, or pass
             fitted probabilities to {.arg .propensity} directly.",
+  problem = "Weights for a binary exposure need a model of the probability of
+             the exposure.",
   call = rlang::caller_env()
 ) {
   family <- model[["family"]]
 
   # A family object is a list. An element that is not one is not a family, and
   # `$` on an atomic vector raises an error of base R's rather than one of this
   # package's, so what the element is gets read before anything is asked of it.
-  #
-  # `isTRUE()` because a family object that carries no name answers the test
-  # with nothing at all, which `&&` can read as neither true nor false. A
-  # family like that is not one of the two, and is refused below by what it is
-  # missing.
+  # A family object that carries no name is not one of the binomial two, and is
+  # refused below by what it is missing.
   family_object <- is.list(family)
-  binomial_family <- family_object &&
-    isTRUE(family$family %in% c("binomial", "quasibinomial"))
 
-  if (binomial_family) {
+  if (is_binomial_family(family)) {
     return(invisible(NULL))
   }
 
@@ -214,8 +224,7 @@ check_binary_model_family.default <- function(
 
   abort(
     c(
-      "Weights for a binary exposure need a model of the probability of the
-       exposure.",
+      problem,
       x = fitted_by,
       i = remedy
     ),
@@ -224,6 +233,15 @@ check_binary_model_family.default <- function(
   )
 }
 
+# Whether a family fits the probability of a binary exposure. A family object is
+# a list, and an element that is not one is not a family. `isTRUE()` because a
+# family object that carries no name answers the test with nothing at all, which
+# `&&` can read as neither true nor false.
+is_binomial_family <- function(family) {
+  is.list(family) &&
+    isTRUE(family$family %in% c("binomial", "quasibinomial"))
+}
+
 # Whether a family object names itself. Every family built by `binomial()` and
 # its neighbors holds its name in `$family`, and a message about one is written
 # from that string; an object that carries a link and no name has nothing to be

R/density.R

@@ -856,7 +856,7 @@ t_sigma_mle <- function(residuals, df, call = rlang::caller_env()) {
       c(
         "The scale of a {.val t} density cannot be estimated by maximum
          likelihood from these residuals.",
-        x = "{n_zero} of the {n} residuals of the propensity score model
+        x = "{n_zero} of the {n} residuals of the model of the exposure
              {?is/are} exactly zero, and with {df} degrees of freedom that
              leaves the likelihood no maximum at a positive scale.",
         i = "Use {.code sigma_method = \"rms\"}, or supply a spread with
@@ -902,7 +902,7 @@ check_sigma_method <- function(.sigma, density, call = rlang::caller_env()) {
       "{.code sigma_method = \"mle\"} cannot be used with {.arg .sigma}.",
       x = "{.code sigma_method = \"mle\"} estimates the scale of the
            {.val {family}} conditional density from the residuals of the
-           propensity score model, and {.arg .sigma} is a spread of your own
+           model of the exposure, and {.arg .sigma} is a spread of your own
            that replaces it.",
       i = "A {.val {family}} density is spread by its own estimator unless you
            ask for another, so {.code sigma_method = \"mle\"} is what it
@@ -984,8 +984,12 @@ continuous_density_ratio <- function(
 
   # The marginal density is the same family read at the exposure's own center
   # and spread, so it is standardized by those rather than by the conditional
-  # spread. A kernel numerator is therefore fit on its own values.
+  # spread. A kernel numerator is therefore fit on its own values, and only on
+  # those of the units the conditional density reads: a unit with no fitted
+  # mean, such as one a trim set aside, is outside the population the moments
+  # were read over and carries no weight, so it has no part in the fit either.
   z_a <- (exposure - mu_a) / sigma_a
+  z_a[is.na(z)] <- NA_real_
   f_num <- density_eval_present(density, z_a, call = call) / sigma_a
 
   f_num / f_den

R/ipw-continuous-models.R

@@ -911,9 +911,10 @@ check_ipw_stabilization_score <- function(
     c(
       headline,
       x = paste("They weight {n} observation{?s} and", held),
-      i = "A per-observation score is one value for each unit, so it does not \\
-      survive the rows being restricted: subsetting the weights drops it, and \\
-      a model frame that drops incomplete rows leaves it at the length the \\
+      i = "A per-observation score is one value for each unit. Subsetting the \\
+      weights with {.code [} keeps it aligned, but a slice that cannot place \\
+      it, such as {.fn vctrs::vec_slice} or {.fn dplyr::arrange}, drops it, \\
+      and a model frame that drops incomplete rows leaves it at the length the \\
       weights were built at.",
       i = "Rebuild the weights on the rows being analyzed, or stabilize on a \\
       single {.arg stabilization_score}, which scales every weight and \\