Discrete Choice Model for Affordable Housing Applications

This notebook implements a Bayesian discrete choice model designed to understand how applicants from different demographic groups make decisions about applying to affordable housing opportunities. The model was developed as part of a collaboration with CHAPA, with the goal of providing a transparent, data-driven way to study equity and access in lottery application patterns.

Rather than simply predicting who applies, the model aims to quantify how different applicant characteristics and property attributes influence the probability of applying, and how these relationships differ across racial and ethnic groups.

---

Notebook Structure

This notebook is organized to balance clarity, transparency, and reproducibility:

• Data Cleaning: Preparation of applicant- and property-level datasets, variable construction, and data diagnostics.
• Final Model (Model 2): The primary Bayesian discrete choice model used for inference. This section includes prior evaluation, convergence diagnostics, posterior predictive checks, and substantive findings.
• Appendix (Models 1A & 1B): Earlier model iterations that document the development process. These versions illustrate how prior choices, model structure, and diagnostic results informed the final specification.

Readers interested in the full modeling workflow may start with data cleaning, review the Appendix models, and then return to the final model for the completed analysis.

---

What the Model Does

At a high level, the model represents each applicant’s decision to apply (or not apply) as a function of:

• Applicant characteristics

  • Age
  • Disability status
  • Number of dependents
  • Household income

• Property characteristics

  • Maximum resale price
  • Distance from the applicant’s home to the property

• Demographic group

  • The model allows each racial or ethnic group to have its own set of coefficients, capturing potential differences in how various factors shape application decisions.

These inputs are combined to form a utility score, which is then passed through a logistic (sigmoid) function to obtain the probability of application for each applicant–property pair.

The model is fully probabilistic: it does not output single estimates, but posterior distributions, giving us uncertainty intervals around all effects.

---

Why This Model Works

Discrete choice models are widely used in economics and policy analysis because they mirror the way real decisions are made. The intuition is simple:

1. Each applicant receives a “utility” from each property, based on its features and their own characteristics.
2. Utility is modeled as a linear combination of effects (e.g., how income interacts with resale price).
3. Random variation is captured through applicant-level uncertainty.
4. The higher the utility, the more likely the applicant is to apply.

This Bayesian formulation offers several advantages:

• Hierarchical race-specific effects: In Model 2, Each racial/ethnic group gets its own set of coefficients, allowing us to examine heterogeneity without overfitting.

• Partial pooling: The model shrinks extreme subgroup estimates toward the overall mean, improving robustness when some groups have fewer observations.

• Uncertainty quantification: Instead of a single point estimate, the model provides a distribution for each effect.

• Applicant-level variability: A random effect (if computationally feasible) absorbs unobserved applicant differences (e.g., motivation, access to information), reducing bias in the fixed effects.

---

What Questions This Model Can Answer

This model can support a wide range of fairness, equity, and program-design questions, including:

1. How do different factors influence the probability of applying?

For example:

• Does household income make applicants more or less likely to apply to higher-priced properties?
• How does distance affect application likelihood?
• How much does disability status increase or reduce the probability of applying?

2. Do these effects differ across racial or ethnic groups?

The hierarchical structure allows us to estimate:

• Which factors matter most for each group
• Where differences between groups are credible vs. driven by noise
• Whether systemic barriers may be impacting certain groups disproportionately

3. What is the probability that a specific type of applicant would apply?

You can simulate:

• A high-income applicant vs. a low-income applicant
• Applicants with vs. without dependents
• Older vs. younger applicants

4. How sensitive are application rates to changes in policy or outreach?

Because the model produces a utility function, you can run counterfactuals like:

• How would applications change if properties had lower resale prices?
• If transportation access improved and distance effectively decreased?
• If outreach increased awareness among groups with lower baseline application probabilities?

5. Can we identify where access gaps exist?

By comparing group-specific coefficients, we can highlight where certain groups face structural disadvantages in the application process.

# import packages

# Data Wrangling
import numpy as np
import pandas as pd
import os

# Bayesian Methods
import pymc as pm
import arviz as az

# visualization
import matplotlib.pyplot as plt
az.style.use("arviz-darkgrid")

# update working directory
os.chdir("..")

# import data
dt = pd.read_parquet("data/processed/applications_clean.parquet")

Creating our Discrete Choice Dataset

# some quick data cleaning

# if 'disabled'/'yes' -> 1; 0 otherwise.
dt['Disability'] = dt['Disability'].isin(["Disabled", "Yes"]).astype(int)

# if NA -> 0
dt['Dependents'] = dt['Dependents'].fillna(0)

# if HH Income missing, drop
dt = dt.dropna(subset = ['HH Income'])

# Impute Age conditional on Income and number of Dependents
def impute_age(df):

    dt['IncomeBin'] = pd.qcut(
        dt['HH Income'],
        q = [0, 0.2, 0.4, 0.6, 0.8, 1],
        labels = ["Q1", "Q2", "Q3", "Q4", "Q5"]
)
    
    # compute group medians
    group_medians = (
        df.groupby(['IncomeBin', 'Dependents'])['Age']
        .median()
        .rename('Age_median')
    )
    
    # merge medians back to df
    df = df.merge(group_medians, on=['IncomeBin', 'Dependents'], how='left')
    
    # impute
    df['Age'] = df['Age'].fillna(df['Age_median'])
    
    # drop helper column
    return df.drop(columns=['Age_median'])

dt = impute_age(dt)

def collapse_race(raw):
    """
    Per current OMB rules.
    """
    r = raw.lower()

    # Hispanic supersedes
    if "hispanic" in r:
        return "Hispanic/Latino"

    # Unknown / prefer not to answer
    if "unknown" in r or "choose" in r or "not" in r:
        return "Unknown"

    # Single-race categories
    if r == "white" or r == "white_mena":
        return "White"
    if "black_african_american" in r and "/" not in r:
        return "Black"
    if r == "asian":
        return "Asian"
    if r == "native_american_alaskan_native":
        return "Native American"

    # Multi-race (non-Hispanic)
    if "/" in r:
        return "Two or More Races"

    # Fallback
    return "Unknown"

