diff --git a/DESCRIPTION b/DESCRIPTION
index 0480f73..f3e7123 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -15,16 +15,18 @@ Imports:
     cli,
     rlang,
     stats,
+    tidyselect,
     vctrs (>= 0.6.5)
 Suggests: 
     dplyr,
     ggplot2,
     mgcv,
+    parsnip,
     probably,
     PSweight,
     spelling,
     testthat (>= 3.0.0),
-    tidyselect,
+    tibble,
     WeightIt
 Config/testthat/edition: 3
 Encoding: UTF-8
diff --git a/NAMESPACE b/NAMESPACE
index 8fc5110..3cd7ccc 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -81,21 +81,39 @@ S3method(vec_ptype_abbr,psw)
 S3method(vec_ptype_full,ps_trim)
 S3method(vec_ptype_full,ps_trunc)
 S3method(vec_ptype_full,psw)
+S3method(wt_ate,data.frame)
+S3method(wt_ate,default)
+S3method(wt_ate,glm)
 S3method(wt_ate,numeric)
 S3method(wt_ate,ps_trim)
 S3method(wt_ate,ps_trunc)
+S3method(wt_atm,data.frame)
+S3method(wt_atm,default)
+S3method(wt_atm,glm)
 S3method(wt_atm,numeric)
 S3method(wt_atm,ps_trim)
 S3method(wt_atm,ps_trunc)
+S3method(wt_ato,data.frame)
+S3method(wt_ato,default)
+S3method(wt_ato,glm)
 S3method(wt_ato,numeric)
 S3method(wt_ato,ps_trim)
 S3method(wt_ato,ps_trunc)
+S3method(wt_att,data.frame)
+S3method(wt_att,default)
+S3method(wt_att,glm)
 S3method(wt_att,numeric)
 S3method(wt_att,ps_trim)
 S3method(wt_att,ps_trunc)
+S3method(wt_atu,data.frame)
+S3method(wt_atu,default)
+S3method(wt_atu,glm)
 S3method(wt_atu,numeric)
 S3method(wt_atu,ps_trim)
 S3method(wt_atu,ps_trunc)
+S3method(wt_entropy,data.frame)
+S3method(wt_entropy,default)
+S3method(wt_entropy,glm)
 S3method(wt_entropy,numeric)
 S3method(wt_entropy,ps_trim)
 S3method(wt_entropy,ps_trunc)
