رگرسیون خطی در پایتون — کدها و برنامه‌ها (بخش دوم)

1import numpy as np
2import pandas as pd
3import matplotlib.pyplot as plt
4
5data = pd.read_csv("Concrete_Data.csv")
6
7
8x = data.iloc[:,0:8]
9y = data.iloc[:,8:]

1from sklearn.cross_validation import train_test_split
2x_train,x_test,y_train,y_test = train_test_split(x,y,test_size = 0.2,random_state = 100)

1from sklearn.linear_model import LinearRegression
2lm = LinearRegression()
3lm = lm.fit(x_train,y_train)   #lm.fit(input,output)

1lm.coef_
2array([[ 0.12415357,  0.10366839,  0.093371  , -0.13429401,  0.28804259,
3         0.02065756,  0.02563037,  0.11461733]])

1coefficients = pd.concat([pd.DataFrame(x_train.columns),pd.DataFrame(np.transpose(lm.coef_))], axis = 1)

10            Cement  0.124154
21             Blast  0.103668
32           Fly Ash  0.093371
43             Water -0.134294
54  Superplasticizer  0.288043
65                CA  0.020658
76                FA  0.025630
87               Age  0.114617

1lm.intercept_
2array([-34.273527])

1y_pred = lm.predict(x_test)
2y_error = y_test - y_pred

1from sklearn.metrics import r2_score
2r2_score(y_test,y_pred)

10.62252008774048395

1import statsmodels.api as sma
2X_train = sma.add_constant(x_train) ## let's add an intercept (beta_0) to our model
3X_test = sma.add_constant(x_test)
4import statsmodels.formula.api as sm
5lm2 = sm.OLS(y_train,X_train).fit()

1lm2.summary()
2"""
3                            OLS Regression Results                            
4==============================================================================
5Dep. Variable:                    CMS   R-squared:                       0.613
6Model:                            OLS   Adj. R-squared:                  0.609
7Method:                 Least Squares   F-statistic:                     161.0
8Date:                Wed, 03 Jan 2018   Prob (F-statistic):          4.37e-162
9Time:                        21:29:10   Log-Likelihood:                -3090.4
10No. Observations:                 824   AIC:                             6199.
11Df Residuals:                     815   BIC:                             6241.
12Df Model:                           8                                         
13Covariance Type:            nonrobust                                         
14====================================================================================
15                       coef    std err          t      P>|t|      [0.025      0.975]
16------------------------------------------------------------------------------------
17const              -34.2735     29.931     -1.145      0.253     -93.025      24.478
18Cement               0.1242      0.010     13.054      0.000       0.105       0.143
19Blast                0.1037      0.011      9.229      0.000       0.082       0.126
20Fly Ash              0.0934      0.014      6.687      0.000       0.066       0.121
21Water               -0.1343      0.046     -2.947      0.003      -0.224      -0.045
22Superplasticizer     0.2880      0.102      2.810      0.005       0.087       0.489
23CA                   0.0207      0.011      1.966      0.050    2.79e-05       0.041
24FA                   0.0256      0.012      2.131      0.033       0.002       0.049
25Age                  0.1146      0.006     19.064      0.000       0.103       0.126
26==============================================================================
27Omnibus:                        3.757   Durbin-Watson:                   2.033
28Prob(Omnibus):                  0.153   Jarque-Bera (JB):                3.762
29Skew:                          -0.165   Prob(JB):                        0.152
30Kurtosis:                       2.974   Cond. No.                     1.07e+05
31==============================================================================
32Warnings:
33[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
34[2] The condition number is large, 1.07e+05. This might indicate that there are
35strong multicollinearity or other numerical problems.
36"""

1y_pred2 = lm2.predict(X_test)

1import numpy as np
2y_test = pd.to_numeric(y_test.CMS, errors='coerce')
3RSS = np.sum((y_pred2 - y_test)**2)
4y_mean = np.mean(y_test)
5TSS = np.sum((y_test - y_mean)**2)
6R2 = 1 - RSS/TSS
7R2
8
9n=X_test.shape[0]
10p=X_test.shape[1] - 1
11
12adj_rsquared = 1 - (1 - R2) * ((n - 1)/(n-p-1))
13adj_rsquared

1influence = lm2.get_influence()  
2resid_student = influence.resid_studentized_external

1 Cement  Blast  Fly Ash  Water  Superplasticizer      CA     FA    Age  \
20   540.0    0.0      0.0  162.0               2.5  1040.0  676.0   28.0   
31   540.0    0.0      0.0  162.0               2.5  1055.0  676.0   28.0   
42   332.5  142.5      0.0  228.0               0.0   932.0  594.0  270.0   
53   332.5  142.5      0.0  228.0               0.0   932.0  594.0  365.0   
64   198.6  132.4      0.0  192.0               0.0   978.4  825.5  360.0   
7   Studentized Residuals  
80               1.559672  
91              -0.917354  
102               1.057443  
113               0.637504  
124              -1.170290

1resid = pd.concat([x_train,pd.Series(resid_student,name = "Studentized Residuals")],axis = 1)
2resid.head()

