How to validate Anova assumption of gaussian residuals in python

Question

I think I have some fundamental misunderstandings when it comes to the assumptions behind ANOVA.

Lets suppose there is a function that maps y to x with the following equation, y=3x+2. Now suppose we don't know that this relationship is true and set about determining the relation between y and x, as well as some irrelevant categorical variable, say - eye colour, which has two values A and B (We don't know that it is irrelevant). To simulate this data run the following code

 from numpy.random import seed
 from numpy.random import randn
 from matplotlib import pyplot
 import pandas as pd
 import scipy.stats as stats

 seed(1)
 err= randn(50)*.5 # random normal distribution error
 err2= randn(50)*.5 # random normal distribution error
 x=list(range(0,50)) # independent variable
 #Generate results based on Category A
 A=[(i*3+2+e) for i,e in zip(x,err)] # adding a random error, because reality is imperfect.
 B=[(i*3+2+e) for i,e in zip(x,err2)]
 data_tuples = list(zip(x,A,B))
 df=pd.DataFrame(data_tuples, columns=['X','A','B'])

Okay - so now we have a nice dataframe with this information. Now someone says "I think the value is correlated to the Categorical variable". We say, "Don't be absurd". They say "Convince me" So we crank out an ANOVA test. I know that typically a t-test is used instead of ANOVA for two categorical variables, but I have my reasons, and ANOVA should boil down to t-test so it shouldn't matter.

We run the following code

 fvalue, pvalue = stats.f_oneway(df['A'], df['B'])
 pvalue

So we see that pvalue is certainly not <=0.05, so we can't rule out the null hypothesis that says that there is no difference between the values for A and B.

But then they say - yeah, but aren't there assumptions for ANOVA? Specifically - don't you need to make sure that the RESIDUALS are approximately normal?

So I scratch my head and say "Ummm ... okay. Not sure how to do that. Here is my best guess":

 d_melt = pd.melt(df.reset_index(), id_vars=['X'], value_vars=['A', 'B']) -- Get the frame into a WIDE format
 model = ols('value ~ C(variable)', data=d_melt).fit() -- Run ordinary least squares
 w, pvalue = stats.shapiro(model.resid) -- Run shapiro test.
 pvalue

Now we see that pvalue is << p, which means that we can reject the null hypothesis that the residuals are approximately normal

Indeed when we graph the residuals we don't get a curve at all

 pyplot.hist(model.resid)
 pyplot.show()

So what am I doing wrong? I suspect it is the ols that is wrong, but how else do you get the residuals in this case?

Answer 1

Okay - I looked into this some more, and think I figured out what my misunderstanding was. The following links helped:

Normality vs normality of residuals

https://www.theanalysisfactor.com/checking-normality-anova-model/#:~:text=So%20you'll%20often%20see,the%20residuals%20are%20normally%20distributed.&text=All%20GLM%20procedures%20have%20an%20option%20to%20save%20residuals.

Basically the main question underlying all of this was "Why do some people refer to the assumption of normality as being about the values of each group, and others refer to the residuals of the model? Are some people just wrong here?"

The answer seems to be - if you are doing a linear regression then you MUST look at the distribution of the residuals. If you are just comparing the mean of one group to another then you can do either. Why? Because if the values in the groups are normally distributed then so are the residuals. How does one calculate the residuals? Just subtract each value by the mean of the group it is in.

Consequently - my code was doomed to fail because the x values I used in the function were just a sequence of 0-49, which is in no way a normal distribution.

Hopefully this is correct and will help someone else puzzling over these ANOVA assumptions.