When you have heteroskedasticity, it doesn't make sense to try to check normality of the entire set of residuals, though you could still check groups individually (with corresponding loss of power of course).
On the other hand, it doesn't really make sense to formally test either normality or heterosckedasticity when checking assumptions, since the hypothesis tests answer the wrong question.
This is because your data aren't actually normal (and it's also very unlikely that your populations have identical variance) - so you already know the answer to the question the hypothesis test checks for. With a nice large sample like you have, the chance that a nice powerful test like the Shapiro-Wilk doesn't pick it up is small - so you'll reject as non-normal data from distributions that will have little impact on the signficance level or the power. That is, you'll tend to reject normality - even at quite small significance levels - when it really doesn't matter. The test is likely to reject when it matters least (i.e. when you have a big sample).
What you actually want to know is the answer to a different question than the test answers - something like "How much does this affect my inference?". The hypothesis test doesn't address that question - it doesn't consider the impact of the non-normality, only your ability to detect it.
Further, when you have a sequence "do this equal-variance normal theory analysis if I don't reject all these tests, otherwise do this analysis if I reject that one, this other analysis if I reject the other one and something else if I reject both" we must consider the properties of the whole scheme. Such programs of testing usually do worse (in terms of test bias, accuracy of significance level and very often, of power) than just assuming you'd rejected both tests.
So you recommend to just stick with the ANOVA and consider no assumptions violated?
Not quite. In fact, if anything my last sentence above suggests that you assume heteroskedasticity and normality are violated from the start. Either one alone being violated is relatively easy to deal with, both together is a little trickier (but still possible). However, in your case I think you're probably okay, since I think you needn't worry about one of the two:
Normality may not be such a problem - the considerations would be what kind of non-normality might you have, and how strongly non-normal, and with how large a sample size?
Your sample size seems reasonably large and the distribution pretty mildly left skew and light-tailed, though that assessment may be confounded with the heteroskedasticity. However, if you had good understanding of the properties of what was being measured - which you may well have - or information from similar studies, you might be able to make an a priori assessment on that basis and so better able to choose an appropriate procedure (though I'd still suggest diagnostic checking).
Since your data are probability estimates, they'll be bounded. In fact the left skewness may simply be caused by some probability assessments getting relatively close to 100. If that's the case, you must also tend to doubt your assumption of linearity, and that will be a likely cause of heteroskedasticity as well. If my guess about getting close to the upper bound is right you'll tend to see lower spread among the groups with higher mean.
You might consider an analysis suited to continuous proportions, perhaps a beta-regression - at least if you have no data exactly on either boundary. (An alternative might be a transformation, but models that deal with the data you have tend to be both easier to defend and more interpretable.)
With your decent-sized sample, you are probably safe enough on non-normality, but heteroskedasticity might be more of an issue - in particular, heteroskedasticity issues don't decrease with sample size.
On the other hand, if your sample sizes are equal (or at least very nearly equal), heteroskedasticity is of little consequence. Your tests will be little impacted in the case of equal sample sizes.
If equal-sample-sizes are not the case, I suggest you:
i) don't assume heteroskedasticity will simply be okay
ii) don't formally test it, for the same reasons outlined above (testing answers the wrong question)
Instead I suggest you start with the assumption that the variances differ - whether it's to use something like the Welch approach (I can't say I know how that works for 2x2 off the top of my head, but it should be quite possible to make it work in that case, since it only affects the calculation of residual variance and its df), or to implement your ANOVA in a regression and move to something like heteroskedasticity-consistent standard errors (which is used more widely in areas like econometrics).
As you know, $W$ is a test statistic. In most cases (all consistent tests), a test statistic is not a suitable effect estimator as the statistic reflects the sample size whereas the effect estimator shall be independent of it. Just think of an asymptotic test to test zero mean under the central limit theorem: The approximate distribution is the same for all $n$, so the test statistic contains even all the information about the sample size. That makes the test statistic unsuitable as effect estimater.
For $W$, it is similar (although the approximate distribution depends on the sample size as well). The lower bound for $W$ is $\frac{a_1^2n}{(n-1)}$, where $a_1$ depends on is the expectation for the smallest order statistic.
So no, it is no suitable effect estimator at all.
In fact, I think you are not yet sure what you are looking for as the term "effect" is a bit more difficult than in the usual parametric world of one-dimensional parameters. Here, the raw effect of a.s. not being normally distributed is infinite dimensional: Each measurable subset of $\mathbb{R}$ can have a different probability from the normal distribution model. For a one-dimensional effect, you need to weight it somehow and be aware of the consequences of various weights to your intended application. This way you would decide if e.g. a certain bimodal distribution with Gaussian tails is more normal than a certain unimodal distribution with heavy tails. In fact trading the tail behaviour against the non-tail behaviour might be the most relevant question to invent a suitable effect.
Then, if will be much easier to find an estimator for this particular effect.
Best Answer
A) What is the best single index of the degree to which the data violates normality?
B) Or is it just better to talk about multiple indices of normality violation (e.g., skewness, kurtosis, outlier prevalence)?
I would vote for B. Different violations have different consequences. For example, unimodal, symmetrical distributions with heavy tails make your CIs very wide and presumably reduce the power to detect any effects. The mean, however, still hits the "typical" value. For very skewed distributions, the mean for example, might not be a very sensible index of "the typical value".
C) How can confidence intervals be calculated (or perhaps a Bayesian approach) for the index?
I don't know about Bayesian statistics, but concerning classical test of normality, I'd like to cite Erceg-Hurn et al. (2008) [2]:
D) What kind of verbal labels could you assign to points on that index to indicate the degree of violation of normality (e.g., mild, moderate, strong, extreme, etc.)?
Micceri (1989) [1] did an analysis of 440 large scale data sets in psychology. He assessed the symmetry and the tail weight and defined criteria and labels. Labels for asymmetry range from 'relatively symmetric' to 'moderate --> extreme --> exponential asymmetry'. Labels for tail weight range from 'Uniform --> less than Gaussian --> About Gaussian --> Moderate --> Extreme --> Double exponential contamination'. Each classification is based on multiple, robust criteria.
He found, that from these 440 data sets only 28% were relatively symmetric, and only 15% were about Gaussian concerning tail weights. Therefore the nice title of the paper:
I wrote an
R
function, that automatically assesses Micceri's criteria and also prints out the labels:Here's a test for the standard normal distribution, a $t$ with 8 df, and a log-normal:
[1] Micceri, T. (1989). The unicorn, the normal curve, and other improbable creatures. Psychological Bulletin, 105, 156-166. doi:10.1037/0033-2909.105.1.156
[2] Erceg-Hurn, D. M., & Mirosevich, V. M. (2008). Modern robust statistical methods: An easy way to maximize the accuracy and power of your research. American Psychologist, 63, 591-601.