I'm conducting a Two-Way ANOVA with my two factors being Sex and Cohort. I have data from two cohorts of subjects, with each cohort consisting of males and females that were measured on one response variable. (Because of some exclusions, there are unequal sample sizes between groups.)

Prior to running the ANOVA, my understanding is that I must test the data for normality and homogeneity of variance (HOV).

  1. Do I test for normality and HOV in each of the four groups separately? (i.e. test for normality in data from cohort 1 males only, then test for normality in data from cohort 1 females only, then cohort 2 males, then cohort 2 females?)

  2. Does the assumption of HOV apply to all four groups, i.e. The null hypothesis is "Cohort 1 male variance = Cohort 1 female variance = Cohort 2 male variance = Cohort 2 female variance?"

  3. I used the Shapiro-Wilk test for normality in each group, and Levene's test of equality of error variances. Unfortunately, in all groups, the data are very non-normal and give highly significant values for Levene's test. I have tried several transformations (square root, log, natural log, square) but nothing has worked to normalize the data so far.

I'm wondering how to proceed? I've read that unlike Welch's test for a one-way ANOVA, there is no good two-way ANOVA equivalent for non-normal data with heterogeneous variances.

Are there any other transformations that could work? If not, would the best option be to simply run the ANOVA, but mention that assumptions were violated that may impact the test results?

EDIT (to add more information):

To clarify, the main issue is lack of homogeneity of variance for the Two-Way ANOVA. I had previously written that the transformations did not work to normalize the data -- I was mistaken (my apologies!). The data were very positively skewed (kurtosis was not really an issue), and the square root transformation successfully normalized the distribution. However, I still have heterogeneous variances. I'm wondering if there's anything I can do to correct this, or if it's acceptable to go ahead with the ANOVA, and explicitly mention the heterogeneous variances in the description of my methods?

EDIT 2 (images added):

Boxplots of untransformed data:

enter image description here

enter image description here

EDIT 3 (raw data added):

**Cohort 1 males (n=12)**: 

**Cohort 1 females (n=12)**: 

**Cohort 2 males (n=11)**: 

**Cohort 2 females (n=11)**:
  • $\begingroup$ You should ideally show us your data if the dataset is not too large; certainly more precise information about the data is needed, e.g. dot plots, box plots or histograms of all four groups, with min, max, mean, SD, median, IQR. Finding a suitable transformation is not shooting in the dark: in particular, if logarithm is a serious candidate, then square cannot be, and vice versa, as they have completely opposite effects. Note that logarithm (base 10, presumably) and natural logarithm are identical in their effects on no-normality.or unequal variances. $\endgroup$
    – Nick Cox
    Feb 12 '16 at 10:11
  • $\begingroup$ In 2) note that the null hypothesis is about the means, not the variances. But indeed it is in the simplest case an assumption (in the usual terminology) that all four variances are the same. (These techniques would be a little easier to understand perhaps if instead of assumptions people talked about ideal conditions.) $\endgroup$
    – Nick Cox
    Feb 12 '16 at 10:14
  • $\begingroup$ You didn't follow up on most of my suggestions. At a wild guess, a stronger transformation such as logarithmic might do more to stabilise variances. I would still rather see data or at least graphs to comment. $\endgroup$
    – Nick Cox
    Feb 12 '16 at 23:36
  • $\begingroup$ Unfortunately I'm not sure how to post the raw data, so I added box plots. It seems like the square root transformation does a better job of stabilizing variances than the log. (The p-value for Levene's test using square-root transformed data is is .04, whereas for log, it is .000 and untransformed, it is .001). $\endgroup$
    – CH.7
    Feb 15 '16 at 0:57
  • $\begingroup$ How big is your data? You should consider non-parametric ANOVA methods. Give Friedman's Test a try. You will lose power, but if you have a lot of data, that shouldn't be a problem and then you won't have to worry about all these assumptions. $\endgroup$ Feb 15 '16 at 6:44

Thanks for posting the data. Posting shows that the box plots concealed, although not intentionally, the sample sizes and important detail too. Whenever I see skewness on a positive response, my first instinct is to reach for logarithms, as they so often work well. Here, however, logarithms drastically over-transform, and plotting everything shows up a small surprise, namely that the two lowest values need care and attention.

