Wojciech Moszczyński - THE DATA SCIENCE LIBRARY https://sigmaquality.pl/tag/wojciech-moszczynski/ Wojciech Moszczyński Fri, 01 May 2020 12:30:14 +0000 pl-PL hourly 1 https://wordpress.org/?v=6.8.3 https://sigmaquality.pl/wp-content/uploads/2019/02/cropped-ryba-32x32.png Wojciech Moszczyński - THE DATA SCIENCE LIBRARY https://sigmaquality.pl/tag/wojciech-moszczynski/ 32 32 Support Vector Regression (SVR) using linear and non-linear kernels https://sigmaquality.pl/uncategorized/support-vector-regression-svr-using-linear-and-non-linear-kernels-en240120201439/ Fri, 24 Jan 2020 18:43:00 +0000 http://sigmaquality.pl/support-vector-regression-svr-using-linear-and-non-linear-kernels-en240120201439/ EN240120201439 What is the kernel method? The kernel method involves bending two-dimensional space. this is perfectly explained here: https://shapeofdata.wordpress.com/2013/05/27/kernels/ Two-dimensional data live in a two-dimensional [...]

Artykuł Support Vector Regression (SVR) using linear and non-linear kernels pochodzi z serwisu THE DATA SCIENCE LIBRARY.

]]> EN240120201439

What is the kernel method?

The kernel method involves bending two-dimensional space. this is perfectly explained here: https://shapeofdata.wordpress.com/2013/05/27/kernels/

Two-dimensional data live in a two-dimensional plane, which we can consider as a piece of paper. The nucleus is a way to place this two-dimensional plane in a space with higher dimensions. In other words, the nucleus is a function from a low-dimensional space to a higher-dimensional space.
So, for example, we can place a plane in three-dimensional space so that it curves along one axis, as in the figure below. Here, the cross-sections of the plane in three-dimensional space are parabolas.

The purpose of the kernel is to make two classes of data points, which can only be separated by a curved line in two-dimensional space, can be separated by a flat plane in three-dimensional space.

By adding more dimensions, we increase the flexibility of lines / planes / hyperplanes to move. For example, in the data set shown below on the left there is clearly no line separating the two blue points from the two green points. However, when we place this data set in three dimensions using the nucleus of the first figure, the two blue points will be raised, as on the right, so there is now another plane (drawn in red) that separates the blue from the green, shown in the lower right corner . This plane intersects the two-dimensional space on the curve shown in the lower left drawing.

In general, the kernel does it so that each plane in three-dimensional space intersects a two-dimensional plane that contains our data in a curve line, not a straight line. If we use a nucleus that places the plane in an even higher dimensional space, then the hyperplanes in these higher spaces intersect the two-dimensional space on potentially more complex curves, which gives us more flexibility in choosing the curve separating the data.

It is worth choosing the simplest possible kernel, is that, as in the case of general regression, the greater the flexibility of the model, the greater the risk of over-fitting.

Epsilon – http://kernelsvm.tripod.com/

The epsilon borders are given in green lines. Blue points represent data instances.

The larger the epsilon parameter, the more skewed the model will be flattened
SVM generalization performance (estimation accuracy) depends on the setting of good meta-parameters C, and kernel parameters.

Parameter C

Parameter C specifies the compromise between the complexity of the model (flatness) and the degree to which deviations greater than are tolerated in the optimization formula, for example, if C is too large (infinity), then the goal is to minimize only the empirical risk, without taking into account part of the model’s complexity optimization formulation.

The parameter controls the width of the sensitivity zone used for training data. The value can affect the number of helper vectors used to construct the regression function. The larger, the less selected auxiliary vectors. On the other hand, higher values ​​result in more „flat” estimates.

RBF kernel

In machine learning, the radial basis function kernel, or RBF kernel, is a popular kernel function used in various kernelized learning algorithms. In particular, it is commonly used in support vector machine classification.

In [1]:
import numpy as np
from sklearn.svm import SVR
import matplotlib.pyplot as plt

Example 1 non-linear kernels SVM SVR(kernel=’rbf’)

In [2]:
x = np.sort(6 * np.random.rand(20, 1), axis=0)
y = np.sin(x).ravel()
In [3]:
plt.figure()
plt.plot(x, y, 'ko', color ='blue', label="Original Noised Data")
plt.legend()
plt.axhline(y=0, color='black', linestyle='--', lw=0.5)
plt.axvline(x= 0, color = 'black', linestyle='--', lw=0.5)
plt.show()

Dokonujemy zniekształcenia wykresu sinusoidalnego.

In [4]:
y = (0.3 * (0.05-np.random.rand(20)))+y  #<<<< distorting element


plt.figure()
plt.plot(x, y, 'ko', color ='blue', label="Original Noised Data")
plt.legend()
plt.axhline(y=0, color='black', linestyle='--', lw=0.5)
plt.axvline(x= 0, color = 'black', linestyle='--', lw=0.5)
plt.show()
In [5]:
print(x.shape)
print(y.shape)

(20, 1)
(20,)

Model SVM non-linear kernels

SVR(kernel=’rbf’)

In [6]:
svr_rbf = SVR(kernel='rbf', C=100, gamma=0.1, epsilon=.1)

# Fit the SVM model according to the given training data.
y_pred = svr_rbf.fit(x, y).predict(x)
In [7]:
import pandas as pd

RBF = pd.DataFrame({'y_actual': y, 'y_predicted': y_pred})
RBF.head(5)
Out[7]:
y_actual y_predicted
0 0.601198 0.577846
1 0.610104 0.668037
2 0.624167 0.724347
3 0.644582 0.736527
4 0.842055 0.770481

Get parameters for this estimator.

In [8]:
from sklearn import metrics

RBF.head(50).plot()
Out[8]:
<matplotlib.axes._subplots.AxesSubplot at 0x285618667b8>
In [9]:
plt.plot(x, y, 'o', color ='blue', label ="data") 
plt.plot(x, y_pred, '--', color ='red', label ="Fitted Curve") 
plt.legend() 
plt.show()

