Source: Beauty/beauty.Rmd
Hamermesh and Parker’s teaching-evaluation data ask a deliberately simple regression question: do instructors rated as more beautiful also receive higher average teaching evaluations? The R original uses rstanarm::stan_glm; the Python version uses the same formula-first Bayesian GLM path through lapylace.
from pathlib import Path
import numpy as np
import pandas as pd
import lapylace as lp
import matplotlib.pyplot as plt
root = Path("../../ROS-Examples")
beauty = pd.read_csv(root / "Beauty/data/beauty.csv")
beauty.head()| eval | beauty | female | age | minority | nonenglish | lower | course_id | |
|---|---|---|---|---|---|---|---|---|
| 0 | 4.3 | 0.201567 | 1 | 36 | 1 | 0 | 0 | 3 |
| 1 | 4.5 | -0.826081 | 0 | 59 | 0 | 0 | 0 | 0 |
| 2 | 3.7 | -0.660333 | 0 | 51 | 0 | 0 | 0 | 4 |
| 3 | 4.3 | -0.766312 | 1 | 40 | 0 | 0 | 0 | 2 |
| 4 | 4.4 | 1.421445 | 1 | 31 | 0 | 0 | 0 | 0 |
def fit_gaussian(formula, seed):
return lp.stan_glm(
formula,
data=beauty,
family=lp.gaussian(),
prior=lp.normal(0, 10),
prior_intercept=lp.normal(0, 10),
prior_aux=lp.exponential(1),
chains=4,
parallel_chains=4,
iter_warmup=500,
iter_sampling=1000,
seed=seed,
refresh=100,
)
def coef_medians(fit):
return pd.Series(
[np.median(fit.alpha_draws()), *np.median(fit.beta_draws(), axis=0)],
index=["Intercept", *fit.columns],
)
def coef_sds(fit):
return pd.Series(
[np.std(fit.alpha_draws()), *np.std(fit.beta_draws(), axis=0)],
index=["Intercept", *fit.columns],
)beauty.describe(include="all")| eval | beauty | female | age | minority | nonenglish | lower | course_id | |
|---|---|---|---|---|---|---|---|---|
| count | 463.000000 | 463.000000 | 463.000000 | 463.000000 | 463.000000 | 463.000000 | 463.000000 | 463.000000 |
| mean | 3.998272 | -0.088349 | 0.421166 | 48.365011 | 0.138229 | 0.060475 | 0.339093 | 4.987041 |
| std | 0.554866 | 0.788648 | 0.494280 | 9.802742 | 0.345513 | 0.238623 | 0.473913 | 8.658995 |
| min | 2.100000 | -1.538843 | 0.000000 | 29.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 |
| 25% | 3.600000 | -0.744618 | 0.000000 | 42.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 |
| 50% | 4.000000 | -0.156363 | 0.000000 | 48.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 |
| 75% | 4.400000 | 0.457253 | 1.000000 | 57.000000 | 0.000000 | 0.000000 | 1.000000 | 6.000000 |
| max | 5.000000 | 1.881674 | 1.000000 | 73.000000 | 1.000000 | 1.000000 | 1.000000 | 30.000000 |
fig, ax = plt.subplots()
ax.scatter(beauty["beauty"], beauty["eval"], alpha=0.75)
ax.set_xlabel("Beauty")
ax.set_ylabel("Average teaching evaluation")Text(0, 0.5, 'Average teaching evaluation')
fit_1 = fit_gaussian("eval ~ beauty", seed=16201)
fit_1.summary(["alpha", "beta", "sigma"]) | Mean | MCSE | StdDev | 5% | 50% | 95% | N_Eff | N_Eff/s | R_hat | |
|---|---|---|---|---|---|---|---|---|---|
| alpha | 4.010180 | 0.000381 | 0.025594 | 3.967090 | 4.010490 | 4.051780 | 4500.960000 | 5295.240000 | 1.000240 |
| beta[1] | 0.133816 | 0.000500 | 0.032340 | 0.081769 | 0.133616 | 0.187690 | 4182.920000 | 4921.080000 | 0.999198 |
| sigma | 0.546355 | 0.000352 | 0.018125 | 0.517281 | 0.545785 | 0.577155 | 2652.029017 | 3120.034138 | 0.999978 |
xs = np.linspace(beauty.beauty.min(), beauty.beauty.max(), 200)
coef_1 = coef_medians(fit_1)
b0, b1 = coef_1["Intercept"], coef_1["beauty"]
sigma = np.median(fit_1.stan_variables()["sigma"])
fig, ax = plt.subplots()
ax.scatter(beauty["beauty"], beauty["eval"], alpha=0.75)
ax.plot(xs, b0 + b1 * xs, color="black")
ax.plot(xs, b0 + b1 * xs + sigma, color="gray", linestyle="--")
ax.plot(xs, b0 + b1 * xs - sigma, color="gray", linestyle="--")
ax.set_xlabel("Beauty")
ax.set_ylabel("Average teaching evaluation")
ax.set_title("Linear fit with +/- one residual SD")Text(0.5, 1.0, 'Linear fit with +/- one residual SD')
The slope is small on the evaluation scale, but positive in this simple regression.
The R page next adds an indicator for women instructors. This gives two parallel regression lines: same beauty slope, different intercepts.
fit_2 = fit_gaussian("eval ~ beauty + female", seed=16202)
fit_2.summary(["alpha", "beta", "sigma"]) | Mean | MCSE | StdDev | 5% | 50% | 95% | N_Eff | N_Eff/s | R_hat | |
|---|---|---|---|---|---|---|---|---|---|
| alpha | 4.094540 | 0.000641 | 0.034279 | 4.036990 | 4.094470 | 4.149790 | 2860.700000 | 3759.130000 | 1.001810 |
| beta[1] | 0.149020 | 0.000569 | 0.032255 | 0.095669 | 0.149517 | 0.201695 | 3211.530000 | 4220.140000 | 0.999814 |
| beta[2] | -0.197132 | 0.000959 | 0.052127 | -0.282705 | -0.197090 | -0.110735 | 2952.240000 | 3879.430000 | 0.999871 |
| sigma | 0.538259 | 0.000297 | 0.017917 | 0.510189 | 0.537229 | 0.568664 | 3634.939584 | 4776.530334 | 1.000758 |
coef_2 = coef_medians(fit_2)
def line_parallel(x, female):
p = coef_2
return p["Intercept"] + p["beauty"] * x + p["female"] * female
fig, axes = plt.subplots(1, 3, figsize=(12, 3.5), sharex=True, sharey=True)
for ax, female, title, color in [(axes[0], 0, "Men", "tab:blue"), (axes[1], 1, "Women", "tab:red")]:
d = beauty[beauty.female == female]
ax.scatter(d["beauty"], d["eval"], alpha=0.75, color=color)
ax.plot(xs, line_parallel(xs, female), color="black")
ax.set_title(title)
ax.set_xlabel("Beauty")
axes[0].set_ylabel("Average teaching evaluation")
axes[2].scatter(beauty.loc[beauty.female == 0, "beauty"], beauty.loc[beauty.female == 0, "eval"], color="tab:blue", alpha=0.65, label="Men")
axes[2].scatter(beauty.loc[beauty.female == 1, "beauty"], beauty.loc[beauty.female == 1, "eval"], color="tab:red", alpha=0.65, label="Women")
axes[2].plot(xs, line_parallel(xs, 0), color="tab:blue")
axes[2].plot(xs, line_parallel(xs, 1), color="tab:red")
axes[2].set_title("Both sexes")
axes[2].set_xlabel("Beauty")
axes[2].legend(frameon=False)
fit_3 = fit_gaussian("eval ~ beauty * female", seed=16203)
fit_3.summary(["alpha", "beta", "sigma"]) | Mean | MCSE | StdDev | 5% | 50% | 95% | N_Eff | N_Eff/s | R_hat | |
|---|---|---|---|---|---|---|---|---|---|
| alpha | 4.104570 | 0.000648 | 0.034140 | 4.047670 | 4.103840 | 4.161350 | 2779.650 | 2988.87000 | 0.999306 |
| beta[1] | 0.201475 | 0.000916 | 0.044307 | 0.128703 | 0.201957 | 0.273947 | 2338.050 | 2514.03000 | 1.000080 |
| beta[2] | -0.205559 | 0.000964 | 0.051030 | -0.291323 | -0.204300 | -0.122626 | 2802.660 | 3013.61000 | 0.999545 |
| beta[3] | -0.113984 | 0.001362 | 0.066373 | -0.224141 | -0.114316 | -0.005010 | 2376.030 | 2554.87000 | 1.000790 |
| sigma | 0.537130 | 0.000280 | 0.018450 | 0.507180 | 0.536560 | 0.568130 | 4304.894 | 4628.91828 | 1.000300 |
coef_3 = coef_medians(fit_3)
def line_interaction(x, female):
p = coef_3
return (
p["Intercept"]
+ p["beauty"] * x
+ p["female"] * female
+ p["beauty:female"] * x * female
)
fig, axes = plt.subplots(1, 2, figsize=(8, 3.5), sharex=True, sharey=True)
for ax, female, title, color in [(axes[0], 0, "Men", "tab:blue"), (axes[1], 1, "Women", "tab:red")]:
d = beauty[beauty.female == female]
ax.scatter(d["beauty"], d["eval"], alpha=0.75, color=color)
ax.plot(xs, line_parallel(xs, female), color="gray", linewidth=1, label="parallel")
ax.plot(xs, line_interaction(xs, female), color="black", linewidth=2, label="interaction")
ax.set_title(title)
ax.set_xlabel("Beauty")
axes[0].set_ylabel("Average teaching evaluation")
axes[0].legend(frameon=False)
The original page fits a sequence of regressions with age, minority status, non-English indicator, lower-division course indicator, and course fixed effects. The same formulas can be passed directly to lapylace.
formulas = {
"beauty + female + age": "eval ~ beauty + female + age",
"beauty + female + minority": "eval ~ beauty + female + minority",
"beauty + female + nonenglish": "eval ~ beauty + female + nonenglish",
"beauty + female + nonenglish + lower": "eval ~ beauty + female + nonenglish + lower",
"beauty + course indicators": "eval ~ beauty + C(course_id)",
}
rows = []
for i, (label, formula) in enumerate(formulas.items(), start=1):
m = fit_gaussian(formula, seed=16210 + i)
coefs = coef_medians(m)
ses = coef_sds(m)
rows.append({
"model": label,
"beauty_coef": coefs.get("beauty", np.nan),
"beauty_sd": ses.get("beauty", np.nan),
"n": len(beauty),
})
pd.DataFrame(rows) | model | beauty_coef | beauty_sd | n | |
|---|---|---|---|---|
| 0 | beauty + female + age | 0.140428 | 0.033940 | 463 |
| 1 | beauty + female + minority | 0.149738 | 0.032677 | 463 |
| 2 | beauty + female + nonenglish | 0.150872 | 0.031848 | 463 |
| 3 | beauty + female + nonenglish + lower | 0.146929 | 0.031100 | 463 |
| 4 | beauty + course indicators | 0.137286 | 0.034266 | 463 |
fit_1.to_arviz()
fit_1.loo()