<jupyter_start><jupyter_text>Chapter 4 Linear regression: before and after fitting the model<jupyter_code>import warnings

import pandas as pd
import proplot as plot
import seaborn as sns
import statsmodels.api as sm
import statsmodels.formula.api as smf
from scipy import stats

warnings.filterwarnings("ignore")
%pylab inline


plt.rcParams["axes.labelweight"] = "bold"
plt.rcParams["font.weight"] = "bold"
earnings_df = pd.read_csv("../data/earnings/heights.tsv.gz", sep="\t")
earnings_df.head()
earnings_df.shape
earnings_df.dropna().shape
earnings_df = earnings_df.dropna()
earnings_df.shape
earnings_df.head()<jupyter_output><empty_output><jupyter_text>Centering and standardizing, especially for models with interactions<jupyter_code>kidiq_df = pd.read_csv("../data/kidiq.tsv.gz", sep="\t")
kidiq_df.head()
model = smf.ols(
    formula="""kid_score ~ mom_hs + mom_iq + mom_hs:mom_iq""", data=kidiq_df
)
results_hs = model.fit()
print(results_hs.summary())
print("Residual SD: {}".format(np.round(results_hs.resid.std(), 2)))<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:              kid_score   R-squared:                       0.230
Model:                            OLS   Adj. R-squared:                  0.225
Method:                 Least Squares   F-statistic:                     42.84
Date:                Sat, 20 Jun 2020   Prob (F-statistic):           3.07e-24
Time:                        23:35:23   Log-Likelihood:                -1867.5
No. Observations:                 434   AIC:                             3743.
Df Residuals:                     430   BIC:                             3759.
Df Model:                           3                                         
Covariance Type:            nonrobust                                         
=================================================================================
                    coef    std err          t   [...]<jupyter_text>Centering by subtracting the mean<jupyter_code>kidiq_df["mom_hs_centered"] = kidiq_df["mom_hs"] - kidiq_df["mom_hs"].mean()
kidiq_df["mom_iq_centered"] = kidiq_df["mom_iq"] - kidiq_df["mom_iq"].mean()


model = smf.ols(
    formula="""kid_score ~ mom_hs_centered + mom_iq_centered + mom_hs_centered:mom_iq_centered""",
    data=kidiq_df,
)
results_hs = model.fit()
print(results_hs.summary())
print("Residual SD: {}".format(np.round(results_hs.resid.std(), 2)))<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:              kid_score   R-squared:                       0.230
Model:                            OLS   Adj. R-squared:                  0.225
Method:                 Least Squares   F-statistic:                     42.84
Date:                Sat, 20 Jun 2020   Prob (F-statistic):           3.07e-24
Time:                        23:35:23   Log-Likelihood:                -1867.5
No. Observations:                 434   AIC:                             3743.
Df Residuals:                     430   BIC:                             3759.
Df Model:                           3                                         
Covariance Type:            nonrobust                                         
===================================================================================================
                               [...]<jupyter_text>The residual SD and $R^2$ do not change as a linear transformation does not affect the fit. Interpretation:        - Intercept (87): Score of children with  mothers who have an average IQ and have gone to a "mean high school"     - Coeff of mom_hs (2.8): Score of children differes by 2.8 for children with mother's mean high school difference of 1 (i.e. went to school vs didn't) as long as they have the same mean IQ.    - Coeff of mom_iq (0.6): Score of children will differ by 0.6 for a 1 point difference in mother's IQ if they are at the same average level of "mean high school"    - Interaction effect: Difference in slope for mom_iq comparing children with mothers who did and did not complete     Coefficient of each main effect corresponds to the predictive difference with the other input at its average value.      Using a conventional centering pointIf we center on an understandable reference - the above values are more interpretable.<jupyter_code>kidiq_df["mom_hs_centered2"] = kidiq_df["mom_hs"] - 0.5
kidiq_df["mom_iq_centered2"] = kidiq_df["mom_iq"] - 100


