<jupyter_start><jupyter_text>Problems**1.Predicting Boston Housing Prices.** The file `BostonHousing.csv` contains information collected by the US Bureau of the Census concerning housing in the area of Boston, Massachusetts. The dataset includes information on 506 census housing tracts in the Boston area. The goal is to predict the median house price in new tracts based on information such as crime rate, pollution, and number of rooms. The dataset contains 13 predictors, and the outcome variable is the median house price (MEDV). The following describes each of the predictors and the outcomevariable:    CRIM: Per capita crime rate by town    ZN: Proportion of residential land zonedfor lots over 25,000 ft squared    INDUS: Proportion of nonretail business acres per town    CHAS: Charles River dummy variable (=1 if tract bounds river; =0 otherwise)    NOX: Nitric oxide concentration (parts per 10 million)    RM: Average number of rooms per dwelling    AGE: Proportion of owner-occupied units built prior to 1940    DIS: Weighted distances to five Boston employment centers    RAD: Index of accessibility to radial highways    TAX: Full-value property-tax rate per $10,000    PTRATIO: Pupil/teacher ratio by town    LSTAT: Percentage lower status of the population    MEDV: Median value of owner-occupied**a.** Why should the data be partitioned into training and validation sets? What will the training set be used for? What will the validation set be used for?When we use the same data both to develop the model and to assess its performance, we introduce an "optimism" bias. This is because when we choose the model that works best with the data, this model's superior performance comes from two sources:- A superior model- Chance aspects of the data that happen to match the chosen model better than they match other modelsThe latter is a particularly serious problem with techniques (such as trees and neural nets) that do not impose linear or other structure on the data, and thus end up overfitting it.To address the overfitting problem, we simply divide (partition) our data and develop our model using only one of the partitions. After we have a model, we try it out on another partition and see how it performs, which we can measure in several ways. In a classification model, we can count the proportion of held-back records that were misclassified. In a prediction model, we can measure the residuals (prediction errors) between the predicted values and the actual values. This evaluation approach in effect mimics the deployment scenario, where our model is applied to data that it hasn't "seen."**Training Partition**: the training partition, typically the largest partition, contains the data used to build the various models we are examining. The same training partition is generally used to develop multiple models.**Validation Partition**: the validation partition (sometimes called the test partition) is used to assess the predictive performance of each model so that you can compare models and choose the best one.**b.** Fit a multiple linear regression model to the median house price (MEDV) as a function of CRIM, CHAS, and RM. Write the equation for predicting the median house price from the predictors in the model.For the above predictors, we will estimate the coefficients of the regression formula from the data using a method called *ordinary least squares* (OLS). This method finds values $\hat{\beta}_0$, $\hat{\beta}_1$, $\hat{\beta}_2$, ..., $\hat{\beta}_p$, where $\hat{\beta_i}$ is the parameter for the predictor *i*, that minimize the sum of squared deviations between the actual outcome values ($Y$) and their predicted values based on that model ($\hat{Y}$).To predict the value of the outcome variable for a record with predictor values $x_1$, $x_2$, ..., $x_p$, we use the equation:            $\hat{Y} = \hat{\beta}_0 + \hat{\beta}_1 x_1 + \hat{\beta}_2 x_2 + ... + \hat{\beta}_p x_p$    One very rough rule of thumb is to have a number of records *n* larger than *5(p + 2)*, where *p* is the number of predictors. Since we have 3 predictors, we need at least *5(3 + 2) = 25* records. We have 506 records, so this is not a concern.Let's estimate the predictors (using all data):<jupyter_code>from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import MinMaxScaler

import seaborn as sns
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.formula.api as sm

from dmutils import adjusted_r2_score, AIC_score
from dmutils import exhaustive_search, backward_elimination
from dmutils import forward_selection, stepwise_selection
from dmutils import regression_summary
from dmutils import lift_chart, gains_chart
# load housing data
housing_df = pd.read_csv("../datasets/BostonHousing.csv")
housing_df
# select predictors
predictors = ["CRIM", "CHAS", "RM"]
outcome = "MEDV"

# partition data
X = housing_df[predictors]
y = housing_df[outcome]

housing_lmr = LinearRegression()
housing_lmr.fit(X, y)

# print coefficients
print(pd.DataFrame({'Predictor': X.columns,
                    'coefficient': housing_lmr.coef_}))
# print intercept
print('\nIntercept: ', housing_lmr.intercept_)
# print performance measures (training data)
regression_summary(y, housing_lmr.predict(X))<jupyter_output>Predictor  coefficient
0      CRIM    -0.260724
1      CHAS     3.763037
2        RM     8.278180

Intercept:  -28.810682506359118

Regression statistics

                      Mean Error (ME) : -0.0000
       Root Mean Squared Error (RMSE) : 6.1453
            Mean Absolute Error (MAE) : 4.1782
          Mean Percentage Error (MPE) : -7.3835
Mean Absolute Percentage Error (MAPE) : 22.9616<jupyter_text>Therefore, we can write the equation for predicting the median house price as:            $\hat{MEDV} = -28.81 - 0.26 x_{CRIM} + 3.76 x_{CHAS} + 8.27 x_{RM}$     **c.** Using the estimated regression model, what median house price is predicted for a tract in the Boston area that does not bound the Charles River, has a crime rate of 0.1, and where the average number of rooms per house is 6?<jupyter_code># the example of a data frame home to be scored using our regression model

new_data = pd.DataFrame({"CHAS": [0], "CRIM": [0.1], "RM": [6]})

print(new_data)
print("\nPredictions: ", housing_lmr.predict(new_data))<jupyter_output>CHAS  CRIM  RM
0     0   0.1   6

Predictions:  [21.23470007]<jupyter_text>Thus the predicted MEDV is \\$21.23 (in \\$1000s) for this house. Since there is no data that is exactly the same as the example given above we are not able to calculate the prediction error (residual). **d.** Reduce the number of predictors:    i. Which predictors are likely to be measuring the same thing among the 13 predictors? Discuss the    relationships among INDUS, NOX, and TAX.    ii. Compute the correlation table for the 12 numerical predictors and search for highly correlated pairs.    These have potential redundancy and can cause multicollinearity. Choose which ones to remove based on this    table.Let's plot the correlation matrix among variables and check for possible multicollinearity (Multicollinearity is the presence of two or more predictors sharing the same linear relationship with the outcome variable).<jupyter_code># simple heatmap of correlations (without values)
corr = housing_df.corr()

