When using Linear Models with random covariates, is it the pearson correlation that determines the reduction of the residual variance? Typically, if you have normally distributed dependent variable Y with variance $\sigma_Y^2$ a treatment indicator and a random covariate that is also normally distributed, then when fitting a linear model with ML-estimates the residual variance is proportional to the factor $(1-r^2)\cdot\sigma_Y^2$. Then, $r^2$ should denote the squared Pearson correlation coefficient between variables $Y$ and variable $X$.
Is that correct? The Pearson correlation would still be used, is if the random variable $X$ is not normally distributed but be exponential or binary distributed?
Would somebody please clear those points up with me and maybe give me some references for further reading?
 A: Sure, you can still use squared correlation as an "$R^2$" statistic under nonnormality, even for non-normal $Y$.  The ML estimate might take a different form, perhaps not involving the sample Pearson correlation, but the ordinary Pearson sample correlation will nevertheless be an asymptotically consistent estimate.
Here is a justification for why the true squared Pearson correlation is an "$R^2$" statistic, even under non-normality of $X$ and $Y$.
First, the law of total variance states that if $(X,Y)$ are jointly distributed with finite variance, then
$$Var(Y) = Var\{f(X)\} + E\{\nu(X)\},$$
where
$$ f(x) = E(Y | X=x)$$
and
$$ \nu(x) = Var(Y | X=x).$$
Since $R^2$ is supposed to be the proportion of variance in $Y$ that is explained by $X$, the true $R^2$ can reasonably be defined as
$$ R^2 = \frac{Var\{f(X)\}}{Var\{f(X)\} + E\{\nu(X)\}} = \frac{Var\{f(X)\}}{Var(Y)}.
$$
Now, under the linearity assumption that $E(Y | X=x) = \beta_0 + \beta_1 x$, we have that $$ \beta_1 = \frac{\sigma_{XY}}{\sigma^2_X} = \rho_{XY}\frac{\sigma_Y}{\sigma_X},$$
where
$$ \sigma_{XY} = E\{(Y-\mu_y)\}\{X-\mu_X)\}$$ is the covariance between $X$ and $Y$, and
$$ \rho_{XY} = \frac{\sigma_{XY}}{\sigma_X\sigma_Y}$$ is the correlation between $X$ and $Y$.
Now, back to using the law of total variance, the true $R^2$ is given by
$$R^2 = \frac{Var\{f(X)\}}{Var(Y)} = \frac{Var\{\beta_0 + \beta_1 X\}}{\sigma^2_Y} = \frac{\beta_1^2 \sigma_X^2}{\sigma^2_Y} = \rho_{XY}^2.$$
The question asked about more than one $X$, in which case the ordinary Pearson correlation does not give you the $R^2$ even under normality. However, the above argument generalizes easily to multiple $X$ variables, with a re-definition of the correlation coefficient that involves the $(X_1,X_2,\dots,Y)$ covariance matrix.  Again, normality is not needed.