1 Cement  Blast  Fly Ash  Water  Superplasticizer     CA     FA  Age  \
2649   166.8  250.2      0.0  203.5               0.0  975.6  692.6  3.0   
3     Studentized Residuals  
4649               3.161183

1resid.loc[np.absolute(resid["Studentized Residuals"]) > 3,:]
2ind = resid.loc[np.absolute(resid["Studentized Residuals"]) > 3,:].index
3ind
4Int64Index([649], dtype='int64')

1y_train.drop(ind,axis = 0,inplace = True)
2x_train.drop(ind,axis = 0,inplace = True)  #Interept column is not there
3X_train.drop(ind,axis = 0,inplace = True)  #Intercept column is there

1from statsmodels.stats.outliers_influence import variance_inflation_factor
2[variance_inflation_factor(x_train.values, j) for j in range(x_train.shape[1])]

1[15.477582601956859,
2 3.2696650121931814,
3 4.1293255012993439,
4 82.210084751631086,
5 5.21853674386234,
6 85.866945489015535,
7 71.816336942930675,
8 1.6861600968467656]

1def calculate_vif(x):
2    thresh = 5.0
3    output = pd.DataFrame()
4    k = x.shape[1]
5    vif = [variance_inflation_factor(x.values, j) for j in range(x.shape[1])]
6    for i in range(1,k):
7        print("Iteration no.")
8        print(i)
9        print(vif)
10        a = np.argmax(vif)
11        print("Max VIF is for variable no.:")
12        print(a)
13        if vif[a] <= thresh :
14            break
15        if i == 1 :          
16            output = x.drop(x.columns[a], axis = 1)
17            vif = [variance_inflation_factor(output.values, j) for j in range(output.shape[1])]
18        elif i > 1 :
19            output = output.drop(output.columns[a],axis = 1)
20            vif = [variance_inflation_factor(output.values, j) for j in range(output.shape[1])]
21    return(output)

1train_out.head()

1   Cement  Blast  Fly Ash  Superplasticizer  Age
2337   275.1    0.0    121.4               9.9   56
3384   516.0    0.0      0.0               8.2   28
4805   393.0    0.0      0.0               0.0   90
5682   183.9  122.6      0.0               0.0   28
6329   246.8    0.0    125.1              12.0    3

1x_test.head()
2x_test.drop(["Water","CA","FA"],axis = 1,inplace = True)
3x_test.head()

1 Cement  Blast  Fly Ash  Superplasticizer  Age
2173   318.8  212.5      0.0              14.3   91
3134   362.6  189.0      0.0              11.6   28
4822   322.0    0.0      0.0               0.0   28
5264   212.0    0.0    124.8               7.8    3
6479   446.0   24.0     79.0              11.6    7

1import statsmodels.api as sma
2import statsmodels.formula.api as sm
3train_out = sma.add_constant(train_out) ## let's add an intercept (beta_0) to our model
4x_test.drop(["Water","CA","FA"],axis = 1,inplace = True)
5X_test = sma.add_constant(x_test)
6lm2 = sm.OLS(y_train,train_out).fit()
7lm2.summary()

1"""
2                            OLS Regression Results                            
3==============================================================================
4Dep. Variable:                    CMS   R-squared:                       0.570
5Model:                            OLS   Adj. R-squared:                  0.567
6Method:                 Least Squares   F-statistic:                     216.3
7Date:                Wed, 10 Jan 2018   Prob (F-statistic):          6.88e-147
8Time:                        15:14:59   Log-Likelihood:                -3128.8
9No. Observations:                 823   AIC:                             6270.
10Df Residuals:                     817   BIC:                             6298.
11Df Model:                           5                                         
12Covariance Type:            nonrobust                                         
13====================================================================================
14                       coef    std err          t      P>|t|      [0.025      0.975]
15------------------------------------------------------------------------------------
16const              -11.1119      1.915     -5.803      0.000     -14.871      -7.353
17Cement               0.1031      0.005     20.941      0.000       0.093       0.113
18Blast                0.0721      0.006     12.622      0.000       0.061       0.083
19Fly Ash              0.0614      0.009      6.749      0.000       0.044       0.079
20Superplasticizer     0.7519      0.077      9.739      0.000       0.600       0.903
21Age                  0.1021      0.006     16.582      0.000       0.090       0.114
22==============================================================================
23Omnibus:                        0.870   Durbin-Watson:                   2.090
24Prob(Omnibus):                  0.647   Jarque-Bera (JB):                0.945
25Skew:                           0.039   Prob(JB):                        0.623
26Kurtosis:                       2.853   Cond. No.                     1.59e+03
27==============================================================================

1from scipy import stats
2stats.shapiro(lm2.resid)

1(0.9983407258987427, 0.6269884705543518)

1import statsmodels.stats.api as sms
2from statsmodels.compat import lzip
3name = ['F statistic', 'p-value']
4test = sms.het_goldfeldquandt(lm2.resid, lm2.model.exog)
5lzip(name, test)

1[('F statistic', 0.9903), ('p-value', 0.539)]