xjtulyc Social Stratification and Intergenerational Mobility
When to Use
Use this skill when you need to:
- Estimate intergenerational income elasticity (IGE) via log-log OLS regression
- Estimate the rank-rank slope (Chetty et al. 2014 approach) for more robust mobility estimates
- Compute upward mobility rates — probability of reaching the top quintile given bottom quintile origins
- Load and use the ISEI (International Socio-Economic Index of Occupational Status) to score occupations on a continuous prestige scale
- Apply the EGP (Erikson-Goldthorpe-Portocarero) class schema to categorize occupations into hierarchical service/intermediate/working class groups
- Build and analyze intergenerational occupational transition matrices
- Perform odds ratio analysis for social fluidity (relative mobility)
- Decompose income variance by education and occupation contributions (Shapley decomposition)
Background
Intergenerational Elasticity (IGE): Regress log child income on log parent income:
log(income_child) = α + β × log(income_parent) + ε
β is the IGE. High β (close to 1) = low mobility (child income closely tracks parent income). The US IGE is approximately 0.45; Nordic countries are closer to 0.15-0.25.
Rank-Rank Slope: Convert income to percentile ranks within cohort, then regress:
rank_child = α + ρ × rank_parent + ε
ρ (rank-rank slope) is less sensitive to outliers than IGE and better suited to censored or top-coded income data. Chetty et al. (2014) found ρ ≈ 0.341 for the United States.
Upward Mobility: P(child in Q5 | parent in Q1) — fraction of children born in the bottom income quintile who reach the top quintile as adults. This measure (sometimes called "Chetty mobility") varies enormously across geographic areas and demographic groups.
ISEI (Hauser & Warren 1997): An occupation scoring system derived from the regression of income and education on the Standard Occupational Classification. Scores range from ~16 (lowest manual) to ~90 (physicians, judges). Assigned from 4-digit ISCO or national occupation codes.
EGP Schema (Erikson, Goldthorpe & Portocarero 1979): Categorical class schema:
- I: Higher professionals, administrators (service class)
- II: Lower professionals, technicians (service class)
- IIIa: Routine non-manual (intermediate)
- IIIb: Personal services (intermediate)
- IVa/IVb/IVc: Self-employed (petty bourgeoisie)
- V/VI: Skilled manual workers
- VIIa/VIIb: Unskilled manual, agricultural workers
Transition matrix and social fluidity: The (origin × destination) mobility table shows probabilities of ending in each class given each origin class. The odds ratio:
OR(a,b;c,d) = (n_ac × n_bd) / (n_ad × n_bc)
where a,b are origin classes and c,d are destination classes. OR = 1 means equal relative odds (perfect fluidity). Deviations from 1 indicate barriers to mobility.
Environment Setup
pip install pandas>=1.5 statsmodels>=0.14 numpy>=1.23 matplotlib>=3.6 scipy>=1.9
PSID (Panel Study of Income Dynamics) data: https://psidonline.isr.umich.edu/ Register for free access. Download cross-year individual/family files.
export PSID_PATH="/data/psid_crossyear.csv"
Core Workflow
import os
import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
from scipy import stats
import warnings
warnings.filterwarnings("ignore")
# ---------------------------------------------------------------------------
# 1. IGE and Rank-Rank Slope
# ---------------------------------------------------------------------------
def intergenerational_elasticity(
parent_income: pd.Series,
child_income: pd.Series,
) -> dict:
"""
Estimate intergenerational income elasticity (log-log OLS).
Parameters
----------
parent_income, child_income : pd.Series
Income series (positive values; zeros will be dropped).
Returns
-------
dict with ige (beta), r_squared, n, intercept.
"""
mask = (parent_income > 0) & (child_income > 0)
lp = np.log(parent_income[mask])
lc = np.log(child_income[mask])
slope, intercept, r, p, se = stats.linregress(lp, lc)
return {
"ige": round(slope, 4),
"intercept": round(intercept, 4),
"r_squared": round(r ** 2, 4),
"se": round(se, 4),
"p_value": round(p, 4),
"n": int(mask.sum()),
}
def rank_rank_slope(
parent_income: pd.Series,
child_income: pd.Series,
bootstrap_n: int = 500,
seed: int = 42,
) -> dict:
"""
Estimate rank-rank slope with bootstrap confidence interval (Chetty method).
Parameters
----------
parent_income, child_income : pd.Series
bootstrap_n : int
Number of bootstrap replications for CI.
seed : int
Returns
-------
dict with rr_slope, intercept, ci_lower, ci_upper, se, n.
"""
mask = parent_income.notna() & child_income.notna()
pr = parent_income[mask].rank(pct=True)
cr = child_income[mask].rank(pct=True)
slope, intercept, r, p, se = stats.linregress(pr, cr)
# Bootstrap CI
rng = np.random.default_rng(seed)
boot_slopes = []
n = len(pr)
pr_arr, cr_arr = pr.values, cr.values
for _ in range(bootstrap_n):
idx = rng.integers(0, n, n)
s, *_ = stats.linregress(pr_arr[idx], cr_arr[idx])
boot_slopes.append(s)
ci_lo, ci_hi = np.percentile(boot_slopes, [2.5, 97.5])
return {
"rr_slope": round(slope, 5),
"intercept": round(intercept, 5),
"ci_lower": round(ci_lo, 5),
"ci_upper": round(ci_hi, 5),
"se": round(np.std(boot_slopes), 5),
"n": int(mask.sum()),
}
def upward_mobility_rate(
parent_income: pd.Series,
child_income: pd.Series,
n_quintiles: int = 5,
from_quintile: int = 1,
to_quintile: int = 5,
) -> dict:
"""
Compute upward mobility: P(child in top quintile | parent in bottom quintile).
Returns
-------
dict with p_upward, n_origin, n_both.
"""
mask = parent_income.notna() & child_income.notna()
pq = pd.qcut(parent_income[mask], n_quintiles, labels=False) + 1
cq = pd.qcut(child_income[mask], n_quintiles, labels=False) + 1
origin_mask = pq == from_quintile
n_origin = origin_mask.sum()
n_both = ((pq == from_quintile) & (cq == to_quintile)).sum()
return {
"p_upward": round(n_both / n_origin, 5) if n_origin > 0 else np.nan,
"n_origin": int(n_origin),
"n_destination": int(n_both),
"from_quintile": from_quintile,
"to_quintile": to_quintile,
}
# ---------------------------------------------------------------------------
# 2. ISEI Occupational Prestige
# ---------------------------------------------------------------------------
# Simplified ISEI lookup by ISCO-08 major group (real application uses 4-digit codes)
ISEI_ISCO_MAJOR = {
1: 68, # Managers
2: 74, # Professionals
3: 56, # Technicians and associate professionals
4: 40, # Clerical support workers
5: 32, # Services and sales workers
6: 23, # Skilled agricultural workers
7: 34, # Craft and related trades
8: 31, # Plant and machine operators
9: 20, # Elementary occupations
0: 47, # Armed forces
}
def assign_isei(occupation_codes: pd.Series, lookup: dict | None = None) -> pd.Series:
"""
Assign ISEI scores from occupation codes.
Parameters
----------
occupation_codes : pd.Series
ISCO major group codes (1-digit integers for simplified lookup).
lookup : dict, optional
Custom {occupation_code: isei_score} mapping.
Returns
-------
pd.Series of ISEI scores.
"""
lkp = lookup or ISEI_ISCO_MAJOR
return occupation_codes.map(lkp)
# EGP schema: mapping from ISCO major group to EGP class (simplified)
EGP_ISCO_MAJOR = {
1: "I", # Managers → Higher service
2: "I", # Professionals → Higher service
3: "II", # Technicians → Lower service
4: "IIIa", # Clerical → Routine non-manual
5: "IIIb", # Service workers → Routine non-manual
6: "IVc", # Agricultural self-employed
7: "V_VI", # Craft → Skilled manual
8: "V_VI", # Operators → Skilled manual
9: "VIIa", # Elementary → Unskilled manual
}
EGP_ORDER = ["I", "II", "IIIa", "IIIb", "IVa", "IVb", "IVc", "V_VI", "VIIa", "VIIb"]
def assign_egp(occupation_codes: pd.Series, lookup: dict | None = None) -> pd.Series:
"""Assign EGP class labels from occupation codes."""
lkp = lookup or EGP_ISCO_MAJOR
return occupation_codes.map(lkp)
# ---------------------------------------------------------------------------
# 3. Intergenerational Transition Matrix
# ---------------------------------------------------------------------------
def mobility_transition_matrix(
origin_class: pd.Series,
destination_class: pd.Series,
class_order: list[str] | None = None,
normalize: str = "origin",
) -> pd.DataFrame:
"""
Compute an intergenerational class transition matrix.
Parameters
----------
origin_class : pd.Series
Parent's (or respondent's father's) class.
destination_class : pd.Series
Respondent's current class.
class_order : list of str, optional
Ordered list of class labels for the matrix.
normalize : str
'origin' (row percentages), 'destination' (column), or 'all'.
Returns
-------
pd.DataFrame — transition proportions (percentages).
"""
classes = class_order or sorted(set(origin_class.dropna()) | set(destination_class.dropna()))
tab = pd.crosstab(
origin_class, destination_class,
values=np.ones(len(origin_class)), aggfunc="sum",
normalize=normalize,
).reindex(index=classes, columns=classes, fill_value=0)
return (tab * 100).round(1)
def compute_odds_ratio(
table: pd.DataFrame,
class_a: str,
class_b: str,
class_c: str,
class_d: str,
) -> float:
"""
Compute odds ratio for social fluidity from a mobility table.
OR = (n_ac × n_bd) / (n_ad × n_bc)
"""
n_ac = table.loc[class_a, class_c]
n_bd = table.loc[class_b, class_d]
n_ad = table.loc[class_a, class_d]
n_bc = table.loc[class_b, class_c]
if n_ad == 0 or n_bc == 0:
return np.nan
return (n_ac * n_bd) / (n_ad * n_bc)
Advanced Usage
Shapley Decomposition of Income Variance
import numpy as np
import pandas as pd
import statsmodels.api as sm
def shapley_r2_decomposition(
outcome: pd.Series,
predictors: dict[str, pd.Series],
) -> pd.DataFrame:
"""
Shapley decomposition of R² across predictor groups.
For each predictor group, compute the average marginal contribution to R²
across all possible orderings (approximated by sequential addition and deletion).
Parameters
----------
outcome : pd.Series
predictors : dict {group_name: pd.Series}
Returns
-------
pd.DataFrame with group, shapley_r2, pct_contribution.
"""
from itertools import combinations
groups = list(predictors.keys())
df_all = pd.concat([outcome] + list(predictors.values()), axis=1).dropna()
y = df_all.iloc[:, 0]
def r2_for_subset(group_subset):
if not group_subset:
return 0.0
X = sm.add_constant(df_all[[g for g in group_subset]])
model = sm.OLS(y, X).fit()
return model.rsquared
n = len(groups)
shapley = {g: 0.0 for g in groups}
weight = 1.0 / n
for size in range(n):
for subset in combinations(groups, size):
base = r2_for_subset(list(subset))
for g in groups:
if g not in subset:
marginal = r2_for_subset(list(subset) + [g]) - base
shapley[g] += marginal / (n * len(list(combinations(
[x for x in groups if x != g], size))))
total = sum(shapley.values())
rows = [{"group": g, "shapley_r2": round(v, 5),
"pct_contribution": round(v / total * 100, 2) if total > 0 else np.nan}
for g, v in shapley.items()]
return pd.DataFrame(rows).sort_values("shapley_r2", ascending=False)
Troubleshooting
| Problem | Cause | Solution |
|---|---|---|
| IGE > 1 or < 0 | Outliers in income distribution | Winsorize at 1st/99th percentile before log transform |
| Rank-rank slope unstable | Small sample size | Bootstrap CI is wide with N < 500; increase sample or report CI |
| ISEI mapping fails | Non-standard occupation codes | Convert national codes to ISCO-08 first using crosswalk table |
| Transition matrix rows don't sum to 100 | Missing values in one class | Check for NaN after assign_egp(); assign "Unknown" category |
| Odds ratio is NaN | Zero cell in the table | Collapse rare classes; add 0.5 Laplace smoothing |
External Resources
- Chetty, R. et al. (2014). Where is the land of opportunity? QJE, 129(4), 1553-1623.
- Erikson, R. & Goldthorpe, J.H. (1992). The Constant Flux. Oxford University Press.
- Hauser, R.M. & Warren, J.R. (1997). Socioeconomic indexes for occupations. Sociological Methodology, 27, 177-298.
- PSID Data Center: https://psidonline.isr.umich.edu/
- Opportunity Insights Data Library: https://opportunityinsights.org/data/
Examples
Example 1: Rank-Rank Regression with Bootstrap CI
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
rng = np.random.default_rng(42)
n = 8000
parent_income = np.exp(rng.normal(10.5, 0.9, n)).clip(5000, 5_000_000)
# True rho = 0.35
noise = rng.normal(0, 0.7, n)
child_income = np.exp(0.35 * np.log(parent_income) + 0.65 * 10.5 + noise).clip(5000, 5_000_000)
df_mob = pd.DataFrame({"parent_income": parent_income, "child_income": child_income})
ige_result = intergenerational_elasticity(df_mob["parent_income"], df_mob["child_income"])
rr_result = rank_rank_slope(df_mob["parent_income"], df_mob["child_income"])
upward = upward_mobility_rate(df_mob["parent_income"], df_mob["child_income"])
print("=== Intergenerational Mobility ===")
print(f" IGE (log-log): {ige_result['ige']} (se={ige_result['se']})")
print(f" Rank-rank slope: {rr_result['rr_slope']} "
f"[{rr_result['ci_lower']}, {rr_result['ci_upper']}]")
print(f" Upward mobility (Q1→Q5): {upward['p_upward']:.3%}")
# Rank-rank scatter
pr = df_mob["parent_income"].rank(pct=True)
cr = df_mob["child_income"].rank(pct=True)
pr_bins = pd.cut(pr, bins=20)
bin_means = pd.DataFrame({"pr": pr, "cr": cr, "bin": pr_bins}).groupby("bin").agg(
pr_mean=("pr", "mean"), cr_mean=("cr", "mean")).reset_index()
fig, ax = plt.subplots(figsize=(8, 6))
ax.scatter(bin_means["pr_mean"], bin_means["cr_mean"], color="#2166ac", s=60, label="Binned means")
slope = rr_result["rr_slope"]
intercept = rr_result["intercept"]
x_line = np.linspace(0, 1, 100)
ax.plot(x_line, intercept + slope * x_line, "r-", linewidth=2,
label=f"Rank-rank slope = {slope:.3f}")
ax.set_xlabel("Parent Income Percentile Rank")
ax.set_ylabel("Child Income Percentile Rank")
ax.set_title("Intergenerational Rank-Rank Mobility")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("rank_rank_mobility.png", dpi=150)
plt.show()
Example 2: Mobility Transition Matrix + Upward Mobility Rate
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
rng = np.random.default_rng(10)
n = 5000
isco_codes = rng.integers(1, 10, n)
df_mob2 = pd.DataFrame({
"parent_isco": isco_codes,
"child_isco": np.clip(isco_codes + rng.integers(-2, 3, n), 1, 9),
})
df_mob2["parent_egp"] = assign_egp(df_mob2["parent_isco"])
df_mob2["child_egp"] = assign_egp(df_mob2["child_isco"])
df_mob2["parent_isei"] = assign_isei(df_mob2["parent_isco"])
df_mob2["child_isei"] = assign_isei(df_mob2["child_isco"])
class_order = [c for c in EGP_ORDER if c in df_mob2["parent_egp"].unique()]
matrix = mobility_transition_matrix(df_mob2["parent_egp"], df_mob2["child_egp"],
class_order=class_order)
print("=== Occupational Transition Matrix (row %) ===")
print(matrix.to_string())
# Heatmap
fig, ax = plt.subplots(figsize=(9, 7))
im = ax.imshow(matrix.values, cmap="YlOrRd", aspect="auto")
ax.set_xticks(range(len(matrix.columns)))
ax.set_yticks(range(len(matrix.index)))
ax.set_xticklabels(matrix.columns.tolist(), rotation=45, ha="right")
ax.set_yticklabels(matrix.index.tolist())
ax.set_xlabel("Destination Class")
ax.set_ylabel("Origin Class")
ax.set_title("Intergenerational Occupational Transition Matrix (%)")
plt.colorbar(im, ax=ax, label="Percentage")
for i in range(len(matrix.index)):
for j in range(len(matrix.columns)):
ax.text(j, i, f"{matrix.iloc[i, j]:.0f}", ha="center", va="center", fontsize=8)
plt.tight_layout()
plt.savefig("mobility_transition_matrix.png", dpi=150)
plt.show()
Example 3: ISEI Distribution by Birth Cohort
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
rng = np.random.default_rng(99)
cohorts = [1950, 1960, 1970, 1980, 1990]
n_per_cohort = 1000
df_isei = pd.concat([
pd.DataFrame({
"cohort": c,
"isco": rng.integers(1, 10, n_per_cohort),
"gender": rng.choice(["Male", "Female"], n_per_cohort),
})
for c in cohorts
], ignore_index=True)
# Later cohorts have higher professional occupations (educational upgrading)
df_isei["isco_adj"] = df_isei.apply(
lambda r: max(1, min(9, int(r["isco"] - (r["cohort"] - 1950) / 20))), axis=1
)
df_isei["isei"] = assign_isei(df_isei["isco_adj"])
fig, ax = plt.subplots(figsize=(10, 5))
for cohort, grp in df_isei.groupby("cohort"):
ax.hist(grp["isei"].dropna(), bins=15, alpha=0.5, label=str(cohort), density=True)
ax.set_xlabel("ISEI Score")
ax.set_ylabel("Density")
ax.set_title("ISEI Occupational Prestige Distribution by Birth Cohort")
ax.legend(title="Birth Cohort")
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("isei_by_cohort.png", dpi=150)
plt.show()
cohort_means = df_isei.groupby("cohort")["isei"].agg(["mean", "std"])
print("=== Mean ISEI by Birth Cohort ===")
print(cohort_means.round(2).to_string())