title: Julia: Multinomial Regression with time series lagged values
body: <p>I need to do some multinomial regression in Julia. In R I get the following result:</p>
  
  <pre><code>library(nnet)
  data &lt;- read.table("Dropbox/scripts/timeseries.txt",header=TRUE)
  multinom(y~X1+X2,data)
  
  # weights:  12 (6 variable)
  initial  value 10985.024274 
  iter  10 value 10438.503738
  final  value 10438.503529 
  converged
  Call:
  multinom(formula = y ~ X1 + X2, data = data)
  
  Coefficients:
    (Intercept)        X1        X2
  2   0.4877087 0.2588725 0.2762119
  3   0.4421524 0.5305649 0.3895339
  
  Residual Deviance: 20877.01 
  AIC: 20889.01 
  </code></pre>
  
  <p>Here is my <a href="https://www.dropbox.com/s/1h7uns55t34c8fr/timeseries.txt?dl=0" rel="nofollow">data</a></p>
  
  <p>My first attempt was using <a href="https://github.com/lindahua/Regression.jl" rel="nofollow">Regression.jl</a>. The documentation is quite sparse for this package so I am not sure which category is used as baseline, which parameters the resulting output corresponds to, etc. I filed an issue to ask about these things <a href="https://github.com/lindahua/Regression.jl" rel="nofollow">here</a>. </p>
  
  <pre><code>using DataFrames
  using Regression
  import Regression: solve, Options, predict
  dat = readtable("timeseries.txt", separator='\t')
  X = convert(Matrix{Float64},dat[:,2:3])
  y = convert(Vector{Int64},dat[:,1])
  
  ret = solve(mlogisticreg(X',y,3), reg=ZeroReg(),  options=Options(verbosity=:iter))
  </code></pre>
  
  <p>the result is </p>
  
  <pre><code>julia&gt; ret.sol
  3x2 Array{Float64,2}:
   -0.573027  -0.531819
    0.173453   0.232029
    0.399575   0.29979 
  </code></pre>
  
  <p>but again, I am not sure what this corresponds to.</p>
  
  <p>Next I tried the Julia wrapper to Python's SciKitLearn:</p>
  
  <pre><code>using ScikitLearn
  @sk_import linear_model: LogisticRegression
  
  model = ScikitLearn.fit!(LogisticRegression(multi_class="multinomial", solver = "lbfgs"), X, y)
  model[:coef_]
  3x2 Array{Float64,2}:
   -0.261902    -0.220771 
   -0.00453731   0.0540354
    0.266439     0.166735 
  </code></pre>
  
  <p><s>but I have not figured out how to extract the coefficients from this model.</s> Updated with coefficients. These also don't look like the R results. </p>
  
  <p>Any help trying to replicate R's results would be appreciate (using whatever package!) </p>
  
  <p>Note the response variables are just the discretized time-lagged response i.e. </p>
  
  <pre><code>julia&gt; dat[1:3,:]
  3x3 DataFrames.DataFrame
  | Row | y | X1 | X2 |
  |-----|---|----|----|
  | 1   | 3 | 1  | 0  |
  | 2   | 3 | 0  | 1  |
  | 3   | 1 | 0  | 1  |
  </code></pre>
  
  <p>for row 2 you can see that the response (0, 1) means the previous observation was a 3. Similarly (1,0) means previous observation was a 2 and (0,0) means previous observation was a 1. </p>
  
  <p>Update:
  For Regression.jl it seems it does not fit an intercept by default (and they call it "bias" instead of an intercept). By adding this term we get results very similar to python (not sure what the third column is though..)</p>
  
  <pre><code>julia&gt; ret = solve(mlogisticreg(X',y,3, bias=1.0), reg=ZeroReg(),  options=Options(verbosity=:iter))
  julia&gt; ret.sol
  3x3 Array{Float64,2}:
   -0.263149    -0.221923   -0.309949
   -0.00427033   0.0543008   0.177753
    0.267419     0.167622    0.132196
  </code></pre>
  
  <p>UPDATE:
  Since the model coefficients are not identifiable I should not be expecting them to be the same acrossed these different implementations. However, the predicted probabilities should be the same, and in fact they are (using R, Regression.jl, or ScikitLearn). </p>