Categorical Data Analysis for Complex Surveys in Python

Design-aware tabulations, cross-tabulations, hypothesis tests, and rank tests

Tutorials

Categorical Analysis

Hypothesis Testing

Python

Analyze categorical survey data in Python with design-adjusted tabulations, cross-tabulations, hypothesis tests, and nonparametric rank tests. Learn to create weighted contingency tables, t-tests, and Wilcoxon / Kruskal-Wallis rank tests using the svy library.

Author

Mamadou S. Diallo, Ph.D.

Published

January 18, 2026

Modified

July 21, 2026

Keywords

categorical data analysis survey Python, survey cross-tabulation Python, weighted contingency table Python, design-adjusted t-test survey Python, weighted frequency table Python, survey tabulation Python, two-sample t-test survey Python, chi-square test complex survey Python, domain estimation survey Python, complex survey hypothesis testing Python, survey rank test Python, Wilcoxon rank test complex survey Python, Kruskal-Wallis test survey Python, design-based nonparametric test Python

Categorical data analysis spans descriptive techniques—contingency tables and cross-tabulations with design-adjusted tests—and model-based approaches for categorical outcomes (logistic, multinomial, loglinear, and mixed-effects GLMs; see Agresti (2013)).

In complex survey applications, ignoring stratification, clustering, or unequal weights can misstate uncertainty. This tutorial shows how to use svy to produce design-aware tabulations, cross-tabulations, t-tests for group differences, and nonparametric rank tests, with standard errors that respect the sample design.

Setting Up the Sample

We’ll use the imaginary country household dataset from World Bank (2023):

import numpy as np
import svy

hld_data = svy.datasets.load(name="hld_sample_wb_2023", source="bundled")

hld_design = svy.Design(stratum=("geo1", "urbrur"), psu="ea", wgt="hhweight")
hld_sample = svy.Sample(data=hld_data, design=hld_design)

# Create household poverty line and binary poverty status
hld_sample = hld_sample.wrangling.mutate(
    {
        "hhpovline": svy.col("hhsize") * 1800,
        "pov_status": svy.when(svy.col("tot_exp") < svy.col("hhpovline")).then(1).otherwise(0),
    }
)

One-Way Tabulation

The tabulate() method produces weighted frequency tables that account for the survey design.

Tabulate the first-level administrative unit (geo1):

hld_admin1_tab = hld_sample.categorical.tabulate(rowvar="geo1")
print(hld_admin1_tab)

╭────────────────────── Table: geo1 ───────────────────────╮
│                                                          │
│  Row      Estimate   Std Err       CV    Lower    Upper  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  geo_01     0.1710    0.0000   0.0000   0.1710   0.1710  │
│  geo_02     0.2679    0.0000   0.0000   0.2679   0.2679  │
│  geo_03     0.2513    0.0000   0.0000   0.2513   0.2513  │
│  geo_04     0.3098    0.0000   0.0000   0.3098   0.3098  │
│                                                          │
╰──────────────────────────────────────────────────────────╯

Changing Output Units

By default, tabulate() produces proportions. Use the units parameter to get counts or percentages:

# Counts
hld_admin1_tab_count = hld_sample.categorical.tabulate(
    rowvar="geo1",
    units="count",
)
print("Table with counts")
print(hld_admin1_tab_count)

Table with counts

╭─────────────────────────── Table: geo1 ────────────────────────────╮
│                                                                    │
│  Row        Estimate   Std Err       CV        Lower        Upper  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  geo_01    9896.0000    0.0000   0.0000    9896.0000    9896.0000  │
│  geo_02   15502.0000    0.0000   0.0000   15502.0000   15502.0000  │
│  geo_03   14541.0000    0.0000   0.0000   14541.0000   14541.0000  │
│  geo_04   17928.0000    0.0000   0.0000   17928.0000   17928.0000  │
│                                                                    │
╰────────────────────────────────────────────────────────────────────╯

# Percentages
hld_admin1_tab_percent = hld_sample.categorical.tabulate(
    rowvar="geo1",
    units="percent",
)
print("Table with percentages:")
print(hld_admin1_tab_percent)

Table with percentages:

╭─────────────────────── Table: geo1 ────────────────────────╮
│                                                            │
│  Row      Estimate   Std Err       CV     Lower     Upper  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  geo_01    17.1013    0.0000   0.0000   17.1013   17.1013  │
│  geo_02    26.7890    0.0000   0.0000   26.7890   26.7890  │
│  geo_03    25.1283    0.0000   0.0000   25.1283   25.1283  │
│  geo_04    30.9814    0.0000   0.0000   30.9814   30.9814  │
│                                                            │
╰────────────────────────────────────────────────────────────╯

Scaling Counts to a Custom Total

Use count_total to express counts on an arbitrary total (useful for scaled headcounts while preserving shares):

# Scale counts so the total sums to 1,000
hld_admin1_tab_n = hld_sample.categorical.tabulate(
    rowvar="geo1",
    count_total=1_000,
)
print(hld_admin1_tab_n)

╭──────────────────────── Table: geo1 ─────────────────────────╮
│                                                              │
│  Row      Estimate   Std Err       CV      Lower      Upper  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  geo_01   171.0128    0.0000   0.0000   171.0128   171.0128  │
│  geo_02   267.8902    0.0000   0.0000   267.8902   267.8902  │
│  geo_03   251.2831    0.0000   0.0000   251.2831   251.2831  │
│  geo_04   309.8139    0.0000   0.0000   309.8139   309.8139  │
│                                                              │
╰──────────────────────────────────────────────────────────────╯

Two-Way Tabulation (Cross-Tabulation)

Cross-tabulations examine the relationship between two categorical variables. Use the colvar parameter to create a two-way table:

# Cross-tabulate urban/rural status by electricity access
urbrur_elec_tab = hld_sample.categorical.tabulate(
    rowvar="urbrur",
    colvar="electricity",
)
print(urbrur_elec_tab)

╭───────────────── Table: urbrur × electricity ─────────────────╮
│                                                               │
│  Row     Col   Estimate   Std Err       CV    Lower    Upper  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  Rural   No      0.0945    0.0186   0.1968   0.0626   0.1403  │
│  Rural   Yes     0.2432    0.0186   0.0765   0.2070   0.2835  │
│  Urban   No      0.0236    0.0145   0.6164   0.0066   0.0813  │
│  Urban   Yes     0.6386    0.0145   0.0228   0.6082   0.6680  │
│                                                               │
╰───────────────────────────────────────────────────────────────╯

Viewing as a Crosstab Matrix

The crosstab() method returns a Polars DataFrame in a familiar matrix format:

# Get crosstab matrix (estimates only)
urbrur_elec_crosstab = urbrur_elec_tab.crosstab()
print(urbrur_elec_crosstab)

shape: (2, 3)
┌────────┬──────────┬──────────┐
│ urbrur ┆ No       ┆ Yes      │
│ ---    ┆ ---      ┆ ---      │
│ str    ┆ f64      ┆ f64      │
╞════════╪══════════╪══════════╡
│ Rural  ┆ 0.094527 ┆ 0.24323  │
│ Urban  ┆ 0.023598 ┆ 0.638645 │
└────────┴──────────┴──────────┘

# Include standard errors with the estimates
urbrur_elec_crosstab_se = urbrur_elec_tab.crosstab(stats=("est", "se"))
print(urbrur_elec_crosstab_se)

shape: (2, 3)
┌────────┬───────────────┬───────────────┐
│ urbrur ┆ No            ┆ Yes           │
│ ---    ┆ ---           ┆ ---           │
│ str    ┆ str           ┆ str           │
╞════════╪═══════════════╪═══════════════╡
│ Rural  ┆ 0.095 ± 0.019 ┆ 0.243 ± 0.019 │
│ Urban  ┆ 0.024 ± 0.015 ┆ 0.639 ± 0.015 │
└────────┴───────────────┴───────────────┘

Two-Way Tables with Different Units

Just like one-way tables, you can change the output units:

# Cross-tabulation with percentages
urbrur_elec_pct = hld_sample.categorical.tabulate(
    rowvar="urbrur",
    colvar="electricity",
    units="percent",
)
print("Cross-tabulation with percentages:")
print(urbrur_elec_pct.crosstab())

Cross-tabulation with percentages:
shape: (2, 3)
┌────────┬──────────┬───────────┐
│ urbrur ┆ No       ┆ Yes       │
│ ---    ┆ ---      ┆ ---       │
│ str    ┆ f64      ┆ f64       │
╞════════╪══════════╪═══════════╡
│ Rural  ┆ 9.452745 ┆ 24.322982 │
│ Urban  ┆ 2.359756 ┆ 63.864517 │
└────────┴──────────┴───────────┘

