Minimum wages and teen employment: recovering partial homogeneity

This vignette works through the paper’s first application: the county-level analysis of minimum wages and teen employment from Callaway and Sant’Anna (2021). It is the case where the cohort-time effects carry genuine, recoverable heterogeneity, and the method finds a partially homogeneous structure that roughly halves the variance of the pooled cells while leaving the headline conclusion intact.

The first stage

phdid does not estimate cohort-time effects itself. It takes them from whichever heterogeneity-robust estimator suits the design, together with their joint sampling covariance. Here that is the Callaway-Sant’Anna estimator, whose influence functions give the exact covariance.

Fit with bstrap = FALSE so that the analytical covariance is computed.

first <- att_gt(
  yname = "lemp", tname = "year", idname = "countyreal",
  gname = "first.treat", control_group = "notyettreated",
  data = mpdta, bstrap = FALSE, cband = FALSE
)

d <- ph_data(first)
d
#> <ph_data>  first-stage cohort-time effects for a partition analysis
#> 
#>   Cells (K)      : 7
#>   Source         : did::att_gt() influence functions
#>   Units          : 500
#>   Covariance     : full (cross-cell correlation retained)
#>   sd(tau) / median(se) : 2.25   <- descriptive signal-to-noise
#> 
#>       cell estimate     se weight
#>  2004:2004  -0.0194 0.0223 0.0687
#>  2004:2005  -0.0783 0.0304 0.0687
#>  2004:2006  -0.1363 0.0354 0.0687
#>  2004:2007  -0.1008 0.0344 0.0687
#>  2006:2006   0.0047 0.0163 0.1375
#>  2006:2007  -0.0412 0.0202 0.1375
#>  2007:2007  -0.0261 0.0167 0.4502
#> 
#>   Use homogeneity_test() to ask whether there is heterogeneity to recover,
#>   then l0_ph() or bayes_ph() to recover it.

ph_data() has kept the seven post-treatment cells, taken \(\hat\Sigma = V_{\text{analytical}} / n\), and weighted the cells by cohort size so that the overall effect is population weighted.

The covariance is dense, and that matters. The cells share counties, so their estimates are correlated. Remark 1 of the paper is explicit that using the exact cross-cell covariance is what makes the reported intervals honest: the diagonal shortcut understates the posterior variance of aggregates and can misstate individual co-clustering probabilities badly.

Is there anything to recover?

Ask before grouping. Both estimators will partition pure noise if asked to.

homogeneity_test(d)
#> <homogeneity_test>  is there heterogeneity to recover?
#> 
#> 1. Common-effect test (pooled GLS deviance)   [package addition]
#>      chi-squared = 25.29 on 6 df,  p = 0.000302
#> 
#> 2. Dispersion decomposition                   [paper 5.2.3, adapted]
#>      observed cross-cell sd    : 0.05014
#>      expected under a common effect : 0.02240
#>      heterogeneity sd (excess) : 0.04485
#>      heterogeneity share       : 80%
#>      signal-to-noise           : 2.24
#>      sd(tau)/median(se)        : 2.25   (descriptive)
#> 
#> --------------------------------------------------------------------
#> Reading: recoverable heterogeneity.
#>   The spread across cells clearly exceeds what sampling noise would produce, so
#>   there is structure for the partition model to find. Proceed to l0_ph() for a
#>   point partition and bayes_ph() for inference, and check the co-clustering
#>   matrix to see which groupings are firm.
#> 
#>   Tags mark provenance: [paper 5.2.3] is the published procedure;
#>   [package addition] and [adapted] are this package's. The reading
#>   above is the package's rule, not the paper's. See ?homogeneity_test.

The common-effect null is decisively rejected, and roughly four fifths of the observed cross-cell dispersion is genuine rather than sampling noise. The descriptive signal-to-noise ratio of 2.25 is the paper’s observation that the cross-cell dispersion (0.050) is more than twice the typical standard error (0.022).

A caution the paper makes explicitly: this descriptive ratio is not the separation \(\delta\) of the simulation section, which is defined on the true adjacent-group gaps and the full joint covariance. It signals that there is recoverable heterogeneity without fixing its exact position on that scale. The co-clustering results below assess that directly.

One partition to report

