6 Analyzing AMPLab data in R

This section provides a practical introduction to analysing AMPLab data in R using the Rrepest package. It is intended for users with a basic working knowledge of R and includes step-by-step examples covering common AMPLab analyses.

Use the contents below to navigate directly to the topics most relevant to your analysis.

6.1 Package installation

To get started, you need to install the required packages once. In the R console, type:

install.packages(c(
  "Rrepest", "haven", "survey", "mitools", "quantreg", "tidyverse", "gt"
))

The haven package is needed to read Stata .dta files.
The Rrepest package automates the handling of plausible values and replicate weights, making it easier to analyse AMPLab data correctly in R.
The survey and mitools packages are used only for the alternative analyses shown in this chapter.
The tidyverse package is used for data wrangling.
The gt package is used to present selected results in readable tables.

You only need to run install.packages() once; there is no need to reinstall each time you open R.

# Load packages at the start of every session

library(haven)
library(Rrepest)
library(survey)
library(mitools)
library(tidyverse)
library(gt)

6.2 Loading data in R

For hands-on examples we’ll start with the student data.

# Paths below are relative to the project root

ampl_dat <- read_dta("data/ampl-student-v1.dta")

It is good practice to run your analysis from an R script rather than the console, as this makes your work reproducible. All examples in this section are written as script code.

6.3 Add variable and value labels

We have supplied an R script that applies variable and value labels to the data. These files are in the data/ folder of this project.

# Apply variable and value labels
source("data/ampl-student-labels.R")
ampl_dat <- apply_ampl_stq_labels(ampl_dat)

After loading the file, let’s check if we have the most important variables that Rrepest uses. The five plausible values for language and mathematics are pv1_read to pv5_read and pv1_math to pv5_math, respectively. The final weight is fwgt. The replicate weights are rwgt1 to rwgt135 and countries are indicated by cnt.

# Confirm key variables exist
ampl_dat |>
  select(matches("^pv[1-5]_(read|math)$"), fwgt, rwgt1:rwgt10, cnt) |>
  glimpse()

## Rows: 22,214
## Columns: 22
## $ pv1_math <dbl> 0.22811525, -2.62520194, -2.76957273, -0.08613674, 0.13072035, 0.47429678, -0.70883948, -0.74843246, 0.32783359, 0.06470735, 1.34257…
## $ pv2_math <dbl> 0.462508053, -1.674260497, -3.285824299, 0.678741813, 0.057586100, -0.041360527, -0.605331600, -0.487574786, -0.000326509, -0.026982…
## $ pv3_math <dbl> 0.89248133, -1.76391876, -3.35175014, 0.95755553, 0.16893062, -0.08748718, -0.18168271, -1.29100907, 0.41965342, 0.42386377, 1.39329…
## $ pv4_math <dbl> 0.60037166, -2.21478987, -3.26386952, -0.47502181, -0.19398750, 0.01470500, -0.34860367, -1.24189603, 0.35222486, 0.07257382, 1.6342…
## $ pv5_math <dbl> 0.76853716, -1.63746548, -3.26517296, 0.27711055, 0.85098785, -0.25193688, 0.37577242, -1.31283343, 0.61149669, 0.60424227, 0.967157…
## $ pv1_read <dbl> 0.768692255, -1.638499618, -2.526624203, 0.606318355, 0.767933547, 0.678709507, 0.351070166, -0.241998076, 1.266379833, 0.826490045,…
## $ pv2_read <dbl> 0.57991463, -0.88148451, -2.35456824, 1.11718643, 0.89872378, 0.64056826, 0.09324962, -0.18746106, 1.17088628, 0.44311303, 1.9112074…
## $ pv3_read <dbl> 1.05408597, -1.38343024, -2.53575850, 1.52389026, 0.83753365, 0.45622772, 0.58972800, -0.49328646, 1.05618250, 1.15197921, 1.9779125…
## $ pv4_read <dbl> 1.198313951, -1.116216302, -2.667570114, 0.841243982, 0.547337115, 0.909113705, -0.142348379, -0.111854941, 0.898920715, 1.050507784…
## $ pv5_read <dbl> 0.94776469, -1.12560523, -2.62032199, 1.02373171, 0.85645103, 0.02454279, 0.43081874, -0.22813433, 1.29445171, 1.11652923, 1.5622578…
## $ fwgt     <dbl> 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 24…
## $ rwgt1    <dbl> 495.3116, 495.3116, 495.3116, 495.3116, 495.3116, 495.3116, 495.3116, 495.3116, 495.3116, 495.3116, 495.3116, 495.3116, 495.3116, 49…
## $ rwgt2    <dbl> 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 24…
## $ rwgt3    <dbl> 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 24…
## $ rwgt4    <dbl> 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 24…
## $ rwgt5    <dbl> 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 24…
## $ rwgt6    <dbl> 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 24…
## $ rwgt7    <dbl> 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 24…
## $ rwgt8    <dbl> 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 24…
## $ rwgt9    <dbl> 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 24…
## $ rwgt10   <dbl> 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 247.6558, 24…
## $ cnt      <chr> "KEN", "KEN", "KEN", "KEN", "KEN", "KEN", "KEN", "KEN", "KEN", "KEN", "KEN", "KEN", "KEN", "KEN", "KEN", "KEN", "KEN", "KEN", "KEN",…

All required variables for Rrepest are available in the data.

6.4 The `Rrepest` package

Rrepest automates the handling of plausible values and replicate weights. Once the survey settings have been specified, you can use it to estimate means, summary statistics, proficiency levels, group differences, percentiles, and regression models while correctly accounting for both plausible values and replicate weights.

6.5 `Rrepest` command syntax

The basic syntax of the Rrepest() function is as follows:

Rrepest(
  data       = <data frame or subset>,
  svy        = <survey name>,
  est        = est(<statistic>, target = <variable(s)>, regressor = <variable(s)>),
  by         = <grouping variable(s)>,  # separate analysis per group
  over       = <subgroup variable(s)>,  # within-group comparison
  cm.weights = c(<final weight>, <replicate weights>),
  var.factor = <variance scaling constant>,
  n.pvs      = <number of plausible values>
)

6.6 Before you begin: Set up `Rrepest` for AMPLab

AMPLab is not one of the studies supported by the Rrepest package, meaning that there aren’t built-in survey specifications. Therefore, before running any analyses, you need to tell Rrepest how the AMPLab data are structured.

You will need to use svy = "SVY" and supply the survey parameters directly. To analyse AMPLab data correctly, Rrepest needs to know the following survey specifications:

AMPLab uses the paired Jackknife method for creating the replicate weights (JK2).
There are 135 replicate weights in the data.
The final weight is given by fwgt.
There are five sets of plausible values for each of mathematics, pv1_math to pv5_math, and reading, pv1_read to pv5_read.

Define the weight vector once and reuse it throughout your analyses:

# Weight helpers — define once, use everywhere
cm <- c("fwgt", paste0("rwgt", 1:135))   # final + 135 replicate weights

We will often present {Rrepest} output with the {gt} package. The helper below reshapes linear-regression results into a tidy table of estimates, standard errors, and \(t\)-ratios.

You do not need to understand this as it is for improving readability of the {Rrepest} output.

# Reshape Rrepest lm output for gt presentation
# b_cols / se_cols: character vectors of column names in the Rrepest result
# terms: labels shown in the table (same order as b_cols / se_cols)

