docs: def _write_comparison_plot_table(spid, models, options, core_results,
                                   fit_results):
      """
      Notes
      -----
      Only applies to analysis using functions from empirical in which models are
      also given.
  
      """
      # TODO: Clean up sorting, may not work if SAR x out of order, e.g.
  
      is_curve = 'x' in core_results[0][1]
      df = core_results[spid][1]
      df.rename(columns={'y': 'empirical'}, inplace=True)
  
      # If distribution, need to sort values so will match sorted rank in fits
      if not is_curve:
          x = np.arange(len(df)) + 1
          df = df.sort(columns='empirical')
          df.insert(0, 'x', x[::-1])
  
      # Add residual column for each model
      for model in models:
          fit_result = fit_results[spid][model]
          df[model] = fit_result[1]
          df[model + "_residual"] = df[model] - df['empirical']
  
      # If curve, sort now for plotting purposes
      if is_curve:
          df = df.sort(columns='x')
  
      # Set up file paths
      f_path = _get_file_path(spid, options, 'data_models.csv')
      p_path = _get_file_path(spid, options, 'data_models.pdf')
  
      # Save table
      df.to_csv(f_path, index=False, float_format='%.4f')  # Table
  
      # Save plot
      fig, (ax1, ax2) = plt.subplots(1, 2)
  
      ax1.scatter(df['x'], df['empirical'], color='k')
      ax1.plot(df['x'], df[models])
      ax1.legend(models + ['empirical'], loc='best')
      ax1.set_xlabel('x')
      ax1.set_ylabel('value')
  
      ax2.hlines(0, np.min(df['x']), np.max(df['x']))
      ax2.plot(df['x'], df[[x + '_residual' for x in models]])
      ax2.legend(models + ['empirical'], loc='best')
      ax2.set_xlabel('x')
      ax2.set_ylabel('residual')
      ax2.set_xlim(ax1.get_xlim())
      ax2.set_ylim(min(ax2.get_ylim()[0], -1), max(ax2.get_ylim()[1], 1))
  
      if options.get('log_y', None):
          ax1.set_yscale('log')
          ax2.set_yscale('symlog', linthreshy=1)
      if options.get('log_x', None):
          ax1.set_xscale('log')
          ax2.set_xscale('log')
  
      if not options.get('log_x', None) and not options.get('log_y', None):
          ax1.set_ylim(bottom=0)
          ax1.set_xlim(left=0)
          ax1 = _pad_plot_frame(ax1)
          ax2 = _pad_plot_frame(ax2)
  
      with warnings.catch_warnings():
          warnings.simplefilter("ignore")
          fig.tight_layout()
      fig.savefig(p_path)
  
      plt.close('all')