model = smf.ols(
    formula="""kid_score ~ mom_hs_centered2 + mom_iq_centered2 + mom_hs_centered2:mom_iq_centered2""",
    data=kidiq_df,
)
results_hs = model.fit()
print(results_hs.summary())
print("Residual SD: {}".format(np.round(results_hs.resid.std(), 2)))<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:              kid_score   R-squared:                       0.230
Model:                            OLS   Adj. R-squared:                  0.225
Method:                 Least Squares   F-statistic:                     42.84
Date:                Sat, 20 Jun 2020   Prob (F-statistic):           3.07e-24
Time:                        23:35:23   Log-Likelihood:                -1867.5
No. Observations:                 434   AIC:                             3743.
Df Residuals:                     430   BIC:                             3759.
Df Model:                           3                                         
Covariance Type:            nonrobust                                         
=====================================================================================================
                             [...]<jupyter_text>- Coefficient of mom_hs_centered2 is the average predictive difference between a child with mom.hs = 1 and mom.hs = 0 for those children with mom.iq = 100.- Coefficient of mom_iq_centered2 is the average predictive difference between a child with mom.iq = 0.5  Standardizing by subtracting the mean and dividing by 2 standard deviations<jupyter_code>kidiq_df["mom_hs_centeredZ"] = (
    (kidiq_df["mom_hs"] - kidiq_df["mom_hs"].mean()) / 2 * kidiq_df["mom_hs"].std()
)
kidiq_df["mom_iq_centeredZ"] = (
    (kidiq_df["mom_iq"] - kidiq_df["mom_iq"].mean()) / 2 * kidiq_df["mom_iq"].std()
)


model = smf.ols(
    formula="""kid_score ~ mom_hs_centeredZ + mom_iq_centeredZ + mom_hs_centeredZ:mom_iq_centeredZ""",
    data=kidiq_df,
)
results_hs = model.fit()
print(results_hs.summary())
print("Residual SD: {}".format(np.round(results_hs.resid.std(), 2)))<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:              kid_score   R-squared:                       0.230
Model:                            OLS   Adj. R-squared:                  0.225
Method:                 Least Squares   F-statistic:                     42.84
Date:                Sat, 20 Jun 2020   Prob (F-statistic):           3.07e-24
Time:                        23:35:23   Log-Likelihood:                -1867.5
No. Observations:                 434   AIC:                             3743.
Df Residuals:                     430   BIC:                             3759.
Df Model:                           3                                         
Covariance Type:            nonrobust                                         
=====================================================================================================
                             [...]<jupyter_text>Why scale by 2 standard deviations?If $x$ was a binary variable them its standard deviation is given by $\sqrt{0.5 \times 0.5}= 0.5 $.If we do not divide by 2, the difference between the extremities of the binary variable  $(0.5/0.5) - (-0.5/0.5) = 2$ and hence the coefficient would correspond half the difference between the two possible values of x. While dividing by 2 will result in a difference of 1.  Logarithmic transformationsA linear model on the logarithmic scale corresponds to a multiplicative model on the original scale.  Height and earnings example<jupyter_code>earnings_df["log_earn"] = np.log(earnings_df["earn"])
earnings_df = earnings_df.replace([np.inf, -np.inf], np.nan)

earnings_df = earnings_df.dropna()
earn_log_model = smf.ols(formula="""log_earn ~ height""", data=earnings_df).fit()
print(earn_log_model.summary())<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:               log_earn   R-squared:                       0.060
Model:                            OLS   Adj. R-squared:                  0.060
Method:                 Least Squares   F-statistic:                     76.44
Date:                Sat, 20 Jun 2020   Prob (F-statistic):           7.62e-18
Time:                        23:35:23   Log-Likelihood:                -1555.6
No. Observations:                1192   AIC:                             3115.
Df Residuals:                    1190   BIC:                             3125.
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|[...]<jupyter_text>InterpretationA coefficient of 0.06 on log scale for height implies a unit change in height will lead to $\exp{0.06}= 1.06$ difference in earnings. Thus a change of 1 unit in height leads to a change of 6% in earnings. Why use a natural log instead of $log_{10}$Coefficients on the natural-log scale are "more" interpretable. A coefficient of $0.06$ for example corresponds toa change of 6%. On the other hand, a change on $\log_{10}$ scale is harder to interpret:<jupyter_code>earnings_df["log10_earn"] = np.log10(earnings_df["earn"])
earnings_df = earnings_df.replace([np.inf, -np.inf], np.nan)

