<jupyter_start><jupyter_text>Start with pure linear regression<jupyter_code>three = pickle.load(open("../data/processed/data2.pickle","rb"))<jupyter_output><empty_output><jupyter_text>Drop all the data with incomplete FA information.<jupyter_code>print(three.shape)<jupyter_output>(2260, 52)<jupyter_text>Before we proceed, we can eliminate some of the columns as they can be defined in terms of the other. $FG = 2P+3P$$TRB = ORB+DRB$$TRB\% = ORB\% + DRB\%$$WS = OWS+DWS$$BPM = OBPM + DBPM$$WS/48 = \frac{WS}{48 min}$We can eliminate at least 6 columns to make things easier.<jupyter_code>columns=['Salary', 'Year', 'Pos', 'Age', 'G', 'GS', 'MP','FGA', 'FG%', '3P', '3PA', '3P%', '2P', 
         '2PA', '2P%', 'eFG%', 'FT','FTA', 'FT%', 'ORB', 'DRB', 'AST', 'STL', 'BLK', 'TOV', 'PF',
       'PTS', 'PER', 'TS%', '3PAr', 'FTr', 'ORB%', 'DRB%', 'AST%','STL%', 'BLK%', 'TOV%', 'USG%', 'OWS', 
         'DWS', 'OBPM','DBPM', 'VORP', 'FA']
df_sm2 = three[columns]
df_sm2.columns = df_sm2.columns.str.replace("%","PCT").str.upper()
df_sm2.columns = df_sm2.columns.str.replace("2","_2").str.replace("3","_3")
del df_sm2["POS"]
del df_sm2["FA"]
df_sm2.columns<jupyter_output><empty_output><jupyter_text>Use Statsmodel<jupyter_code>plt.hist(df_sm2['SALARY']);<jupyter_output><empty_output><jupyter_text>I tried different transforms out to make the target variable more normal. The best one was the cube root.<jupyter_code>plt.hist(np.cbrt(df_sm2['SALARY']));
df_sm2.SALARY = np.cbrt(df_sm2.SALARY)
y, X= patsy.dmatrices('SALARY ~ G + GS + MP + OBPM + FTA + FT + FGA + _3P+_3PA+_2P + _2PA + ORB+ AST + TOV + STL + DRB + PTS + PER + OWS + DWS + VORP + USGPCT+FGPCT+_3PPCT+_2PPCT+EFGPCT+FTPCT+BLK+PF+TSPCT+_3PAR+FTR+ORBPCT+DRBPCT+ASTPCT+STLPCT+BLKPCT+TOVPCT+DBPM', data=df_sm2, return_type = "dataframe")
model = sm.OLS(y,X)

fit = model.fit()

fit.summary()
y, X= patsy.dmatrices('SALARY ~ G + MP + ORB + AST + TOV + PTS + PER + DWS + VORP + USGPCT + PF + FTR + ASTPCT + TOVPCT + DBPM', data=df_sm2, return_type = "dataframe")
model = sm.OLS(y,X)

fit = model.fit()

fit.summary()
plt.figure(figsize=(15,10))

sns.heatmap((df_sm2[['SALARY','G','MP','ORB','AST','TOV','PTS','PER','DWS','VORP','USGPCT','PF','FTR','ASTPCT','TOVPCT','DBPM']]).corr()
            , xticklabels=['SALARY','G','MP','ORB','AST','TOV','PTS','PER','DWS','VORP','USGPCT','PF','FTR','ASTPCT','TOVPCT','DBPM']
            , yticklabels=['SALARY','G','MP','ORB','AST','TOV','PTS','PER','DWS','VORP','USGPCT','PF','FTR','ASTPCT','TOVPCT','DBPM'])

