update (#34)

Author: matthieugomez

Project.toml

@@ -1,10 +1,9 @@
 name = "InteractiveFixedEffectModels"
 uuid = "80307280-efb2-5c5d-af8b-a9c15821677b"
-version = "1.0.1"
+version = "1.1.0"
 
 [deps]
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
-Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 FixedEffectModels = "9d5cd8c9-2029-5cab-9928-427838db53e3"
 FixedEffects = "c8885935-8500-56a7-9867-7708b20db0eb"
@@ -14,22 +13,23 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
 StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 Vcov = "ec2bfdc2-55df-4fc9-b9ae-4958c2cf2486"
 
 [compat]
 DataFrames = "0.21"
-Distributions = "0.23"
 FillArrays = "0.7, 0.8, 0.9"
-FixedEffectModels = "1.2"
-FixedEffects = "1.1"
+FixedEffectModels = "1.3"
+FixedEffects = "2"
 LeastSquaresOptim = "0.7"
 Reexport = "0.2"
 StatsBase = "0.33"
 StatsModels = "0.6"
+StatsFuns = "0.9"
 Tables = "1"
-Vcov = "0.3"
+Vcov = "0.4"
 julia = "1.3"
 
 [extras]

src/InteractiveFixedEffectModels.jl

@@ -5,14 +5,14 @@ module InteractiveFixedEffectModels
 ##
 ##############################################################################
 using DataFrames
-using Distributions
 using FillArrays
 using FixedEffects
 using LeastSquaresOptim
 using LinearAlgebra
 using Printf
 using Statistics
 using StatsBase
+using StatsFuns
 using StatsModels
 using Tables
 using Vcov
@@ -36,6 +36,9 @@ regife
 ##
 ##############################################################################
 include("types.jl")
+include("utils/tss.jl")
+include("utils/formula.jl")
+
 include("methods/gauss_seidel.jl")
 include("methods/ls.jl")
 include("fit.jl")

src/fit.jl

@@ -19,13 +19,12 @@ function regife(
     ##
     ##############################################################################
     df = DataFrame(df; copycols = false)
-    if  (ConstantTerm(0) ∉ FixedEffectModels.eachterm(f.rhs)) & (ConstantTerm(1) ∉ FixedEffectModels.eachterm(f.rhs))
-        formula = FormulaTerm(f.lhs, tuple(ConstantTerm(1), FixedEffectModels.eachterm(f.rhs)...))
+    if  (ConstantTerm(0) ∉ eachterm(f.rhs)) & (ConstantTerm(1) ∉ eachterm(f.rhs))
+        formula = FormulaTerm(f.lhs, tuple(ConstantTerm(1), eachterm(f.rhs)...))
     end
 
-    formula, formula_endo, formula_iv = FixedEffectModels.parse_iv(f)
 
-    m, formula = parse_interactivefixedeffect(df, formula)
+    m, formula = parse_interactivefixedeffect(df, f)
     has_weights = (weights != nothing)
 
 
@@ -60,10 +59,10 @@ function regife(
          weights = uweights(sum(esample))
          sqrtw = ones(length(weights))
      end
-    for a in FixedEffectModels.eachterm(formula.rhs)
+    for a in eachterm(formula.rhs)
        if has_fe(a)
            isa(a, InteractionTerm) && error("Fixed effects cannot be interacted")
-           Symbol(FixedEffectModels.fesymbol(a)) ∉ factor_vars && error("FixedEffect should correspond to id or time dimension of the factor model")
+           Symbol(fesymbol(a)) ∉ factor_vars && error("FixedEffect should correspond to id or time dimension of the factor model")
        end
     end
     fes, ids, formula = FixedEffectModels.parse_fixedeffect(df, formula)
@@ -72,15 +71,15 @@ function regife(
     ## Compute factors, an array of AbtractFixedEffects
     if has_fes
         if any([isa(fe.interaction, Ones) for fe in fes])
-                formula = FormulaTerm(formula.lhs, tuple(ConstantTerm(0), (t for t in FixedEffectModels.eachterm(formula.rhs) if t!= ConstantTerm(1))...))
+                formula = FormulaTerm(formula.lhs, tuple(ConstantTerm(0), (t for t in eachterm(formula.rhs) if t!= ConstantTerm(1))...))
                 has_fes_intercept = true
         end
-        fes = FixedEffect[FixedEffectModels._subset(fe, esample) for fe in fes]
-        feM = FixedEffectModels.AbstractFixedEffectSolver{Float64}(fes, weights, Val{:cpu})
+        fes = FixedEffect[fe[esample] for fe in fes]
+        feM = AbstractFixedEffectSolver{Float64}(fes, weights, Val{:cpu})
     end
 
 
-    has_intercept = ConstantTerm(1) ∈ FixedEffectModels.eachterm(formula.rhs)
+    has_intercept = ConstantTerm(1) ∈ eachterm(formula.rhs)
 
 
     iterations = 0
@@ -100,7 +99,7 @@ function regife(
     formula_schema = apply_schema(formula, schema(formula, subdf, contrasts), StatisticalModel)
 
     y = convert(Vector{Float64}, response(formula_schema, subdf))
-    tss_total = FixedEffectModels.tss(y, has_intercept || has_fes_intercept, weights)
+    tss_total = tss(y, has_intercept || has_fes_intercept, weights)
 
     X = convert(Matrix{Float64}, modelmatrix(formula_schema, subdf))
 
@@ -123,8 +122,8 @@ function regife(
 
     # demean variables
     if has_fes
-        FixedEffectModels.solve_residuals!(y, feM)
-        FixedEffectModels.solve_residuals!(X, feM, progress_bar = false)
+        solve_residuals!(y, feM)
+        solve_residuals!(X, feM, progress_bar = false)
     end
 
 
@@ -136,8 +135,8 @@ function regife(
     ##############################################################################
 
     # initialize factor models at 0.1
-    idpool = fill(0.1, length(id.pool), m.rank)
-    timepool = fill(0.1, length(time.pool), m.rank)
+    idpool = fill(0.1, id.n, m.rank)
+    timepool = fill(0.1, time.n, m.rank)
 
     y .= y .* sqrtw 
     X .= X .* sqrtw
@@ -169,12 +168,12 @@ function regife(
             # y ~ x + γ1 x factors + γ2 x loadings
             # if not, this means fit! ended up on a a local minimum. 
             # restart with randomized coefficients, factors, loadings
-            newfeM = FixedEffectModels.AbstractFixedEffectSolver{Float64}(getfactors(fp, fs), weights, Val{:cpu})
+            newfeM = AbstractFixedEffectSolver{Float64}(getfactors(fp, fs), weights, Val{:cpu})
             ym .= ym ./ sqrtw
-            FixedEffectModels.solve_residuals!(ym, newfeM, tol = tol, maxiter = maxiter)
+            solve_residuals!(ym, newfeM, tol = tol, maxiter = maxiter)
             ym .= ym .* sqrtw
             Xm .= Xm ./ sqrtw
-            FixedEffectModels.solve_residuals!(Xm, newfeM, tol = tol, maxiter = maxiter)
+            solve_residuals!(Xm, newfeM, tol = tol, maxiter = maxiter)
             Xm .= Xm .* sqrtw
             ydiff = Xm  * (fs.b - Xm \ ym)
             if iterations >= maxiter || norm(ydiff)  <= 0.01 * norm(y)
@@ -226,7 +225,7 @@ function regife(
         # compute various r2
         nobs = sum(esample)
         rss = sum(abs2, residualsm)
-        _tss = FixedEffectModels.tss(ym ./ sqrtw, has_intercept || has_fes_intercept, weights)
+        _tss = tss(ym ./ sqrtw, has_intercept || has_fes_intercept, weights)
         r2_within = 1 - rss / _tss 
 
         rss = sum(abs2, residuals)
@@ -266,7 +265,7 @@ function regife(
             subtract_factor!(fp, fs)
             axpy!(-1.0, residuals, oldresiduals)
             # get fixed effect
-            newfes, b, c = FixedEffectModels.solve_coefficients!(oldresiduals, feM; tol = tol, maxiter = maxiter)
+            newfes, b, c = solve_coefficients!(oldresiduals, feM; tol = tol, maxiter = maxiter)
             for j in 1:length(fes)
                 augmentdf[!, ids[j]] = Vector{Union{Float64, Missing}}(missing, length(esample))
                 augmentdf[esample, ids[j]] = newfes[j]

src/types.jl

@@ -281,11 +281,36 @@ StatsBase.islinear(x::InteractiveFixedEffectModel) = false
 StatsBase.rss(x::InteractiveFixedEffectModel) = x.rss
 StatsBase.predict(::InteractiveFixedEffectModel, ::AbstractDataFrame) = error("predict is not defined for linear factor models. Use the option save = true")
 StatsBase.residuals(::InteractiveFixedEffectModel, ::AbstractDataFrame) = error("residuals is not defined for linear factor models. Use the option save = true")
-function StatsBase.confint(x::InteractiveFixedEffectModel)
-    scale = quantile(TDist(dof_residual(x)), 1 - (1-0.95)/2)
+function StatsBase.confint(x::InteractiveFixedEffectModel; level::Real = 0.95)
+    scale = tdistinvcdf(dof_residual(x), 1 - (1 - level) / 2)
     se = stderror(x)
     hcat(x.coef -  scale * se, x.coef + scale * se)
 end
+function StatsBase.coeftable(x::InteractiveFixedEffectModel; level = 0.95)
+    cc = coef(x)
+    se = stderror(x)
+    coefnms = coefnames(x)
+    conf_int = confint(x; level = level)
+    # put (intercept) last
+    if !isempty(coefnms) && ((coefnms[1] == Symbol("(Intercept)")) || (coefnms[1] == "(Intercept)"))
+        newindex = vcat(2:length(cc), 1)
+        cc = cc[newindex]
+        se = se[newindex]
+        conf_int = conf_int[newindex, :]
+        coefnms = coefnms[newindex]
+    end
+    tt = cc ./ se
+    FixedEffectModels.CoefTable(
+        hcat(cc, se, tt, fdistccdf.(Ref(1), Ref(dof_residual(x)), abs2.(tt)), conf_int[:, 1:2]),
+        ["Estimate","Std.Error","t value", "Pr(>|t|)", "Lower 95%", "Upper 95%" ],
+        ["$(coefnms[i])" for i = 1:length(cc)], 4)
+end
+
+##############################################################################
+##
+## Display Result
+##
+##############################################################################
 
 title(::InteractiveFixedEffectModel) = "Interactive Fixed Effect Model"
 top(x::InteractiveFixedEffectModel) = [
@@ -297,10 +322,9 @@ top(x::InteractiveFixedEffectModel) = [
             "Converged" sprint(show, x.converged; context=:compact => true)
             ]
 
+format_scientific(x) = @sprintf("%.3f", x)
+
 function Base.show(io::IO, x::InteractiveFixedEffectModel)
-    show(io, coeftable(x))
-end
-function StatsBase.coeftable(x::InteractiveFixedEffectModel)
     ctitle = title(x)
     ctop = top(x)
     cc = coef(x)
@@ -316,8 +340,73 @@ function StatsBase.coeftable(x::InteractiveFixedEffectModel)
         coefnms = coefnms[newindex]
     end
     tt = cc ./ se
-    FixedEffectModels.CoefTable2(
-        hcat(cc, se, tt, ccdf.(Ref(FDist(1, dof_residual(x))), abs2.(tt)), conf_int[:, 1:2]),
-        ["Estimate","Std.Error","t value", "Pr(>|t|)", "Lower 95%", "Upper 95%" ],
-        ["$(coefnms[i])" for i = 1:length(cc)], 4, ctitle, ctop)
-end
+    mat = hcat(cc, se, tt, fdistccdf.(Ref(1), Ref(dof_residual(x)), abs2.(tt)), conf_int[:, 1:2])
+    nr, nc = size(mat)
+    colnms = ["Estimate","Std.Error","t value", "Pr(>|t|)", "Lower 95%", "Upper 95%"]
+    rownms = ["$(coefnms[i])" for i = 1:length(cc)]
+    pvc = 4
+
+
+    # print
+    if length(rownms) == 0
+        rownms = AbstractString[lpad("[$i]",floor(Integer, log10(nr))+3) for i in 1:nr]
+    end
+    if length(rownms) > 0
+        rnwidth = max(4,maximum([length(nm) for nm in rownms]) + 1)
+        else
+            # if only intercept, rownms is empty collection, so previous would return error
+        rnwidth = 4
+    end
+    rownms = [rpad(nm,rnwidth) for nm in rownms]
+    widths = [length(cn)::Int for cn in colnms]
+    str = [sprint(show, mat[i,j]; context=:compact => true) for i in 1:nr, j in 1:nc]
+    if pvc != 0                         # format the p-values column
+        for i in 1:nr
+            str[i, pvc] = format_scientific(mat[i, pvc])
+        end
+    end
+    for j in 1:nc
+        for i in 1:nr
+            lij = length(str[i, j])
+            if lij > widths[j]
+                widths[j] = lij
+            end
+        end
+    end
+    widths .+= 1
+    totalwidth = sum(widths) + rnwidth
+    if length(ctitle) > 0
+        halfwidth = div(totalwidth - length(ctitle), 2)
+        println(io, " " ^ halfwidth * string(ctitle) * " " ^ halfwidth)
+    end
+    if length(ctop) > 0
+        for i in 1:size(ctop, 1)
+            ctop[i, 1] = ctop[i, 1] * ":"
+        end
+        println(io, "=" ^totalwidth)
+        halfwidth = div(totalwidth, 2) - 1
+        interwidth = 2 +  mod(totalwidth, 2)
+        for i in 1:(div(size(ctop, 1) - 1, 2)+1)
+            print(io, ctop[2*i-1, 1])
+            print(io, lpad(ctop[2*i-1, 2], halfwidth - length(ctop[2*i-1, 1])))
+            print(io, " " ^interwidth)
+            if size(ctop, 1) >= 2*i
+                print(io, ctop[2*i, 1])
+                print(io, lpad(ctop[2*i, 2], halfwidth - length(ctop[2*i, 1])))
+            end
+            println(io)
+        end
+    end
+    println(io,"=" ^totalwidth)
+    println(io," " ^ rnwidth *
+            join([lpad(string(colnms[i]), widths[i]) for i = 1:nc], ""))
+    println(io,"-" ^totalwidth)
+    for i in 1:nr
+        print(io, rownms[i])
+        for j in 1:nc
+            print(io, lpad(str[i,j],widths[j]))
+        end
+        println(io)
+    end
+    println(io,"=" ^totalwidth)
+end

src/utils/formula.jl

@@ -0,0 +1,17 @@
+##############################################################################
+##
+## Iterate on terms
+##
+##############################################################################
+
+eachterm(@nospecialize(x::AbstractTerm)) = (x,)
+eachterm(@nospecialize(x::NTuple{N, AbstractTerm})) where {N} = x
+TermOrTerms = Union{AbstractTerm, NTuple{N, AbstractTerm} where N}
+
+##############################################################################
+##
+## Parse FixedEffect
+##
+##############################################################################
+fesymbol(t::FixedEffectModels.FixedEffectTerm) = t.x
+fesymbol(t::FunctionTerm{typeof(fe)}) = Symbol(t.args_parsed[1])

src/utils/tss.jl

@@ -0,0 +1,8 @@
+function tss(y::AbstractVector, hasintercept::Bool, weights::AbstractWeights)
+    if hasintercept
+        m = mean(y, weights)
+        sum(@inbounds (y[i] - m)^2 * weights[i] for i in eachindex(y))
+    else
+        sum(@inbounds y[i]^2 * weights[i] for i in eachindex(y))
+    end
+end