rrepest_lm_table <- function(result, b_cols, se_cols, terms) {
  df <- as.data.frame(result)
  tibble(
    term = terms,
    estimate = as.numeric(unlist(df[1, b_cols], use.names = FALSE)),
    se = as.numeric(unlist(df[1, se_cols], use.names = FALSE))
  ) |>
    mutate(
      t_value = estimate / se,
      estimate_se = sprintf("%.3f (%.3f)", estimate, se)
    )
}

show_lm_gt <- function(tbl, title, subtitle, source_note = NULL) {
  out <- tbl |>
    select(term, estimate_se, t_value) |>
    gt() |>
    tab_header(title = title, subtitle = subtitle) |>
    cols_label(
      term = "Term",
      estimate_se = "Estimate (SE)",
      t_value = "t"
    ) |>
    fmt_number(columns = t_value, decimals = 2) |>
    tab_style(
      style = cell_text(weight = "bold"),
      locations = cells_column_labels()
    )
  if (!is.null(source_note)) {
    out <- out |> tab_source_note(source_note = source_note)
  }
  out
}

Commands for analysing the AMPLab student data will have the following syntax:

Rrepest(
  data       = subset(ampl_dat, <condition>),
  svy        = "SVY",
  est        = est(<statistic>, target = <variable>),
  cm.weights = cm,
  var.factor = 135,
  n.pvs      = 5
)

6.7 AMPLab Analyses Examples

Before we start analysing proficiency scores, let us get a sense of the countries and grades presented in the data.

ampl_dat |>
  count(cnt, grade_stf)

## # A tibble: 5 × 3
##   cnt   grade_stf     n
##   <chr>     <dbl> <int>
## 1 GMB           3  4106
## 2 KEN           6  5238
## 3 LSO           7  3917
## 4 ZMB           4  4474
## 5 ZMB           7  4479

6.7.1 Calculating Mean Age by Country

Let’s start by calculating the mean age of students in Grade 3 in The Gambia. The age variable is called s_age.

ampl_dat |>
  summarise(
    n = n(),
    min_age = min(s_age, na.rm = TRUE),
    max_age = max(s_age, na.rm = TRUE),
    n_missing_code = sum(s_age == 9999, na.rm = TRUE)
  )

## # A tibble: 1 × 4
##       n min_age   max_age        n_missing_code
##   <int> <dbl+lbl> <dbl+lbl>               <int>
## 1 22214 6         9999 [Missing]            172

Note that ages are not integers, so treating age as a short list of categories would not be appropriate. Missing age values are coded as \(9999\). Because \(9999\) is a placeholder used to indicate missing data rather than a student’s actual age, these observations should be excluded from analyses involving age. Failure to do so may produce misleading results.

To calculate the average age of students in Grade 3 in The Gambia, the syntax is as follows:

# Calculating average age in The Gambia

result_age <- Rrepest(
  data       = subset(ampl_dat, cnt == "GMB" & s_age < 9999),
  svy        = "SVY",
  est        = est("mean", target = "s_age"),
  cm.weights = cm,
  var.factor = 135
)
print(result_age)

## # A tibble: 1 × 3
##   .     b.mean.s_age se.mean.s_age
##   <chr>        <dbl>         <dbl>
## 1 Total         10.7        0.0335

We present the {Rrepest} output above in a reader-friendly format using the {gt} package.

as.data.frame(result_age) |>
  transmute(
    estimate = `b.mean.s_age`,
    se = `se.mean.s_age`
  ) |>
  gt() |>
  tab_header(
    title = "Mean student age",
    subtitle = "The Gambia (Grade 3)"
  ) |>
  cols_label(estimate = "Mean age", se = "SE") |>
  fmt_number(columns = c(estimate, se), decimals = 3) |>
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_column_labels()
  )

Mean age	SE
Mean student age
The Gambia (Grade 3)
10.744	0.034

Because age is an observed variable rather than a plausible-value variable, Rrepest uses only the sampling variance derived from the replicate weights. The reported standard error therefore reflects uncertainty arising from the sample design but not measurement uncertainty.

The estimated mean age of Grade 3 students in The Gambia is \(10.74\) years (\(95\%\) CI: \(10.67–10.81\)). This reflects the substantial grade repetition and late entry common in the region, with many students aged \(11–23\) also enrolled in Grade 3.

Alternative to generate the same result using survey

AMPLab provides replicate weights specifically so that users can reproduce the official variance estimation procedure. The survey package in R can be configured to use the paired jackknife replicate weights.

The syntax to generate the same result is:

options(survey.lonely.psu = "adjust")

des_age <- svrepdesign(
  weights          = ~fwgt,
  repweights       = as.matrix(ampl_dat[, paste0("rwgt", 1:135)]),
  type             = "other",
  scale            = 1,
  rscales          = rep(1, 135),
  combined.weights = TRUE,
  mse              = TRUE,
  data             = ampl_dat
)

svymean(~s_age, subset(des_age, cnt == "GMB" & s_age < 9999), na.rm = TRUE)

##         mean     SE
## s_age 10.744 0.0335

6.7.2 Calculating mean mathematics proficiency

Now we will use \(5\) plausible values to estimate mean mathematics scores in The Gambia.

Plausible values are a set of multiple imputations. The Rrepest package automatically recognises plausible values when the variable name contains the @ symbol.

For example:

est("mean", target = "pv@_read")

tells Rrepest to analyse all five reading plausible values, combine results appropriately, and calculate standard errors that reflect both sampling and measurement uncertainty. This allows researchers to obtain valid estimates without having to implement the multiple-imputation calculations manually.

result_math <- Rrepest(
  data       = subset(ampl_dat, cnt == "GMB"),
  svy        = "SVY",
  est        = est("mean", target = "pv@_math"),
  cm.weights = cm,
  var.factor = 135,
  n.pvs      = 5
)
print(result_math)

## # A tibble: 1 × 3
##   .     `b.mean.pv@_math` `se.mean.pv@_math`
##   <chr>             <dbl>              <dbl>
## 1 Total             -2.34             0.0364

We present the {Rrepest} output above in a reader-friendly format using the {gt} package.

as.data.frame(result_math) |>
  transmute(
    estimate = `b.mean.pv@_math`,
    se = `se.mean.pv@_math`
  ) |>
  gt() |>
  tab_header(
    title = "Mean mathematics proficiency",
    subtitle = "The Gambia (Grade 3) — Learning Progressions Scale"
  ) |>
  cols_label(estimate = "Mean", se = "SE") |>
  fmt_number(columns = c(estimate, se), decimals = 3) |>
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_column_labels()
  )

Mean	SE
Mean mathematics proficiency
The Gambia (Grade 3) — Learning Progressions Scale
−2.339	0.036

The estimated mean mathematics proficiency in The Gambia is \(-2.34\) on the Learning Progressions Scale (\(95\%\) CI: \(-2.41\) to \(-2.27\); SE = \(0.036\)). Scores on this scale are centred near the AMPL proficiency cut scores, so a negative mean indicates that Grade 3 students in The Gambia are, on average, below the end-of-primary MPLb threshold. The standard error reflects both sampling uncertainty and uncertainty arising from the plausible values.

Unlike age, mathematics proficiency is represented by five plausible values. Rrepest estimates the mean separately for each plausible value, combines the five estimates using multiple-imputation formulas, and then incorporates the replicate-weight variance to produce the final standard error.

An alternative approach using the survey and mitools packages reproduces the same result by looping over plausible values manually.

options(survey.lonely.psu = "adjust")