enter image description here

The graph here is a quantile-box plot in which the original data points are plotted in order on scales consistent with the box idea (i.e. about half the points are inside the box and about half outside, the "about" being a side-effect of sample sizes like 11).

A more cautious square root transformation seems about right.

enter image description here

Personally I regard preliminary tests for normality and so forth as over-rated stuff left over from the 1960s. I feel far too queasy about forking paths of the form: pass the test OK, fail the test do something quite different, particularly with small sample sizes. Once you have a scale on which you have approximate symmetry and approximate equality of variances, linear models will work well.

Similarly, skewness and kurtosis from small samples can hardly be trusted. (Actually, skewness and kurtosis from large samples can hardly be trusted.) For some of the reasons see e.g. this paper

Indeed, some fits with generalised linear models with cohort and gender as indicator predictor variables show that results seem consistent over identity, root and log links, even despite the evidence of the first graph. If this were my problem I would push forward with a square root link function. In other words, although transformations are informative about the best scale to work on, you let the link function of a generalised linear model do the work.

Campaign slogan: Conventional box plots with a few groups leave out detail that could easily be interesting or useful and don't make full use of the space available. Use graphs that show more!


Here is token output: predicted values using generalised linear model, root link, normal family, interaction between cohort and females:

  | cohort   females   predicted   Freq. |
  |      1     males       2.056      12 |
  |      1   females       5.024      12 |
  |      2     males      12.712      11 |
  |      2   females      15.348      11 |
  • $\begingroup$ Predicted values using GLM that includes the interaction term should be simply equal to group means? Or wait, is it only true for linear link? (Your numbers do seem to equal group means in this case though.) So just to clarify: what you are advocating here is not exactly the root transformation followed by ANOVA, it is GLM with root link? $\endgroup$
    – amoeba
    Feb 16 '16 at 17:38
  • $\begingroup$ Indeed, but it's a relevant comparison with Frank's result. I don't give the full apparatus of P-values, etc. $\endgroup$
    – Nick Cox
    Feb 16 '16 at 17:42
  • $\begingroup$ GLM with root link: indeed, that's explicit in my answer. $\endgroup$
    – Nick Cox
    Feb 16 '16 at 17:51
  • $\begingroup$ Yes, I see it now. Thanks (+1 by the way). It is just that OP was talking about a transformation followed by ANOVA, that is perhaps why I was thinking along these lines still. $\endgroup$
    – amoeba
    Feb 16 '16 at 17:53
  • $\begingroup$ My stance is that GLMs let you have it both ways, get results for your original measurement scale, yet fit in a better space. The OP quotes up to 10 significant figures, although at a wild guess the numbers are all quotients of some kind, so the implication is that the data are thought of as measurements. $\endgroup$
    – Nick Cox
    Feb 16 '16 at 17:57

When one does not know the transformation in advance, you can hope to get lucky by trying one transformation in addition to the identity transformation (leaving things alone). Trying a total of two transformations is probably not very harmful to statistical inference. But in general we don't know enough about the transformation and all this brings arbitrariness to the analysis and uncertainty in how to control type I error and confidence interval coverage. For these and other reasons I am recommending more and more that semiparametric models be standard choices. For this problem the proportional odds (PO) ordinal logistic model, with no binning of $Y$, is a good choice. This is a generalization of Wilcoxon-Mann-Whitney-Kruskal-Wallis tests. The PO model is transformation-invariant except when using it to estimate $E(Y|X)$. It is robust and competitive with normal-theory methods even if normality holds. The R rms package orm function efficiently handles ordinal models for continuous $Y$.

Here is an ordinal analysis using the R rms package. I have included an interaction between cohort and sex. New: a plot of the estimated underlying conditional distribution of y is added.

