rmv combination (#250)

* rmv combination

* Update README.md

* Update README.md

Author: matthieugomez

Project.toml

@@ -1,6 +1,6 @@
 name = "FixedEffectModels"
 uuid = "9d5cd8c9-2029-5cab-9928-427838db53e3"
-version = "1.9.4"
+version = "1.9.5"
 
 [deps]
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
@@ -19,7 +19,7 @@ Vcov = "ec2bfdc2-55df-4fc9-b9ae-4958c2cf2486"
 
 [compat]
 DataFrames = "0.21, 0.22, 1"
-FixedEffects = "2"
+FixedEffects = "2.3"
 PrecompileTools = "1"
 Reexport = "0.1, 0.2, 1"
 Statistics = "1"

README.md

@@ -32,7 +32,7 @@ reg(df, @formula(Sales ~ NDI + fe(State) + fe(Year)), Vcov.cluster(:State), weig
 ```
 
 
--  A typical formula is composed of one dependent variable, exogeneous variables, endogeneous variables, instrumental variables, and a set of high-dimensional fixed effects.
+-  A typical formula is composed of one dependent variable, exogenous variables, endogenous variables, instrumental variables, and a set of high-dimensional fixed effects.
 
 	```julia
 	dependent variable ~ exogenous variables + (endogenous variables ~ instrumental variables) + fe(fixedeffect variable)
@@ -49,7 +49,7 @@ reg(df, @formula(Sales ~ NDI + fe(State) + fe(Year)), Vcov.cluster(:State), weig
 	reg(df, @formula(Sales ~ (Price ~ Pimin)))
 	```
 
-	To construct formula programatically, use
+	To construct formula programmatically, use
 	```julia
 	reg(df, term(:Sales) ~ term(:NDI) + fe(:State) + fe(:Year))
 	```
@@ -111,7 +111,7 @@ df = dataset("plm", "Cigar")
 reg(df, @formula(Sales ~ NDI + fe(State) + fe(Year)), method = :CUDA)
 ```
 
-The package also supports Apple GPUs with `Metal.jl`, although it does not really improve perfomances
+The package also supports Apple GPUs with `Metal.jl`, although I could not find a way to get better performance
 ```julia
 using Metal, FixedEffectModels
 @assert Metal.functional()

src/fit.jl

@@ -259,21 +259,20 @@ function StatsAPI.fit(::Type{FixedEffectModel},
         # initialize iterations and converged
         iterations = Int[]
         convergeds = Bool[]
+        cols = vcat(eachcol(y), eachcol(Xexo))
         if has_iv
-            Xall = Combination(y, Xexo, Xendo, Z)
-        else
-            Xall = Combination(y, Xexo)
+            append!(cols, eachcol(Xendo), eachcol(Z))
         end
 
         # compute 2-norm (sum of squares) for each variable 
         # (to see if they are collinear with the fixed effects)
-        sumsquares_pre = [sum(abs2, x) for x in eachcol(Xall)]
+        sumsquares_pre = [sum(abs2, x) for x in cols]
 
         # partial out fixed effects
-        _, iterations, convergeds = solve_residuals!(Xall, feM; maxiter = maxiter, tol = tol, progress_bar = progress_bar)
+        _, iterations, convergeds = solve_residuals!(cols, feM; maxiter = maxiter, tol = tol, progress_bar = progress_bar)
 
         # re-compute 2-norm (sum of squares) for each variable
-        sumsquares_post = [sum(abs2, x) for x in eachcol(Xall)]
+        sumsquares_post = [sum(abs2, x) for x in cols]
 
         # mark variables that are likely to be collinear with the fixed effects
         collinear_tol = min(1e-6, tol / 10)

src/partial_out.jl

@@ -112,7 +112,7 @@ function partial_out(
     formula_x_schema = apply_schema(formula_x, schema(formula_x, subdf, contrasts), StatisticalModel)
     X = convert(Matrix{Float64}, modelmatrix(formula_x_schema, subdf))
     if has_fes
-        X, b, c = solve_residuals!(X, feM; maxiter = maxiter, tol = tol, progress_bar = false)
+        _, b, c = solve_residuals!(eachcol(X), feM; maxiter = maxiter, tol = tol, progress_bar = false)
         append!(iterations, b)
         append!(convergeds, c)
     end
@@ -122,7 +122,7 @@ function partial_out(
     end
     # Compute residuals
     if size(X, 2) > 0
-        mul!(Y, X, X\Y, -1.0, 1.0)
+        mul!(Y, X, X \ Y, -1.0, 1.0)
     end
     residuals = Y
 

src/utils/basecol.jl

@@ -1,30 +1,3 @@
-##############################################################################
-##
-## Combination behaves like [A B C ...] without forming it
-## 
-##############################################################################
-
-struct Combination{T} <: AbstractMatrix{T}
-    A::Tuple
-    cumlength::Vector{Int}
-end
-
-function Combination(A::Union{AbstractVector{T}, AbstractMatrix{T}}...) where {T}
-    Combination{T}(A, cumsum([size(x, 2) for x in A]))
-end
-
-Combination() = error("`Combination` requires at least one argument")
-
-Base.size(c::Combination) = (size(c.A[1], 1), c.cumlength[end])
-Base.size(c::Combination, i::Integer) = size(c)[i]
-
-function Base.view(c::Combination, ::Colon, j)
-    index = searchsortedfirst(c.cumlength, j)
-    newj = index == 1 ? j : j - c.cumlength[index-1]
-    view(c.A[index], :, newj)
-end
-
-
 ##############################################################################
 ##
 ## Returns base of [A B C ...]
@@ -54,7 +27,11 @@ end
 
 # generalized 2inverse
 function invsym!(X::AbstractMatrix; has_intercept = false)
-    # The C value adjusts the check to the relative scale of the variable. The C value is equal to the corrected sum of squares for the variable, unless the corrected sum of squares is 0, in which case C is 1. If you specify the NOINT option but not the ABSORB statement, PROC GLM uses the uncorrected sum of squares instead. The default value of the SINGULAR= option, 107, might be too small, but this value is necessary in order to handle the high-degree polynomials used in the literature to compare regression routin
+    # Options from SAS
+    # The C value adjusts the check to the relative scale of the variable. 
+    # The C value is equal to the corrected sum of squares for the variable, unless the corrected sum of squares is 0, in which case C is 1. 
+    # If you specify the NOINT option but not the ABSORB statement, PROC GLM uses the uncorrected sum of squares instead. 
+    # The default value of the SINGULAR= option, 107, might be too small, but this value is necessary in order to handle the high-degree polynomials used in the literature to compare regression routines
     tols = max.(diag(X), 1)
     for j in 1:size(X, 1)
         d = X[j,j]
@@ -82,7 +59,6 @@ function getcols(X::AbstractMatrix,  basecolX::AbstractVector)
     sum(basecolX) == size(X, 2) ? X : X[:, basecolX]
 end
 
-
 ##############################################################################
 # Auxiliary functions to find columns of exogeneous, endogenous and IV variables
 ##############################################################################

test/fit.jl

@@ -661,13 +661,39 @@ end
 	@test coef(x) ≈ [- 0.11752306001586807] atol = 1e-4
 	@test x.iterations <= 50
 
+
+
+	# add tests with missing fixed effects
+	df.Firm_missing = ifelse.(df.Firm .<= 30, missing, df.Firm)
+
+	## test with missing fixed effects
+	m = @formula Wage ~ Emp + fe(Firm_missing)
+	x = reg(df, m)
+	@test coef(x) ≈ [-.1093657] atol = 1e-4
+	@test stderror(x) ≈  [.032949 ] atol = 1e-4
+	@test r2(x) ≈ 0.8703 atol = 1e-2
+	@test adjr2(x) ≈ 0.8502 atol = 1e-2
+	@test x.nobs == 821
+
+	## test with missing interaction
+	df.Year2 = df.Year .>= 1980
+	m = @formula Wage ~ Emp + fe(Firm_missing) & fe(Year2)
+	x = reg(df, m)
+	@test coef(x) ≈ [-0.100863] atol = 1e-4
+	@test stderror(x) ≈  [0.04149] atol = 1e-4
+	@test x.nobs == 821
+end
+
+
+@testset "gpu" begin
 	methods_vec = [:cpu]
 	if CUDA.functional()
 		push!(methods_vec, :CUDA)
 	end
-	if Metal.functional()
-		push!(methods_vec, :Metal)
-	end
+	#if Metal.functional()
+	#	push!(methods_vec, :Metal)
+	#end
+	df = DataFrame(CSV.File(joinpath(dirname(pathof(FixedEffectModels)), "../dataset/EmplUK.csv")))
 	for method in methods_vec
 		# same thing with float32 precision
 		local m = @formula Wage ~ Emp + fe(Firm)
@@ -692,28 +718,6 @@ end
 		local x = reg(df, m, weights = :Output, method = method, double_precision = false)
 		@test coef(x) ≈  [- 0.043475472188120416] atol = 1e-3
 	end
-
-
-	# add tests with missing fixed effects
-	df.Firm_missing = ifelse.(df.Firm .<= 30, missing, df.Firm)
-
-
-	## test with missing fixed effects
-	m = @formula Wage ~ Emp + fe(Firm_missing)
-	x = reg(df, m)
-	@test coef(x) ≈ [-.1093657] atol = 1e-4
-	@test stderror(x) ≈  [.032949 ] atol = 1e-4
-	@test r2(x) ≈ 0.8703 atol = 1e-2
-	@test adjr2(x) ≈ 0.8502 atol = 1e-2
-	@test x.nobs == 821
-
-	## test with missing interaction
-	df.Year2 = df.Year .>= 1980
-	m = @formula Wage ~ Emp + fe(Firm_missing) & fe(Year2)
-	x = reg(df, m)
-	@test coef(x) ≈ [-0.100863] atol = 1e-4
-	@test stderror(x) ≈  [0.04149] atol = 1e-4
-	@test x.nobs == 821
 end
 
 

test/predict.jl

@@ -383,14 +383,18 @@ end
 	model = @formula Sales ~ (Price ~ Pimin) + fe(State) + fe(Year)
 	result = reg(df, model, subset = df.State .<= 30, weights = :Pop, save = true)
 	@test fe(result)[1, :fe_Year] + fe(result)[1, :fe_State] ≈ 158.91798 atol = 1e-4
+end
+
 
+@testset "gpu" begin
+	df = DataFrame(CSV.File(joinpath(dirname(pathof(FixedEffectModels)), "../dataset/Cigar.csv")))
 	methods_vec = [:cpu]
 	if CUDA.functional()
 		push!(methods_vec, :CUDA)
 	end
-	if Metal.functional()
-		push!(methods_vec, :Metal)
-	end
+	#if Metal.functional()
+	#	push!(methods_vec, :Metal)
+	#end
 	for method in methods_vec
 		local model = @formula Sales ~ Price + fe(Year)
 		local result = reg(df, model, save = true, method = method, double_precision = false)