earnings_df = earnings_df.dropna()
earn_log_model = smf.ols(formula="""log10_earn ~ height""", data=earnings_df).fit()
print(earn_log_model.summary())
print("Residual SD: {}".format(np.round(earn_log_model.resid.std(), 2)))<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:             log10_earn   R-squared:                       0.060
Model:                            OLS   Adj. R-squared:                  0.060
Method:                 Least Squares   F-statistic:                     76.44
Date:                Sat, 20 Jun 2020   Prob (F-statistic):           7.62e-18
Time:                        23:35:23   Log-Likelihood:                -561.42
No. Observations:                1192   AIC:                             1127.
Df Residuals:                    1190   BIC:                             1137.
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|[...]<jupyter_text>A coefficient of $0.0255$ is less interpretable here as it still corresponds to a change of $10^{0.0255} = 1.062$   Building a regression model on the log scale<jupyter_code>earnings_df["male"] = 1
# Set rows where sex==2 to female (0)
earnings_df.loc[earnings_df.sex == 2, "male"] = 0
earn_log_model = smf.ols(formula="""log_earn ~ height + male""", data=earnings_df).fit()
print(earn_log_model.summary())
print("Residual SD: {}".format(np.round(earn_log_model.resid.std(), 2)))<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:               log_earn   R-squared:                       0.087
Model:                            OLS   Adj. R-squared:                  0.085
Method:                 Least Squares   F-statistic:                     56.34
Date:                Sat, 20 Jun 2020   Prob (F-statistic):           4.21e-24
Time:                        23:35:23   Log-Likelihood:                -1538.7
No. Observations:                1192   AIC:                             3083.
Df Residuals:                    1189   BIC:                             3099.
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|[...]<jupyter_text>Including an interaction<jupyter_code>earnings_df["male"] = 1
# Set rows where sex==2 to female (0)
earnings_df.loc[earnings_df.sex == 2, "male"] = 0
earn_log_model3 = smf.ols(
    formula="""log_earn ~ height + male + height:male""", data=earnings_df
).fit()
print(earn_log_model3.summary())
print("Residual SD: {}".format(np.round(earn_log_model3.resid.std(), 2)))<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:               log_earn   R-squared:                       0.087
Model:                            OLS   Adj. R-squared:                  0.084
Method:                 Least Squares   F-statistic:                     37.58
Date:                Sat, 20 Jun 2020   Prob (F-statistic):           3.32e-23
Time:                        23:35:23   Log-Likelihood:                -1538.6
No. Observations:                1192   AIC:                             3085.
Df Residuals:                    1188   BIC:                             3106.
Df Model:                           3                                         
Covariance Type:            nonrobust                                         
===============================================================================
                  coef    std err          t      P[...]<jupyter_text>Interpretation:        - Intercept = Predicted log earning when height and male is zero (which is non feasible hence not interpretable)    - Coeff of height (0.017) = Predicted difference in log earnings corresponding to a 1-inch difference in height if male = 0 => The estimated predictive difference in earnings for a 1-inch difference in height is 1.7% for women.    - Coeff of male (-0.0786): Predicted difference in log earnings between women and men if height equals 0. But heights are never zero so this is not interpretable.    - Interaction term height:male (0.0074): Predictive difference in log earnings comparing men to women => An inch of height difference corresponds to 0.7% more of an in increase in earnings among men than among women and the estimated predictive difference in log earnings is 1.7% + 0.7% = 2.4% per inch of height in men.          Linear transformation to make coefficients more interpretable:<jupyter_code>earnings_df["height_z"] = (
    earnings_df["height"] - earnings_df["height"].mean()
) / earnings_df["height"].std()
earn_log_model4 = smf.ols(
    formula="""log_earn ~ height_z + male + height_z:male""", data=earnings_df
).fit()
print(earn_log_model4.summary())
print("Residual SD: {}".format(np.round(earn_log_model4.resid.std(), 2)))
print(earnings_df["height"].mean(), earnings_df["height"].std())<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:               log_earn   R-squared:                       0.087
Model:                            OLS   Adj. R-squared:                  0.084
Method:                 Least Squares   F-statistic:                     37.58
Date:                Sat, 20 Jun 2020   Prob (F-statistic):           3.32e-23
Time:                        23:35:23   Log-Likelihood:                -1538.6
No. Observations:                1192   AIC:                             3085.
Df Residuals:                    1188   BIC:                             3106.
Df Model:                           3                                         
Covariance Type:            nonrobust                                         
=================================================================================
                    coef    std err          t   [...]<jupyter_text>Interpretation:        - Intercept = Predicted log earning when height is average height and male is zero (which is non feasible hence not interpretable) => 66.9 inch women has a log earnings of 9.52    - Coeff of height (0.065) = Predicted difference in log earnings corresponding to a 1-inch standard deviation in height if male = 0 => The estimated predictive difference in earnings for a 3.8 inch difference in height of females is 6% for women.    - Coeff of male (0.4197): Predicted difference in log earnings between women and men if height equals 66.9 inches. Thus, males of height 66.9 inches are expected to have 41.9% higher log earning than women    - Interaction term height:male (0.0286): Predictive difference in log earnings comparing men to women => A 2.8 inch of height difference corresponds to 2.8% more of an in increase in earnings among men than among women and the estimated predictive difference in log earnings is 41.9% + 2.7% = 44.6% per inch of height in men.          Using discrete rather than continous predictors<jupyter_code>kid_score_momwork = smf.ols(formula="""kid_score ~ C(mom_work)""", data=kidiq_df).fit()
