Dichotomous Logit and Probit

The purpose of this session is to show you how to use STATA's "canned" procedures for doing dichotomous Logit and Probit analysis. This includes obtaining predicted probabilities, predictions of the dependent variable, coefficients and marginal effects for the variables, model diagnostics, hypothesis tests, and the heteroskedastic Probit model. In addition, we provide programs that obtain some of these outputs the "hard" way for illustrative purposes only.

/* This program illustrates some of the "canned" procedures in STATA for doing Logit and Probit. In order to use all of this program you will

need to have SPOST and dprobit2/dlogit2 downloaded and installed. The data is a sample of 32 observations from a study by Spector and Mazzeo

on pre- and post-tests for economic literacy. The variables are:

GPA = grade point average

TUCE = test score on teaching college level economics

PSI = program participation variable (school district)

GRADE = response: 1 if grade on post-test is higher than on pre-test, 0 otherwise. */

/* First read the data. */

use "c:\users\Wood\Documents\My Teaching\Maximum Likelihood\Data\AldrichNelson.dta", clear

gen id = _n

keep if (1<=id<=32)

/* Binary choice models: Since the interest is in the PSI variable - this is the treatment effect- we will test the hypothesis

that the coefficient is equal to 0.0. These commands will again illustrate several ways to test hypotheses using STATA. First, the t-statistic

on PSI provides one test.*/

probit grade gpa tuce psi

/* We can also use a Wald test. For a single coefficient the chi-squared is going to be the same as the square of the t-ratio. */

test psi=0

/* For more involved problems, a likelihood ratio test might preferable. To carry out the likelihood ratio test, we require

chi - squared = 2 * (LOGL_with PSI - LOGL_without PSI)

We have the first one from the last estimate. We store the likelihood from the last estimation in A. To get the second,

we just leave PSI out of the model and run the test. */

est store A

probit grade gpa tuce

lrtest A

/* Now let's look at some of the available options on Logit/Probit procedures */

probit grade gpa tuce psi, robust /*Estimate the probit model with robust standard errors.*/

predict probs, p /*Calculate p(y=1) given the model for each y */

predict z, xb /*Calculate the index z for each y */

list grade probs z /*List y, predicted probabilities, and the index z */

/* We can gauge the fit of our model by putting together a prediction success table.*/

gen gradehat = probs >=.5

tab gradehat grade

/* A ROC (receiver operating characteristic) curve can also be used to evaluate the fit of models. In STATA

the best model is the one that comes closest to the left-top corner. In comparing models, the best model

would be the one with the largest area under the curve. It would also be the one closest to the left-top corner.

A completely uninformative model would be one that falls on the reference line along the diagonal of the

plot. */

lroc

/* We can also compare two models using the ROC approach. Note the last line which calls for a

hypothesis test for the first versus second model. */

probit grade gpa tuce psi, robust /*Estimate the first probit model with robust standard errors.*/

predict z1, xb /* Save the predictions */

lroc, nograph /* Estimate the ROC with no graph */

probit grade gpa tuce, robust /*Estimate the second probit model with robust standard errors.*/

predict z2, xb /* Save the predictions */

lroc, nograph /* Estimate the ROC with no graph */

roccomp grade z1 z2, graph summary /* Compare the two ROCs with a hypothesis test. */

/* Now let's use STATA's margins option to obtain some marginal effects. Note that the i. identifies a variable as a factor variable.*/

probit grade gpa tuce i.psi

/* First the marginal effect for a single variable */

margins psi

margins psi, atmeans

margins psi, at(gpa=3.117188)

margins psi, at(gpa=(2(.25)4))

/* Now plot the marginal effect of psi at various gpas */

margins psi, at (gpa=(2(.25)4))

marginsplot

/* Now average marginal effects as derivatives for all variables */

margins, dydx(*)

/* Now the marginal effect as a derivative for gpa conditional on psi */

