LogitModelFit function

Question

I have a problem with LogitModelFit function in Wolfram Mathematica.

Specifically, I am not sure how to include Nominal variables.

This is my approach:

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 -> x]

model["BestFitParameters"]

And these are my estimated parameters:

{-18.5661, -18.5661, 9.28303}

However, when I run the same regression in R programming language, I do not get the same parameters.

The results that I get in R after running the glm function are the following:

> 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

> 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(>|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

So, what is different here? Why the results differ so much?

Am I doing something wrong in Wolfram?

Answer 1

The short answer

Both Mathematica and R are asked to fit 3 parameters associated with the nominal variable: intercept, xA, and xB. But there's a linear dependency among those three parameters and only two can be estimated.

Mathematica chooses to set xB to zero. R chooses xA to be zero. Both approaches give the same predictions (other than potential machine precision and the difference in stopping rules associated with the iterative fitting procedure).

Mathematica (and SAS) choose to set the last value it sees as zero. R chooses the first value it sees as zero. (That is a slightly loose description and can be modified in all.)

The long answer

It is appreciated that you gave a complete example. But unlike regression examples with additive Gaussian errors, 6 data points is inadequate. Binary dependent variables need at least 20 to 100 data points for an instructive example. You'll notice that with your example the estimated standard errors are huge. That should be a red flag.

Here is a larger dataset:

data = {{"a", 0.05, 0.}, {"a", 0.1, 0.}, {"a", 0.15, 1.}, {"a", 0.2, 0.},
{"a", 0.25, 0.}, {"a", 0.3, 0.}, {"a", 0.35, 1.}, {"a", 0.4, 1.},
{"a", 0.45, 1.}, {"a", 0.5, 0.}, {"a", 0.55, 0.}, {"a", 0.6, 0.}, 
{"a", 0.65, 1.}, {"a", 0.7, 1.}, {"a", 0.75, 1.}, {"a", 0.8, 1.}, 
{"a", 0.85, 0.}, {"a", 0.9, 1.}, {"a", 0.95, 0.}, {"a", 1., 1.}, 
{"b", 0.0333333, 1.}, {"b", 0.0666667, 0.}, {"b", 0.1, 1.}, {"b", 0.133333, 1.},
{"b", 0.166667, 1.}, {"b", 0.2, 1.}, {"b", 0.233333, 1.}, {"b", 0.266667, 1.},
{"b", 0.3, 1.}, {"b", 0.333333, 1.}, {"b", 0.366667, 1.}, {"b", 0.4, 1.}, 
{"b", 0.433333, 1.}, {"b", 0.466667, 1.}, {"b", 0.5, 1.}, {"b", 0.533333, 1.},
{"b", 0.566667, 1.}, {"b", 0.6, 1.}, {"b", 0.633333, 1.}, {"b", 0.666667, 1.},
{"b", 0.7, 1.}, {"b", 0.733333, 1.}, {"b", 0.766667, 1.}, {"b", 0.8, 1.}, 
{"b", 0.833333, 1.}, {"b", 0.866667, 1.}, {"b", 0.9, 1.}, {"b", 0.933333, 1.},
{"b", 0.966667, 1.}, {"b", 1., 1.}};

The easiest way to get Mathematica and R to give the same estimates is to set the intercept to zero:

model = LogitModelFit[data, {x, y}, {x, y}, NominalVariables -> x, 
  IncludeConstantBasis -> False];
model["ParameterTable"]

In R one would use

summary(glm(formula = z ~ 0 + x + y, family = "binomial", data = df))

Coefficients:
   Estimate Std. Error z value Pr(>|z|)  
xa   -1.592      1.003  -1.587   0.1126  
xb    2.135      1.151   1.856   0.0635 .
y     3.032      1.671   1.815   0.0696 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1