des_math <- svrepdesign(
  weights          = ~fwgt,
  repweights       = as.matrix(ampl_dat[, paste0("rwgt", 1:135)]),
  type             = "other",
  scale            = 1,
  rscales          = rep(1, 135),
  combined.weights = TRUE,
  mse              = TRUE,
  data             = ampl_dat
)

pv_math <- withPV(
  mapping = list(math ~ pv1_math + pv2_math + pv3_math + pv4_math + pv5_math),
  data    = subset(des_math, cnt == "GMB"),
  action  = function(dsgn) svymean(~math, dsgn, na.rm = TRUE),
  rewrite = TRUE
)
summary(MIcombine(pv_math))

## Multiple imputation results:
##       withPV.svyrep.design(mapping = list(math ~ pv1_math + pv2_math + 
##     pv3_math + pv4_math + pv5_math), data = subset(des_math, 
##     cnt == "GMB"), action = function(dsgn) svymean(~math, dsgn, 
##     na.rm = TRUE), rewrite = TRUE)
##       MIcombine.default(pv_math)
##            results         se    (lower    upper) missInfo
## pv1_math -2.339032 0.03639078 -2.410366 -2.267699      2 %

6.7.3 Proficiency across gender

We can go further and ask how mean proficiency compares between girls and boys within The Gambia. To do this we first create an indicator variable for girls (1 = girl, 0 = boy).

ampl_dat |>
  count(s_gender)

## # A tibble: 2 × 2
##   s_gender       n
##   <dbl+lbl>  <int>
## 1 1 [Female] 11592
## 2 2 [Male]   10622

ampl_dat <- ampl_dat |>
  mutate(
    girl = if_else(s_gender %in% c(1, 2), as.integer(s_gender == 1), NA_integer_)
  )

With the girl indicator created, we can compute mean mathematics and reading proficiency using the by argument.

result_girl <- Rrepest(
  data       = subset(ampl_dat, cnt == "GMB"),
  svy        = "SVY",
  est        = est("mean", target = c("pv@_math", "pv@_read")),
  by         = "girl",
  cm.weights = cm,
  var.factor = 135,
  n.pvs      = 5
)
print(result_girl)

## # A tibble: 2 × 5
##   girl  `b.mean.pv@_math` `se.mean.pv@_math` `b.mean.pv@_read` `se.mean.pv@_read`
##   <chr>             <dbl>              <dbl>             <dbl>              <dbl>
## 1 0                 -2.37             0.0393             -1.40             0.0314
## 2 1                 -2.31             0.0459             -1.33             0.0363

We present the {Rrepest} output above in a reader-friendly format using the {gt} package.

as.data.frame(result_girl) |>
  mutate(
    sex = if_else(as.integer(girl) == 1, "Girls", "Boys"),
    math_b = `b.mean.pv@_math`,
    math_se = `se.mean.pv@_math`,
    read_b = `b.mean.pv@_read`,
    read_se = `se.mean.pv@_read`
  ) |>
  select(sex, math_b, math_se, read_b, read_se) |>
  gt(rowname_col = "sex") |>
  tab_header(
    title = "Mean proficiency by gender",
    subtitle = "The Gambia (Grade 3)"
  ) |>
  tab_spanner(label = "Mathematics", columns = c(math_b, math_se)) |>
  tab_spanner(label = "Reading", columns = c(read_b, read_se)) |>
  cols_label(
    math_b = "Mean",
    math_se = "SE",
    read_b = "Mean",
    read_se = "SE"
  ) |>
  fmt_number(columns = c(math_b, math_se, read_b, read_se), decimals = 3) |>
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_column_labels()
  )

	Mathematics		Reading
Mean proficiency by gender
The Gambia (Grade 3)
	Mean	SE	Mean	SE
Boys	−2.372	0.039	−1.397	0.031
Girls	−2.312	0.046	−1.328	0.036

In The Gambia, girls score slightly higher than boys in both domains. Mean mathematics is \(-2.31\) for girls (SE \(= 0.046\)) versus \(-2.37\) for boys (SE \(= 0.039\)); mean reading is \(-1.33\) for girls (SE \(= 0.036\)) versus \(-1.40\) for boys (SE \(= 0.031\)). The gaps are modest relative to the overall low level of performance on the LPS, but the pattern is consistent across subjects.

6.7.4 Proficiency by Country

Now let us compare mathematics and reading proficiency across countries. Zambia administered AMPLab in two grades. In R we can group by both cnt and grade (or create a combined cnt_grade label) so that Zambia’s Grade 4 and Grade 7 samples are not averaged together.

ampl_dat <- ampl_dat |>
  mutate(
    grade_num = as.integer(haven::zap_labels(grade_stf)),
    cnt_grade = paste(cnt, grade_num)
  )

ampl_dat |>
  count(cnt_grade)

## # A tibble: 5 × 2
##   cnt_grade     n
##   <chr>     <int>
## 1 GMB 3      4106
## 2 KEN 6      5238
## 3 LSO 7      3917
## 4 ZMB 4      4474
## 5 ZMB 7      4479

We can now compute proficiency scores by country (and by grade for Zambia).

result_cntgrade <- Rrepest(
  data       = ampl_dat,
  svy        = "SVY",
  est        = est("mean", target = c("pv@_math", "pv@_read")),
  by         = "cnt_grade",
  cm.weights = cm,
  var.factor = 135,
  n.pvs      = 5
)
print(result_cntgrade)

## # A tibble: 5 × 5
##   cnt_grade `b.mean.pv@_math` `se.mean.pv@_math` `b.mean.pv@_read` `se.mean.pv@_read`
##   <chr>                 <dbl>              <dbl>             <dbl>              <dbl>
## 1 GMB 3                -2.34              0.0364            -1.36              0.0298
## 2 KEN 6                -0.438             0.0339             0.156             0.0390
## 3 LSO 7                -0.853             0.0434            -0.391             0.0494
## 4 ZMB 4                -2.67              0.0352            -1.65              0.0318
## 5 ZMB 7                -1.05              0.0374            -0.523             0.0418

We present the {Rrepest} output above in a reader-friendly format using the {gt} package.

as.data.frame(result_cntgrade) |>
  mutate(
    country_grade = cnt_grade,
    math_b = `b.mean.pv@_math`,
    math_se = `se.mean.pv@_math`,
    read_b = `b.mean.pv@_read`,
    read_se = `se.mean.pv@_read`
  ) |>
  select(country_grade, math_b, math_se, read_b, read_se) |>
  gt(rowname_col = "country_grade") |>
  tab_header(
    title = "Mean proficiency by country and grade",
    subtitle = "Learning Progressions Scale"
  ) |>
  tab_spanner(label = "Mathematics", columns = c(math_b, math_se)) |>
  tab_spanner(label = "Reading", columns = c(read_b, read_se)) |>
  cols_label(
    math_b = "Mean",
    math_se = "SE",
    read_b = "Mean",
    read_se = "SE"
  ) |>
  fmt_number(columns = c(math_b, math_se, read_b, read_se), decimals = 3) |>
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_column_labels()
  )

	Mathematics		Reading
Mean proficiency by country and grade
Learning Progressions Scale
	Mean	SE	Mean	SE
GMB 3	−2.339	0.036	−1.359	0.030
KEN 6	−0.438	0.034	0.156	0.039
LSO 7	−0.853	0.043	−0.391	0.049
ZMB 4	−2.666	0.035	−1.648	0.032
ZMB 7	−1.050	0.037	−0.523	0.042

