What is the probability of infraction given a certain risk score level?

In the article "Risk management systems: using data mining in developing countries’ customs administrations", Bertrand Laporte developed a risk management system that uses statistical scoring techniques [1]. At best, the score should reflect risk of infraction (or even the probability of an infraction occurring). The article also compiled three table related with three ways of computing score: simple average, weighted average, and maximum value. For these three scores it was given a proportion (or can be seen as a probability) of given a score interval what are the cumulative number of declarations with infraction.

These three score values from article is compile and summarized in an Excel file:

https://docs.google.com/spreadsheets/d/1lVEJjvqxdOaDSTpjWu8812vb__ymIG4z/edit?usp=sharing&ouid=106640872116257813737&rtpof=true&sd=true

For a more detailed description about the  context of the problem addressed by the article [1], article [2] could be useful.

From the data, one interesting question that emerges is: "What is the probability of infraction given a certain risk score level". This probability could be estimated employing linear regression as will be done by the next Python commands.

Since the three scores presented a linear correlation with the variable category, with a notable exception in the first category, then a linear regression to predict the probability of a declaration has an infraction, given that it belongs to a certain score interval, could be built.

import pandas as pd

from sklearn import linear_model

def linear_regression_report(x, y):

  regr_ = linear_model.LinearRegression()

  regr_.fit(x, y)

  return regr_

def model_report(model, name='None '):

  print('-------------------------------')

  print(name,' model')

  print('-------------------------------')

  print('Intercept: \n', model.intercept_)

  print('Coefficients: \n', model.coef_)

# Defining independent and dependent variables

x = df[['Category']]

y_sa = df['Simple Average (1)']

y_wa = df['Weighted Average (2)']

y_mv = df['Maximum Value (3)']

# Linear regression with sklearn to predict score using simple average

model_sa = linear_regression_report(x, y_sa)

# Report about the linear regression model using score as an output

model_report(model_sa, 'Score using Simple Average')

# Linear regression with sklearn to predict score using weighted average

model_wa = linear_regression_report(x, y_wa)

# Report about the linear regression model using importation as an output

model_report(model_wa, 'Score using Weighted Average')

# Linear regression with sklearn to predict score using maximum value

model_mv = linear_regression_report(x, y_mv)

# Report about the linear regression model using score as an output

model_report(model_mv, 'Score using Maximum Value')

-------------------------------

Score using Simple Average  model

-------------------------------

Intercept: 

 3.886666666666656

Coefficients: 

 [10.21151515]

-------------------------------

Score using Weighted Average  model

-------------------------------

Intercept: 

 16.146666666666647

Coefficients: 

 [9.10060606]

-------------------------------

Score using Maximum Value  model

-------------------------------

Intercept: 

 33.873333333333335

Coefficients: 

 [7.67030303] 

A more detailed linear regression model could be obtained with the following commands.

import statsmodels.api as sm

def get_model(x, y):

  # Obtaining multiple linear regression using statsmodels: exportation is the output

  xm = sm.add_constant(x) # adding a constant

  model = sm.OLS(y, xm).fit()

  return model, xm

def print_model(model, name='None'):

  summary = model.summary()

  print('--------------------------------------')

  print(name,'model')

  print(summary)

  print('--------------------------------------')

# Obtaining multiple linear regression using statsmodels: sa is the input

model_sa, xm_sa = get_model(x, y_sa)

print_model(model_sa,'Simple Average')

# Obtaining multiple linear regression using statsmodels: wa is the input

model_wa, xm_wa = get_model(x, y_wa)

print_model(model_wa,'Weighted Average')

# Obtaining multiple linear regression using statsmodels: mv is the input

model_mv, xm_mv = get_model(x, y_mv)

print_model(model_mv,'Maximum value')

--------------------------------------

Simple Average model

                            OLS Regression Results                            

==============================================================================

Dep. Variable:     Simple Average (1)   R-squared:                       0.959

Model:                            OLS   Adj. R-squared:                  0.954

Method:                 Least Squares   F-statistic:                     186.8

Date:                Mon, 06 Nov 2023   Prob (F-statistic):           7.91e-07

Time:                        20:29:56   Log-Likelihood:                -32.222

No. Observations:                  10   AIC:                             68.44

Df Residuals:                       8   BIC:                             69.05

Df Model:                           1                                         

Covariance Type:            nonrobust                                         