plt.xticks(np.arange(0.5, 16, step=1), rotation=30);<jupyter_output><empty_output><jupyter_text>Starting DataframeCreate a smaller DF we can start with.<jupyter_code># df_sm2 = three[columns] 
df_sm2 = three[['Salary','Year','G','MP','ORB','AST','TOV','PTS','PER','DWS','VORP','USG%','PF','FTr','AST%','TOV%','DBPM']]
# df_sm2 = three[['Salary','Year','G','ORB','AST','TOV','PTS','DWS','VORP','USG%','PF','AST%','DBPM']]
df_sm2.columns = df_sm2.columns.str.replace("%","PCT").str.upper()
df_sm2.columns = df_sm2.columns.str.replace("2","_2").str.replace("3","_3")
# df_sm2 = df_sm2[df_sm2.PER > 0]
# df_sm2 = df_sm2[df_sm2.PF < 6]
# df_sm2 = df_sm2[(df_sm2["_2PPCT"] < 1.0) & (df_sm2["_2PPCT"] > 0)]
# df_sm2 = df_sm2[(df_sm2["_3PPCT"] < 1.0) & (df_sm2["_3PPCT"] > 0)]
# df_sm2 = df_sm2[(df_sm2["FGPCT"] < 1.0) & (df_sm2["FGPCT"] > 0)]
# df_sm2 = df_sm2[(df_sm2["FTPCT"] < 1.0) & (df_sm2["FTPCT"] > 0)]
# df_sm2 = df_sm2[df_sm2["G"] > 0]
sns.regplot(df_sm2.G,np.cbrt(df_sm2.SALARY));
pp = sns.pairplot(data=df_sm2,
                  y_vars=['SALARY'],
                  x_vars=['AGE', 'G', 'GS', 'MP'])
pp2 = sns.pairplot(data=df_sm2,
                  y_vars=['SALARY'],
                  x_vars=['FGA', 'FGPCT', '_3P', '_3PA', '_3PPCT'])
pp3 = sns.pairplot(data=df_sm2,
                  y_vars=['SALARY'],
                  x_vars=['_2P', '_2PA', '_2PPCT', 'EFGPCT', 'FT'])
pp4 = sns.pairplot(data=df_sm2,
                  y_vars=['SALARY'],
                  x_vars=['FTA', 'FTPCT', 'ORB', 'DRB', 'AST'])

pp5 = sns.pairplot(data=df_sm2,
                  y_vars=['SALARY'],
                  x_vars=['STL', 'BLK', 'TOV', 'PF','PTS'])

pp6 = sns.pairplot(data=df_sm2,
                  y_vars=['SALARY'],
                  x_vars=['PER', 'TSPCT', '_3PAR', 'FTR', 'ORBPCT'])

pp7 = sns.pairplot(data=df_sm2,
                  y_vars=['SALARY'],
                  x_vars=['DRBPCT', 'ASTPCT','STLPCT', 'BLKPCT', 'TOVPCT'])

pp8 = sns.pairplot(data=df_sm2,
                  y_vars=['SALARY'],
                  x_vars=['USGPCT', 'OWS','DWS', 'OBPM','DBPM', 'VORP'])
def lr_test(dataset, target):
    test = dataset[(dataset.YEAR == 2018)]
    train = dataset[dataset.YEAR != 2018]
    test = test.drop('YEAR', axis=1)
    train = train.drop('YEAR', axis=1)
    
    X_test, y_test = pd.get_dummies(test[test.columns[1:]]), test[target]
    X_train, y_train = pd.get_dummies(train[train.columns[1:]]), np.cbrt(train[target])

#     print(X_train.columns)
#     print(X_test.columns)
#     print(X_train.shape, y_train.shape,X_test.shape, y_test.shape)
    
    # Linear Regression
    lr = LinearRegression()
    lr.fit(X_train, y_train)
    
    r2 = lr.score(X_train, y_train)
#     r2_poly = lm_poly.score(X_train_poly, y_train)
    N = X_train.shape[0]
    p = X_train.shape[1]
    
    adjusted_r2 = 1-((1-r2**2)*(N-1)/(N-p-1))
    # Print out the R^2 for the model against the full dataset
    print(f'Linear Regression R2 score: {r2:.6f}')
    print(f'Adjusted R2 Score: {adjusted_r2: .6f}')
#     print(f'\nPolynomial Regression R2 score: {r2_poly:.6f}')
#     print(f'Polynomial Adjusted R2 Score: {adjusted_r2_poly: .6f}')
    
    
    
    # Predictions
    y_pred = lr.predict(X_test)
#     y_pred_poly = lm_poly.predict(X_test_poly)
    
    
    
    print("\nLR Test Data R2 score: ", r2_score(y_test,y_pred**3))
    print("LR Test MAE: ", mean_absolute_error(y_test,y_pred**3))
    print("Median 2018 NBA Salary: ", y_test.median())
    print("Mean 2018 NBA Salary: ", y_test.mean())
    print("Median predicted Salary: ", np.median(y_pred**3))
    
#     plt.figure(figsize=(20,10))
#     plt.subplot(1,2,1)
    
#     res = y_test - y_pred**3
    
