Probit Regression in R, Python, Stata, and SAS

Group by hc

We start by tabulating lfp and hc.

addmargins(table(mroz$lfp, mroz$hc, deparse.level=2))

##         mroz$hc
## mroz$lfp   0   1 Sum
##      0   207 118 325
##      1   251 177 428
##      Sum 458 295 753

Then, we calculate the adjusted predictions of lfp for the two levels of hc variable.

hc_data0 = data.frame(k5 = mean(mroz$k5), k618 = mean(mroz$k618), age=mean(mroz$age), 
                      lwg= mean(mroz$lwg), inc=mean(mroz$inc), hc=0, wc=mean(mroz$wc))
hc_data1 = data.frame(k5 = mean(mroz$k5), k618 = mean(mroz$k618), age=mean(mroz$age), 
                      lwg= mean(mroz$lwg), inc=mean(mroz$inc), hc=1, wc=mean(mroz$wc))
h0 = predict(mroz.probit, hc_data0, type="response", se=TRUE)
h1 = predict(mroz.probit, hc_data1, type="response", se=TRUE)
hc_fit = data.frame(Margin = c(h0$fit[1], h1$fit[1]), se=c(h0$se.fit[1], h1$se.fit[1]))
hc_fit

##      Margin         se
## 1 0.5693807 0.02726260
## 2 0.5917191 0.03463592

The marginal probability of husband being a college graduate is is 0.59, while the marginal probability of husband being a high school graduate is slightly lower at 0.57.

Group by Wife’s College Attendance

Here is the tabulation of lfp and wc. As can be seen below, there is a significant imbalance between labor force participation, when the women are college graduates.

addmargins(table(mroz$lfp, mroz$wc, deparse.level=2))

##         mroz$wc
## mroz$lfp   0   1 Sum
##      0   257  68 325
##      1   284 144 428
##      Sum 541 212 753

We also calculate the adjusted predictions of lfp for the two levels of wc variable, while keeping other variables at mean.

wc_data0 = data.frame(k5 = mean(mroz$k5), k618 = mean(mroz$k618), age=mean(mroz$age), 
                      lwg= mean(mroz$lwg), inc=mean(mroz$inc), hc = mean(mroz$hc), wc=0)
wc_data1 = data.frame(k5 = mean(mroz$k5), k618 = mean(mroz$k618), age=mean(mroz$age), 
                      lwg= mean(mroz$lwg), inc=mean(mroz$inc), hc = mean(mroz$hc), wc=1)
w0 = predict(mroz.probit, wc_data0, type="response", se=TRUE)
w1 = predict(mroz.probit, wc_data1, type="response", se=TRUE)
wc_fit = data.frame(Margin = c(w0$fit[1], w1$fit[1]), se=c(w0$se.fit[1], w1$se.fit[1]))
wc_fit

##      Margin         se
## 1 0.5238094 0.02402095
## 2 0.7081631 0.03847631

The marginal probability increases from 0.524 to 0.708 when women are college graduates vs. only high-school graduates.

Group by Age and Wife’s College Attendance

We can predict the outcome variable by both age and women’s college education. We will use increments of 10 between 30 and 60 as representative ages. As the results show, the older a woman gets, the less are her chances of participating in the labor force. Regardless of age, women’s college education always yields a higher probability of labor force participation.

wc_data0_age=data.frame(k5=rep(mean(mroz$k5), 4), k618=rep(mean(mroz$k618), 4), 
                        age=c(30, 40, 50, 60), lwg=rep(mean(mroz$lwg), 4), 
                        inc=rep(mean(mroz$inc), 4), wc=rep(0, 4), 
                        hc=rep(mean(mroz$hc), 4))
wc_data1_age=data.frame(k5=rep(mean(mroz$k5), 4), k618=rep(mean(mroz$k618), 4), 
                        age=c(30, 40, 50, 60), lwg=rep(mean(mroz$lwg), 4), 
                        inc=rep(mean(mroz$inc), 4), wc=rep(1, 4), 
                        hc=rep(mean(mroz$hc), 4))
m0=predict(mroz.probit, wc_data0_age, type="response", se=TRUE)
m1=predict(mroz.probit, wc_data1_age, type="response", se=TRUE)
wc_fit_age = data.frame(Age=c(30, 40, 50, 60),Margin_wc0=m0$fit, Margin_wc1=m1$fit, se_wc0=m0$se.fit, se_wc1=m1$se.fit)
wc_fit_age