Mean performance differs sharply by country–grade group. On mathematics, Kenya Grade 6 has the highest mean (\(-0.44\)), followed by Lesotho Grade 7 (\(-0.85\)) and Zambia Grade 7 (\(-1.05\)), while Zambia Grade 4 (\(-2.67\)) and The Gambia Grade 3 (\(-2.34\)) are much lower. Reading follows a similar ordering, with Kenya the only group whose mean reading score is positive (\(0.16\)). These contrasts partly reflect grade level as well as system differences, which is why separating Zambia’s Grade 4 and Grade 7 samples matters.

6.7.5 Moving beyond means: Summary statistics with `Rrepest`

Rrepest can estimate a range of statistics beyond the mean, such as percentiles and standard deviations. Below we summarise the mathematics proficiency distribution for Grade 3 students in The Gambia by gender.

result_summ <- Rrepest(
  data       = subset(ampl_dat, cnt == "GMB"),
  svy        = "SVY",
  est        = est(
    c("mean", "std", "quant", 0.05, "quant", 0.25, "quant", 0.5, "quant", 0.75, "quant", 0.95),
    target = "pv@_math"
  ),
  by         = "girl",
  cm.weights = cm,
  var.factor = 135,
  n.pvs      = 5
)
print(as.data.frame(result_summ))

##   girl b.mean.pv@_math se.mean.pv@_math b.std.pv@_math se.std.pv@_math b.quant00.pv@_math se.quant00.pv@_math b.quant025.pv@_math se.quant025.pv@_math
## 1    0       -2.371841       0.03931614       1.065351      0.03348136          -3.850363          0.03842642           -3.121192           0.04034239
## 2    1       -2.311790       0.04594274       1.090902      0.03767271          -3.820493          0.04860671           -3.092192           0.04298113
##   b.quant05.pv@_math se.quant05.pv@_math b.quant075.pv@_math se.quant075.pv@_math b.quant095.pv@_math se.quant095.pv@_math
## 1          -2.522535          0.06196889           -1.790534           0.08457232          -0.3191085            0.1368639
## 2          -2.481752          0.05041328           -1.675401           0.07123222          -0.2163765            0.1640926

We present the {Rrepest} output above in a reader-friendly format using the {gt} package.

as.data.frame(result_summ) |>
  mutate(
    sex = if_else(as.integer(girl) == 1, "Girls", "Boys"),
    mean_b = `b.mean.pv@_math`,
    mean_se = `se.mean.pv@_math`,
    sd_b = `b.std.pv@_math`,
    sd_se = `se.std.pv@_math`,
    p50_b = `b.quant05.pv@_math`,
    p50_se = `se.quant05.pv@_math`,
    p75_b = `b.quant075.pv@_math`,
    p75_se = `se.quant075.pv@_math`,
    p95_b = `b.quant095.pv@_math`,
    p95_se = `se.quant095.pv@_math`
  ) |>
  select(sex, mean_b, mean_se, sd_b, sd_se, p50_b, p50_se, p75_b, p75_se, p95_b, p95_se) |>
  gt(rowname_col = "sex") |>
  tab_header(
    title = "Mathematics distribution by gender",
    subtitle = "The Gambia (Grade 3)"
  ) |>
  tab_spanner(label = "Mean", columns = c(mean_b, mean_se)) |>
  tab_spanner(label = "SD", columns = c(sd_b, sd_se)) |>
  tab_spanner(label = "Median", columns = c(p50_b, p50_se)) |>
  tab_spanner(label = "75th pct", columns = c(p75_b, p75_se)) |>
  tab_spanner(label = "95th pct", columns = c(p95_b, p95_se)) |>
  cols_label(
    mean_b = "Est.", mean_se = "SE",
    sd_b = "Est.", sd_se = "SE",
    p50_b = "Est.", p50_se = "SE",
    p75_b = "Est.", p75_se = "SE",
    p95_b = "Est.", p95_se = "SE"
  ) |>
  fmt_number(decimals = 3) |>
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_column_labels()
  )

	Mean		SD		Median		75th pct		95th pct
Mathematics distribution by gender
The Gambia (Grade 3)
	Est.	SE	Est.	SE	Est.	SE	Est.	SE	Est.	SE
Boys	−2.372	0.039	1.065	0.033	−2.523	0.062	−1.791	0.085	−0.319	0.137
Girls	−2.312	0.046	1.091	0.038	−2.482	0.050	−1.675	0.071	−0.216	0.164

Looking beyond average performance shows that the gender gap is small across the distribution. Boys and girls have similar spread (SD ≈ \(1.07\) vs \(1.09\)). At the median, girls score \(-2.48\) and boys \(-2.52\); at the 75th percentile the gap is a little larger (\(-1.68\) vs \(-1.79\)). Even the strongest performers remain well below zero on the LPS: the 95th percentile is about \(-0.22\) for girls and \(-0.32\) for boys.

To examine the shape of the distribution more closely, we can also estimate selected quantile cut-points (here the 20th, 40th, 60th, and 80th percentiles), which play a similar role to Stata’s quantiletable quintile cuts.

result_quant <- Rrepest(
  data       = subset(ampl_dat, cnt == "GMB"),
  svy        = "SVY",
  est        = est(
    c("quant", 0.2, "quant", 0.4, "quant", 0.6, "quant", 0.8),
    target = c("pv@_math", "pv@_read")
  ),
  by         = "girl",
  cm.weights = cm,
  var.factor = 135,
  n.pvs      = 5
)
print(as.data.frame(result_quant))

##   girl b.quant02.pv@_math se.quant02.pv@_math b.quant04.pv@_math se.quant04.pv@_math b.quant06.pv@_math se.quant06.pv@_math b.quant08.pv@_math
## 1    0          -3.248062          0.04803454          -2.769858          0.04997564          -2.262513          0.04622875          -1.538429
## 2    1          -3.219790          0.04848990          -2.735343          0.05043019          -2.203289          0.06672874          -1.434276
##   se.quant08.pv@_math b.quant02.pv@_read se.quant02.pv@_read b.quant04.pv@_read se.quant04.pv@_read b.quant06.pv@_read se.quant06.pv@_read
## 1          0.08995220          -2.047014          0.03274443          -1.702325          0.02887452          -1.340615          0.03112681
## 2          0.06044208          -2.024973          0.03415830          -1.640170          0.03633984          -1.266828          0.04618761
##   b.quant08.pv@_read se.quant08.pv@_read
## 1         -0.7967462          0.06035399
## 2         -0.6501527          0.06956877

We present the {Rrepest} output above in a reader-friendly format using the {gt} package.

as.data.frame(result_quant) |>
  mutate(
    sex = if_else(as.integer(girl) == 1, "Girls", "Boys"),
    math20 = `b.quant02.pv@_math`,
    math80 = `b.quant08.pv@_math`,
    read20 = `b.quant02.pv@_read`,
    read80 = `b.quant08.pv@_read`
  ) |>
  select(sex, math20, math80, read20, read80) |>
  gt(rowname_col = "sex") |>
  tab_header(
    title = "Selected proficiency quantiles by gender",
    subtitle = "The Gambia (Grade 3) - 20th and 80th percentiles"
  ) |>
  tab_spanner(label = "Mathematics", columns = c(math20, math80)) |>
  tab_spanner(label = "Reading", columns = c(read20, read80)) |>
  cols_label(
    math20 = "20th",
    math80 = "80th",
    read20 = "20th",
    read80 = "80th"
  ) |>
  fmt_number(decimals = 3) |>
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_column_labels()
  )

	Mathematics		Reading