# Cross-tabulation with counts
urbrur_elec_count = hld_sample.categorical.tabulate(
    rowvar="urbrur",
    colvar="electricity",
    units="count",
)
print("Cross-tabulation with counts:")
print(urbrur_elec_count.crosstab())

Cross-tabulation with counts:
shape: (2, 3)
┌────────┬─────────┬──────────┐
│ urbrur ┆ No      ┆ Yes      │
│ ---    ┆ ---     ┆ ---      │
│ str    ┆ f64     ┆ f64      │
╞════════╪═════════╪══════════╡
│ Rural  ┆ 5470.02 ┆ 14074.98 │
│ Urban  ┆ 1365.52 ┆ 36956.48 │
└────────┴─────────┴──────────┘

Crosstab Display Options

The crosstab() method offers several formatting options:

# Show estimates with confidence intervals
print(urbrur_elec_tab.crosstab(stats=("est", "lci", "uci"), precision=4))

shape: (2, 3)
┌────────┬─────────────────────────┬─────────────────────────┐
│ urbrur ┆ No                      ┆ Yes                     │
│ ---    ┆ ---                     ┆ ---                     │
│ str    ┆ str                     ┆ str                     │
╞════════╪═════════════════════════╪═════════════════════════╡
│ Rural  ┆ 0.0945 [0.0626, 0.1403] ┆ 0.2432 [0.2070, 0.2835] │
│ Urban  ┆ 0.0236 [0.0066, 0.0813] ┆ 0.6386 [0.6082, 0.6680] │
└────────┴─────────────────────────┴─────────────────────────┘

# Export full table details to Polars DataFrame
df = urbrur_elec_tab.to_polars()
print(df)

shape: (4, 8)
┌────────┬─────────────┬──────────┬──────────┬──────────┬──────────┬────────────┬───────┐
│ urbrur ┆ electricity ┆ est      ┆ se       ┆ lci      ┆ uci      ┆ table_type ┆ alpha │
│ ---    ┆ ---         ┆ ---      ┆ ---      ┆ ---      ┆ ---      ┆ ---        ┆ ---   │
│ str    ┆ str         ┆ f64      ┆ f64      ┆ f64      ┆ f64      ┆ str        ┆ f64   │
╞════════╪═════════════╪══════════╪══════════╪══════════╪══════════╪════════════╪═══════╡
│ Rural  ┆ No          ┆ 0.094527 ┆ 0.018606 ┆ 0.062597 ┆ 0.140308 ┆ Two-Way    ┆ 0.05  │
│ Rural  ┆ Yes         ┆ 0.24323  ┆ 0.018606 ┆ 0.207045 ┆ 0.283478 ┆ Two-Way    ┆ 0.05  │
│ Urban  ┆ No          ┆ 0.023598 ┆ 0.014546 ┆ 0.006559 ┆ 0.081282 ┆ Two-Way    ┆ 0.05  │
│ Urban  ┆ Yes         ┆ 0.638645 ┆ 0.014546 ┆ 0.608242 ┆ 0.667977 ┆ Two-Way    ┆ 0.05  │
└────────┴─────────────┴──────────┴──────────┴──────────┴──────────┴────────────┴───────┘

T-Tests for Survey Data

The ttest() method performs design-adjusted t-tests that properly account for stratification, clustering, and weighting.

One-Sample T-Test

Test whether a population mean equals a hypothesized value:

# Test: Is the poverty rate different from 25%?
pov_status_mean_h0 = hld_sample.categorical.ttest(
    y="pov_status",
    mean_h0=0.25,
)
print(pov_status_mean_h0)

╭────────────── T-Test: One-sample ───────────────╮
│ Y = 'pov_status'                                │
│ H₀: μ = 0.2500                                  │
│                                                 │
│                                                 │
│  Estimate   Std Err       CV    Lower    Upper  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│    0.2064    0.0389   0.1886   0.1262   0.2866  │
│                                                 │
│                                                 │
│ Test statistic                                  │
│                                                 │
│     diff         t        df   p_value          │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━          │
│  -0.0436   -1.1189   25.0000    0.2738          │
│                                                 │
╰─────────────────────────────────────────────────╯