# include information about values
# (example demonstrate how to control the size of the plot)
fig, ax = plt.subplots()
fig.set_size_inches(11, 7)
sns.heatmap(corr, annot=True, fmt=".1f", cmap="RdBu", center=0, ax=ax)
plt.show()<jupyter_output><empty_output><jupyter_text>**i.** We can see from the above correlation matrix that the correlation between TAX, INDUS and NOX variables are strong, and positive, correlations (knowing the value of one of these variables gives a hint about the others):- TAX and INDUS: 0.7- TAX and NOX: 0.7- NOX and INDUS: 0.8**ii.** Other correlations are:- DIS and INDUS: -0.7- DIS and TAX: -0.7- DIS and AGE: -0.7- DIS and NOX: -0.7These correlation suggest that the weighted distances to five Boston employment centers are inversely proportional to the proportion of nonretail business acres per town, proportion of owner-occupied units built prior to 1940 (old houses), nitric oxide concentration (pollution, parts per 10 million), and full-value property-tax rate per \\$10,000.Therefore, these would be good candidates to be removed from the predictor lists as a way to reduce/remove multicollinearity problem.     iii. Use three subset selection algorithms: backward, forward, and stepwise) to reduce the remaining    predictors. Compute the validation performance for each of the three selected models. Compare RMSE, MAPE,    and mean error, as well as histograms of the errors. Finally, describe the best model.<jupyter_code># select predictors (all variables)
predictors = ["CRIM", "ZN", "INDUS", "CHAS", "NOX",
              "RM", "AGE", "DIS", "RAD", "TAX", "PTRATIO", "LSTAT"]

# partition data, 40% test set
X = housing_df[predictors]
y = housing_df[outcome]
train_X, valid_X, train_y, valid_y = train_test_split(X, y, 
                                                      test_size=0.4,
                                                      random_state=1)<jupyter_output><empty_output><jupyter_text>- Exhaustive search:<jupyter_code>def train_model(variables):
    model = LinearRegression()
    model.fit(train_X[list(variables)], train_y)
    return model

def score_model(model, variables):
    pred_y = model.predict(train_X[list(variables)])
    # we negate as score is optimized to be as low as possible
    return -adjusted_r2_score(train_y, pred_y, model)

all_variables = train_X.columns
results = exhaustive_search(all_variables, train_model, score_model)

data = []
for result in results:
    model = result["model"]
    variables = list(result["variables"])
    AIC = AIC_score(train_y, model.predict(train_X[variables]), model)

    d = {"n": result["n"], "r2adj": -result["score"], "AIC": AIC}
    d.update({var: var in result["variables"] for var in sorted(all_variables)})
    data.append(d)

pd.DataFrame(data, columns=("n", "r2adj", "AIC") + tuple(sorted(all_variables)))<jupyter_output><empty_output><jupyter_text>- Backward Selection<jupyter_code>def train_model(variables):
    model = LinearRegression()
    model.fit(train_X[variables], train_y)
    return model

def score_model(model, variables):
    return AIC_score(train_y, model.predict(train_X[variables]), model)


all_variables = train_X.columns
best_model, best_variables = backward_elimination(all_variables, train_model, score_model,
                                                  verbose=True)
print(best_variables)
regression_summary(valid_y, best_model.predict(valid_X[best_variables]))
all_residuals = valid_y - best_model.predict(valid_X[best_variables])

# Determine the percentage of datapoints with a residual in [-5, 5] = approx. 75%
print("Percentage of datapoints with a residual in [-5, 5]: ",
      len(all_residuals[(all_residuals > -5) & (all_residuals < 5)]) / len(all_residuals))

pd.DataFrame({"Residuals": all_residuals}).hist(bins=25)
plt.show()<jupyter_output>Percentage of datapoints with a residual in [-5, 5]:  0.7684729064039408<jupyter_text>- Forward Selection<jupyter_code># The initial model is the constant model - this requires special handling
# in train_model and score_model
def train_model(variables):
    if len(variables) == 0:
        return None
    model = LinearRegression()
    model.fit(train_X[variables], train_y)
    return model

def score_model(model, variables):
    if len(variables) == 0:
        return AIC_score(train_y, [train_y.mean()] * len(train_y), model, df=1)
    return AIC_score(train_y, model.predict(train_X[variables]), model)

best_model, best_variables = forward_selection(train_X.columns, train_model,
                                               score_model, verbose=True)
print(best_variables)
regression_summary(valid_y, best_model.predict(valid_X[best_variables]))
all_residuals = valid_y - best_model.predict(valid_X[best_variables])

# Determine the percentage of datapoints with a residual in [-5, 5] = approx. 75%
print("Percentage of datapoints with a residual in [-5, 5]: ",
      len(all_residuals[(all_residuals > -5) & (all_residuals < 5)]) / len(all_residuals))

pd.DataFrame({"Residuals": all_residuals}).hist(bins=25)
plt.show()<jupyter_output>Percentage of datapoints with a residual in [-5, 5]:  0.7684729064039408<jupyter_text>- Stepwise Selection<jupyter_code>best_model, best_variables = stepwise_selection(train_X.columns, train_model,
                                                score_model, verbose=True)
print(best_variables)
regression_summary(valid_y, best_model.predict(valid_X[best_variables]))
all_residuals = valid_y - best_model.predict(valid_X[best_variables])

# Determine the percentage of datapoints with a residual in [-5, 5] = approx. 75%
print("Percentage of datapoints with a residual in [-5, 5]: ",
      len(all_residuals[(all_residuals > -5) & (all_residuals < 5)]) / len(all_residuals))

