0
$\begingroup$

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
Adj. R² = -417.1458419 

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?

Improve this question
$\endgroup$

1 Answer 1

Reset to default
2
$\begingroup$

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.

You ideally need to find out what the survey designers recommend as the standard error estimator. It's not as simple as "JKn" because 75 PSUs from 75 strata would give JKn weights with 75-75=0 design degrees of freedom, and that's not what you have.

  • Stata is refusing to give p-values because you've told it the design degrees of freedom are zero, so it's got no reference distribution to compare to. That's correct.

  • R works out the design df from the weights (the rank of the matrix of replicates), and doesn't get zero, so it's happy to give you p-values. That's also correct. I suspect this is closer to what the designers had in mind.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thanks! I wasn't using identical data (!). I fixed that, then added "jkrweight(SRWGT*) multiplier(1))" in Stata and "rscales=1" in R. When I reran the estimates, the coefficients and standard errors came up the same. On the DOF...well, clearly I am out of my depth. I understand what you're saying on Stata (I actually figured out that out after posting my question). But in R, I ran "degf(svdes)" and got 74 for some reason. In my data, 14 countries were divided up into 75 strata in each country, and two schools were sampled in each strata, and students were tested within those schools. Ideas? $\endgroup$
    – lotsofstats
    Commented May 8, 2021 at 2:28
  • $\begingroup$ 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 $\endgroup$
    – Thomas Lumley
    Commented May 8, 2021 at 5:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.