a = svr_rbf.score(x , y , sample_weight = None)
print('Mean Squared Error:    

Mean Squared Error:    0.99

model summary

In [10]:
a = svr_rbf.score(x , y , sample_weight = None)
print('Return the coefficient of determination R^2 of the prediction.')
print('Mean Squared Error:    

Return the coefficient of determination R^2 of the prediction.
Mean Squared Error:    0.99
In [11]:
print('Get parameters for this estimator.')
svr_rbf.get_params(deep=True)

Get parameters for this estimator.
Out[11]:
{'C': 100,
 'cache_size': 200,
 'coef0': 0.0,
 'degree': 3,
 'epsilon': 0.1,
 'gamma': 0.1,
 'kernel': 'rbf',
 'max_iter': -1,
 'shrinking': True,
 'tol': 0.001,
 'verbose': False}
In [12]:
y_pred=svr_rbf.predict(x)

Example 2 non-linear kernels SVM SVR(kernel=’rbf’)

In [13]:
x = np.linspace(0, 10, num = 40)
y = 17.45 * np.sin(0.2 * x) + np.random.normal(size = 40)  
    
y = (0.01 * (0.5-np.random.rand(40)))+y  #<<<< distorting element


plt.figure()
plt.figure(figsize=(17,4))
plt.plot(x, y, 'ko', color ='blue', label="Original Noised Data")
plt.legend()
plt.show()

x.shape

<Figure size 432x288 with 0 Axes>
Out[13]:
(40,)

Adds a dimension to the x vector.

In [14]:
df = pd.DataFrame(x)
x = np.asarray(df)
x.shape
Out[14]:
(40, 1)

Model SVM non-linear kernels

SVR(kernel=’rbf’)

In [15]:
svr_rbf = SVR(kernel='rbf', C=100, gamma=0.1, epsilon=.1)

# Fit the SVM model according to the given training data.
y_pred = svr_rbf.fit(x, y).predict(x)
In [16]:
import pandas as pd

RBF = pd.DataFrame({'y_actual': y, 'y_predicted': y_pred})
RBF.head(5)
Out[16]:
y_actual y_predicted
0 -0.162622 -0.150990
1 0.529202 0.629017
2 2.555832 1.486978
3 1.900103 2.391904
4 3.695815 3.316034
In [17]:
plt.figure(figsize=(17,4))
plt.plot(x, y, 'o', color ='blue', label ="data") 
plt.plot(x, y_pred, '--', color ='red', label ="Fitted Curve") 
plt.legend() 
plt.show()

a = svr_rbf.score(x , y , sample_weight = None)
print('Mean Squared Error:    

Mean Squared Error:    0.98

model summary

In [18]:
y_pred=svr_rbf.predict(x)
a = svr_rbf.score(x , y , sample_weight = None)
print('Return the coefficient of determination R^2 of the prediction.')
print('Mean Squared Error:    
print('Get parameters for this estimator.')
svr_rbf.get_params(deep=True)

Return the coefficient of determination R^2 of the prediction.
Mean Squared Error:    0.98  
Get parameters for this estimator.
Out[18]:
{'C': 100,
 'cache_size': 200,
 'coef0': 0.0,
 'degree': 3,
 'epsilon': 0.1,
 'gamma': 0.1,
 'kernel': 'rbf',
 'max_iter': -1,
 'shrinking': True,
 'tol': 0.001,
 'verbose': False}

Example 3 non-linear kernels SVM SVR(kernel=’rbf’)

In [19]:
x = np.linspace(-10, 10, num = 300) 
y = 1.75 * np.cosh(0.168 * x) + np.random.normal(size = 300)   
    

y = (1.21 * (2.5-np.random.rand(300)))+y  #<<<< distorting element


plt.figure()
plt.figure(figsize=(17,4))
plt.plot(x, y, 'ko', color ='blue', label="Original Noised Data")
plt.legend()
plt.show()

<Figure size 432x288 with 0 Axes>

Adds a dimension to the x vector.

In [20]:
print(x.shape)
df = pd.DataFrame(x)
x = np.asarray(df)
print(x.shape)

(300,)
(300, 1)

Model SVM non-linear kernels

SVR(kernel=’rbf’)

In [21]:
svr_rbf = SVR(kernel='rbf', C=100, gamma=0.1, epsilon=0.1)

# Fit the SVM model according to the given training data.
y_pred = svr_rbf.fit(x, y).predict(x)
In [22]:
import pandas as pd

RBF = pd.DataFrame({'y_actual': y, 'y_predicted': y_pred})
RBF.head(5)
Out[22]:
y_actual y_predicted
0 6.258306 6.525779
1 6.667066 6.567551
2 5.149480 6.603958
3 6.401777 6.634948
4 8.207138 6.660507
In [23]:
plt.figure(figsize=(17,4))
plt.plot(x, y, 'o', color ='blue', label ="data") 
plt.plot(x, y_pred, '--', color ='red', linewidth=5, label ="Fitted Curve")
plt.text(-6, 7.5, r'C=100, gamma=0.1, epsilon=0.1', fontsize=28)
plt.legend() 
plt.show()

a = svr_rbf.score(x , y , sample_weight = None)
print('Mean Squared Error:    

Mean Squared Error:    0.39
In [24]:
y_pred=svr_rbf.predict(x)
a = svr_rbf.score(x , y , sample_weight = None)
print('Return the coefficient of determination R^2 of the prediction.')
print('Mean Squared Error:    
print('Get parameters for this estimator.')
svr_rbf.get_params(deep=True)

Return the coefficient of determination R^2 of the prediction.
Mean Squared Error:    0.39  
Get parameters for this estimator.
Out[24]:
{'C': 100,
 'cache_size': 200,
 'coef0': 0.0,
 'degree': 3,
 'epsilon': 0.1,
 'gamma': 0.1,
 'kernel': 'rbf',
 'max_iter': -1,
 'shrinking': True,
 'tol': 0.001,
 'verbose': False}
In [25]:
svr_rbf = SVR(kernel='rbf', C=100, gamma=6.1, epsilon=.1)
# Fit the SVM model according to the given training data
y_pred = svr_rbf.fit(x, y).predict(x)

plt.figure(figsize=(17,4))
plt.plot(x, y, 'o', color ='blue', label ="data") 
plt.plot(x, y_pred, '--', color ='red', linewidth=5, label ="Fitted Curve") 
plt.text(-6, 7.5, r'C=100, gamma=6.1, epsilon=0.1', fontsize=28)
plt.legend() 
plt.show()

a = svr_rbf.score(x , y , sample_weight = None)
print('Mean Squared Error:    

Mean Squared Error:    0.51
In [26]:
svr_rbf = SVR(kernel='rbf', C=100, gamma=0.1, epsilon=3.4)
# Fit the SVM model according to the given training data
y_pred = svr_rbf.fit(x, y).predict(x)

plt.figure(figsize=(17,4))
plt.plot(x, y, 'o', color ='blue', label ="data") 
plt.plot(x, y_pred, '--', color ='red', linewidth=5,label ="Fitted Curve") 
plt.text(-6, 7.5, r'C=100, gamma=0.1, epsilon=3.4', fontsize=28)
plt.legend() 
plt.show()

a = svr_rbf.score(x , y , sample_weight = None)
print('Mean Squared Error:    

Mean Squared Error:    0.11
In [27]:
svr_rbf = SVR(kernel='rbf', C=6100, gamma=0.1, epsilon=0.1)
# Fit the SVM model according to the given training data
y_pred = svr_rbf.fit(x, y).predict(x)

plt.figure(figsize=(17,4))
plt.plot(x, y, 'o', color ='blue', label ="data") 
plt.plot(x, y_pred, '--', color ='red', linewidth=5, label ="Fitted Curve") 
plt.text(-6, 7.5, r'C=6100, gamma=0.1, epsilon=0.1', fontsize=28)
plt.legend() 
plt.show()

a = svr_rbf.score(x , y , sample_weight = None)
print('Mean Squared Error:    

Mean Squared Error:    0.40

Użycie innego rodzaju jądra

In [28]:
svr_lin = SVR(kernel='linear', C=100, gamma='auto')
svr_poly = SVR(kernel='poly', C=100, gamma='auto', degree=3, epsilon=.1,
               coef0=1)
In [29]:
x = np.linspace(-10, 10, num = 260) 
y = 1.75 * np.exp(0.17 * x) + np.random.normal(size = 260)   

plt.figure()
plt.figure(figsize=(17,4))
plt.plot(x, y, 'ko', color ='blue', label="Original Noised Data")
plt.legend()
plt.axhline(y=0, color='black', linestyle='--', lw=0.5)
plt.axvline(x= 0, color = 'black', linestyle='--', lw=0.5)
plt.show() 

<Figure size 432x288 with 0 Axes>
In [30]:
print(x.shape)
df = pd.DataFrame(x)
x = np.asarray(df)
print(x.shape)

(260,)
(260, 1)
In [33]:
svr_rbf = SVR(kernel='rbf', C=100, gamma=0.1, epsilon=3.4)
svr_linear = SVR(kernel='linear', C=100, gamma='auto')
svr_poly = SVR(kernel='poly', C=100, gamma='auto', degree=3, epsilon=.1,
               coef0=1)
# Fit the SVM model according to the given training data
y_pred_rbf = svr_rbf.fit(x, y).predict(x)
y_pred_linear= svr_linear.fit(x, y).predict(x)
y_pred_poly = svr_poly.fit(x, y).predict(x)
In [42]:
plt.figure(figsize=(17,8))
plt.plot(x, y, 'o', color ='lightgrey', label ="data") 
plt.plot(x, y_pred_rbf, '--', color ='red', linewidth=2,label ="Fitted Curve: rbf") 
plt.plot(x, y_pred_linear,  color ='black', linewidth=2,label ="Fitted Curve: linear") 
plt.plot(x, y_pred_poly, '--', color ='blue', linewidth=2,label ="Fitted Curve: poly") 
#plt.text(-6, 7.5, r'C=100, gamma=0.1, epsilon=3.4', fontsize=28)
plt.legend() 
plt.show()

rbf = svr_rbf.score(x , y , sample_weight = None)
print('Mean Squared Error rbf:        

linear = svr_linear.score(x , y , sample_weight = None)
print('Mean Squared Error svr_linear: 

poly = svr_poly.score(x , y , sample_weight = None)
print('Mean Squared Error svr_poly:   

Mean Squared Error rbf:        0.60  
Mean Squared Error svr_linear: 0.69  
Mean Squared Error svr_poly:   0.40

Artykuł Support Vector Regression (SVR) using linear and non-linear kernels pochodzi z serwisu THE DATA SCIENCE LIBRARY.

]]> Fit curve to data_ scipy.optimize.curve_fit https://sigmaquality.pl/uncategorized/fit-curve-to-data_-scipy-optimize-curve_fit-en220120201529/ Fri, 24 Jan 2020 18:43:00 +0000 http://sigmaquality.pl/fit-curve-to-data_-scipy-optimize-curve_fit-en220120201529/ EN240120201439 EXAMPLE Generates fictitious data In [1]: import numpy as np import matplotlib.pyplot as plt from scipy.optimize import curve_fit from pylab import * x = np.linspace(0, [...]

Artykuł Fit curve to data_ scipy.optimize.curve_fit pochodzi z serwisu THE DATA SCIENCE LIBRARY.

]]> EN240120201439

EXAMPLE

Generates fictitious data

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from pylab import *

x = np.linspace(0, 10, num = 40) 
y = 3.45 * np.sin(1.2 * x) + np.random.normal(size = 40)   

plt.figure()
plt.plot(x, y, 'ko', color ='blue', label="Original Noised Data")
plt.legend()
plt.show()  

<Figure size 640x480 with 1 Axes>

You need to know in outline what equation a given graph can have. We define the equation of the curve.

In [2]:
def test(x, a, b): 
    return a * np.sin(b * x)

The scipy function „scipy.optimize.curve_fit” adopts the type of curve to which you want to fit the data (linear),
– x axis data (x table),
– y axis data (y table),
– guessing parameters (p0).

The function then returns two information:
– popt – Sine function coefficients:
– pcov – estimated parameter covariance

In [3]:
import scipy
from scipy.optimize import curve_fit

parametr, parametr_cov = curve_fit(test, x, y)
In [4]:
parametr
Out[4]:
array([3.40022034, 1.18368923])
In [5]:
print("Sine funcion coefficients:",  parametr) 
print("Covariance of coefficients:", parametr_cov) 

Sine funcion coefficients: [3.40022034 1.18368923]
Covariance of coefficients: [[0.06204388 0.00011354]
 [0.00011354 0.00018356]]

We defined functions as:
a * np.sin (b * x)
Variables a and b are coefficients of a sinusoidal function:

  • parameter [0] = a
  • parameter [1] = b

which contain actual match parameters (popt_linear) and covariance of match parameters (pcov_linear).

In [6]:
print('a: ',parametr[0])
print('b: ',parametr[1])

a:  3.400220342120811
b:  1.1836892329574118

We substitute the formula for the sine function: a e.g. sin (b x)

In [7]:
kot = (parametr[0]*(np.sin(parametr[1]*x)))
In [8]:
plt.plot(x, y, 'o', color ='blue', label ="data") 
plt.plot(x, kot, '--', color ='red', label ="Fitted Curve: 3.6*sin(1.3*x)") 
plt.legend() 
plt.show()

EXAMPLE

Generates fictitious data

In [9]:
x = np.linspace(0, 1, num = 40)  
y = 3.45 * np.exp(2.334 * x) + 2*np.random.normal(size = 40) 

plt.figure()
plt.plot(x, y, 'ko', color ='blue', label="Original Noised Data")
plt.legend()
plt.show()

Trzeba znać w zaryzie jakie równanie może mieć dany wykres. Definiujemy równanie krzywej.

In [10]:
def kot(x, a, b): 
    return a*np.exp(b*x)
In [11]:
import scipy
from scipy.optimize import curve_fit

pa, pa_cov = curve_fit(kot, x, y)
In [12]:
print("Sine funcion coefficients:",  pa) 
print("Covariance of coefficients:")
print(pa_cov)

Sine funcion coefficients: [3.6797288  2.25070778]
Covariance of coefficients:
[[ 0.08865591 -0.02837859]
 [-0.02837859  0.0096322 ]]

Podstawiamy do wzoru: a np.exp(b x)

In [13]:
print('a: 
print('b: 

a: 3.68
b: 2.25
In [14]:
 PZU = pa[0]*np.exp(pa[1]*x) 
#PZU = kot(x, pa[0], pa[1])
In [15]:
plt.plot(x, y, 'o', color ='blue', label ="data") 
plt.plot(x, PZU, '--', color ='red', label ="Fitted Curve") 
plt.legend() 
plt.show()

EXAMPLE

Generates fictitious data

In [16]:
import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(1,10,100)
y = np.linspace(5,200,100)
y_noise = 30*(np.random.ranf(100))
y += y_noise

plt.figure()
plt.plot(x, y, 'ko',color ='blue', label="Original Noised Data")
plt.legend()
plt.show()
In [17]:
def linear(x, a, b):
    return a*x + b
In [18]:
import scipy
from scipy.optimize import curve_fit

par, par_cov = scipy.optimize.curve_fit(linear, x, y)
In [19]:
print("Sine funcion coefficients:",  par) 
print("Covariance of coefficients:")
print(par_cov)

Sine funcion coefficients: [21.97503641 -3.98345159]
Covariance of coefficients:
[[ 0.09584559 -0.52715077]
 [-0.52715077  3.55935696]]

We substitute for the formula:

In [20]:
print('a: 
print('b: 

a: 21.98
b: -3.98
In [21]:
#PKP = par[0]*x + par[1]
PKP = linear(x, par[0], par[1])
In [22]:
plt.plot(x, y, 'o', color ='blue', label ="data") 
plt.plot(x, PKP, '--', color ='red', label ="Fitted Curve") 
plt.legend() 
plt.show()

EXAMPLE

Generates fictitious data

In [23]:
x = np.linspace(0, 10, num = 160) 
y = 1.45 * np.sin(1.2 * x) + np.random.normal(size = 160)   

plt.figure()
plt.plot(x, y, 'ko', color ='blue', label="Original Noised Data")
plt.legend()
plt.show()

You need to know in outline what equation a given graph can have. We define the equation of the curve.

In [24]:
def foka(x, a, b): 
    return a * np.cos(b * x)
In [25]:
import scipy
from scipy.optimize import curve_fit

pac, pac_cov = scipy.optimize.curve_fit(foka, x, y)
In [26]:
print("Sine funcion coefficients:",  pac) 
print("Covariance of coefficients:")
print(pac_cov)

Sine funcion coefficients: [1.07953123 0.96862312]
Covariance of coefficients:
[[ 0.01835107 -0.00023385]
 [-0.00023385  0.00053849]]
In [27]:
print('a: 
print('b: 

a: 1.08
b: 0.97
In [28]:
PKS = foka(x, pac[0], pac[1])
In [29]:
plt.plot(x, y, 'o', color ='blue', label ="data") 
plt.plot(x, PKS, '--', color ='red', label ="Fitted Curve") 
plt.legend() 
plt.show()

EXAMPLE

Generates fictitious data

In [30]:
x = np.linspace(-10, 10, num = 300) 
y = 1.75 * np.cosh(0.168 * x) + np.random.normal(size = 300)   

plt.figure()
plt.plot(x, y, 'ko', color ='blue', label="Original Noised Data")
plt.legend()
plt.axhline(y=0, color='black', linestyle='--', lw=0.5)
plt.axvline(x= 0, color = 'black', linestyle='--', lw=0.5)

plt.show()

Defines the matching function.

In [31]:
def vok(x, a, b):
    return a * np.cosh(x * b)
In [32]:
import scipy
from scipy.optimize import curve_fit

pan, pan_cov = scipy.optimize.curve_fit(vok, x, y)
In [33]:
print("Sine funcion coefficients:",  pan) 
print("Covariance of coefficients:")
print(pan_cov)

Sine funcion coefficients: [1.67013896 0.1780296 ]
Covariance of coefficients:
[[ 7.28891469e-03 -5.93352295e-04]
 [-5.93352295e-04  5.89205781e-05]]
In [34]:
print('a: 
print('b: 

a: 1.67
b: 0.18
In [35]:
POS = vok(x, pan[0], pan[1])
In [36]:
plt.plot(x, y, 'o', color ='blue', label ="data") 
plt.plot(x, POS, '--', color ='orange', linewidth=5.5, label ="Fitted Curve") 
plt.legend() 
plt.axhline(y=0, color='black', linestyle='--', lw=0.5)
plt.axvline(x= 0, color = 'black', linestyle='--', lw=0.5)
plt.show()

Compute tangent

In [37]:
x = np.linspace(-10, 10, num = 60) 
y = 10.45 * np.tan(10.2 * x) + np.random.normal(size = 60)   

plt.figure()
plt.plot(x, y, 'ko', color ='blue', label="Original Noised Data")
plt.legend()
plt.axhline(y=0, color='black', linestyle='--', lw=0.5)
plt.axvline(x= 0, color = 'black', linestyle='--', lw=0.5)
plt.show()

Compute hyperbolic tangent

In [38]:
x = np.linspace(-10, 10, num = 60) 
y = 15.45 * np.tanh(1.2 * x) + np.random.normal(size = 60)   

plt.figure()
plt.plot(x, y, 'ko', color ='blue', label="Original Noised Data")
plt.legend()
plt.axhline(y=0, color='black', linestyle='--', lw=0.5)
plt.axvline(x= 0, color = 'black', linestyle='--', lw=0.5)

plt.show()

Inverse sine

In [39]:
x = np.linspace(-10, 10, num = 260) 
y = 9.45 * np.arcsin(0.7 * x) + np.random.normal(size = 260)   

plt.figure()
plt.plot(x, y, 'ko', color ='blue', label="Original Noised Data")
plt.legend()
plt.axhline(y=0, color='black', linestyle='--', lw=0.5)
plt.axvline(x= 0, color = 'black', linestyle='--', lw=0.5)

plt.show() 

C:ProgramDataAnaconda3libsite-packagesipykernel_launcher.py:2: RuntimeWarning: invalid value encountered in arcsin

Inverse hyperbolic cosine

In [40]:
x = np.linspace(-10, 10, num = 260) 
y = 9.45 * np.arccosh(3.7 * x) + np.random.normal(size = 260)   

plt.figure()
plt.plot(x, y, 'ko', color ='blue', label="Original Noised Data")
plt.legend()
plt.axhline(y=0, color='black', linestyle='--', lw=0.5)
plt.axvline(x= 0, color = 'black', linestyle='--', lw=0.5)
plt.show() 

C:ProgramDataAnaconda3libsite-packagesipykernel_launcher.py:2: RuntimeWarning: invalid value encountered in arccosh

np.cosh

In [41]:
x = np.linspace(-10, 10, num = 260) 
y = 1.45 * np.cosh(0.3 * x) + np.random.normal(size = 260)   

plt.figure()
plt.plot(x, y, 'ko', color ='blue', label="Original Noised Data")
plt.legend()
plt.axhline(y=0, color='black', linestyle='--', lw=0.5)
plt.axvline(x= 0, color = 'black', linestyle='--', lw=0.5)
plt.show()

np.sinh

In [42]:
x = np.linspace(-10, 10, num = 260) 
y = 1.45 * np.sinh(0.3 * x) + np.random.normal(size = 260)   

plt.figure()
plt.plot(x, y, 'ko', color ='blue', label="Original Noised Data")
plt.legend()
plt.axhline(y=0, color='black', linestyle='--', lw=0.5)
plt.axvline(x= 0, color = 'black', linestyle='--', lw=0.5)
plt.show()

np.tanh

In [43]:
x = np.linspace(-10, 10, num = 260) 
y = 7.45 * np.tanh(0.3 * x) + np.random.normal(size = 260)   

plt.figure()
plt.plot(x, y, 'ko', color ='blue', label="Original Noised Data")
plt.legend()
plt.axhline(y=0, color='black', linestyle='--', lw=0.5)
plt.axvline(x= 0, color = 'black', linestyle='--', lw=0.5)
plt.show()

np.arcsinh

In [44]:
x = np.linspace(-10, 10, num = 260) 
y = 7.45 * np.arcsinh(0.3 * x) + np.random.normal(size = 260)   

plt.figure()
plt.plot(x, y, 'ko', color ='blue', label="Original Noised Data")
plt.legend()
plt.axhline(y=0, color='black', linestyle='--', lw=0.5)
plt.axvline(x= 0, color = 'black', linestyle='--', lw=0.5)
plt.show()

np.arccosh

In [45]:
x = np.linspace(-10, 10, num = 260) 
y = 17.45 * np.arccosh(1.3 * x) + np.random.normal(size = 260)   

plt.figure()
plt.plot(x, y, 'ko', color ='blue', label="Original Noised Data")
plt.legend()
plt.axhline(y=0, color='black', linestyle='--', lw=0.5)
plt.axvline(x= 0, color = 'black', linestyle='--', lw=0.5)
plt.show()

C:ProgramDataAnaconda3libsite-packagesipykernel_launcher.py:2: RuntimeWarning: invalid value encountered in arccosh

np.arcsinh

In [46]:
x = np.linspace(-10, 10, num = 260) 
y = 9.45 * np.arcsinh(3.7 * x) + np.random.normal(size = 260)   

plt.figure()
plt.plot(x, y, 'ko', color ='blue', label="Original Noised Data")
plt.legend()
plt.axhline(y=0, color='black', linestyle='--', lw=0.5)
plt.axvline(x= 0, color = 'black', linestyle='--', lw=0.5)
plt.show()

np.arctan

In [47]:
x = np.linspace(-10, 10, num = 260) 
y = 9.45 * np.arctan(3.7 * x) + np.random.normal(size = 260)   

plt.figure()
plt.plot(x, y, 'ko', color ='blue', label="Original Noised Data")
plt.legend()
plt.axhline(y=0, color='black', linestyle='--', lw=0.5)
plt.axvline(x= 0, color = 'black', linestyle='--', lw=0.5)
plt.show()

Inverse hyperbolic tangent

In [48]:
x = np.linspace(-10, 10, num = 260) 
y = 59.45 * np.arctanh(0.7 * x) + np.random.normal(size = 260)   

plt.figure()
plt.plot(x, y, 'ko', color ='blue', label="Original Noised Data")
plt.legend()
plt.axhline(y=0, color='black', linestyle='--', lw=0.5)
plt.axvline(x= 0, color = 'black', linestyle='--', lw=0.5)
plt.show() 

C:ProgramDataAnaconda3libsite-packagesipykernel_launcher.py:2: RuntimeWarning: invalid value encountered in arctanh

rad2deg

In [49]:
x = np.linspace(-10, 10, num = 260) 
y = 100.45 * np.rad2deg(2.7 * x) + np.random.normal(size = 260)   

plt.figure()
plt.plot(x, y, 'ko', color ='blue', label="Original Noised Data")
plt.legend()
plt.axhline(y=0, color='black', linestyle='--', lw=0.5)
plt.axvline(x= 0, color = 'black', linestyle='--', lw=0.5)
plt.show()

Unwrap

In [50]:
x = np.linspace(-10, 10, num = 260) 
y = 0.15 * np.unwrap(2.7 * x) + np.random.normal(size = 260)   

plt.figure()
plt.plot(x, y, 'ko', color ='blue', label="Original Noised Data")
plt.legend()
plt.axhline(y=0, color='black', linestyle='--', lw=0.5)
plt.axvline(x= 0, color = 'black', linestyle='--', lw=0.5)
plt.show()

np.exp

In [51]:
x = np.linspace(-10, 10, num = 260) 
y = 7.45 * np.exp(0.17 * x) + np.random.normal(size = 260)   

plt.figure()
plt.plot(x, y, 'ko', color ='blue', label="Original Noised Data")
plt.legend()
plt.axhline(y=0, color='black', linestyle='--', lw=0.5)
plt.axvline(x= 0, color = 'black', linestyle='--', lw=0.5)
plt.show()

np.log

In [52]:
x = np.linspace(-10, 20, num = 260) 
y = 7.45 * np.log(0.17 * x) + np.random.normal(size = 260)   

plt.figure()
plt.plot(x, y, 'ko', color ='blue', label="Original Noised Data")
plt.legend()
plt.axhline(y=0, color='black', linestyle='--', lw=0.5)
plt.axvline(x= 0, color = 'black', linestyle='--', lw=0.5)
plt.show() 

C:ProgramDataAnaconda3libsite-packagesipykernel_launcher.py:2: RuntimeWarning: invalid value encountered in log

np.log10

In [53]:
x = np.linspace(-10, 30, num = 260) 
y = 7.45 * np.log10(0.17 * x) + np.random.normal(size = 260)   

plt.figure()
plt.plot(x, y, 'ko', color ='blue', label="Original Noised Data")
plt.legend()
plt.axhline(y=0, color='black', linestyle='--', lw=0.5)
plt.axvline(x= 0, color = 'black', linestyle='--', lw=0.5)
plt.show() 

C:ProgramDataAnaconda3libsite-packagesipykernel_launcher.py:2: RuntimeWarning: invalid value encountered in log10
In [54]:
x = np.linspace(-10,30, num = 260) 
y = 7.45 * np.log10(0.17 * x) + np.random.normal(size = 260)   

plt.figure()
plt.plot(x, y, 'ko', color ='blue', label="Original Noised Data")
plt.legend()
plt.axhline(y=0, color='black', linestyle='--', lw=0.5)
plt.axvline(x= 0, color = 'black', linestyle='--', lw=0.5)
plt.show() 

C:ProgramDataAnaconda3libsite-packagesipykernel_launcher.py:2: RuntimeWarning: invalid value encountered in log10

PIERWIASTEK KWADRATOWY

In [55]:
x = np.linspace(-10, 10, num = 260) 
y = 3.45 * np.sqrt(3.17 * x) + np.random.normal(size = 260)   

plt.figure()
plt.plot(x, y, 'ko', color ='blue', label="Original Noised Data")
plt.legend()
plt.axhline(y=0, color='black', linestyle='--', lw=0.5)
plt.axvline(x= 0, color = 'black', linestyle='--', lw=0.5)
plt.show() 

C:ProgramDataAnaconda3libsite-packagesipykernel_launcher.py:2: RuntimeWarning: invalid value encountered in sqrt

PIERWIASTEK PRZESTRZENNY

In [56]:
x = np.linspace(-10, 10, num = 260) 
y = 3.45 * np.cbrt(7.17 * x) + np.random.normal(size = 260)   

plt.figure()
plt.plot(x, y, 'ko', color ='blue', label="Original Noised Data")
plt.legend()
plt.axhline(y=0, color='black', linestyle='--', lw=0.5)
plt.axvline(x= 0, color = 'black', linestyle='--', lw=0.5)
plt.show()

Return the element-wise square of the input.

In [57]:
x = np.linspace(-10, 10, num = 260) 
y = 5.45 * np.square(0.17 * x) + np.random.normal(size = 260)   

plt.figure()
plt.plot(x, y, 'ko', color ='blue', label="Original Noised Data")
plt.legend()
plt.axhline(y=0, color='black', linestyle='--', lw=0.5)
plt.axvline(x= 0, color = 'black', linestyle='--', lw=0.5)
plt.show()

Returns an element-wise indication of the sign of a number

In [58]:
x = np.linspace(-10, 10, num = 260) 
y = 2.30 * np.sign(44.17 * x) + np.random.normal(size = 260)   

plt.figure()
plt.plot(x, y, 'ko', color ='blue', label="Original Noised Data")
plt.legend()
plt.axhline(y=0, color='black', linestyle='--', lw=0.5)
plt.axvline(x= 0, color = 'black', linestyle='--', lw=0.5)
plt.show()
In [78]:
# Generate sample data
x = np.sort(5 * np.random.rand(30, 1), axis=0)
y = np.sin(x).ravel()
In [79]:
plt.figure()
plt.plot(x, y, 'ko', color ='blue', label="Original Noised Data")
plt.legend()
plt.axhline(y=0, color='black', linestyle='--', lw=0.5)
plt.axvline(x= 0, color = 'black', linestyle='--', lw=0.5)
plt.show()
In [76]:
# Generate sample data
x = np.sort(10 * np.random.rand(70, 1), axis=0)
y = np.sin(x).ravel()
In [77]:
plt.figure()
plt.plot(x, y, 'ko', color ='blue', label="Original Noised Data")
plt.legend()
plt.axhline(y=0, color='black', linestyle='--', lw=0.5)
plt.axvline(x= 0, color = 'black', linestyle='--', lw=0.5)
plt.show()

In the next part we will use the menus Support Vector Regression (SVR) using linear and non-linear kernels.

Artykuł Fit curve to data_ scipy.optimize.curve_fit pochodzi z serwisu THE DATA SCIENCE LIBRARY.

]]> ONE-WAY ANOVA analizy wpływuodległości od Paryża na przestępczość we Francji w 1833 roku https://sigmaquality.pl/statistics/one-way-anova-analizy-wplywuodleglosci-od-paryza-na-przestepczosc-we-francji-w-1833-roku-pl190120201944/ Mon, 20 Jan 2020 18:59:00 +0000 http://sigmaquality.pl/one-way-anova-analizy-wplywuodleglosci-od-paryza-na-przestepczosc-we-francji-w-1833-roku-pl190120201944/ PL190120201944 Guerry, Essay on the Moral Statistics of France (1833) Source: http://datavis.ca/gallery/guerry/guerrydat.html 1. dept wydział Num ID działu Standardowe numery dla działów, z wyjątkiem Korsyki [...]

Artykuł ONE-WAY ANOVA analizy wpływuodległości od Paryża na przestępczość we Francji w 1833 roku pochodzi z serwisu THE DATA SCIENCE LIBRARY.

]]>
PL190120201944

Guerry, Essay on the Moral Statistics of France (1833)
Source: http://datavis.ca/gallery/guerry/guerrydat.html

1. dept wydział Num ID działu Standardowe numery dla działów, z wyjątkiem Korsyki (200).
2. Region Region Char (1) Region Region Francji („N” = „Północ”, „S” = „Południe”, „E” = „Wschód”, „W” = „Zachód”, „C” = „Środek”). Korsyka ma kod „” (brak, NA)
3. Department Departament Char (25) Nazwa departamentu Nazwy departamentów są nazywane zgodnie z użyciem w 1830 roku. Zobacz depts.txt, aby uzyskać alternatywne kodowanie nazw departamentów.
4. Crime_pers Crime_pers Num Muzyka pop. za przestępstwa przeciwko osobom Źródło: A2 (Compte général, 1825-1830)
5. Crime_prop Crime_prop Num Muzyka pop. za przestępstwo przeciwko mieniu Źródło: A2 (Compte général, 1825-1830)
6. Literacy Alfabetyzacja Num Procent odczytu i zapisu Procent poborowych wojskowych, którzy potrafią czytać i pisać Źródło: A2
7. Donations Darowizny Num Darowizny dla biednych Źródło: A2 (Bulletin des lois)
8. Infants Niemowlęta Num Muzyka pop. na nielegalne urodzenie Źródło: A2 (Bureaau des Longitudes, 1817-1821)
9. Suicides Samobójstwa Num Muzyka pop. na samobójstwo Źródło: A2 (Compte général, 1827-1830)
10. MainCity Główne Miasto Char (5) Rozmiar głównego miasta Rozmiar głównego miasta („1: Sm”, „2: Med”, „3: Lg”), stosowany jako surogat dla gęstości zalewania. Duże odnosi się do górnej 10, małe do dolnej 10; wszystkie pozostałe są sklasyfikowane jako średnie. Źródło: A1
11. Wealth Bogactwo Ranga Podatek per capita od nieruchomości osobistych Indeks rankingowy oparty na podatkach od nieruchomości osobistych i ruchomych na mieszkańca Źródło: A1
12. Commerce Handel Ranga Handel i przemysł Handel i przemysł mierzony rangą liczby patentów / populacji Źródło: A1
13. Clergy Kler Ranga Dystrybucja duchowieństwa Dystrybucja duchowieństwa mierzona stopniem liczby księży katolickich w czynnej służbie / ludności Źródło: A1 (Almanach officiel du clergy, 1829)
14. Crime_parents Kryminalni rodzice Ranga Przestępstwa przeciwko rodzicom Przestępstwa przeciwko rodzicom, mierzone stopniem stosunku przestępstw przeciwko rodzicom do wszystkich przestępstw – Średnia z lat 1825–1830 Źródło: A1 (Compte général)
15. Infanticide Dzieciobójca Ranga Liczba dzieciobójstw na jednego mieszkańca Wskaźnik wskaźnikowy liczby dzieciobójstw do liczby ludności – Średnia z lat 1825–1830 Źródło: A1 (Compte général)
16. Donation_clergy Duchowieństwo Ranga Darowizny dla duchowieństwa Wskaźniki rankingowe liczby zapisów i darowizn między żywymi dla ludności – Średnia dla lat 1815–1824 Źródło: A1 (Bull. Des lois, ordunn. D’autorization)
17. Lottery Loteria Ranga Zakład per capita w Royal Lottery Ranking rankingowy wpływów z zakładów loterii królewskiej do liczby ludności — Średnia z lat 1822–1826 Źródło: A1 (Compte rendus par le ministre des finances)
18. Desertion Dezercja Ranga Dezercja wojskowa Desercja wojskowa, stosunek liczby młodych żołnierzy oskarżonych o dezercję do siły kontyngentu wojskowego, minus deficyt spowodowany niewystarczalnością dostępnych polan – średnia z lat 1825–1827. Źródło: A1 (Compte du ministere du guerre, 1829 etat V)
19. Instruction Instrukcja Ranga Rankingi instrukcji zapisane z mapy instrukcji Guerry’ego. Uwaga: jest to odwrotnie związane z umiejętnością czytania i pisania (zgodnie z definicją tutaj).
20. Prostitutes Prostytutki Num Prostytutki w Paryżu Liczba prostytutek zarejestrowanych w Paryżu w latach 1816–1834, sklasyfikowanych według departamentu ich urodzenia Źródło: Parent-Duchatelet (1836), De la prostitution en Paris
21. Distance Dystans Num Odległość do Paryża (km) Odległość każdego centroidu departamentu do centroidu Sekwany (Paryż) Źródło: obliczone na podstawie centroidów departamentów
22. Area Powierzchnia Num Obszar (1000 km^2) Źródło: Angeville (1836)
23. Pop1831 Pop1831 Num 1831 populacja ludności w 1831 roku, pochodzi z Angeville (1836), Essai sur la Statistique de la Populacja français w 1000s.

In [2]:
import pandas as pd

df = pd.read_csv('c:/1/Guerry.csv', index_col='Department')
del df['Unnamed: 0']
df.head(4)
Out[2]:
dept Region Crime_pers Crime_prop Literacy Donations Infants Suicides MainCity Wealth … Crime_parents Infanticide Donation_clergy Lottery Desertion Instruction Prostitutes Distance Area Pop1831
Department
Ain 1 E 28870 15890 37 5098 33120 35039 2:Med 73 … 71 60 69 41 55 46 13 218.372 5762 346.03
Aisne 2 N 26226 5521 51 8901 14572 12831 2:Med 22 … 4 82 36 38 82 24 327 65.945 7369 513.00
Allier 3 C 26747 7925 13 10973 17044 114121 2:Med 61 … 46 42 76 66 16 85 34 161.927 7340 298.26
Basses-Alpes 4 E 12935 7289 46 2733 23018 14238 1:Sm 76 … 70 12 37 80 32 29 2 351.399 6925 155.90

4 rows × 22 columns

Celem naszego badanie jest sprawdzenie czy istnieje istotny statystycznie wpływ miasta Paryża na poziom przestępczości przeciwko osobom (Crime_pers). Podzielę wszystkie departamenty na 4 strefy odległości od miasta Paryża (Distance) i porówna poziom przestępczości przeciwko osobom.

– h0: poziom przestępczości (średnia) przeciwko osobom jest taki sam pomiędzy strefami
– h1: poziom przestępczości (średnia) przeciwko osobom jest statycznie inny pomiędzy strefami

In [9]:
df['Class_of_Distance']=pd.qcut(df['Distance'], 4,labels=["do 120 km","120-200 km","200-280 km","280-539 km"])
In [12]:
pd.pivot_table(df, index= ['Class_of_Distance'], values= "Distance", aggfunc= ['min','max', 'count'])
Out[12]:
min max count
Distance Distance Distance
Class_of_Distance
do 120 km 0.000 119.718 22
120-200 km 126.378 199.167 21
200-280 km 202.065 283.810 21
280-539 km 291.624 539.213 22

Nazwałęm te klasy dopiero po tym jak komputer je wskazał. Innymi słowy qcut – podzielil miasteczka na 4 klasy po równo, ja dodałem odpowiednie nazwy klas. Mając klasy Tworzę odpowiednią data frame aby utworzyć One-way ANOVA.

In [14]:
df[['Crime_pers','Class_of_Distance']].head(3)
Out[14]:
Crime_pers Class_of_Distance
Department
Ain 28870 200-280 km
Aisne 26226 do 120 km
Allier 26747 120-200 km

Test ANOVA

In [16]:
import statsmodels.api as sm
from statsmodels.formula.api import ols

model_H = ols('Crime_pers ~ C(Class_of_Distance)', data=df).fit()

anova_table = sm.stats.anova_lm(model_H, typ=2)
print(anova_table)

                            sum_sq    df          F    PR(>F)
C(Class_of_Distance)  1.408371e+09   3.0  11.392976  0.000003
Residual              3.378877e+09  82.0        NaN       NaN

Interpretacja: Wartość P uzyskana z analizy ANOVA jest znacząca (P

ANOVA pokazałe, że istnieje statystycznie istotna różnice w poziomie przestępczości. Niestety ANOVA nie wskazuje, które grupy obszary istotnie różnią się od siebie.
Aby poznać pary obszarów różniące się istotnie poziomem prtzestępczości, należy przeprowadzić analizę wielokrotnego porównania par (porównanie post-hoc) za pomocą testu HSD Tukeya.

Testu HSD Tukeya

In [18]:
from statsmodels.stats.multicomp import pairwise_tukeyhsd

m_comp = pairwise_tukeyhsd(endog=df['Crime_pers'], groups=df['Class_of_Distance'], alpha=0.05)
print(m_comp)

       Multiple Comparison of Means - Tukey HSD,FWER=0.05      
===============================================================
  group1     group2     meandiff     lower      upper    reject
---------------------------------------------------------------
120-200 km 200-280 km   -6203.0   -11398.3884 -1007.6116  True 
120-200 km 280-539 km -10450.2446 -15586.2551 -5314.234   True 
120-200 km do 120 km   -1923.1537  -7059.1642 3212.8569  False 
200-280 km 280-539 km  -4247.2446  -9383.2551  888.766   False 
200-280 km do 120 km   4279.8463   -856.1642  9415.8569  False 
280-539 km do 120 km   8527.0909   3451.1527  13603.0291  True 
---------------------------------------------------------------
In [37]:
m_comp.plot_simultaneous()    # Plot group confidence intervals
plt.vlines(x=20000,ymin=-0.5,ymax=3.5, color="red", alpha=0.8, linestyle='--')
Out[37]:
<matplotlib.collections.LineCollection at 0x14963a6ecc0>

Test Tukey HSD wykazał statystycznie istotną różnicę w przestępczości pomiędzy trzeba parami obszarów0:

– „120-200 km” i „200-280 km”
– „120-200 km” i „280-590 km”
– „280-590 km” i „do 120 km”

Powyższe wyniki z Tukey HSD sugerują, że oprócz wymienionych zestawów, wszystkie inne porównania par zachowują hipotezę zerową czyli w pozostałych obszarach nie istnieją istotne różnice statystyczne w poziomie przestępczości przeciwko osobom.

Sprawdzenie spełnienia warunków ANOVA

Warunki:

  1. wartości rezydualne mają rozkład normalny (test Shapiro Wilksa)
  2. wariancje w grupach są jednorodne (test Levene lub Bartlett)
  3. obserwacje są prowadzone niezależnie od siebie

Test Levene’a Sprawdzenie jednorodności wariancji

Hipoteza zerowa : grupy z populacji mają równe wariancje.

In [20]:
PKS = pd.pivot_table(df, index = 'Department', columns = 'Class_of_Distance', values='Crime_pers')
PKS.head(4)
Out[20]:
Class_of_Distance do 120 km 120-200 km 200-280 km 280-539 km
Department
Ain NaN NaN 28870.0 NaN
Aisne 26226.0 NaN NaN NaN
Allier NaN 26747.0 NaN NaN
Ardeche NaN NaN 9474.0 NaN
Test Levene’a nie toleruje pustych wartości NaN

In [21]:
P01=PKS['do 120 km'].dropna(how='any')
P02=PKS['120-200 km'].dropna(how='any')
P03=PKS['200-280 km'].dropna(how='any')
P04=PKS['280-539 km'].dropna(how='any')

Wykresy przedstawiają statystyki poziomu przestępczości w zadanych czterech obszarach.

In [29]:
PKS.plot.kde()
Out[29]:
<matplotlib.axes._subplots.AxesSubplot at 0x149634c7da0>
In [31]:
import matplotlib.pyplot as plt
PKS.boxplot(column=['do 120 km', '120-200 km', '200-280 km', '280-539 km'], grid=False)
Out[31]:
<matplotlib.axes._subplots.AxesSubplot at 0x149637838d0>
In [30]:
import scipy.stats as stats
w,p = stats.levene(P01,P02,P03,P04)
print("Value:   ",w)
print("p-value: ",p)

Value:    1.2094373013491309
p-value:  0.31156953536619264

Ponieważ wartość P nie jest nieznacząca (p>0.05), nie odrzucamy hipotezy zerowej – czyli grupy wiekowe mają jednorodne wariancje.

Test Shapiro-Wilk Sprawdzenia normalności rozkładu reszt

Hipoteza zerowa : wartości rezydualne mają rozkład normalny.

In [25]:
import scipy.stats as stats
import numpy as np


w, p = stats.shapiro(model_H.resid)
print("Value:   ",w)
print("p-value: ",np.round(p, decimals=2))

Value:    0.9815720915794373
p-value:  0.26

Ponieważ wartość P nie jest znaczące bo większa od współczynnika ufności 0.05 (p>0.05), nie ma podstaw do odrzucenia hipotezy zerowej – wartości rezydualne mają rozkład normalny.

Przyjrzymy się bliżej resztą modelu.

In [26]:
from statsmodels.graphics.gofplots import qqplot
from matplotlib import pyplot

qqplot(model_H.resid, line='s')
pyplot.show()

Faktynie reszty zachowują prawie idealny rozkłąd normalny.

 PL080120201039 

Artykuł ONE-WAY ANOVA analizy wpływuodległości od Paryża na przestępczość we Francji w 1833 roku pochodzi z serwisu THE DATA SCIENCE LIBRARY.

]]> DIAGNOSTYKA REGRESJI LINIOWEJ część 4. Fit Plot https://sigmaquality.pl/uncategorized/diagnostyka-regresji-liniowej-czesc-4-fit-plot-pl170120201440/ Mon, 20 Jan 2020 18:33:00 +0000 http://sigmaquality.pl/diagnostyka-regresji-liniowej-czesc-4-fit-plot-pl170120201440/ PL170120201440 Guerry, Essay on the Moral Statistics of France (1833) Source: http://datavis.ca/gallery/guerry/guerrydat.html dept wydział Num ID działu Standardowe numery dla działów, z wyjątkiem Korsyki (200). [...]

Artykuł DIAGNOSTYKA REGRESJI LINIOWEJ część 4. Fit Plot pochodzi z serwisu THE DATA SCIENCE LIBRARY.

]]>

PL170120201440

Guerry, Essay on the Moral Statistics of France (1833)

Source: http://datavis.ca/gallery/guerry/guerrydat.html

  1. dept wydział Num ID działu Standardowe numery dla działów, z wyjątkiem Korsyki (200).
  2. Region Region Char (1) Region Region Francji („N” = „Północ”, „S” = „Południe”, „E” = „Wschód”, „W” = „Zachód”, „C” = „Środek”). Korsyka ma kod „” (brak, NA)
  3. Department Departament Char (25) Nazwa departamentu Nazwy departamentów są nazywane zgodnie z użyciem w 1830 roku. Zobacz depts.txt, aby uzyskać alternatywne kodowanie nazw departamentów.
  4. Crime_pers Crime_pers Num Muzyka pop. za przestępstwa przeciwko osobom Źródło: A2 (Compte général, 1825-1830)
  5. Crime_prop Crime_prop Num Muzyka pop. za przestępstwo przeciwko mieniu Źródło: A2 (Compte général, 1825-1830)
  6. Literacy Alfabetyzacja Num Procent odczytu i zapisu Procent poborowych wojskowych, którzy potrafią czytać i pisać Źródło: A2
  7. Donations Darowizny Num Darowizny dla biednych Źródło: A2 (Bulletin des lois)
  8. Infants Niemowlęta Num Muzyka pop. na nielegalne urodzenie Źródło: A2 (Bureaau des Longitudes, 1817-1821)
  9. Suicides Samobójstwa Num Muzyka pop. na samobójstwo Źródło: A2 (Compte général, 1827-1830)
  10. MainCity Główne Miasto Char (5) Rozmiar głównego miasta Rozmiar głównego miasta („1: Sm”, „2: Med”, „3: Lg”), stosowany jako surogat dla gęstości zalewania. Duże odnosi się do górnej 10, małe do dolnej 10; wszystkie pozostałe są sklasyfikowane jako średnie. Źródło: A1
  11. Wealth Bogactwo Ranga Podatek per capita od nieruchomości osobistych Indeks rankingowy oparty na podatkach od nieruchomości osobistych i ruchomych na mieszkańca Źródło: A1
  12. Commerce Handel Ranga Handel i przemysł Handel i przemysł mierzony rangą liczby patentów / populacji Źródło: A1
  13. Clergy Kler Ranga Dystrybucja duchowieństwa Dystrybucja duchowieństwa mierzona stopniem liczby księży katolickich w czynnej służbie / ludności Źródło: A1 (Almanach officiel du clergy, 1829)
  14. Crime_parents Kryminalni rodzice Ranga Przestępstwa przeciwko rodzicom Przestępstwa przeciwko rodzicom, mierzone stopniem stosunku przestępstw przeciwko rodzicom do wszystkich przestępstw – Średnia z lat 1825–1830 Źródło: A1 (Compte général)
  15. Infanticide Dzieciobójca Ranga Liczba dzieciobójstw na jednego mieszkańca Wskaźnik wskaźnikowy liczby dzieciobójstw do liczby ludności – Średnia z lat 1825–1830 Źródło: A1 (Compte général)
  16. Donation_clergy Duchowieństwo Ranga Darowizny dla duchowieństwa Wskaźniki rankingowe liczby zapisów i darowizn między żywymi dla ludności – Średnia dla lat 1815–1824 Źródło: A1 (Bull. Des lois, ordunn. D’autorization)
  17. Lottery Loteria Ranga Zakład per capita w Royal Lottery Ranking rankingowy wpływów z zakładów loterii królewskiej do liczby ludności — Średnia z lat 1822–1826 Źródło: A1 (Compte rendus par le ministre des finances)
  18. Desertion Dezercja Ranga Dezercja wojskowa Desercja wojskowa, stosunek liczby młodych żołnierzy oskarżonych o dezercję do siły kontyngentu wojskowego, minus deficyt spowodowany niewystarczalnością dostępnych polan – średnia z lat 1825–1827. Źródło: A1 (Compte du ministere du guerre, 1829 etat V)
  19. Instruction Instrukcja Ranga Rankingi instrukcji zapisane z mapy instrukcji Guerry’ego. Uwaga: jest to odwrotnie związane z umiejętnością czytania i pisania (zgodnie z definicją tutaj).
  20. Prostitutes Prostytutki Num Prostytutki w Paryżu Liczba prostytutek zarejestrowanych w Paryżu w latach 1816–1834, sklasyfikowanych według departamentu ich urodzenia Źródło: Parent-Duchatelet (1836), De la prostitution en Paris
  21. Distance Dystans Num Odległość do Paryża (km) Odległość każdego centroidu departamentu do centroidu Sekwany (Paryż) Źródło: obliczone na podstawie centroidów departamentów
  22. Area Powierzchnia Num Obszar (1000 km^2) Źródło: Angeville (1836)
  23. Pop1831 Pop1831 Num 1831 populacja ludności w 1831 roku, pochodzi z Angeville (1836), Essai sur la Statistique de la Populacja français w 1000s.
In [1]:
import pandas as pd

df = pd.read_csv('c:/1/Guerry.csv', index_col='Department')
del df['Unnamed: 0']
df.head(4)
Out[1]:
dept Region Crime_pers Crime_prop Literacy Donations Infants Suicides MainCity Wealth … Crime_parents Infanticide Donation_clergy Lottery Desertion Instruction Prostitutes Distance Area Pop1831
Department
Ain 1 E 28870 15890 37 5098 33120 35039 2:Med 73 … 71 60 69 41 55 46 13 218.372 5762 346.03
Aisne 2 N 26226 5521 51 8901 14572 12831 2:Med 22 … 4 82 36 38 82 24 327 65.945 7369 513.00
Allier 3 C 26747 7925 13 10973 17044 114121 2:Med 61 … 46 42 76 66 16 85 34 161.927 7340 298.26
Basses-Alpes 4 E 12935 7289 46 2733 23018 14238 1:Sm 76 … 70 12 37 80 32 29 2 351.399 6925 155.90

4 rows × 22 columns

In [2]:
import matplotlib.pyplot as plt

CORREL = df.corr().sort_values('Wealth')

plt.figure(figsize=(10,6))
CORREL['Wealth'].plot(kind='barh', color='red')
plt.title('Korelacja ze zmienną wynikową', fontsize=20)
plt.xlabel('Poziom korelacji')
plt.ylabel('Zmienne nezależne ciągłe')
Out[2]:
Text(0, 0.5, 'Zmienne nezależne ciągłe')

Sprawdziliśmy już w części drugiej kompletność danych oraz format zmiennych dla zbioru: Guerry.csv.

Teraz będę budował regresjię wieloczynnikową.

  • Zmienną zależną będzie poziom lokalnej zamożności mieszkańców (’Wealth’)
  • zmiennymi opisującymi będą: prostytucja (Prostitutes) – wskazująca na istnienie lokalnego nadmiaru siły roboczej i braku zajęcia, poziom darowizn na kościół (’Donation_clergy’), oraz poziom lokalnego przemysłu i handlu (’Commerce’).
In [3]:
from statsmodels.formula.api import ols
import statsmodels.api as sm

X = df[['Commerce','Prostitutes','Donation_clergy']]
y = df['Wealth']

model = sm.OLS(y, sm.add_constant(X))
model_fit = model.fit()

model_fit.summary()

C:ProgramDataAnaconda3libsite-packagesnumpycorefromnumeric.py:2389: FutureWarning: Method .ptp is deprecated and will be removed in a future version. Use numpy.ptp instead.
  return ptp(axis=axis, out=out, **kwargs)
Out[3]:
OLS Regression Results
Dep. Variable: Wealth R-squared: 0.335
Model: OLS Adj. R-squared: 0.311
Method: Least Squares F-statistic: 13.78
Date: Fri, 17 Jan 2020 Prob (F-statistic): 2.33e-07
Time: 14:39:43 Log-Likelihood: -380.69
No. Observations: 86 AIC: 769.4
Df Residuals: 82 BIC: 779.2
Df Model: 3
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
const 17.3567 5.772 3.007 0.003 5.875 28.838
Commerce 0.3788 0.095 4.001 0.000 0.190 0.567
Prostitutes -0.0099 0.004 -2.204 0.030 -0.019 -0.001
Donation_clergy 0.2605 0.092 2.846 0.006 0.078 0.443
Omnibus: 2.237 Durbin-Watson: 1.588
Prob(Omnibus): 0.327 Jarque-Bera (JB): 1.776
Skew: 0.194 Prob(JB): 0.412
Kurtosis: 2.412 Cond. No. 1.39e+03

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.39e+03. This might indicate that there are
strong multicollinearity or other numerical problems.

Przykład dobrej odpowiedzi modelu: zmienna 'Donation_clergy’

In [4]:
fig, ax = plt.subplots(figsize=(12, 8))
fig = sm.graphics.plot_fit(model_fit, "Donation_clergy", ax=ax)
In [5]:
import seaborn as sns

sns.regplot(x='Donation_clergy', y='Wealth', data=df)
Out[5]:
<matplotlib.axes._subplots.AxesSubplot at 0x16115308ac8>

Przykład złej odpowiedzi modelu: zmienna 'Prostitutes’

In [6]:
fig, ax = plt.subplots(figsize=(12, 8))
fig = sm.graphics.plot_fit(model_fit, "Prostitutes", ax=ax)
In [7]:
sns.regplot(x='Prostitutes', y='Wealth', data=df)
Out[7]:
<matplotlib.axes._subplots.AxesSubplot at 0x16115818b00>

Artykuł DIAGNOSTYKA REGRESJI LINIOWEJ część 4. Fit Plot pochodzi z serwisu THE DATA SCIENCE LIBRARY.

]]> DIAGNOSTYKA REGRESJI LINIOWEJ część 3. Partial Regression Plots https://sigmaquality.pl/uncategorized/diagnostyka-regresji-liniowej-czesc-3-partial-regression-plots-pl170120201417/ Mon, 20 Jan 2020 18:28:00 +0000 http://sigmaquality.pl/diagnostyka-regresji-liniowej-czesc-3-partial-regression-plots-pl170120201417/ PL170120201417 Guerry, Essay on the Moral Statistics of France (1833) Source: http://datavis.ca/gallery/guerry/guerrydat.html   dept wydział Num ID działu Standardowe numery dla działów, z wyjątkiem Korsyki [...]

Artykuł DIAGNOSTYKA REGRESJI LINIOWEJ część 3. Partial Regression Plots pochodzi z serwisu THE DATA SCIENCE LIBRARY.

]]> PL170120201417

Guerry, Essay on the Moral Statistics of France (1833)

Source: http://datavis.ca/gallery/guerry/guerrydat.html

 
  1. dept wydział Num ID działu Standardowe numery dla działów, z wyjątkiem Korsyki (200).
  2. Region Region Char (1) Region Region Francji („N” = „Północ”, „S” = „Południe”, „E” = „Wschód”, „W” = „Zachód”, „C” = „Środek”). Korsyka ma kod „” (brak, NA)
  3. Department Departament Char (25) Nazwa departamentu Nazwy departamentów są nazywane zgodnie z użyciem w 1830 roku. Zobacz depts.txt, aby uzyskać alternatywne kodowanie nazw departamentów.
  4. Crime_pers Crime_pers Num Muzyka pop. za przestępstwa przeciwko osobom Źródło: A2 (Compte général, 1825-1830)
  5. Crime_prop Crime_prop Num Muzyka pop. za przestępstwo przeciwko mieniu Źródło: A2 (Compte général, 1825-1830)
  6. Literacy Alfabetyzacja Num Procent odczytu i zapisu Procent poborowych wojskowych, którzy potrafią czytać i pisać Źródło: A2
  7. Donations Darowizny Num Darowizny dla biednych Źródło: A2 (Bulletin des lois)
  8. Infants Niemowlęta Num Muzyka pop. na nielegalne urodzenie Źródło: A2 (Bureaau des Longitudes, 1817-1821)
  9. Suicides Samobójstwa Num Muzyka pop. na samobójstwo Źródło: A2 (Compte général, 1827-1830)
  10. MainCity Główne Miasto Char (5) Rozmiar głównego miasta Rozmiar głównego miasta („1: Sm”, „2: Med”, „3: Lg”), stosowany jako surogat dla gęstości zalewania. Duże odnosi się do górnej 10, małe do dolnej 10; wszystkie pozostałe są sklasyfikowane jako średnie. Źródło: A1
  11. Wealth Bogactwo Ranga Podatek per capita od nieruchomości osobistych Indeks rankingowy oparty na podatkach od nieruchomości osobistych i ruchomych na mieszkańca Źródło: A1
  12. Commerce Handel Ranga Handel i przemysł Handel i przemysł mierzony rangą liczby patentów / populacji Źródło: A1
  13. Clergy Kler Ranga Dystrybucja duchowieństwa Dystrybucja duchowieństwa mierzona stopniem liczby księży katolickich w czynnej służbie / ludności Źródło: A1 (Almanach officiel du clergy, 1829)
  14. Crime_parents Kryminalni rodzice Ranga Przestępstwa przeciwko rodzicom Przestępstwa przeciwko rodzicom, mierzone stopniem stosunku przestępstw przeciwko rodzicom do wszystkich przestępstw – Średnia z lat 1825–1830 Źródło: A1 (Compte général)
  15. Infanticide Dzieciobójca Ranga Liczba dzieciobójstw na jednego mieszkańca Wskaźnik wskaźnikowy liczby dzieciobójstw do liczby ludności – Średnia z lat 1825–1830 Źródło: A1 (Compte général)
  16. Donation_clergy Duchowieństwo Ranga Darowizny dla duchowieństwa Wskaźniki rankingowe liczby zapisów i darowizn między żywymi dla ludności – Średnia dla lat 1815–1824 Źródło: A1 (Bull. Des lois, ordunn. D’autorization)
  17. Lottery Loteria Ranga Zakład per capita w Royal Lottery Ranking rankingowy wpływów z zakładów loterii królewskiej do liczby ludności — Średnia z lat 1822–1826 Źródło: A1 (Compte rendus par le ministre des finances)
  18. Desertion Dezercja Ranga Dezercja wojskowa Desercja wojskowa, stosunek liczby młodych żołnierzy oskarżonych o dezercję do siły kontyngentu wojskowego, minus deficyt spowodowany niewystarczalnością dostępnych polan – średnia z lat 1825–1827. Źródło: A1 (Compte du ministere du guerre, 1829 etat V)
  19. Instruction Instrukcja Ranga Rankingi instrukcji zapisane z mapy instrukcji Guerry’ego. Uwaga: jest to odwrotnie związane z umiejętnością czytania i pisania (zgodnie z definicją tutaj).
  20. Prostitutes Prostytutki Num Prostytutki w Paryżu Liczba prostytutek zarejestrowanych w Paryżu w latach 1816–1834, sklasyfikowanych według departamentu ich urodzenia Źródło: Parent-Duchatelet (1836), De la prostitution en Paris
  21. Distance Dystans Num Odległość do Paryża (km) Odległość każdego centroidu departamentu do centroidu Sekwany (Paryż) Źródło: obliczone na podstawie centroidów departamentów
  22. Area Powierzchnia Num Obszar (1000 km^2) Źródło: Angeville (1836)
  23. Pop1831 Pop1831 Num 1831 populacja ludności w 1831 roku, pochodzi z Angeville (1836), Essai sur la Statistique de la Populacja français w 1000s.
In [1]:
import pandas as pd

df = pd.read_csv('c:/1/Guerry.csv', index_col='Department')
del df['Unnamed: 0']
df.head(4)
Out[1]:
  dept Region Crime_pers Crime_prop Literacy Donations Infants Suicides MainCity Wealth … Crime_parents Infanticide Donation_clergy Lottery Desertion Instruction Prostitutes Distance Area Pop1831
Department                                          
Ain 1 E 28870 15890 37 5098 33120 35039 2:Med 73 … 71 60 69 41 55 46 13 218.372 5762 346.03
Aisne 2 N 26226 5521 51 8901 14572 12831 2:Med 22 … 4 82 36 38 82 24 327 65.945 7369 513.00
Allier 3 C 26747 7925 13 10973 17044 114121 2:Med 61 … 46 42 76 66 16 85 34 161.927 7340 298.26
Basses-Alpes 4 E 12935 7289 46 2733 23018 14238 1:Sm 76 … 70 12 37 80 32 29 2 351.399 6925 155.90

4 rows × 22 columns

In [2]:
import matplotlib.pyplot as plt

CORREL = df.corr().sort_values('Wealth')

plt.figure(figsize=(10,6))
CORREL['Wealth'].plot(kind='barh', color='red')
plt.title('Korelacja ze zmienną wynikową', fontsize=20)
plt.xlabel('Poziom korelacji')
plt.ylabel('Zmienne nezależne ciągłe')
Out[2]:
Text(0, 0.5, 'Zmienne nezależne ciągłe')
 

Sprawdziliśmy już w części drugiej kompletność danych oraz format zmiennych dla zbioru: Guerry.csv.

Teraz będę budował regresjię wieloczynnikową.

  • Zmienną zależną będzie poziom lokalnej zamożności mieszkańców (’Wealth’)
  • zmiennymi opisującymi będą: prostytucja (Prostitutes) – wskazująca na istnienie lokalnego nadmiaru siły roboczej i braku zajęcia, poziom darowizn na kościół (’Donation_clergy’), oraz poziom lokalnego przemysłu i handlu (’Commerce’).
In [3]:
from statsmodels.formula.api import ols
import statsmodels.api as sm

X = df[['Commerce','Prostitutes','Donation_clergy']]
y = df['Wealth']

model = sm.OLS(y, sm.add_constant(X))
model_fit = model.fit()

model_fit.summary()

C:ProgramDataAnaconda3libsite-packagesnumpycorefromnumeric.py:2389: FutureWarning: Method .ptp is deprecated and will be removed in a future version. Use numpy.ptp instead.
  return ptp(axis=axis, out=out, **kwargs)
Out[3]:
OLS Regression Results
Dep. Variable: Wealth R-squared: 0.335
Model: OLS Adj. R-squared: 0.311
Method: Least Squares F-statistic: 13.78
Date: Fri, 17 Jan 2020 Prob (F-statistic): 2.33e-07
Time: 14:15:00 Log-Likelihood: -380.69
No. Observations: 86 AIC: 769.4
Df Residuals: 82 BIC: 779.2
Df Model: 3    
Covariance Type: nonrobust    
  coef std err t P>|t| [0.025 0.975]
const 17.3567 5.772 3.007 0.003 5.875 28.838
Commerce 0.3788 0.095 4.001 0.000 0.190 0.567
Prostitutes -0.0099 0.004 -2.204 0.030 -0.019 -0.001
Donation_clergy 0.2605 0.092 2.846 0.006 0.078 0.443
Omnibus: 2.237 Durbin-Watson: 1.588
Prob(Omnibus): 0.327 Jarque-Bera (JB): 1.776
Skew: 0.194 Prob(JB): 0.412
Kurtosis: 2.412 Cond. No. 1.39e+03

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.39e+03. This might indicate that there are
strong multicollinearity or other numerical problems.

 

Teraz analizujemy wartości odstające:

In [4]:
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(18,13))
fig = sm.graphics.influence_plot(model_fit, ax=ax, criterion="cooks")
plt.axhline(y=0, color='green', linestyle='-', lw=0.6)
plt.axhline(y=2, color='r', linestyle='--', lw=0.6)
plt.axhline(y=-2, color='r', linestyle='--', lw=0.6)
plt.axvline(x= 0.2, color = 'r', linestyle='--', lw=0.6)
plt.axvline(x= 0.014, color = 'r', linestyle='--', lw=0.6)
Out[4]:
<matplotlib.lines.Line2D at 0x135fc79d160>
 

Eliminujemy Seine jako miasto bardzo zniekształacające model. Problem tego miasta można zbadać przy okazji innych badań. Wykres wpływów jest bardzo dobrym narzędziem do wyszukiwania anomalii w danych.
Teraz eliminujemy miasto Seine z danych i robimy model od nowa.

In [5]:
df=df.drop(['Seine','Seine-et-Oise'],axis=0)
X = df[['Commerce','Prostitutes','Donation_clergy']]
y = df['Wealth']

model = sm.OLS(y, sm.add_constant(X))
model_fit = model.fit()

model_fit.summary()
Out[5]:
OLS Regression Results
Dep. Variable: Wealth R-squared: 0.371
Model: OLS Adj. R-squared: 0.348
Method: Least Squares F-statistic: 15.74
Date: Fri, 17 Jan 2020 Prob (F-statistic): 3.92e-08
Time: 14:15:01 Log-Likelihood: -367.73
No. Observations: 84 AIC: 743.5
Df Residuals: 80 BIC: 753.2
Df Model: 3    
Covariance Type: nonrobust    
  coef std err t P>|t| [0.025 0.975]
const 32.2955 7.206 4.482 0.000 17.956 46.635
Commerce 0.2385 0.100 2.379 0.020 0.039 0.438
Prostitutes -0.0808 0.023 -3.497 0.001 -0.127 -0.035
Donation_clergy 0.1857 0.090 2.057 0.043 0.006 0.365
Omnibus: 0.823 Durbin-Watson: 1.630
Prob(Omnibus): 0.663 Jarque-Bera (JB): 0.901
Skew: -0.131 Prob(JB): 0.637
Kurtosis: 2.566 Cond. No. 461.

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

In [10]:
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(18,13))
fig = sm.graphics.influence_plot(model_fit, ax=ax, criterion="cooks")
plt.axhline(y=0, color='green', linestyle='-', lw=0.6)
plt.axhline(y=2, color='r', linestyle='--', lw=0.6)
plt.axhline(y=-2, color='r', linestyle='--', lw=0.6)
plt.axvline(x= 0.2, color = 'r', linestyle='--', lw=0.6)
plt.axvline(x= 0.014, color = 'r', linestyle='--', lw=0.6)
Out[10]:
<matplotlib.lines.Line2D at 0x135fd03e7b8>
 

Partial Regression Plots (Duncan plot)

Ponieważ wykonujemy regresje wielowymiarowe, nie możemy po prostu patrzeć na poszczególne wykresy dwuwymiarowe, aby rozróżnić relacje. Zamiast tego chcemy spojrzeć na zależność zmiennej zależnej i zmiennych niezależnych od innych zmiennych niezależnych. Możemy to zrobić za pomocą wykresów Partial Regression Plots zwanych również added variable plots.

In [7]:
fig, ax = plt.subplots(figsize=(12,8))
fig = sm.graphics.plot_partregress("Wealth", "Commerce", ["Commerce", "Donation_clergy"], data=df, ax=ax)
 

Na powyższym rysunku pokazano jaka jest zależność jakości gospodarki (na osi x) do poziomu bogactwa (na osi y) w kontekście darowizn na kościół.
Zauważmy, że odstają te miejscowości, które mają dużą odległość na wykresie wpływów: Lot-et-Garonne, Bouches-du-Rhone, Ardache, Vosges.

In [8]:
fig, ax = plt.subplots(figsize=(12,8))
fig = sm.graphics.plot_partregress("Wealth", "Commerce", ["Commerce", "Prostitutes"], data=df, ax=ax)
 

Tym razem analizujemy zależność jakości gospodarki (na osi x) do poziomu bogactwa (na osi y) w kontekście prostytucji (czyli poziomu emigracji kobiet do Paryża). Zauważmy, że odstają znowu te miejscowości, które mają dużą odległość na wykresie wpływów Cooka: Lot-et-Garonne, Bouches-du-Rhone (jako biedne), Ardache, Vosges (jako bogate).
Wykres Partial Regression Plots jest doskonałym narzędziem wyszukującym anomalie w danych.

Warto zauważyć, że mający dużą dzwignię na wykresie wpływów: Seine-inferieure na wykresie Duncan plot w kontekście darowizn na kościuł był jako bardzo nisko (biedny) w kontekście prostytucji ma zupełnie inne położenie. Miejscowość ta jest napewno kolejną anomalią.

Jak wyglądał by wykres Partial Regression Plots gdyby analizować tylko poziom 'Prostitutes’ bez kontekstu na 'Commerce’?

In [9]:
fix, ax = plt.subplots(figsize=(12,14))
fig = sm.graphics.plot_partregress("Wealth", "Commerce", ["Prostitutes"], data=df, ax=ax)

Artykuł DIAGNOSTYKA REGRESJI LINIOWEJ część 3. Partial Regression Plots pochodzi z serwisu THE DATA SCIENCE LIBRARY.

]]> DIAGNOSTYKA REGRESJI LINIOWEJ część 2. Wykres wpływu https://sigmaquality.pl/uncategorized/diagnostyka-regresji-liniowej-czesc-2-wykres-wplywu-pl170120201220/ Mon, 20 Jan 2020 18:19:00 +0000 http://sigmaquality.pl/diagnostyka-regresji-liniowej-czesc-2-wykres-wplywu-pl170120201220/ PL170120201220 In [1]: import pandas as pd df = pd.read_csv('c:/1/WorldHappinessReport.csv', index_col='Country') df.head(3) # 19:20 - 15 minut -21:55 Out[1]: Unnamed: 0 Region Happiness Rank Happiness [...]

Artykuł DIAGNOSTYKA REGRESJI LINIOWEJ część 2. Wykres wpływu pochodzi z serwisu THE DATA SCIENCE LIBRARY.

]]>

