Difference in means statistical test for non-normal, heteroscedastic data? Suppose I have a data set with three or more groups.  After some exploratory data analysis, I find that


*

*The groups do not come from a normal distribution

*the variances within the groups are not equal, and violate the rule of thumb that the largest variance is no more than 4 times the smallest variance.
If I'm interested in a statistical test for the difference in means across these groups, what tests can I use?
From what I understand, ANOVA requires normality (although it tends to be robust against deviations from normality) and homoscedasticity. The Kruskal-Wallis test and the Fisher-Pitman permutation test can deal with the non-normality, but I believe both require homoscedasticity.  Using Welch's ANOVA will help with the unequal variances, but it requires normality.
I'm interested in how this community would suggest moving forward with a difference in means analysis.
 A: This is an interesting thread so I thought I would add my own thoughts to it. If I understand things correctly, you have a situation where you are interested in testing differences between k population means (where k > 2) when the populations are normal but have different spreads. 
As a first resort, you could try applying a variance-stabilizing transformation to the data you collected on the outcome variable (e.g., log transformation if these data are strictly positive).  If the transformation works, you can then apply the standard ANOVA to the transformed outcome data (presuming the other assumptions underlying ANOVA still hold). In general, the more complicated the transformation, the more complicated the ensuing interpretation. Also, in some situations, you may be hard-pressed to find an adequate transformation.
A better approach, in my view, would be to use generalized least squares (gls) regression to analyze the untransformed outcome data. In R, the gls approach is implemented via the gls() function in the nlme package:
install.packages("nlme")
require(nlme)

model <- gls(outcome ~ group,
         data = yourdata,
         weights = varIdent(form = ~1|group))

summary(model)

anova(model)

The gls approach will help you relate the outcome variable (expressed on its original scale) to the group variable, while allowing for the possibility that the error variability is different across levels of the group variable.  The approach will estimate the error variability in each group for you, under the assumption of normality of the errors within each group. More importantly, the gls approach will provide you with a flexible framework for testing a priori contrasts or performing various types of post-hoc multiple comparisons. Checking model diagnostics is also straightforward. 
The bootstrapping approach is also a possibility, though I think it's not as easy to implement as the gls approach, especially when you consider the need to test a priori contrasts or perform post-hoc multiple comparisons. 
This question has been asked before and addressed here, for example: Explicitly modelling unequal variances in a linear model to get better conf/pred intervals?.
