[R] Where to find the source code for survey::SE

[Anthony Damico]

hi, are you looking for the code in survey:::svymean.survey.design2   ?
running `debug(svymean)` and `debug(svyrecvar)` might help you step through
the calculations..



On Tue, Jul 5, 2022 at 8:29 AM Iulia Dumitru <iuliadmtru using gmail.com> wrote:

> Hello! I want to port some functionality from the Survey package to Julia
> and I am stuck on standard error calculation. How is the standard error
> calculated in the Survey package?
>
> I searched in the can/survey repo on GitHub and I couldn’t find the source
> code of the SE function. I also looked through Thomas Lumley’s book Complex
> Surveys A Guide to Analysis Using R and from there I understood that the
> standard error is calculated by first calculating the Horvitz-Thompson
> variance estimator and then taking the square root of that, but I get
> completely different results.
>
> This is the R code I used:
>
> > library(survey)
> > data(api)
> > dclus1 <- svydesign(id = ~dnum, weights = ~pw, data = apiclus1, fpc =
> ~fpc)
> > svyby(~api00, by = ~cname, design = dclus1, svymean)
>
> And this is what I obtain:
>
>                   cname    api00           se
> Alameda         Alameda 669.0000 1.264044e-15
> Fresno           Fresno 472.0000 0.000000e+00
> Kern               Kern 452.5000 0.000000e+00
> Los Angeles Los Angeles 647.2667 1.724959e+01
> Mendocino     Mendocino 623.2500 0.000000e+00
> Merced           Merced 519.2500 0.000000e+00
> Orange           Orange 710.5625 0.000000e+00
> Plumas           Plumas 709.5556 1.248655e-13
> San Diego     San Diego 659.4364 2.296800e+00
> San Joaquin San Joaquin 551.1892 1.354175e-13
> Santa Clara Santa Clara 732.0769 1.577292e+01
>
>
> With Julia this is the result I get:
>
> 11×3 DataFrame
>  Row │ cname        mean     se
>      │ String15     Float64  Float64
> ─────┼────────────────────────────────
>    1 │ Alameda      669.0     6469.11
>    2 │ Fresno       472.0     7832.61
>    3 │ Kern         452.5    10703.1
>    4 │ Los Angeles  647.267   5232.24
>    5 │ Mendocino    623.25   10364.8
>    6 │ Merced       519.25    8616.22
>    7 │ Orange       710.563   5534.39
>    8 │ Plumas       709.556   7673.15
>    9 │ San Diego    659.436   1882.06
>   10 │ San Joaquin  551.189   2254.24
>   11 │ Santa Clara  732.077   4000.03
>
> And just to have the full picture, this is my code for the
> Horvitz-Thompson variance estimator:
>
> function hte(y, pw, n)
>         p = 1 .- (1 .- 1 ./ pw) .^ n
>
>         first_sum = sum((1 .- p) ./ (p .^ 2) .* y .^ 2)
>         second_sum = 0
>
>         for i in 1:length(p)
>                 for j in 1:length(p)
>                         if i == j
>                                 continue
>                         end
>
>
>                         p_ij = p[i] + p[j] - (1 - (1 - 1 / pw[i] - 1 /
> pw[j]) ^ n)
>                         second_sum += (p_ij - p[i] * p[j]) / (p[i] * p[j])
> * y[i] * y[j]
>                 end
>         end
>
>         return first_sum + second_sum
> end
>
> The standard error in the resulted data frame above is just the square
> root of the hte result.
>         [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

	[[alternative HTML version deleted]]