All models require a model-specific method for calculating quantities
of interest from the simulated parameters. For a model of class contrib, the appropriate qi() function is qi.contrib(). This function should calculate, at the bare minimum,
the following quantities of interest:
- ev: the expected values, calculated from the analytic
solution for the expected value as a function of the systematic
component and ancillary parameters.
- pr: the predicted values, drawn from a distribution
defined by the predicted values. If R does not have a built-in
random generator for your function, you may take a random draw from
the uniform distribution and use the inverse CDF method to calculate
predicted values.
- fd: first differences in the expected value, calculated
by subtracting the expected values given the specified x from
the expected values given x1.
- ate.ev: the average treatment effect calculated using the
expected values ev. This is simply y - ev, averaged across
simulations for each observation.
- ate.pr: the average treatment effect calculated using
the predicted values pr. This is simply y - pr,
averaged across simulations for each observation.
The required arguments for the qi() function are:
- object: the zelig output object.
- par: the simulated parameters.
- x: the matrix of explanatory variables (created using setx()).
- x1: the optional matrix of alternative values for first
differences (also created using setx()). If first differences
are inappropriate for your model, you should put in a warning() or stop() if x1 is not NULL.
- y: the optional vector or matrix of dependent variables
(for calculating average treatment effects). If average treatment
effects are inappropriate for your model, you should put in a warning() or stop() if conditional prediction has been
selected in the setx() step.
Continuing the Normal regression example from above, the appropriate qi.normal() function is as follows:
qi.normal <- function(object, par, x, x1 = NULL, y = NULL) {
Beta <- parse.par(par, eqn = "mu") # [1]
sigma2 <- parse.par(par, eqn = "sigma2") # [2]
ev <- Beta %*% t(x) # [3a]
pr <- matrix(NA, ncol = ncol(ev), nrow = nrow(ev))
for (i in 1:ncol(ev))
pr[,i] <- rnorm(length(ev[,i]), mean = ev[,i], # [4]
sigma = sd(sigma2[i]))
qi <- list(ev = ev, pr = pr)
qi.name <- list(ev = "Expected Values: E(Y|X)",
pr = "Predicted Values: Y|X")
if (!is.null(x1)){
ev1 <- par %*% t(x1) # [3b]
qi$fd <- ev1 - ev
qi.name$fd <- "First Differences in Expected Values: E(Y|X1)-E(Y|X)"
}
if (!is.null(y)) {
yvar <- matrix(rep(y, nrow(par)), nrow = nrow(par), byrow = TRUE)
tmp.ev <- yvar - qi$ev
tmp.pr <- yvar - qi$pr
qi$ate.ev <- matrix(apply(tmp.ev, 1, mean), nrow = nrow(par))
qi$ate.pr <- matrix(apply(tmp.pr, 1, mean), nrow = nrow(par))
qi.name$ate.ev <- "Average Treatment Effect: Y - EV"
qi.name$ate.pr <- "Average Treatment Effect: Y - PR"
}
list(qi=qi, qi.name=qi.name)
}
There are five lines of code commented above. By changing these five
lines in the following four ways, you can write qi()
function appropriate to almost any model:
- Extract any systematic parameters by substituting the name of
your systematic parameter (defined in describe.mymodel()).
- Extract any ancillary parameters (defined in describe.mymodel()) by substituting their names
here.
- Calculate the expected value using the inverse link function and
141#141
. (For the normal model, this is linear.) You will
need to make this change in two places, at Comment [3a] and [3b].
- Replace rnorm() with a function that takes random draws
from the stochastic component of your model.
Gary King
2010-09-09