pd.DataFrame({"Residuals": all_residuals}).hist(bins=25)
plt.show()<jupyter_output>Percentage of datapoints with a residual in [-5, 5]:  0.7684729064039408<jupyter_text>The chosen ten-predictor models from backward, forward and stepwise selection are identical to the best eight-predictor model chosen by the exhaustive search. The variables are:- CHAS;- CRIM;- DIS;- LSTAT;- NOX;- PTRATIO;- RAD;- RM;- TAX;- ZN;In this example, forward selection, backward selection and stepwise regression end up with this same ten-predictor model. This need not be the case with other datasets. Besides that, the histogram of errors are the same as the resultant models are the same for all variable selector methods.Finally, the validation performance of the ten-predictor models is slightly better than the all-predictor model in terms of validation metrics (see below table). This happens because of the multicollinearity effect that regression coefficients suffer - there are variables in the model that have small coefficients relative to the standard deviation of the noise and also exhibit at least moderate correlation with other variables);                       R_adj       AIC      ten-predictor: 0.731108    1803.571740      all-predictor: 0.729557    1807.231231There is an even smaller model, performing only slightly worse, that is preferable from a parsimony standpoint. **2.Predicting Software Reselling Profits.** Tayko Software is a software catalog firm that sells games and educational software. It started out as a software manufacturer and then added third-party titles to its offerings. It recently revised its collection of items in a new catalog, which it mailed out to its customers. This mailing yielded 2000 purchases. Based on these data, Tayko wants to devise a model for predicting the spending amount that a purchasing customer will yield. The file `Tayko.csv` contains information on 2000 purchases. The following table describes the variables to be used in the problem (the Excel file contains additional variables).    FREQ: Number of transactions in the preceding year    LAST_UPDATE: Number of days since last update to customer record    WEB: Whether customer purchased by Web order at least once    GENDER: Male or female    ADDRESS_RES: Whether it is a residential address    ADDRESS_US: Whether it is a US address    SPENDING: Amount spent by customer in (outcome) test mailing ($)**a.** Explore the spending amount by creating a pivot table for the categorical variables and computing the average and standard deviation of spending in each category.<jupyter_code>columns=["Freq", "last_update_days_ago", "Web order", "Gender=male",
         "Address_is_res", "US", "Spending"]
selling_df = pd.read_csv("../datasets/Tayko.csv")
selling_df = selling_df[columns].copy()
selling_df = selling_df.rename(columns={"Freq": "FREQ",
                                        "last_update_days_ago": "LAST_UPDATE",
                                        "Web order": "WEB",
                                        "Gender=male": "GENDER",
                                        "Address_is_res": "ADDRESS_RES",
                                        "US": "ADDRESS_US",
                                        "Spending": "SPENDING"
                                       })
selling_df.head()<jupyter_output><empty_output><jupyter_text>- Web<jupyter_code>selling_df[["WEB", "SPENDING"]].groupby("WEB").agg({"SPENDING": ["mean", "std"]})<jupyter_output><empty_output><jupyter_text>Web spending are higher on average compared to non-web spending; - Gender<jupyter_code>selling_df[["GENDER", "SPENDING"]].groupby("GENDER").agg({"SPENDING": ["mean", "std"]})<jupyter_output><empty_output><jupyter_text>Male spending are higher on average compared to female spending; - Address Residential<jupyter_code>selling_df[["ADDRESS_RES", "SPENDING"]].groupby("ADDRESS_RES").agg({"SPENDING": ["mean", "std"]})<jupyter_output><empty_output><jupyter_text>Non-address residential purchases spending are higher on average compared to residential addresses spending; - Address US<jupyter_code>selling_df[["ADDRESS_US", "SPENDING"]].groupby("ADDRESS_US").agg({"SPENDING": ["mean", "std"]})<jupyter_output><empty_output><jupyter_text>Address US purchases spending are almost the same (on average) compared to non-us addresses spending; **b.** Explore the relationship between spending and each of the two continuous predictors by creating two scatter plots (Spending vs. Freq, and Spending vs. last_update_days_ago). Does there seem to be a linear relationship?<jupyter_code>plt.scatter(selling_df.SPENDING, selling_df.FREQ, alpha=0.2)
plt.xlabel("SPENDING")
plt.ylabel("FREQ")
plt.show()

plt.scatter(selling_df.SPENDING, selling_df.LAST_UPDATE, alpha=0.2)
plt.xlabel("SPENDING")
plt.ylabel("LAST UPDATE")
plt.show()<jupyter_output><empty_output><jupyter_text>The relationship between spending and frequency of purchase appears appears to be linear, such that when the frequency increases, so does the amount spent. The relationship between spending and last update seems to be less linear. Though it appears to exist a linear relationship between the two variable.Computing the correlation matrix help us to be sure about such relationships:<jupyter_code>selling_df[["FREQ", "LAST_UPDATE", "SPENDING"]].corr()<jupyter_output><empty_output><jupyter_text>From the correlation matrix it is possible to see that both have a linear relationship, being the spending and frequency of transactions in the previous year a positive and strong correlation and the last update (number of days since last update to customer record) negatively, but no strongly, correlated with spending. **c.** To fit a predictive model for Spending:    i. Partition the 2000 records into training and validation sets.<jupyter_code># select predictors
predictors = ["FREQ", "LAST_UPDATE", "WEB", "GENDER",
              "ADDRESS_RES", "ADDRESS_US"]
outcome = "SPENDING"

# partition data, 40% test set
X = selling_df[predictors]
y = selling_df[outcome]
train_X, valid_X, train_y, valid_y = train_test_split(X, y, 
                                                      test_size=0.4,
                                                      random_state=1)<jupyter_output><empty_output><jupyter_text>ii. Run a multiple linear regression model for Spending vs. all six predictors. Give the estimated    predictive equation.<jupyter_code>selling_lmr = LinearRegression()
selling_lmr.fit(train_X, train_y)

# print coefficients
print(pd.DataFrame({'Predictor': train_X.columns,
                    'coefficient': selling_lmr.coef_}))
# print intercept
print('\nIntercept: ', selling_lmr.intercept_)
# print performance measures (training data)
regression_summary(y, selling_lmr.predict(X))<jupyter_output>Predictor  coefficient
0         FREQ    91.274450
1  LAST_UPDATE    -0.010374
2          WEB    18.628731
3       GENDER    -9.111366
4  ADDRESS_RES   -75.815354
5   ADDRESS_US    -4.620293

Intercept:  10.176297413561059

Regression statistics

               Mean Error (ME) : 2.8773