dt["race_norm_final"] = dt["race_norm_final"].apply(collapse_race)

# rename applicant current residence lat/lon
dt = dt.rename({'longitude': 'applicant_longitude', 'latitude': 'applicant_latitude'}, axis=1)

/var/folders/q9/vssjztfs6gd47qy2w2t2f0t80000gn/T/ipykernel_96020/3501486212.py:23: FutureWarning: The default of observed=False is deprecated and will be changed to True in a future version of pandas. Pass observed=False to retain current behavior or observed=True to adopt the future default and silence this warning.
  df.groupby(['IncomeBin', 'Dependents'])['Age']

# unique applicants
app_dt = (dt
          .groupby("ID Number")[['Age', 'Disability', 
                                 'HH Size', 'Dependents', 
                                 'race_norm_final', 
                                 'HH Income', 'HH Assets', 
                                 'applicant_MTSP_Income_Limit',
                                 'applicant_longitude', 'applicant_latitude']].last()).reset_index()

# unique properties
prop_dt = dt[['Application Property', 'Property Maximum Resale Price', 'property_latitude', 'property_longitude']].drop_duplicates().dropna()

# creating cartesian product dataset
choice_dt = pd.merge(app_dt, prop_dt, how='cross')

# create 'applied' binary outcome variable
applied_dt = pd.concat(
    (dt[["ID Number", "Application Property"]], 
     pd.Series( np.ones(len(dt)), name="applied")
     ), axis=1)

# merge using cart-product
choice_dt = choice_dt.merge(applied_dt, on=["ID Number", "Application Property"], how="left")

# find distance in miles
choice_dt["distance_miles"] = 3958.8 * 2 * np.arcsin(
    np.sqrt(
        np.sin(np.radians(choice_dt["property_latitude"] - choice_dt["applicant_latitude"]) / 2) ** 2
        + np.cos(np.radians(choice_dt["applicant_latitude"]))
        * np.cos(np.radians(choice_dt["property_latitude"]))
        * np.sin(np.radians(choice_dt["property_longitude"] - choice_dt["applicant_longitude"]) / 2) ** 2
    )
)

choice_dt["applied"] = choice_dt["applied"].fillna(0)
choice_dt = choice_dt.drop(['applicant_latitude', 'applicant_longitude', 'property_latitude', 'property_longitude'], axis=1)

# TO-DO: NEED BETTER SOLUTION HERE
choice_dt = choice_dt.dropna(subset="distance_miles")

# High household size taming (is this a good fallback value??)
choice_dt.loc[(choice_dt['HH Size']>8), 'HH Size'] = 0
# Alternative: -- accounts for 60% of answers correctly.
# choice_dt.loc[(choice_dt['HH Size']>8), 'HH Size'] = choice_dt.loc[(choice_dt['HH Size']>8), 'Dependents'] + 1

choice_dt.describe()

	ID Number	Age	Disability	HH Size	Dependents	HH Income	HH Assets	applicant_MTSP_Income_Limit	Property Maximum Resale Price	applied	distance_miles
count	138781.000000	138781.000000	138781.000000	138781.000000	138781.000000	1.387810e+05	137161.000000	127760.000000	138781.000000	138781.000000	138781.000000
mean	507439.637386	41.090084	0.066983	2.033578	0.548332	1.203145e+05	45495.000912	60314.898638	240674.258609	0.012343	36.010267
std	289296.522630	13.649327	0.249994	1.202886	0.860182	1.956038e+06	59531.148035	12271.125513	46420.441183	0.110413	21.543833
min	101.000000	0.000000	0.000000	0.000000	0.000000	7.000000e+01	0.000000	40050.000000	148438.000000	0.000000	0.000000
25%	247829.000000	31.000000	0.000000	1.000000	0.000000	5.200000e+04	14000.000000	52550.000000	204948.000000	0.000000	20.584763
50%	521641.000000	38.000000	0.000000	2.000000	0.000000	6.400000e+04	30733.890000	57900.000000	235200.000000	0.000000	32.617685
75%	750684.000000	49.000000	0.000000	3.000000	1.000000	7.700000e+04	59915.430000	66200.000000	273400.000000	0.000000	46.724262
max	998585.000000	86.000000	1.000000	7.000000	7.000000	7.020040e+07	863967.300000	160900.000000	374062.000000	1.000000	168.023139

$$ \text{Final Model Description (Model 2)} $$

---

Hyper-Priors $$ \text{Fixed Effects Mean Hyper-prior: } \mu_\beta \sim \mathcal{N}([0, 0.15, 0.15, 0, -0.2, -1.4],0.1) $$ $$ \text{Fixed Effects Variance Hyper-prior: } \sigma_\beta \sim \text{Exp}(1) $$

$$ \text{Fixed Effects Priors: } \beta_\text{race} \sim \mathcal{N}(\mu_\beta, \sigma_\beta) $$$$ \text{Applicant Slope Mean Prior: } \mu_{\text{app}} \sim \mathcal{N}(-4,1) $$$$ \text{Applicant Slope Variance Prior: } \sigma_{\text{app}}^2 \sim \text{Half-Normal}(0.3) $$

---

Intercept

$$ \gamma \sim \mathcal{N}(\mu_{\text{app}}, \sigma_{\text{app}}^2) $$

---

Utility

$$ \eta_{ij} = \gamma + \mathbf{X_{ij}^T}\beta_\text{race} $$$$ p_{ij} = \frac{1}{1+\text{exp}(-\eta_{ij})} $$$$ Y_{ij} \sim \text{Bernoulli}(p_{ij}) $$

