Does the F-test for multivariable regression work with non-normal residuals but large sample size? Suppose I am building a multiple regression model, perhaps with 5 explanatory variables. Suppose the residuals are not normally distributed (based on a Q-Q plot and a D'Agostino test, due to kurtosis). Suppose I have a sample size in the neighborhood of 100: small enough that I can trust the D'Agostino test to be telling me about a significant departure from normality, but large enough that I can lean on the Central Limit Theorem. 

Can I trust the ANOVA F-statistic and the corresponding p-value given by R? In other words, is there a theoretical reason that, with a sufficiently large sample size (much bigger than 100 if you want), that the ANOVA F-statistic should be F-distributed even if the residuals are not normal? 

This is claimed by Wooldridge on page 176 of his Introductory Econometrics book (4th edition). Wooldridge is usually fairly careful and gives theoretical justifications, but does not in this case. I am a mathematician, so I'd like to understand why the F-statistic should be F-distributed in this setting. I know that the t-statistics for the individual explanatory variables (assuming no multicollinearity) are approximately normally distributed by the Central Limit Theorem. I know that I can do a randomization based test and build an empirical distribution for the F-statistics, and it looks like an F-distribution. I just don't have a theoretical reason to believe it is one, and my students keep asking me this question. Thanks!
 A: I am not entirely sure this answers your question, but here we go. Asymptotically, $q\times F$, with $F$ an F-distributed r.v. with $q$ "numerator degrees of freedom corresponding to the number of restrictions tested in the F-test, converges in distribution to a $\chi^2_q$ r.v. as the denominator degrees of freedom (a function of $n$, the sample size) tend to infinity.
Rearranging this result gives that $F$ itself converges to a $\chi^2_q/q$-distributed r.v., or a $F_{q,\infty}$ r.v. This r.v. can be well approximated by a $F_{q,n}$ r.v. for $n$ "sufficiently large" (as with most asymptotic approximations, there is no general answer to what "sufficiently large" precisely is), just as a $t$-distributed r.v. can be well approximated by a standard normal when degrees of freedom of $t$ are sufficiently large.
The figure below shows that even for 25 denominator degrees of freedom, the approximation is already quite close, for $q=3$.

EDIT:
To (hopefully) address Glenb's valid comment: why would $F$ be approximately $\chi^2_q/q$? We may write $F$ as
$$
F=\frac{(R'\hat{\beta}-r)'\left\{R'(X' X)^{-1}R\right\}^{-1}(R'\hat{\beta}-r)/q}{\hat{\sigma}^{2}},
$$
where $R$ and $r$ specify the null $H_0:R'\beta=r$ to be tested.
Let us take the leading case of $R=I$, so testing restrictions on the coefficients as opposed to, say, testing restrictions on the sum of coefficients. If the null is true, $r=\beta$, so that we may write
$$
qF=\frac{(\hat{\beta}-\beta)'\left\{(X' X)^{-1}\right\}^{-1}(\hat{\beta}-\beta)}{\hat{\sigma}^{2}}=\sqrt{n}(\hat{\beta}-\beta)'\left\{\hat{\sigma}^{2}(X' X/n)^{-1}\right\}^{-1}\sqrt{n}(\hat{\beta}-\beta)
$$
Now, we know that $\sqrt{n}(\hat{\beta}-\beta)$ is asymptotically normal with an asymptotic variance that may be estimated consistently by $\hat{\sigma}^{2}(X' X/n)^{-1}$. So, the middle term $\bigl\{\hat{\sigma}^{2}(X' X/n)^{-1}\bigr\}^{-1}$, loosely speaking, "standardizes" both asymptotic multivariate normal random vectors in the quadratic form by their "standard deviation", so that we (again asymptotically) obtain standard normal vectors. 
In large but finite samples, this vector will only approximately be standard normal. Hence, the resulting sum of squares of the standard normal independent entries of the vector will also only approximately be $\chi^2_q$. 
But for $n$ reasonably large, the approximation will be pretty good. 
