# Simulate probit marginal effects with an interaction term

I have the following probit model:

library(foreign)
library(MASS)

data <- data[data$year > 1959,] model1 <- glm(allons3 ~ confbord + rpc + confbord_rpc + neighlgdp + polity2l + polity2sq + lgdp96l + lnpop + postcoldw + peaceall, data = data, family = "binomial"(link = "probit"))  The data along with a stata do file can be downloaded at: Journal of Peace Research: Bratithwaite The paper can be found at: Braithwaite2010 In the paper from which the data is taken, the author simulates the effects of the interaction term in the model (the confbord_rpc term, the multiplication has been done in advance). "confbord" is a dummy variable and "rpc" ranges from 0 to 7. He first calculates the effects with "confbord" set to 0, and "rpc" varying from 0 to 7 (everything else held constant at their means), and then subtracts the effects from the same models with "confbord" held constant at 1. He obtains a very substantial effect from the "cpr" variable through this procedure. I have been replicating this using the following code:  #simulate coefficients based on the model beta <- coef(model1) covvar <- summary(model1)$cov.unscaled
coefs.sim <- mvrnorm(10000, beta, covvar)

#let rpc vary hold confbord constant at 0, and the other variables at their means
constant <- 1
confbord <- 0
rpc <- seq(0,7,1)
ncivwar_rpc <- ncivwar * rpc
neighlgdp <- 8.077
polity2l <- -0.4691
polity2sq <- 56.51
lgdp96l <- 8.123
lnpop <- 9.028
postcoldw  <- 0
peaceall <- 16.14

#bind them together in a matrix
effects0 <- cbind(constant,ncivwar,rpc, ncivwar_rpc, neighlgdp, polity2l, polity2sq, lgdp96l, lnpop, postcoldw, peaceall)

#multiply the effects matrix with the matrix of simulated coefficients to obtain marginal effects
results0 <- coefs.sim %*% t(effects0)

#evaluate each marginal effect on the standardized normal CDF, to obtain predicted probabilities
results0 <- apply(results0,1:2, function(x) pnorm(x))

#do the same procedure for confbord held constant at 1
constant <- 1
confbord <- 1
rcp <- seq(0,7,1)
ncivwar_rpc <- ncivwar * rpc
neighlgdp <- 8.077
polity2l <- -0.4691
polity2sq <- 56.51
lgdp96l <- 8.123
lnpop <- 9.028
postcoldw  <- 0
peaceall <- 16.14

effects1 <- cbind(constant,ncivwar,rpc, ncivwar_rpc, neighlgdp, polity2l, polity2sq, lgdp96l, lnpop, postcoldw, peaceall)
results1 <- coefs.sim %*% t(effects1)
results1 <- apply(results1,1:2,function(x) pnorm(x))

# subtract the estimated probabilities from each other
results <- results1 - results0

#graph the results
x <- seq(0,7,1)
means <- apply(results,2,mean)
sdUpper <- apply(results,2,function(x) quantile(x, .975))
sdLower <- apply(results,2,function(x) quantile(x, .025))
plot(x,means, ylim = c(-1,.2), type = "l")
lines(x,sdUpper, lty = 2)
lines(x,sdLower, lty = 2)
lines(x,rep(0,8), col = "red")


The problem is that with this procedure I get a much lower effect than the author reports in the paper (figure 1), can anyone spot what I am doing wrong?

You ran the regression for Fig. 2 that uses confbord variable, but you try to get the margins from regression for Fig. 1 that uses nciwvar variable.
• The Stata code looks rather inefficient, even accounting for it being four versions old. These days, all this machinery is bundled in margins and marginsplot, and they produce results that are quite different from what's published in the paper. Approaching this by simulation sounded a bit odd to me. This looks more like a job for Stata's nlcom. Ah well... it's not my paper that we are looking at here :) Commented Aug 25, 2011 at 18:05