features = ['Age', 'Disability', 'Dependents', 'HH Income', 'Property Maximum Resale Price', 'distance_miles']
X_FE = choice_dt[features]
fe_prior_mu = np.array([0, 0.15, 0.15, 0, -0.2, -1.4])

y_applied = choice_dt["applied"]

# Defining sizes
N_app = len(app_dt)
N_FE = len(X_FE.columns)
race_idx = choice_dt["race_norm_final"].astype("category").cat.codes.values
N_race = choice_dt["race_norm_final"].nunique()
N_obs = len(X_FE)

# transform dataframes to arrays
X_FE_arr = X_FE.values
# scale numeric columns
X_FE_arr[:,[0,2,3,4,5]]  = (X_FE_arr[:,[0,2,3,4,5]]  - X_FE_arr[:,[0,2,3,4,5]].mean(0))  / X_FE_arr[:,[0,2,3,4,5]].std(0)
y_arr = y_applied.values

with pm.Model() as model:

    # data
    X_data = pm.Data("X_data", X_FE_arr)
    race_idx_data = pm.Data("race_idx", race_idx)
    y_data = pm.Data("y_data", y_arr)

    # priors
    mu_beta = pm.Normal("mu_beta", fe_prior_mu, 0.1)
    sigma_beta = pm.Exponential("sigma_beta", 1)
    beta_fe = pm.Normal("beta_fe", mu=mu_beta, sigma=sigma_beta, shape=(N_race, N_FE))

    mu_app = pm.Normal("mu_app", mu=-4, sigma=1)
    sigma_app = pm.HalfNormal("sigma_app", sigma=0.3)

    # intercept
    gamma = pm.Normal("gamma", mu=mu_app, sigma=sigma_app)

    # utility calculations
    race_betas = beta_fe[race_idx_data]
    fixed_effects_term = pm.math.sum(race_betas * X_data, axis=1)
    eta = gamma + fixed_effects_term
    eta = eta.flatten() 
    
    p = pm.math.sigmoid(eta)
    p_var = pm.Deterministic("p", p)

    Y = pm.Bernoulli("Y", p=p, observed=y_data)

    # sampling
    prior_predictive = pm.sample_prior_predictive(draws=500, random_seed=42)
    idata = pm.sample(draws=1000, tune=1000, chains=4, random_seed=42)

beta_post = idata.posterior.beta_fe.values
az.summary(idata, var_names=["mu_app", "sigma_app", 'beta_fe', 'mu_beta', 'sigma_beta'])

Sampling: [Y, beta_fe, gamma, mu_app, mu_beta, sigma_app, sigma_beta]
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [mu_beta, sigma_beta, beta_fe, mu_app, sigma_app, gamma]

