Creating Marginal Effect Plots for Linear Regression Models in R

Suppose we estimate the model \(E(y) = \beta_0 + \beta_xx + \beta_zz + \beta_{xz}xz\) and want to calculate the the marginal effect of \(x\) on \(E(y)\) as \(z\) varies \(\left(\dfrac{\partial E(y)}{\partial x}\right)\) and its 90% confidence interval. Brambor, Clark, and Golder describe exactly how to do this and provide Stata code. Below is a bit of R code to do something similar. (Click here to continue reading.)

compactr is now on CRAN

I've been working on a package called compactr  that helps create nice-looking plots in R and it is now up on CRAN.

You can get it by typing

directly into the command line in R. See how it works by typing  example(eplot)  or reading the details here. Below I describe the basic structure and functions. (Click here to continue reading.)

More on What Can Be Learned from Statistical Significance

A while back, I reacted to a post by Justin Esarey, in which he argues that not much can be learned from statistical significance. The basic question he's asking is as follows: For a fixed sample size, what does the posterior probability of the null hypothesis look like if we update based on the result of the hypothesis test rather than the entire data set? His answer is that statistical significance doesn't allow us to update much if our prior is appropriately skeptical.

He makes two points:

  1. A small magnitude but statistically significant result contains virtually no important information
  2. Even a large magnitude, statistically significant result is not especially convincing on its own.

In this post, I explore how we can set up the problem in reasonable way and have hypothesis test be informative. In particular, I'd like to argue that statistical significance can contain a lot of information about the hypothesis. (Click here to continue reading.)

The Problem with Testing for Heteroskedasticity in Probit Models

A friend recently asked whether I trusted the inferences from heteroskedastic probit models. I said no, because the heteroskedastic probit does not allow a researcher to distinguish between non-constant variance and a mis-specified mean function.

In particular, my friend had a hypothesis that the variance of the latent outcome (commonly called "y-star") should increase with an explanatory variable of interest. He was using the heteroskedastic probit model, which looks something like \(Pr(y_i = 1) = \Phi(X_i\beta, e^{Z_i\gamma})\), where \(\Phi()\) is the cumulative normal with mean \(X_i\beta\) and \(e^{Z_i\gamma}\).

He wanted to argue that his explanatory variable increased both the mean function (\(X\beta\)) and the variance function (\(e^{Z\gamma}\)). To do this, he included his variable in both the \(X\) and \(Z\) matrices and tested the statistical significance of the associated coefficients. He found that they were both significant. It would seem that his variance increases the mean and the variance of the latent outcome. He wanted to know if this was good evidence for his theory.

I replied that I did not think so, because a binary outcome variable doesn't contain any direct information about a non-constant variance. Indeed, the variance of a Bernoulli random variable is tied directly to the probability of success. This implies that any inference about changes in \(e^{Z\gamma}\) must come from observed changes in the probability of a success (i.e. changes in the mean function). Because we've assumed a specific (i.e. linear) functional form for the mean function, deviations from this will be attributed to the variance function. Because of this structure, the results are driven totally by our assumption of the linearity of the mean function. Indeed, it would not be hard to find a plausible non-linear mean function (e.g. quadratic specification) that makes the \(\gamma\) parameter no longer significant.

Example One

I thought a good way to illustrate this claim would be to show that for a large but plausible sample size of one million, the heteroskedastic probit will suggest a non-constant variance when the relationship is simply a logit.

To see an illustration of this, start by simulating data from a simple logit model. Then estimate a regular probit and a heteroskedastic probit.

n <- 10^6
x <- runif(n)
y <- rbinom(n, 1, plogis(-4 + 3*x))

r1 <- glm(y~x,family=binomial(link="probit"))

h1 <- hetglm(y ~ x)

Now if the coefficient for x is significant in the model of the scale, then we should conclude there is heteroskedasticity, right? No, because we already know that the latent variance is constant. However, we've barely mis-specified the link function (we're using a probit, the true model is logit). This slight mis-specification causes the results to point toward non-constant variance.

Coefficients (binomial model with probit link):
             Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.090026   0.007269  -287.5   <2e-16 ***
x            1.589261   0.007418   214.2   <2e-16 ***

Latent scale model coefficients (with log link):
  Estimate Std. Error z value Pr(>|z|)
x -0.20780    0.01475  -14.09   <2e-16 ***

Here is a plot of the predicted probabilities from the true, probit, and heteroskedastic probit models. Notice that in the range of the data, the heteroskedastic probit does a great job of representing the relationship. However, that's not because the variance is non-constant as the heteroskedastic probit would suggest. It's because the link function is slightly mis-specified.


Example Two

I think the logit example makes the point powerfully, but let's look at a second example just for kicks. This time, let's say that we believe there's heteroskedasticity that can be accounted for by x, so we estimate a heteroskedastic probit and include x in the mean and variance function. However, again we're wrong. Actually the true model has a constant variance, but a non-linear mean function (\(\beta_0 + \beta_1x^2\)).

If we simulate the data and estimate the model, we see again that our mis-specified mean function leads us to conclude that the variance is non-constant. In this case though, the mis-specification is severe enough that you'll find significant results with much smaller sample. If we made conclusions about the non-constant variance from the statistical significance of coefficients in the model of the variance, then we would be led astray.

Coefficients (binomial model with probit link):
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  -3.6122     0.3726  -9.694   <2e-16 ***
x             3.1493     0.2809  11.210   <2e-16 ***

Latent scale model coefficients (with log link):
  Estimate Std. Error z value Pr(>|z|)
x  -0.8194     0.2055  -3.988 6.66e-05 ***

Again, the plot shows that the heteroskedastic probit does a good job at adjusting for the mis-specified mean function (working much like a non-parametric model).


So What Should You Do?

I think that researchers who have a theory that allows them to speculate about the mean and variance of a latent variable should go ahead and estimate a statistical model that maps cleanly onto their theory (like the heteroskedastic probit). However, these researchers should realize that this model does not allow them to distinguish between non-constant variance and a mis-specified mean function.

Software Signals

This blog post by Sean Taylor generated quite a stir. He discussed the signals one sends by using certain software packages and seems to think that R users are more competent. The reactions ranged from amusement to bashing.



While I don't think this type of post is particularly useful, it is fun (especially the John Myles White line), so I'm writing up my thoughts on the issue.

(Click here to continue reading.)

For more posts, see the Archives.