PL170120201220

In [1]:

import pandas as pd
df = pd.read_csv('c:/1/WorldHappinessReport.csv', index_col='Country')
df.head(3)
# 19:20 - 15 minut -21:55
Out[1]:
Unnamed: 0 Region Happiness Rank Happiness Score Economy (GDP per Capita) Family Health (Life Expectancy) Freedom Trust (Government Corruption) Generosity Dystopia Residual Year
Country
Afghanistan 0 Southern Asia 153.0 3.575 0.31982 0.30285 0.30335 0.23414 0.09719 0.36510 1.95210 2015.0
Albania 1 Central and Eastern Europe 95.0 4.959 0.87867 0.80434 0.81325 0.35733 0.06413 0.14272 1.89894 2015.0
Algeria 2 Middle East and Northern Africa 68.0 5.605 0.93929 1.07772 0.61766 0.28579 0.17383 0.07822 2.43209 2015.0
In [2]:
print(df.index)

Index([         'Afghanistan',              'Albania',              'Algeria',
                     'Angola',            'Argentina',              'Armenia',
                  'Australia',              'Austria',           'Azerbaijan',
                    'Bahrain',
       ...
       'United Arab Emirates',       'United Kingdom',        'United States',
                    'Uruguay',           'Uzbekistan',            'Venezuela',
                    'Vietnam',                'Yemen',               'Zambia',
                   'Zimbabwe'],
      dtype='object', name='Country', length=495)

