<jupyter_start><jupyter_text>There are 50 observations and 5 columns. 4 columns - R&D Spend, Administration and Marketing Spend, and Profile are numeric and one is categorical - State. There is 2 null values in columns R&D Spend feature and and 3 in Marketing Spend. Replace the null values with median for respective state.<jupyter_code>df_null_idx = df[df.isnull().sum(axis = 1) > 0].index
df.iloc[df_null_idx]
median_values = df.groupby("State")[["R&D Spend", "Marketing Spend"]].median()
median_values
df["R&D Spend"] = df.apply(lambda row: median_values.loc[row["State"], "R&D Spend"]  if np.isnan(row["R&D Spend"]) else row["R&D Spend"], axis = 1 )
df["Marketing Spend"] = df.apply(lambda row: median_values.loc[row["State"], "Marketing Spend"]  if np.isnan(row["Marketing Spend"]) else row["Marketing Spend"], axis = 1 )
df.iloc[df_null_idx]
# Check if there are any more null values.
df.isnull().sum()<jupyter_output><empty_output><jupyter_text>Let's see the distribution of the Profit using a histogram plot and see if there is any outliers in the data using bosplot.<jupyter_code>plt.figure(figsize = (8, 6))
plt.subplot(2, 1, 1)
df.Profit.plot.hist(bins = 10, normed = True)
df.Profit.plot.kde(title = "Historgram of Profit")

plt.subplot(2, 1, 2)
df.Profit.plot.box(vert = False, title = "Boxplot of Profit")

plt.tight_layout()<jupyter_output><empty_output><jupyter_text>Profit has one outlier. We can try to take log scale to remove the outlier value before doing any prediction. But for now, let ignore the outlier.Let's plot association between each pair of columns.<jupyter_code>sns.pairplot(df)<jupyter_output><empty_output><jupyter_text>Displays only the numeric column. Let's how the avg Profit plays for each State.<jupyter_code>df.groupby("State").Profit.mean().sort_values().plot.bar(title = "Avg Profit by State")
plt.xlabel("State")
plt.ylabel("Profit")<jupyter_output><empty_output><jupyter_text>Avg Profit is highest in state of Florida and least in California.Let's create the y vector containing the outcome column.<jupyter_code>y = df.Profit.values
y<jupyter_output><empty_output><jupyter_text>Create dummy variables for categorical feature.<jupyter_code>df_features = df.iloc[:, 0:4]
df_dummied = pd.get_dummies(df_features, columns=["State"], drop_first=True)
df_dummied.sample(10)<jupyter_output><empty_output><jupyter_text>State column has been replaced by two additional column - one for Florida and one NY. First value in the categorical values CA has been dropped to avoid collinearity issue.Now, let's create X feature matrix and y outcome vector.<jupyter_code>X = df_dummied.values
X[0, :]<jupyter_output><empty_output><jupyter_text>Let's normalize the feature values to bring them to a similar scale.<jupyter_code>scaler = StandardScaler() 
X_std = scaler.fit_transform(X)
pd.DataFrame(X_std).head()<jupyter_output><empty_output><jupyter_text>Split the X and y into training and test sets.<jupyter_code>X_train, X_test, y_train, y_test = train_test_split(X_std, y, 
                                                    test_size = 0.3, random_state = 100)
print("Training set: ", X_train.shape, y_train.shape)
print("Test set: ", X_test.shape, y_test.shape)<jupyter_output>Test set:  (15, 5) (15,)<jupyter_text>Ratio of the size of the training data<jupyter_code>X_train.shape[0] / df.shape[0]<jupyter_output><empty_output><jupyter_text>Fit linear regression model<jupyter_code>lr = LinearRegression()
lr.fit(X_train, y_train)
lr.intercept_, lr.coef_<jupyter_output><empty_output><jupyter_text>By looking at the cofficients, we can conclude that R&D Spend has the higest influence on the outcome variable.Predict the outcome based on the model<jupyter_code>y_test_pred = lr.predict(X_test)
output = pd.DataFrame({"actual": y_test, "prediction": y_test_pred})
output["error"] = output.actual - output.prediction
output<jupyter_output><empty_output><jupyter_text>A simpliest prediction model could have been the average. Let's how the model did overall against one feature.<jupyter_code>X_test_inv = scaler.inverse_transform(X_test)
plt.scatter(X_test_inv[:, 0], y_test, alpha = 0.3, c = "blue", label = "Actual")
plt.scatter(X_test_inv[:, 0], y_test_pred, c = "red", label = "Predicted")

plt.xlabel("R&D Spend")
plt.ylabel("Profit")
plt.title("Profit Actual vs Estimate")
plt.legend()
np.mean((y_test_pred - y_test) ** 2)
y_train_pred = lr.predict(X_train)<jupyter_output><empty_output><jupyter_text>Compare the root mean squared error (RMSE) of test dataset against the training.<jupyter_code>print("Test rmse: ", sqrt(mean_squared_error(y_test, y_test_pred)), 
      "\nTraining rmse:", sqrt(mean_squared_error(y_train, y_train_pred)))<jupyter_output>Test rmse:  24074.424825793423 
