Estimation of Standard Errors using Replicate Weights

Because TIMSS and other ILSAs sample schools first (the primary sampling unit or PSU) and then classrooms/students within schools using multistage sampling, standard variance formulas that assume independent observations would understate the true sampling error.

Replicate weights are a set of variables used to correctly compute the standard errors of a point estimate (e.g., regression coefficient, mean). Jackknife Repeated Replication (JRR) is a method for generating replicate weights. Conceptually, a statistic of interest is computed as many times as there are replicate weights and the variability of that statistic is used in the computation of the standard error. JRR involves (specifically for a variant referred to as JK2 used in TIMSS):

Creating the weights:

Using the weights:

TIMSS typically uses 75–150+ replicate weights (depending on the number of primary sampling units), and the information regarding the replicate weights are provided directly in the international database— you just apply the JRR weights using software (e.g.,the IEA IDB Analyzer, R) that support this replication method. NOTE: the replication weights are not in the dataset but the JKZONE and JKREP variables are.

Using regular standard errors on TIMSS data without JRR will more likely be statistically significant since it ignores the clustering of students within schools (i.e., standard errors are underestimated).

In statistics, Maurice Quenouille introduced jackknifing (in 1949) and a jackknife (coined by Tukey referring to a multipurpose tool) is referred to as a resampling technique for estimating the variance of a statistic (sometimes referred to as the “delete-1 jackknife”). The standard jackknife works using the following general steps:

  1. Compute the statistic of interest (e.g., mean, regression coefficient; \(\theta_{Full}\)– we don’t use this here though but will use it with TIMSS or PISA.
  2. Recompute the same statistic using n samples but with each sample, leaving out one i observation (\(\theta_{-i}\))
  3. Compute the variability of the statistic (in this case the variability of the mean; used to get the standard error of the mean– not standard deviation of the mean):

\[\frac{(n-1)}{n} \Sigma({\hat{\theta}_{-i}}-\bar{\theta}_.)^2\]

data(mtcars)
n <- nrow(mtcars) #32 observations
full_mean <- mean(mtcars$mpg) #theta
mns <- numeric() #container to store estimated theta
for (i in 1:n){
  mns[i] <- mean(mtcars$mpg[-i]) #leave one out
}
jack_var <- ((n - 1)/n) * sum((mns - mean(mns))^2) #(n - 1) / n is a jackknife correction factor
sqrt(jack_var)
[1] 1.065424
sd(mtcars$mpg) / sqrt(n)
[1] 1.065424

Jackknifing is different from bootstrapping in the sense that it is deterministic– it will always give the same result compared to a bootstrap which is based on random draws. There is also a difference is the original overall mean is used vs. the mean of the means.

An example using a large scale assessment (TIMSS case)

This data file is similar to the LSAs we use (i.e., nested, uses disproportional sampling of different strata, oversampling:)

library(here)
library(dplyr)
library(survey)
library(sjmisc)
# df1 <- read.table(here("05_data/dat1_150.dat"))
# names(df1) <- c('id',  'schwgt', 'y',  'w1',  'x1',
#  'minor',  'one',  'nwt', 'totwgt',  'stuwgt',
#  'strata', 'mos')
df1 <- rio::import('https://github.com/flh3/pubdata/raw/refs/heads/main/ILSA/n150.sav')
n_distinct(df1$id) #school id
[1] 150
n_distinct(df1$strata)
[1] 3

This is just a synthetic dataset of 150 schools, with 25 students per school (25 * 150):

nrow(df1)
[1] 3750
head(df1)
  id schwgt     y     w1     x1 minor one   nwt totwgt stuwgt strata mos
1  8 13.673 1.595 -0.176 -2.011     0   1 0.561 39.653    2.9      2  65
2  8 13.673 0.788 -0.176 -1.525     0   1 0.561 39.653    2.9      2  65
3  8 13.673 0.371 -0.176 -1.735     0   1 0.561 39.653    2.9      2  65
4  8 13.673 2.219 -0.176 -0.627     0   1 0.561 39.653    2.9      2  65
5  8 13.673 1.066 -0.176  1.735     0   1 0.561 39.653    2.9      2  65
6  8 13.673 0.653 -0.176 -0.340     0   1 0.561 39.653    2.9      2  65

The following syntax shows how the zones and rep are created.

First, sort by strata and measure of size (enrollment). Second, a pair of schools is placed in a zone (i.e., 2 schools per zone).

# 1. Build a SCHOOL-level frame, sorted within stratum for pairing
school_frame <- df1 %>%
  distinct(id, strata, mos) %>%
  arrange(strata, mos) %>%
  group_by(strata) %>%
  mutate(
    within_stratum_rank = row_number(),
    n_in_stratum = n(),
    pair_id = ceiling(within_stratum_rank / 2)
  ) %>%
  ungroup() %>%
  mutate(zone = paste(strata, pair_id, sep = "_"))

In cases where are there are no pairs (e.g., there is an extra school in strata 1 since there are 22 pairs + 1 extra since there are 45 schools), a pseudo-pair is created by splitting the last school randomly into two and giving it a pseudo zone (the one school is split in two). See zones 1_23 and 3_38.

Showing the last 6 observations in the 3 strata:

school_frame %>% group_by(strata) %>% slice_tail(n = 6)
# A tibble: 18 × 7
# Groups:   strata [3]
      id strata   mos within_stratum_rank n_in_stratum pair_id zone 
   <dbl>  <dbl> <dbl>               <int>        <int>   <dbl> <chr>
 1  1922      1   137                  40           45      20 1_20 
 2  1286      1   141                  41           45      21 1_21 
 3  1534      1   146                  42           45      21 1_21 
 4  1406      1   158                  43           45      22 1_22 
 5  2256      1   168                  44           45      22 1_22 
 6  1741      1   187                  45           45      23 1_23 
 7   213      2   148                  25           30      13 2_13 
 8   167      2   159                  26           30      13 2_13 
 9   272      2   169                  27           30      14 2_14 
10   123      2   186                  28           30      14 2_14 
11    20      2   214                  29           30      15 2_15 
12   178      2   275                  30           30      15 2_15 
13  1995      3   172                  70           75      35 3_35 
14  1303      3   179                  71           75      36 3_36 
15  2191      3   192                  72           75      36 3_36 
16  2032      3   201                  73           75      37 3_37 
17  2092      3   215                  74           75      37 3_37 
18  1517      3   263                  75           75      38 3_38

Have to identify the unpaired groups– those groups/schools must be split into two “quasi-schools” to complete its zone:

odd_leftover_schools <- school_frame %>%
  group_by(strata) %>%
  filter(n_in_stratum %% 2 == 1, within_stratum_rank == n_in_stratum) %>%
  pull(id)

After, split each leftover school’s students randomly in half; each half becomes a quasi-school for JK2 purposes:

df1$quasi_school <- df1$id

for (sch in odd_leftover_schools) {
  ids  <- which(df1$id == sch)
  half <- sample(ids, size = floor(length(ids) / 2)) #randomly choose
  df1$quasi_school[half] <- paste0(sch, "_a")
  df1$quasi_school[setdiff(ids, half)]   <- paste0(sch, "_b")
}

# Rebuild school_frame at the quasi-school level so split schools count as two separate pairing units
school_frame <- df1 %>%
  distinct(quasi_school, id, strata, mos) %>%
  arrange(strata, mos, quasi_school) %>%
  group_by(strata) %>%
  mutate(pair_id = ceiling(row_number() / 2)) %>%
  ungroup() %>%
  mutate(zone = paste(strata, pair_id, sep = "_")) %>%
  group_by(zone) %>%
  mutate(rep = row_number() - 1) %>%  # 0 or 1 -- which half of the pair
  ungroup()

# Attach zone + member back onto the student-level data
df1 <- df1 %>%
  left_join(
    school_frame %>% select(quasi_school, zone, rep),
    by = "quasi_school"
  )

n_zones <- length(unique(df1$zone))
cat("Number of jackknife zones:", n_zones)  
Number of jackknife zones: 76

There are 76 zones. Within each zone, the schools are split (1 or 0; the jackknife replication code)

library(fhutils)
ctab(df1, zone, rep) %>% head()
 zone/rep          0          1       Total
      1_1 24 (49.0%) 25 (51.0%) 49 (100.0%)
     1_10 25 (50.0%) 25 (50.0%) 50 (100.0%)
     1_11 25 (50.0%) 25 (50.0%) 50 (100.0%)
     1_12 25 (50.0%) 25 (50.0%) 50 (100.0%)
     1_13 25 (50.0%) 25 (50.0%) 50 (100.0%)
     1_14 25 (50.0%) 25 (50.0%) 50 (100.0%)

This is what is seen in TIMSS with the variables: JKZONE (e.g., 75 zones) and JKREP (i.e., 0 vs 1).

Sample TIMSS dataset:

library(BIFIEsurvey)
data(data.timss3)
head(data.timss3)
     IDSTUD   TOTWGT JKZONE JKREP female books lang migrant scsci likesc
1 400010201 17.46818      1     1      1     3    1       0    NA      2
2 400010203 17.46818      1     1      0     3    1       0     2      4
3 400010204 17.46818      1     1      1     5    1       0     2      2
4 400010205 17.46818      1     1      1     3    1       0     1      1
5 400010206 17.46818      1     1      1     3    1       0     2      1
6 400010207 17.46818      1     1      1     2    1       0     3      1
  ASMMAT1 ASSSCI1 ASMMAT2 ASSSCI2 ASMMAT3 ASSSCI3 ASMMAT4 ASSSCI4 ASMMAT5
1 542.768 600.128 557.154 578.214 505.855 570.009 523.688 560.259 577.832
2 521.717 512.277 532.734 519.171 556.618 554.299 510.659 505.586 545.532
3 455.962 496.747 462.111 544.805 445.180 528.147 472.595 549.990 456.757
4 512.242 583.776 510.212 613.559 530.564 569.237 497.112 597.199 527.811
5 505.640 532.747 563.068 567.919 529.584 563.538 488.248 482.694 583.264
6 454.723 483.655 405.011 459.855 408.625 451.032 434.679 464.432 343.622
  ASSSCI5
1 607.334
2 565.149
3 545.888
4 623.008
5 578.385
6 317.674
n_distinct(data.timss3$JKZONE)
[1] 75
n_distinct(data.timss3$JKREP)
[1] 2

Now to use this properly:

We can use the survey package. There are two ways

  1. Make a matrix of replicate weights.

There are 75 columns (or as many zones) of the original weights. The only difference, is that in each zone:

  • One of the schools gets a weight of zero (i.e., zeroed out)
  • The second school in the pair gets double the weight (\(\times\) 2).
  • All other weights are the same (\(\times\) 1)

I have written a function (with help from AI!) that will make the replicate weights. You need to specify the zone, the rep/member, and the overall weight variable.

library(survey)
build_jk2_rw <- function(zone, rep, weight) {
  zn <- unique(zone)
  n_zones <- length(zn)
  n <- length(weight)
  m <- matrix(weight, nrow = n, ncol = n_zones)
  colnames(m) <- paste0("rw", seq_len(n_zones))
  for (h in seq_len(n_zones)) {
    in_zone <- zone == zn[h]
    m[in_zone & rep == 0, h] <- weight[in_zone & rep == 0] * 2
    m[in_zone & rep == 1, h] <- 0
  }
  m
}
rw <- build_jk2_rw(df1$zone, df1$rep, df1$totwgt)
dim(rw)
[1] 3750   76

Once we have the matrix of replicate weights (rw), we can use this in the survey package:

des <- svrepdesign(
  data = df1,
  repweights = rw, #specified here
  weights = ~totwgt, #note the ~,
  mse = TRUE,
  type = 'JK2'
)
summary(des)
Call: svrepdesign.default(data = df1, repweights = rw, weights = ~totwgt, 
    mse = TRUE, type = "JK2")
JK2 jackknife with 76 replicates and MSE variances.
Variables: 
 [1] "id"           "schwgt"       "y"            "w1"           "x1"          
 [6] "minor"        "one"          "nwt"          "totwgt"       "stuwgt"      
[11] "strata"       "mos"          "quasi_school" "zone"         "rep"

After specifying the design, can use the svy functions and instead of the dataframe, pass along the design of the survey.

svymean(~y + w1 + x1 + minor, des)
           mean     SE
y      0.424941 0.0672
w1    -0.022955 0.0852
x1     0.013857 0.0184
minor  0.190580 0.0087

This shows both the mean and the standard error.

NOTE: the standard error is computed as:

\[SE = \sqrt{\Sigma{(\hat{\theta}_j - \bar{\theta})^2}}\] where \(\hat{\theta_j}\) is the statistic of interest computed using \(j\) weights from the \(j\) column. \(\bar{\theta}\) is the overall weighted mean (using \(totwgt\)) (done when mse = TRUE as is recommended in the TIMSS manual). If mse = FALSE then the mean of the \(\hat{\theta_j}\) is used (results can be slightly different).

To compute this manually (for w1):

js <- ncol(rw) #how many replicate weights
wbar <- weighted.mean(df1$w, df1$totwgt) #overall weight
mn <- numeric() #container for j means
for (i in 1:js){
  mn[i] <- weighted.mean(df1$w, rw[,i])
}
hist(mn) #76 means

sqrt(sum((mn - wbar)^2)) #standard error formula
[1] 0.08518286

The standard error is the similar to the one that was done using the function.

  1. The second way to do this is to also use another function in the survey package (acceptable but could be slightly different). As the dataset specifies the ZONE and the REP variables, we can use these:
des2 <- svydesign(
  data = df1,
  weights = ~totwgt,
  strata = ~zone, #JKZONE in TIMSS
  id =~ rep, #JKREP in TIMSS
  nest = TRUE #this is done as rep is 1 or 0; not unique
)

Once that is done, need to specify that as a jackknife replication:

 des2a <- as.svrepdesign(
    des2, 
    type = "JKn"
  )

Note: this is specified as “JKn” (not JK2). Then can use the svy functions as usual:

svymean(~y + w1 + x1 + minor, des2a)
           mean     SE
y      0.424941 0.0673
w1    -0.022955 0.0852
x1     0.013857 0.0185
minor  0.190580 0.0087

The results are similar (the same?) as in the previous example.

NOTE: that the naive standard error (i.e., \(\frac{SD}{\sqrt{n}}\)is too small (especially for w1):

psych::describe(select(df1, y, w1, x1, minor), skew = F, ranges = F)
      vars    n  mean   sd   se
y        1 3750  0.40 1.52 0.02
w1       2 3750 -0.04 0.94 0.02
x1       3 3750  0.01 1.01 0.02
minor    4 3750  0.27 0.44 0.01

Note that this is using unweighted data so the means will be different.

  1. Yet another way to specify this– but this does not really use replicate weights but is similar to using cluster robust standard errors:
des3 <- svydesign(
  data = df1,
  id = ~id, #school id
  weights = ~totwgt
)

There are slight differences– this is still better than the naive standard errors.

svymean(~y + w1 + x1 + minor, des3)
           mean     SE
y      0.424941 0.0658
w1    -0.022955 0.0795
x1     0.013857 0.0196
minor  0.190580 0.0102

There are many svy functions in the survey package. If we want fit a linear regression, can use svyglm:

m1 <- svyglm(y ~ w1 + x1 + minor, des)
m2 <- svyglm(y ~ w1 + x1 + minor, des2a)
m3 <- svyglm(y ~ w1 + x1 + minor, des3)
m4 <- estimatr::lm_robust(y ~ w1 + x1 + minor, data = df1, weight = totwgt) #this is naive OLS; incorrect
library(modelsummary)
modelsummary(list("RW" = m1, "RW2" = m2, "CRSE" = m3, "OLS" = m4), stars = TRUE)
RW RW2 CRSE OLS
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
(Intercept) 0.524*** 0.524*** 0.524*** 0.524***
(0.054) (0.054) (0.056) (0.028)
w1 0.533*** 0.533*** 0.533*** 0.533***
(0.061) (0.062) (0.063) (0.026)
x1 0.489*** 0.489*** 0.489*** 0.489***
(0.022) (0.022) (0.022) (0.024)
minor -0.490*** -0.490*** -0.490*** -0.490***
(0.057) (0.057) (0.053) (0.051)
Num.Obs. 3750 3750 3750 3750
R2 0.234 0.234 0.234 0.234
R2 Adj. -38.887 -38.341 -18.670 0.233
AIC 12671.9 12672.4 12672.4 13206.4
BIC 110409.6 108895.8 71155.6 13237.5
Log.Lik. -55184.202 -54427.325 -35557.225
F 232.274 228.540 215.093
RMSE 1.31 1.31 1.31 1.31

Here is a list of svy functions:

Category Functions
Descriptive Statistics svymean(), svytotal(), svyvar(), svyquantile(), svyby(), svyratio(), svycontrast()
Tabulations & Tests svytable(), svychisq(), svyttest(), svyranktest(), svyhist(), svyboxplot()
Generalized Linear Models svyglm(), svyolr(), svycoxph(), svymle(), svykm(), svysmooth()

Practice this using the TIMSS dataset.

Replication using PISA Replication Weights

The PISA replication weights are constructed differently from the TIMSS weights. This is referred to as using Balanced Repeated Replications (BRR).

data("data.pisaNLD") #from BIFIEsurvey
df2 <- data.pisaNLD
dim(df2)
[1] 3992  405
#has a lot of variables; has 80 replicate weights

The replicate weights in PISA are given and can just be used. The weights are computed differently from what TIMSS does (as shown). The example below specifies the design using these weights. The names of all the replicate weights have to be entered (using the function below using dplyr but can also use the alternative which is using regular expressions).

pisa_des <- svrepdesign(
  data = df2,
  weights = ~W_FSTUWT, #don't forget the ~
  repweights = dplyr::select(df2, starts_with("W_FSTR")), #"W_FSTURWT[0-9]+",
  type = "Fay", # PISA's method: Fay's modified BRR
  rho = 0.5, # Fay's perturbation factor k
  mse = TRUE # PISA default, use original mean
)

Of note: there is 1) type (dictated by PISA) and 2) rho. There is an option combined.weights and the default is fine (TRUE)– just means the the replicate weights are already sampling weights (not just a perturbation factor).