/Users/pavlomysak/opt/anaconda3/envs/chapa_env/lib/python3.10/site-packages/rich/live.py:256: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 1286 seconds.
There were 460 divergences after tuning. Increase `target_accept` or reparameterize.
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details
The effective sample size per chain is smaller than 100 for some parameters.  A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
mu_app	-5.241	0.303	-5.790	-4.602	0.008	0.012	1740.0	1348.0	1.01
sigma_app	0.249	0.179	0.014	0.569	0.006	0.003	269.0	92.0	1.01
beta_fe[0, 0]	-0.082	0.060	-0.190	0.035	0.002	0.001	1354.0	1470.0	1.00
beta_fe[0, 1]	0.111	0.105	-0.076	0.309	0.004	0.002	860.0	2091.0	1.01
beta_fe[0, 2]	0.079	0.055	-0.020	0.189	0.001	0.001	3251.0	1980.0	1.00
beta_fe[0, 3]	-0.016	0.102	-0.206	0.179	0.003	0.002	1578.0	2004.0	1.00
beta_fe[0, 4]	-0.289	0.060	-0.400	-0.166	0.003	0.003	436.0	81.0	1.01
beta_fe[0, 5]	-1.775	0.068	-1.922	-1.666	0.002	0.001	835.0	542.0	1.00
beta_fe[1, 0]	0.012	0.061	-0.087	0.132	0.002	0.001	717.0	779.0	1.00
beta_fe[1, 1]	0.118	0.105	-0.078	0.306	0.004	0.002	638.0	1440.0	1.01
beta_fe[1, 2]	0.079	0.040	-0.002	0.149	0.001	0.001	2376.0	2872.0	1.00
beta_fe[1, 3]	-0.015	0.104	-0.217	0.177	0.003	0.003	1434.0	1839.0	1.01
beta_fe[1, 4]	-0.177	0.062	-0.297	-0.065	0.003	0.001	499.0	1700.0	1.00
beta_fe[1, 5]	-1.779	0.070	-1.916	-1.655	0.003	0.001	649.0	1217.0	1.01
beta_fe[2, 0]	-0.001	0.056	-0.103	0.104	0.002	0.002	726.0	252.0	1.01
beta_fe[2, 1]	0.126	0.114	-0.092	0.352	0.007	0.006	341.0	106.0	1.01
beta_fe[2, 2]	0.060	0.051	-0.032	0.158	0.002	0.001	1072.0	1082.0	1.01
beta_fe[2, 3]	-0.014	0.103	-0.207	0.186	0.003	0.002	1310.0	703.0	1.00
beta_fe[2, 4]	-0.302	0.055	-0.404	-0.200	0.002	0.001	1065.0	2496.0	1.00
beta_fe[2, 5]	-1.650	0.062	-1.778	-1.540	0.002	0.001	954.0	1144.0	1.01
beta_fe[3, 0]	-0.004	0.092	-0.167	0.185	0.002	0.002	2163.0	1911.0	1.00
beta_fe[3, 1]	0.113	0.108	-0.086	0.313	0.004	0.002	897.0	1903.0	1.01
beta_fe[3, 2]	0.129	0.092	-0.036	0.307	0.003	0.002	1127.0	1849.0	1.00
beta_fe[3, 3]	-0.014	0.105	-0.215	0.186	0.002	0.003	1852.0	2162.0	1.01
beta_fe[3, 4]	-0.264	0.093	-0.424	-0.074	0.003	0.002	1184.0	1363.0	1.00
beta_fe[3, 5]	-1.679	0.094	-1.869	-1.519	0.003	0.002	806.0	1404.0	1.00
beta_fe[4, 0]	-0.020	0.086	-0.183	0.146	0.004	0.004	716.0	368.0	1.01
beta_fe[4, 1]	0.116	0.107	-0.086	0.312	0.003	0.002	953.0	1768.0	1.00
beta_fe[4, 2]	0.130	0.088	-0.007	0.331	0.004	0.004	513.0	112.0	1.00
beta_fe[4, 3]	-0.017	0.105	-0.224	0.173	0.003	0.002	1393.0	1742.0	1.00
beta_fe[4, 4]	-0.241	0.088	-0.408	-0.077	0.004	0.002	629.0	1286.0	1.01
beta_fe[4, 5]	-1.648	0.093	-1.824	-1.484	0.003	0.002	762.0	1734.0	1.00
beta_fe[5, 0]	0.008	0.072	-0.134	0.134	0.002	0.001	955.0	1871.0	1.00
beta_fe[5, 1]	0.124	0.109	-0.094	0.323	0.004	0.003	649.0	376.0	1.01
beta_fe[5, 2]	0.081	0.068	-0.047	0.204	0.001	0.002	1972.0	2346.0	1.00
beta_fe[5, 3]	-0.014	0.104	-0.207	0.188	0.003	0.002	1370.0	1449.0	1.00
beta_fe[5, 4]	-0.276	0.071	-0.419	-0.151	0.002	0.002	1011.0	730.0	1.00
beta_fe[5, 5]	-1.658	0.076	-1.794	-1.507	0.002	0.001	1015.0	2479.0	1.00
beta_fe[6, 0]	-0.010	0.030	-0.066	0.045	0.001	0.001	1306.0	453.0	1.00
beta_fe[6, 1]	0.087	0.084	-0.070	0.246	0.002	0.001	1288.0	2098.0	1.00
beta_fe[6, 2]	0.037	0.034	-0.030	0.098	0.001	0.001	1036.0	628.0	1.00
beta_fe[6, 3]	-0.022	0.033	-0.084	0.032	0.001	0.001	1386.0	1681.0	1.01
beta_fe[6, 4]	-0.362	0.037	-0.431	-0.296	0.002	0.001	514.0	1257.0	1.00
beta_fe[6, 5]	-1.705	0.048	-1.792	-1.613	0.001	0.001	1099.0	1767.0	1.00
mu_beta[0]	-0.012	0.043	-0.088	0.075	0.002	0.001	797.0	459.0	1.01
mu_beta[1]	0.117	0.076	-0.026	0.254	0.003	0.001	530.0	1475.0	1.01
mu_beta[2]	0.088	0.040	0.011	0.161	0.001	0.001	878.0	2458.0	1.01
mu_beta[3]	-0.014	0.065	-0.128	0.108	0.002	0.001	1031.0	1597.0	1.00
mu_beta[4]	-0.268	0.046	-0.353	-0.182	0.002	0.001	411.0	104.0	1.01
mu_beta[5]	-1.675	0.055	-1.772	-1.566	0.002	0.001	603.0	455.0	1.00
sigma_beta	0.077	0.027	0.030	0.125	0.002	0.001	248.0	669.0	1.00

The first diagnostic we examine is the Gelman–Rubin statistic ($\hat{R}$), which compares the variability between chains to the variability within chains. Values close to 1 indicate that the chains are sampling from the same target distribution, and a threshold of $\hat{R} \leq 1.01$ is generally considered evidence of satisfactory convergence. Because $\hat{R}$ alone does not guarantee convergence, we also evaluate the effective sample size (ESS).

The table reports two metrics:

• ess_bulk, which reflects the effective sample size for estimates in the bulk of the posterior (e.g., mean, median), and
• ess_tail, which measures the effective sample size in the distribution’s tails.

As a rule of thumb, values above 100 indicate adequate mixing. While a small number of parameters fall slightly below this threshold, the majority exhibit ESS values in the several hundreds or above 1,000, suggesting overall strong sampling efficiency.

We now supplement these numerical diagnostics with trace plots and autocorrelation plots, which provide a visual assessment of chain mixing and potential dependence. Taken together, these diagnostics offer a comprehensive view of the model’s convergence behavior.

# Check trace plots for beta
beta_posterior = beta_post.reshape(4000, 7, 6)

races = ['Asian', 'Black', 'Hispanic/Latino', 'Native American',
         'Two or More Races', 'Unknown', 'White']
features = ['Age', 'Disability', 'Dependents',
            'HH Income', 'Max Resale Price', 'Distance']

fig, axs = plt.subplots(
    nrows=len(races), ncols=len(features),
    figsize=(22, 18),
    sharex=True, sharey=False
)
plt.subplots_adjust(hspace=0.6, wspace=0.4)

for r in range(len(races)):
    for c in range(len(features)):
        vals = beta_posterior[:,r, c]

        axs[r, c].plot(vals)

# row/col labels
for r, race in enumerate(races):
    axs[r, 0].set_ylabel(race, fontsize=12, rotation=0, labelpad=50)
for c, feat in enumerate(features):
    axs[0, c].set_title(feat, fontsize=12)

fig.suptitle('Trace Plots for Beta', fontsize=25)
plt.show()

/var/folders/q9/vssjztfs6gd47qy2w2t2f0t80000gn/T/ipykernel_82210/729612042.py:14: UserWarning: This figure was using a layout engine that is incompatible with subplots_adjust and/or tight_layout; not calling subplots_adjust.
  plt.subplots_adjust(hspace=0.6, wspace=0.4)

# trace plots for gamma
fig, axs = plt.subplots(nrows=2, ncols=2, sharex=True)
axs = axs.flatten()