##   Age Margin_wc0 Margin_wc1     se_wc0     se_wc1
## 1  30  0.7033092  0.8466691 0.03936089 0.03561642
## 2  40  0.5618680  0.7402177 0.02509639 0.03716213
## 3  50  0.4119514  0.6047963 0.03193009 0.04744612
## 4  60  0.2739989  0.4552319 0.04823869 0.06726396

Using ggplot, we can plot the marginal effects to visualize the output. As can be seen from the plot, the variability of predictions increase with the age. Nevertheless, it is clear that a woman’s labor force participation is strongly affected by her age.

ggplot(data=wc_fit_age,aes(y=Age)) +
geom_line(aes( x= Margin_wc0), colour="blue") +
geom_errorbarh(aes(xmin=Margin_wc0-2*se_wc0, xmax=Margin_wc0+2*se_wc0), height=.1, colour="blue")+
geom_line(aes(x = Margin_wc1), colour="red")+
geom_errorbarh(aes(xmin=Margin_wc1-2*se_wc1, xmax=Margin_wc1+2*se_wc1),
               height=.1,colour="red")+scale_color_manual(labels = c("wc=0","wc=1"), values = c("blue", "red"))+
coord_flip()+xlab('Probability(Lfp)') +
  ylab("Age")

Group by Number of Children 5 years or younger

Below is a cross-tab of lfp with k5. The tabulation shows that the more the number of kids aged below 6 years, the less women participate in the labor force.

addmargins(table(mroz$lfp, mroz$k5, deparse.level=2))

##         mroz$k5
## mroz$lfp   0   1   2   3 Sum
##      0   231  72  19   3 325
##      1   375  46   7   0 428
##      Sum 606 118  26   3 753

We can predict lfp at different levels of k5, by keeping other variables at mean. When k5=0, the marginal probability of a woman’s participation in the labor force is 0.657, whereas when they have three or more, the probability is at around 0.013. This range is a good indicator of the significance of k5 variable.

k_data=data.frame(k5=c(0,1,2,3), 
                  k618=rep(mean(mroz$k618), 4), age=rep(mean(mroz$age), 4), 
                  lwg=rep(mean(mroz$lwg), 4), inc=rep(mean(mroz$inc), 4), 
                  wc = rep(mean(mroz$wc), 4), hc=rep(mean(mroz$hc), 4))
k0=predict(mroz.probit, k_data, type="response", se=TRUE)
k_fit = data.frame(k5=c(0,1,2,3),Margin_k=k0$fit, se_k=k0$se.fit)
k_fit

##   k5   Margin_k       se_k
## 1  0 0.65730850 0.02052209
## 2  1 0.31932620 0.03563207
## 3  2 0.08942632 0.03349808
## 4  3 0.01324307 0.01087202

We can also plot the marginal effects of k5. As the plot below shows, there is a significant drop in the probability of a woman being employed with the increase in the number of children below the age of 6.

ggplot(data=k_fit,aes(y=k5)) +
geom_line(aes( x= Margin_k), colour="blue") +
geom_errorbarh(aes(xmin=Margin_k-2*se_k, xmax=Margin_k+2*se_k), height=.1, colour="blue")+
coord_flip()+xlab('Probability(Lfp)') +
  ylab("k5")

Group by Husband’s College Attendance

First we group the data by hc.

To better understand hc, we tabulate lfp and hc. This is done with the pandas crosstab function.

print(pd.crosstab(data["lfp"], data["hc"], margins = True))

hc     False  True  All
lfp                    
False    207   118  325
True     251   177  428
All      458   295  753

Then we calculate adjusted predictions of lfp for two levels of the hc variable. Creating these adjusted predictors is not very clean in Python compared to other languages. We have to build the array from scratch.

hc_data0 = np.column_stack((
    1,
    np.mean(data["k5"]),
    np.mean(data["k618"]),
    np.mean(data["age"]),
    np.mean(data["wc"]),
    0,
    np.mean(data["lwg"]),
    np.mean(data["inc"])
    ))
    
    hc_data1 = np.column_stack((
    1,
    np.mean(data["k5"]),
    np.mean(data["k618"]),
    np.mean(data["age"]),
    np.mean(data["wc"]),
    1,
    np.mean(data["lwg"]),
    np.mean(data["inc"])
    ))

The results of running this prediction is as follows.

print(probit_model.predict(hc_data0))
print(probit_model.predict(hc_data1))

[ 0.56938181]
[ 0.59171968]

We see that the marginal probability of husband being a collage graduate is 0.59, while the marginal probability of husband being a high school graduate is lower at 0.57.

Group by Wife’s College Attendance

