Example Using Replicate Weights with PISA

As the PISA 2025 results were just released this week and I am currently teaching a seminar/class on analyzing international large scale assessments this fall 2026, I thought it would good to show how we actually use these replicate weights when analyzing these datasets. Replicate weights are used to get the correct standard errors when analyzing ILSAs. I don’t think many applied researchers understand what these are– even though a lot has been written about them– and are explained in the documentation of the different ILSAs. It gets confusing also because different ILSAs (e.g., TIMSS vs. PISA) use different versions of these (e.g., JK2 vs Fay’s BRR).

Aside from ILSAs, many other surveys use replicate weights (e.g., ECLS-K). I’ll first show how we can use these replicate weights using base R, then using the survey package, then comparing results from an web app.

We’ll use the PISA 2022 dataset for Albania as an example:

pisa <- rio::import("https://raw.githubusercontent.com/flh3/pubdata/refs/heads/main/ILSA/alb_pisa.csv")
dim(pisa)

[1] 6129   95

names(pisa)

 [1] "CNT"         "CNTRYID"     "W_FSTUWT"    "PV1MATH"     "PV2MATH"     "PV3MATH"    
 [7] "PV4MATH"     "PV5MATH"     "PV6MATH"     "PV7MATH"     "PV8MATH"     "PV9MATH"    
[13] "PV10MATH"    "W_FSTURWT1"  "W_FSTURWT2"  "W_FSTURWT3"  "W_FSTURWT4"  "W_FSTURWT5" 
[19] "W_FSTURWT6"  "W_FSTURWT7"  "W_FSTURWT8"  "W_FSTURWT9"  "W_FSTURWT10" "W_FSTURWT11"
[25] "W_FSTURWT12" "W_FSTURWT13" "W_FSTURWT14" "W_FSTURWT15" "W_FSTURWT16" "W_FSTURWT17"
[31] "W_FSTURWT18" "W_FSTURWT19" "W_FSTURWT20" "W_FSTURWT21" "W_FSTURWT22" "W_FSTURWT23"
[37] "W_FSTURWT24" "W_FSTURWT25" "W_FSTURWT26" "W_FSTURWT27" "W_FSTURWT28" "W_FSTURWT29"
[43] "W_FSTURWT30" "W_FSTURWT31" "W_FSTURWT32" "W_FSTURWT33" "W_FSTURWT34" "W_FSTURWT35"
[49] "W_FSTURWT36" "W_FSTURWT37" "W_FSTURWT38" "W_FSTURWT39" "W_FSTURWT40" "W_FSTURWT41"
[55] "W_FSTURWT42" "W_FSTURWT43" "W_FSTURWT44" "W_FSTURWT45" "W_FSTURWT46" "W_FSTURWT47"
[61] "W_FSTURWT48" "W_FSTURWT49" "W_FSTURWT50" "W_FSTURWT51" "W_FSTURWT52" "W_FSTURWT53"
[67] "W_FSTURWT54" "W_FSTURWT55" "W_FSTURWT56" "W_FSTURWT57" "W_FSTURWT58" "W_FSTURWT59"
[73] "W_FSTURWT60" "W_FSTURWT61" "W_FSTURWT62" "W_FSTURWT63" "W_FSTURWT64" "W_FSTURWT65"
[79] "W_FSTURWT66" "W_FSTURWT67" "W_FSTURWT68" "W_FSTURWT69" "W_FSTURWT70" "W_FSTURWT71"
[85] "W_FSTURWT72" "W_FSTURWT73" "W_FSTURWT74" "W_FSTURWT75" "W_FSTURWT76" "W_FSTURWT77"
[91] "W_FSTURWT78" "W_FSTURWT79" "W_FSTURWT80" "male"        "ESCS"       

There are a lot of weights! The overall weight is W_FSTUWT and the 80 replicate weights are those that begin with W_FSTURWT. I explain elsewhere how these weights are constructed. Know that PISA uses Fay’s balanced repeated replication (BRR; with rho = .50). BRR is a type of resampling technique (along with jackknifing and bootstrapping). These replicate weights use the original sampling weights but are perturbed with some observations having weights increased by a factor of 1.5 and others decreased by a factor of .5. These perturbations will result in different results when we fit the model of interest.

In essence, to use the weights:

1. We fit the main model of interest using the overall student weight (W_FSTUWT) and we get the regression coefficients for this (we don’t need the standard errors here).
2. We then fit 80 models (R = 80) using W_FSTURWT1 to W_FSTURWT80 (instead of W_FSTUWT). These replicate weights will result in 80 sets of different regression coefficients per model fit (i.e., θ₁ to θ₈₀. We save all these coefficients too.
3. We use the 80 regression coefficients to compute the standard errors.

The standard errors are computed using:

$\frac{1}{R(1-k)^2}\Sigma(\theta\_{1...R} - \bar{\theta})^2$

For the case of PISA, k = .50 (or the Fay factor).The first part of the formula is an adjustment factor and reduces to

$\frac{1}{R(1-k)^2}=\frac{1}{80(1-.50)^2}=\frac{1}{20}$

Example 1: Manual method

Let’s say we are interested in seeing how SES and sex predict math outcomes (we’ll just use one plausible value as an example; I explain PVs elsewhere). We fit the model of interest:

NOTE: The lm function in R uses analytic or precision weights– not sampling weights!! The point estimates will be ok but the standard errors are incorrect! However– for this example- we don’t need the standard errors produced anyway! So, for step 1:

orig <- lm(PV1MATH ~ ESCS + male, data = pisa, weight = W_FSTUWT)
cfs <- coef(orig) #save the coefficients
cfs #the coefficients

(Intercept)        ESCS        male 
  392.65420    17.68895   -19.97603 

Now for step 2, we fit 80 models. We can just use a loop (easier to understand for most people):

cfs_i <- matrix(NA, 80, 3) #80 models of 3 coefficients
R <- 80 #number of replicate weights
for (i in 1:R){
  tmp <- lm(PV1MATH ~ ESCS + male, data = pisa, weight = pisa[, 13 + i]) 
  #13 + i are the columns for the replicate weights in the data
  cfs_i[i, ] <- coef(tmp)
}

Now cfs_i has the results from the 80 models. We can inspect this:

head(cfs_i)

         [,1]     [,2]      [,3]
[1,] 393.7959 17.95332 -19.58796
[2,] 391.8748 17.85683 -19.90734
[3,] 394.0433 19.00040 -19.70332
[4,] 394.0977 17.96452 -20.31570
[5,] 393.0942 18.24942 -21.16782
[6,] 393.3589 17.99746 -21.55690

Results will differ per replication. That variability is what allows us to compute the standard errors. We use the formula that was shown above. For example, if we want to get the standard error for ESCS (the second column in the matrix):

sqrt(sum((cfs_i[,2] - cfs[2])^2) * (1 / 20)) #of 1.63

[1] 1.633559

We can do this all at once though instead of one coefficient at a time:

SE <- sqrt(colSums(sweep(cfs_i, 2, cfs, "-")^2) * (1/20))
data.frame(cfs, SE, t = cfs / SE) #putting it all together

                  cfs       SE          t
(Intercept) 392.65420 2.633826 149.081273
ESCS         17.68895 1.633559  10.828471
male        -19.97603 2.595143  -7.697468

Example 2: Using the survey package

We can just use the survey package and compare results:

library(survey)
des <- svrepdesign(
  data = pisa,
  weights = ~W_FSTUWT,
  repweights = dplyr::select(pisa, starts_with("W_FSTURWT")),
  mse = T,
  type = 'Fay',
  rho = .5
)
s1 <- svyglm(PV1MATH ~ ESCS + male, design = des)
summary(s1)


Call:
svyglm(formula = PV1MATH ~ ESCS + male, design = des)

Survey design:
svrepdesign.default(data = pisa, weights = ~W_FSTUWT, repweights = dplyr::select(pisa, 
    starts_with("W_FSTURWT")), mse = T, type = "Fay", rho = 0.5)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  392.654      2.634 149.081  < 2e-16 ***
ESCS          17.689      1.634  10.828  < 2e-16 ***
male         -19.976      2.595  -7.697 3.88e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 6846.107)

Number of Fisher Scoring iterations: 2

Example 3: Use the app here

Try this webapp. After downloading the csv file from github, use that and compare results. They are the same.

NOTE: If you were to do this with TIMSS data, this will be slightly different procedure (the newer TIMSS datasets do not provide the replicate weights per column but the zone and replication info instead).

— END