Selected proficiency quantiles by gender
The Gambia (Grade 3) - 20th and 80th percentiles
	20th	80th	20th	80th
Boys	−3.248	−1.538	−2.047	−0.797
Girls	−3.220	−1.434	−2.025	−0.650

These cut-points summarise ordered parts of the score distribution. For example, among boys the mathematics 20th percentile is about \(-3.25\) and the 80th percentile about \(-1.54\); girls’ corresponding cut-points are slightly higher (about \(-3.22\) and \(-1.43\)), which matches the modest advantage for girls seen in the means and percentiles.

Other useful Rrepest statistics include "freq", "corr", and "lm".

6.7.6 Estimating AMPL Proficiency Levels

AMPL established cut scores that correspond to the minimum proficiency levels required at the end of lower primary school (MPLa) and at the end of primary school (MPLb), for both mathematics and reading.

The cut scores for mathematics are defined by AMPL as follows:

AMPL Proficiency Level	Score
Below MPLa	\(< -1.74\)
Between MPLa & MPLb	\([-1.74, -0.04)\)
Above MPLb	\(\geq -0.04\)

These cut scores for reading and mathematics were established on the Learning Progressions Scale (LPS) with an international standard setting exercise undertaken in 2022. A pairwise comparison method (PCM) study was then used to confirm and validate the locations of those cut scores.

There is a variable for each plausible value that records this cut-score information for both mathematics (pl_pv*_math) and reading (pl_pv*_read). The file also supplies ready-made binary indicators for reaching MPLa (pl_psa*_math, pl_psa*_read) and MPLb (pl_psb*_math, pl_psb*_read).

ampl_dat |>
  count(pl_pv1_read)

## # A tibble: 3 × 2
##   pl_pv1_read                                   n
##   <dbl+lbl>                                 <int>
## 1 1 [Below standard a]                      11986
## 2 2 [Above a standard and below standard b]  8049
## 3 3 [Above standard b]                       2179

ampl_dat |>
  count(pl_pv1_math)

## # A tibble: 3 × 2
##   pl_pv1_math                                   n
##   <dbl+lbl>                                 <int>
## 1 1 [Below standard a]                       9281
## 2 2 [Above a standard and below standard b]  9439
## 3 3 [Above standard b]                       3494

Next, we estimate the proportion of students reaching the minimum proficiency level required at the end of lower primary school (MPLa). The student file already includes binary indicators for reaching MPLa for each plausible value: pl_psa1_math–pl_psa5_math and pl_psa1_read–pl_psa5_read (coded \(1\) = at or above MPLa, \(0\) = below).

ampl_dat |>
  count(pl_psa1_math)

## # A tibble: 2 × 2
##   pl_psa1_math             n
##   <dbl+lbl>            <int>
## 1 0 [Below standard a]  9281
## 2 1 [Above standard a] 12933

Now we estimate the MPLa proportions for reading and mathematics in The Gambia.

result_mpla <- Rrepest(
  data       = subset(ampl_dat, cnt == "GMB"),
  svy        = "SVY",
  est        = est("mean", target = c("pl_psa@_read", "pl_psa@_math")),
  cm.weights = cm,
  var.factor = 135,
  n.pvs      = 5
)
print(result_mpla)

## # A tibble: 1 × 5
##   .     `b.mean.pl_psa@_read` `se.mean.pl_psa@_read` `b.mean.pl_psa@_math` `se.mean.pl_psa@_math`
##   <chr>                 <dbl>                  <dbl>                 <dbl>                  <dbl>
## 1 Total                 0.215                 0.0137                 0.258                 0.0128

We present the {Rrepest} output above in a reader-friendly format using the {gt} package.

as.data.frame(result_mpla) |>
  transmute(
    read_pct = `b.mean.pl_psa@_read` * 100,
    read_se = `se.mean.pl_psa@_read` * 100,
    math_pct = `b.mean.pl_psa@_math` * 100,
    math_se = `se.mean.pl_psa@_math` * 100
  ) |>
  gt() |>
  tab_header(
    title = "Proportion reaching MPLa",
    subtitle = "The Gambia (Grade 3)"
  ) |>
  tab_spanner(label = "Reading", columns = c(read_pct, read_se)) |>
  tab_spanner(label = "Mathematics", columns = c(math_pct, math_se)) |>
  cols_label(
    read_pct = "Percent",
    read_se = "SE",
    math_pct = "Percent",
    math_se = "SE"
  ) |>
  fmt_number(decimals = 1) |>
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_column_labels()
  )

Reading		Mathematics
Proportion reaching MPLa
The Gambia (Grade 3)
Percent	SE	Percent	SE
21.5	1.4	25.8	1.3

In The Gambia, about \(22\%\) of Grade 3 students reach reading MPLa (SE = \(0.014\)) and about \(26\%\) reach mathematics MPLa (SE \(= 0.013\)). In other words, roughly three-quarters of students are still below the end-of-lower-primary minimum proficiency threshold in each domain.

6.7.7 Testing for differences between groups

The table below highlights the distinction between two common types of comparisons in AMPLab data: within-country comparisons (e.g., boys versus girls in the same country) and between-country comparisons (e.g., Lesotho versus Kenya). Because these comparisons involve different survey structures, they require slightly different analytical approaches.

The examples that follow demonstrate how to test for differences in achievement between groups using regression models in Rrepest.

Comparison Type	What it measures	Survey design impact
Within-Country (e.g., Boys vs. Girls in Zambia)	The gap between two demographic subgroups who share the same sampling strata, schools, and teachers.	High covariance. Because groups are clustered together in the same schools, their errors are correlated.
Between-Country (e.g., Lesotho vs Kenya)	The gap between two entirely independent populations with completely separate sampling frames	Zero covariance. Sampling units in Country A have no mathematical relationship to sampling units in Country B

In Rrepest you should use over for within-country comparisons and by for between-country comparisons. You must NOT use over for countries. When using svy = "SVY", the over option is designed for within-country subgroup comparisons only. Using regression achieves the same goal and works correctly with the SVY option.

6.7.7.1 Testing for differences between boys and girls

Are there differences in the proportion of students reaching reading MPLa between girls and boys in Kenya? We can use linear regression to test for this. The girl indicator is already available from earlier in the chapter.

result_kenya <- Rrepest(
  data       = subset(ampl_dat, cnt == "KEN"),
  svy        = "SVY",
  est        = est("lm", target = "pl_psa@_read", regressor = "girl"),
  cm.weights = cm,
  var.factor = 135,
  n.pvs      = 5
)
print(result_kenya)

## # A tibble: 1 × 7
##   .     b.reg_pl_psa@_read.interce…¹ se.reg_pl_psa@_read.…² b.reg_pl_psa@_read.g…³ se.reg_pl_psa@_read.…⁴ b.reg_pl_psa@_read.r…⁵ se.reg_pl_psa@_read.…⁶
##   <chr>                        <dbl>                  <dbl>                  <dbl>                  <dbl>                  <dbl>                  <dbl>
## 1 Total                        0.749                 0.0164                 0.0671                 0.0168                0.00664                0.00327
## # ℹ abbreviated names: ¹`b.reg_pl_psa@_read.intercept`, ²`se.reg_pl_psa@_read.intercept`, ³`b.reg_pl_psa@_read.girl`, ⁴`se.reg_pl_psa@_read.girl`,
## #   ⁵`b.reg_pl_psa@_read.rsqr`, ⁶`se.reg_pl_psa@_read.rsqr`

We present the {Rrepest} output above in a reader-friendly format using the {gt} package.

