I have a network data set on which I am running product diffusion simulations. I have three factors - A (10 levels), B (10 levels) and C (2 levels) producing 200 factor level combinations. For each factor level combination, I am running the diffusion simulation (experiment) 10 times (not sure if we can call it repetition?); each time measuring the total nodes activated ($Y$).

My model is like so:

model1 <- with(data, aov(Y ~ factor(A) * factor(B) * C))
model2 <- with(data, aov(Y^-1 ~ factor(A) * factor(B) * C))

(A and B are numeric, but are really factors in my model)

Model2 is what I am really looking to evaluate. $Y$ is the outcome of the simulation but $\frac{1}{Y}$ is my real interest variable.

The enriching discussion here that presents several alternatives for heteroskedastic, non-normal response variables, with which I am trying to reconcile. The problem as per my diagnosis is in the vast difference in the variance of the response against the factor combinations.

The residual plots from R are as follows:

Model 1: Model1

Model2: Model2

Surely, both models fail the Levene Test for homogeneity of variance.


Levene's Test for Homogeneity of Variance (center = median)
        Df F value    Pr(>F)    
group  199  10.175 < 2.2e-16 ***


Levene's Test for Homogeneity of Variance (center = median)
    Df F value    Pr(>F)    
group  199  10.968 < 2.2e-16 ***

Question 1 Since this robustness check was done after the ANOVA (which was highly significant), I am now unsure if I have a valid ANOVA result or not (given the violation in the assumptions). No matter what transformation I do to $Y$ (log, log(log), sort, cube root), I still cannot clear Levene's Test.

Analysis of Variance Table

Response: Y
                           Df     Sum Sq      Mean Sq   F value    Pr(>F)    
factor(A)                  9      8.7866e+13 9.7628e+12 2498010.3 < 2.2e-16 ***
factor(B)                  9      1.3863e+14 1.5403e+13 3941264.3 < 2.2e-16 ***
C                          1      1.9710e+12 1.9710e+12  504321.5 < 2.2e-16 ***
factor(A):factor(B)        81     5.4508e+13 6.7294e+11  172184.2 < 2.2e-16 ***
factor(A):C                9      2.9917e+11 3.3241e+10    8505.3 < 2.2e-16 ***
factor(B):C                9      6.2798e+11 6.9775e+10   17853.3 < 2.2e-16 ***
factor(A):factor(B):C      81     1.2525e+12 1.5463e+10    3956.6 < 2.2e-16 ***
Residuals                1800     7.0348e+09 3.9082e+06            

However, following the discussion here the Huber-White adjustment produces similar results for model1 but throws an error for model2. I am not sure what is the reason here?

car::Anova(model1, white.adjust = TRUE)
car::Anova(model2, white.adjust = TRUE)

Analysis of Deviance Table (Type II tests)

Response: Y
                           Df      F    Pr(>F)    
factor(A)                  9 1.2607e+07 < 2.2e-16 ***
factor(B)                  9 3.6641e+08 < 2.2e-16 ***
C                          1 2.4293e+05 < 2.2e-16 ***
factor(A):factor(B)       81 5.9477e+06 < 2.2e-16 ***
factor(A):C                9 9.9393e+03 < 2.2e-16 ***
factor(B):C                9 1.9930e+04 < 2.2e-16 ***
factor(A):factor(B):C     81 6.1255e+03 < 2.2e-16 ***
Residuals                1800               

Error in Anova.lm(mod, ...) : residual sum of squares is 0 (within rounding error)

Question 2 Given that both main effects and the interactions are significant, and the ANOVA is not, is it okay to report effect sizes ($\eta^2$ or $\omega^2$) as the primary result? I am worried that the model doesn't really clear the requirements to ascertain the validity of ANOVA, but then the formulae for effect sizes (ratio of sum of squares) do not make any distributional assumptions and to my understanding are more of an accounting of the sum of squares according to the factors.

What I wish for is to find some way to establish the validity of ANOVA (bootstrapping?) and then use this result to present the effect sizes as the main result.

Kindly point out flaws in logic/possible pointers to aid analysis.

  • $\begingroup$ Can you perform a shapiro-wilk test for normality on the residuals and report the p-value? I want to know how significantly non-normal this is model is. $\endgroup$
    – Matt L.
    Commented Jul 28, 2017 at 4:33
  • $\begingroup$ Shapiro-Wilk normality test data: residuals(model1) W = 0.79834, p-value < 2.2e-16 data: residuals(model2) W = 0.59076, p-value < 2.2e-16 $\endgroup$
    – buzaku
    Commented Jul 29, 2017 at 3:14
  • 1
    $\begingroup$ Note that this is no assumption that requires non-normal response variables. The assumption is that you have normal errors. $\endgroup$ Commented Jan 12, 2019 at 18:09

3 Answers 3


As I lay out in this answer, doing hypothesis tests on simulated data doesn't make much sense (other than as a way to evaluate a statistical method). To quote that answer:

The null hypothesis of no difference between treatments is false by design. The p-value is therefore meaningless to begin with.

Note that you can still use the ANOVA framework to partition the variance explained by the different factors - it's the p-values that are useless.

If you were working with numerical predictors and were concerned about nonlinearities in the data, you could use flexible methods like GAMs or random forests to describe variation in the data. Since all your predictors are categorical (if I understand correctly), this would not be very useful to you.


This is a situation for semiparametric regression models. Some resources for those are here.


Okay, for closure on the bootstrapping method for effect sizes (based heavily on the code here (the bootstrapped ANOVA is mentioned in the link)

if (!require("lsr")) {install.packages("lsr"); library("lsr")}
if (!require("boot")) {install.packages("boot"); library("boot")}

ws <- data.frame(A = c(rep(0.02, 10), rep(0.4, 10), rep(0.8, 10)),
                 B = c(rep(2, 10), rep(6, 10), rep(11, 10)),
                 C = c(rep("a", 15), rep("b", 15)),
                 Y = sample(x = 30:9000, size = 30, replace = TRUE))

es_boot <- function(formula, data, i) {

  data_resamp <- data[i,] # Resample rows

  model <- aov(as.formula(formula), data = data_resamp)

  es <- lsr::etaSquared(model, type = 1)[, "eta.sq"] # Extract effect sizes

  return (as.vector(es))

esboot_run <- boot(statistic = es_boot,
                   formula = "Y ~ factor(A) * factor(B) * C", 
                   data = ws, 
                   R = 2000)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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