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

April 18, 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")

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.1370    0.0041   0.0296   0.1292   0.1452  │
│  geo_02     0.0956    0.0038   0.0402   0.0883   0.1034  │
│  geo_03     0.1126    0.0033   0.0296   0.1062   0.1193  │
│  geo_04     0.1481    0.0047   0.0319   0.1391   0.1577  │
│  geo_05     0.0843    0.0031   0.0373   0.0783   0.0907  │
│  geo_06     0.0747    0.0032   0.0425   0.0687   0.0812  │
│  geo_07     0.0727    0.0040   0.0557   0.0651   0.0810  │
│  geo_08     0.0411    0.0039   0.0938   0.0342   0.0494  │
│  geo_09     0.1148    0.0040   0.0351   0.1071   0.1229  │
│  geo_10     0.1191    0.0047   0.0398   0.1101   0.1287  │
│                                                          │
╰──────────────────────────────────────────────────────────╯

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   342733.0000   10680.8501   0.0312   321798.9184   363667.0816  │
│  geo_02   239113.0000   10113.6242   0.0423   219290.6609   258935.3391  │
│  geo_03   281600.0000    8476.1378   0.0301   264987.0752   298212.9248  │
│  geo_04   370596.0000   12852.9351   0.0347   345404.7101   395787.2899  │
│  geo_05   210960.0000    8080.7698   0.0383   195121.9821   226798.0179  │
│  geo_06   186992.0000    8198.4041   0.0438   170923.4232   203060.5768  │
│  geo_07   181766.0000   10639.1062   0.0585   160913.7350   202618.2650  │
│  geo_08   102927.0000    9977.5362   0.0969    83371.3884   122482.6116  │
│  geo_09   287141.0000   10651.1400   0.0371   266265.1492   308016.8508  │
│  geo_10   297927.0000   12817.4855   0.0430   272805.1901   323048.8099  │
│                                                                          │
╰──────────────────────────────────────────────────────────────────────────╯

# 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    13.6997    0.4269   0.0312   12.8629   14.5365  │
│  geo_02     9.5578    0.4043   0.0423    8.7655   10.3501  │
│  geo_03    11.2561    0.3388   0.0301   10.5920   11.9201  │
│  geo_04    14.8134    0.5138   0.0347   13.8065   15.8204  │
│  geo_05     8.4325    0.3230   0.0383    7.7994    9.0656  │
│  geo_06     7.4744    0.3277   0.0438    6.8321    8.1167  │
│  geo_07     7.2655    0.4253   0.0585    6.4320    8.0990  │
│  geo_08     4.1142    0.3988   0.0969    3.3325    4.8959  │
│  geo_09    11.4776    0.4257   0.0371   10.6431   12.3120  │
│  geo_10    11.9087    0.5123   0.0430   10.9046   12.9129  │
│                                                            │
╰────────────────────────────────────────────────────────────╯

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   136.9970    4.2693   0.0312   128.6293   145.3648  │
│  geo_02    95.5781    4.0426   0.0423    87.6547   103.5015  │
│  geo_03   112.5610    3.3881   0.0301   105.9205   119.2015  │
│  geo_04   148.1344    5.1376   0.0347   138.0650   158.2039  │
│  geo_05    84.3248    3.2300   0.0383    77.9940    90.6556  │
│  geo_06    74.7443    3.2771   0.0438    68.3214    81.1673  │
│  geo_07    72.6554    4.2527   0.0585    64.3203    80.9905  │
│  geo_08    41.1419    3.9882   0.0969    33.3252    48.9587  │
│  geo_09   114.7758    4.2575   0.0371   106.4313   123.1203  │
│  geo_10   119.0872    5.1234   0.0430   109.0455   129.1289  │
│                                                              │
╰──────────────────────────────────────────────────────────────╯

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.1520    0.0113   0.0746   0.1310   0.1757  │
│  Rural   Yes     0.2970    0.0120   0.0403   0.2739   0.3211  │
│  Urban   No      0.0186    0.0037   0.2012   0.0125   0.0275  │
│  Urban   Yes     0.5325    0.0074   0.0140   0.5178   0.5471  │
│                                                               │
╰───────────────────────────────────────────────────────────────╯

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.151987 ┆ 0.296956 │
│ Urban  ┆ 0.018563 ┆ 0.532494 │
└────────┴──────────┴──────────┘

# 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.152 ± 0.011 ┆ 0.297 ± 0.012 │
│ Urban  ┆ 0.019 ± 0.004 ┆ 0.532 ± 0.007 │
└────────┴───────────────┴───────────────┘

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  ┆ 15.198691 ┆ 29.695594 │
│ Urban  ┆ 1.856347  ┆ 53.249369 │
└────────┴───────────┴───────────┘