Stawiamy hipotezę zerową, że bogactwo (Economy (GDP per Capita)) jest bezpośrednią przyczyną szczęścia (Happiness Score).

1. Sprawdzamy czy dane są kompletne i jaki mają format

In [3]:
del df['Unnamed: 0']
df = df.dropna(how ='any')
df.isnull().sum()
Out[3]:
Region                           0
Happiness Rank                   0
Happiness Score                  0
Economy (GDP per Capita)         0
Family                           0
Health (Life Expectancy)         0
Freedom                          0
Trust (Government Corruption)    0
Generosity                       0
Dystopia Residual                0
Year                             0
dtype: int64
In [4]:
df.dtypes
Out[4]:
Region                            object
Happiness Rank                   float64
Happiness Score                  float64
Economy (GDP per Capita)         float64
Family                           float64
Health (Life Expectancy)         float64
Freedom                          float64
Trust (Government Corruption)    float64
Generosity                       float64
Dystopia Residual                float64
Year                             float64
dtype: object

2. Sprawdzamy korelację zmiennych niezależnych ze zmienną zależną

In [5]:
CORREL = df.corr().sort_values('Economy (GDP per Capita)')
CORREL['Economy (GDP per Capita)']
Out[5]:
Happiness Rank                  -0.791931
Generosity                      -0.017458
Dystopia Residual                0.039211
Year                             0.133030
Trust (Government Corruption)    0.298720
Freedom                          0.344037
Family                           0.584572
Happiness Score                  0.787068
Health (Life Expectancy)         0.791212
Economy (GDP per Capita)         1.000000
Name: Economy (GDP per Capita), dtype: float64
In [6]:
df.columns = [ 'Region', 'HappinessRank', 'HappinessScore',
       'Economy_GDPperCapita', 'Family', 'Health_LifeExpectancy_',
       'Freedom', 'Trust_GovernmentCorruption_', 'Generosity',
       'Dystopia_Residual', 'Year']

