Chapter 4 Linear regression: before and after fitting the model

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"
Populating the interactive namespace from numpy and matplotlib
earnings_df = pd.read_csv("../data/earnings/heights.tsv.gz", sep="\t")
earnings_df.head()
earn height1 height2 sex race hisp ed yearbn height
0 NaN 5.0 6.0 2 1 2 12 53 66.0
1 NaN 5.0 4.0 1 2 2 12 50 64.0
2 50000.0 6.0 2.0 1 1 2 16 45 74.0
3 60000.0 5.0 6.0 2 1 2 16 32 66.0
4 30000.0 5.0 4.0 2 1 2 16 61 64.0
earnings_df.shape
(2029, 9)
earnings_df.dropna().shape
(1379, 9)
earnings_df = earnings_df.dropna()
earnings_df.shape
(1379, 9)
earnings_df.head()
earn height1 height2 sex race hisp ed yearbn height
2 50000.0 6.0 2.0 1 1 2 16 45 74.0
3 60000.0 5.0 6.0 2 1 2 16 32 66.0
4 30000.0 5.0 4.0 2 1 2 16 61 64.0
6 50000.0 5.0 3.0 2 3 2 16 99 63.0
8 51000.0 5.0 3.0 2 1 2 17 51 63.0

Centering and standardizing, especially for models with interactions

kidiq_df = pd.read_csv("../data/kidiq.tsv.gz", sep="\t")
kidiq_df.head()
kid_score mom_hs mom_iq mom_work mom_age
0 65 1.0 121.117529 4 27
1 98 1.0 89.361882 4 25
2 85 1.0 115.443165 4 27
3 83 1.0 99.449639 3 25
4 115 1.0 92.745710 4 27
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)))
                            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      P>|t|      [0.025      0.975]
---------------------------------------------------------------------------------
Intercept       -11.4820     13.758     -0.835      0.404     -38.523      15.559
mom_hs           51.2682     15.338      3.343      0.001      21.122      81.414
mom_iq            0.9689      0.148      6.531      0.000       0.677       1.260
mom_hs:mom_iq    -0.4843      0.162     -2.985      0.003      -0.803      -0.165
==============================================================================
Omnibus:                        8.014   Durbin-Watson:                   1.660
Prob(Omnibus):                  0.018   Jarque-Bera (JB):                8.258
Skew:                          -0.333   Prob(JB):                       0.0161
Kurtosis:                       2.887   Cond. No.                     3.10e+03
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 3.1e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
Residual SD: 17.91

Centering by subtracting the mean

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)))
                            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      P>|t|      [0.025      0.975]
---------------------------------------------------------------------------------------------------
Intercept                          87.6389      0.908     96.565      0.000      85.855      89.423
mom_hs_centered                     2.8408      2.427      1.171      0.242      -1.929       7.610
mom_iq_centered                     0.5884      0.061      9.712      0.000       0.469       0.707
mom_hs_centered:mom_iq_centered    -0.4843      0.162     -2.985      0.003      -0.803      -0.165
==============================================================================
Omnibus:                        8.014   Durbin-Watson:                   1.660
Prob(Omnibus):                  0.018   Jarque-Bera (JB):                8.258
Skew:                          -0.333   Prob(JB):                       0.0161
Kurtosis:                       2.887   Cond. No.                         42.2
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
Residual SD: 17.91

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 point

If we center on an understandable reference - the above values are more interpretable.

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)))
                            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      P>|t|      [0.025      0.975]
-----------------------------------------------------------------------------------------------------
Intercept                            86.8273      1.213     71.561      0.000      84.442      89.212
mom_hs_centered2                      2.8408      2.427      1.171      0.242      -1.929       7.610
mom_iq_centered2                      0.7268      0.081      8.960      0.000       0.567       0.886
mom_hs_centered2:mom_iq_centered2    -0.4843      0.162     -2.985      0.003      -0.803      -0.165
==============================================================================
Omnibus:                        8.014   Durbin-Watson:                   1.660
Prob(Omnibus):                  0.018   Jarque-Bera (JB):                8.258
Skew:                          -0.333   Prob(JB):                       0.0161
Kurtosis:                       2.887   Cond. No.                         46.9
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
Residual SD: 17.91
  • 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

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)))
                            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      P>|t|      [0.025      0.975]