PISA doesn’t use the JK2 jackknife— it uses Fay’s Balanced Repeated Replication (BRR), where instead of dropping one PSU per replicate (i.e., the original weight x 0 as done in TIMSS), each replicate perturbs roughly half the sample up by a factor of 2−k and half down by a factor of k, with k = .5 (that is fixed or 1.5 vs 0.5).

The halfs are determined by a Hadamard matrix (no need to get into that). In this case, no group ever gets a zero weight. This is also what rho is– it is the perturbation factor k (and given in PISA documentation). We don’t need to construct the weights- they are already given (see the dataset; there are so many of them).

To get the mean of ESCS:

svymean(~ESCS, pisa_des, na.rm = T) #has missing values
         mean     SE
ESCS 0.097788 0.0234

To compute this manually (this helps my understanding of this; you won’t necessarily do this):

rp2 <- dplyr::select(df2, starts_with("W_FSTR")) #matrix of replicate weights
mns <- numeric()
for (x in 1:80){
  mns[x] <- weighted.mean(df2$ESCS, rp2[,x], na.rm = T)
}
omean <- weighted.mean(df2$ESCS, df2$W_FSTUWT, na.rm = TRUE) #the original mean
sqrt((1/20) * sum((mns - omean)^2)) #standard error
[1] 0.02335901