fit_l0 <- l0_ph(d)
fit_l0
#> <l0_ph>  l0-penalised partial-homogeneity estimator
#> 
#>   Cells      : 7
#>   Groups (m) : 4   (selected by BIC, n = 7 in the penalty)
#>   Search     : greedy agglomerative (Appendix B)
#>   Partition  : {2004:2004,2006:2007,2007:2007} {2004:2005,2004:2007} {2004:2006} {2006:2006}
#> 
#>       cell flexible     se group grouped var_ratio
#>  2004:2004  -0.0194 0.0223     1 -0.0285      0.28
#>  2006:2007  -0.0412 0.0202     1 -0.0285      0.33
#>  2007:2007  -0.0261 0.0167     1 -0.0285      0.49
#>  2004:2005  -0.0783 0.0304     2 -0.0871      0.77
#>  2004:2007  -0.1008 0.0344     2 -0.0871      0.61
#>  2004:2006  -0.1363 0.0354     3 -0.1401      0.78
#>  2006:2006   0.0047 0.0163     4  0.0106      0.78
#> 
#>   Mean variance ratio among pooled cells: 0.50
#> 
#>   These intervals condition on the selected partition being correct.
#>   Remark 3: they can under-cover, badly so when the partition is uncertain
#>   (0.57-0.62 in the paper's hard regimes). Use bayes_ph() for inference.

This reproduces Table 6 of the paper. The seven cells collapse into four effect levels: a near-zero group holding the 2004 impact effect and the recent cohorts’ effects, the 2004 cohort’s larger mature effects, its peak, and the lone positive cell.

Pooling roughly halves the variance of the cells that share a group. Notice that even the singleton cells gain precision – in a non-orthogonal design the Gauss-Markov estimator borrows strength across the correlated cells, which is why Proposition 1’s dominance is not confined to groups of size two or more.

plot(fit_l0)

The overall effect is essentially unchanged and slightly more precise:

rbind(
  flexible = aggregate(flex_twfe(d), "overall"),
  l0 = aggregate(fit_l0, "overall")
)
#>             term    estimate  std.error    conf.low   conf.high
#> flexible overall -0.03976363 0.01169526 -0.06268592 -0.01684133
#> l0       overall -0.03882855 0.01153073 -0.06142836 -0.01622874

The substantive conclusion – a roughly 4% teen-employment reduction – survives, while the estimates are summarised by a handful of interpretable levels.

Inference that admits the partition is uncertain

l0_ph() commits to one grouping and its intervals condition on that grouping being correct. With only seven correlated cells that is a strong thing to assume. The Dirichlet Process posterior instead averages over partitions.

fit <- bayes_ph(d, alpha = 1, iters = 6000, burn = 1000, seed = 7,
                progress = FALSE)
fit
#> <bayes_ph>  Dirichlet Process partial-homogeneity posterior
#> 
#>   Cells           : 7
#>   Concentration   : alpha = 1  (prior E[m] = 2.08)
#>   Base measure    : N(-0.04122, 0.2514)
#>   Sampler         : 1 chain(s), 6000 sweeps, 1000 burn-in
#>   Assignment moves: exact covariance
#>   Error variance  : fixed
#> 
#>   Posterior E[# groups] : 2.19
#> 
#>       cell flexible posterior   lower   upper
#>  2004:2004  -0.0194   -0.0107 -0.0393  0.0206
#>  2004:2005  -0.0783   -0.0194 -0.0730  0.0146
#>  2004:2006  -0.1363   -0.0616 -0.1312 -0.0041
#>  2004:2007  -0.1008   -0.0312 -0.1007  0.0071
#>  2006:2006   0.0047   -0.0068 -0.0328  0.0275
#>  2006:2007  -0.0412   -0.0349 -0.0814  0.0013
#>  2007:2007  -0.0261   -0.0216 -0.0553  0.0049
#> 
#>   Intervals are 2.5-97.5% posterior percentiles and already marginalize
#>   over the partition. Use coclustering() to see the grouping structure and
#>   aggregate() for the overall ATT or an event study.

aggregate(fit, "overall")
#>      term    estimate  std.error    conf.low     conf.high
#> 1 overall -0.02391277 0.01190286 -0.04741484 -0.0008181175

Which groupings are firm?

