1
$\begingroup$

Is any one here familiar with an R package called Zelig?

I have a data frame like this:

IQ   AGE
80   50
100  18
90   25

etc.

What I need to do is build a model of IQ given AGE, I am running these commands: z.out <- zelig(IQ~AGE,data=df,model="ls") this runs the what-if given age 110, what would be the IQ x.out <- setx(z.out, AGE=110) This is a simulation model where given the age 110, after running 1 million runs of simulation, what would be the IQ with 95% confidence interval. s.out <- sim(z.out,x.out, num=1000000, level=95)

I have a hard time understanding from what pool of data the sim() function draws the numbers. I read though the docs, but they are written for Ph.D. students, if not more advanced readers. I have asked the Zelig creators this question multiple times but they are directing me to the docs which I read multiple times, with no luck. However, one of the person that works with Zelig sent me this email:

Suppose that you fit $$\text{IQ} = a + \text{Age} * b + e$$ Then you get a table of regression coefficients where a=50, b=2, and their standard errors are something like $\text{s.e.}(a)=\sqrt{10}$ and $\text{s.e.}(b)=1$. These are all hypothetical examples. In maximum likelihood estimation, this regression output is another way of saying that $a$ and $b$ are distributed bivariate normal with means $[50,2]$ and there's a variance-covariance matrix that looks something like this (all numbers are made up): $$\begin{array}{cc} 10 & cov(a,b) \\ cov(a,b) & 1 \\ \end{array}$$ So, the variance of $a$ is 10, the variance of $b$ is 1, and their covariance is $cov(a,b)$. It won't be shown in your regression table, but Zelig remembers it for you. Let's pretend it's 3. This variance-covariance matrix is the inverse of the Hessian I mentioned earlier. Don't worry about it. For this example, you need only remember that $\text{mean}(a,b) = [50,2]$ and $cov(a,b)=\begin{array}{cc}10&3\\3&1\end{array}$. For purposes of plain text email, I'm representing matrices with columns separated by commas and rows by semicolons. In addition, suppose that the error term $e$ is distributed with mean 0 and s.d.=1. Now, one way to predict what IQ you might get for somebody aged 88, based on this regression table, is exactly what you would expect: you simply calculate 50 + 88 * 2 = 226. This is your point estimate. The 95% confidence interval around this point estimate is a function of the standard errors of the coefficient estimates of $a=50$ and $b=2$, and the exact formula for that is in any econometrics textbook. Simulation makes it unnecessary to dig up that textbook. Instead, for 1000 rounds, sim() will come up with 1000 different pairs of $(a,b)$ estimates drawn from the bivariate normal with mean=[50,2] and cov=$\begin{array}{cc}10&3\\3&1\end{array}$. One such pair might be $(47,1.5)$; another might be $(52,3)$; yet another might be $(10,5)$. Whatever they are, sim() plugs them into the formula and gives you 1000 different estimates for the IQ. Their average is your point estimate. If you stack them from lowest to highest, the ends of the 95% confidence interval are the top 25th value and the bottom 25th. That's it. That's all that sim() does.

Given the above explanation, can anybody tell me in lay terms, what numbers sim() is picking? How are those numbers in pool generated? I would greatly appreciate if anyone brings some light into this.

$\endgroup$
  • $\begingroup$ @user777, I dont know what you mean $\endgroup$ – user1471980 Dec 5 '13 at 16:12
4
$\begingroup$

The function sim() randomly selects numbers are randomly chosen from a bivariate normal distribution with the specified vector of means and variance-covariance matrix in order to construct confidence intervals around the parameter estimates.

The vector of means and the covariance matrix are determined by the ML estimates of the model parameters. Therefore, these are simulations of parameter values. Recall that the ML estimates of parameters are just the most likely values. The true value might not actually be the maximum of the likelihood function given the data because of sampling error -- this is the motivation for finding confidence sets: to define a region with high probability of containing the true value.

To understand this topic more fully, you should familiarize yourself with probability density functions generally and the normal distribution and then the bivariate normal specifically. The latter extends the normal from the real line to a real plane, and so is useful in cases when two random variables can vary together.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ your answer is fair. Sim() pulls the numbers from bivariate normal distribution formed from the original values of IQ and AGE? I may be asking too much, but I wished someone run this with a very small data frame to illustrate the whole point. $\endgroup$ – user1471980 Dec 5 '13 at 17:04
  • $\begingroup$ Sim() pulls these values not the distribution of the data, but the distribution of the parameters. Recall that the ML estimates of parameters are just the most likely values. The true value might not actually be the maximum of the likelihood function given the data because of sampling error -- this is the motivation for finding confidence sets: to define a region with high probability of containing the true value. $\endgroup$ – Sycorax Dec 5 '13 at 17:07
  • $\begingroup$ would you be able to run the whole think with a df of 3 rows and 2 columns. That would really help me understand this. $\endgroup$ – user1471980 Dec 5 '13 at 17:17
  • $\begingroup$ You've already written the code to do this in your question; what part is unclear? I'm not sure what you're asking. $\endgroup$ – Sycorax Dec 5 '13 at 17:22
  • $\begingroup$ how do you form a bivariate normal? I think this is the key. Given bivariate normal, you pull numbers from this distrbution with the same mean that is equal to the original coefficients. $\endgroup$ – user1471980 Dec 5 '13 at 18:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.