Software 17.5.1 (Logit by Newton–Raphson). The birthwt data in MASS record \(189\) births at a US hospital, with an indicator of low birth weight (under \(2.5\) kg) and the mother’s weight, smoking, hypertension, uterine irritability and race. Newton–Raphson iterates \(\beta \leftarrow \beta + H^{-1}s\) with the score and Hessian of Theorem 17.1.2:

> library(MASS)
> bw <- birthwt
> bw$race <- factor(bw$race, labels = c("white", "black", "other"))
> X <- model.matrix(~ lwt + smoke + ht + ui + race, bw); y <- bw$low
> b <- rep(0, ncol(X))                          # Newton-Raphson for the logit
> for (it in 1:25) {
+   p <- drop(plogis(X %*% b))
+   H <- crossprod(X * (p * (1 - p)), X)        # minus the Hessian
+   step <- solve(H, crossprod(X, y - p)); b <- b + drop(step)
+   if (max(abs(step)) < 1e-10) break
+ }
> logit <- glm(low ~ lwt + smoke + ht + ui + race, binomial, bw)
> round(cbind(newton = b, glm = coef(logit), se = sqrt(diag(solve(H)))), 4)
             newton     glm     se
(Intercept)  0.0563  0.0563 0.9379
lwt         -0.0167 -0.0167 0.0068
smoke        1.0358  1.0358 0.3926
ht           1.8714  1.8714 0.6909
ui           0.9050  0.9050 0.4476
raceblack    1.3246  1.3246 0.5215
raceother    0.9262  0.9262 0.4304
> c(iterations = it)
iterations
         6

Reading the output. Newton–Raphson converges in six iterations to exactly the glm estimates, as global concavity guarantees, and the inverse of the final \(H\) gives the standard errors. Smoking, hypertension and uterine irritability raise the risk of low birth weight; a heavier mother lowers it.

Stata: logit low lwt smoke ht ui i.race. EViews: Quick \(\to \) Estimate Equation, method BINARY, with the logit option.

Software 17.5.2 (Probit, marginal effects and the linear probability model). The same model by probit and by least squares, with average marginal effects of smoking and the logit odds ratio:

> probit <- glm(low ~ lwt + smoke + ht + ui + race, binomial("probit"), bw)
> lpm <- lm(low ~ lwt + smoke + ht + ui + race, bw)
> ame <- function(fit, dens) mean(dens(predict(fit))) * coef(fit)
> round(rbind(logit = coef(logit), probit = coef(probit),
+             ratio = coef(logit) / coef(probit))[, 2:5], 3)
          lwt smoke    ht    ui
logit  -0.017 1.036 1.871 0.905
probit -0.010 0.633 1.122 0.550
ratio   1.693 1.636 1.668 1.646
> round(c(AME.logit = ame(logit, dlogis)[["smoke"]],
+         AME.probit = ame(probit, dnorm)[["smoke"]],
+         LPM = coef(lpm)[["smoke"]],
+         odds.ratio = exp(coef(logit)[["smoke"]])), 3)
 AME.logit AME.probit        LPM odds.ratio
     0.188      0.193      0.184      2.817
> c(McFadden.R2 = 1 - logLik(logit)[1] / logLik(update(logit, . ~ 1))[1],
+   LR = 2 * (logLik(logit)[1] - logLik(update(logit, . ~ 1))[1]),
+   LPM.outside = sum(fitted(lpm) < 0 | fitted(lpm) > 1))
McFadden.R2          LR LPM.outside
  0.1297785  30.4553763   9.0000000

Reading the output. The logit coefficients are \(1.64\) to \(1.69\) times the probit ones, as the scale normalisation predicts. The three models agree on the effect that matters: smoking raises the probability of a low-weight birth by about \(19\) percentage points (logit \(0.188\), probit \(0.193\), LPM \(0.184\)). In odds terms, smoking multiplies the odds by \(2.8\). The model explains a modest part of the variation (McFadden \(R^2 = 0.13\)), but the likelihood ratio test of all the slopes, \(30.5\) on \(6\) degrees of freedom, is highly significant. Nine of the \(189\) LPM fitted probabilities lie outside \([0, 1]\).

Stata: probit ..., then margins, dydx(smoke); logit ..., or for odds ratios. EViews: the probit option of BINARY, then View \(\to \) Expectation-Prediction Evaluation.

Software 17.5.3 (Tobin’s durable goods data). The tobin data in the survival package are Tobin’s own (1958) sample of \(20\) households: expenditure on durable goods, age of the head, and a liquidity ratio. The Tobit likelihood of Theorem 17.2.1(iii) is maximised directly, and checked against survreg, which fits it as a left-censored normal regression:

> library(survival)
> d <- tobin; c(n = nrow(d), zeros = sum(d$durable == 0))
    n zeros
   20    13
> nll <- function(p) {                          # Tobit log-likelihood, minus
+   xb <- p[1] + p[2] * d$age + p[3] * d$quant; s <- exp(p[4])
+   -sum(ifelse(d$durable > 0, dnorm(d$durable, xb, s, log = TRUE),
+               pnorm(-xb / s, log.p = TRUE)))
+ }
> ols <- coef(lm(durable ~ age + quant, d))
> o <- optim(c(ols, log(sd(d$durable))), nll, method = "BFGS",
+            control = list(reltol = 1e-14, maxit = 1000,
+                           parscale = c(10, 0.1, 0.01, 1)))
> sr <- survreg(Surv(durable, durable > 0, type = "left") ~ age + quant,
+               data = d, dist = "gaussian")
> round(rbind(tobit = c(o$par[1:3], sigma = exp(o$par[4])),
+             survreg = c(coef(sr), sr$scale),
+             OLS.all = c(ols, NA),
+             OLS.pos = c(coef(lm(durable ~ age + quant, d,
+                                 subset = durable > 0)), NA)), 4)
        (Intercept)     age   quant  sigma