Root Mean Squared Error (RMSE) : 130.4020
     Mean Absolute Error (MAE) : 81.1267<jupyter_text>Therefore, we can write the equation for predicting the selling amount as:            $\hat{Y}_{SPENDING} = 10.17 + 91.27 x_{FREQ} - 0.01 x_{LAST\_UPDATE} + 18.63 x_{WEB} - 9.11 x_{GENDER} - 75.81 x_{ADD\_RES} - 4.62 x_{ADD\_US}$         iii. Based on this model, what type of purchaser is most likely to spend a large amount of money?Users that have a higher number of transactions in the past years, done on web, likely to be male and not purchasing for residential address. The other variables do not influence too much the model (`LAST_UPDATE` and `ADDRESS_US`).     iv. If we used backward elimination to reduce the number of predictors, which predictor would be dropped    first from the model?<jupyter_code>def train_model(variables):
    model = LinearRegression()
    model.fit(train_X[variables], train_y)
    return model

def score_model(model, variables):
    return AIC_score(train_y, model.predict(train_X[variables]), model)


all_variables = train_X.columns
best_model, best_variables = backward_elimination(all_variables, train_model, score_model,
                                                  verbose=True)
print(best_variables)
regression_summary(valid_y, best_model.predict(valid_X[best_variables]))<jupyter_output>Variables: FREQ, LAST_UPDATE, WEB, GENDER, ADDRESS_RES, ADDRESS_US
Start: score=15028.53
Step: score=15026.76, remove=ADDRESS_US
Step: score=15026.38, remove=GENDER
Step: score=15026.38, remove=None
['FREQ', 'LAST_UPDATE', 'WEB', 'ADDRESS_RES']

Regression statistics

               Mean Error (ME) : 6.9616
Root Mean Squared Error (RMSE) : 136.5274
     Mean Absolute Error (MAE) : 83.4472<jupyter_text>According to the backward elimitation the first variable would be `ADDRESS_US`, followed by `GENDER`;     v. Show how the prediction and the prediction error are computed for the first purchase in the validation    set.<jupyter_code>valid_X.iloc[0, :]
valid_y.iloc[0]<jupyter_output><empty_output><jupyter_text>The computation would be the following:            $\hat{Y}_{SPENDING} = 10.17 + 91.27 * 2 - 0.01 * 1346 + 18.63 * 0 - 9.11 * 1 - 75.81 * 1 - 4.62 * 1$    Resulting in a value of \\$89.71. The prediction error would be calculated as follows:            $err = Y - \hat{Y}_{SPENDING} = 0 - 89.71 = 89.71$    For the first 20 rows we have the following:<jupyter_code>selling_lmr_pred = selling_lmr.predict(valid_X)
result = pd.DataFrame({"Predicted": selling_lmr_pred,
                       "Actual": valid_y,
                       "Residual": valid_y - selling_lmr_pred})
print(result.head(20))<jupyter_output>Predicted  Actual     Residual
674    89.214915       0   -89.214915
1699  202.231362     184   -18.231362
1282   49.159303       0   -49.159303
1315  824.841659    1289   464.158341
1210    0.121196       0    -0.121196
1636   86.766675       0   -86.766675
613    58.018614       0   -58.018614
447   247.428569    1255  1007.571431
1131   67.036615       0   -67.036615
808    67.825031       0   -67.825031
1496   -7.098168       0     7.098168
1468  194.814024     411   216.185976
1682  -13.480101       0    13.480101
1149  -32.457046       0    32.457046
442    61.247979       0   -61.247979
1813    4.497885     173   168.502115
654   -46.046854       0    46.046854
1264  -32.315195       0    32.315195
858    80.219048       0   -80.219048
1482   51.783900       0   -51.783900<jupyter_text>vi. Evaluate the predictive accuracy of the model by examining its performance on the validation set.<jupyter_code># print performance measures (validation data)
regression_summary(valid_y, selling_lmr_pred)<jupyter_output>Regression statistics

               Mean Error (ME) : 7.1933
Root Mean Squared Error (RMSE) : 136.7397
     Mean Absolute Error (MAE) : 83.6010<jupyter_text>vii. Create a histogram of the model residuals. Do they appear to follow a normal distribution? How does    this affect the predictive performance of the model?<jupyter_code>all_residuals = valid_y - selling_lmr.predict(valid_X)

# Determine the percentage of datapoints with a residual in [-230, 230] = approx. 95%
print("Percentage of datapoints with a residual in [-230, 230]: ",
      len(all_residuals[(all_residuals > -230) & (all_residuals < 230)]) / len(all_residuals))

pd.DataFrame({"Residuals": all_residuals}).hist(bins=25)
plt.show()<jupyter_output>Percentage of datapoints with a residual in [-230, 230]:  0.955<jupyter_text>We can see from the above figure that the residuals follow a normal distribution. Then 95% of them should fall between -230 and 230. If they fall above 230 or below -230, they can be considered unusual. Besides, we can say that the predictions will be unbiased (equal to the true values on average) and will have the smallest mean squared error compared to any unbiased estimates if we make the following assumptions:1. The **noise ε (or equivalently, Y) follows a normal distribution**.2. The choice of predictors and their form is correct (linearity).3. The records are independent of each other.4. The variability in the outcome values for a given set of predictors is the same regardless of the values of the predictors (homoskedasticity).Since the noise ε follows a normal distribution, we know that these are the best predictions possible (unbiased). **3. Predicting Airfare on New Routes.**The following problem takes place in the United States in the late 1990s, when many major US cities were facing issues with airport congestion, partly as a result of the 1978 deregulation of airlines. Both fares and routes were freed from regulation, and low-fare carriers such as Southwest (SW) began competing on existing routes and starting nonstop service on routes that previously lacked it. Building completely new airports is generally not feasible, but sometimes decommissioned military bases or smaller municipal airports can be reconfigured as regional or larger commercial airports. There are numerous players and interests involved in the issue (airlines, city, state and federal authorities, civic groups, the military, airport operators), and an aviation consulting firm is seeking advisory contracts with these players. The firm needs predictive models to support its consulting service. One thing the firm might want to be able to predict is fares, in the event a new airport is brought into service. The firm starts with the file `Airfares.csv`, which contains real data that were collected between Q3-1996 and Q2-1997. The variables in these data are listed below, and are believed to be important in predicting FARE. Some airport-to-airport data are available, but most data are at the city-to-city level. One question that will be of interest in the analysis is the effect that the presence or absence of Southwest has on FARE.    S_CODE: Starting airport’s code    S_CITY: Starting city    E_CODE: Ending airport’s code    E_CITY: Ending city    COUPON: Average number of coupons (a one-coupon flight is a nonstop flight, a two-coupon flight is            a one-stop flight, etc.) for that route     NEW: Number of new carriers entering that route between Q3-1996 and Q2-1997    VACATION: Whether (Yes) or not (No) a vacation route    SW: Whether (Yes) or not (No) Southwest Airlines serves that route    HI: Herfindahl index: measure of market concentration    S_INCOME: Starting city’s average personal income    E_INCOME: Ending city’s average personal income    S_POP: Starting city’s population    E_POP: Ending city’s population    SLOT: Whether or not either endpoint airport is slot-controlled (this is a measure of airport congestion)    GATE: Whether or not either endpoint airport has gate constraints (this is another measure           of airport congestion)    DISTANCE: Distance between two endpoint airports in miles    PAX: Number of passengers on that route during period of data collection     FARE: Average fare on that route**a.** Explore the numerical predictors and outcome (FARE) by creating a correlation table and examining some scatterplots between FARE and those predictors. What seems to be the best single predictor of FARE?<jupyter_code>airfares_df = pd.read_csv("../datasets/Airfares.csv")
airfares_df.head()
# numerical predictors
columns = ["COUPON", "HI", "S_INCOME", "E_INCOME", "S_POP", "E_POP",
           "DISTANCE", "PAX", "NEW", "FARE"]