The adjustment factor (1/20) is based on R (i.e., the number of replicates or 80 for PISA) and the correction factor k (rho or the Fay factor) and is \(\frac{1}{R(1-k)^2}\) or \(\frac{1}{80(1-.5)^2} = \frac{1}{20}\).

We can use the different svy functions since the design is already given.

What about a correlation matrix?

We can use the jtools package for this (survey can do this but without p-values):

library(jtools)
svycor(~PV1MATH + PV2MATH + ESCS, design = pisa_des, 
       na.rm = TRUE, sig.stats = TRUE)
        PV1MATH PV2MATH ESCS 
PV1MATH 1       0.93*   0.43*
PV2MATH 0.93*   1       0.43*
ESCS    0.43*   0.43*   1

NOTE: the na.rm = TRUE is set since there are missing values. If not, you will get an error. Using:

svyvar(~PV1MATH + PV2MATH + ESCS, design = pisa_des,
       na.rm = T) |>
  as.matrix() |>
  cov2cor()

What if there are plausible values (PVs)?

PVs are a form of multiply imputed values for each student. See notes. The estimate (b) is just the average of the estimate (that’s easy). The standard errors are combined using this formula based on Rubin’s rules (\(m\) is the number of plausible values; \(T_B\) is the variance of the estimate):

