<jupyter_start><jupyter_text>Chapter 7 Question 7Further exploration of the `Wage` dataset, using other features, and other non-linear fitting methods.Those other methods:- Polynomial regression- Step functions- Regression splines- Smoothing splines- Local regression- GAMs (the general framework for the above)<jupyter_code>import statsmodels.api as sm
import seaborn as sns
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

import sklearn.linear_model
import sklearn.metrics
import sklearn.pipeline

import statsmodels.stats.anova

sns.set(style="whitegrid")
wage_data = sm.datasets.get_rdataset("Wage", "ISLR").data
wage_data.head()
wage_data.wage.describe()
sns.pairplot(wage_data)<jupyter_output><empty_output><jupyter_text>We have lots of categorical variables here. How many groups does each have?<jupyter_code>for col in ["maritl", "race", "education", "region", "jobclass", "health", "health_ins"]:
    print(col)
    print(wage_data[col].value_counts())
    print()<jupyter_output>maritl
2. Married          2074
1. Never Married     648
4. Divorced          204
5. Separated          55
3. Widowed            19
Name: maritl, dtype: int64

race
1. White    2480
2. Black     293
3. Asian     190
4. Other      37
Name: race, dtype: int64

education
2. HS Grad            971
4. College Grad       685
3. Some College       650
5. Advanced Degree    426
1. < HS Grad          268
Name: education, dtype: int64

region
2. Middle Atlantic    3000
Name: region, dtype: int64

jobclass
1. Industrial     1544
2. Information    1456
Name: jobclass, dtype: int64

health
2. >=Very Good    2142
1. <=Good          858
Name: health, dtype: int64

health_ins
1. Yes    2083
2. No      917
Name: health_ins, dtype: int64<jupyter_text>All these can reasonably be one-hot encoded (i.e. dummy variables created), except "region" which tells us nothing at all<jupyter_code>wage_data = wage_data.drop(columns=["region"])
for col in ["maritl", "race", "education", "jobclass", "health", "health_ins"]:
    dummies = pd.get_dummies(wage_data[col], prefix=col, drop_first=True)
    wage_data = pd.concat([wage_data, dummies], axis=1).drop(columns=[col])
wage_data.head()
X = wage_data.drop(columns=["logwage", "wage"])
y = wage_data.wage
for col in X:    
    if len(X[col].unique()) != 2:
        sns.scatterplot(x=col, y="wage", data=wage_data)
    else:
        sns.violinplot(x=col, y="wage", data=wage_data,orient="v")
    plt.show()
# Route one - a basic OLS regressor with all features
model = sklearn.linear_model.LinearRegression()
scores = sklearn.model_selection.cross_val_score(model, X,y, scoring="neg_mean_squared_error", cv=5)
score = -np.mean(scores)
print(np.sqrt(score))
# Route two - a basic Lasso regressor
lasso = sklearn.linear_model.LassoCV(cv=5)
lasso.fit(X, y)
y_pred = lasso.predict(X)
score = sklearn.metrics.mean_squared_error(y, y_pred)
print(np.sqrt(score))
X.columns[np.where(abs(lasso.coef_) < 1e-7)]<jupyter_output>33.959015003765494<jupyter_text>- Polynomial regression- Step functions- Regression splines- Smoothing splines- Local regression- GAMs (the general framework for the above) I'm using dummy classes for most of the variables, so I don't see how a non-linear function can improve anything! If x is 0, 1, then f(x) can only take two values!Let's try to fit the age and year variables better<jupyter_code>score_list = []
for degree in range(1,20):
    
    pipe = sklearn.pipeline.make_pipeline(sklearn.preprocessing.PolynomialFeatures(degree), sklearn.linear_model.LinearRegression())
    scores = sklearn.model_selection.cross_val_score(pipe, X,y, scoring="neg_mean_squared_error", cv=5)
    score = -np.mean(scores)
    print(f"degree: {degree}, score: {score:.2f}")
    score_list.append(score)                                       
results = []
for degree in range(1,20):
    features =  X[["year", "age"]]
    poly_features = sklearn.preprocessing.PolynomialFeatures(degree=degree).fit_transform(features)
    
    poly_features = np.hstack([X.drop(columns=["year","age"]).to_numpy(), poly_features])
    
    model = sm.OLS(y, poly_features)
    result =  model.fit()
    results.append(result)
    
statsmodels.stats.anova.anova_lm(*results)
# The above suggests degree 3 is useful

degree = 3 
poly = sklearn.preprocessing.PolynomialFeatures(degree=degree)

names = poly.fit(features).get_feature_names()
names = [x.replace("x0", "year").replace("x1", "age") for x in names]
poly_features = poly.transform(features)
poly_features = pd.DataFrame(poly_features, columns=names)
old_features = X.drop(columns=["year","age"])
poly_features = pd.concat([old_features.reset_index(drop=True) ,poly_features.reset_index(drop=True)], axis=1)


model = sklearn.linear_model.LinearRegression()
scores = sklearn.model_selection.cross_val_score(model, poly_features,y, scoring="neg_mean_squared_error", cv=5)
score = -np.mean(scores)
print(np.sqrt(score))

# Route two - a basic Lasso regressor
lasso = sklearn.linear_model.LassoLarsCV(cv=5)
lasso.fit(poly_features, y)
y_pred = lasso.predict(poly_features)
score = sklearn.metrics.mean_squared_error(y, y_pred)
print(np.sqrt(score))
print(lasso.coef_)

print(poly_features.columns[np.where(abs(lasso.coef_) < 1e-7)])
for x,y in zip(poly_features.columns, lasso.coef_):
    print(f"{x}:{y:.5f}")
y = wage_data.wage<jupyter_output>maritl_2. Married:13.45397
maritl_3. Widowed:0.54085
maritl_4. Divorced:0.17927
maritl_5. Separated:7.29461
race_2. Black:-4.71553
race_3. Asian:-2.79674
race_4. Other:-5.79010
education_2. HS Grad:7.54175
education_3. Some College:18.02028
education_4. College Grad:30.53287
education_5. Advanced Degree:53.11785
jobclass_2. Information:3.49839
health_2. >=Very Good:6.25692
health_ins_2. No:-16.40808
1:0.00000
year:1.31583
age:3.41879
year^2:0.00000
year age:0.00000
age^2:0.00000
year^3:0.00000
year^2 age:0.00000
year age^2:-0.00002
age^3:0.00013
{'fraction_non_alphanumeric': 0.12321549966009518, 'fraction_numerical': 0.05387491502379334, 'mean_word_length': 4.631578947368421, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, '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': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '639203', 'license': 'MIT', 'n_tokens_mistral': 2294, 'n_tokens_neox': 2088, 'n_words': 504, 'path': 'Chapter7/question7.ipynb', 'repo_name': 'wllgrnt/islr-python-examples', 'size': '5884'}