# Probit with Equality Constraint

I have a set of data that I analyze using a probit model. However, because of the way that the data are produced, the effect of two variables should be equal. Say, it is a room tilting down from the southheast corner and the x- and y- coordinates of the points in the room should have the same result along the 45 degree northeast sloping lines. How can I make the probit do its analysis under the constraint that the coefficient of the x and y coordinates must be equal?

For example:

 data = {
{0, .9, .1},
{1, .05, .85},
{0, .7, .7},
{1, .3, .3},
{0, .7, .2},
{1, .15, .72}};
pr = ProbitModelFit[data, {x, y}, {x, y}];
pr["ParameterTable"]
<out...>
... Estimate ... P-Value
...
x    -2.8 ...
y    -6.16 ...


The result of the above is different estimated coefficients for x and y, -2.8 and -6.16. I want those two coefficients to be equal. But I want the variables to be different so that they can have different p-values.

This is an extended comment.

I see that there are two issues:

1. There is but a single P-value associated with each parameter. One can't get two P-values when assuming two predictor variables have the same (but unknown) parameter value. Once you assume equality of the parameters, you have one parameter.

2. To clarify what you want you need to specify explicitly the model being fit. Just giving the code to analyze the data isn't sufficient (at least in this case).

The model being fit is $$probit = a_0 + a_x*x + a_y*y$$. If you want to force $$a_x = a_y = a_{xy}$$, then you're fitting the model $$probit = a_0 + a_{xy}*(x+y)$$. You'll need to add a column of the sum of $$x$$ and $$y$$ and you'll get an estimate of one parameter ($$a_{xy}$$) and a single P-value for that parameter.

If what you want is a test of $$a_x = a_y$$, then you can do that with the model you originally fit. (Although I certainly hope that the data provided is just a toy example and not your real data as there's just not enough data to perform a probit analysis.)

(* Get parameter estimates and estimated covariance matrix *)
estimates = pr["BestFitParameters"]
(* {4.20492, -2.80857, -6.1636} *)
covmat = pr["CovarianceMatrix"]
(* {{21.5565, -16.4668, -28.2627}, {-16.4668, 13.2453, 21.3889}, {-28.2627, 21.3889, 38.1694}} *)

(* Test of equality of ax and ay *)
zValue = (estimates[[2]] - estimates[[3]])/Sqrt[covmat[[2, 2]] + covmat[[3, 3]] - 2 covmat[[2, 3]]]
(* 1.14161 *)
pValue = 2 (1 - CDF[NormalDistribution[0, 1], Abs[zValue]])
(* 0.253617 *)

• Thanks Jim. Your first alternative is what I was asking and thanks for the explanation. I need to better specify my question. As for the data, indeed, this is merely an example to illustrate the issue. The real dataset is over 2k observations. Nov 24 '21 at 19:51
• I voted to close this question and I will come back when I have specified it better. Nov 24 '21 at 20:10