#     plt.scatter(y_pred**3,res)
#     plt.title("Residual Plot of 2018 NBA SALARY Prediction")
#     plt.xlabel("Prediction (10 Million $)")
#     plt.ylabel("Residuals (10 Million $)")
    
#     plt.subplot(1,2,2)
    
#     plt.scatter(y_pred**3, y_test)
#     plt.title("Salary: Actual VS Prediction")
#     plt.xlabel("Prediction (10 Million $)")
#     plt.ylabel("Actual (10 Million $)")
    
    return y_pred, lr#, lm_poly<jupyter_output><empty_output><jupyter_text>Polynomial Terms - ^2, sqrt, interaction terms, etc.<jupyter_code>y_pred, lr1 = lr_test(df_sm2, "SALARY")<jupyter_output>Linear Regression R2 score: 0.505045
Adjusted R2 Score:  0.249232

LR Test Data R2 score:  0.30613163743871474
LR Test MAE:  5091888.859591299
Median 2018 NBA Salary:  7000000.0
Mean 2018 NBA Salary:  9558095.951515151
Median predicted Salary:  4667143.632683283<jupyter_text>Add the square terms:<jupyter_code>features = ['G','ORB','AST','TOV','PTS','DWS','VORP','USGPCT','PF','ASTPCT','DBPM']
df_sm2_2 = df_sm2.copy()
for fea in features:
    df_sm2_2[fea+"2"] = df_sm2[fea]**2
y_pred, lr2 = lr_test(df_sm2_2,"SALARY")<jupyter_output>Linear Regression R2 score: 0.517131
Adjusted R2 Score:  0.257415

LR Test Data R2 score:  0.3123044996761797
LR Test MAE:  5070143.299831751
Median 2018 NBA Salary:  7000000.0
Mean 2018 NBA Salary:  9558095.951515151
Median predicted Salary:  4904633.426947447<jupyter_text>Then the Interaction terms:<jupyter_code>df_sm2_3 = df_sm2_2.copy()
for item in itertools.combinations(features,2):
    df_sm2_2[item[0]+item[1]] = df_sm2[item[0]]*df_sm2[item[1]]
y_pred, lr3 = lr_test(df_sm2_3,"SALARY")
df_sm2_3 = df_sm2.copy()
df_sm2_2["PTSG"] = df_sm2["PTS"]*df_sm2["G"]
y_pred, lr2 = lr_test(df_sm2_2, "SALARY")
df_sm2_3 = df_sm2_2.copy()
df_sm2_3["PF2"] = df_sm2["PF"]**2
y_pred, lr3 = lr_test(df_sm2_3, "SALARY")
df_sm2_4 = df_sm2_3.copy()
df_sm2_4["PF2"] = df_sm2["PF"]**2
y_pred, lr4 = lr_test(df_sm2_4, "SALARY")
df_sm2_5 = df_sm2_4.copy()
df_sm2_5["DRBPCT2"] = df_sm2["DRBPCT"]**2
y_pred, lr5 = lr_test(df_sm2_5, "SALARY")
df_sm2_6 = df_sm2_5.copy()
df_sm2_6["G2"] = df_sm2["G"]**2
y_pred, lr6 = lr_test(df_sm2_6, "SALARY")
df_sm2_7 = df_sm2_6.copy()
df_sm2_7["PTSG"] = df_sm2["PTS"]*df_sm2["G"]
y_pred, lr7 = lr_test(df_sm2_7, "SALARY")
df_sm2_8 = df_sm2_7.copy()
df_sm2_8["ASTG"] = df_sm2["AST"]*df_sm2["G"]
df_sm2_8["TOVG"] = df_sm2["TOV"]*df_sm2["G"]
df_sm2_8["VORPG"] = df_sm2["VORP"]*df_sm2["G"]
y_pred, lr8 = lr_test(df_sm2_8, "SALARY")
pickle.dump(lr6, open("../models/LR.pickle", 'wb'))<jupyter_output><empty_output><jupyter_text>Cross-Validation & adding Ridge Regularization<jupyter_code>def cross_val(dataset):
    
    # Split out the y (keep the Year on both for test/train split)
    X,y = dataset.drop('SALARY',axis=1), dataset[['YEAR','SALARY']]
    
    # Test/Train split
    X,X_test = X[X.YEAR != 2018],X[X.YEAR==2018]
    y,y_test = y[y.YEAR != 2018],y[y.YEAR==2018]
    