\[T_B = \overline{U}_b + \left(1 + \frac{1}{m}\right) Bb\]

  • \(\overline{U}_b\) is simply: \(\overline{U}_b = \frac{\Sigma SE_b^2}{m}\). It is the average of the squared standard errors (i.e., the variance)

  • Bb is \(B_b = \frac{\Sigma(b - \bar{b})^2}{m - 1}\) or just the variance of the estimate

If we were to do this manually with 5 PVs:

tmp <- svymean(~PV1MATH + PV2MATH + PV3MATH + PV4MATH + PV5MATH, pisa_des)
coef(tmp) #these are the means
 PV1MATH  PV2MATH  PV3MATH  PV4MATH  PV5MATH 
538.0611 537.7628 537.8091 537.2219 538.2615 
SE(tmp) #these are the standard errors
[1] 3.190609 3.085825 3.074595 3.111499 3.038390

We can just take the average of the estimates to get the estimate. To get \(Bb\), we just take the variance.

est <- mean(coef(tmp)) #the mean of all the imputations
est #our statistics of interest
[1] 537.8233
Bb <- var(coef(tmp)) #the variance of the estimates
Ubar <- mean(SE(tmp)^2) #this is the SE-- don't forget to square it!

Putting this together:

T_b <- Ubar + (1 + 1/5) * Bb

