How do I do an F-test to compare nested linear models in Python?

Viewed 8255

I want to compare two nested linear models, call them m01, and m02 where m01 is the reduced model and m02 is the full model. I want to do a simple F-test to see if the full model adds significant utility over the reduced model.

This is very simple in R. For example:

mtcars <- read.csv("https://raw.githubusercontent.com/focods/WonderfulML/master/data/mtcars.csv")
m01 <- lm(mpg ~ am + wt, mtcars)
m02 <- lm(mpg ~ am + am:wt, mtcars)
anova(m01, m02)

Gives me the following output:

enter image description here

Which tells me that adding the am: wt interaction term significantly improves the model. Is there a way to do something similar to this in Python/sklearn/statsmodels?

Edit: I looked at this question before posting this one and can not figure out how they are the same. The other question is doing an F-test on two vectors. This question is about comparing 2 nested linear models.

I think this is what I need:

http://scikit-learn.org/stable/modules/generated/sklearn.feature_selection.f_regression.html#sklearn.feature_selection.f_regression

but am not sure what exactly to pass this function. If anyone could provide or point to an example, that would be extremely helpful.

2 Answers

Adapting Jeremy's answer in the following way allowed me to get the same result I obtained in R:

import pandas as pd
from statsmodels.formula.api import ols
from statsmodels.stats.anova import anova_lm

cars_df = pd.read_csv("https://raw.githubusercontent.com/focods/WonderfulML/master/data/mtcars.csv")
m01 = ols('mpg ~ am + wt', data=cars_df).fit()
m02 = ols('mpg ~ am + wt + am:wt', data=cars_df).fit()
anovaResults = anova_lm(m01, m02)
print(anovaResults)

This gave me the following results in my jupyter notebook:

enter image description here

I also got these rather cryptic errors:

enter image description here

Anyone have a clue as to what is generating these errors?

I found this book helpful ("An introduction to statistics with python" / Thomas Haslwanter)

Here is the relevant code sample:

import pandas as pd
from statsmodels.formula.api import ols
from statsmodels.stats.anova import anova_lm

df = pd.DataFrame(data, columns=['value', 'treatment'])
model = ols('value ~ C(treatment)', data=df).fit()
anovaResults = anova_lm(model)
print(anovaResults)

Look at the book link above for the print result.

Note: anova_lm sometimes used with the parameter 'typ', try its different possible values to see which suitable for you.

Related