4 Understanding Distributions
This chapter uses {tidyverse} for data manipulation, {srvyr} for survey estimation, {broom} for tidying model-style output, and {scales} for formatting large numbers and percentages in the text.
Every estimate in this chapter comes from the ICAN-ICAR 2025 design object built in Chapter ??. The chunk below reads the data, adds a handful of labelled variables that make the output easier to read, and rebuilds the design object. Documentation for the survey is in the DataFirst metadata record.
options(
survey.lonely.psu = "adjust",
survey.adjust.domain.lonely = TRUE
)
icanicar_2025 <- read_csv("data/ican-icar-2025-v1.csv")
icanicar_des <- icanicar_2025 |>
mutate(
sex = factor(ch03, levels = c(1, 2), labels = c("Female", "Male")),
residence = factor(str_trim(Location), levels = c("Rural", "Urban")),
schooling = factor(
EnrolmentStatus,
levels = c("Currently Enrolled", "Out of School", "Never Enrolled")
),
mpl_math = factor(
MPLMath,
levels = c(0, 1),
labels = c("Below MPL", "At or above MPL")
),
mpl_reading = factor(
MPLReading,
levels = c(0, 1),
labels = c("Below MPL", "At or above MPL")
),
psu = interaction(CountryName, VillageID, drop = TRUE),
stratum = interaction(CountryName, TierOneUnit, drop = TRUE)
) |>
as_survey_design(
ids = c(psu, HHID),
strata = stratum,
weights = HHWeightProvided,
nest = TRUE
)4.1 From a sample of children to a population of learners
The ICAN-ICAR 2025 file holds one row per assessed child: 96,452 children in 11 countries. Those children were not picked one at a time from a national list. They were reached through sampled enumeration areas and sampled households, so some children stand in for far more of their peers than others. This is why a plain count() and a survey-weighted count answer different questions. The first describes who was interviewed; the second describes the population of children those interviews represent.
Descriptive analysis is where that distinction becomes visible for the first time, and it is where most survey analysis begins and ends. Before fitting any model we usually want to know how many children are out of school, what share of them can read a short story or handle two-digit subtraction, how learning outcomes are spread within a country, and whether the gap between rural and urban children is large enough to take seriously. Each of those questions is a statement about a distribution, and each needs a point estimate together with an honest measure of uncertainty.
4.1.1 Matching statistics to the type of variable
Which summary is appropriate depends on the variable. Any variable in the data is made up of two or more categories or attributes. A child’s sex has the categories female and male; the language of assessment has categories such as Swahili, Bangla, and Wolof.
The level of measurement of a variable describes how its categories relate to one another. Three questions settle the level of measurement: are the categories distinct; can they be placed in a meaningful rank order; and are the intervals between them numerically meaningful? Working through those questions gives four types of variable, all four of which appear in ICAN-ICAR 2025:
- Categorical (nominal) variables have categories that cannot be ranked:
CountryName,AssessmentLanguage,Location, and the derivedsexvariable. Summarise these with counts and proportions. - Ordinal variables can be ranked but the gaps between categories have no fixed size. The functional difficulty items such as
ch04a, which runs from no difficulty to cannot do at all, are ordinal. Counts, proportions, and medians are appropriate; means are usually not. - Discrete variables are counted in whole units: household size (
hh06a), the child’s age in completed years (ch02), and the number of correct items on an assessment. - Continuous variables can take any value on an interval: the mathematics and reading ability scores (
MathIRTScore,ReadingIRTScore) and the time each assessment took (IcanAssessTime,IcarAssessTime). Means, quantiles, and measures of spread all apply.
Proficiency indicators such as MPLMath, MPLReading, and MPLBoth sit in a convenient middle ground. They are categorical, coded 0 and 1, so the mean of the indicator is the proportion of children reaching the minimum proficiency level. That equivalence is used repeatedly below.
4.1.2 What this chapter covers
All the estimators used here come from {srvyr}, which wraps the {survey} package in {dplyr} syntax so that survey estimation looks like ordinary data manipulation. They fall into four families:
- Distribution. How often something occurs, and how much of it there is in total:
survey_count(),survey_tally(), andsurvey_total(). - Central tendency. Where the middle of a distribution sits:
survey_mean(),survey_prop(),survey_quantile(), andsurvey_median(). - Relationship. How two variables move together:
survey_ratio()andsurvey_corr(). - Dispersion. How widely values are spread:
survey_var()andsurvey_sd().
The workflow from Chapter ?? applies to all of them:
- Build a
tbl_svyobject withas_survey_design(). - Restrict to a subpopulation with
filter(), if needed. - Define domains of analysis with
group_by(), if needed. - Call the estimator inside
summarize().
Steps 2 and 3 must come after step 1. Filtering or grouping the raw data frame before the design is built hides part of the sample from the variance calculation, and the standard errors that come out the other side will be wrong. Working from icanicar_des rather than icanicar_2025 is the safeguard.
4.1.3 Missing values are not spread evenly
One more thing to check before estimating anything. The assessment outcomes are not observed for every child:
icanicar_2025 |>
summarize(across(
c(MathIRTScore, ReadingIRTScore, MPLMath, MPLBoth, IcanAssessTime),
~ sum(is.na(.x))
))## # A tibble: 1 × 5
## MathIRTScore ReadingIRTScore MPLMath MPLBoth IcanAssessTime
## <int> <int> <int> <int> <int>
## 1 7314 7314 7314 7314 18215In total, 7,314 children have no proficiency classification, and more than twice that many have no recorded ICAN timing. Survey functions in {srvyr} refuse to guess what to do about this: na.rm defaults to FALSE, so an unguarded call returns NA. Setting na.rm = TRUE drops the missing cases and quietly changes the denominator to “children with a valid assessment”. That is usually what we want, but it is a substantive choice rather than a technical one, so this chapter states it every time. Chapter ?? points to the survey documentation on why assessments are missing.
4.2 Counting children and cross-tabulating
survey_count() estimates how many members of the population fall into a category, or into each combination of several categories. The cross-tabulations it produces are the natural first summary of any categorical variable. survey_tally() does the same job for data that have already been grouped.
4.2.1 Calling the functions
survey_count() mirrors dplyr::count(), with extra arguments controlling the uncertainty estimate:
survey_count(
x,
...,
wt = NULL,
sort = FALSE,
name = "n",
.drop = dplyr::group_by_drop_default(x),
vartype = c("se", "ci", "var", "cv")
)Reading the arguments in turn: x is the design object; ... takes the variables to tabulate, exactly as count() would; wt supplies a second weight to apply on top of the survey weights, which we rarely need because the design object already carries them; sort orders the output by size; name sets the column name for the estimate, n by default; .drop controls whether empty combinations disappear from the output; and vartype chooses which uncertainty measures to return. The full argument documentation is in ?survey_count.
survey_tally() takes the same arguments except ... and .drop, because the grouping has already been done by group_by():
The vartype argument appears in nearly every function in this chapter, so it is worth settling now. It accepts any combination of four values, and each one adds columns named after the estimate:
"se"(the default) returns the standard error, suffixed_se. It measures how much the estimate would move around if the survey were repeated, and it is the usual currency for judging precision."ci"returns the bounds of a confidence interval, suffixed_lowand_upp. The default level is 95%, changed through thelevelargument, solevel = 0.9gives a 90% interval."var"returns the variance of the estimate, suffixed_var. It is the square of the standard error."cv"returns the coefficient of variation, suffixed_cv, which is the standard error divided by the estimate. Because it is relative, it is the easiest way to spot estimates that are too imprecise to publish.
Unless a function documents otherwise, confidence intervals are symmetric and based on the t-distribution:
\[\text{estimate} \pm t^*_{df} \times SE\]
The degrees of freedom come from the design rather than the number of children. For a stratified cluster design they are the number of first-stage clusters minus the number of strata, which {survey} reports directly:
## [1] 2654ICAN-ICAR 2025 has 2,927 sampled enumeration areas spread over 273 country-by-region strata, leaving 2,654 degrees of freedom. That is comfortably large, so t-based intervals here are almost identical to normal-based ones. In smaller designs, or in analyses restricted to one country, the difference matters and the df argument lets us override the default.
4.2.2 How many children does the sample represent?
Called with no variables, survey_count() estimates the size of the whole population of interest:
## # A tibble: 1 × 2
## children children_se
## <dbl> <dbl>
## 1 299559314. 4148332.The 96,452 children in the file represent an estimated 299.6 million children across the eleven participating countries, with a standard error of 4.1 million. Two things are worth noticing. First, the estimate is enormous compared with the sample, which is exactly what weights are for. Second, this is a count of children in the age range the survey targets, because the household weight attaches to every child assessed in that household. It is not a count of households, and it is not a count of all children of all ages.
Estimates like this deserve a plausibility check against an external source before they go anywhere near a report. If a weighted count is wildly out of line with the national population of school-age children, the likely explanation is a weighting or subsetting mistake rather than a surprising finding.
4.2.3 Cross-tabulating by country and location
Passing variables to survey_count() breaks the estimate down by every combination of them. Here the estimated number of children is split by country and by rural or urban residence, with confidence intervals instead of standard errors and the count column renamed:
children_by_area <- icanicar_des |>
survey_count(CountryName, residence, name = "children", vartype = "ci")
children_by_area## # A tibble: 22 × 5
## CountryName residence children children_low children_upp
## <chr> <fct> <dbl> <dbl> <dbl>
## 1 Bangladesh Rural 33014359. 29419564. 36609154.
## 2 Bangladesh Urban 18025058. 13614137. 22435979.
## 3 Kenya Rural 9848127. 8015413. 11680840.
## 4 Kenya Urban 10714954. 9014203. 12415705.
## 5 Mali Rural 2797936. 2516464. 3079408.
## 6 Mali Urban 1382335. 1204288. 1560382.
## 7 Mexico Rural 23754324. 19836439. 27672210.
## 8 Mexico Urban 33583260. 29475964. 37690555.
## 9 Mozambique Rural 6697553. 5907171. 7487934.
## 10 Mozambique Urban 3877823. 3079991. 4675655.
## # ℹ 12 more rowsThe table has one row per country and residence combination. In Tanzania, for example, an estimated 17.4 million children live in rural areas and 12.6 million in urban areas. Running this code interactively also produces a stack of warnings about strata containing a single primary sampling unit; the next section explains what they mean and why the chapter suppresses them.
survey_tally() reaches the same answer from a grouped design:
## # A tibble: 2 × 3
## residence children children_se
## <fct> <dbl> <dbl>
## 1 Rural 166999837. 4655677.
## 2 Urban 132559478. 4503887.The two functions are not interchangeable, though. Handing a variable to survey_tally() fails, because it tries to sum that variable rather than group by it:
## Error in `dplyr::summarise()`:
## ℹ In argument: `n = survey_total(CountryName, vartype = vartype, na.rm = TRUE)`.
## Caused by error:
## ! Character vectors not allowed in survey functions, should be used as a grouping variable.The error message is a useful reminder of the division of labour: variables named inside survey_count() are grouping variables, whereas anything passed to a survey_*() estimator is treated as a quantity to be summed or averaged.
4.3 When a stratum contains only one village
The warnings above are specific to this dataset and worth understanding rather than ignoring. Variance estimation for a stratified design works by comparing clusters within each stratum. A stratum that contains a single sampled enumeration area offers nothing to compare, so its contribution to the variance is undefined. A number of strata in ICAN-ICAR 2025 are in exactly that position:
lonely_strata <- icanicar_2025 |>
distinct(CountryName, TierOneUnit, VillageID) |>
count(CountryName, TierOneUnit, name = "villages") |>
filter(villages == 1)
lonely_strata## # A tibble: 63 × 3
## CountryName TierOneUnit villages
## <chr> <chr> <int>
## 1 Mozambique c120 1
## 2 Nepal c000 1
## 3 Tanzania c000 1
## 4 Uganda c106 1
## 5 Uganda c107 1
## 6 Uganda c110 1
## 7 Uganda c113 1
## 8 Uganda c114 1
## 9 Uganda c115 1
## 10 Uganda c117 1
## # ℹ 53 more rowsThere are 63 such strata, spread across 4 countries. Left alone, {survey} raises an error when an estimate touches one of them. The setup chunk therefore sets
which centres those strata at the overall population mean instead of dropping them, giving a conservative (slightly larger) variance rather than no variance at all. The companion option survey.adjust.domain.lonely = TRUE applies the same treatment when a group_by() or filter() leaves a stratum with one cluster in a particular domain, which happens often once estimates are broken down by country and residence together.
Two habits follow from this. Set both options once, near the top of the analysis, so that every estimate is computed the same way. And treat a long list of lonely-stratum warnings as a signal to check whether the domain being estimated is too fine for the design to support, because a conservative variance is still a poor substitute for an adequate sample.
4.4 Adding things up
survey_total() is the survey counterpart of sum(). Applied to a continuous or count variable, it estimates the population total of that quantity. Applied to a 0/1 indicator, it estimates the number of population members for whom the indicator is 1. Unlike survey_count() and survey_tally(), it must be called inside summarize().
4.4.1 Calling the function
survey_total(
x,
na.rm = FALSE,
vartype = c("se", "ci", "var", "cv"),
level = 0.95,
deff = FALSE,
df = NULL
)x is the variable or expression to add up, and it can be left empty to get a population count. na.rm decides whether missing values are dropped. vartype and level behave as described above. deff adds the design effect to the output, discussed in the section on design effects. df overrides the degrees of freedom used for confidence intervals.
4.4.2 How many children reach both minimum proficiency levels?
MPLBoth is coded 1 for children who reach the minimum proficiency level in both reading and mathematics. Totalling it estimates how many children in the population have done so:
mpl_both_total <- icanicar_des |>
summarize(children = survey_total(
MPLBoth,
na.rm = TRUE,
vartype = c("se", "ci")
))
mpl_both_total## # A tibble: 1 × 4
## children children_se children_low children_upp
## <dbl> <dbl> <dbl> <dbl>
## 1 95903033. 1861271. 92253344. 99552723.An estimated 95.9 million children reach both thresholds, with a 95% confidence interval running from 92.3 million to 99.6 million. Set against the 299.6 million children the survey represents, this is the headline result of the whole exercise: only about one child in three clears both bars.
4.4.3 Totals by country
Grouping first gives the same total for each country. Confidence intervals are more informative than standard errors when the point estimates differ by orders of magnitude, and sorting makes the table easier to read:
mpl_by_country <- icanicar_des |>
group_by(CountryName) |>
summarize(children = survey_total(MPLBoth, na.rm = TRUE, vartype = "ci")) |>
arrange(desc(children))
mpl_by_country## # A tibble: 11 × 4
## CountryName children children_low children_upp
## <chr> <dbl> <dbl> <dbl>
## 1 Mexico 29995006. 28134156. 31855856.
## 2 Pakistan 22611627. 20041503. 25181751.
## 3 Bangladesh 15350639. 13932137. 16769142.
## 4 Tanzania 8397215. 7896388. 8898041.
## 5 Kenya 8243830. 7490214. 8997446.
## 6 Nepal 6408482. 5970729. 6846235.
## 7 Uganda 2111516. 1860613. 2362420.
## 8 Mozambique 1049926. 906103. 1193749.
## 9 Senegal 923783. 831009. 1016557.
## 10 Nicaragua 603788. 505058. 702518.
## 11 Mali 207221. 173033. 241409.These totals are driven as much by population size as by learning outcomes: a large country with mediocre proficiency contributes more children above the threshold than a small country doing well. Totals answer “how many”, which is what budgeting and planning need. To compare how well countries are doing, we need proportions, which is the next section.
4.4.4 A total that answers no question
Totals are easy to compute and easy to misuse. hh06a records household size, and every child assessed in a household carries the same value, so households with several assessed children appear several times:
## # A tibble: 1 × 2
## members members_se
## <dbl> <dbl>
## 1 1810518465. 27342661.The result runs into billions, which is far more people than live in the eleven countries combined. Nothing is wrong with the arithmetic; the problem is the unit of analysis. The design object is a sample of children, so summing a household attribute over children counts large households once for each child they contributed. Household-level totals need one row per household, which the chapter on central tendency and dispersion constructs before summarising household size.
The general rule: before totalling a variable, ask which unit it describes, and check that this unit matches the rows of the design object.
4.5 Proportions and means
Proportions and means are the estimates most readers of a report will actually look at. survey_prop() is for categorical variables and survey_mean() for continuous ones, but the two overlap: applied to an indicator coded 0/1, or to a logical variable, survey_mean() returns a proportion.
4.5.1 Calling the functions
survey_mean(
x,
na.rm = FALSE,
vartype = c("se", "ci", "var", "cv"),
level = 0.95,
proportion = FALSE,
prop_method = c("logit", "likelihood", "asin", "beta", "mean"),
deff = FALSE,
df = NULL
)
survey_prop(
na.rm = FALSE,
vartype = c("se", "ci", "var", "cv"),
level = 0.95,
proportion = TRUE,
prop_method =
c("logit", "likelihood", "asin", "beta", "mean", "xlogit"),
deff = FALSE,
df = NULL
)The shared arguments (na.rm, vartype, level, deff, df) work as before. Two differences matter in practice.
The first is where the variable goes. survey_mean() takes it as the argument x. survey_prop() has no such argument: the categories come from a preceding group_by(), and the function reports the share of the population in each. Adding group_by() before survey_mean() does something different — it splits the mean by domain.
The second is the proportion argument, which controls how the confidence interval is built. Wald-type intervals of the form estimate plus or minus a multiple of the standard error can stray outside the range 0 to 1 and can have poor coverage when a proportion sits close to 0 or 1. Setting proportion = TRUE, which is the default in survey_prop(), switches to a method chosen by prop_method:
"logit"(the default) works on the log-odds scale, where an interval cannot run past 0 or 1, and converts the bounds back to proportions afterwards."likelihood"inverts a Rao-Scott scaled chi-squared test on the binomial log-likelihood."asin"builds the interval on the arcsine square-root scale, which holds the binomial variance roughly constant, then transforms back."beta"treats the estimate as binomial with an effective sample size implied by its variance, and reads the bounds from the beta distribution."mean"is the plain Wald interval, and is whatsurvey_mean()gives whenproportion = FALSE."xlogit"transforms the proportion with a logit, builds the interval, and transforms back; this is the default in SUDAAN and SPSS, and is the option to choose when results must line up with output from those packages.
4.5.2 The distribution of a single categorical variable
schooling records whether a child is currently enrolled, has dropped out, or has never been to school. Grouping by it and calling survey_prop() estimates the share of children in each situation:
schooling_p <- icanicar_des |>
group_by(schooling) |>
summarize(p = survey_prop(vartype = "ci"))
schooling_p## # A tibble: 3 × 4
## schooling p p_low p_upp
## <fct> <dbl> <dbl> <dbl>
## 1 Currently Enrolled 0.924 0.918 0.929
## 2 Out of School 0.0246 0.0221 0.0274
## 3 Never Enrolled 0.0519 0.0478 0.0563An estimated 92.4% of children are enrolled, 5.2% have never enrolled, and 2.5% have left school. The proportions add to 1 because the three categories are exhaustive. The interval for the smallest category, 2.2% to 2.7%, is narrow in absolute terms but wide relative to the estimate itself, which is typical of rare outcomes.
Because a proportion is the mean of a set of indicator variables, survey_mean() produces exactly the same numbers when the design is grouped and no variable is named:
## # A tibble: 3 × 4
## schooling p p_low p_upp
## <fct> <dbl> <dbl> <dbl>
## 1 Currently Enrolled 0.924 0.918 0.929
## 2 Out of School 0.0246 0.0220 0.0273
## 3 Never Enrolled 0.0519 0.0476 0.05614.5.3 Conditional and joint proportions
Grouping by two variables changes the meaning of the estimate, and this is the single most common source of confusion in reporting cross-tabulations. With group_by(residence, mpl_math), the proportions are conditional: they sum to 1 within each residence category, so they read as proficiency rates for rural and for urban children.
mpl_by_residence <- icanicar_des |>
filter(!is.na(mpl_math)) |>
group_by(residence, mpl_math) |>
summarize(p = survey_prop(vartype = "ci"))
mpl_by_residence## # A tibble: 4 × 5
## # Groups: residence [2]
## residence mpl_math p p_low p_upp
## <fct> <fct> <dbl> <dbl> <dbl>
## 1 Rural Below MPL 0.536 0.522 0.550
## 2 Rural At or above MPL 0.464 0.450 0.478
## 3 Urban Below MPL 0.439 0.424 0.455
## 4 Urban At or above MPL 0.561 0.545 0.576Read this way, 46.4% of rural children reach the mathematics threshold compared with 56.1% of urban children, a gap of 9.6 percentage points. The confidence intervals do not overlap, so the gap is larger than sampling error alone would explain. The chapter on bivariate analysis shows how to test such a difference formally.
Wrapping the grouping variables in interact() instead gives joint proportions, which sum to 1 across the whole table:
icanicar_des |>
filter(!is.na(mpl_math)) |>
group_by(interact(residence, mpl_math)) |>
summarize(p = survey_prop())## # A tibble: 4 × 4
## residence mpl_math p p_se
## <fct> <fct> <dbl> <dbl>
## 1 Rural Below MPL 0.296 0.00777
## 2 Rural At or above MPL 0.257 0.00771
## 3 Urban Below MPL 0.197 0.00694
## 4 Urban At or above MPL 0.251 0.00819Now each cell is a share of all children rather than of children in that residence category. The urban proficient cell, for instance, is the proportion of children in the population who are both urban and above the threshold, a quantity that depends on how urbanised these countries are as well as on how well urban children perform. Conditional proportions compare groups; joint proportions describe the composition of the population. Choose deliberately, and say in the text which one is being reported.
4.5.4 When the choice of interval method matters
For the estimates above, all six prop_method options return practically the same interval, because the effective sample is large and the proportions are not near a boundary. The methods separate when the estimate is small and the subpopulation is thin. Take the mathematics proficiency rate among five- and six-year-olds in Mali:
young_mali <- icanicar_des |>
filter(CountryName == "Mali", ch02 %in% c(5, 6), !is.na(MPLMath))
prop_methods <- c("logit", "likelihood", "asin", "beta", "mean", "xlogit")
interval_comparison <- prop_methods |>
map(\(method) {
young_mali |>
summarize(p = survey_mean(
MPLMath,
proportion = TRUE,
prop_method = method,
vartype = "ci"
)) |>
mutate(prop_method = method)
}) |>
list_rbind() |>
select(prop_method, p, p_low, p_upp)
interval_comparison## # A tibble: 6 × 4
## prop_method p p_low p_upp
## <chr> <dbl> <dbl> <dbl>
## 1 logit 0.00322 0.00129 0.00800
## 2 likelihood 0.00322 0.00109 0.00710
## 3 asin 0.00322 0.000948 0.00682
## 4 beta 0.00322 0.000981 0.00775
## 5 mean 0.00322 0.000278 0.00615
## 6 xlogit 0.00322 0.00129 0.00800The point estimate is identical in every row, as it must be, but the lower bound varies by a factor of four across methods. The Wald interval ("mean") reaches furthest towards zero and would cross it for a slightly rarer outcome, which is precisely the failure the other methods are designed to avoid. The practical advice is to leave the default "logit" in place, and to report which method was used whenever an estimate is small enough for the choice to change the conclusion.
4.5.5 Means of continuous variables
MathIRTScore is a latent ability score, standardised so that values near zero represent typical performance across the pooled sample. Its population mean, with both a standard error and a confidence interval:
math_mean <- icanicar_des |>
summarize(score = survey_mean(
MathIRTScore,
na.rm = TRUE,
vartype = c("se", "ci")
))
math_mean## # A tibble: 1 × 4
## score score_se score_low score_upp
## <dbl> <dbl> <dbl> <dbl>
## 1 -0.263 0.0128 -0.288 -0.238Grouping splits the mean by domain. Sorting the countries makes the spread of average performance immediately visible:
math_by_country <- icanicar_des |>
group_by(CountryName) |>
summarize(score = survey_mean(MathIRTScore, na.rm = TRUE, vartype = "ci")) |>
arrange(score)
math_by_country## # A tibble: 11 × 4
## CountryName score score_low score_upp
## <chr> <dbl> <dbl> <dbl>
## 1 Mali -1.15 -1.23 -1.07
## 2 Mozambique -0.813 -0.864 -0.762
## 3 Uganda -0.602 -0.639 -0.565
## 4 Senegal -0.437 -0.485 -0.389
## 5 Tanzania -0.385 -0.431 -0.339
## 6 Bangladesh -0.330 -0.366 -0.293
## 7 Nepal -0.186 -0.229 -0.144
## 8 Nicaragua -0.147 -0.218 -0.0756
## 9 Pakistan -0.104 -0.205 -0.00176
## 10 Mexico -0.0501 -0.0857 -0.0144
## 11 Kenya 0.0281 -0.0257 0.0818Average mathematics ability ranges from -1.15 in Mali to 0.03 in Kenya, a spread of more than a standard deviation of the pooled score distribution. Because the score has no natural units, differences like these are best described in relative terms or converted to proficiency rates, which are what most readers will want anyway.
4.6 Quantiles and medians
A mean compresses a distribution into a single number and can be pulled away from the typical case by a few extreme values. Quantiles describe the shape of the distribution instead. survey_quantile() estimates any set of quantiles; survey_median() is a convenience wrapper for the 50th percentile.
4.6.1 Calling the functions
survey_quantile(
x,
quantiles,
na.rm = FALSE,
vartype = c("se", "ci", "var", "cv"),
level = 0.95,
interval_type =
c("mean", "beta", "xlogit", "asin", "score", "quantile"),
qrule = c("math", "school", "shahvaish", "hf1", "hf2", "hf3",
"hf4", "hf5", "hf6", "hf7", "hf8", "hf9"),
df = NULL
)
survey_median(
x,
na.rm = FALSE,
vartype = c("se", "ci", "var", "cv"),
level = 0.95,
interval_type =
c("mean", "beta", "xlogit", "asin", "score", "quantile"),
qrule = c("math", "school", "shahvaish", "hf1", "hf2", "hf3",
"hf4", "hf5", "hf6", "hf7", "hf8", "hf9"),
df = NULL
)The only structural difference is the quantiles argument, a vector of values between 0 and 1 that survey_median() fixes at 0.5. Two further arguments deserve comment. interval_type selects the confidence interval method; besides the options already met for proportions, quantiles offer "score", which inverts a score test and is available for design-based objects only, and "quantile", which uses the replicate distribution and is valid for bootstrap and balanced repeated replication weights but not for jackknife weights. The score method behaves poorly when many observations tie on the same value, which is common for ages, grades, and item counts, so an interval-based method is safer for those variables. qrule decides how a quantile is defined when it falls between two observed values; vignette("qrule", package = "survey") works through the alternatives.
4.6.2 Quartiles of the mathematics score
math_quartiles <- icanicar_des |>
summarize(score = survey_quantile(
MathIRTScore,
quantiles = c(0.25, 0.5, 0.75),
na.rm = TRUE
))
math_quartiles## # A tibble: 1 × 6
## score_q25 score_q50 score_q75 score_q25_se score_q50_se score_q75_se
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 -0.908 -0.230 0.330 0.0165 0.00863 0.0133The quartiles are -0.91, -0.23, and 0.33. The middle half of children therefore sit within a band about 1.24 score points wide. The median, -0.23, is slightly above the mean of -0.26, which points to a mild left skew: a tail of very low-scoring children pulls the average down below the typical child.
Medians can be grouped like any other estimate:
icanicar_des |>
group_by(residence) |>
summarize(score = survey_median(MathIRTScore, na.rm = TRUE, vartype = "ci"))## # A tibble: 2 × 4
## residence score score_low score_upp
## <fct> <dbl> <dbl> <dbl>
## 1 Rural -0.317 -0.343 -0.289
## 2 Urban -0.134 -0.166 -0.1054.6.3 Where quantiles earn their keep
Assessment duration shows why a median is sometimes the only defensible summary. IcanAssessTime records how long the numeracy assessment took in minutes, and a handful of records run to several hours, almost certainly because a tablet was left open rather than because a child worked for that long.
time_summary <- icanicar_des |>
summarize(
mean_minutes = survey_mean(IcanAssessTime, na.rm = TRUE, vartype = NULL),
median_minutes = survey_median(IcanAssessTime, na.rm = TRUE, vartype = NULL),
p90_minutes = survey_quantile(
IcanAssessTime,
quantiles = 0.9,
na.rm = TRUE,
vartype = NULL
),
longest_minutes = unweighted(max(IcanAssessTime, na.rm = TRUE))
)
time_summary## # A tibble: 1 × 4
## mean_minutes median_minutes p90_minutes_q90 longest_minutes
## <dbl> <dbl> <dbl> <dbl>
## 1 8.49 6.8 14.2 713.The median assessment takes 6.8 minutes and nine children in ten finish within 14.2 minutes, yet the mean is 8.5 minutes and the longest record is 713 minutes. Reporting the mean here would misdescribe almost every assessment in the survey. Note also unweighted(), used to pull a sample maximum out of the design object; it is covered properly below.
Quantiles of 0 and 1 are accepted but are a trap:
icanicar_des |>
summarize(score = survey_quantile(
MathIRTScore,
quantiles = c(0, 1),
na.rm = TRUE
))## # A tibble: 1 × 4
## score_q00 score_q100 score_q00_se score_q100_se
## <dbl> <dbl> <dbl> <dbl>
## 1 -3.07 2.14 0.101 0These are the smallest and largest values in the sample. The maximum has a standard error of exactly 0, which is not a claim of perfect precision but a sign that no valid uncertainty can be attached to an extreme order statistic. Extremes do not generalise from a sample to a population, so avoid presenting them as population estimates.
4.7 Ratios
A ratio estimates the total of one variable divided by the total of another:
\[\frac{\sum x_i}{\sum y_i}\]
This is not the same as averaging the individual ratios,
\[\frac{1}{N}\sum \frac{x_i}{y_i},\]
which is what survey_mean() returns if we construct the quotient first. The distinction is easy to state and easy to forget: the first weights each case by the size of its denominator, so big cases count for more; the second treats every case equally.
4.7.1 Calling the function
survey_ratio(
numerator,
denominator,
na.rm = FALSE,
vartype = c("se", "ci", "var", "cv"),
level = 0.95,
deff = FALSE,
df = NULL
)The first two arguments are the numerator and denominator variables; the rest behave as in the other estimators.
4.7.2 How much longer is the numeracy assessment than the reading one?
Every child who completed both assessments has a timing for each. The two ways of expressing “how much longer ICAN takes than ICAR” give different answers:
time_ratio <- icanicar_des |>
summarize(
ican_total = survey_total(IcanAssessTime, na.rm = TRUE, vartype = NULL),
icar_total = survey_total(IcarAssessTime, na.rm = TRUE, vartype = NULL),
ratio_of_totals = survey_ratio(IcanAssessTime, IcarAssessTime, na.rm = TRUE),
mean_of_ratios = survey_mean(
IcanAssessTime / IcarAssessTime,
na.rm = TRUE,
vartype = NULL
)
)
time_ratio## # A tibble: 1 × 5
## ican_total icar_total ratio_of_totals ratio_of_totals_se mean_of_ratios
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2061674526. 1563295302. 1.30 0.0231 1.78The ratio of totals is 1.30: across all the assessment time the survey represents, 1.30 minutes went into numeracy for every minute spent on reading. It is exactly the quotient of the two totals shown alongside it, which is worth confirming by hand once, although only survey_ratio() gives the correct standard error for it. The mean of the individual ratios is 1.78, noticeably larger, because children who raced through the reading assessment produce very large individual quotients and each of them counts as much as a child who took an hour.
Neither number is wrong. The ratio of totals answers “how is total assessment time divided between the two instruments”, which is the question a fieldwork budget asks. The mean of ratios answers “for a typical child, how do the two durations compare”. Ratios of totals also vary by context:
icanicar_des |>
group_by(CountryName) |>
summarize(ratio = survey_ratio(IcanAssessTime, IcarAssessTime, na.rm = TRUE)) |>
arrange(ratio)## # A tibble: 11 × 3
## CountryName ratio ratio_se
## <chr> <dbl> <dbl>
## 1 Bangladesh 1.11 0.0247
## 2 Nepal 1.16 0.0328
## 3 Pakistan 1.20 0.0952
## 4 Mozambique 1.21 0.0893
## 5 Tanzania 1.23 0.0368
## 6 Kenya 1.28 0.0470
## 7 Mexico 1.41 0.0332
## 8 Nicaragua 1.43 0.0765
## 9 Uganda 1.67 0.0425
## 10 Senegal 1.83 0.138
## 11 Mali 2.10 0.2354.8 Correlation
survey_corr() estimates Pearson’s correlation between two continuous variables, applying the survey weights. For a simple random sample the sample correlation is
\[\frac{\sum (x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\sum (x_i - \bar{x})^2}\sqrt{\sum (y_i - \bar{y})^2}}\]
and the survey version replaces the unweighted sums with weighted ones. The result stays between -1 and 1 and measures only the strength of a linear relationship.
4.8.1 Calling the function
x and y are the two variables or expressions; the remaining arguments are the familiar ones.
4.8.2 Do reading and numeracy travel together?
score_corr <- icanicar_des |>
summarize(correlation = survey_corr(
MathIRTScore,
ReadingIRTScore,
na.rm = TRUE,
vartype = "ci"
))
score_corr## # A tibble: 1 × 3
## correlation correlation_low correlation_upp
## <dbl> <dbl> <dbl>
## 1 0.794 0.786 0.803The correlation is 0.79, a strong positive relationship: children who read well tend to do well in mathematics, which is what we would expect if both depend on shared foundations such as language of instruction, years of effective schooling, and home support. It is not close enough to 1 to treat the two scores as interchangeable, and correlation says nothing about which skill supports the other.
Splitting the estimate by domain tests whether the association itself varies:
icanicar_des |>
group_by(residence) |>
summarize(correlation = survey_corr(
MathIRTScore,
ReadingIRTScore,
na.rm = TRUE
))## # A tibble: 2 × 3
## residence correlation correlation_se
## <fct> <dbl> <dbl>
## 1 Rural 0.792 0.00666
## 2 Urban 0.792 0.00503The association is essentially the same in rural and urban areas. Levels of achievement differ sharply between the two, as the earlier proportions showed, but the way reading and numeracy move together does not.
4.9 Variance and standard deviation
Every estimator in this chapter reports the uncertainty of its own estimate, and that is a different thing from the spread of the underlying variable. When the spread itself is the object of interest, survey_var() and survey_sd() estimate the population variance and standard deviation. These come up less often than the other functions, but they are useful for describing inequality in outcomes and for planning the sample size of a future survey.
4.9.1 Calling the functions
survey_var(
x,
na.rm = FALSE,
vartype = c("se", "ci", "var"),
level = 0.95,
df = NULL
)
survey_sd(
x,
na.rm = FALSE
)Since the standard deviation is just the variance put back on the scale of the original variable, the two functions share their arguments, except that survey_sd() accepts no vartype: it is the one estimate in the chapter that comes without a standard error.
4.9.2 How unequal are learning outcomes?
icanicar_des |>
summarize(
variance = survey_var(MathIRTScore, na.rm = TRUE),
std_dev = survey_sd(MathIRTScore, na.rm = TRUE)
)## # A tibble: 1 × 3
## variance variance_se std_dev
## <dbl> <dbl> <dbl>
## 1 1.12 0.0266 1.06The output gives the estimated population variance, its standard error, and the estimated standard deviation. Depending on the versions of R and {survey} installed, survey_var() may also emit a warning about a deprecated recycling operation in its internal calculation; the estimates are unaffected.
Comparing spread across countries is often more revealing than comparing averages, because a country can raise its mean while leaving the weakest children behind:
spread_by_country <- icanicar_des |>
group_by(CountryName) |>
summarize(
variance = survey_var(MathIRTScore, na.rm = TRUE, vartype = NULL),
std_dev = survey_sd(MathIRTScore, na.rm = TRUE)
) |>
arrange(desc(std_dev))
spread_by_country## # A tibble: 11 × 3
## CountryName variance std_dev
## <chr> <dbl> <dbl>
## 1 Pakistan 2.50 1.58
## 2 Tanzania 1.15 1.07
## 3 Senegal 1.04 1.02
## 4 Kenya 0.864 0.930
## 5 Mali 0.806 0.898
## 6 Mozambique 0.668 0.818
## 7 Uganda 0.654 0.809
## 8 Mexico 0.591 0.769
## 9 Nicaragua 0.555 0.745
## 10 Bangladesh 0.538 0.734
## 11 Nepal 0.497 0.705Scores are most dispersed in Pakistan (standard deviation 1.58) and most tightly clustered in Nepal (0.71). A wide distribution and a modest mean together suggest a system serving some children well and others hardly at all, which calls for a different response than a system where nearly everyone performs poorly.
4.10 Rounding out the toolkit
The remaining tools are less about new statistics than about producing analysis that is efficient to write and honest to read.
4.10.1 Comparing weighted and unweighted results
unweighted() evaluates an ordinary R expression on the rows of the design object, ignoring the weights. It describes the children who were assessed rather than the population they represent, and putting the two side by side shows how much the weights are doing:
icanicar_des |>
summarize(
weighted_mean = survey_mean(MathIRTScore, na.rm = TRUE, vartype = NULL),
sample_mean = unweighted(mean(MathIRTScore, na.rm = TRUE)),
children_assessed = unweighted(sum(!is.na(MathIRTScore)))
)## # A tibble: 1 × 3
## weighted_mean sample_mean children_assessed
## <dbl> <dbl> <int>
## 1 -0.263 -0.416 89138The two means differ substantially. The sample mean is the lower of the two because countries and areas with weaker results are over-represented in the sample relative to their share of the population; the weights correct for that. Reporting the unweighted figure as a population estimate would overstate how badly children are doing. unweighted() is for sample descriptions, quality checks, and reporting how many observations an estimate rests on, never for population claims.
4.10.2 Restricting to a subpopulation
Estimates for part of the population are produced by filter(), applied to the design object rather than the data. The design keeps the information about the excluded cases that the variance calculation needs, which is why the order matters.
enrolled_mpl <- icanicar_des |>
filter(EnrolmentStatus == "Currently Enrolled") |>
summarize(p = survey_mean(MPLMath, na.rm = TRUE, vartype = "ci"))
all_mpl <- icanicar_des |>
summarize(p = survey_mean(MPLMath, na.rm = TRUE, vartype = "ci"))
bind_rows(
enrolled_mpl |> mutate(population = "Currently enrolled"),
all_mpl |> mutate(population = "All children")
) |>
relocate(population)## # A tibble: 2 × 4
## population p p_low p_upp
## <chr> <dbl> <dbl> <dbl>
## 1 Currently enrolled 0.530 0.520 0.540
## 2 All children 0.507 0.498 0.517Mathematics proficiency is 53.0% among enrolled children and 50.7% among all children. The difference is modest because most children in the survey are enrolled, but the two figures answer different questions and should never be swapped: one describes what schooling is delivering, the other describes the whole cohort including children school never reached.
4.10.3 Design effects and effective sample size
Clustered samples are cheaper to collect than simple random samples but carry less information per observation, because children in the same village resemble each other. The design effect quantifies that loss: it is the variance of an estimate under the actual design divided by its variance under simple random sampling without replacement. A value above 1 says the design buys less precision than simple random sampling would for the same number of observations; a value below 1, which is uncommon but does happen when the strata line up closely with the outcome, says it buys more. Dividing the sample size by the design effect gives the effective sample size,
\[n_{eff} = \frac{n}{D_{eff}},\]
the size of the simple random sample that would have produced equally precise estimates. Setting deff = TRUE requests it:
math_deff <- icanicar_des |>
filter(!is.na(MathIRTScore)) |>
summarize(
score = survey_mean(MathIRTScore, deff = TRUE),
children = unweighted(length(MathIRTScore))
) |>
mutate(n_effective = children / score_deff)
math_deff## # A tibble: 1 × 5
## score score_se score_deff children n_effective
## <dbl> <dbl> <dbl> <int> <dbl>
## 1 -0.263 0.0128 13.1 89138 6826.The design effect for the mean mathematics score is 13.1, so the 89,138 assessed children carry about as much information as a simple random sample of 6,826. That is a large effect, and it has a straightforward explanation: with roughly 20 households per enumeration area, and learning outcomes strongly clustered by community and school, neighbouring children supply overlapping information. It also explains why ignoring the design would be so misleading here, since standard errors computed as if the sample were random would be far too small.
Design effects are specific to each outcome. Variables that vary within a village more than between villages, such as the child’s sex, show much weaker clustering than variables shaped by the local school and language environment.
4.10.4 Adding summary rows with cascade()
group_by() followed by summarize() gives one row per group and nothing else. Reports usually want a total row alongside the groups. cascade() replaces summarize() and adds it:
.data is the design object and ... takes the same name-value pairs as summarize(). .fill supplies the label for the summary row, NA by default. .fill_level_top moves that row to the top of the output, which works when the grouping variable is a factor.
icanicar_des |>
group_by(residence) |>
cascade(
score = survey_mean(MathIRTScore, na.rm = TRUE),
children = survey_total(MPLBoth, na.rm = TRUE, vartype = NULL),
.fill = "All areas"
)## # A tibble: 3 × 4
## residence score score_se children
## <fct> <dbl> <dbl> <dbl>
## 1 Rural -0.358 0.0190 45176968.
## 2 Urban -0.145 0.0186 50726065.
## 3 All areas -0.263 0.0128 95903033.The final row is the overall estimate, not a third residence category, and giving it a label such as “All areas” through .fill prevents readers from misreading an NA as missing data.
4.10.5 Estimating several outcomes at once
Reporting rarely stops at one indicator. across() applies the same estimator to several variables inside a single summarize(), and .unpack controls how the resulting columns are named:
mpl_wide <- icanicar_des |>
summarize(across(
c(MPLReading, MPLMath, MPLBoth),
~ survey_mean(.x, na.rm = TRUE, vartype = "ci"),
.unpack = "{outer}.{inner}"
))
mpl_wide## # A tibble: 1 × 9
## MPLReading.coef MPLReading._low MPLReading._upp MPLMath.coef MPLMath._low MPLMath._upp MPLBoth.coef MPLBoth._low MPLBoth._upp
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.378 0.366 0.389 0.507 0.498 0.517 0.344 0.333 0.355One wide row is convenient for the computer and awkward for a reader. pivot_longer() turns it into one row per indicator, using a regular expression to split each column name at the dot:
mpl_long <- mpl_wide |>
pivot_longer(
everything(),
names_to = c("indicator", ".value"),
names_pattern = "(.*)\\.(.*)"
) |>
rename(proportion = coef, lower = `_low`, upper = `_upp`) |>
mutate(indicator = recode(
indicator,
MPLReading = "Reading",
MPLMath = "Mathematics",
MPLBoth = "Both subjects"
))
mpl_long## # A tibble: 3 × 4
## indicator proportion lower upper
## <chr> <dbl> <dbl> <dbl>
## 1 Reading 0.378 0.366 0.389
## 2 Mathematics 0.507 0.498 0.517
## 3 Both subjects 0.344 0.333 0.355The table now reads as a set of findings: 37.8% of children reach the reading threshold, 50.7% the mathematics threshold, and 34.4% both. Turning a table like this into publication-ready output is what packages such as {gt} and {gtsummary}, introduced in Chapter ??, are for.
across() does not work with survey_prop(), which takes its variable from group_by() rather than as an argument. When the same breakdown is needed for several categorical variables, write a small function and map over the variable names instead:
mpl_split_by <- function(variable) {
icanicar_des |>
filter(!is.na(.data[[variable]]), !is.na(MPLMath)) |>
group_by(.data[[variable]]) |>
summarize(p = survey_mean(MPLMath, na.rm = TRUE, vartype = "ci")) |>
rename(category = 1) |>
mutate(variable = variable, category = as.character(category)) |>
relocate(variable)
}
c("sex", "residence", "schooling") |>
map(mpl_split_by) |>
list_rbind()## # A tibble: 7 × 5
## variable category p p_low p_upp
## <chr> <chr> <dbl> <dbl> <dbl>
## 1 sex Female 0.499 0.488 0.509
## 2 sex Male 0.515 0.503 0.528
## 3 residence Rural 0.464 0.450 0.478
## 4 residence Urban 0.561 0.545 0.576
## 5 schooling Currently Enrolled 0.530 0.520 0.540
## 6 schooling Out of School 0.335 0.295 0.375
## 7 schooling Never Enrolled 0.0816 0.0597 0.103The function is written once, the mapping extends to as many variables as needed, and list_rbind() stacks the pieces into a single tidy table. Because the grouping variable is supplied as a string, .data[[variable]] is used to look it up, and rename(category = 1) gives the first column a common name so the results can be stacked. This pattern scales to the dozens of background variables a full descriptive report usually needs.
4.11 Exercises
The exercises use icanicar_des from the Prerequisites box.
- How many children in the population have never been enrolled in school? Report a 95% confidence interval. Hint:
schoolingdistinguishes the three enrolment categories. - What proportion of girls and of boys reach the minimum proficiency level in reading? Hint: group by
sexand usempl_reading, or take the mean ofMPLReading. - Estimate the median reading score for rural and urban children, with confidence intervals. Which group has the wider interval, and why might that be?
- Produce a table of the estimated number of children currently enrolled in each country, sorted from largest to smallest, with a summary row for all countries. Hint:
cascade(). - Estimate the design effect attached to the reading proficiency proportion, convert it into an effective sample size, and compare it with the design effect attached to the mean reading score. Which of the two loses more precision to clustering?
- Is the time spent on the numeracy assessment related to the time spent on the reading assessment? Estimate the correlation between
IcanAssessTimeandIcarAssessTime, then repeat it separately for children assessed in their home language and children who were not (AssHomeLang). - Among children who are currently enrolled, estimate the proportion reaching the minimum proficiency level in mathematics for each single year of age from 7 to 12. Describe the pattern in two or three sentences. Hint:
filter()the design first, then group bych02. - Compare the weighted and unweighted proportions of children living in urban areas. What does the difference tell you about how the sample was drawn?
4.12 Further reading
The {srvyr} reference documentation is the authoritative source on the functions used in this chapter, and each help page lists the {survey} function it calls. Zimmer, Powell, and Velásquez-Hernández’s Exploring Complex Survey Data Analysis Using R (tidy-survey-r.github.io/tidy-survey-book) works through the same estimators using United States survey data and goes further into replicate-weight designs. For the sampling design, weights, and questionnaire content behind the estimates above, see the PAL Network documentation collected in the DataFirst metadata record for ICAN-ICAR 2025.