#     print(X_test.shape, y_test.shape)
    
    # Drop the Year column on Test data first
    y_test = y_test.drop('YEAR',axis=1)
    X_test = X_test.drop('YEAR',axis=1)
    
    LR_r2s, LM_REG_r2s = [],[]
    LR_MAES, LM_REG_MAES = [],[]
    
    for i in range(2017,2013,-1):
#         print("Val = Year ", str(i))
        years = [i-j for j in range(1,4)]
#         print(years)
        X_val,y_val = X[X.YEAR == i], y[y.YEAR == i]
        X_val = X_val.drop('YEAR', axis=1)
        y_val = y_val.drop('YEAR', axis=1)
#         print(X_val.shape, y_val.shape)
        X_train, y_train = X[X.YEAR.isin(years)], y[y.YEAR.isin(years)]
        X_train = X_train.drop('YEAR', axis=1)
        y_train = y_train.drop('YEAR', axis=1)
#         print(X_train.shape, y_train.shape)
#         print(X_train.columns)
#         print(y_train.columns)
        
        
        # Set dummies
        X_val,X_test,X_train = pd.get_dummies(X_val), pd.get_dummies(X_test), pd.get_dummies(X_train)

         #simple linear regression
        lm = LinearRegression()
        lm_reg = Ridge(alpha=1)

        lm.fit(X_train, y_train)
        LR_r2s.append(lm.score(X_val, y_val))
        y_pred = lm.predict(X_val)
        LR_MAES.append(mean_absolute_error(y_pred,y_val))

        #ridge with feature scaling
        scaler = StandardScaler()
        X_train_scaled = scaler.fit_transform(X_train)
        X_val_scaled = scaler.transform(X_val)

        lm_reg.fit(X_train_scaled, y_train)
        LM_REG_r2s.append(lm_reg.score(X_val_scaled, y_val))
        y_pred_reg = lm_reg.predict(X_val_scaled)
        LM_REG_MAES.append(mean_absolute_error(y_pred_reg,y_val))
    
#     print(LR_r2s)
#     print(LR_MAES)
    mean_MAE_LR = np.mean(LR_MAES)
    mean_MAE_REG = np.mean(LM_REG_MAES)
    print("Average R2 Score of Linear Regression: ", np.mean(LR_r2s))
    print("Average MAE of Linear Regression: ", mean_MAE_LR)
#     print(LM_REG_r2s)
#     print(LM_REG_MAES)
    print("Average R2 Score of Ridge Regression w/ Scaler: ", np.mean(LM_REG_r2s))
    print("Average MAE of Ridge Regression: ", mean_MAE_REG)
    
    # Test on the whole training set
    y = y.drop('YEAR',axis=1)
    X = X.drop('YEAR',axis=1)
    X = pd.get_dummies(X)
    X_te = scaler.transform(X)
    
    if mean_MAE_REG < mean_MAE_LR:
        print("\nRidge Regression performs better")
        pred = lm_reg.predict(X_te)
    else:
        print("\nLinear Regression performs better")
        pred = lm.predict(X)
    
#     print()
    print("MAE on the whole training set: ", mean_absolute_error(pred,y))
    
    X_test_te = scaler.transform(X_test)
    pred_test = lm_reg.predict(X_test_te)
    print("MAE on the test set: ", mean_absolute_error(pred_test,y_test))
    
    res = pred_test-y_test

#     plt.figure(figsize=(20,10))

#     plt.scatter(pred_test,res)
#     plt.title("Residual Plot of 2018 NBA SALARY Prediction")
#     plt.xlabel("Prediction (10 Million $)")
#     plt.ylabel("Residuals (10 Million $)")    

    
    return lm_reg, lm, pred_test, y_test
lm_reg_cv, lm_cv, pred_test, y_test = cross_val(df_sm2_6)
def ridge_test(dataset):
    
    # Split out the y (keep the Year on both for test/train split)
    X,y = dataset.drop('SALARY',axis=1), dataset[['YEAR','SALARY']]
    
    # Test/Train split
    X,X_test = X[X.YEAR != 2018],X[X.YEAR==2018]
    y,y_test = y[y.YEAR != 2018],y[y.YEAR==2018]
    
#     print(X_test.shape, y_test.shape)
    
    # Drop the Year column on Test data first
    y_test = y_test.drop('YEAR',axis=1)
    X_test = X_test.drop('YEAR',axis=1)
    
    for i in range(2017,2013,-1):
#         print("Val = Year ", str(i))
        years = [i-j for j in range(1,4)]