print(kid_score_momwork.summary())
print("Residual SD: {}".format(np.round(kid_score_momwork.resid.std(), 2)))<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:              kid_score   R-squared:                       0.024
Model:                            OLS   Adj. R-squared:                  0.018
Method:                 Least Squares   F-statistic:                     3.590
Date:                Sat, 20 Jun 2020   Prob (F-statistic):             0.0138
Time:                        23:35:24   Log-Likelihood:                -1918.9
No. Observations:                 434   AIC:                             3846.
Df Residuals:                     430   BIC:                             3862.
Df Model:                           3                                         
Covariance Type:            nonrobust                                         
====================================================================================
                       coef    std err        [...]<jupyter_text>Building regression models for prediction<jupyter_code>mesquite_df = pd.read_csv("../data/mesquite/mesquite.tsv.gz", sep="\t")
mesquite_df.head()
formula = """LeafWt ~ Diam1 + Diam2 + CanHt + TotHt + Dens + Group"""
mesquite_model = smf.ols(formula=formula, data=mesquite_df).fit()
print(mesquite_model.summary())
print("Residual SD: {}".format(np.round(mesquite_model.resid.std(), 2)))<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:                 LeafWt   R-squared:                       0.848
Model:                            OLS   Adj. R-squared:                  0.825
Method:                 Least Squares   F-statistic:                     36.34
Date:                Sat, 20 Jun 2020   Prob (F-statistic):           1.73e-14
Time:                        23:35:24   Log-Likelihood:                -318.82
No. Observations:                  46   AIC:                             651.6
Df Residuals:                      39   BIC:                             664.4
Df Model:                           6                                         
Covariance Type:            nonrobust                                         
================================================================================
                   coef    std err          t     [...]<jupyter_text>Multiplicative Model<jupyter_code>formula = """log(LeafWt) ~ log(Diam1) + log(Diam2) + log(CanHt) + log(TotHt) + log(Dens) + Group"""
