Wrong results in a large data set with one set of FE

[droodman]

I'm running reg() on data set with 9 million rows and nearly 100 non-absorbed regressors, most of which are dummies for certain birth years. I am absorbing only 1 set of FE. Occasionally some of the estimates are clearly wrong. The coefficients on the birth year dummies should vary rather smoothly with birth year, but they sometimes make big jumps. Running the same regression with reghdfe or areg in Stata does not have this problem.

I will paste the output of an example. I cannot share the data publicly but will email it to @matthieugomez

julia> using CSV, DataFrames, Plots, FixedEffectModels

julia> df = CSV.read("c:/users/drood/Downloads/FEbug.csv", DataFrame)
8830997×96 DataFrame
     Row │ part     _Ibirthyr_1907  _Ibirthyr_1908  _Ibirthyr_1909  _Ibirthyr_1910  _Ibirthyr_1911  _Ibirthyr_1912  _I ⋯
         │ Float64  Float64         Float64         Float64         Float64         Float64         Float64         Fl ⋯
─────────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────────
       1 │     0.0             0.0             0.0             0.0             0.0             0.0             0.0     ⋯
       2 │     0.0             0.0             0.0             0.0             0.0             0.0             0.0
       3 │     0.0             0.0             0.0             0.0             0.0             0.0             0.0
       4 │     0.0             0.0             0.0             0.0             0.0             0.0             0.0
       5 │     0.0             0.0             0.0             0.0             0.0             0.0             0.0     ⋯
       6 │     0.0             0.0             0.0             0.0             0.0             0.0             0.0
       7 │     0.0             0.0             0.0             0.0             0.0             0.0             0.0
       8 │     0.0             0.0             0.0             0.0             0.0             0.0             0.0
       9 │     0.0             0.0             0.0             0.0             0.0             0.0             0.0     ⋯
      10 │     0.0             0.0             0.0             0.0             0.0             0.0             0.0
      11 │     0.0             0.0             0.0             0.0             0.0             0.0             0.0
    ⋮    │    ⋮           ⋮               ⋮               ⋮               ⋮               ⋮               ⋮            ⋱
 8830988 │     0.0             0.0             0.0             0.0             0.0             0.0             0.0
 8830989 │     0.0             0.0             0.0             0.0             0.0             0.0             0.0     ⋯
 8830990 │     0.0             0.0             0.0             0.0             0.0             0.0             0.0
 8830991 │     1.0             0.0             0.0             0.0             0.0             0.0             0.0
 8830992 │     0.0             0.0             0.0             0.0             0.0             0.0             0.0
 8830993 │     0.0             0.0             0.0             0.0             0.0             0.0             0.0     ⋯
 8830994 │     1.0             0.0             0.0             0.0             0.0             0.0             0.0
 8830995 │     1.0             0.0             0.0             0.0             0.0             0.0             0.0
 8830996 │     0.0             0.0             0.0             0.0             0.0             0.0             0.0
 8830997 │     0.0             0.0             0.0             0.0             0.0             0.0             0.0     ⋯
                                                                                     89 columns and 8830976 rows omitted
julia> m = reg(df, @formula(part ~ _Ibirthyr_1907+_Ibirthyr_1908+_Ibirthyr_1909+_Ibirthyr_1910+_Ibirthyr_1911+_Ibirthyr_1912+_Ibirthyr_1913+_Ibirthyr_1914+
_Ibirthyr_1915+_Ibirthyr_1916+_Ibirthyr_1917+_Ibirthyr_1918+_Ibirthyr_1919+_Ibirthyr_1920+_Ibirthyr_1921+_Ibirthyr_1922+_Ibirthyr_1923+_Ibirthyr_1924+
_Ibirthyr_1925+_Ibirthyr_1926+_Ibirthyr_1927+_Ibirthyr_1928+_Ibirthyr_1929+_Ibirthyr_1930+_Ibirthyr_1931+_Ibirthyr_1932+_Ibirthyr_1933+_Ibirthyr_1934+
_Ibirthyr_1935+_Ibirthyr_1936+_Ibirthyr_1937+_Ibirthyr_1938+_Ibirthyr_1939+_Ibirthyr_1940+_Ibirthyr_1941+_Ibirthyr_1942+_Ibirthyr_1943+_Ibirthyr_1944+
_Ibirthyr_1945+_Ibirthyr_1946+_Ibirthyr_1947+_Ibirthyr_1948+_Ibirthyr_1949+_Ibirthyr_1950+_Ibirthyr_1951+_Ibirthyr_1952+_Ibirthyr_1953+_Ibirthyr_1954+
_Ibirthyr_1955+_Ibirthyr_1956+_Ibirthyr_1957+_Ibirthyr_1958+_Ibirthyr_1959+_Ibirthyr_1960+__000001+_Ibirthyr_1962+_Ibirthyr_1963+_Ibirthyr_1964+
_Ibirthyr_1965+_Ibirthyr_1966+_Ibirthyr_1967+_Ibirthyr_1968+_Ibirthyr_1969+_Ibirthyr_1970+_Ibirthyr_1971+_Ibirthyr_1972+_Ibirthyr_1973+_Ibirthyr_1974+
_Ibirthyr_1975+_Ibirthyr_1976+_Ibirthyr_1977+_Ibirthyr_1978+_Ibirthyr_1979+_Ibirthyr_1980+_Ibirthyr_1981+_Ibirthyr_1982+_Ibirthyr_1983+_Ibirthyr_1984+
_Ibirthyr_1985+_Ibirthyr_1986+_Ibirthyr_1987+_Ibirthyr_1988+_Ibirthyr_1989+_Ibirthyr_1990+_Ibirthyr_1991+_Ibirthyr_1992+_Ibirthyr_1993+_Ibirthyr_1994+
_Ibirthyr_1995+_Ibirthyr_1996+_Ibirthyr_1997+_Ibirthyr_1998+_Ibirthyr_1999+_Ibirthyr_2000  + fe(survey))    , nthreads = 4 , tol=1.00000000000e-08)

                               Fixed Effect Model
================================================================================
Number of obs:                  8830997  Degrees of freedom:                  92
R2:                               0.024  R2 Adjusted:                      0.024
F-Stat:                         1675.69  p-value:                          0.000
R2 within:                        0.015  Iterations:                           1
================================================================================
part           |    Estimate  Std.Error  t value Pr(>|t|)   Lower 95%  Upper 95%
--------------------------------------------------------------------------------
_Ibirthyr_1907 |   -0.245888  0.0399852 -6.14948    0.000   -0.324257  -0.167518
_Ibirthyr_1908 |   -0.240937  0.0200867 -11.9948    0.000   -0.280306  -0.201568
_Ibirthyr_1909 |   -0.263179  0.0197842 -13.3025    0.000   -0.301955  -0.224402
_Ibirthyr_1910 |   -0.254113  0.0129464 -19.6281    0.000   -0.279487  -0.228738
_Ibirthyr_1911 |   -0.255636  0.0138299 -18.4844    0.000   -0.282742   -0.22853
_Ibirthyr_1912 |   -0.256784  0.0123504 -20.7916    0.000    -0.28099  -0.232578
_Ibirthyr_1913 |   -0.252984  0.0109203 -23.1664    0.000   -0.274388  -0.231581
_Ibirthyr_1914 |   -0.256519 0.00983208   -26.09    0.000    -0.27579  -0.237249
_Ibirthyr_1915 |   -0.258301 0.00771205 -33.4931    0.000   -0.273416  -0.243185
_Ibirthyr_1916 |   -0.251301   0.007272 -34.5574    0.000   -0.265554  -0.237048
_Ibirthyr_1917 |   -0.251159 0.00719553 -34.9049    0.000   -0.265262  -0.237056
_Ibirthyr_1918 |   -0.252904 0.00571224  -44.274    0.000     -0.2641  -0.241708
_Ibirthyr_1919 |   -0.250968 0.00549119 -45.7037    0.000    -0.26173  -0.240205
_Ibirthyr_1920 |   -0.248248 0.00458863 -54.1007    0.000   -0.257241  -0.239254
_Ibirthyr_1921 |   -0.245237 0.00439951 -55.7418    0.000    -0.25386  -0.236614
_Ibirthyr_1922 |   -0.243647 0.00449236 -54.2358    0.000   -0.252452  -0.234842
_Ibirthyr_1923 |   -0.243476 0.00395942 -61.4929    0.000   -0.251237  -0.235716
_Ibirthyr_1924 |   -0.239739 0.00386321  -62.057    0.000   -0.247311  -0.232168
_Ibirthyr_1925 |   -0.240398 0.00329999 -72.8482    0.000   -0.246866   -0.23393
_Ibirthyr_1926 |    -0.23616 0.00329518 -71.6683    0.000   -0.242619  -0.229702
_Ibirthyr_1927 |   -0.234743 0.00311323 -75.4019    0.000   -0.240845  -0.228642
_Ibirthyr_1928 |   -0.229105 0.00275798 -83.0698    0.000   -0.234511    -0.2237
_Ibirthyr_1929 |   -0.228704 0.00270057 -84.6874    0.000   -0.233997  -0.223411
_Ibirthyr_1930 |   -0.225933 0.00233095 -96.9273    0.000   -0.230501  -0.221364
_Ibirthyr_1931 |   -0.228161 0.00233622 -97.6622    0.000    -0.23274  -0.223582
_Ibirthyr_1932 |   -0.228699 0.00226792 -100.841    0.000   -0.233144  -0.224254
_Ibirthyr_1933 |   -0.228606 0.00208562  -109.61    0.000   -0.232694  -0.224518
_Ibirthyr_1934 |   -0.228292 0.00204416  -111.68    0.000   -0.232299  -0.224286
_Ibirthyr_1935 |   -0.226277 0.00187104 -120.936    0.000   -0.229945   -0.22261
_Ibirthyr_1936 |   -0.222548  0.0018566 -119.868    0.000   -0.226187  -0.218909
_Ibirthyr_1937 |   -0.222628 0.00184952  -120.37    0.000   -0.226253  -0.219003
_Ibirthyr_1938 |   -0.220128 0.00173838 -126.628    0.000   -0.223535  -0.216721
_Ibirthyr_1939 |   -0.219995 0.00170748 -128.842    0.000   -0.223342  -0.216648
_Ibirthyr_1940 |   -0.220104 0.00156163 -140.945    0.000   -0.223165  -0.217043
_Ibirthyr_1941 |     -0.2166 0.00160955 -134.572    0.000   -0.219755  -0.213446
_Ibirthyr_1942 |   -0.215136 0.00156292  -137.65    0.000   -0.218199  -0.212072
_Ibirthyr_1943 |   -0.212369   0.001527 -139.076    0.000   -0.215362  -0.209376
_Ibirthyr_1944 |   -0.210578 0.00156345 -134.688    0.000   -0.213642  -0.207513
_Ibirthyr_1945 |   -0.208761 0.00143075  -145.91    0.000   -0.211565  -0.205957
_Ibirthyr_1946 |   -0.135636 0.00145598 -93.1579    0.000    -0.13849  -0.132783
_Ibirthyr_1947 |   -0.134748 0.00147502 -91.3535    0.000   -0.137639  -0.131857
_Ibirthyr_1948 |   -0.132238 0.00141273 -93.6044    0.000   -0.135007  -0.129469
_Ibirthyr_1949 |   -0.132401 0.00140616 -94.1575    0.000   -0.135157  -0.129645
_Ibirthyr_1950 |   -0.127504 0.00130647 -97.5938    0.000   -0.130064  -0.124943
_Ibirthyr_1951 |   -0.125716 0.00133446 -94.2076    0.000   -0.128332  -0.123101
_Ibirthyr_1952 |    -0.12547 0.00132191 -94.9159    0.000   -0.128061   -0.12288
_Ibirthyr_1953 |    -0.12385 0.00127673 -97.0059    0.000   -0.126353  -0.121348
_Ibirthyr_1954 |   -0.122553 0.00126125 -97.1672    0.000   -0.125025  -0.120081
_Ibirthyr_1955 |   -0.173384 0.00121755 -142.404    0.000   -0.175771  -0.170998
_Ibirthyr_1956 |   -0.167128 0.00121466 -137.592    0.000   -0.169509  -0.164748
_Ibirthyr_1957 |   -0.162486 0.00121038 -134.243    0.000   -0.164858  -0.160113
_Ibirthyr_1958 |   -0.156857 0.00117068 -133.988    0.000   -0.159151  -0.154562
_Ibirthyr_1959 |   -0.146719 0.00116607 -125.824    0.000   -0.149004  -0.144433
_Ibirthyr_1960 |   -0.144422  0.0011126 -129.806    0.000   -0.146603  -0.142242
__000001       |         0.0        NaN      NaN      NaN         NaN        NaN
_Ibirthyr_1962 |   -0.122424 0.00113492 -107.871    0.000   -0.124649    -0.1202
_Ibirthyr_1963 |   -0.118468 0.00110999 -106.729    0.000   -0.120643  -0.116292
_Ibirthyr_1964 |   -0.111601 0.00111507 -100.085    0.000   -0.113787  -0.109416
_Ibirthyr_1965 |   -0.110394 0.00108402 -101.837    0.000   -0.112518  -0.108269
_Ibirthyr_1966 |   -0.108949 0.00110209 -98.8573    0.000   -0.111109  -0.106789
_Ibirthyr_1967 |   -0.108906 0.00111755 -97.4501    0.000   -0.111096  -0.106715
_Ibirthyr_1968 |   -0.108402  0.0010982 -98.7087    0.000   -0.110554  -0.106249
_Ibirthyr_1969 |   -0.108043 0.00108709 -99.3873    0.000   -0.110174  -0.105913
_Ibirthyr_1970 |   -0.109195 0.00106885 -102.161    0.000    -0.11129    -0.1071
_Ibirthyr_1971 |   -0.109449 0.00108716 -100.675    0.000    -0.11158  -0.107319
_Ibirthyr_1972 |   -0.107718 0.00108037 -99.7054    0.000   -0.109836  -0.105601
_Ibirthyr_1973 |   -0.110267 0.00108196 -101.915    0.000   -0.112388  -0.108146
_Ibirthyr_1974 |   -0.102156 0.00109414 -93.3658    0.000     -0.1043  -0.100011
_Ibirthyr_1975 |   -0.101302 0.00108169  -93.652    0.000   -0.103422  -0.099182
_Ibirthyr_1976 |  -0.0952319 0.00109914  -86.642    0.000  -0.0973862 -0.0930777
_Ibirthyr_1977 |  -0.0884749 0.00111556   -79.31    0.000  -0.0906614 -0.0862885
_Ibirthyr_1978 |  -0.0823051 0.00112864 -72.9244    0.000  -0.0845171  -0.080093
_Ibirthyr_1979 |   -0.078953 0.00113875 -69.3331    0.000  -0.0811849 -0.0767211
_Ibirthyr_1980 |  -0.0794783 0.00111344 -71.3805    0.000  -0.0816606  -0.077296
_Ibirthyr_1981 |  -0.0677975 0.00114323 -59.3032    0.000  -0.0700382 -0.0655568
_Ibirthyr_1982 |  -0.0657529 0.00113086  -58.144    0.000  -0.0679694 -0.0635365
_Ibirthyr_1983 |   -0.106154  0.0011442 -92.7764    0.000   -0.108397  -0.103912
_Ibirthyr_1984 |   -0.107567 0.00116716 -92.1618    0.000   -0.109855   -0.10528
_Ibirthyr_1985 |   -0.105736 0.00118811 -88.9951    0.000   -0.108065  -0.103407
_Ibirthyr_1986 |   -0.105918 0.00122616 -86.3819    0.000   -0.108322  -0.103515
_Ibirthyr_1987 |   -0.104192 0.00127568 -81.6756    0.000   -0.106693  -0.101692
_Ibirthyr_1988 |   -0.104842 0.00131977 -79.4397    0.000   -0.107429  -0.102255
_Ibirthyr_1989 |   -0.102505 0.00139315 -73.5781    0.000   -0.105236 -0.0997745
_Ibirthyr_1990 |   -0.103979 0.00146044 -71.1967    0.000   -0.106841  -0.101116
_Ibirthyr_1991 |  -0.0979455 0.00159499 -61.4082    0.000   -0.101072 -0.0948194
_Ibirthyr_1992 | 0.000426894 0.00172852 0.246971    0.805 -0.00296094 0.00381473
_Ibirthyr_1993 |  0.00524459 0.00182628  2.87174    0.004  0.00166515 0.00882403
_Ibirthyr_1994 |   0.0174228 0.00187872  9.27377    0.000   0.0137406   0.021105
_Ibirthyr_1995 |   0.0275258 0.00194972  14.1178    0.000   0.0237044  0.0313472
_Ibirthyr_1996 |   0.0341425 0.00205909  16.5813    0.000   0.0301067  0.0381782
_Ibirthyr_1997 |   0.0344863 0.00229145    15.05    0.000   0.0299952  0.0389775
_Ibirthyr_1998 |   0.0286934 0.00274934  10.4364    0.000   0.0233048   0.034082
_Ibirthyr_1999 |     0.01829 0.00367653  4.97479    0.000   0.0110841  0.0254959
_Ibirthyr_2000 |   0.0258094  0.0085951  3.00281    0.003  0.00896333  0.0426555
================================================================================

julia> plot(coef(m))

[matthieugomez]

Thanks for the report. Does this happen without the survey fixed effects too?

[droodman]

Random thought: with WildBootTests.jl, I have struggled with the lack of a function like Stata's invsym()--a matrix inversion function that is precise and gracefully handles rank-deficient, symmetric, positive-semi-definite matrices like X'X when X doesn't have full column rank. It leads to strange numerical behavior in sporadic instances. I've found that pinv() is often imprecise, I think even when I follow the recommendation to set rtol = sqrt(eps(real(float(oneunit(eltype(M)))))).

Perhaps something like that is happening here?

[nilshg]

Don't have anything to contribute on the numerical issues here but just thought I'd point out that when you've got a formula that long it's often useful to constrict it programmatically, here I think your formula is equivalent to:

term("part") ~ sum([term("_IBirthyr_$i") for i ∈ 1907:2000]) + fe(:survey)

[droodman]

Thanks @nilshg. In fact that long command is programmatically generated. I've written a Stata program called reghdfejl that is meant to be a slot-in replacement for reghdfe but calls FixedEffectModels.jl. A distinct component is a Stata-Julia bridge written in C, which so far I've only compiled in Windows. Once I get it compiled in Linux MacOs too, and documented everything, I'll share it all.

[nilshg]

I should have guessed that you're an advanced enough user to know this, but I thought I mention it just in case. Have learned a lot over the years from your blogs, good to see you involved here!

[matthieugomez]

I'm actually using some version of invsym to get the basis of a set of vectors: https://github.com/FixedEffects/FixedEffectModels.jl/blob/master/src/utils/basecol.jl#L56

[droodman]

Oh interesting. I was looking at how GLM does it:. It calls cholesky() with check = false, zeros out certain entries, then divides X'X into X'Y using LAPACK.potrs!(). https://github.com/JuliaStats/GLM.jl/blob/b1ba4c5fdd5030b98a8cf9fe9c46319e5f5eb20e/src/linpred.jl#L186 But I presume your solution works well too.

[droodman]

@matthieugomez, yes, it happens even with lm()! Whereas it does not happen in R and Stata. I found a GLM issue with a long discussion of instability.

One idea in that discussion is that it is more stable to use the QR decomposition of X than the Cholesky decomposition of X'X. This says that's what R does. That may be the root of the matter.

[droodman]

Ah. GLM is gaining a method=:qr option. FixedEffectModels could copy that option, possibly switching to the QR as the default, as seems to be done in R and Stata(?); or use a hybrid approach in which it uses the Cholesky decomposition is used when the condition number is below some threshold.

[matthieugomez]

I cannot reproduce this issue on my computer (Mac). Can you try

using CSV, DataFrames, LinearAlgebra, FixedEffectModels, Plots
df = CSV.read(".../FEbug.csv", DataFrame)
f = @formula(part ~ _Ibirthyr_1907+_Ibirthyr_1908+_Ibirthyr_1909+_Ibirthyr_1910+_Ibirthyr_1911+_Ibirthyr_1912+_Ibirthyr_1913+_Ibirthyr_1914+
     _Ibirthyr_1915+_Ibirthyr_1916+_Ibirthyr_1917+_Ibirthyr_1918+_Ibirthyr_1919+_Ibirthyr_1920+_Ibirthyr_1921+_Ibirthyr_1922+_Ibirthyr_1923+_Ibirthyr_1924+
     _Ibirthyr_1925+_Ibirthyr_1926+_Ibirthyr_1927+_Ibirthyr_1928+_Ibirthyr_1929+_Ibirthyr_1930+_Ibirthyr_1931+_Ibirthyr_1932+_Ibirthyr_1933+_Ibirthyr_1934+
     _Ibirthyr_1935+_Ibirthyr_1936+_Ibirthyr_1937+_Ibirthyr_1938+_Ibirthyr_1939+_Ibirthyr_1940+_Ibirthyr_1941+_Ibirthyr_1942+_Ibirthyr_1943+_Ibirthyr_1944+
     _Ibirthyr_1945+_Ibirthyr_1946+_Ibirthyr_1947+_Ibirthyr_1948+_Ibirthyr_1949+_Ibirthyr_1950+_Ibirthyr_1951+_Ibirthyr_1952+_Ibirthyr_1953+_Ibirthyr_1954+
     _Ibirthyr_1955+_Ibirthyr_1956+_Ibirthyr_1957+_Ibirthyr_1958+_Ibirthyr_1959+_Ibirthyr_1960+_Ibirthyr_1962+_Ibirthyr_1963+_Ibirthyr_1964+
     _Ibirthyr_1965+_Ibirthyr_1966+_Ibirthyr_1967+_Ibirthyr_1968+_Ibirthyr_1969+_Ibirthyr_1970+_Ibirthyr_1971+_Ibirthyr_1972+_Ibirthyr_1973+_Ibirthyr_1974+
     _Ibirthyr_1975+_Ibirthyr_1976+_Ibirthyr_1977+_Ibirthyr_1978+_Ibirthyr_1979+_Ibirthyr_1980+_Ibirthyr_1981+_Ibirthyr_1982+_Ibirthyr_1983+_Ibirthyr_1984+
     _Ibirthyr_1985+_Ibirthyr_1986+_Ibirthyr_1987+_Ibirthyr_1988+_Ibirthyr_1989+_Ibirthyr_1990+_Ibirthyr_1991+_Ibirthyr_1992+_Ibirthyr_1993+_Ibirthyr_1994+
     _Ibirthyr_1995+_Ibirthyr_1996+_Ibirthyr_1997+_Ibirthyr_1998+_Ibirthyr_1999+_Ibirthyr_2000)

X = modelmatrix(f, df) 
y = response(f, df) 

coef_qr = qr(X) \ y
coef_chol = cholesky!(X' * X) \ (X' * y)

plot([coef_qr coef_chol])

And tell me whether coef_qr and coef_chol give correct results?

[droodman]

Yes, they look identical, and they look good.

[matthieugomez]

Thanks. Can you try reg(df, f) then? If it reports correct results, can you try adding fe(survey) and _00001 to the formula, separately, to see which one is creating the incorrect result?

[droodman]

Also, I have tried the original example on a Mac with an M2 and an Intel Mac, and it worked fine on both.

[droodman]

Alas, reg(df,f) is failing:

julia> m=reg(df,f)
                                  FixedEffectModel
=====================================================================================
Number of obs:                    8830997   Converged:                           true
dof (model):                           93   dof (residuals):                  8830902
R²:                                 0.020   R² adjusted:                        0.020
F-statistic:                      2467.09   P-value:                            0.000
=====================================================================================
                  Estimate   Std. Error      t-stat  Pr(>|t|)   Lower 95%   Upper 95%
─────────────────────────────────────────────────────────────────────────────────────
_Ibirthyr_1907  -0.273123   0.0400529      -6.81906    <1e-11  -0.351625   -0.194621
_Ibirthyr_1908  -0.266322   0.0201056     -13.2462     <1e-39  -0.305729   -0.226916
_Ibirthyr_1909  -0.283199   0.0198092     -14.2963     <1e-45  -0.322025   -0.244374
_Ibirthyr_1910  -0.273939   0.0129585     -21.1398     <1e-98  -0.299337   -0.248541
_Ibirthyr_1911  -0.27852    0.0138402     -20.1239     <1e-89  -0.305646   -0.251394
_Ibirthyr_1912  -0.27633    0.0123648     -22.3481     <1e-99  -0.300565   -0.252095
_Ibirthyr_1913  -0.272623   0.0109313     -24.9397     <1e-99  -0.294047   -0.251198
_Ibirthyr_1914  -0.275257   0.00984198    -27.9676     <1e-99  -0.294546   -0.255967
_Ibirthyr_1915  -0.278617   0.00771855    -36.0971     <1e-99  -0.293745   -0.263489
_Ibirthyr_1916  -0.272556   0.00727289    -37.4756     <1e-99  -0.28681    -0.258301
_Ibirthyr_1917  -0.269847   0.00720218    -37.4674     <1e-99  -0.283963   -0.255731
_Ibirthyr_1918  -0.272885   0.00571428    -47.755      <1e-99  -0.284085   -0.261686
_Ibirthyr_1919  -0.269995   0.00549542    -49.1309     <1e-99  -0.280765   -0.259224
_Ibirthyr_1920  -0.267947   0.00459044    -58.3707     <1e-99  -0.276945   -0.25895
_Ibirthyr_1921  -0.26684    0.00440162    -60.623      <1e-99  -0.275467   -0.258212
_Ibirthyr_1922  -0.264173   0.00449576    -58.7605     <1e-99  -0.272984   -0.255361
_Ibirthyr_1923  -0.265196   0.00396114    -66.9493     <1e-99  -0.27296    -0.257432
_Ibirthyr_1924  -0.261272   0.0038648     -67.6032     <1e-99  -0.268847   -0.253698
_Ibirthyr_1925  -0.263561   0.0032994     -79.8814     <1e-99  -0.270028   -0.257094
_Ibirthyr_1926  -0.260826   0.00329492    -79.1601     <1e-99  -0.267284   -0.254368
_Ibirthyr_1927  -0.262224   0.00311456    -84.193      <1e-99  -0.268329   -0.25612
_Ibirthyr_1928  -0.259789   0.00275803    -94.1937     <1e-99  -0.265195   -0.254383
_Ibirthyr_1929  -0.259259   0.00270146    -95.9699     <1e-99  -0.264554   -0.253964
_Ibirthyr_1930  -0.257863   0.00233074   -110.636      <1e-99  -0.262432   -0.253295
_Ibirthyr_1931  -0.257178   0.00233809   -109.995      <1e-99  -0.26176    -0.252595
_Ibirthyr_1932  -0.257484   0.00226968   -113.445      <1e-99  -0.261933   -0.253036
_Ibirthyr_1933  -0.256838   0.00208717   -123.056      <1e-99  -0.260929   -0.252747
_Ibirthyr_1934  -0.25508    0.00204603   -124.671      <1e-99  -0.25909    -0.25107
_Ibirthyr_1935  -0.253795   0.00187253   -135.536      <1e-99  -0.257465   -0.250125
_Ibirthyr_1936  -0.250249   0.00185886   -134.625      <1e-99  -0.253892   -0.246606
_Ibirthyr_1937  -0.250794   0.00185182   -135.431      <1e-99  -0.254424   -0.247165
_Ibirthyr_1938  -0.248538   0.00174025   -142.818      <1e-99  -0.251949   -0.245127
_Ibirthyr_1939  -0.246448   0.00170952   -144.162      <1e-99  -0.249799   -0.243098
_Ibirthyr_1940  -0.246684   0.00156368   -157.758      <1e-99  -0.249749   -0.243619
_Ibirthyr_1941  -0.242315   0.00161241   -150.281      <1e-99  -0.245476   -0.239155
_Ibirthyr_1942  -0.241501   0.00156544   -154.271      <1e-99  -0.244569   -0.238433
_Ibirthyr_1943  -0.239828   0.00152961   -156.79       <1e-99  -0.242826   -0.23683
_Ibirthyr_1944  -0.23724    0.00156606   -151.488      <1e-99  -0.240309   -0.23417
_Ibirthyr_1945  -0.236017   0.00143321   -164.677      <1e-99  -0.238826   -0.233208
_Ibirthyr_1946  -0.231602   0.00145891   -158.75       <1e-99  -0.234461   -0.228742
_Ibirthyr_1947  -0.229845   0.00147762   -155.551      <1e-99  -0.232741   -0.226949
_Ibirthyr_1948  -0.224731   0.00141518   -158.801      <1e-99  -0.227505   -0.221958
_Ibirthyr_1949  -0.221772   0.00140861   -157.44       <1e-99  -0.224533   -0.219011
_Ibirthyr_1950  -0.222545   0.00130881   -170.035      <1e-99  -0.22511    -0.219979
_Ibirthyr_1951  -0.213693   0.00133732   -159.792      <1e-99  -0.216315   -0.211072
_Ibirthyr_1952  -0.211062   0.00132429   -159.378      <1e-99  -0.213657   -0.208466
_Ibirthyr_1953  -0.206814   0.00127914   -161.682      <1e-99  -0.209321   -0.204307
_Ibirthyr_1954  -0.201267   0.00126368   -159.271      <1e-99  -0.203744   -0.198791
_Ibirthyr_1955  -0.148427   0.00121979   -121.683      <1e-99  -0.150818   -0.146036
_Ibirthyr_1956  -0.148154   0.00121726   -121.711      <1e-99  -0.15054    -0.145768
_Ibirthyr_1957  -0.145949   0.00121263   -120.358      <1e-99  -0.148326   -0.143573
_Ibirthyr_1958  -0.146751   0.00117288   -125.12       <1e-99  -0.14905    -0.144452
_Ibirthyr_1959  -0.141585   0.00116838   -121.18       <1e-99  -0.143875   -0.139295
_Ibirthyr_1960  -0.1423     0.00111462   -127.667      <1e-99  -0.144485   -0.140116
_Ibirthyr_1962  -0.137285   0.001137     -120.744      <1e-99  -0.139514   -0.135057
_Ibirthyr_1963  -0.139382   0.00111215   -125.327      <1e-99  -0.141562   -0.137202
_Ibirthyr_1964  -0.13438    0.00111691   -120.313      <1e-99  -0.136569   -0.13219
_Ibirthyr_1965  -0.140646   0.00108565   -129.55       <1e-99  -0.142773   -0.138518
_Ibirthyr_1966  -0.13261    0.00110388   -120.13       <1e-99  -0.134773   -0.130446
_Ibirthyr_1967  -0.131491   0.00111893   -117.515      <1e-99  -0.133684   -0.129298
_Ibirthyr_1968  -0.129631   0.00109942   -117.909      <1e-99  -0.131786   -0.127476
_Ibirthyr_1969  -0.127519   0.00108838   -117.163      <1e-99  -0.129652   -0.125385
_Ibirthyr_1970  -0.131476   0.00106977   -122.901      <1e-99  -0.133573   -0.12938
_Ibirthyr_1971  -0.126526   0.00108828   -116.262      <1e-99  -0.128659   -0.124393
_Ibirthyr_1972  -0.129166   0.00108103   -119.484      <1e-99  -0.131284   -0.127047
_Ibirthyr_1973  -0.133421   0.00108229   -123.277      <1e-99  -0.135543   -0.1313
_Ibirthyr_1974  -0.134938   0.00109439   -123.3        <1e-99  -0.137083   -0.132793
_Ibirthyr_1975  -0.134999   0.00108154   -124.821      <1e-99  -0.137119   -0.132879
_Ibirthyr_1976  -0.134898   0.00109868   -122.781      <1e-99  -0.137051   -0.132744
_Ibirthyr_1977  -0.135823   0.00111433   -121.887      <1e-99  -0.138007   -0.133639
_Ibirthyr_1978  -0.136901   0.00112702   -121.472      <1e-99  -0.139109   -0.134692
_Ibirthyr_1979  -0.137868   0.00113635   -121.325      <1e-99  -0.140095   -0.135641
_Ibirthyr_1980  -0.13839    0.00111066   -124.602      <1e-99  -0.140567   -0.136214
_Ibirthyr_1981  -0.137445   0.00114061   -120.501      <1e-99  -0.13968    -0.135209
_Ibirthyr_1982  -0.135815   0.00112787   -120.417      <1e-99  -0.138025   -0.133604
_Ibirthyr_1983  -0.0870598  0.00114109    -76.2956     <1e-99  -0.0892963  -0.0848233
_Ibirthyr_1984  -0.0800584  0.00116387    -68.7865     <1e-99  -0.0823396  -0.0777773
_Ibirthyr_1985  -0.0726835  0.00118441    -61.3669     <1e-99  -0.0750049  -0.0703621
_Ibirthyr_1986  -0.059753   0.001222      -48.8977     <1e-99  -0.0621481  -0.0573579
_Ibirthyr_1987  -0.0509262  0.00127088    -40.0715     <1e-99  -0.0534171  -0.0484354
_Ibirthyr_1988  -0.0436099  0.0013144     -33.1785     <1e-99  -0.046186   -0.0410337
_Ibirthyr_1989  -0.0421915  0.00138641    -30.4323     <1e-99  -0.0449088  -0.0394742
_Ibirthyr_1990  -0.039882   0.0014523     -27.4613     <1e-99  -0.0427284  -0.0370355
_Ibirthyr_1991  -0.0316633  0.00158544    -19.9713     <1e-88  -0.0347708  -0.0285559
_Ibirthyr_1992  -0.121891   0.00171674    -71.0015     <1e-99  -0.125256   -0.118527
_Ibirthyr_1993  -0.119903   0.00181338    -66.1214     <1e-99  -0.123457   -0.116349
_Ibirthyr_1994  -0.113538   0.00186542    -60.8645     <1e-99  -0.117194   -0.109882
_Ibirthyr_1995  -0.115299   0.00193508    -59.5836     <1e-99  -0.119092   -0.111507
_Ibirthyr_1996  -0.117205   0.00204327    -57.3618     <1e-99  -0.12121    -0.113201
_Ibirthyr_1997  -0.117898   0.00227367    -51.8537     <1e-99  -0.122354   -0.113442
_Ibirthyr_1998  -0.114736   0.00272781    -42.0616     <1e-99  -0.120082   -0.10939
_Ibirthyr_1999  -0.120351   0.00364711    -32.999      <1e-99  -0.127499   -0.113203
_Ibirthyr_2000  -0.0949544  0.00852195    -11.1423     <1e-28  -0.111657   -0.0782517
(Intercept)      0.286457   0.000805943   355.43       <1e-99   0.284877    0.288036
=====================================================================================

julia> plot(coef(m)[2:end])

[matthieugomez]

Sorry just to recap, is this a good summary of what you obtain?

using CSV, DataFrames, LinearAlgebra, FixedEffectModels, Plots
df = CSV.read(".../FEbug.csv", DataFrame)
f = @formula(part ~ _Ibirthyr_1907+_Ibirthyr_1908+_Ibirthyr_1909+_Ibirthyr_1910+_Ibirthyr_1911+_Ibirthyr_1912+_Ibirthyr_1913+_Ibirthyr_1914+
     _Ibirthyr_1915+_Ibirthyr_1916+_Ibirthyr_1917+_Ibirthyr_1918+_Ibirthyr_1919+_Ibirthyr_1920+_Ibirthyr_1921+_Ibirthyr_1922+_Ibirthyr_1923+_Ibirthyr_1924+
     _Ibirthyr_1925+_Ibirthyr_1926+_Ibirthyr_1927+_Ibirthyr_1928+_Ibirthyr_1929+_Ibirthyr_1930+_Ibirthyr_1931+_Ibirthyr_1932+_Ibirthyr_1933+_Ibirthyr_1934+
     _Ibirthyr_1935+_Ibirthyr_1936+_Ibirthyr_1937+_Ibirthyr_1938+_Ibirthyr_1939+_Ibirthyr_1940+_Ibirthyr_1941+_Ibirthyr_1942+_Ibirthyr_1943+_Ibirthyr_1944+
     _Ibirthyr_1945+_Ibirthyr_1946+_Ibirthyr_1947+_Ibirthyr_1948+_Ibirthyr_1949+_Ibirthyr_1950+_Ibirthyr_1951+_Ibirthyr_1952+_Ibirthyr_1953+_Ibirthyr_1954+
     _Ibirthyr_1955+_Ibirthyr_1956+_Ibirthyr_1957+_Ibirthyr_1958+_Ibirthyr_1959+_Ibirthyr_1960+_Ibirthyr_1962+_Ibirthyr_1963+_Ibirthyr_1964+
     _Ibirthyr_1965+_Ibirthyr_1966+_Ibirthyr_1967+_Ibirthyr_1968+_Ibirthyr_1969+_Ibirthyr_1970+_Ibirthyr_1971+_Ibirthyr_1972+_Ibirthyr_1973+_Ibirthyr_1974+
     _Ibirthyr_1975+_Ibirthyr_1976+_Ibirthyr_1977+_Ibirthyr_1978+_Ibirthyr_1979+_Ibirthyr_1980+_Ibirthyr_1981+_Ibirthyr_1982+_Ibirthyr_1983+_Ibirthyr_1984+
     _Ibirthyr_1985+_Ibirthyr_1986+_Ibirthyr_1987+_Ibirthyr_1988+_Ibirthyr_1989+_Ibirthyr_1990+_Ibirthyr_1991+_Ibirthyr_1992+_Ibirthyr_1993+_Ibirthyr_1994+
     _Ibirthyr_1995+_Ibirthyr_1996+_Ibirthyr_1997+_Ibirthyr_1998+_Ibirthyr_1999+_Ibirthyr_2000)