Accessing Results Programmatically

# Test statistics
print(f"t-statistic: {pov_status_mean_h0.stats.t:.4f}")
print(f"Degrees of freedom: {pov_status_mean_h0.stats.df}")
print(f"p-value: {pov_status_mean_h0.stats.p_value:.4f}")

t-statistic: -1.1189
Degrees of freedom: 25.0
p-value: 0.2738

# Difference from hypothesized mean
diff = pov_status_mean_h0.diff[0]
print(f"Difference: {diff.diff:.4f}")
print(f"95% CI: [{diff.lci:.4f}, {diff.uci:.4f}]")

Difference: -0.0436
95% CI: [-0.1238, 0.0366]

Export to DataFrame

# Combined test results (default)
pov_status_mean_h0.to_polars()

shape: (1, 8)

y	diff	se	lci	uci	t	df	p_value
str	f64	f64	f64	f64	f64	f64	f64
"pov_status"	-0.043567	0.038938	-0.123761	0.036627	-1.118891	25.0	0.273823

# Raw estimates only
pov_status_mean_h0.to_polars("estimates")

shape: (1, 5)

est	se	cv	lci	uci
f64	f64	f64	f64	f64
0.206433	0.038938	0.188623	0.126239	0.286627

Two-Sample T-Test

Compare means between two groups:

# Test: Does total expenditure differ between urban and rural areas?
exp_by_urbrur = hld_sample.categorical.ttest(
    y="tot_exp",
    group="urbrur",
)
print(exp_by_urbrur)

╭─────────────────────── T-Test: Two-sample (unpaired) ────────────────────────╮
│ Y = 'tot_exp'                                                                │
│ Groups: urbrur = ['Rural' vs 'Urban']                                        │
│                                                                              │
│                                                                              │
│  Group    Level     Estimate     Std Err       CV        Lower        Upper  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  urbrur   Rural    8059.2760    962.8024   0.1195    6076.3473   10042.2048  │
│  urbrur   Urban   14573.2217   1011.0039   0.0694   12491.0201   16655.4233  │
│                                                                              │
│                                                                              │
│ Test statistic                                                               │
│                                                                              │
│       diff        t        df   p_value                                      │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━                                      │
│  6513.9457   4.6658   25.0000   <0.0001                                      │
│                                                                              │
╰──────────────────────────────────────────────────────────────────────────────╯

Two-Sample Results

# Group difference
diff = exp_by_urbrur.diff[0]
print(f"Difference: {diff.diff:.2f}")
print(f"95% CI: [{diff.lci:.2f}, {diff.uci:.2f}]")

Difference: 6513.95
95% CI: [3638.61, 9389.28]

# Individual group estimates
for est in exp_by_urbrur.estimates:
    print(f"{est.group_level}: {est.est:.2f} (SE: {est.se:.2f})")

Rural: 8059.28 (SE: 962.80)
Urban: 14573.22 (SE: 1011.00)

# Export to DataFrame
exp_by_urbrur.to_polars()

shape: (1, 10)

y	group_var	paired	diff	se	lci	uci	t	df	p_value
str	str	bool	f64	f64	f64	f64	f64	f64	f64
"tot_exp"	"urbrur"	false	6513.945685	1396.107996	3638.607443	9389.283927	4.665789	25.0	0.000089

Alternative Hypotheses

Specify one-sided tests with the alternative parameter:

# Test: Is urban expenditure GREATER than rural expenditure?
exp_by_urbrur_greater = hld_sample.categorical.ttest(
    y="tot_exp",
    group="urbrur",
    alternative="greater",
)
print(f"One-sided p-value: {exp_by_urbrur_greater.stats.p_value:.6f}")
print(f"Alternative: {exp_by_urbrur_greater.alternative}")

One-sided p-value: 0.000044
Alternative: greater

Options:

`alternative`	Null Hypothesis	Alternative Hypothesis
`"two-sided"`	μ₁ = μ₂	μ₁ ≠ μ₂
`"less"`	μ₁ ≥ μ₂	μ₁ < μ₂
`"greater"`	μ₁ ≤ μ₂	μ₁ > μ₂

For one-sample tests, the comparison is against mean_h0.

Domain Estimation with `by`