#         print(years)
        X_val,y_val = X[X.YEAR == i], y[y.YEAR == i]
        X_val = X_val.drop('YEAR', axis=1)
        y_val = y_val.drop('YEAR', axis=1)
#         print(X_val.shape, y_val.shape)
        X_train, y_train = X[X.YEAR.isin(years)], y[y.YEAR.isin(years)]
        X_train = X_train.drop('YEAR', axis=1)
        y_train = y_train.drop('YEAR', axis=1)
#         print(X_train.shape, y_train.shape)
#         print(X_train.columns)
#         print(y_train.columns)
        
        
        # Set dummies
        X_val,X_test,X_train = pd.get_dummies(X_val), pd.get_dummies(X_test), pd.get_dummies(X_train)
        values = []
        
        #
        alphalist = 10**(np.linspace(-2,2,200))
        err_vec_val = np.zeros(len(alphalist))
        err_vec_train = np.zeros(len(alphalist))

        for k,curr_alpha in enumerate(alphalist):

            # note the use of a new sklearn utility: Pipeline to pack
            # multiple modeling steps into one fitting process 
            steps = [('standardize', StandardScaler()), 
                     ('ridge', Ridge(alpha = curr_alpha))]

            pipe = Pipeline(steps)
            pipe.fit(X_train.values, y_train)

            val_set_pred = pipe.predict(X_val.values)
            err_vec_val[k] = mean_absolute_error(y_val, val_set_pred)
        
        print("Minimum MAE: ", np.min(err_vec_val))
        print("When alpha = ", alphalist[np.argmin(err_vec_val)])
#             plt.plot(np.log10(alphalist), err_vec_val)
        
#         # ridge with feature scaling
#         scaler = StandardScaler()
#         X_train_scaled = scaler.fit_transform(X_train)
#         X_val_scaled = scaler.transform(X_val)
        
#         # fit and predict
#         lm_reg.fit(X_train_scaled, y_train)
#         LM_REG_r2s.append(lm_reg.score(X_val_scaled, y_val))
#         y_pred_reg = lm_reg.predict(X_val_scaled)
#         LM_REG_MAES.append(mean_absolute_error(y_pred_reg,y_val))
    
# #     print(LR_r2s)
# #     print(LR_MAES)
#     mean_MAE_LR = np.mean(LR_MAES)
#     mean_MAE_REG = np.mean(LM_REG_MAES)
#     print("Average R2 Score of Linear Regression: ", np.mean(LR_r2s))
#     print("Average MAE of Linear Regression: ", mean_MAE_LR)
# #     print(LM_REG_r2s)
# #     print(LM_REG_MAES)
#     print("Average R2 Score of Ridge Regression w/ Scaler: ", np.mean(LM_REG_r2s))
#     print("Average MAE of Ridge Regression: ", mean_MAE_REG)
    
#     # Test on the whole training set
#     y = y.drop('YEAR',axis=1)
#     X = X.drop('YEAR',axis=1)
#     X = pd.get_dummies(X)
#     X_te = scaler.transform(X)
    
#     if mean_MAE_REG < mean_MAE_LR:
#         print("\nRidge Regression performs better")
#         pred = lm_reg.predict(X_te)
#     else:
#         print("\nLinear Regression performs better")
#         pred = lm.predict(X)
    
# #     print()
#     print("MAE on the whole training set: ", mean_absolute_error(pred,y))
    
#     X_test_te = scaler.transform(X_test)
#     pred_test = lm_reg.predict(X_test_te)
#     print("MAE on the test set: ", mean_absolute_error(pred_test,y_test))
    
#     res = pred_test-y_test

#     plt.figure(figsize=(20,10))

#     plt.scatter(pred_test,res)
#     plt.title("Residual Plot of 2018 NBA SALARY Prediction")
#     plt.xlabel("Prediction (10 Million $)")
#     plt.ylabel("Residuals (10 Million $)")    

    
    return steps, pred_test, y_test
lm_reg_cv, pred_test, y_test = ridge_test(df_sm2_6)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17144438976377951, 'fraction_numerical': 0.04023129921259842, 'mean_word_length': 3.4859271523178808, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 44, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, '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': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '22796', 'license': 'MIT', 'n_tokens_mistral': 7225, 'n_tokens_neox': 6852, 'n_words': 1105, 'path': 'notebooks/Modeling-1yrdata.ipynb', 'repo_name': 'youngjeong46/nba-fa-salary-pred', 'size': '16256'}