Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Directly store replicate weights #168

Merged
merged 11 commits into from
Jan 10, 2023

Conversation

ayushpatnaikgit
Copy link
Member

  1. Replicate weights are directly stored instead of storing the scaling factor and asking functions to do the multiplication later.

Tests and doctests are failing. Need to be changed.

@smishr
Copy link
Contributor

smishr commented Jan 4, 2023

you want to fix before merging?

Copy link
Contributor

@smishr smishr left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

looks good. the small mistake in the total.jl also resolved.

Copy link
Contributor

@smishr smishr left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

adding some missing checks? with disallowmissing

@iuliadmtru
Copy link
Contributor

I am working on fixing the tests and adding more tests for this functionality. I will convert this to a draft PR until that is done.

@iuliadmtru iuliadmtru marked this pull request as draft January 9, 2023 07:33
@iuliadmtru
Copy link
Contributor

iuliadmtru commented Jan 9, 2023

Does bootstrapping take into account fpc? In the code fpc is not used, but I understand that bootstrapping might account for fpc implicitly. Is that true?

This is relevant for testing. I am comparing our results to the ones in R. What is the correct way to compare?

Let's take total and mean for a simple random sample using the apisrs dataset for example. The Julia code for generating the design and the estimates is:

julia> apisrs = load_data("apisrs");

julia> srs = SurveyDesign(apisrs; weights = :pw) |> bootweights;

julia> tot = total(:api00, srs)
1×2 DataFrame
 Row │ total      SE
     │ Float64    Float64
─────┼────────────────────
   14.06689e6  57644.8

julia> mn = mean(:api00, srs)
1×2 DataFrame
 Row │ mean     SE
     │ Float64  Float64
─────┼──────────────────
   1656.585  9.30656

For R, I believe the corresponding code is:

> srs <- svydesign(data=apisrs, id=~1, weights=~pw)
> srsrep <- as.svrepdesign(srs, type="bootstrap", replicates=4000)
> svytotal(~api00, srsrep)
        total    SE
api00 4066888 56295
> svymean(~api00, srsrep)
        mean     SE
api00 656.58 9.0886

However, specifying fpc gives results closer to ours:

> srs <- svydesign(data=apisrs, id=~1, weights=~pw, fpc=~fpc)
> srsrep <- as.svrepdesign(srs, type="bootstrap", replicates=4000)
> svytotal(~api00, srsrep)
        total    SE
api00 4066888 57549
> svymean(~api00, srsrep)
        mean     SE
api00 656.58 9.2912

While this is not an indication that the latter is the actual corresponding code, it does make me wonder. Which one is the correct approach that I should use for testing?

@iuliadmtru
Copy link
Contributor

iuliadmtru commented Jan 9, 2023

Also, regarding testing, so far we've been using absolute tolerance to check correctness as compared to the equivalent R results. However, we should consider using relative tolerance instead.

For estimated standard errors we should use rtol and for estimated statistics, we should use atol (in theory) since we expect the results to be exactly the same as those in R. Hence, say, atol = 1e-3 should be a good tolerance for testing. However, in some cases R and Julia don't agree on rounding methods and the results end up within an absolute tolerance of 1. An example is:

Julia:

julia> apisrs = load_data("apisrs");

julia> srs = SurveyDesign(apisrs; weights = :pw) |> bootweights

julia> tot = total(:api00, srs)
1×2 DataFrame
 Row │ total      SE
     │ Float64    Float64
─────┼────────────────────
   14.06689e6  57644.8

R:

> srs <- svydesign(data=apisrs, id=~1, weights=~pw)
> srsrep <- as.svrepdesign(srs, type="bootstrap", replicates=4000)
> svytotal(~api00, srsrep)
        total    SE
api00 4066888 56295

An absolute tolerance of 1 is too large for most cases. Hence, I suggest we use relative tolerance for estimated statistics, as well as for estimated standard errors. This has the advantage of consistency.

There is also the case when the result is close to 0. We might get this especially for standard errors. In this case, the relative tolerance would be too high and we should use atol.

@iuliadmtru
Copy link
Contributor

ReplicateDesign should be supported by ratio and quantile as well. @ayushpatnaikgit should we do this in a separate PR after this? And add tests for ratio and quantile in that PR.

@iuliadmtru
Copy link
Contributor

We should also add back the CategoricalArray support.

@ayushpatnaikgit
Copy link
Member Author

Does bootstrapping take into account fpc? In the code fpc is not used, but I understand that bootstrapping might account for fpc implicitly. Is that true?

This is relevant for testing. I am comparing our results to the ones in R. What is the correct way to compare?

