# Regression: fixed effects + quantiles

**URL:** <https://discourse.julialang.org/t/regression-fixed-effects-quantiles/45427>\
**Category:** Statistics\
**Created:** [August 23, 2020, 7:39pm UTC](https://discourse.julialang.org/t/regression-fixed-effects-quantiles/45427 "2020-08-23T19:39:49Z")\
**Posts on this page:** 17\
**Page:** 1

<div class="post-metadata">

**Author:** ![danicaratelli](https://avatars.discourse-cdn.com/v4/letter/d/90db22/32.png) [@danicaratelli](https://discourse.julialang.org/u/danicaratelli)\
**Post date:** [August 23, 2020, 7:39pm UTC](https://discourse.julialang.org/t/regression-fixed-effects-quantiles/45427/1 "2020-08-23T19:39:49Z")

</div>

I am trying to run a **quantile regression** with **fixed effects**. Does anyone know if a package for this exists?

Thanks!

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [August 23, 2020, 8:16pm UTC](https://discourse.julialang.org/t/regression-fixed-effects-quantiles/45427/2 "2020-08-23T20:16:21Z")

</div>

We have one package for fixed effects:  
[https://github.com/FixedEffects](https://github.com/FixedEffects)

A separate for quantile regression:  
[https://github.com/pkofod/QuantileRegressions.jl](https://github.com/pkofod/QuantileRegressions.jl)

Currently no package that automates both together.

I would create fixed effects by hand & then run a quantile regression.

@pkofod  
@matthieu

---

<div class="post-metadata">

**Author:** ![danicaratelli](https://avatars.discourse-cdn.com/v4/letter/d/90db22/32.png) [@danicaratelli](https://discourse.julialang.org/u/danicaratelli)\
**Post date:** [August 23, 2020, 8:17pm UTC](https://discourse.julialang.org/t/regression-fixed-effects-quantiles/45427/3 "2020-08-23T20:17:22Z")

</div>

Yeah, that’s what I was doing. I guess better than nothing. Thank you!

---

<div class="post-metadata">

**Author:** ![tlienart](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tlienart/32/7640_2.png) [@tlienart](https://discourse.julialang.org/u/tlienart)\
**Post date:** [August 23, 2020, 9:08pm UTC](https://discourse.julialang.org/t/regression-fixed-effects-quantiles/45427/4 "2020-08-23T21:08:00Z")

</div>

btw [MLJLinearModels](https://github.com/alan-turing-institute/MLJLinearModels.jl) also supports QuantileRegression an is generally faster than the package cited above (actually it’s also faster than the specialised R package 🤷‍♂️) that being said it’s probably more barebones, depends what you want.

---

<div class="post-metadata">

**Author:** ![pkofod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pkofod/32/2179_2.png) [@pkofod](https://discourse.julialang.org/u/pkofod)\
**Post date:** [August 24, 2020, 7:13am UTC](https://discourse.julialang.org/t/regression-fixed-effects-quantiles/45427/5 "2020-08-24T07:13:55Z")

</div>

Do you have the benchmarks against quantreg? I’m curious about the design matrices you used.

---

<div class="post-metadata">

**Author:** ![pkofod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pkofod/32/2179_2.png) [@pkofod](https://discourse.julialang.org/u/pkofod)\
**Post date:** [August 24, 2020, 7:16am UTC](https://discourse.julialang.org/t/regression-fixed-effects-quantiles/45427/6 "2020-08-24T07:16:57Z")

</div>

Without using a sparse designmatrix that seems very expensive 😬 In what sense are you including a fixed effect here? How are you keeping it constant across fixed effects in QuantileRegressions.jl?

---

<div class="post-metadata">

**Author:** ![Paul\_Soderlind](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paul_soderlind/32/1753_2.png) [@Paul\_Soderlind](https://discourse.julialang.org/u/Paul_Soderlind)\
**Post date:** [August 24, 2020, 7:27am UTC](https://discourse.julialang.org/t/regression-fixed-effects-quantiles/45427/7 "2020-08-24T07:27:12Z")

</div>

> I would create fixed effects by hand & then run a quantile regression.

is a good way forward.  
The first part is just creating the usual (ystar,xstar) by individual de-meaning. Then run a quantile regression on that. For instance, here is a working code

```julia
"""
    QuantRegrIRLSPs(y,x::Matrix,q=0.5;prec=1e-8,epsu=1e-6,maxiter=1000)

Quantile regression by using an IRLS algorithm.

Estimate a quantile regression, quantile q. The outputs are the point estimates
and three different variance-covariance matrices of the estimates.

# Input
- `y::Vector`: T vector, dependent variable
- `x::Matrix`: TXK, regressors (including any constant)
- `q::Number`: quantile to estimate at, 0<q<1
- `prec::Float64`: convergence criterion, 1e-8
- `epsu::Float64`: lower bound on 1/weight, 1e-6
- `maxiter::Int`: maximum number of iterations, 1000

# Output
- `theta::Vector`: K vector, estimated coefficients
- `vcv::Matrix`: KxK, traditional covariance matrix
- `vcv2::Matrix`: KxK, Powell (1991) covariance matrix
- `vcv3::Matrix`: KxK, Powell (1991) covariance matrix, uniform

"""
function QuantRegrIRLSPs(y,x::Matrix,q=0.5;prec=1e-8,epsu=1e-6,maxiter=1000)

  (T,K) = size(x)
  xw = copy(x)

  (b_old,b,u,iter) = (zeros(K),fill(1e+6,K) .+ prec,zeros(T),0)

  while maximum(abs,b - b_old) > prec
    copyto!(b_old, b)
    b .= (xw'*x)\(xw'*y)
    u .= y - x*b
    #u .= ifelse.(u.>0,1-q,q).*abs.(u) #as in Python code, divide loss fn by q(1-q) to get it
    u .= ifelse.(u.>0,1/q,1/(1-q)).*abs.(u) #abs(u)/q if u>0, abs(u)/(1-q) if u<0
    u .= max.(u,epsu) #not smaller than epsu
    xw .= x./u
    iter = iter + 1
    if iter > maxiter
      @warn("$iter > maxiter")
      b = NaN
      break
    end
  end

  res = y - x*b

  D = x'x/T
  h = 1.06*std(res)/T^0.2 #Silverman (1986) recommendation
  fx = exp.(-0.5*((res/h).^2))./(h*sqrt(2*pi)) #Green 7th ed, using N() kernel
  f0 = mean(fx)
  C = f0*x'x/T
  C_1 = inv(C)
  vcv = q*(1-q)*C_1*D*C_1/T #variance-covariance matrix

  C = (fx.*x)'x/T #Wooldrige 2dn ed, Powell 1991
  C_1 = inv(C) #but with Gaussian kernel
  vcv2 = q*(1-q)*C_1*D*C_1/T #caputures (x,res) dependence

  fx = (abs.(res) .<= h/2)/h #Wooldrige 2nd ed, Powell 1991
  C = (fx.*x)'x/T #uniform kernel
  C_1 = inv(C)
  vcv3 = q*(1-q)*C_1*D*C_1/T

  return b, vcv, vcv2, vcv3

end

```

---

<div class="post-metadata">

**Author:** ![danicaratelli](https://avatars.discourse-cdn.com/v4/letter/d/90db22/32.png) [@danicaratelli](https://discourse.julialang.org/u/danicaratelli)\
**Post date:** [August 24, 2020, 2:32pm UTC](https://discourse.julialang.org/t/regression-fixed-effects-quantiles/45427/8 "2020-08-24T14:32:28Z")

</div>

I first modify the data by manually including FEs.  
Then I write the formula corresponding to the regression with FE.  
This is included here:

```julia
#this function runs a quantile regression with fixed effects
function QregFE(Xdata::DataFrame,dep_var::String,indep_vars::Array{String,1},FEs::Array{String,1},icon=true)
    Xdata = copy(Xdata);
    nr = size(Xdata,1);
    indep_vars
    #including FE
    for ife = 1:length(FEs)
        fe_list = sort(unique(Xdata[!,Symbol(FEs[ife])]));
        nFE = length(fe_list);
        for i=1:nFE-1
            Xdata[!,Symbol(FEs[ife]*"_"*string(i))] = Bool.(zeros(nr));
            Xdata[findall(Xdata[!,Symbol(FEs[ife])].==fe_list[i]),Symbol(FEs[ife]*"_"*string(i))] .=true;
            indep_vars = vcat(indep_vars,FEs[ife]*"_"*string(i));
        end
    end
    #creating formula using the strings
    reg_string = "@formula(";
    reg_string = reg_string*" "*dep_var*" ~ ";
    if icon; reg_string = reg_string*" 1 + "; end;
    reg_string = reg_string*join(indep_vars, '+');
    reg_string = reg_string*")";
    fm = eval(Meta.parse(reg_string));
    return fm, Xdata;
end

```

Finally, I use your function `qreg`. This is my approach.

---

<div class="post-metadata">

**Author:** ![pdeffebach](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pdeffebach/32/10320_2.png) [@pdeffebach](https://discourse.julialang.org/u/pdeffebach)\
**Post date:** [August 24, 2020, 2:34pm UTC](https://discourse.julialang.org/t/regression-fixed-effects-quantiles/45427/9 "2020-08-24T14:34:36Z")

</div>

The formula construction part is much more complicated than it needs to be. You can construct a formula programatically pretty easily.

[https://juliastats.org/StatsModels.jl/stable/formula/#Constructing-a-formula-programmatically-1](https://juliastats.org/StatsModels.jl/stable/formula/#Constructing-a-formula-programmatically-1)

---

<div class="post-metadata">

**Author:** ![tlienart](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tlienart/32/7640_2.png) [@tlienart](https://discourse.julialang.org/u/tlienart)\
**Post date:** [August 24, 2020, 3:05pm UTC](https://discourse.julialang.org/t/regression-fixed-effects-quantiles/45427/10 "2020-08-24T15:05:15Z")

</div>

I used the stuff that’s used in the R package quantreg; but I’d be interested in building a more general test suite for generalised regression; btw speed is usually not a key element for that kind of stuff unless people have “big” data but iirc your package makes use of RWLS (I’m typing on my phone here and may be remembering wrong!) whereas the package I linked to offers a few solvers, some which can be faster.

Anyway this is maybe for another discussion; I’m generally interested in building a generalised regression toolbox in Julia that makes good use of tools that are already available like Optim and IterativeSolvers and can offer all the standard penalised regressions out there that would beat standard regression packages in R and Matlab etc. MLJLM is one (reasonably serious) attempt at that but I’d be happy to do something with more people. (GLMs is a bit too constrained in its design so it would have to be something else)

There were a few other attempts like the defunct Regressions.jl and mostly defunct SparseRegressions and there’s probably others out there…

---

<div class="post-metadata">

**Author:** ![pkofod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pkofod/32/2179_2.png) [@pkofod](https://discourse.julialang.org/u/pkofod)\
**Post date:** [August 25, 2020, 11:26am UTC](https://discourse.julialang.org/t/regression-fixed-effects-quantiles/45427/11 "2020-08-25T11:26:31Z")

</div>

> [@tlienart](#):
>
> btw speed is usually not a key element for that kind of stuff unless people have “big” data but iirc your package makes use of RWLS (I’m typing on my phone here and may be remembering wrong!) whereas the package I linked to offers a few solvers, some which can be faster.

it should default to the interior point solver but does also have iteratively reweighted least squares. The first hopefully matches quantregs default as it was programmed up according to the article (The _Gaussian Hare_ and the Laplacian Tortoise).

Anyway yeah, it’s all about use-cases. The usecase in OP is quite performance sensitive though, and does require sparse liinear algebra to be acceptable. I know because I’ve estimated these models quite a few times with quite a few individuals. Using the penalized version of Koenker and also the plug-in estimator by (striking out on the author here).

---

<div class="post-metadata">

**Author:** ![tlienart](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tlienart/32/7640_2.png) [@tlienart](https://discourse.julialang.org/u/tlienart)\
**Post date:** [August 25, 2020, 4:47pm UTC](https://discourse.julialang.org/t/regression-fixed-effects-quantiles/45427/12 "2020-08-25T16:47:43Z")

</div>

So it turns out that in some case you can do a projected gradient descent and it’s significantly faster than the IP (let alone IRLS) to get to essentially the same coefs; but I’m aware these are always tricky claims to make and I’d be happy to investigate more in depth though maybe a different topic would be more appropriate 🙂

---

<div class="post-metadata">

**Author:** ![pkofod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pkofod/32/2179_2.png) [@pkofod](https://discourse.julialang.org/u/pkofod)\
**Post date:** [August 26, 2020, 9:02am UTC](https://discourse.julialang.org/t/regression-fixed-effects-quantiles/45427/13 "2020-08-26T09:02:18Z")

</div>

Interesting, yeah. Happy to add it directly to QuantileRegressions.jl. You could also consider benchmarking against some of the other non-default solvers in quantreg such as Barrodal and Roberts. I don’t know of the proximal versions are in quantreg but there are some variations out there. Anyway, this is getting a bit off-topic 🙂

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [October 12, 2020, 5:34pm UTC](https://discourse.julialang.org/t/regression-fixed-effects-quantiles/45427/14 "2020-10-12T17:34:07Z")

</div>

My experience (& this discussion) makes me feel that Julia stats packages can work better together (be more modular). Currently if I want fixed effects (outside of [FixedEffectModels](https://github.com/FixedEffects/FixedEffectModels.jl).jl) I have to create them manually (by hand).

It would be great if a general package such as `StatsModels.jl` could automatically create fixed effects. Then these can be used w/ other stats packages in the ecosystem (such as `QuantileRegressions.jl`).

---

<div class="post-metadata">

**Author:** ![pdeffebach](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pdeffebach/32/10320_2.png) [@pdeffebach](https://discourse.julialang.org/u/pdeffebach)\
**Post date:** [October 12, 2020, 6:16pm UTC](https://discourse.julialang.org/t/regression-fixed-effects-quantiles/45427/15 "2020-10-12T18:16:42Z")

</div>

I’m not sure what you mean. Non-`Real` variables are automatically treated as fixed effects in `glm`.

```julia
julia> df = DataFrame(income = rand(1000), edu = rand(1000), state = string.(rand(1:50,1000)));

julia> lm(@formula(income ~ edu + state), df)
StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Array{Float64,1}},GLM.DensePredChol{Float64,LinearAlgebra.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}

income ~ 1 + edu + state

Coefficients:
──────────────────────────────────────────────────────────────────────────────
                    Coef. Std. Error t Pr(>|t|) Lower 95% Upper 95%
──────────────────────────────────────────────────────────────────────────────
(Intercept) 0.507432 0.0681312 7.45 <1e-12 0.373726 0.641137
edu 0.0327742 0.0324312 1.01 0.3125 -0.0308709 0.0964193

```

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [October 12, 2020, 7:13pm UTC](https://discourse.julialang.org/t/regression-fixed-effects-quantiles/45427/16 "2020-10-12T19:13:07Z")

</div>

@pdeffebach  
Imagine I have a numeric code for `state` (like a FIPS code).  
Compare: [GLM](https://github.com/JuliaStats/GLM.jl).jl vs [QuantileRegressions](https://github.com/pkofod/QuantileRegressions.jl).jl vs [FixedEffectModels](https://github.com/FixedEffects/FixedEffectModels.jl).jl

```julia
using DataFrames, GLM, QuantileRegressions, FixedEffectModels;
N=1000;
df=DataFrame(income = rand(N), edu = rand(N), state = rand(1:50,N));
f1=@formula(income ~ edu + state)
f2=@formula(income ~ edu + fe(state))
#
a1 = GLM.lm(f1, df)
b1 = qreg(f1, df, .5)
c1 = reg(df, f1)
#
a2 = GLM.lm(f2, df) #Error
b2 = qreg(f2, df, .5) #Error
c2 = reg(df, f2)

```

Formula `f1` treats `state` as a continuous variable & works w/ all packages that accept such formulae.

Formula `f2` treats `state` as a fixed effect, and only works with [FixedEffectModels](https://github.com/FixedEffects/FixedEffectModels.jl).jl but not with [GLM](https://github.com/JuliaStats/GLM.jl).jl & [QuantileRegressions](https://github.com/pkofod/QuantileRegressions.jl).jl.  
It would be nice if `fe()` was in `StatsModels.jl` (or somewhere) allowing fe to automatically work w/ all formulas.

Similarly, @matthieu’s [Vcov](https://github.com/matthieugomez/Vcov.jl).jl currently only works w/ [FixedEffectModels](https://github.com/FixedEffects/FixedEffectModels.jl).jl but not yet with [GLM](https://github.com/JuliaStats/GLM.jl).jl etc.

```julia
using Vcov
GLM.lm(f1, df, Vcov.cluster(:state)) #Error
reg(df, f1, Vcov.cluster(:state)) #no Error ;)

```

I hope we can increase modularity in the Julia stats ecosystem…

Just to be clear, in the trivial example above formula `f2` can be easily done by hand.  
My note is more generally about putting certain features that many different users want like fe into a more general packages available to everybody.  
In STATA: `reg y F(-1 0 1 2).D i.x1#(c.x2 i.x3)` works across methods

---

<div class="post-metadata">

**Author:** ![pdeffebach](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pdeffebach/32/10320_2.png) [@pdeffebach](https://discourse.julialang.org/u/pdeffebach)\
**Post date:** [October 12, 2020, 7:23pm UTC](https://discourse.julialang.org/t/regression-fixed-effects-quantiles/45427/17 "2020-10-12T19:23:00Z")

</div>

I absolutely think we need a better way of using fixed effects in GLM and statsmodels mode generally. I think either `fe(x)` or `categorical(x)` or something should be used to refer to fixed effects. For reference, the way to do this in StatsModels currently is

```julia
julia> lm(@formula(income ~ edu + state), df, contrasts = Dict(:state => EffectsCoding()))

```

Also, you have a typo in your post

```julia
reg(df, f1, Vcov.cluster(:state)) #Error

```

This works fine, does not error.

I definitely think that

1. StatsModels should support `fe`
2. FixedEffectModels should support `glm` and the general ecosystem better.