diff --git a/.gitignore b/.gitignore index 47d3741a..e4441b23 100644 --- a/.gitignore +++ b/.gitignore @@ -5,3 +5,7 @@ /docs/Manifest.toml /docs/build/ /test/Manifest.toml +/dev/* +.gitignore +.DS_Store +*.json diff --git a/Project.toml b/Project.toml index 1cda24c2..70c896c9 100644 --- a/Project.toml +++ b/Project.toml @@ -7,8 +7,8 @@ version = "0.11.1" AlgebraOfGraphics = "cbdf2221-f076-402e-a563-3d30da359d67" CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" +CategoricalArrays = "324d7699-5711-5eae-9e2f-1d82baa6b597" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" -GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" @@ -18,6 +18,5 @@ AlgebraOfGraphics = "0.6" CSV = "0.10" CairoMakie = "0.8, 0.9, 0.10" DataFrames = "1" -GLM = "1" StatsBase = "0.33" julia = "1" diff --git a/README.md b/README.md index 4865a73e..ac9cdb69 100644 --- a/README.md +++ b/README.md @@ -7,121 +7,104 @@ [![Milestones](https://img.shields.io/badge/-milestones-brightgreen)](https://github.com/xKDR/Survey.jl/milestones) -This package is used to study complex survey data. It is the Julia implementation of the [Survey package in R](https://cran.r-project.org/web/packages/survey/index.html) developed by [Professor Thomas Lumley](https://www.stat.auckland.ac.nz/people/tlum005). +This package is used to study complex survey data. It aims to be a fast alternative to the [Survey package in R](https://cran.r-project.org/web/packages/survey/index.html) developed by [Professor Thomas Lumley](https://www.stat.auckland.ac.nz/people/tlum005). -As the size of survey datasets have become larger, processing the records can take hours or days in R. We endeavour to solve this problem by implementing the Survey package in Julia. +This package currently supports simple random sample and stratified sample. In future releases, it will support multistage sampling as well. ## How to install - - add "https://github.com/xKDR/Survey.jl.git" - -## Basic usage - -In the following example, we will load the Academic Performance Index dataset for Californian schools and produce the weighted mean for each county. ```julia -using Survey - -apiclus1 = load_data("apiclus1") -## This function loads a commonly used dataset, Academic Performance Index (API), as an example. -## Any DataFrame object can be used with this package. - -dclus1 = svydesign(id = :1, weights = :pw, data = apiclus1) - -svyby(:api00, :cname, dclus1, svymean) -11×3 DataFrame - Row │ cname mean SE - │ String15 Float64 Float64 -─────┼──────────────────────────────── - 1 │ Alameda 669.0 16.2135 - 2 │ Fresno 472.0 9.85278 - 3 │ Kern 452.5 29.5049 - 4 │ Los Angeles 647.267 23.5116 - 5 │ Mendocino 623.25 24.216 - 6 │ Merced 519.25 10.4925 - 7 │ Orange 710.562 28.9123 - 8 │ Plumas 709.556 13.2174 - 9 │ San Diego 659.436 12.2082 - 10 │ San Joaquin 551.189 11.578 - 11 │ Santa Clara 732.077 12.2291 -``` - -This example is from the Survey package in R. The [examples section of the documentation](https://xkdr.github.io/Survey.jl/dev/examples/) shows the R and the Julia code side by side for this and a few other examples. - -## Performance -We will measure the performance of the R and Julia for the example shown above. - -**R** - -```R -library(survey) -library(microbenchmark) -data(api) -dclus1 <- svydesign(id = ~1, weights = ~pw, data = apiclus1) -microbenchmark(svyby(~api00, by = ~cname, design = dclus1, svymean), units = "us") +] add "https://github.com/xKDR/Survey.jl.git" ``` +## Basic usage -```R - expr min lq - svyby(~api00, by = ~cname, design = dclus1, svymean) 10180.47 12102.61 - mean median uq max neval - 12734.43 12421.93 12788.55 17242.35 100 -``` +### Simple Random Sample -**Julia** +In the following example, we will load a simple random sample of the Academic Performance Index dataset for Californian schools and do basic analysis. ```julia -using Survey, BenchmarkTools -apiclus1 = load_data("apiclus1") -dclus1 = svydesign(id=:1, weights=:pw, data = apiclus1) -@benchmark svyby(:api00, :cname, dclus1, svymean) -``` +using Survey -```julia -BenchmarkTools.Trial: 10000 samples with 1 evaluation. - Range (min … max): 54.464 μs … 6.070 ms ┊ GC (min … max): 0.00% … 94.01% - Time (median): 72.468 μs ┊ GC (median): 0.00% - Time (mean ± σ): 81.833 μs ± 190.657 μs ┊ GC (mean ± σ): 7.62% ± 3.23% - ``` - -The Julia code is about 171 times faster than the R code. - -We increase the complexity by grouping the data by two variables and then performing the same operations. -**R** - -```R -library(survey) -library(microbenchmark) -data(api) -dclus1 <- svydesign(id = ~1, weights = ~pw, data = apiclus1) -microbenchmark(svyby(~api00, by = ~cname+meals, design = dclus1, svymean, keep.var = FALSE), units = "us") +srs = load_data("apisrs") + +dsrs = SimpleRandomSample(srs; weights = :pw) + +mean(:api00, dsrs) +1×2 DataFrame + Row │ mean SE + │ Float64 Float64 +─────┼────────────────── + 1 │ 656.585 9.24972 + +total(:enroll, dsrs) +1×2 DataFrame + Row │ total SE + │ Float64 Float64 +─────┼───────────────────── + 1 │ 3.62107e6 1.6952e5 + +mean(:api00, :cname, dsrs) +38×3 DataFrame + Row │ cname mean SE + │ String15 Float64 Float64 +─────┼──────────────────────────────────── + 1 │ Kern 573.6 42.8026 + 2 │ Los Angeles 658.156 21.0728 + 3 │ Orange 749.333 27.0613 + ⋮ │ ⋮ ⋮ ⋮ + 36 │ Napa 727.0 46.722 + 37 │ Lake 804.0 NaN + 38 │ Merced 595.0 NaN + +quantile(:enroll,dsrs,[0.1,0.2,0.5,0.75,0.95]) +5×2 DataFrame + Row │ probability quantile + │ Float64 Float64 +─────┼─────────────────────── + 1 │ 0.1 245.5 + 2 │ 0.2 317.6 + 3 │ 0.5 453.0 + 4 │ 0.75 668.5 + 5 │ 0.95 1473.1 ``` -```R -Unit: microseconds - expr min lq - svyby(~api00, by = ~cname + meals, design = dclus1, svymean) 132468.1 149914 - mean median uq max neval - 166121.9 160571.3 172301.6 304979.2 100 -``` +### Stratified Sample -**Julia** -```julia -using Survey, BenchmarkTools -apiclus1 = load_data("apiclus1") -dclus1 = svydesign(id=:1, weights=:pw, data = apiclus1) -@benchmark svyby(:api00, [:cname, :meals], dclus1, svymean) -``` +In the following example, we will load a stratified sample of the Academic Performance Index dataset for Californian schools and do basic analysis. ```julia -BenchmarkTools.Trial: 10000 samples with 1 evaluation. - Range (min … max): 219.387 μs … 8.284 ms ┊ GC (min … max): 0.00% … 90.94% - Time (median): 265.214 μs ┊ GC (median): 0.00% - Time (mean ± σ): 325.100 μs ± 513.020 μs ┊ GC (mean ± σ): 14.23% ± 8.58% - ``` +using Survey -The Julia code is about 605 times faster than the R code. +strat = load_data("apistrat") + +dstrat = StratifiedSample(strat, :stype; weights = :pw, popsize = :fpc) + +mean(:api00, dstrat) +1×2 DataFrame + Row │ mean SE + │ Float64 Float64 +─────┼────────────────── + 1 │ 662.287 9.40894 + +total(:api00, dstrat) +1×2 DataFrame + Row │ total SE + │ Float64 Float64 +─────┼──────────────────── + 1 │ 4.10221e6 58279.0 + +mean(:api00, :cname, dstrat) +40×3 DataFrame + Row │ cname mean SE + │ String15 Float64 Float64 +─────┼─────────────────────────────────────── + 1 │ Los Angeles 633.511 21.3912 + 2 │ Ventura 707.172 31.6856 + 3 │ Kern 678.235 53.1337 + ⋮ │ ⋮ ⋮ ⋮ + 39 │ Mendocino 632.018 1.04942 + 40 │ Butte 627.0 0.0 +``` ## Strategic goals - We want to implement all the features provided by the [Survey package in R](https://cran.r-project.org/web/packages/survey/index.html) The [milestones](https://github.com/xKDR/Survey.jl/milestones) sections of the repository contains a list of features that contributors can implement in the short-term. diff --git a/docs/make.jl b/docs/make.jl index 8fef3fa2..df10b794 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -15,10 +15,10 @@ makedocs(; ), pages=[ "Home" => "index.md", - "Examples" => "examples.md", - "Comparison with R" => "R_comparison.md", - "Performance" => "performance.md", + "Moving from R" => "R_comparison.md", + "API reference" => "api.md" ], + checkdocs=:exports, ) deploydocs(; diff --git a/docs/src/R_comparison.md b/docs/src/R_comparison.md index 3fe15a6c..1d1a4b2d 100644 --- a/docs/src/R_comparison.md +++ b/docs/src/R_comparison.md @@ -1,98 +1,131 @@ -# Comparison with R +# Moving from R to Julia +This section presents examples to help move from R to Julia. Examples show R and Julia code for common operations in survey analysis.
+For the same operation, first the R and then the Julia code is presented. -In the following examples, we'll compare Julia's performance to R's on the same set of operations. +## Simple random sample -## Installing and loading the package -**R** +The `apisrs` data, which is provided in both `survey` and `Survey.jl`, is used as an example. It's a simple random sample of the Academic Performance Index of Californian schools. -```r -install.package("survey") +### 1. Creating a survey design +Instantiating a simple random sample survey design. + +```R library(survey) +data(api) +dsrs = svydesign(id = ~1, data = apisrs, weights = ~pw, fpc = ~fpc) ``` -**Julia** ```julia -using Pkg -Pkg.add(url = "https://github.com/xKDR/Survey.jl.git") using Survey +srs = load_data("apisrs") +dsrs = SimpleRandomSample(srs; popsize = :fpc) ``` -The following command in the Pkg REPL may also be used to install the package. +### 2. Mean +In the following example the mean of the variable `api00` is calculated. + +```R +svymean(~api00, dsrs) ``` -add "https://github.com/xKDR/Survey.jl.git" +```julia +mean(:api00, dsrs) ``` -## API data +### 3. Total +In the following example the sum of the variable `api00` is calculated. -The Academic Performance Index is computed for all California schools based on standardised testing of students. The [data sets](https://cran.r-project.org/web/packages/survey/survey.pdf) contain information for all schools with at least 100 students and for various probability samples of the data. apiclus1 is a cluster sample of school districts, apistrat is a sample stratified by stype. +```R +svytotal(~api00, dsrs) +``` +```julia +total(:api00, dsrs) +``` -In the following examples, we'll use the apiclus1 data from the api dataset. +### 4. Quantile +In the following example the median of the variable `api00` is calculated. +```R +svyquantile(~api00, dsrs, 0.5) +``` +```julia +quantile(:api00, dsrs, 0.5) +``` -The api dataset can be loaded using the following command: +### 5. Domain estimation +In the following example the mean of the variable `api00` is calculated grouped by the variable `cname`. -**R** -```r -data(api) +```R +svyby(~api00, ~cname, dsrs, svymean) ``` -**Julia** ```julia -apiclus1 = load_data("apiclus1") +mean(:api00, :cname, dsrs) ``` -## svydesign -[The ```svydesign``` object combines a data frame and all the survey design information needed to analyse it.](https://www.rdocumentation.org/packages/survey/versions/4.1-1/topics/svydesign) - -A ```svydesign``` object can be constructed with the following command: +In the following example the total of the variable `api00` is calculated grouped by the variable `cname`. -**R** -```r -dclus1 <-svydesign(id = ~1, weights = ~pw, data = apiclus1, fpc = ~fpc) +```R +svyby(~api00, ~cname, dsrs, svytotal) ``` -**Julia** ```julia -dclus1 = svydesign(id = :1, weights = :pw, data = apiclus1, fpc = :fpc) +total(:api00, :cname, dsrs) ``` -## svyby -The `svyby` function can be used to generate stratified estimates. +## Stratified sample + +The `apistrat` data, which is provided in both `survey` and `Survey`, is used as an example. It's a stratified sample of the Academic Performance Index of Californian schools. -### Mean -Weighted mean of a variable by strata can be computed using the following command: +### 1. Creating a design object +The following example shows how to construct a design object for a stratified sample. -**R** -```r -svyby(~api00, by = ~cname, design = dclus1, svymean) +```R +library(survey) +data(api) +dstrat = svydesign(id = ~1, data = apistrat, strata = ~stype, weights = ~pw, fpc = ~fpc) ``` -**Julia** ```julia -svyby(:api00, :cname, dclus1, svymean) +using Survey +strat = load_data("apistrat") +dstrat = StratifiedSample(strat, :stype; popsize = :fpc) ``` -### Sum -Weighted sum of a variable by strata can be computed using the following command: +### 2. Mean +In the following example the mean of the variable `api00` is calculated. -**R** -```r -svyby(~api00, by = ~cname, design = dclus1, svytotal) +```R +svymean(~api00, dstrat) ``` - -**Julia** ```julia -svyby(:api00, :cname, dclus1, svytotal) +mean(:api00, dstrat) ``` -### Quantile -Weighted quantile of a variable by strata can be computed using the following command: +### 3. Total +In the following example the sum of the variable `api00` is calculated. -**R** -```r -svyby(~api00, by = ~cname, design = dclus1, svyquantile, quantile = 0.63) +```R +svytotal(~api00, dstrat) +``` +```julia +total(:api00, dstrat) ``` -**Julia** +### 4. Quantile +In the following example the median of the variable `api00` is calculated. +```R +svyquantile(~api00, dstrat, 0.5) +``` ```julia -svyby(:api00, :cname, dclus1, svyquantile, 0.63) +quantile(:api00, dstrat, 0.5) +``` + +### 5. Domain estimation +In the following example the mean of the variable `api00` is calculated grouped by the variable `cname`. + +```R +svyby(~api00, ~cname, dstrat, svymean) ``` + +```julia +mean(:api00, :cname, dstrat) +``` \ No newline at end of file diff --git a/docs/src/api.md b/docs/src/api.md new file mode 100644 index 00000000..ea9f9bfe --- /dev/null +++ b/docs/src/api.md @@ -0,0 +1,40 @@ +# API + +## Index + +```@index +Module = [Survey] +Order = [:type, :function] +Private = false +``` +Survey data can be loaded from a `DataFrame` into a survey design. The package currently supports simple random sample and stratified sample designs. +```@docs +AbstractSurveyDesign +SimpleRandomSample +StratifiedSample +``` + +```@docs +load_data +Survey.mean(x::Symbol, design::SimpleRandomSample) +total(x::Symbol, design::SimpleRandomSample) +quantile +``` + +It is often required to estimate population parameters for sub-populations of interest. For example, you may have a sample of heights, but you want the average heights of males and females separately. +```@docs +mean(x::Symbol, by::Symbol, design::SimpleRandomSample) +total(x::Symbol, by::Symbol, design::SimpleRandomSample) +``` +```@docs +plot(design::AbstractSurveyDesign, x::Symbol, y::Symbol; kwargs...) +boxplot(design::AbstractSurveyDesign, x::Symbol, y::Symbol; kwargs...) +hist(design::AbstractSurveyDesign, var::Symbol, + bins::Union{Integer, AbstractVector} = freedman_diaconis(design, var); + normalization = :density, + kwargs... + ) +dim(design::AbstractSurveyDesign) +dimnames(design::AbstractSurveyDesign) +colnames(design::AbstractSurveyDesign) +``` diff --git a/docs/src/assets/boxplot.png b/docs/src/assets/boxplot.png index 0fc249da..b7f44543 100644 Binary files a/docs/src/assets/boxplot.png and b/docs/src/assets/boxplot.png differ diff --git a/docs/src/assets/hist.png b/docs/src/assets/hist.png index f8ef5881..966d3ea2 100644 Binary files a/docs/src/assets/hist.png and b/docs/src/assets/hist.png differ diff --git a/docs/src/assets/scatter.png b/docs/src/assets/scatter.png index 1c1304ec..fff128d7 100644 Binary files a/docs/src/assets/scatter.png and b/docs/src/assets/scatter.png differ diff --git a/docs/src/examples.md b/docs/src/examples.md deleted file mode 100644 index 1776ce2d..00000000 --- a/docs/src/examples.md +++ /dev/null @@ -1,7 +0,0 @@ -# Examples - -The following examples use the Academic Performance Index (API) dataset for Californian schools. - -```@docs -svyby(formula::Symbol, by, design::svydesign, func::Function, params = []) -``` diff --git a/docs/src/index.md b/docs/src/index.md index 72f72983..a099c73c 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -4,14 +4,59 @@ CurrentModule = Survey # Survey -This package is the Julia implementation of the [Survey package in R](https://cran.r-project.org/web/packages/survey/index.html) developed by [Professor Thomas Lumley](https://www.stat.auckland.ac.nz/people/tlum005). +This package is used to study complex survey data. It aims to be a fast alternative to the [Survey package in R](https://cran.r-project.org/web/packages/survey/index.html) developed by [Professor Thomas Lumley](https://www.stat.auckland.ac.nz/people/tlum005). -## The need for moving the code to Julia. +This package currently supports simple random sample and stratified sample. In future releases, it will support multistage sampling as well. -At [xKDR](https://xkdr.org/) we processed millions of records from household surveys using the survey package in R. This process took hours of computing time. By implementing the code in Julia, we are able to do the processing in seconds. In this package we have implemented the functions `svymean`, `svyquantile` and `svysum`. We have kept the syntax between the two packages similar so that we can easily move our existing code to the new language. +## Basic demo -Documentation for [Survey](https://github.com/Survey.jl). +The following demo uses the +[Academic Performance Index](https://r-survey.r-forge.r-project.org/survey/html/api.html) +(API) dataset for Californian schools. The data sets contain information for all schools +with at least 100 students and for various probability samples of the data. -```@autodocs -Modules = [Survey] +The API program has been discontinued at the end of 2018. Information is archived at +[https://www.cde.ca.gov/re/pr/api.asp](https://www.cde.ca.gov/re/pr/api.asp) + +Firstly, a survey design needs a dataset from which to gather information. + + +The sample datasets provided with the package can be loaded as `DataFrames` using the `load_data` function: + +```julia +julia> apisrs = load_data("apisrs"); +``` +`apisrs` is a simple random sample of the Academic Performance Index of Californian schools. + +Next, we can build a design. The design corresponding to a simple random sample is [`SimpleRandomSample`](@ref), which can be instantiated by calling the constructor: + +```julia +julia> srs = SimpleRandomSample(apisrs; weights = :pw) +SimpleRandomSample: +data: 200x42 DataFrame +weights: 31.0, 31.0, 31.0, ..., 31.0 +probs: 0.0323, 0.0323, 0.0323, ..., 0.0323 +fpc: 6194, 6194, 6194, ..., 6194 +popsize: 6194 +sampsize: 200 +sampfraction: 0.0323 +ignorefpc: false +``` + +With a `SimpleRandomSample` (as well as with any subtype of [`AbstractSurveyDesign`](@ref)) it is possible to calculate estimates of the mean, population total, etc., for a given variable, along with the corresponding standard errors. + +```julia +julia> mean(:api00, srs) +1×2 DataFrame + Row │ mean sem + │ Float64 Float64 +─────┼────────────────── + 1 │ 656.585 9.24972 + +julia> total(:api00, srs) +1×2 DataFrame + Row │ total se_total + │ Float64 Float64 +─────┼───────────────────── + 1 │ 4.06689e6 57292.8 ``` diff --git a/docs/src/performance.md b/docs/src/performance.md deleted file mode 100644 index 478b8f73..00000000 --- a/docs/src/performance.md +++ /dev/null @@ -1,74 +0,0 @@ -# Performance - -## Grouping by a single column -**R** - -```R -> library(survey) -> library(microbenchmark) -> data(api) -> dclus1 <- svydesign(id = ~dnum, weights = ~pw, data = apiclus1, fpc = ~fpc) -> microbenchmark(svyby(~api00, by = ~cname, design = dclus1, svymean, keep.var = FALSE), units = "us") -``` - -```R -Unit: microseconds - expr - svyby(~api00, by = ~cname, design = dclus1, svymean, keep.var = FALSE) - min lq mean median uq max neval - 9427.043 10587.81 11269.22 10938.55 11219.24 17620.25 100 -``` - -**Julia** -```julia -using Survey, BenchmarkTools -apiclus1 = load_data("apiclus1") -dclus1 = svydesign(id=:dnum, weights=:pw, data = apiclus1, fpc=:fpc) -@benchmark svyby(:api00, :cname, dclus1, svymean) -``` - -```julia -BenchmarkTools.Trial: 10000 samples with 1 evaluation. - Range (min … max): 43.567 μs … 5.905 ms ┊ GC (min … max): 0.00% … 90.27% - Time (median): 53.680 μs ┊ GC (median): 0.00% - Time (mean ± σ): 58.090 μs ± 125.671 μs ┊ GC (mean ± σ): 4.36% ± 2.00% -``` - -**The median time is about 198 times lower in Julia as compared to R.** - -## Grouping by two columns. - -**R** - -```R -> library(survey) -> library(microbenchmark) -> data(api) -> dclus1 <- svydesign(id = ~dnum, weights = ~pw, data = apiclus1, fpc = ~fpc) -> microbenchmark(svyby(~api00, by = ~cname+meals, design = dclus1, svymean, keep.var = FALSE), units = "us") -``` - -```R -Unit: microseconds - expr - svyby(~api00, by = ~cname + meals, design = dclus1, svymean, keep.var = FALSE) - min lq mean median uq max neval - 120823.6 131472.8 141797.3 134375.8 140818.3 263964.3 100 -``` - -**Julia** -```julia -using Survey, BenchmarkTools -apiclus1 = load_data("apiclus1") -dclus1 = svydesign(id=:dnum, weights=:pw, data = apiclus1, fpc=:fpc) -@benchmark svyby(:api00, [:cname, :meals], dclus1, svymean) -``` - -```julia -BenchmarkTools.Trial: 10000 samples with 1 evaluation. - Range (min … max): 64.591 μs … 6.559 ms ┊ GC (min … max): 0.00% … 77.46% - Time (median): 78.204 μs ┊ GC (median): 0.00% - Time (mean ± σ): 89.447 μs ± 235.344 μs ┊ GC (mean ± σ): 8.48% ± 3.19% -``` - - **The median time is about 1718 times lower in Julia as compared to R.** diff --git a/src/Survey.jl b/src/Survey.jl index 243343ff..f902bdc5 100644 --- a/src/Survey.jl +++ b/src/Survey.jl @@ -2,43 +2,36 @@ module Survey using DataFrames using Statistics +import Statistics: quantile using StatsBase +import StatsBase: mean,quantile using CSV -using GLM using LinearAlgebra using CairoMakie using AlgebraOfGraphics +using CategoricalArrays -include("svydesign.jl") -include("svyby.jl") -include("example.jl") -include("svyglm.jl") -include("svyhist.jl") -include("svyplot.jl") +include("SurveyDesign.jl") +include("mean.jl") +include("quantile.jl") +include("total.jl") +include("load_data.jl") +include("hist.jl") +include("plot.jl") include("dimnames.jl") -include("svyboxplot.jl") +include("boxplot.jl") +include("show.jl") -export svydesign, svyby, svyglm export load_data -export svymean, svytotal, svyquantile -export @formula -export svyhist, sturges, freedman_diaconis -export svyplot +export AbstractSurveyDesign, SimpleRandomSample, StratifiedSample +export SurveyDesign +export by +export ht_calc export dim, colnames, dimnames -export svyboxplot -export - #families - Normal, - Binomial, - Gamma, - Poisson, - #links - IdentityLink, - InverseLink, - LogLink, - LogitLink, - ProbitLink, - CauchitLink, - CloglogLink, - SqrtLink -end +export mean, total, quantile +export plot +export hist, sturges, freedman_diaconis +export boxplot +export ht_total, ht_mean + +end \ No newline at end of file diff --git a/src/SurveyDesign.jl b/src/SurveyDesign.jl new file mode 100644 index 00000000..13003458 --- /dev/null +++ b/src/SurveyDesign.jl @@ -0,0 +1,323 @@ +""" + AbstractSurveyDesign + +Supertype for every survey design type. + +!!! note + + The data passed to a survey constructor is modified. To avoid this pass a copy of the data + instead of the original. +""" +abstract type AbstractSurveyDesign end + +""" + SimpleRandomSample <: AbstractSurveyDesign + + +Survey design sampled by simple random sampling. + +# Arguments: +`data::AbstractDataFrame`: the survey dataset (!this gets modified by the constructor). +`sampsize::Union{Nothing,Symbol,<:Unsigned,Vector{<:Real}}=UInt(nrow(data))`: the survey sample size. +`popsize::Union{Nothing,Symbol,<:Unsigned,Vector{<:Real}}=nothing`: the (expected) survey population size. +`weights::Union{Nothing,Symbol,Vector{<:Real}}=nothing`: the sampling weights. +`probs::Union{Nothing,Symbol,Vector{<:Real}}=nothing: the sampling probabilities. +`ignorefpc::Bool=false`: choose to ignore finite population correction and assume all weights equal to 1.0 + +The precedence order of using `popsize`, `weights` and `probs` is `popsize` > `weights` > `probs`. +E.g. If `popsize` is given then it is assumed to be the ground truth over `weights` or `probs`. + +If `popsize` is not given `weights` or `probs` must be given. `popsize` is then calculated +using the weights and the sample size. + +```jldoctest +julia> apisrs = load_data("apisrs"); + +julia> srs = SimpleRandomSample(apisrs; popsize=:fpc) +SimpleRandomSample: +data: 200x42 DataFrame +weights: 31.0, 31.0, 31.0, ..., 31.0 +probs: 0.0323, 0.0323, 0.0323, ..., 0.0323 +fpc: 6194, 6194, 6194, ..., 6194 +popsize: 6194 +sampsize: 200 +sampfraction: 0.0323 +ignorefpc: false +``` +""" +struct SimpleRandomSample <: AbstractSurveyDesign + data::AbstractDataFrame + sampsize::Union{Unsigned,Nothing} + popsize::Union{Unsigned,Nothing} + sampfraction::Float64 + fpc::Float64 + ignorefpc::Bool + function SimpleRandomSample(data::AbstractDataFrame; + popsize::Union{Nothing,Symbol,Unsigned,Vector{<:Real}}=nothing, + sampsize::Union{Nothing,Symbol,Unsigned,Vector{<:Real}}=nrow(data) |> UInt, + weights::Union{Nothing,Symbol,Vector{<:Real}}=nothing, + probs::Union{Nothing,Symbol,Vector{<:Real}}=nothing, + ignorefpc::Bool=false + ) + # If any of weights or probs given as Symbol, + # find the corresponding column in `data` + if isa(weights, Symbol) + weights = data[!, weights] + end + if isa(probs, Symbol) + probs = data[!, probs] + end + # If weights/probs vector not numeric/real, ie. string column passed for weights, then raise error + if !isa(weights, Union{Nothing,Vector{<:Real}}) + error("weights should be Vector{<:Real}. You passed $(typeof(weights))") + elseif !isa(probs, Union{Nothing,Vector{<:Real}}) + error("sampling probabilities should be Vector{<:Real}. You passed $(typeof(probs))") + end + # If popsize given as Symbol or Vector, check all records equal + if isa(popsize, Symbol) + if !all(w -> w == first(data[!, popsize]), data[!, popsize]) + error("popsize must be same for all observations in Simple Random Sample") + end + popsize = first(data[!, popsize]) |> UInt + elseif isa(popsize, Vector{<:Real}) + if !all(w -> w == first(popsize), popsize) + error("popsize must be same for all observations in Simple Random Sample") + end + popsize = first(popsize) |> UInt + end + # If sampsize given as Symbol or Vector, check all records equal + if isa(sampsize, Symbol) + if !all(w -> w == first(data[!, sampsize]), data[!, sampsize]) + error("sampsize must be same for all observations in Simple Random Sample") + end + sampsize = first(data[!, sampsize]) |> UInt + elseif isa(sampsize, Vector{<:Real}) + if !all(w -> w == first(sampsize), sampsize) + error("sampsize must be same for all observations in Simple Random Sample") + end + sampsize = first(sampsize) |> UInt + end + # If both `weights` and `probs` given, then `weights` is assumed to be ground truth for probs. + if !isnothing(weights) && !isnothing(probs) + probs = 1 ./ weights + data[!, :probs] = probs + end + # popsize must be nothing or <:Unsigned by now + if isnothing(popsize) + # If popsize not given, fallback to weights, probs and sampsize to estimate `popsize` + @warn "popsize not given. using weights/probs and sampsize to estimate `popsize`" + # Check that all weights (or probs if weights not given) are equal, as SRS is by definition equi-weighted + if typeof(weights) <: Vector{<:Real} + if !all(w -> w == first(weights), weights) + error("all frequency weights must be equal for Simple Random Sample") + end + elseif typeof(probs) <: Vector{<:Real} + if !all(p -> p == first(probs), probs) + error("all probability weights must be equal for Simple Random Sample") + end + weights = 1 ./ probs + else + error("either weights or probs must be given if `popsize` not given") + end + # Estimate population size + popsize = round(sampsize * first(weights)) |> UInt + elseif typeof(popsize) <: Unsigned + weights = fill(popsize / sampsize, nrow(data)) # If popsize is given, weights vector is made concordant with popsize and sampsize, regardless of given weights argument + probs = 1 ./ weights + else + error("something went wrong, please check validity of inputs.") + end + # If sampsize greater than popsize than illogical arguments specified. + if sampsize > popsize + error("population size was estimated to be less than given sampsize. Please check input arguments.") + end + # If ignorefpc then set weights to 1 ?? + # TODO: This works under some cases, but should find better way to process ignoring fpc + if ignorefpc + @warn "assuming all weights are equal to 1.0" + weights = ones(nrow(data)) + probs = 1 ./ weights + end + # sum of weights must equal to `popsize` for SRS + if !isnothing(weights) && !(isapprox(sum(weights), popsize; atol=1e-4)) + if ignorefpc && !(isapprox(sum(weights), sampsize; atol=1e-4)) # Change if ignorefpc functionality changes + error("sum of sampling weights should be equal to `sampsize` for `SimpleRandomSample` with `ignorefpc`") + elseif !ignorefpc + error("sum of sampling weights must be equal to `popsize` for `SimpleRandomSample`") + end + end + # sum of probs must equal popsize for SRS + if !isnothing(probs) && !(isapprox(sum(1 ./ probs), popsize; atol=1e-4)) + if ignorefpc && !(isapprox(sum(1 ./ probs), sampsize; atol=1e-4)) # Change if ignorefpc functionality changes + error("sum of inverse sampling probabilities should be equal to `sampsize` for `SimpleRandomSample` with `ignorefpc`") + elseif !ignorefpc + error("sum of inverse of sampling probabilities must be equal to `popsize` for Simple Random Sample") + end + end + ## Set remaining parts of data structure + # set sampling fraction + sampfraction = sampsize / popsize + # set fpc + fpc = ignorefpc ? 1 : 1 - (sampsize / popsize) + # add columns for frequency and probability weights to `data` + data[!, :weights] = weights + if isnothing(probs) + probs = 1 ./ data[!, :weights] + end + data[!, :probs] = probs + # Initialise the structure + new(data, sampsize, popsize, sampfraction, fpc, ignorefpc) + end +end + +""" + StratifiedSample <: AbstractSurveyDesign + +Survey design sampled by stratification. + +`strata` must be specified as a Symbol name of a column in `data`. + +# Arguments: +`data::AbstractDataFrame`: the survey dataset (!this gets modified by the constructor). +`strata::Symbol`: the stratification variable - must be given as a column in `data`. +`sampsize::Union{Nothing,Symbol,<:Unsigned,Vector{<:Real}}=UInt(nrow(data))`: the survey sample size. +`popsize::Union{Nothing,Symbol,<:Unsigned,Vector{<:Real}}=nothing`: the (expected) survey population size. +`weights::Union{Nothing,Symbol,Vector{<:Real}}=nothing`: the sampling weights. +`probs::Union{Nothing,Symbol,Vector{<:Real}}=nothing: the sampling probabilities. +`ignorefpc::Bool=false`: choose to ignore finite population correction and assume all weights equal to 1.0 + +The `popsize`, `weights` and `probs` parameters follow the same rules as for [`SimpleRandomSample`](@ref). + +```jldoctest +julia> apistrat = load_data("apistrat"); + +julia> dstrat = StratifiedSample(apistrat, :stype; popsize=:fpc) +StratifiedSample: +data: 200x45 DataFrame +strata: stype +weights: 44.2, 44.2, 44.2, ..., 15.1 +probs: 0.0226, 0.0226, 0.0226, ..., 0.0662 +fpc: 0.977, 0.977, 0.977, ..., 0.934 +popsize: 4421, 4421, 4421, ..., 755 +sampsize: 100, 100, 100, ..., 50 +sampfraction: 0.0226, 0.0226, 0.0226, ..., 0.0662 +ignorefpc: false +``` +""" +struct StratifiedSample <: AbstractSurveyDesign + data::AbstractDataFrame + strata::Symbol + ignorefpc::Bool + function StratifiedSample(data::AbstractDataFrame, strata::Symbol; + popsize::Union{Nothing,Symbol}=nothing, + sampsize::Union{Nothing,Symbol}=nothing, + weights::Union{Nothing,Symbol,Vector{<:Real}}=nothing, + probs::Union{Nothing,Symbol,Vector{<:Real}}=nothing, + ignorefpc::Bool=false + ) + # Store the iterator over each strata, as used multiple times + data_groupedby_strata = groupby(data, strata) + # If any of weights or probs given as Symbol, find the corresponding column in `data` + if isa(weights, Symbol) + for each_strata in keys(data_groupedby_strata) + if !all(w -> w == first(data_groupedby_strata[each_strata][!, weights]), data_groupedby_strata[each_strata][!, weights]) + error("sampling weights within each strata must be equal in StratifiedSample") + end + end + # original_weights_colname = copy(weights) + weights = data[!, weights] # If all good with weights column, then store it as Vector + end + if isa(probs, Symbol) + for each_strata in keys(data_groupedby_strata) + if !all(p -> p == first(data_groupedby_strata[each_strata][!, probs]), data_groupedby_strata[each_strata][!, probs]) + error("sampling probabilities within each strata must be equal in StratifiedSample") + end + end + # original_probs_colname = copy(probs) + probs = data[!, probs] # If all good with probs column, then store it as Vector + end + # If weights/probs vector not numeric/real, ie. string column passed for weights, then raise error + if !isa(weights, Union{Nothing,Vector{<:Real}}) + error("weights should be Vector{<:Real}. You passed $(typeof(weights))") + elseif !isa(probs, Union{Nothing,Vector{<:Real}}) + error("sampling probabilities should be Vector{<:Real}. You passed $(typeof(probs))") + end + # If popsize given as Symbol or Vector, check all records equal in each strata + if isa(popsize, Symbol) + for each_strata in keys(data_groupedby_strata) + if !all(w -> w == first(data_groupedby_strata[each_strata][!, popsize]), data_groupedby_strata[each_strata][!, popsize]) + error("popsize must be same for all observations within each strata in StratifiedSample") + end + end + # original_popsize_colname = copy(popsize) + popsize = data[!, popsize] + end + # If sampsize given as Symbol or Vector, check all records equal + if isa(sampsize, Symbol) + if isnothing(popsize) && isnothing(weights) && isnothing(probs) + error("if sampsize given, and popsize not given, then weights or probs must given to calculate popsize") + end + for each_strata in keys(data_groupedby_strata) + if !all(w -> w == first(data_groupedby_strata[each_strata][!, sampsize]), data_groupedby_strata[each_strata][!, sampsize]) + error("sampsize must be same for all observations within each strata in StratifiedSample") + end + end + # original_sampsize_colname = copy(sampsize) + sampsize = data[!, sampsize] + # If sampsize column not provided in constructor call, set it as nrow of strata + elseif isnothing(sampsize) + sampsize = transform(groupby(data, strata), nrow => :counts).counts + end + # If both `weights` and `probs` given, then `weights` is assumed to be ground truth for probs. + if !isnothing(weights) && !isnothing(probs) + probs = 1 ./ weights + data[!, :probs] = probs + end + # `popsize` is either nothing or a Vector{<:Real} by now + if isnothing(popsize) + # If popsize not given, fallback to weights, probs and sampsize to estimate `popsize` + @warn "popsize not given. using weights/probs and sampsize to estimate `popsize` for StratifiedSample" + # Check that all weights (or probs if weights not given) are equal, as SRS is by definition equi-weighted + if typeof(probs) <: Vector{<:Real} + weights = 1 ./ probs + elseif !(typeof(weights) <: Vector{<:Real}) + error("either weights or probs must be given if `popsize` not given") + end + # Estimate population size + popsize = sampsize .* weights + elseif typeof(popsize) <: Vector{<:Real} # Still need to check if the provided Column is of <:Real + # If popsize is given, weights and probs made concordant with popsize and sampsize, regardless of supplied arguments + weights = popsize ./ sampsize + probs = 1 ./ weights + else + error("something went wrong. Please check validity of inputs.") + end + # If sampsize greater than popsize than illogical arguments specified. + if any(sampsize .> popsize) + error("population sizes were estimated to be less than sampsize. please check input arguments.") + end + # If ignorefpc then set weights to 1 ?? + # TODO: This works under some cases, but should find better way to process ignoring fpc + if ignorefpc + @warn "assuming all weights are equal to 1.0" + weights = ones(nrow(data)) + probs = 1 ./ weights + end + ## Set remaining parts of data structure + # set sampling fraction + sampfraction = sampsize ./ popsize + # set fpc + fpc = ignorefpc ? fill(1, size(data, 1)) : 1 .- (sampsize ./ popsize) + # add columns for frequency and probability weights to `data` + data[!, :weights] = weights + if isnothing(probs) + probs = 1 ./ data[!, :weights] + end + data[!, :probs] = probs + data[!, :sampsize] = sampsize + data[!, :popsize] = popsize + data[!, :fpc] = fpc + data[!, :sampfraction] = sampfraction + new(data, strata, ignorefpc) + end +end \ No newline at end of file diff --git a/src/boxplot.jl b/src/boxplot.jl new file mode 100644 index 00000000..46f6958a --- /dev/null +++ b/src/boxplot.jl @@ -0,0 +1,25 @@ +""" + boxplot(design, x, y; kwargs...) + +Box plot of survey design variable `y` grouped by column `x`. + +Weights can be specified by a `Symbol` using the keyword argument `weights`. + +The keyword arguments are all the arguments that can be passed to `mapping` in +[AlgebraOfGraphics](https://docs.juliahub.com/AlgebraOfGraphics/CHIaw/0.4.7/). + +```@example boxplot +apisrs = load_data("apisrs"); +srs = SimpleRandomSample(apisrs; weights = :pw); +bp = boxplot(srs, :stype, :enroll; weights = :pw) +save("boxplot.png", bp); nothing # hide +``` + +![](assets/boxplot.png) +""" +function boxplot(design::AbstractSurveyDesign, x::Symbol, y::Symbol; kwargs...) + map = mapping(x, y; kwargs...) + data = AlgebraOfGraphics.data(design.data) + + data * visual(BoxPlot) * map |> draw +end diff --git a/src/dimnames.jl b/src/dimnames.jl index 09746613..91bd473e 100644 --- a/src/dimnames.jl +++ b/src/dimnames.jl @@ -1,33 +1,31 @@ """ dim(design) -Get the dimensions of a survey design. -```jldoctest -julia> using Survey +Get the dimensions of a `SurveyDesign`. -julia> apistrat = load_data("apistrat"); +```jldoctest +julia> apisrs = load_data("apisrs"); -julia> dstrat = svydesign(data = apistrat, id = :1, strata = :stype, weights = :pw, fpc = :fpc); +julia> srs = SimpleRandomSample(apisrs; popsize =:fpc); -julia> dim(dstrat) -(200, 45) +julia> dim(srs) +(200, 42) ``` """ -dim(design::svydesign) = size(design.variables) +dim(design::AbstractSurveyDesign) = size(design.data) """ colnames(design) -Get the column names of a survey design. -```jldoctest -julia> using Survey +Get the column names of a `SurveyDesign`. -julia> apistrat = load_data("apistrat"); +```jldoctest +julia> apisrs = load_data("apisrs"); -julia> dstrat = svydesign(data = apistrat, id = :1, strata = :stype, weights = :pw, fpc = :fpc); +julia> srs = SimpleRandomSample(apisrs; popsize=:fpc); -julia> colnames(dstrat) -45-element Vector{String}: +julia> colnames(srs) +42-element Vector{String}: "Column1" "cds" "stype" @@ -39,34 +37,33 @@ julia> colnames(dstrat) "cname" "cnum" ⋮ + "avg.ed" + "full" + "emer" "enroll" "api.stu" "pw" "fpc" - "probs" "weights" - "popsize" - "sampsize" - "strata" + "probs" ``` """ -colnames(design::svydesign) = names(design.variables) +colnames(design::AbstractSurveyDesign) = names(design.data) """ dimnames(design) -Get the names of the rows and columns of a survey design. -```jldoctest -julia> using Survey +Get the names of the rows and columns of a `SurveyDesign`. -julia> apistrat = load_data("apistrat"); +```jldoctest +julia> apisrs = load_data("apisrs"); -julia> dstrat = svydesign(data = apistrat, id = :1, strata = :stype, weights = :pw, fpc = :fpc); +julia> srs = SimpleRandomSample(apisrs;popsize=:fpc); -julia> dimnames(dstrat) +julia> dimnames(srs) 2-element Vector{Vector{String}}: ["1", "2", "3", "4", "5", "6", "7", "8", "9", "10" … "191", "192", "193", "194", "195", "196", "197", "198", "199", "200"] - ["Column1", "cds", "stype", "name", "sname", "snum", "dname", "dnum", "cname", "cnum" … "emer", "enroll", "api.stu", "pw", "fpc", "probs", "weights", "popsize", "sampsize", "strata"] + ["Column1", "cds", "stype", "name", "sname", "snum", "dname", "dnum", "cname", "cnum" … "grad.sch", "avg.ed", "full", "emer", "enroll", "api.stu", "pw", "fpc", "weights", "probs"] ``` """ -dimnames(design::svydesign) = [string.(1:size(design.variables, 1)), names(design.variables)] +dimnames(design::AbstractSurveyDesign) = [string.(1:size(design.data, 1)), names(design.data)] diff --git a/src/example.jl b/src/example.jl deleted file mode 100644 index 0c49e0d4..00000000 --- a/src/example.jl +++ /dev/null @@ -1,26 +0,0 @@ -const PKG_DIR = joinpath(pathof(Survey), "..", "..") |> normpath -asset_path(args...) = joinpath(PKG_DIR, "assets", args...) - -""" -The Academic Performance Index is computed for all California schools based on standardised testing of students. -The data sets contain information for all schools with at least 100 students and for various probability samples of the data. - -Use `load_data(name)` to load API data, with -`name ∈ ["apiclus1", "apiclus2", "apipop", "apistrat", "apisrs"]` -being the name of the dataset. - -```julia -df = load_data("apiclus1") -``` - -Details about the columns of the dataset can be found here: https://r-survey.r-forge.r-project.org/survey/html/api.html - -The API program has been discontinued at the end of 2018. Information is archived at https: -//www.cde.ca.gov/re/pr/api.asp -""" -function load_data(name) - name = name * ".csv" - @assert name ∈ readdir(asset_path()) - - CSV.read(asset_path(name), DataFrame, missingstring="NA") -end diff --git a/src/hist.jl b/src/hist.jl new file mode 100644 index 00000000..c140e59e --- /dev/null +++ b/src/hist.jl @@ -0,0 +1,85 @@ +sturges(n::Integer) = ceil(Int, log2(n)) + 1 +sturges(vec::AbstractVector) = ceil(Int, log2(length(vec))) + 1 +sturges(df::DataFrame, var::Symbol) = ceil(Int, log2(size(df[!, var], 1))) + 1 + +""" + sturges(design::SurveyDesign, var::Symbol) + +Calculate the number of bins to use in a histogram using the Sturges rule. + +# Examples +```jldoctest +julia> apisrs = load_data("apisrs"); + +julia> srs = SimpleRandomSample(apisrs;popsize=:fpc); + +julia> sturges(srs, :enroll) +9 +``` +""" +sturges(design::AbstractSurveyDesign, var::Symbol) = sturges(design.data, var) + +freedman_diaconis(v::AbstractVector) = round(Int, length(v)^(1 / 3) * (maximum(v) - minimum(v)) / (2 * iqr(v))) +freedman_diaconis(df::DataFrame, var::Symbol) = freedman_diaconis(df[!, var]) + +""" + freedman_diaconis(design::SurveyDesign, var::Symbol) + +Calculate the number of bins to use in a histogram using the Freedman-Diaconis rule. + +# Examples +```jldoctest +julia> apisrs = load_data("apisrs"); + +julia> srs = SimpleRandomSample(apisrs;popsize=:fpc); + +julia> freedman_diaconis(srs, :enroll) +18 +``` +""" +freedman_diaconis(design::AbstractSurveyDesign, var::Symbol) = freedman_diaconis(design.data[!, var]) + +""" + hist(design, var, bins = freedman_diaconis; normalization = :density, kwargs...) + +Histogram plot of a survey design variable given by `var`. + +`bins` can be an `Integer` specifying the number of equal-width bins, +an `AbstractVector` specifying the bins intervals, or a `Function` specifying +the function used for calculating the number of bins. The possible functions +are `sturges` and `freedman_diaconis`. + +The normalization can be set to `:none`, `:density`, `:probability` or `:pdf`. +See [AlgebraOfGraphics.histogram](https://docs.juliahub.com/AlgebraOfGraphics/CHIaw/0.4.9/generated/datatransformations/#AlgebraOfGraphics.histogram) +for more information. + +For the complete argument list see [Makie.hist](https://makie.juliaplots.org/stable/examples/plotting_functions/hist/). + +!!! note + + The `weights` argument should be a `Symbol` specifying a design variable. + +```@example histogram +apisrs = load_data("apisrs"); +srs = SimpleRandomSample(apisrs;popsize=:fpc); +h = hist(srs, :enroll) +save("hist.png", h); nothing # hide +``` + +![](assets/hist.png) +""" +function hist(design::AbstractSurveyDesign, var::Symbol, + bins::Union{Integer, AbstractVector} = freedman_diaconis(design, var); + normalization = :density, + kwargs... + ) + hist = histogram(bins = bins, normalization = normalization, kwargs...) + data(design.data) * mapping(var, weights = :weights) * hist |> draw +end + +function hist(design::AbstractSurveyDesign, var::Symbol, + bins::Function; + kwargs... + ) + hist(design, var, bins(design, var); kwargs...) +end diff --git a/src/load_data.jl b/src/load_data.jl new file mode 100644 index 00000000..66344051 --- /dev/null +++ b/src/load_data.jl @@ -0,0 +1,41 @@ +const PKG_DIR = joinpath(pathof(Survey), "..", "..") |> normpath +asset_path(args...) = joinpath(PKG_DIR, "assets", args...) + +""" + load_data(name) + +Load a sample dataset as a `DataFrame`. + +All available datasets can be found [here](https://github.com/xKDR/Survey.jl/tree/main/assets). + +```jldoctest +julia> apisrs = load_data("apisrs") +200×40 DataFrame + Row │ Column1 cds stype name sname ⋯ + │ Int64 Int64 String1 String15 String ⋯ +─────┼────────────────────────────────────────────────────────────────────────── + 1 │ 1039 15739081534155 H McFarland High McFarland High ⋯ + 2 │ 1124 19642126066716 E Stowers (Cecil Stowers (Cecil B.) E + 3 │ 2868 30664493030640 H Brea-Olinda Hig Brea-Olinda High + 4 │ 1273 19644516012744 E Alameda Element Alameda Elementary + 5 │ 4926 40688096043293 E Sunnyside Eleme Sunnyside Elementary ⋯ + 6 │ 2463 19734456014278 E Los Molinos Ele Los Molinos Elementa + 7 │ 2031 19647336058200 M Northridge Midd Northridge Middle + 8 │ 1736 19647336017271 E Glassell Park E Glassell Park Elemen + ⋮ │ ⋮ ⋮ ⋮ ⋮ ⋮ ⋱ + 194 │ 4880 39686766042782 E Tyler Skills El Tyler Skills Element ⋯ + 195 │ 993 15636851531987 H Desert Junior/S Desert Junior/Senior + 196 │ 969 15635291534775 H North High North High + 197 │ 1752 19647336017446 E Hammel Street E Hammel Street Elemen + 198 │ 4480 37683386039143 E Audubon Element Audubon Elementary ⋯ + 199 │ 4062 36678196036222 E Edison Elementa Edison Elementary + 200 │ 2683 24657716025621 E Franklin Elemen Franklin Elementary + 36 columns and 185 rows omitted +``` +""" +function load_data(name) + name = name * ".csv" + @assert name ∈ readdir(asset_path()) + + CSV.read(asset_path(name), DataFrame, missingstring="NA") +end diff --git a/src/mean.jl b/src/mean.jl new file mode 100644 index 00000000..95869857 --- /dev/null +++ b/src/mean.jl @@ -0,0 +1,167 @@ +""" + mean(x, design) + +Estimate the population mean of a variable of a simple random sample, and the corresponding standard error. + +The calculations were done according to the book [Sampling Techniques](https://www.academia.edu/29684662/Cochran_1977_Sampling_Techniques_Third_Edition) +by William Cochran. + +```jldoctest +julia> apisrs = load_data("apisrs"); + +julia> srs = SimpleRandomSample(apisrs;popsize=:fpc); + +julia> mean(:enroll, srs) +1×2 DataFrame + Row │ mean SE + │ Float64 Float64 +─────┼────────────────── + 1 │ 584.61 27.3684 + +julia> mean([:api00, :api99], srs) +2×3 DataFrame + Row │ names mean SE + │ String Float64 Float64 +─────┼────────────────────────── + 1 │ api00 656.585 9.24972 + 2 │ api99 624.685 9.5003 + +julia> strat = load_data("apistrat"); + +julia> dstrat = StratifiedSample(strat, :stype; popsize = :fpc); + +julia> mean(:api00, dstrat) +1×2 DataFrame + Row │ mean SE + │ Float64 Float64 +─────┼────────────────── + 1 │ 662.287 9.40894 +``` +""" +function mean(x::Symbol, design::SimpleRandomSample) + function se(x::Symbol, design::SimpleRandomSample) + variance = design.fpc * Statistics.var(design.data[!, x]) / design.sampsize + return sqrt(variance) + end + if isa(design.data[!, x], CategoricalArray) + gdf = groupby(design.data, x) + p = combine(gdf, nrow => :counts) + p.mean = p.counts ./ sum(p.counts) + # variance of proportion + p.var = design.fpc .* p.mean .* (1 .- p.mean) ./ (design.sampsize - 1) + p.se = sqrt.(p.var) + return select(p, Not([:counts, :var])) + end + return DataFrame(mean=mean(design.data[!, x]), SE=se(x, design)) +end + +function mean(x::Vector{Symbol}, design::SimpleRandomSample) + df = reduce(vcat, [mean(i, design) for i in x]) + insertcols!(df, 1, :names => String.(x)) + return df +end + +function mean(x::Symbol, design::StratifiedSample) + if x == design.strata + gdf = groupby(design.data, x) + p = combine(gdf, :weights => sum => :Nₕ) + p.Wₕ = p.Nₕ ./ sum(p.Nₕ) + p = select!(p, Not(:Nₕ)) + return p + elseif isa(design.data[!, x], CategoricalArray) + gdf = groupby(design.data, x) + p = combine(gdf, nrow => :counts) + p.proportion = p.counts ./ sum(p.counts) + # variance of proportion + p.var = design.fpc .* p.proportion .* (1 .- p.proportion) ./ (design.sampsize - 1) + p.SE = sqrt.(p.var) + return p + end + gdf = groupby(design.data, design.strata) + ȳₕ = combine(gdf, x => mean => :mean).mean + Nₕ = combine(gdf, :weights => sum => :Nₕ).Nₕ + nₕ = combine(gdf, nrow => :nₕ).nₕ + fₕ = nₕ ./ Nₕ + Wₕ = Nₕ ./ sum(Nₕ) + Ȳ̂ = sum(Wₕ .* ȳₕ) + s²ₕ = combine(gdf, x => var => :s²h).s²h + V̂Ȳ̂ = sum((Wₕ .^ 2) .* (1 .- fₕ) .* s²ₕ ./ nₕ) + SE = sqrt(V̂Ȳ̂) + return DataFrame(mean=Ȳ̂, SE=SE) +end + +""" + mean(x, by, design) + +Estimate the subpopulation mean of a variable `x`. + +The calculations were done according to the book [Model-Assisted Survey Sampling](https://link.springer.com/book/9780387406206) +by Carl-Erik Sarndal, Bengt Swensson, Jan Wretman, section 3.3 and Chap 10. Assumes popsize is known and subpopulation size is unknown. + +```jldoctest +julia> using Survey; + +julia> srs = load_data("apisrs"); + +julia> srs = SimpleRandomSample(srs; popsize = :fpc); + +julia> mean(:api00, :cname, srs) |> first +DataFrameRow + Row │ cname mean SE + │ String15 Float64 Float64 +─────┼──────────────────────────── + 1 │ Kern 573.6 42.8026 + +julia> strat = load_data("apistrat"); + +julia> dstrat = StratifiedSample(strat, :stype; popsize = :fpc); + +julia> mean(:api00, :cname, dstrat) |> first +DataFrameRow + Row │ cname mean SE + │ String15 Float64 Float64 +─────┼─────────────────────────────── + 1 │ Los Angeles 633.511 21.3912 +``` +""" +function mean(x::Symbol, by::Symbol, design::SimpleRandomSample) + function domain_mean(x::AbstractVector, design::SimpleRandomSample, weights) + function se(x::AbstractVector, design::SimpleRandomSample) + nd = length(x) # domain size + n = design.sampsize + fpc = design.fpc + variance = (nd / n)^(-2) / n * fpc * ((nd - 1) / (n - 1)) * var(x) + return sqrt(variance) + end + return DataFrame(mean=Statistics.mean(x), SE=se(x, design)) + end + gdf = groupby(design.data, by) + combine(gdf, [x, :weights] => ((a, b) -> domain_mean(a, design, b)) => AsTable) +end + +function mean(x::Symbol, by::Symbol, design::StratifiedSample) + function domain_mean(x::AbstractVector, popsize::AbstractVector, sampsize::AbstractVector, sampfraction::AbstractVector, strata::AbstractVector) + df = DataFrame(x=x, popsize=popsize, sampsize=sampsize, sampfraction=sampfraction, strata=strata) + function calculate_components(x, popsize, sampsize, sampfraction) + return DataFrame( + nsdh = length(x), + nsh = length(x), + substrata_domain_totals = sum(x), + ȳsdh = mean(x), + Nh = first(popsize), + nh = first(sampsize), + fh = first(sampfraction), + sigma_ȳsh_squares = sum((x .- mean(x)).^2) + ) + end + components = combine(groupby(df, :strata), [:x, :popsize, :sampsize, :sampfraction] => calculate_components => AsTable) + domain_mean = sum(components.Nh .* components.substrata_domain_totals ./ components.nh) / sum(components.Nh .* components.nsdh ./ components.nh) + pdh = components.nsdh ./ components.nh + N̂d = sum(components.Nh .* pdh) + domain_var = sum(components.Nh .^ 2 .* (1 .- components.fh) .* (components.sigma_ȳsh_squares .+ (components.nsdh .* (1 .- pdh) .* (components.ȳsdh .- domain_mean) .^ 2)) ./ (components.nh .* (components.nh .- 1))) ./ N̂d .^ 2 + domain_mean_se = sqrt(domain_var) + return DataFrame(mean=domain_mean, SE=domain_mean_se) + end + gdf_domain = groupby(design.data, by) + combine(gdf_domain, [x, :popsize,:sampsize,:sampfraction, design.strata] => domain_mean => AsTable) +end diff --git a/src/plot.jl b/src/plot.jl new file mode 100644 index 00000000..cb9792d1 --- /dev/null +++ b/src/plot.jl @@ -0,0 +1,20 @@ +""" + plot(design, x, y; kwargs...) + +Scatter plot of survey design variables `x` and `y`. + +The plot takes into account the frequency weights specified by the user +in the design. + +```@example plot +apisrs = load_data("apisrs"); +srs = SimpleRandomSample(apisrs; weights = :pw); +s = plot(srs, :api99, :api00) +save("scatter.png", s); nothing # hide +``` + +![](assets/scatter.png) +""" +function plot(design::AbstractSurveyDesign, x::Symbol, y::Symbol; kwargs...) + data(design.data) * mapping(x, y, markersize = :weights) * visual(Scatter, marker = '￮') |> draw +end diff --git a/src/quantile.jl b/src/quantile.jl new file mode 100644 index 00000000..8ee6000c --- /dev/null +++ b/src/quantile.jl @@ -0,0 +1,53 @@ +""" + quantile(var, design, p; kwargs...) +Estimate quantiles for a complex survey. + +Hyndman and Fan compiled a taxonomy of nine algorithms to estimate quantiles. These are implemented in Statistics.quantile, which this function calls. +The Julia, R and Python-numpy use the same defaults + +# References: +- Hyndman, R.J and Fan, Y. (1996) ["Sample Quantiles in Statistical Packages"](https://www.amherst.edu/media/view/129116/original/Sample+Quantiles.pdf), The American Statistician, Vol. 50, No. 4, pp. 361-365. +- [Quantiles](https://en.m.wikipedia.org/wiki/Quantile) on wikipedia +- [Complex Surveys: a guide to analysis using R](https://r-survey.r-forge.r-project.org/svybook/), Section 2.4.1 and Appendix C.4. + +```jldoctest +julia> apisrs = load_data("apisrs"); + +julia> srs = SimpleRandomSample(apisrs;popsize=:fpc); + +julia> quantile(:api00,srs,0.5) +1×2 DataFrame + Row │ probability quantile + │ Float64 Float64 +─────┼─────────────────────── + 1 │ 0.5 659.0 + +julia> quantile(:enroll,srs,[0.1,0.2,0.5,0.75,0.95]) +5×2 DataFrame + Row │ probability quantile + │ Float64 Float64 +─────┼─────────────────────── + 1 │ 0.1 245.5 + 2 │ 0.2 317.6 + 3 │ 0.5 453.0 + 4 │ 0.75 668.5 + 5 │ 0.95 1473.1 +``` +""" +function quantile(var::Symbol, design::SimpleRandomSample, p::Union{<:Real,Vector{<:Real}}; + alpha::Float64=0.05, ci::Bool=false, se::Bool=false, qrule="hf7",kwargs...) + v = design.data[!, var] + probs = design.data[!, :probs] + df = DataFrame(probability=p, quantile=Statistics.quantile(v, ProbabilityWeights(probs), p)) + # TODO: Add CI and SE of the quantile + return df +end + +function quantile(var::Symbol, design::StratifiedSample, p::Union{<:Real,Vector{<:Real}}; + alpha::Float64=0.05, ci::Bool=false, se::Bool=false, qrule="hf7",kwargs...) + v = design.data[!, var] + probs = design.data[!, :probs] + df = DataFrame(probability=p, quantile=Statistics.quantile(v, ProbabilityWeights(probs), p)) # Not sure which quantile defintion this returns + # TODO: Add CI and SE of the quantile + return df +end \ No newline at end of file diff --git a/src/show.jl b/src/show.jl new file mode 100644 index 00000000..0daffe24 --- /dev/null +++ b/src/show.jl @@ -0,0 +1,49 @@ +""" +Helper function that transforms a given `Number` or `Vector` into a short-form string. +""" +function makeshort(x) + if isa(x[1], Float64) + x = round.(x, sigdigits=3) + end + # print short vectors or single values as they are, compress otherwise + x = length(x) < 3 ? join(x, ", ") : join(x[1:3], ", ") * ", ..., " * string(last(x)) +end + +""" +Print information in the form: + **name:** content[\n] +""" +function printinfo(io::IO, name::String, content::String; newline::Bool=true) + printstyled(io, name, ": "; bold=true) + newline ? println(io, content) : print(io, content) +end + +"Print information about a survey design." +function Base.show(io::IO, ::MIME"text/plain", design::AbstractSurveyDesign) + type = typeof(design) + printstyled(io, "$type:\n"; bold=true) + printstyled(io, "data: "; bold=true) + println(io, size(design.data, 1), "x", size(design.data, 2), " DataFrame") + printinfo(io, "weights", makeshort(design.data.weights)) + printinfo(io, "probs", makeshort(design.data.probs)) + printinfo(io, "fpc", makeshort(design.data.fpc)) + printinfo(io, "popsize", makeshort(design.popsize)) + printinfo(io, "sampsize", makeshort(design.sampsize)) + printinfo(io, "sampfraction", makeshort(design.sampfraction)) + printinfo(io, "ignorefpc", string(design.ignorefpc); newline=false) +end + +function Base.show(io::IO, ::MIME"text/plain", design::StratifiedSample) + type = typeof(design) + printstyled(io, "$type:\n"; bold=true) + printstyled(io, "data: "; bold=true) + println(io, size(design.data, 1), "x", size(design.data, 2), " DataFrame") + printinfo(io, "strata", string(design.strata); newline=true) + printinfo(io, "weights", makeshort(design.data.weights)) + printinfo(io, "probs", makeshort(design.data.probs)) + printinfo(io, "fpc", makeshort(design.data.fpc)) + printinfo(io, "popsize", makeshort(design.data.popsize)) + printinfo(io, "sampsize", makeshort(design.data.sampsize)) + printinfo(io, "sampfraction", makeshort(design.data.sampfraction)) + printinfo(io, "ignorefpc", string(design.ignorefpc); newline=false) +end \ No newline at end of file diff --git a/src/svyboxplot.jl b/src/svyboxplot.jl deleted file mode 100644 index 9e999144..00000000 --- a/src/svyboxplot.jl +++ /dev/null @@ -1,29 +0,0 @@ -""" -``` - svyboxplot(design, x, y; kwargs...) -``` -Box plot of survey design variable `y` grouped by column `x`. - -Weights can be specified by a Symbol using the keyword argument `weights`. - -The keyword arguments are all the arguments that can be passed to `mapping` in -[AlgebraOfGraphics](https://docs.juliahub.com/AlgebraOfGraphics/CHIaw/0.4.7/). - -```@example svyboxplot -julia> using survey - -julia> data(api); - -julia> dstrat = svydesign(data = apistrat, id = :1, strata = :stype, weights = :pw, fpc = :fpc); - -julia> bp = svyboxplot(dstrat, :stype, :enroll; weights = :pw) -``` - -![](./assets/boxplot.png) -""" -function svyboxplot(design::svydesign, x::Symbol, y::Symbol; kwargs...) - map = mapping(x, y; kwargs...) - data = AlgebraOfGraphics.data(design.variables) - - data * visual(BoxPlot) * map |> draw -end diff --git a/src/svyby.jl b/src/svyby.jl deleted file mode 100644 index 33c683ad..00000000 --- a/src/svyby.jl +++ /dev/null @@ -1,85 +0,0 @@ -function svymean(y, design::svydesign) - # popsize correction isn't implemented yet - ss = maximum(design.variables.sampsize) - w = design.variables.probs - x = design.variables[!, y] - function SE(x, w, ss) - f = sqrt(1 - 1 / length(x) + 1 / ss) - var = sum(w .* (x .- sum(w .* x) / sum(w)).^2) / sum(w) - sd = sqrt(var / (length(x) - 1)) - return f * sd - end - - return DataFrame(mean = mean(x, weights(w)), SE = SE(x, w, ss)) -end - -function svymean(x, w, popsize, sampsize) - # popsize correction isn't implemented yet - ss = maximum(sampsize) - function SE(x, w, ss) - f = sqrt(1 - 1/length(x) + 1/ss) - var = sum(w.*(x.- sum(w.*x)/sum(w)).^2)/sum(w) - sd = sqrt(var / (length(x) -1)) - return f * sd - end - return DataFrame(mean = mean(x, weights(w)), SE = SE(x, w, ss)) -end - -function svyquantile(var, design::svydesign, q) - x = design.variables[!, var] - w = design.variables.probs - df = DataFrame(tmp = quantile(Float32.(x), weights(w), q)) - rename!(df,:tmp => Symbol(string(q) .* "th percentile")) - return df -end - -function svyquantile(x, w, popsize, sampsize, q) - df = DataFrame(tmp = quantile(Float32.(x), weights(w), q)) - rename!(df,:tmp => Symbol(string(q) .* "th percentile")) - return df -end - -function svytotal(var, design::svydesign) - w = design.variables.probs; - x = design.variables[!, var]; - df = DataFrame(total = wsum(Float32.(x), weights(1 ./ w))); - return df -end - -function svytotal(var, w, popsize, sampsize) - df = DataFrame(total = wsum(Float32.(var), weights(1 ./ w))) - return df -end - -""" -The `svyby` function can be used to generate subsets of a survey design. - -```jldoctest -julia> using Survey - -julia> apiclus1 = load_data("apiclus1"); - -julia> dclus1 = svydesign(id=:1, weights=:pw, data = apiclus1); - -julia> svyby(:api00, :cname, dclus1, svymean) -11×3 DataFrame - Row │ cname mean SE - │ String15 Float64 Float64 -─────┼──────────────────────────────── - 1 │ Alameda 669.0 16.2135 - 2 │ Fresno 472.0 9.85278 - 3 │ Kern 452.5 29.5049 - 4 │ Los Angeles 647.267 23.5116 - 5 │ Mendocino 623.25 24.216 - 6 │ Merced 519.25 10.4925 - 7 │ Orange 710.562 28.9123 - 8 │ Plumas 709.556 13.2174 - 9 │ San Diego 659.436 12.2082 - 10 │ San Joaquin 551.189 11.578 - 11 │ Santa Clara 732.077 12.2291 -``` -""" -function svyby(formula::Symbol, by, design::svydesign, func::Function, params = []) - gdf = groupby(design.variables, by) - return combine(gdf, [formula, :probs, :popsize, :sampsize] => ((a, b, c, d) -> func(a, b, c, d, params...)) => AsTable) -end diff --git a/src/svydesign.jl b/src/svydesign.jl deleted file mode 100644 index 73321d74..00000000 --- a/src/svydesign.jl +++ /dev/null @@ -1,121 +0,0 @@ -function print_short(x) - if length(x) < 3 - print(x) - else - print( x[1], ", ", x[2], ", ", x[3], " ...", " (length = ", length(x), ")") - end -end -""" -The `svydesign` object combines a data frame and all the survey design information needed to analyse it. - -```jldoctest -julia> using Survey; - -julia> apiclus1 = load_data("apiclus1"); - -julia> dclus1 = svydesign(id= :dnum, weights= :pw, data = apiclus1, fpc= :fpc) |> print -Survey Design: -variables: 183x45 DataFrame -id: dnum -strata: 1, 1, 1 ... (length = 183) -probs: 0.029544719150814778, 0.029544719150814778, 0.029544719150814778 ... (length = 183) -fpc: - popsize: 757, 757, 757 ... (length = 183) - sampsize: 183, 183, 183 ... (length = 183) -nest: false -check_strat: true -``` -""" -struct svydesign - id - variables::DataFrame - nest::Bool - check_strat::Bool -end - -function get_weights(data, wt::Vector) - if nrow(data) == length(wt) - return wt - else - @error "length of the weights vector is not equal to the number of rows in the dataset" - end -end - -function get_weights(data, wt::Nothing) - return ones(nrow(data)) -end - -function get_weights(data, wt::Symbol) - wt = data[!, String(wt)] -end - -function get_probs(data, wt, probs::Symbol) - return data[!, String(probs)] -end - -function get_probs(data, wt, probs::Nothing) - return 1 ./ wt -end - -function get_fpc(data, fpc::Symbol) - return data[!, String(fpc)] -end - -function get_fpc(data, fpc::Nothing) - return repeat([nothing], nrow(data)) -end - -function get_fpc(data, fpc::Vector) - return fpc -end - -function get_fpc(data, fpc::Real) - return repeat([fpc], nrow(data)) -end - -function get_strata(data, strata::Symbol) - return data[!, String(strata)] -end - -function get_strata(data, strata::Vector) - return strata -end - -function get_strata(data, strata::Nothing) - return repeat([1], nrow(data)) -end - -function svydesign(; data = DataFrame(), id = Symbol(), probs = nothing, strata = nothing, fpc = nothing, nest = false, check_strat = !nest, weights = nothing) - wt = get_weights(data, weights) - if isnothing(probs) & isnothing(weights) - @warn "No weights or probabilities supplied, assuming equal probability" - end - df = data - df.probs = ProbabilityWeights(get_probs(data, wt, probs)) - df.weights = FrequencyWeights(get_weights(data, weights)) - df.popsize = get_fpc(data, fpc) - df.sampsize = repeat([nrow(data)], nrow(data)) - df.strata = get_strata(data, strata) - return svydesign(id, df, nest, check_strat) -end - -function Base.show(io::IO, design::svydesign) - printstyled("Survey Design:\n") - printstyled("variables: "; bold = true) - print(size(design.variables)[1], "x", size(design.variables)[2], " DataFrame") - printstyled("\nid: "; bold = true) - print(design.id) - printstyled("\nstrata: "; bold = true) - print_short(design.variables.strata) - printstyled("\nprobs: "; bold = true) - print_short(design.variables.probs) - printstyled("\nfpc: "; bold = true) - print("\n popsize: ") - print_short(design.variables.popsize) - print("\n sampsize: ") - print_short(design.variables.sampsize) - printstyled("\nnest: "; bold = true) - print(design.nest) - printstyled("\ncheck_strat: "; bold = true) - print(design.check_strat) -end \ No newline at end of file diff --git a/src/svyglm.jl b/src/svyglm.jl deleted file mode 100644 index 962bd0c6..00000000 --- a/src/svyglm.jl +++ /dev/null @@ -1,106 +0,0 @@ -function nullformula(f) - FormulaTerm(f.lhs,ConstantTerm(1)) -end - -mutable struct control - rtol - atol - maxiter -end - -""" -```julia -svyglm(formula, design, dist, link) -``` - -The `svyglm` function can be used to fit glms on svydesign. - -```jldoctest -julia> using Survey - -julia> apiclus1 = load_data("apiclus1"); - -julia> dclus1 = svydesign(id=:dnum, weights=:pw, data = apiclus1); - -julia> svyglm(@formula(ell~meals),dclus1,Normal(),IdentityLink()) - -``` -""" -mutable struct svyglm - glm - coefficients - data - weights #wrkwts - aic - family - rank - formula - model #subset of data as described in formula - deviance - offset - y - linear_predictors - fitted_values - prior_weights #initial weights - residuals - converged - control - terms - contrasts - naive_cov - df_null - null_deviance - df_residual - function svyglm_cons(glm, nullglm, data, weights,rtol,atol,maxiter) - out = new() - out.glm = glm - out.coefficients = coef(glm) - out.data = data - out.weights = glm.model.rr.wrkwt - out.aic = aic(glm) - out.family = glm.model.rr.d - out.rank = rank(glm.model.pp.X) - out.formula = formula(glm) - out.model = glm.mf.data - out.deviance = deviance(glm) - out.offset = glm.model.rr.offset - out.y = glm.model.rr.y - out.linear_predictors = predict(glm) - out.fitted_values = fitted(glm) - out.prior_weights = weights - out.residuals = glm.model.rr.wrkresid - out.converged = true - out.control = control(rtol,atol,maxiter) - out.terms = glm.mf.f - out.contrasts = [] - out.naive_cov = vcov(glm)/GLM.dispersion(glm.model,true) - out.df_null = dof_residual(nullglm) - out.null_deviance = deviance(nullglm) - out.df_residual = dof_residual(glm) - out - end - - function svyglm(formula, design, dist=Normal(), link=canonicallink(dist)) - data = design.variables - rtol = 1e-8 - atol = 1e-8 - maxiter = 30 - weights = 1 ./ data.probs - - glmout = glm(formula, data, dist, link, wts = weights, rtol = rtol, atol = atol, maxiter = maxiter) - nullglm = glm(nullformula(formula), data, dist, link, wts = weights, rtol = rtol, atol = atol, maxiter = maxiter) - svyglm_cons(glmout, nullglm, data, weights, rtol, atol, maxiter) - end - - svyglm(;formula, design, dist=Normal(), link=canonicallink(dist)) = svyglm(formula,design,dist,link) - -end - -function Base.show(io::IO, g::svyglm) - print(g.glm) - println("") - println("Degrees of Freedom: $(g.df_null) (i.e. Null); $(g.df_residual) Residual") - println("Null Deviance: $(g.null_deviance)") - println("Residual Deviance: $(g.deviance)") - println("AIC: $(g.aic)") -end diff --git a/src/svyhist.jl b/src/svyhist.jl deleted file mode 100644 index 693c0b3f..00000000 --- a/src/svyhist.jl +++ /dev/null @@ -1,154 +0,0 @@ -""" - sturges(v) - -Calculate the number of bins to use in a histogram using the Sturges rule. - -# Examples -```jldoctest -julia> sturges(10) -5 - -julia> sturges([10, 20, 30, 40, 50]) -4 -``` -""" -sturges(n::Integer) = ceil(Int, log2(n)) + 1 -sturges(vec::AbstractVector) = ceil(Int, log2(length(vec))) + 1 - -""" - sturges(df::DataFrame, var::Symbol) - -Calculate the number of bins for a `DataFrame` variable. - -# Examples -```jldoctest -julia> using DataFrames - -julia> df = DataFrame((a=[1, 2, 3, 4, 5], b=[10, 20, 30, 40, 50])) -5×2 DataFrame - Row │ a b - │ Int64 Int64 -─────┼────────────── - 1 │ 1 10 - 2 │ 2 20 - 3 │ 3 30 - 4 │ 4 40 - 5 │ 5 50 - -julia> sturges(df, :b) -4 -``` -""" -sturges(df::DataFrame, var::Symbol) = ceil(Int, log2(size(df[!, var], 1))) + 1 - -""" - sturges(design::svydesign, var::Symbol) - -Calculate the number of bins for a survey design variable. - -# Examples -```jldoctest -julia> apistrat = load_data("apistrat"); - -julia> dstrat = svydesign(data = apistrat, id = :1, strata = :stype, weights = :pw, fpc = :fpc); - -julia> sturges(dstrat, :enroll) -9 -``` -""" -sturges(design::svydesign, var::Symbol) = sturges(design.variables, var) - -""" - freedman_diaconis(v::AbstractVector) - -Calculate the number of bins to use in a histogram using the Freedman-Diaconis rule. - -# Examples -```jldoctest -julia> freedman_diaconis([10, 20, 30, 40, 50]) -2 -``` -""" -freedman_diaconis(v::AbstractVector) = round(Int, length(v)^(1 / 3) * (maximum(v) - minimum(v)) / (2 * iqr(v))) - -""" - freedman_diaconis(df::DataFrame, var::Symbol) - -Calculate the number of bins for a `DataFrame` variable. - -# Examples -```jldoctest -julia> using DataFrames - -julia> df = DataFrame((a=[1, 2, 3, 4, 5], b=[10, 20, 30, 40, 50])); - -julia> freedman_diaconis(df, :b) -2 -``` -""" -freedman_diaconis(df::DataFrame, var::Symbol) = freedman_diaconis(df[!, var]) - -""" - freedman_diaconis(design::svydesign, var::Symbol) - -Calculate the number of bins for a survey design variable. - -# Examples -```jldoctest -julia> apistrat = load_data("apistrat"); - -julia> dstrat = svydesign(data = apistrat, id = :1, strata = :stype, weights = :pw, fpc = :fpc); - -julia> freedman_diaconis(dstrat, :enroll) -15 -``` -""" -freedman_diaconis(design::svydesign, var::Symbol) = freedman_diaconis(design.variables[!, var]) - -""" -```julia -svyhist(design, var, bins = freedman_diaconis; normalization = :density, weights = ones(size(design.variables, 1), ...) -``` -Histogram plot of a survey design variable given by `var`. - -`bins` can be an `Integer` specifying the number of equal-width bins, -an `AbstractVector` specifying the bins intervals, or a `Function` specifying -the function used for calculating the number of bins. The possible functions -are `sturges` and `freedman_diaconis`. - -The normalization can be set to `:none`, `:density`, `:probability` or `:pdf`. -See [Makie.hist](https://makie.juliaplots.org/stable/examples/plotting_functions/hist/) -for more information. - -The `weights` argument can be an `AbstractVector` or a `Symbol` specifying a -design variable. - -For the complete argument list see [Makie.hist](https://makie.juliaplots.org/stable/examples/plotting_functions/hist/). - -```@example e1 -julia> using Survey - -julia> apistrat = load_data("apistrat"); - -julia> dstrat = svydesign(data = apistrat, id = :1, strata = :stype, weights = :pw, fpc = :fpc); - -julia> h = svyhist(dstrat, :enroll) -``` - -![](./assets/hist.png) -""" -function svyhist(design::svydesign, var::Symbol, - bins::Union{Integer, AbstractVector} = freedman_diaconis(design, var); - normalization = :density, - kwargs... - ) - hist = histogram(bins = bins, normalization = normalization, kwargs...) - data(design.variables) * mapping(var, weights = :weights) * hist |> draw -end - -function svyhist(design::svydesign, var::Symbol, - bins::Function; - kwargs... - ) - svyhist(design, var, bins(design, var); kwargs...) -end diff --git a/src/svyplot.jl b/src/svyplot.jl deleted file mode 100644 index 654cec0e..00000000 --- a/src/svyplot.jl +++ /dev/null @@ -1,24 +0,0 @@ -""" -``` - svyplot(design, x, y; kwargs...) -``` -Scatter plot of survey design variables `x` and `y`. - -The plot takes into account the frequency weights specified by the user -in the design. - -```@example svyplot -julia> using Survey - -julia> apistrat = load_data("apistrat"); - -julia> dstrat = svydesign(data = apistrat, id = :1, strata = :stype, weights = :pw, fpc = :fpc); - -julia> s = svyplot(dstrat, :api99, :api00) -``` - -![](./assets/scatter.png) -""" -function svyplot(design::svydesign, x::Symbol, y::Symbol; kwargs...) - data(design.variables) * mapping(x, y, markersize = :weights) * visual(Scatter, marker = '￮') |> draw -end diff --git a/src/total.jl b/src/total.jl new file mode 100644 index 00000000..be8879cb --- /dev/null +++ b/src/total.jl @@ -0,0 +1,123 @@ +""" + total(x, design) + +Estimate the population total for the variable specified by `x`. + +```jldoctest +julia> using Survey; + +julia> apisrs = load_data("apisrs"); + +julia> srs = SimpleRandomSample(apisrs; popsize=:fpc); + +julia> total(:enroll, srs) +1×2 DataFrame + Row │ total SE + │ Float64 Float64 +─────┼───────────────────── + 1 │ 3.62107e6 1.6952e5 + +julia> strat = load_data("apistrat"); + +julia> dstrat = StratifiedSample(strat, :stype; popsize=:fpc); + +julia> total(:api00, dstrat) +1×2 DataFrame + Row │ total SE + │ Float64 Float64 +─────┼──────────────────── + 1 │ 4.10221e6 58279.0 + +julia> total([:api00, :enroll], dstrat) +2×3 DataFrame + Row │ names total SE + │ String Float64 Float64 +─────┼────────────────────────────────── + 1 │ api00 4.10221e6 58279.0 + 2 │ enroll 3.68718e6 1.14642e5 +``` +""" +function total(x::Symbol, design::SimpleRandomSample) + function se(x::Symbol, design::SimpleRandomSample) + function variance(x::Symbol, design::SimpleRandomSample) + return design.popsize^2 * design.fpc * var(design.data[!, x]) / design.sampsize + end + return sqrt(variance(x, design)) + end + if isa(design.data[!, x], CategoricalArray) + gdf = groupby(design.data, x) + p = combine(gdf, nrow => :count) + p.total = design.popsize .* p.count ./ sum(p.count) + p.proportion = p.total ./ design.popsize + p = select!(p, Not(:count)) # count column is not necessary for `total` + p.var = design.popsize^2 .* design.fpc .* p.proportion .* + (1 .- p.proportion) ./ (design.sampsize - 1) # N^2 .* variance of proportion + p.SE = sqrt.(p.var) + return select(p, Not([:proportion, :var])) + end + m = mean(x,design) + total = design.popsize * m.mean[1] + return DataFrame(total=total, SE=se(x, design)) +end + +function total(x::Symbol, design::StratifiedSample) + # TODO: check if statement + if x == design.strata + gdf = groupby(design.data, x) + return combine(gdf, :weights => sum => :Nₕ) + end + gdf = groupby(design.data, design.strata) + grand_total = sum(combine(gdf, [x, :weights] => ((a, b) -> wsum(a, b)) => :total).total) + # variance estimation using closed-form formula + Nₕ = combine(gdf, :weights => sum => :Nₕ).Nₕ + nₕ = combine(gdf, nrow => :nₕ).nₕ + fₕ = nₕ ./ Nₕ + + s²ₕ = combine(gdf, x => var => :s²h).s²h + # the only difference between total and mean variance is the Nₕ instead of Wₕ + V̂Ȳ̂ = sum((Nₕ .^ 2) .* (1 .- fₕ) .* s²ₕ ./ nₕ) + SE = sqrt(V̂Ȳ̂) + return DataFrame(total=grand_total, SE=SE) +end + +function total(x::Vector{Symbol}, design::AbstractSurveyDesign) + df = reduce(vcat, [total(i, design) for i in x]) + insertcols!(df, 1, :names => String.(x)) + return df +end + +""" + total(x, by, design) + +Estimate the subpopulation total of a variable `x`. + +```jldoctest +julia> using Survey; + +julia> apisrs = load_data("apisrs"); + +julia> srs = SimpleRandomSample(apisrs;popsize=:fpc); + +julia> total(:api00, :cname, srs) |> first +DataFrameRow + Row │ cname total SE + │ String15 Float64 Float64 +─────┼────────────────────────────── + 1 │ Kern 1.77644e5 55600.8 + +``` +""" +function total(x::Symbol, by::Symbol, design::SimpleRandomSample) + function domain_total(x::AbstractVector, design::SimpleRandomSample, weights) + function se(x::AbstractVector, design::SimpleRandomSample) + # vector of length equal to `sampsize` containing `x` and zeros + z = cat(zeros(design.sampsize - length(x)), x; dims=1) + variance = design.popsize^2 / design.sampsize * design.fpc * var(z) + return sqrt(variance) + end + total = wsum(x, weights) + return DataFrame(total=total, SE=se(x, design::SimpleRandomSample)) + end + gdf = groupby(design.data, by) + combine(gdf, [x, :weights] => ((a, b) -> domain_total(a, design, b)) => AsTable) +end diff --git a/test/Project.toml b/test/Project.toml index 0c363327..62dd5dde 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,2 +1,3 @@ [deps] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +CategoricalArrays = "324d7699-5711-5eae-9e2f-1d82baa6b597" diff --git a/test/SurveyDesign.jl b/test/SurveyDesign.jl new file mode 100644 index 00000000..fa882a31 --- /dev/null +++ b/test/SurveyDesign.jl @@ -0,0 +1,166 @@ +# Work on copies, keep original +@testset "SimpleRandomSample" begin + ##### SimpleRandomSample tests + # Load API datasets + apisrs_original = load_data("apisrs") + apisrs_original[!, :derived_probs] = 1 ./ apisrs_original.pw + apisrs_original[!, :derived_sampsize] = fill(200.0, size(apisrs_original, 1)) + ############################## + ### Valid type checking tests + apisrs = copy(apisrs_original) + @test_throws TypeError SimpleRandomSample(apisrs, popsize=-2.83, ignorefpc=true) + @test_throws TypeError SimpleRandomSample(apisrs, sampsize=-300) + @test_throws TypeError SimpleRandomSample(apisrs, sampsize=-2.8, ignorefpc=true) + @test_throws TypeError SimpleRandomSample(apisrs, weights=50) + @test_throws TypeError SimpleRandomSample(apisrs, probs=1) + ############################## + ### weights or probs as Symbol + apisrs = copy(apisrs_original) + srs_weights = SimpleRandomSample(apisrs; weights=:pw) + @test srs_weights.data.weights[1] ≈ 30.97 atol = 1e-4 + @test srs_weights.data.weights == 1 ./ srs_weights.data.probs + ### probs as Symbol + apisrs = copy(apisrs_original) + srs_probs_sym = SimpleRandomSample(apisrs; probs=:derived_probs) + @test srs_probs_sym.data.probs[1] ≈ 0.032289 atol = 1e-4 + @test srs_probs_sym.data.probs == 1 ./ srs_probs_sym.data.weights + ############################## + ### Weights or probs as non-numeric error + apisrs = copy(apisrs_original) + @test_throws ErrorException SimpleRandomSample(apisrs, weights=:stype) + @test_throws ErrorException SimpleRandomSample(apisrs, probs=:cname) + ############################## + ### popsize given as Symbol + apisrs = copy(apisrs_original) + srs_popsize_sym = SimpleRandomSample(apisrs; popsize=:fpc) + @test srs_popsize_sym.data.weights == 1 ./ srs_popsize_sym.data.probs # weights should be inverse of probs + @test srs_popsize_sym.sampsize > 0 + ### popsize given as Vector + apisrs = copy(apisrs_original) + srs_popsize_vec = SimpleRandomSample(apisrs; popsize=apisrs.fpc) + @test srs_popsize_vec.data.weights == 1 ./ srs_popsize_vec.data.probs # weights should be inverse of probs + @test srs_popsize_vec.sampsize > 0 + ############################## + ### sampsize given as Symbol + apisrs = copy(apisrs_original) + srs_sampsize_sym = SimpleRandomSample(apisrs; sampsize=:derived_sampsize, weights=:pw) + @test srs_sampsize_sym.data.weights == 1 ./ srs_sampsize_sym.data.probs # weights should be inverse of probs + @test srs_sampsize_sym.sampsize > 0 + ### sampsize given as Vector + apisrs = copy(apisrs_original) + srs_sampsize_vec = SimpleRandomSample(apisrs; sampsize=apisrs.derived_sampsize, probs=:derived_probs) + @test srs_sampsize_vec.data.weights == 1 ./ srs_sampsize_vec.data.probs # weights should be inverse of probs + @test srs_sampsize_vec.sampsize > 0 + ############################## + ### both weights and probs given + # If weights given, probs is superfluous + apisrs = copy(apisrs_original) + srs_weights_probs = SimpleRandomSample(apisrs; weights=:pw, probs=:derived_probs) + srs_weights_probs = SimpleRandomSample(apisrs; weights=:pw, probs=:pw) + ############################## + ### sum of weights and probs condition check + apisrs = copy(apisrs_original) + @test_throws ErrorException SimpleRandomSample(apisrs, weights=fill(0.3, size(apisrs_original, 1))) + apisrs = copy(apisrs_original) + @test_throws ErrorException SimpleRandomSample(apisrs, probs=fill(0.3, size(apisrs_original, 1))) + ############################## + ### weights only as Vector + apisrs = copy(apisrs_original) + srs_weights = SimpleRandomSample(apisrs; weights=apisrs.pw) + @test srs_weights.data.weights[1] == 30.97 + @test srs_weights.data.weights == 1 ./ srs_weights.data.probs + ### probs only as Vector + apisrs = copy(apisrs_original) + srs_freq = SimpleRandomSample(apisrs; probs=apisrs.derived_probs) + @test srs_freq.data.weights[1] == 30.97 + @test srs_freq.data.weights == 1 ./ srs_freq.data.probs + ############################## + ### ignorefpc tests. TODO: change if ignorefpc functionality changed + apisrs = copy(apisrs_original) + srs_ignorefpc = SimpleRandomSample(apisrs; popsize=:fpc, ignorefpc=true) + @test srs_ignorefpc.data.weights == 1 ./ srs_ignorefpc.data.probs # weights should be inverse of probs + @test srs_ignorefpc.sampsize > 0 + ### incorrect probs with correct popsize, ignorefpc = true + apisrs = copy(apisrs_original) + srs_w_p = SimpleRandomSample(apisrs, popsize=:fpc, probs=fill(0.3, size(apisrs_original, 1)), ignorefpc=true) + @test srs_w_p.data.probs == 1 ./ srs_w_p.data.weights + ### ingorefpc = true with probs given + apisrs = copy(apisrs_original) + srs = SimpleRandomSample(apisrs, ignorefpc=true, probs=:derived_probs) + @test srs.data.probs == 1 ./ srs.data.weights + ############################## + ### probs as vector declared on-the-fly + apisrs = copy(apisrs_original) + srs_prob = SimpleRandomSample(apisrs; probs=1 ./ apisrs.pw) + @test srs_prob.data.weights[1] == 30.97 + @test srs_prob.data.weights == 1 ./ srs_prob.data.probs +end + +@testset "StratifiedSample" begin + ### StratifiedSample tests + # Load API datasets + apistrat_original = load_data("apistrat") + apistrat_original[!, :derived_probs] = 1 ./ apistrat_original.pw + apistrat_original[!, :derived_sampsize] = apistrat_original.fpc ./ apistrat_original.pw + ############################## + ### Valid type checking tests + apistrat = copy(apistrat_original) + @test_throws TypeError StratifiedSample(apistrat,:stype; popsize=-2.83, ignorefpc=true) + @test_throws TypeError StratifiedSample(apistrat,:stype; sampsize=-300) + @test_throws TypeError StratifiedSample(apistrat,:stype; sampsize=-2.8, ignorefpc=true) + @test_throws TypeError StratifiedSample(apistrat,:stype; weights=50) + @test_throws TypeError StratifiedSample(apistrat,:stype; probs=1) + ############################## + ### weights as Symbol + apistrat = copy(apistrat_original) + strat_wt = StratifiedSample(apistrat, :stype; weights=:pw) + @test strat_wt.data.probs == 1 ./ strat_wt.data.weights + ### probs as Symbol + apistrat = copy(apistrat_original) + strat_probs = StratifiedSample(apistrat, :stype; probs=:derived_probs) + @test strat_probs.data.probs == 1 ./ strat_probs.data.weights + ### weights as Vector{<:Real} + apistrat = copy(apistrat_original) + strat_wt = StratifiedSample(apistrat, :stype; weights=apistrat.pw) + @test strat_wt.data.probs == 1 ./ strat_wt.data.weights + ### probs as Vector{<:Real} + apistrat = copy(apistrat_original) + strat_probs = StratifiedSample(apistrat, :stype; probs=apistrat.derived_probs) + @test strat_probs.data.probs == 1 ./ strat_probs.data.weights + ############################## + ### popsize as Symbol + apistrat = copy(apistrat_original) + strat_pop = StratifiedSample(apistrat, :stype; popsize=:fpc) + @test strat_pop.data.probs == 1 ./ strat_pop.data.weights + ### popsize given as Vector (should give error for now, not implemented Vector input directly for popsize) + apistrat = copy(apistrat_original) + @test_throws TypeError StratifiedSample(apistrat,:stype; popsize=apistrat.fpc) + ############################## + ### sampsize given as Symbol + apistrat = copy(apistrat_original) + strat_sampsize_sym = StratifiedSample(apistrat,:stype; sampsize=:derived_sampsize, weights=:pw) + @test strat_sampsize_sym.data.weights == 1 ./ strat_sampsize_sym.data.probs # weights should be inverse of probs + ### sampsize given as symbol without weights or probs, and popsize not given - raise error + apistrat = copy(apistrat_original) + @test_throws ErrorException StratifiedSample(apistrat,:stype; sampsize=:derived_sampsize) + ############################## + ### both weights and probs given + # If weights given, probs is superfluous + apistrat = copy(apistrat_original) + strat_weights_probs = StratifiedSample(apistrat,:stype; weights=:pw, probs=:derived_probs) + strat_weights_probs = StratifiedSample(apistrat,:stype; weights=:pw, probs=:pw) + ############################## + ### ignorefpc test (Modify if ignorefpc changed) + apistrat = copy(apistrat_original) + strat_ignorefpc=StratifiedSample(apistrat,:stype; popsize=:fpc, ignorefpc=true) + @test strat_ignorefpc.data.probs == 1 ./ strat_ignorefpc.data.weights + ############################## + # For now, no sum checks on probs and weights for StratifiedSample (unlike SRS) + apistrat = copy(apistrat_original) + strat_probs1 = StratifiedSample(apistrat, :stype; probs=fill(0.3, size(apistrat, 1))) + @test strat_probs1.data.probs == 1 ./ strat_probs1.data.weights + ############################## + #should throw error because sampsize > popsize + apistrat = copy(apistrat_original) + @test_throws ErrorException StratifiedSample(apistrat, :stype; popsize= :pw, sampsize=:fpc) +end \ No newline at end of file diff --git a/test/boxplot.jl b/test/boxplot.jl new file mode 100644 index 00000000..f40dd87f --- /dev/null +++ b/test/boxplot.jl @@ -0,0 +1,11 @@ +@testset "boxplot.jl" begin + # SimpleRandomSample + apisrs = load_data("apisrs") + srs = SimpleRandomSample(apisrs,popsize = apisrs.fpc) + bp = boxplot(srs, :stype, :enroll; weights = :pw) + + @test bp.grid[1].entries[1].positional[2] == srs.data[!, :enroll] + @test bp.grid[1].entries[1].named[:weights] == srs.data[!, :pw] + + # StratifiedSample +end diff --git a/test/dimnames.jl b/test/dimnames.jl index 3c631419..6241bd64 100644 --- a/test/dimnames.jl +++ b/test/dimnames.jl @@ -1,15 +1,16 @@ -using Survey -using Test - @testset "dimnames.jl" begin - apistrat = load_data("apistrat") - dstrat = svydesign(data = apistrat, id = :1, strata = :stype, weights = :pw, fpc = :fpc) - - @test dim(dstrat)[1] == 200 - @test dim(dstrat)[2] == size(dstrat.variables, 2) - - @test length(colnames(dstrat)) == dim(dstrat)[2] + # Simple random sampling tests + apisrs = load_data("apisrs") + # make a copy to not modify the original dataset + apisrs_copy = copy(apisrs) + srs = SimpleRandomSample(apisrs_copy,popsize=:fpc,ignorefpc = true) + # `dim` + @test dim(srs)[2] == 42 + # `colnames` + @test length(colnames(srs)) == dim(srs)[2] + # `dimnames` + @test length(dimnames(srs)[1]) == parse(Int, last(dimnames(srs)[1])) + @test dimnames(srs)[2] == colnames(srs) - @test length(dimnames(dstrat)[1]) == parse(Int, last(dimnames(dstrat)[1])) - @test dimnames(dstrat)[2] == colnames(dstrat) + # Stratified sampling tests end diff --git a/test/hist.jl b/test/hist.jl new file mode 100644 index 00000000..64e260b9 --- /dev/null +++ b/test/hist.jl @@ -0,0 +1,22 @@ +@testset "hist.jl" begin + @test Survey.sturges(10) == 5 + @test Survey.sturges([1, 2, 5, 10, 15, 17, 20]) == 4 + + # SimpleRandomSample + apisrs = load_data("apisrs") + srs = SimpleRandomSample(apisrs,popsize=:fpc) + + h = hist(srs, :enroll) + @test h.grid[1].entries[1].positional[2] |> length == 21 + + h = hist(srs, :enroll, 9) + @test h.grid[1].entries[1].positional[2] |> length == 11 + + h = hist(srs, :enroll, Survey.sturges) + @test h.grid[1].entries[1].positional[2] |> length == 11 + + h = hist(srs, :enroll, [0, 1000, 2000, 3000]) + @test h.grid[1].entries[1].positional[2] |> length == 3 + + # StratifiedSample +end diff --git a/test/mean.jl b/test/mean.jl new file mode 100644 index 00000000..1221a0aa --- /dev/null +++ b/test/mean.jl @@ -0,0 +1,53 @@ +@testset "mean_SimpleRandomSample" begin + ##### SimpleRandomSample tests + # Load API datasets + apisrs_original = load_data("apisrs") + apisrs_original[!, :derived_probs] = 1 ./ apisrs_original.pw + apisrs_original[!, :derived_sampsize] = fill(200.0, size(apisrs_original, 1)) + ############################## + ### Basic functionality + apisrs = copy(apisrs_original) + srs = SimpleRandomSample(apisrs, popsize = :fpc) + @test mean(:api00, srs).mean[1] ≈ 656.585 atol = 1e-4 + @test mean(:api00, srs).SE[1] ≈ 9.249722039282807 atol = 1e-4 + @test mean(:enroll, srs).mean[1] ≈ 584.61 atol = 1e-4 + @test mean(:enroll, srs).SE[1] ≈ 27.36836524766856 atol = 1e-4 + # ignorefpc = true + apisrs = copy(apisrs_original) + srs = SimpleRandomSample(apisrs, popsize=:fpc,ignorefpc = true) + @test mean(:api00, srs).mean[1] ≈ 656.585 atol = 1e-4 + @test mean(:api00, srs).SE[1] ≈ 9.402772170880636 atol = 1e-4 + @test mean(:enroll, srs).mean[1] ≈ 584.61 atol = 1e-4 + @test mean(:enroll, srs).SE[1] ≈ 27.821214737089324 atol = 1e-4 + ############################## + ### Vector of Symbols + apisrs = copy(apisrs_original) + srs = SimpleRandomSample(apisrs, popsize = :fpc) + mean_vec_sym = mean([:api00,:enroll], srs) + @test mean_vec_sym.mean[1] ≈ 656.585 atol = 1e-4 + @test mean_vec_sym.SE[1] ≈ 9.249722039282807 atol = 1e-4 + @test mean_vec_sym.mean[2] ≈ 584.61 atol = 1e-4 + @test mean_vec_sym.SE[2] ≈ 27.36836524766856 atol = 1e-4 + ############################## + ### Categorical Array - estimating proportions + apisrs_categ = copy(apisrs_original) + apisrs_categ.stype = CategoricalArray(apisrs_categ.stype) # Convert a column to CategoricalArray + srs_design_categ = SimpleRandomSample(apisrs_categ, popsize = :fpc) + #>>>>>>>>> complete this suite + mean_categ = mean(:stype,srs_design_categ) + # complete this +end + +@testset "mean_Stratified" begin + ## Add tests +end + +@testset "mean_svyby_SimpleRandomSample" begin + ## Add tests +end + +@testset "mean_svyby_Stratified" begin + ## Add tests +end + + diff --git a/test/plot.jl b/test/plot.jl new file mode 100644 index 00000000..e9e59c36 --- /dev/null +++ b/test/plot.jl @@ -0,0 +1,10 @@ +@testset "plot.jl" begin + # SimpleRandomSample + apisrs = load_data("apisrs") + srs = SimpleRandomSample(apisrs,popsize=:fpc) + s = plot(srs, :api99, :api00) + @test s.grid[1].entries[1].named[:markersize] == srs.data.weights + @test s.grid[1].entries[1].positional[1] == srs.data.api99 + @test s.grid[1].entries[1].positional[2] == srs.data.api00 + # StratifiedSample +end diff --git a/test/quantile.jl b/test/quantile.jl new file mode 100644 index 00000000..cab18fdb --- /dev/null +++ b/test/quantile.jl @@ -0,0 +1,34 @@ +@testset "quantile_SimpleRandomSample" begin + ##### SimpleRandomSample tests + # Load API datasets + apisrs_original = load_data("apisrs") + apisrs_original[!, :derived_probs] = 1 ./ apisrs_original.pw + apisrs_original[!, :derived_sampsize] = fill(200.0, size(apisrs_original, 1)) + ############################## + apisrs = copy(apisrs_original) + srs_design = SimpleRandomSample(apisrs,popsize=:fpc) + @test quantile(:api00,srs_design,0.5)[!,2][1] ≈ 659.0 atol=1e-4 + @test quantile(:api00,srs_design,[0.1753,0.25,0.5,0.75,0.975])[!,2] ≈ [512.8847,544,659,752.5,905] atol = 1e-4 + @test quantile(:enroll,srs_design,[0.1,0.2,0.5,0.75,0.95])[!,2] ≈ [245.5,317.6,453.0,668.5,1473.1] atol = 1e-4 +end + +@testset "quantile_Stratified" begin + ##### StratifiedSample tests + # Load API datasets + apistrat_original = load_data("apistrat") + apistrat_original[!, :derived_probs] = 1 ./ apistrat_original.pw + apistrat_original[!, :derived_sampsize] = apistrat_original.fpc ./ apistrat_original.pw + # base functionality + apistrat = copy(apistrat_original) + dstrat = StratifiedSample(apistrat, :stype; popsize = :fpc) + # Check which definition of quantile for StratifiedSample + # @test quantile(:enroll,dstrat,[0.1,0.2,0.5,0.75,0.95])[!,2] ≈ [262,309.3366,446.4103,658.8764,1589.7881] atol = 1e-4 +end + +@testset "quantile_by_SimpleRandomSample" begin + ## Add tests +end + +@testset "quantile_by_Stratified" begin + ## Add tests +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 59ab86db..00e850b5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,24 +1,19 @@ using Survey using Test +using CategoricalArrays @testset "Survey.jl" begin apiclus1 = load_data("apiclus1") @test size(apiclus1) == (183, 40) @test size(load_data("apiclus2")) == (126, 41) @test size(load_data("apipop")) == ((6194, 38)) - dclus1 = svydesign(id=:1, strata=:stype, weights=:pw, data = apiclus1, fpc=:fpc) - @test dclus1.variables.strata[1] == "H" - @test length(dclus1.variables.probs) == 183 - @test dclus1.id == 1 - api00_by_cname = svyby(:api00, :cname, dclus1, svymean).mean - @test api00_by_cname ≈ [669.0000000000001, 472.00000000000006, 452.5, 647.2666666666668, 623.25, 519.25, 710.5625000000001, 709.5555555555557, 659.4363636363635, 551.1891891891892, 732.0769230769226] - api00_by_cname_meals = svyby(:api00, [:cname, :meals], dclus1, svymean) - @test api00_by_cname_meals[1,3] ≈ 608.0 end -include("svyglm.jl") -include("svyhist.jl") -include("svyplot.jl") +include("SurveyDesign.jl") +include("total.jl") +include("quantile.jl") +include("mean.jl") include("dimnames.jl") -include("svyplot.jl") -include("svyboxplot.jl") +include("plot.jl") +include("hist.jl") +include("boxplot.jl") \ No newline at end of file diff --git a/test/sampling.jl b/test/sampling.jl new file mode 100644 index 00000000..276ed8a1 --- /dev/null +++ b/test/sampling.jl @@ -0,0 +1,3 @@ +""" + Testing suite for sampling functions +""" \ No newline at end of file diff --git a/test/svyboxplot.jl b/test/svyboxplot.jl deleted file mode 100644 index bfbeed48..00000000 --- a/test/svyboxplot.jl +++ /dev/null @@ -1,11 +0,0 @@ -using Survey -using Test - -@testset "svyplot.jl" begin - apistrat = load_data("apistrat") - dstrat = svydesign(data = apistrat, id = :1, strata = :stype, weights = :pw, fpc = :fpc) - bp = svyboxplot(dstrat, :stype, :enroll; weights = :pw) - - @test bp.grid[1].entries[1].positional[2] == dstrat.variables[!, :enroll] - @test bp.grid[1].entries[1].named[:weights] == dstrat.variables[!, :pw] -end diff --git a/test/svyglm.jl b/test/svyglm.jl deleted file mode 100644 index a8e871c4..00000000 --- a/test/svyglm.jl +++ /dev/null @@ -1,2411 +0,0 @@ -using Survey -using Test - -function genWts(data,chaos) - if chaos==1 - wts = ones(size(data)[1]) - else - wts = ones(size(data)[1]) - wts[1] = 200 - end - wts -end - -@testset "svyglm.jl" begin - rtol=0.1 #rtol which test uses is diffrent from - apiclus1 = load_data("apiclus1") - apiclus2 = load_data("apiclus2") - apipop = load_data("apipop") - apistrat = load_data("apistrat") - apisrs = load_data("apisrs") - - rcoef = [19.852539307789417, 0.016802996766905515] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [13.821691894538569, 0.020932345630754703] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.04873273798101878, -2.545900442141623e-5] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.06602717760486693, -4.837774447484219e-5] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.000475022138835, 0.0006674250431904819] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [2.626328291498775, 0.001128251554586315] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.4675561050493435, 0.0016794583729629272] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.6919682339147015, 0.0024858792643755187] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [20.81754723634115, 0.014574128759780829] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [16.9504183227201, 0.013414574139003429] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.047326430818283034, -2.273224677320802e-5] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.06092784931279272, -3.790260657717406e-5] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.0427930460495958, 0.0005781143228868981] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [2.8081808689622867, 0.0007223757899528744] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.5709143099027445, 0.001451428755866715] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.093801084880702, 0.001554601904188512] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [20.396962685204315, 0.015578775639845489] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [15.757897032885943, 0.016354051972422335] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.04802942131078968, -2.4136270761976737e-5] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.06396051126631712, -4.4334324326067683e-5] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.023674047116632, 0.0006198204624428633] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [2.7306451504600457, 0.0009013786283701581] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.525317648001372, 0.0015554611944026753] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.937974903181009, 0.001926878861268805] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.476009259711357, 0.45896567043916364, -0.2189808315090155, 0.016762348305265986] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(growth~meals+ell+dnum),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [7.405836020305114, 0.4072566213513245, -0.15502327824484397, 0.014823594631887456] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(growth~meals+ell+dnum),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [17.499309726273815, 1.0362247060905792, 0.20689279530841123] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [10.618862171412498, 0.9470888770809175, 0.2554170143650083] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.03105097698342063, -0.00025951964795808497, -4.980649281479226e-5] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.03564902053538176, -0.00029846548118673334, -7.127665001634578e-5] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.200988225159865, 0.01853539304483344, 0.004401765720794799] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.016183666106752, 0.019353379360543335, 0.005674410063167833] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.4050112100743, 0.07570792572126114, 0.01816200488597775] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.972520325284209, 0.07053682209109562, 0.02109203630823282] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [96.58554019967326, 0.09528706473480851, 0.04746149077760165, -0.29581091687660227] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [102.56903957327081, 0.2994671678196024, 0.05694530556083325, -0.45227877815926576] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.010298668498595407, -1.1475337581433332e-5, -6.476114427094932e-6, 3.866960426539555e-5] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Normal(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.00971042834048784, -4.293936615201201e-5, -6.686320985047049e-6, 5.7285543560066516e-5] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Normal(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.573015337880906, 0.0010496565303077222, 0.0005597422989190293, -0.0033886072141109664] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Normal(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.632906410122461, 0.0036298370322293233, 0.0006277035681771951, -0.005127700303215906] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Normal(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.834055808096808, 0.0050047147493627675, 0.00258380795582847, -0.015837488847978982] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Normal(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [10.13399218244361, 0.016534390287256482, 0.0030050612704961583, -0.024121159310513514] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Normal(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [97.04322155899567, 0.11444814943941939, 0.04894444703603661, -0.3154905816230446] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Gamma(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [103.36576358179042, 0.27958665616726947, 0.06210745944693044, -0.4640788209061506] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Gamma(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.010259348226255835, -1.3551287902882181e-5, -6.804004558645566e-6, 4.085883946939977e-5] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Gamma(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.009620561658963203, -4.1340773001206705e-5, -7.77958780972283e-6, 5.957478962653988e-5] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Gamma(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.57722584780227, 0.0012472563878204615, 0.0005829761882771319, -0.0035939835822592836] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Gamma(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.641357580037837, 0.0034312167113636086, 0.00071388149921673, -0.005291340773091543] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Gamma(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.855962779482573, 0.005975538095906088, 0.0026782307142671764, -0.016840058593652295] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Gamma(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [10.174989337839179, 0.015526010287809732, 0.003356942476353457, -0.02481863615015127] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Gamma(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [96.78174537092866, 0.10403982721444952, 0.04811675254370905, -0.30465358460709974] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Poisson(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [102.93852488351948, 0.28896639160064413, 0.0597468338664144, -0.4574349763919181] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Poisson(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.01028153721826924, -1.2443211144164435e-5, -6.630761733850919e-6, 3.967635761261705e-5] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.009668490487342555, -4.2193491493239046e-5, -7.206862262535231e-6, 5.8378528534779286e-5] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.574836900133671, 0.0011409175094760585, 0.0005704953953500102, -0.0034820528248527775] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.636835034133741, 0.0035296478514372305, 0.0006704849256610509, -0.005202532036672366] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.843493008377969, 0.005450655556148467, 0.0026267322072815956, -0.016290879285155466] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [10.153024510485578, 0.016012306833782283, 0.003185833666927896, -0.024432788452762772] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [93.82000643799564, 0.032714209742078004, -0.027322045119543895] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [93.82846194140612, 0.032704317447096384, -0.027327105792050974] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.01063501615215185, -3.3548983370761818e-6, 2.7996866875279524e-6] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Normal(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.01063573587989556, -3.355628580324203e-6, 2.7991420970934927e-6] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Normal(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.542522453794568, 0.0003313304964353771, -0.0002766121742660158] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Normal(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.542530169818557, 0.000331322066236489, -0.00027661739942401377] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Normal(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.68888483545306, 0.0016461974399421373, -0.001374601593992681] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Normal(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.689117241108823, 0.0016459348302332616, -0.0013747501452265492] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Normal(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [93.7681632462419, 0.03309794966873599, -0.02763990805626656] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Gamma(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [93.7948170030615, 0.03306397401308171, -0.0276529878883247] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Gamma(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.010640396277748213, -3.3949562818952374e-6, 2.8329786089132253e-6] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Gamma(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.010639296824781657, -3.393730062080016e-6, 2.833696372778119e-6] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Gamma(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.54199424095025, 0.0003352523414203465, -0.00027986631489690383] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Gamma(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.542184017118389, 0.00033502548614321057, -0.000279974732243069] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Gamma(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.686268493077382, 0.0016655938946646084, -0.0013906821665282774] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Gamma(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.687411021387943, 0.0016641834373892195, -0.0013912894932434584] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Gamma(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [93.79425168307411, 0.03290474796528548, -0.02747990117196297] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Poisson(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [93.81169317983408, 0.032883440833357924, -0.027489412238860728] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Poisson(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.010637689518962693, -3.374790155020553e-6, 2.8162208220632363e-6] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.010637511053698747, -3.3746001369769957e-6, 2.8163466313745392e-6] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.542260056838586, 0.00033327761172990305, -0.00027822803717126405] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.54235764872683, 0.00033316608351502614, -0.0002782890764011935] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.687584823332436, 0.001655829977966438, -0.0013825887368813059] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.688266708309156, 0.0016550238100658532, -0.0013829878528847806] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [13.810831277596188, 0.31523877591740934] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(growth~meals),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [13.649256827507362, 0.31664824141234876] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(growth~meals),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [12.313991611568767, 0.09081027644029858] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.054897182273514, 0.119240397293647] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.07516798436631947, -0.0002901194904550819] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Normal(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.09425628209020485, -0.0004983617414145051] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Normal(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [2.5564115357446133, 0.005121203237218301] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Normal(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [2.2952503238004067, 0.007883880687574145] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Normal(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.553454726833889, 0.010778623695115919] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Normal(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.0839429944953802, 0.015466993740961069] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Normal(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [11.613082768340933, 0.10746041110244206] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Gamma(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.354958120506199, 0.11118815542811016] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Gamma(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.07698785070208947, -0.00032269454273518557] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Gamma(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.0950437919713874, -0.0005140366435911466] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Gamma(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [2.520096325684129, 0.005861604838146036] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Gamma(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [2.2959818816454685, 0.007865038609518699] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Gamma(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.473276551733995, 0.012536527872094224] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Gamma(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.1043143405178815, 0.01495305289832234] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Gamma(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [12.009772837768539, 0.09757144487679642] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Poisson(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.204694341305, 0.1154006782615505] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Poisson(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.07596610583667521, -0.00030382549601011566] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.09462468372609342, -0.0005052680185240707] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [2.540525387358834, 0.005428890150162429] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [2.2950326950962068, 0.007888508617627793] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.5185108949759534, 0.011501587664729895] - apistrat.wts = genWts(apistrat,1) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.092836130186039, 0.015254388159506728] - apistrat.wts = genWts(apistrat,2) - sd = svydesign(data=apistrat,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [33.6084325542773, -0.005375331118991334] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [39.679599355676565, -0.04561035202548833] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.02945276542256612, 5.995955704110545e-6] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.015021044227906446, 0.00011440794587404663] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.5193936176789187, -0.0001790606328083504] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.8681535659383326, -0.0023837748207087782] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [5.803681258038523, -0.0004902638734413207] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [6.607322649924698, -0.0053450503387416125] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.029657127683424198, 5.464318207618596e-6] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.012285606750201688, 0.0001271015752112013] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.5135990865026363, -0.00016458932419128794] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.736710402810912, -0.0021066177395622355] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [5.788166924530174, -0.00045224932753613217] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.029560537101786844, 5.711668214418187e-6] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.013651766694489979, 0.00012206005677876978] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.516356929247808, -0.00017138229820086258] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.8560252928224816, -0.002352665272706395] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [5.795565626257508, -0.00047014410019493165] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [6.3575094132992565, -0.0047912177647261986] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [12.241096093599847, -0.019335518503428153, 0.6511483736914527, 0.019335321551248532] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(growth~meals+ell+dnum),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [11.917944512214367, -0.15179706694825687, 0.7855459974961514, 0.029589449113631185] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(growth~meals+ell+dnum),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [21.005257946639674, 0.9827280154221911, 0.06440098939709134] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [6.352198003579708, 1.3318472822744785, -0.12269369401616546] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.027792901997291015, -0.00025440519280463784, -1.7283397138803052e-6] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Normal(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.040600875904590064, -0.00047563577795723067, 3.469068821071962e-5] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Normal(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.39544398317539, 0.016829838807992666, 0.0006057859458406191] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Normal(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [2.8943948034882707, 0.02838907699404031, -0.002922129069810155] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Normal(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [5.06194201193598, 0.06626668368659366, 0.0034032563958906505] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Normal(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.5209221799040256, 0.10630810917708457, -0.010756449992144744] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Normal(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [23.951994049494704, 0.8476808238116741, 0.07432666726600255] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Gamma(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.03049724650374688, -0.00031040801568211357, -8.65111043794249e-6] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Gamma(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.04444131725553145, -0.0005992844618031765, 6.729349738539915e-5] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Gamma(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.294386557501264, 0.018892903837350317, 0.0013626494243970657] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Gamma(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [2.909071624307819, 0.026532177379337444, -0.002421939822045831] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Gamma(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.997536119654526, 0.06528667192785456, 0.005541577438365458] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Gamma(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - - rcoef = [22.27893587644691, 0.9249351756659369, 0.07306160622542461] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Poisson(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.02900153868999018, -0.0002787483789302613, -3.714046173347455e-6] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.04250716415868932, -0.0005313757155334976, 5.074691632759612e-5] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.3384521314092517, 0.018019245870671954, 0.0009269057356260197] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [2.8501735319279895, 0.029977873166426503, -0.0031733863598370888] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.9805779782576725, 0.06727363929756108, 0.004599642547288176] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.839737043387107, 0.09072877337608513, -0.008027613943817407] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [96.03716474761991, -0.09308529164420191, -0.036018740662903614, -0.08783495038760757] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [90.3532717624607, -0.16995872954241745, -0.0794215124999823, 0.055299585889525125] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.01033909445120728, 1.4027187847251898e-5, 5.118455742927665e-6, 1.011366679594371e-5] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Normal(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.010945924318483602, 2.6240415050236824e-5, 1.1519196624599343e-5, -7.582469261208434e-6] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Normal(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.568200842570717, -0.001150346080770433, -0.000431346173217862, -0.0009431197751283893] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Normal(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.5091882292109515, -0.002118374127034952, -0.0009608541532906671, 0.0006488825394376682] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Normal(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.808247179639459, -0.005183079743361119, -0.001973243115231818, -0.004551504760315564] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Normal(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.518283909884522, -0.009493008999780896, -0.004372438890290817, 0.002996492388453739] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Normal(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [96.19997349022383, -0.08112540464795365, -0.03337353092988485, -0.09950248367709875] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Gamma(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [89.80111693127198, -0.15240190328661482, -0.07479753872974587, 0.05302799980394516] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Gamma(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.010318813114430597, 1.2636763157210844e-5, 4.828607750973035e-6, 1.1505682120285014e-5] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Gamma(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.01102274615789953, 2.353694021052259e-5, 1.0981916918685278e-5, -7.2900996490785434e-6] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Gamma(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.569981640922017, -0.001020531709309143, -0.0004035288855108989, -0.0010703663583560355] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Gamma(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.502545685860094, -0.0018968912882356078, -0.000909411252982136, 0.000622692170254524] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Gamma(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.816713667155282, -0.00455912936908498, -0.0018374552022377537, -0.005160446307074326] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Gamma(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.487874891964788, -0.0085036465505177, -0.004126936781971756, 0.002874125349800617] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Gamma(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [96.11073916966183, -0.08712890881687134, -0.03470691705860483, -0.09351077144199814] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Poisson(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [90.0662131882656, -0.16069925408764113, -0.07706125540266325, 0.05412225524672907] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Poisson(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.010329599340290442, 1.3335833800504176e-5, 4.975046630228947e-6, 1.0794485882605418e-5] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.010985138174732451, 2.4834231191175263e-5, 1.1253811080152784e-5, -7.432930068597907e-6] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.569020069648979, -0.001085736387074946, -0.0004175679706983154, -0.0010051986301528165] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.505759774633123, -0.0020021728348470834, -0.0009349398788967742, 0.0006353847329247284] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.81210659920837, -0.00487240223582593, -0.0019059456220664262, -0.004848172760120932] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.502525731572922, -0.008972392904489234, -0.004247868830253232, 0.0029332273499584416] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [95.98893837597574, 0.0020960735021797266, 0.0027002245697935108] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [93.79511451018823, -0.0028859724945897707, 0.010208612864784947] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.010407216183973248, -2.189097324293306e-7, -2.680253702441079e-7] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Normal(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.010644049752293425, 2.9653800869084056e-7, -1.0616672483584035e-6] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Normal(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.564754639383343, 2.1426512833542692e-5, 2.6901137258989773e-5] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Normal(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.541949724453297, -2.9255050298881498e-5, 0.0001041169630217841] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Normal(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.798685241433185, 0.00010596871257151579, 0.00013475723968354983] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Normal(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.686834861073187, -0.0001452845945416358, 0.0005154955040020731] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Normal(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [95.97505170833615, 0.0020414618237609003, 0.0027810754686418315] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Gamma(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [93.83095148212792, -0.002931568660409474, 0.010196038392819766] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Gamma(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.01040864645309545, -2.1360808510942383e-7, -2.760009268015381e-7] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Gamma(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.010640528012209534, 3.0125000182228103e-7, -1.0606927910520759e-6] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Gamma(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.56461366896614, 2.0888338015041562e-5, 2.7704246128490972e-5] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Gamma(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.542305473848398, -2.9718669103659766e-5, 0.00010400459642094935] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Gamma(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.797985475185868, 0.00010325737037441503, 0.00013878722330251187] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Gamma(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.688620894219156, -0.00014758393432827202, 0.000514899387580319] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Gamma(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [95.98200837110346, 0.002069019624359683, 0.002740353337861943] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Poisson(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [93.81299654257235, -0.002908827511031285, 0.010202460545207983] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Poisson(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.01040792954438476, -2.1628463658097162e-7, -2.7198221219006557e-7] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.010642294973494804, 2.9889919723739873e-7, -1.0611971339654083e-6] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.564684299216584, 2.1159935740659737e-5, 2.7299716458315883e-5] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.54212705622947, -2.9487270939920538e-5, 0.00010406226398345519] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.798336223309624, 0.00010462626847812745, 0.00013675640549820373] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.68772560477788, -0.00014643673068323988, 0.0005152045531879398] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [21.167610614103, 0.31707690934463073] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(growth~meals),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [39.58046593197649, 0.06449167349644262] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(growth~meals),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [12.090620308879862, 0.1132587027979] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [12.832173614643434, 0.1030861677186473] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [12.021198719138594, 0.1146324215395372] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Poisson(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [12.822689686481139, 0.10336422900945573] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Poisson(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.07378016352719716, -0.00032235346434947635] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.07290466475616456, -0.00031153430201177396] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [2.550081370417678, 0.006252315172727608] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [2.5875951193039626, 0.005752547055895223] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.523547778919546, 0.01349352049543981] - apiclus1.wts = genWts(apiclus1,1) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.615139565189028, 0.012234478433057711] - apiclus1.wts = genWts(apiclus1,2) - sd = svydesign(data=apiclus1,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [27.011269086057517, 0.003306192641207121] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [29.734049229078398, -0.0014070978500682226] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.036938393932331164, -4.063172858912875e-6] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.03362313700719904, 1.6781725044485508e-6] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.2973536965077845, 0.00011602267513052464] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.3924146112829914, -4.8583837257764576e-5] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [5.1986478770345155, 0.00030975374140493827] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [5.45305834407665, -0.00013072506464489688] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.036952194300148095, -4.091077081292087e-6] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.03362979906986281, 1.6426568649014016e-6] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.297120670062287, 0.0001165008822296672] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.392222217346555, -4.7584480751816696e-5] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [5.198235306324976, 0.00031060672832945165] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [5.4525416725551095, -0.00012807537447841126] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [27.00956101107512, 0.003309740103588972] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [29.731254655169913, -0.0013928292921025396] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.0369460820311807, -4.0786483407123045e-6] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.03362649700530502, 1.6601034700652453e-6] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.2972152168020687, 0.0001163056530614922] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.3923176684279333, -4.807598830577685e-5] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [5.198383722058313, 0.0003102978033010925] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [5.45279759934814, -0.00012937662394039547] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [17.60045300824547, 0.4401307103518129, -0.2963239543592062, -0.0064756189409267105] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(growth~meals+ell+dnum),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [23.366592068220267, 0.4579682134496839, -0.38234614425540525, -0.015537773084009322] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(growth~meals+ell+dnum),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [12.56943704549076, 1.2488227303735902, 0.22499718656631085] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [15.782434939170027, 1.109970069168465, 0.27018542019039005] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.03765607369805905, -0.0004569825787645028, -1.3387955820354623e-5] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.037581405798499494, -0.0004552900675052336, -1.3784259156516248e-5] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [2.9888346810284623, 0.025631126141300413, 0.004237770493224647] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.0775582657514304, 0.022813907586810494, 0.004942217095911617] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.9291890330013164, 0.09461176914137939, 0.01881611879880435] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.248852659908276, 0.08141677047777836, 0.02261652011414011] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [95.45020244242633, 0.0005588104719308948, 0.07184620304870884, -0.11719767549448347] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [97.50830278066785, -0.1501915695578266, 0.0991286642643912, -0.0827072162975507] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.010471398184637227, -2.287779697832461e-7, -8.420925532842145e-6, 1.4029146029205811e-5] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Normal(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.010251518986302837, 1.854053299246571e-5, -1.066542888306889e-5, 9.049022040323355e-6] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Normal(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.558880071814485, 1.5518260652862364e-5, 0.0007804112370919204, -0.0012855797672350405] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Normal(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.580125691479207, -0.0016634598760997912, 0.0010334112706571997, -0.0008715288619721561] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Normal(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.770566179984149, 5.416592024853234e-5, 0.0037471027054722106, -0.006141369936511209] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Normal(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.875087326434127, -0.007895923900889997, 0.005067185630205803, -0.004253232806394067] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Normal(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [95.3892306817184, -0.004280987804747864, 0.07030851863109519, -0.11254932145244756] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Gamma(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [97.42516077240988, -0.14545429897550416, 0.10006778624810236, -0.081908878573277] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Gamma(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.010476821402501613, 2.2773503010888788e-7, -8.35416858377594e-6, 1.364078865343635e-5] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Gamma(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.010260143273310084, 1.7561848408119722e-5, -1.1021537892053484e-5, 9.310316299235217e-6] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Gamma(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.558306030733681, -3.1860165672387765e-5, 0.0007688343373653318, -0.0012422786842810297] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Gamma(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.579285034632418, -0.0015912476819282507, 0.0010557053716828647, -0.0008803886849389655] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Gamma(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.76760761267514, -0.00018578256887032257, 0.003678952020250707, -0.005916020059623148] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Gamma(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.870913679440493, -0.007597670411100468, 0.005145852613271459, -0.0042547667579689455] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Gamma(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [95.4207615435587, -0.001796603812703806, 0.07111531283749667, -0.11493754633496521] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Poisson(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [97.4676223518129, -0.14745592552744205, 0.0997382122143591, -0.08253357296678293] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Poisson(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.010474021950952017, -9.28166853800264e-9, -8.393704840108051e-6, 1.3844484817938591e-5] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.010255782012722644, 1.8014865770872542e-5, -1.0856010322072613e-5, 9.199778637433041e-6] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.558603013451985, -7.402002325582733e-6, 0.0007751389155923243, -0.001264742977847314] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.579712287270685, -0.0016236520821931214, 0.0010459610650708719, -0.0008781884977022112] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.769138499888868, -6.224448665522241e-5, 0.0037152964599643182, -0.006032367083667221] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.87304041592248, -0.007728297155601642, 0.005113619975757302, -0.004265337911943751] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(full~ell+growth+meals),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [14.53576103415235, 0.3202608744513976] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(growth~meals),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [25.020794198515205, 0.23479025486014074] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(growth~meals),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [12.76310479339429, 0.09708826730389312] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [11.630532315251829, 0.10632063404240416] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.07441320727925439, -0.0003184272913034482] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Normal(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.07953385436246366, -0.00037495300789860676] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Normal(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [2.5706849654989345, 0.005672599554745491] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Normal(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [2.4934549317488672, 0.006499278608911069] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Normal(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.592865175451691, 0.011798948393309843] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Normal(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.4435763394989913, 0.013274331461837475] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Normal(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [12.881903718693145, 0.09362179374541775] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Gamma(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [11.991554948704405, 0.09418865315916253] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Gamma(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.07517164786272039, -0.00033352515401263743] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Gamma(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.07984623288093695, -0.00038294500546953535] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Gamma(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [2.5690688537524875, 0.005707914917914804] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Gamma(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [2.501637438666087, 0.006249588711716856] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Gamma(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.599745092745747, 0.011612791584109484] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Gamma(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.4755712426357372, 0.012250438779389793] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Gamma(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [12.818133337538974, 0.09559845282984288] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Poisson(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [11.804674184014628, 0.1006762631091293] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Poisson(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.07484165355326194, -0.00032655750027481777] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.07974585105130189, -0.0003801983903780943] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [2.5690292767409866, 0.00570933814907725] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [2.496142121256949, 0.006421623745290465] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.595079837431775, 0.01174392442459022] - apiclus2.wts = genWts(apiclus2,1) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.4573083909377567, 0.012853608051332578] - apiclus2.wts = genWts(apiclus2,2) - sd = svydesign(data=apiclus2,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [22.23439009133479, 0.012149721254178248] - apipop.wts = genWts(apipop,1) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [19.379529133476442, 0.01742166488113855] - apipop.wts = genWts(apipop,2) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.044591975674265094, -1.800163600556982e-5] - apipop.wts = genWts(apipop,1) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.04856811075892835, -2.4546945309033293e-5] - apipop.wts = genWts(apipop,2) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.104137912347123, 0.00047214256447071723] - apipop.wts = genWts(apipop,1) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [2.995316601369728, 0.0006621119744598153] - apipop.wts = genWts(apipop,2) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.717518227886106, 0.0011993560271005886] - apipop.wts = genWts(apipop,1) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.437798064124778, 0.001701603911033396] - apipop.wts = genWts(apipop,2) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [22.762355853199104, 0.010890816534124902] - apipop.wts = genWts(apipop,1) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [19.184153968994853, 0.01792487652396429] - apipop.wts = genWts(apipop,2) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.04379312473822179, -1.6349701448620263e-5] - apipop.wts = genWts(apipop,1) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.04886103152540363, -2.513777150904968e-5] - apipop.wts = genWts(apipop,2) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.126199258968172, 0.00042322960070233684] - apipop.wts = genWts(apipop,1) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [2.9911201104493093, 0.0006717156962124961] - apipop.wts = genWts(apipop,2) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.772150587890555, 0.001073798269964286] - apipop.wts = genWts(apipop,1) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.424993204844972, 0.0017327856012052977] - apipop.wts = genWts(apipop,2) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [22.52319890293024, 0.011476311849907966] - apipop.wts = genWts(apipop,1) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [19.331830023124613, 0.017536405545714867] - apipop.wts = genWts(apipop,2) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.04417974224219384, -1.7169152638988327e-5] - apipop.wts = genWts(apipop,1) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.04870474977209929, -2.4813549674945644e-5] - apipop.wts = genWts(apipop,2) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.1160082436544827, 0.0004464195161463967] - apipop.wts = genWts(apipop,1) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [2.994332284456228, 0.0006642503312315288] - apipop.wts = genWts(apipop,2) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.7472404972029, 0.0011325722413844626] - apipop.wts = genWts(apipop,1) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.43527072494797, 0.001707409990271073] - apipop.wts = genWts(apipop,2) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [23.27674238741806, 0.2265038812633155, -0.00983094170459761, -0.002640596086200574] - apipop.wts = genWts(apipop,1) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(growth~meals+ell+dnum),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [25.370353249623857, 0.20538405354664227, 0.007899341624985498, -0.005378728200668248] - apipop.wts = genWts(apipop,2) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(growth~meals+ell+dnum),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [20.221785634008462, 1.0685971201271238, 0.1027513982394961] - apipop.wts = genWts(apipop,1) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [19.27922518039677, 1.0828137394397237, 0.09586409553223135] - apipop.wts = genWts(apipop,2) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.02921238961457218, -0.00024953358486822603, -2.0100348861915888e-5] - apipop.wts = genWts(apipop,1) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.029860332287800304, -0.0002584426611375959, -1.895714110208761e-5] - apipop.wts = genWts(apipop,2) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.2861593680350025, 0.018434020710060492, 0.002082857777291423] - apipop.wts = genWts(apipop,1) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.2578407478391203, 0.018896914028401737, 0.0019930414448396305] - apipop.wts = genWts(apipop,2) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.710235887248018, 0.07597655522651227, 0.008326425900809984] - apipop.wts = genWts(apipop,1) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.6368022388746715, 0.07711454809242572, 0.007855745686824159] - apipop.wts = genWts(apipop,2) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [22.166081452544613, 0.22136802130499902] - apipop.wts = genWts(apipop,1) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(growth~meals),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [23.217489630436187, 0.20742473434125808] - apipop.wts = genWts(apipop,2) - sd = svydesign(data=apipop,weights=:wts) - jcoef = svyglm(@formula(growth~meals),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [21.650497884802157, 0.010062258530350853] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [15.35079371648053, 0.010789008933764766] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.04635672714581871, -1.734629735040647e-5] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.0756971167146116, -5.52124732894222e-5] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.071106349570877, 0.00042181523025128627] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [2.605024332066764, 0.0008751857555967217] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.647488946818263, 0.0010317473465799375] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.80601783995445, 0.001522967962829505] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Normal(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [22.27052208954838, 0.008619250571327296] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.04525938933001788, -1.508817945336351e-5] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.06199168147294855, -2.702923670413523e-5] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.0985569829286574, 0.0003617521980159552] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [2.82683644686873, 0.0003875904288854227] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.713276960829454, 0.0008832981950201349] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.146048442631385, 0.0007531548099478192] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Gamma(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [21.98820569449091, 0.00929075174026864] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [16.714541439893186, 0.0076530171542340646] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.04578484118130555, -1.6194183092118965e-5] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.06826295663324547, -4.045740695109782e-5] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.085876554263612, 0.00039013243354418185] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [2.758193027328452, 0.0005431593616227393] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.68319697148335, 0.0009527027756396931] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.03634850583009, 0.0010069046509206749] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(cnum~dnum),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [24.82758069617155, 0.29286700952579436, -0.03691707470178347, -0.015295946199160597] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(growth~meals+ell+dnum),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [15.918845413841463, 0.47359734194085656, -0.315499444817528, -0.015108840993373984] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(growth~meals+ell+dnum),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [22.66672036065673, 0.935692198094081, 0.1592001186738126] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [18.024577819737416, 0.8510129277921251, 0.24725334660838746] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [18.099336966787398, 0.8609168623565143, 0.23219125882826197] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Poisson(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.028499566477846162, -0.00021378711667711803, -5.356345419722975e-5] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.029760993395621924, -0.00022181936917489264, -6.328125974334514e-5] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.358908660414955, 0.015906123249569144, 0.003126673407781379] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.2833974839712683, 0.01581238650376551, 0.004084526918153274] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.962994355580012, 0.06610488008628344, 0.011662470868456358] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.713235440661724, 0.061550100222985245, 0.01637522740421099] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(meals~ell+growth),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [95.21202843265861, 0.017003718835692697, -0.012776703721702535] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [96.50687250475366, 0.012448891628779435, -0.009830114718549651] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.010491164829391124, -1.7415002720766583e-6, 1.306050744614016e-6] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Normal(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.010358454770769356, -1.277530313276944e-6, 1.0066062453698471e-6] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Normal(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.5566742056722935, 0.0001721078418155973, -0.00012920286803554638] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Normal(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.569785116214157, 0.00012611837038504933, -9.948212900254918e-5] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Normal(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.759063360836693, 0.0008553805085390379, -0.0006424458516566352] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Normal(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.82421324184559, 0.0006265148456543873, -0.0004944595873635184] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Normal(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [95.19897220467752, 0.017106586862824846, -0.012863906853071563] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Gamma(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [96.50754129204527, 0.012484997800712633, -0.009869011930619832] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Gamma(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.01049258733947373, -1.753113256391414e-6, 1.315989194899257e-6] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Gamma(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.010358391909760878, -1.281581353760541e-6, 1.010948644379981e-6] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Gamma(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.556537767977723, 0.00017320330673554853, -0.0001301362310643327] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Gamma(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.56979160769325, 0.00012650126864536293, -9.98935468578143e-5] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Gamma(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.758395738260907, 0.0008606921127718869, -0.0006469604135616767] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Gamma(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.824246196624872, 0.0006283744981996406, -0.0004964603113298763] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Gamma(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [95.20548474291016, 0.01705560394852365, -0.012820763202976856] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Poisson(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [96.50720669487956, 0.012467110352926788, -0.009849736292381363] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Poisson(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.010491877069819897, -1.7473443873234547e-6, 1.3110589060458006e-6] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.010358423358912553, -1.2795719372383925e-6, 1.0087942838682328e-6] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.5566058616573315, 0.0001726597035805286, -0.00012967378159337847] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [4.56978835984881, 0.00012631145729649884, -9.968954901681832e-5] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.758728938646206, 0.0008580571492429526, -0.0006447245061143383] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [9.82422970863238, 0.0006274529139253154, -0.0004954685576497072] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(pcttest~api00+api99),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [17.618057537052124, 0.2855817329123753] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(growth~meals),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [7.055586336420611, 0.3385660012670978] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(growth~meals),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [14.452546232206487, 0.06013704794628077] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [7.208916670262022, 0.09647309140449355] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Normal(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [14.551824595499689, 0.05815187771446325] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Poisson(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [8.71612535007701, 0.06441336951270263] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Poisson(),IdentityLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.06818079188235583, -0.0002048934351749668] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [0.12917734532987823, -0.0008410752426298812] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Poisson(),InverseLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [2.682504168784929, 0.003441095906950557] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [2.078774738912212, 0.0078697423218258] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Poisson(),LogLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [3.8197222674497087, 0.0070671224153743325] - apisrs.wts = genWts(apisrs,1) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - - rcoef = [2.892120154959868, 0.011198654714211375] - apisrs.wts = genWts(apisrs,2) - sd = svydesign(data=apisrs,weights=:wts) - jcoef = svyglm(@formula(mobility~meals),sd,Poisson(),SqrtLink()).coefficients - @test maximum(abs.((jcoef-rcoef)./rcoef))<=rtol - -end diff --git a/test/svyhist.jl b/test/svyhist.jl deleted file mode 100644 index cc2e8246..00000000 --- a/test/svyhist.jl +++ /dev/null @@ -1,22 +0,0 @@ -using Survey -using Test - -@testset "svyhist.jl" begin - apistrat = load_data("apistrat") - dstrat = svydesign(data = apistrat, id = :1, strata = :stype, weights = :pw, fpc = :fpc) - - @test Survey.sturges(10) == 5 - @test Survey.sturges([1, 2, 5, 10, 15, 17, 20]) == 4 - - h = svyhist(dstrat, :enroll) - @test h.grid[1].entries[1].positional[2] |> length == 16 - - h = svyhist(dstrat, :enroll, 9) - @test h.grid[1].entries[1].positional[2] |> length == 7 - - h = svyhist(dstrat, :enroll, Survey.sturges) - @test h.grid[1].entries[1].positional[2] |> length == 7 - - h = svyhist(dstrat, :enroll, [0, 1000, 2000, 3000]) - @test h.grid[1].entries[1].positional[2] |> length == 3 -end diff --git a/test/svyplot.jl b/test/svyplot.jl deleted file mode 100644 index e9fb4415..00000000 --- a/test/svyplot.jl +++ /dev/null @@ -1,12 +0,0 @@ -using Survey -using Test - -@testset "svyplot.jl" begin - apistrat = load_data("apistrat") - dstrat = svydesign(data = apistrat, id = :1, strata = :stype, weights = :pw, fpc = :fpc) - s = svyplot(dstrat, :api99, :api00) - - @test s.grid[1].entries[1].named[:markersize] == dstrat.variables.weights - @test s.grid[1].entries[1].positional[1] == dstrat.variables.api99 - @test s.grid[1].entries[1].positional[2] == dstrat.variables.api00 -end diff --git a/test/total.jl b/test/total.jl new file mode 100644 index 00000000..80e36cba --- /dev/null +++ b/test/total.jl @@ -0,0 +1,94 @@ +@testset "total_SimpleRandomSample" begin + apisrs_original = load_data("apisrs") + + # base functionality + apisrs = copy(apisrs_original) + srs = SimpleRandomSample(apisrs; popsize = :fpc) + tot = total(:api00, srs) + @test tot.total[1] ≈ 4.06688749e6 atol = 1e-4 + @test tot.SE[1] ≈ 57292.7783113177 atol = 1e-4 + # without fpc + srs_ignorefpc = SimpleRandomSample(apisrs; popsize = :fpc, ignorefpc = true) + tot = total(:api00, srs_ignorefpc) + # TODO: uncomment after correcting `total` function + # @test tot.total[1] ≈ 131317 atol = 1 + # @test tot.SE[1] ≈ 1880.6 atol = 1e-1 + + # CategoricalArray + apisrs = copy(apisrs_original) + apisrs[!, :cname] = CategoricalArrays.categorical(apisrs.cname) + srs = SimpleRandomSample(apisrs; popsize = :fpc) + tot = total(:cname, srs) + @test size(tot)[1] == apisrs.cname |> unique |> length + @test filter(:cname => ==("Alameda"), tot).total[1] ≈ 340.67 atol = 1e-2 + @test filter(:cname => ==("Alameda"), tot).SE[1] ≈ 98.472 atol = 1e-3 + @test filter(:cname => ==("Los Angeles"), tot).total[1] ≈ 1393.65 atol = 1e-2 + @test filter(:cname => ==("Los Angeles"), tot).SE[1] ≈ 180.368 atol = 1e-3 + + # Vector{Symbol} + apisrs = copy(apisrs_original) + srs = SimpleRandomSample(apisrs; popsize = :fpc) + tot = total([:api00, :enroll], srs) + ## :api00 + @test tot.total[1] ≈ 4066888 atol = 1 + @test tot.SE[1] ≈ 57293 atol = 1 + ## :enroll + @test tot.total[2] ≈ 3621074 atol = 1 + @test tot.SE[2] ≈ 169520 atol = 1 + + # subpopulation + apisrs = copy(apisrs_original) + srs = SimpleRandomSample(apisrs; popsize = :fpc) + tot = total(:api00, :cname, srs) + @test size(tot)[1] == apisrs.cname |> unique |> length + @test filter(:cname => ==("Los Angeles"), tot).total[1] ≈ 917238.49 atol = 1e-2 + @test filter(:cname => ==("Los Angeles"), tot).SE[1] ≈ 122289.00 atol = 1e-2 + @test filter(:cname => ==("Monterey"), tot).total[1] ≈ 74947.40 atol = 1e-2 + @test filter(:cname => ==("Monterey"), tot).SE[1] ≈ 37616.17 atol = 1e-2 +end + +@testset "total_Stratified" begin + apistrat_original = load_data("apistrat") + + # base functionality + apistrat = copy(apistrat_original) + strat = StratifiedSample(apistrat, :stype; popsize = :fpc) + tot = total(:api00, strat) + @test tot.total[1] ≈ 4102208 atol = 1e-1 + @test tot.SE[1] ≈ 58279 atol = 1e-1 + # without fpc + apistrat = copy(apistrat_original) + strat_ignorefpc = StratifiedSample(apistrat, :stype; popsize = :fpc, ignorefpc = true) + tot = total(:api00, strat_ignorefpc) + @test tot.total[1] ≈ 130564 atol = 1e-4 + # TODO: uncomment after correcting `total` function + # @test tot.SE[1] ≈ 1690.4 atol = 1e-1 + + # CategoricalArray + apistrat = copy(apistrat_original) + apistrat[!, :cname] = CategoricalArrays.categorical(apistrat.cname) + strat = StratifiedSample(apistrat, :stype; popsize = :fpc) + # TODO: uncomment after adding `CategoricalArray` support + # @test tot.SE[1] ≈ 1690.4 atol = 1e-1 + # tot = total(:cname, strat) + # @test size(tot)[1] == apistrat.cname |> unique |> length + # @test filter(:cname => ==("Kern"), tot).total[1] ≈ 291.97 atol = 1e-2 + # @test filter(:cname => ==("Kern"), tot).SE[1] ≈ 101.760 atol = 1e-3 + # @test filter(:cname => ==("Los Angeles"), tot).total[1] ≈ 1373.15 atol = 1e-2 + # @test filter(:cname => ==("Los Angeles"), tot).SE[1] ≈ 199.635 atol = 1e-3 + + # Vector{Symbol} + apistrat = copy(apistrat_original) + strat = StratifiedSample(apistrat, :stype; popsize = :fpc) + tot = total([:api00, :enroll], strat) + ## :api00 + @test tot.total[1] ≈ 4102208 atol = 1 + @test tot.SE[1] ≈ 58279 atol = 1 + ## :enroll + @test tot.total[2] ≈ 3687178 atol = 1 + @test tot.SE[2] ≈ 114642 atol = 1 + + # subpopulation + # TODO: add functionality in `src/total.jl` + # TODO: add tests +end