update stderr -> stderror

Author: matthieugomez

REQUIRE

@@ -1,6 +1,6 @@
 julia 0.6.0
 Distributions 0.4.6
-StatsBase 0.7.1
+StatsBase 0.22.0
 StatsModels 0.2.4
 DataArrays 0.3.4
 CategoricalArrays

benchmark/benchmark.jl

@@ -5,21 +5,21 @@ id1 = rand(1:(N/K), N)
 id2 = rand(1:K, N)
 x1 =  randn(N)
 x2 =  randn(N)
-w = cos(id1)
-y= 3 .* x1 .+ 5 .* x2 .+ cos(id1) .+ cos(id2).^2 .+ randn(N)
+w = cos.(id1)
+y= 3 .* x1 .+ 5 .* x2 .+ cos.(id1) .+ cos.(id2).^2 .+ randn(N)
 df = DataFrame(id1 = pool(id1), id2 = pool(id2), x1 = x1, x2 = x2, w = w, y = y)
-@time reg(y ~ x1 + x2, df)
+@time reg(df, @model(y ~ x1 + x2))
 # 1.258554 seconds (723 allocations: 1.205 GB, 18.70% gc time)
 #852MB to obtain the matrix etc, and then 1.205 for Regressions part
-@time reg(y ~ x1 + x2, df, VcovCluster(:id2))
+@time reg(df, @model(y ~ x1 + x2, vcov = cluster(id2)))
 #  1.569679 seconds (843 allocations: 1.293 GB, 25.21% gc time)
-@time reg(y ~ x1 + x2 |> id1, df)
+@time reg(df, @model(y ~ x1 + x2, fe = id1))
 # 1.476390 seconds (890 allocations: 1.175 GB, 20.15% gc time)
-@time reg(y ~ x1 + x2 |> id1, df, VcovCluster(:id1))
+@time reg(df, @model(y ~ x1 + x2, fe = id1, vcov = cluster(id1)))
 # 1.974738 seconds (1.04 k allocations: 1.255 GB, 15.85% gc time)
-@time reg(y ~ x1 + x2 |> id1 + id2, df)
+@time reg(df, @model(y ~ x1 + x2, fe = id1 + id2))
 # 4.554836 seconds (1.01 k allocations: 1.188 GB, 9.56% gc time)
-@time reg(y ~ x1 + x2 |> id1, df, weight = :w)
+@time reg(df, @model(y ~ x1 + x2, fe = id1, weights = w))
 # 2.000850 seconds (20.00 M allocations: 1.010 GB, 18.82% gc time)
-@time reg(y ~ x1 + x2 |> id1 + id2, df, weight = :w)
+@time reg(df, @model(y ~ x1 + x2, fe = id1 + id2, weights = w))
 # 4.118974 seconds (20.00 M allocations: 1.018 GB, 9.60% gc time)

src/FixedEffectModels.jl

@@ -15,7 +15,7 @@ import CategoricalArrays: CategoricalArray, CategoricalVector, compress, categor
 import DataFrames: DataFrame, AbstractDataFrame, completecases, names!, ismissing
 import StatsModels: ModelMatrix, ModelFrame, Terms, coefnames, Formula, completecases, names!,  @formula
 using Reexport
-import StatsBase: coef, nobs, coeftable, vcov, predict, residuals, var, RegressionModel, model_response, stderr, confint, fit, CoefTable, df_residual
+import StatsBase: coef, nobs, coeftable, vcov, predict, residuals, var, RegressionModel, model_response, stderror, confint, fit, CoefTable, df_residual
 @reexport using StatsBase
 ##############################################################################
 ##

src/RegressionResult.jl

@@ -16,7 +16,7 @@ df_residual(x::AbstractRegressionResult) = x.df_residual
 
 function confint(x::AbstractRegressionResult) 
     scale = quantile(TDist(x.df_residual), 1 - (1-0.95)/2)
-    se = stderr(x)
+    se = stderror(x)
     hcat(x.coef -  scale * se, x.coef + scale * se)
 end
 
@@ -66,7 +66,7 @@ function coeftable(x::AbstractRegressionResult)
     ctitle = title(x)
     ctop = top(x)
     cc = coef(x)
-    se = stderr(x)
+    se = stderror(x)
     coefnms = coefnames(x)
     conf_int = confint(x)
     # put (intercept) last

src/fixedeffect/FixedEffectProblem_LSMR.jl

@@ -113,13 +113,13 @@ size(fem::FixedEffectMatrix, dim::Integer) = (dim == 1) ? fem.m :
 function A_mul_B_helper!(α::Number, fe::FixedEffect, 
                         x::Vector{Float64}, y::AbstractVector{Float64}, cache::Vector{Float64})
     for (i, j) in zip(1:length(y), eachindex(y))
-        y[j] += α * x[fe.refs[i]] * cache[i]
+        @inbounds y[j] += α * x[fe.refs[i]] * cache[i]
     end
 end
 function A_mul_B!(α::Number, fem::FixedEffectMatrix, fev::FixedEffectVector, 
                 β::Number, y::AbstractVector{Float64})
     safe_scale!(y, β)
