regression-model-diagnosis/diagnosis.py
def diagnostic_plots(X, y, model_fit=None):
  """
  Function to reproduce the 4 base plots of an OLS model in R.

  ---
  Inputs:

  X: A numpy array or pandas dataframe of the features to use in building the linear regression model

  y: A numpy array or pandas series/dataframe of the target variable of the linear regression model

  model_fit [optional]: a statsmodel.api.OLS model after regressing y on X. If not provided, will be
                        generated from X, y
  """

  if not model_fit:
      model_fit = sm.OLS(y, sm.add_constant(X)).fit()

  # create dataframe from X, y for easier plot handling
  dataframe = pd.concat([X, y], axis=1)

  # model values
  model_fitted_y = model_fit.fittedvalues
  # model residuals
  model_residuals = model_fit.resid
  # normalized residuals
  model_norm_residuals = model_fit.get_influence().resid_studentized_internal
  # absolute squared normalized residuals
  model_norm_residuals_abs_sqrt = np.sqrt(np.abs(model_norm_residuals))
  # absolute residuals
  model_abs_resid = np.abs(model_residuals)
  # leverage, from statsmodels internals
  model_leverage = model_fit.get_influence().hat_matrix_diag
  # cook's distance, from statsmodels internals
  model_cooks = model_fit.get_influence().cooks_distance[0]

  plot_lm_1 = plt.figure()
  plot_lm_1.axes[0] = sns.residplot(model_fitted_y, dataframe.columns[-1], data=dataframe,
                            lowess=True,
                            scatter_kws={'alpha': 0.5},
                            line_kws={'color': 'red', 'lw': 1, 'alpha': 0.8})

  plot_lm_1.axes[0].set_title('Residuals vs Fitted')
  plot_lm_1.axes[0].set_xlabel('Fitted values')
  plot_lm_1.axes[0].set_ylabel('Residuals');

  # annotations
  abs_resid = model_abs_resid.sort_values(ascending=False)
  abs_resid_top_3 = abs_resid[:3]
  for i in abs_resid_top_3.index:
      plot_lm_1.axes[0].annotate(i,
                                 xy=(model_fitted_y[i],
                                     model_residuals[i]));

  QQ = ProbPlot(model_norm_residuals)
  plot_lm_2 = QQ.qqplot(line='45', alpha=0.5, color='#4C72B0', lw=1)
  plot_lm_2.axes[0].set_title('Normal Q-Q')
  plot_lm_2.axes[0].set_xlabel('Theoretical Quantiles')
  plot_lm_2.axes[0].set_ylabel('Standardized Residuals');
  # annotations
  abs_norm_resid = np.flip(np.argsort(np.abs(model_norm_residuals)), 0)
  abs_norm_resid_top_3 = abs_norm_resid[:3]
  for r, i in enumerate(abs_norm_resid_top_3):
      plot_lm_2.axes[0].annotate(i,
                                 xy=(np.flip(QQ.theoretical_quantiles, 0)[r],
                                     model_norm_residuals[i]));

  plot_lm_3 = plt.figure()
  plt.scatter(model_fitted_y, model_norm_residuals_abs_sqrt, alpha=0.5);
  sns.regplot(model_fitted_y, model_norm_residuals_abs_sqrt,
              scatter=False,
              ci=False,
              lowess=True,
              line_kws={'color': 'red', 'lw': 1, 'alpha': 0.8});
  plot_lm_3.axes[0].set_title('Scale-Location')
  plot_lm_3.axes[0].set_xlabel('Fitted values')
  plot_lm_3.axes[0].set_ylabel('$\sqrt{|Standardized Residuals|}$');

  # annotations
  abs_sq_norm_resid = np.flip(np.argsort(model_norm_residuals_abs_sqrt), 0)
  abs_sq_norm_resid_top_3 = abs_sq_norm_resid[:3]
  for i in abs_norm_resid_top_3:
      plot_lm_3.axes[0].annotate(i,
                                 xy=(model_fitted_y[i],
                                     model_norm_residuals_abs_sqrt[i]));


  plot_lm_4 = plt.figure();
  plt.scatter(model_leverage, model_norm_residuals, alpha=0.5);
  sns.regplot(model_leverage, model_norm_residuals,
              scatter=False,
              ci=False,
              lowess=True,
              line_kws={'color': 'red', 'lw': 1, 'alpha': 0.8});
  plot_lm_4.axes[0].set_xlim(0, max(model_leverage)+0.01)
  plot_lm_4.axes[0].set_ylim(-3, 5)
  plot_lm_4.axes[0].set_title('Residuals vs Leverage')
  plot_lm_4.axes[0].set_xlabel('Leverage')
  plot_lm_4.axes[0].set_ylabel('Standardized Residuals');

  # annotations
  leverage_top_3 = np.flip(np.argsort(model_cooks), 0)[:3]
  for i in leverage_top_3:
      plot_lm_4.axes[0].annotate(i,
                                 xy=(model_leverage[i],
                                     model_norm_residuals[i]));

  p = len(model_fit.params) # number of model parameters
  graph(lambda x: np.sqrt((0.5 * p * (1 - x)) / x),
        np.linspace(0.001, max(model_leverage), 50),
        'Cook\'s distance') # 0.5 line
  graph(lambda x: np.sqrt((1 * p * (1 - x)) / x),
        np.linspace(0.001, max(model_leverage), 50)) # 1 line
  plot_lm_4.legend(loc='upper right');

  Function to reproduce the 4 base plots of an OLS model in R.

  ---
  Inputs:

  X: A numpy array or pandas dataframe of the features to use in building the linear regression model

  y: A numpy array or pandas series/dataframe of the target variable of the linear regression model

  model_fit [optional]: a statsmodel.api.OLS model after regressing y on X. If not provided, will be
                        generated from X, y
  