rrepest_lm_table(
  result_kenya,
  b_cols = c(
    "b.reg_pl_psa@_read.intercept",
    "b.reg_pl_psa@_read.girl",
    "b.reg_pl_psa@_read.rsqr"
  ),
  se_cols = c(
    "se.reg_pl_psa@_read.intercept",
    "se.reg_pl_psa@_read.girl",
    "se.reg_pl_psa@_read.rsqr"
  ),
  terms = c("(Intercept)", "Girl (vs Boy)", "R-squared")
) |>
  show_lm_gt(
    title = "Gender difference in reading MPLa",
    subtitle = "Kenya — weighted linear regression",
    source_note = "Dependent variable: pl_psa*_read (1 = at or above MPLa). Girl = 1, Boy = 0."
  )

Term	Estimate (SE)	t
Gender difference in reading MPLa
Kenya — weighted linear regression
(Intercept)	0.749 (0.016)	45.52
Girl (vs Boy)	0.067 (0.017)	3.98
R-squared	0.007 (0.003)	2.03
Dependent variable: pl_psa*_read (1 = at or above MPLa). Girl = 1, Boy = 0.

In Kenya, girls are about \(6.7\) percentage points more likely than boys to reach reading MPLa (coefficient \(= 0.067\); SE \(= 0.017\)). The intercept implies that about \(75\%\) of boys meet reading MPLa, so the girl estimate corresponds to roughly \(82\%\) among girls. The \(95\%\) confidence interval for the girl coefficient is \(0.034\) to \(0.100\).

Rrepest does not always print p-values in the default table, but you can compute a two-sided p-value from the coefficient and standard error.

# Two-sided p-value from z (large-sample normal approximation)
2 * pnorm(-abs(0.0670775 / 0.0168496))

## [1] 6.863893e-05

A p-value near \(0\) means a gap of this size would be highly unusual if there were no true gender difference in the population. Substantively, Kenya shows a clear within-country advantage for girls on reading MPLa, even though a large share of both girls and boys already meet the threshold.

6.7.7.2 Testing for differences between countries

Next, we ask whether mean mathematics scores differ between students enrolled in Grade 7 in Zambia and Lesotho. Create a dummy variable that equals \(1\) for Grade 7 Lesotho, \(0\) for Grade 7 Zambia, and missing for everyone else.

ampl_dat <- ampl_dat |>
  mutate(
    grade_num = as.integer(haven::zap_labels(grade_stf)),
    LSO_7 = case_when(
      cnt == "LSO" & grade_num == 7 ~ 1L,
      cnt == "ZMB" & grade_num == 7 ~ 0L,
      TRUE ~ NA_integer_
    )
  )

result_lso_zmb <- Rrepest(
  data       = subset(ampl_dat, !is.na(LSO_7)),
  svy        = "SVY",
  est        = est("lm", target = "pv@_math", regressor = "LSO_7"),
  cm.weights = cm,
  var.factor = 135,
  n.pvs      = 5
)
print(result_lso_zmb)

## # A tibble: 1 × 7
##   .     `b.reg_pv@_math.intercept` se.reg_pv@_math.interce…¹ `b.reg_pv@_math.lso_7` se.reg_pv@_math.lso_…² `b.reg_pv@_math.rsqr` `se.reg_pv@_math.rsqr`
##   <chr>                      <dbl>                     <dbl>                  <dbl>                  <dbl>                 <dbl>                  <dbl>
## 1 Total                      -1.05                    0.0374                  0.197                 0.0582               0.00240                0.00147
## # ℹ abbreviated names: ¹`se.reg_pv@_math.intercept`, ²`se.reg_pv@_math.lso_7`

We present the {Rrepest} output above in a reader-friendly format using the {gt} package.

rrepest_lm_table(
  result_lso_zmb,
  b_cols = c(
    "b.reg_pv@_math.intercept",
    "b.reg_pv@_math.lso_7",
    "b.reg_pv@_math.rsqr"
  ),
  se_cols = c(
    "se.reg_pv@_math.intercept",
    "se.reg_pv@_math.lso_7",
    "se.reg_pv@_math.rsqr"
  ),
  terms = c("(Intercept)", "Lesotho Grade 7 (vs Zambia Grade 7)", "R-squared")
) |>
  show_lm_gt(
    title = "Mathematics difference: Lesotho vs Zambia (Grade 7)",
    subtitle = "Weighted linear regression on plausible values",
    source_note = "Dependent variable: mathematics plausible values. LSO_7 = 1 for Lesotho Grade 7."
  )

Term	Estimate (SE)	t
Mathematics difference: Lesotho vs Zambia (Grade 7)
Weighted linear regression on plausible values
(Intercept)	-1.050 (0.037)	−28.07
Lesotho Grade 7 (vs Zambia Grade 7)	0.197 (0.058)	3.39
R-squared	0.002 (0.001)	1.63
Dependent variable: mathematics plausible values. LSO_7 = 1 for Lesotho Grade 7.

Lesotho Grade 7 students score about \(0.20\) LPS points higher in mathematics than Zambia Grade 7 students (SE \(= 0.058\); \(95\%\) CI: \(0.08\) to \(0.31\)). The intercept of \(-1.05\) is the Zambia Grade 7 mean. The gap is statistically clear and matches the earlier country–grade means.

We can also compare a lower part of the distribution using the 25th percentile by country group. (Stata’s qreg estimates a quantile regression coefficient directly; here we report the 25th percentiles for each Grade 7 sample and interpret their difference.)

result_lso_zmb_q25 <- Rrepest(
  data       = subset(ampl_dat, !is.na(LSO_7)),
  svy        = "SVY",
  est        = est(c("quant", 0.25), target = "pv@_math"),
  by         = "LSO_7",
  cm.weights = cm,
  var.factor = 135,
  n.pvs      = 5
)
print(result_lso_zmb_q25)

## # A tibble: 2 × 3
##   LSO_7 `b.quant025.pv@_math` `se.quant025.pv@_math`
##   <chr>                 <dbl>                  <dbl>
## 1 0                     -1.72                 0.0324
## 2 1                     -1.50                 0.0420

We present the {Rrepest} output above in a reader-friendly format using the {gt} package.

as.data.frame(result_lso_zmb_q25) |>
  mutate(
    group = if_else(as.integer(LSO_7) == 1, "Lesotho Grade 7", "Zambia Grade 7"),
    q25 = `b.quant025.pv@_math`,
    se = `se.quant025.pv@_math`
  ) |>
  select(group, q25, se) |>
  gt(rowname_col = "group") |>
  tab_header(
    title = "25th percentile of mathematics proficiency",
    subtitle = "Lesotho vs Zambia (Grade 7)"
  ) |>
  cols_label(q25 = "25th percentile", se = "SE") |>
  fmt_number(decimals = 3) |>
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_column_labels()
  )

	25th percentile	SE
25th percentile of mathematics proficiency
Lesotho vs Zambia (Grade 7)
Zambia Grade 7	−1.723	0.032
Lesotho Grade 7	−1.498	0.042

At the 25th percentile, Zambia Grade 7 scores about \(-1.72\) and Lesotho Grade 7 about \(-1.50\), a gap of roughly \(0.22\) LPS points. The gap is similar in size to the mean difference, which suggests that Lesotho’s higher mathematics performance is not confined to students near the centre of the distribution.

At Grade 7, the relevant minimum proficiency level is MPLb. The student file already includes binary indicators for reaching MPLb for each plausible value: pl_psb1_math–pl_psb5_math and pl_psb1_read–pl_psb5_read (coded \(1\) = at or above MPLb, \(0\) = below).