3. Robimy model OLS

In [7]:
from statsmodels.formula.api import ols
import statsmodels.api as sm

X = df['Economy_GDPperCapita']
y = df['HappinessScore']

model = sm.OLS(y, sm.add_constant(X))
model_fit = model.fit()

print(model_fit.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:         HappinessScore   R-squared:                       0.619
Model:                            OLS   Adj. R-squared:                  0.619
Method:                 Least Squares   F-statistic:                     760.3
Date:                Fri, 17 Jan 2020   Prob (F-statistic):          4.90e-100
Time:                        11:56:55   Log-Likelihood:                -499.12
No. Observations:                 469   AIC:                             1002.
Df Residuals:                     467   BIC:                             1011.
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
========================================================================================
                           coef    std err          t      P>|t|      [0.025      0.975]
----------------------------------------------------------------------------------------
const                    3.3706      0.079     42.419      0.000       3.214       3.527
Economy_GDPperCapita     2.1585      0.078     27.573      0.000       2.005       2.312
==============================================================================
Omnibus:                        1.916   Durbin-Watson:                   2.080
Prob(Omnibus):                  0.384   Jarque-Bera (JB):                1.718
Skew:                           0.037   Prob(JB):                        0.423
Kurtosis:                       2.713   Cond. No.                         4.68
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

C:ProgramDataAnaconda3libsite-packagesnumpycorefromnumeric.py:2389: FutureWarning: Method .ptp is deprecated and will be removed in a future version. Use numpy.ptp instead.
  return ptp(axis=axis, out=out, **kwargs)
In [8]:
df.index
Out[8]:
Index(['Afghanistan', 'Albania', 'Algeria', 'Angola', 'Argentina', 'Armenia',
       'Australia', 'Austria', 'Azerbaijan', 'Bahrain',
       ...
       'United Arab Emirates', 'United Kingdom', 'United States', 'Uruguay',
       'Uzbekistan', 'Venezuela', 'Vietnam', 'Yemen', 'Zambia', 'Zimbabwe'],
      dtype='object', name='Country', length=469)
In [9]:
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(18,13))
fig = sm.graphics.influence_plot(model_fit, ax=ax, criterion="cooks")
plt.axhline(y=0, color='green', linestyle='-', lw=0.6)
plt.axhline(y=2, color='r', linestyle='--', lw=0.6)
plt.axhline(y=-2, color='r', linestyle='--', lw=0.6)
plt.axvline(x= 0.01, color = 'r', linestyle='--', lw=0.6)
plt.axvline(x= 0.014, color = 'r', linestyle='--', lw=0.6)
Out[9]:
<matplotlib.lines.Line2D at 0x21ee1dd8588>

Interpretacja:
To jest wykres resztek modelu (wartości rezydualnych). Wielkość bąbli oznacza odległość Cooka.

Odległość Cooka

Jest powszechnie stosowanym oszacowaniem wpływu punktu danych podczas przeprowadzania analizy regresji metodą najmniejszych kwadratów. [1] W praktycznej analizie zwykłych metodą najmniejszych kwadratów odległość Cooka można wykorzystać na kilka sposobów:

  • do identyfikacji wpływowych punktów danych, które są szczególnie warte sprawdzenia pod kątem ważności, lub
  • wskazania obszaru przestrzeni projektowej, w których dobrze byłoby uzyskać więcej punktów danych.

Oznacza to, że odległość Cooksa mierzy wpływ każdej obserwacji w modelu lub odpowiada na pytanie: „co by się stało, gdyby określonej obserwacji nie było w modelu”, jest to jeden ze sposobów wykrywania wartości odstających, które szczególnie wpływają na linię regresji.
Wartości odstające w modelu mogą nie być najbardziej reprezentatywne lub odpowiednie, co może prowadzić do nieprawidłowych wniosków.

H Leverage oraz dystanst Cooka mierzy wpływ poszczególnych obserwacji na model OLS.
Wpływ obserwacji można rozpatrywać w kategoriach, jak bardzo przewidywane wyniki dla innych obserwacji byłyby różne, gdyby nie uwzględniono danej obserwacji. Jeśli przewidywania są takie same z daną obserwacją lub bez niej, obserwacja nie ma wpływu na model regresji.
Powszechną powszechną zasadą jest to, że obserwacja o wartości D Cooka powyżej 1,0 ma zbyt duży wpływ.

Wpływ obserwacji jest funkcją dwóch czynników:

  1. dźwignią obserwacji (leverage) jak bardzo wartość obserwacji na zmiennej predyktora różni się od średniej zmiennej predyktora i
  2. odległością obserwacji różnica między przewidywanym wynikiem dla obserwacji a jej rzeczywistym wynikiem.
    http://onlinestatbook.com/2/regression/influential.html

INTERPRETACJA: Somalia, Mozambique oraz SomaliandRegion to bardzo wpływowe dane ponieważ bomble są wielki (czym większe ty większe odstawanie od danych wg Cooka). Na wykresie widać, że wartości rezydualne dla tych krajów bardzo odstają od średniej.

Wymienione kraje mają również wysoką dźwignię. Czyli można wysnuć wniosek że kraje te destabilizują model. To tak zwane wyjątki potwierdzające regułę.

Obserwacja z dużej odległości nie będzie miała tak dużego wpływu, jeśli jej dźwignia będzie niska. To połączenie dźwigni obserwacji i odległości determinuje jej wpływ.

Leverage

W statystyce, a zwłaszcza w analizie regresji, dźwignia jest miarą tego, jak daleko wartości zmiennych niezależnych w obserwacji są od wartości innych obserwacji.
Punkty o dużej dźwigni to ewentualne obserwacje dokonywane przy skrajnych lub odległych wartościach zmiennych niezależnych, tak że brak obserwacji sąsiednich oznacza, że ​​dopasowany model regresji przejdzie blisko tej konkretnej obserwacji.
https://en.wikipedia.org/wiki/Leverage_(statistics)

Innymi słowy mówiąc, punkt o wysokiej dźwigni odsuwa linię regresji OLS od jej w łaściwego dopasowania. Czym Wyższy poziom levarage dla danej obserwacji tym obserwacja ta bardziej destabilizuje model OLS.

In [10]:
model_fit.params[0]
Out[10]:
3.370636826209402
In [11]:
X = df['Economy_GDPperCapita']
y = df['HappinessScore']

## wykres--------------------------------------###
plt.figure()
plt.scatter(df.Economy_GDPperCapita, df.HappinessScore, c = 'grey')
plt.plot(df.Economy_GDPperCapita, model_fit.params[0] + model_fit.params[1] * df.Economy_GDPperCapita, c = 'r')
plt.xlabel('Economy_GDPperCapita')
plt.ylabel('HappinessScore')
plt.title("Linear Regression Plot")
Out[11]:
Text(0.5, 1.0, 'Linear Regression Plot')
In [12]:
plt.hist(model_fit.resid_pearson)
plt.ylabel('Count')
plt.xlabel('Normalized residuals')
Out[12]:
Text(0.5, 0, 'Normalized residuals')
In [13]:
from statsmodels.graphics.regressionplots import *

plot_leverage_resid2(model_fit)
plt.axvline(x= 4, color = 'r', linestyle='--', lw=0.9)
Out[13]:
<matplotlib.lines.Line2D at 0x21ee1eb6da0>

Guerry, Essay on the Moral Statistics of France (1833)

Source: http://datavis.ca/gallery/guerry/guerrydat.html

Num Zmienna Rodzaj Etykieta / opis / źródło

  1. dept wydział Num ID działu Standardowe numery dla działów, z wyjątkiem Korsyki (200).
  2. Region Region Char (1) Region Region Francji („N” = „Północ”, „S” = „Południe”, „E” = „Wschód”, „W” = „Zachód”, „C” = „Środek”). Korsyka ma kod „” (brak, NA)
  3. Department Departament Char (25) Nazwa departamentu Nazwy departamentów są nazywane zgodnie z użyciem w 1830 roku. Zobacz depts.txt, aby uzyskać alternatywne kodowanie nazw departamentów.
  4. Crime_pers Crime_pers Num Muzyka pop. za przestępstwa przeciwko osobom Źródło: A2 (Compte général, 1825-1830)
  5. Crime_prop Crime_prop Num Muzyka pop. za przestępstwo przeciwko mieniu Źródło: A2 (Compte général, 1825-1830)
  6. Literacy Alfabetyzacja Num Procent odczytu i zapisu Procent poborowych wojskowych, którzy potrafią czytać i pisać Źródło: A2
  7. Donations Darowizny Num Darowizny dla biednych Źródło: A2 (Bulletin des lois)
  8. Infants Niemowlęta Num Muzyka pop. na nielegalne urodzenie Źródło: A2 (Bureaau des Longitudes, 1817-1821)
  9. Suicides Samobójstwa Num Muzyka pop. na samobójstwo Źródło: A2 (Compte général, 1827-1830)
  10. MainCity Główne Miasto Char (5) Rozmiar głównego miasta Rozmiar głównego miasta („1: Sm”, „2: Med”, „3: Lg”), stosowany jako surogat dla gęstości zalewania. Duże odnosi się do górnej 10, małe do dolnej 10; wszystkie pozostałe są sklasyfikowane jako średnie. Źródło: A1
  11. Wealth Bogactwo Ranga Podatek per capita od nieruchomości osobistych Indeks rankingowy oparty na podatkach od nieruchomości osobistych i ruchomych na mieszkańca Źródło: A1
  12. Commerce Handel Ranga Handel i przemysł Handel i przemysł mierzony rangą liczby patentów / populacji Źródło: A1
  13. Clergy Kler Ranga Dystrybucja duchowieństwa Dystrybucja duchowieństwa mierzona stopniem liczby księży katolickich w czynnej służbie / ludności Źródło: A1 (Almanach officiel du clergy, 1829)
  14. Crime_parents Kryminalni rodzice Ranga Przestępstwa przeciwko rodzicom Przestępstwa przeciwko rodzicom, mierzone stopniem stosunku przestępstw przeciwko rodzicom do wszystkich przestępstw – Średnia z lat 1825–1830 Źródło: A1 (Compte général)
  15. Infanticide Dzieciobójca Ranga Liczba dzieciobójstw na jednego mieszkańca Wskaźnik wskaźnikowy liczby dzieciobójstw do liczby ludności – Średnia z lat 1825–1830 Źródło: A1 (Compte général)
  16. Donation_clergy Duchowieństwo Ranga Darowizny dla duchowieństwa Wskaźniki rankingowe liczby zapisów i darowizn między żywymi dla ludności – Średnia dla lat 1815–1824 Źródło: A1 (Bull. Des lois, ordunn. D’autorization)
  17. Lottery Loteria Ranga Zakład per capita w Royal Lottery Ranking rankingowy wpływów z zakładów loterii królewskiej do liczby ludności — Średnia z lat 1822–1826 Źródło: A1 (Compte rendus par le ministre des finances)
  18. Desertion Dezercja Ranga Dezercja wojskowa Desercja wojskowa, stosunek liczby młodych żołnierzy oskarżonych o dezercję do siły kontyngentu wojskowego, minus deficyt spowodowany niewystarczalnością dostępnych polan – średnia z lat 1825–1827. Źródło: A1 (Compte du ministere du guerre, 1829 etat V)
  19. Instruction Instrukcja Ranga Rankingi instrukcji zapisane z mapy instrukcji Guerry’ego. Uwaga: jest to odwrotnie związane z umiejętnością czytania i pisania (zgodnie z definicją tutaj).
  20. Prostitutes Prostytutki Num Prostytutki w Paryżu Liczba prostytutek zarejestrowanych w Paryżu w latach 1816–1834, sklasyfikowanych według departamentu ich urodzenia Źródło: Parent-Duchatelet (1836), De la prostitution en Paris
  21. Distance Dystans Num Odległość do Paryża (km) Odległość każdego centroidu departamentu do centroidu Sekwany (Paryż) Źródło: obliczone na podstawie centroidów departamentów
  22. Area Powierzchnia Num Obszar (1000 km^2) Źródło: Angeville (1836)
  23. Pop1831 Pop1831 Num 1831 populacja ludności w 1831 roku, pochodzi z Angeville (1836), Essai sur la Statistique de la Populacja français w 1000s.

Tworzymy model analizujący ilość czynów kriminalnych (Crime_pers + Crime_prop) na podstawie poziomu wykształcenia (Literacy) oraz liczbą księży katolickich w czynnej służbie (Clergy) oraz poziomu prostytucji (Prostitutes).

In [14]:
df2 = pd.read_csv('c:/1/Guerry.csv', index_col='Department')
del df2['Unnamed: 0']
print(df2.index)
df2.head(4)

Index(['Ain', 'Aisne', 'Allier', 'Basses-Alpes', 'Hautes-Alpes', 'Ardeche',
       'Ardennes', 'Ariege', 'Aube', 'Aude', 'Aveyron', 'Bouches-du-Rhone',
       'Calvados', 'Cantal', 'Charente', 'Charente-Inferieure', 'Cher',
       'Correze', 'Cote-d'Or', 'Cotes-du-Nord', 'Creuse', 'Dordogne', 'Doubs',
       'Drome', 'Eure', 'Eure-et-Loir', 'Finistere', 'Gard', 'Haute-Garonne',
       'Gers', 'Gironde', 'Herault', 'Ille-et-Vilaine', 'Indre',
       'Indre-et-Loire', 'Isere', 'Jura', 'Landes', 'Loir-et-Cher', 'Loire',
       'Haute-Loire', 'Loire-Inferieure', 'Loiret', 'Lot', 'Lot-et-Garonne',
       'Lozere', 'Maine-et-Loire', 'Manche', 'Marne', 'Haute-Marne', 'Mayenne',
       'Meurthe', 'Meuse', 'Morbihan', 'Moselle', 'Nievre', 'Nord', 'Oise',
       'Orne', 'Pas-de-Calais', 'Puy-de-Dome', 'Basses-Pyrenees',
       'Hautes-Pyrenees', 'Pyrenees-Orientales', 'Bas-Rhin', 'Haut-Rhin',
       'Rhone', 'Haute-Saone', 'Saone-et-Loire', 'Sarthe', 'Seine',
       'Seine-Inferieure', 'Seine-et-Marne', 'Seine-et-Oise', 'Deux-Sevres',
       'Somme', 'Tarn', 'Tarn-et-Garonne', 'Var', 'Vaucluse', 'Vendee',
       'Vienne', 'Haute-Vienne', 'Vosges', 'Yonne', 'Corse'],
      dtype='object', name='Department')
Out[14]:
dept Region Crime_pers Crime_prop Literacy Donations Infants Suicides MainCity Wealth … Crime_parents Infanticide Donation_clergy Lottery Desertion Instruction Prostitutes Distance Area Pop1831
Department
Ain 1 E 28870 15890 37 5098 33120 35039 2:Med 73 … 71 60 69 41 55 46 13 218.372 5762 346.03
Aisne 2 N 26226 5521 51 8901 14572 12831 2:Med 22 … 4 82 36 38 82 24 327 65.945 7369 513.00
Allier 3 C 26747 7925 13 10973 17044 114121 2:Med 61 … 46 42 76 66 16 85 34 161.927 7340 298.26
Basses-Alpes 4 E 12935 7289 46 2733 23018 14238 1:Sm 76 … 70 12 37 80 32 29 2 351.399 6925 155.90

4 rows × 22 columns

1. Sprawdzamy kompletność i format danych

In [15]:
df2.isnull().sum()
Out[15]:
dept               0
Region             1
Crime_pers         0
Crime_prop         0
Literacy           0
Donations          0
Infants            0
Suicides           0
MainCity           0
Wealth             0
Commerce           0
Clergy             0
Crime_parents      0
Infanticide        0
Donation_clergy    0
Lottery            0
Desertion          0
Instruction        0
Prostitutes        0
Distance           0
Area               0
Pop1831            0
dtype: int64
In [16]:
df2.dtypes
Out[16]:
dept                 int64
Region              object
Crime_pers           int64
Crime_prop           int64
Literacy             int64
Donations            int64
Infants              int64
Suicides             int64
MainCity            object
Wealth               int64
Commerce             int64
Clergy               int64
Crime_parents        int64
Infanticide          int64
Donation_clergy      int64
Lottery              int64
Desertion            int64
Instruction          int64
Prostitutes          int64
Distance           float64
Area                 int64
Pop1831            float64
dtype: object

Dane nadają się do wykorzystania w procesie budowy modelu regresji liniowej.

2. Tworzymy połączoną zmienną wynikową Crime = Crime_pers + Crime_prop jako wskaźnik skorygowany wielkością lokalnej populacji

In [17]:
df2['Crime'] = ((df2['Crime_pers'] + df2['Crime_prop'])/df2['Pop1831'])
df2[['Crime','Crime_pers', 'Crime_prop','Pop1831']].sample(3)
Out[17]:
Crime Crime_pers Crime_prop Pop1831
Department
Correze 95.685649 15262 12949 294.83
Hautes-Alpes 198.776143 17488 8174 129.10
Lozere 97.613110 7710 5990 140.35