-----------------------------------------------------------------------------------------------------
Intercept                            87.6389      0.908     96.565      0.000      85.855      89.423
mom_hs_centeredZ                     13.8304     11.814      1.171      0.242      -9.391      37.051
mom_iq_centeredZ                      0.0785      0.008      9.712      0.000       0.063       0.094
mom_hs_centeredZ:mom_iq_centeredZ    -0.3144      0.105     -2.985      0.003      -0.521      -0.107
==============================================================================
Omnibus:                        8.014   Durbin-Watson:                   1.660
Prob(Omnibus):                  0.018   Jarque-Bera (JB):                8.258
Skew:                          -0.333   Prob(JB):                       0.0161
Kurtosis:                       2.887   Cond. No.                     1.54e+03
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.54e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
Residual SD: 17.91

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 transformations

A linear model on the logarithmic scale corresponds to a multiplicative model on the original scale.

Height and earnings example

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())
                            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>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept      5.7785      0.451     12.815      0.000       4.894       6.663
height         0.0588      0.007      8.743      0.000       0.046       0.072
==============================================================================
Omnibus:                      223.552   Durbin-Watson:                   1.904
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              457.878
Skew:                          -1.079   Prob(JB):                    3.74e-100
Kurtosis:                       5.136   Cond. No.                     1.17e+03
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.17e+03. This might indicate that there are
strong multicollinearity or other numerical problems.

Interpretation

A 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 to a change of 6%. On the other hand, a change on \(\log_{10}\) scale is harder to interpret:

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)))
                            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>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept      2.5096      0.196     12.815      0.000       2.125       2.894
height         0.0255      0.003      8.743      0.000       0.020       0.031
==============================================================================
Omnibus:                      223.552   Durbin-Watson:                   1.904
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              457.878
Skew:                          -1.079   Prob(JB):                    3.74e-100
Kurtosis:                       5.136   Cond. No.                     1.17e+03
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.17e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
Residual SD: 0.39

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

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)))
                            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>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept      8.1527      0.603     13.530      0.000       6.970       9.335
height         0.0207      0.009      2.218      0.027       0.002       0.039
male           0.4232      0.072      5.841      0.000       0.281       0.565
==============================================================================
Omnibus:                      225.577   Durbin-Watson:                   1.889
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              462.608
Skew:                          -1.088   Prob(JB):                    3.52e-101
Kurtosis:                       5.141   Cond. No.                     1.59e+03
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.59e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
Residual SD: 0.88

Including an interaction

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)))
                            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>|t|      [0.025      0.975]
-------------------------------------------------------------------------------
Intercept       8.3885      0.844      9.945      0.000       6.734      10.043
height          0.0170      0.013      1.304      0.193      -0.009       0.043
male           -0.0786      1.258     -0.062      0.950      -2.546       2.389
height:male     0.0074      0.019      0.400      0.690      -0.029       0.044
==============================================================================
Omnibus:                      225.961   Durbin-Watson:                   1.889
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              463.655
Skew:                          -1.089   Prob(JB):                    2.08e-101
Kurtosis:                       5.143   Cond. No.                     4.16e+03
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 4.16e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
Residual SD: 0.88

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:

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())
                            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>|t|      [0.025      0.975]
---------------------------------------------------------------------------------
Intercept         9.5266      0.045    210.878      0.000       9.438       9.615
height_z          0.0654      0.050      1.304      0.193      -0.033       0.164
male              0.4197      0.073      5.748      0.000       0.276       0.563
height_z:male     0.0286      0.072      0.400      0.690      -0.112       0.169
==============================================================================
Omnibus:                      225.961   Durbin-Watson:                   1.889
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              463.655
Skew:                          -1.089   Prob(JB):                    2.08e-101
Kurtosis:                       5.143   Cond. No.                         4.64
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
Residual SD: 0.88
66.91694630872483 3.846637848956869

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

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)))
                            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          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------------
Intercept           82.0000      2.305     35.568      0.000      77.469      86.531
C(mom_work)[T.2]     3.8542      3.095      1.245      0.214      -2.229       9.937
C(mom_work)[T.3]    11.5000      3.553      3.237      0.001       4.517      18.483
C(mom_work)[T.4]     5.2098      2.704      1.927      0.055      -0.105      10.524
==============================================================================
Omnibus:                       12.371   Durbin-Watson:                   1.477
Prob(Omnibus):                  0.002   Jarque-Bera (JB):               12.950
Skew:                          -0.411   Prob(JB):                      0.00154
Kurtosis:                       2.796   Cond. No.                         5.92
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
Residual SD: 20.16

Building regression models for prediction

mesquite_df = pd.read_csv("../data/mesquite/mesquite.tsv.gz", sep="\t")
mesquite_df.head()
Obs Group Diam1 Diam2 TotHt CanHt Dens LeafWt
0 1 MCD 1.8 1.15 1.30 1.00 1 401.3
1 2 MCD 1.7 1.35 1.35 1.33 1 513.7
2 3 MCD 2.8 2.55 2.16 0.60 1 1179.2
3 4 MCD 1.3 0.85 1.80 1.20 1 308.0
4 5 MCD 3.3 1.90 1.55 1.05 1 855.2
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)))
                            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      P>|t|      [0.025      0.975]