Similarly to hc, we tabulate our response variable lfp and predictor wc. We can see that there is a significant imbalance between labor force participation for women that are college graduates.

print(pd.crosstab(data["lfp"], data["wc"], margins = True))

wc     False  True  All
lfp                    
False    257    68  325
True     284   144  428
All      541   212  753

We also calculate the adjusted predictions of lfp for the two levels of wc variable, while keeping other variables at the mean.

wc_data0 = np.column_stack((
    1,
    np.mean(data["k5"]),
    np.mean(data["k618"]),
    np.mean(data["age"]),
    0,
    np.mean(data["hc"]),
    np.mean(data["lwg"]),
    np.mean(data["inc"])
    ))
wc_data1 = np.column_stack((
    1,
    np.mean(data["k5"]),
    np.mean(data["k618"]),
    np.mean(data["age"]),
    1,
    np.mean(data["hc"]),
    np.mean(data["lwg"]),
    np.mean(data["inc"])
    ))

The result of running this prediction is as follows:

print(probit_model.predict(wc_data0))
print(probit_model.predict(wc_data1))

[ 0.52380974]
[ 0.70816505]

It appears that the marginal probability increases from 0.524 to 0.708 when women are college graduates compared to highschool graduates.

Group by Age and Wife’s College Attendance

We can predict the outcome variable by both age and women’s college education. We use increments of 10 between 30 and 60 as our representative ages.

wc_data0 = np.column_stack((
    np.repeat(1,4),
    np.repeat(np.mean(data["k5"]),4),
    np.repeat(np.mean(data["k618"]), 4),
    (30,40,50,60),
    np.repeat(0,4),
    np.repeat(np.mean(data["hc"]),4),
    np.repeat(np.mean(data["lwg"]),4),
    np.repeat(np.mean(data["inc"]),4)
    ))
wc_data1 = np.column_stack((
    np.repeat(1,4),
    np.repeat(np.mean(data["k5"]),4),
    np.repeat(np.mean(data["k618"]), 4),
    (30,40,50,60),
    np.repeat(1,4),
    np.repeat(np.mean(data["hc"]),4),
    np.repeat(np.mean(data["lwg"]),4),
    np.repeat(np.mean(data["inc"]),4)
))

The results of running the prediction is as follows:

print(probit_model.predict(wc_data0))
print(probit_model.predict(wc_data1))

[ 0.70330952  0.5618684   0.41195181  0.27399923]
[ 0.84667045  0.74021954  0.6047985   0.45523423]

As we can see, the older a women gets, the less chance of her participating in the labor force. Regardless of age, women’s college education always yields a higher probability of labor force participation.

Group by Number of Children 5 years or younger (K5)

Finally we look at the effects of number of children 5 years or younger. As before, here is a cross tabulation of this variable.

print(pd.crosstab(data["lfp"], data["k5"], margins = True))

k5       0    1   2  3  All
lfp                        
False  231   72  19  3  325
True   375   46   7  0  428
All    606  118  26  3  753

Then we predict lfp at different levels of k5, while keeping other variables at their means.

k5_data = np.column_stack((
    np.repeat(1,4),
    (0,1,2,3),
    np.repeat(np.mean(data["k618"]), 4),
    np.repeat(np.mean(data["age"]),4),
    np.repeat(np.mean(data["wc"]),4),
    np.repeat(np.mean(data["hc"]),4),
    np.repeat(np.mean(data["lwg"]),4),
    np.repeat(np.mean(data["inc"]),4)
    ))
    
    print(probit_model.predict(k5_data))

[ 0.65730924  0.31932735  0.08942703  0.01324326]

We see that if there are three or more children under the age of 5, there is a 0.013 chance of a woman being in the work force, but when there is no children there is a 0.6573 chance. This is a good indication of the significance of this variable.

Extension

Finally, an overall marginal effect can be observed by calling the get_margeff() method of the probit model class. This is unique to Python. Results are as follows.

mfx = probit_model.get_margeff()
print(mfx.summary())

 Probit Marginal Effects   
 
=====================================
Dep. Variable:                    lfp
Method:                          dydx
At:                           overall
==============================================================================
                dy/dx    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
k5            -0.2997      0.034     -8.726      0.000      -0.367      -0.232
k618          -0.0132      0.014     -0.955      0.340      -0.040       0.014
age           -0.0130      0.002     -5.219      0.000      -0.018      -0.008
wc             0.1673      0.045      3.696      0.000       0.079       0.256
hc             0.0196      0.042      0.461      0.645      -0.064       0.103
lwg            0.1253      0.029      4.311      0.000       0.068       0.182
inc           -0.0070      0.002     -4.451      0.000      -0.010      -0.004
==============================================================================

