<jupyter_start><jupyter_text>Multiple Linear RegressionAnalysis of basic UK Climate Data using Multiple Regression (with [statsmodels](https://www.statsmodels.org/stable/index.html)) SetupLoad required packages<jupyter_code>import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt
import seaborn as sns<jupyter_output><empty_output><jupyter_text>The data are for 20 weather stations across the UK. Read the data file into a Pandas DataFrame and check<jupyter_code>df = pd.read_csv('../../data/climate/UKJanClimate.csv', header=0)
df.info()
df.head(10)<jupyter_output><empty_output><jupyter_text>Checking CollinearityIndependent variables must not be collinear (i.e. correlated amongst themselves). Generally variables need to be r  -0.5Check using a correlation matrix on our dataframe<jupyter_code>df.corr().round(3)
corrmat = df.corr()  #save as an object
sns.heatmap(corrmat)  #then pass to heatmap
#with our own palette
sns.heatmap(df.corr(method='spearman'), cmap=sns.diverging_palette(20, 220, n=200))<jupyter_output><empty_output><jupyter_text>Fitting & Selecting Models  Fitting For a single variable we could do:<jupyter_code>#regplot adds a regression line https://seaborn.pydata.org/tutorial/regression.html
sns.regplot(x="Lat",
           y="min_temp",
           data=df)
plt.show()
#require import statsmodels.formula.api as smf
lat_mod = smf.ols("min_temp ~ Lat", data=df).fit()
print (lat_mod.summary())<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:               min_temp   R-squared:                       0.414
Model:                            OLS   Adj. R-squared:                  0.382
Method:                 Least Squares   F-statistic:                     12.72
Date:                Fri, 25 Feb 2022   Prob (F-statistic):            0.00220
Time:                        17:58:25   Log-Likelihood:                -23.402
No. Observations:                  20   AIC:                             50.80
Df Residuals:                      18   BIC:                             52.79
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|[...]<jupyter_text>For multiple predictor (independent) variables:<jupyter_code>#require import statsmodels.formula.api as smf
lat_elev_mod = smf.ols("min_temp ~ Lat + Elev", data=df).fit()
print (lat_elev_mod.summary())
#require import statsmodels.formula.api as smf
lat_elev_lon_mod = smf.ols("min_temp ~ Lat + Elev + Long", data=df).fit()
print (lat_elev_lon_mod.summary())<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:               min_temp   R-squared:                       0.716
Model:                            OLS   Adj. R-squared:                  0.662
Method:                 Least Squares   F-statistic:                     13.43
Date:                Fri, 25 Feb 2022   Prob (F-statistic):           0.000123
Time:                        17:58:25   Log-Likelihood:                -16.171
No. Observations:                  20   AIC:                             40.34
Df Residuals:                      16   BIC:                             44.32
Df Model:                           3                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|[...]<jupyter_text>Model SelectionComparing the last few models:- which has the best $r^2$?- which independent variables are statistically significant at 95% confidence? - which model has the smallest AIC?- which models have differences in AIC (from the min) >2?- how do the Adjusted $r^2$ of the models look?The three-variable model has the largest $r^2$ (and Adjusted $r^2$) so explains most overall variation. However, the `Long` variable is not statistically significant (at 95% confidence). And the difference in AIC between two- and three-variables models is <2.   So... - For explanation, the two-variable model would be best as it is the most parsimonious and all independent variables are stat. sig. at 95%- For prediction, use three-variable model as it explains the most overall variance (i.e. is the most predictive)  Model Diagnostics As we created models above we really should have check model diagnostics to check assumptions of regression had not been violated. Let's just do it for the two-variable model here.  Normality of Residuals<jupyter_code>#n=20 so important to check normality!
plt.hist(lat_elev_mod.resid_pearson, 
         edgecolor='black',
         bins=[-2,-1.5,-1,-0.5,0,0.5,1,1.5,2])
plt.xlabel('Residuals')
#better to make our own plot using standardised residuals (resid_pearson from the model object)
plt.plot(lat_elev_mod.fittedvalues, lat_elev_mod.resid_pearson,
         'oC0',  
         markeredgecolor='black')
plt.axhline(c = 'red')
plt.xlabel('Predicted y Values')
plt.ylabel('Standardized Residuals (z)')  
#lots of other ways to check (see https://www.statsmodels.org/dev/generated/statsmodels.graphics.gofplots.qqplot.html)
qqfig = sm.qqplot(lat_elev_mod.resid_pearson, line='r')<jupyter_output>/home/james/anaconda3/envs/sds2021/lib/python3.8/site-packages/statsmodels/graphics/gofplots.py:993: UserWarning: marker is redundantly defined by the 'marker' keyword argument and the fmt string "bo" (-> marker='o'). The keyword argument will take precedence.
  ax.plot(x, y, fmt, **plot_style)<jupyter_text>Leverage and Influence*Leverage* is an indicator of how extreme our independent variable (i.e. elevation, latitude) values are for each observation (i.e. each location). It's important to identify observations with high leverage because in certain situations they may have an undue *influence* on the estimated regression function (and parameters). Points with undue influence should be considered for removal from the data used for a regression model (because there is something about them that means that are 'different' from the remaining data).Note, as [<NAME>](https://www.real-statistics.com/multiple-regression/outliers-and-influencers/) points out:> Keep in mind that since we are dealing with a multi-dimensional model, there may be data points that look perfectly fine in any single dimension but are multivariate outliers. E.g. for the general population, there is nothing unusual about a 6-foot man or a 125-pound man, but a 6-foot man that weighs 125 pounds is unusual.There's a nice intuitive blog on [Leverage and Influence in a Nutshell](https://www.r-bloggers.com/2016/06/leverage-and-influence-in-a-nutshell/) with an interactive plot that might help to understand the different issues discussed here (also see the other sources linked to here, and readings).    The Influence PlotWe can use [the `influence_plot` from statsmodels](https://www.statsmodels.org/dev/generated/statsmodels.graphics.regressionplots.influence_plot.html) to visualise leverage and influence:<jupyter_code>fig, ax = plt.subplots(figsize=(10,6))
fig = sm.graphics.influence_plot(lat_elev_mod, alpha=0.05, 
                                 criterion="cooks", ax=ax)<jupyter_output><empty_output><jupyter_text>The y-axis on this plot is the 'studentized' residuals - this is the same as 'standardized' residuals (z). This axis helps to identify outliers:> In linear regression, an outlier is an observation with large residual. In other words, it is an observation whose dependent-variable value is unusual given its values on the predictor variables. An outlier may indicate a sample peculiarity or may indicate a data entry error or other problem. [[source](https://songhuiming.github.io/pages/2016/11/27/linear-regression-in-python-outliers-leverage-detect/)]The x-axis on this plot is the measure of leverage. Leverage is:>  a measure of the distance between the $x$ value for the data point [i.e., observation] and the mean of the $x$ values for all $n$ data points. The leverage is a number between 0 and 1, inclusive. The sum of the equals $p$, the number of parameters (regression coefficients including the intercept). [[source](https://online.stat.psu.edu/stat501/lesson/11/11.2)]High leverage points [can be assumed](https://online.stat.psu.edu/stat501/lesson/11/11.2) to be those with leverage more than 3 times larger than the mean leverage.The size of the point indicates the value of [Cook's Distance](https://www.statsmodels.org/stable/generated/statsmodels.graphics.regressionplots.influence_plot.html) (we specified that with the `criterion` argument to `influence_plot`).> If we remove a point from the sample, then the equation for regression line changes. Points that have the most influence produce the largest change in the equation of the regression line. A measure of this influence is called *Cook’s Distance*. [[source](https://www.real-statistics.com/multiple-regression/outliers-and-influencers/)]Different people argue for different thresholds of Cook's Distance ($D$) greater than which a point should be considered as having undue influence. Some argue ($D=1$) is an appropriate threshold (as Cook himself did), while [others argue](https://stats.stackexchange.com/questions/87962/cooks-distance-cut-off-value) the threshold should be $4/n$, where $n$ is the number of data points. A problem with the statsmodels `influence_plot` is that it is difficult to know what value of $D$ a point has from the size of the point (especially as no legend is provided; we'll see how to output the actual $D$ values below).  Note that three points are labelled on the plot. Points that are labelled are:- outliers on the y-axis. - high leverage points on the x-axis. The definition of an outlier for statsmodels is given by the `alpha` argument. Here we used `0.05` so an outlier is any point with a standardized residual 1.96 [[source](https://songhuiming.github.io/pages/2016/11/27/linear-regression-in-python-outliers-leverage-detect/)]. Statsmodels defines a high leverage point as $2\frac{(df + 1)}{n}$ where $df$ is the degrees of freedom of the model.  However, just because a point is indicated in the `influence_plot` with a label, does not mean these points actually have *undue* influence. Ultimately, because influence can be thought of as the product of leverage and outlierness, the points we want to be most worried about are:- those with high $D$ (i.e the largest points)- those points with both high leverage and large (absolute) residuals (i.e. those in the top right or bottom right corners of the plot).  As noted above it is difficult to assess from the point size what the value of Cook's Distance $D$ is. We can access the values directly as follows:<jupyter_code>lat_elev_mod.get_influence().cooks_distance[0]<jupyter_output><empty_output><jupyter_text>From this we can see we have no points with $D>1$ but using the $D > 4/n$ (0.2) we find that two points could be considered to have undue influence. However, > It’s important to note that Cook’s Distance is often used as a way to identify influential data points. Just because a data point is influential doesn’t mean it should necessarily be deleted – first you should check to see if the data point has simply been incorrectly recorded or if there is something strange about the data point that may point to an interesting finding. [[source](https://www.statology.org/how-to-identify-influential-data-points-using-cooks-distance/)]In this case, points 8 and 9 are ones we should consider investigating further (8 due to highest $D$ value and 9 because it's an outlier with some leverage, but not too extreme) and all other points are likely fine.   Partial Regression PlotsAnother way to visualise leverage and influence that some people use is [partial regression plots](https://en.wikipedia.org/wiki/Partial_regression_plot). The partial regression plot is useful to detect outliers and check that the relationships between independent and dependent variables make sense. In univariate regression identifying these things would be straight-forward with a scatter plot of the dependent variable against the (single) independent variable. However, it is not so striaght-forward to do this for multivariate regression models. Remember,>The estimates of the beta values measure the expected change in the Y value, given a one-unit change in the X variable, *holding all other X variables constant*. (Harris 2016, p.197, emphasis added)Simply regressing multiple independent variables against the dependent variable do not visualise the results of the multiple regression model, because the other independent variables are not being held constant (and as we saw above, all independent variables are correlated to a degree). However, partial regression plots *are* able to to control for this variation. A partial regression plot is a scatter plot of residuals from regression models of:- the dependent variable against all the independent variables except the one of interest- the independent variable of interest against all the other independent variablesThis means we get a scatter plot for *every* independent variable (as well as the intercept). We can use [the `plot_partregress_grid` function](https://www.statsmodels.org/dev/generated/statsmodels.graphics.regressionplots.plot_partregress_grid.html) in statsmodels to produce a partial regression plot for our `lat_elev_mod` fit model object:<jupyter_code>fig = plt.figure(figsize=(12,8))
fig = sm.graphics.plot_partregress_grid(lat_elev_mod, fig=fig)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.098705966930266, 'fraction_numerical': 0.017253774263120056, 'mean_word_length': 3.5926048200726313, 'pattern_counts': {'":': 0, '<': 36, '<?xml version=': 0, '>': 47, 'https://': 16, 'lorem ipsum': 0, 'www.': 9, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '231251', 'license': 'MIT', 'n_tokens_mistral': 3932, 'n_tokens_neox': 3583, 'n_words': 1534, 'path': 'code/python/5SSG2060_MultiRegression-Climate.ipynb', 'repo_name': 'jamesdamillington/sds-principles-part2-2022', 'size': '13919'}