Add a type RegressionBasedDIDResult

Author: junyuan-chen

src/InteractionWeightedDIDs.jl

@@ -13,15 +13,16 @@ using SplitApplyCombine: group, mapview
 using StatsBase: AbstractWeights, CovarianceEstimator, UnitWeights, NoQuote, vcov
 using StatsFuns: fdistccdf
 using StatsModels: coefnames
-using Tables: getcolumn
+using Tables: getcolumn, columnnames
 using TypedTables: Table
 using Vcov
 @reexport using DiffinDiffsBase
-using DiffinDiffsBase: termvars, hasintercept, omitsintercept, isintercept, parse_intercept
+using DiffinDiffsBase: termvars, hasintercept, omitsintercept, isintercept, parse_intercept,
+    _getsubcolumns, _treatnames
 
 import Base: show
-import DiffinDiffsBase: required, default, transformed, combinedargs, _getsubcolumns,
-    valid_didargs, result
+import DiffinDiffsBase: required, default, transformed, combinedargs,
+    _get_default, valid_didargs, result
 import FixedEffectModels: has_fe
 
 export Vcov,

src/did.jl

@@ -4,8 +4,129 @@
 Estimation procedure for regression-based difference-in-differences.
 """
 const RegressionBasedDID = DiffinDiffsEstimator{:RegressionBasedDID,
-    Tuple{CheckData, CheckVcov, CheckVars, MakeWeights, CheckFEs, MakeFESolver,
-    MakeYXCols, MakeTreatCols}}
+    Tuple{CheckData, CheckVcov, CheckVars, CheckFEs, MakeWeights, MakeFESolver,
+    MakeYXCols, MakeTreatCols, SolveLeastSquares, EstVcov}}
 
 const Reg = RegressionBasedDID
 
+function _get_default(::Reg, @nospecialize(ntargs::NamedTuple))
+    defaults = (subset=nothing, weightname=nothing, vce=Vcov.RobustCovariance(),
+        treatintterms=(), xterms=(), drop_singletons=true, nfethreads=6,
+        contrasts=Dict{Symbol,Any}(), fetol=1.0e-8, femaxiter=10000, cohortinteracted=true)
+    return merge(defaults, ntargs)
+end
+
+function valid_didargs(d::Type{Reg}, ::DynamicTreatment{SharpDesign},
+        ::TrendParallel{Unconditional, Exact}, @nospecialize(ntargs::NamedTuple))
+    ntargs = _get_default(d(), ntargs)
+    name = haskey(ntargs, :name) ? ntargs.name : ""
+    return name, d, ntargs
+end
+
+struct RegressionBasedDIDResult{CohortInteracted} <: DIDResult
+    coef::Vector{Float64}
+    vcov::Matrix{Float64}
+    vce::CovarianceEstimator
+    tr::AbstractTreatment
+    pr::AbstractParallel
+    cellweights::Vector{Float64}
+    cellcounts::Vector{Int}
+    esample::BitVector
+    nobs::Int
+    dof_residual::Int
+    F::Float64
+    p::Float64
+    yname::String
+    coefnames::Vector{String}
+    coefinds::Dict{String, Int}
+    treatinds::Table
+    yxterms::Dict{AbstractTerm, AbstractTerm}
+    treatname::Symbol
+    contrasts::Dict{Symbol, Any}
+    weightname::Union{Symbol, Nothing}
+    fenames::Vector{Symbol}
+    nfeiterations::Union{Int, Nothing}
+    feconverged::Union{Bool, Nothing}
+    nfesingledropped::Int
+end
+
+function result(::Type{Reg}, @nospecialize(nt::NamedTuple))
+    cellweights = [nt.cellweights[k] for k in nt.treatinds]
+    cellcounts = [nt.cellcounts[k] for k in nt.treatinds]
+    yname = coefnames(nt.yxterms[nt.yterm])
+    tnames = _treatnames(nt.treatinds)
+    cnames = [coefnames(nt.yxterms[x]) for x in nt.xterms if width(nt.yxterms[x])>0]
+    cnames = vcat(tnames, cnames)[nt.basecols]
+    coefinds = Dict(cnames .=> 1:length(cnames))
+    didresult = RegressionBasedDIDResult{nt.cohortinteracted}(
+        nt.coef, nt.vcov_mat, nt.vce, nt.tr, nt.pr, cellweights, cellcounts,
+        nt.esample, sum(nt.esample), nt.dof_resid, nt.F, nt.p,
+        yname, cnames, coefinds, nt.treatinds, nt.yxterms, nt.treatname, nt.contrasts,
+        nt.weightname, nt.fenames, nt.nfeiterations, nt.feconverged, nt.nsingle)
+    return merge(nt, (result=didresult,))
+end
+
+has_fe(r::RegressionBasedDIDResult) = !isempty(r.fenames)
+
+format_scientific(f::Number) = @sprintf("%.3f", f)
+
+function _summary(r::RegressionBasedDIDResult)
+    top = ["Number of obs" sprint(show, nobs(r), context=:compact=>true);
+        "Degrees of freedom" sprint(show, nobs(r) - dof_residual(r), context=:compact=>true);
+        "F-statistic" sprint(show, r.F, context = :compact => true);
+        "p-value" format_scientific(r.p)]
+    fe_info = has_fe(r) ?  
+        ["Converged" sprint(show, r.feconverged, context=:compact=>true);
+        "Singletons dropped" sprint(show, r.nfesingledropped, context=:compact=>true);
+        ] : nothing
+    return top, fe_info
+end
+
+_nunique(t::Table, s::Symbol) = length(unique(getproperty(t, s)))
+
+_excluded_rel_str(tr::DynamicTreatment) = 
+    isempty(tr.exc) ? "None" : join(string.(tr.exc), " ")
+
+_treat_info(tr::DynamicTreatment, trinds::Table, treatname::Symbol) =
+    ["Number of cohorts" sprint(show, _nunique(trinds, treatname), context=:compact=>true);
+    "Interactions within cohorts" sprint(show, length(columnnames(trinds))-2, context=:compact=>true);
+    "Relative time periods" sprint(show, _nunique(trinds, :rel), context=:compact=>true);
+    "Excluded periods" sprint(show, NoQuote(_excluded_rel_str(tr)), context=:compact=>true)]
+
+_treat_spec(r::RegressionBasedDIDResult{true}, tr::DynamicTreatment{SharpDesign}) =
+    "Cohort-interacted sharp dynamic specification"
+
+_treat_spec(r::RegressionBasedDIDResult{false}, tr::DynamicTreatment{SharpDesign}) =
+    "Sharp dynamic specification"
+    
+function show(io::IO, r::RegressionBasedDIDResult;
+    totalwidth::Int=70, interwidth::Int=4+mod(totalwidth,2))
+    halfwidth = div(totalwidth-interwidth, 2)
+    top, fe_info = _summary(r)
+    tr_info = _treat_info(r.tr, r.treatinds, r.treatname)
+    blocks = [top, tr_info, fe_info]
+    fes = has_fe(r) ? join(string.(r.fenames), " ") : "none"
+    fetitle = string("Fixed effects: ", fes)
+    blocktitles = ["Summary of results",
+                   _treat_spec(r, r.tr),
+                   fetitle[1:min(totalwidth,length(fetitle))]]
+    for (ib, b) in enumerate(blocks)
+        for i in 1:size(b, 1)
+            b[i, 1] = b[i, 1] * ":"
+        end
+        println(io, "─" ^totalwidth)
+        println(io, blocktitles[ib])
+        println(io, "─" ^totalwidth)
+        for i in 1:(div(size(b, 1) - 1, 2)+1)
+            print(io, b[2*i-1, 1])
+            print(io, lpad(b[2*i-1, 2], halfwidth - length(b[2*i-1, 1])))
+            print(io, " " ^interwidth)
+            if size(b, 1) >= 2*i
+                print(io, b[2*i, 1])
+                print(io, lpad(b[2*i, 2], halfwidth - length(b[2*i, 1])))
+            end
+            println(io)
+        end
+    end
+    println(io, "─" ^totalwidth)
+end

src/procedures.jl

@@ -67,10 +67,10 @@ Construct `FixedEffects.AbstractFixedEffectSolver`.
 See also [`MakeFESolver`](@ref).
 """
 function makefesolver(fenames::Vector{Symbol}, weights::AbstractWeights, esample::BitVector,
-        fes::Vector{FixedEffect})
+        nfethreads::Int, fes::Vector{FixedEffect})
     if !isempty(fes)
         fes = FixedEffect[fe[esample] for fe in fes]
