<jupyter_start><jupyter_text>BIOEE 4940 : **Introduction to Quantitative Analysis in Ecology** ***Spring 2021*** Instructor: **<NAME>** ( ✉️ <EMAIL>) Teaching Assistant: **Yanqiu (Autumn) Zhou** (✉️ <EMAIL>)--- Lecture 6 *Regression III: Model Selection**Partly adapted from [How to be a quantitative ecologist](https://www.researchgate.net/publication/310239832_How_to_be_a_Quantitative_Ecologist_The_'A_to_R'_of_Green_Mathematics_and_Statistics) and [All of Statistics](https://www.stat.cmu.edu/~larry/all-of-statistics/)*<jupyter_code># import packages and prepare data to be used
# import packages and read the data

import pandas as pd
import numpy as np


baad_data_url = 'https://raw.githubusercontent.com/xiangtaoxu/QuantitativeEcology/main/Lab1/baad_data.csv'
baad_dictionary_url = 'https://raw.githubusercontent.com/xiangtaoxu/QuantitativeEcology/main/Lab1/baad_dictionary.csv'

df_data = pd.read_csv(baad_data_url, encoding='latin_1') # can also read local files
df_dict = pd.read_csv(baad_dictionary_url, encoding='latin_1')

df_ms = df_data[['a.lf','d.bh','h.t','r.st','ma.ilf','species','family','map','mat']].dropna()

# rename columns to avoid mis-interpreation for formula-based regressions
df_ms.rename(columns={'d.bh' : 'dbh','h.t' : 'h','r.st':'wd','ma.ilf':'lma','a.lf' : 'la'},inplace=True)

df_ms.dropna(inplace=True)

df_ms.shape<jupyter_output><empty_output><jupyter_text>1. Goal of model selection (Tredennick et al. 2021):      * Exploration:        What are possible candidate covariates to include?            Hypotheses generation based on inductive reasoning (e.g. Is $X_i$ influencing Y).            Trade-off: thorough vs false discoveries (type I error)<jupyter_code>pd.plotting.scatter_matrix(df_ms[['la','dbh','h','wd','lma','map','mat']])
plt.show()
# which variables should be included to predict la?
# example of type I error for casting a wide net for explanatory variables
# recall what is type I error

import numpy as np
from scipy import stats

from statsmodels.stats.multicomp import MultiComparison


import matplotlib.pyplot as plt
%matplotlib inline


# generate 10 random samples
data = []
group = []
for i in range(20):
    data.extend(np.random.rand(20))
    group.extend(np.ones((20,)) * i)
    
rand_df = pd.DataFrame({'data' : data, 'group' : group})
                 
# conduct a pairwise correlation analysis
mcomp = MultiComparison(rand_df['data'],rand_df['group'])

res = mcomp.allpairtest(stats.pearsonr,alpha=0.05,method='Holm')
print(res[0])

pvals = np.array(res[0].data)[1::,3].astype(float)

print(np.sum(pvals <= 0.05) / len(pvals))<jupyter_output><empty_output><jupyter_text>* Inference:          Test detailed hypotheses in more rigorous ways.            Is $X_i$ and the underlying natural process an important explanatory factor for Y?            What is the senstivity of Y to $X_i$?            Usually involve comparing alternative models.            Need replication and validation across a range of conditions before being accepted as scientific fact            Collinearity can be a major challenge.<jupyter_code># example of wd effect
import statsmodels.formula.api as smf
import statsmodels.api as sm

# no wd
res0 = smf.ols('h ~ dbh',data=df_ms).fit()
print(res0.summary())


# wd influences intercept
res1 = smf.ols('h ~ dbh + wd',data=df_ms).fit()
print(res1.summary())

# wd influences slope
res2 = smf.ols('h ~ dbh + dbh:wd',data=df_ms).fit()
print(res2.summary())

# wd influences both
res3 = smf.ols('h ~ dbh * wd',data=df_ms).fit()
print(res3.summary())

# how to tell including one extra variable is statistically meaningful?
from statsmodels.stats.anova import anova_lm
print(anova_lm(res0,res1,res2,res3))
# likelihood test
# res.compare_lr_test

res3.compare_lr_test(res2)<jupyter_output><empty_output><jupyter_text>* Prediction:        One of the most common goal of statistical models (e.g. statistical down-scaling, ecological forecasting)            Given the equation $\hat{Y} = \hat{\beta}X$ acquired from regression, prediction is more about $\hat{Y}$ while inference is more about $\hat{\beta}$            Confronting model predictions with new data is the ultimate test of our understanding --> the importance of 'out of sample' predictions.            The optimal model for prediction may not be suitable for inference<jupyter_code># a simple experiment for different bewteen inference and predictions

x = np.random.rand(100)

y = 2. * x.copy() + 1. # y is linearly related with x

z = x ** 2 + 0.1 # z is linearly related with x^2

fig = plt.figure()
plt.scatter(x,y,c='k')
plt.scatter(x,z,c='r')

plt.show()
# now consider in the real world, observations of x,y,z all come with errors
# let's assume these errors are normal with zero means
sigma_x = 0.1
sigma_y = 0.1
sigma_z = 0.1
x_obs = x + np.random.randn(100) * sigma_x
y_obs = y + np.random.randn(100) * sigma_y
z_obs = z + np.random.randn(100) * sigma_z

fig = plt.figure()
plt.scatter(x_obs,y_obs,c='k')
plt.scatter(x_obs,z_obs,c='r')

plt.show()
# let's try some regressions
df_test = pd.DataFrame({'x' : x_obs,'y' : y_obs, 'z' : z_obs})

# we use the first 50 data points to train the model
# we then use the last 50 data to validate the model (out-of-sample test)

df_train = df_test.iloc[0:50]
df_valid = df_test.iloc[50::]

res1 = smf.ols('y ~ x',data=df_train).fit()
res2 = smf.ols('y ~ x + z',data=df_train).fit()
print(res1.summary())
print(res2.summary())
print(anova_lm(res1,res2))
# predictions
Y1_predict = res1.predict(df_valid)
Y2_predict = res2.predict(df_valid)

# compare RMSE
print('RMSE:')
print('Model 1:',np.sqrt(np.mean(np.power(Y1_predict - df_valid['y'],2))))
print('Model 2:',np.sqrt(np.mean(np.power(Y2_predict - df_valid['y'],2))))<jupyter_output><empty_output><jupyter_text>2. How to 'score' a model      In order to rank/select models, we first need a method to score each model.      Obviously, we have $R^2$, which denotes the fraction of variance explained by the independent variables. However, this metric has two potential problems: (1) R2 tends to always increase with new variables, adding to the risk of type I error; (2) R2 is ultimately a measure of 'explanatory' power of the model and has nothing to do with prediction, especially out-of-sample prediction.    We briefly talked about $R^2_{adj}$ before, which includes some penalty of model complexity (number of independent variables) by considering the changes in degree of freedom when calculating $R^2$.    Let's see a simple example below.<jupyter_code># examine how R2 and R2adj changes with increasing number of polynomial terms of dbh when explaining la

# create the data columns
df_ms['log_la'] = np.log(df_ms['la'])
for i in range(1,10+1):
    df_ms[f'log_dbh{i}'] = np.power(np.log(df_ms['dbh'].values),i)
    
# create regressions

reg_ress = [] # store regression results
reg_r2 = [] # store regression r2
reg_r2adj = []# stored regression r2 adj

for i in range(1,10+1):
    x_strs = [f'log_dbh{j}' for j in range(1,i+1)]
    
    reg_str = 'log_la ~ ' + ' + '.join(x_strs)
    
    print(reg_str)
    
    res = smf.ols(reg_str,df_ms).fit()
    
    reg_ress.append(res)
    reg_r2.append(res.rsquared)
    reg_r2adj.append(res.rsquared_adj)
    





fig = plt.figure()
x_val = range(1,10+1)
plt.plot(x_val,reg_r2,'k-o')
plt.plot(x_val,reg_r2adj,'r-s')

# the model with maximum r2 and r2adj
print(np.argmax(reg_r2))
print(np.argmax(reg_r2adj))
plt.show()

# visualize the last model
dbh_val = np.log(np.arange(0.001,1.,0.0001))
data_predict = {}
for i in range(1,10+1):
    data_predict[f'log_dbh{i}'] = np.power(dbh_val,i)
    
Y_pred = reg_ress[3].predict(pd.DataFrame(data_predict))

fig = plt.figure()
plt.scatter(df_ms['log_dbh1'],df_ms['log_la'])
plt.plot(dbh_val,Y_pred,'r-')
plt.show()

# r2_adj helps but is not perfect at excluding false positives
print(reg_ress[3].summary())<jupyter_output><empty_output><jupyter_text>* Akaike Information Criterion (AIC)    Althought $R^2_{adj}$ does not work perfectly, the underlying motivation makes sense, which separates the model score into **goodness of fit** and **penalty of complexity**.        AIC uses log-likelihood (recall what is likelihood - P(Y | model, X)) and number of model parameters to construct the score        AIC = - ($log_e(L)$ - 2k), where k is the number of parameters. The negative sign upfront is made so that models with smaller AIC will be 'better'.        AICc corrects for small sample size (e.g. < 20). AICc = AIC + $\frac{2k^2+2k}{n-k-1}$<jupyter_code>reg_aic = [] # store regression AIC

for i in range(1,10+1):
    x_strs = [f'log_dbh{j}' for j in range(1,i+1)]
    
    reg_str = 'log_la ~ ' + ' + '.join(x_strs)
    
    print(reg_str)
    
    res = smf.ols(reg_str,df_ms).fit()
    
    reg_aic.append(res.aic)
# compare with another model
res_h = smf.ols('log_la ~ log_dbh1 + np.log(h)',df_ms).fit()
print(res_h.summary())
fig = plt.figure()
x_val = range(1,10+1)
plt.plot(x_val,reg_aic,'b-d')
plt.show()

# the results look similar to R2_adj but they can differ in certain cases
# AIC is more recommended
# Check Burnham & Anderson 2004 for more details
# 1. <NAME>, <NAME>, Multimodel inference - understanding AIC and BIC in model selection. Sociol. Methods Res. 33, 261–304 (2004).<jupyter_output><empty_output><jupyter_text>* Cross-validation    If we are mostly interested in predictions, we can use **cross-validation** to estimate the predictive risk even without new data.        Cross-validation separates the whole data set into *training* and *validation* data sets. A general method is called **K-fold cross validation**, we divide the data into k groups; often 5-10. We omit one group of data and fit the models to the remaining data. We use the fitted model to predict the data in the groups that was omitted (out-of-sample). Note that when k is equal to n, it becomes **leave-one-out cross validation**        We can use RMSE to assess the performance of the model.<jupyter_code># use scikit learn package
# a powerful package for machine learning and data analysis
from sklearn.model_selection import KFold

kf = KFold(5)
kf_indexes = kf.split(df_ms)
# create data
kf_data = []
for train_idx, test_idx in kf_indexes:
    df_train = df_ms.iloc[train_idx]
    df_test  = df_ms.iloc[test_idx]
    
    kf_data.append((df_train,df_test))

# loop over models
mod_cv_rmse = []
for i in range(1,10+1):
    x_strs = [f'log_dbh{j}' for j in range(1,i+1)]
    
    reg_str = 'log_la ~ ' + ' + '.join(x_strs)
    
    print(reg_str)
    
    y_pre = []
    y_obs = []
    
    # loop over K-Fold groups
    for df_train,df_test in kf_data:
        
        res = smf.ols(reg_str,df_train).fit()
        y_obs.extend(df_test['log_la'].values)
        y_pre.extend(res.predict(df_test))
    
    y_obs = np.array(y_obs)
    y_pre = np.array(y_pre)
    
    
    
    mod_cv_rmse.append(np.sqrt(np.mean(np.power(y_obs - y_pre,2))))


fig = plt.figure()
x_val = range(1,10+1)
plt.plot(x_val,mod_cv_rmse,'b-d')
plt.show()
print(np.argmin(mod_cv_rmse))<jupyter_output><empty_output><jupyter_text>3. How to search for the bestNow that we know a few metrics to score each statistical model, the next question is how to search through all the models to find the best one.The most naive option is to conduct an **exhaustive** search. This method can ensure we find the optimal model from all the candidate explanatory variables. However, If there are *k* covariates, there are $2^k$ possible models. Even if we only have 10 candidate variables (which can be easily exceeded for modeling ecological processes which integrate various physical, chemical, and biological processes), there will be 1024 models to evaluate. This approach becomes computationally consuming and even infeasible when k is very large.Naturally, we will consider trimming the unnecessary model evaluations. One common set of methods is **stepwise regression**. It can be both *forward* or *backward*. In forward stepwise regression, we start with no covariates and then add the one variable that leads to the best score (requires evaluation of *k* different models). We continue adding variables one at a time (always the one leading to the best score, which requires k-1, k-2, ..., evaluations) until the score does not improve. Backwards stepwise regression is the same except that we start with all possible covariates and drop one variable at a time. Both are 'greedy' searches (i.e. maximize the gain at the current time step) and thus neither is guaranteed to find the model with the best score because greedy searches ensure convergence to local optimal solutions but not necessarily global optimal solutions. Despite its limitation, these are probably the most common search method.<jupyter_code># example for model covariate/feature selection to predict la
# use linear_regression from sklearn
from sklearn.linear_model import LinearRegression
# use Recursive Feature Elimination, largely equiavlent to backward stepwise selection
# R has a similar step function as well
from sklearn.feature_selection import RFECV

# Let's consider predict log_la
# We will use 10 dbh variables (1-> 10th power), 10 h variables (1 -> 10th power),
# 4 wd variables (wd, wd:dbh, wd:h, wd:dbh:h)
# 2 climate variables, MAT and MAP

df_rfe = df_ms[['log_la','mat','map','wd']]

# create dbh and h vars
for i in range(1,10+1):
    df_rfe[f'log_dbh{i}'] = np.power(np.log(df_ms['dbh'].values),i)
    df_rfe[f'log_h{i}'] = np.power(np.log(df_ms['h'].values),i)
    
# create wd interaction terms
df_rfe['log_wd'] = np.log(df_rfe['wd'])
df_rfe['log_wd_dbh'] = df_rfe['log_wd'] * df_rfe['log_dbh1']
df_rfe['log_wd_h'] = df_rfe['log_wd'] * df_rfe['log_h1']
df_rfe['log_wd_dbh_h'] = df_rfe['log_wd'] * df_rfe['log_dbh1'] *df_rfe['log_h1']

print(df_rfe.shape)
# have 27 covariates/features in total

# use rfe to select 10 predictors
rfe = RFECV(LinearRegression(),verbose=1)
# get X
x_str = df_rfe.columns.tolist()
x_str.remove('log_la')

rfe = rfe.fit(df_rfe[x_str].values,df_rfe['log_la'].values)
print(rfe.support_)
print(rfe.ranking_)
print(np.array(x_str)[rfe.support_])
print(rfe.grid_scores_) # psuedo R2, 1- RSS/TSS<jupyter_output><empty_output><jupyter_text>* RegularizationAn alternative to stepwise model selection is to apply regression with **statistical regularization**. It includes *ridge regression, LASSO (least absolute shrinkage and selection operator), and the elastic net*.The idea behind regularization is similar to model scoring with consideration of both goodness of fit and model complexity. Instead of analyzing the score after fitting the model, statisticians developed method to directly include the model complexity penalty during model fitting so that we can the model fitting process can directly 'select features'.Some more details of regularization can be found in the supporting materials of Tredennick et al. (2021). In brief, for the case of OLS regression, instead of minimizing the squared residuals (good-ness of fit), we can also minimize the penalty associated with model complexity$loss function = \sum_{i=1}{n} (y_i - \beta_0 - \sum_{j=1}^p x_{ij}\beta_j)^2 + \gamma_1\sum_{j=1}^{p} |\beta_j|^{\gamma_2}$.Here $\gamma_1$ is referred to as the 'penalization' parameter determining the strenght of penalty towards more complex models.When $\gamma_2$ is equal to 2, it becomes **Ridge regression**. This will penalize on models with too large coefficients (which is usually the case for model overfitting), which will lead to "shrinkage" of model parametersWhen $\gamma_2$ is equal to 1, it becomes **LASSO**, which can shrink the coefficient estimates all the way to **zero**. This will lead to automatic feature selection during model fitting. The  of final selected covariates depend on the value of $\gamma_1$**Elastic net** is a combination of Ridge regression and LASSO that includes both kinds of regulators with adjustable weights for each regulator.<jupyter_code># Ridge/LASSO regression for polynomial fit between la and dbh
x_str = [f'log_dbh{i}' for i in range(1,10+1)]

# first compare it with a linear model
from sklearn.linear_model import LinearRegression
# it is recommended to standarize all variables for Ridge/Lasso/ElasticNet
# we do the same for linear regression
from sklearn import preprocessing

scaler = preprocessing.StandardScaler().fit(df_ms[x_str + ['log_la']].values)
data_scaled = scaler.transform(df_ms[x_str + ['log_la']].values)

X_scaled = data_scaled[:,0:-1]
y_scaled = data_scaled[:,-1]

res_lin = LinearRegression().fit(X_scaled,y_scaled)
print(res_lin.score(X_scaled,y_scaled)) # get R2
print(res_lin.coef_) # get coef
from sklearn.linear_model import RidgeCV

# test different alpha/gamma_1 values
res_ridge = RidgeCV(alphas=np.logspace(-4,1,10)).fit(X_scaled, y_scaled)
print(res_ridge.score(X_scaled, y_scaled)) # R2
print(res_ridge.coef_)
print(res_ridge.alpha_) # best alpha

# test to use different alphas
print(np.logspace(-4,1,10))
from sklearn.linear_model import LassoCV
# test different alpha/gamma_1 values
res_lasso = LassoCV(alphas=np.logspace(-4,1,10)).fit(X_scaled, y_scaled)
print(res_lasso.score(X_scaled, y_scaled)) # R2
print(res_lasso.coef_)
print(res_lasso.alpha_) # best alpha

# test using different alphas
# compare results
# visualize three different models
dbh_val = np.log(np.arange(0.001,1.,0.0001))
# the last one is dummy variable
X = np.array([np.power(dbh_val,i) for i in range(1,11+1)]).T
print(X.shape)

X_scaled = scaler.transform(X)[:,0:10] # only the first 10
    
y_lin = res_lin.predict(X_scaled)
y_ridge = res_ridge.predict(X_scaled)
y_lasso = res_lasso.predict(X_scaled)

fig = plt.figure()
plt.scatter(scaler.transform(df_ms[x_str + ['log_la']].values)[:,0],y_scaled)
plt.plot(X_scaled[:,0],y_lin,'k-',label='OLS')
plt.plot(X_scaled[:,0],y_ridge,'r-',label='Ridge')
plt.plot(X_scaled[:,0],y_lasso,'g-',label='LASSO')
plt.legend()
plt.show()


# in practice, try different values of alpha to get the best Cross Validation results
# we will use elastic net as an example to select the best model to predict la
from sklearn.linear_model import ElasticNetCV
from sklearn import preprocessing

# get X
x_str = df_rfe.columns.tolist()
x_str.remove('log_la')

data_str = x_str + ['log_la'] # the last column is log_la

scaler = preprocessing.StandardScaler().fit(df_rfe[data_str].values)
data_scaled = scaler.transform(df_rfe[data_str].values)

X_scaled = data_scaled[:,0:-1]
y_scaled = data_scaled[:,-1]

res_ee = ElasticNetCV(l1_ratio=np.logspace(-4,0,20)).fit(X_scaled,y_scaled) # half way between ridge and lasso

print(np.logspace(-4,0,20))
print(res_ee.alpha_) # the best alpha
print(res_ee.l1_ratio_) # the best l1_ratio
print(res_ee.score(X_scaled,y_scaled)) # best model score
print(pd.DataFrame({'covariate': x_str,'coef': res_ee.coef_})) # best model coefficient<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11623831775700935, 'fraction_numerical': 0.017470263381478335, 'mean_word_length': 4.81804139635465, 'pattern_counts': {'":': 0, '<': 44, '<?xml version=': 0, '>': 45, 'https://': 4, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '141881', 'license': 'MIT', 'n_tokens_mistral': 6604, 'n_tokens_neox': 6309, 'n_words': 1971, 'path': 'Lecture7/Lecture7_model_selection.ipynb', 'repo_name': 'xiangtaoxu/QuantitativeEcology', 'size': '18872'}