round(coclustering(fit), 2)
#>           2004:2004 2004:2005 2004:2006 2004:2007 2006:2006 2006:2007 2007:2007
#> 2004:2004      1.00      0.78      0.13      0.57      0.86      0.45      0.71
#> 2004:2005      0.78      1.00      0.20      0.70      0.71      0.55      0.72
#> 2004:2006      0.13      0.20      1.00      0.44      0.09      0.54      0.33
#> 2004:2007      0.57      0.70      0.44      1.00      0.50      0.71      0.69
#> 2006:2006      0.86      0.71      0.09      0.50      1.00      0.37      0.64
#> 2006:2007      0.45      0.55      0.54      0.71      0.37      1.00      0.67
#> 2007:2007      0.71      0.72      0.33      0.69      0.64      0.67      1.00
plot(fit)

The small-effect cells recur together with moderate-to-high probability: the 2004 impact effect co-clusters with the lone positive 2006:2006 cell, and with the 2004:2005 effect, at around 0.85, so a low-effect cluster is a robust feature. The large mature 2004 effects are more distinct, with the peak 2004:2006 cell co-clustering with every other cell at roughly 0.5 or below.

(This vignette uses short chains so that it builds quickly. The reported probabilities are therefore a little noisier than the paper’s, which uses 20,000 sweeps; the exact-enumeration comparison below is the check that the sampler is nonetheless landing in the right place.)

The matrix is more diffuse than any single partition. That is not a failure – it is the correct report that with only seven correlated cells the finer grouping is genuinely uncertain. This is the honest reading the co-clustering is meant to deliver: a recurring coarse structure with calibrated uncertainty about the finer detail.

Checking the sampler

With seven cells there are only \(B_7 = 877\) partitions, so the posterior can be computed exactly and the sampler checked against it.

exact <- enumerate_partitions(d, alpha = 1)
exact
#> <ph_enumeration>  exact partition posterior
#> 
#>   Cells           : 7
#>   Partitions       : 877 (enumerated in full)
#>   Concentration   : alpha = 1
#>   Base measure    : N(-0.04122, 0.2514)
#> 
#>   Posterior E[# groups] : 2.195
#>   Overall effect        : -0.0240  [-0.0477, -0.0003]
#> 
#>   Most probable partitions:
#>     0.227  {2004:2004,2004:2005,2004:2007,2006:2006,2006:2007,2007:2007} {2004:2006}
#>     0.102  {2004:2004,2004:2005,2004:2007,2006:2006,2007:2007} {2004:2006,2006:2007}
#>     0.083  {2004:2004,2004:2005,2006:2006} {2004:2006,2004:2007,2006:2007,2007:2007}
#>     0.080  {2004:2004,2004:2005,2004:2006,2004:2007,2006:2006,2006:2007,2007:2007}
#>     0.075  {2004:2004,2004:2005,2006:2006,2007:2007} {2004:2006,2004:2007,2006:2007}

data.frame(
  quantity = c("E[# groups]", "overall effect", "interval lower",
               "interval upper"),
  gibbs = c(fit$m_mean, aggregate(fit, "overall")$estimate,
            aggregate(fit, "overall")$conf.low,
            aggregate(fit, "overall")$conf.high),
  exact = c(exact$m_mean, exact$overall$mean, exact$overall$lower,
            exact$overall$upper)
)
#>         quantity         gibbs         exact
#> 1    E[# groups]  2.1950000000  2.1953344255
#> 2 overall effect -0.0239127733 -0.0239750744
#> 3 interval lower -0.0474148440 -0.0477445245
#> 4 interval upper -0.0008181175 -0.0003249729

cat("largest co-clustering discrepancy:",
    round(max(abs(fit$coclust - exact$coclust)), 4), "\n")
#> largest co-clustering discrepancy: 0.0241

The two agree to within Monte Carlo error, which is the Appendix E benchmark.

It is worth seeing what makes that agreement possible:

cmp <- covariance_check(d, alpha = 1, iters = 1500, burn = 400, seed = 1)
cmp$summary
#>                quantity       exact    diagonal         gap
#> 1        overall effect -0.02439166 -0.02687496 0.002483299
#> 2      overall CI width  0.05116939  0.04538664 0.005782752
#> 3           E[# groups]  2.18363636  2.62545455 0.441818182
#> 4 max co-clustering gap          NA          NA 0.690909091

The diagonal shortcut reproduces the aggregates but moves individual co-clustering probabilities, which is why the exact form is the default.

How much does the prior matter?

With only seven cells, quite a lot for the aggregate. Rather than fix alpha, report the path.

