title: Difference between glm and LogitModelFit
body: <p>I have a problem with <strong>glm</strong> function in R.</p>
  
  <p>Specifically, I am not sure how to include <strong>nominal variables</strong>.</p>
  
  <p>The results that I get in R after running the glm function are the following:</p>
  
  <pre><code>&gt; df
  
     x1 x2 y
  1  a  2  0
  2  b  4  1
  3  a  4  0
  4  b  2  1
  5  a  4  1
  6  b  2  0
  
  &gt; str(df)
  'data.frame':   6 obs. of  3 variables:
   $ x1: Factor w/ 2 levels "a","b": 1 2 1 2 1 2
   $ x2: num  2 4 4 2 4 2
   $ y: Factor w/ 2 levels "0","1": 1 2 1 2 2 1
  
  Call:
  glm(formula = y ~ x1 + x2, family = "binomial", data = df)
  
  Coefficients:
               Estimate Std. Error z value Pr(&gt;|z|)
  (Intercept)   -39.132  15208.471  -0.003    0.998
  x1b            19.566   7604.236   0.003    0.998
  x2              9.783   3802.118   0.003    0.998
  </code></pre>
  
  <p>However, when I run the <strong>LogitModelFit</strong> function in <strong>Wolfram Mathematica</strong> I get different parameters.</p>
  
  <p>The code in Wolfram is provided below:</p>
  
  <pre><code>data = {{a, 2, 0}, {b, 4, 1}, {a, 4, 0}, {b, 2, 1}, {a, 4, 1}, {b, 2, 0}};
  
  model = LogitModelFit[data, {x, y}, {x, y}, NominalVariables -&gt; x]
  
  model["BestFitParameters"]
  </code></pre>
  
  <p>And these are my estimated parameters:</p>
  
  <pre><code>{-18.5661, -18.5661, 9.28303}
  
  model // Normal
  
  1/(1 + E^(18.5661 - 9.28303 y + 18.5661 DiscreteIndicator[x, a, {a, b}]))
  </code></pre>
  
  <p>So, what is different here? Why the results differ so much?</p>
  
  <p>Am I doing something wrong in R or in Wolfram?</p>