The standard error is the square root of this value or \(\sqrt{T_b}\):

sqrt(T_b) #this is the correct pooled SE
[1] 3.130174

We won’t be doing this manualy. This is easy to specify once you have the design. The syntax just looks clunky!

library(mitools) #for combining the values properly

# Define the mapping of the analysis variable name to the dataset columns
pv_map <- list(math ~ PV1MATH + PV2MATH + PV3MATH + PV4MATH + PV5MATH)

# Run svy function across all plausible values
# VERY IMPORTANT::: only use one variable --> DO NOT DO ~math + ESCS

pv_models <- withPV(
  mapping = pv_map,
  data = pisa_des, #indicate the design
  action = function(x) svymean(~math, design = x) #x is the design
)

# Pool the estimates and standard errors using Rubin's rules
pooled_res <- MIcombine(pv_models)

# Inspect pooled results
summary(pooled_res)
Multiple imputation results:
      withPV.svyrep.design(mapping = pv_map, data = pisa_des, action = function(x) svymean(~math, 
    design = x))
      MIcombine.default(pv_models)
         results       se   (lower  upper) missInfo
PV1MATH 537.8233 3.130174 531.6876 543.959      2 %

Can do this for a regression model (the syntax just looks strange but what’s important is to change the regression model):

# pv_reg <- withPV(
#   mapping = pv_map,
#   data = pisa_des, #indicate the design
#   rewrite = F,
#   action = quote(svyglm(math ~ ESCS + factor(LANG) +
#           factor(ST03Q01), design = .DESIGN))
# )