Function to reproduce the 4 base plots of an OLS model in R.
Inputs.

A numpy array or pandas dataframe of the features to use in building the linear regression model

A numpy array or pandas series/dataframe of the target variable of the linear regression model

model_fit [optional]: a statsmodel.api.OLS model after regressing y on X. If not provided, will be
generated from X, y
Function to reproduce the 4 base plots of an OLS model in R . Inputs . A numpy array or pandas dataframe of the features to use in building the linear regression model A numpy array or pandas series / dataframe of the target variable of the linear regression model model_fit [ optional ] : a statsmodel . api . OLS model after regressing y on X . If not provided will be generated from X y
def diagnostic_plots(X, y, model_fit=None):
  if not model_fit:
      model_fit = sm.OLS(y, sm.add_constant(X)).fit()
  dataframe = pd.concat([X, y], axis=1)
  model_fitted_y = model_fit.fittedvalues
  model_residuals = model_fit.resid
  model_norm_residuals = model_fit.get_influence().resid_studentized_internal
  model_norm_residuals_abs_sqrt = np.sqrt(np.abs(model_norm_residuals))
  model_abs_resid = np.abs(model_residuals)
  model_leverage = model_fit.get_influence().hat_matrix_diag
  model_cooks = model_fit.get_influence().cooks_distance[0]
  plot_lm_1 = plt.figure()
  plot_lm_1.axes[0] = sns.residplot(model_fitted_y, dataframe.columns[-1], data=dataframe,
                            lowess=True,
                            scatter_kws={'alpha': 0.5},
                            line_kws={'color': 'red', 'lw': 1, 'alpha': 0.8})
  plot_lm_1.axes[0].set_title('Residuals vs Fitted')
  plot_lm_1.axes[0].set_xlabel('Fitted values')
  plot_lm_1.axes[0].set_ylabel('Residuals');
  abs_resid = model_abs_resid.sort_values(ascending=False)
  abs_resid_top_3 = abs_resid[:3]
  for i in abs_resid_top_3.index:
      plot_lm_1.axes[0].annotate(i,
                                 xy=(model_fitted_y[i],
                                     model_residuals[i]));
  QQ = ProbPlot(model_norm_residuals)
  plot_lm_2 = QQ.qqplot(line='45', alpha=0.5, color='#4C72B0', lw=1)
  plot_lm_2.axes[0].set_title('Normal Q-Q')
  plot_lm_2.axes[0].set_xlabel('Theoretical Quantiles')
  plot_lm_2.axes[0].set_ylabel('Standardized Residuals');
  abs_norm_resid = np.flip(np.argsort(np.abs(model_norm_residuals)), 0)
  abs_norm_resid_top_3 = abs_norm_resid[:3]
  for r, i in enumerate(abs_norm_resid_top_3):
      plot_lm_2.axes[0].annotate(i,
                                 xy=(np.flip(QQ.theoretical_quantiles, 0)[r],
                                     model_norm_residuals[i]));
  plot_lm_3 = plt.figure()
  plt.scatter(model_fitted_y, model_norm_residuals_abs_sqrt, alpha=0.5);
  sns.regplot(model_fitted_y, model_norm_residuals_abs_sqrt,
              scatter=False,
              ci=False,
              lowess=True,
              line_kws={'color': 'red', 'lw': 1, 'alpha': 0.8});
  plot_lm_3.axes[0].set_title('Scale-Location')
  plot_lm_3.axes[0].set_xlabel('Fitted values')
  plot_lm_3.axes[0].set_ylabel('$\sqrt{|Standardized Residuals|}$');
  abs_sq_norm_resid = np.flip(np.argsort(model_norm_residuals_abs_sqrt), 0)
  abs_sq_norm_resid_top_3 = abs_sq_norm_resid[:3]
  for i in abs_norm_resid_top_3:
      plot_lm_3.axes[0].annotate(i,
                                 xy=(model_fitted_y[i],
                                     model_norm_residuals_abs_sqrt[i]));
  plot_lm_4 = plt.figure();
  plt.scatter(model_leverage, model_norm_residuals, alpha=0.5);
  sns.regplot(model_leverage, model_norm_residuals,
              scatter=False,
              ci=False,
              lowess=True,
              line_kws={'color': 'red', 'lw': 1, 'alpha': 0.8});
  plot_lm_4.axes[0].set_xlim(0, max(model_leverage)+0.01)
  plot_lm_4.axes[0].set_ylim(-3, 5)
  plot_lm_4.axes[0].set_title('Residuals vs Leverage')
  plot_lm_4.axes[0].set_xlabel('Leverage')
  plot_lm_4.axes[0].set_ylabel('Standardized Residuals');
  leverage_top_3 = np.flip(np.argsort(model_cooks), 0)[:3]
  for i in leverage_top_3:
      plot_lm_4.axes[0].annotate(i,
                                 xy=(model_leverage[i],
                                     model_norm_residuals[i]));
  p = len(model_fit.params) 
  graph(lambda x: np.sqrt((0.5 * p * (1 - x)) / x),
        np.linspace(0.001, max(model_leverage), 50),
        'Cook\'s distance') 
  graph(lambda x: np.sqrt((1 * p * (1 - x)) / x),
        np.linspace(0.001, max(model_leverage), 50)) 
  plot_lm_4.legend(loc='upper right');