Are there differences in the percentage reaching the end-of-primary minimum proficiency level (MPLb) for mathematics?

result_mplb <- Rrepest(
  data       = subset(ampl_dat, !is.na(LSO_7)),
  svy        = "SVY",
  est        = est("lm", target = "pl_psb@_math", regressor = "LSO_7"),
  cm.weights = cm,
  var.factor = 135,
  n.pvs      = 5
)
print(result_mplb)

## # A tibble: 1 × 7
##   .     b.reg_pl_psb@_math.interce…¹ se.reg_pl_psb@_math.…² b.reg_pl_psb@_math.l…³ se.reg_pl_psb@_math.…⁴ b.reg_pl_psb@_math.r…⁵ se.reg_pl_psb@_math.…⁶
##   <chr>                        <dbl>                  <dbl>                  <dbl>                  <dbl>                  <dbl>                  <dbl>
## 1 Total                        0.161                 0.0139                 0.0359                 0.0211               0.000608               0.000725
## # ℹ abbreviated names: ¹`b.reg_pl_psb@_math.intercept`, ²`se.reg_pl_psb@_math.intercept`, ³`b.reg_pl_psb@_math.lso_7`, ⁴`se.reg_pl_psb@_math.lso_7`,
## #   ⁵`b.reg_pl_psb@_math.rsqr`, ⁶`se.reg_pl_psb@_math.rsqr`

We present the {Rrepest} output above in a reader-friendly format using the {gt} package.

rrepest_lm_table(
  result_mplb,
  b_cols = c(
    "b.reg_pl_psb@_math.intercept",
    "b.reg_pl_psb@_math.lso_7",
    "b.reg_pl_psb@_math.rsqr"
  ),
  se_cols = c(
    "se.reg_pl_psb@_math.intercept",
    "se.reg_pl_psb@_math.lso_7",
    "se.reg_pl_psb@_math.rsqr"
  ),
  terms = c("(Intercept)", "Lesotho Grade 7 (vs Zambia Grade 7)", "R-squared")
) |>
  show_lm_gt(
    title = "Mathematics MPLb difference: Lesotho vs Zambia (Grade 7)",
    subtitle = "Weighted linear regression",
    source_note = "Dependent variable: pl_psb*_math (1 = at or above MPLb). Intercept = Zambia Grade 7 rate."
  )

Term	Estimate (SE)	t
Mathematics MPLb difference: Lesotho vs Zambia (Grade 7)
Weighted linear regression
(Intercept)	0.161 (0.014)	11.56
Lesotho Grade 7 (vs Zambia Grade 7)	0.036 (0.021)	1.70
R-squared	0.001 (0.001)	0.84
Dependent variable: pl_psb*_math (1 = at or above MPLb). Intercept = Zambia Grade 7 rate.

About \(16\%\) of Zambia Grade 7 students reach mathematics MPLb (intercept \(= 0.161\)). The Lesotho coefficient is \(0.036\) (SE \(= 0.021\)), so Lesotho’s estimated rate is about \(20\%\). The \(95\%\) CI for the difference (\(-0.006\) to \(0.077\)) includes zero, so unlike the mean score gap, the MPLb gap is not statistically clear at conventional levels.

If you prefer an explicit p-value:

# Two-sided p-value from z (large-sample normal approximation)
2 * pnorm(-abs(0.0359181 / 0.0211479))

## [1] 0.08942777

A p-value around \(0.09\) means this MPLb difference would not be unusual if the two Grade 7 populations had the same proficiency rate. The mean score contrast is clearer than the binary MPLb contrast because many students in both countries remain below the MPLb cut score.

6.7.8 Multivariate regression

In this example we look at student and family factors associated with learning in Lesotho. We start with the association between reading scores (pv@_read), gender, and the household wealth index (wealth).

ampl_dat |>
  summarise(
    n = n(),
    min_wealth = min(wealth, na.rm = TRUE),
    max_wealth = max(wealth, na.rm = TRUE),
    n_missing_code = sum(wealth == 999, na.rm = TRUE)
  )

## # A tibble: 1 × 4
##       n min_wealth max_wealth    n_missing_code
##   <int> <dbl+lbl>  <dbl+lbl>              <int>
## 1 22214 -3.79      999 [Missing]            346

Always exclude observations with missing wealth codes (wealth = 999). Including the placeholder value as if it were a real score would distort the regression.

result_multi_read <- Rrepest(
  data       = subset(ampl_dat, cnt == "LSO" & wealth < 999),
  svy        = "SVY",
  est        = est("lm", target = "pv@_read", regressor = c("girl", "wealth")),
  cm.weights = cm,
  var.factor = 135,
  n.pvs      = 5
)
print(result_multi_read)

## # A tibble: 1 × 9
##   .     `b.reg_pv@_read.intercept` se.reg_pv@_read.interce…¹ `b.reg_pv@_read.girl` `se.reg_pv@_read.girl` b.reg_pv@_read.wealt…² se.reg_pv@_read.weal…³
##   <chr>                      <dbl>                     <dbl>                 <dbl>                  <dbl>                  <dbl>                  <dbl>
## 1 Total                     -0.447                    0.0517                 0.260                 0.0397                  0.159                 0.0287
## # ℹ abbreviated names: ¹`se.reg_pv@_read.intercept`, ²`b.reg_pv@_read.wealth`, ³`se.reg_pv@_read.wealth`
## # ℹ 2 more variables: `b.reg_pv@_read.rsqr` <dbl>, `se.reg_pv@_read.rsqr` <dbl>

We present the {Rrepest} output above in a reader-friendly format using the {gt} package.

rrepest_lm_table(
  result_multi_read,
  b_cols = c(
    "b.reg_pv@_read.intercept",
    "b.reg_pv@_read.girl",
    "b.reg_pv@_read.wealth",
    "b.reg_pv@_read.rsqr"
  ),
  se_cols = c(
    "se.reg_pv@_read.intercept",
    "se.reg_pv@_read.girl",
    "se.reg_pv@_read.wealth",
    "se.reg_pv@_read.rsqr"
  ),
  terms = c("(Intercept)", "Girl (vs Boy)", "Household wealth", "R-squared")
) |>
  show_lm_gt(
    title = "Predictors of reading proficiency",
    subtitle = "Lesotho — weighted linear regression",
    source_note = "Dependent variable: reading plausible values. Restricted to wealth < 999."
  )

Term	Estimate (SE)	t
Predictors of reading proficiency
Lesotho — weighted linear regression
(Intercept)	-0.447 (0.052)	−8.63
Girl (vs Boy)	0.260 (0.040)	6.55
Household wealth	0.159 (0.029)	5.53
R-squared	0.073 (0.019)	3.76
Dependent variable: reading plausible values. Restricted to wealth < 999.

In Lesotho, girls score about \(0.26\) LPS points higher in reading than boys after controlling for wealth (SE \(= 0.040\); \(95\%\) CI: \(0.18\) to \(0.34\)). A one-unit increase in the wealth index is associated with about \(0.16\) points higher reading (SE \(= 0.029\); \(95\%\) CI: \(0.10\) to \(0.22\)). Both associations are statistically clear. The intercept (\(-0.45\)) is the expected reading score for boys with wealth equal to zero on the index.

Next, we add parents’ highest education level and switch the outcome to mathematics MPLb, so the coefficients are on the probability scale. Because Rrepest treats named regressors as continuous unless you create separate indicator variables, we build ISCED dummies with “did not complete ISCED level 1” (hisced == 0) as the reference group.

