# Differing standard errors with jackknife dataset - R compared with Stata

I am working with a stratified international assessment (14 countries from IEA's ICCS 2016 assessment) and am trying to understand why R's survey package and Stata's survey module differ slightly in their results. As you can see below, the difference in minor, but they are different.

The dataset I'm working with has a set of 75 jackknife leave-one-out replicate weights which appear as separate columns in the dataset (SRWGT1-SRWGT75). There is one replicate for each of the 75 strata (which are provided in a variable called JKZONES). I am comparing Stata and R by running a simple linear regression using the jackknife replicates. My point estimates are the same between the two packages, but SE and t-values are slightly different.

Stata commands:

svyset IDSCHOOL [pweight=TOTWGTS], strata(JKZONES) vce(jackknife) jkrweight(SRWGT*) mse
svy: regress S_INTRUST S_NISB


Stata results (no p-values or confidence intervals are produced):

Survey: Linear regression

Number of strata = 75                              Number of obs   =    30,716
Population size = 1,048,437
Replications    =        75
Design df       =         0
F(1, 0)         =         .
Prob > F        =         .
R-squared       =    0.0062

------------------------------------------------------------------------------
|              Jknife *
S_INTRUST | Coefficient  std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
S_NISB |   .6750581   .0714951     9.44       .            .           .
_cons |   52.20439   .1077827   484.35       .            .           .
------------------------------------------------------------------------------


R commands (I didn't see a place for identifying the variables for the unique PSU IDs or the strata in the svrepdesign help section):

svdes <- svrepdesign(
data = svdata,
type = "JKn" ,
repweights = "SRWGT[0-9]" ,
weights = ~TOTWGTS,
rscales=rep(1, 75),
mse = TRUE)
summ(svyglm(S_INTRUST ~ S_NISB, design=svdes), digits=7)


R survey results (produces p-values, unlike my Stata commands):

MODEL INFO:
Observations: 30716
Dependent Variable: S_INTRUST
Type: Survey-weighted linear regression

MODEL FIT:
R² = 0.0061974

Standard errors: Robust
--------------------------------------------------------------------
Est.        S.E.        t val.           p
----------------- ------------ ----------- ------------- -----------
(Intercept)         52.2038545   0.1077395   484.5378497   0.0000000
S_NISB               0.6769446   0.0721005     9.3889068   0.0000000
--------------------------------------------------------------------

Estimated dispersion parameter = 2280595


Does anyone have any suggestions on improving the way I am entering the commands? Or does anyone know why the two sets of results don't match?

You say the point estimates are the same but the SE and $$t$$-values are different. Actually, the point estimates agree about as well as the standard errors do. Which is strange -- I'd expect the point estimates to be identical. I'd expect the SEs and t-values to be slightly different, because you have specified rscales=1 in R but not multiplier(1) in Stata's jkrw. I would check very carefully that you have identical datasets in R and Stata, given that the point estimates are not identical and this isn't an iterative computation.
• degf is actually defined as one less than the rank of the matrix of replicates, so if you have 75 strata you might well get 74 as the answer. There isn't a perfect analogue of the PSUs-strata rule for replicate-weight designs, but R's rule gives reasonable answers across a wide range of settings -- and you can specify df.resid to summary() if you don't like the defaults May 8, 2021 at 5:33