# 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  ┆ 380234.000686 ┆ 742910.999314 │
│ Urban  ┆ 46441.251273  ┆ 1.3322e6      │
└────────┴───────────────┴───────────────┘

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.1520 [0.1310, 0.1757] ┆ 0.2970 [0.2739, 0.3211] │
│ Urban  ┆ 0.0186 [0.0125, 0.0275] ┆ 0.5325 [0.5178, 0.5471] │
└────────┴─────────────────────────┴─────────────────────────┘

# 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.151987 ┆ 0.011336 ┆ 0.130996 ┆ 0.175662 ┆ Two-Way    ┆ 0.05  │
│ Rural  ┆ Yes         ┆ 0.296956 ┆ 0.011976 ┆ 0.27394  ┆ 0.321051 ┆ Two-Way    ┆ 0.05  │
│ Urban  ┆ No          ┆ 0.018563 ┆ 0.003736 ┆ 0.012477 ┆ 0.027537 ┆ Two-Way    ┆ 0.05  │
│ Urban  ┆ Yes         ┆ 0.532494 ┆ 0.007446 ┆ 0.517817 ┆ 0.547115 ┆ 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.2288    0.0140   0.0613   0.2012   0.2564  │
│                                                 │
│                                                 │
│ Test statistic                                  │
│                                                 │
│     diff         t         df   p_value         │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━         │
│  -0.0212   -1.5141   300.0000    0.1310         │
│                                                 │
╰─────────────────────────────────────────────────╯

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.5141
Degrees of freedom: 300.0
p-value: 0.1310

# 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.0212
95% CI: [-0.0488, 0.0064]

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.021237	0.014026	-0.048839	0.006365	-1.514116	300.0	0.131049

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

shape: (1, 5)

est	se	cv	lci	uci
f64	f64	f64	f64	f64
0.228763	0.014026	0.061313	0.201161	0.256365

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    9116.6293   305.5200   0.0335    8515.3957    9717.8630  │
│  urbrur   Urban   14437.9184   326.4021   0.0226   13795.5907   15080.2461  │
│                                                                             │
│                                                                             │
│ Test statistic                                                              │
│                                                                             │
│       diff         t         df   p_value                                   │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━                                   │
│  5321.2891   11.9023   300.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: 5321.29
95% CI: [4441.48, 6201.10]

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

Rural: 9116.63 (SE: 305.52)
Urban: 14437.92 (SE: 326.40)

# 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	5321.289092	447.080293	4441.478438	6201.099746	11.902312	300.0	5.1425e-27

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.000000
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.1641, p=0.0000
Region geo_02: diff=0.0916, p=0.1097
Region geo_03: diff=-0.0314, p=0.4071
Region geo_04: diff=-0.0785, p=0.0275
Region geo_05: diff=0.0732, p=0.2347
Region geo_06: diff=-0.0313, p=0.5456
Region geo_07: diff=-0.0223, p=0.6639
Region geo_08: diff=0.0028, p=0.9730
Region geo_09: diff=0.0701, p=0.1176
Region geo_10: diff=-0.0229, p=0.5351

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 (lumley2013?), 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.2357   11.4681   300.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: 11.4681
Degrees of freedom: 300.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.2357
95% CI: [0.1952, 0.2761]

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.235669	0.02055	0.195229	0.276109	11.468117	300.0	1.7227e-25

# 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 (10 levels)                  │
│                                           │
│                                           │
│  Estimate   Std Err   CV   Lower   Upper  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│                                           │
│                                           │
│ Test statistic                            │
│                                           │
│  df     Chisq        F   p_value          │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━          │
│   9   24.6464   2.7385    0.0044          │
│                                           │
╰───────────────────────────────────────────╯

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: 2.7385
Numerator df: 9
Denominator df: 292
p-value: 0.0044

# 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"	9	292.0	24.64644	2.738493	0.004394

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    │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━    │
│  0.7858   10.9563   300.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"	0.785846	0.071725	0.644698	0.926995	10.956315	300.0	1.0138e-23

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.3653   11.9190   300.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.365278	0.030647	0.304968	0.425587	11.919049	300.0	4.4873e-27

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.2664   11.8310   300.0000   <0.0001    │
│                                           │
╰───────────────────────────────────────────╯

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.26645	0.022521	0.22213	0.31077	11.831011	300.0	9.1831e-27

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.000000
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.

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`