Using Replicate Weights in R

Posted: Sep 8, 2026

I was preparing for a class and I had AI come up with a small app that calculates standard errors using replicate weights. The app runs entirely in the browser, so your data never leaves your device. It supports BRR/Fay, JK1, JK2, and bootstrap replicate variance estimation, as well as using plausible values.

I haven’t tested the app thoroughly, so double check your results!

Here’s an example of how to do this using R (using a PISA data)

Example 1:

> library(survey) #to fit the models
> library(dplyr)
> pisa2 <- read.csv("https://francish.net/post/2026_repwgt_app/example.csv")
> 
> des <- svrepdesign(
  data = pisa2,
  repweights = select(pisa2, starts_with("W_FSTR")),
  weights =~W_FSTUWT,
  rho = .5,
  type = 'Fay'
)
> 
> library(mitools) #to combine PV results using Rubin's rules
> 
> pv_map <- list(math ~ PV1MATH + PV2MATH + PV3MATH + PV4MATH + PV5MATH)
> 
> ## just get the mean
> pv_models <- withPV(
  mapping = pv_map,
  data = des, #indicate the design
  action = function(x) svymean(~math, design = x) #x is the design
)
> MIcombine(pv_models) %>% summary()
Multiple imputation results:
      withPV.svyrep.design(mapping = pv_map, data = des, action = function(x) svymean(~math, 
    design = x))
      MIcombine.default(pv_models)
         results       se  (lower   upper) missInfo
PV1MATH 537.8233 3.127436 531.693 543.9536      2 %
> 
> ## try a regression
> pv_models <- withPV(
  mapping = pv_map,
  data = des, 
  action = function(x) svyglm(math ~ ESCS, design = x) 
)
> 
> MIcombine(pv_models) %>% summary()
Multiple imputation results:
      withPV.svyrep.design(mapping = pv_map, data = des, action = function(x) svyglm(math ~ 
    ESCS, design = x))
      MIcombine.default(pv_models)
              results       se    (lower    upper) missInfo
(Intercept) 539.00347 2.454642 534.19024 543.81669      4 %
ESCS         44.73243 2.362189  40.09589  49.36896      7 %

For example 2:


> library(survey)
> df2 <- read.csv("https://francish.net/post/2026_repwgt_app/timss_ex2.csv")
> # des <- svrepdesign(
> #   data = df2,
> #   repweights = select(df2, starts_with("rw")),
> #   weights =~totwgt,
> #   type = 'JK2'
> # ) #alternative using generated repweights
> 
> des <- svydesign(
   data = df2,
   strata = ~zone, #JKZONE
   ids = ~member, #JKREP
   weights = ~totwgt,
   mse = TRUE,
   nest = TRUE
+ )
> mod1 <- svyglm(y ~ w1 + x1 + minor, des)
> summary(mod1)

Call:
svyglm(formula = y ~ w1 + x1 + minor, design = des)

Survey design:
svydesign(data = df2, strata = ~zone, ids = ~member, weights = ~totwgt, 
    mse = TRUE, nest = TRUE)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.47118    0.05552   8.486 1.71e-12 ***
w1           0.59031    0.05613  10.517 2.85e-16 ***
x1           0.49160    0.02196  22.391  < 2e-16 ***
minor       -0.42578    0.06302  -6.756 2.94e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 1.700741)

Number of Fisher Scoring iterations: 2

Cite as:

Huang, F. L. (2026, September 8). Using Replicate Weights in R. FLH Website. https://francish.net/post/2026_rw2/

Using Replicate Weights in R

Related