Perform separate t-tests for each level of a domain variable:

# Test poverty rate vs 25% separately by region
pov_by_region = hld_sample.categorical.ttest(
    y="pov_status",
    mean_h0=0.25,
    by="geo1",
)

# Returns a list of TTestOneGroup objects
for r in pov_by_region:
    diff = r.diff[0]
    print(f"Region {diff.by_level}: diff={diff.diff:.4f}, p={r.stats.p_value:.4f}")

Region geo_01: diff=-0.0980, p=0.3163
Region geo_02: diff=-0.0771, p=0.2264
Region geo_03: diff=0.0124, p=0.8993
Region geo_04: diff=-0.0300, p=0.6347

Rank Tests for Survey Data

When comparing distributions across groups, the classic Wilcoxon rank-sum and Kruskal–Wallis tests assume independent observations with equal selection probabilities—assumptions that break down under complex sampling. The ranktest() method implements the design-based rank tests of Lumley and Scott (2013), which replace raw ranks with estimated population mid-ranks computed from the survey weights and then apply a rank-score transformation. Standard errors come from the same Taylor-linearization machinery used elsewhere in svy, so stratification, clustering, and unequal weights are handled automatically.

Wilcoxon vs. Kruskal–Wallis

Under the hood, both tests use the same score function ("kruskal-wallis"). The distinction is purely about the number of groups: with exactly two groups svy runs a Wilcoxon-style t-test on the difference in mean rank scores; with three or more groups it runs a Kruskal–Wallis-style Wald / F-test on the full set of contrasts. You do not need to choose—ranktest() detects the number of groups and picks the right test automatically.

Two-Sample Rank Test (Wilcoxon)

Compare the total expenditure distribution between urban and rural households:

# Two-sample Wilcoxon rank test: expenditure by urban/rural
exp_rank_urbrur = hld_sample.categorical.ranktest(
    "tot_exp",
    group="urbrur",
    method="kruskal-wallis",
    drop_nulls=True,
)
print(exp_rank_urbrur)

╭──── Rank Test: Two-sample (Wilcoxon) ─────╮
│ Y = 'tot_exp'                             │
│ Groups: urbrur = ['Rural' vs 'Urban']     │
│                                           │
│                                           │
│  Estimate   Std Err   CV   Lower   Upper  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│                                           │
│                                           │
│ Test statistic                            │
│                                           │
│    diff        t        df   p_value      │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━      │
│  0.2991   5.0586   25.0000   <0.0001      │
│                                           │
╰───────────────────────────────────────────╯

Even though we passed "kruskal-wallis", the result reports itself as “Wilcoxon” because there are only two groups:

print(f"Method displayed: {exp_rank_urbrur.method_name}")

Method displayed: Wilcoxon

Accessing Two-Sample Results

# Test statistics
print(f"t-statistic: {exp_rank_urbrur.stats.value:.4f}")
print(f"Degrees of freedom: {exp_rank_urbrur.stats.df}")
print(f"p-value: {exp_rank_urbrur.stats.p_value:.4f}")

t-statistic: 5.0586
Degrees of freedom: 25.0
p-value: 0.0000

# Difference in mean rank score
diff = exp_rank_urbrur.diff[0]
print(f"Difference: {diff.diff:.4f}")
print(f"95% CI: [{diff.lci:.4f}, {diff.uci:.4f}]")

Difference: 0.2991
95% CI: [0.1773, 0.4209]

Export to DataFrame

# Combined test output (default)
exp_rank_urbrur.to_polars()

shape: (1, 10)

y	group_var	method	diff	se	lci	uci	t	df	p_value
str	str	str	f64	f64	f64	f64	f64	f64	f64
"tot_exp"	"urbrur"	"Wilcoxon"	0.299125	0.059132	0.177341	0.420909	5.058626	25.0	0.000032

# Group-level estimates only
exp_rank_urbrur.to_polars("estimates")

shape: (0, 8)

y	group	group_level	est	se	cv	lci	uci
str	str	str	f64	f64	f64	f64	f64

K-Sample Rank Test (Kruskal–Wallis)

When the grouping variable has three or more levels, ranktest() automatically performs a Kruskal–Wallis F-test:

# K-sample rank test: total expenditure by region
exp_rank_region = hld_sample.categorical.ranktest(
    "tot_exp",
    group="geo1",
    method="kruskal-wallis",
    drop_nulls=True,
)
print(exp_rank_region)