for chain in range(4):
    axs[chain].plot(idata.posterior.gamma.values[chain, :])
    axs[chain].set_title(f"Chain {chain} Trace")
    axs[chain].set_yticks([])
fig.suptitle('Trace Plots for Gamma', fontsize=15)
plt.show()

# autocorrelation for gamma
az.plot_autocorr(idata.posterior.gamma)
plt.show()

The trace plots show no evidence of persistent stickiness or slow exploration across any of the beta coefficients or the gamma parameter, indicating that the sampler mixed well across chains.

A brief inspection of the autocorrelation function (illustrated here for the gamma parameter) shows a rapid decay, approaching zero by approximately lag 15. Given the complexity and noisiness typical of real-world data, this level of autocorrelation is entirely acceptable and further supports the conclusion that the sampler has converged.

With convergence diagnostics complete, we now proceed to the final step before inference: the posterior predictive check. With a dataset like ours, where our target value is extremely unbalanced, a typical visual posterior predictive check won't be incredibly insightful. We revert to PPP-Values (Posterior Predictive P-Values) to help us.

with model:
    pm.sample_posterior_predictive(idata, extend_inferencedata=True, random_seed=42)

Sampling: [Y]

/Users/pavlomysak/opt/anaconda3/envs/chapa_env/lib/python3.10/site-packages/rich/live.py:256: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

az.plot_bpv(idata, kind="t_stat", t_stat="mean")
plt.show()

A low posterior predictive p-value suggests poor model fit, while a value near 0.5 indicates that the model generates data consistent with what was observed. Our posterior predictive p-value, which is approximately 0.5, reflects strong agreement between the model and the data and provides additional reassurance that the model is appropriate for inference.

Many of our key questions of interest are addressed through the beta coefficients. The inspection of these coefficients follows.

fig, axs = plt.subplots(nrows=1, ncols=beta_posterior.shape[2], figsize=(15,6))
axs = axs.flatten()

for feat in range(beta_posterior.shape[2]):

    vals = beta_posterior[:, :, feat]
    means = vals.mean(axis=0)
    ci50 = np.percentile(vals, [25, 75], axis=0)
    ci95 = np.percentile(vals, [2.5, 97.5], axis=0)
    y = np.arange(len(races))

    axs[feat].hlines(y, ci95[0], ci95[1], color='black')
    axs[feat].hlines(y, ci50[0], ci50[1], color='black', linewidth=4)
    axs[feat].plot(means, y, 'o', color='black')
    axs[feat].vlines(means.mean(), ymin=0, ymax=6, color="red", ls="--", alpha=0.2)
    axs[feat].set_title(f"{features[feat]}")

axs[0].set_yticks(y)
axs[0].set_yticklabels(races)
for feat in range(beta_posterior.shape[2]):
    axs[feat].invert_yaxis()

plt.figure(dpi=600)
plt.show()

<Figure size 4320x2880 with 0 Axes>

In the forest plot above, we see each covariates impact on the probability of an applicant applying. The dotted red line represents the between-group mean, helping us visually assess which groups might deviate from the group norm.

Broadly, the model finds that most applicant-level factors exert similar effects across racial and ethnic groups, with meaningful differences emerging primarily in responses to property attributes.

Across all groups, the effects associated with age and household income were near zero, indicating that these characteristics do not materially alter application behavior. A slight deviation appears among Asian applicants, for whom the associated probability of applying decreases as age increases.

Disability status and the number of dependents showed consistent positive associations with application probability. Applicants reporting a disability were more likely to apply, and the likelihood of applying increased incrementally with each additional dependent. Importantly, these effects did not differ meaningfully across demographic groups, suggesting that these factors operate similarly regardless of racial or ethnic background.

Group-specific variation becomes more evident when examining property-level characteristics. As the distance between an applicant’s residence and a property increased, Asian and Black applicants exhibited a sharper decline in the probability of applying than other groups. Hispanic/Latino applicants, by contrast, showed comparatively little sensitivity to distance, implying that geographic proximity plays a smaller role in their decision-making. Differences also emerged in responses to maximum resale price. Higher resale prices were associated with lower application likelihood for White applicants, while Black applicants appeared less sensitive to price. Although these effects were modest in magnitude, they were directionally consistent and suggest that price and location considerations are evaluated differently across groups.

The results indicate that while baseline applicant characteristics generally influence all groups in similar ways, heterogeneity arises in how groups respond to key property features. These findings highlight that distance and pricing structures contribute to observable differences in application behavior, even as other determinants remain stable across demographic categories.

Below is an alternate look at these posterior distributions.

fig, axs = plt.subplots(
    nrows=len(races), ncols=len(features),
    figsize=(22, 18),
    sharex=False, sharey=False, dpi=600
)

plt.subplots_adjust(hspace=0.6, wspace=0.4)

for r in range(len(races)):
    for c in range(len(features)):
        vals = beta_posterior[:, r, c]

        axs[r, c].hist(vals, density=True, bins=30)

        axs[r, c].axvline(0, color="red")

        # remove extra ticks
        axs[r, c].set_yticks([])
        
        if r != (len(races)-1):
            axs[r, c].set_xticks([])


        axs[r, c].tick_params(labelsize=8)

        if c == (len(features)-1):
            axs[r, c].set_xlim(-2.2, -1.2)
        else:
            axs[r, c].set_xlim(-0.45, 0.45)

# row/col labels
for r, race in enumerate(races):
    axs[r, 0].set_ylabel(race, fontsize=12, rotation=0, labelpad=50)
for c, feat in enumerate(features):
    axs[0, c].set_title(feat, fontsize=12)

plt.show()

/var/folders/q9/vssjztfs6gd47qy2w2t2f0t80000gn/T/ipykernel_82210/2044758254.py:7: UserWarning: This figure was using a layout engine that is incompatible with subplots_adjust and/or tight_layout; not calling subplots_adjust.
  plt.subplots_adjust(hspace=0.6, wspace=0.4)

fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(15,6))
axs = axs.flatten()

for feat in range(2):

    vals = beta_posterior[:, :, feat+4]
    means = vals.mean(axis=0)
    ci50 = np.percentile(vals, [25, 75], axis=0)
    ci95 = np.percentile(vals, [2.5, 97.5], axis=0)
    y = np.arange(len(races))

    axs[feat].hlines(y, ci95[0], ci95[1], color='black')
    axs[feat].hlines(y, ci50[0], ci50[1], color='black', linewidth=4)
    axs[feat].plot(means, y, 'o', color='black')
    axs[feat].vlines(means.mean(), ymin=0, ymax=6, color="red", ls="--", alpha=0.2)
    axs[feat].set_title(f"{features[feat+4]}")

axs[0].set_yticks(y)
axs[0].set_yticklabels(races)
for feat in range(2):
    axs[feat].invert_yaxis()

plt.figure(dpi=600)
plt.show()

<Figure size 4320x2880 with 0 Axes>

Above is a closer look at the two covariates showing the greatest heterogeneity across groups. While we observe notable differences, it is worth noting that in many cases the 95% credible intervals still overlap with the overall group means.

Modeling Assumptions and Identifiability

Structural Assumptions of the Discrete-Choice Framework

We assume that an applicant’s decision to apply can be represented through a latent utility model where the log-odds of application are a linear function of applicant characteristics, property attributes, and group-specific effects. This implies that:

• covariate effects are additive,
• applicants consider each property independently, and
• unobserved factors influencing application behavior are absorbed into the error term and group-level intercepts.

Hierarchical Structure and Partial Pooling

The model includes group-specific coefficients for racial and ethnic groups. Partial pooling assumes:

• groups share an underlying population-level distribution of effects,
• differences between groups arise as deviations from a common mean, and
• shrinkage appropriately stabilizes estimates for groups with limited data.

Prior Specification and Weak Identifiability

Bayesian models with hierarchical random effects can exhibit weak identifiability without appropriate priors. To address this:

• fixed-effect priors provide regularization to avoid extreme log-odds estimates,
• priors on the group-level variance constrain the degree of heterogeneity, and
• the intercept prior incorporates known class imbalance to anchor the baseline probability.

Multicollinearity and Parameter Interpretability

Highly correlated covariates (e.g., household size and dependents) can lead to non-identifiability of individual coefficients. To preserve interpretability and identifiability:

• one of the collinear variables was removed in later model iterations, and
• priors were adjusted to reflect plausible effect sizes.

Exchangeability Assumptions

Within-group applicants are assumed exchangeable after conditioning on covariates. Similarly, group-level effects are assumed exchangeable under the hierarchical prior. These assumptions enable borrowing strength across groups while still allowing meaningful group-level differences.

---

Appendix - Iterative Model Refinement Process

$$ \text{Model 1 Description} $$

Prior Structure and Rationale:

In Model 1A, I place a weakly informative regularizing prior of $\mathcal{N}(0, 0.5)$ on the coefficients. This scale is small enough to shrink implausibly large effects toward zero, while still allowing the data substantial influence. It is an appropriate starting point for a logistic model where we do not expect extreme log-odds shifts from most covariates.

For the applicant-level intercept, the data exhibit substantial class imbalance, with far fewer applications than non-applications. It is therefore reasonable to expect the baseline log-odds of applying to be negative. To encode this, I use a prior of $\mathcal{N}(-2, 1)$, which reflects a conservative belief in a negative baseline probability while remaining broad enough to let the likelihood meaningfully update the intercept.

Finally, for the applicant-level slope variance, I use a half-normal prior with scale 1. This provides a weakly informative constraint that keeps group-level variability within a plausible range and stabilizes sampling, while not imposing strong assumptions about the exact magnitude of heterogeneity across applicants.

Notice I make use of the normal distribution often due to the nice property of parameter interpretability. See model outline below.

---

Priors $$ \text{Fixed Effects Priors: } \beta \sim \mathcal{N}(0,0.5) $$ $$ \text{Applicant Mean Prior: } \mu_{\text{app}} \sim \mathcal{N}(-2,1) $$ $$ \text{Applicant Variance Prior: } \sigma_{\text{app}}^2 \sim \text{Half-Normal}(1) $$

---

Applicant-Level Random Effects

$$ \gamma_i \sim \mathcal{N}(\mu_{\text{app}}, \sigma_{\text{app}}^2) $$

---

Utility

$$ \eta_{ij} = \gamma_i + \mathbf{X_{ij}^T}\beta $$$$ p_{ij} = \frac{1}{1+\text{exp}(-\eta_{ij})} $$$$ Y_{ij} \sim \text{Bernoulli}(p_{ij}) $$

$\text{Model 1A: A First Look}$

features = ['Age', 'Disability', 'HH Size', 'Dependents', 'HH Income', 'Property Maximum Resale Price', 'distance_miles']
X_FE = choice_dt[features]
y_applied = choice_dt["applied"]

# Defining sizes
N_app = len(app_dt)
N_FE = len(X_FE.columns)
N_int = 1 #len(X_int.columns)
N_obs = len(X_FE)

app_id_to_idx = {id_num: idx for idx, id_num in enumerate(app_dt['ID Number'].unique())}
app_idx_sequential = np.array([app_id_to_idx[id_num] for id_num in choice_dt['ID Number']])

# transform dataframes to arrays
X_FE_arr = X_FE.values
# scale numeric columns
X_FE_arr[:,[0,2,3,4,5,6]]  = (X_FE_arr[:,[0,2,3,4,5,6]]  - X_FE_arr[:,[0,2,3,4,5,6]].mean(0))  / X_FE_arr[:,[0,2,3,4,5,6]].std(0)
y_arr = y_applied.values

