7
$\begingroup$

I am comparing linear regression with and without intercept for the general sampling case. For this, I have $n$ samples of two correlated random variables $X \sim N\left(0,\sigma_X^2\right)$ and $Y \sim N\left(0, \sigma_Y^2\right)$ with correlation $\rho$.

For the random samples, I calculate the linear regression models with and without intercept
(1) $y_i=\alpha_0+\alpha_1x_i+\epsilon_1$ and
(2) $y_i=\beta_1x_i+\epsilon_2$

Using numerical experiments, I have found that $E[\hat\alpha_1] = E[\hat\beta_1]$, which seems logical to me. However, I have also found thtat $\text{Var}(\hat\alpha_1) \neq \text{Var}(\hat\beta_1)$, which I am currently trying to understand.

In another question of mine, I have found that for the general sampling case $\text{Var}(\hat \alpha_1) = \frac{\sigma_Y^2}{\sigma_X^2} \frac{1-\rho^2}{N-3}$ for the model with intercept and am trying to find $\text{Var}(\hat\beta_1)$.

Overall, I am therefore trying to find $\text{Var}(\hat\beta_1)=\text{Var}\left(\frac{\sum x_iy_i}{\sum x_i^2}\right)$.

The denominator is clearly gamma distributed. However, the distribution of the numerator as a sum of products of normal distributed random variables is tough, not to mention the ratio.

Calculating $\text{Var}(\hat\beta_1)=E[\hat\beta_1^2] - E[\hat\beta_1]^2$ isn't much easier, I think.

After spending hours in the local university library and searching research papers, I am turning to CrossValidated for help (again).

Does somebody know a way to calculate the variance in question?

$\endgroup$
  • 1
    $\begingroup$ Did you switch the definitions of $\alpha$ and $\beta$ partway through your post? $\endgroup$ – Dougal Aug 7 '14 at 3:57
  • $\begingroup$ Thank you, your comment is indeed correct. I have corrected the initial post. $\endgroup$ – sebastianb Oct 2 '14 at 13:49
  • 3
    $\begingroup$ Just as a comment: hope that you are aware that regression without intercept is generally a bad idea, see: stats.stackexchange.com/questions/7948/… $\endgroup$ – Tim Nov 25 '15 at 9:24
  • 4
    $\begingroup$ Is it purposeful that the true data generating process doesn't include a constant? In other words, that there's a zero in $y_i = 0 + \beta x_i + \epsilon_i$? If so, this question in some sense is, "What happens when I include a completely unrelated variable in a regression vs. when I don't include it?" (In this case, the completely unrelated variable would be the constant.) $\endgroup$ – Matthew Gunn Jun 12 '16 at 2:02
2
$\begingroup$

In your problem, assuming joint-normality of the variables, you can write the joint distribution of a single data point $(X_i, Y_i)$ as:

$$\begin{bmatrix} X_i \\ Y_i \end{bmatrix} \text{ ~ N} \left( \begin{bmatrix} 0 \\ 0 \end{bmatrix}, \begin{bmatrix} \sigma_X^2 & \rho \sigma_X \sigma_y \\ \rho \sigma_X \sigma_y & \sigma_Y^2 \end{bmatrix} \right).$$

With a bit of algebra, the value $Y_i$ can be shown to be equivalent to:

$$Y_i = \rho \frac{\sigma_Y}{\sigma_X} \cdot X_i + \sqrt{1 - \rho^2} \sigma_Y \cdot \varepsilon_i,$$

where $\varepsilon_i$ is an independent standard normal error term. Hence, the true regression model is:

$$\begin{matrix} Y_i = \beta_0 + \beta_1 \cdot X_i + \sigma \cdot \varepsilon_i & & & \beta_0 = 0 & \beta_1 = \rho \frac{\sigma_Y}{\sigma_X} & \sigma = \sqrt{1 - \rho^2} \sigma_Y \end{matrix}.$$