Training rmse: 12413.672826747377<jupyter_text>r2 score can have a max value 1, negative values of R2 means suboptimal model<jupyter_code>r2_score(y_test, y_test_pred), r2_score(y_train, y_train_pred)<jupyter_output><empty_output><jupyter_text>On the training the both RMSE and R2 scores perform natually better than those on the test dataset. Let's calculate R2 score manually.<jupyter_code>SSR = np.sum((y_train - y_train_pred) ** 2)  # Sum of squared residuals
SST = np.sum((y_train - np.mean(y_train_pred)) ** 2) # Sum of squared totals
R2 = 1 - SSR/SST
R2<jupyter_output><empty_output><jupyter_text>R2 can be viewed as (1 - mse/variance(y))  Significance Scores for feature selection<jupyter_code>from sklearn.feature_selection import f_regression
_, p_vals = f_regression(X_train, y_train)
p_vals
pd.DataFrame({"feature": df_dummied.columns, "p_value": p_vals})<jupyter_output><empty_output><jupyter_text>p-value indicates the significant scores for each feature. p-value < 0.05 indicates, the corresponding feature is statistically significant. We can rebuild the model excluding the non-significant features one by one until all remaining features are significant.    Power Plant DatasetLet's look at another dataset<jupyter_code>df = pd.read_csv("/data/Combined_Cycle_Power_Plant.csv")
df.head()
X = df.iloc[:, 0:4].values
y = df.PE.values
sns.pairplot(df)
scaler = StandardScaler()
X_std = scaler.fit_transform(X)
X_train, X_test, y_train, y_test = train_test_split(X_std, y, test_size = 0.3, random_state = 1)
def rmse(y_true, y_pred):
    return sqrt(mean_squared_error(y_true, y_pred))
lr = LinearRegression(normalize=False)
lr.fit(X_train, y_train)
y_train_pred = lr.predict(X_train)
y_test_pred = lr.predict(X_test)
rmse(y_test, y_test_pred)
from scipy import stats
residuals = y_test - y_test_pred

plt.figure(figsize=(15, 6))
plt.subplot(1, 2, 1)
plt.scatter(y_test, residuals)
plt.xlabel("y_test")
plt.ylabel("Residuals")
plt.hlines([0], xmin = 420, xmax = 500, linestyles = "dashed")

plt.subplot(1, 2, 2)
stats.probplot(residuals, plot=plt)<jupyter_output><empty_output><jupyter_text>Residual plots show there are outliers in the lower end of the y_test values. qqPlot shows that residuals do not exhibit normaality, indicating non linearity in the model.<jupyter_code>poly = PolynomialFeatures(degree=2)

X = df.iloc[:, 0:4].values
X_poly = poly.fit_transform(X)
X_poly_train, X_poly_test, y_train, y_test = train_test_split(X_poly, y, test_size = 0.3, random_state = 100)
X_poly_train_std = scaler.fit_transform(X_poly_train)
X_poly_test_std = scaler.transform(X_poly_test)

pd.DataFrame(X_poly_train_std).head()
lr.fit(X_poly_train_std, y_train)
print("Train rmse: ", rmse(y_train, lr.predict(X_poly_train_std)))
print("Test rmse: ", rmse(y_test, lr.predict(X_poly_test_std)))
print(lr.intercept_, lr.coef_)<jupyter_output>454.309651218 [  0.     -39.9701 -42.9304  92.6406  61.3886   4.9224   7.0201  20.0357
  -2.6156  -1.2008  36.9954   0.0829 -91.48   -57.9414  -3.7064]<jupyter_text>Polynomial regression generally sufferes from overfitting. Let's regularize the model using Lasso.<jupyter_code>lasso = Lasso(alpha=0.03, max_iter=10000, normalize=False, random_state=100)
lasso.fit(X_poly_train_std, y_train)
print("Train rmse: ", rmse(y_train, lasso.predict(X_poly_train_std)))
print("Test rmse: ", rmse(y_test, lasso.predict(X_poly_test_std)))
print(lasso.intercept_, lasso.coef_)<jupyter_output>Train rmse:  4.330404075518093
Test rmse:  4.35914090833783
454.309651218 [  0.     -17.3624  -1.4063   0.4679   0.       4.471    0.      -0.
  -1.0568  -0.      -0.435   -1.6751   0.       0.      -0.    ]<jupyter_text>Let's find cross validation score that accuracy score is more reliable in a sense that it incorporates every piece of is incorporated in both training and testing.<jupyter_code>X_poly_std = scaler.fit_transform(X_poly)
