Previous: Algorithms for setx

Up: Formulae - A Peek

Next: References

Algorithms for simqi

simqi simulates quantities of interest based on the parameters that were generated by estsimp and the $x$-values that were chosen with setx. The program obtains simulations of the dependent variable and uses them to calculate expected values, probabilities, first differences, and other quantities of interest. This procedure works in all cases but involves some approximation error, which users can make arbitrarily small by choosing a sufficient number of simulations. In many cases, though, shortcuts exist that can curtail both computation time and approximation error. simqi employs such shortcuts whenever possible. Here, we sketch the algorithms for each model that ${\mathfrak{C}}$larify supports.

regress: The exact algorithm in simqi depends on whether the user has transformed the dependent variable (e.g., taken the log of $y$) prior to estimation. If no such transformation has occurred, the program generates one predicted value according to the formula $\tilde{y} = X_c\tilde{\beta} + \tilde{\epsilon}$, where $\tilde\beta$ is a vector of simulated effect coefficients and $\tilde{\epsilon}$ is one draw from $N(0,\tilde{\sigma}^2)$. Likewise, the program simulates one expected value as $\tilde{E}(y) = X_c\tilde{\beta}$. The algorithm becomes a bit more complicated if the user transformed the dependent variable prior to estimation, and would like to reverse the transformation when interpreting the results. Let $f$ represent a function, as identified by the tfunc() option, that reverses the transformation. If $f$ has been specified, the program simulates one predicted value according to the formula $f(X_c\tilde{\beta} + \tilde{\epsilon})$. For an expected value, the program draws $m$ values of $\tilde{\epsilon}_d$ $(d=1,2, \ldots ,m)$ from $N(0,\tilde{\sigma}^2)$ and then computes $(1/m)\sum_{d=1}^{m} f(X_c\tilde{\beta}+\tilde{\epsilon}_d)$, which is the average of $m$ predicted values.

logit: The formula for $\tilde{\pi}$, the simulated probability that the dependent variable $y$ takes on a value of 1, is $1/(1+e^{-X_c\tilde{\beta}})$. To obtain one simulation of $y$, the program draws a number from the Bernoulli distribution with parameter $\tilde{\pi}$.

probit: The formula for $\tilde{\pi}$, the simulated probability that the dependent variable $y$ takes on a value of 1, is $\Phi(X_c\tilde{\beta})$ where $\Phi$ is the c.d.f. of the standard normal distribution. To obtain one simulation of $y$, the program draws a number from the Bernoulli distribution with parameter $\tilde{\pi}$.

