diff --git a/.Rbuildignore b/.Rbuildignore index 9bff1d1..e617e8d 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -20,3 +20,7 @@ ^bml-brms-alignment-plan\.md$ ^micro_macro_framework_slides\.qmd$ ^data/.*\.feather$ +^vignettes/\.fits +^vignettes/examples\.html$ +^vignettes/\.fits +^vignettes/examples\.html$ diff --git a/.gitignore b/.gitignore index 39cdd5e..1adf676 100644 --- a/.gitignore +++ b/.gitignore @@ -15,3 +15,9 @@ inst/doc /.quarto/ /backup/ data/*.feather +vignettes/.fits/ +vignettes/examples.html +.Renviron +vignettes/.fits/ +vignettes/examples.html +.Renviron diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 0000000..8403fcb --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,5 @@ +{ + "files.associations": { + "*.Rmd.orig": "rmd" + } +} diff --git a/vignettes/examples-README.md b/vignettes/examples-README.md new file mode 100644 index 0000000..61a0278 --- /dev/null +++ b/vignettes/examples-README.md @@ -0,0 +1,308 @@ +# Building the `examples` vignette + +This note explains how the Examples vignette is built, what data it needs, and where that +data has to live. It is written for someone picking the repository up fresh. + +The short version: the vignette is **precompiled**. The committed `examples.Rmd` contains no +executable code, so it renders for anyone with no data and no JAGS. You only need the data +and the setup below if you want to change a model and rebuild. + +--- + +## 1. Why the vignette is precompiled + +Fitting the fourteen models in this vignette takes about half an hour of JAGS time, and two +of the datasets cannot be redistributed. If the vignette were an ordinary `.Rmd`, every +`R CMD build`, every pkgdown run and every CRAN check would try to refit everything and fail +without the data. + +Precompiling solves both problems. The models are fit once, cached, and the results are baked +into a committed artifact. + +``` +fit_models.R → vignettes/.fits/*.rds → examples.Rmd.orig → precompile.R → examples.Rmd + (JAGS, ~30 min) (git-ignored) (the source) (seconds) (committed) +``` + +Three files do the work, all in `vignettes/`: + +| File | Role | Needs data? | Needs JAGS? | +|---|---|---|---| +| `fit_models.R` | fits every model once, caches to `.fits/` | yes | yes | +| `examples.Rmd.orig` | the **source** vignette; loads cached fits, never fits anything | no | no | +| `precompile.R` | knits `.orig` into the committed `examples.Rmd` | no | no | + +**The committed `examples.Rmd` is a build artifact. Never edit it by hand.** Every code block +in it is inert display text, so an edit there is silently discarded the next time anyone runs +`precompile.R`. Edit `examples.Rmd.orig` instead. + +--- + +## 2. What data the vignette uses + +Four datasets, two of which ship and two of which do not. + +| Example | Data | Ships? | How it is obtained | +|---|---|---|---| +| 1 | CILS4EU Germany, Wave 1 and 2 | **no** | restricted access, see §3 | +| 2 | Boston 1970 census tracts | yes | `spData` package | +| 3, 7 | `coalgov` | yes | this package | +| 4, 5, 6 | O\*NET occupation-activity structure | **no** | built by a Python script, see §4 | + +Examples 2, 3 and 7 therefore rebuild with no setup at all. Examples 1 and 4 through 6 need +the two external datasets. + +--- + +## 3. CILS4EU (Example 1) + +Three Feather files, exactly these names: + +``` +nodedat.feather Wave 1 student records (smoking, gender) +edgedat.feather friendship nominations with closeness rank +nodedat-w2.feather Wave 2 student records (the outcome) +``` + +Put them in one directory **outside this repository** and point at it with an environment +variable. See §5. + +CILS4EU is restricted-access data and must not be committed. `fit_models.R` has no default +path for it, by design. An earlier version fell back to `"data"`, which resolves to the +package's own `data/` directory when run from the repository root, and that directory is +tracked by git and shipped to users. The fallback was removed so a missing environment +variable fails loudly rather than quietly writing restricted data somewhere public. + +--- + +## 4. O\*NET (Examples 4, 5 and 6) + +These examples pair O\*NET's occupation-activity structure with activity-level AI exposure +scores and BLS wage and employment series. The dataset is assembled by +`build_onet_ai_dataset.py`. + +### 4.1 What the vignette reads + +Two Parquet files: + +``` +long_occupation_member_dwa.parquet one row per occupation-activity pair +model_occupation_year_dwa.parquet one row per occupation +``` + +`fit_models.R` reads both and joins them. The long table carries the membership structure, +the importance weights and each activity's AI exposure. The occupation-level outcomes and +moderators live only in the model table, so both files are required. + +**Use the `_dwa` build, not `_task`.** The script writes both. At the task level each task +belongs to exactly one occupation, so there is no multiple membership and the `mm()` term +would add nothing over an ordinary occupation-level regression. At the DWA level the same +work activity recurs across occupations, which is the structure the model needs. The script's +own diagnostic reports this: `multiple_membership_present` is `true` for `_dwa` and `false` +for `_task`, and 97% of activities appear in more than one occupation. + +### 4.2 Regenerating the Parquet files + +Only needed if the underlying sources change. The script expects these raw files: + +``` +onet/Task Ratings.txt task importance ratings, the weights +onet/Tasks to DWAs.txt links each task to its detailed work activities +onet/Job Zones.txt occupation preparation band +onet/Work Activities.txt activity names, used for logging only +OpenAI/full_labelset.tsv Eloundou et al. task-level AI exposure scores +BLS/oesm22nat/national_M2022_dl.xlsx base-year wages and employment +BLS/oesm25nat/national_M2025_dl.xlsx later-year wages and employment +BLS/occupation.xlsx employment projections, sheet "Table 1.2" +``` + +Set the two paths at the top of the script and run it: + +```python +"RAW_DIR": "/abs/path/to/onet/raw_data", +"OUT_DIR": "/abs/path/to/onet/processed_data", +``` + +```bash +python build_onet_ai_dataset.py +``` + +It writes eight files to `OUT_DIR`, including `data_checks_dwa.json`. That file is worth +opening: it reports the counts the vignette's prose asserts, so the text can be checked +against the build without running any models. + +| Vignette says | JSON key | +|---|---| +| 17,537 occupation-activity pairs | `n_member_rows` | +| 894 occupations | `n_occupations` | +| 2,080 activities | `multiple_membership.n_distinct_members` | +| about 20 activities per occupation | `multiple_membership.members_per_occupation_mean` | +| about 8 occupations per activity | `multiple_membership.occupations_per_member_mean` | +| up to 91 | `multiple_membership.occupations_per_member_max` | + +`occupations_matched_to_oews` is 854, which is the N reported in Examples 4 and 5. Nothing is +lost beyond the OEWS wage match. + +--- + +## 5. Where to put the data and how to point at it + +Keep both datasets outside the repository. Nothing then depends on an ignore rule holding. + +``` +~/research-data/ +├── cils/ +│ ├── nodedat.feather +│ ├── edgedat.feather +│ └── nodedat-w2.feather +└── onet/ + ├── raw_data/ # only if regenerating the parquets + └── processed_data/ # the two _dwa parquets live here +``` + +Then set two environment variables in `~/.Renviron`, which R reads at startup: + +``` +BML_CILS_DIR=/abs/path/to/research-data/cils +ONET_DATA_DIR=/abs/path/to/research-data/onet/processed_data +``` + +No quotes, no spaces around `=`, and absolute paths rather than `~`, which `arrow` does not +expand. Restart R, then confirm: + +```r +Sys.getenv(c("BML_CILS_DIR", "ONET_DATA_DIR")) +``` + +Both are required. `fit_models.R` stops with a named message if either is missing. + +Two things about that file are worth knowing in advance, because both fail silently. + +`~/.Renviron` means the file in your home directory. If you create it from a shell with the +path quoted, as in `"~/.Renviron"`, the tilde does not expand and you get a literal directory +named `~` inside whatever directory you happened to be in, with the file inside it. R never +looks there. The symptom is simply that both variables come back empty, with nothing to +indicate the file exists in the wrong place. `cat /Users/yourname/.Renviron` confirms it is +where you think it is. The easiest way to avoid the problem is +`Rscript -e 'usethis::edit_r_environ()'`, which opens the right file. + +R reads `.Renviron` once, at startup. Editing it in a running session has no effect, so +restart R before checking `Sys.getenv()`. + +--- + +## 6. Requirements + +**JAGS** must be installed at system level, plus the `rjags` R package. + +For `fit_models.R`: `dplyr`, `spData`, `sf`, `spdep`, `bml`, `arrow`, `loo`. + +For `precompile.R`, additionally: `spatialreg`, `ggplot2`, `bayesplot`, `tidyr`, `tibble`, +`posterior`, `knitr`. + +`bayesplot` is easy to miss. Without it the posterior predictive figure in Example 1 renders +as an error message rather than a plot, and nothing else fails, so the build looks clean. + +--- + +## 7. Rebuilding + +```bash +cd /path/to/bml + +# 1. clear stale figures. REQUIRED if any chunk was added or removed. +# Figures are named by chunk position, so adding a chunk renames every +# downstream image. Old files keep their old names, the new tags point +# at new ones, and you get broken images with nothing in the build log. +rm -rf vignettes/examples-figures/ + +# 2. fit. About 30 minutes. Needs JAGS and both environment variables. +Rscript vignettes/fit_models.R + +# 3. precompile. Seconds. No JAGS, no data. +Rscript -e 'source("vignettes/precompile.R")' + +# 4. render, to check it +Rscript -e 'rmarkdown::render("vignettes/examples.Rmd")' +``` + +Step 2 writes 24 files into `vignettes/.fits/`: fourteen fitted models, seven `loo` objects, +one posterior predictive sample, and two raw-data caches. It stops early with a specific +message if a required column is missing from the O\*NET tables, or if the `smax` shape +parameter did not survive caching. + +If step 4 fails on pandoc from a plain shell, RStudio bundles a copy: + +```bash +export RSTUDIO_PANDOC="/Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools" +``` + +Knitting in RStudio with Ctrl+Shift+K works too. Note that `.Rmd` files should be knitted +rather than previewed through Quarto, since this is an R Markdown vignette and `quarto +preview` rejects it. + +Only steps 1, 3 and 4 are needed for a prose-only change. Step 2 is required only when a +model's specification, data or seed changes. + +--- + +## 8. Checking a rebuild + +The prose quotes numbers that the rebuild has to reproduce. If any of these move, something +in the setup differs. + +| Check | Expected | +|---|---| +| Ex 1 DIC | 18,700 / 17,577 / 16,964 | +| Ex 2 DIC | 292 / 229 | +| Ex 2 social dissimilarity `b1` | 0.592 [0.403, 0.797] | +| Ex 3 DIC | 3864 / 3868 / 3854 | +| Ex 4 exposure spread `V_ai_exposure` | -0.096 [-0.165, -0.027] | +| Ex 5 `fn[kappa]` | -1.57 [-2.73, -0.40] | +| Ex 6 `Aexp:education` | -0.009 [-0.025, 0.007] | + +A few figures in the prose appear in no table, mainly random-effect standard deviations. +`vignettes/check_numbers.R` reads them straight out of the cached fits and prints them, so +they can be verified without refitting anything. + +Open the rendered HTML and confirm four figures appear. A missing figure is the failure mode +that step 1 prevents and it is silent in the log. + +--- + +## 9. What is committed and what is not + +Committed: `examples.Rmd.orig`, `fit_models.R`, `precompile.R`, `examples.Rmd`, and +`examples-figures/`. + +Not committed: `vignettes/.fits/`. Two independent rules keep it out, and both are needed. + +- `.gitignore` has `vignettes/.fits/` +- `.Rbuildignore` has `^vignettes/\.fits` + +`R CMD build` copies the whole `vignettes/` directory into the source tarball and does not +read `.gitignore`, so without the second rule the cache ships inside the tarball. That matters +because the cache is not only fitted parameters. `cils_raw.rds` and `onet_raw.rds` are full +copies of the source tables, and `ppc_lim.rds` holds the observed outcome vector for 4,002 +students. + +Worth running once after any change to those rules: + +```bash +R CMD build . && tar -tzf bml_*.tar.gz | grep -Ei 'fits|parquet|feather|raw_data' +``` + +It should return nothing. + +--- + +## 10. A known limitation + +`coalgov` ships, so Examples 3 and 7 are fully reproducible by anyone. Examples 1 and 4 +through 6 are not, because their data cannot be redistributed. The rendered vignette is +complete for readers, but a user cannot rerun `fit_models.R` themselves. + +For the O\*NET examples this is fixable if desired. A trimmed occupation-activity table with +only the columns the models use compresses to well under 1 MB as an `.rda`, which would make +Examples 4 through 6 as reproducible as Example 3. The CILS4EU restriction is not something +the package can work around. diff --git a/vignettes/examples-figures/monet-kappa-1.png b/vignettes/examples-figures/monet-kappa-1.png new file mode 100644 index 0000000..4feac63 Binary files /dev/null and b/vignettes/examples-figures/monet-kappa-1.png differ diff --git a/vignettes/examples-figures/ppc-lim-1.png b/vignettes/examples-figures/ppc-lim-1.png new file mode 100644 index 0000000..e437fcd Binary files /dev/null and b/vignettes/examples-figures/ppc-lim-1.png differ diff --git a/vignettes/examples-figures/unnamed-chunk-12-1.png b/vignettes/examples-figures/unnamed-chunk-12-1.png index 3875260..e437fcd 100644 Binary files a/vignettes/examples-figures/unnamed-chunk-12-1.png and b/vignettes/examples-figures/unnamed-chunk-12-1.png differ diff --git a/vignettes/examples-figures/unnamed-chunk-19-1.png b/vignettes/examples-figures/unnamed-chunk-19-1.png deleted file mode 100644 index fb871f8..0000000 Binary files a/vignettes/examples-figures/unnamed-chunk-19-1.png and /dev/null differ diff --git a/vignettes/examples-figures/unnamed-chunk-33-1.png b/vignettes/examples-figures/unnamed-chunk-33-1.png index 6226c5a..3c64c61 100644 Binary files a/vignettes/examples-figures/unnamed-chunk-33-1.png and b/vignettes/examples-figures/unnamed-chunk-33-1.png differ diff --git a/vignettes/examples-figures/unnamed-chunk-41-1.png b/vignettes/examples-figures/unnamed-chunk-41-1.png deleted file mode 100644 index 22df8e0..0000000 Binary files a/vignettes/examples-figures/unnamed-chunk-41-1.png and /dev/null differ diff --git a/vignettes/examples-figures/weight-dissim-1.png b/vignettes/examples-figures/weight-dissim-1.png new file mode 100644 index 0000000..4101ecc Binary files /dev/null and b/vignettes/examples-figures/weight-dissim-1.png differ diff --git a/vignettes/examples.Rmd b/vignettes/examples.Rmd index 0e6f1be..c5eaebb 100644 --- a/vignettes/examples.Rmd +++ b/vignettes/examples.Rmd @@ -20,7 +20,7 @@ editor_options: This vignette presents seven examples that show how the multiple-membership multilevel model (MMMM) can be applied to network, spatial, and aggregation problems. The first three cover the additive micro-macro link: peer influence in a friendship network, neighbourhood effects on home values, and how coalition parties shape government survival. -The last four continue the coalition example to showcase aggregation beyond the weighted mean: distributional features (`fn("var")`), estimating the aggregation function itself (`fn("smax", kappa = est())`), cross-level interactions through named blocks, and heterogeneous member effects (`re(1 + x)`). +The last four go beyond the weighted mean to showcase distributional features (`fn("var")`), estimating the aggregation function itself (`fn("smax", kappa = est())`), cross-level interactions through named blocks, and heterogeneous member effects (`re(1 + x)`). # 1. All friends, or just your best friend? (network regression) @@ -45,50 +45,53 @@ The model uses the following variables: | `rank` | weight | Friendship closeness (1 = best friend, 5 = fifth-closest) | | `n` | weight | Number of friends nominated by the ego | -Two models are fit and compared. -The economics standard for this kind of network problem is the linear-in-means model: every friend carries equal weight, so what matters is the average across all of them, implying that each friend exerts the same influence. -The alternative is that influence is concentrated: a single close friend drives most of the effect while looser ties matter little. -`bml` fits both and compares them directly, changing only one line of the weight function. +Three models are fit and compared, tracing a spectrum from the most concentrated weighting rule to the most diffuse. +At one end, influence is concentrated: a single close friend drives the effect while looser ties matter little. +At the other is the economics standard for this kind of network problem, the linear-in-means model: every friend carries equal weight, so what matters is the average across all of them. +Between the two sits a graded rule in which closeness matters but no friend is discarded. +All three are fixed weighting rules, nothing inside `w()` is estimated, so the comparison isolates the shape of the aggregation, changing only one line of the weight function each time. -## Model 1: linear-in-means +## Model 1: best friend only -Every friend receives equal weight, `w ~ 1/n`, so the peer term is the simple average of all nominated friends' Wave 1 smoking. +Only the closest-ranked friend contributes; all others are dropped. +`rank == min(rank)` evaluates to `TRUE` for the closest-ranked friend within each ego's friend set (rank 1), placing all weight on that friend. +With `scale = TRUE` the weight is renormalised to 1 within each ego, so the peer term becomes that single friend's smoking. ``` r -mod_lim <- +mod_bf <- bml( - smoking_w2 ~ smoking_ego + gender_ego + # outcome ~ ego-level predictors - mm( # multiple-membership term (the peers) - id = id(alter, ego), # member = alter (friend), group = ego - vars = vars(smoking_alter), # member variable to aggregate - w = w(~ 1/n, scale = TRUE), # equal weights, renormalised to sum to 1 - fn = fn("sum"), # aggregation function: the weighted mean - RE = TRUE # member-level random effect + smoking_w2 ~ smoking_ego + gender_ego + + mm( + id = id(alter, ego), + vars = vars(smoking_alter), + w = w(~ rank == min(rank), scale = TRUE), # all weight on the closest friend + fn = fn("sum"), # aggregation function: the weighted mean + RE = TRUE ), - family = gaussian(), # continuous outcome + family = gaussian(), data = cils_long ) ``` -## Model 2: best friend only +## Model 2: graded decay -Only the closest-ranked friend contributes; all others are dropped. -`rank == min(rank)` evaluates to `TRUE` for the closest-ranked friend within each ego's friend set (rank 1), placing all weight on that friend. -With `scale = TRUE` the weight is renormalised to 1 within each ego, so the peer term becomes that single friend's smoking. +The best-friend rule is an extreme as it throws away four of the five nominations. +A softer version keeps every friend but lets closeness matter (`w ~ 1/rank`) +gives the best friend weight 1, the second-closest 1/2, and so on down to 1/5. ``` r -mod_bf <- +mod_decay <- bml( smoking_w2 ~ smoking_ego + gender_ego + mm( id = id(alter, ego), vars = vars(smoking_alter), - w = w(~ rank == min(rank), scale = TRUE), # all weight on the closest friend - fn = fn("sum"), # same aggregation function as Model 1 + w = w(~ 1/rank, scale = TRUE), # closeness-graded weights + fn = fn("sum"), RE = TRUE ), family = gaussian(), @@ -98,15 +101,39 @@ mod_bf <- -## Comparing the two models +## Model 3: linear-in-means + +Every friend receives equal weight, `w ~ 1/n`, so the peer term is the simple average of all nominated friends' Wave 1 smoking. + + +``` r +mod_lim <- + bml( + smoking_w2 ~ smoking_ego + gender_ego + # outcome ~ ego-level predictors + mm( # multiple-membership term (the peers) + id = id(alter, ego), # member = alter (friend), group = ego + vars = vars(smoking_alter), # member variable to aggregate + w = w(~ 1/n, scale = TRUE), # equal weights, renormalised to sum to 1 + fn = fn("sum"), # aggregation function: the weighted mean + RE = TRUE # member-level random effect + ), + family = gaussian(), # continuous outcome + data = cils_long + ) +``` + + + +## Comparing the three models `bmlCompare()` lays out the coefficient estimates and fit statistics of several models side by side, one column per model, with readable labels supplied through `labels`. ``` r bmlCompare( - "Linear-in-means" = mod_lim, "Best friend only" = mod_bf, + "Graded decay" = mod_decay, + "Linear-in-means" = mod_lim, terms = c("(Intercept)", "smoking_ego", "gender_ego", "A_smoking_alter"), labels = c("Intercept", "Own smoking (W1)", "Gender", "Peers' smoking (W1)") ) @@ -114,23 +141,28 @@ bmlCompare( -|Term |Linear-in-means |Best friend only | -|:-------------------|:---------------------|:---------------------| -|Intercept |2.556 (2.345, 2.770) |2.623 (2.448, 2.793) | -|Own smoking (W1) |0.365 (0.337, 0.390) |0.360 (0.334, 0.386) | -|Gender |0.047 (-0.042, 0.133) |0.059 (-0.025, 0.144) | -|Peers' smoking (W1) |0.042 (0.001, 0.084) |0.026 (-0.000, 0.052) | -|N |4002 |4002 | -|DIC |17029 |18463 | +|Term |Best friend only |Graded decay |Linear-in-means | +|:-------------------|:---------------------|:---------------------|:---------------------| +|Intercept |2.622 (2.448, 2.791) |2.567 (2.362, 2.768) |2.557 (2.346, 2.770) | +|Own smoking (W1) |0.360 (0.334, 0.386) |0.362 (0.336, 0.389) |0.365 (0.337, 0.390) | +|Gender |0.059 (-0.025, 0.144) |0.048 (-0.044, 0.136) |0.047 (-0.042, 0.133) | +|Peers' smoking (W1) |0.026 (-0.000, 0.052) |0.042 (0.002, 0.079) |0.042 (0.001, 0.083) | +|N |4002 |4002 |4002 | +|DIC |18700 |17577 |16964 | -Both models find a positive peer effect: students whose friends smoke more tend to smoke more themselves one year later, over and above their own Wave 1 smoking. -The own-smoking coefficient (about 0.36 in both models) dominates, which is expected since past behaviour is the strongest predictor of future behaviour. +All three models find a positive peer effect: students whose friends smoke more tend to smoke more themselves one year later, over and above their own Wave 1 smoking. +The own-smoking coefficient (about 0.36 across the board) dominates, which is expected since past behaviour is the strongest predictor of future behaviour. -The peer effect is present but modest. -Under linear-in-means, averaging across all nominated friends gives a coefficient of about 0.04 (CI [0.00, 0.08]), which just clears zero. -Limiting to the best friend gives about 0.03 (CI [-0.00, 0.05]): a smaller point estimate with a tighter interval that now just touches zero. +The peer effect itself depends on how the friend set is aggregated. +Restricting attention to the best friend gives about 0.03 with an interval that touches zero ([-0.00, 0.05]). +Both of the rules that keep every friend give about 0.04 and clear zero (graded decay [0.00, 0.08], linear-in-means [0.00, 0.08]). +Discarding the looser ties does not sharpen the estimate; it loses signal. -On fit, the linear-in-means model is clearly preferred (DIC 17,029 vs 18,463, a gap of over 1,400). +The fit statistics fall monotonically as the weighting becomes more diffuse: DIC drops from 18,700 for the best-friend rule to 17,577 for graded decay and 16,964 for linear-in-means. +Each relaxation of the concentration assumption buys a substantial improvement, and the largest single gain comes from not throwing the other four friends away. + +The member-level standard deviation tells the same story from another angle, rising from 0.40 under best-friend weighting to 0.83 and then 0.93. +Under the best-friend rule most nominated friends receive zero weight in every group they belong to, so the data carry little information about their individual effects and the estimated spread shrinks toward the prior. Peer influence on smoking appears to be spread across the friend set rather than concentrated in the single closest tie. ## Cross-validation: which aggregation does the data prefer? @@ -140,24 +172,27 @@ For gaussian outcomes `bml` monitors the pointwise log-likelihood in the generat ``` r -loo_lim <- loo(mod_lim) -loo_bf <- loo(mod_bf) -loo::loo_compare(loo_lim, loo_bf) +loo_bf <- loo(mod_bf) +loo_decay <- loo(mod_decay) +loo_lim <- loo(mod_lim) +loo::loo_compare(loo_bf, loo_decay, loo_lim) ``` ``` -#> elpd_diff se_diff -#> model1 0.0 0.0 -#> model2 -37.1 9.8 +#> model elpd_diff se_diff p_worse diag_diff diag_elpd +#> model3 0.0 0.0 NA 8 k_psis > 0.7 +#> model2 -6.3 5.1 0.89 3 k_psis > 0.7 +#> model1 -37.1 9.7 1.00 ``` -The linear-in-means model comes out on top: the best-friend model loses about 37 points of expected log predictive density (standard error about 10), agreeing with the DIC ranking. -Both criteria point the same way — influence is spread over the whole friend set rather than concentrated in the closest tie. +model1 is best-friend, model2 graded decay, model3 linear-in-means, in the order passed. +PSIS-LOO ranks them as DIC does: the best-friend model loses about 37 points of expected log predictive density (standard error 10, p_worse 1.00), while graded decay sits only 6 points behind linear-in-means (standard error 5, p_worse 0.89), decisive at the bottom of the ladder, close at the top. +The Pareto-k flags bias elpd_diff toward whichever model carries more of them (8 for linear-in-means, 3 for graded decay), so that narrow margin is best read as an upper bound. +Both criteria point the same way: influence is spread over the friend set rather than concentrated in the closest tie. -A posterior predictive check makes sure the preferred model reproduces the shape of the outcome at all. `pp_check()` overlays replicated outcome distributions on the observed one: @@ -165,7 +200,7 @@ A posterior predictive check makes sure the preferred model reproduces the shape pp_check(mod_lim) ``` -Posterior predictive check: replicated outcome densities overlaid on the observed smoking distribution +Posterior predictive check: replicated outcome densities overlaid on the observed smoking distribution The replicated densities track the observed distribution's spread and skew; the discreteness of the smoking scale shows up as bumps the gaussian model smooths over, which is worth knowing but does not affect the aggregation comparison. @@ -220,84 +255,96 @@ mod_bml <- RE = TRUE # neighbour-level random effect ), family = gaussian(), + iter = 5000, warmup = 500, chains = 3, seed = 1, data = boston_df ) ``` -## Parameterised weights: similarity across covariates +## Parameterised weights: social similarity -Equal weighting assumes every neighbour matters the same regardless of how similar it is to the focal tract. -A natural alternative is that tracts are more strongly influenced by neighbours that resemble them, a spatial homophily effect. -Rather than collapse similarity into a single number, the weight is allowed to depend on the dissimilarity in each covariate separately, using the functional form from Rosche (2026): +Equal weighting assumes every adjacent tract matters the same, regardless of how +much it resembles the focal tract. +A natural alternative is that spillovers run more strongly between tracts that +are alike socially, social distance layered on top of physical adjacency. +One parameter is enough to test it: $$ -w_{ij} = \frac{1}{1 + (n_i - 1)\exp\!\bigl(-(b_0 + b_1 \cdot d^{\text{CRIM}}_{ij} + b_2 \cdot d^{\text{AGE}}_{ij})\bigr)} +w_{ij} \propto \exp\!\bigl(-b_1 \cdot d^{\text{LSTAT}}_{ij}\bigr) $$ -Each `d` is the standardised absolute difference between a focal tract and its neighbour on that covariate. -When all the `b`'s are zero this collapses to `1/n_i`, the equal-weight baseline. -A negative coefficient means neighbours that are similar on that covariate (low `d`) receive more weight; a positive coefficient means dissimilar neighbours dominate. -Letting the two dissimilarities enter on their own lets the model decide which kind of similarity drives the aggregation, instead of imposing a single combined measure. +where $d^{\text{LSTAT}}_{ij}$ is the standardised absolute difference between +tract $i$ and neighbour $j$ in the share of lower-status population. +At $b_1 = 0$ the weights collapse to $1/n_i$ and the equal-weight baseline is +recovered exactly, so the baseline is nested inside this model. +Positive $b_1$ means socially similar neighbours carry more weight. +`LSTAT` appears nowhere else in the model — not among the own-tract covariates, +not in `vars()`, so it acts purely as an aggregation variable. ``` r -mod_bml_w <- +mod_bml_lstat <- bml( lnCMEDV ~ NOX + CRIM + RM + DIS + AGE + mm( id = id(tid_nb, tid), vars = vars(NOX_nb + CRIM_nb + RM_nb + DIS_nb + AGE_nb), - # weight as a function of covariate dissimilarity; nests 1/n when all b's are 0 - w = w(~ 1 / (1 + (n - 1) * exp(-(b0 + b1 * d_CRIM + b2 * d_AGE))), scale = TRUE), + # weight falls with social dissimilarity; nests 1/n at b1 = 0 + w = w(~ exp(-b1 * d_LSTAT), scale = TRUE), fn = fn("sum"), RE = TRUE ), family = gaussian(), - prior = prior(normal(0, 1), class = "w"), # weakly informative; prevents exp() overflow - data = boston_df2 + prior = prior(normal(0, 1), class = "w"), # weakly informative; keeps exp() stable + iter = 5000, warmup = 500, chains = 3, seed = 1, + data = boston_df ) ``` -The estimated weight-function parameters are: - ``` r bmlCompare( - "Similarity weights" = mod_bml_w, + "Social similarity" = mod_bml_lstat, component = "weights", - terms = c("w[b0] (mm.1)", "w[b1] (mm.1)", "w[b2] (mm.1)"), - labels = c("Weight intercept (b0)", "Crime dissimilarity (b1)", "Age dissimilarity (b2)") + terms = c("w[d_LSTAT] (mm.1)"), + labels = c("Social dissimilarity (b1)") ) ``` -|Term |Similarity weights | -|:------------------------|:-----------------------| -|Weight intercept (b0) |1.561 (-0.282, 2.787) | -|Crime dissimilarity (b1) |-1.758 (-2.588, -1.165) | -|Age dissimilarity (b2) |0.713 (0.074, 1.523) | -|N |506 | -|DIC |295 | +|Term |Social similarity | +|:-------------------------|:--------------------| +|Social dissimilarity (b1) |0.592 (0.403, 0.797) | +|N |506 | +|DIC |229 | -The crime-dissimilarity coefficient is about −1.76 (CI [−2.59, −1.16]), clearly negative, so neighbours with similar crime levels carry more weight. -The age-dissimilarity coefficient is about 0.71 (CI [0.07, 1.52]), pointing the other way: neighbours with different age profiles carry more weight. -The model picks up two kinds of similarity at once, in opposite directions. +The estimate is about 0.59 with a 95% interval of [0.40, 0.80], comfortably +away from zero. +Because `d_LSTAT` is standardised, a one-standard-deviation increase in social +dissimilarity multiplies a neighbour's unnormalised weight by +$\exp(-0.59) \approx 0.55$, roughly halving its influence. +Adjacency alone does not determine who counts: among neighbours that all share +a border, the socially similar ones dominate the aggregate. ## Do the weights actually vary with similarity? -A direct check plots the estimated weight each neighbour received against the dissimilarity between that neighbour and its focal tract. -If crime homophily drives the weights, low-`d_CRIM` pairs should cluster toward high weights and high-`d_CRIM` pairs toward low weights. +A direct check plots the estimated weight each neighbour received against its +social dissimilarity from the focal tract. +If social homophily drives the weights, low-`d_LSTAT` pairs should sit toward +high weights and high-`d_LSTAT` pairs toward low weights. -Estimated neighbour weight plotted against covariate dissimilarity, faceted by crime and age dissimilarity +Estimated neighbour weight plotted against social dissimilarity -The two panels show the similarity effects acting in opposite directions. -The weight falls as crime dissimilarity rises (the homophily the negative crime coefficient pointed to), while it rises gently with age dissimilarity, matching the positive age coefficient. -Each panel holds only a partial view, since every neighbour's weight reflects both dissimilarities at once, but the opposing slopes are visible. +The estimated weight declines with social dissimilarity: neighbours close to the +focal tract on lower-status share carry the most weight, and the weight falls +away as dissimilarity grows. +The scatter around the trend is renormalisation at work, a pair's final weight +depends on how many neighbours it competes with within its tract, but the +downward slope, the homophily the positive `b1` pointed to, is unmistakable. ## Benchmark: spatial Durbin error model @@ -320,35 +367,39 @@ The two bml models and the SDEM benchmark are placed in one table. The SDEM's neighbour-covariate (`lag.`) terms are matched to the bml neighbour effects, and all rows are given readable labels. -|Term |Equal weights (1/n) |Similarity weights |SDEM (benchmark) | -|:------------------------|:-----------------------|:-----------------------|:-----------------------| -|Intercept |3.033 (2.997, 3.070) |3.022 (2.987, 3.057) |3.015 (2.950, 3.080) | -|Air quality (NOX) |-0.106 (-0.161, -0.048) |-0.075 (-0.131, -0.015) |-0.112 (-0.160, -0.064) | -|Crime |-0.061 (-0.082, -0.041) |-0.055 (-0.080, -0.030) |-0.067 (-0.086, -0.047) | -|Rooms |0.157 (0.137, 0.177) |0.165 (0.147, 0.182) |0.169 (0.150, 0.187) | -|Distance to employment |-0.096 (-0.208, 0.015) |-0.071 (-0.164, 0.027) |-0.086 (-0.169, -0.003) | -|Pre-1940 share |-0.095 (-0.129, -0.062) |-0.103 (-0.131, -0.075) |-0.088 (-0.118, -0.057) | -|Air quality (neighbours) |0.020 (-0.072, 0.111) |-0.001 (-0.090, 0.085) |0.012 (-0.076, 0.101) | -|Crime (neighbours) |-0.132 (-0.184, -0.082) |-0.143 (-0.212, -0.074) |-0.103 (-0.156, -0.050) | -|Rooms (neighbours) |0.073 (0.025, 0.121) |0.092 (0.049, 0.136) |0.040 (-0.006, 0.085) | -|Distance (neighbours) |-0.012 (-0.145, 0.117) |-0.026 (-0.144, 0.086) |-0.038 (-0.148, 0.072) | -|Pre-1940 (neighbours) |0.007 (-0.066, 0.080) |0.014 (-0.051, 0.079) |-0.020 (-0.094, 0.055) | -|Weight intercept (b0) |— |1.561 (-0.282, 2.787) |— | -|Crime dissimilarity (b1) |— |-1.758 (-2.588, -1.165) |— | -|Age dissimilarity (b2) |— |0.713 (0.074, 1.523) |— | +|Term |Equal weights (1/n) |Social similarity |SDEM (benchmark) | +|:-------------------------|:-----------------------|:-----------------------|:-----------------------| +|Intercept |3.033 (2.997, 3.070) |3.033 (2.999, 3.067) |3.015 (2.950, 3.080) | +|Air quality (NOX) |-0.105 (-0.161, -0.049) |-0.059 (-0.115, -0.004) |-0.112 (-0.160, -0.064) | +|Crime |-0.061 (-0.082, -0.040) |-0.052 (-0.072, -0.032) |-0.067 (-0.086, -0.047) | +|Rooms |0.157 (0.136, 0.178) |0.130 (0.110, 0.151) |0.169 (0.150, 0.187) | +|Distance to employment |-0.099 (-0.205, 0.011) |-0.080 (-0.186, 0.024) |-0.086 (-0.169, -0.003) | +|Pre-1940 share |-0.095 (-0.128, -0.062) |-0.071 (-0.103, -0.038) |-0.088 (-0.118, -0.057) | +|Air quality (neighbours) |0.019 (-0.072, 0.108) |-0.022 (-0.107, 0.060) |0.012 (-0.076, 0.101) | +|Crime (neighbours) |-0.132 (-0.185, -0.080) |-0.105 (-0.153, -0.057) |-0.103 (-0.156, -0.050) | +|Rooms (neighbours) |0.072 (0.024, 0.121) |0.124 (0.078, 0.172) |0.040 (-0.006, 0.085) | +|Distance (neighbours) |-0.010 (-0.138, 0.120) |-0.025 (-0.149, 0.099) |-0.038 (-0.148, 0.072) | +|Pre-1940 (neighbours) |0.008 (-0.067, 0.082) |-0.024 (-0.096, 0.046) |-0.020 (-0.094, 0.055) | +|Social dissimilarity (b1) |— |0.592 (0.403, 0.797) |— | Model fit (the SDEM is a maximum-likelihood model and has no DIC, so only the two bml models are compared on fit): |Model | DIC| |:-------------------|---:| -|Equal weights (1/n) | 278| -|Similarity weights | 295| +|Equal weights (1/n) | 292| +|Social similarity | 229| + +The equal-weight bml model lines up closely with the SDEM, air quality at -0.105 against -0.112, crime at -0.061 against -0.067, rooms at 0.157 against 0.169. If equal-weight neighbour effects were all that was needed, errorsarlm would serve, and it would be faster. + +Estimating the aggregation changes the picture. +Letting weights fall with social dissimilarity lowers the DIC from 292 to 229, and the residual standard deviation from 0.128 to 0.114 which is an improvement in fit, not merely a reshuffling of the complexity penalty. +The member-level standard deviation is essentially unchanged (0.386 to 0.369), so the multiple-membership structure is still doing its work; what has changed is which neighbours it reads. -The own and neighbour coefficients from the equal-weight bml model line up with the SDEM, the reassurance worth having: the bml baseline reproduces the classical benchmark. -If equal-weight neighbour effects were all that was needed, `errorsarlm` would serve, and it would be faster. -What the parameterised weight function adds is the aggregation itself as something to estimate, which the SDEM cannot do. -On fit, the equal-weight model keeps the lower DIC (278 vs 295): the similarity structure in the weights is real — both dissimilarity coefficients sit away from zero — but the added flexibility does not pay for itself in overall fit here. The SDEM, a maximum-likelihood model, has no DIC, so it serves only as a coefficient yardstick. +The coefficients shift accordingly: the own-tract air-quality effect halves, from -0.105 to -0.059, while the neighbour rooms effect rises from 0.07 to 0.12. +This is the point of estimating the weights rather than imposing them. +Under equal weighting, part of what looks like an own-tract air-quality effect is really the influence of socially similar neighbours, misattributed because the aggregate averaged over every adjacent tract alike. +The departure from the SDEM here is expected: the benchmark can only ever apply the fixed weights it is given, whereas the weight function is exactly what this model estimates. [↑ Back to top](#top) @@ -582,7 +633,8 @@ monetPlot(mod_pm_n, "b.w.1", label = "b1: 1 = all weight on PM party, 0 = equal Posterior density and MCMC trace plot for the PM-party weight parameter b1 The posterior density is a single hump centred at about 0.23, just right of the zero line, with most of its mass to the right but a clear slice left of zero, so zero stays inside the 95% interval. -In the trace panel the three chains overlap and cover the same range with no separation — the visual counterpart to the clean diagnostics above. +The figures printed along the bottom axis are the 5th, 50th and 95th percentiles, so they sit inside the 95% interval reported above. +In the trace panel the three chains overlap and cover the same range with no separation, the visual counterpart to the clean diagnostics above. ## brms-style accessors @@ -619,38 +671,95 @@ posterior::summarise_draws(posterior::subset_draws(draws, variable = "b.w.1")) # 4. Does the average matter, or the spread? (emergent features) -The three examples so far all aggregated member attributes with a weighted *mean*: the additive micro-macro link. -But a coalition whose parties are uniformly moderate on some trait and a coalition mixing one extreme party with several counterbalancing ones can have the *same* mean. -If what breaks governments is internal contrast rather than the average level, the mean is silent about the mechanism. +The first three examples asked how friends, neighbours, and coalition partners combine. +This example and the two that follow ask the same of work: an occupation is not a single job but a bundle of activities, and the same activity turns up in many occupations at once. +The question is whether an occupation's fortunes follow the average of what its activities expose it to, or the shape of that bundle. + +## The data + +The data pair O*NET's occupation–activity structure with activity-level AI exposure scores and BLS wage series. +Each row is one activity's place in one occupation: 17,537 occupation–activity pairs spanning 894 occupations and 2,080 activities. +An occupation draws on about 20 activities, and an activity is used by about 8 occupations, up to 91 for the most general ones. + +| Variable | Level | Meaning | +|------------------|------------|-------------------------------------------------------| +| `occupation_id` | id | Occupation identifier (the group) | +| `member_id` | id | Work activity identifier (the member) | +| `wage_level` | outcome | Log median annual wage | +| `wage_growth` | outcome | Log change in median annual wage | +| `employment_growth` | outcome | Log change in occupational employment | +| `baseline_wage` | occupation | Median annual wage at the start of the growth window | +| `education` | occupation | Typical education needed for entry, ordinal 0–6 | +| `job_zone` | occupation | O*NET Job Zone, 1–5, a preparation-level band | +| `ai_exposure` | activity | The activity's exposure to AI (standardised) | +| `importance_raw` | activity | The activity's O*NET importance within the occupation | + +Unlike the coalition models, the weights here are observed rather than imposed. +`w(~ importance_raw, scale = TRUE)` lets each activity count in proportion to how important O*NET records it as being for that occupation, renormalised to sum to one within each occupation. +Because activities are shared across as many as 91 occupations, the activity random effects couple occupations that an occupation-level regression would treat as independent. + +## Mean and spread + +Two occupations can share an average exposure and be built very differently. +One spreads moderate exposure evenly across every activity it draws on; the other combines a few heavily exposed activities with several that are barely exposed at all. +If what matters is that contrast rather than the average level, the mean is silent about it. `fn()` selects the aggregation function a block applies to the weighted member records. -`fn("sum")` (used in every example so far) produces the weighted mean `A_finance`; `fn("var")` produces the weighted variance `V_finance`, a feature of the whole set. -Blocks stack, so both features can enter one model, each with its own main-model coefficient: +`fn("sum")` (used in every example so far) produces the weighted mean `A_ai_exposure`; `fn("var")` produces the weighted variance `V_ai_exposure`, a feature of the whole set. +Blocks stack, so both features can enter one model, each with its own main-model coefficient. +Occupation-level controls, required education and O*NET Job Zone are held throughout. + +The baseline is the mean on its own: ``` r -mod_var <- +mod_mean <- bml( - Surv(dur_wkb, event_wkb) ~ 1 + majority + mwc + + wage_level ~ 1 + education + job_zone + mm( - id = id(pid, gid), - vars = vars(finance), - w = w(~ 1/n, scale = TRUE), + id = id(member_id, occupation_id), + vars = vars(ai_exposure), + w = w(~ importance_raw, scale = TRUE), fn = fn("sum"), - RE = TRUE # mean block keeps the party random effects + RE = TRUE + ), + family = gaussian(), + iter = 2000, + warmup = 500, + chains = 3, + seed = 1, + data = onet + ) +``` + + + +Adding a second block, identical but for `fn("var")`, gives the spread its own coefficient: + + +``` r +mod_spread <- + bml( + wage_level ~ 1 + education + job_zone + + mm( + id = id(member_id, occupation_id), + vars = vars(ai_exposure), + w = w(~ importance_raw, scale = TRUE), + fn = fn("sum"), + RE = TRUE # mean block keeps the activity random effects ) + mm( - id = id(pid, gid), - vars = vars(finance), - w = w(~ 1/n, scale = TRUE), - fn = fn("var") # spread block: V_finance + id = id(member_id, occupation_id), + vars = vars(ai_exposure), + w = w(~ importance_raw, scale = TRUE), + fn = fn("var") # spread block: V_ai_exposure ), - family = weibull(), - iter = 5000, + family = gaussian(), + iter = 2000, warmup = 500, chains = 3, - seed = 1, - data = coalgov + seed = 1, + data = onet ) ``` @@ -659,29 +768,48 @@ mod_var <- ``` r bmlCompare( - "Mean only" = mod_eq, - "Mean + spread" = mod_var, - terms = c("A_finance", "V_finance"), - labels = c("Finance (mean)", "Finance (spread)") + "Mean only" = mod_mean, + "Mean + spread" = mod_spread, + terms = c("A_ai_exposure", "V_ai_exposure"), + labels = c("Exposure (mean)", "Exposure (spread)") ) ``` -|Term |Mean only |Mean + spread | -|:----------------|:-----------------------|:-----------------------| -|Finance (mean) |-0.313 (-0.562, -0.064) |-0.306 (-0.600, -0.006) | -|Finance (spread) |— |-0.015 (-0.268, 0.234) | -|N |628 |628 | -|DIC |3864 |3868 | +|Term |Mean only |Mean + spread | +|:-----------------|:--------------------|:-----------------------| +|Exposure (mean) |0.064 (0.007, 0.120) |0.070 (0.015, 0.126) | +|Exposure (spread) |— |-0.096 (-0.165, -0.027) | +|N |854 |854 | +|DIC |2633 |2316 | -The mean effect `A_finance` is essentially unchanged by adding the spread term (about -0.31 in both models), which is the first thing to check: the two features answer different questions rather than competing for the same variation. -The spread coefficient `V_finance` is about -0.02 with a 95% interval of [-0.27, 0.23] — a null. -For these data, what matters about a coalition's financial exposure is its level, not its internal contrast, and the DIC agrees (3864 for the mean-only model vs 3868 with the spread term). -A null on the spread is itself informative: it is the empirical license behind the mean-only models of Example 3, checked rather than assumed. + +``` r +loo_mean <- loo(mod_mean) +loo_spread <- loo(mod_spread) +loo::loo_compare(loo_mean, loo_spread) +``` + + + + +``` +#> model elpd_diff se_diff p_worse diag_diff diag_elpd +#> model2 0.0 0.0 NA 278 k_psis > 0.7 +#> model1 -2.2 3.3 0.74 |elpd_diff| < 4 297 k_psis > 0.7 +``` + +The mean effect `A_ai_exposure` is barely moved by adding the spread term (0.064 against 0.070), which is the first thing to check: the two features answer different questions rather than competing for the same variation. +The spread coefficient `V_ai_exposure` is -0.096 with a 95% interval of [-0.165, -0.027], comfortably clear of zero. +The fit statistics are less decisive than the coefficient itself. +DIC falls sharply, from 2633 to 2316, but the cross-validation comparison separates the two models by only 2.2 points of expected log predictive density against a standard error of 3.3, and flags the gap as too small to read. +Roughly a third of the observations exceed the Pareto-k threshold in that comparison, so loo is straining here; the coefficient carries this result rather than the model ranking. +The residual standard deviation holds at 0.134 and the activity-level standard deviation barely moves (0.844 to 0.837), so the spread term is not absorbing variation the mean block was already carrying. +Two occupations with the same average exposure are not equivalent: the one whose activities are unevenly exposed pays less. Both aggregation functions read the *same* weighted member records; only the reduction differs. -That is the sense in which the variance term is emergent: it is a property of the set that no single party's contribution can carry. +That is the sense in which the variance term is emergent: it is a property of the set that no single activity's contribution can carry. [↑ Back to top](#top) @@ -696,31 +824,32 @@ $$ $$ which runs from the *minimum* of the member attributes (as $\kappa \to -\infty$) through the weighted *mean* (at $\kappa \to 0$) to the *maximum* (as $\kappa \to +\infty$). -With `kappa = est()` the data choose the point on that path: is government survival driven by the coalition's average financial exposure, by its least exposed party (weakest link), or by its most exposed one? +With `kappa = est()` the data choose the point on that path: does an occupation's wage follow the average exposure of its activities, its most exposed activity, or its least exposed one? ``` r mod_smax <- bml( - Surv(dur_wkb, event_wkb) ~ 1 + majority + mwc + + wage_level ~ 1 + education + mm( - id = id(pid, gid), - vars = vars(finance), - w = w(~ 1/n, scale = TRUE), + id = id(member_id, occupation_id), + vars = vars(ai_exposure), + w = w(~ importance_raw, scale = TRUE), fn = fn("smax", kappa = est()) # kappa < 0: min-like, kappa > 0: max-like ) + mm( - id = id(pid, gid), - w = w(~ 1/n, scale = TRUE), + id = id(member_id, occupation_id), + w = w(~ importance_raw, scale = TRUE), fn = fn("sum"), - RE = TRUE # party random effects in their own block + RE = TRUE # activity random effects in their own block ), - family = weibull(), - iter = 5000, + family = gaussian(), + prior = prior(normal(0, 1), class = "fn"), + iter = 2000, warmup = 500, chains = 3, - seed = 1, - data = coalgov + seed = 1, + data = onet ) ``` @@ -728,26 +857,70 @@ mod_smax <- Shape parameters like `kappa` are not coefficients: they are reported in their own component (`fn[...]`), never in the coefficient table, and take priors through `prior(..., class = "fn")`. `get_prior()` lists what is settable. +The baseline column below is a mean-only fit carrying the same control as the smooth-max model, education alone. + + +``` r +bmlCompare( + "Mean" = mod_mean_edu, + "Smooth max" = mod_smax, + terms = c("A_ai_exposure", "smax_ai_exposure"), + labels = c("Exposure (mean)", "Exposure (smooth max)") +) +``` + + + +|Term |Mean |Smooth max | +|:---------------------|:--------------------|:--------------------| +|Exposure (mean) |0.095 (0.037, 0.153) |— | +|Exposure (smooth max) |— |0.126 (0.067, 0.184) | +|N |854 |854 | +|DIC |2643 |2528 | + + +``` r +loo_mean_edu <- loo(mod_mean_edu) +loo_smax <- loo(mod_smax) +loo::loo_compare(loo_mean_edu, loo_smax) +``` + + + + +``` +#> model elpd_diff se_diff p_worse diag_diff diag_elpd +#> model2 0.0 0.0 NA 298 k_psis > 0.7 +#> model1 -8.9 3.2 1.00 300 k_psis > 0.7 +``` + The posterior of `kappa` answers the question directly: ``` r -monetPlot(mod_smax, "fn.kappa.1", label = "kappa: < 0 weakest link, 0 mean, > 0 strongest link") +monetPlot(mod_smax, "fn.kappa.1", label = "kappa: < 0 least exposed activity, 0 mean, > 0 most exposed") ``` -Posterior density and MCMC trace for the smax shape parameter kappa +Posterior density and MCMC trace for the smax shape parameter kappa + +The feature's coefficient is 0.126 [0.067, 0.184], and `kappa` is -1.57 with a 95% interval of [-2.73, -0.40]. +Because `smax` collapses to the weighted mean at $\kappa = 0$, an interval clear of zero is a direct test rejecting the mean as the aggregation rule for these data, and 99.8% of the posterior mass sits below zero. +The negative sign places the aggregate toward the minimum: what tracks an occupation's wage is closer to its *least* exposed activity than to its average one. +Cross-validation agrees. +The smooth max beats the weighted mean by 8.9 points of expected log predictive density against a standard error of 3.2, better than two standard errors, and DIC moves the same way, from 2643 to 2528. +Here too about a third of the observations exceed the Pareto-k threshold, so the exact margin should be read loosely; the direction is not in doubt. -The feature's coefficient is clearly negative (`smax_finance` about -0.24 [-0.52, -0.06]), so the exposure feature matters — but the posterior for `kappa` is centred near 2.5 with a 95% interval of roughly [-6, 11], spanning weakest-link through strongest-link readings. -Two things follow. -First, the caveat that belongs with this estimand: `kappa` reaches the outcome only through the feature, so it is identified only insofar as differently-shaped aggregates fit differently — and with small coalitions (most have two to four parties), the min, mean, and max of finance are highly correlated, so the likelihood is nearly flat in `kappa`. -Second, a diffuse posterior here is not a failure of the sampler; it is the honest statement that these data cannot pin down the aggregation function, which is exactly what making `f` an estimand is for. +`kappa` reaches the outcome only through the feature, so it is identified only insofar as differently-shaped aggregates fit differently. +That is a real constraint: where members are few, the minimum, mean, and maximum of an attribute are nearly the same number and the likelihood is close to flat. +Occupations here draw on about 20 activities each, which is what gives the likelihood enough curvature to locate the shape. [↑ Back to top](#top) # 6. Does context change the aggregate's effect? (cross-level interaction) -Coalition majority status is a government-level (macro) attribute. -Because it is constant across the parties of a government, interacting it with the aggregated finance profile collapses to a product of two group-level quantities: a cross-level interaction. +The education an occupation typically requires is an occupation-level (macro) attribute. +Because it is constant across the activities that make up the occupation, interacting it with the aggregated exposure profile collapses to a product of two group-level quantities: a cross-level interaction. +A natural hypothesis is that required education buffers exposure, that AI exposure weighs less heavily on employment where more schooling is needed. In `bml` this is written by giving the block a `name` and referencing that name in the main formula. The macro variable must also appear as a main effect: @@ -756,21 +929,21 @@ The macro variable must also appear as a main effect: ``` r mod_xl <- bml( - Surv(dur_wkb, event_wkb) ~ 1 + majority + mwc + Afin:majority + + employment_growth ~ 1 + education + baseline_wage + Aexp:education + mm( - name = Afin, # named block: referencable feature - id = id(pid, gid), - vars = vars(finance), - w = w(~ 1/n, scale = TRUE), + name = Aexp, # named block: referencable feature + id = id(member_id, occupation_id), + vars = vars(ai_exposure), + w = w(~ importance_raw, scale = TRUE), fn = fn("sum"), RE = TRUE ), - family = weibull(), + family = gaussian(), iter = 5000, warmup = 500, chains = 3, - seed = 1, - data = coalgov + seed = 1, + data = onet_wg ) ``` @@ -778,21 +951,31 @@ mod_xl <- ``` r -fixef(mod_xl) +fixef(mod_xl)[c("Aexp", "Aexp:education"), ] ``` ``` -#> Estimate Est.Error Q2.5 Q97.5 -#> (Intercept) 6.8099757 0.1533719 6.51490635 7.1262164 -#> majority 0.2362897 0.1603613 -0.08576551 0.5481584 -#> mwc 0.4907232 0.1503586 0.19334021 0.7821277 -#> Afin -0.2281652 0.2198792 -0.65680427 0.2096516 -#> Afin:majority -0.1072704 0.2371408 -0.57573029 0.3637120 +#> Estimate Est.Error Q2.5 Q97.5 +#> Aexp 0.002506506 0.02279086 -0.04304832 0.045714862 +#> Aexp:education -0.009236561 0.00824449 -0.02499706 0.006929471 ``` -The interaction row `Afin:majority` asks whether the financial-exposure penalty differs between majority and minority governments. -The estimate is about -0.11 [-0.58, 0.36]: no evidence that majority status moderates the exposure penalty. -The name refers to the block's *fixed* feature only — the interaction multiplies `majority` with the aggregate `A_finance`, not with the party random effects, which continue to enter additively. +The `Aexp:education` interaction is -0.009 with a 95% interval of [-0.025, 0.007]. +The point estimate sits in the buffering direction, so occupations needing more schooling do show a slightly weaker exposure effect, but the interval covers zero and the estimate holds at that value across three seeds. +Required education does not measurably change how an occupation's exposure profile bears on its employment growth. + +A single moderator is thin evidence for a null, so the same specification was refit with Job Zone in place of education, O*NET's 1 to 5 banding of how much preparation an occupation needs. +That interaction comes back at 0.001 with a 95% interval of [-0.024, 0.024], centred on zero rather than merely overlapping it. +The two are close cousins as measures of required preparation, so their agreement is a weaker check than two independent moderators would be, but it is the check these data allow. + +Within these data, occupation-level context does not moderate the aggregate exposure effect. +The effect of an occupation's exposure profile is much the same whether the occupation demands a little schooling or a lot. +An ordinary occupation-level regression reaches the same conclusion, and here the multiple-membership model and the naive regression agree, which is worth noting precisely because they need not. +Example 2 showed the two diverging once the aggregation was estimated rather than imposed. + +The example still illustrates the mechanics that matter. +The `name` refers to the block's *fixed* feature only: the interaction multiplies the macro variable with the aggregate `A_ai_exposure`, not with the activity random effects, which continue to enter additively. +A well-specified cross-level interaction that returns a null is itself informative — it says the aggregate's effect is stable across this dimension of context, rather than leaving the question unasked. [↑ Back to top](#top) diff --git a/vignettes/examples.Rmd.orig b/vignettes/examples.Rmd.orig index 7103cd7..948f3d1 100644 --- a/vignettes/examples.Rmd.orig +++ b/vignettes/examples.Rmd.orig @@ -34,7 +34,7 @@ knitr::opts_chunk$set( This vignette presents seven examples that show how the multiple-membership multilevel model (MMMM) can be applied to network, spatial, and aggregation problems. The first three cover the additive micro-macro link: peer influence in a friendship network, neighbourhood effects on home values, and how coalition parties shape government survival. -The last four continue the coalition example to showcase aggregation beyond the weighted mean: distributional features (`fn("var")`), estimating the aggregation function itself (`fn("smax", kappa = est())`), cross-level interactions through named blocks, and heterogeneous member effects (`re(1 + x)`). +The last four go beyond the weighted mean to showcase distributional features (`fn("var")`), estimating the aggregation function itself (`fn("smax", kappa = est())`), cross-level interactions through named blocks, and heterogeneous member effects (`re(1 + x)`). # 1. All friends, or just your best friend? (network regression) @@ -62,50 +62,53 @@ The model uses the following variables: | `rank` | weight | Friendship closeness (1 = best friend, 5 = fifth-closest) | | `n` | weight | Number of friends nominated by the ego | -Two models are fit and compared. -The economics standard for this kind of network problem is the linear-in-means model: every friend carries equal weight, so what matters is the average across all of them, implying that each friend exerts the same influence. -The alternative is that influence is concentrated: a single close friend drives most of the effect while looser ties matter little. -`bml` fits both and compares them directly, changing only one line of the weight function. +Three models are fit and compared, tracing a spectrum from the most concentrated weighting rule to the most diffuse. +At one end, influence is concentrated: a single close friend drives the effect while looser ties matter little. +At the other is the economics standard for this kind of network problem, the linear-in-means model: every friend carries equal weight, so what matters is the average across all of them. +Between the two sits a graded rule in which closeness matters but no friend is discarded. +All three are fixed weighting rules, nothing inside `w()` is estimated, so the comparison isolates the shape of the aggregation, changing only one line of the weight function each time. -## Model 1: linear-in-means +## Model 1: best friend only -Every friend receives equal weight, `w ~ 1/n`, so the peer term is the simple average of all nominated friends' Wave 1 smoking. +Only the closest-ranked friend contributes; all others are dropped. +`rank == min(rank)` evaluates to `TRUE` for the closest-ranked friend within each ego's friend set (rank 1), placing all weight on that friend. +With `scale = TRUE` the weight is renormalised to 1 within each ego, so the peer term becomes that single friend's smoking. ```{r, eval = FALSE} -mod_lim <- +mod_bf <- bml( - smoking_w2 ~ smoking_ego + gender_ego + # outcome ~ ego-level predictors - mm( # multiple-membership term (the peers) - id = id(alter, ego), # member = alter (friend), group = ego - vars = vars(smoking_alter), # member variable to aggregate - w = w(~ 1/n, scale = TRUE), # equal weights, renormalised to sum to 1 - fn = fn("sum"), # aggregation function: the weighted mean - RE = TRUE # member-level random effect + smoking_w2 ~ smoking_ego + gender_ego + + mm( + id = id(alter, ego), + vars = vars(smoking_alter), + w = w(~ rank == min(rank), scale = TRUE), # all weight on the closest friend + fn = fn("sum"), # aggregation function: the weighted mean + RE = TRUE ), - family = gaussian(), # continuous outcome + family = gaussian(), data = cils_long ) ``` ```{r, echo = FALSE} -mod_lim <- .load("mod_lim") +mod_bf <- .load("mod_bf") ``` -## Model 2: best friend only +## Model 2: graded decay -Only the closest-ranked friend contributes; all others are dropped. -`rank == min(rank)` evaluates to `TRUE` for the closest-ranked friend within each ego's friend set (rank 1), placing all weight on that friend. -With `scale = TRUE` the weight is renormalised to 1 within each ego, so the peer term becomes that single friend's smoking. +The best-friend rule is an extreme as it throws away four of the five nominations. +A softer version keeps every friend but lets closeness matter (`w ~ 1/rank`) +gives the best friend weight 1, the second-closest 1/2, and so on down to 1/5. ```{r, eval = FALSE} -mod_bf <- +mod_decay <- bml( smoking_w2 ~ smoking_ego + gender_ego + mm( id = id(alter, ego), vars = vars(smoking_alter), - w = w(~ rank == min(rank), scale = TRUE), # all weight on the closest friend - fn = fn("sum"), # same aggregation function as Model 1 + w = w(~ 1/rank, scale = TRUE), # closeness-graded weights + fn = fn("sum"), RE = TRUE ), family = gaussian(), @@ -114,30 +117,60 @@ mod_bf <- ``` ```{r, echo = FALSE} -mod_bf <- .load("mod_bf") +mod_decay <- .load("mod_decay") ``` -## Comparing the two models +## Model 3: linear-in-means + +Every friend receives equal weight, `w ~ 1/n`, so the peer term is the simple average of all nominated friends' Wave 1 smoking. + +```{r, eval = FALSE} +mod_lim <- + bml( + smoking_w2 ~ smoking_ego + gender_ego + # outcome ~ ego-level predictors + mm( # multiple-membership term (the peers) + id = id(alter, ego), # member = alter (friend), group = ego + vars = vars(smoking_alter), # member variable to aggregate + w = w(~ 1/n, scale = TRUE), # equal weights, renormalised to sum to 1 + fn = fn("sum"), # aggregation function: the weighted mean + RE = TRUE # member-level random effect + ), + family = gaussian(), # continuous outcome + data = cils_long + ) +``` + +```{r, echo = FALSE} +mod_lim <- .load("mod_lim") +``` + +## Comparing the three models `bmlCompare()` lays out the coefficient estimates and fit statistics of several models side by side, one column per model, with readable labels supplied through `labels`. ```{r} bmlCompare( - "Linear-in-means" = mod_lim, "Best friend only" = mod_bf, + "Graded decay" = mod_decay, + "Linear-in-means" = mod_lim, terms = c("(Intercept)", "smoking_ego", "gender_ego", "A_smoking_alter"), labels = c("Intercept", "Own smoking (W1)", "Gender", "Peers' smoking (W1)") ) ``` -Both models find a positive peer effect: students whose friends smoke more tend to smoke more themselves one year later, over and above their own Wave 1 smoking. -The own-smoking coefficient (about 0.36 in both models) dominates, which is expected since past behaviour is the strongest predictor of future behaviour. +All three models find a positive peer effect: students whose friends smoke more tend to smoke more themselves one year later, over and above their own Wave 1 smoking. +The own-smoking coefficient (about 0.36 across the board) dominates, which is expected since past behaviour is the strongest predictor of future behaviour. + +The peer effect itself depends on how the friend set is aggregated. +Restricting attention to the best friend gives about 0.03 with an interval that touches zero ([-0.00, 0.05]). +Both of the rules that keep every friend give about 0.04 and clear zero (graded decay [0.00, 0.08], linear-in-means [0.00, 0.08]). +Discarding the looser ties does not sharpen the estimate; it loses signal. -The peer effect is present but modest. -Under linear-in-means, averaging across all nominated friends gives a coefficient of about 0.04 (CI [0.00, 0.08]), which just clears zero. -Limiting to the best friend gives about 0.03 (CI [-0.00, 0.05]): a smaller point estimate with a tighter interval that now just touches zero. +The fit statistics fall monotonically as the weighting becomes more diffuse: DIC drops from 18,700 for the best-friend rule to 17,577 for graded decay and 16,964 for linear-in-means. +Each relaxation of the concentration assumption buys a substantial improvement, and the largest single gain comes from not throwing the other four friends away. -On fit, the linear-in-means model is clearly preferred (DIC 17,029 vs 18,463, a gap of over 1,400). +The member-level standard deviation tells the same story from another angle, rising from 0.40 under best-friend weighting to 0.83 and then 0.93. +Under the best-friend rule most nominated friends receive zero weight in every group they belong to, so the data carry little information about their individual effects and the estimated spread shrinks toward the prior. Peer influence on smoking appears to be spread across the friend set rather than concentrated in the single closest tie. ## Cross-validation: which aggregation does the data prefer? @@ -146,31 +179,34 @@ DIC is convenient but crude; PSIS-LOO cross-validation compares the models on th For gaussian outcomes `bml` monitors the pointwise log-likelihood in the generated JAGS model, so `loo()` works directly on the fitted objects, and `loo::loo_compare()` ranks them: ```{r, eval = FALSE} -loo_lim <- loo(mod_lim) -loo_bf <- loo(mod_bf) -loo::loo_compare(loo_lim, loo_bf) +loo_bf <- loo(mod_bf) +loo_decay <- loo(mod_decay) +loo_lim <- loo(mod_lim) +loo::loo_compare(loo_bf, loo_decay, loo_lim) ``` ```{r, echo = FALSE} -loo_lim <- readRDS(file.path(".fits", "loo_lim.rds")) -loo_bf <- readRDS(file.path(".fits", "loo_bf.rds")) +loo_bf <- readRDS(file.path(".fits", "loo_bf.rds")) +loo_decay <- readRDS(file.path(".fits", "loo_decay.rds")) +loo_lim <- readRDS(file.path(".fits", "loo_lim.rds")) ``` ```{r, echo = FALSE} -loo::loo_compare(loo_lim, loo_bf) +loo::loo_compare(loo_bf, loo_decay, loo_lim) ``` -The linear-in-means model comes out on top: the best-friend model loses about 37 points of expected log predictive density (standard error about 10), agreeing with the DIC ranking. -Both criteria point the same way — influence is spread over the whole friend set rather than concentrated in the closest tie. +model1 is best-friend, model2 graded decay, model3 linear-in-means, in the order passed. +PSIS-LOO ranks them as DIC does: the best-friend model loses about 37 points of expected log predictive density (standard error 10, p_worse 1.00), while graded decay sits only 6 points behind linear-in-means (standard error 5, p_worse 0.89), decisive at the bottom of the ladder, close at the top. +The Pareto-k flags bias elpd_diff toward whichever model carries more of them (8 for linear-in-means, 3 for graded decay), so that narrow margin is best read as an upper bound. +Both criteria point the same way: influence is spread over the friend set rather than concentrated in the closest tie. -A posterior predictive check makes sure the preferred model reproduces the shape of the outcome at all. `pp_check()` overlays replicated outcome distributions on the observed one: ```{r, eval = FALSE} pp_check(mod_lim) ``` -```{r, echo = FALSE, fig.width = 6, fig.height = 3.5, fig.alt = "Posterior predictive check: replicated outcome densities overlaid on the observed smoking distribution"} +```{r ppc-lim, echo = FALSE, fig.width = 6, fig.height = 3.5, fig.alt = "Posterior predictive check: replicated outcome densities overlaid on the observed smoking distribution"} ppc <- readRDS(file.path(".fits", "ppc_lim.rds")) bayesplot::ppc_dens_overlay(ppc$y, ppc$yrep) ``` @@ -202,7 +238,7 @@ library(bml) boston <- read_sf(system.file("shapes/boston_tracts.gpkg", package = "spData")) |> - select(CMEDV, NOX, CRIM, RM, DIS, AGE, geom) |> + select(CMEDV, NOX, CRIM, RM, DIS, AGE, LSTAT, geom) |> st_transform(crs = 5070) |> mutate( tid = row_number(), @@ -229,13 +265,9 @@ boston_df <- select(-CMEDV, -lnCMEDV, -geom) |> rename_with(~ paste0(.x, "_nb")), by = c("tid_nb" = "tid_nb") - ) - -boston_df2 <- boston_df |> - mutate( - d_CRIM = as.numeric(scale(log1p(abs(CRIM - CRIM_nb)))), # log first: crime is right-skewed - d_AGE = as.numeric(scale(abs(AGE - AGE_nb))) - ) + ) |> + # aggregation variable only: LSTAT enters no formula, just the weight function + mutate(d_LSTAT = as.numeric(scale(abs(LSTAT - LSTAT_nb)))) ``` The data are the 1970 Boston Standard Metropolitan Statistical Area census tracts shipped with `spData`. @@ -272,6 +304,7 @@ mod_bml <- RE = TRUE # neighbour-level random effect ), family = gaussian(), + iter = 5000, warmup = 500, chains = 3, seed = 1, data = boston_df ) ``` @@ -280,87 +313,96 @@ mod_bml <- mod_bml <- .load("mod_bml") ``` -## Parameterised weights: similarity across covariates +## Parameterised weights: social similarity -Equal weighting assumes every neighbour matters the same regardless of how similar it is to the focal tract. -A natural alternative is that tracts are more strongly influenced by neighbours that resemble them, a spatial homophily effect. -Rather than collapse similarity into a single number, the weight is allowed to depend on the dissimilarity in each covariate separately, using the functional form from Rosche (2026): +Equal weighting assumes every adjacent tract matters the same, regardless of how +much it resembles the focal tract. +A natural alternative is that spillovers run more strongly between tracts that +are alike socially, social distance layered on top of physical adjacency. +One parameter is enough to test it: $$ -w_{ij} = \frac{1}{1 + (n_i - 1)\exp\!\bigl(-(b_0 + b_1 \cdot d^{\text{CRIM}}_{ij} + b_2 \cdot d^{\text{AGE}}_{ij})\bigr)} +w_{ij} \propto \exp\!\bigl(-b_1 \cdot d^{\text{LSTAT}}_{ij}\bigr) $$ -Each `d` is the standardised absolute difference between a focal tract and its neighbour on that covariate. -When all the `b`'s are zero this collapses to `1/n_i`, the equal-weight baseline. -A negative coefficient means neighbours that are similar on that covariate (low `d`) receive more weight; a positive coefficient means dissimilar neighbours dominate. -Letting the two dissimilarities enter on their own lets the model decide which kind of similarity drives the aggregation, instead of imposing a single combined measure. +where $d^{\text{LSTAT}}_{ij}$ is the standardised absolute difference between +tract $i$ and neighbour $j$ in the share of lower-status population. +At $b_1 = 0$ the weights collapse to $1/n_i$ and the equal-weight baseline is +recovered exactly, so the baseline is nested inside this model. +Positive $b_1$ means socially similar neighbours carry more weight. +`LSTAT` appears nowhere else in the model — not among the own-tract covariates, +not in `vars()`, so it acts purely as an aggregation variable. ```{r, eval = FALSE} -mod_bml_w <- +mod_bml_lstat <- bml( lnCMEDV ~ NOX + CRIM + RM + DIS + AGE + mm( id = id(tid_nb, tid), vars = vars(NOX_nb + CRIM_nb + RM_nb + DIS_nb + AGE_nb), - # weight as a function of covariate dissimilarity; nests 1/n when all b's are 0 - w = w(~ 1 / (1 + (n - 1) * exp(-(b0 + b1 * d_CRIM + b2 * d_AGE))), scale = TRUE), + # weight falls with social dissimilarity; nests 1/n at b1 = 0 + w = w(~ exp(-b1 * d_LSTAT), scale = TRUE), fn = fn("sum"), RE = TRUE ), family = gaussian(), - prior = prior(normal(0, 1), class = "w"), # weakly informative; prevents exp() overflow - data = boston_df2 + prior = prior(normal(0, 1), class = "w"), # weakly informative; keeps exp() stable + iter = 5000, warmup = 500, chains = 3, seed = 1, + data = boston_df ) ``` ```{r, echo = FALSE} -mod_bml_w <- .load("mod_bml_w") +mod_bml_lstat <- .load("mod_bml_lstat") ``` -The estimated weight-function parameters are: - ```{r} bmlCompare( - "Similarity weights" = mod_bml_w, + "Social similarity" = mod_bml_lstat, component = "weights", - terms = c("w[b0] (mm.1)", "w[b1] (mm.1)", "w[b2] (mm.1)"), - labels = c("Weight intercept (b0)", "Crime dissimilarity (b1)", "Age dissimilarity (b2)") + terms = c("w[d_LSTAT] (mm.1)"), + labels = c("Social dissimilarity (b1)") ) ``` -The crime-dissimilarity coefficient is about −1.76 (CI [−2.59, −1.16]), clearly negative, so neighbours with similar crime levels carry more weight. -The age-dissimilarity coefficient is about 0.71 (CI [0.07, 1.52]), pointing the other way: neighbours with different age profiles carry more weight. -The model picks up two kinds of similarity at once, in opposite directions. +The estimate is about 0.59 with a 95% interval of [0.40, 0.80], comfortably +away from zero. +Because `d_LSTAT` is standardised, a one-standard-deviation increase in social +dissimilarity multiplies a neighbour's unnormalised weight by +$\exp(-0.59) \approx 0.55$, roughly halving its influence. +Adjacency alone does not determine who counts: among neighbours that all share +a border, the socially similar ones dominate the aggregate. ## Do the weights actually vary with similarity? -A direct check plots the estimated weight each neighbour received against the dissimilarity between that neighbour and its focal tract. -If crime homophily drives the weights, low-`d_CRIM` pairs should cluster toward high weights and high-`d_CRIM` pairs toward low weights. +A direct check plots the estimated weight each neighbour received against its +social dissimilarity from the focal tract. +If social homophily drives the weights, low-`d_LSTAT` pairs should sit toward +high weights and high-`d_LSTAT` pairs toward low weights. -```{r, fig.width = 6, fig.height = 4, echo = FALSE, fig.alt = "Estimated neighbour weight plotted against covariate dissimilarity, faceted by crime and age dissimilarity"} -w_mat <- mod_bml_w$w[[1]] # 506 x 15 matrix, NA padded +```{r weight-dissim, fig.width = 6, fig.height = 4, echo = FALSE, fig.alt = "Estimated neighbour weight plotted against social dissimilarity"} +w_mat <- mod_bml_lstat$w[[1]] # tracts x max-neighbours, NA padded -weight_df <- boston_df2 |> +weight_df <- boston_df |> group_by(tid) |> mutate(pos = row_number()) |> ungroup() |> mutate(weight = w_mat[cbind(tid, pos)]) |> filter(!is.na(weight)) -weight_long <- weight_df |> - tidyr::pivot_longer(c(d_CRIM, d_AGE), names_to = "covariate", values_to = "dissimilarity") - -ggplot(weight_long, aes(x = dissimilarity, y = weight)) + +ggplot(weight_df, aes(x = d_LSTAT, y = weight)) + geom_point(alpha = 0.15, size = 0.8) + - geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = FALSE, color = "#2c7bb6") + - facet_wrap(~ covariate, scales = "free_x") + - labs(x = "Dissimilarity", y = "Estimated neighbour weight") + + geom_smooth(method = "loess", se = FALSE, color = "#2c7bb6") + + labs(x = "Social dissimilarity (d_LSTAT)", y = "Estimated neighbour weight") + theme_minimal() ``` -The two panels show the similarity effects acting in opposite directions. -The weight falls as crime dissimilarity rises (the homophily the negative crime coefficient pointed to), while it rises gently with age dissimilarity, matching the positive age coefficient. -Each panel holds only a partial view, since every neighbour's weight reflects both dissimilarities at once, but the opposing slopes are visible. +The estimated weight declines with social dissimilarity: neighbours close to the +focal tract on lower-status share carry the most weight, and the weight falls +away as dissimilarity grows. +The scatter around the trend is renormalisation at work, a pair's final weight +depends on how many neighbours it competes with within its tract, but the +downward slope, the homophily the positive `b1` pointed to, is unmistakable. ## Benchmark: spatial Durbin error model @@ -395,8 +437,8 @@ relabel <- function(t) { # Two bml models: fixed effects + weight parameters tab_bml <- bind_rows( - tidy(mod_bml, component = "all") |> mutate(model = "Equal weights (1/n)"), - tidy(mod_bml_w, component = "all") |> mutate(model = "Similarity weights") + tidy(mod_bml, component = "all") |> mutate(model = "Equal weights (1/n)"), + tidy(mod_bml_lstat, component = "all") |> mutate(model = "Social similarity") ) |> filter(component %in% c("fixed", "weights")) |> transmute(model, term = relabel(term), estimate, conf.low, conf.high) @@ -419,7 +461,7 @@ tab_sdem <- tibble::tibble( term_levels <- c("Intercept", "NOX", "CRIM", "RM", "DIS", "AGE", "NOX (nb)", "CRIM (nb)", "RM (nb)", "DIS (nb)", "AGE (nb)", - "w[b0]", "w[b1]", "w[b2]") + "w[d_LSTAT]") label_map <- c( "Intercept" = "Intercept", @@ -433,9 +475,7 @@ label_map <- c( "RM (nb)" = "Rooms (neighbours)", "DIS (nb)" = "Distance (neighbours)", "AGE (nb)" = "Pre-1940 (neighbours)", - "w[b0]" = "Weight intercept (b0)", - "w[b1]" = "Crime dissimilarity (b1)", - "w[b2]" = "Age dissimilarity (b2)" + "w[d_LSTAT]" = "Social dissimilarity (b1)" ) bind_rows(tab_bml, tab_sdem) |> @@ -453,16 +493,22 @@ Model fit (the SDEM is a maximum-likelihood model and has no DIC, so only the tw ```{r, echo = FALSE} tibble::tibble( - Model = c("Equal weights (1/n)", "Similarity weights"), - DIC = c(glance(mod_bml)$DIC, glance(mod_bml_w)$DIC) + Model = c("Equal weights (1/n)", "Social similarity"), + DIC = c(glance(mod_bml)$DIC, glance(mod_bml_lstat)$DIC) ) |> knitr::kable(digits = 0) ``` -The own and neighbour coefficients from the equal-weight bml model line up with the SDEM, the reassurance worth having: the bml baseline reproduces the classical benchmark. -If equal-weight neighbour effects were all that was needed, `errorsarlm` would serve, and it would be faster. -What the parameterised weight function adds is the aggregation itself as something to estimate, which the SDEM cannot do. -On fit, the equal-weight model keeps the lower DIC (278 vs 295): the similarity structure in the weights is real — both dissimilarity coefficients sit away from zero — but the added flexibility does not pay for itself in overall fit here. The SDEM, a maximum-likelihood model, has no DIC, so it serves only as a coefficient yardstick. +The equal-weight bml model lines up closely with the SDEM, air quality at -0.105 against -0.112, crime at -0.061 against -0.067, rooms at 0.157 against 0.169. If equal-weight neighbour effects were all that was needed, errorsarlm would serve, and it would be faster. + +Estimating the aggregation changes the picture. +Letting weights fall with social dissimilarity lowers the DIC from 292 to 229, and the residual standard deviation from 0.128 to 0.114 which is an improvement in fit, not merely a reshuffling of the complexity penalty. +The member-level standard deviation is essentially unchanged (0.386 to 0.369), so the multiple-membership structure is still doing its work; what has changed is which neighbours it reads. + +The coefficients shift accordingly: the own-tract air-quality effect halves, from -0.105 to -0.059, while the neighbour rooms effect rises from 0.07 to 0.12. +This is the point of estimating the weights rather than imposing them. +Under equal weighting, part of what looks like an own-tract air-quality effect is really the influence of socially similar neighbours, misattributed because the aggregate averaged over every adjacent tract alike. +The departure from the SDEM here is expected: the benchmark can only ever apply the fixed weights it is given, whereas the weight function is exactly what this model estimates. [↑ Back to top](#top) @@ -670,7 +716,8 @@ monetPlot(mod_pm_n, "b.w.1", label = "b1: 1 = all weight on PM party, 0 = equal ``` The posterior density is a single hump centred at about 0.23, just right of the zero line, with most of its mass to the right but a clear slice left of zero, so zero stays inside the 95% interval. -In the trace panel the three chains overlap and cover the same range with no separation — the visual counterpart to the clean diagnostics above. +The figures printed along the bottom axis are the 5th, 50th and 95th percentiles, so they sit inside the 95% interval reported above. +In the trace panel the three chains overlap and cover the same range with no separation, the visual counterpart to the clean diagnostics above. ## brms-style accessors @@ -690,60 +737,136 @@ posterior::summarise_draws(posterior::subset_draws(draws, variable = "b.w.1")) # 4. Does the average matter, or the spread? (emergent features) -The three examples so far all aggregated member attributes with a weighted *mean*: the additive micro-macro link. -But a coalition whose parties are uniformly moderate on some trait and a coalition mixing one extreme party with several counterbalancing ones can have the *same* mean. -If what breaks governments is internal contrast rather than the average level, the mean is silent about the mechanism. +The first three examples asked how friends, neighbours, and coalition partners combine. +This example and the two that follow ask the same of work: an occupation is not a single job but a bundle of activities, and the same activity turns up in many occupations at once. +The question is whether an occupation's fortunes follow the average of what its activities expose it to, or the shape of that bundle. + +## The data + +The data pair O*NET's occupation–activity structure with activity-level AI exposure scores and BLS wage series. +Each row is one activity's place in one occupation: 17,537 occupation–activity pairs spanning 894 occupations and 2,080 activities. +An occupation draws on about 20 activities, and an activity is used by about 8 occupations, up to 91 for the most general ones. + +| Variable | Level | Meaning | +|------------------|------------|-------------------------------------------------------| +| `occupation_id` | id | Occupation identifier (the group) | +| `member_id` | id | Work activity identifier (the member) | +| `wage_level` | outcome | Log median annual wage | +| `wage_growth` | outcome | Log change in median annual wage | +| `employment_growth` | outcome | Log change in occupational employment | +| `baseline_wage` | occupation | Median annual wage at the start of the growth window | +| `education` | occupation | Typical education needed for entry, ordinal 0–6 | +| `job_zone` | occupation | O*NET Job Zone, 1–5, a preparation-level band | +| `ai_exposure` | activity | The activity's exposure to AI (standardised) | +| `importance_raw` | activity | The activity's O*NET importance within the occupation | + +Unlike the coalition models, the weights here are observed rather than imposed. +`w(~ importance_raw, scale = TRUE)` lets each activity count in proportion to how important O*NET records it as being for that occupation, renormalised to sum to one within each occupation. +Because activities are shared across as many as 91 occupations, the activity random effects couple occupations that an occupation-level regression would treat as independent. + +## Mean and spread + +Two occupations can share an average exposure and be built very differently. +One spreads moderate exposure evenly across every activity it draws on; the other combines a few heavily exposed activities with several that are barely exposed at all. +If what matters is that contrast rather than the average level, the mean is silent about it. `fn()` selects the aggregation function a block applies to the weighted member records. -`fn("sum")` (used in every example so far) produces the weighted mean `A_finance`; `fn("var")` produces the weighted variance `V_finance`, a feature of the whole set. -Blocks stack, so both features can enter one model, each with its own main-model coefficient: +`fn("sum")` (used in every example so far) produces the weighted mean `A_ai_exposure`; `fn("var")` produces the weighted variance `V_ai_exposure`, a feature of the whole set. +Blocks stack, so both features can enter one model, each with its own main-model coefficient. +Occupation-level controls, required education and O*NET Job Zone are held throughout. + +The baseline is the mean on its own: ```{r, eval = FALSE} -mod_var <- +mod_mean <- bml( - Surv(dur_wkb, event_wkb) ~ 1 + majority + mwc + + wage_level ~ 1 + education + job_zone + mm( - id = id(pid, gid), - vars = vars(finance), - w = w(~ 1/n, scale = TRUE), + id = id(member_id, occupation_id), + vars = vars(ai_exposure), + w = w(~ importance_raw, scale = TRUE), + fn = fn("sum"), + RE = TRUE + ), + family = gaussian(), + iter = 2000, + warmup = 500, + chains = 3, + seed = 1, + data = onet + ) +``` + +```{r, echo = FALSE} +mod_mean <- .load("mod_mean") +``` + +Adding a second block, identical but for `fn("var")`, gives the spread its own coefficient: + +```{r, eval = FALSE} +mod_spread <- + bml( + wage_level ~ 1 + education + job_zone + + mm( + id = id(member_id, occupation_id), + vars = vars(ai_exposure), + w = w(~ importance_raw, scale = TRUE), fn = fn("sum"), - RE = TRUE # mean block keeps the party random effects + RE = TRUE # mean block keeps the activity random effects ) + mm( - id = id(pid, gid), - vars = vars(finance), - w = w(~ 1/n, scale = TRUE), - fn = fn("var") # spread block: V_finance + id = id(member_id, occupation_id), + vars = vars(ai_exposure), + w = w(~ importance_raw, scale = TRUE), + fn = fn("var") # spread block: V_ai_exposure ), - family = weibull(), - iter = 5000, + family = gaussian(), + iter = 2000, warmup = 500, chains = 3, - seed = 1, - data = coalgov + seed = 1, + data = onet ) ``` ```{r, echo = FALSE} -mod_var <- .load("mod_var") +mod_spread <- .load("mod_spread") ``` ```{r} bmlCompare( - "Mean only" = mod_eq, - "Mean + spread" = mod_var, - terms = c("A_finance", "V_finance"), - labels = c("Finance (mean)", "Finance (spread)") + "Mean only" = mod_mean, + "Mean + spread" = mod_spread, + terms = c("A_ai_exposure", "V_ai_exposure"), + labels = c("Exposure (mean)", "Exposure (spread)") ) ``` -The mean effect `A_finance` is essentially unchanged by adding the spread term (about -0.31 in both models), which is the first thing to check: the two features answer different questions rather than competing for the same variation. -The spread coefficient `V_finance` is about -0.02 with a 95% interval of [-0.27, 0.23] — a null. -For these data, what matters about a coalition's financial exposure is its level, not its internal contrast, and the DIC agrees (3864 for the mean-only model vs 3868 with the spread term). -A null on the spread is itself informative: it is the empirical license behind the mean-only models of Example 3, checked rather than assumed. +```{r, eval = FALSE} +loo_mean <- loo(mod_mean) +loo_spread <- loo(mod_spread) +loo::loo_compare(loo_mean, loo_spread) +``` + +```{r, echo = FALSE} +loo_mean <- readRDS(file.path(".fits", "loo_mean.rds")) +loo_spread <- readRDS(file.path(".fits", "loo_spread.rds")) +``` + +```{r, echo = FALSE} +loo::loo_compare(loo_mean, loo_spread) +``` + +The mean effect `A_ai_exposure` is barely moved by adding the spread term (0.064 against 0.070), which is the first thing to check: the two features answer different questions rather than competing for the same variation. +The spread coefficient `V_ai_exposure` is -0.096 with a 95% interval of [-0.165, -0.027], comfortably clear of zero. +The fit statistics are less decisive than the coefficient itself. +DIC falls sharply, from 2633 to 2316, but the cross-validation comparison separates the two models by only 2.2 points of expected log predictive density against a standard error of 3.3, and flags the gap as too small to read. +Roughly a third of the observations exceed the Pareto-k threshold in that comparison, so loo is straining here; the coefficient carries this result rather than the model ranking. +The residual standard deviation holds at 0.134 and the activity-level standard deviation barely moves (0.844 to 0.837), so the spread term is not absorbing variation the mean block was already carrying. +Two occupations with the same average exposure are not equivalent: the one whose activities are unevenly exposed pays less. Both aggregation functions read the *same* weighted member records; only the reduction differs. -That is the sense in which the variance term is emergent: it is a property of the set that no single party's contribution can carry. +That is the sense in which the variance term is emergent: it is a property of the set that no single activity's contribution can carry. [↑ Back to top](#top) @@ -758,56 +881,91 @@ $$ $$ which runs from the *minimum* of the member attributes (as $\kappa \to -\infty$) through the weighted *mean* (at $\kappa \to 0$) to the *maximum* (as $\kappa \to +\infty$). -With `kappa = est()` the data choose the point on that path: is government survival driven by the coalition's average financial exposure, by its least exposed party (weakest link), or by its most exposed one? +With `kappa = est()` the data choose the point on that path: does an occupation's wage follow the average exposure of its activities, its most exposed activity, or its least exposed one? ```{r, eval = FALSE} mod_smax <- bml( - Surv(dur_wkb, event_wkb) ~ 1 + majority + mwc + + wage_level ~ 1 + education + mm( - id = id(pid, gid), - vars = vars(finance), - w = w(~ 1/n, scale = TRUE), + id = id(member_id, occupation_id), + vars = vars(ai_exposure), + w = w(~ importance_raw, scale = TRUE), fn = fn("smax", kappa = est()) # kappa < 0: min-like, kappa > 0: max-like ) + mm( - id = id(pid, gid), - w = w(~ 1/n, scale = TRUE), + id = id(member_id, occupation_id), + w = w(~ importance_raw, scale = TRUE), fn = fn("sum"), - RE = TRUE # party random effects in their own block + RE = TRUE # activity random effects in their own block ), - family = weibull(), - iter = 5000, + family = gaussian(), + prior = prior(normal(0, 1), class = "fn"), + iter = 2000, warmup = 500, chains = 3, - seed = 1, - data = coalgov + seed = 1, + data = onet ) ``` ```{r, echo = FALSE} -mod_smax <- .load("mod_smax") +mod_mean_edu <- .load("mod_mean_edu") +mod_smax <- .load("mod_smax_onet") ``` Shape parameters like `kappa` are not coefficients: they are reported in their own component (`fn[...]`), never in the coefficient table, and take priors through `prior(..., class = "fn")`. `get_prior()` lists what is settable. +The baseline column below is a mean-only fit carrying the same control as the smooth-max model, education alone. + +```{r} +bmlCompare( + "Mean" = mod_mean_edu, + "Smooth max" = mod_smax, + terms = c("A_ai_exposure", "smax_ai_exposure"), + labels = c("Exposure (mean)", "Exposure (smooth max)") +) +``` + +```{r, eval = FALSE} +loo_mean_edu <- loo(mod_mean_edu) +loo_smax <- loo(mod_smax) +loo::loo_compare(loo_mean_edu, loo_smax) +``` + +```{r, echo = FALSE} +loo_mean_edu <- readRDS(file.path(".fits", "loo_mean_edu.rds")) +loo_smax <- readRDS(file.path(".fits", "loo_smax_onet.rds")) +``` + +```{r, echo = FALSE} +loo::loo_compare(loo_mean_edu, loo_smax) +``` + The posterior of `kappa` answers the question directly: -```{r, fig.alt = "Posterior density and MCMC trace for the smax shape parameter kappa"} -monetPlot(mod_smax, "fn.kappa.1", label = "kappa: < 0 weakest link, 0 mean, > 0 strongest link") +```{r monet-kappa, fig.alt = "Posterior density and MCMC trace for the smax shape parameter kappa"} +monetPlot(mod_smax, "fn.kappa.1", label = "kappa: < 0 least exposed activity, 0 mean, > 0 most exposed") ``` -The feature's coefficient is clearly negative (`smax_finance` about -0.24 [-0.52, -0.06]), so the exposure feature matters — but the posterior for `kappa` is centred near 2.5 with a 95% interval of roughly [-6, 11], spanning weakest-link through strongest-link readings. -Two things follow. -First, the caveat that belongs with this estimand: `kappa` reaches the outcome only through the feature, so it is identified only insofar as differently-shaped aggregates fit differently — and with small coalitions (most have two to four parties), the min, mean, and max of finance are highly correlated, so the likelihood is nearly flat in `kappa`. -Second, a diffuse posterior here is not a failure of the sampler; it is the honest statement that these data cannot pin down the aggregation function, which is exactly what making `f` an estimand is for. +The feature's coefficient is 0.126 [0.067, 0.184], and `kappa` is -1.57 with a 95% interval of [-2.73, -0.40]. +Because `smax` collapses to the weighted mean at $\kappa = 0$, an interval clear of zero is a direct test rejecting the mean as the aggregation rule for these data, and 99.8% of the posterior mass sits below zero. +The negative sign places the aggregate toward the minimum: what tracks an occupation's wage is closer to its *least* exposed activity than to its average one. +Cross-validation agrees. +The smooth max beats the weighted mean by 8.9 points of expected log predictive density against a standard error of 3.2, better than two standard errors, and DIC moves the same way, from 2643 to 2528. +Here too about a third of the observations exceed the Pareto-k threshold, so the exact margin should be read loosely; the direction is not in doubt. + +`kappa` reaches the outcome only through the feature, so it is identified only insofar as differently-shaped aggregates fit differently. +That is a real constraint: where members are few, the minimum, mean, and maximum of an attribute are nearly the same number and the likelihood is close to flat. +Occupations here draw on about 20 activities each, which is what gives the likelihood enough curvature to locate the shape. [↑ Back to top](#top) # 6. Does context change the aggregate's effect? (cross-level interaction) -Coalition majority status is a government-level (macro) attribute. -Because it is constant across the parties of a government, interacting it with the aggregated finance profile collapses to a product of two group-level quantities: a cross-level interaction. +The education an occupation typically requires is an occupation-level (macro) attribute. +Because it is constant across the activities that make up the occupation, interacting it with the aggregated exposure profile collapses to a product of two group-level quantities: a cross-level interaction. +A natural hypothesis is that required education buffers exposure, that AI exposure weighs less heavily on employment where more schooling is needed. In `bml` this is written by giving the block a `name` and referencing that name in the main formula. The macro variable must also appear as a main effect: @@ -815,35 +973,48 @@ The macro variable must also appear as a main effect: ```{r, eval = FALSE} mod_xl <- bml( - Surv(dur_wkb, event_wkb) ~ 1 + majority + mwc + Afin:majority + + employment_growth ~ 1 + education + baseline_wage + Aexp:education + mm( - name = Afin, # named block: referencable feature - id = id(pid, gid), - vars = vars(finance), - w = w(~ 1/n, scale = TRUE), + name = Aexp, # named block: referencable feature + id = id(member_id, occupation_id), + vars = vars(ai_exposure), + w = w(~ importance_raw, scale = TRUE), fn = fn("sum"), RE = TRUE ), - family = weibull(), + family = gaussian(), iter = 5000, warmup = 500, chains = 3, - seed = 1, - data = coalgov + seed = 1, + data = onet_wg ) ``` ```{r, echo = FALSE} -mod_xl <- .load("mod_xl") +mod_xl <- .load("mod_xl_onet") ``` ```{r} -fixef(mod_xl) +fixef(mod_xl)[c("Aexp", "Aexp:education"), ] ``` -The interaction row `Afin:majority` asks whether the financial-exposure penalty differs between majority and minority governments. -The estimate is about -0.11 [-0.58, 0.36]: no evidence that majority status moderates the exposure penalty. -The name refers to the block's *fixed* feature only — the interaction multiplies `majority` with the aggregate `A_finance`, not with the party random effects, which continue to enter additively. +The `Aexp:education` interaction is -0.009 with a 95% interval of [-0.025, 0.007]. +The point estimate sits in the buffering direction, so occupations needing more schooling do show a slightly weaker exposure effect, but the interval covers zero and the estimate holds at that value across three seeds. +Required education does not measurably change how an occupation's exposure profile bears on its employment growth. + +A single moderator is thin evidence for a null, so the same specification was refit with Job Zone in place of education, O*NET's 1 to 5 banding of how much preparation an occupation needs. +That interaction comes back at 0.001 with a 95% interval of [-0.024, 0.024], centred on zero rather than merely overlapping it. +The two are close cousins as measures of required preparation, so their agreement is a weaker check than two independent moderators would be, but it is the check these data allow. + +Within these data, occupation-level context does not moderate the aggregate exposure effect. +The effect of an occupation's exposure profile is much the same whether the occupation demands a little schooling or a lot. +An ordinary occupation-level regression reaches the same conclusion, and here the multiple-membership model and the naive regression agree, which is worth noting precisely because they need not. +Example 2 showed the two diverging once the aggregation was estimated rather than imposed. + +The example still illustrates the mechanics that matter. +The `name` refers to the block's *fixed* feature only: the interaction multiplies the macro variable with the aggregate `A_ai_exposure`, not with the activity random effects, which continue to enter additively. +A well-specified cross-level interaction that returns a null is itself informative — it says the aggregate's effect is stable across this dimension of context, rather than leaving the question unasked. [↑ Back to top](#top) diff --git a/vignettes/fit_models.R b/vignettes/fit_models.R index 4e8de2f..b6845c8 100644 --- a/vignettes/fit_models.R +++ b/vignettes/fit_models.R @@ -7,9 +7,14 @@ # examples.Rmd.orig then just readRDS() them, so precompile.R runs in seconds. # # Re-run this script only when a model's specification, data, or seed changes. -# Requires JAGS + rjags and the CILS4EU feathers used by Example 1 (not shipped): -# point at them via BML_CILS_DIR, e.g. -# Sys.setenv(BML_CILS_DIR = "C:/path/to/cils/data") +# Requires JAGS + rjags, plus two datasets that do not ship with the package and +# should live OUTSIDE this repository: +# * the CILS4EU feathers used by Example 1 +# * the "_dwa" parquets written by build_onet_ai_dataset.py, used by Examples 4-6 +# Both paths are required env vars (no defaults, by design -- see the notes at +# each read). Set them once in ~/.Renviron: +# BML_CILS_DIR=/abs/path/to/cils +# ONET_DATA_DIR=/abs/path/to/data/processed_data # # Run from the repository root: Rscript vignettes/fit_models.R # The .fits/ directory is git-ignored; the committed artdefact is examples.Rmd. @@ -61,8 +66,17 @@ save_fit <- function(obj, name) saveRDS(slim_fit(obj), file.path(fits_dir, paste # the arrow native library can segfault when loaded in the same session as # sf/GEOS. Delete cils_raw.rds (or the whole .fits/ dir) to force a re-read. # ----------------------------------------------------------------------------- -cils_dir <- Sys.getenv("BML_CILS_DIR", unset = "data") + +# No default: the old fallback was "data", which resolves to the package's own +# data/ directory when this is run from the repo root. That is a tracked, shipped +# directory and restricted microdata must never land there, so require the env var. +cils_dir <- Sys.getenv("BML_CILS_DIR", unset = "") +if (!nzchar(cils_dir)) { + stop("BML_CILS_DIR is not set. Point it at the folder holding nodedat.feather, ", + "edgedat.feather and nodedat-w2.feather (outside this repository).") +} cils_rds <- file.path(fits_dir, "cils_raw.rds") + if (!file.exists(cils_rds)) { conv <- tempfile(fileext = ".R") writeLines(c( @@ -79,6 +93,7 @@ if (!file.exists(cils_rds)) { stop("Failed to read CILS feathers in an isolated process. Check BML_CILS_DIR.") } } + cils <- readRDS(cils_rds) nodedat <- cils$nodedat edgedat <- cils$edgedat @@ -112,7 +127,7 @@ mod_lim <- ), family = gaussian(), iter = 5000, - warmup = 1000, + warmup = 500, chains = 3, seed = 1, data = cils_long @@ -131,18 +146,40 @@ mod_bf <- ), family = gaussian(), iter = 5000, - warmup = 1000, + warmup = 500, chains = 3, seed = 1, data = cils_long ) save_fit(mod_bf, "mod_bf") +# Third weighting rule: closeness-graded, keeps every friend (vignette Model 2). +mod_decay <- + bml( + smoking_w2 ~ smoking_ego + gender_ego + + mm( + id = id(alter, ego), + vars = bml::vars(smoking_alter), + w = w(~ 1/rank, scale = TRUE), + fn = fn("sum"), + RE = TRUE + ), + family = gaussian(), + iter = 5000, + warmup = 500, + chains = 3, + seed = 1, + data = cils_long + ) +save_fit(mod_decay, "mod_decay") + # Cross-validation and posterior-predictive artifacts for the vignette. These need # the full draws, which slim_fit strips, so compute them here and cache the small # results: loo objects (~n x 3 pointwise) and a 30-draw yrep matrix. -saveRDS(loo(mod_lim), file.path(fits_dir, "loo_lim.rds")) -saveRDS(loo(mod_bf), file.path(fits_dir, "loo_bf.rds")) +# (save_fit() does not mutate its argument, so the in-memory fits are still full.) +saveRDS(loo(mod_lim), file.path(fits_dir, "loo_lim.rds")) +saveRDS(loo(mod_bf), file.path(fits_dir, "loo_bf.rds")) +saveRDS(loo(mod_decay), file.path(fits_dir, "loo_decay.rds")) saveRDS( list(y = mod_lim$input$y, yrep = posterior_predict(mod_lim, ndraws = 30)), file.path(fits_dir, "ppc_lim.rds") @@ -150,10 +187,14 @@ saveRDS( # ----------------------------------------------------------------------------- # Example 2: Boston tracts (spData; reconstructed live in the vignette too) +# +# NOTE: examples.Rmd.orig rebuilds this exact frame in a hidden chunk. Any change +# here must be mirrored there or the two drift apart. # ----------------------------------------------------------------------------- + boston <- read_sf(system.file("shapes/boston_tracts.gpkg", package = "spData")) |> - select(CMEDV, NOX, CRIM, RM, DIS, AGE, geom) |> + select(CMEDV, NOX, CRIM, RM, DIS, AGE, LSTAT, geom) |> st_transform(crs = 5070) |> mutate( tid = row_number(), @@ -179,13 +220,9 @@ boston_df <- select(-CMEDV, -lnCMEDV, -geom) |> rename_with(~ paste0(.x, "_nb")), by = c("tid_nb" = "tid_nb") - ) - -boston_df2 <- boston_df |> - mutate( - d_CRIM = as.numeric(scale(log1p(abs(CRIM - CRIM_nb)))), - d_AGE = as.numeric(scale(abs(AGE - AGE_nb))) - ) + ) |> + # aggregation variable only: LSTAT enters no formula, just the weight function + mutate(d_LSTAT = as.numeric(scale(abs(LSTAT - LSTAT_nb)))) mod_bml <- bml( @@ -201,18 +238,18 @@ mod_bml <- iter = 5000, warmup = 500, chains = 3, - seed = 42, + seed = 1, data = boston_df ) save_fit(mod_bml, "mod_bml") -mod_bml_w <- +mod_bml_lstat <- bml( lnCMEDV ~ NOX + CRIM + RM + DIS + AGE + mm( id = id(tid_nb, tid), vars = bml::vars(NOX_nb + CRIM_nb + RM_nb + DIS_nb + AGE_nb), - w = w(~ 1 / (1 + (n - 1) * exp(-(b0 + b1 * d_CRIM + b2 * d_AGE))), scale = TRUE), + w = w(~ exp(-b1 * d_LSTAT), scale = TRUE), # nests 1/n at b1 = 0 fn = fn("sum"), RE = TRUE ), @@ -221,16 +258,18 @@ mod_bml_w <- iter = 5000, warmup = 500, chains = 3, - seed = 42, - data = boston_df2 + seed = 1, + data = boston_df ) -save_fit(mod_bml_w, "mod_bml_w") +save_fit(mod_bml_lstat, "mod_bml_lstat") # ----------------------------------------------------------------------------- # Example 3: coalition governments (coalgov ships with the package) # ----------------------------------------------------------------------------- + data(coalgov) coalgov$prime <- as.integer(coalgov$prime) + coalgov <- coalgov |> group_by(gid) |> mutate(is_smallest = as.integer(pseat == min(pseat) & pid == pid[which.min(pseat)])) |> @@ -299,89 +338,242 @@ mod_pm_p <- save_fit(mod_pm_p, "mod_pm_p") # ----------------------------------------------------------------------------- -# Example 4: mean vs spread (sum + var features, stacked) +# Examples 4-6: O*NET occupation-activity structure (data not shipped) +# +# Written by build_onet_ai_dataset.py into its CONFIG$OUT_DIR. We use the "_dwa" +# build: members are detailed work activities, which recur across occupations, so +# the multiple-membership structure is present. The "_task" build puts each task +# in exactly one occupation -- no multiple membership, and mm() would add nothing +# over an ordinary occupation-level regression. +# +# TWO tables are needed. The long table carries the occupation-activity structure +# (importance weights, member exposure); the occupation-level outcomes and +# moderators live only in the model table. They are joined below. +# +# Both are read in an ISOLATED R process and cached as onet_raw.rds, for the same +# reason as the CILS feathers above: the arrow native library can segfault when +# loaded alongside sf/GEOS, and this script loads both. +# Delete onet_raw.rds (or the whole .fits/ dir) to force a re-read. # ----------------------------------------------------------------------------- -mod_var <- + +onet_dir <- Sys.getenv("ONET_DATA_DIR", unset = "") +if (!nzchar(onet_dir)) { + stop("ONET_DATA_DIR is not set. Point it at the folder build_onet_ai_dataset.py ", + "writes to (its CONFIG$OUT_DIR, e.g. /abs/path/to/data/processed_data).") +} + +onet_long_file <- "long_occupation_member_dwa.parquet" +onet_model_file <- "model_occupation_year_dwa.parquet" +onet_rds <- file.path(fits_dir, "onet_raw.rds") + +if (!file.exists(onet_rds)) { + conv <- tempfile(fileext = ".R") + writeLines(c( + sprintf("onet_dir <- %s", deparse(onet_dir)), + "saveRDS(list(", + sprintf(" long = arrow::read_parquet(file.path(onet_dir, %s)),", deparse(onet_long_file)), + sprintf(" model = arrow::read_parquet(file.path(onet_dir, %s))", deparse(onet_model_file)), + sprintf("), %s)", deparse(onet_rds)) + ), conv) + rscript <- file.path(R.home("bin"), + if (.Platform$OS.type == "windows") "Rscript.exe" else "Rscript") + if (system2(rscript, shQuote(conv)) != 0) { + stop("Failed to read O*NET parquets in an isolated process. Check ONET_DATA_DIR ", + "and that both _dwa parquets are present.") + } +} + +onet_parts <- readRDS(onet_rds) + +# The model table is one row per occupation, so this is a clean 1-to-many join +# onto the long rows. Only the columns the vignette models are carried over. +onet_raw <- + onet_parts$long |> + inner_join( + onet_parts$model |> + select(occupation_id, wage_level, employment_growth, + education, job_zone, baseline_wage), + by = "occupation_id" + ) + +# bml's id() needs the member and group identifiers. If it rejects the character +# O*NET codes ("15-1252.00", "4.A.2.a.4.i.3"), uncomment these two lines -- +# relabelling identifiers does not change the model, only how they are indexed. +# onet_raw$occupation_id <- as.integer(factor(onet_raw$occupation_id)) +# onet_raw$member_id <- as.integer(factor(onet_raw$member_id)) + +# Fail loudly and specifically if the column names differ from what the vignette +# prose documents, rather than erroring somewhere inside bml(). +.onet_req <- c("occupation_id", "member_id", "wage_level", "employment_growth", + "education", "job_zone", "baseline_wage", + "ai_exposure", "importance_raw") +.onet_missing <- setdiff(.onet_req, names(onet_raw)) +if (length(.onet_missing)) { + stop("O*NET frame is missing: ", paste(.onet_missing, collapse = ", "), + ". Check onet_file and the column names against onet_runs.Rmd.") +} + +# Two complete-case frames, because Examples 4-5 and Example 6 sit on different +# outcomes. +# +# The vignette's w(~ importance_raw, scale = TRUE) reproduces the pipeline's own +# `w` column exactly: build_onet_ai_dataset.py drops members lacking exposure and +# then renormalises, and both filters below drop whole occupations rather than +# individual members, so the within-occupation normalisation is unchanged. +# +# onet -- wage_level. ONE frame for both Example 4 and Example 5: Example 5 +# controls on education alone, but it must be estimated on the same rows +# as Example 4 or the two sections report different N and the tables stop +# being comparable. Target: N = 854 occupations. +# onet_wg -- employment_growth, for Example 6. Complete-cased on exactly the +# variables mod_xl_onet uses, so the frame matches the confirmed run. +onet <- onet_raw |> + filter(!is.na(wage_level), !is.na(education), !is.na(job_zone), + !is.na(ai_exposure), !is.na(importance_raw)) + +onet_wg <- onet_raw |> + filter(!is.na(employment_growth), !is.na(education), !is.na(baseline_wage), + !is.na(ai_exposure), !is.na(importance_raw)) + +# Example 4 baseline: weighted mean only, education + job_zone +mod_mean <- bml( - Surv(dur_wkb, event_wkb) ~ 1 + majority + mwc + + wage_level ~ 1 + education + job_zone + mm( - id = id(pid, gid), - vars = vars(finance), - w = w(~ 1/n, scale = TRUE), + id = id(member_id, occupation_id), + vars = bml::vars(ai_exposure), + w = w(~ importance_raw, scale = TRUE), + fn = fn("sum"), + RE = TRUE + ), + family = gaussian(), + iter = 2000, + warmup = 500, + chains = 3, + seed = 1, + data = onet + ) +save_fit(mod_mean, "mod_mean") + +# Example 4: mean + spread (sum + var features, stacked) +mod_spread <- + bml( + wage_level ~ 1 + education + job_zone + + mm( + id = id(member_id, occupation_id), + vars = bml::vars(ai_exposure), + w = w(~ importance_raw, scale = TRUE), fn = fn("sum"), RE = TRUE ) + mm( - id = id(pid, gid), - vars = vars(finance), - w = w(~ 1/n, scale = TRUE), + id = id(member_id, occupation_id), + vars = bml::vars(ai_exposure), + w = w(~ importance_raw, scale = TRUE), fn = fn("var") ), - family = weibull(), - monitor = TRUE, - iter = 5000, + family = gaussian(), + iter = 2000, warmup = 500, chains = 3, seed = 1, - data = coalgov + data = onet ) -save_fit(mod_var, "mod_var") +save_fit(mod_spread, "mod_spread") + +# Example 5 baseline: weighted mean only, education alone -- matches the controls +# of mod_smax_onet, so the two columns of the Example 5 table are comparable. +mod_mean_edu <- + bml( + wage_level ~ 1 + education + + mm( + id = id(member_id, occupation_id), + vars = bml::vars(ai_exposure), + w = w(~ importance_raw, scale = TRUE), + fn = fn("sum"), + RE = TRUE + ), + family = gaussian(), + iter = 2000, + warmup = 500, + chains = 3, + seed = 1, + data = onet + ) +save_fit(mod_mean_edu, "mod_mean_edu") -# ----------------------------------------------------------------------------- # Example 5: what is the aggregation function? (smax with estimated kappa) -# ----------------------------------------------------------------------------- -mod_smax <- +# monitor = TRUE is required: the vignette calls monetPlot() on fn.kappa.1. +mod_smax_onet <- bml( - Surv(dur_wkb, event_wkb) ~ 1 + majority + mwc + + wage_level ~ 1 + education + mm( - id = id(pid, gid), - vars = vars(finance), - w = w(~ 1/n, scale = TRUE), + id = id(member_id, occupation_id), + vars = bml::vars(ai_exposure), + w = w(~ importance_raw, scale = TRUE), fn = fn("smax", kappa = est()) ) + mm( - id = id(pid, gid), - w = w(~ 1/n, scale = TRUE), + id = id(member_id, occupation_id), + w = w(~ importance_raw, scale = TRUE), fn = fn("sum"), RE = TRUE ), - family = weibull(), + family = gaussian(), + prior = prior(normal(0, 1), class = "fn"), monitor = TRUE, - iter = 5000, + iter = 2000, warmup = 500, chains = 3, seed = 1, - data = coalgov + data = onet ) -save_fit(mod_smax, "mod_smax") +save_fit(mod_smax_onet, "mod_smax_onet") -# ----------------------------------------------------------------------------- -# Example 6: cross-level interaction via a named block -# ----------------------------------------------------------------------------- -mod_xl <- +# Example 6: cross-level interaction via a named block, on employment growth. +# Kept at 5000/500 to match the three-seed confirmation the pasted numbers come +# from. monitor = TRUE because the vignette calls fixef() on this fit. +mod_xl_onet <- bml( - Surv(dur_wkb, event_wkb) ~ 1 + majority + mwc + Afin:majority + + employment_growth ~ 1 + education + baseline_wage + Aexp:education + mm( - name = Afin, - id = id(pid, gid), - vars = vars(finance), - w = w(~ 1/n, scale = TRUE), + name = Aexp, + id = id(member_id, occupation_id), + vars = bml::vars(ai_exposure), + w = w(~ importance_raw, scale = TRUE), fn = fn("sum"), RE = TRUE ), - family = weibull(), + family = gaussian(), monitor = TRUE, iter = 5000, warmup = 500, chains = 3, seed = 1, - data = coalgov + data = onet_wg ) -save_fit(mod_xl, "mod_xl") +save_fit(mod_xl_onet, "mod_xl_onet") + +# loo objects for the Example 4 and 5 comparisons. As above, these need the full +# draws, so they are computed here from the un-slimmed in-memory fits. +saveRDS(loo(mod_mean), file.path(fits_dir, "loo_mean.rds")) +saveRDS(loo(mod_spread), file.path(fits_dir, "loo_spread.rds")) +saveRDS(loo(mod_mean_edu), file.path(fits_dir, "loo_mean_edu.rds")) +saveRDS(loo(mod_smax_onet), file.path(fits_dir, "loo_smax_onet.rds")) + +# slim_fit()'s keep-list is rownames(reg.table) + "deviance". kappa is reported in +# the fn[...] component, which may not be in reg.table, so confirm it survived the +# round-trip: monetPlot() in Example 5 reads it back off the cached object. +.m <- readRDS(file.path(fits_dir, "mod_smax_onet.rds")) +if (!"fn.kappa.1" %in% dimnames(.m$jags.out$BUGSoutput$sims.array)[[3]]) { + stop("fn.kappa.1 was dropped by slim_fit(); extend its keep-list before precompiling.") +} +rm(.m) # ----------------------------------------------------------------------------- # Example 7: heterogeneous member effects (explained + residual in one block) # ----------------------------------------------------------------------------- + mod_het <- bml( Surv(dur_wkb, event_wkb) ~ 1 + majority + mwc +