with pm.Model() as model:

    # creating data objs
    X_data = pm.Data("X_data", X_FE_arr)
    app_idx = pm.Data("app_idx", app_idx_sequential)
    y_data = pm.Data("y_data", y_arr)

    # priors
    beta_fe = pm.Normal("beta_fe", mu=0, sigma=0.5, shape=N_FE)
    
    mu_app = pm.Normal("mu_app", mu=-2, sigma=1)
    sigma_app = pm.HalfNormal("sigma_app", sigma=1)

    # Applicant-level RE
    gamma_i = pm.Normal("gamma_i", mu=mu_app, sigma=sigma_app, shape=N_app)

    # Utility calculations
    gamma_for_obs = gamma_i[app_idx].reshape((-1,1))
    fixed_effects_term = X_data @ beta_fe.reshape((-1,1))
    eta = gamma_for_obs + fixed_effects_term
    eta = eta.flatten() 
    
    p = pm.invlogit(eta)
    p_var = pm.Deterministic("p", p)

    Y = pm.Bernoulli("Y", p=p, observed=y_data)

    # sampling
    prior_predictive = pm.sample_prior_predictive(draws=500, random_seed=42)
    idata = pm.sample(draws=1000, tune=1000, chains=4, random_seed=42)

beta_post = idata.posterior.beta_fe.values
az.summary(idata, var_names=["mu_app", "sigma_app", 'beta_fe'])

Sampling: [Y, beta_fe, gamma_i, mu_app, sigma_app]
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [beta_fe, mu_app, sigma_app, gamma_i]

/Users/pavlomysak/opt/anaconda3/envs/chapa_env/lib/python3.10/site-packages/rich/live.py:256: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 925 seconds.
There were 2116 divergences after tuning. Increase `target_accept` or reparameterize.
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details
The effective sample size per chain is smaller than 100 for some parameters.  A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
mu_app	-5.375	0.033	-5.424	-5.316	0.014	0.005	6.0	54.0	1.78
sigma_app	0.065	0.018	0.039	0.093	0.006	0.002	8.0	12.0	1.49
beta_fe[0]	-0.012	0.022	-0.056	0.025	0.004	0.002	28.0	252.0	1.11
beta_fe[1]	0.040	0.084	-0.121	0.203	0.009	0.005	73.0	362.0	1.05
beta_fe[2]	0.006	0.043	-0.061	0.087	0.014	0.002	10.0	96.0	1.32
beta_fe[3]	0.050	0.042	-0.024	0.120	0.011	0.002	14.0	108.0	1.21
beta_fe[4]	-0.018	0.033	-0.078	0.031	0.004	0.001	56.0	435.0	1.06
beta_fe[5]	-0.317	0.032	-0.373	-0.254	0.011	0.004	9.0	34.0	1.40
beta_fe[6]	-1.756	0.034	-1.832	-1.699	0.010	0.004	10.0	26.0	1.30

Notice our $\hat{R}$ and effective sample sizes (ess) are quite poor. Let's continue to investigate and see if we get any ideas on how to fix this.

fig, axs = plt.subplots(ncols = beta_post.shape[0], nrows = beta_post.shape[2], figsize=(11,7))

for chain in range(beta_post.shape[0]):
    for col in range(beta_post.shape[2]):
        axs[col, chain].plot(beta_post[chain,:,col])
        axs[col, chain].set_yticks([])
        axs[col, chain].set_xticks([])

for chain in range(beta_post.shape[0]):
    axs[0, chain].set_title(f"chain {chain+1}")
plt.show()

fig, axs = plt.subplots(ncols = 3, nrows = 3, figsize=(11,7))
axs = axs.flatten()

for col in range(beta_post.shape[2]):
    for chain in range(beta_post.shape[0]):
        axs[col].hist(beta_post[chain,:,col], density=True, label=f"chain {chain+1}", alpha = 0.4)
    axs[col].set_title(f"{features[col]}", fontsize=8)
    
axs[2].legend()
plt.show()

Notice some stickiness in our trace and disagreement in our posterior distributions (especially in chains 3 & 4). This can be caused by a few things, but it looks like our model is having a hard time sampling. Next things to check: priors & multicolinearity.

# it also looks like our priors can be wayyy more informative...
az.plot_forest(prior_predictive.prior, var_names=["beta_fe"], combined=True)

array([<Axes: title={'center': '94.0% HDI'}>], dtype=object)

Our priors are relatively diffuse at this stage, allowing parameters to take on a wide range of values. These can be tightened using domain knowledge without introducing problematic bias. For instance, we know that higher maximum resale prices should reduce the probability of application, and the priors can be constructed to reflect this structural expectation.

We also observe some disagreement between chains for the HH Size and Dependents coefficients. Given the strong correlation between these variables, this likely reflects multicollinearity rather than sampling issues, and suggests that both should not be included simultaneously.

In the next iteration, we address this multicollinearity and refine the prior specifications accordingly.

$ \text{Model 1B: More informative priors, reducing multicollinearity} $

Updated Prior Structure & Rationale:

The fixed‐effect priors have been updated to incorporate domain knowledge and improve regularization. Each covariate now has its own prior mean to reflect directional expectations:

• Age: 0
• Disability: 0.15
• Dependents: 0.15
• Household Income: 0
• Maximum Resale Price: –0.2
• Distance: –1.4

All fixed‐effect priors now share a tightened variance of 0.1, which places more mass around reasonable effect sizes while still allowing the data to update these parameters meaningfully.

Based on diagnostics from the initial model, I also revised the prior on the applicant-level intercept. Given the substantial class imbalance and very low baseline probability of application, a more negative prior mean is justified. I therefore updated the intercept prior to $\mathcal{N}(-4, 1)$.

Finally, the variance of the applicant-level random intercept is now modeled with a HalfNormal(0.3) prior. This tighter prior helps stabilize estimation by preventing unrealistically large between-applicant variability, which was evident in the earlier iteration.