ologit: The exact formula depends on the number of categories in the dependent variable. Suppose there are three categories. Let $\tilde\beta$ represent one simulated vector of effect coefficients and let $\tilde\tau_{lo}$ and $\tilde\tau_{hi}$ stand for draws of the cutpoints. To obtain one simulation of the probabilities for each category $(y=0, y=1, y=2)$, the program calculates: $\tilde{\pi}_0 \equiv \tilde{Pr}(y=0) = \frac{1}{ 1 + e^ {\left( X_c\tilde{\beta} -\tilde\tau_{lo} \right) }}$, $\tilde{\pi}_1 \equiv \tilde{Pr}(y=1) = \frac{1} {1 + e^ {\left( X_c\tilde{\bet... ...ght) }} - \frac{1}{1 + e^ {\left( X_c\tilde{\beta} - \tilde\tau_{lo} \right) }}$, and $\tilde{\pi}_2 \equiv \tilde{Pr}(y=2) = 1 - \frac{1}{1 + e^ {\left( X_c\tilde{\beta} - \tilde\tau_{hi} \right)} }$. With these results, the program can draw a predicted value, $\tilde{y}$, from a multinomial distribution with parameters $\tilde{\pi}_0$, $\tilde{\pi}_1$, $\tilde{\pi}_2$, and $n=1$.

oprobit: The exact formula depends on the number of categories in the dependent variable. Suppose there are three categories. Let $\tilde\beta$ represent one simulated vector of effect coefficients and let $\tilde\tau_{lo}$ and $\tilde\tau_{hi}$ stand for draws of the cutpoints. To obtain one simulation of the probabilities for each category $(y=0, y=1, y=2)$, the program calculates $\tilde{\pi}_0 \equiv \tilde{Pr}(y=0) = \Phi(\tilde\tau_{lo} - X_c\tilde{\beta})$, $\tilde{\pi}_1 \equiv \tilde{Pr}(y=1) = \Phi(\tilde\tau_{hi} - X_c\tilde{\beta}) - \Phi(\tilde\tau_{lo} - X_c\tilde{\beta})$, and $\tilde{\pi}_2 \equiv \tilde{Pr}(y=2) = \Phi(X_c\tilde{\beta} - \tilde\tau_{hi})$. With these results, the program can draw a predicted value, $\tilde y$, from a multinomial distribution with parameters $\tilde{\pi}_0$, $\tilde{\pi}_1$, $\tilde{\pi}_2$, and $n=1$.

mlogit: The probability equation for the $K$ nominal outcomes of the multinomial logit is $\tilde{\pi}_j \equiv \tilde{Pr}(y=j) = \frac{e^{X_c\tilde{\beta}_j}}{\sum_{k=1}^{K}e^{X_c\tilde{\beta}_k}}$, where one of the $J$ outcomes is the base category, such that the effect coefficients $\tilde\beta$ for that category are set to zero. With these results, the program can draw a predicted value, $\tilde y$, from a multinomial distribution with parameters equal to the $\tilde\pi$s and $n=1$.

poisson: The formula for the expected value $\tilde\mu$ is $e^{X_c\tilde\beta}$, and the probability that the dependent variable takes on the integer value $j$ is $\tilde{Pr}(y=j) = \frac{e^{-\tilde\mu}\tilde\mu^j}{j!}$. To obtain one predicted value, the program draws $\tilde{y}$ from a Poisson distribution with parameter $\tilde\mu$. The Poisson simulator is adapted from Press, et al. (1992), pp. 293-95.

nbreg: The formula for the expected value $\tilde\mu$ is $e^{X_c\tilde\beta}$, just as in the Poisson regression model (Long 1997, pp. 230-33). The probability that the dependent value takes on the integer value $j$ can be simulated as $\tilde{Pr}(y=j)=\frac{\Gamma(j+\tilde\alpha^{-1})}{j!\Gamma(\tilde\alpha^{-1})... ...\alpha^{-1}} \left( \frac{\tilde\mu}{\tilde\alpha^{-1} + \tilde\mu} \right) ^ j$. simqi obtains $\tilde\alpha$, the ``overdispersion'' parameter, by drawing simulations of $ln(\alpha)$ and the other parameters from the multivariate normal distribution and then calculating $e^{ln(\alpha)}$. To obtain a predicted value $\tilde{y}$, the program draws one number from a poisson distribution with mean $e^{X_c\tilde{\beta} + \tilde\epsilon}$, where $e^{\tilde\epsilon}$ is simulated from a gamma distribution with shape parameter $\tilde\alpha^{-1}$ and scale parameter $\tilde\alpha$. When $\tilde\alpha^{-1} < 1$, the gamma simulator is based on the algorithm developed by Ahrens and Dieter, as described in Ripley (1987, p. 88). For other values of $\tilde\alpha^{-1}$, the gamma simulator is based on the procedure by Best, as described in Devroye (1986, p. 410).

sureg: As with regress, the algorithm for interpreting the results of a sureg depends on whether the user transformed the dependent variable. If the user estimated the model without transforming the dependent variable, the program generates one predicted value for equation $k$ according to the formula $\tilde{y}_k = X_c\tilde{\beta}_k + \tilde{\epsilon}_k$, where $\tilde\beta_k$ is a vector of simulated effect coefficients for equation $k$ and $\tilde{\epsilon}_k$ is a simulated disturbance term for that equation. Disturbances for all equations are drawn simultaneously from a multivariate normal distribution with mean 0 and variance matrix $\tilde\Sigma$, as obtained from the inverse Wishart. Likewise, the program simulates one expected value for equation $k$ as $\tilde{E}(y_k) = X_c\tilde{\beta}_k$. If the user has transformed the dependent variable, let $f$ represent the function that reverses the transformation. The program simulates one predicted value for equation $k$ according to the formula $f(X_c\tilde{\beta}_k + \tilde{\epsilon}_k)$. For an expected value, the program draws $m$ sets of disturbance terms from $N(0,\tilde\Sigma)$ and indexes them as $\tilde{\epsilon}_{k,d}$, where $k$ marks the equation and $d=1,2, \ldots ,m$. Then, for each equation $k$ the program computes $(1/m)\sum_{d=1}^{m} f(X_c\tilde{\beta}_k+\tilde{\epsilon}_{k,d})$, which is the average of $m$ predicted values.

weibull: The algorithm depends on which metric, proportional hazard (PH) or accelerated failure-time (AFT) metric, was used at the estsimp stage. The expected value is defined as $\tilde\lambda^{1/\tilde{p}}\Gamma \left( 1 + 1/{\tilde{p}} \right)$. In the AFT metric, $\tilde\lambda = e^{-X_c \tilde\beta \tilde{p}}$; in the PH metric, $\tilde\lambda = e^{X_c \tilde\beta}$. The program obtains simulations of the ancillary shape parameter $p$ drawing by $ln(p)$ and the other parameters from a multivariate normal distribution and then calculating $e^{ln(p)}$. To obtain a predicted value, the program draws one number from the Weibull distribution with parameters $\tilde\lambda$ and $\tilde{p}$.