--------------------------------------------------------------------------------
Intercept    -1091.8880    176.456     -6.188      0.000   -1448.804    -734.972
Group[T.MCD]   363.2951    100.184      3.626      0.001     160.654     565.936
Diam1          189.6690    112.760      1.682      0.101     -38.410     417.748
Diam2          371.4621    124.378      2.987      0.005     119.883     623.041
CanHt          355.6653    209.843      1.695      0.098     -68.782     780.113
TotHt         -101.7325    185.574     -0.548      0.587    -477.091     273.626
Dens           131.2542     34.355      3.820      0.000      61.764     200.744
==============================================================================
Omnibus:                        9.865   Durbin-Watson:                   2.137
Prob(Omnibus):                  0.007   Jarque-Bera (JB):               19.723
Skew:                           0.366   Prob(JB):                     5.21e-05
Kurtosis:                       6.123   Cond. No.                         26.9
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
Residual SD: 250.39

Multiplicative Model

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)))
                            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      P>|t|      [0.025      0.975]
--------------------------------------------------------------------------------
Intercept        4.7680      0.155     30.747      0.000       4.454       5.082
Group[T.MCD]     0.5834      0.129      4.534      0.000       0.323       0.844
log(Diam1)       0.3938      0.282      1.397      0.170      -0.177       0.964
log(Diam2)       1.1512      0.210      5.477      0.000       0.726       1.576
log(CanHt)       0.3732      0.281      1.330      0.191      -0.194       0.941
log(TotHt)       0.3943      0.313      1.260      0.215      -0.239       1.027
log(Dens)        0.1093      0.122      0.896      0.376      -0.137       0.356
==============================================================================
Omnibus:                        1.935   Durbin-Watson:                   1.362
Prob(Omnibus):                  0.380   Jarque-Bera (JB):                1.047
Skew:                          -0.283   Prob(JB):                        0.593
Kurtosis:                       3.474   Cond. No.                         12.6
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
Residual SD: 0.31

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.

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)))
                            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>|t|      [0.025      0.975]
-------------------------------------------------------------------------------
Intercept       5.1697      0.083     62.066      0.000       5.002       5.338
log(CanVol)     0.7224      0.055     13.234      0.000       0.612       0.832
==============================================================================
Omnibus:                        2.765   Durbin-Watson:                   0.996
Prob(Omnibus):                  0.251   Jarque-Bera (JB):                2.066
Skew:                          -0.515   Prob(JB):                        0.356
Kurtosis:                       3.125   Cond. No.                         2.59
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
Residual SD: 0.41

So leaf weight is approximately proportional to 0.72th power to Canopy Volume. It has a impressive \(R^2\) of 0.8

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)))
                            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      P>|t|      [0.025      0.975]
--------------------------------------------------------------------------------
Intercept        4.7680      0.155     30.747      0.000       4.454       5.082
Group[T.MCD]     0.5834      0.129      4.534      0.000       0.323       0.844
log(CanVol)      0.3732      0.281      1.330      0.191      -0.194       0.941
log(CanArea)     0.3993      0.294      1.357      0.183      -0.196       0.994
log(Shape)      -0.3787      0.231     -1.642      0.109      -0.845       0.088
log(TotHt)       0.3943      0.313      1.260      0.215      -0.239       1.027
log(Dens)        0.1093      0.122      0.896      0.376      -0.137       0.356
==============================================================================
Omnibus:                        1.935   Durbin-Watson:                   1.362
Prob(Omnibus):                  0.380   Jarque-Bera (JB):                1.047
Skew:                          -0.283   Prob(JB):                        0.593
Kurtosis:                       3.474   Cond. No.                         22.2
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
Residual SD: 0.31
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)))
                            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      P>|t|      [0.025      0.975]
--------------------------------------------------------------------------------
Intercept        4.6975      0.118     39.812      0.000       4.459       4.936
Group[T.MCD]     0.5269      0.116      4.558      0.000       0.294       0.760
log(CanVol)      0.6148      0.191      3.215      0.003       0.229       1.001
log(CanArea)     0.2889      0.238      1.215      0.231      -0.191       0.768
==============================================================================
Omnibus:                        1.398   Durbin-Watson:                   1.282
Prob(Omnibus):                  0.497   Jarque-Bera (JB):                0.694
Skew:                          -0.261   Prob(JB):                        0.707
Kurtosis:                       3.300   Cond. No.                         13.3
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
Residual SD: 0.33