Let's take total and mean for a simple random sample using the apisrs dataset for example. The Julia code for generating the design and the estimates is:

julia> apisrs = load_data("apisrs");

julia> srs = SurveyDesign(apisrs; weights = :pw) |> bootweights;

julia> tot = total(:api00, srs)
1×2 DataFrame
 Row │ total      SE
     │ Float64    Float64
─────┼────────────────────
   14.06689e6  57644.8

julia> mn = mean(:api00, srs)
1×2 DataFrame
 Row │ mean     SE
     │ Float64  Float64
─────┼──────────────────
   1656.585  9.30656

For R, I believe the corresponding code is:

> srs <- svydesign(data=apisrs, id=~1, weights=~pw)
> srsrep <- as.svrepdesign(srs, type="bootstrap", replicates=4000)
> svytotal(~api00, srsrep)
        total    SE
api00 4066888 56295
> svymean(~api00, srsrep)
        mean     SE
api00 656.58 9.0886

However, specifying fpc gives results closer to ours:

> srs <- svydesign(data=apisrs, id=~1, weights=~pw, fpc=~fpc)
> srsrep <- as.svrepdesign(srs, type="bootstrap", replicates=4000)
> svytotal(~api00, srsrep)
        total    SE
api00 4066888 57549
> svymean(~api00, srsrep)
        mean     SE
api00 656.58 9.2912

While this is not an indication that the latter is the actual corresponding code, it does make me wonder. Which one is the correct approach that I should use for testing?

svymean(~api00, srsrep)

I tried it with 100,000 replicates. The answers are very close to the analytical solutions of with and without fpc. So fpc matters.

@ayushpatnaikgit
Copy link
Member Author

We can merge the PR once the tests and doctests are complete.

@smishr
Copy link
Contributor

smishr commented Jan 9, 2023

We should also add back the CategoricalArray support.

In separate PR, not this one

@smishr
Copy link
Contributor

smishr commented Jan 9, 2023

ReplicateDesign should be supported by ratio and quantile as well. @ayushpatnaikgit should we do this in a separate PR after this? And add tests for ratio and quantile in that PR.

Yes this is higher priority than CategoricalArray.

@codecov-commenter
Copy link

codecov-commenter commented Jan 9, 2023

Codecov Report

Merging #168 (76275cd) into singledesign (eb290eb) will increase coverage by 0.18%.
The diff coverage is 96.00%.

@@               Coverage Diff                @@
##           singledesign     #168      +/-   ##
================================================
+ Coverage         64.61%   64.79%   +0.18%     
================================================
  Files                13       13              
  Lines               195      196       +1     
================================================
+ Hits                126      127       +1     
  Misses               69       69              
Impacted Files Coverage Δ
src/bootstrap.jl 95.65% <90.90%> (+0.19%) ⬆️
src/SurveyDesign.jl 92.30% <100.00%> (ø)
src/by.jl 100.00% <100.00%> (ø)
src/mean.jl 100.00% <100.00%> (ø)
src/total.jl 100.00% <100.00%> (ø)

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

@iuliadmtru
Copy link
Contributor

Tests now pass. I will soon get into documentation and doctests.

We should consider designing more R-independent tests.

@smishr
Copy link
Contributor

smishr commented Jan 10, 2023

documentation still failing, as not all references to analytical solutions removed in this branch. that and better docs with current branch.

merge this now doc pr separate

@ayushpatnaikgit
Copy link
Member Author

@iuliadmtru can you update README.md?

@iuliadmtru
Copy link
Contributor

iuliadmtru commented Jan 10, 2023

All tests and doctests should pass now. I made some modifications:

  • All text from api.md is gone. I believe that text is better fit for index.md in the tutorial/demo section. We should also have more complete docstrings where we write all the relevant information. Right now our dosctrings are lacking substance.
  • There are no more references to the old designs (SimpleRandomSample, StratifiedSample) in the documentation.
  • The sturges and freedman_diaconis functions are exported again. They are necessary for hist. However, I don't think they are relevant enough to our package to include them in the API reference. Therefore I didn't include them. Documenter now gives a warning because of this and I don't think we can avoid it.
  • The files that I changed now have a coherent style.

There are still a few problems left:

I think we can now close this PR and have a clean source code once again so that the tests for our following PRs don't fail because of this one.

@iuliadmtru iuliadmtru marked this pull request as ready for review January 10, 2023 12:08
@iuliadmtru iuliadmtru requested a review from smishr January 10, 2023 12:08
Copy link
Contributor

@smishr smishr left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

looks good to merge

@smishr smishr merged commit 7e946cc into xKDR:singledesign Jan 10, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants