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

Thanks!

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

Thanks!

1 Like

We have one package for fixed effects:

A separate for quantile regression:

Currently no package that automates both together.

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

3 Likes

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

btw MLJLinearModels 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.

3 Likes

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

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?

1 Like

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

```
"""
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
```

2 Likes

I first modify the data by manually including FEs.

Then I write the formula corresponding to the regression with FE.

This is included here:

```
#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.

1 Like

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

2 Likes

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â€¦

3 Likes

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).

2 Likes

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

1 Like

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

1 Like