sens <- alpha_sensitivity(
  d, alpha_grid = c(0.1, 0.5, 1, 5, 14, 100),
  iters = 1500, burn = 400, seed = 3
)
sens
#>   alpha    term    estimate    conf.low     conf.high        m  prior_m
#> 1   0.1 overall -0.01703034 -0.04266925  0.0074815216 1.497273 0.426268
#> 2   0.5 overall -0.02181250 -0.04462486  0.0018062038 1.934545 1.354025
#> 3   1.0 overall -0.02467258 -0.04971142 -0.0003261453 2.172727 2.079442
#> 4   5.0 overall -0.02875588 -0.05219768 -0.0055716366 3.070909 4.377344
#> 5  14.0 overall -0.03333214 -0.05827615 -0.0099963448 4.098182 5.676512
#> 6 100.0 overall -0.03823008 -0.05971888 -0.0154249135 6.153636 6.765865
plot_sensitivity(sens)

The effect slides monotonically from about -0.017 under heavy pooling toward the flexible estimate as alpha grows, and is bounded away from zero only for alpha above roughly 1. As a reference point, the BIC pass above selected four groups, which the posterior expected number of groups matches at alpha around 14.

The grid here is coarse and the chains short, again to keep the build fast; widen alpha_grid and raise iters for a publication figure.

That is offered as a reference, not a recipe. Choosing alpha from the data turns the prior into a data-dependent object, with the usual consequences of empirical Bayes – understated posterior uncertainty and a double use of the data – so showing the sensitivity is preferable to committing to a single value.

Watching the individual effects, not just the aggregate

The path above is the overall effect. It is reassuringly stable, but that stability is partly a property of aggregation: the paper notes that the ATT tolerates a wide range of penalties while the cohort-time effects do not. To see what the regularisation is actually doing, ask for the same path per cell.

sens_cells <- alpha_sensitivity(
  d, alpha_grid = c(0.1, 0.5, 1, 5, 14, 100),
  type = "cells", iters = 1500, burn = 400, seed = 5
)
plot_sensitivity(sens_cells)

Every cell starts at the fully pooled value on the left and fans out toward its flexible estimate as alpha grows. Cells whose lines stay together are the ones the model keeps merging; the peak 2004:2006 effect separates earliest, which is the same message the co-clustering matrix gave.

The l0 estimator has the same structure, indexed by its penalty instead. Its path is a step function, so it is easier to read indexed by the number of groups, which visits every partition on the agglomeration path:

lam <- lambda_sensitivity(d, type = "cells", by = "m")
plot_sensitivity(lam)

Reading right to left, this is the agglomeration in reverse: at m = 7 the lines are the seven flexible estimates, and each step merges one more pair until everything meets at the pooled value.

One caution specific to the l0 path. Indexing by lambda rather than by m can skip group counts entirely: the stopping rule compares each merge cost against lambda * |A| * |B|, and that factor grows as groups absorb one another, so the effective threshold is not monotone in the merge order. On this data the thresholds run 0.079, 0.686, 0.397, 1.350, 2.765, 3.571 – so no value of lambda selects five groups. Use by = "m" when you want the complete path.

Note also that these l0 bands narrow as the penalty tightens, because they condition on the selected partition being correct. The Bayesian bands above do not have this problem: they average over partitions.

The event-study reading

aggregate(fit, "dynamic")
#>   term    estimate  std.error    conf.low     conf.high
#> 1  e=0 -0.01736403 0.01155017 -0.04077529  0.0048477245
#> 2  e=1 -0.02976200 0.01775717 -0.06710489  0.0006997113
#> 3  e=2 -0.06163171 0.03335810 -0.13120197 -0.0041376141
#> 4  e=3 -0.03118668 0.02719331 -0.10074599  0.0071489075

The effect is small at impact and accumulates over the horizon, the familiar pattern of a minimum-wage effect that builds over several years.

One caveat on interpretation: the composition of the event-time aggregate changes across the horizon. Event times 2 and 3 are identified off the 2004 cohort alone, so the longer-run trajectory is that cohort’s, and the design cannot say whether later cohorts would follow it. The impact-year effects, by contrast, are small and similar across cohorts. The partial homogeneity the method recovers here is largely the homogeneity of the short-run response across cohorts, with the accumulating longer-run effects carried by the 2004 cohort’s later cells kept separate.