3. Analizujemy poziom korelacji zmiennych

In [18]:
CORREL = df2.corr().sort_values('Crime')

plt.figure(figsize=(10,6))
CORREL['Crime'].plot(kind='barh', color='red')
plt.title('Korelacja ze zmienną wynikową', fontsize=20)
plt.xlabel('Poziom korelacji')
plt.ylabel('Zmienne nezależne ciągłe')
Out[18]:
Text(0, 0.5, 'Zmienne nezależne ciągłe')

4. Budujemy model regresji liniowej

In [19]:
X = df2['Pop1831']
y = df2['Crime']

model = sm.OLS(y, sm.add_constant(X))
model_fit = model.fit()

print(model_fit.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                  Crime   R-squared:                       0.348
Model:                            OLS   Adj. R-squared:                  0.340
Method:                 Least Squares   F-statistic:                     44.83
Date:                Fri, 17 Jan 2020   Prob (F-statistic):           2.28e-09
Time:                        11:56:57   Log-Likelihood:                -409.35
No. Observations:                  86   AIC:                             822.7
Df Residuals:                      84   BIC:                             827.6
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const        133.7714      8.470     15.794      0.000     116.928     150.614
Pop1831       -0.1395      0.021     -6.696      0.000      -0.181      -0.098
==============================================================================
Omnibus:                       25.716   Durbin-Watson:                   1.734
Prob(Omnibus):                  0.000   Jarque-Bera (JB):               53.099
Skew:                           1.080   Prob(JB):                     2.95e-12
Kurtosis:                       6.187   Cond. No.                     1.12e+03
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.12e+03. This might indicate that there are
strong multicollinearity or other numerical problems.

C:ProgramDataAnaconda3libsite-packagesnumpycorefromnumeric.py:2389: FutureWarning: Method .ptp is deprecated and will be removed in a future version. Use numpy.ptp instead.
  return ptp(axis=axis, out=out, **kwargs)

5.

In [20]:
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(18,13))
fig = sm.graphics.influence_plot(model_fit, ax=ax, criterion="cooks")
plt.axhline(y=0, color='green', linestyle='-', lw=0.6)
plt.axhline(y=2, color='r', linestyle='--', lw=0.6)
plt.axhline(y=-2, color='r', linestyle='--', lw=0.6)
plt.axvline(x= 0.05, color = 'r', linestyle='--', lw=0.6)
plt.axvline(x= 0.014, color = 'r', linestyle='--', lw=0.6)
Out[20]:
<matplotlib.lines.Line2D at 0x21ee1ff1668>
In [21]:
from statsmodels.graphics.regressionplots import *

plot_leverage_resid2(model_fit)
plt.axvline(x= 3.5, color = 'r', linestyle='--', lw=0.9)
plt.axhline(y=0.15, color='r', linestyle='--', lw=0.6)
Out[21]:
<matplotlib.lines.Line2D at 0x21ee203f390>
Interpretacja: na wykresach widać, że kilka miast mocno destabilizuje model OLS (wysoka dźwignia i wysoka odległość Cooka).

5. Usuwamy trzy obserwacje destabilizujące proces: 'Nord’, 'Creuse’, 'Hautes-Alpes’, 'Corse’, 'Seine’, 'Indre’, 'Ardennes’

In [22]:
df3=df2.drop(['Nord', 'Creuse', 'Hautes-Alpes', 'Corse', 'Seine', 'Indre', 'Ardennes'],axis=0)
In [23]:
X = df3['Pop1831']
y = df3['Crime']

model2 = sm.OLS(y, sm.add_constant(X))
model_fit2 = model2.fit()

print(model_fit.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                  Crime   R-squared:                       0.348
Model:                            OLS   Adj. R-squared:                  0.340
Method:                 Least Squares   F-statistic:                     44.83
Date:                Fri, 17 Jan 2020   Prob (F-statistic):           2.28e-09
Time:                        11:56:58   Log-Likelihood:                -409.35
No. Observations:                  86   AIC:                             822.7
Df Residuals:                      84   BIC:                             827.6
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const        133.7714      8.470     15.794      0.000     116.928     150.614
Pop1831       -0.1395      0.021     -6.696      0.000      -0.181      -0.098
==============================================================================
Omnibus:                       25.716   Durbin-Watson:                   1.734
Prob(Omnibus):                  0.000   Jarque-Bera (JB):               53.099
Skew:                           1.080   Prob(JB):                     2.95e-12
Kurtosis:                       6.187   Cond. No.                     1.12e+03
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.12e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
In [24]:
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(18,13))
fig = sm.graphics.influence_plot(model_fit2, ax=ax, criterion="cooks")
plt.axhline(y=0, color='green', linestyle='-', lw=0.6)
plt.axhline(y=2, color='r', linestyle='--', lw=0.6)
plt.axhline(y=-2, color='r', linestyle='--', lw=0.6)
plt.axvline(x= 0.05, color = 'r', linestyle='--', lw=0.6)
plt.axvline(x= 0.014, color = 'r', linestyle='--', lw=0.6)
Out[24]:
<matplotlib.lines.Line2D at 0x21ee20d0710>

Model w zakresie R^2 nieznacznie się poprawił.

In [34]:
X = df2['Pop1831']
y = df2['Crime']

## wykres--------------------------------------###
plt.figure()
plt.scatter(df2.Pop1831, df2.Crime, c = 'r')
plt.plot(df2.Pop1831, model_fit.params[0] + model_fit.params[1] * df2.Pop1831, c = 'r')
plt.scatter(df3.Pop1831, df3.Crime, c = 'green')
plt.plot(df3.Pop1831, model_fit2.params[0] + model_fit2.params[1] * df3.Pop1831, c = 'green')
plt.xlabel('Local population')
plt.ylabel('Crime')
plt.title("Regression")
Out[34]:
Text(0.5, 1.0, 'Regression')

Na wykresie

  • czerwona linia i czerwone punkty symbolizują zbiór danych przed eliminacją obserwacji odstających.
  • punkty i linia zielona pokazuje zbiór danych i linie regresji po eliminacji odstających punktów.

Artykuł DIAGNOSTYKA REGRESJI LINIOWEJ część 2. Wykres wpływu pochodzi z serwisu THE DATA SCIENCE LIBRARY.

]]> Przykład klasyfikacji wykonanej za pomocą regresji logistycznej. Eliminacja współliniowości zmiennych niezależnych za pomocą VIF https://sigmaquality.pl/uncategorized/przyklad-klasyfikacji-wykonanej-za-pomoca-regresji-logistycznej-eliminacja-wspolliniowosci-zmiennych-niezaleznych-za-pomoca-vif-pl140120202018/ Tue, 14 Jan 2020 19:21:24 +0000 http://sigmaquality.pl/przyklad-klasyfikacji-wykonanej-za-pomoca-regresji-logistycznej-eliminacja-wspolliniowosci-zmiennych-niezaleznych-za-pomoca-vif-pl140120202018/ PL140120202018 Rozważania teoretyczne na temat znaczenia (szkodliwości) wspołiniowości zmiennych niezależnych (multicollinearity) zostały wyjaśnione w poprzedniej części. Problem występowania współliniowości zmiennych niezależnych w regresji logistycznej. Przykład [...]

Artykuł Przykład klasyfikacji wykonanej za pomocą regresji logistycznej. Eliminacja współliniowości zmiennych niezależnych za pomocą VIF pochodzi z serwisu THE DATA SCIENCE LIBRARY.

]]> PL140120202018

Rozważania teoretyczne na temat znaczenia (szkodliwości) wspołiniowości zmiennych niezależnych (multicollinearity) zostały wyjaśnione w poprzedniej części.

Problem występowania współliniowości zmiennych niezależnych w regresji logistycznej. Przykład użycia testu VIF.

Tym razem będziemy analizowali zagadnienie zanieczyszczenia miasta.
Źródło danych można znaleźć tutaj.

In [1]:
import pandas as pd
df = pd.read_csv('c:/TF/AirQ_filled.csv')
df.head(3)
Out[1]:
Unnamed: 0 Date Time CO(GT) PT08.S1(CO) C6H6(GT) PT08.S2(NMHC) NOx(GT) PT08.S3(NOx) NO2(GT) PT08.S4(NO2) PT08.S5(O3) T RH AH
0 0 10/03/2004 18.00.00 2.6 1360.0 11.9 1046.0 166.0 1056.0 113.0 1692.0 1268.0 13.6 48.9 0.7578
1 1 10/03/2004 19.00.00 2.0 1292.0 9.4 955.0 103.0 1174.0 92.0 1559.0 972.0 13.3 47.7 0.7255
2 2 10/03/2004 20.00.00 2.2 1402.0 9.0 939.0 131.0 1140.0 114.0 1555.0 1074.0 11.9 54.0 0.7502

Na początku sprawdzimy kompletność danych oraz ich format.

In [2]:
del df['Unnamed: 0']
df.isnull().sum()
Out[2]:
Date             0
Time             0
CO(GT)           0
PT08.S1(CO)      0
C6H6(GT)         0
PT08.S2(NMHC)    0
NOx(GT)          0
PT08.S3(NOx)     0
NO2(GT)          0
PT08.S4(NO2)     0
PT08.S5(O3)      0
T                0
RH               0
AH               0
dtype: int64
In [3]:
df.dtypes
Out[3]:
Date              object
Time              object
CO(GT)           float64
PT08.S1(CO)      float64
C6H6(GT)         float64
PT08.S2(NMHC)    float64
NOx(GT)          float64
PT08.S3(NOx)     float64
NO2(GT)          float64
PT08.S4(NO2)     float64
PT08.S5(O3)      float64
T                float64
RH               float64
AH               float64
dtype: object

Dane są kompletne i mają właściwy format do prowadzenia dalszej analizy.

Zadanie

Naszym zadaniem jest podział na trzy klasy zanieczyszczenia oznaczonego jako PT08.S5(03). Następnie musimy zbudować model klasyfikacji logistycznej oparty na pozostałych zmiennych w postaci poziomów zanieczyszczeń.

Korelacja zmiennych niezależnych ze zmienną zależną.

In [4]:
CORREL = df.corr()

import seaborn as sns
sns.heatmap(CORREL, annot=True, cbar=False, cmap="coolwarm")
Out[4]:
<matplotlib.axes._subplots.AxesSubplot at 0x243da06fb38>

Powyższa macierz wskazuje, że istnieje współliniowość, która negatywnie wpłynie na poziom klasyfikacji. Nie zważając na to kontynujemy realizację naszego zadania. Poniżej przeprowadziłem analizę korelacji zmiennej zależnej: 'PT08.S5(O3)’ ze zmiennymi niezależnymi. W wiekszości występuje wysoki poziom korelacji co jest zjawiskiem pozytywnym.

In [5]:
CORREL2 = df.corr().sort_values('PT08.S5(O3)')
CORREL2['PT08.S5(O3)'].plot(kind='bar')
Out[5]:
<matplotlib.axes._subplots.AxesSubplot at 0x243da5f1400>
In [6]:
df['Categores_PT08.S5(O3)'] = pd.qcut(df['PT08.S5(O3)'],2)
df['Categores_PT08.S5(O3)'].value_counts().to_frame()
Out[6]:
Categores_PT08.S5(O3)
(220.999, 964.0] 4683
(964.0, 2523.0] 4674

Badana kategoria PT08.S5(O3) została podzielona na równe przedziały. Mniejsza o to czy ten podział jest uzasadniony chemicznie. Najważniejsze że teraz możemy utworzyć model klasyfikacji oparty na regresji logistycznej.
Równy podział zbiorów oddala nas od negatywnego zjawiska niezbilansowania zbiorów i koniczności tworzenia oversampling. Dzisiaj nie to jest najważniejsze lecz przećwiczenie radzenia sobie ze zjawiskiej multicollinearity.

Zbudjemy model klasyfikacji tak jakby nie dotyczył nas problem multicollinearity.

Wskazujemy zmienne niezależne:
In [7]:
KOT = df[['CO(GT)', 'PT08.S1(CO)', 'C6H6(GT)', 'PT08.S2(NMHC)',
       'NOx(GT)', 'PT08.S3(NOx)', 'NO2(GT)', 'PT08.S4(NO2)', 
       'T', 'RH', 'AH','Categores_PT08.S5(O3)']]
KOT.columns = ['CO_GT', 'PT08_S1_CO', 'C6H6_GT', 'PT08_S2_NMHC',
       'NOx_GT', 'PT08_S3_NOx', 'NO2_GT', 'PT08_S4_NO2', 
       'T', 'RH', 'AH','Categores_PT08_S5_O3']

Przekształcamy zmienną zależną: Categores_PT08.S5(O3) na zmienną ciągłą.

In [8]:
KOT['Categores_PT08_S5_O3'] = KOT['Categores_PT08_S5_O3'].astype(str)
KOT['Categores_PT08_S5_O3'].dtypes

