Independent samples t-test: Do data really need to be normally distributed for large sample sizes? Let's say I want to test if two independent samples have different means. I know the underlying distribution is not normal.
If I understand correctly, my test statistic is the mean, and for large enough sample sizes, the mean should become normally distributed even if the samples are not. So a parametric significance test should be valid in this case, right? I have read conflicting and confusing information about this so I would appreciate some confirmation (or explanation why I'm wrong).
Also, I've read that for large sample sizes, I should use the z-statistic instead of the t-statistic. But in practice, the t-distribution will just converge to the normal distribution and the two statistics should be the same, no?
Edit: Below are some sources describing the z-test. They both state that the populations must be normally distributed:
Here, it says "Irrespective of the type of Z-test used it is assumed that the populations from which the samples are drawn are normal."
And here, the requirements for the z-test are listed as "Two normally distributed but independent populations, σ is known".
 A: I think this is a common misunderstanding of the CLT. Not only does the CLT have nothing to do with preserving type II error (which no one has mentioned here) but it is often not applicable when you must estimate the population variance. The sample variance can be very far from a scaled chi-squared distribution when the data are non-Gaussian, so the CLT may not apply even when the sample size exceeds tens of thousands. For many distributions the SD is not even a good measure of dispersion.
To really use the CLT, one of two things must be true:  (1) the sample standard deviation works as a measure of dispersion for the true unknown distribution or (2) the true population standard deviation is known. That is very often not the case. And an example of n=20,000 being far too small for the CLT to "work" comes from drawing samples from the lognormal distribution as discussed elsewhere on this site.
The sample standard deviation "works" as a dispersion measure if for example the distribution is symmetric and does not have tails that are heavier than the Gaussian distribution.
I do not want to rely on the CLT for any of my analyses.
A: I am leaving this paragraph for the comments to make sense: Probably the assumption of normality in the original populations is too restrictive, and can be forgone focusing on the sampling distribution, and thanks to the central limit theorem, especially for large samples. 
Applying the $t$ test is probably a good idea if (as is usually the case) you don't know the population variance, and you are instead using the sample variances as estimators. Note that the assumption of identical variances may need to be tested with an F test of variances or a Lavene test before applying a pooled variance - I have some notes on GitHub here.
As you mention, the t-distribution does converge to the normal distribution as the sample increases, as this quick R plot demonstrates:

In red is the pdf of a normal distribution, and in purple, you can see the progressive change in the "fat tails" (or heavier tails) of the pdf of the $t$ distribution as the degrees of freedom increase until it finally blends with the normal plot.
So applying a z-test would likely be OK with large samples.

Addressing the issues with my initial answer. Thank you, Glen_b for your help with the OP (the likely new mistakes in interpretation are entirely mine).


*

*THE T STATISTIC FOLLOWS A T DISTRIBUTION UNDER NORMALITY ASSUMPTION:
Leaving aside complexities in the formulas for one-sample v. two-sample (paired and unpaired), the general t statistic focusing on the case of comparing a sample mean to a population mean is:
$\text{t-test}= \Large \frac{\bar X-\mu}{\frac{s}{\sqrt{n}}}=\large\frac{\frac{\bar{X}-\mu}{\sigma/\sqrt{n}}}{\sqrt{\frac{s^2}{\sigma^2}}} =\displaystyle \large\frac{\frac{\bar{X}-\mu}{\sigma/\sqrt{n}}}{\sqrt{\frac{\frac{\sum_{x=1}^n(X - \bar{X})^2}{n-1}}{\sigma^2}}} \tag1$
If $X$ follows a normal distribution with mean $\mu$ and variance $\sigma^2$:


*

*The numerator of $(1)$ $\sim N(1,0)$.

*The denominator of $(1)$ will be the square root of $\frac{s^2/\sigma^2}{n-1}\sim\frac{1}{n-1}\,\,\chi^2_{n-1}$ (scaled chi squared), since $(n-1)s^2/\sigma^2\sim\chi^2_{n-1}$ as derived here.

*Numerator and denominator should be independent.


Under these conditons the $\text{t-statistic} \sim t(df=n-1)$.


*CENTRAL LIMIT THEOREM:
The tendency towards normality of the sampling distribution of the sample means
as the sample size increases can justify assuming a normal distribution of the numerator even if the population is not normal. However, it does not influence the other two conditions (chi square distribution of the denominator and independence of the numerator from the denominator).
But not all is lost, in this post it is discussed how Slutzky theorem supports the asymptotic convergence towards a normal distribution even if the chi distribution of the denominator is not met.


*ROBUSTNESS:
On the paper "A  More Realistic Look at the  Robustness and  Type II Error Propertiesof the  t Test to Departures From  Population Normality" by Sawilowsky SS and Blair RC in Psychological  Bulletin, 1992, Vol. 111, No. 2, 352-360, where they tested less ideal or  more "real world" (less normal) distributions for power and for type I errors, the following assertions can be found: "Despite the conservative nature with regard to Type I error of the  t test  for some of these real distributions, there was little effect  on the  power levels for the  variety of treatment conditions and  sample sizes studied. Researchers  may easily compensate for the slight loss in power by selecting a slightly larger sample size".
"The prevailing view seems to be that the independent-samples t test is reasonably robust, insofar as Type I errors are concerned,  to non-Gaussian population shape so long as (a) sample  sizes are  equal  or  nearly so, (b) sample sizes are  fairly large (Boneau, 1960, mentions sample sizes of 25 to 30), and (c) tests are two-tailed rather than one-tailed. Note also that when these conditions are met and  differences  between nominal alpha  and  actual alpha  do  occur, discrepancies are usually of a conservative  rather than  of a liberal nature."
The authors do stress the controversial aspects of the topic, and I look forward to working on some simulations based on the lognormal distribution as mentioned by Professor Harrell. I would also like to come up with some Monte Carlo comparisons with non-parametric methods (e.g. Mann–Whitney U test). So it's a work in progress...

SIMULATIONS:
Disclaimer: What follows is one of these exercises in "proving it myself" one way or another. The results cannot be used to make generalizations (at least not by me), but I guess I can say that these two (probably flawed) MC simulations don't seem to be too discouraging as to the use of the t test in the circumstances described.
Type I error:
On the issue of type I errors, I ran a Monte Carlo simulation using the Lognormal distribution. Extracting what would be considered larger samples ($n=50$) many times from a lognormal distribution with parameters $\mu=0$ and $\sigma=1$, I calculated the t-values and p-values that would result if we were to compare the means of these samples, all of them arising from the same population, and all of the same size. The lognormal was chosen based on the comments and the marked skewness of the distribution to the right:

Setting a significance level of $5\%$ the actual type I error rate would have been $4.5\%$, not too bad...
In fact the plot of the density of the t tests obtained seemed to overlap the actual pdf of the t-distribution:

The most interesting part was looking at the "denominator" of the t test, the part that was supposed to follow a chi-squared distribution:
$$(n-1)s^2/\sigma^2=98\,\frac{(49 \,  (\text{SD}_A^2 + \text{SD}_A^2))/98} {(e^{\sigma^2}-1) \, e^{2\mu+\sigma^2}}$$.
Here we are using the common standard deviation, as in this Wikipedia entry:
$$S_{X_1X_2}=\sqrt{\frac{(n_1 -1)\,S_{X_1}^2 + (n_2 -1)\,S_{X_2}^2}{n_1+n_2-2}}$$
And, surprisingly (or not) the plot was extremely unlike the superimposed chi-squared pdf:

Type II Error and Power:
The distribution of blood pressure is possible log-normal, which comes extremely handy to set up a synthetic scenario in which the comparison groups are separate in average values by a distance of clinical relevance, say in a clinical study testing the effect of a blood pressure drug focusing on the diastolic BP, a significant effect could be considered an average drop of $10$ mmHg (a SD of approximately $9$ mmHg was chosen):

Running comparison t-tests on an otherwise similar Monte Carlo simulation as for type I errors between these fictitious groups, and with a significance level of $5\%$ we end up with $0.024\%$ type II errors, and a power of only $99\%$.
The code is here.