# we believe there to be multicollinearity between 'HH Size' and 'Dependents'. Which shall we choose to go into the model?
choice_dt[['HH Size', 'Dependents', 'applied']].corr()

	HH Size	Dependents	applied
HH Size	1.000000	0.776442	0.003878
Dependents	0.776442	1.000000	0.004909
applied	0.003878	0.004909	1.000000

# let's also add some stronger priors!!
features = ['Age', 'Disability', 'Dependents', 'HH Income', 'Property Maximum Resale Price', 'distance_miles']
X_FE = choice_dt[features]
fe_prior_mu = np.array([0, 0.15, 0.15, 0, -0.2, -1.4])

# Defining sizes
N_app = len(app_dt)
N_FE = len(X_FE.columns)
N_obs = len(X_FE)
app_id_to_idx = {id_num: idx for idx, id_num in enumerate(app_dt['ID Number'].unique())}
app_idx_sequential = np.array([app_id_to_idx[id_num] for id_num in choice_dt['ID Number']])

# transform dataframes to arrays
X_FE_arr = X_FE.values
# scale numeric columns
X_FE_arr[:,[0,2,3,4,5]]  = (X_FE_arr[:,[0,2,3,4,5]]  - X_FE_arr[:,[0,2,3,4,5]].mean(0))  / X_FE_arr[:,[0,2,3,4,5]].std(0)
y_arr = y_applied.values

with pm.Model() as model:

    # creating data objs
    X_data = pm.Data("X_data", X_FE_arr)
    app_idx = pm.Data("app_idx", app_idx_sequential)
    y_data = pm.Data("y_data", y_arr)

    # priors
    beta_fe = pm.Normal("beta_fe", mu=fe_prior_mu, sigma=0.1)
    
    mu_app = pm.Normal("mu_app", mu=-4, sigma=1)
    sigma_app = pm.HalfNormal("sigma_app", sigma=0.3)

    # Applicant-level RE
    gamma_i = pm.Normal("gamma_i", mu=mu_app, sigma=sigma_app, shape=N_app)

    # Utility calculations
    gamma_for_obs = gamma_i[app_idx].reshape((-1,1))
    fixed_effects_term = X_data @ beta_fe.reshape((-1,1))
    eta = gamma_for_obs + fixed_effects_term
    eta = eta.flatten() 
    
    p = pm.invlogit(eta)
    p_var = pm.Deterministic("p", p)

    Y = pm.Bernoulli("Y", p=p, observed=y_data)

    # sampling
    prior_predictive = pm.sample_prior_predictive(draws=500, random_seed=42)
    idata = pm.sample(draws=1000, tune=1500, chains=4, random_seed=42)

beta_post = idata.posterior.beta_fe.values
az.summary(idata, var_names=["mu_app", "sigma_app", 'beta_fe'])

Sampling: [Y, beta_fe, gamma_i, mu_app, sigma_app]
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [beta_fe, mu_app, sigma_app, gamma_i]

/Users/pavlomysak/opt/anaconda3/envs/chapa_env/lib/python3.10/site-packages/rich/live.py:256: UserWarning: install 
"ipywidgets" for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

Sampling 4 chains for 1_500 tune and 1_000 draw iterations (6_000 + 4_000 draws total) took 1056 seconds.
There were 625 divergences after tuning. Increase `target_accept` or reparameterize.
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details
The effective sample size per chain is smaller than 100 for some parameters.  A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
mu_app	-5.357	0.036	-5.429	-5.289	0.008	0.010	22.0	15.0	1.41
sigma_app	0.057	0.027	0.025	0.108	0.012	0.005	6.0	31.0	1.70
beta_fe[0]	-0.015	0.024	-0.056	0.031	0.001	0.000	585.0	1932.0	1.01
beta_fe[1]	0.109	0.068	-0.016	0.236	0.002	0.002	973.0	1502.0	1.00
beta_fe[2]	0.060	0.023	0.019	0.105	0.001	0.001	572.0	700.0	1.01
beta_fe[3]	-0.022	0.033	-0.084	0.033	0.002	0.002	472.0	289.0	1.01
beta_fe[4]	-0.311	0.028	-0.362	-0.258	0.002	0.001	168.0	185.0	1.01
beta_fe[5]	-1.728	0.036	-1.798	-1.662	0.006	0.006	41.0	41.0	1.10

Notice, our $\hat{R}$ is starting to improve a bit, but it's not great yet...

fig, axs = plt.subplots(ncols = beta_post.shape[0], nrows = beta_post.shape[2], figsize=(11,7))

for chain in range(beta_post.shape[0]):
    for col in range(beta_post.shape[2]):
        axs[col, chain].plot(beta_post[chain,:,col])
        axs[col, chain].set_yticks([])
        axs[col, chain].set_xticks([])

for chain in range(beta_post.shape[0]):
    axs[0, chain].set_title(f"chain {chain+1}")
plt.show()

There's some stickiness happening in chain 2 and chain 3. Will this come across in the distributions?

fig, axs = plt.subplots(ncols = 3, nrows = 2, figsize=(11,7))
axs = axs.flatten()

for col in range(beta_post.shape[2]):
    for chain in range(beta_post.shape[0]):
        axs[col].hist(beta_post[chain,:,col], density=True, label=f"chain {chain+1}", alpha = 0.4)
    axs[col].set_title(f"{features[col]}", fontsize=8)
    
axs[2].legend()
plt.show()

# these priors look way better.
az.plot_forest(prior_predictive.prior, var_names=["beta_fe"], combined=True)

array([<Axes: title={'center': '94.0% HDI'}>], dtype=object)

Priors look better and model fit has improved (convergence diagnostics)... but we can do better. Let's add race/ethnicity as a group-level effect, drop the applicant-level random effects (too computationally expensive), and retain our current informative priors.

See the Model 2 section to view the final model results.