C:ProgramDataAnaconda3libsite-packagesipykernel_launcher.py:1: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  """Entry point for launching an IPython kernel.
Out[8]:
dtype('O')
In [9]:
from sklearn.model_selection import train_test_split 


y = KOT['Categores_PT08_S5_O3']
X = KOT.drop('Categores_PT08_S5_O3', axis=1)


Xtrain, Xtest, ytrain, ytest = train_test_split(X,y, test_size=0.33, stratify = y, random_state = 148)
In [10]:
KOT.dtypes
Out[10]:
CO_GT                   float64
PT08_S1_CO              float64
C6H6_GT                 float64
PT08_S2_NMHC            float64
NOx_GT                  float64
PT08_S3_NOx             float64
NO2_GT                  float64
PT08_S4_NO2             float64
T                       float64
RH                      float64
AH                      float64
Categores_PT08_S5_O3     object
dtype: object
In [11]:
print ('Zbiór X treningowy: ',Xtrain.shape)
print ('Zbiór X testowy:    ', Xtest.shape)
print ('Zbiór y treningowy: ', ytrain.shape)
print ('Zbiór y testowy:    ', ytest.shape)

Zbiór X treningowy:  (6269, 11)
Zbiór X testowy:     (3088, 11)
Zbiór y treningowy:  (6269,)
Zbiór y testowy:     (3088,)
In [12]:
import numpy as np
from sklearn import model_selection
from sklearn.pipeline import make_pipeline
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import GridSearchCV

Parameteres = {'C': np.power(10.0, np.arange(-3, 3))}
LR = LogisticRegression(warm_start = True)
LR_Grid = GridSearchCV(LR, param_grid = Parameteres, scoring = 'roc_auc', n_jobs = 5, cv=2)

LR_Grid.fit(Xtrain, ytrain)

C:ProgramDataAnaconda3libsite-packagessklearnlinear_modellogistic.py:432: FutureWarning: Default solver will be changed to 'lbfgs' in 0.22. Specify a solver to silence this warning.
  FutureWarning)
Out[12]:
GridSearchCV(cv=2, error_score='raise-deprecating',
             estimator=LogisticRegression(C=1.0, class_weight=None, dual=False,
                                          fit_intercept=True,
                                          intercept_scaling=1, l1_ratio=None,
                                          max_iter=100, multi_class='warn',
                                          n_jobs=None, penalty='l2',
                                          random_state=None, solver='warn',
                                          tol=0.0001, verbose=0,
                                          warm_start=True),
             iid='warn', n_jobs=5,
             param_grid={'C': array([1.e-03, 1.e-02, 1.e-01, 1.e+00, 1.e+01, 1.e+02])},
             pre_dispatch='2*n_jobs', refit=True, return_train_score=False,
             scoring='roc_auc', verbose=0)

Ocena modelu regresji logistycznej bez redukcji zjawiska multicollinearity

In [13]:
ypred = LR_Grid.predict(Xtest)

from sklearn.metrics import classification_report, confusion_matrix
from sklearn import metrics

co_matrix = metrics.confusion_matrix(ytest, ypred)
co_matrix
Out[13]:
array([[1414,  131],
       [ 176, 1367]], dtype=int64)
In [14]:
print(classification_report(ytest, ypred))

                  precision    recall  f1-score   support

(220.999, 964.0]       0.89      0.92      0.90      1545
 (964.0, 2523.0]       0.91      0.89      0.90      1543

        accuracy                           0.90      3088
       macro avg       0.90      0.90      0.90      3088
    weighted avg       0.90      0.90      0.90      3088

Model wykazuje się doskonałymi zdolnościmi prognostycznymi w klasyfikacji zanieczyszczenia. Nie oznacza to jednak, że wszystko jest doskonale. Wśród zmiennych opisujących występuje wysoka multicollinearity. Jest to niekorzystne zjawiski, które teraz wyelinimujemy.

Eliniminacja zjawiska multicollinearity

Variance Inflation Factor (VIF)

Variance Inflation Factor (VIF) to miara multicollinearity między zmiennymi predykcyjnymi w regresji wielokrotnej. Określa ilościowo nasilenie multicollinearity w zwykłej analizie regresji metodą najmniejszych kwadratów. Zapewnia wskaźnik, który mierzy, o ile wariancja (kwadrat odchylenia standardowego oszacowania) szacowanego współczynnika regresji jest zwiększona z powodu kolinearności.

Kroki wdrażania VIF

  1. Uruchom regresję wieloraką.
  2. Oblicz współczynniki VIF.
  3. Sprawdź współczynniki dla każdej zmiennej predykcyjnej, jeśli VIF wynosi między 5-10, prawdopodobnie występuje multicollinearity i powinieneś rozważyć usunięcie tej zmiennej.

Czyli przed rozpoczęciem pracy nad modelem regresji logistycznej należy utworzyć model regresji wielorakiej i przez VIF wyselekcjonować, wybrać dane do modelu regresji logistycznej.

Jak wspominają tutaj: https://www.ibm.com/support/pages/multicollinearity-diagnostics-logistic-regression-nomreg-or-plum

Procedury regresji dla zmiennych zależnych kategorycznie nie mają diagnostyki kolinearności. W tym celu można jednak użyć procedury regresji liniowej. Statystyka kolinearności w regresji dotyczy relacji między predyktorami, ignorując zmienną zależną. Możesz więc uruchomić REGRESSION z tą samą listą predyktorów i zmiennej zależnej, jakiej chcesz użyć w REGRESJI LOGISTYCZNEJ (na przykład) i zażądać diagnostyki kolinearności. Uruchom regresję logistyczną, aby uzyskać właściwe współczynniki, przewidywane prawdopodobieństwa itp. Po podjęciu niezbędnych decyzji (porzucenie predyktorów itp.) Wynikających z analizy kolinearności.

In [16]:
from patsy import dmatrices
from statsmodels.stats.outliers_influence import variance_inflation_factor
import statsmodels.formula.api as smf

lm = smf.ols(formula = 'Categores_PT08_S5_O3 ~ CO_GT+PT08_S1_CO+C6H6_GT+PT08_S2_NMHC+NOx_GT+PT08_S3_NOx+NO2_GT+PT08_S4_NO2+T+RH+AH', data = KOT).fit()
y, X = dmatrices('Categores_PT08_S5_O3 ~ CO_GT+PT08_S1_CO+C6H6_GT+PT08_S2_NMHC+NOx_GT+PT08_S3_NOx+NO2_GT+PT08_S4_NO2+T+RH+AH', data = KOT, return_type = "dataframe")
vif = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
print(vif)

[648.230997747899, 11.43831319456862, 7.939972023928728, 46.44356640618427, 66.47134191824114, 5.220873142678388, 5.43267527445646, 3.714154563031008, 14.52279194616348, 14.359778937836424, 8.176237173633366, 10.619321015188527]
In [17]:
vif  =  np.round(vif, decimals=2) 
vif = list(map(float, vif))
name = list(X)

s1=pd.Series(name,name='name')
s2=pd.Series( vif,name='vif')

RFE_list = pd.concat([s1,s2], axis=1)

RFE_list
Out[17]:
name vif
0 Intercept 648.23
1 CO_GT 11.44
2 PT08_S1_CO 7.94
3 C6H6_GT 46.44
4 PT08_S2_NMHC 66.47
5 NOx_GT 5.22
6 PT08_S3_NOx 5.43
7 NO2_GT 3.71
8 PT08_S4_NO2 14.52
9 T 14.36
10 RH 8.18
11 AH 10.62

Interpretacja: wynik w postaci wektora reprezentuje zmienną w określonej kolejności jak w modelu. W zaleceniu VIF wskazano, że jeśli współczynnik przypisany do zmiennej jest większy niż 5, zmienna ta jest wysoce skorelowana z innymi zmiennymi i powinna zostać wyeliminowana z modelu.

Test wykazał, że zmienne: density oraz factorA powinny zostać usunięta z modelu.
jeszcze raz oglądamy macierz korelacji

Jeszcze raz tworzymy macierz korelacji vif i sprawdzamy zmienne po wyeliminowaniu:

  • C6H6_GT
  • PT08_S2_NMHC
  • PT08_S4_NO2Eliminuje też zmienne gdzie korelacja ze zmienną zależną było bardzo małe a zmienne miały VIF powyżej 5:
  • T
  • RH
  • AH
In [18]:
from patsy import dmatrices
from statsmodels.stats.outliers_influence import variance_inflation_factor
import statsmodels.formula.api as smf

lm = smf.ols(formula = 'Categores_PT08_S5_O3 ~ CO_GT+PT08_S1_CO+NOx_GT+PT08_S3_NOx+NO2_GT', data = KOT).fit()
y, X = dmatrices('Categores_PT08_S5_O3 ~ CO_GT+PT08_S1_CO+NOx_GT+PT08_S3_NOx+NO2_GT', data = KOT, return_type = "dataframe")
vif = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
print(vif)

[200.1579775895654, 5.895905610386816, 5.807909737315139, 3.1029777937073924, 2.7926706710126528, 2.4675790616762594]
In [19]:
vif  =  np.round(vif, decimals=2) 
vif = list(map(float, vif))
name = list(X)

s1=pd.Series(name,name='name')
s2=pd.Series( vif,name='vif')

RFE_list = pd.concat([s1,s2], axis=1)

RFE_list
Out[19]:
name vif
0 Intercept 200.16
1 CO_GT 5.90
2 PT08_S1_CO 5.81
3 NOx_GT 3.10
4 PT08_S3_NOx 2.79
5 NO2_GT 2.47

Teraz zmienne, pomimo wysokich stanów przekraczających 5, to są: CO_GT i PT08_S1_CO są akceptowane w zakresie poziomu multicollinearity. Tworzę nowy model klasyfikacji na okrojonym zbiorze zmienych opisujących.

In [20]:
KOT2 = KOT[['Categores_PT08_S5_O3', 'CO_GT', 'PT08_S1_CO','NOx_GT', 'PT08_S3_NOx', 'NO2_GT']]

KOT2.sample(3)
Out[20]:
Categores_PT08_S5_O3 CO_GT PT08_S1_CO NOx_GT PT08_S3_NOx NO2_GT
790 (220.999, 964.0] 1.4 1035.0 96.0 1222.0 82.0
2484 (220.999, 964.0] 0.9 890.0 93.0 922.0 67.0
4306 (220.999, 964.0] 0.5 828.0 72.0 1281.0 48.0

Budujemy model klasyfikacji po VIF

In [22]:
#from sklearn.model_selection import train_test_split 


y = KOT2['Categores_PT08_S5_O3']
X = KOT2.drop('Categores_PT08_S5_O3', axis=1)


Xtrain, Xtest, ytrain, ytest = train_test_split(X,y, test_size=0.33, stratify = y, random_state = 148)
In [23]:
#import numpy as np
#from sklearn import model_selection
#from sklearn.pipeline import make_pipeline
#from sklearn.linear_model import LogisticRegression
#from sklearn.model_selection import GridSearchCV

Parameteres = {'C': np.power(10.0, np.arange(-3, 3))}
LR = LogisticRegression(warm_start = True)
LR_VIF = GridSearchCV(LR, param_grid = Parameteres, scoring = 'roc_auc', n_jobs = 5, cv=2)

LR_VIF.fit(Xtrain, ytrain)

C:ProgramDataAnaconda3libsite-packagessklearnlinear_modellogistic.py:432: FutureWarning: Default solver will be changed to 'lbfgs' in 0.22. Specify a solver to silence this warning.
  FutureWarning)
Out[23]:
GridSearchCV(cv=2, error_score='raise-deprecating',
             estimator=LogisticRegression(C=1.0, class_weight=None, dual=False,
                                          fit_intercept=True,
                                          intercept_scaling=1, l1_ratio=None,
                                          max_iter=100, multi_class='warn',
                                          n_jobs=None, penalty='l2',
                                          random_state=None, solver='warn',
                                          tol=0.0001, verbose=0,
                                          warm_start=True),
             iid='warn', n_jobs=5,
             param_grid={'C': array([1.e-03, 1.e-02, 1.e-01, 1.e+00, 1.e+01, 1.e+02])},
             pre_dispatch='2*n_jobs', refit=True, return_train_score=False,
             scoring='roc_auc', verbose=0)

Ocena modelu regresji logistycznej bez redukcji zjawiska multicollinearity

In [24]:
ypred = LR_VIF.predict(Xtest)

from sklearn.metrics import classification_report, confusion_matrix
from sklearn import metrics

co_matrix = metrics.confusion_matrix(ytest, ypred)
co_matrix
Out[24]:
array([[1387,  158],
       [ 182, 1361]], dtype=int64)
In [26]:
print(classification_report(ytest, ypred))

                  precision    recall  f1-score   support

(220.999, 964.0]       0.88      0.90      0.89      1545
 (964.0, 2523.0]       0.90      0.88      0.89      1543

        accuracy                           0.89      3088
       macro avg       0.89      0.89      0.89      3088
    weighted avg       0.89      0.89      0.89      3088

Nowy model zbudowany na pięciu niezwiązanych przez multicollinearity jest nieznacznie gorszy (na poziomie 1-2

Artykuł Przykład klasyfikacji wykonanej za pomocą regresji logistycznej. Eliminacja współliniowości zmiennych niezależnych za pomocą VIF pochodzi z serwisu THE DATA SCIENCE LIBRARY.

]]> Praktyczny przykład wykorzystania ANOVA two-way w Pythonie https://sigmaquality.pl/uncategorized/praktyczny-przyklad-wykorzystania-anova-two-way-w-pythonie-pl080120201039/ Wed, 08 Jan 2020 19:43:00 +0000 http://sigmaquality.pl/praktyczny-przyklad-wykorzystania-anova-two-way-w-pythonie-pl080120201039/ Source of data: http://faculty.webster.edu/woolflm/8canswer.html Model ANOVA two-way (two-factor) W tej analizie robimy porównie dwóch czyników w tych samych populacjach. Porównuje się na przykład temperaturę i [...]

Artykuł Praktyczny przykład wykorzystania ANOVA two-way w Pythonie pochodzi z serwisu THE DATA SCIENCE LIBRARY.

]]>

 PL080120201039 

Source of data: http://faculty.webster.edu/woolflm/8canswer.html

Model ANOVA two-way (two-factor)

W tej analizie robimy porównie dwóch czyników w tych samych populacjach.
Porównuje się na przykład temperaturę i ciśnienie w procesie w kilku różnych okresach.
Tworzy się trzy hipotezy zerowe:

  • h01: średnia temperatury w procesach jest równa
  • h02: średnie ciśnienie w procesach jest równe
  • h03: nie ma powiązania między ciśnieniem a temperaturą

Mówimy, że występuje interakcja pomiędzy dwoma czynnikami, jeżeli efekt uzyskany przy danym poziomie jednego czynnika zależy od poziomu drugiego czynnika. Brak interakcji oznacza, że czynniki są addytywne.

Przeprowadzając two-way ANOVA powinniśmy najpierw sprawdzić trzecią hipotezę zerową.

Analiza ANOVA jest dla H01 i H02 dwiema ANOWA jednoczynnikowymi jednak przy H03 jest już połączeniem tych dwóch ANOVA one-way czyli ANOVA two-way.

ANOVA trzy czynnikowa

Polega na trzech hipotezach zerowych ANOVA one-way o równości średnich dla czynników A, B, C oraz kilku hipotez interakcji: AB, AC, BC, ABC.

Case Study ANOVA

Przeprowadzono badania, których celem było określenie poziomu radości życia. Postawiono następujące hipotezy zerowe:

  • h0a – średnia radość życia jest taka sama dla kobiet i mężczyzn
  • H0b – średnia radość życia jest taka sama wśród trzech grup: 'Young Adult’, 'Middle Adult’, 'Older Adult’.
  • H0c – nie istnieje interakcji pomiędzy grupami: 'Young Adult’, 'Middle Adult’, 'Older Adult’ a płcią badanych osób mającej wpływ na poziom radość życia.
In [1]:
import pandas as pd

cc =['Male', 'Male', 'Male', 'Male','Male','Female','Female','Female','Female','Female']
aaa	 = [4,2,3,4,2,7,4,3,6,5]
bbb = [7,5,7,5,6,8,10,7,7,8]
ccc = [10,7,9,8,11,10,9,12,11,13]
df = pd.DataFrame({'Group': cc, 'Young Adult':aaa, 'Middle Adult':bbb, 'Older Adult':ccc })
df
Out[1]:
Group Young Adult Middle Adult Older Adult
0 Male 4 7 10
1 Male 2 5 7
2 Male 3 7 9
3 Male 4 5 8
4 Male 2 6 11
5 Female 7 8 10
6 Female 4 10 9
7 Female 3 7 12
8 Female 6 7 11
9 Female 5 8 13

Warunki ANOVA

Analiza ANOVA to analiza stwierdzająca istnienie różnic pomiędzy średnimi w kilku populacjch. ANOVA weryfikuje zgodność średnich w kilku populachach anlizując ich wariancję. Celem AVONA jest wykrycie różnic między średnimi jednak nazywa się ją analizą wariancji. Założenia stosowania anailzy wariancji ANOVA:

  • Próby zostały pobrane niezależnie od siebie z każdej z r populacji
  • W każdej z r badanych populacji rozkłąd jest normalny o tej samej wariancji. Średnie w r populacjach mogą być równę, lecz nie muszą*

*źródło: Amir D. Aczel, „Statystyka w zarządzaniu” PWN str.391

Spełnienie warunków 1 i 2 jest niezbędne do zastosowania analizy ANOVA. Dopuszcza się rozkłady zbliżone do normalnych. Jeżel rozkłady są silnie skośne itp. znacznie odbiegają od rozkłądów normalnych lub gdy wariancje w przybliżeniu nie są jednakowe nie należy stosować ANOVA. W takim wypadku należy zastosować test nieparametryczny Kruskala-Wallisa.

Sprawdzenie czy zmienne niezależne mają rozkłąd normalny. Wyszły nieistotne wartości p-value (tzn. p>0,05). Wszystkie zmienne mają rozkład normalny.

In [2]:
from scipy.stats import normaltest

stat, p = normaltest(df['Young Adult'])
print("Young Adult:  ",'p=

stat, p = normaltest(df['Middle Adult'])
print("Middle Adult: ",'p=

stat, p = normaltest(df['Older Adult'])
print("Older Adult:  ",'p=

Young Adult:   p=0.691
Middle Adult:  p=0.564
Older Adult:   p=0.582

C:ProgramDataAnaconda3libsite-packagesscipystatsstats.py:1416: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=10
  "anyway, n=
In [3]:
import matplotlib.pyplot as plt
import seaborn as sns


fig, ax = plt.subplots()
df.plot.kde(ax=ax, legend=True, title='Histogram: A vs. B')
ax.set_ylabel('Probability')
ax.grid(axis='y')

Test wykazał, że grupy mają rozkłady normalne.

In [4]:
df.boxplot(column=['Young Adult', 'Middle Adult', 'Older Adult'], grid=False)
Out[4]:
<matplotlib.axes._subplots.AxesSubplot at 0x22ff1f3ecc0>
In [5]:
df.plot.kde()
Out[5]:
<matplotlib.axes._subplots.AxesSubplot at 0x22ff1f3e518>

Sprawdzamy czy wariancje w tych grupach są podobne.

In [6]:
k = df['Young Adult'].var(ddof=0)
print("'Young Adult:  ",'k=

k = df['Middle Adult'].var(ddof=0)
print("Middle Adult:  ",'k=

k = df['Older Adult'].var(ddof=0)
print("Older Adult:   ",'k=

'Young Adult:   k=2.4
Middle Adult:   k=2.0
Older Adult:    k=3.0

Wariancje mają zbliżoną wartoś, można zastosować test ANOVA.

Przekształcamy dataframe do postaci wymaganej przez statmodels.

In [7]:
df_melt = pd.melt(df.reset_index(),  id_vars=['Group'], value_vars=['Young Adult', 'Middle Adult', 'Older Adult'])
df_melt.sample(5)
Out[7]:
Group variable value
24 Male Older Adult 11
20 Male Older Adult 10
3 Male Young Adult 4
16 Female Middle Adult 10
2 Male Young Adult 3
In [8]:
import statsmodels.api as sm
from statsmodels.formula.api import ols


model = ols('value ~ C(Group) + C(variable) + C(Group):C(variable)', data=df_melt).fit()
anova_table = sm.stats.anova_lm(model, typ=2)
anova_table
Out[8]:
sum_sq df F PR(>F)
C(Group) 3.000000e+01 1.0 1.636364e+01 4.700077e-04
C(variable) 1.800000e+02 2.0 4.909091e+01 3.299559e-09
C(Group):C(variable) 1.183291e-29 2.0 3.227158e-30 1.000000e+00
Residual 4.400000e+01 24.0 NaN NaN
In [9]:
import matplotlib.pyplot as plt
import seaborn as sns

kot = ["#c0c7ce", "#e40c2b"]

sns.set(font_scale=1.5)
sns.factorplot('variable','value',hue='Group',data=df_melt , palette=kot, height=6, aspect=1.2)
plt.title('Agitation by type of music', fontsize=20)

C:ProgramDataAnaconda3libsite-packagesseaborncategorical.py:3666: UserWarning: The `factorplot` function has been renamed to `catplot`. The original name will be removed in a future release. Please update your code. Note that the default `kind` in `factorplot` (`'point'`) has changed `'strip'` in `catplot`.
  warnings.warn(msg)
Out[9]:
Text(0.5, 1, 'Agitation by type of music')

Odrzucamy hipotezę H0a: Istnieją znaczące różnice pomiędzy płcią (F = 16,36, p=0,0004700077, p <0,01)

Odrzucamy hipotezę H0b: Istnieją znaczące różnice pomiędzy grupami wiekowymi (F = 49,09, p=0,000000003299559, p <0,01)

Zachowujemy hipotezę H0c: Brak efektu interakcji (F = 0,00, p=1 nieistotny).

Z danych wynika, że starsi dorośli mają najwyższy poziom zadowolenia z życia, a młodsi dorośli mają najniższy poziom zadowolenia z życia. Kobiety mają także znacznie wyższą satysfakcję z życia niż mężczyźni.

Dzielimy grupę na mężczyzn i kobiety i oddzielnie testujemy Tukey HSD

In [10]:
Male = df_melt[df_melt['Group']=='Male']
Female = df_melt[df_melt['Group']=='Female']

Tukey HSD Test dla mężczyzn

In [11]:
import scipy.stats as stats
from statsmodels.stats.multicomp import (pairwise_tukeyhsd,MultiComparison)

# Tworzę specjalną wartość
KOT = MultiComparison(Male['value'],Male['variable'])
AA = KOT.tukeyhsd()
print(AA.summary())

   Multiple Comparison of Means - Tukey HSD,FWER=0.05   
========================================================
   group1       group2   meandiff  lower   upper  reject
--------------------------------------------------------
Middle Adult Older Adult   3.0     0.9345  5.0655  True 
Middle Adult Young Adult   -3.0   -5.0655 -0.9345  True 
Older Adult  Young Adult   -6.0   -8.0655 -3.9345  True 
--------------------------------------------------------
In [12]:
AA.plot_simultaneous()    
plt.vlines(x=160,ymin=-0.5,ymax=3.5, color="red")
Out[12]:
<matplotlib.collections.LineCollection at 0x22ff2f43278>

Istnieją duże różnice w zadowoleniu z życia pomiędzy grupami mężczyzn: 'Young Adult’, 'Middle Adult’, 'Older Adult’.

Test Bonferroniego

Jest to test substytucyjny do testu tukey HDS.

In [13]:
import scipy.stats as stats
from statsmodels.stats.multicomp import (pairwise_tukeyhsd,MultiComparison)

# Tworzę specjalną wartość
KOT = MultiComparison(Male['value'],Male['variable'])
comp = KOT.allpairtest(stats.ttest_rel, method='Holm')
print (comp[0])

Test Multiple Comparison ttest_rel 
FWER=0.05 method=Holm
alphacSidak=0.02, alphacBonf=0.017
========================================================
   group1       group2     stat   pval  pval_corr reject
--------------------------------------------------------
Middle Adult Older Adult -5.4772 0.0054   0.0108   True 
Middle Adult Young Adult  5.4772 0.0054   0.0108   True 
Older Adult  Young Adult  7.1714 0.002    0.006    True 
--------------------------------------------------------

Tukey HSD Test dla kobiet

In [14]:
import scipy.stats as stats
from statsmodels.stats.multicomp import (pairwise_tukeyhsd,MultiComparison)

# Tworzę specjalną wartość
KOT = MultiComparison(Female['value'],Female['variable'])
FF = KOT.tukeyhsd()
print(AA.summary())

   Multiple Comparison of Means - Tukey HSD,FWER=0.05   
========================================================
   group1       group2   meandiff  lower   upper  reject
--------------------------------------------------------
Middle Adult Older Adult   3.0     0.9345  5.0655  True 
Middle Adult Young Adult   -3.0   -5.0655 -0.9345  True 
Older Adult  Young Adult   -6.0   -8.0655 -3.9345  True 
--------------------------------------------------------
In [15]:
FF.plot_simultaneous()    
plt.vlines(x=160,ymin=-0.5,ymax=3.5, color="red")
Out[15]:
<matplotlib.collections.LineCollection at 0x22ff2f9b438>

Zarówno wśród kobiet jak i wśród mężczyzn istnieje wysoki poziom różnic pomiędzy grupami wiekowymi.

 PL080120201039 

Artykuł Praktyczny przykład wykorzystania ANOVA two-way w Pythonie pochodzi z serwisu THE DATA SCIENCE LIBRARY.

]]> Problem występowania współliniowości zmiennych niezależnych w regresji logistycznej. Przykład użycia testu VIF. https://sigmaquality.pl/uncategorized/multicollinearity-w-regresji-logistycznej-vif-pl202001081905/ Wed, 08 Jan 2020 18:08:51 +0000 http://sigmaquality.pl/multicollinearity-w-regresji-logistycznej-vif-pl202001081905/ PL202001081905 Jednym z najważniejszych warunków w modelu regresji logistycznej jest brak współliniowości (multicollinearity) zmiennych niezależnych. Oznacza to, że model powinien mieć niewielką współliniowość zmiennych opisujących [...]

Artykuł Problem występowania współliniowości zmiennych niezależnych w regresji logistycznej. Przykład użycia testu VIF. pochodzi z serwisu THE DATA SCIENCE LIBRARY.

]]>
PL202001081905

Jednym z najważniejszych warunków w modelu regresji logistycznej jest brak współliniowości (multicollinearity) zmiennych niezależnych. Oznacza to, że model powinien mieć niewielką współliniowość zmiennych opisujących lub jej brak.

Multicollinearity zmiennych niezależnych prowadzi do niewiarygodnych i niestabilnych oszacowań. Model nie jest stabilny, niewielkie zmiany w zmiennych niezależnych mogą diametralnie zmieniać wyniki klasyfikacji modelu.

Mimo że regresja logistyczna nie jest regresją lecz klasyfikacją, obowiązują tu podobne zasady. Ważna jest zmiana zmiennej zależnej dla każdej 1 zmiany jednostki w zmiennej niezależnej, gdy wszystkie pozostałe zmienne niezależne są stałe.

Chodzi o to, żeby zmienna jednej zmiennej nie wiązała się ze zmianą innych zmiennych niezależnych. Gdy zmienne niezależne są skorelowane, oznacza to, że zmiany w jednej zmiennej są powiązane ze zmianami w innej zmiennej. Im silniejsza korelacja, tym trudniej jest zmienić jedną zmienną bez zmiany drugiej. Modelowi trudno jest oszacować związek między każdą zmienną niezależną a zmienną zależną niezależnie, ponieważ wtedy gdy zmienne niezależne są skorelowane zmieniają się zgodnie.

Dowodem na istnienie współliniowości (multicollinearity) zmiennych niezależnych w modelu regresji logistycznej jest ich skorelowanie. Ta korelacja stanowi problem, ponieważ zmienne niezależne powinny być niezależne. Jeśli stopień korelacji między zmiennymi jest wystarczająco wysoki, może to powodować problemy w dopasowaniu modelu i interpretacji wyników.

Jak radzić sobie z problemem multicollinearity w modelu regresji logistycznej?

Multicolinearity można sprawdzić robiąc prostą macierz regresji dla zmiennych niezależnych.
Korfelacja zmiennych niezależnych to nie to samo co multicollinearity ponieważ multicollinearity może się pojawić w modelu, nawet gdy izolowane pary zmiennych nie są współliniowe (skorelowane), pojawia się interakcja 3 zmiennych itd.
Dlatego multicollinearity jest czasem bardziej trudne do wykrycia.

Variance Inflation Factor (VIF)

Variance Inflation Factor (VIF) to miara multicollinearity między zmiennymi predykcyjnymi w regresji wielokrotnej. Określa ilościowo nasilenie multicollinearity w zwykłej analizie regresji metodą najmniejszych kwadratów . Zapewnia wskaźnik, który mierzy, o ile wariancja (kwadrat odchylenia standardowego oszacowania ) szacowanego współczynnika regresji jest zwiększona z powodu kolinearności.

Pierwiastek kwadratowy współczynnika inflacji wariancji wskazuje, o ile większy błąd standardowy wzrasta w porównaniu z tym, czy zmienna ta ma korelację 0 z innymi zmiennymi predykcyjnymi w modelu.

*Przykład
Jeśli współczynnik inflacji wariancji zmiennej predykcyjnej wynosił 5,27 (√5,27 = 2,3), oznacza to, że błąd standardowy współczynnika tej zmiennej predykcyjnej jest 2,3 razy większy niż w przypadku, gdy zmienna predyktorowa ma 0 korelacji z innymi zmiennymi predykcyjnymi.*

(źródło:https://en.wikipedia.org/wiki/Variance_inflation_factor)

Kroki wdrażania VIF

  1. Uruchom regresję wieloraką.
  2. Oblicz współczynniki VIF.
  3. Sprawdź współczynniki dla każdej zmiennej predykcyjnej, jeśli VIF wynosi między 5-10, prawdopodobnie występuje multicollinearity i powinieneś rozważyć usunięcie tej zmiennej.

źródło: https://etav.github.io/python/vif_factor_python.html

Czyli przed rozpoczęciem pracy nad modelem regresji logistycznej należy utworzyć model regresji wielorakiej i przez VIF wyselekcjonować, wybrać dane do modelu regresji logistycznej.

Jak wspominają tutaj: https://www.ibm.com/support/pages/multicollinearity-diagnostics-logistic-regression-nomreg-or-plum

Procedury regresji dla zmiennych zależnych kategorycznie nie mają diagnostyki kolinearności. W tym celu można jednak użyć procedury regresji liniowej. Statystyka kolinearności w regresji dotyczy relacji między predyktorami, ignorując zmienną zależną. Możesz więc uruchomić REGRESSION z tą samą listą predyktorów i zmiennej zależnej, jakiej chcesz użyć w REGRESJI LOGISTYCZNEJ (na przykład) i zażądać diagnostyki kolinearności. Uruchom regresję logistyczną, aby uzyskać właściwe współczynniki, przewidywane prawdopodobieństwa itp. Po podjęciu niezbędnych decyzji (porzucenie predyktorów itp.) Wynikających z analizy kolinearności.

In [1]:
import pandas as pd
import numpy as np

df = pd.read_csv('c:/2/poliaxid.csv', index_col=0)
#del df['Unnamed: 0.1']
df.head(5)
Out[1]:
nr. factorA factorB citric catoda residual butanol caroton stable nodinol sulfur in nodinol density pH noracid lacapon quality class
0 0 4.933333 0.466667 0.000000 1.266667 0.050667 7.333333 22.666667 0.665200 2.340000 0.373333 6.266667 1
1 1 5.200000 0.586667 0.000000 1.733333 0.065333 16.666667 44.666667 0.664533 2.133333 0.453333 6.533333 1
2 2 5.200000 0.506667 0.026667 1.533333 0.061333 10.000000 36.000000 0.664667 2.173333 0.433333 6.533333 1
3 3 7.466667 0.186667 0.373333 1.266667 0.050000 11.333333 40.000000 0.665333 2.106667 0.386667 6.533333 2
4 4 4.933333 0.466667 0.000000 1.266667 0.050667 7.333333 22.666667 0.665200 2.340000 0.373333 6.266667 1

Z ciekawości sprawdzamy jaki jest poziom korelacji pomiędzy zmiennymi niezależnymi.

In [2]:
import seaborn as sns
import matplotlib.pyplot as plt
plt.figure(figsize=(10,10))


CORREL = df.corr()
sns.heatmap(CORREL, cmap="YlGnBu", annot=True, cbar=False)
plt.show
Out[2]:
<function matplotlib.pyplot.show(*args, **kw)>

Widać wyraźnie, że niektóre zmienne są ze sobą silnie skorelowane. Nie możemy ich ślepo eliminować z modelu. Następnym krokiem jest stworzenie prostego modelu regresji wielorakiej OLS.

Rutynowo sprawdzamy kompletność i format danych.

In [3]:
df.isnull().sum()
Out[3]:
nr.                  0
factorA              0
factorB              0
citric catoda        0
residual butanol     0
caroton              0
stable nodinol       0
sulfur in nodinol    0
density              0
pH                   0
noracid              0
lacapon              0
quality class        0
dtype: int64
In [4]:
df.dtypes
Out[4]:
nr.                    int64
factorA              float64
factorB              float64
citric catoda        float64
residual butanol     float64
caroton              float64
stable nodinol       float64
sulfur in nodinol    float64
density              float64
pH                   float64
noracid              float64
lacapon              float64
quality class          int64
dtype: object

Budujemy regresję liniową w bibliotece statmodel.

Tworzymy model regresji liniowej

In [5]:
df.columns
Out[5]:
Index(['nr.', 'factorA', 'factorB', 'citric catoda', 'residual butanol',
       'caroton', 'stable nodinol', 'sulfur in nodinol', 'density', 'pH',
       'noracid', 'lacapon', 'quality class'],
      dtype='object')
In [6]:
df.columns=['nr.', 'factorA', 'factorB', 'citric_catoda', 'residual_butanol',
       'caroton', 'stable_nodinol', 'sulfur_in_nodinol', 'density', 'pH',
       'noracid', 'lacapon', 'quality_class']
In [7]:
#import statsmodels.api as sm
#from statsmodels.formula.api import ols

VIF (Variation Inflation Factor)

Dzięki temu czynnikowi jesteśmy w stanie wskazać, która zmienna powinna zostać wyeliminowana z modelu.

In [8]:
df.columns
Out[8]:
Index(['nr.', 'factorA', 'factorB', 'citric_catoda', 'residual_butanol',
       'caroton', 'stable_nodinol', 'sulfur_in_nodinol', 'density', 'pH',
       'noracid', 'lacapon', 'quality_class'],
      dtype='object')
In [9]:
from patsy import dmatrices
from statsmodels.stats.outliers_influence import variance_inflation_factor
import statsmodels.formula.api as smf

lm = smf.ols(formula = 'quality_class ~ factorA+factorB+citric_catoda+residual_butanol+caroton+stable_nodinol+density+sulfur_in_nodinol+pH+noracid+lacapon', data = df).fit()
y, X = dmatrices('quality_class ~ factorA+factorB+citric_catoda+residual_butanol+caroton+stable_nodinol+density+sulfur_in_nodinol+pH+noracid+lacapon', data = df, return_type = "dataframe")
vif = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
print(vif)

[1709157.8270185539, 7.777537721773227, 1.7885509640263009, 3.1329923025205733, 1.7030403470293642, 1.4834333245817521, 1.9734391365941484, 6.3474498618969335, 2.20192092692558, 3.323580122673945, 1.429921652882204, 3.034070351762669]

Ustawiamy macierz wyników.

In [10]:
X.columns
Out[10]:
Index(['Intercept', 'factorA', 'factorB', 'citric_catoda', 'residual_butanol',
       'caroton', 'stable_nodinol', 'density', 'sulfur_in_nodinol', 'pH',
       'noracid', 'lacapon'],
      dtype='object')
In [11]:
vif  =  np.round(vif, decimals=2)
In [12]:
vif = list(map(float, vif))
name = list(X)
In [13]:
s1=pd.Series(name,name='name')
s2=pd.Series( vif,name='vif')

RFE_list = pd.concat([s1,s2], axis=1)

RFE_list
Out[13]:
name vif
0 Intercept 1709157.83
1 factorA 7.78
2 factorB 1.79
3 citric_catoda 3.13
4 residual_butanol 1.70
5 caroton 1.48
6 stable_nodinol 1.97
7 density 6.35
8 sulfur_in_nodinol 2.20
9 pH 3.32
10 noracid 1.43
11 lacapon 3.03

Interpretacja: wynik w postaci wektora reprezentuje zmienną w określonej kolejności jak w modelu. W zaleceniu VIF wskazano, że jeśli współczynnik przypisany do zmiennej jest większy niż 5, zmienna ta jest wysoce skorelowana z innymi zmiennymi i powinna zostać wyeliminowana z modelu.

Test wykazał, że zmienne: density oraz factorA powinny zostać usunięta z modelu.
jeszcze raz oglądamy macierz korelacji

Jeszcze raz tworzymy macierz korelacji vif i sprawdzamy zmienne po wyeliminowaniu: factorA i density.

In [14]:
lm = smf.ols(formula = 'quality_class ~ factorB+citric_catoda+residual_butanol+caroton+stable_nodinol+sulfur_in_nodinol+pH+noracid+lacapon', data = df).fit()
y, X = dmatrices('quality_class ~ factorB+citric_catoda+residual_butanol+caroton+stable_nodinol+sulfur_in_nodinol+pH+noracid+lacapon', data = df, return_type = "dataframe")
vif2 = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
print(vif2)

[809.2230489081707, 1.6530213053949512, 2.1806661580635818, 1.093561268844909, 1.3702675584858117, 1.9475860748759193, 2.0181057450324356, 1.6077595870352759, 1.3354874384986997, 1.2793132834168612]
In [15]:
X2 = 'Intercept', 'factorB', 'citric_catoda', 'residual_butanol','caroton', 'stable_nodinol', 'sulfur_in_nodinol', 'pH','noracid', 'lacapon'


vif  =  np.round(vif2, decimals=2) 
vif = list(map(float, vif2))
name = list(X2)

s1=pd.Series(name,name='name')
s2=pd.Series( vif2,name='vif')

RFE_list = pd.concat([s1,s2], axis=1)

RFE_list
Out[15]:
name vif
0 Intercept 809.223049
1 factorB 1.653021
2 citric_catoda 2.180666
3 residual_butanol 1.093561
4 caroton 1.370268
5 stable_nodinol 1.947586
6 sulfur_in_nodinol 2.018106
7 pH 1.607760
8 noracid 1.335487
9 lacapon 1.279313

Teraz wyraźnie zmniejszył się poziom vif, oznacza to że zmniejszył się stopień multicollinearity.

ciąg dalszy tutaj:

Przykład klasyfikacji wykonanej za pomocą regresji logistycznej. Eliminacja współliniowości zmiennych niezależnych za pomocą VIF

Artykuł Problem występowania współliniowości zmiennych niezależnych w regresji logistycznej. Przykład użycia testu VIF. pochodzi z serwisu THE DATA SCIENCE LIBRARY.

]]> Graficzny przykład ANOVA two-way – interaction_plot https://sigmaquality.pl/anova/graficzny-przyklad-anova-two-way-interaction_plot-pl202001081306/ Wed, 08 Jan 2020 18:08:00 +0000 http://sigmaquality.pl/graficzny-przyklad-anova-two-way-interaction_plot-pl202001081306/ PL202001081306 Source of data: http://faculty.webster.edu/woolflm/8canswer.html Model ANOVA two-way (two-factor) W tej analizie robimy porównie dwóch czyników w tych samych populacjach. Analiza ANOVA jest dla H01 [...]

Artykuł Graficzny przykład ANOVA two-way – interaction_plot pochodzi z serwisu THE DATA SCIENCE LIBRARY.

]]> PL202001081306

Source of data: http://faculty.webster.edu/woolflm/8canswer.html

Model ANOVA two-way (two-factor)

W tej analizie robimy porównie dwóch czyników w tych samych populacjach.

Analiza ANOVA jest dla H01 i H02 dwiema ANOVA jednoczynnikowymi jednak przy H03 jest już połączeniem tych dwóch ANOVA one-way czyli ANOVA two-way.

Przeprowadzono badania, których celem było określenie poziomu radości życia. Postawiono następujące hipotezy zerowe:

  • h0a – średnia radość życia jest taka sama dla kobiet i mężczyzn
  • H0b – średnia radość życia jest taka sama wśród trzech grup: 'Young Adult’, 'Middle Adult’, 'Older Adult’.
  • H0c – nie istnieje interakcji pomiędzy grupami: 'Young Adult’, 'Middle Adult’, 'Older Adult’ a płcią badanych osób mającej wpływ na poziom radość życia.

Kodujemy bazwy grup:

  • ’Young Adult’ :1
  • ’Middle Adult’:2
  • ’Older Adult’ :3
In [5]:
import pandas as pd

cc =['Male', 'Male', 'Male', 'Male','Male','Female','Female','Female','Female','Female']
aaa	 = [4,2,3,4,2,7,4,3,6,5]
bbb = [7,5,7,5,6,8,10,7,7,8]
ccc = [10,7,9,8,11,10,9,12,11,13]
df = pd.DataFrame({'Group': cc, 1:aaa, 2:bbb, 3:ccc })
df
Out[5]:
Group 1 2 3
0 Male 4 7 10
1 Male 2 5 7
2 Male 3 7 9
3 Male 4 5 8
4 Male 2 6 11
5 Female 7 8 10
6 Female 4 10 9
7 Female 3 7 12
8 Female 6 7 11
9 Female 5 8 13

Badaliśy już ten zbiór danych pod kontem możliwości wykorzystania go w analizie ANOVA.

  • Próby zostały pobrane niezależnie od siebie z każdej z r populacji
  • W każdej z r badanych populacji rozkłąd jest normalny o tej samej wariancji. Średnie w r populacjach mogą być równę, lecz nie muszą*

Analiza ANOVA tej próby jest opisana w linku poniżej:
http://sigmaquality.pl/python/praktyczny-przyklad-wykorzystania-anova-two-way-w-pythonie-pl080120201039/

Przekształacamy dataframe do postaci zgodniej ze statmodels.

In [12]:
df_melt = pd.melt(df.reset_index(),  id_vars=['Group'], value_vars=[1, 2, 3])
df_melt.columns = ['Sex','Age','value']
df_melt.sample(5)
Out[12]:
Sex Age value
9 Female 1 5
3 Male 1 4
0 Male 1 4
20 Male 3 10
28 Female 3 11
In [34]:
df_melt['Age'] = df_melt['Age'].astype('int64')
df_melt.dtypes
Out[34]:
Sex      object
Age       int64
value     int64
dtype: object

Teraz tworzymy dwa modele regresji wielorakiej OLS.

In [49]:
from statsmodels.formula.api import ols
from statsmodels.stats.anova import anova_lm

Reg1 = ols(formula = "value ~ Age + Sex", data = df_melt)
Fit1 = Reg1.fit()

Reg2 = ols(formula = "value ~ Age*Sex", data = df_melt)
Fit2 = Reg2.fit()
In [50]:
print(df_melt['value'].dtype)
print(df_melt['Age'].dtype)
print(df_melt['Sex'].dtype)

int64
int64
object
In [63]:
from statsmodels.graphics.factorplots import interaction_plot


fig = interaction_plot(df_melt['Age'], df_melt['Sex'],df_melt['value'],
             colors=['black','darkgreen'], markers=['D','x'], linestyles=["-", "--"], ylabel='radość życia', xlabel='Grupa wiekowa')
fig = interaction_plot(df_melt['Age'],df_melt['Sex'], Fit1.fittedvalues,
            colors=['black','darkgreen'], markers=['D','x'],linestyles=["-", "--"], ylabel='radość życia - OLS', xlabel='Grupa wiekowa')
fig = interaction_plot(df_melt['Age'],df_melt['Sex'], Fit2.fittedvalues,
             colors=['black','darkgreen'], markers=['D','x'], linestyles=["-", "--"], ylabel='radość życia - OLS interakcja', xlabel='Grupa wiekowa')

import matplotlib.pyplot as plt
plt.show()

Interpretacja wykresu interakcji

Wykres interakcji pokazuje w jaki sposób związek między jednym czynnikiem kategorycznym (płeć) a ciągłą reakcją zależy od wartości drugiego czynnika kategorycznego (grupa wiekowa). Ten wykres wyświetla średnie dla poziomów jednego czynnika na osi y oraz osobną linię dla każdego czynnika.

  • Równoległe linie oznaczają brak interakcji.
  • Linie nierównoległe – występuje interakcja.

Im bardziej nierównoległe są linie, tym większa jest siła interakcji.
Chociaż można użyć tego wykresu do wyświetlenia efektów, należy wykonać odpowiedni test ANOVA i ocenić istotność statystyczną efektów. Jeśli efekty interakcji są znaczące, nie można interpretować głównych efektów bez uwzględnienia efektów interakcji.

 PL080120201039 

Artykuł Graficzny przykład ANOVA two-way – interaction_plot pochodzi z serwisu THE DATA SCIENCE LIBRARY.

]]>