corr = airfares_df[columns].corr()

# include information about values
fig, ax = plt.subplots()
fig.set_size_inches(11, 7)
sns.heatmap(corr, annot=True, fmt=".1f", cmap="RdBu", center=0, ax=ax)
plt.show()<jupyter_output><empty_output><jupyter_text>Looking at the correlation matrix, the `COUPON` and `DISTANCE` exhibit interesting linear relationships with `FARE` - approx. 0.5 and 0.7, respectively. We can also see from the matrix that `DISTANCE` has the strongest correlation and therefore is the best single predictor of `FARE`. The `S_INCOME`, `E_INCOME` and `E_POP` also showed a positive correlation, though these correlations are not strong. Let's plot these variables:<jupyter_code>columns = ["COUPON", "S_INCOME", "E_INCOME", "E_POP", "DISTANCE", "FARE"]
plot_df = airfares_df[columns]

fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(8, 8),
                         tight_layout=True, sharey=True)
fig.delaxes(axes[2, 1])

axes[0][0].scatter(plot_df.COUPON, plot_df.FARE, s=3)
axes[0][1].scatter(plot_df.S_INCOME, plot_df.FARE, s=3)
axes[1][0].scatter(plot_df.E_INCOME, plot_df.FARE, s=3)
axes[1][1].scatter(plot_df.E_POP, plot_df.FARE, s=3)
axes[2][0].scatter(plot_df.DISTANCE, plot_df.FARE, s=3)

axes[0][0].set_ylabel("FARE")
axes[0][0].set_xlabel("COUPON")
axes[0][1].set_xlabel("S_INCOME")
axes[1][0].set_ylabel("FARE")
axes[1][0].set_xlabel("E_INCOME")
axes[1][1].set_xlabel("E_POP")
axes[2][0].set_ylabel("FARE")
axes[2][0].set_xlabel("DISTANCE")

plt.show()<jupyter_output><empty_output><jupyter_text>**b.** Explore the categorical predictors (excluding the first four) by computing the percentage of flights in each category. Create a pivot table with the average fare in each category. Which categorical predictor seems best for predicting FARE? - Vacation<jupyter_code>airfares_df[["VACATION", "FARE"]].groupby("VACATION").agg({"FARE": ["mean", "std"]})<jupyter_output><empty_output><jupyter_text>- SW<jupyter_code>airfares_df[["SW", "FARE"]].groupby("SW").agg({"FARE": ["mean", "std"]})<jupyter_output><empty_output><jupyter_text>- Slot<jupyter_code>airfares_df[["SLOT", "FARE"]].groupby("SLOT").agg({"FARE": ["mean", "std"]})<jupyter_output><empty_output><jupyter_text>- GATE<jupyter_code>airfares_df[["GATE", "FARE"]].groupby("GATE").agg({"FARE": ["mean", "std"]})<jupyter_output><empty_output><jupyter_text>All categorical values present different means for each of its categorical values (True/False), though two of them appear to have more predictive power over the others (`VACATION` and `SW`). The `SW` mean/std for each group seem to be particular interesting as a predictor of `FARE`.<jupyter_code># categorical predictors
columns = ["VACATION", "SW", "SLOT", "GATE", "FARE"]
corr = airfares_df[columns].copy()

# tranforming categorical predictors to bool
# another way to do this is using pd.get_dummies(airfares_df[columns])
corr["VACATION"] = corr.VACATION.map({"Yes": 1, "No": 0})
corr["SW"] = corr.SW.map({"Yes": 1, "No": 0})
corr["SLOT"] = corr.SLOT.map({"Controlled": 1, "Free": 0})
corr["GATE"] = corr.GATE.map({"Constrained": 1, "Free": 0})
corr = corr.corr()
# include information about values
fig, ax = plt.subplots()
fig.set_size_inches(11, 7)
sns.heatmap(corr, annot=True, fmt=".1f", cmap="RdBu", center=0, ax=ax)
plt.show()<jupyter_output><empty_output><jupyter_text>**c.** Find a model for predicting the average fare on a new route:    i. Convert categorical variables into dummy variables. Then, partition the data into training and       validation  sets. The model will be fit to the training data and evaluated on the validation set.<jupyter_code># categorical predictors
columns = ["VACATION", "SW", "SLOT", "GATE"]
data_df = airfares_df[columns].copy()

# tranforming categorical predictors to bool
# another way to do this is using pd.get_dummies(airfares_df[columns])
data_df["VACATION"] = data_df.VACATION.map({"Yes": 1, "No": 0})
data_df["SW"] = data_df.SW.map({"Yes": 1, "No": 0})
data_df["SLOT"] = data_df.SLOT.map({"Controlled": 1, "Free": 0})
data_df["GATE"] = data_df.GATE.map({"Constrained": 1, "Free": 0})

# numerical predictors
columns = ["COUPON", "HI", "S_INCOME", "E_INCOME", "S_POP", "E_POP",
           "DISTANCE", "PAX", "NEW", "FARE"]

data_df = data_df.join(airfares_df[columns]).copy()
data_df
predictors = data_df.columns[:-1]
outcome = data_df.columns[-1]