X = modelmatrix(f, df) 
y = response(f, df) 
# works everywhere
coef_qr = qr(X) \ y
coef_chol = cholesky!(X' * X) \ (X' * y)
# works correctly with Mac Intel or M2 but fails on Linux
reg(df, f)

[droodman]

Yes! Would it make sense for me to try the QR decomposition here ? https://github.com/FixedEffects/FixedEffectModels.jl/blob/eda1b0874d615e10493ab29ee62a2f358f6a5bc8/src/fit.jl#L397

On Windows...

[matthieugomez]

Sure that would be helful. Could you also try this (note that I add an explicit intercept in the formula, which matters when using modelmatrix directly)

using CSV, DataFrames, LinearAlgebra, FixedEffectModels, Plots
df = CSV.read(".../FEbug.csv", DataFrame)
f = @formula(part ~ 1 + _Ibirthyr_1907+_Ibirthyr_1908+_Ibirthyr_1909+_Ibirthyr_1910+_Ibirthyr_1911+_Ibirthyr_1912+_Ibirthyr_1913+_Ibirthyr_1914+
     _Ibirthyr_1915+_Ibirthyr_1916+_Ibirthyr_1917+_Ibirthyr_1918+_Ibirthyr_1919+_Ibirthyr_1920+_Ibirthyr_1921+_Ibirthyr_1922+_Ibirthyr_1923+_Ibirthyr_1924+
     _Ibirthyr_1925+_Ibirthyr_1926+_Ibirthyr_1927+_Ibirthyr_1928+_Ibirthyr_1929+_Ibirthyr_1930+_Ibirthyr_1931+_Ibirthyr_1932+_Ibirthyr_1933+_Ibirthyr_1934+
     _Ibirthyr_1935+_Ibirthyr_1936+_Ibirthyr_1937+_Ibirthyr_1938+_Ibirthyr_1939+_Ibirthyr_1940+_Ibirthyr_1941+_Ibirthyr_1942+_Ibirthyr_1943+_Ibirthyr_1944+
     _Ibirthyr_1945+_Ibirthyr_1946+_Ibirthyr_1947+_Ibirthyr_1948+_Ibirthyr_1949+_Ibirthyr_1950+_Ibirthyr_1951+_Ibirthyr_1952+_Ibirthyr_1953+_Ibirthyr_1954+
     _Ibirthyr_1955+_Ibirthyr_1956+_Ibirthyr_1957+_Ibirthyr_1958+_Ibirthyr_1959+_Ibirthyr_1960+_Ibirthyr_1962+_Ibirthyr_1963+_Ibirthyr_1964+
     _Ibirthyr_1965+_Ibirthyr_1966+_Ibirthyr_1967+_Ibirthyr_1968+_Ibirthyr_1969+_Ibirthyr_1970+_Ibirthyr_1971+_Ibirthyr_1972+_Ibirthyr_1973+_Ibirthyr_1974+
     _Ibirthyr_1975+_Ibirthyr_1976+_Ibirthyr_1977+_Ibirthyr_1978+_Ibirthyr_1979+_Ibirthyr_1980+_Ibirthyr_1981+_Ibirthyr_1982+_Ibirthyr_1983+_Ibirthyr_1984+
     _Ibirthyr_1985+_Ibirthyr_1986+_Ibirthyr_1987+_Ibirthyr_1988+_Ibirthyr_1989+_Ibirthyr_1990+_Ibirthyr_1991+_Ibirthyr_1992+_Ibirthyr_1993+_Ibirthyr_1994+
     _Ibirthyr_1995+_Ibirthyr_1996+_Ibirthyr_1997+_Ibirthyr_1998+_Ibirthyr_1999+_Ibirthyr_2000)