==============================================================================

                 coef    std err          t      P>|t|      [0.025      0.975]

------------------------------------------------------------------------------

const          3.8867      4.636      0.838      0.426      -6.804      14.577

Category      10.2115      0.747     13.668      0.000       8.489      11.934

==============================================================================

Omnibus:                        4.516   Durbin-Watson:                   1.548

Prob(Omnibus):                  0.105   Jarque-Bera (JB):                1.623

Skew:                          -0.958   Prob(JB):                        0.444

Kurtosis:                       3.471   Cond. No.                         13.7

==============================================================================

Notes:

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

--------------------------------------

--------------------------------------

Weighted Average model

                             OLS Regression Results                             

================================================================================

Dep. Variable:     Weighted Average (2)   R-squared:                       0.937

Model:                              OLS   Adj. R-squared:                  0.930

Method:                   Least Squares   F-statistic:                     119.8

Date:                  Mon, 06 Nov 2023   Prob (F-statistic):           4.31e-06

Time:                          20:29:56   Log-Likelihood:                -33.291

No. Observations:                    10   AIC:                             70.58

Df Residuals:                         8   BIC:                             71.19

Df Model:                             1                                         

Covariance Type:              nonrobust                                         

==============================================================================

                 coef    std err          t      P>|t|      [0.025      0.975]

------------------------------------------------------------------------------

const         16.1467      5.159      3.130      0.014       4.251      28.043

Category       9.1006      0.831     10.946      0.000       7.183      11.018

==============================================================================

Omnibus:                        8.987   Durbin-Watson:                   1.240

Prob(Omnibus):                  0.011   Jarque-Bera (JB):                3.858

Skew:                          -1.413   Prob(JB):                        0.145

Kurtosis:                       4.126   Cond. No.                         13.7

==============================================================================

Notes:

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

--------------------------------------

--------------------------------------

Maximum value model

                            OLS Regression Results                            

==============================================================================

Dep. Variable:      Maximum Value (3)   R-squared:                       0.771

Model:                            OLS   Adj. R-squared:                  0.742

Method:                 Least Squares   F-statistic:                     26.92

Date:                Mon, 06 Nov 2023   Prob (F-statistic):           0.000834

Time:                        20:29:56   Log-Likelihood:                -39.046

No. Observations:                  10   AIC:                             82.09

Df Residuals:                       8   BIC:                             82.70

Df Model:                           1                                         

Covariance Type:            nonrobust                                         

==============================================================================

                 coef    std err          t      P>|t|      [0.025      0.975]

------------------------------------------------------------------------------

const         33.8733      9.173      3.693      0.006      12.721      55.025

Category       7.6703      1.478      5.189      0.001       4.261      11.079

==============================================================================

Omnibus:                        9.301   Durbin-Watson:                   1.114

Prob(Omnibus):                  0.010   Jarque-Bera (JB):                3.984

Skew:                          -1.416   Prob(JB):                        0.136

Kurtosis:                       4.242   Cond. No.                         13.7

==============================================================================

It is also possible to obtain a graphical comparison between the values from the dataset and from the model's predictions.

import matplotlib.pyplot as plt

def model_predict(model, xm):

   y_hat = model.predict(xm)

   return y_hat

def draw_model(model, y, yname, y_hat):

  xp = list(range(1,len(x)+1))

  plt.plot(xp,y,'ob',label='Data')

  plt.plot(xp,y,'-r')

  plt.plot(xp,y_hat,'--g',label='Prediction')

  plt.title('Prediction on ' + str(yname))

  plt.xlabel('Time')

  plt.ylabel(yname)

  plt.legend()

  plt.grid()

  plt.show()

# Predict the values used to build the sa model coefficients.

y_hat_sa = model_predict(model_sa, xm_sa)

#print(y_hat_export)

# Drawing to compare data x prediction in score data.

draw_model(model_sa, y_sa,  'SA', y_hat_sa)

# Predict the values used to build the wa model coefficients.

y_hat_wa = model_predict(model_wa, xm_wa)

#print(y_hat_export)

# Drawing to compare data x prediction in score data.

draw_model(model_wa, y_wa,  'WA', y_hat_wa)

# Predict the values used to build the mv model coefficients.

y_hat_mv = model_predict(model_mv, xm_mv)

#print(y_hat_export)

# Drawing to compare data x prediction in score data.

draw_model(model_mv, y_mv,  'MV', y_hat_mv)