def diagnostic_plots ( X , y , model_fit = None ) : if not model_fit : model_fit = sm . OLS ( y , sm . add_constant ( X ) ) . fit ( ) dataframe = pd . concat ( [ X , y ] , axis = 1 ) model_fitted_y = model_fit . fittedvalues model_residuals = model_fit . resid model_norm_residuals = model_fit . get_influence ( ) . resid_studentized_internal model_norm_residuals_abs_sqrt = np . sqrt ( np . abs ( model_norm_residuals ) ) model_abs_resid = np . abs ( model_residuals ) model_leverage = model_fit . get_influence ( ) . hat_matrix_diag model_cooks = model_fit . get_influence ( ) . cooks_distance [ 0 ] plot_lm_1 = plt . figure ( ) plot_lm_1 . axes [ 0 ] = sns . residplot ( model_fitted_y , dataframe . columns [ - 1 ] , data = dataframe , lowess = True , scatter_kws = { 'alpha' : 0.5 } , line_kws = { 'color' : 'red' , 'lw' : 1 , 'alpha' : 0.8 } ) plot_lm_1 . axes [ 0 ] . set_title ( 'Residuals vs Fitted' ) plot_lm_1 . axes [ 0 ] . set_xlabel ( 'Fitted values' ) plot_lm_1 . axes [ 0 ] . set_ylabel ( 'Residuals' ) ; abs_resid = model_abs_resid . sort_values ( ascending = False ) abs_resid_top_3 = abs_resid [ : 3 ] for i in abs_resid_top_3 . index : plot_lm_1 . axes [ 0 ] . annotate ( i , xy = ( model_fitted_y [ i ] , model_residuals [ i ] ) ) ; QQ = ProbPlot ( model_norm_residuals ) plot_lm_2 = QQ . qqplot ( line = '45' , alpha = 0.5 , color = '#4C72B0' , lw = 1 ) plot_lm_2 . axes [ 0 ] . set_title ( 'Normal Q-Q' ) plot_lm_2 . axes [ 0 ] . set_xlabel ( 'Theoretical Quantiles' ) plot_lm_2 . axes [ 0 ] . set_ylabel ( 'Standardized Residuals' ) ; abs_norm_resid = np . flip ( np . argsort ( np . abs ( model_norm_residuals ) ) , 0 ) abs_norm_resid_top_3 = abs_norm_resid [ : 3 ] for r , i in enumerate ( abs_norm_resid_top_3 ) : plot_lm_2 . axes [ 0 ] . annotate ( i , xy = ( np . flip ( QQ . theoretical_quantiles , 0 ) [ r ] , model_norm_residuals [ i ] ) ) ; plot_lm_3 = plt . figure ( ) plt . scatter ( model_fitted_y , model_norm_residuals_abs_sqrt , alpha = 0.5 ) ; sns . regplot ( model_fitted_y , model_norm_residuals_abs_sqrt , scatter = False , ci = False , lowess = True , line_kws = { 'color' : 'red' , 'lw' : 1 , 'alpha' : 0.8 } ) ; plot_lm_3 . axes [ 0 ] . set_title ( 'Scale-Location' ) plot_lm_3 . axes [ 0 ] . set_xlabel ( 'Fitted values' ) plot_lm_3 . axes [ 0 ] . set_ylabel ( '$\sqrt{|Standardized Residuals|}$' ) ; abs_sq_norm_resid = np . flip ( np . argsort ( model_norm_residuals_abs_sqrt ) , 0 ) abs_sq_norm_resid_top_3 = abs_sq_norm_resid [ : 3 ] for i in abs_norm_resid_top_3 : plot_lm_3 . axes [ 0 ] . annotate ( i , xy = ( model_fitted_y [ i ] , model_norm_residuals_abs_sqrt [ i ] ) ) ; plot_lm_4 = plt . figure ( ) ; plt . scatter ( model_leverage , model_norm_residuals , alpha = 0.5 ) ; sns . regplot ( model_leverage , model_norm_residuals , scatter = False , ci = False , lowess = True , line_kws = { 'color' : 'red' , 'lw' : 1 , 'alpha' : 0.8 } ) ; plot_lm_4 . axes [ 0 ] . set_xlim ( 0 , max ( model_leverage ) + 0.01 ) plot_lm_4 . axes [ 0 ] . set_ylim ( - 3 , 5 ) plot_lm_4 . axes [ 0 ] . set_title ( 'Residuals vs Leverage' ) plot_lm_4 . axes [ 0 ] . set_xlabel ( 'Leverage' ) plot_lm_4 . axes [ 0 ] . set_ylabel ( 'Standardized Residuals' ) ; leverage_top_3 = np . flip ( np . argsort ( model_cooks ) , 0 ) [ : 3 ] for i in leverage_top_3 : plot_lm_4 . axes [ 0 ] . annotate ( i , xy = ( model_leverage [ i ] , model_norm_residuals [ i ] ) ) ; p = len ( model_fit . params ) graph ( lambda x : np . sqrt ( ( 0.5 * p * ( 1 - x ) ) / x ) , np . linspace ( 0.001 , max ( model_leverage ) , 50 ) , 'Cook\'s distance' ) graph ( lambda x : np . sqrt ( ( 1 * p * ( 1 - x ) ) / x ) , np . linspace ( 0.001 , max ( model_leverage ) , 50 ) ) plot_lm_4 . legend ( loc = 'upper right' ) ;
Function to reproduce the 4 base plots of an OLS model in R.
Inputs:
Function to reproduce the 4 base plots of an OLS model in R . Inputs :
"""
  Function to reproduce the 4 base plots of an OLS model in R.

  ---
  Inputs:

  X: A numpy array or pandas dataframe of the features to use in building the linear regression model

  y: A numpy array or pandas series/dataframe of the target variable of the linear regression model

  model_fit [optional]: a statsmodel.api.OLS model after regressing y on X. If not provided, will be
                        generated from X, y
  """ # create dataframe from X, y for easier plot handling # model values # model residuals # normalized residuals # absolute squared normalized residuals # absolute residuals # leverage, from statsmodels internals # cook's distance, from statsmodels internals # annotations # annotations # annotations # annotations # number of model parameters # 0.5 line # 1 line
X
y
model_fit