-    for i in 1:length(fev._)
+    Threads.@threads for i in 1:length(fev._)
         A_mul_B_helper!(α, fem._[i], fev._[i], y, fem.cache[i])
     end
     return y
@@ -129,14 +129,13 @@ end
 function Ac_mul_B_helper!(α::Number, fe::FixedEffect, 
                         y::AbstractVector{Float64}, x::Vector{Float64}, cache::Vector{Float64})
     for (i, j) in zip(1:length(y), eachindex(y))
-        x[fe.refs[i]] += α * y[j] * cache[i]
+        @inbounds x[fe.refs[i]] += α * y[j] * cache[i]
     end
 end
 function Ac_mul_B!(α::Number, fem::FixedEffectMatrix, 
                 y::AbstractVector{Float64}, β::Number, fev::FixedEffectVector)
-   safe_scale!(fev, β)
-    for i in 1:length(fev._)
-       # @code_warntype Ac_mul_B_helper!(α, fem._[i], y, fev._[i])
+    safe_scale!(fev, β)
+    Threads.@threads for i in 1:length(fev._)
         Ac_mul_B_helper!(α, fem._[i], y, fev._[i], fem.cache[i])
     end
     return fev

src/reg.jl

@@ -125,7 +125,11 @@ function reg(df::AbstractDataFrame, f::Formula;
     sqrtw = get_weights(df, esample, weights)
 
     # remove unusused levels
-    subdf = df[esample, all_vars]
+    if any(esample)
+        subdf = df[esample, all_vars]
+    else
+        subdf = df[all_vars]
+    end
     main_vars = unique(convert(Vector{Symbol}, vcat(vars, endo_vars, iv_vars)))
     for v in main_vars
         # in case subdataframe, don't construct subdf[v] if you dont need to do it

test/reg.jl

@@ -186,55 +186,55 @@ x = reg(df, m)
 # Simple
 m = @model y ~ x1
 x = reg(df, m)
-@test stderr(x) ≈ [1.52126, 0.01889] atol = 1e-4
+@test stderror(x) ≈ [1.52126, 0.01889] atol = 1e-4
 # Stata ivreg
 m = @model y ~ (x1 ~ z1)
 x = reg(df, m)
-@test stderr(x) ≈ [1.53661, 0.01915] atol = 1e-4
+@test stderror(x) ≈ [1.53661, 0.01915] atol = 1e-4
 # Stata areg
 m = @model y ~ x1 fe = pid1
 x = reg(df, m)
-@test stderr(x) ≈  [0.00980] atol = 1e-4
+@test stderror(x) ≈  [0.00980] atol = 1e-4
 
 # White
 # Stata reg
 m = @model y ~ x1 vcov = robust
 x = reg(df, m)
-@test stderr(x) ≈ [1.68679, 0.01670] atol = 1e-4
+@test stderror(x) ≈ [1.68679, 0.01670] atol = 1e-4
 # Stata ivreg
 m = @model y ~ (x1 ~ z1) vcov = robust
 x = reg(df, m)
-@test stderr(x) ≈ [1.63305, 0.01674] atol = 1e-4
+@test stderror(x) ≈ [1.63305, 0.01674] atol = 1e-4
 # Stata areg
 m = @model y ~ x1 fe = pid1 vcov = robust
 x = reg(df, m)
-@test stderr(x) ≈ [0.01100] atol = 1e-4
+@test stderror(x) ≈ [0.01100] atol = 1e-4
 
 # cluster
 m = @model y ~ x1 vcov = cluster(pid1)
 x = reg(df, m)
-@test stderr(x)[2] ≈ 0.03792 atol = 1e-4
+@test stderror(x)[2] ≈ 0.03792 atol = 1e-4
 m = @model y ~ x1 fe = pid1 vcov = cluster(pid2)
 x = reg(df, m)
-@test stderr(x) ≈ [0.02205] atol = 1e-4
+@test stderror(x) ≈ [0.02205] atol = 1e-4
 # stata reghxe
 m = @model y ~ x1 fe = pid1 vcov = cluster(pid1)
 x = reg(df, m)
-@test stderr(x) ≈ [0.03573] atol = 1e-4
+@test stderror(x) ≈ [0.03573] atol = 1e-4
 
 # no reference
 m = @model y ~ x1 vcov = cluster(pid1 + pid2)
 x = reg(df, m)
-@test stderr(x)[1] ≈ 6.17025 atol = 1e-4
+@test stderror(x)[1] ≈ 6.17025 atol = 1e-4
 # no reference
 m = @model y ~ x1 fe = pid1 vcov = cluster(pid1 + pid2)
 x = reg(df, m)
-@test stderr(x)[1] ≈ 0.04037 atol = 1e-4
+@test stderror(x)[1] ≈ 0.04037 atol = 1e-4
 
 # TO CHECK WITH STATA
 m = @model y ~ x1 fe = pid1 vcov = cluster(pid1&pid2)
 x = reg(df, m)
-@test stderr(x) ≈ [0.0110005] atol = 1e-4
+@test stderror(x) ≈ [0.0110005] atol = 1e-4
 
 
 @test_throws ErrorException reg(df, @model(y ~ x1, vcov = cluster(State)))