Simplifying the variance problem: When you estimate the model coefficients in the model with an intercept term, you would expect to get an intercept estimate close to the true value of zero, which means you would also expect the estimated slope coefficients to be similar with or without the inclusion of an intercept term in the model (as you have pointed out). However, the inclusion of an intercept term will tend to reduce the variance of the estimated slope coefficient. You have:

$$X_i Y_i = \rho \frac{\sigma_Y}{\sigma_X} \cdot X_i^2 + \sqrt{1 - \rho^2} \sigma_Y \cdot X_i \varepsilon_i.$$

Defining $\boldsymbol{Z} \equiv \boldsymbol{X} / \sigma_X$ and $\boldsymbol{U} \equiv \| \boldsymbol{Z} \|^2$ allows us to write:

$$\hat{\beta}_1 = \frac{\sum_{i=1}^n X_i Y_i}{\sum_{i=1}^n X_i^2} = \rho \frac{\sigma_Y}{\sigma_X} + \sqrt{1 - \rho^2} \frac{\sigma_Y}{\sigma_X} \cdot \frac{\boldsymbol{Z} \cdot \boldsymbol{\varepsilon}}{\boldsymbol{U}}.$$

Taking the variance gives:

$$\mathbb{V}{(\hat{\beta}_1)} = (1 - \rho^2) \frac{\sigma_Y^2}{\sigma_X^2} \mathbb{V}\left(\frac{\boldsymbol{Z} \cdot \boldsymbol{\varepsilon}}{\boldsymbol{U}} \right).$$

Now, the three vectors in the variance operator are independent (using Cochran's theorem). The numerator is a sum of products of independent standard normal random variables, and the denominator is a chi-squared random variable.

I must confess that I am not sure where to go from here. I do not recognise the variance expression as any simple form, though maybe others will. In any case, I hope that gives you some progress towards what you want.

$\endgroup$
0
$\begingroup$

I don't know how you found that $\sigma_{\hat{\alpha_1}}^2 \neq \sigma_{\hat{\beta_1}}^2$. Using this R code

alpha <- c()
beta <- c()
for(i in 1:1000){
    dat<-matrix(rnorm(200000,0,1),nrow=100000,byrow=T)
    cov <- chol(matrix(c(5,2,2,5),byrow=T,nrow=2))
    dat <- dat%*%cov
    X1 <- matrix(0,nrow=100000,ncol=2)
    X1[,1] <- 1
    X1[,2] <- dat[,2]
    alpha <- cbind(alpha,solve(t(X1)%*%(X1))%*%t(X1)%*%dat[,1])
    beta <- cbind(beta,solve(dat[,2]%*%dat[,2])%*%dat[,2]%*%dat[,1])
}

I got

> var(alpha[2,])
[1] 8.310618e-06
> var(beta[1,])
[1] 8.309166e-06

However, your model is equal to $y = Xb+e$, where $\sigma_{b}^2 = (X^{T}X)^{-1}\sigma_e^2$. Since $\sigma_e^2 = \sigma_y^2 - \sigma_{xy} \frac{\sigma_{xy}}{\sigma_x^2}$, which is equal to $\frac{\sigma_y^2}{\sigma_x^2}(1-r^2)$, $\sigma_{b_1}^2 = \frac{\sigma_e^2}{\sum_i^n x_i^2}$, and $\sigma_{b_0}^2 = \frac{\sigma_e^2}{n}$

If you expand the above script for calculating $e$, and storing $e^{T}e$, $y^{T}y$, $x^{T}x$ and $x^{T}y$, you can calculate all these values.

Cheers

$\endgroup$
  • 2
    $\begingroup$ Thank you for your answer! $var(\hat \alpha_1) = var(\hat \beta_1)$ is indeed true for big sample sizes. You used a bis sample with n=100,000 in your example. When using a tiny sample (for instance $n=20$) and a higher number of iterations (for instance 20,000) the difference is observable: > var(alpha[2,]) [1] 0.0502148 > var(beta[1,]) [1] 0.04727873 $\endgroup$ – sebastianb Jul 1 '14 at 17:31

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.