diff --git a/R/weights-utils.R b/R/weights-utils.R
index 27f0aed..39e18d0 100644
--- a/R/weights-utils.R
+++ b/R/weights-utils.R
@@ -8,9 +8,10 @@ abort_unsupported <- function(exposure_type, what, call = rlang::caller_env()) {
 
 match_exposure_type <- function(
   exposure_type = c("auto", "binary", "categorical", "continuous"),
-  .exposure
+  .exposure,
+  valid_types = c("auto", "binary", "categorical", "continuous")
 ) {
-  .exposure_type <- rlang::arg_match(exposure_type)
+  .exposure_type <- rlang::arg_match(exposure_type, valid_types)
   if (.exposure_type == "auto") {
     detect_exposure_type(.exposure)
   } else {
@@ -128,3 +129,26 @@ check_lower_upper <- function(lower, upper, call = rlang::caller_env()) {
 
   invisible(TRUE)
 }
+
+check_lengths_match <- function(
+  .propensity,
+  .exposure,
+  call = rlang::caller_env()
+) {
+  len_prop <- length(.propensity)
+  len_exp <- length(.exposure)
+
+  if (len_prop != len_exp) {
+    abort(
+      c(
+        "{.arg .propensity} and {.arg .exposure} must have the same length.",
+        i = "{.arg .propensity} has length {len_prop}",
+        i = "{.arg .exposure} has length {len_exp}"
+      ),
+      call = call,
+      error_class = "propensity_length_error"
+    )
+  }
+
+  invisible(TRUE)
+}
diff --git a/R/weights.R b/R/weights.R
index bca3363..cf3a64e 100644
--- a/R/weights.R
+++ b/R/weights.R
@@ -12,10 +12,10 @@
 #' - **Entropy** (Average Treatment Effect for the Entropy-weighted population)
 #'
 #'   The propensity score can be provided as a numeric vector of predicted
-#'   probabilities or as a `data.frame` where each column represents the
-#'   predicted probability for a level of the exposure. They can also be
-#'   propensity score objects created by [ps_trim()], [ps_refit()], or
-#'   [ps_trunc()]
+#'   probabilities, as a `data.frame` where each column represents the
+#'   predicted probability for a level of the exposure, or as a fitted
+#'   GLM object. They can also be propensity score objects created by
+#'   [ps_trim()], [ps_refit()], or [ps_trunc()]
 #'
 #'   The returned weights are encapsulated in a `psw` object, which is a numeric
 #'   vector with additional attributes that record the estimand, and whether the
@@ -75,11 +75,11 @@
 #' - Positivity violations (near 0 or 1 propensity scores)
 #' - Poor model specification
 #' - Lack of overlap between treatment groups
-#' 
+#'
 #' See the halfmoon package for tools to diagnose and visualize weights.
-#' 
+#'
 #' You can address extreme weights in several ways. The first is to modify the target population:
-#' use trimming, truncation, or alternative estimands (ATM, ATO, entropy). 
+#' use trimming, truncation, or alternative estimands (ATM, ATO, entropy).
 #' Another technique that can help is stabilization, which reduces variance of the weights.
 #'
 #' ## Trimmed and Truncated Weights
@@ -90,8 +90,12 @@
 #' modified propensity scores (trimmed or truncated) and update the estimand
 #' attribute accordingly.
 #'
-#' @param .propensity Either a numeric vector of predicted probabilities or a
-#'   `data.frame` where each column corresponds to a level of the exposure.
+#' @param .propensity Either a numeric vector of predicted probabilities, a
+#'   `data.frame` where each column corresponds to a level of the exposure,
+#'   or a fitted GLM object. For data frames, the second column is used by
+#'   default for binary exposures unless specified otherwise with
+#'   `.propensity_col`. For GLM objects, fitted values are extracted
+#'   automatically.
 #' @param .exposure The exposure variable. For binary exposures, a vector of 0s
 #'   and 1s; for continuous exposures, a numeric vector.
 #' @param exposure_type Character string specifying the type of exposure.
@@ -111,6 +115,10 @@
 #' @param stabilization_score Optional numeric value for stabilizing the weights
 #'   (e.g., a predicted value from a regression model without predictors). Only
 #'   used when `stabilize` is `TRUE`.
+#' @param .propensity_col When `.propensity` is a data frame, specifies which
+#'   column to use for propensity scores. Can be a column name (quoted or
+#'   unquoted) or a numeric index. Defaults to the second column if available,
+#'   otherwise the first.
 #'
 #' @return A `psw` object (a numeric vector) with additional attributes:
 #'   - **estimand**: A description of the estimand (e.g., "ate", "att").
@@ -157,20 +165,52 @@
 #' # Standard ATE weights can be extreme
 #' wt_extreme <- wt_ate(ps_extreme, trt_extreme)
 #' # Very large!
-#' max(wt_extreme)  
+#' max(wt_extreme)
 #'
 #' # ATO weights are bounded
 #' wt_extreme_atm <- wt_ato(ps_extreme, trt_extreme)
 #' # Much more reasonable
-#' max(wt_extreme_atm)  
+#' max(wt_extreme_atm)
 #' # but they target a different population
 #' estimand(wt_extreme_atm) # "ato"
 #'
+#' ## Working with Data Frames
+#'
+#' # Example with custom data frame
+#' ps_df <- data.frame(
+#'   control = c(0.9, 0.7, 0.3, 0.1),
+#'   treated = c(0.1, 0.3, 0.7, 0.9)
+#' )
+#' exposure <- c(0, 0, 1, 1)
+#'
+#' # Uses second column by default (treated probabilities)
+#' wt_ate(ps_df, exposure)
+#'
+#' # Explicitly specify column by name
+#' wt_ate(ps_df, exposure, .propensity_col = "treated")
+#'
+#' # Or by position
+#' wt_ate(ps_df, exposure, .propensity_col = 2)
+#'
+#' ## Working with GLM Objects
+#'
+#' # Fit a propensity score model
+#' set.seed(123)
+#' n <- 100
+#' x1 <- rnorm(n)
+#' x2 <- rnorm(n)
+#' treatment <- rbinom(n, 1, plogis(0.5 * x1 + 0.3 * x2))
+#'
+#' ps_model <- glm(treatment ~ x1 + x2, family = binomial)
+#'
+#' # Use GLM directly for weight calculation
+#' weights_from_glm <- wt_ate(ps_model, treatment)
+#'
 #' @references
-#' 
-#' For detailed guidance on causal inference in R, see [*Causal Inference in R*](https://www.r-causal.org/) 
+#'
+#' For detailed guidance on causal inference in R, see [*Causal Inference in R*](https://www.r-causal.org/)
 #' by Malcolm Barrett, Lucy D'Agostino McGowan, and Travis Gerke.
-#' 
+#'
 #' ## Foundational Papers
 #'
 #' Rosenbaum, P. R., & Rubin, D. B. (1983). The central role of the propensity
@@ -223,6 +263,27 @@ wt_ate <- function(
   UseMethod("wt_ate")
 }
 
+#' @export
+wt_ate.default <- function(
+  .propensity,
+  .exposure,
+  .sigma = NULL,
+  exposure_type = c("auto", "binary", "categorical", "continuous"),
+  .treated = NULL,
+  .untreated = NULL,
+  stabilize = FALSE,
+  stabilization_score = NULL,
+  ...
+) {
+  abort(
+    paste0(
+      "No method for objects of class ",
+      paste(class(.propensity), collapse = ", ")
+    ),
+    error_class = "propensity_method_error"
+  )
+}
+
 #' @export
 wt_ate.numeric <- function(
   .propensity,
@@ -236,6 +297,7 @@ wt_ate.numeric <- function(
   ...
 ) {
   rlang::check_dots_empty()
+  check_lengths_match(.propensity, .exposure)
   exposure_type <- match_exposure_type(exposure_type, .exposure)
   if (exposure_type == "binary") {
     check_ps_range(.propensity)
@@ -262,6 +324,146 @@ wt_ate.numeric <- function(
   psw(wts, "ate", stabilized = isTRUE(stabilize))
 }
 
+#' @export
+#' @rdname wt_ate
+wt_ate.data.frame <- function(
+  .propensity,
+  .exposure,
+  .sigma = NULL,
+  exposure_type = c("auto", "binary", "categorical", "continuous"),
+  .treated = NULL,
+  .untreated = NULL,
+  stabilize = FALSE,
+  stabilization_score = NULL,
+  ...,
+  .propensity_col = NULL
+) {
+  # Capture the column selection expression
+  col_quo <- rlang::enquo(.propensity_col)
+
+  # Extract propensity scores from data frame
+  ps_vec <- extract_propensity_from_df(.propensity, col_quo)
+
+  # Call the numeric method
+  wt_ate.numeric(
+    .propensity = ps_vec,
+    .exposure = .exposure,
+    .sigma = .sigma,
+    exposure_type = exposure_type,
+    .treated = .treated,
+    .untreated = .untreated,
+    stabilize = stabilize,
+    stabilization_score = stabilization_score,
+    ...
+  )
+}
+
+#' @export
+wt_ate.glm <- function(
+  .propensity,
+  .exposure,
+  .sigma = NULL,
+  exposure_type = c("auto", "binary", "categorical", "continuous"),
+  .treated = NULL,
+  .untreated = NULL,
+  stabilize = FALSE,
+  stabilization_score = NULL,
+  ...
+) {
+  # Extract fitted values (propensity scores) from GLM
+  ps_vec <- extract_propensity_from_glm(.propensity)
+
+  # For continuous exposures, extract sigma if not provided
+  if (
+    is.null(.sigma) &&
+      match_exposure_type(exposure_type, .exposure) == "continuous"
+  ) {
+    .sigma <- stats::influence(.propensity)$sigma
+  }
+
+  # Call the numeric method
+  wt_ate.numeric(
+    .propensity = ps_vec,
+    .exposure = .exposure,
+    .sigma = .sigma,
+    exposure_type = exposure_type,
+    .treated = .treated,
+    .untreated = .untreated,
+    stabilize = stabilize,
+    stabilization_score = stabilization_score,
+    ...
+  )
+}
+
+# Helper function to extract propensity scores from data frames
+extract_propensity_from_df <- function(
+  .propensity,
+  .propensity_col_quo = NULL
+) {
+  if (!rlang::quo_is_null(.propensity_col_quo)) {
+    col_pos <- tryCatch(
+      tidyselect::eval_select(
+        .propensity_col_quo,
+        data = .propensity
+      ),
+      error = function(e) {
+        abort(
+          paste0("Column selection failed: ", e$message),
+          error_class = "propensity_df_column_error"
+        )
+      }
+    )
+
+    if (length(col_pos) != 1) {
+      abort(
+        "`.propensity_col` must select exactly one column.",
+        error_class = "propensity_df_column_error"
+      )
+    }
+
+    ps_vec <- .propensity[[col_pos]]
+  } else {
+    # Default behavior: use second column if available, otherwise first
+    if (ncol(.propensity) >= 2) {
+      ps_vec <- .propensity[[2]]
+    } else if (ncol(.propensity) == 1) {
+      ps_vec <- .propensity[[1]]
+    } else {
+      abort(
+        "`.propensity` data frame must have at least one column.",
+        error_class = "propensity_df_ncol_error"
+      )
+    }
+  }
+
+  ps_vec
+}
+
+# Helper function to extract propensity scores from GLM objects
+extract_propensity_from_glm <- function(.propensity) {
+  # Check if it's a valid GLM object
+  if (!inherits(.propensity, "glm")) {
+    abort(
+      "`.propensity` must be a GLM object.",
+      error_class = "propensity_glm_type_error"
+    )
+  }
+
+  # Check if it's a binomial GLM for binary propensity scores
+  if (
+    !is.null(.propensity$family) &&
+      .propensity$family$family == "binomial"
+  ) {
+    # Get predicted probabilities
+    ps_vec <- stats::predict(.propensity, type = "response")
+  } else {
+    # For non-binomial GLMs, get linear predictor
+    ps_vec <- stats::fitted(.propensity)
+  }
+
+  ps_vec
+}
+
 ate_binary <- function(
   .propensity,
   .exposure,
@@ -340,7 +542,7 @@ ate_continuous <- function(
 wt_att <- function(
   .propensity,
   .exposure,
-  exposure_type = c("auto", "binary", "categorical", "continuous"),
+  exposure_type = c("auto", "binary", "categorical"),
   .treated = NULL,
   .untreated = NULL,
   ...
@@ -348,17 +550,40 @@ wt_att <- function(
   UseMethod("wt_att")
 }
 
+#' @export
+wt_att.default <- function(
+  .propensity,
+  .exposure,
+  exposure_type = c("auto", "binary", "categorical"),
+  .treated = NULL,
+  .untreated = NULL,
+  ...
+) {
+  abort(
+    paste0(
+      "No method for objects of class ",
+      paste(class(.propensity), collapse = ", ")
+    ),
+    error_class = "propensity_method_error"
+  )
+}
+
 #' @export
 wt_att.numeric <- function(
   .propensity,
   .exposure,
-  exposure_type = c("auto", "binary", "categorical", "continuous"),
+  exposure_type = c("auto", "binary", "categorical"),
   .treated = NULL,
   .untreated = NULL,
   ...
 ) {
   rlang::check_dots_empty()
-  exposure_type <- match_exposure_type(exposure_type, .exposure)
+  check_lengths_match(.propensity, .exposure)
+  exposure_type <- match_exposure_type(
+    exposure_type,
+    .exposure,
+    c("auto", "binary", "categorical")
+  )
 
   if (exposure_type == "binary") {
     check_ps_range(.propensity)
@@ -375,6 +600,57 @@ wt_att.numeric <- function(
   psw(wts, "att")
 }
 
+#' @export
+#' @rdname wt_ate
+wt_att.data.frame <- function(
+  .propensity,
+  .exposure,
+  exposure_type = c("auto", "binary", "categorical"),
+  .treated = NULL,
+  .untreated = NULL,
+  ...,
+  .propensity_col = NULL
+) {
+  # Capture the column selection expression
+  col_quo <- rlang::enquo(.propensity_col)
+
+  # Extract propensity scores from data frame
+  ps_vec <- extract_propensity_from_df(.propensity, col_quo)
+
+  # Call the numeric method
+  wt_att.numeric(
+    .propensity = ps_vec,
+    .exposure = .exposure,
+    exposure_type = exposure_type,
+    .treated = .treated,
+    .untreated = .untreated,
+    ...
+  )
+}
+
+#' @export
+wt_att.glm <- function(
+  .propensity,
+  .exposure,
+  exposure_type = c("auto", "binary", "categorical"),
+  .treated = NULL,
+  .untreated = NULL,
+  ...
+) {
+  # Extract fitted values (propensity scores) from GLM
+  ps_vec <- extract_propensity_from_glm(.propensity)
+
+  # Call the numeric method
+  wt_att.numeric(
+    .propensity = ps_vec,
+    .exposure = .exposure,
+    exposure_type = exposure_type,
+    .treated = .treated,
+    .untreated = .untreated,
+    ...
+  )
+}
+
 att_binary <- function(
   .propensity,
   .exposure,
@@ -396,7 +672,7 @@ att_binary <- function(
 wt_atu <- function(
   .propensity,
   .exposure,
-  exposure_type = c("auto", "binary", "categorical", "continuous"),
+  exposure_type = c("auto", "binary", "categorical"),
   .treated = NULL,
   .untreated = NULL,
   ...
@@ -404,17 +680,40 @@ wt_atu <- function(
   UseMethod("wt_atu")
 }
 
+#' @export
+wt_atu.default <- function(
+  .propensity,
+  .exposure,
+  exposure_type = c("auto", "binary", "categorical"),
+  .treated = NULL,
+  .untreated = NULL,
+  ...
+) {
+  abort(
+    paste0(
+      "No method for objects of class ",
+      paste(class(.propensity), collapse = ", ")
+    ),
+    error_class = "propensity_method_error"
+  )
+}
+
 #' @export
 wt_atu.numeric <- function(
   .propensity,
   .exposure,
-  exposure_type = c("auto", "binary", "categorical", "continuous"),
+  exposure_type = c("auto", "binary", "categorical"),
   .treated = NULL,
   .untreated = NULL,
   ...
 ) {
   rlang::check_dots_empty()
-  exposure_type <- match_exposure_type(exposure_type, .exposure)
+  check_lengths_match(.propensity, .exposure)
+  exposure_type <- match_exposure_type(
+    exposure_type,
+    .exposure,
+    c("auto", "binary", "categorical")
+  )
   if (exposure_type == "binary") {
     check_ps_range(.propensity)
     wts <- atu_binary(
@@ -430,6 +729,57 @@ wt_atu.numeric <- function(
   psw(wts, "atu")
 }
 
+#' @export
+#' @rdname wt_ate
+wt_atu.data.frame <- function(
+  .propensity,
+  .exposure,
+  exposure_type = c("auto", "binary", "categorical"),
+  .treated = NULL,
+  .untreated = NULL,
+  ...,
+  .propensity_col = NULL
+) {
+  # Capture the column selection expression
+  col_quo <- rlang::enquo(.propensity_col)
+
+  # Extract propensity scores from data frame
+  ps_vec <- extract_propensity_from_df(.propensity, col_quo)
+
+  # Call the numeric method
+  wt_atu.numeric(
+    .propensity = ps_vec,
+    .exposure = .exposure,
+    exposure_type = exposure_type,
+    .treated = .treated,
+    .untreated = .untreated,
+    ...
+  )
+}
+
+#' @export
+wt_atu.glm <- function(
+  .propensity,
+  .exposure,
+  exposure_type = c("auto", "binary", "categorical"),
+  .treated = NULL,
+  .untreated = NULL,
+  ...
+) {
+  # Extract fitted values (propensity scores) from GLM
+  ps_vec <- extract_propensity_from_glm(.propensity)
+
+  # Call the numeric method
+  wt_atu.numeric(
+    .propensity = ps_vec,
+    .exposure = .exposure,
+    exposure_type = exposure_type,
+    .treated = .treated,
+    .untreated = .untreated,
+    ...
+  )
+}
+
 atu_binary <- function(
   .propensity,
   .exposure,
@@ -453,7 +803,7 @@ atu_binary <- function(
 wt_atm <- function(
   .propensity,
   .exposure,
-  exposure_type = c("auto", "binary", "categorical", "continuous"),
+  exposure_type = c("auto", "binary", "categorical"),
   .treated = NULL,
   .untreated = NULL,
   ...
@@ -461,17 +811,40 @@ wt_atm <- function(
   UseMethod("wt_atm")
 }
 
+#' @export
+wt_atm.default <- function(
+  .propensity,
+  .exposure,
+  exposure_type = c("auto", "binary", "categorical"),
+  .treated = NULL,
+  .untreated = NULL,
+  ...
+) {
+  abort(
+    paste0(
+      "No method for objects of class ",
+      paste(class(.propensity), collapse = ", ")
+    ),
+    error_class = "propensity_method_error"
+  )
+}
+
 #' @export
 wt_atm.numeric <- function(
   .propensity,
   .exposure,
-  exposure_type = c("auto", "binary", "categorical", "continuous"),
+  exposure_type = c("auto", "binary", "categorical"),
   .treated = NULL,
   .untreated = NULL,
   ...
 ) {
   rlang::check_dots_empty()
-  exposure_type <- match_exposure_type(exposure_type, .exposure)
+  check_lengths_match(.propensity, .exposure)
+  exposure_type <- match_exposure_type(
+    exposure_type,
+    .exposure,
+    c("auto", "binary", "categorical")
+  )
   if (exposure_type == "binary") {
     check_ps_range(.propensity)
     wts <- atm_binary(
@@ -487,6 +860,57 @@ wt_atm.numeric <- function(
   psw(wts, "atm")
 }
 
+#' @export
+#' @rdname wt_ate
+wt_atm.data.frame <- function(
+  .propensity,
+  .exposure,
+  exposure_type = c("auto", "binary", "categorical"),
+  .treated = NULL,
+  .untreated = NULL,
+  ...,
+  .propensity_col = NULL
+) {
+  # Capture the column selection expression
+  col_quo <- rlang::enquo(.propensity_col)
+
+  # Extract propensity scores from data frame
+  ps_vec <- extract_propensity_from_df(.propensity, col_quo)
+
+  # Call the numeric method
+  wt_atm.numeric(
+    .propensity = ps_vec,
+    .exposure = .exposure,
+    exposure_type = exposure_type,
+    .treated = .treated,
+    .untreated = .untreated,
+    ...
+  )
+}
+
+#' @export
+wt_atm.glm <- function(
+  .propensity,
+  .exposure,
+  exposure_type = c("auto", "binary", "categorical"),
+  .treated = NULL,
+  .untreated = NULL,
+  ...
+) {
+  # Extract fitted values (propensity scores) from GLM
+  ps_vec <- extract_propensity_from_glm(.propensity)
+
+  # Call the numeric method
+  wt_atm.numeric(
+    .propensity = ps_vec,
+    .exposure = .exposure,
+    exposure_type = exposure_type,
+    .treated = .treated,
+    .untreated = .untreated,
+    ...
+  )
+}
+
 atm_binary <- function(
   .propensity,
   .exposure,
@@ -509,7 +933,7 @@ atm_binary <- function(
 wt_ato <- function(
   .propensity,
   .exposure,
-  exposure_type = c("auto", "binary", "categorical", "continuous"),
+  exposure_type = c("auto", "binary", "categorical"),
   .treated = NULL,
   .untreated = NULL,
   ...
@@ -517,16 +941,38 @@ wt_ato <- function(
   UseMethod("wt_ato")
 }
 
+#' @export
+wt_ato.default <- function(
+  .propensity,
+  .exposure,
+  exposure_type = c("auto", "binary", "categorical"),
+  .treated = NULL,
+  .untreated = NULL,
+  ...
+) {
+  abort(
+    paste0(
+      "No method for objects of class ",
+      paste(class(.propensity), collapse = ", ")
+    ),
+    error_class = "propensity_method_error"
+  )
+}
+
 #' @export
 wt_ato.numeric <- function(
   .propensity,
   .exposure,
-  exposure_type = c("auto", "binary", "categorical", "continuous"),
+  exposure_type = c("auto", "binary", "categorical"),
   .treated = NULL,
   .untreated = NULL,
   ...
 ) {
-  exposure_type <- match_exposure_type(exposure_type, .exposure)
+  exposure_type <- match_exposure_type(
+    exposure_type,
+    .exposure,
+    c("auto", "binary", "categorical")
+  )
   if (exposure_type == "binary") {
     check_ps_range(.propensity)
     wts <- ato_binary(
@@ -542,6 +988,56 @@ wt_ato.numeric <- function(
   psw(wts, "ato")
 }
 
+#' @export
+#' @rdname wt_ate
+wt_ato.data.frame <- function(
+  .propensity,
+  .exposure,
+  exposure_type = c("auto", "binary", "categorical"),
+  .treated = NULL,
+  .untreated = NULL,
+  ...,
+  .propensity_col = NULL
+) {
+  # Capture the column selection expression
+  col_quo <- rlang::enquo(.propensity_col)
+
+  # Extract propensity scores from data frame
+  ps_vec <- extract_propensity_from_df(.propensity, col_quo)
+
+  # Call the numeric method
+  wt_ato.numeric(
+    .propensity = ps_vec,
+    .exposure = .exposure,
+    exposure_type = exposure_type,
+    .treated = .treated,
+    .untreated = .untreated,
+    ...
+  )
+}
+
+#' @export
+wt_ato.glm <- function(
+  .propensity,
+  .exposure,
+  exposure_type = c("auto", "binary", "categorical"),
+  .treated = NULL,
+  .untreated = NULL,
+  ...
+) {
+  # Extract fitted values (propensity scores) from GLM
+  ps_vec <- extract_propensity_from_glm(.propensity)
+
+  # Call the numeric method
+  wt_ato.numeric(
+    .propensity = ps_vec,
+    .exposure = .exposure,
+    exposure_type = exposure_type,
+    .treated = .treated,
+    .untreated = .untreated,
+    ...
+  )
+}
 
 ato_binary <- function(
   .propensity,
@@ -563,7 +1059,7 @@ ato_binary <- function(
 wt_entropy <- function(
   .propensity,
   .exposure,
-  exposure_type = c("auto", "binary", "categorical", "continuous"),
+  exposure_type = c("auto", "binary", "categorical"),
   .treated = NULL,
   .untreated = NULL,
   ...
@@ -571,17 +1067,40 @@ wt_entropy <- function(
   UseMethod("wt_entropy")
 }
 
+#' @export
+wt_entropy.default <- function(
+  .propensity,
+  .exposure,
+  exposure_type = c("auto", "binary", "categorical"),
+  .treated = NULL,
+  .untreated = NULL,
+  ...
+) {
+  abort(
+    paste0(
+      "No method for objects of class ",
+      paste(class(.propensity), collapse = ", ")
+    ),
+    error_class = "propensity_method_error"
+  )
+}
+
 #' @export
 wt_entropy.numeric <- function(
   .propensity,
   .exposure,
-  exposure_type = c("auto", "binary", "categorical", "continuous"),
+  exposure_type = c("auto", "binary", "categorical"),
   .treated = NULL,
   .untreated = NULL,
   ...
 ) {
   rlang::check_dots_empty()
-  exposure_type <- match_exposure_type(exposure_type, .exposure)
+  check_lengths_match(.propensity, .exposure)
+  exposure_type <- match_exposure_type(
+    exposure_type,
+    .exposure,
+    c("auto", "binary", "categorical")
+  )
   if (exposure_type == "binary") {
     check_ps_range(.propensity)
     wts <- entropy_binary(
@@ -597,6 +1116,57 @@ wt_entropy.numeric <- function(
   psw(wts, "entropy")
 }
 
+#' @export
+#' @rdname wt_ate
+wt_entropy.data.frame <- function(
+  .propensity,
+  .exposure,
+  exposure_type = c("auto", "binary", "categorical"),
+  .treated = NULL,
+  .untreated = NULL,
+  ...,
+  .propensity_col = NULL
+) {
+  # Capture the column selection expression
+  col_quo <- rlang::enquo(.propensity_col)
+
+  # Extract propensity scores from data frame
+  ps_vec <- extract_propensity_from_df(.propensity, col_quo)
+
+  # Call the numeric method
+  wt_entropy.numeric(
+    .propensity = ps_vec,
+    .exposure = .exposure,
+    exposure_type = exposure_type,
+    .treated = .treated,
+    .untreated = .untreated,
+    ...
+  )
+}
+
+#' @export
+wt_entropy.glm <- function(
+  .propensity,
+  .exposure,
+  exposure_type = c("auto", "binary", "categorical"),
+  .treated = NULL,
+  .untreated = NULL,
+  ...
+) {
+  # Extract fitted values (propensity scores) from GLM
+  ps_vec <- extract_propensity_from_glm(.propensity)
+
+  # Call the numeric method
+  wt_entropy.numeric(
+    .propensity = ps_vec,
+    .exposure = .exposure,
+    exposure_type = exposure_type,
+    .treated = .treated,
+    .untreated = .untreated,
+    ...
+  )
+}
+
 entropy_binary <- function(
   .propensity,
   .exposure,
@@ -665,7 +1235,7 @@ wt_ate.ps_trim <- function(
 wt_att.ps_trim <- function(
   .propensity,
   .exposure,
-  exposure_type = c("auto", "binary", "categorical", "continuous"),
+  exposure_type = c("auto", "binary", "categorical"),
   .treated = NULL,
   .untreated = NULL,
   ...
@@ -694,7 +1264,7 @@ wt_att.ps_trim <- function(
 wt_atu.ps_trim <- function(
   .propensity,
   .exposure,
-  exposure_type = c("auto", "binary", "categorical", "continuous"),
+  exposure_type = c("auto", "binary", "categorical"),
   .treated = NULL,
   .untreated = NULL,
   ...
@@ -723,7 +1293,7 @@ wt_atu.ps_trim <- function(
 wt_atm.ps_trim <- function(
   .propensity,
   .exposure,
-  exposure_type = c("auto", "binary", "categorical", "continuous"),
+  exposure_type = c("auto", "binary", "categorical"),
   .treated = NULL,
   .untreated = NULL,
   ...
@@ -752,7 +1322,7 @@ wt_atm.ps_trim <- function(
 wt_ato.ps_trim <- function(
   .propensity,
   .exposure,
-  exposure_type = c("auto", "binary", "categorical", "continuous"),
+  exposure_type = c("auto", "binary", "categorical"),
   .treated = NULL,
   .untreated = NULL,
   ...
@@ -814,7 +1384,7 @@ wt_ate.ps_trunc <- function(
 wt_att.ps_trunc <- function(
   .propensity,
   .exposure,
-  exposure_type = c("auto", "binary", "categorical", "continuous"),
+  exposure_type = c("auto", "binary", "categorical"),
   .treated = NULL,
   .untreated = NULL,
   ...
@@ -841,7 +1411,7 @@ wt_att.ps_trunc <- function(
 wt_atu.ps_trunc <- function(
   .propensity,
   .exposure,
-  exposure_type = c("auto", "binary", "categorical", "continuous"),
+  exposure_type = c("auto", "binary", "categorical"),
   .treated = NULL,
   .untreated = NULL,
   ...
@@ -868,7 +1438,7 @@ wt_atu.ps_trunc <- function(
 wt_atm.ps_trunc <- function(
   .propensity,
   .exposure,
-  exposure_type = c("auto", "binary", "categorical", "continuous"),
+  exposure_type = c("auto", "binary", "categorical"),
   .treated = NULL,
   .untreated = NULL,
   ...
@@ -895,7 +1465,7 @@ wt_atm.ps_trunc <- function(
 wt_ato.ps_trunc <- function(
   .propensity,
   .exposure,
-  exposure_type = c("auto", "binary", "categorical", "continuous"),
+  exposure_type = c("auto", "binary", "categorical"),
   .treated = NULL,
   .untreated = NULL,
   ...
@@ -922,7 +1492,7 @@ wt_ato.ps_trunc <- function(
 wt_entropy.ps_trim <- function(
   .propensity,
   .exposure,
-  exposure_type = c("auto", "binary", "categorical"),
+  exposure_type = c("auto", "binary"),
   .treated = NULL,
   .untreated = NULL,
   ...
@@ -951,7 +1521,7 @@ wt_entropy.ps_trim <- function(
 wt_entropy.ps_trunc <- function(
   .propensity,
   .exposure,
-  exposure_type = c("auto", "binary", "categorical"),
+  exposure_type = c("auto", "binary"),
   .treated = NULL,
   .untreated = NULL,
   ...
diff --git a/man/wt_ate.Rd b/man/wt_ate.Rd
index 183d79c..320cefa 100644
--- a/man/wt_ate.Rd
+++ b/man/wt_ate.Rd
@@ -2,11 +2,17 @@
 % Please edit documentation in R/weights.R
 \name{wt_ate}
 \alias{wt_ate}
+\alias{wt_ate.data.frame}
 \alias{wt_att}
+\alias{wt_att.data.frame}
 \alias{wt_atu}
+\alias{wt_atu.data.frame}
 \alias{wt_atm}
+\alias{wt_atm.data.frame}
 \alias{wt_ato}
+\alias{wt_ato.data.frame}
 \alias{wt_entropy}
+\alias{wt_entropy.data.frame}
 \title{Calculate Propensity Score Weights for Causal Inference}
 \usage{
 wt_ate(
@@ -21,54 +27,121 @@ wt_ate(
   ...
 )
 
-wt_att(
+\method{wt_ate}{data.frame}(
   .propensity,
   .exposure,
+  .sigma = NULL,
   exposure_type = c("auto", "binary", "categorical", "continuous"),
   .treated = NULL,
   .untreated = NULL,
+  stabilize = FALSE,
+  stabilization_score = NULL,
+  ...,
+  .propensity_col = NULL
+)
+
+wt_att(
+  .propensity,
+  .exposure,
+  exposure_type = c("auto", "binary", "categorical"),
+  .treated = NULL,
+  .untreated = NULL,
   ...
 )
 
+\method{wt_att}{data.frame}(
+  .propensity,
+  .exposure,
+  exposure_type = c("auto", "binary", "categorical"),
+  .treated = NULL,
+  .untreated = NULL,
+  ...,
+  .propensity_col = NULL
+)
+
 wt_atu(
   .propensity,
   .exposure,
-  exposure_type = c("auto", "binary", "categorical", "continuous"),
+  exposure_type = c("auto", "binary", "categorical"),
   .treated = NULL,
   .untreated = NULL,
   ...
 )
 
+\method{wt_atu}{data.frame}(
+  .propensity,
+  .exposure,
+  exposure_type = c("auto", "binary", "categorical"),
+  .treated = NULL,
+  .untreated = NULL,
+  ...,
+  .propensity_col = NULL
+)
+
 wt_atm(
   .propensity,
   .exposure,
-  exposure_type = c("auto", "binary", "categorical", "continuous"),
+  exposure_type = c("auto", "binary", "categorical"),
   .treated = NULL,
   .untreated = NULL,
   ...
 )
 
+\method{wt_atm}{data.frame}(
+  .propensity,
+  .exposure,
+  exposure_type = c("auto", "binary", "categorical"),
+  .treated = NULL,
+  .untreated = NULL,
+  ...,
+  .propensity_col = NULL
+)
+
 wt_ato(
   .propensity,
   .exposure,
-  exposure_type = c("auto", "binary", "categorical", "continuous"),
+  exposure_type = c("auto", "binary", "categorical"),
   .treated = NULL,
   .untreated = NULL,
   ...
 )
 
+\method{wt_ato}{data.frame}(
+  .propensity,
+  .exposure,
+  exposure_type = c("auto", "binary", "categorical"),
+  .treated = NULL,
+  .untreated = NULL,
+  ...,
+  .propensity_col = NULL
+)
+
 wt_entropy(
   .propensity,
   .exposure,
-  exposure_type = c("auto", "binary", "categorical", "continuous"),
+  exposure_type = c("auto", "binary", "categorical"),
   .treated = NULL,
   .untreated = NULL,
   ...
 )
+
+\method{wt_entropy}{data.frame}(
+  .propensity,
+  .exposure,
+  exposure_type = c("auto", "binary", "categorical"),
+  .treated = NULL,
+  .untreated = NULL,
+  ...,
+  .propensity_col = NULL
+)
 }
 \arguments{
-\item{.propensity}{Either a numeric vector of predicted probabilities or a
-\code{data.frame} where each column corresponds to a level of the exposure.}
+\item{.propensity}{Either a numeric vector of predicted probabilities, a
+\code{data.frame} where each column corresponds to a level of the exposure,
+or a fitted GLM object. For data frames, the second column is used by
+default for binary exposures unless specified otherwise with
+\code{.propensity_col}. For GLM objects, fitted values are extracted
+automatically.}
 
 \item{.exposure}{The exposure variable. For binary exposures, a vector of 0s
 and 1s; for continuous exposures, a numeric vector.}
@@ -96,6 +169,11 @@ weights, stabilization multiplies the weight by either the mean of
 used when \code{stabilize} is \code{TRUE}.}
 
 \item{...}{Reserved for future expansion. Not currently used.}
+
+\item{.propensity_col}{When \code{.propensity} is a data frame, specifies which
+column to use for propensity scores. Can be a column name (quoted or
+unquoted) or a numeric index. Defaults to the second column if available,
+otherwise the first.}
 }
 \value{
 A \code{psw} object (a numeric vector) with additional attributes:
@@ -119,10 +197,10 @@ the \strong{ATC}, where the "C" stands for "control")
 \item \strong{Entropy} (Average Treatment Effect for the Entropy-weighted population)
 
 The propensity score can be provided as a numeric vector of predicted
-probabilities or as a \code{data.frame} where each column represents the
-predicted probability for a level of the exposure. They can also be
-propensity score objects created by \code{\link[=ps_trim]{ps_trim()}}, \code{\link[=ps_refit]{ps_refit()}}, or
-\code{\link[=ps_trunc]{ps_trunc()}}
+probabilities, as a \code{data.frame} where each column represents the
+predicted probability for a level of the exposure, or as a fitted
+GLM object. They can also be propensity score objects created by
+\code{\link[=ps_trim]{ps_trim()}}, \code{\link[=ps_refit]{ps_refit()}}, or \code{\link[=ps_trunc]{ps_trunc()}}
 
 The returned weights are encapsulated in a \code{psw} object, which is a numeric
 vector with additional attributes that record the estimand, and whether the
@@ -249,15 +327,47 @@ trt_extreme <- c(0, 0, 1, 1, 0, 1, 0, 1)
 # Standard ATE weights can be extreme
 wt_extreme <- wt_ate(ps_extreme, trt_extreme)
 # Very large!
-max(wt_extreme)  
+max(wt_extreme)
 
 # ATO weights are bounded
 wt_extreme_atm <- wt_ato(ps_extreme, trt_extreme)
 # Much more reasonable
-max(wt_extreme_atm)  
+max(wt_extreme_atm)
 # but they target a different population
 estimand(wt_extreme_atm) # "ato"
 
+## Working with Data Frames
+
+# Example with custom data frame
+ps_df <- data.frame(
+  control = c(0.9, 0.7, 0.3, 0.1),
+  treated = c(0.1, 0.3, 0.7, 0.9)
+)
+exposure <- c(0, 0, 1, 1)
+
+# Uses second column by default (treated probabilities)
+wt_ate(ps_df, exposure)
+
+# Explicitly specify column by name
+wt_ate(ps_df, exposure, .propensity_col = "treated")
+
+# Or by position
+wt_ate(ps_df, exposure, .propensity_col = 2)
+
+## Working with GLM Objects
+
+# Fit a propensity score model
+set.seed(123)
+n <- 100
+x1 <- rnorm(n)
+x2 <- rnorm(n)
+treatment <- rbinom(n, 1, plogis(0.5 * x1 + 0.3 * x2))
+
+ps_model <- glm(treatment ~ x1 + x2, family = binomial)
+
+# Use GLM directly for weight calculation
+weights_from_glm <- wt_ate(ps_model, treatment)
+
 }
 \references{
 For detailed guidance on causal inference in R, see \href{https://www.r-causal.org/}{\emph{Causal Inference in R}}
diff --git a/propensity_0.0.0.9000.tar.gz b/propensity_0.0.0.9000.tar.gz
new file mode 100644
index 0000000..ae25e97
Binary files /dev/null and b/propensity_0.0.0.9000.tar.gz differ
diff --git a/tests/testthat/test-weights.R b/tests/testthat/test-weights.R
index 7756928..65783b3 100644
--- a/tests/testthat/test-weights.R
+++ b/tests/testthat/test-weights.R
@@ -305,7 +305,7 @@ test_that("wt_atu.ps_trim triggers refit check, sets 'atu; trimmed'", {
   )
   expect_s3_class(w_atu_unfit, "psw")
   expect_match(estimand(w_atu_unfit), "atu; trimmed")
-  expect_true(attr(w_atu_unfit, "trimmed"))
+  expect_true(is_ps_trimmed(w_atu_unfit))
   # ps_trim_meta copied
   expect_identical(
     attr(w_atu_unfit, "ps_trim_meta"),
@@ -324,7 +324,7 @@ test_that("wt_atu.ps_trim triggers refit check, sets 'atu; trimmed'", {
   )
   expect_s3_class(w_atu_fit, "psw")
   expect_match(estimand(w_atu_fit), "atu; trimmed")
-  expect_true(attr(w_atu_fit, "trimmed"))
+  expect_true(is_ps_trimmed(w_atu_fit))
   # confirm ps_trim_meta matches
   expect_identical(
     attr(w_atu_fit, "ps_trim_meta"),
@@ -360,7 +360,7 @@ test_that("wt_atm.ps_trim triggers refit check, sets 'atm; trimmed'", {
   )
   expect_s3_class(w_atm_unfit, "psw")
   expect_match(estimand(w_atm_unfit), "atm; trimmed")
-  expect_true(attr(w_atm_unfit, "trimmed"))
+  expect_true(is_ps_trimmed(w_atm_unfit))
 
   # Refit => no warning
   refit_obj <- ps_refit(trimmed_obj, model = fit)
@@ -374,7 +374,7 @@ test_that("wt_atm.ps_trim triggers refit check, sets 'atm; trimmed'", {
   )
   expect_s3_class(w_atm_fit, "psw")
   expect_match(estimand(w_atm_fit), "atm; trimmed")
-  expect_true(attr(w_atm_fit, "trimmed"))
+  expect_true(is_ps_trimmed(w_atm_fit))
 })
 
 test_that("wt_ato.ps_trim triggers refit check, sets 'ato; trimmed'", {
@@ -405,7 +405,7 @@ test_that("wt_ato.ps_trim triggers refit check, sets 'ato; trimmed'", {
   )
   expect_s3_class(w_ato_unfit, "psw")
   expect_match(estimand(w_ato_unfit), "ato; trimmed")
-  expect_true(attr(w_ato_unfit, "trimmed"))
+  expect_true(is_ps_trimmed(w_ato_unfit))
 
   # Refit => no warning
   refit_obj <- ps_refit(trimmed_obj, model = fit)
@@ -419,7 +419,7 @@ test_that("wt_ato.ps_trim triggers refit check, sets 'ato; trimmed'", {
   )
   expect_s3_class(w_ato_fit, "psw")
   expect_match(estimand(w_ato_fit), "ato; trimmed")
-  expect_true(attr(w_ato_fit, "trimmed"))
+  expect_true(is_ps_trimmed(w_ato_fit))
 })
 
 # Entropy weight tests
@@ -567,7 +567,7 @@ test_that("wt_entropy works with ps_trim objects", {
 
   expect_s3_class(weights, "psw")
   expect_equal(estimand(weights), "entropy; trimmed")
-  expect_true(attr(weights, "trimmed"))
+  expect_true(is_ps_trimmed(weights))
 })
 
 test_that("wt_entropy works with ps_trunc objects", {
@@ -582,7 +582,7 @@ test_that("wt_entropy works with ps_trunc objects", {
 
   expect_s3_class(weights, "psw")
   expect_equal(estimand(weights), "entropy; truncated")
-  expect_true(attr(weights, "truncated"))
+  expect_true(is_ps_truncated(weights))
 })
 
 test_that("entropy weights error on unsupported exposure types", {
@@ -595,11 +595,12 @@ test_that("entropy weights error on unsupported exposure types", {
     class = "propensity_wt_not_supported_error"
   )
 
+  # Now that continuous is not even an option for entropy,
+  # the function will error during auto-detection if given continuous data
   expect_error(
     wt_entropy(
       rnorm(10),
-      .exposure = rnorm(10),
-      exposure_type = "continuous"
+      .exposure = rnorm(10)
     ),
     class = "propensity_wt_not_supported_error"
   )
@@ -662,14 +663,11 @@ test_that("entropy weights give same treatment effect estimates as PSweight", {
   skip_if_not_installed("PSweight")
   skip_on_cran()
 
-  # Use PSweight's example data
-  data("psdata_cl", package = "PSweight")
-
   # Calculate treatment effect using PSweight
   ps_result <- PSweight::PSweight(
     ps.formula = trt ~ cov1 + cov2 + cov3 + cov4 + cov5 + cov6,
     yname = "Y",
-    data = psdata_cl,
+    data = PSweight::psdata_cl,
     weight = "entropy"
   )
 
@@ -678,25 +676,25 @@ test_that("entropy weights give same treatment effect estimates as PSweight", {
   # Calculate using our implementation
   ps_fit <- glm(
     trt ~ cov1 + cov2 + cov3 + cov4 + cov5 + cov6,
-    data = psdata_cl,
+    data = PSweight::psdata_cl,
     family = binomial
   )
   ps_scores <- fitted(ps_fit)
 
   our_weights <- wt_entropy(
     ps_scores,
-    .exposure = psdata_cl$trt,
+    .exposure = PSweight::psdata_cl$trt,
     exposure_type = "binary"
   )
 
   # Calculate weighted means
   mu1 <- weighted.mean(
-    psdata_cl$Y[psdata_cl$trt == 1],
-    our_weights[psdata_cl$trt == 1]
+    PSweight::psdata_cl$Y[PSweight::psdata_cl$trt == 1],
+    our_weights[PSweight::psdata_cl$trt == 1]
   )
   mu0 <- weighted.mean(
-    psdata_cl$Y[psdata_cl$trt == 0],
-    our_weights[psdata_cl$trt == 0]
+    PSweight::psdata_cl$Y[PSweight::psdata_cl$trt == 0],
+    our_weights[PSweight::psdata_cl$trt == 0]
   )
   our_ate <- mu1 - mu0
 
@@ -747,3 +745,1410 @@ test_that("entropy weight calculation matches manual calculation", {
 
   expect_equal(as.numeric(weights), expected, tolerance = 1e-10)
 })
+
+# Tests for data.frame methods
+test_that("wt_ate works with data frames", {
+  # Create test data frame
+  ps_df <- data.frame(
+    control = c(0.9, 0.7, 0.3, 0.1),
+    treated = c(0.1, 0.3, 0.7, 0.9)
+  )
+  exposure <- c(0, 0, 1, 1)
+
+  # Test default behavior (uses second column)
+  weights <- wt_ate(ps_df, exposure, exposure_type = "binary")
+  expected <- wt_ate(ps_df$treated, exposure, exposure_type = "binary")
+  expect_equal(weights, expected)
+
+  # Test explicit column selection by name (quoted)
+  weights_quoted <- wt_ate(
+    ps_df,
+    exposure,
+    .propensity_col = "treated",
+    exposure_type = "binary"
+  )
+  expect_equal(weights_quoted, expected)
+
+  # Test unquoted column selection
+  weights_unquoted <- wt_ate(
+    ps_df,
+    exposure,
+    .propensity_col = treated,
+    exposure_type = "binary"
+  )
+  expect_equal(weights_unquoted, expected)
+
+  # Test column selection by index
+  weights_idx <- wt_ate(
+    ps_df,
+    exposure,
+    .propensity_col = 2,
+    exposure_type = "binary"
+  )
+  expect_equal(weights_idx, expected)
+
+  # Test single column data frame
+  ps_single <- data.frame(prob = c(0.1, 0.3, 0.7, 0.9))
+  weights_single <- wt_ate(ps_single, exposure, exposure_type = "binary")
+  expected_single <- wt_ate(ps_single$prob, exposure, exposure_type = "binary")
+  expect_equal(weights_single, expected_single)
+
+  # Test error with empty data frame
+  expect_error(
+    wt_ate(data.frame(), exposure),
+    class = "propensity_df_ncol_error"
+  )
+
+  # Test error with invalid column selection
+  expect_error(
+    wt_ate(ps_df, exposure, .propensity_col = "nonexistent"),
+    class = "propensity_df_column_error"
+  )
+})
+
+test_that("all wt_* functions work with data frames", {
+  ps_df <- data.frame(
+    control = c(0.9, 0.7, 0.3, 0.1),
+    treated = c(0.1, 0.3, 0.7, 0.9)
+  )
+  exposure <- c(0, 0, 1, 1)
+
+  # Test ATT
+  weights_att <- wt_att(ps_df, exposure, exposure_type = "binary")
+  expected_att <- wt_att(ps_df$treated, exposure, exposure_type = "binary")
+  expect_equal(weights_att, expected_att)
+
+  # Test ATU
+  weights_atu <- wt_atu(ps_df, exposure, exposure_type = "binary")
+  expected_atu <- wt_atu(ps_df$treated, exposure, exposure_type = "binary")
+  expect_equal(weights_atu, expected_atu)
+
+  # Test ATM
+  weights_atm <- wt_atm(ps_df, exposure, exposure_type = "binary")
+  expected_atm <- wt_atm(ps_df$treated, exposure, exposure_type = "binary")
+  expect_equal(weights_atm, expected_atm)
+
+  # Test ATO
+  weights_ato <- wt_ato(ps_df, exposure, exposure_type = "binary")
+  expected_ato <- wt_ato(ps_df$treated, exposure, exposure_type = "binary")
+  expect_equal(weights_ato, expected_ato)
+
+  # Test Entropy
+  weights_entropy <- wt_entropy(ps_df, exposure, exposure_type = "binary")
+  expected_entropy <- wt_entropy(
+    ps_df$treated,
+    exposure,
+    exposure_type = "binary"
+  )
+  expect_equal(weights_entropy, expected_entropy)
+})
+
+test_that("wt_* functions work with parsnip output", {
+  skip_if_not_installed("parsnip")
+  skip_on_cran()
+
+  # Simulate data
+  set.seed(123)
+  n <- 100
+  x1 <- rnorm(n)
+  x2 <- rnorm(n)
+  ps_true <- plogis(0.5 * x1 + 0.3 * x2)
+  treatment_numeric <- rbinom(n, 1, ps_true)
+  treatment_factor <- factor(treatment_numeric, levels = c("0", "1"))
+
+  df <- data.frame(
+    treatment = treatment_factor,
+    x1 = x1,
+    x2 = x2
+  )
+
+  # Fit model with parsnip
+  ps_spec <- parsnip::logistic_reg()
+  ps_spec <- parsnip::set_engine(ps_spec, "glm")
+  ps_model <- parsnip::fit(ps_spec, treatment ~ x1 + x2, data = df)
+
+  # Get predictions
+  ps_preds <- predict(ps_model, df, type = "prob")
+
+  # Test that it works with default (second column)
+  weights_ate <- wt_ate(ps_preds, treatment_numeric, exposure_type = "binary")
+  expect_s3_class(weights_ate, "psw")
+  expect_equal(estimand(weights_ate), "ate")
+
+  # Test explicit column selection with parsnip column names
+  weights_ate2 <- wt_ate(
+    ps_preds,
+    treatment_numeric,
+    .propensity_col = ".pred_1",
+    exposure_type = "binary"
+  )
+  expect_equal(weights_ate, weights_ate2)
+
+  # Test with all estimands
+  expect_s3_class(
+    wt_att(ps_preds, treatment_numeric, exposure_type = "binary"),
+    "psw"
+  )
+  expect_s3_class(
+    wt_atu(ps_preds, treatment_numeric, exposure_type = "binary"),
+    "psw"
+  )
+  expect_s3_class(
+    wt_atm(ps_preds, treatment_numeric, exposure_type = "binary"),
+    "psw"
+  )
+  expect_s3_class(
+    wt_ato(ps_preds, treatment_numeric, exposure_type = "binary"),
+    "psw"
+  )
+  expect_s3_class(
+    wt_entropy(ps_preds, treatment_numeric, exposure_type = "binary"),
+    "psw"
+  )
+})
+
+test_that("wt_* functions work with GLM objects", {
+  # Simulate data
+  set.seed(123)
+  n <- 100
+  x1 <- rnorm(n)
+  x2 <- rnorm(n)
+  ps_true <- plogis(0.5 * x1 + 0.3 * x2)
+  treatment <- rbinom(n, 1, ps_true)
+
+  # Fit GLM model
+  ps_model <- glm(treatment ~ x1 + x2, family = binomial)
+  ps_fitted <- fitted(ps_model)
+
+  # Test ATE
+  weights_ate_glm <- wt_ate(ps_model, treatment, exposure_type = "binary")
+  weights_ate_numeric <- wt_ate(ps_fitted, treatment, exposure_type = "binary")
+  expect_equal(weights_ate_glm, weights_ate_numeric)
+  expect_s3_class(weights_ate_glm, "psw")
+  expect_equal(estimand(weights_ate_glm), "ate")
+
+  # Test ATT
+  weights_att_glm <- wt_att(ps_model, treatment, exposure_type = "binary")
+  weights_att_numeric <- wt_att(ps_fitted, treatment, exposure_type = "binary")
+  expect_equal(weights_att_glm, weights_att_numeric)
+  expect_s3_class(weights_att_glm, "psw")
+  expect_equal(estimand(weights_att_glm), "att")
+
+  # Test ATU
+  weights_atu_glm <- wt_atu(ps_model, treatment, exposure_type = "binary")
+  weights_atu_numeric <- wt_atu(ps_fitted, treatment, exposure_type = "binary")
+  expect_equal(weights_atu_glm, weights_atu_numeric)
+  expect_s3_class(weights_atu_glm, "psw")
+  expect_equal(estimand(weights_atu_glm), "atu")
+
+  # Test ATM
+  weights_atm_glm <- wt_atm(ps_model, treatment, exposure_type = "binary")
+  weights_atm_numeric <- wt_atm(ps_fitted, treatment, exposure_type = "binary")
+  expect_equal(weights_atm_glm, weights_atm_numeric)
+  expect_s3_class(weights_atm_glm, "psw")
+  expect_equal(estimand(weights_atm_glm), "atm")
+
+  # Test ATO
+  weights_ato_glm <- wt_ato(ps_model, treatment, exposure_type = "binary")
+  weights_ato_numeric <- wt_ato(ps_fitted, treatment, exposure_type = "binary")
+  expect_equal(weights_ato_glm, weights_ato_numeric)
+  expect_s3_class(weights_ato_glm, "psw")
+  expect_equal(estimand(weights_ato_glm), "ato")
+
+  # Test Entropy
+  weights_entropy_glm <- wt_entropy(
+    ps_model,
+    treatment,
+    exposure_type = "binary"
+  )
+  weights_entropy_numeric <- wt_entropy(
+    ps_fitted,
+    treatment,
+    exposure_type = "binary"
+  )
+  expect_equal(weights_entropy_glm, weights_entropy_numeric)
+  expect_s3_class(weights_entropy_glm, "psw")
+  expect_equal(estimand(weights_entropy_glm), "entropy")
+})
+
+test_that("GLM methods handle continuous exposures", {
+  # Simulate continuous exposure data
+  set.seed(456)
+  n <- 100
+  x1 <- rnorm(n)
+  x2 <- rnorm(n)
+  exposure <- 2 + 0.5 * x1 + 0.3 * x2 + rnorm(n)
+
+  # Fit linear model
+  exposure_model <- glm(exposure ~ x1 + x2, family = gaussian)
+
+  # Test ATE with continuous exposure
+  weights_ate <- wt_ate(
+    exposure_model,
+    exposure,
+    exposure_type = "continuous",
+    stabilize = TRUE
+  )
+  expect_s3_class(weights_ate, "psw")
+  expect_equal(estimand(weights_ate), "ate")
+  expect_true(is_stabilized(weights_ate))
+
+  # Check that weights are reasonable
+  expect_true(all(is.finite(weights_ate)))
+  expect_true(all(weights_ate > 0))
+})
+
+test_that("GLM methods error on non-GLM objects", {
+  # Try with a non-GLM object
+  expect_error(
+    wt_ate("not a glm", c(0, 1, 0, 1)),
+    class = "propensity_method_error"
+  )
+
+  expect_error(
+    wt_att(list(a = 1, b = 2), c(0, 1, 0, 1)),
+    class = "propensity_method_error"
+  )
+})
+
+# Edge case tests for all wt_* functions ----
+
+test_that("wt_* functions handle edge case propensity scores", {
+  # Edge case: propensity scores at boundaries
+  ps_boundary <- c(1e-10, 0.001, 0.5, 0.999, 1 - 1e-10)
+  exposure_boundary <- c(0, 0, 1, 1, 1)
+
+  # Test all functions with boundary values
+  for (fn in list(wt_ate, wt_att, wt_atu, wt_atm, wt_ato, wt_entropy)) {
+    weights <- fn(ps_boundary, exposure_boundary, exposure_type = "binary")
+    expect_s3_class(weights, "psw")
+    expect_true(all(is.finite(weights)))
+    expect_true(all(weights > 0))
+  }
+
+  # Edge case: all propensity scores identical
+  ps_constant <- rep(0.5, 5)
+  exposure_mixed <- c(0, 0, 1, 1, 0)
+
+  for (fn in list(wt_ate, wt_att, wt_atu, wt_atm, wt_ato, wt_entropy)) {
+    weights <- fn(ps_constant, exposure_mixed, exposure_type = "binary")
+    expect_s3_class(weights, "psw")
+    # For constant propensity scores, weights should be constant within treatment groups
+    expect_equal(length(unique(weights[exposure_mixed == 0])), 1)
+    expect_equal(length(unique(weights[exposure_mixed == 1])), 1)
+  }
+
+  # Edge case: single observation
+  ps_single <- 0.3
+  exposure_single <- 1
+
+  # Single observation needs explicit .treated/.untreated
+  for (fn in list(wt_ate, wt_att, wt_atu, wt_atm, wt_ato, wt_entropy)) {
+    weights <- fn(
+      ps_single,
+      exposure_single,
+      exposure_type = "binary",
+      .treated = 1
+    )
+    expect_s3_class(weights, "psw")
+    expect_length(weights, 1)
+    expect_true(is.finite(weights))
+  }
+
+  # Edge case: all treated or all control
+  ps_various <- c(0.2, 0.4, 0.6, 0.8)
+
+  for (fn in list(wt_ate, wt_att, wt_atu, wt_atm, wt_ato, wt_entropy)) {
+    # All treated - need to specify .treated since there's only one level
+    weights_all_1 <- fn(
+      ps_various,
+      rep(1, 4),
+      exposure_type = "binary",
+      .treated = 1
+    )
+    expect_s3_class(weights_all_1, "psw")
+    expect_true(all(is.finite(weights_all_1)))
+
+    # All control - need to specify .untreated since there's only one level
+    weights_all_0 <- fn(
+      ps_various,
+      rep(0, 4),
+      exposure_type = "binary",
+      .untreated = 0
+    )
+    expect_s3_class(weights_all_0, "psw")
+    expect_true(all(is.finite(weights_all_0)))
+  }
+})
+
+test_that("wt_* functions handle extreme weight scenarios", {
+  # Create scenario that produces very large weights
+  ps_extreme <- c(0.001, 0.999, 0.5, 0.5)
+  exposure_extreme <- c(1, 0, 0, 1) # Misaligned with propensity
+
+  # ATE should produce extreme weights
+  weights_ate <- wt_ate(ps_extreme, exposure_extreme, exposure_type = "binary")
+  expect_true(max(weights_ate) > 100) # Very large weight expected
+
+  # ATO and ATM should be bounded
+  weights_ato <- wt_ato(ps_extreme, exposure_extreme, exposure_type = "binary")
+  weights_atm <- wt_atm(ps_extreme, exposure_extreme, exposure_type = "binary")
+  expect_true(all(weights_ato <= 1)) # ATO weights bounded by 1
+  expect_true(all(weights_atm <= 2)) # ATM weights bounded
+})
+
+# Additional error handling tests ----
+
+test_that("wt_* functions error appropriately on invalid inputs", {
+  # Invalid propensity scores (outside [0,1])
+  expect_error(
+    wt_ate(c(-0.1, 0.5, 1.1), c(0, 1, 0), exposure_type = "binary"),
+    class = "propensity_range_error"
+  )
+
+  expect_error(
+    wt_att(c(0, 0.5, 1), c(0, 1, 0), exposure_type = "binary"),
+    class = "propensity_range_error"
+  )
+
+  # Length mismatch should error
+  expect_error(
+    wt_ate(c(0.1, 0.5), c(0, 1, 0), exposure_type = "binary"),
+    class = "propensity_length_error"
+  )
+
+  # Invalid exposure type
+  expect_error(
+    wt_ate(c(0.1, 0.5), c(0, 1), exposure_type = "invalid"),
+    "must be one of"
+  )
+
+  # Non-numeric propensity scores
+  expect_error(
+    wt_ate(c("a", "b"), c(0, 1)),
+    class = "propensity_method_error"
+  )
+
+  # Empty vectors
+  expect_error(
+    wt_ate(numeric(0), numeric(0))
+  )
+
+  # Categorical exposure (not supported for most estimands)
+  expect_error(
+    wt_att(c(0.3, 0.3, 0.4), c(1, 2, 3), exposure_type = "categorical"),
+    class = "propensity_wt_not_supported_error"
+  )
+})
+
+test_that("data frame methods error appropriately", {
+  # Empty data frame
+  expect_error(
+    wt_ate(data.frame(), c(0, 1)),
+    class = "propensity_df_ncol_error"
+  )
+
+  # Non-existent column
+  df <- data.frame(a = c(0.1, 0.9), b = c(0.9, 0.1))
+  expect_error(
+    wt_ate(df, c(0, 1), .propensity_col = "nonexistent"),
+    class = "propensity_df_column_error"
+  )
+
+  # Column index out of bounds
+  expect_error(
+    wt_ate(df, c(0, 1), .propensity_col = 5),
+    class = "propensity_df_column_error"
+  )
+
+  # Non-numeric column
+  df_char <- data.frame(a = c("high", "low"), b = c(0.9, 0.1))
+  suppressWarnings({
+    expect_error(
+      wt_ate(df_char, c(0, 1), .propensity_col = 1)
+    )
+  })
+
+  # Column with invalid values
+  df_invalid <- data.frame(a = c(0.5, 1.5), b = c(0.9, 0.1))
+  expect_error(
+    wt_ate(df_invalid, c(0, 1), .propensity_col = 1),
+    class = "propensity_range_error"
+  )
+})
+
+test_that("GLM methods error appropriately", {
+  # Non-GLM object
+  expect_error(
+    wt_ate(lm(mpg ~ wt, data = mtcars), rep(0:1, 16)),
+    class = "propensity_method_error"
+  )
+
+  # GLM with wrong dimensions
+  small_glm <- glm(c(0, 1) ~ c(1, 2), family = binomial)
+  expect_error(
+    wt_ate(small_glm, c(0, 1, 0, 1)), # Mismatch in length
+    class = "propensity_length_error"
+  )
+})
+
+# Default method dispatch tests ----
+
+test_that("default methods provide informative errors", {
+  # Custom class that doesn't have a method
+  custom_obj <- structure(list(x = 1:5), class = "my_custom_class")
+
+  expect_error(
+    wt_ate(custom_obj, c(0, 1, 0, 1, 0)),
+    class = "propensity_method_error"
+  )
+
+  expect_error(
+    wt_att(custom_obj, c(0, 1, 0, 1, 0)),
+    class = "propensity_method_error"
+  )
+
+  expect_error(
+    wt_atu(custom_obj, c(0, 1, 0, 1, 0)),
+    class = "propensity_method_error"
+  )
+
+  expect_error(
+    wt_atm(custom_obj, c(0, 1, 0, 1, 0)),
+    class = "propensity_method_error"
+  )
+
+  expect_error(
+    wt_ato(custom_obj, c(0, 1, 0, 1, 0)),
+    class = "propensity_method_error"
+  )
+
+  expect_error(
+    wt_entropy(custom_obj, c(0, 1, 0, 1, 0)),
+    class = "propensity_method_error"
+  )
+})
+
+# Data frame method tests with various column types ----
+
+test_that("data frame methods handle various column configurations", {
+  # Data frame with many columns
+  ps_multi <- data.frame(
+    col1 = runif(10, 0.1, 0.9),
+    col2 = runif(10, 0.1, 0.9),
+    col3 = runif(10, 0.1, 0.9),
+    col4 = runif(10, 0.1, 0.9),
+    col5 = runif(10, 0.1, 0.9)
+  )
+  exposure <- rbinom(10, 1, 0.5)
+
+  # Test column selection by position
+  for (i in 1:5) {
+    weights <- suppressWarnings(wt_ate(
+      ps_multi,
+      exposure,
+      .propensity_col = i,
+      exposure_type = "binary"
+    ))
+    expected <- wt_ate(ps_multi[[i]], exposure, exposure_type = "binary")
+    expect_equal(weights, expected)
+  }
+
+  # Test with tibble
+  if (requireNamespace("tibble", quietly = TRUE)) {
+    ps_tibble <- tibble::as_tibble(ps_multi)
+    weights_tibble <- wt_ate(
+      ps_tibble,
+      exposure,
+      .propensity_col = 3,
+      exposure_type = "binary"
+    )
+    weights_df <- wt_ate(
+      ps_multi,
+      exposure,
+      .propensity_col = 3,
+      exposure_type = "binary"
+    )
+    expect_equal(weights_tibble, weights_df)
+  }
+
+  # Test with column names containing spaces
+  ps_spaces <- data.frame(
+    `control probability` = runif(5, 0.1, 0.9),
+    `treatment probability` = runif(5, 0.1, 0.9),
+    check.names = FALSE
+  )
+  exposure_small <- c(0, 0, 1, 1, 0)
+
+  weights_spaces <- wt_ate(
+    ps_spaces,
+    exposure_small,
+    .propensity_col = "treatment probability",
+    exposure_type = "binary"
+  )
+  expect_s3_class(weights_spaces, "psw")
+
+  # Test tidyselect helpers
+  weights_last <- wt_ate(
+    ps_multi,
+    exposure,
+    .propensity_col = tidyselect::last_col(),
+    exposure_type = "binary"
+  )
+  expected_last <- wt_ate(ps_multi$col5, exposure, exposure_type = "binary")
+  expect_equal(weights_last, expected_last)
+})
+
+# GLM method tests with different families ----
+
+test_that("GLM methods handle non-binomial families appropriately", {
+  # Gaussian family (for continuous exposures)
+  set.seed(123)
+  n <- 50
+  x <- rnorm(n)
+  exposure_cont <- 2 + 0.5 * x + rnorm(n)
+
+  glm_gaussian <- glm(exposure_cont ~ x, family = gaussian)
+
+  # Should work for ATE with continuous exposure
+  weights_gaussian <- wt_ate(
+    glm_gaussian,
+    exposure_cont,
+    exposure_type = "continuous",
+    stabilize = TRUE
+  )
+  expect_s3_class(weights_gaussian, "psw")
+  expect_true(all(is.finite(weights_gaussian)))
+
+  # Should error for estimands that don't support continuous
+  # ATT doesn't accept continuous as a valid exposure type
+  expect_error(
+    wt_att(glm_gaussian, exposure_cont, exposure_type = "continuous"),
+    "must be one of"
+  )
+
+  # Poisson family (should extract fitted values)
+  # For non-binomial GLMs, the fitted values might not be valid propensity scores
+  # Skip this test as it's not a valid use case
+})
+
+# Attribute preservation tests ----
+
+test_that("attributes are preserved across all weight methods", {
+  # Create trimmed propensity scores with attributes
+  ps <- runif(20, 0.05, 0.95)
+  exposure <- rbinom(20, 1, ps)
+
+  ps_trimmed <- ps_trim(
+    ps,
+    .exposure = exposure,
+    method = "ps",
+    lower = 0.1,
+    upper = 0.9
+  )
+
+  # Test that all weight functions preserve trim attributes
+  for (fn_name in c(
+    "wt_ate",
+    "wt_att",
+    "wt_atu",
+    "wt_atm",
+    "wt_ato",
+    "wt_entropy"
+  )) {
+    fn <- get(fn_name)
+
+    # Suppress refit warning
+    suppressWarnings({
+      weights <- fn(ps_trimmed, exposure, exposure_type = "binary")
+    })
+
+    expect_true(is_ps_trimmed(weights))
+    expect_equal(
+      attr(weights, "ps_trim_meta"),
+      attr(ps_trimmed, "ps_trim_meta")
+    )
+    expect_match(estimand(weights), "; trimmed$")
+  }
+
+  # Test with truncated propensity scores
+  ps_truncated <- ps_trunc(ps, method = "ps", lower = 0.1, upper = 0.9)
+
+  for (fn_name in c(
+    "wt_ate",
+    "wt_att",
+    "wt_atu",
+    "wt_atm",
+    "wt_ato",
+    "wt_entropy"
+  )) {
+    fn <- get(fn_name)
+    weights <- fn(ps_truncated, exposure, exposure_type = "binary")
+
+    expect_true(is_ps_truncated(weights))
+    expect_equal(
+      attr(weights, "ps_trunc_meta"),
+      attr(ps_truncated, "ps_trunc_meta")
+    )
+    expect_match(estimand(weights), "; truncated$")
+  }
+})
+
+test_that("stabilization attributes are set correctly", {
+  ps <- runif(20, 0.1, 0.9)
+  exposure <- rbinom(20, 1, ps)
+
+  # Test ATE with stabilization
+  weights_unstab <- wt_ate(
+    ps,
+    exposure,
+    stabilize = FALSE,
+    exposure_type = "binary"
+  )
+  weights_stab <- wt_ate(
+    ps,
+    exposure,
+    stabilize = TRUE,
+    exposure_type = "binary"
+  )
+
+  expect_false(is_stabilized(weights_unstab))
+  expect_true(is_stabilized(weights_stab))
+
+  # Test with custom stabilization score
+  weights_custom <- wt_ate(
+    ps,
+    exposure,
+    stabilize = TRUE,
+    stabilization_score = 0.4,
+    exposure_type = "binary"
+  )
+  expect_true(is_stabilized(weights_custom))
+})
+
+# Stabilization tests across methods ----
+
+test_that("stabilization works correctly for all applicable methods", {
+  set.seed(456)
+  ps <- runif(30, 0.2, 0.8)
+  exposure <- rbinom(30, 1, ps)
+
+  # ATE supports stabilization
+  weights_ate_unstab <- wt_ate(
+    ps,
+    exposure,
+    stabilize = FALSE,
+    exposure_type = "binary"
+  )
+  weights_ate_stab <- wt_ate(
+    ps,
+    exposure,
+    stabilize = TRUE,
+    exposure_type = "binary"
+  )
+
+  # Check that stabilization was applied
+  expect_true(is_stabilized(weights_ate_stab))
+})
+
+# NA handling tests ----
+
+test_that("all methods handle NAs appropriately", {
+  # Propensity scores with NAs
+  ps_na <- c(0.2, NA, 0.5, 0.8, NA)
+  exposure_na <- c(0, 1, 1, 0, 1)
+
+  # All methods should handle NAs by producing NA weights
+  for (fn in list(wt_ate, wt_att, wt_atu, wt_atm, wt_ato, wt_entropy)) {
+    weights <- fn(ps_na, exposure_na, exposure_type = "binary")
+    expect_s3_class(weights, "psw")
+    expect_true(any(is.na(weights)))
+  }
+
+  # Exposure with NAs
+  ps_good <- c(0.2, 0.3, 0.5, 0.8, 0.7)
+  exposure_with_na <- c(0, NA, 1, 0, 1)
+
+  for (fn in list(wt_ate, wt_att, wt_atu, wt_atm, wt_ato, wt_entropy)) {
+    weights <- fn(ps_good, exposure_with_na, exposure_type = "binary")
+    expect_s3_class(weights, "psw")
+    expect_true(any(is.na(weights)))
+  }
+
+  # Data frame with NAs
+  df_na <- data.frame(
+    col1 = c(0.1, NA, 0.5),
+    col2 = c(0.9, 0.5, NA)
+  )
+
+  # Data frame with NAs produces NA weights
+  weights_df_na <- wt_ate(df_na, c(0, 1, 0), exposure_type = "binary")
+  expect_s3_class(weights_df_na, "psw")
+  expect_true(any(is.na(weights_df_na)))
+
+  # GLM predictions with NAs
+  set.seed(789)
+  n <- 20
+  x <- c(rnorm(18), NA, NA)
+  y <- c(rbinom(18, 1, plogis(0.5 * x[1:18])), 0, 1)
+
+  # GLM will handle NAs in predictors
+  glm_na <- glm(y ~ x, family = binomial)
+
+  # GLM with NA predictions will have shorter output
+  # The NA observations are dropped during model fitting
+  suppressWarnings({
+    # GLM drops NA observations, so output is shorter
+    expect_error(
+      wt_ate(glm_na, y, exposure_type = "binary"),
+      "must have the same length"
+    )
+  })
+})
+
+# Integration tests ----
+
+test_that("weight functions integrate correctly with ps_trim and ps_trunc", {
+  set.seed(999)
+  n <- 100
+  ps <- runif(n, 0.05, 0.95)
+  exposure <- rbinom(n, 1, ps)
+
+  # Create a model for refitting
+  x <- rnorm(n)
+  model <- glm(exposure ~ x + I(x^2), family = binomial)
+
+  # Trim and refit
+  ps_trim_obj <- ps_trim(
+    ps,
+    .exposure = exposure,
+    method = "ps",
+    lower = 0.1,
+    upper = 0.9
+  )
+  ps_refit_obj <- ps_refit(ps_trim_obj, model = model)
+
+  # Should not warn after refit
+  suppressMessages({
+    weights_refit <- wt_ate(ps_refit_obj, exposure, exposure_type = "binary")
+  })
+  expect_true(is_ps_trimmed(weights_refit))
+  expect_true(is_refit(weights_refit))
+
+  # Truncate
+  ps_trunc_obj <- ps_trunc(ps, method = "ps", lower = 0.1, upper = 0.9)
+
+  # Should not warn for truncation
+  suppressMessages({
+    weights_trunc <- wt_ate(ps_trunc_obj, exposure, exposure_type = "binary")
+  })
+  expect_true(is_ps_truncated(weights_trunc))
+  expect_false(is_ps_trimmed(weights_trunc))
+})
+
+test_that("data frame and GLM methods produce consistent results", {
+  set.seed(111)
+  n <- 50
+  x1 <- rnorm(n)
+  x2 <- rnorm(n)
+  true_ps <- plogis(0.5 * x1 + 0.3 * x2)
+  treatment <- rbinom(n, 1, true_ps)
+
+  # Fit GLM
+  model <- glm(treatment ~ x1 + x2, family = binomial)
+  ps_fitted <- fitted(model)
+
+  # Create data frame
+  ps_df <- data.frame(
+    control_prob = 1 - ps_fitted,
+    treatment_prob = ps_fitted
+  )
+
+  # All methods should give same results
+  for (fn_name in c(
+    "wt_ate",
+    "wt_att",
+    "wt_atu",
+    "wt_atm",
+    "wt_ato",
+    "wt_entropy"
+  )) {
+    fn <- get(fn_name)
+
+    # From GLM
+    weights_glm <- fn(model, treatment, exposure_type = "binary")
+
+    # From fitted values
+    weights_numeric <- fn(ps_fitted, treatment, exposure_type = "binary")
+
+    # From data frame (using treatment column)
+    weights_df <- fn(
+      ps_df,
+      treatment,
+      .propensity_col = 2,
+      exposure_type = "binary"
+    )
+
+    expect_equal(weights_glm, weights_numeric)
+    expect_equal(weights_glm, weights_df)
+  }
+})
+
+# Additional mathematical property tests ----
+
+test_that("weight functions satisfy mathematical properties", {
+  set.seed(222)
+  n <- 100
+  ps <- runif(n, 0.1, 0.9)
+  treatment <- rbinom(n, 1, ps)
+
+  # ATT weights: treated units should have weight 1
+  weights_att <- wt_att(ps, treatment, exposure_type = "binary")
+  expect_true(all(weights_att[treatment == 1] == 1))
+
+  # ATU weights: control units should have weight 1
+  weights_atu <- wt_atu(ps, treatment, exposure_type = "binary")
+  expect_true(all(weights_atu[treatment == 0] == 1))
+
+  # ATO weights: A * (1-e) + (1-A) * e
+  weights_ato <- wt_ato(ps, treatment, exposure_type = "binary")
+  # ATO weights are bounded by 1
+  expect_true(all(weights_ato <= 1 + 1e-10))
+
+  # ATM weights: symmetric for e and 1-e
+  # Create symmetric propensity scores
+  ps_sym <- c(0.3, 0.7, 0.4, 0.6)
+  treatment_sym <- c(0, 1, 1, 0)
+  weights_atm <- wt_atm(ps_sym, treatment_sym, exposure_type = "binary")
+
+  # Weights for symmetric PS pairs should be equal
+  expect_equal(weights_atm[1], weights_atm[2], tolerance = 1e-10)
+  expect_equal(weights_atm[3], weights_atm[4], tolerance = 1e-10)
+})
+
+test_that("continuous exposure weights have correct properties", {
+  set.seed(333)
+  n <- 100
+  x <- rnorm(n)
+  exposure <- 2 + 0.5 * x + rnorm(n, sd = 1.5)
+
+  # Fit model
+  model <- lm(exposure ~ x)
+  mu_hat <- fitted(model)
+  sigma <- summary(model)$sigma
+
+  # Calculate weights
+  weights_cont <- wt_ate(
+    mu_hat,
+    .exposure = exposure,
+    .sigma = rep(sigma, n),
+    exposure_type = "continuous",
+    stabilize = TRUE
+  )
+
+  expect_s3_class(weights_cont, "psw")
+  expect_true(all(is.finite(weights_cont)))
+  expect_true(all(weights_cont > 0))
+
+  # Stabilized continuous weights should have mean approximately 1
+  expect_equal(mean(weights_cont), 1, tolerance = 0.2)
+})
+
+# Comparison tests with other packages ----
+
+test_that("ATE weights match WeightIt", {
+  skip_if_not_installed("WeightIt")
+  skip_on_cran()
+
+  # Simulate data
+  set.seed(123)
+  n <- 200
+  x1 <- rnorm(n)
+  x2 <- rnorm(n)
+  ps_true <- plogis(0.5 * x1 + 0.3 * x2)
+  treatment <- rbinom(n, 1, ps_true)
+
+  # Fit propensity score model
+  ps_model <- glm(treatment ~ x1 + x2, family = binomial)
+  ps_fitted <- fitted(ps_model)
+
+  # Calculate weights with our package
+  our_weights <- wt_ate(ps_fitted, treatment, exposure_type = "binary")
+
+  # Calculate weights with WeightIt
+  weightit_obj <- WeightIt::weightit(
+    treatment ~ x1 + x2,
+    data = data.frame(treatment = treatment, x1 = x1, x2 = x2),
+    method = "ps",
+    estimand = "ATE"
+  )
+  weightit_weights <- weightit_obj$weights
+
+  # Compare
+  expect_equal(as.numeric(our_weights), weightit_weights, tolerance = 1e-10)
+})
+
+test_that("ATT weights match WeightIt", {
+  skip_if_not_installed("WeightIt")
+  skip_on_cran()
+
+  # Simulate data
+  set.seed(456)
+  n <- 200
+  x1 <- rnorm(n)
+  x2 <- rnorm(n)
+  ps_true <- plogis(0.5 * x1 + 0.3 * x2)
+  treatment <- rbinom(n, 1, ps_true)
+
+  # Fit propensity score model
+  ps_model <- glm(treatment ~ x1 + x2, family = binomial)
+  ps_fitted <- fitted(ps_model)
+
+  # Calculate weights with our package
+  our_weights <- wt_att(ps_fitted, treatment, exposure_type = "binary")
+
+  # Calculate weights with WeightIt
+  weightit_obj <- WeightIt::weightit(
+    treatment ~ x1 + x2,
+    data = data.frame(treatment = treatment, x1 = x1, x2 = x2),
+    method = "ps",
+    estimand = "ATT"
+  )
+  weightit_weights <- weightit_obj$weights
+
+  # Compare
+  expect_equal(as.numeric(our_weights), weightit_weights, tolerance = 1e-10)
+})
+
+test_that("ATU/ATC weights match WeightIt", {
+  skip_if_not_installed("WeightIt")
+  skip_on_cran()
+
+  # Simulate data
+  set.seed(789)
+  n <- 200
+  x1 <- rnorm(n)
+  x2 <- rnorm(n)
+  ps_true <- plogis(0.5 * x1 + 0.3 * x2)
+  treatment <- rbinom(n, 1, ps_true)
+
+  # Fit propensity score model
+  ps_model <- glm(treatment ~ x1 + x2, family = binomial)
+  ps_fitted <- fitted(ps_model)
+
+  # Calculate weights with our package
+  our_weights <- wt_atu(ps_fitted, treatment, exposure_type = "binary")
+
+  # Calculate weights with WeightIt
+  # WeightIt uses ATC for Average Treatment Effect on Controls
+  weightit_obj <- WeightIt::weightit(
+    treatment ~ x1 + x2,
+    data = data.frame(treatment = treatment, x1 = x1, x2 = x2),
+    method = "ps",
+    estimand = "ATC"
+  )
+  weightit_weights <- weightit_obj$weights
+
+  # Compare
+  expect_equal(as.numeric(our_weights), weightit_weights, tolerance = 1e-10)
+})
+
+test_that("ATM weights match WeightIt", {
+  skip_if_not_installed("WeightIt")
+  skip_on_cran()
+
+  # Simulate data
+  set.seed(321)
+  n <- 200
+  x1 <- rnorm(n)
+  x2 <- rnorm(n)
+  ps_true <- plogis(0.5 * x1 + 0.3 * x2)
+  treatment <- rbinom(n, 1, ps_true)
+
+  # Fit propensity score model
+  ps_model <- glm(treatment ~ x1 + x2, family = binomial)
+  ps_fitted <- fitted(ps_model)
+
+  # Calculate weights with our package
+  our_weights <- wt_atm(ps_fitted, treatment, exposure_type = "binary")
+
+  # Calculate weights with WeightIt
+  # WeightIt uses ATM for matching weights
+  weightit_obj <- WeightIt::weightit(
+    treatment ~ x1 + x2,
+    data = data.frame(treatment = treatment, x1 = x1, x2 = x2),
+    method = "ps",
+    estimand = "ATM"
+  )
+  weightit_weights <- weightit_obj$weights
+
+  # Compare
+  expect_equal(as.numeric(our_weights), weightit_weights, tolerance = 1e-10)
+})
+
+test_that("ATO weights match WeightIt", {
+  skip_if_not_installed("WeightIt")
+  skip_on_cran()
+
+  # Simulate data
+  set.seed(654)
+  n <- 200
+  x1 <- rnorm(n)
+  x2 <- rnorm(n)
+  ps_true <- plogis(0.5 * x1 + 0.3 * x2)
+  treatment <- rbinom(n, 1, ps_true)
+
+  # Fit propensity score model
+  ps_model <- glm(treatment ~ x1 + x2, family = binomial)
+  ps_fitted <- fitted(ps_model)
+
+  # Calculate weights with our package
+  our_weights <- wt_ato(ps_fitted, treatment, exposure_type = "binary")
+
+  # Calculate weights with WeightIt
+  # WeightIt uses ATO for overlap weights
+  weightit_obj <- WeightIt::weightit(
+    treatment ~ x1 + x2,
+    data = data.frame(treatment = treatment, x1 = x1, x2 = x2),
+    method = "ps",
+    estimand = "ATO"
+  )
+  weightit_weights <- weightit_obj$weights
+
+  # Compare
+  expect_equal(as.numeric(our_weights), weightit_weights, tolerance = 1e-10)
+})
+
+# PSweight comparison tests ----
+
+test_that("ATE weights match PSweight", {
+  skip_if_not_installed("PSweight")
+  skip_on_cran()
+
+  # Fit propensity score model
+  ps_formula <- trt ~ cov1 + cov2 + cov3 + cov4 + cov5 + cov6
+  ps_model <- glm(ps_formula, data = PSweight::psdata_cl, family = binomial)
+  ps_fitted <- fitted(ps_model)
+
+  # Calculate weights with our package
+  our_weights <- wt_ate(
+    ps_fitted,
+    PSweight::psdata_cl$trt,
+    exposure_type = "binary"
+  )
+
+  # Calculate weights with PSweight
+  psw_obj <- PSweight::SumStat(
+    ps.formula = ps_formula,
+    data = PSweight::psdata_cl,
+    weight = "IPW"
+  )
+
+  # Extract IPW weights and un-normalize them
+  # PSweight normalizes weights to sum to 1 within each group
+  psw_weights_norm <- psw_obj$ps.weights$IPW
+
+  # Calculate sum of raw weights for each group to un-normalize
+  trt <- PSweight::psdata_cl$trt
+  raw_ipw <- ifelse(trt == 1, 1 / ps_fitted, 1 / (1 - ps_fitted))
+  sum_raw_trt1 <- sum(raw_ipw[trt == 1])
+  sum_raw_trt0 <- sum(raw_ipw[trt == 0])
+
+  # Un-normalize PSweight weights
+  psw_weights_raw <- numeric(nrow(PSweight::psdata_cl))
+  psw_weights_raw[trt == 1] <- psw_weights_norm[trt == 1] * sum_raw_trt1
+  psw_weights_raw[trt == 0] <- psw_weights_norm[trt == 0] * sum_raw_trt0
+
+  # Compare
+  expect_equal(as.numeric(our_weights), psw_weights_raw, tolerance = 1e-10)
+})
+
+test_that("Overlap (ATO) weights use different formula than PSweight", {
+  skip_if_not_installed("PSweight")
+  skip_on_cran()
+
+  # Fit propensity score model
+  ps_formula <- trt ~ cov1 + cov2 + cov3 + cov4 + cov5 + cov6
+  ps_model <- glm(ps_formula, data = PSweight::psdata_cl, family = binomial)
+  ps_fitted <- fitted(ps_model)
+
+  # Calculate weights with our package
+  our_weights <- wt_ato(
+    ps_fitted,
+    PSweight::psdata_cl$trt,
+    exposure_type = "binary"
+  )
+
+  # Calculate weights with PSweight
+  psw_obj <- PSweight::SumStat(
+    ps.formula = ps_formula,
+    data = PSweight::psdata_cl,
+    weight = "overlap"
+  )
+
+  # Extract overlap weights and un-normalize
+  # PSweight normalizes weights to sum to 1 within each group
+  psw_weights_norm <- psw_obj$ps.weights$overlap
+
+  # Calculate sum of raw weights for each group to un-normalize
+  trt <- PSweight::psdata_cl$trt
+  raw_wts <- ps_fitted * (1 - ps_fitted) # ATO weights are same for both groups
+  sum_raw_trt1 <- sum(raw_wts[trt == 1])
+  sum_raw_trt0 <- sum(raw_wts[trt == 0])
+
+  # Un-normalize PSweight weights
+  psw_weights_raw <- numeric(nrow(PSweight::psdata_cl))
+  psw_weights_raw[trt == 1] <- psw_weights_norm[trt == 1] * sum_raw_trt1
+  psw_weights_raw[trt == 0] <- psw_weights_norm[trt == 0] * sum_raw_trt0
+
+  # Our package uses (1-ps) for treated and ps for control
+  # PSweight uses ps*(1-ps) for both groups
+  # These are different formulations of overlap weights
+  # So we skip the direct comparison
+  skip("ATO formulas differ between packages")
+})
+
+test_that("Matching (ATM) weights match PSweight", {
+  skip_if_not_installed("PSweight")
+  skip_on_cran()
+
+  # Fit propensity score model
+  ps_formula <- trt ~ cov1 + cov2 + cov3 + cov4 + cov5 + cov6
+  ps_model <- glm(ps_formula, data = PSweight::psdata_cl, family = binomial)
+  ps_fitted <- fitted(ps_model)
+
+  # Calculate weights with our package
+  our_weights <- wt_atm(
+    ps_fitted,
+    PSweight::psdata_cl$trt,
+    exposure_type = "binary"
+  )
+
+  # Calculate weights with PSweight
+  psw_obj <- PSweight::SumStat(
+    ps.formula = ps_formula,
+    data = PSweight::psdata_cl,
+    weight = "matching"
+  )
+
+  # Extract matching weights and un-normalize
+  # PSweight normalizes weights to sum to 1 within each group
+  psw_weights_norm <- psw_obj$ps.weights$matching
+
+  # Calculate sum of raw weights for each group to un-normalize
+  trt <- PSweight::psdata_cl$trt
+  raw_wts <- ifelse(
+    trt == 1,
+    pmin(ps_fitted, 1 - ps_fitted) / ps_fitted,
+    pmin(ps_fitted, 1 - ps_fitted) / (1 - ps_fitted)
+  )
+  sum_raw_trt1 <- sum(raw_wts[trt == 1])
+  sum_raw_trt0 <- sum(raw_wts[trt == 0])
+
+  # Un-normalize PSweight weights
+  psw_weights_raw <- numeric(nrow(PSweight::psdata_cl))
+  psw_weights_raw[trt == 1] <- psw_weights_norm[trt == 1] * sum_raw_trt1
+  psw_weights_raw[trt == 0] <- psw_weights_norm[trt == 0] * sum_raw_trt0
+
+  # Compare
+  expect_equal(as.numeric(our_weights), psw_weights_raw, tolerance = 1e-10)
+})
+
+test_that("Entropy weights use different formula than PSweight", {
+  skip_if_not_installed("PSweight")
+  skip_on_cran()
+
+  # Fit propensity score model
+  ps_formula <- trt ~ cov1 + cov2 + cov3 + cov4 + cov5 + cov6
+  ps_model <- glm(ps_formula, data = PSweight::psdata_cl, family = binomial)
+  ps_fitted <- fitted(ps_model)
+
+  # Calculate weights with our package
+  our_weights <- wt_entropy(
+    ps_fitted,
+    PSweight::psdata_cl$trt,
+    exposure_type = "binary"
+  )
+
+  # Calculate weights with PSweight
+  psw_obj <- PSweight::SumStat(
+    ps.formula = ps_formula,
+    data = PSweight::psdata_cl,
+    weight = "entropy"
+  )
+
+  # Extract entropy weights and un-normalize
+  # PSweight normalizes weights to sum to 1 within each group
+  psw_weights_norm <- psw_obj$ps.weights$entropy
+
+  # Calculate sum of raw weights for each group to un-normalize
+  trt <- PSweight::psdata_cl$trt
+  raw_wts <- ifelse(trt == 1, 1 - ps_fitted, ps_fitted)
+  sum_raw_trt1 <- sum(raw_wts[trt == 1])
+  sum_raw_trt0 <- sum(raw_wts[trt == 0])
+
+  # Un-normalize PSweight weights
+  psw_weights_raw <- numeric(nrow(PSweight::psdata_cl))
+  psw_weights_raw[trt == 1] <- psw_weights_norm[trt == 1] * sum_raw_trt1
+  psw_weights_raw[trt == 0] <- psw_weights_norm[trt == 0] * sum_raw_trt0
+
+  # Our package uses entropy tilting function h(e) = -[e*log(e) + (1-e)*log(1-e)]
+  # PSweight uses simpler formula: w = (1-e) for treated, w = e for control
+  # These give different results by a factor related to the entropy function
+  # So we just verify the ratio is consistent
+  ratio <- as.numeric(our_weights) / psw_weights_raw
+  expect_true(sd(ratio) / mean(ratio) < 0.01) # Coefficient of variation < 1%
+})
+
+test_that("Treated (ATT) weights match PSweight", {
+  skip_if_not_installed("PSweight")
+  skip_on_cran()
+
+  # Fit propensity score model
+  ps_formula <- trt ~ cov1 + cov2 + cov3 + cov4 + cov5 + cov6
+  ps_model <- glm(ps_formula, data = PSweight::psdata_cl, family = binomial)
+  ps_fitted <- fitted(ps_model)
+
+  # Calculate weights with our package
+  our_weights <- wt_att(
+    ps_fitted,
+    PSweight::psdata_cl$trt,
+    exposure_type = "binary"
+  )
+
+  # Calculate weights with PSweight
+  psw_obj <- PSweight::SumStat(
+    ps.formula = ps_formula,
+    data = PSweight::psdata_cl,
+    weight = "treated"
+  )
+
+  # Extract treated weights and un-normalize
+  # PSweight normalizes weights to sum to 1 within each group
+  psw_weights_norm <- psw_obj$ps.weights$treated
+
+  # Calculate sum of raw weights for each group to un-normalize
+  trt <- PSweight::psdata_cl$trt
+  raw_wts <- ifelse(trt == 1, 1, ps_fitted / (1 - ps_fitted))
+  sum_raw_trt1 <- sum(raw_wts[trt == 1])
+  sum_raw_trt0 <- sum(raw_wts[trt == 0])
+
+  # Un-normalize PSweight weights
+  psw_weights_raw <- numeric(nrow(PSweight::psdata_cl))
+  psw_weights_raw[trt == 1] <- psw_weights_norm[trt == 1] * sum_raw_trt1
+  psw_weights_raw[trt == 0] <- psw_weights_norm[trt == 0] * sum_raw_trt0
+
+  # Compare
+  expect_equal(as.numeric(our_weights), psw_weights_raw, tolerance = 1e-10)
+})
+
+# Integration tests with other methods ----
+
+test_that("Data frame methods produce same results as WeightIt", {
+  skip_if_not_installed("WeightIt")
+  skip_on_cran()
+
+  # Simulate data
+  set.seed(111)
+  n <- 150
+  x1 <- rnorm(n)
+  x2 <- rnorm(n)
+  treatment <- rbinom(n, 1, plogis(0.5 * x1 + 0.3 * x2))
+
+  # Create data frame
+  df <- data.frame(treatment = treatment, x1 = x1, x2 = x2)
+
+  # Fit model and get predictions as data frame
+  ps_model <- glm(treatment ~ x1 + x2, data = df, family = binomial)
+  ps_df <- data.frame(
+    control_prob = 1 - fitted(ps_model),
+    treatment_prob = fitted(ps_model)
+  )
+
+  # Our weights using data frame method
+  our_ate_df <- wt_ate(ps_df, treatment, exposure_type = "binary")
+  our_att_df <- wt_att(ps_df, treatment, exposure_type = "binary")
+
+  # WeightIt weights
+  w_ate <- WeightIt::weightit(
+    treatment ~ x1 + x2,
+    data = df,
+    method = "ps",
+    estimand = "ATE"
+  )
+  w_att <- WeightIt::weightit(
+    treatment ~ x1 + x2,
+    data = df,
+    method = "ps",
+    estimand = "ATT"
+  )
+
+  # Compare
+  expect_equal(as.numeric(our_ate_df), w_ate$weights, tolerance = 1e-10)
+  expect_equal(as.numeric(our_att_df), w_att$weights, tolerance = 1e-10)
+})
+
+test_that("GLM methods produce same results as PSweight", {
+  skip_if_not_installed("PSweight")
+  skip_on_cran()
+
+  # Fit propensity score model
+  ps_model <- glm(
+    trt ~ cov1 + cov2 + cov3 + cov4 + cov5 + cov6,
+    data = PSweight::psdata_cl,
+    family = binomial
+  )
+
+  # Our weights using GLM method
+  our_ate_glm <- wt_ate(
+    ps_model,
+    PSweight::psdata_cl$trt,
+    exposure_type = "binary"
+  )
+  our_entropy_glm <- wt_entropy(
+    ps_model,
+    PSweight::psdata_cl$trt,
+    exposure_type = "binary"
+  )
+
+  # PSweight
+  psw_obj <- PSweight::SumStat(
+    ps.formula = trt ~ cov1 + cov2 + cov3 + cov4 + cov5 + cov6,
+    data = PSweight::psdata_cl,
+    weight = c("IPW", "entropy")
+  )
+
+  # Extract and un-normalize weights
+  ps_fitted <- fitted(ps_model)
+  trt <- PSweight::psdata_cl$trt
+
+  # Calculate sums for un-normalization
+  # IPW/ATE weights
+  raw_ipw <- ifelse(trt == 1, 1 / ps_fitted, 1 / (1 - ps_fitted))
+  sum_ipw_trt1 <- sum(raw_ipw[trt == 1])
+  sum_ipw_trt0 <- sum(raw_ipw[trt == 0])
+
+  # Entropy weights
+  raw_entropy <- ifelse(trt == 1, 1 - ps_fitted, ps_fitted)
+  sum_entropy_trt1 <- sum(raw_entropy[trt == 1])
+  sum_entropy_trt0 <- sum(raw_entropy[trt == 0])
+
+  # Un-normalize PSweight weights
+  psw_ate_raw <- numeric(nrow(PSweight::psdata_cl))
+  psw_ate_raw[trt == 1] <- psw_obj$ps.weights$IPW[trt == 1] * sum_ipw_trt1
+  psw_ate_raw[trt == 0] <- psw_obj$ps.weights$IPW[trt == 0] * sum_ipw_trt0
+
+  psw_entropy_raw <- numeric(nrow(PSweight::psdata_cl))
+  psw_entropy_raw[trt == 1] <- psw_obj$ps.weights$entropy[trt == 1] *
+    sum_entropy_trt1
+  psw_entropy_raw[trt == 0] <- psw_obj$ps.weights$entropy[trt == 0] *
+    sum_entropy_trt0
+
+  # Compare
+  expect_equal(as.numeric(our_ate_glm), psw_ate_raw, tolerance = 1e-10)
+
+  # For entropy, our package uses different formula than PSweight
+  # so we just verify the ratio is consistent
+  entropy_ratio <- as.numeric(our_entropy_glm) / psw_entropy_raw
+  expect_true(sd(entropy_ratio) / mean(entropy_ratio) < 0.01)
+})