lasso = Lasso(alpha=0.03, max_iter=10000, random_state=100)
scores = cross_val_score(lasso, X_poly_std, y, cv = 10, scoring="neg_mean_squared_error")
scores = np.sqrt(-scores)
print("RMSE scores", scores)
print("Mean rmse: ", np.mean(scores))<jupyter_output>RMSE scores [ 4.307   4.371   4.3509  4.3279  4.3506  4.3315  4.4165  4.2604  4.3743
  4.3033]
Mean rmse:  4.33935182329<jupyter_text>Encapsulate the steps in a pipeline<jupyter_code>from sklearn.pipeline import Pipeline
pipeline = Pipeline(steps = [
    ("poly", PolynomialFeatures(degree=2, include_bias=False)),
    ("scaler", StandardScaler()),
    ("lasso", Lasso(alpha=0.03, max_iter=10000, normalize=False, random_state=1))
])

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.3, random_state = 1)
pipeline.fit(X_train, y_train)
rmse(y_test, pipeline.predict(X_test))<jupyter_output><empty_output><jupyter_text>LassoCV helps find the best alpha. We could also use model tuning techqniues to find best alpha as well.<jupyter_code># Find best alpha
lassocv = LassoCV(cv = 10, max_iter=10000, tol=1e-5)
lassocv.fit(X_poly_std, y)
print("Lassocv alpha: ", lassocv.alpha_)

# Apply the best alpha to find cross validation score
lasso = Lasso(alpha = lassocv.alpha_, max_iter=10000, random_state=100)
scores = cross_val_score(lasso, X_poly_std, y, cv = 10, scoring="neg_mean_squared_error")
print("Mean rmse: ", np.mean(np.sqrt(-scores)))<jupyter_output>Lassocv alpha:  0.0161806999848
Mean rmse:  4.32254217026<jupyter_text>Look at the cofficients values. Many of the features are not zero making the model parsimonious hence more robust - that is less prone to overfitting. Let's plot how coefficient reached 0 values by varying the alpha valuess.<jupyter_code>coefs = []
alphas = 10 ** np.linspace(-5, 5, 20)
for alpha in alphas:
    lasso = Lasso(alpha=alpha, max_iter=10000, tol=1e-5,random_state=100)
    lasso.fit(X_poly_std, y)
    coefs.append(lasso.coef_)

plt.plot(alphas, coefs)
plt.xscale("log")
plt.xlabel("Alpha (penalty term on the coefficients)")
plt.ylabel("Coefficients of the features")<jupyter_output>/Users/abulbasar/anaconda/lib/python3.6/site-packages/sklearn/linear_model/coordinate_descent.py:491: ConvergenceWarning: Objective did not converge. You might want to increase the number of iterations. Fitting data with very small alpha may cause precision problems.
  ConvergenceWarning)<jupyter_text>From this graph, which alpha values should we select. That question can be answered by looking which alpha values gives the best performance (rmse for example). lassocv function does that for us, or we can use model tuning techniques using grid search - that will be explained later.  Xgboost<jupyter_code>poly = PolynomialFeatures(degree=2)

X = df.iloc[:, 0:4].values
X_poly = poly.fit_transform(X)
X_poly_train, X_poly_test, y_train, y_test = train_test_split(X_poly, y, test_size = 0.3, random_state = 100)
X_poly_train_std = scaler.fit_transform(X_poly_train)
X_poly_test_std = scaler.transform(X_poly_test)
gbm = xgb.XGBRegressor(max_depth=10, learning_rate=0.1, n_estimators=100, 
                     objective='reg:linear', booster='gbtree', 
                     reg_alpha=0.01, reg_lambda=1, random_state=0)
gbm.fit(X_poly_train_std, y_train)
print("rmse:", rmse(y_test, gbm.predict(X_poly_test_std)))
param = {'silent':1, 
         'objective':'reg:linear', 
         'booster':'gbtree',
         'alpha': 0.01, 
         'lambda': 1
        }

dtrain = xgb.DMatrix(X_poly_train_std, label=y_train)
dtest = xgb.DMatrix(X_poly_test_std, label=y_test)
watchlist  = [(dtrain,'eval'), (dtest, 'train')]
num_round = 100
bst = xgb.train(param, dtrain, num_round, watchlist, verbose_eval=False)
print("rmse:", rmse(y_test, bst.predict(dtest)))

plt.figure(figsize=(8, 10))
xgb.plot_importance(bst)<jupyter_output>rmse: 2.019509243562749
{'fraction_non_alphanumeric': 0.11649782923299566, 'fraction_numerical': 0.04076218041485769, 'mean_word_length': 5.749321758003256, 'pattern_counts': {'":': 4, '<': 94, '<?xml version=': 0, '>': 94, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '25', 'dir': 'jupyter-structured-clean-dedup', 'id': '343109', 'license': 'Apache-2.0', 'n_tokens_mistral': 4921, 'n_tokens_neox': 4497, 'n_words': 1104, 'path': 'Scikit - 03 Linear Regression.ipynb', 'repo_name': 'ashishpatel26/Machine-Learning-and-Deep-Learning-with-sklearn-and-sparkml', 'size': '12438'}