margins psi, dydx(gpa)

/* Now let's use SPOST to obtain various fit and interpretation statistics. */

probit grade gpa tuce psi

fitstat /*Obtain various fit statistics on the probit */

listcoef, help /*List the coefficients and standardized coefficients*/

prvalue, x(psi=0) rest(mean) /*Compute p(grade=1) when psi=0, and rest of variables at means */

prvalue, x(psi=1) rest(mean) /*Compute p(grade=1) when psi=1, and rest of variables at means */

prchange, help /*Compute various first differences */

/*Now let's show how to graph the probit probabilities for psi at 0 and 1 */

gen ps1_0=normal(-7.45+1.62*gpa+.052*21.938)

gen ps1_1=normal(-7.45+1.62*gpa+.052*21.938+1.4263)

sort gpa

graph twoway (connected ps1_0 ps1_1 gpa if gpa>=2, s(..) c(l))

/*Now let's use dprobit to obtain various marginal effects from the probit model. */

dprobit grade gpa tuce psi /*The model in derivatives for continuous variables from means and with shift from 0-1 for dummies*/

dprobit2 grade gpa tuce psi /*The model in derivatives with all variables at their means*/

dprobit2 grade gpa tuce psi, at(psi==0) /*The model in derivatives when psi=0*/

dprobit2 grade gpa tuce psi, at(psi==1) /*The model in derivatives when psi=1*/

drop probs z

/*We can also use most of the preceding options with an estimated logit. */

logit grade gpa tuce psi /*Estimate the logit model */

predict probs, p /*Calculate p(y=1) given the model for each y */

predict z, xb /*Calculate the index z for each y */

list grade probs z /*List y, predicted probabilities, and the index z */

/*In addition, logit has an option to report odds ratios, rather than coefficients */

logistic grade gpa tuce psi

/* Now let's use STATA's margins option to obtain some marginal effects. Note that the i. identifies a variable as a factor variable.*/

logit grade gpa tuce i.psi

/* First the marginal effect for a single variable */

margins psi

margins psi, atmeans

margins psi, at(gpa=3.117188)

margins psi, at(gpa=(2(.25)4))

/* Now plot the marginal effect of psi at various gpas */

margins psi, at (gpa=(2(.25)4))

marginsplot

/* Now average marginal effects as derivatives for all variables */

margins, dydx(*)

/* Now the marginal effect as a derivative for gpa conditional on psi */

margins psi, dydx(gpa)

/*SPOST also works with logit */

logit grade gpa tuce psi

fitstat /*Obtain various fit statistics on the logit */

listcoef, help /*List the coefficients and standardized coefficients*/

prvalue, x(psi=0) rest(mean) /*Compute p(grade=1) when psi=0, and rest at mean */

prvalue, x(psi=1) rest(mean) /*Compute p(grade=1) when psi=1, and rest at mean */

prchange, help /*Compute various first differences */

/*Now let's use dlogit2 to obtain various marginal effects from the probit model. */

dlogit2 grade gpa tuce psi /*The model in derivatives with all variables at means*/

dlogit2 grade gpa tuce psi, at(psi==0) /*The model in derivatives when psi=0*/

dlogit2 grade gpa tuce psi, at(psi==1) /*The model in derivatives when psi=1*/

/*Now let's estimate the heteroskedastic probit, with the heteroskedasticity a function of the log of gpa */

gen loggpa=log(gpa)

hetprob grade gpa tuce psi, het(loggpa)

/*Now let's do a Wald test for heteroskedasticity from the previous model */

test loggpa

/*The preceding makes use of STATA's built-in functions. However, the curious among you may not be

satisfied with this, given your earlier experience with STATA's programming functions. Here is a program that

illustrates how to do Logit/Probit using your own log-likelihood functions.*/

clear

/* Let's create some simulated data to use in this demonstration. x and e are random normal numbers with

mean zero and standard deviation one. From x and e we create the index z. Then from z we create y the dependent

variable for logit and probit.*/

