<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Linear regression with statsmodels (OLS)
Step1:  We load the boston house-prices dataset and X are our features and y is the target variable medv (Median value of owner-occupied homes in $1000s).
Step2:  Let's split the data in a test and training set.
Step3:  Fitting models - the standard way
Full model without an intercept
Step4:  Full model with an intercept
Step5:  Fitting models using R-style formulas
We can also fit a model with the R syntax y ~ x_1 + x_2 and build some complexer models.
Step6:  Full model with an intercept
Step7:  Model with a polynomial and the target variable log transformed
Step8:  Let's plot the QQ-Plot for the residuals
<ASSISTANT_TASK:>
Python Code:
import statsmodels.api as sm
import statsmodels.formula.api as smf
import statsmodels.tools.eval_measures as eval_measures
import seaborn as sns
import scipy.stats as stats
import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
from sklearn import datasets, model_selection, metrics
Explanation: Linear regression with statsmodels (OLS)
End of explanation
boston = datasets.load_boston()
print(boston.DESCR)
Explanation: We load the boston house-prices dataset and X are our features and y is the target variable medv (Median value of owner-occupied homes in $1000s).
End of explanation
X = pd.DataFrame(boston.data, columns=boston.feature_names)
y = boston.target
X_train, X_test, y_train, y_test = model_selection.train_test_split(X, y, train_size=0.7)
df_train = pd.DataFrame(y_train, columns=['target'])
df_train['type'] = 'train'
df_test = pd.DataFrame(y_test, columns=['target'])
df_test['type'] = 'test'
df_set = df_train.append(df_test)
_ = sns.displot(df_set, x="target" ,hue="type", kind="kde", log_scale=False)
Explanation: Let's split the data in a test and training set.
End of explanation
model = sm.OLS(y_train, X_train)
result = model.fit()
print(result.summary())
_ = sm.qqplot(result.resid, fit=True, line="s")
result.pvalues < 0.05
predicted = result.predict(X_test)
print("r2 score: {}".format(metrics.r2_score(y_test, predicted)))
print("mse: {}".format(metrics.mean_squared_error(y_test, predicted)))
print("rmse: {}".format(np.sqrt(metrics.mean_squared_error(y_test, predicted))))
print("mae: {}".format(metrics.mean_absolute_error(y_test, predicted)))
Explanation: Fitting models - the standard way
Full model without an intercept
End of explanation
model = sm.OLS(y_train, sm.add_constant(X_train))
result = model.fit()
print(result.summary())
_ = sm.qqplot(result.resid, fit=True, line="s")
predicted = result.predict(sm.add_constant(X_test))
print("r2 score: {}".format(metrics.r2_score(y_test, predicted)))
print("mse: {}".format(metrics.mean_squared_error(y_test, predicted)))
print("rmse: {}".format(np.sqrt(metrics.mean_squared_error(y_test, predicted))))
print("mae: {}".format(metrics.mean_absolute_error(y_test, predicted)))
Explanation: Full model with an intercept
End of explanation
dat = X_train.copy()
dat['MEDV'] = y_train
dat.head()
Explanation: Fitting models using R-style formulas
We can also fit a model with the R syntax y ~ x_1 + x_2 and build some complexer models.
End of explanation
result = smf.ols('MEDV ~ CRIM + ZN + INDUS + CHAS + NOX + RM + AGE + DIS + RAD + TAX + PTRATIO + B', data=dat).fit()
print(result.summary())
_ = sm.qqplot(result.resid, fit=True, line="s")
predicted = result.predict(X_test)
print("r2 score: {}".format(metrics.r2_score(y_test, predicted)))
print("mse: {}".format(metrics.mean_squared_error(y_test, predicted)))
print("rmse: {}".format(np.sqrt(metrics.mean_squared_error(y_test, predicted))))
print("mae: {}".format(metrics.mean_absolute_error(y_test, predicted)))
Explanation: Full model with an intercept
End of explanation
result = smf.ols('np.log(MEDV) ~ CRIM + CHAS + NOX + RM + DIS + RAD + TAX + PTRATIO + B + pow(AGE, 2)', data=dat).fit()
print(result.summary())
Explanation: Model with a polynomial and the target variable log transformed
End of explanation
result.pvalues < 0.05
predicted = np.exp(result.predict(X_test))
print("r2 score: {}".format(metrics.r2_score(y_test, predicted)))
print("mse: {}".format(metrics.mean_squared_error(y_test, predicted)))
print("rmse: {}".format(np.sqrt(metrics.mean_squared_error(y_test, predicted))))
print("mae: {}".format(metrics.mean_absolute_error(y_test, predicted)))
_ = sm.qqplot(result.resid, fit=True, line="q")
Explanation: Let's plot the QQ-Plot for the residuals
End of explanation
<END_TASK>