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

Calculating Quantities of Interest

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:

The required arguments for the qi() function are:

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:
  1. Extract any systematic parameters by substituting the name of your systematic parameter (defined in describe.mymodel()).
  2. Extract any ancillary parameters (defined in describe.mymodel()) by substituting their names here.
  3. 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].
  4. Replace rnorm() with a function that takes random draws from the stochastic component of your model.



Gary King 2010-09-09