# partition data, 40% test set
X = data_df[predictors]
y = data_df[outcome]
train_X, valid_X, train_y, valid_y = train_test_split(X, y, 
                                                      test_size=0.4,
                                                      random_state=1)<jupyter_output><empty_output><jupyter_text>ii. Use stepwise regression to reduce the number of predictors. You can ignore the first four predictors        (S_CODE, S_CITY, E_CODE, E_CITY). Report the estimated model selected.<jupyter_code># The initial model is the constant model - this requires special handling
# in train_model and score_model
def train_model(variables):
    if len(variables) == 0:
        return None
    model = LinearRegression()
    model.fit(train_X[variables], train_y)
    return model

def score_model(model, variables):
    if len(variables) == 0:
        return AIC_score(train_y, [train_y.mean()] * len(train_y), model, df=1)
    return AIC_score(train_y, model.predict(train_X[variables]), model)

best_model, best_variables = stepwise_selection(train_X.columns, train_model,
                                                score_model, verbose=True)
print(best_variables)
regression_summary(valid_y, best_model.predict(valid_X[best_variables]))<jupyter_output>Variables: VACATION, SW, SLOT, GATE, COUPON, HI, S_INCOME, E_INCOME, S_POP, E_POP, DISTANCE, PAX, NEW
Start: score=4379.42, constant
Step: score=4149.88, add DISTANCE
Step: score=4025.89, add SW
Step: score=3913.59, add VACATION
Step: score=3890.27, add HI
Step: score=3873.33, add GATE
Step: score=3852.81, add SLOT
Step: score=3850.30, add PAX
Step: score=3844.58, add E_POP
Step: score=3826.43, add S_POP
Step: score=3821.88, add E_INCOME
Step: score=3821.88, unchanged None
['DISTANCE', 'SW', 'VACATION', 'HI', 'GATE', 'SLOT', 'PAX', 'E_POP', 'S_POP', 'E_INCOME']

Regression statistics

                      Mean Error (ME) : -1.8591
       Root Mean Squared Error (RMSE) : 36.1129
            Mean Absolute Error (MAE) : 28.5252
          Mean Percentage Error (MPE) : -6.7084
Mean Absolute Percentage Error (MAPE) : 21.7198<jupyter_text>iii. Repeat (ii) using exhaustive search instead of stepwise regression. Compare the resulting best         model to the one you obtained in (ii) in terms of the predictors that are in the model.<jupyter_code>def train_model(variables):
    model = LinearRegression()
    model.fit(train_X[list(variables)], train_y)
    return model

def score_model(model, variables):
    pred_y = model.predict(train_X[list(variables)])
    # we negate as score is optimized to be as low as possible
    return -adjusted_r2_score(train_y, pred_y, model)

all_variables = train_X.columns
results = exhaustive_search(all_variables, train_model, score_model)

data = []
for result in results:
    model = result["model"]
    variables = list(result["variables"])
    AIC = AIC_score(train_y, model.predict(train_X[variables]), model)

    d = {"n": result["n"], "r2adj": -result["score"], "AIC": AIC}
    d.update({var: var in result["variables"] for var in sorted(all_variables)})
    data.append(d)

pd.DataFrame(data, columns=("n", "r2adj", "AIC") + tuple(sorted(all_variables)))<jupyter_output><empty_output><jupyter_text>The chosen ten-predictor model chosen by the exhaustive model are: `DISTANCE`, `E_INCOME`, `E_POP`, `GATE`, `HI`, `PAX`, `SLOT`, `SW`, `S_POP`, `VACATION`. The result of stepwise regression is identical to the best ten-predictor model chosen by the exhaustive model.     iv. Compare the predictive accuracy of both models (ii) and (iii) using measures such as RMSE and        average error and lift charts.Once the models were the same, the following measures and graph apply to both of them. **Important note:**Cumulative gains and lift charts for a continuous response are relevant only when we are searching for a set of records that gives the highest cumulative predicted values. Such charts are not relevant if we are interested in predicting the outcome value for each new record.<jupyter_code>pred_y = pd.Series(best_model.predict(valid_X[best_variables]))

regression_summary(valid_y, pred_y)
print("\nr2adj: ", adjusted_r2_score(valid_y, pred_y, best_model))

fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(8, 4))
ax = gains_chart(pred_y.sort_values(ascending=False), ax=axes[0])
ax.set_ylabel("Cumulative Fares")
ax.set_title("Cumulative Gains Chart")

ax = lift_chart(pred_y.sort_values(ascending=False), ax=axes[1])
ax.set_ylabel("Lift")
plt.tight_layout()
plt.show()<jupyter_output>Regression statistics

                      Mean Error (ME) : -1.8591
       Root Mean Squared Error (RMSE) : 36.1129
            Mean Absolute Error (MAE) : 28.5252
          Mean Percentage Error (MPE) : -6.7084
Mean Absolute Percentage Error (MAPE) : 21.7198

r2adj:  0.7790766452927753<jupyter_text>The above figure shows a cumulative gains chart and decile lift chart based on fitting a linear regression model to the `Airfares.csv`. The charts are based on the validation data of 256 fares and routes. It can be seen that themodel's predictive performance in terms of gains is better than the baseline model (average), since its cumulative gains curve is higher than that of the baseline model.The overall measures of predictive accuracy show the following:                            Mean Error (ME) : $-1.8591         Root Mean Squared Error (RMSE) : $36.1129                  Mean Absolute Error (MAE) : $28.5252            Mean Percentage Error (MPE) : $-6.7084      Mean Absolute Percentage Error (MAPE) : $21.7198Note that, in terms of the performance, this model resulted pretty good.     v. Using model (iii), predict the average fare on a route with the following characteristics:       COUPON = 1.202, NEW = 3, VACATION = No, SW = No, HI = 4442.141, S_INCOME = $28,760, E_INCOME = $27,664,       S_POP = 4,557,004, E_POP = 3,195,503, SLOT = Free,       GATE = Free, PAX = 12,782, DISTANCE = 1976 miles.<jupyter_code># the example of a data frame home to be scored using our regression model
new_data = pd.DataFrame({"VACATION": [0], "SW": [0], "SLOT": [0], "GATE": [0],
                         "COUPON": [1.202], "HI": [4442.141], "S_INCOME": [28760], "E_INCOME": [27664],
                         "S_POP": [4557004], "E_POP": [3195503], "DISTANCE": [1976], "PAX": [12782], "NEW": [3]})