d1 <- data.frame(cohort='one', sex='male', y=c(.476,
d2 <- data.frame(cohort='one', sex='female', y=c(4.548333,
d3 <- data.frame(cohort='two', sex='male', y=c(7.9645,
d4 <- data.frame(cohort='two', sex='female', y=c(11.2865,
d <- rbind(d1, d2, d3, d4)
dd <- datadist(d); options(datadist='dd')

# Fit the default ordinal model (prop. odds)
f <- orm(y ~ cohort * sex, data=d)

Logistic (Proportional Odds) Ordinal Regression Model

orm(formula = y ~ cohort * sex, data = d)
                      Model Likelihood          Discrimination          Rank Discrim.    
                         Ratio Test                 Indexes                Indexes       
Obs            46    LR chi2      58.46    R2                  0.720    rho     0.854    
Unique Y       46    d.f.             3    g                   3.502                     
Median Y  6.68525    Pr(> chi2) <0.0001    gr                 33.176                     
max |deriv| 0.002    Score chi2   52.40    |Pr(Y>=median)-0.5| 0.410                     
                     Pr(> chi2) <0.0001                                                  

                        Coef    S.E.   Wald Z Pr(>|Z|)
cohort=two               6.8607 1.3333  5.15  <0.0001 
sex=female               2.6922 0.8680  3.10  0.0019  
cohort=two * sex=female -1.8481 1.1579 -1.60  0.1105  

            Wald Statistics          Response: y 

 Factor                                      Chi-Square d.f. P     
 cohort  (Factor+Higher Order Factors)       28.92      2    <.0001
  All Interactions                            2.55      1    0.1105
 sex  (Factor+Higher Order Factors)          10.82      2    0.0045
  All Interactions                            2.55      1    0.1105
 cohort * sex  (Factor+Higher Order Factors)  2.55      1    0.1105
 TOTAL                                       32.59      3    <.0001

# Show intercepts as a function of y to estimate the underlying
# conditional distribution.  Result: more uniform than Gaussian
alphas <- coef(f)[1 : num.intercepts(f)]
yunique <- f$yunique[-1]
plot(yunique, alphas)
# Compare to distribution of residuals
plot(ecdf(resid(ols(y ~ cohort * sex, data=d))), main='')

Intercepts vs. y and residuals from OLS

M <- Mean(f)
# Confidence intervals for means are approximate
# Confidence intervals for odds ratios or exceedance probabilities
# are correct for ordinal models
Predict(f, cohort, sex, fun=M)

  cohort    sex      yhat      lower     upper
1    one   male  2.051195  0.7412913  4.029275
2    two   male 13.089852  8.7310555 17.054696
3    one female  5.261155  3.7446728  7.000745
4    two female 14.884409 10.3247910 18.616770

Response variable (y):  

Limits are 0.95 confidence limits

# Ordinary sample means with t- confidence limits:
with(d, summarize(y, llist(cohort, sex), smean.cl.normal))
  cohort    sex         y      Lower     Upper
2    one   male  2.055708  0.8934179  3.217999
1    one female  5.024132  3.7586617  6.289602
4    two   male 12.711545  9.6236006 15.799490
3    two female 15.348114 13.1603031 17.535924
  • 1
    $\begingroup$ +1 Frank. I Wholehearted, agree. This this is, by far the most straight-forward approach. Iterating through these transformations in hopes of getting the "right" one, without carefully thinking though the impact on errors, make me nervous. This approach remedies this. $\endgroup$ Feb 15 '16 at 20:37
  • $\begingroup$ This is without doubt an intriguing approach. But arguably it solves questions about measurement scale by ignoring them. From some points of view that's not a weakness. I wonder how easily its conclusions map back to the science from which these data arise, whatever they are (just anonymous "RV", so far as we are told, apart from there being two genders). $\endgroup$
    – Nick Cox
    Feb 15 '16 at 20:47
  • 1
    $\begingroup$ @amoeba Indeed, but Frank's proposal, as I understand it, is to treat it as ordered with as many distinct levels as are present in the response. Disclosure: I have yet to read the 2/e of his book Regression modeling strategies. $\endgroup$
    – Nick Cox
    Feb 15 '16 at 21:54
  • 1
    $\begingroup$ All continuous interval-level variables are ordinal. Think of how we use the Wilcoxon 2-sample test. It only assumes ordinality, not interval level, and works perfectly well even with no ties in the data. $\endgroup$ Feb 15 '16 at 22:44
  • 1
    $\begingroup$ This is a fantastic example of ordinal logistic regression. Something similar can be achieved in Stata using the ologit command, though it lacks some of the extensive conveniences the R packages rms and Hmisc bring. $\endgroup$ Feb 16 '16 at 21:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.