<jupyter_start><jupyter_text>Comparing Regularization Techniques---Advanced Applied Machine Learning- Lab 2By <NAME> ![variableselections.png](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAfIAAAFsCAYAAAA3/TvjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeVxU1fvA8c8AoiIooIBpaGpqKpqmqIWouATue4KKe2XuaaVkbrmUWWZhaqGpqT81Nbfc1xRLEjQLJM0lZVEWZRt25P7+mK+TxOIAA8PE8369eDVz7r3nPnOYeDznnnuuSlEUBSGEEEIYJRNDByCEEEKIopNELoQQQhgxSeRCCCGEEZNELoQQQhgxSeRCCCGEEZNELoQQQhgxSeSiTPDy8mLr1q25ynfs2IGXl5fO9fz++++MGzfuqft16dKFwMDAXOUBAQF0795d5/M99t1339GnTx88PDzo2rUrs2bNQq1WF3hMeHg4TZs2LfS5Hrty5Qp//vknAFu2bGHlypVFruvfVqxYQYcOHdi9e3ex6ilKuxTk3LlzREZGAuDr68ucOXMKXcf333+vfT1q1ChCQkJ0Pnb27Nm0b98eDw8P3Nzc6Nu3L5s2bSI7O/upx2ZkZLB3795Cx6uLW7ducfHixRKpW5R9kshFmTBw4EAOHDiQq3zfvn0MHDhQ53patGjB+vXr9R<KEY>//<KEY>5EEIIYcKkkAshhBAmTAq5EEIIYcKkkAshhBAm7P8BD8fvG+gNS20AAAAASUVORK5CYII=) This lab will compare multiple different regularization techniques on both real and simulated data. Each regressor will be run with multiple different hyperperameter settings to better compare that regressor on each dataset. Regressors Compared:1.   [Ridge Regression](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Ridge.html)2.   [Lasso Regression](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Lasso.html)3.   [Elastic Net](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.ElasticNet.html) 4.   [SCAD (Smoothly Clipped Absolute Deviation)](https://hmjianggatech.github.io/picasso/pycasso.html)5.   [Square Root Lasso](https://www.statsmodels.org/stable/generated/statsmodels.regression.linear_model.OLS.fit_regularized.html)<jupyter_code>import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import keras
from keras.models import Sequential
from keras.layers import Dense
from keras.layers import Dropout
from sklearn.metrics import r2_score
from yellowbrick.regressor import ResidualsPlot
import sklearn.datasets
import sklearn
import sklearn.datasets as dt
from sklearn.model_selection import train_test_split as tts
from sklearn.linear_model import Lasso
from sklearn.preprocessing import StandardScaler
def scad_penalty(beta_hat, lambda_val, a_val):
    is_linear = (np.abs(beta_hat) <= lambda_val)
    is_quadratic = np.logical_and(lambda_val < np.abs(beta_hat), np.abs(beta_hat) <= a_val * lambda_val)
    is_constant = (a_val * lambda_val) < np.abs(beta_hat)
    
    linear_part = lambda_val * np.abs(beta_hat) * is_linear
    quadratic_part = (2 * a_val * lambda_val * np.abs(beta_hat) - beta_hat**2 - lambda_val**2) / (2 * (a_val - 1)) * is_quadratic
    constant_part = (lambda_val**2 * (a_val + 1)) / 2 * is_constant
    return linear_part + quadratic_part + constant_part
    
def scad_derivative(beta_hat, lambda_val, a_val):
    return lambda_val * ((beta_hat <= lambda_val) + (a_val * lambda_val - beta_hat)*((a_val * lambda_val - beta_hat) > 0) / ((a_val - 1) * lambda_val) * (beta_hat > lambda_val))
def plotModelVariables(modeltype, model1, alpha1, model2, alpha2, model3, alpha3, data_name):
  if data_name == "Synthetic":
    XTrain, XTest, YTrain, YTest = synthXTrain, synthXTest, synthYTrain, synthYTest
    ticks = list(range(0, 6))
    xticks = synth_xticks
  elif data_name == "Diabetes":
    XTrain, XTest, YTrain, YTest = diabetesXTrain, diabetesXTest, diabetesYTrain, diabetesYTest
    ticks = list(range(0, 10))
    xticks = diabetes_xticks
  elif data_name == "Boston Housing":
    XTrain, XTest, YTrain, YTest = bostonXTrain, bostonXTest, bostonYTrain, bostonYTest
    ticks = list(range(0, 12))
    xticks = boston_xticks
  plt.plot(np.reshape(model1.coef_, (len(xticks), 1)),alpha=0.7,linestyle='none',marker='*',markersize=8,color='red',label=r'%s; parameter = %s' %(modeltype, str(alpha1)), zorder=7) 
  plt.plot(np.reshape(model2.coef_, (len(xticks), 1)),alpha=0.5,linestyle='none',marker='d',markersize=7,color='blue',label=r'%s; parameter = %s' %(modeltype, str(alpha2)))
  plt.plot(np.reshape(model3.coef_, (len(xticks), 1)),alpha=0.5,linestyle='none',marker='o',markersize=7,color='green',label=r'%s; parameter = %s' %(modeltype, str(alpha3)))  
  plt.xticks(ticks, xticks)
  plt.title("%s Regression Variable Selection for %s Dataset" %(modeltype, data_name))
  plt.xlabel('Variable Index',fontsize=16)
  plt.ylabel('Coefficient Magnitude',fontsize=16)
  plt.legend(fontsize=13,loc=0)
  plt.show()
def fitModel(model, dataset, ModelName, best = False):
  newmodel = model
  units = ''
  if dataset == "Synthetic":
    XTrain, XTest, YTrain, YTest = synthXTrain, synthXTest, synthYTrain, synthYTest
  elif dataset == "Diabetes":
    XTrain, XTest, YTrain, YTest = diabetesXTrain, diabetesXTest, diabetesYTrain, diabetesYTest
  elif dataset == "Boston Housing":
    XTrain, XTest, YTrain, YTest = bostonXTrain, bostonXTest, bostonYTrain, bostonYTest
    units = "$"
  else:
    print("Incorrect Dataset, try again")
    return
  newmodel.fit(XTrain, YTrain)
  yhat = newmodel.predict(XTest)
  mse = sklearn.metrics.mean_squared_error(YTest, yhat)
  if dataset == "Boston Housing":
      mse = mse*1000
  model_train_score=newmodel.score(XTrain,YTrain)
  model_test_score=newmodel.score(XTest,YTest)
  print("%s R^2 Accuracy on %s dataset:" %(ModelName, dataset))
  print("Train set : %f" %(model_train_score))
  print("Test set: %f" %(model_test_score))
  print("Mean Square Error: %s %.2f" %(units, mse))
  if best == True:
    return newmodel, mse, (model_test_score)
  return newmodel

def printImportances(model, modelName, data_name):
  if data_name == "Synthetic":
    ticks = list(range(0, 6))
    xticks = synth_xticks
  elif data_name == "Diabetes":
    ticks = list(range(0, 10))
    xticks = diabetes_xticks
  elif data_name == "Boston Housing":
    ticks = list(range(0, 12))
    xticks = boston_xticks
  importances = model.coef_
  print(importances)
  for i,v in enumerate(importances):
	  print('Feature: %0d, Score: %.5f' % (i,v))
#  plot feature importance
  plt.bar([x for x in range(len(importances))], importances)
  plt.xticks(ticks, xticks)
  plt.title("%s Regression Variable Importance for %s Dataset" %(modelName, data_name))
  plt.xlabel('Variable Index',fontsize=16)
  plt.ylabel('Coefficient Magnitude',fontsize=16)
  plt.show()<jupyter_output><empty_output><jupyter_text>Datasets---  [Boston Housing](https://www.cs.toronto.edu/~delve/data/boston/bostonDetail.html)The Boston Housing Dataset is a dataset of 506 observations and 12 X variables. Each observation represents the average characteristics of a neighborhood or block in Boston, Massachusetts. The y variable that we will try to predict is the Median Price of a house in that neighborhood, in units of $1000. The X variables for this dataset are described by the data host as -1.   CRIM - per capita crime rate by town2.   ZN - proportion of residential land zoned for lots over 25,000 sq.ft.3.   INDUS - proportion of non-retail business acres per town.4.   CHAS - Charles River dummy variable (1 if tract bounds river; 0 otherwise)5.   NOX - nitric oxides concentration (parts per 10 million)6.   RM - average number of rooms per dwelling7.   AGE - proportion of owner-occupied units built prior to 19408.   DIS - weighted distances to five Boston employment centres9.   RAD - index of accessibility to radial highways10.  TAX - full-value property-tax rate per \$10,00011.  PTRATIO - pupil-teacher ratio by town12.  B - 1000(Bk - 0.63)^2 where Bk is the proportion of African-American residents 13.  LSTAT - % of residents considered low-income or impoverishedThe dataset appears to have been "censored" at $50,000, with an unusual number of datapoints recording that exact price. The data will be centered and normalized by Z-score before analysis is begun, to make it appropriate for Ridge and Lasso regressions.<jupyter_code>boston = sklearn.datasets.load_boston()
bostonY = boston['target']
bostonDF = pd.DataFrame(boston.data, columns = boston.feature_names)
#show data before scaling
bostonDF['median price'] = bostonY
bostonX = bostonDF[['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD', 'PTRATIO', 'B', 'LSTAT']].values
scaler = StandardScaler()
scaler.fit_transform(bostonX)
boston_xticks = ['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD', 'PTRATIO', 'B', 'LSTAT']
bostonDF.head()


bostonXTrain, bostonXTest, bostonYTrain, bostonYTest = tts(bostonX, bostonY, test_size = .2, random_state = 1693, shuffle = True)<jupyter_output><empty_output><jupyter_text>[Diabetes](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.load_diabetes.htmlsklearn.datasets.load_diabetes)The Diabetes dataset represents blood measurements taken from 442 diabetes patients. The dataset has 10 X variables. The y variable is a 'quantitative measure of disease progression.' The X variables of this dataset represent-1.   'age'- age in years2.   'sex'- sex of patient3.   'bmi'- body mass index of patient4.   'bp'-   average blood pressure of patient5.    "s1" - TC, total serum cholesterol6.    "s2" - LDL, low-density lipoproteins7.    "s3"-  HDL, high-density lipoproteins8.    "s4" -  TCH, total cholesterol / HDL9.    "s5" -  LTG, log of serum triglycerides level10.   "s6" -  GLU, blood sugar levelThe data has already been standardized and centered, and will not be further pre-processed before analysis.<jupyter_code>diabetes = sklearn.datasets.load_diabetes()
diabetesY = diabetes['target']
diabetesDF = pd.DataFrame(diabetes.data, columns = diabetes.feature_names)
diabetesDF['disease progression'] = diabetesY
diabetesX = diabetesDF[['age', 'sex', 'bmi', 'bp', 's1', 's2', 's3', 's4', 's5', 's6']]
diabetes_xticks = ['age', 'sex', 'bmi', 'bp', 's1', 's2', 's3', 's4', 's5', 's6'] 
diabetesDF.head()

diabetesXTrain, diabetesXTest, diabetesYTrain, diabetesYTest = tts(diabetesX, diabetesY, test_size = .2, random_state = 1693, shuffle = True)<jupyter_output><empty_output><jupyter_text>[Synthetic Data](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.make_regression.html)A set of sythetic data will be created with a prescribed covariance matrix and 6 input features. Some noise is added to the data. There are 1000 observations and 6 X variables.<jupyter_code>num_samples =1000
#define strength of noise
sigma = 4

# The desired mean values of the sample.
mu = np.array([5.0, 0.0, 10.0, 7, 1.0, 25.0])
#Matrix size
matrixSize = 6 
#creating random positive-semidefinite square matrix 
A = np.random.rand(matrixSize, matrixSize)
B = np.dot(A, A.transpose())

# Generate the random samples.
synthX = np.random.multivariate_normal(mu, B, size=num_samples)
beta = np.array([2,0,-3, 0, 5, 0]).reshape(-1,1)
synthY = np.matmul(synthX,beta) + sigma*np.random.normal(0,1,[num_samples,1])
synth_xticks = ['x1', 'x2', 'x3', 'x4', 'x5', 'x6'] 
synthX = pd.DataFrame(synthX, columns = synth_xticks)
synthX
betaDF = pd.DataFrame(beta, index = synth_xticks, columns= ["True Values"])
betaDF
synthXTrain, synthXTest, synthYTrain, synthYTest = tts(synthX, synthY, test_size = .2, random_state = 1693, shuffle = True)<jupyter_output><empty_output><jupyter_text>[Ridge Regression](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Ridge.html)Ridge Regression is a linear least squares regression that uses L2 regularization.Parameter: *Alpha*---Ridge's *alpha* parameter controls the strength of the regularization- higher values correspond to a stronger regularization. We'll demonstrate 3 different Ridge fits, utilizing *alpha* 1, 10, and 100.<jupyter_code>from sklearn.linear_model import Ridge

ridgeRegressor1 = Ridge(alpha = 1.0, fit_intercept = False)

ridgeRegressor10 = Ridge(alpha = 10.0, fit_intercept = False)

ridgeRegressor100 = Ridge(alpha = 100, fit_intercept = False)<jupyter_output><empty_output><jupyter_text>Ridge - BostonOn the Boston dataset, we can see that many variable coefficients are shrunk towards zero. The models with higher *alpha* settings shrink variable coefficients closer to 0 due to the stronger regression and higher error penalties.\ Of Rooms is the most-predictive variable and is given the highest positive coefficient by all 3 models. Proximity to the river also appears to be a positively correlated trait. Nitrous Oxide concentration is given a negative coefficient by most models, suggesting that this is an undesirable trait for homeowners.The Ridge Regression with *alpha* = 1 has the lowest Mean Square Error and highest R2 on the test set, and will be compared with the other models in the conclusion of this project.<jupyter_code>Ridge1Boston, Ridge1BostonMSE, Ridge1BostonR2  = fitModel(ridgeRegressor1, 'Boston Housing', 'Ridge 1', best = True)
Ridge10Boston = fitModel(ridgeRegressor10, 'Boston Housing', 'Ridge 10')
Ridge100Boston = fitModel(ridgeRegressor100, 'Boston Housing', 'Ridge 100')
plotModelVariables("Ridge", ridgeRegressor1, "1", ridgeRegressor10, '10', ridgeRegressor100, '100', "Boston Housing")<jupyter_output><empty_output><jupyter_text>Ridge - Diabetes DataThe Ridge regression appears to be poorly suited to predicting the Diabetes data. All 3 models have very poor R2 values, suggesting that a linear regression is not an appropriate fit for the Diabetes dataset.The *alpha* = 100 regressor reduced almost all coefficients very close to 0. The messy R2 values make analysis of the other variable coefficients unreliable.The *alpha* = 1 regressor had the lowest Mean Square Error and will be compared with the other models at the conclusion of this paper.<jupyter_code>Ridge1Diabetes, Ridge1DiabetesMSE, Ridge1DiabetesR2 = fitModel(ridgeRegressor1, 'Diabetes', 'Ridge 1', best = True)
Ridge10Diabetes = fitModel(ridgeRegressor10, 'Diabetes', 'Ridge 10')
Ridge100Diabetes = fitModel(ridgeRegressor100, 'Diabetes', 'Ridge 100')
plotModelVariables("Ridge", ridgeRegressor1, "1", ridgeRegressor10, '10', ridgeRegressor100, '100', "Diabetes")<jupyter_output><empty_output><jupyter_text>Ridge - Synthetic DataAll three of the Ridge models roughly selected the correct variables for the synthetic data.They each identified X1 and X5 as having the strong positive coefficients (2 and 5, respectively), while X3 had a negative coefficient.Ridge with *alpha* = 1 had the highest R2 on the test set and lowest MSE, and will be compared to the other models in the conclusion of this project.<jupyter_code>betaDF
Ridge1Synthetic, Ridge1SyntheticMSE, Ridge1SyntheticR2 = fitModel(ridgeRegressor1, 'Synthetic', 'Ridge 1', best = True)
Ridge10Synthetic = fitModel(ridgeRegressor10, "Synthetic", 'Ridge10')
Ridge100Synthetic = fitModel(ridgeRegressor100, 'Synthetic', 'Ridge 100')
plotModelVariables("Ridge", ridgeRegressor1, "1", ridgeRegressor10, '10', ridgeRegressor100, '100', "Synthetic")<jupyter_output><empty_output><jupyter_text>[Lasso Regression](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Lasso.html)Lasso, an acronym for Least Absolute Shrinkage and Selection Operator, is a linear regression method that uses the L1 regularization technique. L1 regularization will force some variable coefficient weights to zero, eliminating them from use in predictions. Parameter: *Alpha*---Lasso's *alpha* parameter controls the strength of the L1 regularization. I demonstrate Lasso regressions using *alpha* values of 1, 10, and 100.<jupyter_code>lasso1 = Lasso(alpha = 1)
lasso10 = Lasso(alpha = 10)
lasso100 = Lasso(alpha = 100)<jupyter_output><empty_output><jupyter_text>Lasso - Boston Housing DataThe Lasso regressors struggle with the Boston Housing dataset. All three of them have high Mean Square Error and low R2 values. Both Lasso with *alpha* 10 and 100 reduce almost all variables to 0, selecting  of Rooms as an important positively correlated variable, and LSTAT (% of residents impoverished) as a negatively correlated variable. Lasso with *alpha* = 1 has the highest R2 and lowest MSE, and will be compared to the other regressors in the conclusion. Interestingly, this regressor mostly selects DIS (distance from employment centers), PTRATIO (Pupil/Teacher Ratio in Nearby Schools), and LSTAT (% of residents in poverty) as negative corellations with stronger coefficients than  of Rooms.<jupyter_code>lasso1Boston, lasso1BostonMSE, lasso1BostonR2 = fitModel(lasso1, 'Boston Housing', 'Lasso 1', best = True)
lasso10Boston = fitModel(lasso10, 'Boston Housing', 'Lasso 10')
lasso100Boston = fitModel(lasso100, 'Boston Housing', 'Lasso 100')
plotModelVariables("Lasso", lasso1, "1", lasso10, '10', lasso100, '100', "Boston Housing")<jupyter_output><empty_output><jupyter_text>Lasso - Diabetes DataLasso regression struggles with the Diabetes dataset as well. Both the *alpha* = 10 and 100 regressors set all variables to zero, showing that this is too high of an *alpha* for this dataset. The *alpha* = 1 Lasso has the best performance of the three, though it still has very low R2 and high MSE. It reduces all variables to 0 other than BMI, s5,  and Blood Pressure.<jupyter_code>lasso1Diabetes, lasso1DiabetesMSE, lasso1DiabetesR2 = fitModel(lasso1, 'Diabetes', 'Lasso 1', best = True)
lasso10Diabetes = fitModel(lasso10, 'Diabetes', 'Lasso 10')

lasso100Diabetes = fitModel(lasso100, 'Diabetes', 'Lasso 100')
plotModelVariables("Lasso", lasso1, "1",lasso10, '10', lasso100, "100", "Diabetes")<jupyter_output><empty_output><jupyter_text>Lasso - Synthetic DataThe Lasso regressors perform poorly on the Synthetically generated dataset. The *alpha* =10 and 100 regressors have too high of a regularization penalty and reduce all variables to 0. The *alpha* = 1 regressor selects x1 and x5 as variables with high positive correlation, but reduces all other variables to 0, missing the negative correlation in x3.The *alpha* = 1 regressor shows the highest R2s and lowest MSE and will be compared to the other models in the conclusion of this paper.<jupyter_code>betaDF
lasso1Synthetic, lasso1SyntheticMSE, lasso1SyntheticR2 = fitModel(lasso1, 'Synthetic', 'Lasso 1', best = True)
lasso10Synthetic = fitModel(lasso10, 'Synthetic', 'Lasso 10')
lasso100Synthetic = fitModel(lasso100, 'Synthetic', 'Lasso 100')
plotModelVariables("Lasso", lasso1, "1", lasso10, '10', lasso100, "100", "Synthetic")<jupyter_output><empty_output><jupyter_text>[Elastic Net](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.ElasticNet.html)Elastic Net is another method of linear regression. Where Ridge and Lasso Regressions use L2 and L1 regularizations respectively, the Elastic Net regressor uses a combination of both regularization methods.  *L1 Ratio* Parameter---Elastic Net's *L1 Ratio* Parameter determines how the regressor combines the L2 and L1 regularization methods. When l1_ratio = 1, the regression is equivalent to a Lasso regression. l1_ratio = 0 corresponds to a Ridge regression. Values between 0 and 1 present combine both regularization methods.I demonstrate Elastic Net regressors using *l1_ratio* of  0.2, 0.5, and 0.8<jupyter_code>from sklearn.linear_model import ElasticNet
elastic20 = ElasticNet(alpha=1.0, l1_ratio=0.2)
elastic50 = ElasticNet(alpha = 1.0, l1_ratio = 0.5)
elastic80 = ElasticNet(alpha=1.0, l1_ratio = 0.8)<jupyter_output><empty_output><jupyter_text>Elastic Net - Boston Housing DataThe Elastic Net regressors all perform similarly well on the Boston Housing Dataset. They all identify  of Rooms as a positively correlated variable, and DIS, Lstat and PtRatio as negatively correlated variables- similar to the Lasso regressor's variable selection.The Elastic Net with L1 Ratio = 20% has the highest R2 and lowest MSE, and will be compared to other models on this dataset in the conclusion.<jupyter_code>ENet20Boston, ENet20BostonMSE, ENet20BostonR2 = fitModel(elastic20, 'Boston Housing', 'Elastic L1 = 20', best = True)
ENet50Boston = fitModel(elastic50, 'Boston Housing', 'Elastic L1 = 50')
ENet80Boston = fitModel(elastic80, 'Boston Housing', 'Elastic L1 = 80')
plotModelVariables("Elastic Net", ENet20Boston, "20", ENet50Boston, '50', ENet80Boston, "80", "Boston Housing")<jupyter_output><empty_output><jupyter_text>Elastic Net - Diabetes DataLike the Ridge and Lasso regressions, Elastic Net struggles with the Diabetes dataset. All three models have extremely low R2 and high MSE. The Elastic Net regressor with L1_Ratio = 80% has the highest R2 and lowest MSE (even though both metrics are pretty poor).<jupyter_code>ENet20Diabetes = fitModel(elastic20, 'Diabetes', 'Elastic L1 = 20')
ENet50Diabetes = fitModel(elastic50, 'Diabetes', 'Elastic L1 = 50')
ENet80Diabetes, ENet80DiabetesMSE, ENet80DiabetesR2 = fitModel(elastic80, 'Diabetes', 'Elastic L1 = 80', best = True)
plotModelVariables("Elastic Net", ENet20Diabetes, "20", ENet50Diabetes, '50', ENet80Diabetes, "80", "Diabetes")<jupyter_output><empty_output><jupyter_text>Elastic Net - Synthetic DatasetThe Elastic Net regressors perform better on the Synthetic data than they did on the more complicated Diabetes dataset, though they still perform worse than the Ridge and Lasso regressions. The regressors identify X1 and X5 as the most positively correlated variables, but give them equivalent coefficients when X5 is truly more than double X1. The regressors also identify X3's negative correlation, but give it a much smaller coefficient than it should have in reality. The Elastic Net with L1_Ratio = 20% has the best performance on this dataset and will be compared to the other models at the end of this paper.<jupyter_code>betaDF
ENet20Synthetic, ENet20SyntheticMSE, ENet20SyntheticR2 = fitModel(elastic20, 'Synthetic', 'Elastic L1 = 20%', best = True)
ENet50Synthetic = fitModel(elastic50, 'Synthetic', 'Elastic L1 = 50%')
ENet80Synthetic = fitModel(elastic80, 'Synthetic', 'Elastic L1 = 80%')
plotModelVariables("Elastic Net", ENet20Synthetic, "20", ENet50Synthetic, '50', ENet80Synthetic, "80", "Synthetic")<jupyter_output><empty_output><jupyter_text>[SCAD (Smoothly Clipped Absolute Deviation)](https://fan.princeton.edu/papers/01/penlike.pdf)The SCAD Estimator, proposed by Fan and Li in the above-linked paper, is a modification of Lasso regression. SCAD maintains Lasso's desirable qualities such as handling orthoganal cases and providing continuous solutions, but minimizes some of Lasso's problems when dealing with large regression coefficients.I use Professor Vasiliu's SCAD implementation demonstrated in lecture, which is partially adapted from [this](https://andrewcharlesjones.github.io/journal/scad.html) article by <NAME>. I won't show any hyperparameter options for this regressor.<jupyter_code>def scad_penalty(beta_hat, lambda_val, a_val):
    is_linear = (np.abs(beta_hat) <= lambda_val)
    is_quadratic = np.logical_and(lambda_val < np.abs(beta_hat), np.abs(beta_hat) <= a_val * lambda_val)
    is_constant = (a_val * lambda_val) < np.abs(beta_hat)
    
    linear_part = lambda_val * np.abs(beta_hat) * is_linear
    quadratic_part = (2 * a_val * lambda_val * np.abs(beta_hat) - beta_hat**2 - lambda_val**2) / (2 * (a_val - 1)) * is_quadratic
    constant_part = (lambda_val**2 * (a_val + 1)) / 2 * is_constant
    return linear_part + quadratic_part + constant_part
    
def scad_derivative(beta_hat, lambda_val, a_val):
    return lambda_val * ((beta_hat <= lambda_val) + (a_val * lambda_val - beta_hat)*((a_val * lambda_val - beta_hat) > 0) / ((a_val - 1) * lambda_val) * (beta_hat > lambda_val))
def scad(beta):
  beta = beta.flatten()
  beta = beta.reshape(-1,1)
  n = len(y)
  return 1/n*np.sum((y-X.dot(beta))**2) + np.sum(scad_penalty(beta,lam,a))
  
def dscad(beta):
  beta = beta.flatten()
  beta = beta.reshape(-1,1)
  n = len(y)
  return np.array(-2/n*np.transpose(X).dot(y-X.dot(beta))+scad_derivative(beta,lam,a)).flatten()<jupyter_output><empty_output><jupyter_text>SCAD - Boston Housing DataSCAD identifies that  of Rooms is the most positively correlated variable with housing price, but it puzzlingly identifies Nitrous Oxide concentration as a positive factor, when most other models identified a negative correlation between this variable or reduced this variable to 0.It reduces most other variables close to 0 other than Distance from Employment Centers and Poor % of Neighborhood Population, which matches the results from most other models.<jupyter_code>import scipy
lam = 0.5
a = 1.001
p = bostonXTrain.shape[1]
b0 = np.random.normal(1,1,p)
y = bostonYTrain
X = bostonXTrain
y = y.reshape((y.shape[0]), 1)
bostonSCAD = scipy.optimize.minimize(scad, b0 , method='L-BFGS-B', jac=dscad,options={'gtol': 1e-8, 'maxiter': 50000,'maxls': 25,'disp': True})
bostonSCAD.x
importances = bostonSCAD.x
for i,v in enumerate(importances):
  print('Feature: %0d, Score: %.5f' % (i,v))
plt.bar([x for x in range(len(importances))], importances)
plt.xticks(list(range(0,12)), boston_xticks)
plt.title("%s Regression Variable Importance for %s Dataset" %("SCAD", "Boston Housing"))
plt.xlabel('Variable Index',fontsize=16)
plt.ylabel('Coefficient Magnitude',fontsize=16)
plt.show()<jupyter_output>Feature: 0, Score: -0.12910
Feature: 1, Score: 0.07186
Feature: 2, Score: -0.19533
Feature: 3, Score: 0.36914
Feature: 4, Score: 1.59761
Feature: 5, Score: 4.65989
Feature: 6, Score: -0.00560
Feature: 7, Score: -1.47938
Feature: 8, Score: 0.04484
Feature: 9, Score: -0.07000
Feature: 10, Score: 0.02003
Feature: 11, Score: -0.50303<jupyter_text>SCAD - Diabetes DatasetSCAD's variable selection for the diabetes dataset provides strange results that do not align with most of the other regressors. SCAD identifies BMI, BP, and S5 as having positive correlations, which is similar to other models. However, SCAD does not reduce any variables to 0, and instead provides extremely negative coefficients to most of the other x variables, in contrast to other regressors such as Square Root Lasso, which reduce most of these variables to 0.<jupyter_code>lam = 0.5
a = 1.001
p = diabetesXTrain.shape[1]
b0 = np.random.normal(1,1,p)
y = diabetesYTrain
X = diabetesXTrain
y = y.reshape((y.shape[0]), 1)
diabetesSCAD = scipy.optimize.minimize(scad, b0 , method='L-BFGS-B', jac=dscad,options={'gtol': 1e-8, 'maxiter': 50000,'maxls': 25,'disp': True})
diabetesSCAD.x
importances = diabetesSCAD.x
for i,v in enumerate(importances):
  print('Feature: %0d, Score: %.5f' % (i,v))
plt.bar([x for x in range(len(importances))], importances)
plt.xticks(list(range(0,10)), diabetes_xticks)
plt.title("%s Regression Variable Importance for %s Dataset" %("SCAD", "Diabetes"))
plt.xlabel('Variable Index',fontsize=16)
plt.ylabel('Coefficient Magnitude',fontsize=16)
plt.show()<jupyter_output>Feature: 0, Score: 213.69666
Feature: 1, Score: -288.79951
Feature: 2, Score: 363.87068
Feature: 3, Score: 407.88279
Feature: 4, Score: -16.01311
Feature: 5, Score: -144.52672
Feature: 6, Score: -385.47365
Feature: 7, Score: 81.63260
Feature: 8, Score: 383.55577
Feature: 9, Score: 113.10442<jupyter_text>SCAD - Synthetic DatasetSCAD performs some puzzling variable selection on the synthetically generated data.It identifies X1 and X5 as having strong positive correlations, and X3 as having a negative correlation. However, it leaves X1 with a stronger coefficient than X5 even though the true relationship is reversed. It also does not reduce X2 or X4 to their true values of 0.<jupyter_code>betaDF

lam = 0.5
a = 1.001
p = synthXTrain.shape[1]
b0 = np.random.normal(1,1,p)
y = synthYTrain
X = synthXTrain
y = y.reshape((y.shape[0]), 1)
syntheticSCAD = scipy.optimize.minimize(scad, b0 , method='L-BFGS-B', jac=dscad,options={'gtol': 1e-8, 'maxiter': 50000,'maxls': 25,'disp': True})
syntheticSCAD.x
importances = syntheticSCAD.x
for i,v in enumerate(importances):
  print('Feature: %0d, Score: %.5f' % (i,v))
plt.bar([x for x in range(len(importances))], importances)
plt.xticks(list(range(0,6)), synth_xticks)
plt.title("%s Regression Variable Importance for %s Dataset" %("SCAD", "Synthetic"))
plt.xlabel('Variable Index',fontsize=16)
plt.ylabel('Coefficient Magnitude',fontsize=16)
plt.show()<jupyter_output>Feature: 0, Score: 4.08512
Feature: 1, Score: 1.10038
Feature: 2, Score: -3.00377
Feature: 3, Score: -1.25970
Feature: 4, Score: 2.64774
Feature: 5, Score: 0.03298<jupyter_text>[Square Root Lasso](https://www.jstor.org/stable/23076172?seq=1)I will use the [StatsModels implementation of Square Root Lasso](https://www.statsmodels.org/devel/generated/statsmodels.regression.linear_model.OLS.fit_regularized.html) to demonstrate this regression.I will not demonstrate any hyperparameter choices for this model. Instead of calculating MSE or R2, we'll just examine the variable selection of this model.<jupyter_code>#! pip install --upgrade Cython
#! pip install --upgrade git+https://github.com/statsmodels/statsmodels
import statsmodels.api as sm<jupyter_output>/usr/local/lib/python3.7/dist-packages/statsmodels/tools/_testing.py:19: FutureWarning: pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead.
  import pandas.util.testing as tm<jupyter_text>Square Root Lasso - Boston Housing DatasetSquare Root Lasso selects similar variables to the other models above. It chooses  of Rooms and CHAS as the most correlated X features, and assigns NOX a negative coefficient.<jupyter_code>sqrLBoston = sm.OLS(bostonYTrain,bostonXTrain)
bostonresult = sqrLBoston.fit_regularized(method='sqrt_lasso', alpha=0.5, l1_wt = 1.0)
importances = bostonresult.params
for i,v in enumerate(importances):
  print('Feature: %0d, Score: %.5f' % (i,v))
plt.bar([x for x in range(len(importances))], importances)
plt.xticks(list(range(0,12)), boston_xticks)
plt.title("%s Regression Variable Importance for %s Dataset" %("Square Root Lasso", "Boston Housing"))
plt.xlabel('Variable Index',fontsize=16)
plt.ylabel('Coefficient Magnitude',fontsize=16)
plt.show()<jupyter_output>Feature: 0, Score: -0.09428
Feature: 1, Score: 0.04247
Feature: 2, Score: -0.06469
Feature: 3, Score: 2.49946
Feature: 4, Score: -2.08656
Feature: 5, Score: 5.39102
Feature: 6, Score: 0.00019
Feature: 7, Score: -0.86776
Feature: 8, Score: 0.04230
Feature: 9, Score: -0.37381
Feature: 10, Score: 0.01619
Feature: 11, Score: -0.47415<jupyter_text>Square Root Lasso - Diabetes DatasetLike the other models, Square Root Lasso assigns BMI and Blood Pressure high coefficients. It also assigns s5 a high coefficient, not seen in the other models.<jupyter_code>sqrLDiabetes = sm.OLS(diabetesYTrain,diabetesXTrain)
diabetesresult = sqrLDiabetes.fit_regularized(method='sqrt_lasso', alpha=0.5)
importances = diabetesresult.params
for i,v in enumerate(importances):
  print('Feature: %0d, Score: %.5f' % (i,v))
plt.bar([x for x in range(len(importances))], importances)
plt.xticks(list(range(0,10)), diabetes_xticks)
plt.title("%s Regression Variable Importance for %s Dataset" %("Square Root Lasso", "Diabetes"))
plt.xlabel('Variable Index',fontsize=16)
plt.ylabel('Coefficient Magnitude',fontsize=16)
plt.show()<jupyter_output>Feature: 0, Score: 117.60309
Feature: 1, Score: -29.11350
Feature: 2, Score: 361.07091
Feature: 3, Score: 383.94198
Feature: 4, Score: -0.00009
Feature: 5, Score: 0.00000
Feature: 6, Score: -201.16819
Feature: 7, Score: 0.00090
Feature: 8, Score: 493.31043
Feature: 9, Score: 0.00419<jupyter_text>Square Root Lasso - Synthetic DatasetSquare Root Lasso performs very well on the synthetic dataset, identifying X1 and X5 as positive correlations and giving x5 a higher coefficient. X3 is correctly identified as approximately -3 while the other variables are all reduced to close to zero. Overall this method performs excellently for variable selection on this dataset.<jupyter_code>sqrLSynthetic = sm.OLS(synthYTrain,synthXTrain)
syntheticresult = sqrLSynthetic.fit_regularized(method='sqrt_lasso', alpha=0.5)
betaDF
importances = syntheticresult.params
for i,v in enumerate(importances):
  print('Feature: %0d, Score: %.5f' % (i,v))
plt.bar([x for x in range(len(importances))], importances)
plt.xticks(list(range(0,6)), synth_xticks)
plt.title("%s Regression Variable Importance for %s Dataset" %("Square Root Lasso", "Diabetes"))
plt.xlabel('Variable Index',fontsize=16)
plt.ylabel('Coefficient Magnitude',fontsize=16)
plt.show()<jupyter_output>Feature: 0, Score: 2.37764
Feature: 1, Score: 0.27677
Feature: 2, Score: -3.15225
Feature: 3, Score: 0.16484
Feature: 4, Score: 4.31172
Feature: 5, Score: -0.03129<jupyter_text>Final ComparisonsNow that I've compared the regressors' hyperparameters, I'll compare the different regressors accuaracy and variable selection.I'll create graphs comparing the regressors' importancies and error on each dataset. I'll use the version of each regressor that we had the lowest error on each dataset above.<jupyter_code>column_names = ["Model", "Mean Squared Error", "Test R2"]
model_names = ['Ridge', 'Lasso', "Elastic Net"]<jupyter_output><empty_output><jupyter_text>Boston Housing DatasetAll of the models identified that  of rooms is positively correlated with housing price. Models also identified proximity to the river as being a positive trait, while Distance from Employment Centers, Student/Teacher Ratio, and Portion of Poor Residents (LStat) were all negatively correlated with housing price.Ridge Regression had the lowest MSE and the best R2 on the test data.Interestingly, the Square Root Lasso variable selection mirrors the Ridge variable selection, where I expected it to mirror the Lasso selection. I'm not sure if this points to an error in implementing the regressor, or if I misunderstood how the Square Root Lasso functions.<jupyter_code>boston_mses = [Ridge1BostonMSE, lasso1BostonMSE,ENet20BostonMSE]
boston_r2s = [Ridge1BostonR2, lasso1BostonR2, ENet20BostonR2]
bostonacc = pd.DataFrame()
bostonacc["Model"] = model_names
bostonacc["Mean Squared Error ($)"] = boston_mses
bostonacc["Test R2 "] = boston_r2s
bostonacc.round(decimals = 2)

ticks = list(range(0, 12))
xticks = boston_xticks
plt.plot(np.reshape(Ridge1Boston.coef_, (len(xticks), 1)),alpha=0.7,linestyle='none',marker='*',markersize=8,color='red',label=r'%s; $\alpha = %s $' %("Ridge", "1")) 
plt.plot(lasso1Boston.coef_,alpha=0.5,linestyle='none',marker='d',markersize=7,color='blue',label=r'%s; $\alpha = %s$' %('Lasso', "1"))
plt.plot(ENet20Boston.coef_,alpha=0.5,linestyle='none',marker='o',markersize=7,color='green',label=r'%s; L1 Ratio = %s' %("Elastic Net", '20%')) 
plt.plot(bostonSCAD.x,alpha = 0.5, linestyle = 'none', marker = 'H', markersize = 7, color = "red", label = "SCAD" )
plt.plot(bostonresult.params,alpha=0.5,linestyle='none',marker='v',markersize=7,color='orange',label=r'%s' %('Square Root Lasso'))
plt.xticks(ticks, xticks)
plt.title("Variable Selection for Boston Housing Dataset")
plt.xlabel('Variable Index',fontsize=16)
plt.ylabel('Coefficient Magnitude',fontsize=16)
plt.legend(fontsize=13,loc=1)
plt.show()<jupyter_output><empty_output><jupyter_text>Diabetes DatasetAll of the models performed very poorly on the Diabetes datasets. Their terrible R2 values and high MSE makes me question if examining their variable selections is worthwhile.Lasso had the highest R2 and lowest MSE.Four out of five models selected BMI and Blood pressure as having a high correlation with disease progression. S5, one of the blood tests, was also identified as a positive correlation by Ridge, Lasso, SCAD, and Square Root Lasso.<jupyter_code>diabetes_mses = [Ridge1DiabetesMSE, lasso1DiabetesMSE,ENet80DiabetesMSE]
diabetes_r2s = [Ridge1DiabetesR2, lasso1DiabetesR2, ENet80DiabetesR2]
diabetesacc = pd.DataFrame()
diabetesacc["Model"] = model_names
diabetesacc["Mean Squared Error"] = diabetes_mses
diabetesacc["Test R2 "] = diabetes_r2s
diabetesacc.round(decimals = 2)

ticks = list(range(0, 10))
xticks = diabetes_xticks
plt.plot(Ridge1Diabetes.coef_,alpha=0.5,linestyle='none',marker='*',markersize=8,color='red',label=r'%s; $\alpha = %s $' %("Ridge", "1")) 
plt.plot(lasso1Diabetes.coef_,alpha=0.5,linestyle='none',marker='d',markersize=7,color='blue',label=r'%s; $\alpha = %s$' %('Lasso', "1"))
plt.plot(ENet80Diabetes.coef_,alpha=0.5,linestyle='none',marker='o',markersize=7,color='green',label=r'%s; L1 Ratio = %s' %("Elastic Net", '80%')) 
plt.plot(diabetesresult.params,alpha=0.5,linestyle='none',marker='v',markersize=7,color='orange',label=r'%s' %('Square Root Lasso'))
plt.plot(diabetesSCAD.x,alpha = 0.5, linestyle = 'none', marker = 'H', markersize = 7, color = "red", label = "SCAD" )
plt.xticks(ticks, xticks)
plt.title("Variable Selection for Diabetes Dataset")
plt.xlabel('Variable Index',fontsize=16)
plt.ylabel('Coefficient Magnitude',fontsize=16)
plt.legend(fontsize=13,loc=0)
plt.show()<jupyter_output><empty_output><jupyter_text>Synthetic DataThe models identified similar trends in the synthetic variables, but had varying degrees of success at locating the ground truth. Square Root Lasso was one of the most successful approaches, successfully reducing X4 and X6 to zero, X2 to near zero, and identifying the correct correlations with X1, X3, and X5.<jupyter_code>synth_mses = [Ridge1SyntheticMSE, lasso1SyntheticMSE,ENet20SyntheticMSE]
synth_r2s = [Ridge10SyntheticR2, lasso1SyntheticR2, ENet80SyntheticR2]
synthacc = pd.DataFrame()
synthacc["Model"] = model_names
synthacc["Mean Squared Error"] = synth_mses
synthacc["Test R2 "] = synth_r2s
synthacc.round(decimals = 2)
betaDF

ticks = list(range(0, 6))
xticks = synth_xticks
plt.plot(np.reshape(Ridge10Synthetic.coef_, (len(xticks), 1)),alpha=0.7,linestyle='none',marker='*',markersize=8,color='red',label=r'%s; $\alpha = %s $' %("Ridge", "1")) 
plt.plot(lasso1Synthetic.coef_,alpha=0.5,linestyle='none',marker='d',markersize=7,color='blue',label=r'%s; $\alpha = %s$' %('Lasso', "1"))
plt.plot(ENet20Synthetic.coef_,alpha=0.5,linestyle='none',marker='o',markersize=7,color='green',label=r'%s; L1 Ratio = %s' %("Elastic Net", '20%')) 
plt.plot(syntheticresult.params,alpha=0.5,linestyle='none',marker='v',markersize=7,color='orange',label=r'%s' %('Square Root Lasso'))
plt.plot(syntheticSCAD.x,alpha = 0.5, linestyle = 'none', marker = 'H', markersize = 7, color = "red", label = "SCAD" )
plt.plot(beta,alpha = 0.3, linestyle = 'none', marker = 's', markersize = 7, color = "black", label = "True Values" )
plt.xticks(ticks, xticks)
plt.title("Variable Selection for Synthetic Dataset")
plt.xlabel('Variable Index',fontsize=16)
plt.ylabel('Coefficient Magnitude',fontsize=16)
plt.legend(fontsize=10,loc=0)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11377026956567078, 'fraction_numerical': 0.04103446477918557, 'mean_word_length': 6.118018185192058, 'pattern_counts': {'":': 10, '<': 120, '<?xml version=': 0, '>': 114, 'https://': 16, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '28', 'dir': 'jupyter-structured-clean-dedup', 'id': '313944', 'license': 'MIT', 'n_tokens_mistral': 14313, 'n_tokens_neox': 12743, 'n_words': 3765, 'path': 'AAMLP2.ipynb', 'repo_name': 'NewFigton/newfigton.github.io', 'size': '70324'}