ampl_dat |>
  count(hisced)

## # A tibble: 7 × 2
##   hisced                                 n
##   <dbl+lbl>                          <int>
## 1 0 [Did not complete ISCED level 1]  1480
## 2 1 [ISCED level 1]                   3434
## 3 2 [ISCED level 2]                   3461
## 4 3 [ISCED level 3]                   5400
## 5 4 [ISCED level 4 or 5]              3328
## 6 5 [ISCED level 6 or higher]         4304
## 7 9 [Missing/Unknown]                  807

ampl_dat_lso <- ampl_dat |>
  filter(cnt == "LSO", wealth < 999, hisced < 9) |>
  mutate(
    hisced1 = as.integer(hisced == 1),
    hisced2 = as.integer(hisced == 2),
    hisced3 = as.integer(hisced == 3),
    hisced4 = as.integer(hisced == 4),
    hisced5 = as.integer(hisced == 5)
  )

result_multi_mplb <- Rrepest(
  data       = ampl_dat_lso,
  svy        = "SVY",
  est        = est(
    "lm",
    target = "pl_psb@_math",
    regressor = c("girl", "wealth", "hisced1", "hisced2", "hisced3", "hisced4", "hisced5")
  ),
  cm.weights = cm,
  var.factor = 135,
  n.pvs      = 5
)
print(result_multi_mplb)

## # A tibble: 1 × 19
##   .     b.reg_pl_psb@_math.interce…¹ se.reg_pl_psb@_math.…² b.reg_pl_psb@_math.g…³ se.reg_pl_psb@_math.…⁴ b.reg_pl_psb@_math.w…⁵ se.reg_pl_psb@_math.…⁶
##   <chr>                        <dbl>                  <dbl>                  <dbl>                  <dbl>                  <dbl>                  <dbl>
## 1 Total                       0.0914                 0.0320                 0.0624                 0.0148                 0.0313                0.00840
## # ℹ abbreviated names: ¹`b.reg_pl_psb@_math.intercept`, ²`se.reg_pl_psb@_math.intercept`, ³`b.reg_pl_psb@_math.girl`, ⁴`se.reg_pl_psb@_math.girl`,
## #   ⁵`b.reg_pl_psb@_math.wealth`, ⁶`se.reg_pl_psb@_math.wealth`
## # ℹ 12 more variables: `b.reg_pl_psb@_math.hisced1` <dbl>, `se.reg_pl_psb@_math.hisced1` <dbl>, `b.reg_pl_psb@_math.hisced2` <dbl>,
## #   `se.reg_pl_psb@_math.hisced2` <dbl>, `b.reg_pl_psb@_math.hisced3` <dbl>, `se.reg_pl_psb@_math.hisced3` <dbl>, `b.reg_pl_psb@_math.hisced4` <dbl>,
## #   `se.reg_pl_psb@_math.hisced4` <dbl>, `b.reg_pl_psb@_math.hisced5` <dbl>, `se.reg_pl_psb@_math.hisced5` <dbl>, `b.reg_pl_psb@_math.rsqr` <dbl>,
## #   `se.reg_pl_psb@_math.rsqr` <dbl>

We present the {Rrepest} output above in a reader-friendly format using the {gt} package.

rrepest_lm_table(
  result_multi_mplb,
  b_cols = c(
    "b.reg_pl_psb@_math.intercept",
    "b.reg_pl_psb@_math.girl",
    "b.reg_pl_psb@_math.wealth",
    "b.reg_pl_psb@_math.hisced1",
    "b.reg_pl_psb@_math.hisced2",
    "b.reg_pl_psb@_math.hisced3",
    "b.reg_pl_psb@_math.hisced4",
    "b.reg_pl_psb@_math.hisced5",
    "b.reg_pl_psb@_math.rsqr"
  ),
  se_cols = c(
    "se.reg_pl_psb@_math.intercept",
    "se.reg_pl_psb@_math.girl",
    "se.reg_pl_psb@_math.wealth",
    "se.reg_pl_psb@_math.hisced1",
    "se.reg_pl_psb@_math.hisced2",
    "se.reg_pl_psb@_math.hisced3",
    "se.reg_pl_psb@_math.hisced4",
    "se.reg_pl_psb@_math.hisced5",
    "se.reg_pl_psb@_math.rsqr"
  ),
  terms = c(
    "(Intercept)",
    "Girl (vs Boy)",
    "Household wealth",
    "ISCED 1 (vs below ISCED 1)",
    "ISCED 2 (vs below ISCED 1)",
    "ISCED 3 (vs below ISCED 1)",
    "ISCED 4 (vs below ISCED 1)",
    "ISCED 5 (vs below ISCED 1)",
    "R-squared"
  )
) |>
  show_lm_gt(
    title = "Predictors of mathematics MPLb",
    subtitle = "Lesotho - weighted linear regression",
    source_note = "Dependent variable: pl_psb*_math (1 = at or above MPLb). Restricted to wealth < 999 and hisced < 9."
  )

Term	Estimate (SE)	t
Predictors of mathematics MPLb
Lesotho - weighted linear regression
(Intercept)	0.091 (0.032)	2.86
Girl (vs Boy)	0.062 (0.015)	4.21
Household wealth	0.031 (0.008)	3.72
ISCED 1 (vs below ISCED 1)	0.012 (0.032)	0.37
ISCED 2 (vs below ISCED 1)	0.062 (0.032)	1.91
ISCED 3 (vs below ISCED 1)	0.114 (0.041)	2.77
ISCED 4 (vs below ISCED 1)	0.226 (0.043)	5.31
ISCED 5 (vs below ISCED 1)	0.185 (0.050)	3.70
R-squared	0.073 (0.018)	3.98
Dependent variable: pl_psb*_math (1 = at or above MPLb). Restricted to wealth < 999 and hisced < 9.

After controlling for wealth and parental education, girls remain about \(6.2\) percentage points more likely than boys to reach mathematics MPLb (SE \(= 0.015\)). Wealth retains a positive association (\(0.031\) per unit; SE \(= 0.008\)). Relative to parents who did not complete ISCED level 1, higher parental education is associated with higher MPLb rates, especially at ISCED levels 3–5 (coefficients from about \(0.11\) to \(0.23\)). The lowest education contrast (ISCED 1) is not statistically distinguishable from the reference group. Restricting to hisced < 9 drops missing/unknown parental education codes, which is important for valid inference.

6 Analyzing AMPLab data in R

6.1 Package installation

6.2 Loading data in R

6.3 Add variable and value labels

6.4 The Rrepest package

6.5 Rrepest command syntax

6.6 Before you begin: Set up Rrepest for AMPLab

6.7 AMPLab Analyses Examples

6.7.1 Calculating Mean Age by Country

Alternative to generate the same result using survey

6.7.2 Calculating mean mathematics proficiency

6.7.3 Proficiency across gender

6.7.4 Proficiency by Country

6.7.5 Moving beyond means: Summary statistics with Rrepest

6.7.6 Estimating AMPL Proficiency Levels

6.7.7 Testing for differences between groups

6.7.7.1 Testing for differences between boys and girls

6.7.7.2 Testing for differences between countries

6.7.8 Multivariate regression

6.4 The `Rrepest` package

6.5 `Rrepest` command syntax

6.6 Before you begin: Set up `Rrepest` for AMPLab

6.7.5 Moving beyond means: Summary statistics with `Rrepest`