drawnorm x, n(200)

drawnorm e, n(200)

gen z=x+e

gen y=z>0

summarize

/*First, let's do logit on the simulated data */

program define mylogit

version 7

args lnf theta

quietly replace `lnf'=ln(1/(1+exp(-`theta'))) if $ML_y1==1

quietly replace `lnf'=ln(1/(1+exp(`theta'))) if $ML_y1==0

end

ml model lf mylogit (y=x)

ml maximize

program drop mylogit

/*Verify we've done it right using the canned procedure.*/

logit y x

/*Now, let's do probit on the simulated data */

program define myprobit

version 7

args lnf theta

quietly replace `lnf'=ln(normprob(`theta')) if $ML_y1==1

quietly replace `lnf'=ln(normprob(-`theta')) if $ML_y1==0

end

ml model lf myprobit (y=x)

ml maximize

program drop myprobit

/*Verify we've done it right using the canned procedure.*/

probit y x

/*Some of you may be curious concerning how STATA calculates the derivatives or marginal effects. Let's demonstrate how this can be done.

We return to the original dataset so that we can compare.*/

use "c:\users\Wood\My Documents\My Teaching\Maximum Likelihood\Data\AldrichNelson.dta", clear

gen id = _n

keep if (1<=id<=32)

gen one=1

mkmat gpa tuce psi one, matrix(x) /*create the x matrix */

mkmat one, matrix(i) /*create a vector of means of the x matrix*/

matrix ii=i'*i

matrix xbar=inv(ii)*x'*i

matrix list xbar

/*First let's do the probit derivatives.*/

probit grade gpa tuce psi /*Get the vector of estimated probit coefficients */

matrix betap=e(b)

matrix list betap

matrix zp=betap*xbar /*Retrieve the probit index z=xb calculated at the variable means */

matrix list zp

svmat zp, name(zp) /*Get the scale factor for probit */

scalar probitscale=normalden(zp)

display probitscale

matrix dprobit=probitscale*betap /*Calculate probit derivatives (or marginal effects) with all variables at the mean. */

matrix list dprobit

/*Now let's calculate the logit derivatives */

logit grade gpa tuce psi /*Get the vector of estimated logit coefficients */

matrix betal=e(b)

matrix list betal

matrix zl=betal*xbar /*Create the logit index z=xb at the means */

matrix list zl

svmat zl, name(zl) /*Get the scale factor for logit */

scalar logitscale=(1/(1+exp(-zl)))*(1-(1/(1+exp(-zl))))

display logitscale

matrix dlogit=logitscale*betal /*Calculate logit derivatives (or marginal effects) with all variables at the mean. */

matrix list dlogit

/* Some of you may also be interested in how to calculate the probit and logit probabilities. Let's calculate both

just for fun.*/

matrix xbp=x*betap' /* First the probit probabilities*/

svmat xbp, name(zip)

gen pp=normal(zip)

matrix xbl=x*betal' /* Now the logit probabilities */

svmat xbl, name(zil)

gen lp=(1/(1+exp(-zil)))

list grade zip zil pp lp /*List grade, the probit and logit indices, and the probit and logit probabilities side by side */

/*Now, let's do heteroskedastic probit using STATA's programming functions. We use starting values from

the regular probit above. */

program define heteroprob

version 7

args lnf theta1 theta2

quietly replace `lnf'=ln(normprob(`theta1'/exp(`theta2'))) if $ML_y1==1

quietly replace `lnf'=ln(normprob(-`theta1'/exp(`theta2'))) if $ML_y1==0

end

ml model lf heteroprob (eq1:grade=gpa tuce psi) (eq2:grade=psi,noconstant)

ml init 1.6 .05 1.4 -7.45 1, copy

ml maximize

program drop heteroprob

exit

/*That's the end of this lesson */