<jupyter_start><jupyter_text>Linear regression with Professor Mittens, a.k.a. recipe for linear regression. OverviewIn this notebook we will learn how to use regression to study the factors that affect the number of pats cats will recieve. This will start with a visual inspection of the data, followed by the development of a linear model to explain the data. Along the way we will answer a few questions such as: does coat colour influence the number of pats, is a long coat better than a short coat, and how important is the volume of a meow. Specifying regression modelsA very popular way to describe regression models is with "formulae" as popularised by R. The [R documentation on formulae](https://cran.r-project.org/doc/manuals/R-intro.htmlFormulae-for-statistical-models) is a good place to learn how to use these properly. For example, here is the syntax we will use today,- `y ~ x1 + x2` will make a linear model with the predictors $x_1$ and $x_2$.- `y ~ x1 * x2` includes the terms $x_1 + x_2 + x_1x_2$- `y ~ x1 : x2` includes *just* the interaction term $x_1x_2$- `y ~ 1 + x` explicitly states that there should be an intercept term- `y ~ 0 + x` explicitly states that there should _not_ be an intercept term- `y ~ .` in R would mean all the variables but not in statsmodels- `y ~ C(x)` specifies that $x$ is a catagorical variable **NOTE** this is not necessary in R.<jupyter_code>%matplotlib inline
import pandas as pd
import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf
import altair as alt
from functools import reduce<jupyter_output><empty_output><jupyter_text>Helping cats get more patsProfessor Mittens in interested in helping cats optimise the number of pats they can get. To learn more about this, he has interviewed 1000 cats and taken measurements of their behaviour and appearance. The data in `cat-pats.csv` contains measurments of the following:- `time_outdoors` is the number of hours that the cat is out of their primary dwelling,- `coat_colour` is either tortoiseshell, white, or "other" encoded as integers 1, 2, and 3 respectively,- `weight` is the weight of the cat in kilograms,- `height` is their height in centimeters,- `loudness` is a measure of how loud their meow is, the units are not known,- `whisker_length` is the length of their whiskers in centimeters,- `is_longhaired` is a Boolean variable equal to 1 if the cat is of a longhaired breed and 0 if it is of a shorthaired breed,- `coat_length` is the length of their fur in centimeters,- and `num_pats` is the number of pats they received on the day they were interviewed.The variable we are interested in explaining is `num_pats`. Although this is a discrete variable, we will ignore this aspect of the data and consider it as a continuous value. This is a useful simplifying assumption, as you learn more about regression, in particular generalized linear models, you will see additional ways to handle this. For this example, you can consider it a continuous variable though.The types of questions that Professor Mittens is interested in answering are as follows:1. Do any of the variables correlate with the number of pats that the cats recieve?2. Under a naive model, how much of the variability in pats can they explain? Do all the variables need to be included?3. Does the coat colour matter?4. Among short-haired breeds they say longer hair is better, among long-haired breeds they say short hair is better, who is correct?5. **If a cat can choose to spend more time outdoors, or practise meowing louder, which will get them more pats?**  Read in the data and generate some scatter plots to see if there are any good predictors of the number of patsThe data is in the file `cat-pats.csv` so read this into a data frame using `pd.read_csv` and go from there. I have used altair to generate my scatter plots based on [this example](https://altair-viz.github.io/gallery/scatter_matrix.html) but you can use whatever you feel most comfortable with. It might be useful to use colour to see if `coat_colour` and `is_longhaired` are important.  QuestionBased on these figures, what variables appear to relate to the number of pats? What do you notice about the catagorical variables `coat_colour` and `is_longhaired`?<jupyter_code>cat_df = pd.read_csv("cat-pats.csv")
col_names = cat_df.columns.tolist()
predictor_names = col_names.copy()
predictor_names.remove("num_pats")

alt.Chart(cat_df).mark_circle().encode(
    alt.X(alt.repeat("column"), type='quantitative', scale=alt.Scale(zero=False)),
    alt.Y(alt.repeat("row"), type='quantitative', scale=alt.Scale(zero=False)),
    color = "is_longhaired:O"
).properties(
    width=100,
    height=100
).repeat(
    row=["num_pats"],
    column=predictor_names
)
alt.Chart(cat_df).mark_circle().encode(
    alt.X(alt.repeat("column"), type='quantitative', scale=alt.Scale(zero=False)),
    alt.Y(alt.repeat("row"), type='quantitative', scale=alt.Scale(zero=False)),
    color = "coat_colour:O"
).properties(
    width=100,
    height=100
).repeat(
    row=["num_pats"],
    column=predictor_names
)<jupyter_output><empty_output><jupyter_text>Partial answer- Based on the plots, it looks like the time the cat spends outdoors and their coat length (both whether they are a longhaired breed and the actual length of their fur) are important.- Also note that when you colour based on whether the breed is a long hair some of the data seems to cluster while this is not the case in for the colour of the coat.  Compute the correlation between each variable and the number of pats, what looks important QuestionDoes the the correlation matrix raise any further questions? Does it handle the catagorical variables correctly?<jupyter_code>cat_df.corr()<jupyter_output><empty_output><jupyter_text>Partial answer- There is a strong relationship between whisker length and height which we had not seen in the plots above.- Since the catagorical variables were encoded as numbers they are being handled as qualitative variables so the reported values are not correct for those rows/columns.  What is $R^2$?Sometimes called the *coefficient of determination*, this statistic measures the proportion of the variance in the response variable that is explained by the regression model. In the case of simple linear regression it is just the correlation squared, it can also be calculated as the ratio of the regression sum of squares and the total sum of squares. $$R^2 = \frac{\text{RegSS}}{\text{TSS}}$$It can be thought of as the proportion of the total variance that is explained by the regression model. What is an *adjusted* $R^2$?For a fixed number of observations, as the number of covariates increases you can get explain as much of the variability as you want! The adjusted $R^2$ is a way to penalise using too many covariates. The adjusted $R^2$ for a model with $n$ observations and $p$ coefficients is given by the following:$$\tilde{R}^2 = 1 - \frac{n - 1}{n - p}\left(1 - R^2\right)$$ Under a naive model, how much of the variability in pats can they explain?Run an ordinary linear regression with all of the variables and see what percentage of the variability in the number of pats is explained. Make sure that you have used the catagorical variables correctly. Can be be confident in rejecting the null hypothesis that none of these variables is associated with the number of pats received?<jupyter_code>lm_1 = smf.ols("num_pats ~ time_outdoors + C(coat_colour) + weight + height + loudness + whisker_length + C(is_longhaired) + coat_length", cat_df).fit()
print(lm_1.summary())
sm.graphics.influence_plot(lm_1, criterion="cooks")
sm.graphics.plot_fit(lm_1, "coat_length")
cat_df[543:546]
cat_df.describe()<jupyter_output><empty_output><jupyter_text>Question: Is colinearity an issue in this model? Do all of the variables need to be included?Compute the VIF to see if there is a concerning amount of colinearity between any of the covariates.<jupyter_code>col_names = cat_df.copy().columns.tolist()
col_names.remove("num_pats")

def join_strings(xs, sep):
    return reduce(lambda a, b: a + sep + b, xs)
    
for v in col_names:
    cns = col_names.copy()
    cns.remove(v)
    formula = v + " ~ " + join_strings(cns, " + ")
    coef_det = smf.ols(formula, data = cat_df).fit()
    vif = 1 / (1 - coef_det.rsquared)
    if vif > 3:
        print("\n" + 80 * "=")
        print(v)
        print(vif)
lm_3 = smf.ols("num_pats ~ time_outdoors + C(coat_colour) + weight + height + loudness + C(is_longhaired) + coat_length", cat_df).fit()
print(lm_3.summary())<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:               num_pats   R-squared:                       0.571
Model:                            OLS   Adj. R-squared:                  0.567
Method:                 Least Squares   F-statistic:                     164.7
Date:                Thu, 12 Nov 2020   Prob (F-statistic):          3.82e-176
Time:                        15:16:59   Log-Likelihood:                -3240.3
No. Observations:                1000   AIC:                             6499.
Df Residuals:                     991   BIC:                             6543.
Df Model:                           8                                         
Covariance Type:            nonrobust                                         
=========================================================================================
                            coef    std e[...]<jupyter_text>Partial answerThe VIF of whisker length is >5 and removing it does not reduce the $R^2$ and means that loudness is now closer to being significant. Together this suggests that there was some colinearity there, so removing the whisker length variable seems sensible.  Does coat colour matter?1. Make a box plot of the number of pats by coat colour to see this pattern.2. Fit an additional linear model without the coat colour as a covariate to see how much of the explained variability comes from the inclusion of coat colour in the model.<jupyter_code>coat_df = cat_df.loc[:,["coat_colour", "num_pats"]].copy()
alt.Chart(coat_df).mark_boxplot().encode(
    x='coat_colour:O',
    y='num_pats:Q'
).properties(
    width=600,
    height=300
)

lm_with_colour = smf.ols("num_pats ~ time_outdoors + C(coat_colour) + weight + height + loudness + C(is_longhaired) + coat_length", cat_df).fit()
lm_without_colour = smf.ols("num_pats ~ time_outdoors + weight + height + loudness + C(is_longhaired) + coat_length", cat_df).fit()

a = (lm_with_colour.resid ** 2).sum()
b = (lm_without_colour.resid ** 2).sum()

(b - a) / b<jupyter_output><empty_output><jupyter_text>Partial answerThe box plots show that tortoiseshell cats tend to get fewer pats than white cats which get fewer than any other colour. The coat colour is an important predictor to include, it accounts for about 20% of the residual sum of squares without coat colour  Among short-haired breeds they say longer hair is better, among long-haired breeds they say short hair is better, who is correct?Since in the figures above we saw that the breed longhaired/shorthaired appears to separate the data, it may be useful to consider different models on each subset. Fit a linear model to each subset of the data and see that the effect of the coat length is in each case.<jupyter_code>lm_4a = smf.ols("num_pats ~ time_outdoors + C(coat_colour) + weight + height + loudness + coat_length", cat_df[cat_df["is_longhaired"] == 1]).fit()
lm_4b = smf.ols("num_pats ~ time_outdoors + C(coat_colour) + weight + height + loudness + coat_length", cat_df[cat_df["is_longhaired"] == 0]).fit()
print(lm_4a.summary())<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:               num_pats   R-squared:                       0.444
Model:                            OLS   Adj. R-squared:                  0.436
Method:                 Least Squares   F-statistic:                     55.53
Date:                Thu, 12 Nov 2020   Prob (F-statistic):           3.57e-58
Time:                        15:16:59   Log-Likelihood:                -1563.8
No. Observations:                 494   AIC:                             3144.
Df Residuals:                     486   BIC:                             3177.
Df Model:                           7                                         
Covariance Type:            nonrobust                                         
=======================================================================================
                          coef    std err  [...]<jupyter_text>Partial answerWhen we fit the model to each subset of the data we see that the coefficient changes direction suggesting there is an interaction here.  Fit a model with an interacion term between the coat length and the long/shorthaired breedWhat does this tell us about the age old debate about cat hair length?<jupyter_code>lm_4c = smf.ols("num_pats ~ time_outdoors + C(coat_colour) + weight + height + loudness + C(is_longhaired) : coat_length", cat_df).fit()
print(lm_4c.summary())<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:               num_pats   R-squared:                       0.579
Model:                            OLS   Adj. R-squared:                  0.576
Method:                 Least Squares   F-statistic:                     170.5
Date:                Thu, 12 Nov 2020   Prob (F-statistic):          2.10e-180
Time:                        15:16:59   Log-Likelihood:                -3230.3
No. Observations:                1000   AIC:                             6479.
Df Residuals:                     991   BIC:                             6523.
Df Model:                           8                                         
Covariance Type:            nonrobust                                         
===================================================================================================
                               [...]<jupyter_text>the principle of marginality tells us in general we should include the main effects... but here we showed that we can remove the main effects and just include the interaction termi.e. use the * and not the :  Partial answerFitting the model with an interaction term suggests that both are correct: if you are a shorthaired breed, increasing your coat lenght tends to increase the number of pats and if you are a longhaired breed it will decrease the number of pats on average.  How else could we handle coat length?We could instead have included quadratic terms for coat length to see if this was a better way to explain the non-linear effect.  Shouldn't we check for influential points?We can generate a plot of the studentized residuals and the leverage to check if there are any influential points.If there is a potential outlier, does removing it change anything?<jupyter_code>sm.graphics.influence_plot(lm_4c, criterion="cooks")
cat_df[630:635]
cat_df.describe()
sm.graphics.plot_fit(lm_4c, "time_outdoors")
# Warning!
# ========
# We have to be very careful removing data or only work with copies of the
# dataframe.
cat_df_2 = cat_df.copy()
print(cat_df_2[632:635])
cat_df_2.drop([633], axis=0, inplace=True)
print(80 *"=")
print(cat_df_2[633:635])
lm_5 = smf.ols("num_pats ~ time_outdoors + C(coat_colour) + weight + height + loudness + C(is_longhaired) : coat_length", cat_df_2).fit()
print(lm_5.summary())
print(lm_4c.summary())<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:               num_pats   R-squared:                       0.574
Model:                            OLS   Adj. R-squared:                  0.570
Method:                 Least Squares   F-statistic:                     166.4
Date:                Thu, 12 Nov 2020   Prob (F-statistic):          2.43e-177
Time:                        15:17:01   Log-Likelihood:                -3223.8
No. Observations:                 999   AIC:                             6466.
Df Residuals:                     990   BIC:                             6510.
Df Model:                           8                                         
Covariance Type:            nonrobust                                         
===================================================================================================
                               [...]
{'fraction_non_alphanumeric': 0.1122068550811786, 'fraction_numerical': 0.019542994588093807, 'mean_word_length': 2.8803079794680353, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 35, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '27648', 'license': 'MIT', 'n_tokens_mistral': 4732, 'n_tokens_neox': 4307, 'n_words': 1864, 'path': 'example-2/example-2-answers.ipynb', 'repo_name': 'xw-song/aas-extended-examples', 'size': '16630'}