docs: def nstar(self, image, star_groups):
          """
          Fit, as appropriate, a compound or single model to the given
          ``star_groups``. Groups are fitted sequentially from the
          smallest to the biggest. In each iteration, ``image`` is
          subtracted by the previous fitted group.
  
          Parameters
          ----------
          image : numpy.ndarray
              Background-subtracted image.
          star_groups : `~astropy.table.Table`
              This table must contain the following columns: ``id``,
              ``group_id``, ``x_0``, ``y_0``, ``flux_0``.  ``x_0`` and
              ``y_0`` are initial estimates of the centroids and
              ``flux_0`` is an initial estimate of the flux. Additionally,
              columns named as ``<param_name>_0`` are required if any
              other parameter in the psf model is free (i.e., the
              ``fixed`` attribute of that parameter is ``False``).
  
          Returns
          -------
          result_tab : `~astropy.table.Table`
              Astropy table that contains photometry results.
          image : numpy.ndarray
              Residual image.
          """
  
          result_tab = Table()
          for param_tab_name in self._pars_to_output.keys():
              result_tab.add_column(Column(name=param_tab_name))
  
          unc_tab = Table()
          for param, isfixed in self.psf_model.fixed.items():
              if not isfixed:
                  unc_tab.add_column(Column(name=param + "_unc"))
  
          y, x = np.indices(image.shape)
  
          star_groups = star_groups.group_by('group_id')
          for n in range(len(star_groups.groups)):
              group_psf = get_grouped_psf_model(self.psf_model,
                                                star_groups.groups[n],
                                                self._pars_to_set)
              usepixel = np.zeros_like(image, dtype=np.bool)
  
              for row in star_groups.groups[n]:
                  usepixel[overlap_slices(large_array_shape=image.shape,
                                          small_array_shape=self.fitshape,
                                          position=(row['y_0'], row['x_0']),
                                          mode='trim')[0]] = True
  
              fit_model = self.fitter(group_psf, x[usepixel], y[usepixel],
                                      image[usepixel])
              param_table = self._model_params2table(fit_model,
                                                     len(star_groups.groups[n]))
              result_tab = vstack([result_tab, param_table])
  
              if 'param_cov' in self.fitter.fit_info.keys():
                  unc_tab = vstack([unc_tab,
                                    self._get_uncertainties(
                                        len(star_groups.groups[n]))])
              try:
                  from astropy.nddata.utils import NoOverlapError
              except ImportError:
                  raise ImportError("astropy 1.1 or greater is required in "
                                    "order to use this class.")
              # do not subtract if the fitting did not go well
              try:
                  image = subtract_psf(image, self.psf_model, param_table,
                                       subshape=self.fitshape)
              except NoOverlapError:
                  pass
  
          if 'param_cov' in self.fitter.fit_info.keys():
              result_tab = hstack([result_tab, unc_tab])
  
          return result_tab, image