Previous: Simulating Parameters | Up: Compatibility with sim() | Next: Getting Ready for the |

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.

`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

- 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