# using the best variables only
print("\nPredictions: ", best_model.predict(new_data[best_variables]))
new_data[[col for col in new_data.columns if col in best_variables]]<jupyter_output>Predictions:  [251.72153114]<jupyter_text>Let's compare the new data with others with similar `FARE`:<jupyter_code>mask = (data_df.FARE >= 245) & (data_df.FARE <= 255)
data_df[mask]<jupyter_output><empty_output><jupyter_text>We can say they are reasonably similar in terms of predictors, specially the important ones (e.g. `SW`, `VACATION`)     vi. Predict the reduction in average fare on the route in (v) if Southwest decides to cover this route        [using model (iii)].<jupyter_code># the example of a data frame home to be scored using our regression model
new_data = pd.DataFrame({"VACATION": [0], "SW": [1], "SLOT": [0], "GATE": [0],
                         "COUPON": [1.202], "HI": [4442.141], "S_INCOME": [28760], "E_INCOME": [27664],
                         "S_POP": [4557004], "E_POP": [3195503], "DISTANCE": [1976], "PAX": [12782], "NEW": [3]})

# using the best variables only
print("\nPredictions: ", best_model.predict(new_data[best_variables]))
new_data[[col for col in new_data.columns if col in best_variables]]<jupyter_output>Predictions:  [208.69025908]<jupyter_text>Almost \\$50 dollars!     vii. In reality, which of the factors will not be available for predicting the average fare         from a new airport (i.e., before flights start operating on those routes)?         Which ones can be estimated? How?The variables that correlates with the number of passengers/carriers/market would be difficult to estimate before flights start operating on a specific route (`PAX`, `NEW` and `HI`). `COUPON`, `SW`, Incomes, and Populations can all be predicted with a high level of confidence since those factors are independent the proposed flight. `VACATION`, `SLOT`, and `GATE` can be predicted based on the proposed route and knowledge of the airport's policies but are highly susceptible to changes.    viii. Select a model that includes only factors that are available before flights begin to operate          on the new route. Use an exhaustive search to find such a model.This means we will need to remove `PAX`, `NEW` and `HI` variables. Looking at the previous computation by the exhaustive search table, the best model would be a three-predictive model with these variables:- DISTANCE- SW- VACATION     ix. Use the model in (viii) to predict the average fare on a route with characteristics:        COUPON = 1.202, NEW = 3, VACATION = No, SW = No, HI = 4442.141, S_INCOME = 28,760,        E_INCOME = $27,664, S_ POP = 4,557,004, E_POP = 3,195,503, SLOT = Free,        GATE = Free, PAX = 12,782, DISTANCE = 1976 miles.<jupyter_code>predictors = ["DISTANCE", "SW", "VACATION"]
outcome = "FARE"

airfare_lmr = LinearRegression()
airfare_lmr.fit(train_X[predictors], train_y)

# the example of a data frame home to be scored using our regression model
new_data = pd.DataFrame({"VACATION": [0], "SW": [1], "SLOT": [0], "GATE": [0],
                         "COUPON": [1.202], "HI": [4442.141], "S_INCOME": [28760], "E_INCOME": [27664],
                         "S_POP": [4557004], "E_POP": [3195503], "DISTANCE": [1976], "PAX": [12782], "NEW": [3]})

# using the best variables only
print("\nPredictions: ", airfare_lmr.predict(new_data[predictors]))
new_data[[col for col in new_data.columns if col in predictors]]<jupyter_output>Predictions:  [204.90881318]<jupyter_text>Not bad, \\$50 less compared to the model using ten variables.     x. Compare the predictive accuracy of this model with model (iii). Is this model good       enough, or is it worthwhile reevaluating the model once flights begin on the new route?<jupyter_code>pred_y = pd.Series(airfare_lmr.predict(valid_X[predictors]))

regression_summary(valid_y, pred_y)
print("\nr2adj: ", adjusted_r2_score(valid_y, pred_y, airfare_lmr))

fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(8, 4))
ax = gains_chart(pred_y.sort_values(ascending=False), ax=axes[0])
ax.set_ylabel("Cumulative Fares")
ax.set_title("Cumulative Gains Chart")

ax = lift_chart(pred_y.sort_values(ascending=False), ax=axes[1])
ax.set_ylabel("Lift")
plt.tight_layout()
plt.show()<jupyter_output>Regression statistics

                      Mean Error (ME) : -1.2464
       Root Mean Squared Error (RMSE) : 42.5705
            Mean Absolute Error (MAE) : 33.7933
          Mean Percentage Error (MPE) : -8.0557
Mean Absolute Percentage Error (MAPE) : 25.1156

