Support ScaledArray and update dependency (#7)

Author: junyuan-chen

Project.toml

@@ -4,6 +4,7 @@ authors = ["Junyuan Chen <norman.chen.git@gmail.com>"]
 version = "0.2.0"
 
 [deps]
+DataAPI = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
 DiffinDiffsBase = "7fc23ce2-6e83-4e3c-822f-a79085ccb3e8"
 FixedEffectModels = "9d5cd8c9-2029-5cab-9928-427838db53e3"
 FixedEffects = "c8885935-8500-56a7-9867-7708b20db0eb"
@@ -16,9 +17,10 @@ Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 Vcov = "ec2bfdc2-55df-4fc9-b9ae-4958c2cf2486"
 
 [compat]
-DataFrames = "0.22"
-DiffinDiffsBase = "0.3"
-FixedEffectModels = "1.4"
+DataAPI = "1.5"
+DataFrames = "1"
+DiffinDiffsBase = "0.3.2"
+FixedEffectModels = "1.5"
 FixedEffects = "2"
 Reexport = "0.2, 1"
 StatsBase = "0.33"
@@ -30,7 +32,8 @@ julia = "1.3"
 
 [extras]
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["DataFrames", "Test"]
+test = ["DataFrames", "Dates", "Test"]

src/InteractionWeightedDIDs.jl

@@ -1,10 +1,11 @@
 module InteractionWeightedDIDs
 
 using Base: Callable
+using DataAPI: refarray
 using FixedEffectModels: FixedEffectTerm, Combination,
-    fe, _parse_fixedeffect, basecol, isnested, nunique
+    fe, _parse_fixedeffect, invsym!, isnested, nunique
 using FixedEffects
-using LinearAlgebra: Cholesky, Factorization, Symmetric, cholesky!
+using LinearAlgebra: Cholesky, Factorization, Symmetric, cholesky!, diag
 using Reexport
 using StatsBase: AbstractWeights, CovarianceEstimator, UnitWeights, PValue, TestStat, NoQuote
 using StatsFuns: fdistccdf
@@ -13,7 +14,7 @@ using Tables
 using Tables: getcolumn, columnnames
 using Vcov
 @reexport using DiffinDiffsBase
-using DiffinDiffsBase: TimeType, termvars, isintercept, parse_intercept!,
+using DiffinDiffsBase: ValidTimeType, termvars, isintercept, parse_intercept!,
     _count!, _groupfind, _treatnames, _parse_bycells!, _parse_subset
 
 import Base: show

src/did.jl

@@ -78,7 +78,7 @@ function result(::Type{Reg}, @nospecialize(nt::NamedTuple))
     yterm = nt.yxterms[nt.yterm]
     yname = coefnames(yterm)
     cnames = _treatnames(nt.treatcells)
-    cnames = append!(cnames, coefnames.(nt.xterms))[nt.basecols]
+    cnames = append!(cnames, coefnames.(nt.xterms))[nt.basiscols]
     coefinds = Dict(cnames .=> 1:length(cnames))
     didresult = RegressionBasedDIDResult{typeof(nt.tr), nt.cohortinteracted}(
         nt.coef, nt.vcov_mat, nt.vce, nt.tr, nt.pr, nt.cellweights, nt.cellcounts,

src/procedures.jl

@@ -4,11 +4,17 @@
 Exclude rows that are invalid for variance-covariance estimator.
 See also [`CheckVcov`](@ref).
 """
-checkvcov!(data, esample::BitVector,
+checkvcov!(data, esample::BitVector, aux::BitVector,
     vce::Union{Vcov.SimpleCovariance, Vcov.RobustCovariance}) = NamedTuple()
 
-function checkvcov!(data, esample::BitVector, vce::Vcov.ClusterCovariance)
-    esample .&= Vcov.completecases(data, vce)
+function checkvcov!(data, esample::BitVector, aux::BitVector, vce::Vcov.ClusterCovariance)
+    for name in Vcov.names(vce)
+        col = getcolumn(data, name)
+        if Missing <: eltype(col)
+            aux .= .!ismissing.(col)
+            esample .&= aux
+        end
+    end
     return (esample=esample,)
 end
 
@@ -20,7 +26,7 @@ exclude rows that are invalid for variance-covariance estimator.
 """
 const CheckVcov = StatsStep{:CheckVcov, typeof(checkvcov!), true}
 
-required(::CheckVcov) = (:data, :esample)
+required(::CheckVcov) = (:data, :esample, :aux)
 default(::CheckVcov) = (vce=Vcov.robust(),)
 copyargs(::CheckVcov) = (2,)
 
@@ -163,7 +169,7 @@ default(::MakeYXCols) = (contrasts=nothing, fetol=1e-8, femaxiter=10000)
 
 function combinedargs(::MakeYXCols, allntargs)
     yx = TermSet()
-    @inbounds for nt in allntargs
+    for nt in allntargs
         push!(yx, nt.yterm)
         foreach(x->push!(yx, x), nt.xterms)
     end
@@ -182,15 +188,15 @@ function maketreatcols(data, treatname::Symbol, treatintterms::TermSet,
         weights::AbstractWeights, esample::BitVector,
         cohortinteracted::Bool, fetol::Real, femaxiter::Int,
         ::Type{DynamicTreatment{SharpDesign}}, time::Symbol,
-        exc::IdDict{Int,Int}, notreat::IdDict{TimeType,Int})
+        exc::Dict{Int,Int}, notreat::IdDict{ValidTimeType,Int})
 
     nobs = sum(esample)
     # Putting treatname before time avoids sorting twice if cohortinteracted
     cellnames = Symbol[treatname, time, sort!(termvars(treatintterms))...]
     cols = subcolumns(data, cellnames, esample)
     cells, rows = cellrows(cols, findcell(cols))
 
-    rel = cells[2] .- cells[1]
+    rel = refarray(cells[2]) .- refarray(cells[1])
     kept = .!haskey.(Ref(exc), rel) .& .!haskey.(Ref(notreat), cells[1])
     treatrows = rows[kept]
     # Construct cells needed for treatment indicators
@@ -283,17 +289,17 @@ combinedargs(step::MakeTreatCols, allntargs) =
 
 # Obtain the relative time periods excluded by all tr in allntargs
 function combinedargs(::MakeTreatCols, allntargs, ::Type{DynamicTreatment{SharpDesign}})
-    exc = IdDict{Int,Int}()
-    notreat = IdDict{TimeType,Int}()
+    exc = Dict{Int,Int}()
+    notreat = IdDict{ValidTimeType,Int}()
     @inbounds for nt in allntargs
         foreach(x->_count!(exc, x), nt.tr.exc)
         foreach(x->_count!(notreat, x), nt.pr.e)
     end
     nnt = length(allntargs)
-    @inbounds for (k, v) in exc
+    for (k, v) in exc
         v == nnt || delete!(exc, k)
     end
-    @inbounds for (k, v) in notreat
+    for (k, v) in notreat
         v == nnt || delete!(notreat, k)
     end
     return (exc, notreat)
@@ -339,14 +345,14 @@ function solveleastsquares!(tr::DynamicTreatment{SharpDesign}, pr::TrendParallel
     X = hcat(tcols..., (yxcols[x] for x in xs)...)
 
     ntcols = length(tcols)
-    basecols = trues(size(X,2))
+    basiscols = trues(size(X,2))
     if size(X, 2) > ntcols
-        basecols = basecol(X)
+        basiscols = diag(invsym!(X'X)) .> 0
         # Do not drop any treatment indicator
-        sum(basecols[1:ntcols]) == ntcols ||
+        sum(basiscols[1:ntcols]) == ntcols ||
             error("covariates are collinear with treatment indicators")
-        sum(basecols) < size(X, 2) &&
-            (X = X[:, basecols])
+        sum(basiscols) < size(X, 2) &&
+            (X = X[:, basiscols])
     end
 
     crossx = cholesky!(Symmetric(X'X))
@@ -356,7 +362,7 @@ function solveleastsquares!(tr::DynamicTreatment{SharpDesign}, pr::TrendParallel
     return (coef=coef::Vector{Float64}, X=X::Matrix{Float64},
         crossx=crossx::Cholesky{Float64,Matrix{Float64}},
         residuals=residuals::Vector{Float64}, treatcells=treatcells::VecColumnTable,
-        xterms=xs::Vector{AbstractTerm}, basecols=basecols::BitVector)
+        xterms=xs::Vector{AbstractTerm}, basiscols=basiscols::BitVector)
 end
 
 """

test/did.jl

@@ -34,7 +34,7 @@
 
     @test r.ycellweights == r.ycellcounts
     @test r.ycellcounts == repeat([252, 176, 163, 65], inner=4)
-    @test all(i->r.coef[i]≈sum(r.lsweights[:,i].*r.ycellmeans), 1:ntreatcoef(r))
+    @test all(i->r.coef[i] ≈ sum(r.lsweights[:,i].*r.ycellmeans), 1:ntreatcoef(r))
 
     @test sprint(show, r) == "Regression-based DID result"
     pv = VERSION < v"1.6.0" ? " <1e-7" : "<1e-07"
@@ -55,10 +55,21 @@
         Converged:                   true    Singletons dropped:             0
         ──────────────────────────────────────────────────────────────────────"""
 
+    df = DataFrame(hrs)
+    df.wave = Date.(df.wave)
+    df.wave_hosp = Date.(df.wave_hosp)
+    df.wave = settime(df, :wave, step=Year(1))
+    df.wave_hosp = settime(df, :wave_hosp, start=Date(7), step=Year(1))
+    r1 = @did(Reg, data=df, dynamic(:wave, -1), notyettreated(Date(11)),
+        vce=Vcov.cluster(:hhidpn), yterm=term(:oop_spend), treatname=:wave_hosp,
+        treatintterms=(), xterms=(fe(:wave)+fe(:hhidpn)), solvelsweights=false)
+    @test coef(r1) ≈ coef(r)
+    @test r1.coefnames[1] == "wave_hosp: 0008-01-01 & rel: 0"
+
     r = @did(Reg, data=hrs, dynamic(:wave, -1), notyettreated([11]),
         vce=Vcov.cluster(:hhidpn), yterm=term(:oop_spend), treatname=:wave_hosp,
         treatintterms=(), cohortinteracted=false, lswtnames=(:wave_hosp, :wave))
-    @test all(i->r.coef[i]≈sum(r.lsweights[:,i].*r.ycellmeans), 1:ntreatcoef(r))
+    @test all(i->r.coef[i] ≈ sum(r.lsweights[:,i].*r.ycellmeans), 1:ntreatcoef(r))
 
     @test sprint(show, MIME("text/plain"), r) == """
         ──────────────────────────────────────────────────────────────────────
@@ -75,8 +86,8 @@
         ──────────────────────────────────────────────────────────────────────"""
 
     sr = view(r, 1:3)
-    @test coef(sr)[1] ≈ 1560.2174484214308
-    @test vcov(sr)[1] ≈ 1.1248572269146878e6
+    @test coef(sr)[1] == r.coef[1]
+    @test vcov(sr)[1] == r.vcov[1]
     @test parent(sr) === r
 
     tr = rescale(r, fill(2, 5), 1:5)
@@ -137,7 +148,7 @@ end
 
     @test a.cellweights == a.cellcounts
     @test a.cellcounts == [163, 339, 591, 428, 252]
-    @test all(i->a.coef[i]≈sum(a.lsweights[:,i].*r.ycellmeans), 1:ntreatcoef(a))
+    @test all(i->a.coef[i] ≈ sum(a.lsweights[:,i].*r.ycellmeans), 1:ntreatcoef(a))
 
     a1 = agg(r, (:rel,), subset=:rel=>isodd)
     @test length(coef(a1)) == 2