# Use of StatsModels?

**URL:** https://discourse.julialang.org/t/use-of-statsmodels/16979
**Category:** Statistics
**Created:** [October 30, 2018, 6:12pm UTC](https://discourse.julialang.org/t/use-of-statsmodels/16979 "2018-10-30T18:12:36Z")
**Posts on this page:** 8
**Page:** 1

<div class="post-metadata">

### Author: ![mcreel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcreel/32/30088_2.png) [@mcreel](https://discourse.julialang.org/u/mcreel)
#### Post date: [October 30, 2018, 6:12pm UTC](https://discourse.julialang.org/t/use-of-statsmodels/16979/1 "2018-10-30T18:12:36Z")

</div>

I’m interested in making some of my modeling functions accept dataframes and fomulae as inputs, using the StatsModels package. Is the following the proper way to achieve this?

My plan would be to do something like

```julia
using CSV, StatsModels, DataFrames
nerlove = CSV.read("nerlove.csv")
nerlove[:lnC] = log.(nerlove[:cost])
nerlove[:lnQ] = log.(nerlove[:output])
nerlove[:lnPL] = log.(nerlove[:labor])
nerlove[:lnPF] = log.(nerlove[:fuel])
nerlove[:lnPK] = log.(nerlove[:capital])
f = @formula(lnC ~ 1 + lnQ + lnPL + lnPF + lnPK)

```

At this point, I have a formula and a dataframe. Then, I would call a fitting function, e.g., a linear regression function that accepts a formula and a dataframe, like this:

```julia
ols(f,nerlove)

```

Inside ols(), I would get the dependent variables and regressor matrix doing something like

```julia
m = ModelFrame(f, nerlove)
mm = ModelMatrix(m)
y = Float64.(nerlove[f.lhs])
x = mm.m
b = x\y
etc.

```

I have tried this, and it works. My doubts are whether or not the last block is the best way to create ordinary arrays for computing regression coefficients, etc. I also don’t know how to deal with the possibility that the dataframe has missings in it, which is why I do the type conversion, assuming that there aren’t any. Pointers to examples would be very welcome.

---

<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 30, 2018, 6:53pm UTC](https://discourse.julialang.org/t/use-of-statsmodels/16979/2 "2018-10-30T18:53:10Z")

</div>

Is there a reason why you don’t want to use `GLM.jl`? If you have some custom estimators in mind, then you _might_ want to explore your own regression output type, but `GLM` is very flexible, and its easy to write a closure to store the outputs you want.

---

<div class="post-metadata">

### Author: ![mcreel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcreel/32/30088_2.png) [@mcreel](https://discourse.julialang.org/u/mcreel)
#### Post date: [October 30, 2018, 7:31pm UTC](https://discourse.julialang.org/t/use-of-statsmodels/16979/3 "2018-10-30T19:31:32Z")

</div>

Yes, the idea is to use it with more general estimators, OLS is just an example. I eventually would also like to use the table formatting, and other features in GLM, and related packages, instead of my home-cooked stuff.

---

<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: [October 30, 2018, 7:38pm UTC](https://discourse.julialang.org/t/use-of-statsmodels/16979/4 "2018-10-30T19:38:54Z")

</div>

For users, packages like GLM can be great, but for researchers and teachers there can be good reasons to avoid monolithic packages that end up working against you for your specific intents and purposes

---

<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 30, 2018, 8:23pm UTC](https://discourse.julialang.org/t/use-of-statsmodels/16979/5 "2018-10-30T20:23:11Z")

</div>

If you want to go all the way for full inter-operability with the rest of the ecosystem, then you can create your own model object `<: StatisticalModel` that implements the methods described [here](http://juliastats.github.io/StatsBase.jl/stable/statmodels.html).

---

<div class="post-metadata">

### Author: ![nalimilan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nalimilan/32/147_2.png) [@nalimilan](https://discourse.julialang.org/u/nalimilan)
#### Post date: [October 31, 2018, 7:10am UTC](https://discourse.julialang.org/t/use-of-statsmodels/16979/6 "2018-10-31T07:10:54Z")

</div>

IIRC observations with missing values are dropped automatically. Maybe @dave.f.kleinschmidt can confirm. You can find examples in GLM, Econometrics, Microeconometrics and FixedEffectsModels.

BTW, soon it should be possible to put the `log` calls directly in the formula (thanks to [this PR](https://github.com/JuliaStats/StatsModels.jl/pull/71)).

---

<div class="post-metadata">

### Author: ![mcreel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcreel/32/30088_2.png) [@mcreel](https://discourse.julialang.org/u/mcreel)
#### Post date: [October 31, 2018, 10:20am UTC](https://discourse.julialang.org/t/use-of-statsmodels/16979/7 "2018-10-31T10:20:46Z")

</div>

Thanks to all for the suggestions. I now have the basic stuff working, and I think I see where to look to find examples.

---

<div class="post-metadata">

### Author: ![dave.f.kleinschmidt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dave.f.kleinschmidt/32/55_2.png) [@dave.f.kleinschmidt](https://discourse.julialang.org/u/dave.f.kleinschmidt)
#### Post date: [October 31, 2018, 12:44pm UTC](https://discourse.julialang.org/t/use-of-statsmodels/16979/8 "2018-10-31T12:44:05Z")

</div>

Yup, any row with a missing in a column involved in the formula are dropped. That happens (currently) in the [`ModelFrame` constructor](https://github.com/JuliaStats/StatsModels.jl/blob/master/src/modelframe.jl#L145).

Note that you can get “automatic” support for formula/dataframe argument by specifying `RegressionModel` or `StatisticalModel` as the supertype of your model types. Then [this `fit` method will apply](https://github.com/JuliaStats/StatsModels.jl/blob/master/src/statsmodel.jl#L64) automatically. That’s the approach that GLM takes. However, the downside is that what is returned is actually a `DataFrameRegressionModel` wrapper, _not_ an instance of whatever you asked to be fit in the first place. For that reason, there’s discussion about removing this automatic fallback in [this issue](https://github.com/JuliaStats/StatsModels.jl/issues/32), which AFAIK everyone is in favor of but still not clear on how best to handle fitting models with formulae (or other kinds of data transformations).