-        feM = AbstractFixedEffectSolver{Float64}(fes, weights, Val{:cpu}, Threads.nthreads())
+        feM = AbstractFixedEffectSolver{Float64}(fes, weights, Val{:cpu}, nfethreads)
         return (feM=feM,), true
     else
         return (feM=nothing,), true
@@ -86,6 +86,7 @@ The returned object named `feM` may be shared across multiple specifications.
 const MakeFESolver = StatsStep{:MakeFESolver, typeof(makefesolver)}
 
 required(::MakeFESolver) = (:fenames, :weights, :esample)
+default(::MakeFESolver) = (nfethreads=Threads.nthreads(),)
 # Determine equality of fes by fenames
 combinedargs(::MakeFESolver, allntargs) = (allntargs[1].fes,)
 
@@ -214,7 +215,7 @@ function maketreatcols(data, treatname::Symbol, treatintterms::Terms,
         weightname::Union{Symbol, Nothing}, weights::AbstractWeights,
         esample::BitVector, tr_rows::BitVector,
         cohortinteracted::Bool, fetol::Real, femaxiter::Int,
-        ::Type{<:DynamicTreatment{SharpDesign}}, time::Symbol, exc::Set{<:Integer})
+        ::Type{DynamicTreatment{SharpDesign}}, time::Symbol, exc::Set{<:Integer})
 
     nobs = sum(esample)
     tnames = (time, treatname, termvars(treatintterms)...)
