library(ggplot2)
te_pal <- list(forest = "#275139", green = "#2f8f63", sage = "#93a87f",
clay = "#b5534e", gold = "#cda23f", line = "#dad9ca",
ink = "#16241d", paper = "#f5f4ee")
theme_te <- function() {
theme_minimal(base_size = 12) +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_line(colour = "#e7e6dc"),
plot.background = element_rect(fill = "#f5f4ee", colour = NA),
panel.background = element_rect(fill = "#f5f4ee", colour = NA),
plot.title = element_text(face = "bold", colour = te_pal$ink),
axis.title = element_text(colour = "#2c3a31"),
axis.text = element_text(colour = "#2c3a31"))
}Transporting an effect to a new region
R
causal inference
transportability
external validity
ecology tutorial
An effect measured in one valley is wanted for a whole region. Measuring the sign flip, the transport formula, the precision it costs and the two failures.
Corrected 14 October 2026: the section on an unmeasured modifier set the error that transport left behind against the gap between the two truths of the one valley analysed first; it now sets it, there and in The honest limit, against how far transport moved the estimate averaged over the same replicate valleys, and the share it reports is a little smaller.
Corrected 1 October 2026: the coverage figure drew its dashed valley truth line at the truth of the one valley analysed first, to the left of the replicate estimates; it now marks the valley truth averaged over the simulated valleys, the quantity the bias and coverage are measured against.
Corrected 27 September 2026: the post said its two effect functions were identical inside the valley’s range. They agree over the dense range but part company above thirty per cent moisture, which the valley’s wettest plots still reach, and the sentence now says so.
Updated 25 September 2026: the post named the augmented transport estimator as the natural next step. A new section, The augmented transport estimator, builds it and measures it on both targets and both truths, and Where to go next now points back to it.
A grazing exclusion trial in one valley: three hundred and twenty permanent plots on semi-natural grassland, each either fenced against sheep or left open, with vascular plant richness recorded in a four square metre quadrat after six seasons. The plots were not fenced at random. The wetter ground at the valley floor was easier to secure and got fenced more often than the dry upper slopes, so soil moisture drives both the treatment and the outcome, and the raw comparison is confounded. Soil moisture was measured on every plot, so the confounding is the kind that adjustment handles.
The result goes into a report, and the report is read by an agency that regulates grazing across the whole upland district: forty thousand mapped grassland parcels, moisture known for every one of them from a soil survey and a wetness index, and no outcome data anywhere except in that one valley. The question the agency asks is not the one the trial answered. The trial answered what fencing did to the plots in the valley; the agency wants to know what a district-wide prescription would do to the district.
Those two numbers are not the same, and the gap is not sampling error. Herbivore removal does different things on different ground: Bakker et al (2006) found across grasslands that grazing raises plant diversity where habitats are productive and lowers it where they are not, so the sign of a grazing effect depends on where you stand. The study valley is dry, and the district includes a great deal of wet, productive ground that the valley barely represents. Averaging one effect function over two different moisture distributions gives two different answers, and the difference can be larger than the effect.
This post is positioned directly against G-computation and standardisation. That post fits an outcome model, predicts the potential outcome for every unit under each treatment level, averages, and takes the difference. Its own code names the target in a comment: the true SAMPLE-marginal effect. The averaging runs over the covariate distribution of the sample in hand, which is right if that sample is the population you care about, and nothing there asks whether it is. Here the machinery is nearly identical and one thing changes: the averaging runs over a different population’s covariate distribution. The arithmetic is a line of code apart. The assumption behind it is not, and that assumption has no test.
The other neighbour is Checking causal assumptions, which stresses the internal-validity premise with negative controls and a robustness value. Everything it checks is about whether the estimate is right for the units that were studied. Nothing in it, and nothing else in this blog’s causal cluster, asks whether those units are the units the decision is about. That second question is what follows.
Seven things get measured: the effect in the valley, what the valley’s interval covers when the district is the target, the transport formula run two ways, the precision it costs, a diagnostic that warns early, two ways the whole thing fails while looking healthy, and what the augmented estimator that combines both routes is worth on each target.
Two populations, one effect function
The generating model has one covariate doing three jobs. Soil moisture confounds, because it drives fencing; it is prognostic, because wetter grassland carries more species anyway; and it modifies, because the sign of the fencing effect depends on it. Write \(W\) for moisture as per cent volumetric water content, \(A\) for the fencing indicator and \(Y\) for richness in the quadrat:
\[E[Y \mid A, W] = \mu(W) + \tau(W)\,A, \qquad \tau(W) = \tau_0 + \tau_1 W\]
The single decision that makes transport possible is that \(\tau(\cdot)\) is the same function in the valley and in the district: two plots with the same soil moisture respond to fencing the same way, wherever they are. That is the transportability assumption in full, and the simulation is built so that it holds. What differs between the two places is only the distribution of \(W\).
w_lo <- 4; w_hi <- 46; n_study <- 320; n_target <- 40000; sd_noise <- 1.6
tau_0 <- 2.6; tau_1 <- -0.10; zero_at <- -tau_0 / tau_1
tau_of <- function(w) tau_0 + tau_1 * w
mu_of <- function(w) 9 + 0.22 * w
prop_of <- function(w) 1 / (1 + exp(-(-0.4 + 0.08 * (w - 18))))
set.seed(20260802)
w_target <- w_lo + (w_hi - w_lo) * rbeta(n_target, 7.235, 6.515)
study_of <- function(seed) {
set.seed(seed)
w <- w_lo + (w_hi - w_lo) * rbeta(n_study, 2.2, 4.4)
a <- rbinom(n_study, 1, prop_of(w))
data.frame(y = mu_of(w) + tau_of(w) * a + rnorm(n_study, 0, sd_noise),
a = a, w = w)
}
d_val <- study_of(20260803)
truth_val <- mean(tau_of(d_val$w))
truth_tgt <- mean(tau_of(w_target))
print(round(c(plots = n_study, fenced = sum(d_val$a), parcels = n_target,
valley_mean = mean(d_val$w), valley_sd = sd(d_val$w),
valley_min = min(d_val$w), valley_max = max(d_val$w)), 3)) plots fenced parcels valley_mean valley_sd valley_min
320.000 137.000 40000.000 18.468 7.115 4.980
valley_max
41.145 print(round(c(district_mean = mean(w_target), district_sd = sd(w_target),
district_min = min(w_target), district_max = max(w_target),
propensity_min = min(prop_of(d_val$w)),
propensity_max = max(prop_of(d_val$w))), 3)) district_mean district_sd district_min district_max propensity_min
26.091 5.450 8.817 43.155 0.191
propensity_max
0.810 print(round(c(zero_effect_moisture = zero_at, truth_valley = truth_val,
truth_district = truth_tgt), 4))zero_effect_moisture truth_valley truth_district
26.0000 0.7532 -0.0091 The valley averages 18.47 per cent moisture with a standard deviation of 7.11; the district averages 26.09 per cent with a standard deviation of 5.45. Both sit inside the same physical range, so this is not yet a story about the district holding conditions the valley has never seen; that case comes later. It is a story about the same conditions mixed in different proportions.
The effect function crosses zero at 26 per cent moisture. Below that, fencing adds species; above it, fencing removes them, the ungrazed sward closing over and the small forbs going out. Average \(\tau(W)\) over the valley’s 320 plots and the answer is 0.7532 species per quadrat. Average the same function over the district’s 40000 parcels and the answer is -0.0091. One number says fencing gains most of a species per quadrat; the other says fencing does nothing at all. Neither is wrong. They are averages of one function over two different populations.
The fencing propensity runs from 0.191 on the driest plot to 0.81 on the wettest, so both treatments are possible everywhere in the valley and internal positivity is comfortable. The failures measured later are not positivity failures inside the study; they are failures of a different condition that looks similar.
The valley answer, done properly
Standardisation on the study data is the g-computation post’s procedure without a comma changed. Fit an outcome model that allows the effect to depend on moisture, predict every plot’s richness twice, once with the fence and once without, and take the mean difference.
gcomp <- function(fit, newdat) {
mean(predict(fit, transform(newdat, a = 1)) -
predict(fit, transform(newdat, a = 0)))
}
se_marg <- function(fit, wbar, extra = 0) {
kv <- setNames(rep(0, length(coef(fit))), names(coef(fit)))
kv["a"] <- 1; kv["a:w"] <- wbar
sqrt(drop(t(kv) %*% vcov(fit) %*% kv) + extra)
}
fit_val <- lm(y ~ a * w, data = d_val)
naive_val <- mean(d_val$y[d_val$a == 1]) - mean(d_val$y[d_val$a == 0])
est_val <- gcomp(fit_val, d_val)
se_val <- se_marg(fit_val, mean(d_val$w),
unname(coef(fit_val)["a:w"])^2 * var(d_val$w) / n_study)
ci_val <- est_val + c(-1, 1) * qt(0.975, df.residual(fit_val)) * se_val
print(round(summary(fit_val)$coefficients, 4)) Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.0379 0.3183 28.3964 0
a 3.1083 0.5014 6.1993 0
w 0.2151 0.0178 12.1123 0
a:w -0.1233 0.0250 -4.9357 0print(round(c(naive_difference = naive_val, standardised_valley = est_val,
se = se_val, lower = ci_val[1], upper = ci_val[2],
valley_truth = truth_val), 4)) naive_difference standardised_valley se lower
1.4069 0.8302 0.1861 0.4641
upper valley_truth
1.1964 0.7532 The unadjusted difference between fenced and open plots is 1.4069 species, nearly twice the truth. The confounding runs through the main effect: fenced plots are the wet ones, wet grassland carries more species whatever is done to it, and the raw comparison charges that baseline difference to the fence. Standardisation returns 0.8302 species with an interval of 0.4641 to 1.1964, covering the valley truth of 0.7532.
Two details in the standard error recur below. The estimator is a linear combination of the fitted coefficients, \(\hat\beta_A + \hat\beta_{AW}\bar{W}\), so its variance follows from the coefficient covariance without a bootstrap; and because the valley’s own mean moisture is estimated from 320 plots, the valley version carries an extra term \(\hat\beta_{AW}^2 \operatorname{Var}(\bar W)\) that the district does not, the district’s mean being a census figure over 40000 parcels.
One dataset settles nothing, so the procedure runs over many valleys. Each replicate redraws the plots, the fencing and the noise, refits, and records the estimate, its interval and the valley truth for that draw.
n_rep <- 500
rep_tab <- matrix(NA_real_, n_rep, 5)
for (b in seq_len(n_rep)) {
db <- study_of(41000 + b)
fb <- lm(y ~ a * w, data = db)
rep_tab[b, ] <- c(
gcomp(fb, db),
se_marg(fb, mean(db$w), unname(coef(fb)["a:w"])^2 * var(db$w) / n_study),
gcomp(fb, data.frame(w = w_target)),
se_marg(fb, mean(w_target)), mean(tau_of(db$w)))
}
colnames(rep_tab) <- c("est_v", "se_v", "est_t", "se_t", "truth_v")
rep_tab <- as.data.frame(rep_tab)
t_crit <- qt(0.975, n_study - 4)
cov_own <- mean(abs(rep_tab$est_v - rep_tab$truth_v) < t_crit * rep_tab$se_v)
cov_wrong <- mean(abs(rep_tab$est_v - truth_tgt) < t_crit * rep_tab$se_v)
bias_val <- mean(rep_tab$est_v - rep_tab$truth_v)
mc_val <- sd(rep_tab$est_v - rep_tab$truth_v) / sqrt(n_rep)
print(round(c(replicates = n_rep, mean_estimate = mean(rep_tab$est_v),
mean_valley_truth = mean(rep_tab$truth_v),
bias = bias_val, mc_se_of_bias = mc_val,
coverage_of_valley_truth = cov_own,
coverage_of_district_truth = cov_wrong,
mean_se = mean(rep_tab$se_v),
sd_of_estimates = sd(rep_tab$est_v)), 5)) replicates mean_estimate
500.00000 0.81262
mean_valley_truth bias
0.79984 0.01278
mc_se_of_bias coverage_of_valley_truth
0.00878 0.94000
coverage_of_district_truth mean_se
0.01400 0.19511
sd_of_estimates
0.19759 Over 500 replicates the standardised estimate averages 0.8126 against a mean valley truth of 0.7998, a bias of 0.01278 with a Monte Carlo standard error of 0.00878, and the interval covers the valley truth 94 per cent of the time. The model form is correct, the confounder is measured, positivity holds, and the answer is unbiased with nominal coverage.
What that interval covers, and what it does not
The same numbers read against the other target tell a different story. Nothing about the estimator changes; only the question does.
gap_val <- truth_val - truth_tgt
gap_in_se <- (mean(rep_tab$est_v) - truth_tgt) / mean(rep_tab$se_v)
nearest_edge <- min(abs(rep_tab$est_v - t_crit * rep_tab$se_v - truth_tgt))
print(round(c(valley_truth = truth_val, district_truth = truth_tgt,
gap = gap_val, gap_in_standard_errors = gap_in_se,
coverage_of_district_truth = cov_wrong,
replicates_covering = round(cov_wrong * n_rep),
closest_approach_of_lower_limit = nearest_edge), 4)) valley_truth district_truth
0.7532 -0.0091
gap gap_in_standard_errors
0.7623 4.2117
coverage_of_district_truth replicates_covering
0.0140 7.0000
closest_approach_of_lower_limit
0.0021 The valley interval contains the district truth in 1.4 per cent of replicates, 7 out of 500. The estimate sits 4.212 standard errors from the quantity the agency will act on, and the smallest gap in either direction between that quantity and a replicate’s lower limit was 0.0021 species.
This is a confident, well calibrated interval around the wrong estimand. Every diagnostic an analyst would reach for is clean, because nothing is wrong with the model: residuals fine, coverage nominal for the quantity estimated, confounding handled, a negative control null. The failure is in the sentence written under the fit, which silently swaps the population the average was taken over for the population the reader has in mind. No statistic sees that swap, because it happens in the prose.
The transport formula
The repair is one substitution. Standardisation averages the predicted contrast over a covariate distribution; supply the district’s distribution instead of the valley’s. Written out, with \(P\) the study population and \(Q\) the target,
\[\hat\tau_Q \;=\; \frac{1}{n_Q}\sum_{i \in Q}\Big[\hat{m}(1, W_i) - \hat{m}(0, W_i)\Big]\]
where \(\hat{m}\) is the outcome model fitted on the study data alone. It never sees a district outcome, because there are none; it sees 40000 district moisture values and predicts into them. This is the outcome-model route, and it is the object Pearl and Bareinboim (2014) obtain from their selection diagrams when the difference between populations is confined to the distribution of a measured modifier.
A second route never predicts anything. Reweight the study plots so that their moisture distribution matches the district’s, then estimate the effect in the reweighted sample. The weight is the density ratio \(f_Q(W)/f_P(W)\), estimated by pooling the study plots with a sample of district parcels, fitting a model for which population a unit came from, and taking the odds. Cole and Stuart (2010) introduced this for taking a clinical trial to a national target population, and Westreich et al (2017) set out the same weight as inverse odds of sampling weights; Stuart et al (2011) use the participation model as a diagnostic for how far apart two populations are before any estimate is attempted.
est_out <- gcomp(fit_val, data.frame(w = w_target))
se_out <- se_marg(fit_val, mean(w_target))
ci_out <- est_out + c(-1, 1) * qt(0.975, df.residual(fit_val)) * se_out
n_pool <- 2000
set.seed(77001)
w_pool <- sample(w_target, n_pool)
ratio_of <- function(dstudy, wpop) {
pooled <- data.frame(w = c(dstudy$w, wpop),
s = c(rep(1L, nrow(dstudy)), rep(0L, length(wpop))))
m_part <- suppressWarnings(glm(s ~ w + I(w^2), binomial, pooled))
p_in <- predict(m_part, data.frame(w = dstudy$w), type = "response")
((1 - p_in) / p_in) * (nrow(dstudy) / length(wpop))
}
hajek <- function(dstudy, e_hat, wt) {
a1 <- dstudy$a == 1; a0 <- !a1
sum(wt[a1] * dstudy$y[a1] / e_hat[a1]) / sum(wt[a1] / e_hat[a1]) -
sum(wt[a0] * dstudy$y[a0] / (1 - e_hat[a0])) /
sum(wt[a0] / (1 - e_hat[a0]))
}
dens_ratio <- ratio_of(d_val, w_pool)
e_val <- fitted(glm(a ~ w, family = binomial, data = d_val))
est_wt <- hajek(d_val, e_val, dens_ratio)
set.seed(77002)
n_boot <- 600
boot_wt <- numeric(n_boot)
for (b in seq_len(n_boot)) {
db <- d_val[sample.int(n_study, n_study, replace = TRUE), ]
eb <- fitted(suppressWarnings(glm(a ~ w, binomial, db)))
boot_wt[b] <- hajek(db, eb, ratio_of(db, w_pool))
}
ci_wt <- unname(quantile(boot_wt, c(0.025, 0.975)))
print(round(c(valley_standardised = est_val,
transported_outcome_model = est_out,
transported_weighting = est_wt, district_truth = truth_tgt), 4)) valley_standardised transported_outcome_model transported_weighting
0.8302 -0.1101 -0.1309
district_truth
-0.0091 print(round(c(se_outcome = se_out, lower_outcome = ci_out[1],
upper_outcome = ci_out[2], se_weighting = sd(boot_wt),
lower_weighting = ci_wt[1], upper_weighting = ci_wt[2]), 4)) se_outcome lower_outcome upper_outcome se_weighting lower_weighting
0.2568 -0.6154 0.3953 0.2896 -0.6701
upper_weighting
0.4366 On this valley the outcome-model transport returns -0.1101 species with an interval of -0.6154 to 0.3953, and the weighting route returns -0.1309 with a bootstrap interval of -0.6701 to 0.4366. Both cover the district truth of -0.0091, and both sit a long way from the 0.8302 the same data gave for the valley. Two estimators, one built on the outcome model and one on a pair of propensity-style models, land in the same place.
The weighting route is wider here, its standard error 0.2896 against 0.2568, which is the usual ordering: it throws away the parametric structure the outcome model exploits. Dahabreh et al (2020) set out both kinds together with the augmented combination that the AIPW tutorial builds in the internal-validity setting; whether it buys the same insurance here from the same two models is measured further down, and the answer depends on the target. Whether transport is unbiased, rather than lucky on one valley, is the question already asked of standardisation, and the replicate loop stored the answer.
cov_trans <- mean(abs(rep_tab$est_t - truth_tgt) < t_crit * rep_tab$se_t)
bias_trans <- mean(rep_tab$est_t) - truth_tgt
mc_trans <- sd(rep_tab$est_t) / sqrt(n_rep)
se_price <- mean(rep_tab$se_t) / mean(rep_tab$se_v)
print(c(mean_transported = formatC(mean(rep_tab$est_t), format = "f", digits = 5),
district_truth = formatC(truth_tgt, format = "f", digits = 5),
bias = formatC(bias_trans, format = "f", digits = 5),
mc_se_of_mean = formatC(mc_trans, format = "f", digits = 5)))mean_transported district_truth bias mc_se_of_mean
"-0.00005" "-0.00913" "0.00908" "0.01285" print(round(c(coverage_of_district_truth = cov_trans,
mean_se_valley = mean(rep_tab$se_v),
mean_se_transported = mean(rep_tab$se_t), se_ratio = se_price,
sd_ratio = sd(rep_tab$est_t) / sd(rep_tab$est_v)), 4))coverage_of_district_truth mean_se_valley
0.9320 0.1951
mean_se_transported se_ratio
0.2755 1.4120
sd_ratio
1.4545 Over the same 500 replicates the transported estimate averages -0.00005 against a district truth of -0.00913, a bias of 0.00908 with a Monte Carlo standard error of 0.01285, which is to say no measurable bias at this number of replicates. Coverage of the district truth is 93.2 per cent, against the 1.4 per cent the untransported interval managed.
The price
Transport is not free, and the bill arrives in the standard error.
lev_gap <- (mean(w_target) - mean(d_val$w)) / sd(d_val$w)
width_gain <- 100 * (se_price - 1)
n_equiv <- n_study / se_price^2
print(round(c(valley_mean_moisture = mean(d_val$w),
district_mean_moisture = mean(w_target),
gap_in_valley_sd_units = lev_gap, se_ratio = se_price,
per_cent_wider = width_gain, plots_equivalent_to = n_equiv,
plots_actual = n_study,
plots_thrown_away = n_study - n_equiv), 3)) valley_mean_moisture district_mean_moisture gap_in_valley_sd_units
18.468 26.091 1.071
se_ratio per_cent_wider plots_equivalent_to
1.412 41.198 160.506
plots_actual plots_thrown_away
320.000 159.494 The transported standard error is 1.412 times the valley one, so the interval is 41.2 per cent wider. The reason is geometric. The estimator is \(\hat\beta_A + \hat\beta_{AW}\bar{W}\), and the variance of that combination is smallest when \(\bar{W}\) sits at the centre of the study’s own moisture distribution and grows as it moves away. The district’s mean sits 1.071 valley standard deviations from the valley’s mean, and the study data is thin out there: the plots that carry information about the effect on wet ground are the handful of plots that happened to be wet.
Read in units the field crew would recognise, an interval 1.412 times wider is what 161 plots would have given at the valley’s own target. The trial has 320. Transporting the answer costs the equivalent of 159 plots, spent not on measuring anything but on aiming the estimate somewhere the design was not centred.
A diagnostic that fires before the interval does
The density ratio is worth looking at for its own sake, because it is the only part of the procedure that speaks about the two populations rather than about one model. A weighted average with unequal weights behaves like a smaller sample, and the standard measure of how much smaller is the effective sample size,
\[n_{\text{eff}} = \frac{\left(\sum_i w_i\right)^2}{\sum_i w_i^2}\]
which equals the sample size when all weights are equal and falls towards one as a single plot takes over.
ess_of <- function(wt) sum(wt)^2 / sum(wt^2)
# a second, wetter candidate target, used from here on as a contrast
set.seed(20260806)
w_wet <- w_lo + (w_hi - w_lo) * rbeta(n_target, 12, 4)
set.seed(77003)
w_pool_wet <- sample(w_wet, n_pool)
ess_main <- ess_of(dens_ratio)
share_main <- max(dens_ratio) / sum(dens_ratio)
top10_main <- sum(sort(dens_ratio, decreasing = TRUE)[1:10]) / sum(dens_ratio)
print(round(c(plots = n_study, effective_sample_size = ess_main,
per_cent_of_plots = 100 * ess_main / n_study,
largest_weight = max(dens_ratio),
largest_weight_share_pct = 100 * share_main,
top_ten_share_pct = 100 * top10_main,
se_ratio_measured = se_price), 3)) plots effective_sample_size per_cent_of_plots
320.000 143.950 44.985
largest_weight largest_weight_share_pct top_ten_share_pct
3.537 1.116 11.145
se_ratio_measured
1.412 The transport weights have an effective sample size of 144 plots out of 320, 45 per cent of the sample. The largest single weight takes 1.12 per cent of the total and the ten largest take 11.14 per cent, so no one plot is running the show. That agrees with what the standard error did: an effective sample size of 45 per cent and a standard error ratio of 1.412 are two views of one loss of information. The diagnostic can be computed before any outcome is used, from the study covariates and the target covariates alone, so a team that runs it at the design stage learns how much of its sample will survive the journey and can still move plots.
An effect modifier nobody measured
Everything so far worked because soil moisture was the whole story: the only thing that modified the effect, and measured in both places. Drop either half of that and the transport formula returns a number with the same shape and none of the meaning.
Suppose a second modifier exists. Soil phosphorus varies within the valley, and the fencing effect falls as phosphorus rises, the extra fertility feeding the tall grasses that shade the forbs out once grazing stops. Phosphorus is not a confounder here: it has nothing to do with which plots got fenced, so leaving it out of the model costs the valley estimate nothing. But the district lies on a different parent material and is markedly richer in phosphorus, and nobody has mapped it.
p_shift <- 1.2
p_slope <- -0.55
study_p <- function(seed) {
set.seed(seed)
w <- w_lo + (w_hi - w_lo) * rbeta(n_study, 2.2, 4.4)
phos <- rnorm(n_study, 0, 1)
a <- rbinom(n_study, 1, prop_of(w))
y <- mu_of(w) + 0.5 * phos + (tau_of(w) + p_slope * phos) * a +
rnorm(n_study, 0, sd_noise)
data.frame(y = y, a = a, w = w, phos = phos)
}
truth_tgt_p <- truth_tgt + p_slope * p_shift
n_rep_p <- 400
rep_p <- matrix(NA_real_, n_rep_p, 4)
for (b in seq_len(n_rep_p)) {
db <- study_p(52000 + b)
fb <- lm(y ~ a * w, data = db)
rep_p[b, ] <- c(gcomp(fb, data.frame(w = w_target)),
se_marg(fb, mean(w_target)), gcomp(fb, db),
mean(tau_of(db$w) + p_slope * db$phos))
}
colnames(rep_p) <- c("est_t", "se_t", "est_v", "truth_v")
rep_p <- as.data.frame(rep_p)
cov_p <- mean(abs(rep_p$est_t - truth_tgt_p) < t_crit * rep_p$se_t)
bias_p <- mean(rep_p$est_t) - truth_tgt_p
half_p <- t_crit * mean(rep_p$se_t)
correction <- mean(rep_p$est_v) - mean(rep_p$est_t)
print(round(c(district_truth_with_phosphorus = truth_tgt_p,
transported_estimate = mean(rep_p$est_t), bias = bias_p,
mc_se = sd(rep_p$est_t) / sqrt(n_rep_p),
coverage_of_district_truth = cov_p,
mean_interval_half_width = half_p,
bias_over_half_width = bias_p / half_p), 4))district_truth_with_phosphorus transported_estimate
-0.6691 0.0130
bias mc_se
0.6822 0.0137
coverage_of_district_truth mean_interval_half_width
0.3150 0.5592
bias_over_half_width
1.2198 print(round(c(valley_estimate = mean(rep_p$est_v),
valley_truth = mean(rep_p$truth_v),
correction_transport_made = correction,
error_it_left_behind = bias_p,
error_as_share_of_correction = bias_p / correction), 4)) valley_estimate valley_truth
0.8147 0.8001
correction_transport_made error_it_left_behind
0.8017 0.6822
error_as_share_of_correction
0.8509 The district runs 1.2 standard deviations higher in phosphorus than the valley, and the fencing effect falls by 0.55 species per standard deviation of it. The transported estimate averages 0.013 species over 400 replicates, against a district truth that is now -0.6691. The bias is 0.6822, the interval covers the truth 31.5 per cent of the time, and the interval itself has not changed shape at all: its half width is 0.5592, essentially what it was in the working case. The valley estimate is still fine, averaging 0.8147 against a valley truth of 0.8001, because phosphorus is balanced across fencing within the valley and its omission costs nothing there.
Put the two numbers side by side. Averaged over the same 400 replicates, transporting moved the estimate 0.8017 species, from the valley average to the district average, and the error left behind by the unmeasured modifier is 0.6822 species, 85.1 per cent of the correction just made. The procedure travelled most of the distance, stopped, and gave no indication of where. Catford et al (2022) make the general version of this point: context dependence is not noise to be averaged away but a set of specific moderators, and a claim that travels is a claim about which moderators were captured.
There is no test for this. The valley data hold no information about a variable never measured, and the district data hold no outcomes at all. The only defence is subject knowledge: list what plausibly changes the sign or size of the effect, and check that each item is measured in both populations. A modifier that fails the check does not make transport impossible; it makes the transported number conditional on that modifier being similarly distributed, a claim to be argued rather than computed.
Where the target leaves the study behind
The second failure is different in kind. So far both populations occupied the same range of soil moisture and only their proportions differed. Now the target is a wetter district, averaging 35.5 per cent moisture, with a substantial share of parcels above anything the valley contains.
hinge_at <- 30
sup_hi <- unname(quantile(d_val$w, 0.99))
frac_in_main <- mean(w_target <= sup_hi)
frac_in_wet <- mean(w_wet <= sup_hi)
frac_above_max <- mean(w_wet > max(d_val$w))
dens_wet <- ratio_of(d_val, w_pool_wet)
ess_wet <- ess_of(dens_wet)
share_wet <- max(dens_wet) / sum(dens_wet)
print(round(c(valley_max = max(d_val$w), valley_99th = sup_hi,
wet_target_mean = mean(w_wet), wet_target_max = max(w_wet),
share_of_wet_target_inside = frac_in_wet,
share_of_wet_above_valley_max = frac_above_max,
share_of_district_inside = frac_in_main), 4)) valley_max valley_99th
41.1445 37.0883
wet_target_mean wet_target_max
35.4947 45.4707
share_of_wet_target_inside share_of_wet_above_valley_max
0.6031 0.0846
share_of_district_inside
0.9828 print(round(c(ess_district = ess_main, ess_wet_target = ess_wet,
ess_wet_as_pct_of_plots = 100 * ess_wet / n_study,
largest_weight_district = max(dens_ratio),
largest_weight_wet = max(dens_wet),
largest_share_wet_pct = 100 * share_wet), 3)) ess_district ess_wet_target ess_wet_as_pct_of_plots
143.950 9.377 2.930
largest_weight_district largest_weight_wet largest_share_wet_pct
3.537 91.172 23.419 The diagnostics fire loudly. Only 60.3 per cent of the wetter target lies below the valley’s ninety-ninth percentile, against 98.3 per cent for the district used earlier, and 8.5 per cent of it lies above the wettest plot the trial ever measured. The effective sample size of the transport weights collapses from 144 plots to 9.4, or 2.9 per cent of the sample, and one plot now carries 23.4 per cent of the total weight against 1.12 per cent before. A weighting estimator built on that is effectively an experiment with 9 plots in it.
The outcome-model route does not collapse in the same visible way, and that is what makes it dangerous. It will happily predict at 45.5 per cent moisture, the wettest parcel in the target, because a fitted line has a value everywhere. Whether that value is right depends on whether the effect function really is the straight line the model assumed, out where no plot ever was, and the study data cannot settle it. Two truths, identical over the valley’s dense range and different in its sparse wet tail and beyond it: the straight line used so far, and one where the fencing penalty accelerates once the ground is waterlogged past 30 per cent moisture.
tau_hinge <- function(w) tau_of(w) - 0.02 * pmax(w - hinge_at, 0)^2
study_h <- function(seed, tfun) {
set.seed(seed)
w <- w_lo + (w_hi - w_lo) * rbeta(n_study, 2.2, 4.4)
a <- rbinom(n_study, 1, prop_of(w))
data.frame(y = mu_of(w) + tfun(w) * a + rnorm(n_study, 0, sd_noise),
a = a, w = w)
}
n_rep_o <- 400
rep_o <- matrix(NA_real_, n_rep_o, 9)
for (b in seq_len(n_rep_o)) {
d_lin <- study_h(63000 + b, tau_of)
d_hin <- study_h(63000 + b, tau_hinge)
f_lin <- lm(y ~ a * w, data = d_lin)
f_hin <- lm(y ~ a * w, data = d_hin)
keep <- w_wet <= unname(quantile(d_hin$w, 0.99))
rep_o[b, ] <- c(gcomp(f_lin, data.frame(w = w_wet)),
se_marg(f_lin, mean(w_wet)),
gcomp(f_hin, data.frame(w = w_wet)),
se_marg(f_hin, mean(w_wet)),
gcomp(f_hin, data.frame(w = w_wet[keep])),
se_marg(f_hin, mean(w_wet[keep])),
mean(tau_hinge(w_wet[keep])), mean(keep),
gcomp(f_hin, d_hin) - mean(tau_hinge(d_hin$w)))
}
colnames(rep_o) <- c("lin", "lin_se", "hin", "hin_se", "res", "res_se",
"res_truth", "covered", "in_sample_bias")
rep_o <- as.data.frame(rep_o)
truth_wet_lin <- mean(tau_of(w_wet))
truth_wet_hin <- mean(tau_hinge(w_wet))
cov_lin <- mean(abs(rep_o$lin - truth_wet_lin) < t_crit * rep_o$lin_se)
cov_hin <- mean(abs(rep_o$hin - truth_wet_hin) < t_crit * rep_o$hin_se)
bias_hin <- mean(rep_o$hin) - truth_wet_hin
print(round(c(wet_truth_straight_line = truth_wet_lin,
estimate_straight_line = mean(rep_o$lin),
bias_straight_line = mean(rep_o$lin) - truth_wet_lin,
coverage_straight_line = cov_lin), 4))wet_truth_straight_line estimate_straight_line bias_straight_line
-0.9495 -0.9438 0.0057
coverage_straight_line
0.9550 print(round(c(wet_truth_accelerating = truth_wet_hin,
estimate_accelerating = mean(rep_o$hin),
bias_accelerating = bias_hin,
coverage_accelerating = cov_hin,
mean_se_straight = mean(rep_o$lin_se),
mean_se_accelerating = mean(rep_o$hin_se),
in_sample_bias = mean(rep_o$in_sample_bias)), 4))wet_truth_accelerating estimate_accelerating bias_accelerating
-1.9111 -1.1588 0.7523
coverage_accelerating mean_se_straight mean_se_accelerating
0.6525 0.4843 0.4855
in_sample_bias
0.0065 Under the straight-line truth, transport to the wetter target is unbiased: -0.9438 against -0.9495, coverage 95.5 per cent, even though 8.5 per cent of the target is outside the study’s observed range. Extrapolation into a region the data never saw is free when the assumed functional form happens to be correct there.
Under the accelerating truth it is not free. The estimate averages -1.1588 against a truth of -1.9111, a bias of 0.7523 species, and coverage drops to 65.25 per cent. The two cases have standard errors of 0.4843 and 0.4855, effectively the same, and identical overlap diagnostics, those being functions of the covariates only. Within the valley the curvature is invisible: the standardised valley estimate under it is biased by 0.0065 species, only a handful of plots being wet enough to show it.
The honest response is to stop extrapolating and say so: restrict the target to the region of common support, report the effect there, and report what fraction of the target it covers.
cov_res <- mean(abs(rep_o$res - rep_o$res_truth) < t_crit * rep_o$res_se)
bias_res <- mean(rep_o$res - rep_o$res_truth)
print(round(c(share_of_wet_target_covered = mean(rep_o$covered),
sd_of_that_share = sd(rep_o$covered),
share_using_this_valley = frac_in_wet,
restricted_estimate = mean(rep_o$res),
restricted_truth = mean(rep_o$res_truth),
bias = bias_res, coverage = cov_res), 4))share_of_wet_target_covered sd_of_that_share
0.4635 0.1036
share_using_this_valley restricted_estimate
0.6031 -0.7290
restricted_truth bias
-0.7546 0.0257
coverage
0.9525 print(round(c(whole_target_truth = truth_wet_hin,
whole_target_estimate = mean(rep_o$hin),
restricted_truth = mean(rep_o$res_truth),
understatement = mean(rep_o$res_truth) - truth_wet_hin), 4)) whole_target_truth whole_target_estimate restricted_truth
-1.9111 -1.1588 -0.7546
understatement
1.1565 Restricted to the parcels below the valley’s ninety-ninth percentile of moisture, which is 46.4 per cent of the wetter target, the transported estimate is -0.729 against a restricted truth of -0.7546, a bias of 0.0257 with 95.25 per cent coverage. The estimator works again, because the question shrank to fit the data.
How far it shrank is itself uncertain. The covered share averages 46.4 per cent across replicate valleys with a standard deviation of 10.4 percentage points, against 60.3 per cent for the valley used above. The common-support boundary is a ninety-ninth percentile of 320 plots, so the size of the region a study can honestly speak for is itself a statistic with a standard error.
That shrinkage is the cost, and it should be reported rather than buried. The restricted answer, -0.7546 species, understates the whole-target effect of -1.9111 by 1.1565 species, because the parcels left out are exactly the ones where fencing does most harm. A regulator reading only the restricted number would underestimate the damage; a regulator reading only the extrapolated number would get a figure 39.4 per cent off with an interval that does not admit it. The choice between them is not statistical. It is a choice about which kind of error the decision can tolerate, and it can only be made if both numbers and the coverage fraction are on the page.
The augmented transport estimator
Each transport route so far leans on one kind of model. The outcome-model route needs the line to be right wherever the target has parcels; the weighting route needs the participation model and the fencing propensity to be right. The AIPW tutorial combines the two kinds of model for the internal-validity question, and Dahabreh et al (2020) give the transport version. Their normalised form, written with this post’s pieces, is
\[\hat\tau_{\text{aug}} = \frac{1}{n_Q}\sum_{j \in Q}\Big[\hat{m}(1, W_j) - \hat{m}(0, W_j)\Big] + \frac{\sum_i v_{1i}\{Y_i - \hat{m}(1, W_i)\}}{\sum_i v_{1i}} - \frac{\sum_i v_{0i}\{Y_i - \hat{m}(0, W_i)\}}{\sum_i v_{0i}}\]
where the sums over \(i\) run over the study plots, \(v_{1i} = \hat{r}(W_i)A_i/\hat{e}(W_i)\) and \(v_{0i} = \hat{r}(W_i)(1 - A_i)/\{1 - \hat{e}(W_i)\}\), \(\hat{r}\) is the density ratio from ratio_of() and \(\hat{e}\) the fitted fencing propensity. The first term is the outcome-model transport. The other two are the weighting route applied to the outcome model’s residuals, so they add back what the model missed, in the proportions the target would have. In Dahabreh et al the first term averages over a sample from the target; here it runs over the census of parcels, and the target sample is used only to fit the participation model.
The outcome model is y ~ a * w throughout: right under the straight-line truth, wrong under the accelerating one. The propensity model is the true one. The participation model s ~ w + I(w^2) is only an approximation, because the log ratio of two beta densities is linear in the logarithm of rescaled moisture and of its complement, not in moisture and its square. So the weighting route is not guaranteed to be unbiased either, and it is measured here rather than assumed. The chunk runs all three estimators on both targets under both truths, on fresh seeds, and each replicate draws its own sample of 2000 parcels from the target to fit the participation model, as a new study would. Two variants ride along as checks: the single parcel sample of the weight diagnostics reused in every replicate, and the exact ratio of the two beta densities in place of a fitted one. For this outcome model the standardised contrast is exactly \(\hat\beta_A + \hat\beta_{AW}\bar{W}\), which the loop uses in place of gcomp() to save time; the chunk checks that the two agree on the valley data. The augmented estimator’s interval uses a plug-in standard error, the standard deviation of the per-plot correction terms over \(\sqrt{n}\). That is the same kind of influence-function estimator the AIPW tutorial uses, and it treats all three fitted models as known. The weighting route gets no per-replicate interval here.
eff_mean <- function(fit, w) unname(coef(fit)["a"] + coef(fit)["a:w"] * mean(w))
aug_corr <- function(dstudy, fit, e_hat, wt) {
a1 <- dstudy$a == 1
r1 <- dstudy$y - predict(fit, transform(dstudy, a = 1))
r0 <- dstudy$y - predict(fit, transform(dstudy, a = 0))
v1 <- wt * a1 / e_hat
v0 <- wt * (!a1) / (1 - e_hat)
c1 <- sum(v1 * r1) / sum(v1)
c0 <- sum(v0 * r0) / sum(v0)
phi <- v1 / mean(v1) * (r1 - c1) - v0 / mean(v0) * (r0 - c0)
c(c1 - c0, sd(phi) / sqrt(nrow(dstudy)))
}
arm_means <- function(dstudy, e_hat, wt) {
a1 <- dstudy$a == 1
c(sum(wt[a1] * dstudy$w[a1] / e_hat[a1]) / sum(wt[a1] / e_hat[a1]),
sum(wt[!a1] * dstudy$w[!a1] / (1 - e_hat[!a1])) / sum(wt[!a1] / (1 - e_hat[!a1])))
}
gcomp_gap <- gcomp(fit_val, data.frame(w = w_wet)) - eff_mean(fit_val, w_wet)
case_a <- c("district, straight line", "district, accelerating",
"wetter, straight line", "wetter, accelerating")
truth_a <- c(mean(tau_of(w_target)), mean(tau_hinge(w_target)),
mean(tau_of(w_wet)), mean(tau_hinge(w_wet)))
n_rep_a <- 1000
ratio_exact <- function(w, a, b) {
x <- (w - w_lo) / (w_hi - w_lo)
dbeta(x, a, b) / dbeta(x, 2.2, 4.4)
}
rep_a <- array(NA_real_, c(n_rep_a, 4, 10), dimnames = list(
NULL, case_a, c("om", "om_se", "wt", "aug", "aug_se", "beyond",
"wt_fixed", "aug_fixed", "wt_exact", "aug_exact")))
arm_w <- matrix(NA_real_, n_rep_a, 6, dimnames = list(NULL, c(
"fenced_d", "open_d", "fenced_wet", "open_wet", "fenced_wet_exact", "open_wet_exact")))
for (b in seq_len(n_rep_a)) {
d_l <- study_h(64000 + b, tau_of)
d_h <- study_h(64000 + b, tau_hinge)
e_b <- fitted(glm(a ~ w, family = binomial, data = d_l))
set.seed(65000 + b)
pool_b <- list(sample(w_target, n_pool), sample(w_wet, n_pool))
wt_b <- list(ratio_of(d_l, pool_b[[1]]), ratio_of(d_l, pool_b[[2]]))
wt_f <- list(ratio_of(d_l, w_pool), ratio_of(d_l, w_pool_wet))
wt_x <- list(ratio_exact(d_l$w, 7.235, 6.515), ratio_exact(d_l$w, 12, 4))
dat_b <- list(d_l, d_h)
fit_b <- list(lm(y ~ a * w, data = d_l), lm(y ~ a * w, data = d_h))
arm_w[b, ] <- c(arm_means(d_l, e_b, wt_b[[1]]), arm_means(d_l, e_b, wt_b[[2]]),
arm_means(d_l, e_b, wt_x[[2]]))
for (k in 1:4) {
j <- 2 - k %% 2; g <- 1 + (k > 2)
w_to <- if (g == 1) w_target else w_wet
tfun <- if (j == 1) tau_of else tau_hinge
cf <- coef(fit_b[[j]])
past <- w_to > max(d_l$w)
om <- eff_mean(fit_b[[j]], w_to)
ac <- aug_corr(dat_b[[j]], fit_b[[j]], e_b, wt_b[[g]])
rep_a[b, k, ] <- c(om, se_marg(fit_b[[j]], mean(w_to)),
hajek(dat_b[[j]], e_b, wt_b[[g]]), om + ac[1], ac[2],
sum(cf["a"] + cf["a:w"] * w_to[past] - tfun(w_to[past])) /
length(w_to),
hajek(dat_b[[j]], e_b, wt_f[[g]]),
om + aug_corr(dat_b[[j]], fit_b[[j]], e_b, wt_f[[g]])[1],
hajek(dat_b[[j]], e_b, wt_x[[g]]),
om + aug_corr(dat_b[[j]], fit_b[[j]], e_b, wt_x[[g]])[1])
}
}
err_a <- sweep(rep_a[, , c("om", "wt", "aug")], 2, truth_a)
mc_of <- function(x) sd(x) / sqrt(length(x))
aug_tab <- data.frame(
truth = truth_a,
bias_om = colMeans(err_a[, , "om"]), mc_om = apply(err_a[, , "om"], 2, mc_of),
bias_wt = colMeans(err_a[, , "wt"]), mc_wt = apply(err_a[, , "wt"], 2, mc_of),
bias_aug = colMeans(err_a[, , "aug"]), mc_aug = apply(err_a[, , "aug"], 2, mc_of),
sd_om = apply(err_a[, , "om"], 2, sd), sd_wt = apply(err_a[, , "wt"], 2, sd),
sd_aug = apply(err_a[, , "aug"], 2, sd),
rmse_om = sqrt(colMeans(err_a[, , "om"]^2)),
rmse_wt = sqrt(colMeans(err_a[, , "wt"]^2)),
rmse_aug = sqrt(colMeans(err_a[, , "aug"]^2)),
cov_om = colMeans(abs(err_a[, , "om"]) < t_crit * rep_a[, , "om_se"]),
cov_aug = colMeans(abs(err_a[, , "aug"]) < t_crit * rep_a[, , "aug_se"]),
se_aug = colMeans(rep_a[, , "aug_se"]),
om_error_beyond_valley = colMeans(rep_a[, , "beyond"]),
row.names = case_a)
dmse_a <- colMeans(err_a[, , "aug"]^2 - err_a[, , "om"]^2)
mc_dmse_a <- apply(err_a[, , "aug"]^2 - err_a[, , "om"]^2, 2, mc_of)
mc_cov_a <- sqrt(0.95 * 0.05 / n_rep_a)
slope_mu <- mu_of(1) - mu_of(0)
wt_pred <- c(tau_1 * (mean(arm_w[, "fenced_d"]) - mean(w_target)) +
slope_mu * (mean(arm_w[, "fenced_d"]) - mean(arm_w[, "open_d"])),
tau_1 * (mean(arm_w[, "fenced_wet"]) - mean(w_wet)) +
slope_mu * (mean(arm_w[, "fenced_wet"]) - mean(arm_w[, "open_wet"])))
err_x <- sweep(rep_a[, , c("wt_fixed", "aug_fixed", "wt_exact", "aug_exact")], 2, truth_a)
alt_tab <- data.frame(
bias_wt_fixed = colMeans(err_x[, , "wt_fixed"]),
mc_wt_fixed = apply(err_x[, , "wt_fixed"], 2, mc_of),
bias_aug_fixed = colMeans(err_x[, , "aug_fixed"]),
mc_aug_fixed = apply(err_x[, , "aug_fixed"], 2, mc_of),
bias_wt_exact = colMeans(err_x[, , "wt_exact"]),
mc_wt_exact = apply(err_x[, , "wt_exact"], 2, mc_of),
bias_aug_exact = colMeans(err_x[, , "aug_exact"]),
mc_aug_exact = apply(err_x[, , "aug_exact"], 2, mc_of),
fixed_minus_fresh_aug = colMeans(rep_a[, , "aug_fixed"] - rep_a[, , "aug"]),
mc_paired_aug = apply(rep_a[, , "aug_fixed"] - rep_a[, , "aug"], 2, mc_of),
row.names = case_a)
mc_cov_obs <- sqrt(aug_tab$cov_aug * (1 - aug_tab$cov_aug) / n_rep_a)
se_to_sd <- apply(rep_a[, , "aug_se"], 2, median) / aug_tab$sd_aug
print(c(replicates = n_rep_a, gcomp_agrees_to_1e_10 = as.integer(abs(gcomp_gap) < 1e-10))) replicates gcomp_agrees_to_1e_10
1000 1 print(round(aug_tab[, 1:7], 4)) truth bias_om mc_om bias_wt mc_wt bias_aug mc_aug
district, straight line -0.0091 -0.0134 0.0085 0.0103 0.0104 -0.0072 0.0093
district, accelerating -0.0821 -0.0489 0.0086 0.0227 0.0103 0.0041 0.0093
wetter, straight line -0.9495 -0.0175 0.0150 0.2924 0.0352 0.0059 0.0309
wetter, accelerating -1.9111 0.7225 0.0154 0.5051 0.0349 0.2148 0.0328print(round(aug_tab[, 8:17], 4)) sd_om sd_wt sd_aug rmse_om rmse_wt rmse_aug cov_om
district, straight line 0.2687 0.3278 0.2926 0.2689 0.3278 0.2925 0.956
district, accelerating 0.2735 0.3253 0.2939 0.2777 0.3260 0.2938 0.943
wetter, straight line 0.4753 1.1116 0.9759 0.4754 1.1489 0.9754 0.956
wetter, accelerating 0.4868 1.1028 1.0369 0.8711 1.2125 1.0584 0.676
cov_aug se_aug om_error_beyond_valley
district, straight line 0.945 0.2919 -0.0002
district, accelerating 0.942 0.2925 0.0102
wetter, straight line 0.791 0.6530 -0.0037
wetter, accelerating 0.784 0.7013 0.4575print(round(rbind(mse_aug_minus_mse_om = dmse_a, mc_se = mc_dmse_a), 4)) district, straight line district, accelerating
mse_aug_minus_mse_om 0.0133 0.0092
mc_se 0.0020 0.0023
wetter, straight line wetter, accelerating
mse_aug_minus_mse_om 0.7254 0.3615
mc_se 0.0441 0.0533print(round(c(district_mean = mean(w_target), fenced_reweighted_d = mean(arm_w[, 1]),
open_reweighted_d = mean(arm_w[, 2]), wet_mean = mean(w_wet),
fenced_reweighted_wet = mean(arm_w[, 3]),
open_reweighted_wet = mean(arm_w[, 4]),
weighting_bias_implied_d = wt_pred[1],
weighting_bias_implied_wet = wt_pred[2]), 4)) district_mean fenced_reweighted_d
26.0913 26.0254
open_reweighted_d wet_mean
25.9743 35.4947
fenced_reweighted_wet open_reweighted_wet
34.7034 33.7694
weighting_bias_implied_d weighting_bias_implied_wet
0.0178 0.2846 print(round(alt_tab, 4)) bias_wt_fixed mc_wt_fixed bias_aug_fixed mc_aug_fixed
district, straight line 0.0303 0.0103 -0.0076 0.0092
district, accelerating 0.0370 0.0102 -0.0043 0.0092
wetter, straight line 0.2999 0.0383 0.0057 0.0332
wetter, accelerating 0.4040 0.0380 0.1130 0.0358
bias_wt_exact mc_wt_exact bias_aug_exact mc_aug_exact
district, straight line 0.0017 0.0105 -0.0089 0.0092
district, accelerating 0.0001 0.0105 -0.0097 0.0092
wetter, straight line 0.3053 0.0359 0.0056 0.0310
wetter, accelerating 0.4488 0.0369 0.1509 0.0344
fixed_minus_fresh_aug mc_paired_aug
district, straight line -0.0003 0.0005
district, accelerating -0.0084 0.0006
wetter, straight line -0.0003 0.0042
wetter, accelerating -0.1017 0.0054print(round(c(fixed_pool_mean_d = mean(w_pool), fixed_pool_mean_wet = mean(w_pool_wet),
fenced_exact_wet = mean(arm_w[, "fenced_wet_exact"]),
open_exact_wet = mean(arm_w[, "open_wet_exact"])), 4)) fixed_pool_mean_d fixed_pool_mean_wet fenced_exact_wet open_exact_wet
25.8859 35.4477 34.8793 33.8040 print(round(rbind(coverage_mc_se = mc_cov_obs, median_plugin_se_over_sd = se_to_sd), 4)) district, straight line district, accelerating
coverage_mc_se 0.0072 0.0074
median_plugin_se_over_sd 0.9848 0.9838
wetter, straight line wetter, accelerating
coverage_mc_se 0.0129 0.0130
median_plugin_se_over_sd 0.6261 0.6337The linear combination and gcomp() agree on the valley data to rounding error, which is what the 1 in the first printed line records, so the shortcut changes nothing. The district comes first. Under the straight-line truth the outcome model is right, and its bias over 1000 replicate valleys is -0.0134 (Monte Carlo standard error 0.0085), with the augmented estimator at -0.0072 (0.0093). Under the accelerating truth the outcome model is wrong, though only a little over the district’s range, and its bias is -0.0489 (0.0086). The augmentation takes that to 0.0041 (0.0093), which is no measurable bias.
The weighting route is close to centred on the district as well. Its bias is 0.0103 (0.0104) under the straight line, no measurable bias, and 0.0227 (0.0103) under the accelerating truth, a little over two Monte Carlo standard errors from zero. Under the straight-line truth its expectation has a closed form. The Hajek contrast of weighted means is \(\tau(\bar W_1) + \gamma(\bar W_1 - \bar W_0)\) plus noise with mean zero, where \(\bar W_1\) and \(\bar W_0\) are the reweighted mean moistures of the fenced and the open plots and \(\gamma\) = 0.22 is the slope of baseline richness on moisture. After weighting, the fenced plots average 26.03 per cent moisture and the open ones 25.97, against the district’s 26.09, and the formula puts the weighting bias at 0.0178, within one Monte Carlo standard error of the measured one. With the exact density ratio the bias is 0.0017 (0.0105).
The parcel sample behind the participation model is part of the estimator. The single draw of 2000 parcels used in the weight diagnostics averages 25.89 per cent moisture against the census’s 26.09, and reused in every replicate it puts the weighting bias at 0.0303 (0.0103) and 0.0370 (0.0102) under the two truths, because every replicate then inherits that one draw’s shortfall. The same swap moves the augmented estimator by -0.0003 (0.0005) species under the straight line, because its weights only multiply the outcome model’s residuals, and under the straight-line truth those average zero at every moisture, whatever the weights. Under the accelerating truth both of its models are somewhat wrong and the residuals of the wettest plots no longer average zero, but the product of two small errors is smaller still, and the swap moves it by -0.0084 (0.0006).
The protection has a premium. The augmented estimator’s standard deviation is 0.2926 under the straight line and 0.2939 under the accelerating truth, against 0.2687 and 0.2735 for the outcome model, with the weighting route wider again at 0.3278 and 0.3253. The bias removed under the accelerating truth is too small to pay for the extra spread: the augmented estimator’s mean squared error exceeds the outcome model’s by 0.0092 (Monte Carlo standard error 0.0023), and by 0.0133 (0.0020) under the straight line, where there was nothing to remove. Its plug-in interval covers the district truth 94.5 and 94.2 per cent of the time, against 95.6 and 94.3 for the outcome model’s, each with a Monte Carlo standard error of 0.7 points.
On the wetter target the reasoning in advance was that the correction has no plots wetter than the valley’s wettest, so it cannot touch the part of the target beyond it, and there it would reduce to the outcome model’s extrapolation. The measurement says otherwise. Under the accelerating truth the outcome model’s bias is 0.7225 (0.0154), and 0.4575 of it comes from parcels wetter than the wettest plot of each replicate valley. The augmented estimator’s bias is 0.2148 (0.0328), 70.3 per cent less than the outcome model’s and smaller than the out-of-range part on its own, so the correction reaches past the valley’s edge. It can only do that through the plots it has. The normalised weights spend the whole target’s mass on the study plots, the density ratio climbs steeply with moisture (the wetter curve in the weight figure above), so the wettest plots carry the share of parcels wetter than any of them, and under the accelerating truth their residuals already show the bend. They show it only as far as those plots reach, and how much of the rest the correction recovers depends on how the weights behave at the valley’s edge. With the exact density ratio the augmented bias is 0.1509 (0.0344); with the single parcel sample of the weight diagnostics reused in every replicate it is 0.1130 (0.0358), a paired shift of -0.1017 (0.0054) that comes from that one draw alone.
The price is spread. Under the accelerating truth the augmented estimator’s standard deviation is 1.0369 against the outcome model’s 0.4868, and its root mean squared error is 1.0584 against 0.8711: the mean squared error difference of 0.3615 (Monte Carlo standard error 0.0533) says the corrected estimate is further from the truth, on average, than the biased extrapolation it corrects. Under the straight line both are centred, the augmented estimator at 0.0059 (0.0309), but its standard deviation is 0.9759 against 0.4753. The plug-in interval covers 79.1 and 78.4 per cent under the two truths (Monte Carlo standard error 1.3 points), because its mean standard error, 0.6530 and 0.7013, falls well short of the actual spread: with the valley’s weights down to an effective 9.4 plots, the plug-in variance is itself estimated from those few heavily weighted plots, and under both truths the median plug-in standard error is about 0.63 of the actual standard deviation. The outcome model’s interval covers 67.6 per cent under the accelerating truth. The weighting route is biased on this target under either truth, by 0.2924 (0.0352) under the straight line and 0.5051 (0.0349) under the accelerating truth, and for the straight line the closed form above puts it at 0.2846. Most of that is imbalance between the arms rather than a shortfall against the target: after weighting the fenced plots average 34.70 per cent moisture and the open ones 33.77, and the baseline slope turns that gap into 0.2055 species of spurious effect. The participation model is not the cause: with the exact density ratio the bias is 0.3053 (0.0359) and the arms average 34.88 and 33.80. The fencing propensity has the true form as well, so in the population the two arms would reweight to the same mean; the gap is what normalised weighted means do on so few effective plots.
So the question left open earlier, whether the augmented combination buys the same insurance here from the same two models, has a split answer. On the district, which the valley’s range covers almost entirely, it behaves as the AIPW tutorial’s version does and stays centred when a model is wrong, but here the bias it guards against is small and the plain outcome model has the smaller error. Beyond the valley’s range the AIPW tutorial’s own limit applies, that double robustness “does not repair poor overlap”: the augmentation removed most of the extrapolation bias, by an amount that moved with the parcel sample behind the weights, and the heavy weights that let it do so made the estimate more variable and its interval too short. The point estimate shows neither.
The honest limit
Transportability is an assumption that a function is the same in two places, and it is not testable from the data at hand. The study valley identifies \(\tau(W)\) over the moisture values it contains; the district supplies a distribution of \(W\) and nothing else. If the conditional effect genuinely differs, because the sheep are a different breed, the soils hold water differently, or the seed bank has been depleted for longer, then no reweighting recovers the district effect and the output looks exactly as it does when transport works. Bareinboim and Pearl (2016) put this formally: which quantities transport, and by what formula, is decided by a graph encoding what is assumed to differ between populations, and the graph is an input rather than an output.
Two of the failures above are limits rather than results. The unmeasured modifier left a bias of 0.6822 species, 85.1 per cent of the correction transport had just made, with an interval of unchanged width and 31.5 per cent coverage. Avoiding it requires the modifier to be measured in both populations, which is what a mapped region and a small experiment usually do not share: the region has whatever the soil survey recorded, the trial has whatever the fieldwork recorded, and the intersection is often one or two variables.
The overlap failure is the one with a diagnostic, and even there the diagnostic is only half of what is needed. Effective sample size and the range comparison fired identically under the straight-line truth, where transport was unbiased, and under the accelerating truth, where it was off by 0.7523 species with 65.25 per cent coverage. Those diagnostics tell you where the data are thin. They cannot tell you whether the model’s shape is right where the data are thin, and that second question is the one that determines the bias.
One thing did not go the way I expected before running it. The precision cost of transport came out modest, a standard error ratio of 1.412 and an effective sample size of 45 per cent, for a shift of 1.07 valley standard deviations in the modifier mean. That is a smaller penalty than the size of the estimand change suggests: that change is 0.7623 species, several times the valley standard error, while the interval widened by 41.2 per cent. Precision loss is a poor alarm here. By the time the weights are bad enough to widen the interval noticeably, as in the wetter target where the effective sample size fell to 9.4 plots, the estimate has usually left the region where it can be trusted at all.
Where to go next
The augmented transport estimator above combines the outcome model with the participation and propensity models, the construction Doubly robust estimation with AIPW uses for the internal-validity question. On this post’s two targets its double protection was not the same bargain as there: on the district it removed a bias that was already small and cost a little precision, and on the wetter target it removed most of a large bias and gave back more than that in variance. Two things it was not given are the obvious next tries. A participation model in the logarithms of rescaled moisture and its complement would estimate the exact density ratio that the section used only as a check, and a bootstrap interval would replace the plug-in one where that fell short. Neither changes that on the wetter target the correction rests on a handful of edge plots standing in for parcels wetter than any of them: even the exact ratio left 0.1509 species of the bias there.
The other direction is design. Everything measured here was a property of a trial whose plots were already in the ground. The question asked before that is more interesting: if the point of a trial is to inform a district, how should its plots be spread over the district’s moisture gradient? The transport standard error has an explicit form, so it can be minimised over the design rather than reported after the fact, and that is not the calculation that maximises power for the valley.
References
Pearl J, Bareinboim E 2014 Statistical Science 29(4):579-595 (10.1214/14-STS486)
Bareinboim E, Pearl J 2016 Proceedings of the National Academy of Sciences 113(27):7345-7352 (10.1073/pnas.1510507113)
Cole SR, Stuart EA 2010 American Journal of Epidemiology 172(1):107-115 (10.1093/aje/kwq084)
Stuart EA, Cole SR, Bradshaw CP, Leaf PJ 2011 Journal of the Royal Statistical Society Series A 174(2):369-386 (10.1111/j.1467-985X.2010.00673.x)
Westreich D, Edwards JK, Lesko CR, Stuart EA, Cole SR 2017 American Journal of Epidemiology 186(8):1010-1014 (10.1093/aje/kwx164)
Dahabreh IJ, Robertson SE, Steingrimsson JA, Stuart EA, Hernan MA 2020 Statistics in Medicine 39(14):1999-2014 (10.1002/sim.8426)
Bakker ES, Ritchie ME, Olff H, Milchunas DG, Knops JMH 2006 Ecology Letters 9(7):780-788 (10.1111/j.1461-0248.2006.00925.x)
Catford JA, Wilson JRU, Pysek P, Hulme PE, Duncan RP 2022 Trends in Ecology and Evolution 37(2):158-170 (10.1016/j.tree.2021.09.007)