Group by hc

First, we make a table of frequency counts of hc and lfp and we predict the lfp for two groups: hc=0 and hc=1, while keeping other variables at mean.

tab lfp hc

*use margins for each level of hc
margins hc, atmeans

The marginal probability of hc=1 (husband has attained college) is 0.59 and it slightly higher than the marginal probability of hc=0 (husband has not attained college), which is 0.57. There is not an obivious differnce. It is reasonable because the p-value of hc is very high.

Group by wc

The table of frequency shows that when wc=0, the proportion of lfp is average, which is close to 0.5. However, when wc=1, the proportion of lfp=1 is much higher.

tab lfp wc

We predict the lfp for two groups: wc=0 and wc=1, and we keep other variables at mean.

*use margins for each level of wc
margins wc, atmeans

The result shows that the marginal probability is 0.71 when wc=1, and the marginal probability is 0.52 when wc=0. The probability of participating in the labor-force is higher when wife has attended college. We can say that wife’s college attendance is an important predictor.

We can go deeper on the predictor wc. We predict lfp for group by age and wc. Age is at every 10 years of age from 30 to 60. Since the output of marginal function is long, we make a plot to visualize the output and it is easier to interpert.

*use margins for each level of wc and age
margins, at(age=(30(10)60) wc=(0 1)) atmeans vsquish

marginsplot

Looking at the table, we can see the predicted probability of lfp is 0.7033095 if age is 30 and wc=0.It will increase to 0.8466704 if age is 30 and wc=1. From the marginal plot, we can conclude that when age is increasing, the probability is decreasing. Also, The probability of wc=1 is always higher than wc=0. At age 60, the variablity is the highest because the 95% confidence interval is the widest.

Group by k5

The table of frequency shows that the proportion of lfp is decreasing when k5 is increasing.

tab lfp k5

We predict the lfp by k5= 0 1 2 3, and we keep other variables at mean. Also, we make a plot to visualize the data.

*use margins for each level of k5
margins, at(k5=(0 1 2 3)) atmeans

marginsplot

The output shows that when women do not have any children 5 years old or younger, the probability of participating in the labor-force is 0.66 which is higher than the average. However, after they had children, the probability of participating in the labor-force is decreasing. Therefore, we can conclude that k5 is a significant predictor.

Display the summary table of the dataset

(Seperate the 3 binary variables and the others).

/*For the other non-binary variables*/
proc means data=mroz;
 var k5 k618 age lwg inc;
run;

/*For those 3 binary variables, consider individuals and interactions between reponse and wc/hc.*/
proc freq data=mroz;
 tables lfp wc hc lfp*wc lfp*hc;
run;

Marginal effects plot with 95%CI, for wc, at means:

/*Marginal effects plot with 95%CI, for wc, at means*/
* k5=0.238 k618=1.353 age=42.54 lwg=1.097 inc=20.13 hc=0.392*;
*In this plot, actually we only need the 2 points on the curve: where wc=0 and wc=1*;
proc logistic data=mroz descending plots=EFFECT;
  class lfp / param=ref ;
  model lfp = k5 k618 age lwg inc wc hc /link=probit;
  output out=estimated predicted=estprob l=lower95 u=upper95;
run;

In this plot, we only need the 2 points on the curve: where wc=0 and wc=1.

Marginal effects plot with 95%CI, for hc, at means:

/*Marginal effects plot with 95%CI, for hc, at means*/
* k5=0.238 k618=1.353 age=42.56 lwg=1.097 inc=20.13 wc=0.282*;
*In this plot, actually we only need the 2 points on the curve: where hc=0 and hc=1*;
proc logistic data=mroz descending plots=EFFECT;
  class lfp / param=ref ;
  model lfp = k5 k618 age lwg inc wc hc /link=probit;
  output out=estimated predicted=estprob l=lower95 u=upper95;
run;

In this plot, we only need the 2 points on the curve: where hc=0 and hc=1.

Marginal effects plot for wc*hc, at means

/*Marginal effects plot for wc*hc, at means*/
proc logistic data=mroz descending plots=EFFECT;
  class lfp wc hc / param=ref ;
  model lfp = k5 k618 age lwg inc wc hc /link=probit;
  output out=estimated predicted=estprob l=lower95 u=upper95;
run;

Diagnosis

The code above also creates several plots for model diagnosis:

We can see from these diagonosis plots, the probit regression model’s performance is not bad.