╭── Rank Test: K-sample (Kruskal-Wallis) ───╮
│ Y = 'tot_exp'                             │
│ Groups: geo1 (4 levels)                   │
│                                           │
│                                           │
│  Estimate   Std Err   CV   Lower   Upper  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│                                           │
│                                           │
│ Test statistic                            │
│                                           │
│  df     Chisq        F   p_value          │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━          │
│   3   27.1816   9.0605    0.0004          │
│                                           │
╰───────────────────────────────────────────╯

print(f"Method displayed: {exp_rank_region.method_name}")
print(f"F-statistic: {exp_rank_region.stats.value:.4f}")
print(f"Numerator df: {int(exp_rank_region.stats.df_num)}")
print(f"Denominator df: {exp_rank_region.stats.df_den:.0f}")
print(f"p-value: {exp_rank_region.stats.p_value:.4f}")

Method displayed: Kruskal-Wallis
F-statistic: 9.0605
Numerator df: 3
Denominator df: 23
p-value: 0.0004

# Export to DataFrame
exp_rank_region.to_polars()

shape: (1, 8)

y	group_var	method	ndf	ddf	chisq	f_stat	p_value
str	str	str	i64	f64	f64	f64	f64
"tot_exp"	"geo1"	"Kruskal-Wallis"	3	23.0	27.181565	9.060522	0.000381

Alternative Score Functions

Besides the default Wilcoxon / Kruskal–Wallis scores, svy supports two additional built-in score transformations. Each applies a different function g(r) to the estimated population mid-ranks r, with N denoting the estimated population total:

Method	Score	Sensitive to
`"kruskal-wallis"`	g(r) = r / N	General location shift
`"vander-waerden"`	g(r) = Φ⁻¹(r / N)	Location shift (emphasises tails)
`"median"`	g(r) = I(r > N / 2)	Difference in medians

Van der Waerden Scores

Van der Waerden scores transform the proportional ranks through the standard-normal quantile function. This gives more weight to observations in the tails and tends to be more powerful than Wilcoxon scores when the underlying distribution is close to normal:

# Van der Waerden rank test
exp_vdw = hld_sample.categorical.ranktest(
    "tot_exp",
    group="urbrur",
    method="vander-waerden",
    drop_nulls=True,
)
print(exp_vdw)

╭─ Rank Test: Two-sample (van der Waerden) ─╮
│ Y = 'tot_exp'                             │
│ Groups: urbrur = ['Rural' vs 'Urban']     │
│                                           │
│                                           │
│  Estimate   Std Err   CV   Lower   Upper  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│                                           │
│                                           │
│ Test statistic                            │
│                                           │
│    diff        t        df   p_value      │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━      │
│  1.0241   4.7891   25.0000   <0.0001      │
│                                           │
╰───────────────────────────────────────────╯

exp_vdw.to_polars()

shape: (1, 10)

y	group_var	method	diff	se	lci	uci	t	df	p_value
str	str	str	f64	f64	f64	f64	f64	f64	f64
"tot_exp"	"urbrur"	"van der Waerden"	1.024072	0.213832	0.583677	1.464467	4.789145	25.0	0.000064

Mood’s Median Test

Mood’s median test reduces each observation to a binary indicator of whether its rank exceeds the population median. It is robust to outliers but less powerful than Wilcoxon or van der Waerden for detecting location shifts:

# Mood's median rank test
exp_median = hld_sample.categorical.ranktest(
    "tot_exp",
    group="urbrur",
    method="median",
    drop_nulls=True,
)
print(exp_median)

╭───── Rank Test: Two-sample (Median) ──────╮
│ Y = 'tot_exp'                             │
│ Groups: urbrur = ['Rural' vs 'Urban']     │
│                                           │
│                                           │
│  Estimate   Std Err   CV   Lower   Upper  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│                                           │
│                                           │
│ Test statistic                            │
│                                           │
│    diff        t        df   p_value      │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━      │
│  0.4313   5.1862   25.0000   <0.0001      │
│                                           │
╰───────────────────────────────────────────╯

exp_median.to_polars()

shape: (1, 10)

