Use this Skill for political survey analysis: complex sampling with weights, ANES/CCES/ESS data loading, weighted logit/ordered logit, and cross-national equivalence testing.
Use this Skill for political survey analysis: complex sampling with weights, ANES/CCES/ESS data loading, weighted logit/ordered logit, and cross-national equivalence testing.
Load and clean ANES (American National Election Studies), CCES (Cooperative Congressional
Election Study), or ESS (European Social Survey) data files
Compute weighted frequency tables, weighted means, and weighted cross-tabulations
Estimate logit or ordered logit models that account for complex survey designs (stratification,
clustering, probability weights)
Calibrate survey weights through post-stratification or iterative raking
Compare survey measurements across countries and test for cross-national measurement equivalence
Perform Rao-Scott chi-square adjustments for design-based inference
This skill is not a replacement for dedicated survey software (Stata svy, R survey package).
It provides Python implementations suitable for reproducible research workflows.
Background
Political surveys rarely use simple random sampling. The ANES, for example, uses a stratified
multi-stage area probability sample. Ignoring the complex design produces understated standard
errors and invalid inference. Three design features matter:
Feature
Effect if ignored
Probability weights
Biased point estimates
Stratification
Overestimated standard errors
Clustering (PSU)
Underestimated standard errors
Design-based vs. model-based SE: Design-based inference treats the finite population as fixed
and the sample selection as random. Model-based inference conditions on the sample and assumes a
data-generating process. For descriptive inference about populations, design-based SE is preferred.
ANES structure: Each respondent has a weight variable (e.g., V201617x in 2020 ANES). The
pre-election and post-election waves have separate weights. Weights sum to the target population
(eligible voters or adult citizens).
ESS structure: Multi-country survey with a design weight (dweight) correcting for unequal
selection probabilities within countries, and a post-stratification weight (pspwght). For
cross-national analysis, use pweight (population size weight) to make country samples
proportional to national populations.
Raking (iterative proportional fitting): When post-stratification requires simultaneous
calibration on multiple marginal distributions (age × gender × education), raking iterates through
each marginal until convergence. The resulting weights satisfy all marginal totals simultaneously.
Ordered logit for Likert outcomes: Survey items often use 5- or 7-point scales. OLS treats
the ordinal scale as metric; ordered logit respects the ordinal nature and estimates cut-points
between categories.
Measurement equivalence across countries proceeds in steps:
Configural invariance: same factor structure across groups
Metric invariance: equal factor loadings
Scalar invariance: equal item intercepts (required for mean comparison)
import os
import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
from scipy import stats
import warnings
warnings.filterwarnings("ignore", category=FutureWarning)
# ---------------------------------------------------------------------------# 1. Weighted Frequency Tables# ---------------------------------------------------------------------------defweighted_crosstab(
df: pd.DataFrame,
row_var: str,
col_var: str,
weight_var: str,
normalize: str = "row",
) -> pd.DataFrame:
"""
Compute a weighted cross-tabulation.
Parameters
----------
df : pd.DataFrame
row_var : str
Row variable name.
col_var : str
Column variable name.
weight_var : str
Survey weight column.
normalize : str
'row', 'col', or 'all' — passed to pd.crosstab.
Returns
-------
pd.DataFrame
Weighted percentage table.
"""
tab = pd.crosstab(
df[row_var],
df[col_var],
values=df[weight_var],
aggfunc="sum",
normalize=normalize,
)
return (tab * 100).round(1)
defweighted_mean(series: pd.Series, weights: pd.Series) -> float:
"""Compute weighted mean, ignoring NaN in either series."""
mask = series.notna() & weights.notna()
return np.average(series[mask], weights=weights[mask])
defweighted_summary(
df: pd.DataFrame, var: str, weight_var: str, group_var: str | None = None) -> pd.DataFrame:
"""
Weighted mean and std by optional group.
Returns
-------
pd.DataFrame with columns: group (optional), mean, std, n_eff
"""def_stats(sub: pd.DataFrame) -> dict:
w = sub[weight_var].fillna(0)
y = sub[var]
mask = y.notna() & (w > 0)
w, y = w[mask], y[mask]
iflen(w) == 0:
return {"mean": np.nan, "std": np.nan, "n_eff": 0}
mu = np.average(y, weights=w)
var_w = np.average((y - mu) ** 2, weights=w)
n_eff = w.sum() ** 2 / (w ** 2).sum()
return {"mean": round(mu, 4), "std": round(var_w ** 0.5, 4), "n_eff": round(n_eff)}
if group_var isNone:
return pd.DataFrame([_stats(df)])
return df.groupby(group_var).apply(_stats).apply(pd.Series).reset_index()
# ---------------------------------------------------------------------------# 2. Weighted Logit with Survey Weights# ---------------------------------------------------------------------------defweighted_logit(
df: pd.DataFrame,
outcome: str,
predictors: list[str],
weight_var: str,
add_constant: bool = True,
) -> sm.regression.linear_model.RegressionResultsWrapper:
"""
Estimate a logit model using frequency weights as an approximation
to probability-weighted MLE.
Parameters
----------
df : pd.DataFrame
outcome : str
Binary (0/1) dependent variable.
predictors : list of str
weight_var : str
Survey weight column. Weights are scaled to sum to N (sample size)
to preserve degrees of freedom.
add_constant : bool
Whether to add an intercept.
Returns
-------
statsmodels GLMResultsWrapper
"""
sub = df[[outcome] + predictors + [weight_var]].dropna()
y = sub[outcome]
X = sub[predictors]
if add_constant:
X = sm.add_constant(X)
# Scale weights to sum to sample size
w = sub[weight_var]
w_scaled = w / w.mean()
model = sm.GLM(
y,
X,
family=sm.families.Binomial(),
freq_weights=w_scaled,
)
result = model.fit()
return result
deflogit_coeff_table(result) -> pd.DataFrame:
"""
Extract a clean coefficient table with odds ratios.
Returns
-------
pd.DataFrame with columns: coef, se, z, p, OR, OR_lower, OR_upper
"""
tbl = pd.DataFrame({
"coef": result.params,
"se": result.bse,
"z": result.tvalues,
"p": result.pvalues,
})
tbl["OR"] = np.exp(tbl["coef"])
tbl["OR_lower"] = np.exp(tbl["coef"] - 1.96 * tbl["se"])
tbl["OR_upper"] = np.exp(tbl["coef"] + 1.96 * tbl["se"])
return tbl.round(4)
# ---------------------------------------------------------------------------# 3. Ordered Logit for Likert Outcomes# ---------------------------------------------------------------------------defordered_logit(
df: pd.DataFrame,
outcome: str,
predictors: list[str],
weight_var: str | None = None,
) -> object:
"""
Fit an ordered logit (proportional odds) model via statsmodels.
Parameters
----------
outcome : str
Ordinal outcome (integer-coded Likert scale).
predictors : list of str
weight_var : str, optional
If provided, use freq_weights.
Returns
-------
statsmodels OrderedModel result
"""from statsmodels.miscmodels.ordinal_model import OrderedModel
sub = df[[outcome] + predictors + ([weight_var] if weight_var else [])].dropna()
y = sub[outcome].astype(int)
X = sub[predictors]
freq_w = sub[weight_var] / sub[weight_var].mean() if weight_var elseNone
om = OrderedModel(y, X, distr="logit")
result = om.fit(method="bfgs", disp=False)
return result
# ---------------------------------------------------------------------------# 4. Post-stratification Raking# ---------------------------------------------------------------------------defrake_weights(
df: pd.DataFrame,
initial_weight_col: str,
targets: dict[str, dict],
max_iter: int = 50,
tol: float = 1e-6,
) -> pd.Series:
"""
Iterative proportional fitting (raking) to calibrate survey weights.
Parameters
----------
df : pd.DataFrame
initial_weight_col : str
Starting weights (e.g., design weights).
targets : dict
{variable_name: {category_value: target_proportion, ...}}
Example: {'age_group': {1: 0.20, 2: 0.35, 3: 0.30, 4: 0.15}}
max_iter : int
Maximum raking iterations.
tol : float
Convergence tolerance (max relative change in weights).
Returns
-------
pd.Series
Calibrated weights, same index as df.
"""
weights = df[initial_weight_col].copy().astype(float)
for iteration inrange(max_iter):
max_change = 0.0for var, target_props in targets.items():
for cat, target_prop in target_props.items():
mask = df[var] == cat
current_share = weights[mask].sum() / weights.sum()
if current_share > 0:
adjustment = target_prop / current_share
old_w = weights[mask].copy()
weights[mask] *= adjustment
change = np.abs(weights[mask] - old_w).max() / (old_w.max() + 1e-12)
max_change = max(max_change, change)
if max_change < tol:
print(f"Raking converged in {iteration + 1} iterations.")
breakelse:
print(f"Warning: raking did not converge after {max_iter} iterations.")
# Normalize to original total
weights *= df[initial_weight_col].sum() / weights.sum()
return weights
# ---------------------------------------------------------------------------# 5. Rao-Scott Chi-Square Adjustment# ---------------------------------------------------------------------------defrao_scott_chisq(observed: np.ndarray, weights: np.ndarray) -> dict:
"""
First-order Rao-Scott chi-square adjustment for design effect.
Parameters
----------
observed : np.ndarray
2D contingency table (raw counts).
weights : np.ndarray
1D array of weights for each respondent in the table.
Returns
-------
dict with keys: chisq_rs, df, pvalue, deff
"""from scipy.stats import chi2
# Unweighted Pearson chi-sq
row_totals = observed.sum(axis=1)
col_totals = observed.sum(axis=0)
total = observed.sum()
expected = np.outer(row_totals, col_totals) / total
chisq_pearson = ((observed - expected) ** 2 / expected).sum()
# Design effect approximation
n = weights.sum()
deff = (weights ** 2).sum() * n / (weights.sum() ** 2)
chisq_rs = chisq_pearson / deff
df = (observed.shape[0] - 1) * (observed.shape[1] - 1)
pvalue = 1 - chi2.cdf(chisq_rs, df)
return {"chisq_rs": round(chisq_rs, 4), "df": df, "pvalue": round(pvalue, 4), "deff": round(deff, 4)}
Advanced Usage
Cross-National ESS Analysis
The ESS runs every two years across 20+ European countries. Country-level weights (pspwght)
correct for within-country stratification and non-response. The cross-national weight (pweight)
makes country sample sizes proportional to population size, enabling continent-wide estimates.
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
defload_ess(path: str, variables: list[str], countries: list[str] | None = None) -> pd.DataFrame:
"""
Load ESS data from Stata .dta file.
Parameters
----------
path : str
Path to ESS .dta file.
variables : list of str
Variables to retain plus essential columns.
countries : list of str, optional
Filter by cntry (ISO2 country code).
Returns
-------
pd.DataFrame
"""
essential = ["cntry", "idno", "dweight", "pspwght", "pweight"]
keep = list(set(essential + variables))
df = pd.read_stata(path, columns=[c for c in keep], convert_categoricals=False)
if countries:
df = df[df["cntry"].isin(countries)]
# Combined analysis weight = pspwght * pweight
df["analysis_weight"] = df["pspwght"] * df["pweight"]
return df
defess_country_means(
df: pd.DataFrame, var: str, weight_col: str = "pspwght") -> pd.DataFrame:
"""Compute weighted country means for a variable."""
results = []
for country, grp in df.groupby("cntry"):
w = grp[weight_col]
y = grp[var]
mask = y.notna() & w.notna() & (w > 0)
if mask.sum() < 10:
continue
mu = np.average(y[mask], weights=w[mask])
n = mask.sum()
se = y[mask].std() / np.sqrt(n)
results.append({"country": country, "mean": mu, "se": se, "n": n})
return pd.DataFrame(results).sort_values("mean", ascending=False)
defplot_country_means(means_df: pd.DataFrame, var_label: str, save_path: str | None = None):
"""Horizontal bar chart of country-level weighted means with 95% CI."""
df = means_df.sort_values("mean")
fig, ax = plt.subplots(figsize=(9, max(5, len(df) * 0.35)))
y_pos = range(len(df))
ax.barh(y_pos, df["mean"], xerr=1.96 * df["se"], align="center",
color="#4c72b0", ecolor="#c44e52", capsize=3, alpha=0.85)
ax.set_yticks(list(y_pos))
ax.set_yticklabels(df["country"].tolist())
ax.set_xlabel(var_label)
ax.set_title(f"Country-Level Weighted Means: {var_label}")
ax.grid(True, axis="x", alpha=0.3)
plt.tight_layout()
if save_path:
fig.savefig(save_path, dpi=150)
return fig
# Example: Trust in parliament across ESS Round 10 countries
ESS_PATH = os.environ.get("ESS_PATH", "ESS10.dta")
# df_ess = load_ess(ESS_PATH, variables=["trstprl", "trstplt", "age", "eduyrs", "gndr"])# means = ess_country_means(df_ess, "trstprl")# print(means.head(10).to_string(index=False))# plot_country_means(means, "Trust in Parliament (0-10)", save_path="trust_parliament.png")# Raking example: calibrate ANES weights to Census targets
rake_targets = {
"age_group": {1: 0.15, 2: 0.20, 3: 0.25, 4: 0.22, 5: 0.18},
"gender": {1: 0.49, 2: 0.51},
"educ3": {1: 0.28, 2: 0.38, 3: 0.34},
}
# Assuming df_anes has columns: age_group, gender, educ3, base_weight# df_anes["raked_weight"] = rake_weights(df_anes, "base_weight", rake_targets)# Verify margins after raking:# for var, targets in rake_targets.items():# w = df_anes["raked_weight"]# for cat, tgt in targets.items():# actual = w[df_anes[var] == cat].sum() / w.sum()# print(f"{var}={cat}: target={tgt:.3f}, actual={actual:.3f}")
Troubleshooting
Problem
Cause
Solution
KeyError on weight variable
Different weight names across ANES years
Inspect codebook; 2020 ANES pre-election weight is V201617x
Raking does not converge
Conflicting marginal targets or zero cells
Check that targets sum to 1.0 per variable; increase max_iter
OrderedModel fails
Outcome not integer-coded
Cast with .astype(int) after mapping categories
Design effect >> 3
High clustering in PSUs
Consider explicit cluster SE using cov_kwds={'groups': psu_col}
ESS pweight missing
Country not in cross-national file
Download the integrated ESS file, not country-specific files
Weighted logit perfect separation
Sparse cells after weighting
Regularize with alpha in fit_regularized() or collapse categories
Lumley, T. (2010). Complex Surveys: A Guide to Analysis Using R. Wiley.
Pasek, J. (2018). anesrake: ANES Raking Implementation. CRAN.
Lehtonen, R. & Pahkinen, E. (2004). Practical Methods for Design and Analysis of Complex Surveys.
Rao, J.N.K. & Scott, A.J. (1981). The analysis of categorical data from complex sample surveys.
Journal of the American Statistical Association, 76(374), 221-230.
Examples
Example 1: Weighted Logit — Vote Choice in ANES 2020