# I like this way better
pv_reg <- withPV(
  mapping = pv_map,
  data = pisa_des, #indicate the design
  action = \(x) svyglm(math ~ ESCS + factor(LANG) +
          factor(ST03Q01), design = x)
)

res2 <- MIcombine(pv_reg)
summary(res2)
Multiple imputation results:
      withPV.svyrep.design(mapping = pv_map, data = pisa_des, action = function(x) svyglm(math ~ 
    ESCS + factor(LANG) + factor(ST03Q01), design = x))
      MIcombine.default(pv_reg)
                    results       se     (lower    upper) missInfo
(Intercept)      540.737455 3.134911 534.591710 546.88320      3 %
ESCS              42.589146 2.291321  38.094239  47.08405      6 %
factor(LANG)1    -47.885371 8.717825 -65.004928 -30.76581      8 %
factor(ST03Q01)2   4.815302 3.412780  -1.874784  11.50539      2 %

Manually running a regression and then combining (note this is just using cluster robust variance estimation, not using the replicate weights).

Need some data management. Create a list of imputed datasets, one dataset per plausible value.

library(tidyr)
sm <- select(df2, PV1MATH:PV5MATH, ESCS, LANG, ST03Q01, SCHOOLID, W_FSTUWT)
tall <- pivot_longer(sm, PV1MATH:PV5MATH, values_to = "MATH")
tall$pv <- readr::parse_number(tall$name)
datlist <- split(tall, tall$pv)
length(datlist)
[1] 5