y	group_var	method	diff	se	lci	uci	t	df	p_value
str	str	str	f64	f64	f64	f64	f64	f64	f64
"tot_exp"	"urbrur"	"Median"	0.431315	0.083166	0.260032	0.602599	5.186194	25.0	0.000023

Custom Score Functions

You can supply your own rank-score transformation via the score_fn parameter. The function must accept (r, N) where r is the array of estimated population mid-ranks and N is the estimated population total, and return an array of scores. When using score_fn, omit the method parameter—providing both raises an error.

import numpy as np

# Custom: upper-quartile indicator
def upper_quartile(r, N):
    return (r > 0.75 * N).astype(np.float64)

exp_custom = hld_sample.categorical.ranktest(
    "tot_exp",
    group="urbrur",
    score_fn=upper_quartile,
    drop_nulls=True,
)
print(exp_custom)

╭─ Rank Test: Two-sample (upper_quartile) ──╮
│ Y = 'tot_exp'                             │
│ Groups: urbrur = ['Rural' vs 'Urban']     │
│                                           │
│                                           │
│  Estimate   Std Err   CV   Lower   Upper  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│                                           │
│                                           │
│ Test statistic                            │
│                                           │
│    diff        t        df   p_value      │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━      │
│  0.2566   3.8704   25.0000    0.0007      │
│                                           │
╰───────────────────────────────────────────╯

The method name is inferred from the function’s __name__:

print(f"Method displayed: {exp_custom.method_name}")

Method displayed: upper_quartile

exp_custom.to_polars()

shape: (1, 10)

y	group_var	method	diff	se	lci	uci	t	df	p_value
str	str	str	f64	f64	f64	f64	f64	f64	f64
"tot_exp"	"urbrur"	"upper_quartile"	0.256637	0.066307	0.120075	0.393199	3.870438	25.0	0.000691

One-Sided Rank Tests

For two-sample tests, you can specify alternative="less" or alternative="greater":

# One-sided: is urban rank-score distribution shifted above rural?
exp_rank_greater = hld_sample.categorical.ranktest(
    "tot_exp",
    group="urbrur",
    method="kruskal-wallis",
    alternative="greater",
    drop_nulls=True,
)
print(f"One-sided p-value: {exp_rank_greater.stats.p_value:.6f}")
print(f"Alternative: {exp_rank_greater.alternative}")

One-sided p-value: 0.000016
Alternative: greater

Degrees of Freedom

For complex surveys, degrees of freedom depend on the design:

Design	Degrees of Freedom
Weights only (SRS)	n - 1
Stratified	n_strata - 1
Clustered	n_psu - 1
Stratified + Clustered	n_psu - n_strata - 1

Where:

n = number of observations
n_strata = number of strata
n_psu = number of primary sampling units

Next Steps

Continue to Generalized Linear Models to learn how to fit linear and logistic regression models with design-adjusted standard errors.

Ready for regression modeling?
Learn GLMs in Generalized Linear Models →

References

Agresti, Alan. 2013. Categorical Data Analysis, 3rd edn. John Wiley & Sons, Hoboken, New Jersey.

Lumley, Thomas, and Alastair J. Scott. 2013. “Two-Sample Rank Tests Under Complex Sampling.” Biometrika 100 (4): 831–42. https://doi.org/10.1093/biomet/ast027.

World Bank. 2023. “Synthetic Data for an Imaginary Country, Sample, 2023.” World Bank, Development Data Group. https://doi.org/10.48529/MC1F-QH23.

Setting Up the Sample

One-Way Tabulation

Changing Output Units

Scaling Counts to a Custom Total

Two-Way Tabulation (Cross-Tabulation)

Viewing as a Crosstab Matrix

Two-Way Tables with Different Units

Crosstab Display Options

T-Tests for Survey Data

One-Sample T-Test

Accessing Results Programmatically

Export to DataFrame

Two-Sample T-Test

Two-Sample Results

Alternative Hypotheses

Domain Estimation with by

Rank Tests for Survey Data

Wilcoxon vs. Kruskal–Wallis

Two-Sample Rank Test (Wilcoxon)

Accessing Two-Sample Results

Export to DataFrame

K-Sample Rank Test (Kruskal–Wallis)

Alternative Score Functions

Van der Waerden Scores

Mood’s Median Test

Custom Score Functions

One-Sided Rank Tests

Degrees of Freedom

Next Steps

References

Domain Estimation with `by`