@@ -279,7 +280,7 @@ transformed(::MakeTreatCols, @nospecialize(nt::NamedTuple)) =
 combinedargs(step::MakeTreatCols, allntargs) =
     combinedargs(step, allntargs, typeof(allntargs[1].tr))
 
-combinedargs(::MakeTreatCols, allntargs, ::Type{<:DynamicTreatment{SharpDesign}}) =
+combinedargs(::MakeTreatCols, allntargs, ::Type{DynamicTreatment{SharpDesign}}) =
     (Set(intersect((nt.tr.exc for nt in allntargs)...)),)
 
 """
@@ -351,7 +352,8 @@ function _vce(data, esample::BitVector, vce::Vcov.ClusterCovariance,
     concrete_vce = Vcov.materialize(cludata, vce)
     dof_absorb = 0
     for fe in fes
-        any(c->isnested(fe, c.refs), concrete_vce.clusters) && (dof_absorb += 1)
+        # ! To be fixed
+        dof_absorb += any(c->isnested(fe, c.refs), concrete_vce.clusters) ? 1 : nunique(fe)
     end
     return concrete_vce, dof_absorb
 end

test/did.jl

@@ -0,0 +1,17 @@
+@testset "RegressionBasedDID" begin
+    hrs = exampledata("hrs")
+    r = @did(Reg, data=hrs, dynamic(:wave, -1), notyettreated([11]), vce=Vcov.cluster(:hhidpn),
+        yterm=term(:oop_spend), treatname=:wave_hosp, treatintterms=(),
+        xterms=(fe(:wave)+fe(:hhidpn)))
+    @test coef(r, "rel: -3 & wave_hosp: 10") ≈ 591.04639 atol=1e-5
+    @test coef(r, "rel: -2 & wave_hosp: 9") ≈ 298.97735 atol=1e-5
+    @test coef(r, "rel: -2 & wave_hosp: 10") ≈ 410.58102 atol=1e-5
+    @test coef(r, "rel: 0 & wave_hosp: 8") ≈ 2825.5659 atol=1e-4
+    @test coef(r, "rel: 0 & wave_hosp: 9") ≈ 3030.8408 atol=1e-4
+    @test coef(r, "rel: 0 & wave_hosp: 10") ≈ 3091.5084 atol=1e-4
+    @test coef(r, "rel: 1 & wave_hosp: 8") ≈ 825.14585 atol=1e-5
+    @test coef(r, "rel: 1 & wave_hosp: 9") ≈ 106.83785 atol=1e-5
+    @test coef(r, "rel: 2 & wave_hosp: 8") ≈ 800.10647 atol=1e-5
+    @test nobs(r) == 2624
+    
+end

test/procedures.jl

@@ -43,7 +43,8 @@ end
     nobs = size(hrs, 1)
     fes = FixedEffect[FixedEffect(hrs.hhidpn)]
     fenames = [:fe_hhidpn]
-    nt = (fenames=fenames, weights=uweights(nobs), esample=trues(nobs), fes=fes)
+    nt = (fenames=fenames, weights=uweights(nobs), esample=trues(nobs),
+        default(MakeFESolver())..., fes=fes)
     ret, share = makefesolver(nt...)
     @test ret.feM isa FixedEffects.FixedEffectSolverCPU{Float64}
     nt = merge(nt, (fenames=Symbol[], fes=FixedEffect[]))