r2adj:  0.7015303109350637<jupyter_text>The above table shows the validation performance of the 3-predictor model as well as the gain and lift charts, which turns out to be only slightly worse than the 10-predictor model in terms of validation metrics. Since this model performs only slightly worse, it would be preferable from a **parsimony** standpoint. **d.** In competitive industries, a new entrant with a novel business plan can have a disruptive effect on existing firms. If a new entrant's business model is sustainable, other players are forced to respond by changing their business practices. If the goal of the analysis was to evaluate the effect of Southwest Airlines' presence on the airline industry rather than predicting fares on new routes, how would the analysis be different? Describe technical and conceptual aspects.In particular, the difference in the objective behind fitting a regression model would be:    Explaining or quantifying the average effect of inputs on an outcome (explanatory or descriptive task,    respectively) rather than predicting the outcome value for new records, given their input values    (predictive task)    In that scenario, the data would be treated as a random sample from a larger population of interest. The regression model estimated from this sample would be an attempt to capture the *average* relationship in the larger population when comparing the group of data that has `SW` on new routes and those that has not. This model is then used in decision-making to generate statements such as "the `SW` presence ($X_{SW}$) is associated with an average increase of \\$50 dollars in the average fair ($Y$), all other factors ($X_{VACATION}$, $X_{GATE}$, ..., $X_{POPULATION}$) being equal." If $X_{SW}$ is known to *cause* $Y$, then such a statement indicates actionable policy changes - this is called explanatory modeling. When the causal structure is unknown, then this model quantifies the degree of *association* between the inputs and outcome variable, and the approach is called descriptive modeling.Besides that, in explanatory and descriptive modeling, where the focus is on modeling the average record, we try to fit the best model to the data in an attempt to learn about the underlying relationship in the population. In contrast, in predictive modeling (data mining), the goal is to find a regression model that best predicts new individual records. A regression model that fits the existing data too well is not likely to perform well with new data. **4. Predicting Prices of Used Cars.** The file `ToyotaCorolla.csv` contains data on used cars (Toyota Corolla) on sale during late summer of 2004 in the Netherlands. It has 1,436 records containing details on 38 attributes, including Price, Age, Kilometers, HP, and other specifications. The goal is to predict the price of a used Toyota Corolla based on its specifications.Split the data into training (50%), validation (30%), and test (20%) datasets. Run a multiple linear regression with the outcome variable `Price` and predictor variables `Age_08_04`, `KM`, `Fuel_Type`, `HP`, `Automatic`, `Doors`, `Quarterly_Tax`, `Mfr_Guarantee`, `Guarantee_Period`, `Airco`, `Automatic_airco`, `CD_Player`, `Powered_Windows`, `Sport_Model`, and `Tow_Bar`.**a.** What appear to be the three or four most important car specifications for predicting the car’s price?Some considerations before start this exercise. The importance of a variable in a regression models depends of some attributes of the model. In summary, they are:- Larger coefficients don't necessarily identify more important predictor variables. Why? Regression coefficients use different scales, so we can't compare them directly. However, if we standardize/normalize the regression coefficients, so they're based on the same scale, we can compare them.- Low p-values don't necessarily identify predictor variables that are practically important. Why? A very low p-value can reflect properties other than importance, such as a very precise estimate and a large sample size- Combining these two statistics we can determine variable importance. Therefore, let's use the `statsmodels` here. Let's also normalize the data.From [here](https://blog.minitab.com/en/adventures-in-statistics-2/how-to-identify-the-most-important-predictor-variables-in-regression-models).<jupyter_code>car_df = pd.read_csv('../datasets/ToyotaCorolla.csv')

# Fuel_Type - a categorical variable with 3 different values
# We will use only N-1 dummy variables when we have categorical variables
car_df.Fuel_Type = car_df.Fuel_Type.astype("category")
car_df = pd.get_dummies(car_df,
                        columns=["Fuel_Type"], prefix_sep="_",
                        drop_first=False)

predictors = ['Age_08_04', 'KM', 'Fuel_Type_CNG', 'Fuel_Type_Diesel', 'HP', 'Automatic',
              'Doors', 'Quarterly_Tax', 'Mfr_Guarantee', 'Guarantee_Period', 'Airco',
              'Automatic_airco', 'CD_Player', 'Powered_Windows', 'Sport_Model', 'Tow_Bar']
outcome = 'Price'

# normalizing data
scaler = MinMaxScaler()
car_df_norm = pd.DataFrame(scaler.fit_transform(car_df[predictors + [outcome]]),
                           columns=predictors + [outcome])


# partition data
# random_state is set to a defined value to get the same partitions when re-running the code
# training: 50 %
# validation: 30% => 0.5 * 0.6
# test: 20% => 0.5 * 0.4
train_data, temp = train_test_split(car_df_norm[predictors + [outcome]], test_size=0.50, random_state=1)
valid_data, test_data = train_test_split(temp, test_size=0.40, random_state=1)

print("Training:   ", train_data.shape)
print("Validation: ", valid_data.shape)
print("Test:       ", test_data.shape)

train_data.head()
# run a linear regression of Price on the remaining 11 predictors in the
# training set
train_df = train_data[predictors].join(train_data[outcome])
formula = 'Price ~ ' + ' + '.join(predictors)
car_lm = sm.ols(formula=formula, data=train_df).fit()
car_lm.summary()<jupyter_output><empty_output><jupyter_text>Most important variables are:- `Age_08_04`- `KM`- `HP`- `Quarterly_Tax`- `Automatic_airco` **b.** Using metrics you consider useful, assess the performance of the model in predicting prices. Let's keep redo the model maintaining the `Price` scale, so we can assess the errors in their true dimensions:<jupyter_code># partition data
# random_state is set to a defined value to get the same partitions when re-running the code
# training: 50 %
# validation: 30% => 0.5 * 0.6
# test: 20% => 0.5 * 0.4
train_data, temp = train_test_split(car_df[predictors + [outcome]], test_size=0.50, random_state=1)
valid_data, test_data = train_test_split(temp, test_size=0.40, random_state=1)

car_lmr = LinearRegression()
car_lmr.fit(train_data[predictors], train_data[outcome])

# print coefficients
print(pd.DataFrame({'Predictor': predictors,
                    'coefficient': car_lmr.coef_}))
# print performance measures (training data)
regression_summary(train_data[outcome],
                   car_lmr.predict(train_data[predictors]))
# print performance measures (valid data)
pred_y = pd.Series(car_lmr.predict(valid_data[predictors]))

regression_summary(valid_data[outcome], pred_y)
print("\nr2adj - validation: ", adjusted_r2_score(valid_data[outcome],
                                                  pred_y, car_lmr))

fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(8, 4))
ax = gains_chart(pred_y.sort_values(ascending=False), ax=axes[0])
ax.set_ylabel("Cumulative Prices")
ax.set_title("Cumulative Gains Chart")

ax = lift_chart(pred_y.sort_values(ascending=False), ax=axes[1])
ax.set_ylabel("Lift")
plt.tight_layout()

plt.show()
# print performance measures (test data)
pred_y = pd.Series(car_lmr.predict(test_data[predictors]))

regression_summary(test_data[outcome], pred_y)
print("\nr2adj - test: ", adjusted_r2_score(test_data[outcome],
                                                  pred_y, car_lmr))

fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(8, 4))
ax = gains_chart(pred_y.sort_values(ascending=False), ax=axes[0])
ax.set_ylabel("Cumulative Prices")
ax.set_title("Cumulative Gains Chart")

ax = lift_chart(pred_y.sort_values(ascending=False), ax=axes[1])
ax.set_ylabel("Lift")
plt.tight_layout()<jupyter_output>Regression statistics

                      Mean Error (ME) : 91.2764
       Root Mean Squared Error (RMSE) : 1271.8760
            Mean Absolute Error (MAE) : 958.7654
          Mean Percentage Error (MPE) : 0.3914
Mean Absolute Percentage Error (MAPE) : 9.6015

r2adj - test:  0.8672175994750433
StarcoderdataJupyterStructuredCleanDedup