tobit       15.1448 -0.1291 -0.0455 5.5725
survreg     15.1449 -0.1291 -0.0455 5.5725
OLS.all     11.0743 -0.0261 -0.0346     NA
OLS.pos     13.7781  0.3392 -0.1045     NA
> xb <- drop(cbind(1, d$age, d$quant) %*% o$par[1:3]); s <- exp(o$par[4])
> round(c(P.positive = mean(pnorm(xb / s)),
+         dE.dquant = mean(pnorm(xb / s)) * o$par[[3]]), 4)
P.positive  dE.dquant
    0.3604    -0.0164

Reading the output. Thirteen of the twenty households bought nothing. The direct maximisation and survreg agree. Least squares on all the data attenuates both coefficients, the age coefficient by a factor of five; least squares on the seven buyers even reverses the sign of the age effect. The average probability of buying is \(0.36\), so by the McDonald–Moffitt decomposition a unit increase in the liquidity ratio lowers expected expenditure by \(0.36 \times 0.0455 = 0.016\), much less than the latent coefficient. With twenty observations none of the coefficients is precisely estimated, but the example shows clearly why the zeros cannot be ignored.

Stata: tobit durable age quant, ll(0). EViews: method CENSORED, with left censoring at zero.

Software 17.5.4 (Heckman’s two-step estimator). A simulated labour market with \(2000\) people: participation depends on \(z\) (for instance, non-labour income) and \(x\), wage offers depend on \(x\) only, and the two disturbances have correlation \(0.6\):

> set.seed(10); n <- 2000
> z <- rnorm(n); x <- rnorm(n)                  # z: not in the wage equation
> e <- MASS::mvrnorm(n, c(0, 0), matrix(c(1, 0.6, 0.6, 1), 2))
> work <- 0.5 + 1 * z + 0.5 * x + e[, 1] > 0    # participation
> w <- 1 + 0.5 * x + e[, 2]                     # wage offer, seen if work
> sel <- glm(work ~ z + x, binomial("probit"))
> zg <- predict(sel); mills <- dnorm(zg) / pnorm(zg)
> heck <- lm(w ~ x + mills, subset = work)
> round(rbind(OLS.workers = coef(lm(w ~ x, subset = work)),
+             Heckman = coef(heck)[1:2], truth = c(1, 0.5)), 3)
            (Intercept)     x
OLS.workers       1.270 0.405
Heckman           0.977 0.487
truth             1.000 0.500
> round(c(rho.sigma = coef(heck)[["mills"]], share.working = mean(work)), 3)
    rho.sigma share.working
        0.734         0.642

Reading the output. About \(64\%\) work. Least squares on the workers overstates the intercept and understates the slope: workers with low \(x\) are in work only if their wage disturbance is high, which flattens the observed relation. The two-step estimator recovers \(0.977\) and \(0.487\), close to the true \(1\) and \(0.5\). The coefficient of the inverse Mills ratio, \(0.73\), estimates \(\rho \sigma _u = 0.6\); it is significantly positive, which is evidence of selection.

Stata: heckman w x, select(work = z x) twostep. EViews: method HECKIT, choosing the two-step option.

Software 17.5.5 (Poisson and negative binomial regression). The quine data in MASS give the days absent from school in a year for \(146\) children in New South Wales, by ethnic background (Aboriginal or not), sex, age group and learner status. The robust standard errors are those of Theorem 17.4.1:

> q <- quine; c(mean = mean(q$Days), var = var(q$Days))
    mean      var
 16.4589 264.1673
> pois <- glm(Days ~ Eth + Sex + Age + Lrn, poisson, q)
> Xq <- model.matrix(pois); mu <- fitted(pois)
> B <- solve(crossprod(Xq * mu, Xq))           # inverse information
> robust <- sqrt(diag(B %*% crossprod(Xq * (q$Days - mu)) %*% B))
> nb <- glm.nb(Days ~ Eth + Sex + Age + Lrn, q)
> round(rbind(poisson = coef(pois), se.poisson = sqrt(diag(vcov(pois))),
+             se.robust = robust, negbin = coef(nb),
+             se.negbin = sqrt(diag(vcov(nb))))[, 1:3], 3)
           (Intercept)   EthN  SexM
poisson          2.715 -0.534 0.162
se.poisson       0.065  0.042 0.043
se.robust        0.235  0.153 0.155
negbin           2.895 -0.569 0.082
se.negbin        0.228  0.153 0.160
> c(dispersion = sum(residuals(pois, "pearson")^2) / df.residual(pois),
+   theta.nb = nb$theta)
dispersion   theta.nb
 13.166843   1.274893

Reading the output. The variance of days absent is \(16\) times the mean, so the data are heavily overdispersed. The Poisson coefficient of non-Aboriginal background, \(-0.534\), means \(100(1 - e^{-0.534}) = 41\%\) fewer days absent, other things equal. Its Poisson standard error, \(0.042\), is far too small: the robust standard error is \(0.153\), \(3.6\) times larger, which is close to \(\sqrt {13.2}\), the square root of the estimated dispersion. With the robust errors the sex difference is no longer significant. The negative binomial model gives similar coefficients with standard errors close to the robust ones, and an estimated \(\theta = 1.27\), implying substantial extra-Poisson variation.

Stata: poisson days i.eth i.sex i.age i.lrn, vce(robust) and nbreg. EViews: method COUNT, with the Poisson or negative binomial option and Huber/White standard errors.