different results grouped survey standard deviation in Stata and R

Question

I have a question concerning the calculation of the grouped variance or standard deviation in R (survey-packge by Thomas Lumley) and Stata (using svyset and svy prefix). I'd like to know if there ara different methods used for the calculation as the grouped variances/sds are different in Stata and R.
It seems the difference is limited to the estimation of the grouped standard deviation and variance. Because grouped means are the same in both cases and also the non grouped standard deviation/variance.

I would be thankful for every idea concerning the difference in estimation. Here is my R and Stata code.

R

#libraries
library(srvyr)
library(survey)

#load data
mtcars <- read.table("https://forge.scilab.org/index.php/p/rdataset/source/file/master/csv/datasets/mtcars.csv",
              sep=".", 
              header=TRUE)

mtcars_cplx <- mtcars %>% as_survey_design(id = cyl, weights = qsec)

# ungrouped standard deviation
var_mpg <- svyvar(~mpg, mtcars_cplx) #variance
sd_mpg<-sqrt(var_mpg) #standard deviation
sd_mpg

# grouped standard deviation
var_grouped_mpg <- svyby(~mpg, ~vs+am, mtcars_cplx, svyvar, na.rm=TRUE, na.rm.all = TRUE) #variance
sd_grouped_mpg <- sqrt(var_grouped_mpg) # standard deviation
sd_grouped_mpg

#grouped mean
mean_grouped_mpg <- svyby(~mpg, ~vs+am, mtcars_cplx, svymean, na.rm=TRUE, na.rm.all = TRUE)
mean_grouped_mpg

Stata

* load data
import delim using "https://forge.scilab.org/index.php/p/rdataset/source/file/master/csv/datasets/mtcars.csv"
svyset cyl [pweight=qsec]

* mean and standard deviation
svy: mean mpg
estat sd

* grouped mean and standard deviation
svy: mean mpg, over(vs am)
estat sd

R results

grouped mean

    vs am      mpg           se
0.0  0  0 15.03757 1.665335e-16
1.0  1  0 20.82073 1.617214e+00
0.1  0  1 19.94862 2.190127e+00
1.1  1  1 28.42264 1.776357e-15

grouped sd

    vs am      mpg           se
0.0  0  0 2.799371 1.053671e-08
1.0  1  0 2.471672 6.100837e-01
0.1  0  1 3.988702 3.429033e+00
1.1  1  1 4.788752 5.161914e-08

Stata results

grouped mean and sd

-------------------------------------
        Over |       Mean   Std. Dev.
-------------+-----------------------
mpg          |
   _subpop_1 |   15.03757    2.778608
   _subpop_2 |   20.82073    2.198151
   _subpop_3 |   19.94862    3.932388
   _subpop_4 |   28.42264    4.400745
-------------------------------------

Thank your for your help Stephan

Answer 1

What's happening is that R and Stata use quite different approaches to estimating the population standard deviation. R goes directly via the mean of $(X-\bar X)^2$ and Stata goes via an estimated standard error of the mean (under simple random sampling).

As a result, Stata's estimate depends on values outside the subpopulation (as is unarguably correct for the standard error of the mean) and R's doesn't (as is at least arguably correct for the subpopulation standard deviation).

For example, you can get the 2.198 in the Stata output as

> sqrt(2.47162^2*571.16/139.77*6/31)
[1] 2.198105

where 2.47162 is the R estimate, 571.16 is the estimated population size (sum of weights), 139.77 is the estimated subpopulation size (sum of weights), 6 is 7-1 and 31 is 32-1.

As you can see, it isn't easy to convert; it's probably easier to just implement the formula from the Stata manual if you want the Stata version.