title: Normal probability plot for multiple factors
body: <p>I'm trying to run diagnostics for normality in a 2^4 factorial problem with two replicates. Here is my code: </p>
  
  <pre><code>n = 2
  
  A &lt;- factor(c(rep("-", 1*n), rep("+", 1*n)))
  B &lt;- factor(c(rep("-", 2*n), rep("+", 2*n)))
  C &lt;- factor(c(rep("-", 4*n), rep("+", 4*n)))
  D &lt;- factor(c(rep("-", 8*n), rep("+", 8*n)))
  
  obs &lt;- c(90, 93,
           74, 78,
           81, 85,
           83, 80,
           77, 78,
           81, 80,
           88, 82,
           73, 70,
           98, 95,
           72, 76,
           87, 83,
           85, 86,
           99, 90,
           79, 75,
           87, 84,
           80, 80)
  
  df &lt;- data.frame(A, B, C, D, obs)
  
  model &lt;- aov(obs ~ A*B*C*D, data = df)
  
  summary(model)
  
  
  par(mfrow=c(1,2))
  
  qqnorm(resid(model), ylab = "Residuals", xlab = "Quantiles", pch = 16)
  qqline(resid(model))
  
  plot(resid(model) ~ fitted(model), ylab = "Residual", xlab = "Predicted",            pch = 16)
  abline(0,0)
  </code></pre>
  
  <p>The ANOVA table is giving me the correct values, but when I analyze the normality conditions using a Normal Q-Q plot, it incorrectly gives me symmetric residuals. I have noticed that I only run into this issue when I am analyzing four or more interactions. All the residual plots for three interaction or less has the correct expected output with the same code.</p>
  
  <p>Any help would be greatly appreciated </p>