mesquite_model = smf.ols(formula=formula, data=mesquite_df).fit()
print(mesquite_model.summary())
print("Residual SD: {}".format(np.round(mesquite_model.resid.std(), 2)))<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:            log(LeafWt)   R-squared:                       0.887
Model:                            OLS   Adj. R-squared:                  0.870
Method:                 Least Squares   F-statistic:                     51.17
Date:                Sat, 20 Jun 2020   Prob (F-statistic):           5.73e-17
Time:                        23:35:24   Log-Likelihood:                -10.406
No. Observations:                  46   AIC:                             34.81
Df Residuals:                      39   BIC:                             47.61
Df Model:                           6                                         
Covariance Type:            nonrobust                                         
================================================================================
                   coef    std err          t     [...]<jupyter_text>Interpretation => A difference of x% in canopy height leads to a difference of 0.37x% difference in leaf weight.The above model is a complicated one. Maybe we can start with something very simple.<jupyter_code>mesquite_df["CanVol"] = mesquite_df['Diam1'] * mesquite_df['Diam2'] * mesquite_df['CanHt'] 
formula = """log(LeafWt) ~ log(CanVol)"""
mesquite_model = smf.ols(formula=formula, data=mesquite_df).fit()
print(mesquite_model.summary())
print("Residual SD: {}".format(np.round(mesquite_model.resid.std(), 2)))<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:            log(LeafWt)   R-squared:                       0.799
Model:                            OLS   Adj. R-squared:                  0.795
Method:                 Least Squares   F-statistic:                     175.1
Date:                Sat, 20 Jun 2020   Prob (F-statistic):           6.09e-17
Time:                        23:35:24   Log-Likelihood:                -23.687
No. Observations:                  46   AIC:                             51.37
Df Residuals:                      44   BIC:                             55.03
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
===============================================================================
                  coef    std err          t      P[...]<jupyter_text>So leaf weight is approximately proportional to 0.72th power to Canopy Volume. It has a impressive $R^2$ of 0.8<jupyter_code>mesquite_df["CanVol"] = mesquite_df['Diam1'] * mesquite_df['Diam2'] * mesquite_df['CanHt'] 
mesquite_df["CanArea"] = mesquite_df['Diam1'] * mesquite_df['Diam2']
mesquite_df["Shape"] = mesquite_df['Diam1'] / mesquite_df['Diam2']

formula = """log(LeafWt) ~ log(CanVol) + log(CanArea) + log(Shape) + log(TotHt) + log(Dens) + Group"""

mesquite_model = smf.ols(formula=formula, data=mesquite_df).fit()
print(mesquite_model.summary())
print("Residual SD: {}".format(np.round(mesquite_model.resid.std(), 2)))
formula = """log(LeafWt) ~ log(CanVol) + log(CanArea) + Group"""

mesquite_model = smf.ols(formula=formula, data=mesquite_df).fit()
print(mesquite_model.summary())
print("Residual SD: {}".format(np.round(mesquite_model.resid.std(), 2)))<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:            log(LeafWt)   R-squared:                       0.873
Model:                            OLS   Adj. R-squared:                  0.864
Method:                 Least Squares   F-statistic:                     96.01
Date:                Sat, 20 Jun 2020   Prob (F-statistic):           7.79e-19
Time:                        23:35:24   Log-Likelihood:                -13.198
No. Observations:                  46   AIC:                             34.40
Df Residuals:                      42   BIC:                             41.71
Df Model:                           3                                         
Covariance Type:            nonrobust                                         
================================================================================
                   coef    std err          t     [...]
{'fraction_non_alphanumeric': 0.17266677357280308, 'fraction_numerical': 0.036561898652982684, 'mean_word_length': 1.7211737754990728, 'pattern_counts': {'":': 0, '<': 47, '<?xml version=': 0, '>': 56, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 50, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '430907', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 7307, 'n_tokens_neox': 6306, 'n_words': 1760, 'path': 'Chapters/04_Chapter04_Linear_regresion_before_and_after_fitting_the_model.ipynb', 'repo_name': 'saketkc/ARMpy', 'size': '24944'}