datlist is now a list of the five datasets, one for each plausible value.

library(estimatr)
allmods <- lapply(datlist, \(x)
    lm_robust(MATH ~ ESCS + factor(LANG) + factor(ST03Q01), weights = W_FSTUWT, cluster = SCHOOLID, data = x, se_type = 'stata'))

pooled <- MIcombine(allmods)
summary(pooled)
Multiple imputation results:
      MIcombine.default(allmods)
                    results       se     (lower    upper) missInfo
(Intercept)      540.737455 4.779036 531.370312 550.10460      1 %
ESCS              42.589146 2.748725  37.199428  47.97886      4 %
factor(LANG)1    -47.885371 8.428535 -64.441448 -31.32929      9 %
factor(ST03Q01)2   4.815302 3.510376  -2.065975  11.69658      2 %

Results are similar. Can also use:

library(mitml)
testEstimates(allmods)

Call:

testEstimates(model = allmods)

Final parameter estimates and inferences obtained from 5 imputed data sets.

                  Estimate Std.Error   t.value        df   P(>|t|)       RIV       FMI 
(Intercept)        540.737     4.779   113.148 28050.506     0.000     0.012     0.012 
ESCS                42.589     2.749    15.494  2815.845     0.000     0.039     0.038 
factor(LANG)1      -47.885     8.429    -5.681   549.699     0.000     0.093     0.089 
factor(ST03Q01)2     4.815     3.510     1.372  7816.348     0.170     0.023     0.023 

Unadjusted hypothesis test as appropriate in larger samples.

— END (2026.09.03)