It’s been a year since my last response and unfortunately, I need to repeat: the whole problem is that the LogitModelFit code requests the fit of a model that is overparameterized. This is not an opinion in that I think the model is too complex and needs simplification to make it more easily understood. Stated another way: The given code results in a singular design matrix and the blame resides in the request to attempt to fit an overparameterized model and not with Mathematica. (Although Mathematica could certainly put out a more direct error message.)
The model being requested to fit with LogitModelFit[data, {1, x, a}, {x, a}, NominalVariables -> {a}]
is
$$\log\left({{p_i}\over{1-p_i}}\right)=\mu+a_i \alpha_1+(1-a_i)\alpha_0 +\beta x_i$$
where
$p_i$ is the probability a success,
$\mu$ is the intercept (requested by the 1 in {1, x, a}),
$\alpha_0$ and
$\alpha_1$ are the effects of the nominal variable a (requested by the a in {1, x, a} and NominalVariables -> {a}), and
$\beta$ is the coefficient multiplied by
$x_i$ (requested by the x in {1, x, a}).
This model is overparameterized and LogitModelFit doesn’t handle that very well.
Why is it overparameterized? Let’s rewrite the equation as
$$\log\left({{p_i}\over{1-p_i}}\right)=\mu+\alpha_0+a_i(\alpha_1-\alpha_0)\beta x_i$$
We see that there are really on 3 parameters
$\mu+\alpha_0$,
$\alpha_1-\alpha_0$, and
$\beta$ that can be estimated. You might think that request LogitModelFit[data, {1, x, a}, {x, a}, NominalVariables -> {a}] should be only fitting 3 parameters but that construction tells LogitModelFit to attempt to fit 4 parameters and it chokes.
Not convinced? Consider the design matrix given the data.
(* Get the columns of the design matrix ready *)
n = Length[data];
ones = ConstantArray[1, n]; (* For the intercept mu *)
covariate = data[[All, 1]]; (* x *)
a = {#[[2]], 1 - #[[2]]} & /@ data;
(* Design matrix *)
m = Transpose[{ones, covariate, a[[All, 1]], a[[All, 2]]}];
MatrixRank[m]
(* 3 *)
The rank is 3 (and not 4). The rank of 3 indicates that only 3 unique parameters can be fit (or put another way "uniquely identified"). So it is not logical to expect nice results from an overparameterized model. As long as only 3 parameters are requested to be estimated everything is fine. And there are several ways to obtain identical model predictions.
So how to get the model that one thinks they should be getting? There are several ways. Here are 4 equivalent ways to obtain the identical predictions.
(* Equivalent models *)
m1 = Transpose[{ones, data[[All, 1]], data[[All, 2]]}];
lm1 = LogitModelFit[{m1, data[[All, 3]]}]
m2 = Transpose[{data[[All, 1]], 1 - data[[All, 2]], data[[All, 2]]}];
lm2 = LogitModelFit[{m2, data[[All, 3]]}]
lm3 = LogitModelFit[data, {x, a}, {x, a}, NominalVariables -> a,
IncludeConstantBasis -> False]
lm4 = LogitModelFit[data, {x, a}, {x, a}]
(* Equivalent predictions *)
lm1["PredictedResponse"]
(* {0.971958, 0.779499, 0.297112, 0.514681, 0.95155, 0.951966, 0.527294,
0.753359, 0.737949, 0.55053, 0.955198, 0.825052, 0.875761, 0.535885,
0.187753, 0.998126, 0.98149, 0.696078, 0.997495, 0.911263} *)
lm2["PredictedResponse"]
(* {0.971958, 0.779499, 0.297112, 0.514681, 0.95155, 0.951966, 0.527294,
0.753359, 0.737949, 0.55053, 0.955198, 0.825052, 0.875761, 0.535885,
0.187753, 0.998126, 0.98149, 0.696078, 0.997495, 0.911263} *)
lm3["PredictedResponse"]
(* {0.971958, 0.779499, 0.297112, 0.514681, 0.95155, 0.951966, 0.527294,
0.753359, 0.737949, 0.55053, 0.955198, 0.825052, 0.875761, 0.535885,
0.187753, 0.998126, 0.98149, 0.696078, 0.997495, 0.911263} *)
lm4["PredictedResponse"]
(* {0.971958, 0.779499, 0.297112, 0.514681, 0.95155, 0.951966, 0.527294,
0.753359, 0.737949, 0.55053, 0.955198, 0.825052, 0.875761, 0.535885,
0.187753, 0.998126, 0.98149, 0.696078, 0.997495, 0.911263} *)
Until better error messages are added and the documentation deals with more detail about NominalVariables one should construct their own nominal/dummy variables.