X = modelmatrix(f, df) 
y = response(f, df) 
coef_qr = qr(X) \ y
coef_chol = cholesky!(X' * X) \ (X' * y)
coef_chol2 = cholesky!(Symmetric(X' * X)) \ (X' * y)

[droodman]

Those all look good in Windows:

plot([coef_qr coef_chol.+.01 coef_chol2.+.02][2:end,:])

[droodman]

Separately, yes, switching to the QR decomposition fixes the problem! I just replaced ldiv!(crossx, Xhat'y) with qr(Xhat) \ y. This confirms my initial guess that this is a numerical issue related to matrix inversion.

Sources I looked at a few days ago were nearly unanimous in saying that software should use the QR, that it's probably what's used in Stata and core R packages. E.g., "QR decomposition is indeed the standard way of solving least square problems effectively with much simplicity"

So I think this package should have an option for QR at least in the post-FE-removal stage, and I would have reghdfejl call it by default.

Relatedly, see this discussion of argument naming. method is getting used different ways in different places.

[matthieugomez]

I don't understand why coef_chol2 gives you different results than reg — they should literally give you the same thing. I am also not convinced that QR can matter relative to Cholesky for economically meaningful regressions (especially because your coefficients are estimated with small SEs). I think something else may be happening.

Just before https://github.com/FixedEffects/FixedEffectModels.jl/blob/master/src/fit.jl#L397 I have added the lines

    f2 = @formula(part ~ 1 + _Ibirthyr_1907+_Ibirthyr_1908+_Ibirthyr_1909+_Ibirthyr_1910+_Ibirthyr_1911+_Ibirthyr_1912+_Ibirthyr_1913+_Ibirthyr_1914+
    _Ibirthyr_1915+_Ibirthyr_1916+_Ibirthyr_1917+_Ibirthyr_1918+_Ibirthyr_1919+_Ibirthyr_1920+_Ibirthyr_1921+_Ibirthyr_1922+_Ibirthyr_1923+_Ibirthyr_1924+
    _Ibirthyr_1925+_Ibirthyr_1926+_Ibirthyr_1927+_Ibirthyr_1928+_Ibirthyr_1929+_Ibirthyr_1930+_Ibirthyr_1931+_Ibirthyr_1932+_Ibirthyr_1933+_Ibirthyr_1934+
    _Ibirthyr_1935+_Ibirthyr_1936+_Ibirthyr_1937+_Ibirthyr_1938+_Ibirthyr_1939+_Ibirthyr_1940+_Ibirthyr_1941+_Ibirthyr_1942+_Ibirthyr_1943+_Ibirthyr_1944+
    _Ibirthyr_1945+_Ibirthyr_1946+_Ibirthyr_1947+_Ibirthyr_1948+_Ibirthyr_1949+_Ibirthyr_1950+_Ibirthyr_1951+_Ibirthyr_1952+_Ibirthyr_1953+_Ibirthyr_1954+
    _Ibirthyr_1955+_Ibirthyr_1956+_Ibirthyr_1957+_Ibirthyr_1958+_Ibirthyr_1959+_Ibirthyr_1960+_Ibirthyr_1962+_Ibirthyr_1963+_Ibirthyr_1964+
    _Ibirthyr_1965+_Ibirthyr_1966+_Ibirthyr_1967+_Ibirthyr_1968+_Ibirthyr_1969+_Ibirthyr_1970+_Ibirthyr_1971+_Ibirthyr_1972+_Ibirthyr_1973+_Ibirthyr_1974+
    _Ibirthyr_1975+_Ibirthyr_1976+_Ibirthyr_1977+_Ibirthyr_1978+_Ibirthyr_1979+_Ibirthyr_1980+_Ibirthyr_1981+_Ibirthyr_1982+_Ibirthyr_1983+_Ibirthyr_1984+
    _Ibirthyr_1985+_Ibirthyr_1986+_Ibirthyr_1987+_Ibirthyr_1988+_Ibirthyr_1989+_Ibirthyr_1990+_Ibirthyr_1991+_Ibirthyr_1992+_Ibirthyr_1993+_Ibirthyr_1994+
    _Ibirthyr_1995+_Ibirthyr_1996+_Ibirthyr_1997+_Ibirthyr_1998+_Ibirthyr_1999+_Ibirthyr_2000)
   Xhat2 = modelmatrix(f2, df)
   y2 = response(f2, df)
   coef2 = cholesky!(Xhat2' * Xhat2) \ (Xhat2' * y2)

   crossx = cholesky!(Symmetric(Xhat'Xhat))
   coef = ldiv!(crossx, Xhat'y)

   @assert sum(basis_coef) == size(Xhat2, 2)
   @assert Xhat == Xhat2
   @assert y == y2
   @assert coef == coef2
   @show "success"

When I run reg(df, f), the four assertions pass for me and "success" gets printed. They should not pass for you. Could you try to run it? Understanding which one fails would help us pinpoint the issue.

[droodman]

I get:

ERROR: AssertionError: coef == coef2

[matthieugomez]

Thanks. Can you try again after replacing the line coef = ldiv!(crossx, Xhat'y) by coef = crossx \ (Xhat' * y)? If that still fails, I'm out of idea on what is driving the issue.

[droodman]

That still fails.

So I drilled down some with the debugger. It seems to have to do with multiplication by a SentinalArrays.ChainedVector. Look at the items in the Watch pane on the left.

[droodman]

The upshot here is that replacing y2 in your test code with collect(y2) in order to convert it to Vector{Float64} lets it pass all the @asserts on my machine. I'm not sure what's going on with that. I'm also not sure it's relevant because the actual variable in the code, y, is already Vector{Float64}.

Meanwhile, using collect(y2) doesn't solve my original problem. What does solve the problem is qr().

Here's a complicated screenshot. The last item in the Watch list pulls together coef2 as defined in what is now line 408, involving collect(y2), the computation of coef in the actual code as ldiv!(crossx, Xhat'y), and the QR-based computation qr(Xhat) \ y. The call to show() in the last watch item dumps the full results into the debugger output window below. I copied the resulting 3-column matrix into another Julia session and plotted the three columns, as shown on the right. Only the QR-based result is correct.

In my experience, this is how it is with numerical issues in linear algebra routines. You drill and drill to find where the computation goes off the rails, and it turns out it's not really a logical bug, just imprecision in the matrix inversion/factorization that you previously didn't need to think about. Then you have to confront a speed-precision tradeoff.

[matthieugomez]

Thanks a lot. I am open to using QR in some cases but, if you look at the conditioning number of X, it is not particularly different than the conditioning number of any other matrix (which, I think, is the same as saying that the system is not close to being collinear). So I'm unsure why this happens and how one could predict this issue.

[matthieugomez]

I think that's what SAS does btw (and Stata is probably similar): https://blogs.sas.com/content/iml/2018/04/18/sweep-operator-sas.html, https://blogs.sas.com/content/iml/2021/07/14/performance-ls-regression.html

[droodman]

Yes, this other SAS blog post about the QR says "most SAS regression procedures use the sweep operator to construct least-squares solutions." OTOH, it starts out with "A SAS programmer recently mentioned that some open-source software uses the QR algorithm to solve least-squares regression problems" and we can guess what that refers to. And indeed it looks like R uses the QR (this, this).

I don't know what Stata does, but evidently it sometimes uses something more precise than Cholesky. It looks like Numpy uses the SVD, which I read is even slower and more precise than the QR (this, this).

I'm getting a condition number of 346 for Xhat, thus 346^2 for crossx (which is faster to compute especially since crossx will still be computed anyway, right?). Isn't that pretty big?

I can see a few options:

1. Leave it up to the user, which seems to be where lm() is headed.
2. Always use the QR, apparently like R.
3. Use some threshold condition number on crossx to switch between Cholesky and QR. Matlab does something like this (see the flow chart).

[s3alfisc]

Hi both, just for reference - in R, lm() calls .lm.fit(), which uses the QR decomposition. Here's the documentation: link.

lm.fit (x, y,    offset = NULL, method = "qr", tol = 1e-7,
       singular.ok = TRUE, ...)

method: currently, only method = "qr" is supported.
tol: tolerance for the [qr](https://stat.ethz.ch/R-manual/R-devel/library/Matrix/html/qr-methods.html) decomposition. Default is 1e-7.

[ericqu]

Apologies as it is a bit off-topic, I was wondering if the data could be shared or tested with the LinearRegressionKit.jl, I would like to know out of curiosity if this data set is an issue with the sweep operator.

Personally, I am not a big fan of a threshold on the condition number as it only considers the Xs, in my view, this leaves too much uncertainty to leave the Y out.

My preference would be to have a QR by default, that could be changed by the user.

[droodman]

@ericqu I'd be happy to send you the data if you provide a more private way to do so.

[matthieugomez]

I am still not convinced that QR is needed in a setup where one obtains small SEs for all coefficients. Relatedly, to compute standard errors, one needs to compute inv(X'X), even though the computer science / numerical methods people would say it's a horrible thing to do. In any case, having QR as a default would be bad as it would slow down everything (this is one reason why R is much slower than Stata for standard OLS regressions with big data).

I have pushed a version that solves for coefficients using sweep operations, using the same algorithm I was already using to test for collinearity. Hopefully that solves your problem. https://github.com/FixedEffects/FixedEffectModels.jl/pull/253

[droodman]

Thank you @matthieugomez. I upgraded to 1.10.0 but unfortunately, the results have not changed.

If you are getting tired of this, I can just fork the project, tweak it, and point reghdfejl to that.

[matthieugomez]

Well I would prefer to understand the difference, but it's hard to me to solve something that does not happen on my computer (I have double checked and I get the exact same coefficients in Stata and Julia). If you time Stata, you can see that it returns a result faster than what a QR factorization would take, so I really don't think they use QR.

[droodman]

reg(df,f; nthreads=6) takes 48.640588 seconds. In Stata, after set processors 6, the regression takes 4.85 seconds. (Not a typo.)

[matthieugomez]

Part of it is https://github.com/JuliaStats/StatsModels.jl/issues/222. You should be able to avoid it by passing factors directly rather than the continuous variables.

In any case, the point is that 4.85 seconds is less than what my computer takes to compute a QR factorization.

[droodman]

Interesting. Yeah on my machine the QR takes ~8.1 seconds.

[matthieugomez]

That still fails.

So I drilled down some with the debugger. It seems to have to do with multiplication by a SentinalArrays.ChainedVector. Look at the items in the Watch pane on the left.

I am looking at this screenshot. Am I correct in understanding that you are getting (Xhat' y)[end] = 320.0 but (Xhat[:, end]' y) = 477.0? If it is the case (which seems wild to me), can you get it with this simpler code (i.e. not involving reg):

f = @formula(part ~ 1 + _Ibirthyr_1907+_Ibirthyr_1908+_Ibirthyr_1909+_Ibirthyr_1910+_Ibirthyr_1911+_Ibirthyr_1912+_Ibirthyr_1913+_Ibirthyr_1914+
    _Ibirthyr_1915+_Ibirthyr_1916+_Ibirthyr_1917+_Ibirthyr_1918+_Ibirthyr_1919+_Ibirthyr_1920+_Ibirthyr_1921+_Ibirthyr_1922+_Ibirthyr_1923+_Ibirthyr_1924+
    _Ibirthyr_1925+_Ibirthyr_1926+_Ibirthyr_1927+_Ibirthyr_1928+_Ibirthyr_1929+_Ibirthyr_1930+_Ibirthyr_1931+_Ibirthyr_1932+_Ibirthyr_1933+_Ibirthyr_1934+
    _Ibirthyr_1935+_Ibirthyr_1936+_Ibirthyr_1937+_Ibirthyr_1938+_Ibirthyr_1939+_Ibirthyr_1940+_Ibirthyr_1941+_Ibirthyr_1942+_Ibirthyr_1943+_Ibirthyr_1944+
    _Ibirthyr_1945+_Ibirthyr_1946+_Ibirthyr_1947+_Ibirthyr_1948+_Ibirthyr_1949+_Ibirthyr_1950+_Ibirthyr_1951+_Ibirthyr_1952+_Ibirthyr_1953+_Ibirthyr_1954+
    _Ibirthyr_1955+_Ibirthyr_1956+_Ibirthyr_1957+_Ibirthyr_1958+_Ibirthyr_1959+_Ibirthyr_1960+_Ibirthyr_1962+_Ibirthyr_1963+_Ibirthyr_1964+
    _Ibirthyr_1965+_Ibirthyr_1966+_Ibirthyr_1967+_Ibirthyr_1968+_Ibirthyr_1969+_Ibirthyr_1970+_Ibirthyr_1971+_Ibirthyr_1972+_Ibirthyr_1973+_Ibirthyr_1974+
    _Ibirthyr_1975+_Ibirthyr_1976+_Ibirthyr_1977+_Ibirthyr_1978+_Ibirthyr_1979+_Ibirthyr_1980+_Ibirthyr_1981+_Ibirthyr_1982+_Ibirthyr_1983+_Ibirthyr_1984+
    _Ibirthyr_1985+_Ibirthyr_1986+_Ibirthyr_1987+_Ibirthyr_1988+_Ibirthyr_1989+_Ibirthyr_1990+_Ibirthyr_1991+_Ibirthyr_1992+_Ibirthyr_1993+_Ibirthyr_1994+
    _Ibirthyr_1995+_Ibirthyr_1996+_Ibirthyr_1997+_Ibirthyr_1998+_Ibirthyr_1999+_Ibirthyr_2000)
Xhat = modelmatrix(f, df)
y_sentinel = response(f, df)
y_vector = convert(Vector{Float64}, y_sentinel)
@show (Xhat' * y_vector)[end]
@show Xhat[:, end]' * y_vector

[droodman]

Yes you are correct. I don't know if calls for a new issue on that package.

julia> @show (Xhat' * y_vector)[end]
(Xhat' * y_vector)[end] = 320.0
320.0

julia> @show Xhat[:, end]' * y_vector
(Xhat[:, end])' * y_vector = 477.0
477.0

[matthieugomez]

Well so that explains the problem — there is an issue with your version of Julia or with your underlying BLAS / MKL. Both Xhat and y_vectors are vectors/matrices of 0 and 1 so there should not be any difference between the two numbers. Have you tried with Julia 1.10? Can you try to find a version of this issue with simpler arrays so that you can open an issue on the Julia repository? (maybe Xhat= Xhat[:, end:end] and y = ones(length(y_vector))

[grantmcdermott]

A few unsolicited, but hopefully still helpful thoughts.

1. The potential instability of Cholesky is well documented and, IMO, is good cause for defaulting to QR. Luke Tierney (R Core member) has a nice lecture series on this sort of thing here: http://homepage.divms.uiowa.edu/~luke/classes/STAT7400-2023/
2. Speaking of options, in R you can invoke different decomposition strategies through different packages. The (very fast) collapse package offers a flm convenience frontend to a bunch of these. All of which is to say that, if you can, doing something similar and offering your users options isn't a bad compromise.
3. Optimised BLAS introduces a whole new can of worms, due to competing parallelism / race conditions, among other things. Relying on MKL, in particular, is asking for trouble. Having been bit by wicked esoteric bugs in the past, I can only pass on the advice that was given to me: Friends don't let friends use MKL! (Didn't MacKinnon have a working paper demonstrating the instability of Stata results under different BLAS regimes, or am I misremembering?)

[matthieugomez]

@grantmcdermott I could be convinced to add an option for QR but the first step would be for someone to come up with an example where the difference matters relative to the standard errors of the coefficients. The issue right now is a bug, not of instability of Cholesky versus QR.

[droodman]

@matthieugomez I'm not sure what you mean by "the difference matters relative to the standard errors." In the motivating example here, the discontinuity is very large relative to the standard errors.

[matthieugomez]

Your example has nothing to do with the instability of Cholesky versus QR. It's about the fact that your computer returns a wrong result for X'y (X and y are just arrays of 0-1, so it is not related to numerical instability). You should make sure you are on the latest version of Julia (1.9.3), simplify the example so that you can share it, and open an issue mentioning your OS and the result of LinearAlgebra.BLAS.get_config() on https://github.com/JuliaLang/julia.

Here is some code that might help you get a reproductible example (it produces a vector y and a matrix X similar to the ones obtained from your data).

N = 10_000_000
X = hcat(ones(N), rand(N, 100) .>= 0.95)
y = Vector{Float64}(rand(N) .>= 0.9)
@assert (X'y)[end] ≈ X[:, end]'y

[droodman]

Submitted. I had thought that conversion of that Sentinel vector to a regular vector had eliminated the effect on the computations in this package. It did so in the debugging watches, but that's not quite the same thing.

Ironically @grantmcdermott, switching from OpenBLAS to MKL eliminates the problem(!).

What I don't understand is why sticking in a simple coef = qr(X) \ y also did. I guess the ldiv operator is not having this weird math issue.

[matthieugomez]

Thanks. If you replace y by view(y, :), the issue does not happen anymore right? If so, I will do that until this OpenBlas issue gets fixed.

[droodman]

I tried adding y = view(y,:) in fit.jl after the line constructing crossx. That did not help. But maybe that's not what you mean.

Following this comment I added using BLISBLAS in reghdfejl before calling FixedEffectModels. It works!

[matthieugomez]

Great, but I'd like a solution even for people that don't use reghdfejl. What if you replace X' y by X' reshape(y, length(y), 1) in your reproductible example?

[droodman]

Yes, that works too! But now I don't trust OpenBLAS. What if you put using BLISBLAS in FixedEffectModels, at least for x86_64 architecture?

[matthieugomez]

Well BLISBLAS can probably also create other types of bugs in the future... Out of curiosity, were you able to replicate this bug with other matrices? I'm wondering how common it is in the wild

[matthieugomez]

Solved by https://github.com/FixedEffects/FixedEffectModels.jl/pull/255