From 5f8604ed5b0e5d4ac21571a4347b8591b14dd9f7 Mon Sep 17 00:00:00 2001
From: Malcolm Barrett <malcolmbarrett@gmail.com>
Date: Wed, 16 Sep 2026 10:21:45 -0400
Subject: [PATCH 1/7] add failing tests for the Laplace scale estimator and the
 shared scale seam

---
 tests/testthat/helper-density.R          |  25 ++
 tests/testthat/test-density.R            | 216 +++++++++++++++-
 tests/testthat/test-ipw-psi.R            | 129 ++++++++++
 tests/testthat/test-weights-continuous.R | 299 ++++++++++++++++++++++-
 4 files changed, 657 insertions(+), 12 deletions(-)

diff --git a/tests/testthat/helper-density.R b/tests/testthat/helper-density.R
index 5f0d2b0d..01b4d6fb 100644
--- a/tests/testthat/helper-density.R
+++ b/tests/testthat/helper-density.R
@@ -50,3 +50,28 @@ t_scale_mle <- function(residuals, df) {
 
   exp(log_scale)
 }
+
+# The maximum likelihood scale of a Laplace density fit to `residuals`, for the
+# tests that hold the package's closed form against an oracle.
+#
+# The likelihood is maximized directly, over the log of the scale, rather than
+# evaluated at the mean absolute residual the closed form returns, so that the
+# closed form is the thing under test rather than the thing assumed. The search
+# is bracketed off the root mean square, which is the estimator the closed form
+# is not, and reaches a factor of e^6 either side of it.
+laplace_scale_mle <- function(residuals) {
+  residuals <- residuals[!is.na(residuals)]
+
+  loglik <- function(log_scale) {
+    sum(-log(2) - log_scale - abs(residuals) / exp(log_scale))
+  }
+
+  exp(
+    stats::optimize(
+      loglik,
+      interval = log(sqrt(mean(residuals^2))) + c(-6, 6),
+      maximum = TRUE,
+      tol = .Machine$double.eps^0.5
+    )$maximum
+  )
+}
diff --git a/tests/testthat/test-density.R b/tests/testthat/test-density.R
index efc331ea..ec00f01e 100644
--- a/tests/testthat/test-density.R
+++ b/tests/testthat/test-density.R
@@ -63,7 +63,7 @@ test_that("dens_laplace() is the standard Laplace density", {
   spec <- dens_laplace()
 
   expect_s3_class(spec, "propensity_density")
-  expect_named(spec, c("family", "params", "fn"))
+  expect_named(spec, c("family", "params", "fn", "sigma_method"))
   expect_identical(spec$family, "laplace")
   expect_identical(spec$params, list())
   expect_equal(spec$fn(z), exp(-abs(z)) / 2)
@@ -103,6 +103,37 @@ test_that("dens_t() refuses a scale estimator it does not have", {
   )
 })
 
+test_that("dens_laplace() records how its scale is estimated", {
+  # The Laplace has a scale of its own, and the same choice of estimator for it
+  # that the t has: the root mean square that spreads every other family, or
+  # maximum likelihood under the Laplace itself, which is the default. It is
+  # recorded where the t records it, beside the parameters rather than among
+  # them, because `laplace` is the same density however its scale was arrived
+  # at.
+  spec <- dens_laplace()
+
+  expect_named(spec, c("family", "params", "fn", "sigma_method"))
+  expect_identical(spec$sigma_method, "mle")
+  expect_identical(dens_laplace(sigma_method = "mle")$sigma_method, "mle")
+  expect_identical(dens_laplace(sigma_method = "rms")$sigma_method, "rms")
+
+  # The estimator is not a parameter, so it reaches neither the parameters nor
+  # the printed density.
+  expect_identical(dens_laplace(sigma_method = "rms")$params, list())
+  expect_identical(format(dens_laplace(sigma_method = "rms")), "laplace")
+  expect_equal(
+    dens_laplace(sigma_method = "rms")$fn(density_z()),
+    exp(-abs(density_z())) / 2
+  )
+})
+
+test_that("dens_laplace() refuses a scale estimator it does not have", {
+  expect_error(
+    dens_laplace(sigma_method = "median"),
+    class = "rlang_error"
+  )
+})
+
 test_that("dens_kernel() records its stats::density() controls", {
   spec <- dens_kernel()
 
@@ -642,6 +673,189 @@ test_that("density_eval() drops the shape of a matrix result", {
   )
 })
 
+# ---- the scale a density is read at -----------------------------------------
+
+# `density_scale_estimate()` is the one estimator of the spread both densities
+# of a continuous exposure's ratio are read at, and `laplace_sigma_mle()` is the
+# closed form it reaches for the Laplace. The equation the estimate is the root
+# of lives beside the stacked system that solves it, and is tested in
+# tests/testthat/test-ipw-psi.R. All three are internal, and are called here by
+# name the way the rest of the suite calls internals.
+
+# Residuals from a law with tails heavier than a normal's, so that the mean
+# absolute residual and the root mean square stand about a factor of sqrt(2)
+# apart and no test below can pass on one where it meant the other.
+density_scale_residuals <- withr::with_seed(20260902, {
+  u <- stats::runif(200) - 0.5
+
+  -sign(u) * log(1 - 2 * abs(u))
+})
+
+test_that("laplace_sigma_mle() is the mean absolute residual", {
+  residuals <- density_scale_residuals
+
+  # The oracle maximizes the Laplace likelihood itself rather than evaluating
+  # the closed form, so what is under test is that the closed form is where the
+  # likelihood is maximized. The tolerance is the search's rather than the
+  # estimator's: a maximum found by search settles to about the square root of
+  # the machine epsilon.
+  expect_equal(
+    laplace_sigma_mle(residuals),
+    laplace_scale_mle(residuals),
+    tolerance = 1e-6
+  )
+  expect_equal(laplace_sigma_mle(residuals), mean(abs(residuals)))
+
+  # It is the estimator the root mean square is not.
+  expect_false(isTRUE(all.equal(
+    laplace_sigma_mle(residuals),
+    sqrt(mean(residuals^2)),
+    tolerance = 0.1
+  )))
+
+  # A residual that is not there is one the likelihood cannot read, exactly as
+  # it is one the root mean square does not average.
+  gappy <- residuals
+  gappy[c(3, 17)] <- NA
+
+  expect_equal(laplace_sigma_mle(gappy), mean(abs(residuals[-c(3, 17)])))
+})
+
+test_that("density_scale_estimate() takes the root mean square by default", {
+  residuals <- density_scale_residuals
+  rms <- sqrt(mean(residuals^2))
+
+  # A density that names no estimator, and every family that carries no
+  # `sigma_method` at all, is spread by the moment estimator, which is the
+  # maximum likelihood scale of the normal and the spread every family has
+  # always been read at.
+  expect_equal(density_scale_estimate(residuals, NULL), rms)
+  expect_equal(density_scale_estimate(residuals, dens_normal()), rms)
+  expect_equal(density_scale_estimate(residuals, dens_kernel()), rms)
+  expect_equal(
+    density_scale_estimate(residuals, dens_fn(function(z) stats::dnorm(z))),
+    rms
+  )
+
+  # A family that has an estimator of its own and was asked not to use it is
+  # read at the moment estimator too. The estimator follows what the density
+  # asks for, not the family it names.
+  expect_equal(
+    density_scale_estimate(residuals, dens_laplace(sigma_method = "rms")),
+    rms
+  )
+  expect_equal(
+    density_scale_estimate(residuals, dens_t(4, sigma_method = "rms")),
+    rms
+  )
+})
+
+test_that("density_scale_estimate() reads the family a density is fit under", {
+  residuals <- density_scale_residuals
+
+  expect_equal(
+    density_scale_estimate(residuals, dens_laplace()),
+    laplace_scale_mle(residuals),
+    tolerance = 1e-6
+  )
+  expect_equal(
+    density_scale_estimate(residuals, dens_laplace(sigma_method = "mle")),
+    mean(abs(residuals))
+  )
+  expect_equal(
+    density_scale_estimate(residuals, dens_t(4, sigma_method = "mle")),
+    t_scale_mle(residuals, df = 4),
+    tolerance = 1e-8
+  )
+
+  # The Laplace answer is the mean absolute residual rather than the moment
+  # estimator, which these residuals put about 40 percent away from it, so an
+  # estimate that reached the moment estimator could not pass the assertions
+  # above. The t answer is not separated here and is not asked to be: a t with
+  # four degrees of freedom fit to Laplace residuals lands within a thousandth
+  # of their mean absolute deviation, so what tells the two estimators apart is
+  # the oracle each is held against rather than a gap between their answers.
+  expect_false(isTRUE(all.equal(
+    density_scale_estimate(residuals, dens_laplace()),
+    density_scale_estimate(residuals, dens_normal()),
+    tolerance = 0.05
+  )))
+})
+
+test_that("density_scale_estimate() reads the residuals that are there", {
+  residuals <- density_scale_residuals
+  gappy <- residuals
+  gappy[c(3, 17)] <- NA
+  present <- residuals[-c(3, 17)]
+
+  expect_equal(
+    density_scale_estimate(gappy, dens_laplace()),
+    mean(abs(present))
+  )
+  expect_equal(
+    density_scale_estimate(gappy, dens_normal()),
+    sqrt(mean(present^2))
+  )
+})
+
+test_that("a scale estimator the package does not have is refused", {
+  # Only the t and the Laplace can ask to be fit under themselves, because no
+  # other constructor takes `sigma_method`. A specification that asks anyway is
+  # one the package cannot build, and it is refused rather than passed over for
+  # a spread of nothing, which would leave every weight missing and the stacked
+  # system solving an equation no weight was built at.
+  spec <- new_density_spec(
+    "normal",
+    fn = function(z) stats::dnorm(z),
+    sigma_method = "mle"
+  )
+
+  expect_error(
+    density_scale_estimate(density_scale_residuals, spec),
+    class = "propensity_density_error"
+  )
+
+  # The refusal reads against the function the user called rather than against
+  # the estimator, which is internal.
+  cnd <- rlang::catch_cnd(
+    density_scale_estimate(
+      density_scale_residuals,
+      spec,
+      call = rlang::call2("wt_ate")
+    ),
+    classes = "propensity_density_error"
+  )
+
+  expect_identical(conditionCall(cnd), quote(wt_ate()))
+})
+
+test_that("density_specs_agree() tells the scale estimators apart", {
+  # Two Laplace weights built with different scale estimators are ratios of
+  # densities of different widths, so the records they carry do not agree and
+  # the weights do not combine silently.
+  expect_true(density_specs_agree(dens_laplace(), dens_laplace()))
+  expect_true(density_specs_agree(
+    dens_laplace(sigma_method = "rms"),
+    dens_laplace(sigma_method = "rms")
+  ))
+  expect_true(density_specs_agree(
+    dens_laplace(sigma_method = "mle"),
+    dens_laplace()
+  ))
+
+  expect_false(density_specs_agree(
+    dens_laplace(sigma_method = "mle"),
+    dens_laplace(sigma_method = "rms")
+  ))
+
+  # The family still tells them apart on its own, whatever the estimator.
+  expect_false(density_specs_agree(dens_laplace(), dens_normal()))
+  expect_false(density_specs_agree(
+    dens_laplace(sigma_method = "rms"),
+    dens_t(4, sigma_method = "rms")
+  ))
+})
+
 # ---- print and format -------------------------------------------------------
 
 # The printed forms are pinned as exact text here rather than by snapshot: the
diff --git a/tests/testthat/test-ipw-psi.R b/tests/testthat/test-ipw-psi.R
index 3324cd02..c96d0a3b 100644
--- a/tests/testthat/test-ipw-psi.R
+++ b/tests/testthat/test-ipw-psi.R
@@ -1632,3 +1632,132 @@ test_that("a degenerate maximum likelihood scale is refused in the caller's name
 
   expect_identical(conditionCall(cnd), quote(ipw()))
 })
+
+# ---- the equation a density's scale is the root of ---------------------------
+
+# The residuals the scale rows are read at, heavier tailed than a normal's so
+# that the estimators the rows belong to answer differently.
+density_scale_row_residuals <- withr::with_seed(20260903, {
+  u <- stats::runif(150) - 0.5
+
+  -sign(u) * log(1 - 2 * abs(u))
+})
+
+test_that("density_scale_row() is the moment row unless the density asks otherwise", {
+  residuals <- c(-2.1, 0.4, 1.7, -0.3, 3.2)
+  sigma2 <- 1.6
+  moment <- matrix(residuals^2 - sigma2, nrow = 1)
+
+  # The uncentered second moment, which is the equation the root mean square is
+  # the root of and the row the continuous block has always carried.
+  expect_equal(density_scale_row(residuals, sigma2, NULL), moment)
+  expect_equal(density_scale_row(residuals, sigma2, dens_normal()), moment)
+  expect_equal(density_scale_row(residuals, sigma2, dens_kernel()), moment)
+
+  # A family with an estimator of its own that was asked not to use it emits the
+  # moment row too: the row follows what the density asks for, not the family it
+  # names, which is what keeps it the equation the weights were built at.
+  expect_equal(
+    density_scale_row(residuals, sigma2, dens_laplace(sigma_method = "rms")),
+    moment
+  )
+  expect_equal(
+    density_scale_row(residuals, sigma2, dens_t(4, sigma_method = "rms")),
+    moment
+  )
+})
+
+test_that("density_scale_row() is the score of the family a density is fit under", {
+  residuals <- c(-2.1, 0.4, 1.7, -0.3, 3.2)
+  sigma2 <- 1.6
+
+  # The Laplace score for the scale, multiplied through by the scale, is the
+  # absolute residual against the scale itself, so it is written against the
+  # square root of the parameter the block carries.
+  expect_equal(
+    density_scale_row(residuals, sigma2, dens_laplace()),
+    matrix(abs(residuals) - sqrt(sigma2), nrow = 1)
+  )
+  expect_equal(
+    density_scale_row(residuals, sigma2, dens_t(4, sigma_method = "mle")),
+    matrix(5 * residuals^2 / (4 * sigma2 + residuals^2) - 1, nrow = 1)
+  )
+
+  # One row for the one parameter it estimates, whatever the family.
+  expect_identical(
+    dim(density_scale_row(residuals, sigma2, dens_laplace())),
+    c(1L, length(residuals))
+  )
+})
+
+test_that("the scale a density is estimated at is the root of the row that estimates it", {
+  # The seam the whole stacked system rests on: one estimator and one equation,
+  # so the spread the weights were built at is exactly where the row averages to
+  # zero. The preflight that rebuilds the weights at the seed holds only while
+  # that is so, and it holds only because the two are keyed on the same
+  # question.
+  residuals <- density_scale_row_residuals
+
+  # Named by the estimator as well as by the family, because the printed density
+  # is the same whichever estimator it asked for and a failure has to say which
+  # of the five came apart.
+  densities <- list(
+    "normal" = dens_normal(),
+    "laplace, mle" = dens_laplace(),
+    "laplace, rms" = dens_laplace(sigma_method = "rms"),
+    "t(df = 4), mle" = dens_t(4, sigma_method = "mle"),
+    "t(df = 4), rms" = dens_t(4, sigma_method = "rms")
+  )
+
+  for (name in names(densities)) {
+    scale <- density_scale_estimate(residuals, densities[[name]])
+    row <- density_scale_row(residuals, scale^2, densities[[name]])
+
+    expect_lt(abs(mean(row)), 1e-8, label = name)
+  }
+})
+
+test_that("density_scale_row() leaves a missing residual missing", {
+  residuals <- c(-2.1, NA, 1.7, -0.3, 3.2)
+  sigma2 <- 1.6
+
+  laplace <- density_scale_row(residuals, sigma2, dens_laplace())
+  moment <- density_scale_row(residuals, sigma2, dens_normal())
+
+  # A unit whose exposure is missing leaves no residual for the equation to
+  # read, and the row it would have contributed is missing rather than zero,
+  # which would be a unit the equation is solved at. The row keeps its width
+  # either way, since the block it sits in is one row for each parameter and one
+  # column for each observation.
+  expect_equal(laplace, matrix(abs(residuals) - sqrt(sigma2), nrow = 1))
+  expect_equal(moment, matrix(residuals^2 - sigma2, nrow = 1))
+
+  expect_true(is.na(laplace[1, 2]))
+  expect_false(anyNA(laplace[1, -2]))
+  expect_identical(dim(laplace), c(1L, length(residuals)))
+})
+
+test_that("a scale row the package cannot write is refused", {
+  # The default arm aborts rather than falling through to nothing. Only the t
+  # and the Laplace can ask to be fit under themselves, so a specification that
+  # asks anyway is one the package cannot write a row for, and a block of
+  # nothing would leave the stacked system solving for a parameter no equation
+  # names.
+  spec <- new_density_spec(
+    "normal",
+    fn = function(z) stats::dnorm(z),
+    sigma_method = "mle"
+  )
+
+  expect_error(
+    density_scale_row(c(-1, 0.5, 2), 1.4, spec),
+    class = "propensity_density_error"
+  )
+
+  cnd <- rlang::catch_cnd(
+    density_scale_row(c(-1, 0.5, 2), 1.4, spec, call = rlang::call2("ipw")),
+    classes = "propensity_density_error"
+  )
+
+  expect_identical(conditionCall(cnd), quote(ipw()))
+})
diff --git a/tests/testthat/test-weights-continuous.R b/tests/testthat/test-weights-continuous.R
index 9e7fcdf3..add56aa7 100644
--- a/tests/testthat/test-weights-continuous.R
+++ b/tests/testthat/test-weights-continuous.R
@@ -676,14 +676,23 @@ test_that("the refusal of a matrix of conditional means reads plainly", {
 
 test_that("the weights record the density family they were built from", {
   families <- list(
-    normal = list(input = "normal", recorded = "normal"),
-    laplace = list(input = dens_laplace(), recorded = "laplace"),
-    t = list(input = dens_t(df = 4), recorded = "t(df = 4)"),
+    normal = list(input = "normal", recorded = "normal", sigma = "pooled"),
+    laplace = list(
+      input = dens_laplace(),
+      recorded = "laplace",
+      sigma = "mle"
+    ),
+    t = list(input = dens_t(df = 4), recorded = "t(df = 4)", sigma = "pooled"),
     kernel = list(
       input = dens_kernel(adjust = 1.5),
-      recorded = 'kernel(bw = "nrd0", adjust = 1.5, kernel = "gaussian", n = 512)'
+      recorded = 'kernel(bw = "nrd0", adjust = 1.5, kernel = "gaussian", n = 512)',
+      sigma = "pooled"
     ),
-    fn = list(input = function(z) stats::dlogis(z), recorded = "function")
+    fn = list(
+      input = function(z) stats::dlogis(z),
+      recorded = "function",
+      sigma = "pooled"
+    )
   )
 
   for (family in families) {
@@ -695,10 +704,12 @@ test_that("the weights record the density family they were built from", {
     expect_identical(format(record$density), family$recorded)
     expect_identical(exposure_type(weights), "continuous")
 
-    # The rest of the record is what it was: the family says nothing about
-    # what stabilized the weights or where the spread came from.
+    # The family says nothing about what stabilized the weights, and says where
+    # the spread came from only when it has an estimator of its own: a family
+    # read at the moment estimator records the pooled spread whether it took
+    # that estimator for want of another or asked for it.
     expect_identical(record$numerator, "marginal")
-    expect_identical(record$sigma, "pooled")
+    expect_identical(record$sigma, family$sigma)
   }
 })
 
@@ -2178,15 +2189,26 @@ test_that("an interpolated numerator that dips below zero is refused", {
 # ---- what the integrated weights record -------------------------------------
 
 test_that("integrated weights record the numerator they were built from", {
-  for (family in list("normal", dens_t(df = 4), "laplace", "kernel")) {
+  families <- list(
+    list(input = "normal", sigma = "pooled"),
+    list(input = dens_t(df = 4), sigma = "pooled"),
+    list(input = "laplace", sigma = "mle"),
+    list(input = "kernel", sigma = "pooled")
+  )
+
+  for (family in families) {
     record <- density_meta(continuous_density_wt(
-      .density = family,
+      .density = family$input,
       numerator = "integrated"
     ))
 
     expect_s3_class(record, "propensity_density_meta")
     expect_identical(record$numerator, "integrated")
-    expect_identical(record$sigma, "pooled")
+
+    # An integrated numerator reads the conditional density in both places, so
+    # it carries the conditional spread and records the estimator that spread
+    # was reached by.
+    expect_identical(record$sigma, family$sigma)
   }
 
   # And the marginal numerator still records itself, written out or not.
@@ -2757,6 +2779,261 @@ test_that("a likelihood with no maximum at a positive scale is refused", {
   )
 })
 
+# ---- the Laplace scale by maximum likelihood --------------------------------
+
+# A continuous-exposure problem whose residuals follow the family the weights
+# assert. The scale of a Laplace is its mean absolute deviation and its standard
+# deviation is sqrt(2) times that, so the two estimators of the spread stand far
+# enough apart here that no test below can pass on one where it meant the other.
+continuous_laplace_data <- local({
+  set.seed(20260904)
+
+  n <- 300
+  x <- rnorm(n)
+  u <- runif(n) - 0.5
+  exposure <- 1 + 0.7 * x - sign(u) * log(1 - 2 * abs(u))
+
+  list(
+    x = x,
+    exposure = exposure,
+    mu = as.numeric(fitted(lm(exposure ~ x)))
+  )
+})
+
+continuous_laplace_residuals <- continuous_laplace_data$exposure -
+  continuous_laplace_data$mu
+
+# Weights for that problem. They are unstabilized unless the test asks for the
+# marginal numerator, so that the conditional density is on its own wherever
+# what is under test is the conditional spread.
+continuous_laplace_wt <- function(
+  ...,
+  exposure = continuous_laplace_data$exposure,
+  stabilize = FALSE
+) {
+  wt_ate(
+    continuous_laplace_data$mu,
+    exposure,
+    exposure_type = "continuous",
+    stabilize = stabilize,
+    ...
+  )
+}
+
+# The Laplace density on the exposure's own scale: the standard density read at
+# the residual standardized by `scale`, divided by the scale that standardized
+# it.
+laplace_density <- function(residuals, scale) {
+  exp(-abs(residuals / scale)) / 2 / scale
+}
+
+test_that("a Laplace denominator is read at the scale of the Laplace", {
+  residuals <- continuous_laplace_residuals
+  scale <- mean(abs(residuals))
+
+  # The closed form is where the Laplace likelihood is maximized, held against
+  # an oracle that maximizes it by search.
+  expect_equal(scale, laplace_scale_mle(residuals), tolerance = 1e-6)
+
+  weights <- continuous_laplace_wt(.density = dens_laplace())
+
+  expect_equal(
+    as.numeric(weights),
+    1 / laplace_density(residuals, scale),
+    tolerance = 1e-12
+  )
+
+  # The spread the weights were built at is the scale rather than the standard
+  # deviation, and these residuals tell the two apart.
+  rms <- sqrt(mean(residuals^2))
+
+  expect_false(isTRUE(all.equal(
+    as.numeric(weights),
+    1 / laplace_density(residuals, rms),
+    tolerance = 1e-2
+  )))
+})
+
+test_that("a Laplace denominator asked for the root mean square is read at it", {
+  residuals <- continuous_laplace_residuals
+  rms <- sqrt(mean(residuals^2))
+
+  weights <- continuous_laplace_wt(
+    .density = dens_laplace(sigma_method = "rms")
+  )
+
+  expect_equal(
+    as.numeric(weights),
+    1 / laplace_density(residuals, rms),
+    tolerance = 1e-12
+  )
+
+  # Which is the same ratio as the one a caller who supplied that number by hand
+  # would be given, since the estimator is the only thing the setting changes.
+  supplied <- continuous_laplace_wt(
+    .density = dens_laplace(sigma_method = "rms"),
+    .sigma = rms
+  )
+
+  expect_equal(as.numeric(weights), as.numeric(supplied), tolerance = 1e-12)
+})
+
+test_that("the marginal numerator of a Laplace ratio is read at the same estimator", {
+  exposure <- continuous_laplace_data$exposure
+  residuals <- continuous_laplace_residuals
+  scale <- mean(abs(residuals))
+  mu_a <- mean(exposure)
+  scale_a <- mean(abs(exposure - mu_a))
+
+  weights <- continuous_laplace_wt(
+    .density = dens_laplace(),
+    stabilize = TRUE
+  )
+
+  expect_equal(
+    as.numeric(weights),
+    laplace_density(exposure - mu_a, scale_a) /
+      laplace_density(residuals, scale),
+    tolerance = 1e-12
+  )
+
+  # The numerator's location is the exposure's own mean, as it has always been;
+  # its spread alone is what the family decides. A numerator spread by the root
+  # mean square while the denominator is spread by the scale is a ratio of two
+  # densities of different widths, and is the ratio these weights are not.
+  sd_a <- sqrt(mean((exposure - mu_a)^2))
+
+  expect_false(isTRUE(all.equal(
+    as.numeric(weights),
+    laplace_density(exposure - mu_a, sd_a) /
+      laplace_density(residuals, scale),
+    tolerance = 1e-3
+  )))
+})
+
+test_that("a Laplace ratio asked for the root mean square is read at it on both halves", {
+  # The estimator follows what the density asks for rather than the family it
+  # names, on the numerator as well as on the denominator. An estimator keyed on
+  # the family alone would spread the denominator by the root mean square, which
+  # is what was asked for, and the numerator by the mean absolute deviation,
+  # which is what the family implies, leaving the weights the ratio of neither
+  # request.
+  exposure <- continuous_laplace_data$exposure
+  residuals <- continuous_laplace_residuals
+  rms <- sqrt(mean(residuals^2))
+  mu_a <- mean(exposure)
+  sd_a <- sqrt(mean((exposure - mu_a)^2))
+
+  weights <- continuous_laplace_wt(
+    .density = dens_laplace(sigma_method = "rms"),
+    stabilize = TRUE
+  )
+
+  expect_equal(
+    as.numeric(weights),
+    laplace_density(exposure - mu_a, sd_a) / laplace_density(residuals, rms),
+    tolerance = 1e-12
+  )
+
+  expect_false(isTRUE(all.equal(
+    as.numeric(weights),
+    laplace_density(exposure - mu_a, mean(abs(exposure - mu_a))) /
+      laplace_density(residuals, rms),
+    tolerance = 1e-3
+  )))
+
+  # The two settings are two different sets of weights, so neither test above
+  # could be passing on the other's numbers.
+  under_mle <- continuous_laplace_wt(
+    .density = dens_laplace(),
+    stabilize = TRUE
+  )
+
+  expect_false(isTRUE(all.equal(
+    as.numeric(weights),
+    as.numeric(under_mle),
+    tolerance = 1e-3
+  )))
+})
+
+test_that("a Laplace scale is fit on the residuals that are there", {
+  exposure <- continuous_laplace_data$exposure
+  exposure[c(5, 22)] <- NA
+  residuals <- exposure - continuous_laplace_data$mu
+  scale <- mean(abs(residuals), na.rm = TRUE)
+
+  weights <- continuous_laplace_wt(
+    exposure = exposure,
+    .density = dens_laplace()
+  )
+
+  # A unit with no exposure has no residual for the likelihood to read and no
+  # weight to be given, and it is the only unit that comes back missing.
+  expect_equal(
+    as.numeric(weights)[-c(5, 22)],
+    (1 / laplace_density(residuals, scale))[-c(5, 22)],
+    tolerance = 1e-12
+  )
+  expect_true(all(is.na(as.numeric(weights)[c(5, 22)])))
+})
+
+test_that("weights record the Laplace scale estimator they were built at", {
+  under_mle <- continuous_laplace_wt(.density = dens_laplace())
+  under_rms <- continuous_laplace_wt(
+    .density = dens_laplace(sigma_method = "rms")
+  )
+
+  # A scale estimated under the family is a source of its own, the one `ipw()`
+  # has an equation of its own to solve. It is estimated from the residuals, as
+  # the pooled spread is, so there is no constant to keep either.
+  expect_identical(density_meta(under_mle)$sigma, "mle")
+  expect_null(density_meta(under_mle)$sigma_value)
+  expect_identical(density_meta(under_rms)$sigma, "pooled")
+
+  # The string names the family at its default, which is the constructor at its
+  # default.
+  by_string <- continuous_laplace_wt(.density = "laplace")
+
+  expect_identical(density_meta(by_string)$sigma, "mle")
+  expect_equal(
+    as.numeric(by_string),
+    as.numeric(under_mle),
+    tolerance = 1e-12
+  )
+})
+
+test_that("a spread supplied and a Laplace scale estimated under itself are refused together", {
+  # `.sigma` says the spread is a number of the caller's own, and a family that
+  # estimates its scale under itself says it is read off the residuals: two
+  # instructions about the same quantity, one of which would go unread whichever
+  # way the pairing was resolved. It is the pairing that is refused rather than
+  # the family, so the moment estimator takes a supplied spread as it always
+  # has.
+  expect_error(
+    continuous_laplace_wt(.density = dens_laplace(), .sigma = 0.9),
+    class = "propensity_density_error"
+  )
+  expect_error(
+    continuous_laplace_wt(.density = "laplace", .sigma = 0.9),
+    class = "propensity_density_error"
+  )
+  expect_no_error(
+    continuous_laplace_wt(
+      .density = dens_laplace(sigma_method = "rms"),
+      .sigma = 0.9
+    )
+  )
+
+  # The check keys on the field rather than on the family, so it reads a Laplace
+  # exactly as it reads a t.
+  expect_error(
+    check_sigma_method(0.9, dens_laplace()),
+    class = "propensity_density_error"
+  )
+  expect_null(check_sigma_method(0.9, dens_laplace(sigma_method = "rms")))
+  expect_null(check_sigma_method(NULL, dens_laplace()))
+})
+
 # ---- a numerator model supplied to stabilize --------------------------------
 
 # A conditional numerator for the same problem: the exposure fit on a baseline

From cf13d68b47ff964bced0240ff7f1e2db668aa8ce Mon Sep 17 00:00:00 2001
From: Malcolm Barrett <malcolmbarrett@gmail.com>
Date: Wed, 16 Sep 2026 11:10:50 -0400
Subject: [PATCH 2/7] estimate the Laplace scale under the Laplace through a
 shared scale seam

---
 NEWS.md                                       |  56 +++---
 R/continuous-models.R                         |   9 +-
 R/density.R                                   | 168 +++++++++++++-----
 R/ipw-continuous-models.R                     |  15 +-
 R/ipw-joint-models.R                          |  20 +--
 R/ipw-mestimation.R                           |   8 +-
 R/ipw-psi.R                                   |  72 +++++---
 R/ipw.R                                       |  27 +--
 R/psw.R                                       |   7 +-
 R/weights.R                                   |  82 +++++----
 man/dens_normal.Rd                            |  65 ++++---
 man/exposure_type.Rd                          |   7 +-
 man/ipw-methods.Rd                            |  27 +--
 man/wt_ate.Rd                                 |  72 ++++----
 tests/testthat/_snaps/ipw-joint-continuous.md |   4 +-
 tests/testthat/_snaps/weights-continuous.md   |   6 +-
 tests/testthat/test-ipw-continuous.R          |  27 ++-
 tests/testthat/test-ipw-psi.R                 |  10 +-
 tests/testthat/test-weights-continuous.R      | 110 +++++++++---
 19 files changed, 501 insertions(+), 291 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index 752f29cd..93291099 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -290,20 +290,20 @@
   freedom, and the root mean square is pulled outward by the large residuals a
   heavy tail produces, which is what the family was chosen to accommodate.
   `sigma_method = "mle"` now estimates the scale under the t itself, by the
-  score each residual enters through a bounded term of, so a residual far out in
-  the tail moves the estimate by less than it moves the root mean square. It
-  spreads the conditional density alone; the marginal density that stabilizes
-  the weights is still read at the exposure's own mean and root mean square.
-  Weights built that way record `sigma = "mle"` in `density_meta()`, and `ipw()`
-  stacks the score of the t for the scale in place of the moment equation the
-  pooled spread is the root of, so the sandwich accounts for having estimated
-  it. The default is unchanged, and a `.sigma` supplied alongside
-  `sigma_method = "mle"` is refused with an error of class
-  `propensity_density_error`, being a second instruction about the same
-  quantity. Residuals a model reproduced exactly say nothing about the spread
-  of the density around it, so a fit with enough of them that the likelihood
-  has no maximum at a positive scale is refused with that same class, naming
-  how many they are.
+  score each residual enters through a bounded term of, so a residual far out
+  in the tail moves the estimate by less than it moves the root mean square. It
+  spreads both densities of the ratio: the marginal density that stabilizes the
+  weights is read at the exposure's own mean and at the scale the same
+  estimator gives for it. Weights built that way record `sigma = "mle"` in
+  `density_meta()`, and `ipw()` stacks that score for each of the two scales in
+  place of the moment equation the pooled spread is the root of, so the
+  sandwich accounts for having estimated them. The default is unchanged, and a
+  `.sigma` supplied alongside `sigma_method = "mle"` is refused with an error
+  of class `propensity_density_error`, being a second instruction about the
+  same quantity. Residuals a model reproduced exactly say nothing about the
+  spread of the density around it, so a fit with enough of them that the
+  likelihood has no maximum at a positive scale is refused with that same
+  class, naming how many they are.
 
 * Weights for a continuous exposure now print the record of the density ratio
   they are under their values, as the three lines `density_meta()` renders. Two
@@ -686,23 +686,25 @@
 
   Both densities in the ratio are now evaluated on a standardized residual,
   `(A - mu) / sigma` for the conditional density and the exposure standardized
-  by its own mean and standard deviation for the marginal one, and each is
-  divided by the spread that standardized it. That factor is the Jacobian of
-  the change of variable, and it returns both densities to the exposure's own
-  units, so every family is read on one scale and the normal family returns the
-  weights the package returned before, to within a rounding error in the last
-  binary digit. Whatever the density gives back is checked before it becomes a
-  weight: one finite, non-negative value for each standardized residual, and
-  not zero at every one of them. Anything else is an error of class
-  `propensity_density_error` whose message reports the standardized residuals
-  it failed at. What the ratio was built from is recorded on the weights and
-  read back with `density_meta()`.
+  by its own mean and the spread its family estimates for the marginal one, and
+  each is divided by the spread that standardized it. That factor is the
+  Jacobian of the change of variable, and it returns both densities to the
+  exposure's own units, so every family is read on one scale and the normal
+  family returns the weights the package returned before, to within a rounding
+  error in the last binary digit. Whatever the density gives back is checked
+  before it becomes a weight: one finite, non-negative value for each
+  standardized residual, and not zero at every one of them. Anything else is an
+  error of class `propensity_density_error` whose message reports the
+  standardized residuals it failed at. What the ratio was built from is
+  recorded on the weights and read back with `density_meta()`.
 
 * `wt_ate()` and `wt_cens()` gain a `numerator` argument, which chooses how the
   marginal density that stabilizes a continuous exposure's weights is arrived
   at. `"marginal"`, the default and the behavior of every earlier version,
-  reads the family `.density` names at the population mean and standard
-  deviation of the exposure; those two moments are parameters of the weights,
+  reads the family `.density` names at the population mean of the exposure and
+  at the spread the same estimator gives for it, which is the standard
+  deviation unless the family estimates a scale of its own. That mean and that
+  spread are parameters of the weights,
   and `ipw()` estimates them alongside the rest of its parameter vector.
   `"integrated"` marginalizes the conditional density numerically instead,
   averaging it over the units at each of 50 points spanning the exposure and
diff --git a/R/continuous-models.R b/R/continuous-models.R
index 1ad61985..4dc0d040 100644
--- a/R/continuous-models.R
+++ b/R/continuous-models.R
@@ -4,8 +4,9 @@
 # center of the conditional density: one fitted conditional mean for each unit,
 # on the scale of the exposure, around which a single spread describes the
 # residuals. Every class here reports that mean through `fitted()`, whatever
-# link it was fit under, and every one of them is spread by the pooled residual
-# root mean square unless the caller supplies `.sigma`.
+# link it was fit under, and every one of them is spread by the estimator its
+# density's family asks for, the pooled residual root mean square unless the
+# caller supplies `.sigma` or the family estimates a scale of its own.
 #
 # Methods are resolved in class order, so a `gam` is read as a `gam` rather than
 # as the `glm` it inherits from, and an `rlm` as an `rlm` rather than as the
@@ -66,8 +67,8 @@ check_continuous_model_response <- function(
 # iteration, so its fitted values are read the way an `lm`'s are. Its spread is
 # not: `rlm` reports a robust scale estimate in `fit$s`, which resists the
 # extreme residuals rather than pooling all of them, while the density these
-# weights are a ratio of is spread by the pooled residual root mean square, as
-# it is for every other class. That is deliberate, and the robust scale is
+# weights are a ratio of is spread by the estimator its family asks for, as it
+# is for every other class. That is deliberate, and the robust scale is
 # available by passing `.sigma = fit$s`.
 #' @export
 extract_continuous_ps.rlm <- function(model, call = rlang::caller_env()) {
diff --git a/R/density.R b/R/density.R
index 30700d20..444ad0a9 100644
--- a/R/density.R
+++ b/R/density.R
@@ -51,10 +51,13 @@ new_density_spec <- function(
 #' * `dens_normal()` is the standard normal density, the default.
 #' * `dens_laplace()` is the standard Laplace density,
 #'   \eqn{\exp(-|z|) / 2}, which puts more mass in the tails than the normal.
+#'   Its scale is estimated under the Laplace itself, where it is the mean
+#'   absolute residual, or by the root mean square of the residuals, which is
+#'   the spread every family without an estimator of its own is read at.
 #' * `dens_t()` is Student's t density with `df` degrees of freedom, heavier
 #'   tailed still, and heavier the smaller `df` is. Its scale is estimated by
-#'   the root mean square of the residuals, as every other family's spread is,
-#'   or by maximum likelihood under the t itself.
+#'   the root mean square of the residuals, or by maximum likelihood under the
+#'   t itself.
 #' * `dens_kernel()` is a kernel density estimate of the standardized
 #'   residuals, fit with [stats::density()] and interpolated to each
 #'   observation. It assumes no family at all, at the cost of a density that
@@ -86,39 +89,49 @@ new_density_spec <- function(
 #'   `family`, the name of the family; `params`, the parameters that identify
 #'   it, which is an empty list for a family that takes none; and `fn`, the
 #'   function that evaluates the density, which is `NULL` for `dens_kernel()`.
-#'   `dens_t()` carries a fourth element, `sigma_method`, the name of the
-#'   estimator its scale is read with.
+#'   `dens_t()` and `dens_laplace()` carry a fourth element, `sigma_method`,
+#'   the name of the estimator their scale is read with.
 #'
-#' @section The spread of a t density:
+#' @section The spread of the conditional density:
 #'
 #' Both densities are evaluated on a residual standardized by a spread, and for
 #' the conditional density that spread is estimated from the residuals of the
-#' propensity score model. Every family takes the root mean square of those
-#' residuals, which is `sigma_method = "rms"`, the default, and is the maximum
-#' likelihood estimator of the spread of a normal density.
+#' propensity score model. A family with no scale estimator of its own takes
+#' the root mean square of those residuals, which is `sigma_method = "rms"` and
+#' is the maximum likelihood estimator of the spread of a normal density.
 #'
-#' It is not the maximum likelihood estimator of the scale of a t. The scale of
-#' a t is smaller than its standard deviation by a factor that depends on `df`,
+#' It is not the maximum likelihood estimator of the scale of a heavier-tailed
+#' family. The scale of such a family is smaller than its standard deviation,
 #' and the root mean square is pulled outward by the large residuals a heavy
 #' tail produces, which is what the family was chosen to accommodate.
-#' `sigma_method = "mle"` estimates the scale under the t itself, as the root of
+#' `sigma_method = "mle"` estimates the scale under the family itself.
+#' `dens_laplace()` takes it by default, where it is the mean absolute
+#' residual. For `dens_t()` it is the root of
 #' \eqn{\sum_i \left[(\nu + 1) r_i^2 / (\nu \sigma^2 + r_i^2)\right] = n},
 #' where \eqn{r_i} is the residual and \eqn{\nu} is `df`. Each residual enters
 #' that sum through a bounded term, so a residual far out in the tail moves the
 #' estimate by less than it moves the root mean square. Prefer it when the
-#' residuals are heavy tailed, which is the case the t family is for; the two
+#' residuals are heavy tailed, which is the case these families are for; the two
 #' estimators answer the same question, and answer it alike, as `df` grows and
 #' the t approaches the normal.
 #'
-#' The choice describes the conditional density alone. The marginal density
-#' that stabilizes the weights is the exposure's own, read at the exposure's
-#' mean and root mean square, whatever the conditional spread was estimated by.
-#' A scale estimated by maximum likelihood is recorded by [density_meta()] as
-#' `sigma = "mle"`, and [ipw()] estimates it alongside the propensity score
-#' model's coefficients, solving the equation above as part of its stacked
-#' system so that the standard errors account for it. Supplying a `.sigma` says
-#' the spread is a number of your own rather than one estimated from the
-#' residuals, so the two cannot be given together.
+#' `dens_normal()` offers no such choice, because the root mean square is
+#' already its maximum likelihood estimator. Neither do `dens_kernel()` and
+#' `dens_fn()`: a kernel is fit to the standardized residuals and divided by the
+#' same spread, so the spread cancels up to the grid the estimate is
+#' interpolated on, and a density you write yourself names no family an
+#' estimator could be derived from.
+#'
+#' The choice describes both densities of the ratio. The marginal density that
+#' stabilizes the weights is the exposure's own, read at the exposure's mean and
+#' at the spread the same estimator gives for it, so the two halves of the ratio
+#' are densities of the same width. A scale estimated by maximum likelihood is
+#' recorded by [density_meta()] as `sigma = "mle"`, and [ipw()] estimates it
+#' alongside the propensity score model's coefficients, solving the equation it
+#' is the root of as part of its stacked system so that the standard errors
+#' account for it. Supplying a `.sigma` says the spread is a number of your own
+#' rather than one estimated from the residuals, so the two cannot be given
+#' together.
 #'
 #' @examples
 #' dens_normal()
@@ -127,6 +140,8 @@ new_density_spec <- function(
 #'
 #' dens_t(df = 4, sigma_method = "mle")
 #'
+#' dens_laplace(sigma_method = "rms")
+#'
 #' dens_kernel(adjust = 1.5)
 #'
 #' dens_fn(function(z) stats::dt(z, df = 4))
@@ -139,8 +154,14 @@ dens_normal <- function() {
 
 #' @rdname dens_normal
 #' @export
-dens_laplace <- function() {
-  new_density_spec("laplace", fn = function(z) exp(-abs(z)) / 2)
+dens_laplace <- function(sigma_method = c("mle", "rms")) {
+  sigma_method <- rlang::arg_match(sigma_method)
+
+  new_density_spec(
+    "laplace",
+    fn = function(z) exp(-abs(z)) / 2,
+    sigma_method = sigma_method
+  )
 }
 
 #' @rdname dens_normal
@@ -418,7 +439,7 @@ check_numerator <- function(
              model estimated at every unit's fitted mean, and {.arg .sigma}
              replaces the spread of that density with one of your own.",
         i = "Leave {.arg .sigma} unset to spread the conditional density by the
-             pooled residual root mean square, or use {.arg numerator} =
+             estimator its family asks for, or use {.arg numerator} =
              {.val marginal}, which takes a spread you supply."
       ),
       error_class = "propensity_numerator_error",
@@ -601,9 +622,9 @@ check_numerator_model_not_multinom <- function(
 
 # What the numerator model contributes to the ratio: the conditional mean it
 # fits for each unit and the spread of its residuals around them. The spread is
-# the uncentered root mean square, the estimator the conditional density is read
-# at, so the two densities of the ratio are spread by the same estimator applied
-# to each model's own residuals.
+# the uncentered root mean square of that model's own residuals, taken about
+# zero rather than about their mean because the fitted mean is what the ratio
+# reads the density at.
 numerator_model_moments <- function(numerator_model) {
   mu <- as.numeric(stats::fitted(numerator_model))
   residuals <- as.numeric(stats::residuals(numerator_model, type = "response"))
@@ -669,12 +690,11 @@ check_continuous_finite <- function(x, arg, call = rlang::caller_env()) {
   )
 }
 
-# The spread of the conditional density: the one the caller supplied, the scale
-# of a t density estimated by maximum likelihood, or the pooled uncentered root
-# mean square of the residuals. The root mean square is uncentered because the
-# residuals of a fitted model already average to zero, and because the
-# estimating equation `ipw()` solves for the same quantity is the uncentered
-# moment.
+# The spread of the conditional density: the one the caller supplied, or the one
+# the residuals leave under the family that reads them. The supplied spread and
+# the estimator are the whole of the resolution, so nothing here can answer a
+# different question from the estimator every other reader of the spread comes
+# through.
 continuous_sigma <- function(
   exposure,
   mu,
@@ -686,18 +706,75 @@ continuous_sigma <- function(
     return(.sigma)
   }
 
-  residuals <- exposure - mu
+  density_scale_estimate(exposure - mu, density, call = call)
+}
 
-  if (density_sigma_is_mle(density)) {
-    return(t_sigma_mle(residuals, density$params$df, call = call))
+# The spread a density is read at, estimated from the residuals it is read
+# against: the uncentered root mean square unless the density asks to be fit
+# under itself, and then that family's own estimator. The root mean square is
+# uncentered because the residuals of a fitted model already average to zero,
+# and because the estimating equation `ipw()` solves for the same quantity is
+# the uncentered moment.
+#
+# Every spread estimated from residuals comes through here, the conditional
+# density's and the marginal one's alike, so the two halves of a ratio cannot be
+# spread by different estimators. What is asked is what the density asked for
+# rather than what family it names: a family with an estimator of its own that
+# was asked for the root mean square is spread by the root mean square on both
+# halves.
+density_scale_estimate <- function(
+  residuals,
+  density,
+  call = rlang::caller_env()
+) {
+  if (!density_sigma_is_mle(density)) {
+    return(sqrt(mean(residuals^2, na.rm = TRUE)))
   }
 
-  sqrt(mean(residuals^2, na.rm = TRUE))
+  switch(
+    density$family,
+    t = t_sigma_mle(residuals, density$params$df, call = call),
+    laplace = laplace_sigma_mle(residuals),
+    abort_no_scale_estimator(density, call = call)
+  )
+}
+
+# The maximum likelihood scale of a Laplace density fit to the residuals that
+# are there, which is their mean absolute value in closed form, so it needs
+# neither a search nor a refusal. It is zero only where every residual is zero,
+# which is where the root mean square is zero too, so the two estimators agree
+# on the degenerate case.
+laplace_sigma_mle <- function(residuals) {
+  mean(abs(residuals), na.rm = TRUE)
+}
+
+# A specification that asks to be fit under itself while naming a family with no
+# estimator of its own. Only `dens_t()` and `dens_laplace()` take
+# `sigma_method`, so nothing a user writes reaches this; a specification built
+# by hand does, and it is refused rather than spread by nothing, which would
+# leave every weight missing and the stacked system solving an equation no
+# weight was built at.
+abort_no_scale_estimator <- function(density, call = rlang::caller_env()) {
+  family <- density$family
+
+  abort(
+    c(
+      "The scale of a {.val {family}} density cannot be estimated by maximum
+       likelihood.",
+      x = "{.code sigma_method = \"mle\"} asks for an estimator this package
+           has only for {.fun dens_t} and {.fun dens_laplace}.",
+      i = "Build the density with {.fun dens_t} or {.fun dens_laplace}, or
+           leave its scale to the root mean square of the residuals."
+    ),
+    error_class = "propensity_density_error",
+    call = call
+  )
 }
 
 # Whether a density asks for its scale to be estimated under itself rather than
-# by the root mean square every family otherwise takes. Only `dens_t()` records
-# the field, so every other specification answers no.
+# by the root mean square every family otherwise takes. Only `dens_t()` and
+# `dens_laplace()` record the field, so every other specification answers no,
+# and so does no density at all.
 density_sigma_is_mle <- function(density) {
   identical(density$sigma_method, "mle")
 }
@@ -805,15 +882,18 @@ check_sigma_method <- function(.sigma, density, call = rlang::caller_env()) {
     return(invisible(NULL))
   }
 
+  family <- density$family
+
   abort(
     c(
       "{.code sigma_method = \"mle\"} cannot be used with {.arg .sigma}.",
       x = "{.code sigma_method = \"mle\"} estimates the scale of the
-           conditional density from the residuals of the propensity score
-           model, and {.arg .sigma} is a spread of your own that replaces it.",
-      i = "Drop {.arg .sigma} to estimate the scale under the t density, or
-           build the density with {.code sigma_method = \"rms\"} to spread the
-           one you supplied."
+           {.val {family}} conditional density from the residuals of the
+           propensity score model, and {.arg .sigma} is a spread of your own
+           that replaces it.",
+      i = "Drop {.arg .sigma} to estimate the scale under the {.val {family}}
+           density, or build the density with {.code sigma_method = \"rms\"}
+           to spread the one you supplied."
     ),
     error_class = "propensity_density_error",
     call = call
diff --git a/R/ipw-continuous-models.R b/R/ipw-continuous-models.R
index 5c350dc4..2525039a 100644
--- a/R/ipw-continuous-models.R
+++ b/R/ipw-continuous-models.R
@@ -1258,13 +1258,14 @@ check_ipw_numerator_model <- function(
 }
 
 # The spread the stacked system reads the conditional density at. A pooled
-# spread is the residual moment the system estimates alongside the coefficients,
-# and a scale fit by maximum likelihood is the score of the t estimated in the
-# same place; a single spread the caller supplied is a constant it holds fixed,
-# carrying none of its uncertainty, which is what fixing it says. A spread
-# supplied for each observation is neither: it is a function of the data that no
-# parameter value here reproduces, so it is refused before anything is solved
-# rather than reported afterwards as two vectors that disagree.
+# spread is the residual moment the system estimates alongside the
+# coefficients, and a scale fit by maximum likelihood is the score of the
+# family it names, estimated in the same place; a single spread the caller
+# supplied is a constant it holds fixed, carrying none of its uncertainty,
+# which is what fixing it says. A spread supplied for each observation is
+# neither: it is a function of the data that no parameter value here
+# reproduces, so it is refused before anything is solved rather than reported
+# afterwards as two vectors that disagree.
 ipw_continuous_spread <- function(meta, call = rlang::caller_env()) {
   if (identical(meta$sigma, "pooled")) {
     return(list(kind = "pooled", value = NULL))
diff --git a/R/ipw-joint-models.R b/R/ipw-joint-models.R
index 112038ea..e40cbc00 100644
--- a/R/ipw-joint-models.R
+++ b/R/ipw-joint-models.R
@@ -2272,7 +2272,6 @@ ipw_joint_models_score_rows <- function(
   link,
   ps_fns = NULL,
   sigma_row = TRUE,
-  sigma = NULL,
   density = NULL
 ) {
   if (identical(block$type, "categorical")) {
@@ -2297,7 +2296,7 @@ ipw_joint_models_score_rows <- function(
 
   rbind(
     score,
-    ipw_continuous_sigma_row(sigma, y - block$ps, block$sigma2_d, density)
+    density_scale_row(y - block$ps, block$sigma2_d, density)
   )
 }
 
@@ -2324,7 +2323,12 @@ ipw_joint_models_stab_rows <- function(spec, th_stab, exposures) {
 
     if (identical(stab$numerator, "marginal")) {
       if (identical(stab$type, "continuous")) {
-        return(deli::ee_mean_variance(th, y = a))
+        r_a <- a - th[[1]]
+
+        return(rbind(
+          matrix(r_a, nrow = 1),
+          density_scale_row(r_a, th[[2]], spec$density)
+        ))
       }
 
       # A categorical component's exposure is its reference-first indicator
@@ -2439,12 +2443,7 @@ ipw_init_joint_models_ps <- function(spec, call = rlang::caller_env()) {
     c(
       alpha,
       stats::setNames(
-        ipw_continuous_sigma2_seed(
-          spec$sigma,
-          resid,
-          spec$density,
-          call = call
-        ),
+        ipw_continuous_sigma2_seed(resid, spec$density, call = call),
         paste0("sigma2_", spec$names[[i]])
       )
     )
@@ -2486,7 +2485,7 @@ ipw_init_joint_models_stab <- function(spec) {
       if (identical(stab$type, "continuous")) {
         mu_a <- mean(a_fit)
         return(stats::setNames(
-          c(mu_a, mean((a_fit - mu_a)^2)),
+          c(mu_a, density_scale_estimate(a_fit - mu_a, spec$density)^2),
           paste0(prefix, c("mu", "sigma2"))
         ))
       }
@@ -2627,7 +2626,6 @@ ipw_psi_joint_models <- function(
           ps_link[[i]],
           ps_fns = ps_fns,
           sigma_row = sigma_row,
-          sigma = spec$sigma,
           density = spec$density
         )
       })
diff --git a/R/ipw-mestimation.R b/R/ipw-mestimation.R
index 3bb84404..189123b9 100644
--- a/R/ipw-mestimation.R
+++ b/R/ipw-mestimation.R
@@ -2865,15 +2865,15 @@ ipw_compare_weights <- function(
     )
     # A dose is the only exposure whose weights carry a spread, and it is a dose
     # whether it is the exposure or the second component of a joint one. The
-    # stacked system estimates one pooled residual variance, so weights built on
-    # observation-level standard deviations are a different function of the data
-    # and cannot be reproduced here at any parameter value.
+    # stacked system estimates a single spread, the one its family estimates, so
+    # weights built on observation-level standard deviations are a different
+    # function of the data and cannot be reproduced here at any parameter value.
     dose <- identical(exposure_type, "continuous") ||
       "continuous" %in% components
     sigma_hint <- if (dose) {
       "Weights built with an observation-level {.arg .sigma}, such as \\
       {.code influence(model)$sigma}, are one cause: {.fun ipw} models the \\
-      conditional density with a single pooled residual root mean square, \\
+      conditional density with a single spread, the one its family estimates, \\
       which is what {.fun wt_ate} uses when no {.arg .sigma} is given."
     } else {
       NULL
diff --git a/R/ipw-psi.R b/R/ipw-psi.R
index 7f6a4991..93e42b79 100644
--- a/R/ipw-psi.R
+++ b/R/ipw-psi.R
@@ -744,7 +744,6 @@ ipw_init_continuous <- function(spec, call = rlang::caller_env()) {
     c(
       alpha,
       sigma2_d = ipw_continuous_sigma2_seed(
-        spec$sigma,
         ipw_continuous_fit_residuals(spec),
         spec$density,
         call = call
@@ -779,7 +778,7 @@ ipw_init_continuous <- function(spec, call = rlang::caller_env()) {
   } else if (identical(spec$numerator, "marginal")) {
     a_fit <- ipw_continuous_fit_exposure(spec)
     mu_a <- mean(a_fit)
-    sigma2_a <- mean((a_fit - mu_a)^2)
+    sigma2_a <- density_scale_estimate(a_fit - mu_a, spec$density)^2
     stab_block <- c(mu_a = mu_a, sigma2_a = sigma2_a)
   } else {
     stab_block <- numeric(0)
@@ -1248,11 +1247,11 @@ ipw_psi_continuous <- function(
     } else {
       rbind(
         ps_score,
-        ipw_continuous_sigma_row(
-          spec$sigma,
+        density_scale_row(
           a - inputs$mu,
           inputs$extras$sigma2_d,
-          spec$density
+          spec$density,
+          call = call
         )
       )
     }
@@ -1276,7 +1275,12 @@ ipw_psi_continuous <- function(
         matrix((a - inputs$extras$mu_n)^2 - th_stab[[p_n + 1L]], nrow = 1)
       )
     } else {
-      deli::ee_mean_variance(th_stab, y = a)
+      r_a <- a - th_stab[[1]]
+
+      rbind(
+        matrix(r_a, nrow = 1),
+        density_scale_row(r_a, th_stab[[2]], spec$density, call = call)
+      )
     }
 
     w <- weight_fn(inputs$mu, a, inputs$extras)
@@ -1374,24 +1378,41 @@ ipw_numerator_grid <- function(exposure, numerator) {
   seq(min(present), max(present), length.out = continuous_grid_n)
 }
 
-# The equation the stacked system estimates the conditional variance of the
-# density by, which is the equation the spread the weights were built with is
-# the root of. The pooled spread is the root mean square of the residuals, so
-# its row is the uncentered second moment; a scale fit by maximum likelihood is
-# the root of the score of the t itself for the scale, so its row is that score,
-# multiplied through by the scale so that each residual enters through a bounded
-# term.
-ipw_continuous_sigma_row <- function(sigma, residuals, sigma2_d, density) {
-  if (identical(sigma$kind, "mle")) {
-    df <- density$params$df
-
-    return(matrix(
-      (df + 1) * residuals^2 / (df * sigma2_d + residuals^2) - 1,
-      nrow = 1
-    ))
+# The equation the stacked system estimates a density's squared scale by, which
+# is the equation the spread the weights were built at is the root of. The
+# pooled spread is the root mean square of the residuals, so its row is the
+# uncentered second moment; a scale fit under the family itself is the root of
+# that family's own score for the scale, multiplied through by the scale so that
+# each residual enters through a bounded term.
+#
+# It asks the question `density_scale_estimate()` asks, so the estimator and the
+# equation it is the root of cannot come apart: the row follows what the density
+# asked for rather than the family it names, and a family asked for the root
+# mean square emits the moment row. The conditional density's residuals are
+# `a - mu`, a marginal numerator's are `a - mu_a`, and a numerator model's are
+# `a - mu_n`.
+density_scale_row <- function(
+  residuals,
+  sigma2,
+  density,
+  call = rlang::caller_env()
+) {
+  if (!density_sigma_is_mle(density)) {
+    return(matrix(residuals^2 - sigma2, nrow = 1))
   }
 
-  matrix(residuals^2 - sigma2_d, nrow = 1)
+  switch(
+    density$family,
+    t = matrix(
+      (density$params$df + 1) *
+        residuals^2 /
+        (density$params$df * sigma2 + residuals^2) -
+        1,
+      nrow = 1
+    ),
+    laplace = matrix(abs(residuals) - sqrt(sigma2), nrow = 1),
+    abort_no_scale_estimator(density, call = call)
+  )
 }
 
 # The seed for that parameter: the exact root of whichever row estimates it, so
@@ -1401,14 +1422,9 @@ ipw_continuous_sigma_row <- function(sigma, residuals, sigma2_d, density) {
 # over; the rows the system goes on to solve over are a separate question the
 # row itself answers.
 ipw_continuous_sigma2_seed <- function(
-  sigma,
   residuals,
   density,
   call = rlang::caller_env()
 ) {
-  if (identical(sigma$kind, "mle")) {
-    return(t_sigma_mle(residuals, density$params$df, call = call)^2)
-  }
-
-  mean(residuals^2)
+  density_scale_estimate(residuals, density, call = call)^2
 }
diff --git a/R/ipw.R b/R/ipw.R
index 734fe1e9..d7c38c79 100644
--- a/R/ipw.R
+++ b/R/ipw.R
@@ -986,8 +986,8 @@
 #' as a known constant. It solves an equation of its own that the stacked system
 #' does not write, so none of its sampling variability is propagated, and the
 #' standard errors are those of a psi score read at a scale treated as fixed.
-#' The spread of the conditional density is a separate quantity, and is the
-#' pooled residual spread there as it is for every other class.
+#' The spread of the conditional density is a separate quantity, and comes from
+#' the estimator its family asks for, there as for every other class.
 #'
 #' [MASS::psi.bisquare()] and [MASS::psi.hampel()] redescend, and that is worth
 #' knowing about the standard errors they get. A redescending psi's estimating
@@ -1341,17 +1341,18 @@
 #' propensity score model; a mismatch errors, on both standard error methods.
 #'
 #' For a continuous exposure that requirement bears on the spread of the
-#' conditional density. `ipw()` stacks a single pooled residual variance
-#' alongside the propensity score coefficients, so the weights it rebuilds are
-#' the ones [wt_ate()] produces with its pooled default. A single `.sigma` is
-#' taken instead as a known constant: the weights are rebuilt at the number that
-#' was supplied, and the stacked system carries none of that number's
-#' uncertainty, which is what fixing a spread says. Weights built with an
-#' observation-level `.sigma`, such as `influence(model)$sigma`, are a different
-#' function of the data with no counterpart in the stacked system, and are
-#' refused before anything is solved, with an error of class
-#' `propensity_ipw_sigma_error`. Rebuild the weights with the pooled default or
-#' with one number to use `ipw()`.
+#' conditional density. `ipw()` stacks a single conditional spread alongside
+#' the propensity score coefficients, by the moment the pooled spread is the
+#' root of or, for a density built with `sigma_method = "mle"`, by the score of
+#' the family it names, so the weights it rebuilds are the ones [wt_ate()]
+#' produces when no `.sigma` is given. A single `.sigma` is taken instead as a
+#' known constant: the weights are rebuilt at the number that was supplied, and
+#' the stacked system carries none of that number's uncertainty, which is what
+#' fixing a spread says. Weights built with an observation-level `.sigma`, such
+#' as `influence(model)$sigma`, are a different function of the data with no
+#' counterpart in the stacked system, and are refused before anything is
+#' solved, with an error of class `propensity_ipw_sigma_error`. Rebuild the
+#' weights with the pooled default or with one number to use `ipw()`.
 #'
 #' It bears on a `stabilization_score` in the same way. A score written per
 #' observation is one value per unit, so it does not survive the rows being
diff --git a/R/psw.R b/R/psw.R
index 1c91e139..19022874 100644
--- a/R/psw.R
+++ b/R/psw.R
@@ -335,9 +335,10 @@ stabilization_score <- function(wt) {
 #'   to different values of the same terms, record different numerators and are
 #'   combined the way any other disagreeing records are.
 #' * `sigma`, where the residual spread of the conditional density came from:
-#'   `"pooled"` for the pooled residual root mean square, `"mle"` for the scale
-#'   a [dens_t()] built with `sigma_method = "mle"` estimates under the t
-#'   itself, and `"supplied"` for a `.sigma` the caller gave.
+#'   `"pooled"` for the pooled residual root mean square, `"mle"` for a scale
+#'   estimated under the family that reads it, which [dens_t()] and
+#'   [dens_laplace()] take with `sigma_method = "mle"`, and `"supplied"` for a
+#'   `.sigma` the caller gave.
 #' * `sigma_value`, the single spread the caller supplied, and `NULL` for a
 #'   spread estimated from the residuals, by either estimator, and for one
 #'   supplied per observation. A spread that is one
diff --git a/R/weights.R b/R/weights.R
index e5e484c5..ea67b6ee 100644
--- a/R/weights.R
+++ b/R/weights.R
@@ -63,9 +63,9 @@
 #' model of the exposure reads, reducing variance. When no `stabilization_score` is
 #' supplied, that estimate is the marginal mean of `.exposure` for a binary or
 #' categorical exposure, and for a continuous exposure the marginal density of
-#' the family `.density` names, evaluated at the population mean and standard
-#' deviation of `.exposure`. Stabilization is supported for ATE and censoring
-#' weights (`wt_ate()` and `wt_cens()`) alone.
+#' the family `.density` names, evaluated at the population mean of `.exposure`
+#' and at the spread the same estimator gives for it. Stabilization is
+#' supported for ATE and censoring weights (`wt_ate()` and `wt_cens()`) alone.
 #'
 #' Whether it is applied is read from the exposure type unless `stabilize` says
 #' so outright. A continuous exposure is stabilized, since the unstabilized
@@ -75,17 +75,17 @@
 #'
 #' For a continuous exposure, `numerator` chooses how that marginal density is
 #' arrived at. `"marginal"`, the default, reads the family `.density` names at
-#' the mean and standard deviation of `.exposure` over the rows the propensity
-#' model kept, which are the rows the weights describe: a model that dropped
-#' rows leaves `.propensity` missing at them, and those rows carry no weight.
-#' Those two moments are parameters of the weights, and [ipw()] estimates them
-#' alongside the rest
-#' of its parameter vector, so the standard errors account for the numerator
-#' having been estimated. `"integrated"` marginalizes the conditional density
-#' numerically instead: it averages \eqn{f_{A|X}(t \mid X_i)} over the units at
-#' each of 50 points spanning `.exposure`, then interpolates that average back
-#' to each observed exposure with a cubic spline. It estimates no parameters of
-#' its own, and it is what WeightIt has done since version 2.0.0.
+#' the mean of `.exposure` and the spread the same estimator gives for it, both
+#' over the rows the propensity model kept, which are the rows the weights
+#' describe: a model that dropped rows leaves `.propensity` missing at them,
+#' and those rows carry no weight. Those two moments are parameters of the
+#' weights, and [ipw()] estimates them alongside the rest of its parameter
+#' vector, so the standard errors account for the numerator having been
+#' estimated. `"integrated"` marginalizes the conditional density numerically
+#' instead: it averages \eqn{f_{A|X}(t \mid X_i)} over the units at each of 50
+#' points spanning `.exposure`, then interpolates that average back to each
+#' observed exposure with a cubic spline. It estimates no parameters of its
+#' own, and it is what WeightIt has done since version 2.0.0.
 #'
 #' The two agree whenever the family is normal and the fitted conditional means
 #' are themselves normal, since an average of normal densities centered on
@@ -294,14 +294,15 @@
 #' is evaluated at \eqn{z_i = (A_i - \hat{A}_i) / \sigma}, where
 #' \eqn{\hat{A}_i} is the fitted conditional mean in `.propensity` and
 #' \eqn{\sigma} is the residual spread; the marginal density is evaluated at
-#' \eqn{z^A_i = (A_i - \bar{A}) / s_A}, where \eqn{\bar{A}} and \eqn{s_A} are
-#' the mean and standard deviation of `.exposure` over the rows the propensity
-#' model kept. Both moments are read over those rows whether the spread was
-#' pooled from the model's own residuals or supplied through `.sigma`, so the
-#' two halves of the ratio describe one set of units. Each density is
-#' then divided by the spread that standardized it, the Jacobian of that change
-#' of variable, which returns both to the exposure's own units so that each
-#' integrates to one:
+#' \eqn{z^A_i = (A_i - \bar{A}) / s_A}, where \eqn{\bar{A}} is the mean of
+#' `.exposure` over the rows the propensity model kept and \eqn{s_A} is the
+#' spread the same estimator gives for it, the standard deviation unless the
+#' family estimates a scale of its own. Both moments are read over those rows
+#' whether the spread was pooled from the model's own residuals or supplied
+#' through `.sigma`, so the two halves of the ratio describe one set of units.
+#' Each density is then divided by the spread that standardized it, the
+#' Jacobian of that change of variable, which returns both to the exposure's
+#' own units so that each integrates to one:
 #'
 #' \deqn{w_i = \frac{g(z^A_i) / s_A}{g(z_i) / \sigma}}
 #'
@@ -355,13 +356,14 @@
 #' on, a linear probability model included.
 #'
 #' \eqn{\sigma} is the pooled residual spread
-#' \eqn{\sqrt{\mathrm{mean}((A - \hat{A})^2)}} unless `.sigma` supplies a single
-#' standard deviation, or one for each unit, which the model methods do not do
-#' on their own. [dens_t()] offers a second estimator, the maximum likelihood
-#' scale of the t itself, through `sigma_method = "mle"`; see **The spread of a
-#' t density** in its documentation. The marginal density is spread by
-#' \eqn{s_A} whatever the conditional spread was arrived at, the spread of the
-#' exposure itself being no business of the conditional model's.
+#' \eqn{\sqrt{\mathrm{mean}((A - \hat{A})^2)}} unless `.sigma` supplies a
+#' single standard deviation, or one for each unit, which the model methods do
+#' not do on their own. [dens_t()] and [dens_laplace()] offer a second
+#' estimator, the maximum likelihood scale of the family itself, through
+#' `sigma_method = "mle"`; see **The spread of the conditional density** in
+#' [dens_t()]'s documentation. A family that estimates its own scale spreads
+#' both densities of the ratio, so the marginal density is read at that scale
+#' too and the two halves are densities of one width rather than two.
 #'
 #' Every model class is spread that same way, [MASS::rlm()] included. `rlm`
 #' reports a robust scale estimate of its own in `fit$s`, which resists the
@@ -376,7 +378,7 @@
 #' [ipw()] models the conditional density with a single conditional variance,
 #' estimated jointly with the rest of the parameter vector, by the moment the
 #' pooled spread is the root of or, for a density built with
-#' `sigma_method = "mle"`, by the score of the t. A single
+#' `sigma_method = "mle"`, by the score of the family it names. A single
 #' `.sigma` is taken instead as a known constant: the weights are rebuilt at the
 #' number that was supplied, and the stacked system carries none of that
 #' number's uncertainty. An observation-level `.sigma` is a different function
@@ -491,8 +493,9 @@
 #'   exposures: a single standard deviation applied to every unit, or one for
 #'   each observation (e.g., `influence(model)$sigma`). Optional: with none
 #'   supplied, including when `.propensity` is a fitted model, the conditional
-#'   density uses the pooled residual spread of `.exposure` around
-#'   `.propensity`.
+#'   density uses the spread its family estimates from the residuals of
+#'   `.exposure` around `.propensity`, the pooled root mean square unless the
+#'   family estimates a scale of its own.
 #'
 #'   The two shapes differ downstream. A single spread is a constant the weights
 #'   can be rebuilt from, so it is recorded in `density_meta()` as `sigma_value`
@@ -532,8 +535,9 @@
 #' @param numerator How the marginal density that stabilizes a continuous
 #'   exposure's weights is obtained, described under **Stabilization** in
 #'   Details. Either `"marginal"`, the default, which reads the family
-#'   `.density` names at the mean and standard deviation of `.exposure` over
-#'   the rows the propensity model kept, or `"integrated"`, which averages the
+#'   `.density` names at the mean of `.exposure` and the spread the same
+#'   estimator gives for it, both over the rows the propensity model kept, or
+#'   `"integrated"`, which averages the
 #'   conditional density
 #'   over the units on a grid spanning `.exposure` and interpolates the result
 #'   back to each observed exposure.
@@ -1376,13 +1380,21 @@ ate_continuous <- function(
   # population the weights describe. Reading the marginal moments over them
   # while the conditional spread above is pooled over the rows that remain would
   # leave the two halves of the ratio describing different sets of units.
+  #
+  # Its spread is read by the estimator the family asks for, which is the one
+  # the conditional density above was spread by, so the two halves of the ratio
+  # are densities of the same width.
   mu_a <- NULL
   sigma_a <- NULL
   if (identical(numerator, "marginal")) {
     kept <- !is.na(.exposure) & !is.na(.propensity)
     exposure_kept <- .exposure[kept]
     mu_a <- mean(exposure_kept)
-    sigma_a <- sqrt(mean((exposure_kept - mu_a)^2))
+    sigma_a <- density_scale_estimate(
+      exposure_kept - mu_a,
+      .density,
+      call = call
+    )
   }
 
   # Whether the weights this call is about to build have a finite second moment.
diff --git a/man/dens_normal.Rd b/man/dens_normal.Rd
index c1a08540..47833c6c 100644
--- a/man/dens_normal.Rd
+++ b/man/dens_normal.Rd
@@ -10,7 +10,7 @@
 \usage{
 dens_normal()
 
-dens_laplace()
+dens_laplace(sigma_method = c("mle", "rms"))
 
 dens_t(df, sigma_method = c("rms", "mle"))
 
@@ -19,13 +19,13 @@ dens_kernel(bw = "nrd0", adjust = 1, kernel = "gaussian", n = 512)
 dens_fn(f)
 }
 \arguments{
-\item{df}{Degrees of freedom for Student's t, a single positive, finite
-number.}
-
 \item{sigma_method}{How the spread of the conditional density is estimated
 from the residuals of the propensity score model, either \code{"rms"} or
 \code{"mle"}. See the section below.}
 
+\item{df}{Degrees of freedom for Student's t, a single positive, finite
+number.}
+
 \item{bw}{The bandwidth passed to \code{\link[stats:density]{stats::density()}}: a single positive
 number, or the name of one of its selection rules
 (\code{"nrd0"}, \code{"nrd"}, \code{"ucv"}, \code{"bcv"}, \code{"SJ"}, \code{"SJ-ste"}, or \code{"SJ-dpi"}).}
@@ -48,8 +48,8 @@ An object of class \code{propensity_density}: a list with the elements
 \code{family}, the name of the family; \code{params}, the parameters that identify
 it, which is an empty list for a family that takes none; and \code{fn}, the
 function that evaluates the density, which is \code{NULL} for \code{dens_kernel()}.
-\code{dens_t()} carries a fourth element, \code{sigma_method}, the name of the
-estimator its scale is read with.
+\code{dens_t()} and \code{dens_laplace()} carry a fourth element, \code{sigma_method},
+the name of the estimator their scale is read with.
 }
 \description{
 Weights for a continuous exposure are a ratio of densities. Both densities
@@ -68,10 +68,13 @@ is wrapped with \code{dens_fn()}.
 \item \code{dens_normal()} is the standard normal density, the default.
 \item \code{dens_laplace()} is the standard Laplace density,
 \eqn{\exp(-|z|) / 2}, which puts more mass in the tails than the normal.
+Its scale is estimated under the Laplace itself, where it is the mean
+absolute residual, or by the root mean square of the residuals, which is
+the spread every family without an estimator of its own is read at.
 \item \code{dens_t()} is Student's t density with \code{df} degrees of freedom, heavier
 tailed still, and heavier the smaller \code{df} is. Its scale is estimated by
-the root mean square of the residuals, as every other family's spread is,
-or by maximum likelihood under the t itself.
+the root mean square of the residuals, or by maximum likelihood under the
+t itself.
 \item \code{dens_kernel()} is a kernel density estimate of the standardized
 residuals, fit with \code{\link[stats:density]{stats::density()}} and interpolated to each
 observation. It assumes no family at all, at the cost of a density that
@@ -82,37 +85,47 @@ estimate.
 \item \code{dens_fn()} is a density you write yourself.
 }
 }
-\section{The spread of a t density}{
+\section{The spread of the conditional density}{
 
 
 Both densities are evaluated on a residual standardized by a spread, and for
 the conditional density that spread is estimated from the residuals of the
-propensity score model. Every family takes the root mean square of those
-residuals, which is \code{sigma_method = "rms"}, the default, and is the maximum
-likelihood estimator of the spread of a normal density.
+propensity score model. A family with no scale estimator of its own takes
+the root mean square of those residuals, which is \code{sigma_method = "rms"} and
+is the maximum likelihood estimator of the spread of a normal density.
 
-It is not the maximum likelihood estimator of the scale of a t. The scale of
-a t is smaller than its standard deviation by a factor that depends on \code{df},
+It is not the maximum likelihood estimator of the scale of a heavier-tailed
+family. The scale of such a family is smaller than its standard deviation,
 and the root mean square is pulled outward by the large residuals a heavy
 tail produces, which is what the family was chosen to accommodate.
-\code{sigma_method = "mle"} estimates the scale under the t itself, as the root of
+\code{sigma_method = "mle"} estimates the scale under the family itself.
+\code{dens_laplace()} takes it by default, where it is the mean absolute
+residual. For \code{dens_t()} it is the root of
 \eqn{\sum_i \left[(\nu + 1) r_i^2 / (\nu \sigma^2 + r_i^2)\right] = n},
 where \eqn{r_i} is the residual and \eqn{\nu} is \code{df}. Each residual enters
 that sum through a bounded term, so a residual far out in the tail moves the
 estimate by less than it moves the root mean square. Prefer it when the
-residuals are heavy tailed, which is the case the t family is for; the two
+residuals are heavy tailed, which is the case these families are for; the two
 estimators answer the same question, and answer it alike, as \code{df} grows and
 the t approaches the normal.
 
-The choice describes the conditional density alone. The marginal density
-that stabilizes the weights is the exposure's own, read at the exposure's
-mean and root mean square, whatever the conditional spread was estimated by.
-A scale estimated by maximum likelihood is recorded by \code{\link[=density_meta]{density_meta()}} as
-\code{sigma = "mle"}, and \code{\link[=ipw]{ipw()}} estimates it alongside the propensity score
-model's coefficients, solving the equation above as part of its stacked
-system so that the standard errors account for it. Supplying a \code{.sigma} says
-the spread is a number of your own rather than one estimated from the
-residuals, so the two cannot be given together.
+\code{dens_normal()} offers no such choice, because the root mean square is
+already its maximum likelihood estimator. Neither do \code{dens_kernel()} and
+\code{dens_fn()}: a kernel is fit to the standardized residuals and divided by the
+same spread, so the spread cancels up to the grid the estimate is
+interpolated on, and a density you write yourself names no family an
+estimator could be derived from.
+
+The choice describes both densities of the ratio. The marginal density that
+stabilizes the weights is the exposure's own, read at the exposure's mean and
+at the spread the same estimator gives for it, so the two halves of the ratio
+are densities of the same width. A scale estimated by maximum likelihood is
+recorded by \code{\link[=density_meta]{density_meta()}} as \code{sigma = "mle"}, and \code{\link[=ipw]{ipw()}} estimates it
+alongside the propensity score model's coefficients, solving the equation it
+is the root of as part of its stacked system so that the standard errors
+account for it. Supplying a \code{.sigma} says the spread is a number of your own
+rather than one estimated from the residuals, so the two cannot be given
+together.
 }
 
 \examples{
@@ -122,6 +135,8 @@ dens_t(df = 4)
 
 dens_t(df = 4, sigma_method = "mle")
 
+dens_laplace(sigma_method = "rms")
+
 dens_kernel(adjust = 1.5)
 
 dens_fn(function(z) stats::dt(z, df = 4))
diff --git a/man/exposure_type.Rd b/man/exposure_type.Rd
index 9a30ae34..9517410b 100644
--- a/man/exposure_type.Rd
+++ b/man/exposure_type.Rd
@@ -71,9 +71,10 @@ weights stabilized on models that read different terms, or that were fit
 to different values of the same terms, record different numerators and are
 combined the way any other disagreeing records are.
 \item \code{sigma}, where the residual spread of the conditional density came from:
-\code{"pooled"} for the pooled residual root mean square, \code{"mle"} for the scale
-a \code{\link[=dens_t]{dens_t()}} built with \code{sigma_method = "mle"} estimates under the t
-itself, and \code{"supplied"} for a \code{.sigma} the caller gave.
+\code{"pooled"} for the pooled residual root mean square, \code{"mle"} for a scale
+estimated under the family that reads it, which \code{\link[=dens_t]{dens_t()}} and
+\code{\link[=dens_laplace]{dens_laplace()}} take with \code{sigma_method = "mle"}, and \code{"supplied"} for a
+\code{.sigma} the caller gave.
 \item \code{sigma_value}, the single spread the caller supplied, and \code{NULL} for a
 spread estimated from the residuals, by either estimator, and for one
 supplied per observation. A spread that is one
diff --git a/man/ipw-methods.Rd b/man/ipw-methods.Rd
index 9cbdd350..07a0a942 100644
--- a/man/ipw-methods.Rd
+++ b/man/ipw-methods.Rd
@@ -1251,8 +1251,8 @@ stacked at the values it was retuned to. That scale enters the location score
 as a known constant. It solves an equation of its own that the stacked system
 does not write, so none of its sampling variability is propagated, and the
 standard errors are those of a psi score read at a scale treated as fixed.
-The spread of the conditional density is a separate quantity, and is the
-pooled residual spread there as it is for every other class.
+The spread of the conditional density is a separate quantity, and comes from
+the estimator its family asks for, there as for every other class.
 
 \code{\link[MASS:psi.bisquare]{MASS::psi.bisquare()}} and \code{\link[MASS:psi.hampel]{MASS::psi.hampel()}} redescend, and that is worth
 knowing about the standard errors they get. A redescending psi's estimating
@@ -1609,17 +1609,18 @@ model. The outcome model weights must match the values implied by the
 propensity score model; a mismatch errors, on both standard error methods.
 
 For a continuous exposure that requirement bears on the spread of the
-conditional density. \code{ipw()} stacks a single pooled residual variance
-alongside the propensity score coefficients, so the weights it rebuilds are
-the ones \code{\link[=wt_ate]{wt_ate()}} produces with its pooled default. A single \code{.sigma} is
-taken instead as a known constant: the weights are rebuilt at the number that
-was supplied, and the stacked system carries none of that number's
-uncertainty, which is what fixing a spread says. Weights built with an
-observation-level \code{.sigma}, such as \code{influence(model)$sigma}, are a different
-function of the data with no counterpart in the stacked system, and are
-refused before anything is solved, with an error of class
-\code{propensity_ipw_sigma_error}. Rebuild the weights with the pooled default or
-with one number to use \code{ipw()}.
+conditional density. \code{ipw()} stacks a single conditional spread alongside
+the propensity score coefficients, by the moment the pooled spread is the
+root of or, for a density built with \code{sigma_method = "mle"}, by the score of
+the family it names, so the weights it rebuilds are the ones \code{\link[=wt_ate]{wt_ate()}}
+produces when no \code{.sigma} is given. A single \code{.sigma} is taken instead as a
+known constant: the weights are rebuilt at the number that was supplied, and
+the stacked system carries none of that number's uncertainty, which is what
+fixing a spread says. Weights built with an observation-level \code{.sigma}, such
+as \code{influence(model)$sigma}, are a different function of the data with no
+counterpart in the stacked system, and are refused before anything is
+solved, with an error of class \code{propensity_ipw_sigma_error}. Rebuild the
+weights with the pooled default or with one number to use \code{ipw()}.
 
 It bears on a \code{stabilization_score} in the same way. A score written per
 observation is one value per unit, so it does not survive the rows being
diff --git a/man/wt_ate.Rd b/man/wt_ate.Rd
index ac650f87..46016c8b 100644
--- a/man/wt_ate.Rd
+++ b/man/wt_ate.Rd
@@ -275,8 +275,9 @@ and are carried through to the weights as missing.}
 exposures: a single standard deviation applied to every unit, or one for
 each observation (e.g., \code{influence(model)$sigma}). Optional: with none
 supplied, including when \code{.propensity} is a fitted model, the conditional
-density uses the pooled residual spread of \code{.exposure} around
-\code{.propensity}.
+density uses the spread its family estimates from the residuals of
+\code{.exposure} around \code{.propensity}, the pooled root mean square unless the
+family estimates a scale of its own.
 
 The two shapes differ downstream. A single spread is a constant the weights
 can be rebuilt from, so it is recorded in \code{density_meta()} as \code{sigma_value}
@@ -418,8 +419,9 @@ the weights were built from is recorded on the result and read back with
 \item{numerator}{How the marginal density that stabilizes a continuous
 exposure's weights is obtained, described under \strong{Stabilization} in
 Details. Either \code{"marginal"}, the default, which reads the family
-\code{.density} names at the mean and standard deviation of \code{.exposure} over
-the rows the propensity model kept, or \code{"integrated"}, which averages the
+\code{.density} names at the mean of \code{.exposure} and the spread the same
+estimator gives for it, both over the rows the propensity model kept, or
+\code{"integrated"}, which averages the
 conditional density
 over the units on a grid spanning \code{.exposure} and interpolates the result
 back to each observed exposure.
@@ -551,9 +553,9 @@ f_A(A) (continuous), or by an estimate of either conditional on variables a
 model of the exposure reads, reducing variance. When no \code{stabilization_score} is
 supplied, that estimate is the marginal mean of \code{.exposure} for a binary or
 categorical exposure, and for a continuous exposure the marginal density of
-the family \code{.density} names, evaluated at the population mean and standard
-deviation of \code{.exposure}. Stabilization is supported for ATE and censoring
-weights (\code{wt_ate()} and \code{wt_cens()}) alone.
+the family \code{.density} names, evaluated at the population mean of \code{.exposure}
+and at the spread the same estimator gives for it. Stabilization is
+supported for ATE and censoring weights (\code{wt_ate()} and \code{wt_cens()}) alone.
 
 Whether it is applied is read from the exposure type unless \code{stabilize} says
 so outright. A continuous exposure is stabilized, since the unstabilized
@@ -563,17 +565,17 @@ and \code{stabilize = TRUE} asks for it there.
 
 For a continuous exposure, \code{numerator} chooses how that marginal density is
 arrived at. \code{"marginal"}, the default, reads the family \code{.density} names at
-the mean and standard deviation of \code{.exposure} over the rows the propensity
-model kept, which are the rows the weights describe: a model that dropped
-rows leaves \code{.propensity} missing at them, and those rows carry no weight.
-Those two moments are parameters of the weights, and \code{\link[=ipw]{ipw()}} estimates them
-alongside the rest
-of its parameter vector, so the standard errors account for the numerator
-having been estimated. \code{"integrated"} marginalizes the conditional density
-numerically instead: it averages \eqn{f_{A|X}(t \mid X_i)} over the units at
-each of 50 points spanning \code{.exposure}, then interpolates that average back
-to each observed exposure with a cubic spline. It estimates no parameters of
-its own, and it is what WeightIt has done since version 2.0.0.
+the mean of \code{.exposure} and the spread the same estimator gives for it, both
+over the rows the propensity model kept, which are the rows the weights
+describe: a model that dropped rows leaves \code{.propensity} missing at them,
+and those rows carry no weight. Those two moments are parameters of the
+weights, and \code{\link[=ipw]{ipw()}} estimates them alongside the rest of its parameter
+vector, so the standard errors account for the numerator having been
+estimated. \code{"integrated"} marginalizes the conditional density numerically
+instead: it averages \eqn{f_{A|X}(t \mid X_i)} over the units at each of 50
+points spanning \code{.exposure}, then interpolates that average back to each
+observed exposure with a cubic spline. It estimates no parameters of its
+own, and it is what WeightIt has done since version 2.0.0.
 
 The two agree whenever the family is normal and the fitted conditional means
 are themselves normal, since an average of normal densities centered on
@@ -785,14 +787,15 @@ Both densities are read on a standardized residual. The conditional density
 is evaluated at \eqn{z_i = (A_i - \hat{A}_i) / \sigma}, where
 \eqn{\hat{A}_i} is the fitted conditional mean in \code{.propensity} and
 \eqn{\sigma} is the residual spread; the marginal density is evaluated at
-\eqn{z^A_i = (A_i - \bar{A}) / s_A}, where \eqn{\bar{A}} and \eqn{s_A} are
-the mean and standard deviation of \code{.exposure} over the rows the propensity
-model kept. Both moments are read over those rows whether the spread was
-pooled from the model's own residuals or supplied through \code{.sigma}, so the
-two halves of the ratio describe one set of units. Each density is
-then divided by the spread that standardized it, the Jacobian of that change
-of variable, which returns both to the exposure's own units so that each
-integrates to one:
+\eqn{z^A_i = (A_i - \bar{A}) / s_A}, where \eqn{\bar{A}} is the mean of
+\code{.exposure} over the rows the propensity model kept and \eqn{s_A} is the
+spread the same estimator gives for it, the standard deviation unless the
+family estimates a scale of its own. Both moments are read over those rows
+whether the spread was pooled from the model's own residuals or supplied
+through \code{.sigma}, so the two halves of the ratio describe one set of units.
+Each density is then divided by the spread that standardized it, the
+Jacobian of that change of variable, which returns both to the exposure's
+own units so that each integrates to one:
 
 \deqn{w_i = \frac{g(z^A_i) / s_A}{g(z_i) / \sigma}}
 
@@ -846,13 +849,14 @@ a conditional mean is refused there whatever its fitted values happen to fall
 on, a linear probability model included.
 
 \eqn{\sigma} is the pooled residual spread
-\eqn{\sqrt{\mathrm{mean}((A - \hat{A})^2)}} unless \code{.sigma} supplies a single
-standard deviation, or one for each unit, which the model methods do not do
-on their own. \code{\link[=dens_t]{dens_t()}} offers a second estimator, the maximum likelihood
-scale of the t itself, through \code{sigma_method = "mle"}; see \strong{The spread of a
-t density} in its documentation. The marginal density is spread by
-\eqn{s_A} whatever the conditional spread was arrived at, the spread of the
-exposure itself being no business of the conditional model's.
+\eqn{\sqrt{\mathrm{mean}((A - \hat{A})^2)}} unless \code{.sigma} supplies a
+single standard deviation, or one for each unit, which the model methods do
+not do on their own. \code{\link[=dens_t]{dens_t()}} and \code{\link[=dens_laplace]{dens_laplace()}} offer a second
+estimator, the maximum likelihood scale of the family itself, through
+\code{sigma_method = "mle"}; see \strong{The spread of the conditional density} in
+\code{\link[=dens_t]{dens_t()}}'s documentation. A family that estimates its own scale spreads
+both densities of the ratio, so the marginal density is read at that scale
+too and the two halves are densities of one width rather than two.
 
 Every model class is spread that same way, \code{\link[MASS:rlm]{MASS::rlm()}} included. \code{rlm}
 reports a robust scale estimate of its own in \code{fit$s}, which resists the
@@ -867,7 +871,7 @@ propagating the uncertainty of neither.
 \code{\link[=ipw]{ipw()}} models the conditional density with a single conditional variance,
 estimated jointly with the rest of the parameter vector, by the moment the
 pooled spread is the root of or, for a density built with
-\code{sigma_method = "mle"}, by the score of the t. A single
+\code{sigma_method = "mle"}, by the score of the family it names. A single
 \code{.sigma} is taken instead as a known constant: the weights are rebuilt at the
 number that was supplied, and the stacked system carries none of that
 number's uncertainty. An observation-level \code{.sigma} is a different function
diff --git a/tests/testthat/_snaps/ipw-joint-continuous.md b/tests/testthat/_snaps/ipw-joint-continuous.md
index 8eafc3b6..9fb3ed6e 100644
--- a/tests/testthat/_snaps/ipw-joint-continuous.md
+++ b/tests/testthat/_snaps/ipw-joint-continuous.md
@@ -52,7 +52,7 @@
       ! The "ate" weights recomputed from `wt_mod` differ from the weights supplied to `outcome_mod` (compared at relative tolerance 1e-6).
       i The estimand the weights were built for may differ from the one `ipw()` resolved.
       i `ipw()` rebuilt these weights as a "t(df = 4)" density with a "marginal" numerator.
-      i Weights built with an observation-level `.sigma`, such as `influence(model)$sigma`, are one cause: `ipw()` models the conditional density with a single pooled residual root mean square, which is what `wt_ate()` uses when no `.sigma` is given.
+      i Weights built with an observation-level `.sigma`, such as `influence(model)$sigma`, are one cause: `ipw()` models the conditional density with a single spread, the one its family estimates, which is what `wt_ate()` uses when no `.sigma` is given.
       i Weights trimmed, truncated, or normalized after `wt_mod` was fit differ from the ones rebuilt here, which come from that model alone.
       i `.data` values that differ from the data the models were fit to move the recomputed weights on their own and leave the supplied weights exactly right.
       i Refit `outcome_mod` with weights from the two treatment models `wt_mod` holds, and this estimand, if the weights are the cause.
@@ -112,7 +112,7 @@
       ! The "ate" weights recomputed from `wt_mod` differ from the weights supplied to `outcome_mod` (compared at relative tolerance 1e-6).
       i The estimand the weights were built for may differ from the one `ipw()` resolved.
       i `ipw()` rebuilt these weights as a "normal" density with a "marginal" numerator.
-      i Weights built with an observation-level `.sigma`, such as `influence(model)$sigma`, are one cause: `ipw()` models the conditional density with a single pooled residual root mean square, which is what `wt_ate()` uses when no `.sigma` is given.
+      i Weights built with an observation-level `.sigma`, such as `influence(model)$sigma`, are one cause: `ipw()` models the conditional density with a single spread, the one its family estimates, which is what `wt_ate()` uses when no `.sigma` is given.
       i A component built with a fixed `stabilization_score` is one cause: the record on `e` keeps no score, so `ipw()` rebuilt its numerator from the exposure's own marginal distribution instead. A product assembled by hand, built by a version of this package that recorded no score, or subset after it was built, which drops a score held per observation, carries such a record.
       i Weights trimmed, truncated, or normalized after `wt_mod` was fit differ from the ones rebuilt here, which come from that model alone.
       i `.data` values that differ from the data the models were fit to move the recomputed weights on their own and leave the supplied weights exactly right.
diff --git a/tests/testthat/_snaps/weights-continuous.md b/tests/testthat/_snaps/weights-continuous.md
index 862c493f..5460af6e 100644
--- a/tests/testthat/_snaps/weights-continuous.md
+++ b/tests/testthat/_snaps/weights-continuous.md
@@ -157,7 +157,7 @@
       Error in `wt_ate()`:
       ! `numerator` = "integrated" cannot be used with `.sigma`.
       x The marginalization reads the conditional density the propensity model estimated at every unit's fitted mean, and `.sigma` replaces the spread of that density with one of your own.
-      i Leave `.sigma` unset to spread the conditional density by the pooled residual root mean square, or use `numerator` = "marginal", which takes a spread you supply.
+      i Leave `.sigma` unset to spread the conditional density by the estimator its family asks for, or use `numerator` = "marginal", which takes a spread you supply.
 
 ---
 
@@ -187,8 +187,8 @@
     Condition <propensity_density_error>
       Error in `wt_ate()`:
       ! `sigma_method = "mle"` cannot be used with `.sigma`.
-      x `sigma_method = "mle"` estimates the scale of the conditional density from the residuals of the propensity score model, and `.sigma` is a spread of your own that replaces it.
-      i Drop `.sigma` to estimate the scale under the t density, or build the density with `sigma_method = "rms"` to spread the one you supplied.
+      x `sigma_method = "mle"` estimates the scale of the "t" conditional density from the residuals of the propensity score model, and `.sigma` is a spread of your own that replaces it.
+      i Drop `.sigma` to estimate the scale under the "t" density, or build the density with `sigma_method = "rms"` to spread the one you supplied.
 
 # a likelihood with no maximum at a positive scale is refused
 
diff --git a/tests/testthat/test-ipw-continuous.R b/tests/testthat/test-ipw-continuous.R
index 8497e71d..c8669862 100644
--- a/tests/testthat/test-ipw-continuous.R
+++ b/tests/testthat/test-ipw-continuous.R
@@ -2067,8 +2067,8 @@ test_that("ipw() stacks a density the user wrote", {
 # ---- densities and numerators end to end ------------------------------------
 #
 # The combinations the sandwich has to rebuild weights for. Each is a family and
-# a numerator; the spread is pooled throughout, which is the case a fixed
-# `.sigma` and an observation-level one are held against below.
+# a numerator, each spread by the estimator its own family asks for, which is
+# the case a fixed `.sigma` and an observation-level one are held against below.
 continuous_density_cases <- list(
   list(.density = dens_t(4), numerator = "marginal"),
   list(.density = "laplace", numerator = "marginal"),
@@ -2388,7 +2388,10 @@ test_that("ipw() continuous m-estimation standard errors match WeightIt", {
   densities <- list(
     list(ours = "normal", theirs = NULL),
     list(ours = dens_t(df = 4), theirs = "dt_4"),
-    list(ours = "laplace", theirs = "dlaplace")
+    # WeightIt standardizes by the root mean square whatever the family, so the
+    # parity is written at that estimator rather than at the scale of the
+    # Laplace, which is what the family asks for on its own.
+    list(ours = dens_laplace(sigma_method = "rms"), theirs = "dlaplace")
   )
 
   for (dens in densities) {
@@ -2452,8 +2455,22 @@ test_that("ipw() solves the scale equation of a t spread by maximum likelihood",
   # the number it estimates supplied as a constant. The weights are equal, so
   # the point estimates are, and what differs between the two fits is what the
   # sandwich accounts for.
-  mle <- fit_continuous_models(dat, .density = dens_t(6, sigma_method = "mle"))
-  fixed <- fit_continuous_models(dat, .density = dens_t(6), .sigma = scale)
+  #
+  # Unstabilized, so that the conditional density is the whole of the ratio. A
+  # marginal numerator is read at the estimator the family it names asks for,
+  # which a spread supplied for the conditional density does not replace, so the
+  # two arms would divide by the same number and multiply by different ones.
+  mle <- fit_continuous_models(
+    dat,
+    .density = dens_t(6, sigma_method = "mle"),
+    stabilize = FALSE
+  )
+  fixed <- fit_continuous_models(
+    dat,
+    .density = dens_t(6),
+    .sigma = scale,
+    stabilize = FALSE
+  )
   expect_equal(as.double(mle$wts), as.double(fixed$wts), tolerance = 1e-8)
 
   res_mle <- ipw(mle$ps_mod, mle$outcome_mod)
diff --git a/tests/testthat/test-ipw-psi.R b/tests/testthat/test-ipw-psi.R
index c96d0a3b..9db9c4db 100644
--- a/tests/testthat/test-ipw-psi.R
+++ b/tests/testthat/test-ipw-psi.R
@@ -432,11 +432,10 @@ continuous_spec <- function(
     ),
     density = meta$density,
     numerator = meta$numerator,
-    sigma = if (is.null(.sigma)) {
-      list(kind = "pooled", value = NULL)
-    } else {
-      list(kind = "fixed", value = .sigma)
-    },
+    # Read off the record the weights carry, the way `ipw_continuous_spread()`
+    # reads it, rather than off `.sigma` alone: a family that estimates its
+    # scale under itself records a spread of its own with no `.sigma` in sight.
+    sigma = ipw_continuous_spread(meta),
     outcome = list(
       X = model.matrix(msm),
       y = dat[[outcome_var]],
@@ -1622,7 +1621,6 @@ test_that("a degenerate maximum likelihood scale is refused in the caller's name
   # than against the seed.
   cnd <- rlang::catch_cnd(
     ipw_continuous_sigma2_seed(
-      list(kind = "mle"),
       c(rep(0, 9), 1),
       dens_t(4, sigma_method = "mle"),
       call = rlang::call2("ipw")
diff --git a/tests/testthat/test-weights-continuous.R b/tests/testthat/test-weights-continuous.R
index add56aa7..300a462b 100644
--- a/tests/testthat/test-weights-continuous.R
+++ b/tests/testthat/test-weights-continuous.R
@@ -192,12 +192,16 @@ test_that("Student's t weights are the t density ratio computed by hand", {
 })
 
 test_that("Laplace weights are the Laplace density ratio computed by hand", {
-  z <- continuous_density_z()
-  sigma <- continuous_density_pooled()
-  sd_a <- continuous_density_sd_a()
+  exposure <- continuous_density_data$exposure
+  mu <- continuous_density_data$mu
+
+  # The Laplace is read at its own scale, the mean absolute residual, on both
+  # halves of the ratio.
+  scale <- mean(abs(exposure - mu))
+  scale_a <- mean(abs(exposure - mean(exposure)))
 
-  f_den <- exp(-abs(z)) / 2 / sigma
-  f_num <- exp(-abs(continuous_density_z_a())) / 2 / sd_a
+  f_den <- exp(-abs((exposure - mu) / scale)) / 2 / scale
+  f_num <- exp(-abs((exposure - mean(exposure)) / scale_a)) / 2 / scale_a
 
   by_string <- continuous_density_wt(.density = "laplace")
   by_spec <- continuous_density_wt(.density = dens_laplace())
@@ -1468,7 +1472,10 @@ test_that("the denominator is the one WeightIt divides by", {
       g = stats::dt(z, df = 4)
     ),
     list(
-      ours = "laplace",
+      # WeightIt standardizes by the root mean square whatever the family, so
+      # the parity is written at that estimator rather than at the scale of the
+      # Laplace, which is what the family asks for on its own.
+      ours = dens_laplace(sigma_method = "rms"),
       theirs = "dlaplace",
       g = exp(-abs(z)) / 2
     )
@@ -1548,9 +1555,9 @@ continuous_integrated_numerator <- function(
 continuous_integrated_wt <- function(
   g,
   exposure = continuous_density_data$exposure,
-  mu = continuous_density_data$mu
+  mu = continuous_density_data$mu,
+  sigma = sqrt(mean((exposure - mu)^2))
 ) {
-  sigma <- sqrt(mean((exposure - mu)^2))
   f_den <- g((exposure - mu) / sigma) / sigma
   f_num <- continuous_integrated_numerator(
     g,
@@ -1562,6 +1569,21 @@ continuous_integrated_wt <- function(
   f_num / f_den
 }
 
+# The same ratio for the Laplace, which is read at the scale of the Laplace
+# rather than at the root mean square. An integrated numerator reads the
+# conditional density in both places, so the one scale spreads the whole ratio.
+continuous_integrated_laplace_wt <- function(
+  exposure = continuous_density_data$exposure,
+  mu = continuous_density_data$mu
+) {
+  continuous_integrated_wt(
+    function(z) exp(-abs(z)) / 2,
+    exposure = exposure,
+    mu = mu,
+    sigma = mean(abs(exposure - mu))
+  )
+}
+
 # The integrated ratio for a kernel, which is one estimate rather than two: the
 # same fit answers for the standardized residuals and for the whole standardized
 # grid, so it is fit over a range covering both.
@@ -1603,13 +1625,29 @@ continuous_integrated_kernel_wt <- function(
 }
 
 test_that("the integrated numerator is the grid marginalization by hand", {
+  # An integrated numerator reads the conditional density in both places, so the
+  # whole ratio is spread by the estimator the family it names asks for.
+  pooled <- continuous_density_pooled()
+  laplace_scale <- mean(abs(
+    continuous_density_data$exposure - continuous_density_data$mu
+  ))
+
   families <- list(
-    list(input = "normal", g = function(z) stats::dnorm(z)),
-    list(input = dens_t(df = 4), g = function(z) stats::dt(z, df = 4)),
-    list(input = "laplace", g = function(z) exp(-abs(z)) / 2),
+    list(input = "normal", g = function(z) stats::dnorm(z), sigma = pooled),
+    list(
+      input = dens_t(df = 4),
+      g = function(z) stats::dt(z, df = 4),
+      sigma = pooled
+    ),
+    list(
+      input = "laplace",
+      g = function(z) exp(-abs(z)) / 2,
+      sigma = laplace_scale
+    ),
     list(
       input = function(z) stats::dlogis(z),
-      g = function(z) stats::dlogis(z)
+      g = function(z) stats::dlogis(z),
+      sigma = pooled
     )
   )
 
@@ -1621,7 +1659,7 @@ test_that("the integrated numerator is the grid marginalization by hand", {
 
     expect_equal(
       as.numeric(weights),
-      continuous_integrated_wt(family$g),
+      continuous_integrated_wt(family$g, sigma = family$sigma),
       tolerance = 1e-12
     )
   }
@@ -2297,7 +2335,7 @@ test_that("the numerator reaches the weights through a data frame", {
 
   expect_equal(
     as.numeric(from_frame),
-    continuous_integrated_wt(function(z) exp(-abs(z)) / 2),
+    continuous_integrated_laplace_wt(),
     tolerance = 1e-12
   )
   expect_identical(density_meta(from_frame)$numerator, "integrated")
@@ -2398,7 +2436,7 @@ test_that("censoring weights carry a numerator through their other methods", {
 
   expect_equal(
     as.numeric(from_frame),
-    continuous_integrated_wt(function(z) exp(-abs(z)) / 2),
+    continuous_integrated_laplace_wt(),
     tolerance = 1e-12
   )
   expect_identical(estimand(from_frame), "uncensored")
@@ -2489,7 +2527,10 @@ test_that("integrated weights are the weights WeightIt gives", {
   families <- list(
     list(ours = "normal", theirs = NULL),
     list(ours = dens_t(df = 4), theirs = "dt_4"),
-    list(ours = "laplace", theirs = "dlaplace"),
+    # WeightIt standardizes by the root mean square whatever the family, so the
+    # parity is written at that estimator rather than at the scale of the
+    # Laplace, which is what the family asks for on its own.
+    list(ours = dens_laplace(sigma_method = "rms"), theirs = "dlaplace"),
     list(ours = "kernel", theirs = "kernel")
   )
 
@@ -2640,18 +2681,17 @@ test_that("weights record a maximum likelihood spread as a source of its own", {
   expect_null(density_meta(mle)$sigma_value)
 })
 
-test_that("a maximum likelihood scale spreads the conditional density alone", {
+test_that("a maximum likelihood scale spreads the density it is read under", {
   # `.sigma` and `sigma_method` describe the same thing: the spread of the
-  # conditional density the weights divide by. The marginal density that
-  # stabilizes them is the exposure's own, read at the exposure's own mean and
-  # root mean square, which neither of them has ever changed.
+  # conditional density the weights divide by, which is the whole of an
+  # unstabilized ratio.
   scale <- t_scale_mle(continuous_t_residuals, df = 6)
 
-  stabilized <- wt_ate(
+  mle <- wt_ate(
     continuous_t_data$mu,
     continuous_t_data$exposure,
     exposure_type = "continuous",
-    stabilize = TRUE,
+    stabilize = FALSE,
     .density = dens_t(6, sigma_method = "mle")
   )
   oracle <- wt_ate(
@@ -2659,11 +2699,33 @@ test_that("a maximum likelihood scale spreads the conditional density alone", {
     continuous_t_data$exposure,
     .sigma = scale,
     exposure_type = "continuous",
-    stabilize = TRUE,
+    stabilize = FALSE,
     .density = dens_t(6)
   )
 
-  expect_equal(as.numeric(stabilized), as.numeric(oracle), tolerance = 1e-8)
+  expect_equal(as.numeric(mle), as.numeric(oracle), tolerance = 1e-8)
+
+  # The marginal density that stabilizes the weights is the exposure's own, read
+  # at its own mean and at the spread the same estimator gives for it, so the
+  # two halves of the ratio are densities of one width rather than two.
+  exposure <- continuous_t_data$exposure
+  scale_a <- t_scale_mle(exposure - mean(exposure), df = 6)
+
+  stabilized <- wt_ate(
+    continuous_t_data$mu,
+    exposure,
+    exposure_type = "continuous",
+    stabilize = TRUE,
+    .density = dens_t(6, sigma_method = "mle")
+  )
+
+  expect_equal(
+    as.numeric(stabilized),
+    stats::dt((exposure - mean(exposure)) / scale_a, df = 6) /
+      scale_a *
+      as.numeric(mle),
+    tolerance = 1e-8
+  )
 })
 
 test_that("a maximum likelihood scale is fit on the residuals that are there", {

From dc96bc8cb976afe11df2293ad497860b58a8362e Mon Sep 17 00:00:00 2001
From: Malcolm Barrett <malcolmbarrett@gmail.com>
Date: Wed, 16 Sep 2026 12:20:39 -0400
Subject: [PATCH 3/7] add failing tests for the t scale default and the
 numerator model spread

---
 tests/testthat/test-density.R            |  12 +-
 tests/testthat/test-ipw-continuous.R     |  55 ++++++-
 tests/testthat/test-weights-continuous.R | 190 ++++++++++++++++++++++-
 3 files changed, 243 insertions(+), 14 deletions(-)

diff --git a/tests/testthat/test-density.R b/tests/testthat/test-density.R
index ec00f01e..1f53295c 100644
--- a/tests/testthat/test-density.R
+++ b/tests/testthat/test-density.R
@@ -83,12 +83,12 @@ test_that("dens_t() records its degrees of freedom", {
 })
 
 test_that("dens_t() records how its scale is estimated", {
-  # The scale of the standardized residuals is estimated by the root mean square
-  # that spreads every other family, or by maximum likelihood under the t
-  # itself. It is not a parameter of the density: `t(df = 4)` is the same
-  # density however its scale was arrived at, so the choice sits beside the
-  # parameters rather than among them, and the printed density is what it was.
-  expect_identical(dens_t(df = 4)$sigma_method, "rms")
+  # The scale of the standardized residuals is estimated by maximum likelihood
+  # under the t itself, which is the default, or by the root mean square. It is
+  # not a parameter of the density: `t(df = 4)` is the same density however its
+  # scale was arrived at, so the choice sits beside the parameters rather than
+  # among them, and the printed density is what it was.
+  expect_identical(dens_t(df = 4)$sigma_method, "mle")
   expect_identical(dens_t(df = 4, sigma_method = "rms")$sigma_method, "rms")
   expect_identical(dens_t(df = 4, sigma_method = "mle")$sigma_method, "mle")
 
diff --git a/tests/testthat/test-ipw-continuous.R b/tests/testthat/test-ipw-continuous.R
index c8669862..0519a447 100644
--- a/tests/testthat/test-ipw-continuous.R
+++ b/tests/testthat/test-ipw-continuous.R
@@ -2195,9 +2195,18 @@ test_that("ipw_spec_continuous reads the density, numerator, and spread off the
   spec <- ipw_spec_continuous(mods$ps_mod, mods$outcome_mod)
   expect_true(density_specs_agree(spec$density, dens_t(3)))
   expect_equal(spec$numerator, "integrated")
-  expect_equal(spec$sigma$kind, "pooled")
+
+  # A t is spread under itself and a normal by the root mean square, which is
+  # the normal's own scale parameter, so both estimators a fit can be made under
+  # are read back off the weights.
+  expect_equal(spec$sigma$kind, "mle")
   expect_null(spec$sigma$value)
 
+  normal <- fit_continuous_models(dat, .density = dens_normal())
+  normal_spec <- ipw_spec_continuous(normal$ps_mod, normal$outcome_mod)
+  expect_equal(normal_spec$sigma$kind, "pooled")
+  expect_null(normal_spec$sigma$value)
+
   fixed <- fit_continuous_models(dat, .sigma = 1.25)
   spec_fixed <- ipw_spec_continuous(fixed$ps_mod, fixed$outcome_mod)
   expect_equal(spec_fixed$sigma$kind, "fixed")
@@ -2783,6 +2792,50 @@ test_that("a stacked continuous numerator block solves over the rows .data keeps
   )))
 })
 
+test_that("a numerator model's spread is stacked at the scale its family asks for", {
+  dat <- sim_continuous()
+  num_mod <- lm(A ~ x1, data = dat)
+  fits <- fit_continuous_models(
+    dat,
+    stabilize = num_mod,
+    .density = dens_laplace()
+  )
+
+  # Both halves of the ratio are Laplace densities, so the numerator's own block
+  # carries the scale of a Laplace fit to its residuals rather than the second
+  # moment a normal's block carries. The weights and the block's seed are two
+  # computations of that one quantity, and `ipw_weights_at_init()` refuses at
+  # 1e-6 when they part, so a route that corrected the weights alone would not
+  # reach this line.
+  res <- muffle_coverage_warning(ipw(
+    fits$ps_mod,
+    fits$outcome_mod,
+    .data = dat
+  ))
+
+  expect_s3_class(res, "ipw")
+
+  theta <- coef(res$fit)
+  expect_true("sigma2_n" %in% names(theta))
+
+  # The row the block solves is the Laplace score, whose root is the mean
+  # absolute residual, squared because the parameter carries a squared scale.
+  expect_equal(
+    unname(theta[["sigma2_n"]]),
+    mean(abs(residuals(num_mod)))^2,
+    tolerance = 1e-8
+  )
+
+  # The second moment is what that parameter would hold if the numerator were
+  # still read at the estimator the normal asks for, and these residuals tell
+  # the two apart.
+  expect_false(isTRUE(all.equal(
+    unname(theta[["sigma2_n"]]),
+    mean(residuals(num_mod)^2),
+    tolerance = 1e-3
+  )))
+})
+
 # ---- a dose numerator covariate `.data` supplies as another type ------------
 #
 # A numerator model is rebuilt from `.data` the way the other models are, so a
diff --git a/tests/testthat/test-weights-continuous.R b/tests/testthat/test-weights-continuous.R
index 300a462b..d38a9e43 100644
--- a/tests/testthat/test-weights-continuous.R
+++ b/tests/testthat/test-weights-continuous.R
@@ -686,7 +686,7 @@ test_that("the weights record the density family they were built from", {
       recorded = "laplace",
       sigma = "mle"
     ),
-    t = list(input = dens_t(df = 4), recorded = "t(df = 4)", sigma = "pooled"),
+    t = list(input = dens_t(df = 4), recorded = "t(df = 4)", sigma = "mle"),
     kernel = list(
       input = dens_kernel(adjust = 1.5),
       recorded = 'kernel(bw = "nrd0", adjust = 1.5, kernel = "gaussian", n = 512)',
@@ -2229,7 +2229,7 @@ test_that("an interpolated numerator that dips below zero is refused", {
 test_that("integrated weights record the numerator they were built from", {
   families <- list(
     list(input = "normal", sigma = "pooled"),
-    list(input = dens_t(df = 4), sigma = "pooled"),
+    list(input = dens_t(df = 4), sigma = "mle"),
     list(input = "laplace", sigma = "mle"),
     list(input = "kernel", sigma = "pooled")
   )
@@ -2598,6 +2598,12 @@ continuous_t_data <- local({
 
 continuous_t_residuals <- continuous_t_data$exposure - continuous_t_data$mu
 
+# The t density on the exposure's own scale: the standard density read at the
+# residual standardized by `scale`, divided by the scale that standardized it.
+t_density <- function(residuals, scale, df) {
+  stats::dt(residuals / scale, df = df) / scale
+}
+
 # Unstabilized weights for that problem, which are the conditional density and
 # nothing else. The marginal density that stabilizes weights is read at the
 # exposure's own moments, and what is under test here is the spread of the
@@ -2656,15 +2662,52 @@ test_that("the two spreads meet as the t approaches the normal", {
   expect_lt(max(abs(as.numeric(mle) / as.numeric(rms) - 1)), 0.01)
 })
 
-test_that("the t family is spread by the root mean square unless asked otherwise", {
-  # The default is the estimator every family takes, so weights written the way
-  # they have always been written are the weights they have always been.
+test_that("a t denominator is read at the scale of the t", {
+  residuals <- continuous_t_residuals
+  scale <- t_scale_mle(residuals, df = 4)
+
+  # The whole of an unstabilized ratio is the conditional density, so the
+  # weights say on their own what spread that density was read at, held against
+  # an oracle that maximizes the likelihood by search rather than by solving the
+  # equation the package solves.
+  weights <- continuous_t_wt(.density = dens_t(df = 4))
+
+  expect_equal(
+    as.numeric(weights),
+    1 / t_density(residuals, scale, df = 4),
+    tolerance = 1e-8
+  )
+
+  # The spread the weights were built at is the scale parameter of the t rather
+  # than the root mean square, and these residuals tell the two apart: a t with
+  # four degrees of freedom has a variance of twice its squared scale.
+  rms <- sqrt(mean(residuals^2))
+
+  expect_false(isTRUE(all.equal(
+    as.numeric(weights),
+    1 / t_density(residuals, rms, df = 4),
+    tolerance = 1e-2
+  )))
+})
+
+test_that("the t family is spread under itself unless asked otherwise", {
+  # The density is read at a residual standardized by the spread and divided by
+  # that spread, so the spread is the scale parameter of the family reading it.
+  # The root mean square is the scale parameter of the normal alone, and a
+  # caller who wants a t spread by it says so.
   default <- continuous_t_wt(.density = dens_t(6))
   asked <- continuous_t_wt(.density = dens_t(6, sigma_method = "rms"))
+  under_mle <- continuous_t_wt(.density = dens_t(6, sigma_method = "mle"))
 
-  expect_identical(as.numeric(default), as.numeric(asked))
-  expect_identical(density_meta(default)$sigma, "pooled")
+  expect_identical(density_meta(default)$sigma, "mle")
   expect_identical(density_meta(asked)$sigma, "pooled")
+
+  expect_identical(as.numeric(default), as.numeric(under_mle))
+  expect_false(isTRUE(all.equal(
+    as.numeric(default),
+    as.numeric(asked),
+    tolerance = 1e-3
+  )))
 })
 
 test_that("weights record a maximum likelihood spread as a source of its own", {
@@ -3472,6 +3515,139 @@ test_that("censoring weights reach the numerator model route too", {
   expect_identical(density_meta(weights)$numerator_model, fit)
 })
 
+# ---- the spread a numerator model's density is read at ----------------------
+
+# Conditional numerators for the two problems whose residuals are not normal:
+# the exposure fit on a baseline covariate the propensity score model does not
+# read, so that the numerator model leaves residuals of its own and a spread of
+# its own for its half of the ratio to be read at.
+continuous_laplace_numerator_model <- local({
+  v <- withr::with_seed(
+    20260906,
+    stats::rnorm(length(continuous_laplace_data$exposure))
+  )
+  exposure <- continuous_laplace_data$exposure
+
+  stats::lm(exposure ~ v)
+})
+
+continuous_t_numerator_model <- local({
+  v <- withr::with_seed(
+    20260907,
+    stats::rnorm(length(continuous_t_data$exposure))
+  )
+  exposure <- continuous_t_data$exposure
+
+  stats::lm(exposure ~ v)
+})
+
+test_that("a Laplace numerator model is read at the scale of the Laplace", {
+  exposure <- continuous_laplace_data$exposure
+  fit <- continuous_laplace_numerator_model
+  numerator_residuals <- as.numeric(stats::residuals(fit))
+
+  denominator_scale <- mean(abs(continuous_laplace_residuals))
+  numerator_scale <- mean(abs(numerator_residuals))
+
+  weights <- wt_ate(
+    continuous_laplace_data$mu,
+    exposure,
+    exposure_type = "continuous",
+    stabilize = fit,
+    .density = dens_laplace()
+  )
+
+  # Both halves of the ratio are Laplace densities, so both are spread by the
+  # scale of a Laplace: the numerator by the scale of its own model's residuals,
+  # the denominator by the scale of the propensity score model's.
+  expect_equal(
+    as.numeric(weights),
+    laplace_density(numerator_residuals, numerator_scale) /
+      laplace_density(continuous_laplace_residuals, denominator_scale),
+    tolerance = 1e-12
+  )
+
+  # A numerator spread by the root mean square while the denominator is spread
+  # by the scale is a ratio of two densities of different widths, and is the
+  # ratio these weights are not.
+  expect_false(isTRUE(all.equal(
+    as.numeric(weights),
+    laplace_density(
+      numerator_residuals,
+      sqrt(mean(numerator_residuals^2))
+    ) /
+      laplace_density(continuous_laplace_residuals, denominator_scale),
+    tolerance = 1e-3
+  )))
+})
+
+test_that("a numerator model asked for the root mean square is read at it", {
+  # The estimator follows what the density asks for rather than the family it
+  # names, on a fitted numerator as well as on a marginal one.
+  exposure <- continuous_laplace_data$exposure
+  fit <- continuous_laplace_numerator_model
+  numerator_residuals <- as.numeric(stats::residuals(fit))
+
+  weights <- wt_ate(
+    continuous_laplace_data$mu,
+    exposure,
+    exposure_type = "continuous",
+    stabilize = fit,
+    .density = dens_laplace(sigma_method = "rms")
+  )
+
+  expect_equal(
+    as.numeric(weights),
+    laplace_density(
+      numerator_residuals,
+      sqrt(mean(numerator_residuals^2))
+    ) /
+      laplace_density(
+        continuous_laplace_residuals,
+        sqrt(mean(continuous_laplace_residuals^2))
+      ),
+    tolerance = 1e-12
+  )
+})
+
+test_that("a t numerator model is read at the scale of the t", {
+  exposure <- continuous_t_data$exposure
+  fit <- continuous_t_numerator_model
+  numerator_residuals <- as.numeric(stats::residuals(fit))
+
+  denominator_scale <- t_scale_mle(continuous_t_residuals, df = 4)
+  numerator_scale <- t_scale_mle(numerator_residuals, df = 4)
+
+  weights <- wt_ate(
+    continuous_t_data$mu,
+    exposure,
+    exposure_type = "continuous",
+    stabilize = fit,
+    .density = dens_t(df = 4)
+  )
+
+  expect_equal(
+    as.numeric(weights),
+    t_density(numerator_residuals, numerator_scale, df = 4) /
+      t_density(continuous_t_residuals, denominator_scale, df = 4),
+    tolerance = 1e-8
+  )
+
+  # The root mean square of residuals this heavy-tailed is far from the scale
+  # the likelihood reaches, so a numerator read at it is a numerator of a
+  # different width from the denominator it is divided by.
+  expect_false(isTRUE(all.equal(
+    as.numeric(weights),
+    t_density(
+      numerator_residuals,
+      sqrt(mean(numerator_residuals^2)),
+      df = 4
+    ) /
+      t_density(continuous_t_residuals, denominator_scale, df = 4),
+    tolerance = 1e-3
+  )))
+})
+
 # ---- the rows a gappy fit's numerator is read over --------------------------
 #
 # A propensity score model fit with `na.action = na.exclude` reports its fitted

From 7ddd698dca1efd400f4717dfab13fff2819c7187 Mon Sep 17 00:00:00 2001
From: Malcolm Barrett <malcolmbarrett@gmail.com>
Date: Wed, 16 Sep 2026 12:53:40 -0400
Subject: [PATCH 4/7] estimate the t scale under the t by default and read a
 numerator model at its family's spread

---
 NEWS.md                                     | 34 +++++++--
 R/density.R                                 | 38 ++++++----
 R/ipw-continuous-models.R                   | 17 +++--
 R/ipw-joint-models.R                        | 14 +++-
 R/ipw-mestimation.R                         |  1 +
 R/ipw-psi.R                                 | 30 ++++++--
 R/weights.R                                 |  9 +--
 man/dens_normal.Rd                          | 16 +++--
 man/wt_ate.Rd                               |  7 +-
 tests/testthat/_snaps/psw-records.md        |  4 +-
 tests/testthat/_snaps/weights-continuous.md |  1 +
 tests/testthat/test-ipw-continuous.R        |  8 +--
 tests/testthat/test-weights-continuous.R    | 79 +++++++++++++--------
 13 files changed, 174 insertions(+), 84 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index 93291099..8ba78d86 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,30 @@
 # propensity 0.1.0.9000 (development version)
 
+* `dens_t()` now estimates the scale of the t under the t itself by default,
+  as `dens_laplace()` estimates the scale of the Laplace. Both densities of a
+  continuous exposure's weights are read at a residual standardized by a spread
+  and divided by that spread, so the spread is the scale parameter of the family
+  reading it, and the root mean square is the scale parameter of the normal
+  alone. Weights built from `dens_t()` without a `sigma_method` change: they
+  record `sigma = "mle"` in `density_meta()`, and `ipw()` stacks the score the
+  t itself gives for the scale in place of the moment equation the root mean
+  square is the root of. Write `dens_t(df, sigma_method = "rms")` for the
+  spread the family was read at before, which is also the spread software that
+  standardizes every family alike reads it at. A `.sigma` supplied alongside a
+  `dens_t()` written without a `sigma_method` is now refused, being a second
+  instruction about the same quantity, and the refusal names
+  `sigma_method = "rms"` as the way to supply a spread of your own.
+
+* A numerator model supplied to `stabilize` for a continuous exposure is now
+  read at the spread its family asks for rather than always at the root mean
+  square of its residuals. Under `dens_laplace()`, or a `dens_t()` fit under
+  itself, the numerator was the one density of the ratio still spread by the
+  root mean square, so the two halves of the ratio were densities of different
+  widths. The weights, the value `ipw()` seeds that spread at, and the
+  estimating equation its stacked system solves for it now all read the
+  estimator the family asks for, applied to the numerator model's own
+  residuals, on the single-treatment route and the joint one alike.
+
 * A joint weight's record now keeps one copy of what stabilized each
   component. A component weighting a dose kept its numerator model in the
   `density` record `joint_wt_meta()` reports and in that record's
@@ -297,10 +322,11 @@
   estimator gives for it. Weights built that way record `sigma = "mle"` in
   `density_meta()`, and `ipw()` stacks that score for each of the two scales in
   place of the moment equation the pooled spread is the root of, so the
-  sandwich accounts for having estimated them. The default is unchanged, and a
-  `.sigma` supplied alongside `sigma_method = "mle"` is refused with an error
-  of class `propensity_density_error`, being a second instruction about the
-  same quantity. Residuals a model reproduced exactly say nothing about the
+  sandwich accounts for having estimated them. It is what `dens_t()` takes by
+  default, as the entry on that default describes, and a `.sigma` supplied
+  alongside `sigma_method = "mle"` is refused with an error of class
+  `propensity_density_error`, being a second instruction about the same
+  quantity. Residuals a model reproduced exactly say nothing about the
   spread of the density around it, so a fit with enough of them that the
   likelihood has no maximum at a positive scale is refused with that same
   class, naming how many they are.
diff --git a/R/density.R b/R/density.R
index 444ad0a9..7f3a9cf8 100644
--- a/R/density.R
+++ b/R/density.R
@@ -55,9 +55,9 @@ new_density_spec <- function(
 #'   absolute residual, or by the root mean square of the residuals, which is
 #'   the spread every family without an estimator of its own is read at.
 #' * `dens_t()` is Student's t density with `df` degrees of freedom, heavier
-#'   tailed still, and heavier the smaller `df` is. Its scale is estimated by
-#'   the root mean square of the residuals, or by maximum likelihood under the
-#'   t itself.
+#'   tailed still, and heavier the smaller `df` is. Its scale is estimated
+#'   under the t itself, by maximum likelihood, or by the root mean square of
+#'   the residuals.
 #' * `dens_kernel()` is a kernel density estimate of the standardized
 #'   residuals, fit with [stats::density()] and interpolated to each
 #'   observation. It assumes no family at all, at the cost of a density that
@@ -104,9 +104,11 @@ new_density_spec <- function(
 #' family. The scale of such a family is smaller than its standard deviation,
 #' and the root mean square is pulled outward by the large residuals a heavy
 #' tail produces, which is what the family was chosen to accommodate.
-#' `sigma_method = "mle"` estimates the scale under the family itself.
-#' `dens_laplace()` takes it by default, where it is the mean absolute
-#' residual. For `dens_t()` it is the root of
+#' `sigma_method = "mle"` estimates the scale under the family itself, and both
+#' `dens_laplace()` and `dens_t()` take it by default, since the spread the
+#' density is divided by is the scale parameter of the family reading it. For
+#' `dens_laplace()` it is the mean absolute residual. For `dens_t()` it is the
+#' root of
 #' \eqn{\sum_i \left[(\nu + 1) r_i^2 / (\nu \sigma^2 + r_i^2)\right] = n},
 #' where \eqn{r_i} is the residual and \eqn{\nu} is `df`. Each residual enters
 #' that sum through a bounded term, so a residual far out in the tail moves the
@@ -166,7 +168,7 @@ dens_laplace <- function(sigma_method = c("mle", "rms")) {
 
 #' @rdname dens_normal
 #' @export
-dens_t <- function(df, sigma_method = c("rms", "mle")) {
+dens_t <- function(df, sigma_method = c("mle", "rms")) {
   check_density_number(df, "df")
   sigma_method <- rlang::arg_match(sigma_method)
   force(df)
@@ -621,15 +623,24 @@ check_numerator_model_not_multinom <- function(
 }
 
 # What the numerator model contributes to the ratio: the conditional mean it
-# fits for each unit and the spread of its residuals around them. The spread is
-# the uncentered root mean square of that model's own residuals, taken about
+# fits for each unit and the spread of its residuals around them, taken about
 # zero rather than about their mean because the fitted mean is what the ratio
-# reads the density at.
-numerator_model_moments <- function(numerator_model) {
+# reads the density at. The spread is read by the estimator the family asks
+# for, which is the estimator the conditional density was spread by, so both
+# densities of the ratio are spread by one estimator applied to each model's
+# own residuals and are densities of the same width.
+numerator_model_moments <- function(
+  numerator_model,
+  density,
+  call = rlang::caller_env()
+) {
   mu <- as.numeric(stats::fitted(numerator_model))
   residuals <- as.numeric(stats::residuals(numerator_model, type = "response"))
 
-  list(mu = mu, sigma = sqrt(mean(residuals^2, na.rm = TRUE)))
+  list(
+    mu = mu,
+    sigma = density_scale_estimate(residuals, density, call = call)
+  )
 }
 
 # What a binary exposure's numerator model contributes to the ratio: the fitted
@@ -891,6 +902,9 @@ check_sigma_method <- function(.sigma, density, call = rlang::caller_env()) {
            {.val {family}} conditional density from the residuals of the
            propensity score model, 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
+           carries whether or not you wrote it.",
       i = "Drop {.arg .sigma} to estimate the scale under the {.val {family}}
            density, or build the density with {.code sigma_method = \"rms\"}
            to spread the one you supplied."
diff --git a/R/ipw-continuous-models.R b/R/ipw-continuous-models.R
index 2525039a..aeea2ee2 100644
--- a/R/ipw-continuous-models.R
+++ b/R/ipw-continuous-models.R
@@ -974,6 +974,7 @@ ipw_continuous_ratio_meta <- function(
 # stabilized components unable to tell which fit is refused.
 ipw_numerator_model_block <- function(
   numerator_model,
+  density,
   .data = NULL,
   component = NULL,
   call = rlang::caller_env()
@@ -1026,13 +1027,15 @@ ipw_numerator_model_block <- function(
     penalty = entry$penalty,
     coefs = stats::coef(numerator_model),
     # The spread the numerator's density was read at, taken the way
-    # `numerator_model_moments()` takes it: the mean square of the fit's own
-    # response-scale residuals, over the fit's own rows. The design above is
-    # rebuilt over the rows `.data` leaves, and the same moment read over those
-    # rows is a different spread from the one the weights carry.
-    sigma2_fit = mean(
-      as.numeric(stats::residuals(numerator_model, type = "response"))^2,
-      na.rm = TRUE
+    # `numerator_model_moments()` takes it: the estimator the family asks for,
+    # read on the fit's own response-scale residuals, over the fit's own rows.
+    # The design above is rebuilt over the rows `.data` leaves, and the same
+    # estimator read over those rows is a different spread from the one the
+    # weights carry.
+    sigma2_fit = ipw_continuous_sigma2_seed(
+      as.numeric(stats::residuals(numerator_model, type = "response")),
+      density,
+      call = call
     )
   )
 }
diff --git a/R/ipw-joint-models.R b/R/ipw-joint-models.R
index e40cbc00..f3217aeb 100644
--- a/R/ipw-joint-models.R
+++ b/R/ipw-joint-models.R
@@ -1917,6 +1917,7 @@ ipw_joint_models_stab_components <- function(
         scores[[i]],
         names[[i]],
         n,
+        joint_wt_dose_density(wts, dose_idx),
         .data = .data,
         call = call
       ))
@@ -2034,6 +2035,7 @@ ipw_joint_models_dose_stab <- function(
   score,
   name,
   n,
+  density,
   .data = NULL,
   call = rlang::caller_env()
 ) {
@@ -2077,6 +2079,7 @@ ipw_joint_models_dose_stab <- function(
 
   model <- ipw_numerator_model_block(
     numerator_model,
+    density,
     .data = .data,
     component = name,
     call = call
@@ -2465,7 +2468,7 @@ ipw_init_joint_models_ps <- function(spec, call = rlang::caller_env()) {
 # that fit's rows, since those are the moments the weights were built at. The
 # rows the spec analyzes are what the equations solve the moments over, which is
 # the answer rather than the starting value.
-ipw_init_joint_models_stab <- function(spec) {
+ipw_init_joint_models_stab <- function(spec, call = rlang::caller_env()) {
   blocks <- lapply(seq_along(spec$stab$components), function(i) {
     stab <- spec$stab$components[[i]]
 
@@ -2527,7 +2530,12 @@ ipw_init_joint_models_stab <- function(spec) {
     c(
       coefs,
       stats::setNames(
-        ipw_numerator_fit_sigma2(model, spec$exposure[[i]]),
+        ipw_numerator_fit_sigma2(
+          model,
+          spec$exposure[[i]],
+          spec$density,
+          call = call
+        ),
         paste0(prefix, "sigma2")
       )
     )
@@ -2541,7 +2549,7 @@ ipw_init_joint_models_stab <- function(spec) {
 ipw_init_joint_models <- function(spec, call = rlang::caller_env()) {
   beta <- spec$outcome$coefs
   ps_block <- ipw_init_joint_models_ps(spec, call = call)
-  stab_block <- ipw_init_joint_models_stab(spec)
+  stab_block <- ipw_init_joint_models_stab(spec, call = call)
 
   # A dose has no cells, so the surface reports the outcome block directly and
   # the system carries no marginal means and no contrasts over them.
diff --git a/R/ipw-mestimation.R b/R/ipw-mestimation.R
index 189123b9..51db2745 100644
--- a/R/ipw-mestimation.R
+++ b/R/ipw-mestimation.R
@@ -2396,6 +2396,7 @@ ipw_spec_continuous <- function(
   numerator_block <- if (stacked) {
     ipw_numerator_model_block(
       numerator_mod,
+      ratio$density,
       .data = if (!is.null(.data)) mm_data,
       call = call
     )
diff --git a/R/ipw-psi.R b/R/ipw-psi.R
index 93e42b79..85a65180 100644
--- a/R/ipw-psi.R
+++ b/R/ipw-psi.R
@@ -714,17 +714,23 @@ ipw_continuous_fit_design <- function(spec) {
 }
 
 # The same reading for a numerator model, whose block carries the spread its own
-# fit's residuals came to. A block assembled without that record is one whose
-# design is its own fit's, so the dose it is stacked against, passed here as
-# `exposure`, answers the same question.
-ipw_numerator_fit_sigma2 <- function(model, exposure) {
+# fit's residuals came to under the estimator the family asks for. A block
+# assembled without that record is one whose design is its own fit's, so the
+# dose it is stacked against, passed here as `exposure`, answers the same
+# question, read by that same estimator.
+ipw_numerator_fit_sigma2 <- function(
+  model,
+  exposure,
+  density,
+  call = rlang::caller_env()
+) {
   if (!is.null(model$sigma2_fit)) {
     return(model$sigma2_fit)
   }
 
   fitted_n <- ipw_numerator_model_fns(model)$mean(model$X, model$coefs)
 
-  mean((exposure - fitted_n)^2)
+  ipw_continuous_sigma2_seed(exposure - fitted_n, density, call = call)
 }
 
 ipw_init_continuous <- function(spec, call = rlang::caller_env()) {
@@ -773,7 +779,12 @@ ipw_init_continuous <- function(spec, call = rlang::caller_env()) {
     coefs <- model$coefs
     stab_block <- c(
       stats::setNames(coefs, paste0("stab_", colnames(model$X))),
-      sigma2_n = ipw_numerator_fit_sigma2(model, spec$exposure)
+      sigma2_n = ipw_numerator_fit_sigma2(
+        model,
+        spec$exposure,
+        spec$density,
+        call = call
+      )
     )
   } else if (identical(spec$numerator, "marginal")) {
     a_fit <- ipw_continuous_fit_exposure(spec)
@@ -1272,7 +1283,12 @@ ipw_psi_continuous <- function(
           model$X,
           a
         ),
-        matrix((a - inputs$extras$mu_n)^2 - th_stab[[p_n + 1L]], nrow = 1)
+        density_scale_row(
+          a - inputs$extras$mu_n,
+          th_stab[[p_n + 1L]],
+          spec$density,
+          call = call
+        )
       )
     } else {
       r_a <- a - th_stab[[1]]
diff --git a/R/weights.R b/R/weights.R
index ea67b6ee..c679300d 100644
--- a/R/weights.R
+++ b/R/weights.R
@@ -154,8 +154,8 @@
 #'
 #' For a continuous exposure they are \eqn{f(A \mid V) / f(A \mid X)}: the
 #' family `.density` names, read at the numerator model's fitted mean and the
-#' root mean square of its residuals, over the same family read at the
-#' propensity score model's.
+#' spread its own residuals leave under that family, over the same family read
+#' at the propensity score model's.
 #'
 #' For a categorical exposure they are
 #' \eqn{P(Z = z_i \mid V_i) / f(z_i \mid X_i)}: the multinomial model's fitted
@@ -589,7 +589,8 @@
 #'     are \eqn{P(A = a_i \mid V_i) / f(a_i \mid X_i)} for the variables
 #'     \eqn{V} the numerator model reads. For a continuous exposure the
 #'     numerator is the family `.density` names, read at the model's fitted mean
-#'     and the root mean square of its residuals, so the weights are
+#'     and the spread its own residuals leave under that family, so the weights
+#'     are
 #'     \eqn{f(A \mid V) / f(A \mid X)}. For a categorical exposure the model is
 #'     an [nnet::multinom()] fit and the numerator is its fitted probability of
 #'     the level each unit took, read from the column named for that level, so
@@ -1368,7 +1369,7 @@ ate_continuous <- function(
   mu_n <- NULL
   sigma_n <- NULL
   if (identical(numerator, "model")) {
-    moments <- numerator_model_moments(numerator_model)
+    moments <- numerator_model_moments(numerator_model, .density, call = call)
     mu_n <- moments$mu
     sigma_n <- moments$sigma
   }
diff --git a/man/dens_normal.Rd b/man/dens_normal.Rd
index 47833c6c..1016cc60 100644
--- a/man/dens_normal.Rd
+++ b/man/dens_normal.Rd
@@ -12,7 +12,7 @@ dens_normal()
 
 dens_laplace(sigma_method = c("mle", "rms"))
 
-dens_t(df, sigma_method = c("rms", "mle"))
+dens_t(df, sigma_method = c("mle", "rms"))
 
 dens_kernel(bw = "nrd0", adjust = 1, kernel = "gaussian", n = 512)
 
@@ -72,9 +72,9 @@ Its scale is estimated under the Laplace itself, where it is the mean
 absolute residual, or by the root mean square of the residuals, which is
 the spread every family without an estimator of its own is read at.
 \item \code{dens_t()} is Student's t density with \code{df} degrees of freedom, heavier
-tailed still, and heavier the smaller \code{df} is. Its scale is estimated by
-the root mean square of the residuals, or by maximum likelihood under the
-t itself.
+tailed still, and heavier the smaller \code{df} is. Its scale is estimated
+under the t itself, by maximum likelihood, or by the root mean square of
+the residuals.
 \item \code{dens_kernel()} is a kernel density estimate of the standardized
 residuals, fit with \code{\link[stats:density]{stats::density()}} and interpolated to each
 observation. It assumes no family at all, at the cost of a density that
@@ -98,9 +98,11 @@ It is not the maximum likelihood estimator of the scale of a heavier-tailed
 family. The scale of such a family is smaller than its standard deviation,
 and the root mean square is pulled outward by the large residuals a heavy
 tail produces, which is what the family was chosen to accommodate.
-\code{sigma_method = "mle"} estimates the scale under the family itself.
-\code{dens_laplace()} takes it by default, where it is the mean absolute
-residual. For \code{dens_t()} it is the root of
+\code{sigma_method = "mle"} estimates the scale under the family itself, and both
+\code{dens_laplace()} and \code{dens_t()} take it by default, since the spread the
+density is divided by is the scale parameter of the family reading it. For
+\code{dens_laplace()} it is the mean absolute residual. For \code{dens_t()} it is the
+root of
 \eqn{\sum_i \left[(\nu + 1) r_i^2 / (\nu \sigma^2 + r_i^2)\right] = n},
 where \eqn{r_i} is the residual and \eqn{\nu} is \code{df}. Each residual enters
 that sum through a bounded term, so a residual far out in the tail moves the
diff --git a/man/wt_ate.Rd b/man/wt_ate.Rd
index 46016c8b..2c5dc0a4 100644
--- a/man/wt_ate.Rd
+++ b/man/wt_ate.Rd
@@ -338,7 +338,8 @@ is its fitted probability of the level each unit took, so the weights
 are \eqn{P(A = a_i \mid V_i) / f(a_i \mid X_i)} for the variables
 \eqn{V} the numerator model reads. For a continuous exposure the
 numerator is the family \code{.density} names, read at the model's fitted mean
-and the root mean square of its residuals, so the weights are
+and the spread its own residuals leave under that family, so the weights
+are
 \eqn{f(A \mid V) / f(A \mid X)}. For a categorical exposure the model is
 an \code{\link[nnet:multinom]{nnet::multinom()}} fit and the numerator is its fitted probability of
 the level each unit took, read from the column named for that level, so
@@ -642,8 +643,8 @@ carrying \eqn{P(A = 1 \mid V)}, which is the mistake the
 
 For a continuous exposure they are \eqn{f(A \mid V) / f(A \mid X)}: the
 family \code{.density} names, read at the numerator model's fitted mean and the
-root mean square of its residuals, over the same family read at the
-propensity score model's.
+spread its own residuals leave under that family, over the same family read
+at the propensity score model's.
 
 For a categorical exposure they are
 \eqn{P(Z = z_i \mid V_i) / f(z_i \mid X_i)}: the multinomial model's fitted
diff --git a/tests/testthat/_snaps/psw-records.md b/tests/testthat/_snaps/psw-records.md
index e049d418..1b611d72 100644
--- a/tests/testthat/_snaps/psw-records.md
+++ b/tests/testthat/_snaps/psw-records.md
@@ -12,10 +12,10 @@
       print(continuous_records_psw(.density = dens_t(4)))
     Output
       <psw{estimand = ate; stabilized}[6]>
-      [1] 0.13555778 0.10683962 0.15808856 0.12701939 0.19685396 0.07671633
+      [1] 0.13395369 0.10437312 0.16618639 0.12879842 0.20832037 0.07119641
       density:   t(df = 4)
       numerator: marginal
-      sigma:     pooled
+      sigma:     mle
     Code
       print(continuous_records_psw(stabilize = records_numerator_model()))
     Output
diff --git a/tests/testthat/_snaps/weights-continuous.md b/tests/testthat/_snaps/weights-continuous.md
index 5460af6e..9d9705af 100644
--- a/tests/testthat/_snaps/weights-continuous.md
+++ b/tests/testthat/_snaps/weights-continuous.md
@@ -188,6 +188,7 @@
       Error in `wt_ate()`:
       ! `sigma_method = "mle"` cannot be used with `.sigma`.
       x `sigma_method = "mle"` estimates the scale of the "t" conditional density from the residuals of the propensity score model, and `.sigma` is a spread of your own that replaces it.
+      i A "t" density is spread by its own estimator unless you ask for another, so `sigma_method = "mle"` is what it carries whether or not you wrote it.
       i Drop `.sigma` to estimate the scale under the "t" density, or build the density with `sigma_method = "rms"` to spread the one you supplied.
 
 # a likelihood with no maximum at a positive scale is refused
diff --git a/tests/testthat/test-ipw-continuous.R b/tests/testthat/test-ipw-continuous.R
index 0519a447..06419d25 100644
--- a/tests/testthat/test-ipw-continuous.R
+++ b/tests/testthat/test-ipw-continuous.R
@@ -2396,10 +2396,10 @@ test_that("ipw() continuous m-estimation standard errors match WeightIt", {
 
   densities <- list(
     list(ours = "normal", theirs = NULL),
-    list(ours = dens_t(df = 4), theirs = "dt_4"),
     # WeightIt standardizes by the root mean square whatever the family, so the
-    # parity is written at that estimator rather than at the scale of the
-    # Laplace, which is what the family asks for on its own.
+    # parity is written at that estimator rather than at the scale each family
+    # asks for on its own.
+    list(ours = dens_t(df = 4, sigma_method = "rms"), theirs = "dt_4"),
     list(ours = dens_laplace(sigma_method = "rms"), theirs = "dlaplace")
   )
 
@@ -2476,7 +2476,7 @@ test_that("ipw() solves the scale equation of a t spread by maximum likelihood",
   )
   fixed <- fit_continuous_models(
     dat,
-    .density = dens_t(6),
+    .density = dens_t(6, sigma_method = "rms"),
     .sigma = scale,
     stabilize = FALSE
   )
diff --git a/tests/testthat/test-weights-continuous.R b/tests/testthat/test-weights-continuous.R
index d38a9e43..ff17ac5a 100644
--- a/tests/testthat/test-weights-continuous.R
+++ b/tests/testthat/test-weights-continuous.R
@@ -172,11 +172,16 @@ test_that("Student's t weights are the t density ratio computed by hand", {
   f_den <- stats::dt(z, df = 4) / sigma
   f_num <- stats::dt(continuous_density_z_a(), df = 4) / sd_a
 
-  stabilized <- continuous_density_wt(.density = dens_t(df = 4))
+  # Both halves of the oracle are spread by the root mean square, so the family
+  # is written at that estimator rather than at the scale a t asks for on its
+  # own. What is under test here is the density, not the spread it is read at.
+  stabilized <- continuous_density_wt(
+    .density = dens_t(df = 4, sigma_method = "rms")
+  )
   expect_equal(as.numeric(stabilized), f_num / f_den, tolerance = 1e-12)
 
   unstabilized <- continuous_density_wt(
-    .density = dens_t(df = 4),
+    .density = dens_t(df = 4, sigma_method = "rms"),
     stabilize = FALSE
   )
   expect_equal(as.numeric(unstabilized), 1 / f_den, tolerance = 1e-12)
@@ -238,7 +243,9 @@ test_that("a density the user writes is evaluated as written", {
   as_function <- continuous_density_wt(
     .density = function(z) stats::dt(z, df = 4)
   )
-  as_spec <- continuous_density_wt(.density = dens_t(df = 4))
+  as_spec <- continuous_density_wt(
+    .density = dens_t(df = 4, sigma_method = "rms")
+  )
 
   expect_equal(
     as.numeric(as_function),
@@ -305,7 +312,7 @@ test_that("a supplied spread scales the density it is given to", {
   expect_equal(
     as.numeric(continuous_density_wt(
       .sigma = scalar,
-      .density = dens_t(df = 4)
+      .density = dens_t(df = 4, sigma_method = "rms")
     )),
     f_num / f_den_scalar,
     tolerance = 1e-12
@@ -317,7 +324,7 @@ test_that("a supplied spread scales the density it is given to", {
   expect_equal(
     as.numeric(continuous_density_wt(
       .sigma = per_observation,
-      .density = dens_t(df = 4)
+      .density = dens_t(df = 4, sigma_method = "rms")
     )),
     f_num / f_den_each,
     tolerance = 1e-12
@@ -328,7 +335,7 @@ test_that("a supplied spread scales the density it is given to", {
   expect_equal(
     as.numeric(continuous_density_wt(
       .sigma = per_observation,
-      .density = dens_t(df = 4),
+      .density = dens_t(df = 4, sigma_method = "rms"),
       stabilize = FALSE
     )),
     1 / f_den_each,
@@ -342,7 +349,7 @@ test_that("a stabilization score replaces the numerator of any family", {
   f_den <- stats::dt(continuous_density_z(), df = 4) / sigma
 
   weights <- continuous_density_wt(
-    .density = dens_t(df = 4),
+    .density = dens_t(df = 4, sigma_method = "rms"),
     stabilization_score = score
   )
 
@@ -1067,7 +1074,7 @@ test_that("a trimmed propensity score leaves the units it set aside missing", {
 
   families <- list(
     list(
-      input = dens_t(df = 4),
+      input = dens_t(df = 4, sigma_method = "rms"),
       g = function(values) stats::dt(values, df = 4)
     ),
     list(input = "kernel", g = continuous_density_kde)
@@ -1467,14 +1474,14 @@ test_that("the denominator is the one WeightIt divides by", {
       g = stats::dnorm(z)
     ),
     list(
-      ours = dens_t(df = 4),
+      # WeightIt standardizes by the root mean square whatever the family, so
+      # the parity is written at that estimator rather than at the scale each
+      # family asks for on its own.
+      ours = dens_t(df = 4, sigma_method = "rms"),
       theirs = "dt_4",
       g = stats::dt(z, df = 4)
     ),
     list(
-      # WeightIt standardizes by the root mean square whatever the family, so
-      # the parity is written at that estimator rather than at the scale of the
-      # Laplace, which is what the family asks for on its own.
       ours = dens_laplace(sigma_method = "rms"),
       theirs = "dlaplace",
       g = exp(-abs(z)) / 2
@@ -1635,7 +1642,7 @@ test_that("the integrated numerator is the grid marginalization by hand", {
   families <- list(
     list(input = "normal", g = function(z) stats::dnorm(z), sigma = pooled),
     list(
-      input = dens_t(df = 4),
+      input = dens_t(df = 4, sigma_method = "rms"),
       g = function(z) stats::dt(z, df = 4),
       sigma = pooled
     ),
@@ -1908,7 +1915,7 @@ test_that("the integrated numerator reads only the units with a residual", {
   # those units make on their own, and a unit with no residual has no weight.
   families <- list(
     list(
-      input = dens_t(df = 4),
+      input = dens_t(df = 4, sigma_method = "rms"),
       oracle = function(exposure, mu) {
         continuous_integrated_wt(
           function(z) stats::dt(z, df = 4),
@@ -2311,7 +2318,7 @@ test_that("the numerator reaches the weights through a fitted model", {
     exposure,
     exposure_type = "continuous",
     stabilize = TRUE,
-    .density = dens_t(df = 4),
+    .density = dens_t(df = 4, sigma_method = "rms"),
     numerator = "integrated"
   )
 
@@ -2347,7 +2354,7 @@ test_that("censoring weights take an integrated numerator of their own", {
     continuous_density_data$exposure,
     exposure_type = "continuous",
     stabilize = TRUE,
-    .density = dens_t(df = 4),
+    .density = dens_t(df = 4, sigma_method = "rms"),
     numerator = "integrated"
   )
 
@@ -2407,7 +2414,7 @@ test_that("the numerator reaches the weights through a modified score", {
       problem$exposure,
       exposure_type = "continuous",
       stabilize = TRUE,
-      .density = dens_t(df = 4),
+      .density = dens_t(df = 4, sigma_method = "rms"),
       numerator = "integrated"
     )
 
@@ -2453,7 +2460,7 @@ test_that("censoring weights carry a numerator through their other methods", {
     continuous_density_data$exposure,
     exposure_type = "continuous",
     stabilize = TRUE,
-    .density = dens_t(df = 4),
+    .density = dens_t(df = 4, sigma_method = "rms"),
     numerator = "integrated"
   )
 
@@ -2486,7 +2493,7 @@ test_that("censoring weights carry a numerator through a modified score", {
       problem$exposure,
       exposure_type = "continuous",
       stabilize = TRUE,
-      .density = dens_t(df = 4),
+      .density = dens_t(df = 4, sigma_method = "rms"),
       numerator = "integrated"
     )
 
@@ -2526,10 +2533,10 @@ test_that("integrated weights are the weights WeightIt gives", {
   # the scale ever drifts.
   families <- list(
     list(ours = "normal", theirs = NULL),
-    list(ours = dens_t(df = 4), theirs = "dt_4"),
     # WeightIt standardizes by the root mean square whatever the family, so the
-    # parity is written at that estimator rather than at the scale of the
-    # Laplace, which is what the family asks for on its own.
+    # parity is written at that estimator rather than at the scale each family
+    # asks for on its own.
+    list(ours = dens_t(df = 4, sigma_method = "rms"), theirs = "dt_4"),
     list(ours = dens_laplace(sigma_method = "rms"), theirs = "dlaplace"),
     list(ours = "kernel", theirs = "kernel")
   )
@@ -2626,13 +2633,16 @@ test_that("the t scale by maximum likelihood is the root of the t score", {
   scale <- t_scale_mle(continuous_t_residuals, df = 3)
 
   mle <- continuous_t_wt(.density = dens_t(3, sigma_method = "mle"))
-  oracle <- continuous_t_wt(.density = dens_t(3), .sigma = scale)
+  oracle <- continuous_t_wt(
+    .density = dens_t(3, sigma_method = "rms"),
+    .sigma = scale
+  )
 
   expect_equal(as.numeric(mle), as.numeric(oracle), tolerance = 1e-8)
 })
 
 test_that("a maximum likelihood scale parts from the root mean square in heavy tails", {
-  rms <- continuous_t_wt(.density = dens_t(3))
+  rms <- continuous_t_wt(.density = dens_t(3, sigma_method = "rms"))
   mle <- continuous_t_wt(.density = dens_t(3, sigma_method = "mle"))
 
   relative <- abs(as.numeric(mle) / as.numeric(rms) - 1)
@@ -2653,7 +2663,10 @@ test_that("the two spreads meet as the t approaches the normal", {
 
   expect_equal(scale, continuous_density_pooled(), tolerance = 0.005)
 
-  rms <- continuous_density_wt(stabilize = FALSE, .density = dens_t(1000))
+  rms <- continuous_density_wt(
+    stabilize = FALSE,
+    .density = dens_t(1000, sigma_method = "rms")
+  )
   mle <- continuous_density_wt(
     stabilize = FALSE,
     .density = dens_t(1000, sigma_method = "mle")
@@ -2743,7 +2756,7 @@ test_that("a maximum likelihood scale spreads the density it is read under", {
     .sigma = scale,
     exposure_type = "continuous",
     stabilize = FALSE,
-    .density = dens_t(6)
+    .density = dens_t(6, sigma_method = "rms")
   )
 
   expect_equal(as.numeric(mle), as.numeric(oracle), tolerance = 1e-8)
@@ -2783,7 +2796,7 @@ test_that("a maximum likelihood scale is fit on the residuals that are there", {
   )
   oracle <- continuous_t_wt(
     exposure = exposure,
-    .density = dens_t(4),
+    .density = dens_t(4, sigma_method = "rms"),
     .sigma = scale
   )
 
@@ -2803,9 +2816,13 @@ test_that("a spread supplied and a spread estimated under the density are refuse
     continuous_t_wt(.density = dens_t(4, sigma_method = "mle"), .sigma = 0.9)
   )
 
-  # The default estimator takes a supplied spread as every other family does, so
-  # what is refused is the pairing rather than the family.
-  expect_no_error(continuous_t_wt(.density = dens_t(4), .sigma = 0.9))
+  # A family asked for the root mean square takes a supplied spread as every
+  # family without an estimator of its own does, so what is refused is the
+  # pairing rather than the family.
+  expect_no_error(continuous_t_wt(
+    .density = dens_t(4, sigma_method = "rms"),
+    .sigma = 0.9
+  ))
 })
 
 # A propensity score model that all but interpolates its exposure: every unit
@@ -2843,7 +2860,7 @@ test_that("a maximum likelihood scale far below the root mean square is found",
     exposure,
     exposure_type = "continuous",
     stabilize = FALSE,
-    .density = dens_t(4),
+    .density = dens_t(4, sigma_method = "rms"),
     .sigma = scale
   )
 

From 0df85221cab1726c3b37af67a468853cab9f3a2b Mon Sep 17 00:00:00 2001
From: Malcolm Barrett <malcolmbarrett@gmail.com>
Date: Wed, 16 Sep 2026 13:25:12 -0400
Subject: [PATCH 5/7] add a Laplace bootstrap arm and name the estimator the
 heavy-tailed block tests

---
 tests/testthat/test-ipw-continuous.R | 47 ++++++++++++++++++++++++++--
 1 file changed, 45 insertions(+), 2 deletions(-)

diff --git a/tests/testthat/test-ipw-continuous.R b/tests/testthat/test-ipw-continuous.R
index 06419d25..762212fa 100644
--- a/tests/testthat/test-ipw-continuous.R
+++ b/tests/testthat/test-ipw-continuous.R
@@ -2244,9 +2244,14 @@ test_that("the sandwich SE for a heavy-tailed integrated ratio tracks a bootstra
   skip_on_cran()
 
   dat <- sim_continuous(seed = 2024, n = 600)
+  # Written against the root mean square, which was the t's default when this
+  # block was first set down. The estimator is named here rather than left to
+  # the default so the block says which row it holds against the bootstrap: at
+  # `sigma_method = "mle"` the stacked system solves the t's own score for the
+  # scale, not the moment equation.
   mods <- fit_continuous_models(
     dat,
-    .density = dens_t(4),
+    .density = dens_t(4, sigma_method = "mle"),
     numerator = "integrated"
   )
   mest_se <- ipw(mods$ps_mod, mods$outcome_mod)$estimates$std.err
@@ -2256,7 +2261,7 @@ test_that("the sandwich SE for a heavy-tailed integrated ratio tracks a bootstra
     w <- continuous_weights(
       as.double(fitted(ps)),
       d$A,
-      .density = dens_t(4),
+      .density = dens_t(4, sigma_method = "mle"),
       numerator = "integrated"
     )
     msm <- lm(yc ~ A, data = d, weights = as.double(w))
@@ -2278,6 +2283,44 @@ test_that("the sandwich SE for a heavy-tailed integrated ratio tracks a bootstra
   expect_lt(abs(mest_se - boot_se) / boot_se, 0.15)
 })
 
+test_that("the sandwich SE for a Laplace ratio tracks a bootstrap", {
+  skip_on_cran()
+
+  # The Laplace scale is the one estimator no bootstrap held the sandwich
+  # against, and it is the arm the simulation found the largest coverage error
+  # on when the scale was read by the moment instead. Its row is a different
+  # shape from the t's, `abs(r) - sqrt(sigma2)` rather than a score, and it is
+  # stacked twice, once for the conditional scale and once for the marginal.
+  dat <- sim_continuous(seed = 2024, n = 600)
+  mods <- fit_continuous_models(dat, .density = dens_laplace())
+  mest_se <- ipw(mods$ps_mod, mods$outcome_mod)$estimates$std.err
+
+  boot_slope <- function(d) {
+    ps <- lm(A ~ x1 + x2, data = d)
+    w <- continuous_weights(
+      as.double(fitted(ps)),
+      d$A,
+      .density = dens_laplace()
+    )
+    msm <- lm(yc ~ A, data = d, weights = as.double(w))
+    unname(coef(msm)[["A"]])
+  }
+
+  withr::local_seed(918)
+  reps <- 400L
+  n <- nrow(dat)
+  boot <- vapply(
+    seq_len(reps),
+    function(i) {
+      boot_slope(dat[sample.int(n, n, replace = TRUE), , drop = FALSE])
+    },
+    numeric(1)
+  )
+
+  boot_se <- stats::sd(boot)
+  expect_lt(abs(mest_se - boot_se) / boot_se, 0.15)
+})
+
 test_that("the sandwich SE for a log-link propensity model tracks a bootstrap", {
   skip_on_cran()
 

From d4c6a64ee93bbcda5818f9b97358565547f41d0f Mon Sep 17 00:00:00 2001
From: Malcolm Barrett <malcolmbarrett@gmail.com>
Date: Wed, 16 Sep 2026 13:56:45 -0400
Subject: [PATCH 6/7] document the scale estimator each density family reads
 its spread at

---
 NEWS.md                  | 60 +++++++++++++++++++++-------------------
 R/density.R              |  6 ++--
 R/ipw.R                  | 35 ++++++++++++-----------
 R/weights.R              | 23 +++++++++------
 man/dens_normal.Rd       |  6 ++--
 man/ipw-methods.Rd       | 35 ++++++++++++-----------
 man/wt_ate.Rd            | 23 +++++++++------
 vignettes/propensity.Rmd |  5 +++-
 8 files changed, 109 insertions(+), 84 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index 8ba78d86..8a4ccfef 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -572,7 +572,7 @@
   reports when an outcome is constant in one exposure group only; the
   weights-mismatch report no longer offers the focal level as a cause on the
   routes that resolve none; and the spread of a conditional density is named as
-  the pooled residual root mean square, which is what the package computes.
+  the spread its family estimates, which is what the package computes.
 
 * Weights for a continuous exposure now read a one-column matrix and a
   one-dimensional array of conditional means as the vector of one mean per unit
@@ -764,16 +764,16 @@
   `lm` method, which a `MASS::rlm()` reaches by inheritance, and their `glm`
   method reads a `gaussian()` fit under any of its links as well as an
   `mgcv::gam()` fit with it. Each of these classes reports its conditional mean
-  on the scale of the exposure, so a log, inverse, or square root link never has
-  to be undone, and every one of them uses the same pooled residual spread.
-  `rlm` reports a robust scale estimate of its own in `fit$s`, which resists the
-  extreme residuals rather than pooling all of them, and is used only when it is
-  asked for with `.sigma = fit$s`. A family whose spread changes with its fitted
-  values, such as `poisson()` or `quasi(variance = "mu")`, describes a different
-  density for every unit, which a single spread cannot stand in for, and is
-  refused with an error of class `propensity_model_family_error`.
-  `quasi(variance = "constant")` is the gaussian variance under another name and
-  is accepted.
+  on the scale of the exposure, so a log, inverse, or square root link never
+  has to be undone, and every one of them is spread by the same estimator.
+  `rlm` reports a robust scale estimate of its own in `fit$s`, which resists
+  the extreme residuals rather than pooling all of them, and is used only when
+  it is asked for with `.sigma = fit$s`. A family whose spread changes with its
+  fitted values, such as `poisson()` or `quasi(variance = "mu")`, describes a
+  different density for every unit, which a single spread cannot stand in for,
+  and is refused with an error of class `propensity_model_family_error`.
+  `quasi(variance = "constant")` is the gaussian variance under another name
+  and is accepted.
 
   The binary path is now held to the same rule from the other side, which
   changes what some calls return. Only `binomial()` and `quasibinomial()` fit
@@ -792,18 +792,20 @@
   `NULL` for weights built without an exposure to describe, such as those
   written with `psw()` directly. `density_meta()` describes the ratio a
   continuous exposure's weights are, and returns `NULL` for weights that are no
-  such ratio. Its elements are `density`, the specification of the
-  conditional density family; `numerator`, what stabilized the weights, which is
+  such ratio. Its elements are `density`, the specification of the conditional
+  density family; `numerator`, what stabilized the weights, which is
   `"marginal"` for the marginal density of the exposure, `"integrated"` for the
   conditional density marginalized over the units, `"score"` for a
   `stabilization_score` the caller supplied, and `"none"` for weights that were
-  not stabilized; `sigma`, where the residual spread came from, either
-  `"pooled"` or `"supplied"`; and `sigma_value`, the number a single supplied
-  spread was, which is `NULL` for a pooled spread and for one supplied per
-  observation. A spread that is one number is a constant the weights can be
-  rebuilt from, which is what `ipw()` needs of it; a spread that changes with
-  the observation is not, so the record holds where it came from and nothing
-  more.
+  not stabilized; `sigma`, where the residual spread came from, which is
+  `"pooled"` for the pooled residual root mean square, `"mle"` for a scale
+  estimated under the family that reads it, and `"supplied"` for one the caller
+  gave; and `sigma_value`, the number a single supplied spread was, which is
+  `NULL` for a spread estimated from the residuals, by either estimator, and
+  for one supplied per observation. A spread that is one number is a constant
+  the weights can be rebuilt from, which is what `ipw()` needs of it; a spread
+  that changes with the observation is not, so the record holds where it came
+  from and nothing more.
 
   Both records describe the exposure rather than the units, so neither holds a
   length of its own and both survive subsetting, arithmetic, and anything else
@@ -851,15 +853,15 @@
   than the sum of squares, so the stacked system carries the Huber score its
   coefficients are the root of, clipped where the fit itself clipped: at the
   psi's own constant, including one a caller passed as `k`, times the scale the
-  fit settled on. That scale enters the score as a known constant whose sampling
-  variability is not propagated, and the spread of the conditional density is
-  the pooled residual spread, as it is for every other class. A fit with a psi
-  other than Huber, and one fit with `method = "MM"`, which finishes on a
-  redescending psi, are refused with an error of class
-  `propensity_ipw_robust_psi_error`; one whose iteration stopped short of its
-  own tolerance is refused with `propensity_ipw_convergence_error`. The first
-  two point to a Huber refit or to a bootstrap of the whole fit written by hand;
-  the last points to a larger `maxit` or a looser `acc`.
+  fit settled on. That scale enters the score as a known constant whose
+  sampling variability is not propagated, and the spread of the conditional
+  density comes from the estimator its family asks for, there as for every
+  other class. A fit with a psi other than Huber, and one fit with `method =
+  "MM"`, which finishes on a redescending psi, are refused with an error of
+  class `propensity_ipw_robust_psi_error`; one whose iteration stopped short of
+  its own tolerance is refused with `propensity_ipw_convergence_error`. The
+  first two point to a Huber refit or to a bootstrap of the whole fit written
+  by hand; the last points to a larger `maxit` or a looser `acc`.
 
 * `wt_joint()` now reads each component's exposure type off the component
   rather than being told it. A weight function records the type of exposure it
diff --git a/R/density.R b/R/density.R
index 7f3a9cf8..108b4f99 100644
--- a/R/density.R
+++ b/R/density.R
@@ -70,8 +70,10 @@ new_density_spec <- function(
 #' @param df Degrees of freedom for Student's t, a single positive, finite
 #'   number.
 #' @param sigma_method How the spread of the conditional density is estimated
-#'   from the residuals of the propensity score model, either `"rms"` or
-#'   `"mle"`. See the section below.
+#'   from the residuals of the propensity score model, either `"rms"`, the root
+#'   mean square, or `"mle"`, the scale estimated under the family itself, which
+#'   both `dens_laplace()` and `dens_t()` take by default. See the section
+#'   below.
 #' @param bw The bandwidth passed to [stats::density()]: a single positive
 #'   number, or the name of one of its selection rules
 #'   (`"nrd0"`, `"nrd"`, `"ucv"`, `"bcv"`, `"SJ"`, `"SJ-ste"`, or `"SJ-dpi"`).
diff --git a/R/ipw.R b/R/ipw.R
index d7c38c79..0baf7ac0 100644
--- a/R/ipw.R
+++ b/R/ipw.R
@@ -1083,22 +1083,25 @@
 #' existed.
 #'
 #' A numerator estimated by a model the caller passed to [wt_ate()]'s
-#' `stabilize` contributes that model: one parameter per coefficient, named with
-#' a `stab_` prefix so the terms it shares with the propensity score model stay
-#' apart from it, and one for the spread its density is read at, `sigma2_n`. Its
-#' score and the moment its spread is the root of are stacked alongside the
-#' propensity score model's, so the standard errors account for the numerator
-#' having been fitted, which is what separates this from a
-#' `stabilization_score`: a score the caller computed is carried as a known
-#' constant and contributes no parameter at all. The numerator model is read
-#' through the registry above, so a class, family, or link whose score this
-#' system cannot write is refused there in the terms that registry refuses a
-#' propensity score model in, naming `stabilize`. A model of a response other
-#' than the exposure, or one fit to a different set of observations, errors with
-#' class `propensity_ipw_numerator_error`. One fit with case weights errors with
-#' class `propensity_ipw_ps_weights_error`, for the reason a weighted propensity
-#' score model does: the block written for it is its unweighted score, and its
-#' coefficients are not at the root of that.
+#' `stabilize` contributes that model: one parameter per coefficient, named
+#' with a `stab_` prefix so the terms it shares with the propensity score model
+#' stay apart from it, and one for the spread its density is read at,
+#' `sigma2_n`. That parameter holds a squared scale rather than a variance,
+#' which are the same number only for a family whose scale is its standard
+#' deviation. Its score, and the equation that spread is the root of, which is
+#' the residual moment or the family's own score depending on what the density
+#' asks for, are stacked alongside the propensity score model's, so the
+#' standard errors account for the numerator having been fitted, which is what
+#' separates this from a `stabilization_score`: a score the caller computed is
+#' carried as a known constant and contributes no parameter at all. The
+#' numerator model is read through the registry above, so a class, family, or
+#' link whose score this system cannot write is refused there in the terms that
+#' registry refuses a propensity score model in, naming `stabilize`. A model of
+#' a response other than the exposure, or one fit to a different set of
+#' observations, errors with class `propensity_ipw_numerator_error`. One fit
+#' with case weights errors with class `propensity_ipw_ps_weights_error`, for
+#' the reason a weighted propensity score model does: the block written for it
+#' is its unweighted score, and its coefficients are not at the root of that.
 #'
 #' A dose numerator model's design is one of the designs `ipw()` rebuilds when
 #' `.data` is supplied, on the terms a discrete exposure's is: it is rebuilt
diff --git a/R/weights.R b/R/weights.R
index c679300d..e92f716f 100644
--- a/R/weights.R
+++ b/R/weights.R
@@ -355,15 +355,20 @@
 #' a conditional mean is refused there whatever its fitted values happen to fall
 #' on, a linear probability model included.
 #'
-#' \eqn{\sigma} is the pooled residual spread
-#' \eqn{\sqrt{\mathrm{mean}((A - \hat{A})^2)}} unless `.sigma` supplies a
-#' single standard deviation, or one for each unit, which the model methods do
-#' not do on their own. [dens_t()] and [dens_laplace()] offer a second
-#' estimator, the maximum likelihood scale of the family itself, through
-#' `sigma_method = "mle"`; see **The spread of the conditional density** in
-#' [dens_t()]'s documentation. A family that estimates its own scale spreads
-#' both densities of the ratio, so the marginal density is read at that scale
-#' too and the two halves are densities of one width rather than two.
+#' \eqn{\sigma} is the spread the family `.density` names asks for, read from
+#' the residuals of the propensity score model. For a family with no scale
+#' estimator of its own that is the pooled root mean square
+#' \eqn{\sqrt{\mathrm{mean}((A - \hat{A})^2)}}; for [dens_t()] and
+#' [dens_laplace()] it is the maximum likelihood scale of the family itself,
+#' which both take by default. `sigma_method = "rms"` asks either of them for
+#' the root mean square instead; see **The spread of the conditional density**
+#' in [dens_t()]'s documentation. `.sigma` supplies a spread of your own, a
+#' single standard deviation or one for each unit, which the model methods do
+#' not do on their own, and is refused alongside a family that estimates its
+#' own scale, being a second instruction about the same quantity. A family that
+#' estimates its own scale spreads both densities of the ratio, so the marginal
+#' density is read at that scale too and the two halves are densities of one
+#' width rather than two.
 #'
 #' Every model class is spread that same way, [MASS::rlm()] included. `rlm`
 #' reports a robust scale estimate of its own in `fit$s`, which resists the
diff --git a/man/dens_normal.Rd b/man/dens_normal.Rd
index 1016cc60..73735796 100644
--- a/man/dens_normal.Rd
+++ b/man/dens_normal.Rd
@@ -20,8 +20,10 @@ dens_fn(f)
 }
 \arguments{
 \item{sigma_method}{How the spread of the conditional density is estimated
-from the residuals of the propensity score model, either \code{"rms"} or
-\code{"mle"}. See the section below.}
+from the residuals of the propensity score model, either \code{"rms"}, the root
+mean square, or \code{"mle"}, the scale estimated under the family itself, which
+both \code{dens_laplace()} and \code{dens_t()} take by default. See the section
+below.}
 
 \item{df}{Degrees of freedom for Student's t, a single positive, finite
 number.}
diff --git a/man/ipw-methods.Rd b/man/ipw-methods.Rd
index 07a0a942..3bdf43f7 100644
--- a/man/ipw-methods.Rd
+++ b/man/ipw-methods.Rd
@@ -1349,22 +1349,25 @@ by hand, are read as the normal ratio the package built before the record
 existed.
 
 A numerator estimated by a model the caller passed to \code{\link[=wt_ate]{wt_ate()}}'s
-\code{stabilize} contributes that model: one parameter per coefficient, named with
-a \code{stab_} prefix so the terms it shares with the propensity score model stay
-apart from it, and one for the spread its density is read at, \code{sigma2_n}. Its
-score and the moment its spread is the root of are stacked alongside the
-propensity score model's, so the standard errors account for the numerator
-having been fitted, which is what separates this from a
-\code{stabilization_score}: a score the caller computed is carried as a known
-constant and contributes no parameter at all. The numerator model is read
-through the registry above, so a class, family, or link whose score this
-system cannot write is refused there in the terms that registry refuses a
-propensity score model in, naming \code{stabilize}. A model of a response other
-than the exposure, or one fit to a different set of observations, errors with
-class \code{propensity_ipw_numerator_error}. One fit with case weights errors with
-class \code{propensity_ipw_ps_weights_error}, for the reason a weighted propensity
-score model does: the block written for it is its unweighted score, and its
-coefficients are not at the root of that.
+\code{stabilize} contributes that model: one parameter per coefficient, named
+with a \code{stab_} prefix so the terms it shares with the propensity score model
+stay apart from it, and one for the spread its density is read at,
+\code{sigma2_n}. That parameter holds a squared scale rather than a variance,
+which are the same number only for a family whose scale is its standard
+deviation. Its score, and the equation that spread is the root of, which is
+the residual moment or the family's own score depending on what the density
+asks for, are stacked alongside the propensity score model's, so the
+standard errors account for the numerator having been fitted, which is what
+separates this from a \code{stabilization_score}: a score the caller computed is
+carried as a known constant and contributes no parameter at all. The
+numerator model is read through the registry above, so a class, family, or
+link whose score this system cannot write is refused there in the terms that
+registry refuses a propensity score model in, naming \code{stabilize}. A model of
+a response other than the exposure, or one fit to a different set of
+observations, errors with class \code{propensity_ipw_numerator_error}. One fit
+with case weights errors with class \code{propensity_ipw_ps_weights_error}, for
+the reason a weighted propensity score model does: the block written for it
+is its unweighted score, and its coefficients are not at the root of that.
 
 A dose numerator model's design is one of the designs \code{ipw()} rebuilds when
 \code{.data} is supplied, on the terms a discrete exposure's is: it is rebuilt
diff --git a/man/wt_ate.Rd b/man/wt_ate.Rd
index 2c5dc0a4..20669486 100644
--- a/man/wt_ate.Rd
+++ b/man/wt_ate.Rd
@@ -849,15 +849,20 @@ exposure is held to the same rule from the other side: only \code{\link[=binomia
 a conditional mean is refused there whatever its fitted values happen to fall
 on, a linear probability model included.
 
-\eqn{\sigma} is the pooled residual spread
-\eqn{\sqrt{\mathrm{mean}((A - \hat{A})^2)}} unless \code{.sigma} supplies a
-single standard deviation, or one for each unit, which the model methods do
-not do on their own. \code{\link[=dens_t]{dens_t()}} and \code{\link[=dens_laplace]{dens_laplace()}} offer a second
-estimator, the maximum likelihood scale of the family itself, through
-\code{sigma_method = "mle"}; see \strong{The spread of the conditional density} in
-\code{\link[=dens_t]{dens_t()}}'s documentation. A family that estimates its own scale spreads
-both densities of the ratio, so the marginal density is read at that scale
-too and the two halves are densities of one width rather than two.
+\eqn{\sigma} is the spread the family \code{.density} names asks for, read from
+the residuals of the propensity score model. For a family with no scale
+estimator of its own that is the pooled root mean square
+\eqn{\sqrt{\mathrm{mean}((A - \hat{A})^2)}}; for \code{\link[=dens_t]{dens_t()}} and
+\code{\link[=dens_laplace]{dens_laplace()}} it is the maximum likelihood scale of the family itself,
+which both take by default. \code{sigma_method = "rms"} asks either of them for
+the root mean square instead; see \strong{The spread of the conditional density}
+in \code{\link[=dens_t]{dens_t()}}'s documentation. \code{.sigma} supplies a spread of your own, a
+single standard deviation or one for each unit, which the model methods do
+not do on their own, and is refused alongside a family that estimates its
+own scale, being a second instruction about the same quantity. A family that
+estimates its own scale spreads both densities of the ratio, so the marginal
+density is read at that scale too and the two halves are densities of one
+width rather than two.
 
 Every model class is spread that same way, \code{\link[MASS:rlm]{MASS::rlm()}} included. \code{rlm}
 reports a robust scale estimate of its own in \code{fit$s}, which resists the
diff --git a/vignettes/propensity.Rmd b/vignettes/propensity.Rmd
index 3e41720e..1b2b7495 100644
--- a/vignettes/propensity.Rmd
+++ b/vignettes/propensity.Rmd
@@ -300,7 +300,10 @@ density_meta(wts_dose)
 
 The default reads that ratio in the normal family, which is a strong claim
 about the residuals of a dose model. `.density` chooses another family, and a
-heavier tail holds down the weight of a unit whose dose the model fits poorly:
+heavier tail holds down the weight of a unit whose dose the model fits poorly.
+A family is a model of the residual distribution, so its scale is estimated
+under that family rather than by the root mean square, which is the scale of a
+normal alone; `dens_t()` and `dens_laplace()` do that by default:
 
 ```{r continuous-density}
 wts_t <- wt_ate(ps_dose, .density = dens_t(df = 4))

From f6413a1ed1f32187681a7372b33475cdd889ba3e Mon Sep 17 00:00:00 2001
From: Malcolm Barrett <malcolmbarrett@gmail.com>
Date: Wed, 16 Sep 2026 15:00:21 -0400
Subject: [PATCH 7/7] record the outcome model formula itself so the printed
 call is stable across R versions

---
 tests/testthat/_snaps/ipw-categorical.md | 4 ++--
 tests/testthat/_snaps/ipw-continuous.md  | 2 +-
 tests/testthat/test-ipw-categorical.R    | 8 +++++++-
 tests/testthat/test-ipw-continuous.R     | 3 +++
 4 files changed, 13 insertions(+), 4 deletions(-)

diff --git a/tests/testthat/_snaps/ipw-categorical.md b/tests/testthat/_snaps/ipw-categorical.md
index 221e625d..bea52c39 100644
--- a/tests/testthat/_snaps/ipw-categorical.md
+++ b/tests/testthat/_snaps/ipw-categorical.md
@@ -21,8 +21,8 @@
           reltol = 1e-14, maxit = 2000) 
       
       Outcome Model:
-        Call: glm(formula = fmla, family = quasibinomial(), data = dat, weights = wts, 
-          control = glm.control(epsilon = 1e-14, maxit = 200)) 
+        Call: glm(formula = y ~ a + x1, family = quasibinomial(), data = dat, 
+          weights = wts, control = glm.control(epsilon = 1e-14, maxit = 200)) 
       
       Marginal estimates:
                      estimate std.err     z ci.lower ci.upper conf.level p.value    
diff --git a/tests/testthat/_snaps/ipw-continuous.md b/tests/testthat/_snaps/ipw-continuous.md
index 1820b836..95ee327e 100644
--- a/tests/testthat/_snaps/ipw-continuous.md
+++ b/tests/testthat/_snaps/ipw-continuous.md
@@ -11,7 +11,7 @@
         Call: lm(formula = A ~ x1 + x2, data = dat) 
       
       Outcome Model:
-        Call: lm(formula = msm_fmla, data = dat, weights = wts) 
+        Call: lm(formula = yc ~ A, data = dat, weights = wts) 
       
       Marginal estimates:
             estimate  std.err      z ci.lower ci.upper conf.level   p.value    
diff --git a/tests/testthat/test-ipw-categorical.R b/tests/testthat/test-ipw-categorical.R
index 05f5f133..d1b0c99f 100644
--- a/tests/testthat/test-ipw-categorical.R
+++ b/tests/testthat/test-ipw-categorical.R
@@ -108,7 +108,7 @@ fit_outcome <- function(
   outcome_var <- if (outcome_family == "binomial") "y" else "yc"
   rhs <- if (covariates) c("a", "x1") else "a"
   fmla <- stats::reformulate(rhs, response = outcome_var)
-  if (outcome_family == "binomial") {
+  mod <- if (outcome_family == "binomial") {
     glm(
       fmla,
       data = dat,
@@ -119,6 +119,12 @@ fit_outcome <- function(
   } else {
     lm(fmla, data = dat, weights = wts)
   }
+  # Record the formula itself rather than the name that held it. R-devel
+  # substitutes the formula into the call it stores, where released versions
+  # keep the symbol, so a printed call built from a variable is not the same
+  # text on every R version. Normalizing it here makes the snapshot stable.
+  mod$call$formula <- fmla
+  mod
 }
 
 # Build the propensity score model, the categorical weights, and the weighted
diff --git a/tests/testthat/test-ipw-continuous.R b/tests/testthat/test-ipw-continuous.R
index 762212fa..6405c865 100644
--- a/tests/testthat/test-ipw-continuous.R
+++ b/tests/testthat/test-ipw-continuous.R
@@ -136,6 +136,9 @@ fit_continuous_models <- function(
   } else {
     outcome_mod <- lm(msm_fmla, data = dat, weights = wts)
   }
+  # As in fit_outcome() in test-ipw-categorical.R: record the formula itself so
+  # the printed call reads the same on every R version.
+  outcome_mod$call$formula <- msm_fmla
 
   list(ps_mod = ps_mod, outcome_mod = outcome_mod, wts = wts)
 }