task: metadata: install: pip install -e .[all]
patch: diff --git a/CHANGES.rst b/CHANGES.rst
index 2a699374..ea0fea0b 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -27,6 +27,9 @@ New Features
 
   - Added a progress bar for fitting PSF photometry [#1426]
 
+  - Added a ``subshape`` keyword to the PSF-fitting classes to define
+    the shape over which the PSF is subtracted. [#1477]
+
 - ``photutils.segmentation``
 
   - Added the ability to slice ``SegmentationImage`` objects. [#1413]
diff --git a/docs/whats_new/1.6.rst b/docs/whats_new/1.6.rst
index ad1aaccf..accea92c 100644
--- a/docs/whats_new/1.6.rst
+++ b/docs/whats_new/1.6.rst
@@ -81,6 +81,17 @@ The progress bars require installation of the `tqdm
 <https://tqdm.github.io/>`_ optional dependency.
 
 
+New subshape keyword in PSF fitting
+===================================
+
+A new ``subshape`` keyword was added to the PSF-fitting classes to
+define the shape over which the PSF is subtracted when computing the
+residual image. Previously, the PSF-subtraction region was always
+defined by the ``fitshape`` keyword. By default (and for backwards
+compatibility), ``subshape`` is set to `None`, which means the
+``fitshape`` value will be used.
+
+
 Other changes
 =============
 
diff --git a/photutils/psf/photometry.py b/photutils/psf/photometry.py
index 60646e39..49fc29ed 100644
--- a/photutils/psf/photometry.py
+++ b/photutils/psf/photometry.py
@@ -68,12 +68,11 @@ class BasicPSFPhotometry:
         ``y_0`` and ``flux``. `~photutils.psf.prepare_psf_model` can be
         used to prepare any 2D model to match this assumption.
     fitshape : int or length-2 array-like
-        Rectangular shape around the center of a star which will be
-        used to collect the data to do the fitting. Can be an integer
-        to be the same along both axes. For example, 5 is the same as
-        (5, 5), which means to fit only at the following relative pixel
-        positions: [-2, -1, 0, 1, 2]. Each element of ``fitshape`` must
-        be an odd number.
+        Rectangular shape around the center of a star that will be
+        used to define the PSF-fitting region. If ``fitshape`` is a
+        scalar then a square shape of size ``fitshape`` will be used.
+        If ``fitshape`` has two elements, they must be in ``(ny, nx)``
+        order. Each element of ``fitshape`` must be an odd number.
     finder : callable or instance of any \
             `~photutils.detection.StarFinderBase` subclasses or None
         ``finder`` should be able to identify stars, i.e., compute a
@@ -103,6 +102,13 @@ class BasicPSFPhotometry:
         List of additional columns for parameters derived by any of the
         intermediate fitting steps (e.g., ``finder``), such as roundness
         or sharpness.
+    subshape : `None`, int, or length-2 array-like
+        Rectangular shape around the center of a star that will be
+        used to define the PSF-subtraction region. If `None`, then
+        ``fitshape`` will be used. If ``subshape`` is a scalar then a
+        square shape of size ``subshape`` will be used. If ``subshape``
+        has two elements, they must be in ``(ny, nx)`` order. Each
+        element of ``subshape`` must be an odd number.
 
     Notes
     -----
@@ -128,7 +134,7 @@ class BasicPSFPhotometry:
 
     def __init__(self, group_maker, bkg_estimator, psf_model, fitshape, *,
                  finder=None, fitter=LevMarLSQFitter(), aperture_radius=None,
-                 extra_output_cols=None):
+                 extra_output_cols=None, subshape=None):
         self.group_maker = group_maker
         self.bkg_estimator = bkg_estimator
         self.psf_model = psf_model
@@ -140,6 +146,7 @@ class BasicPSFPhotometry:
         self._pars_to_output = None
         self._residual_image = None
         self._extra_output_cols = extra_output_cols
+        self.subshape = subshape
 
     @property
     def fitshape(self):
@@ -167,6 +174,36 @@ class BasicPSFPhotometry:
             raise ValueError('fitshape must have two dimensions, '
                              f'received fitshape={value}')
 
+    @property
+    def subshape(self):
+        return self._subshape
+
+    @subshape.setter
+    def subshape(self, value):
+        if value is None:
+            self._subshape = self._fitshape
+            return
+
+        value = np.asarray(value)
+
+        # assume a lone value should mean both axes
+        if value.shape == ():
+            value = np.array((value, value))
+
+        if value.size == 2:
+            if np.all(value) > 0:
+                if np.all(value % 2) == 1:
+                    self._subshape = tuple(value)
+                else:
+                    raise ValueError('subshape must be odd integer-valued, '
+                                     f'received subshape={value}')
+            else:
+                raise ValueError('subshape must have positive elements, '
+                                 f'received subshape={value}')
+        else:
+            raise ValueError('subshape must have two dimensions, '
+                             f'received subshape={value}')
+
     @property
     def aperture_radius(self):
         return self._aperture_radius
@@ -498,7 +535,7 @@ class BasicPSFPhotometry:
             # do not subtract if the fitting did not go well
             try:
                 image = subtract_psf(image, self.psf_model, param_table,
-                                     subshape=self.fitshape)
+                                     subshape=self.subshape)
             except NoOverlapError:
                 pass
 
@@ -665,12 +702,11 @@ class IterativelySubtractedPSFPhotometry(BasicPSFPhotometry):
         ``y_0`` and ``flux``. `~photutils.psf.prepare_psf_model` can be
         used to prepare any 2D model to match this assumption.
     fitshape : int or length-2 array-like
-        Rectangular shape around the center of a star which will be
-        used to collect the data to do the fitting. Can be an integer
-        to be the same along both axes. For example, 5 is the same as
-        (5, 5), which means to fit only at the following relative pixel
-        positions: [-2, -1, 0, 1, 2]. Each element of ``fitshape`` must
-        be an odd number.
+        Rectangular shape around the center of a star that will be
+        used to define the PSF-fitting region. If ``fitshape`` is a
+        scalar then a square shape of size ``fitshape`` will be used.
+        If ``fitshape`` has two elements, they must be in ``(ny, nx)``
+        order. Each element of ``fitshape`` must be an odd number.
     finder : callable or instance of any \
             `~photutils.detection.StarFinderBase` subclasses
         ``finder`` should be able to identify stars, i.e., compute a
@@ -701,6 +737,13 @@ class IterativelySubtractedPSFPhotometry(BasicPSFPhotometry):
         List of additional columns for parameters derived by any of the
         intermediate fitting steps (e.g., ``finder``), such as roundness
         or sharpness.
+    subshape : `None`, int, or length-2 array-like
+        Rectangular shape around the center of a star that will be
+        used to define the PSF-subtraction region. If `None`, then
+        ``fitshape`` will be used. If ``subshape`` is a scalar then a
+        square shape of size ``subshape`` will be used. If ``subshape``
+        has two elements, they must be in ``(ny, nx)`` order. Each
+        element of ``subshape`` must be an odd number.
 
     Notes
     -----
@@ -720,12 +763,13 @@ class IterativelySubtractedPSFPhotometry(BasicPSFPhotometry):
 
     def __init__(self, group_maker, bkg_estimator, psf_model, fitshape,
                  finder, *, fitter=LevMarLSQFitter(), niters=3,
-                 aperture_radius=None, extra_output_cols=None):
+                 aperture_radius=None, extra_output_cols=None, subshape=None):
 
         super().__init__(group_maker, bkg_estimator, psf_model, fitshape,
                          finder=finder, fitter=fitter,
                          aperture_radius=aperture_radius,
-                         extra_output_cols=extra_output_cols)
+                         extra_output_cols=extra_output_cols,
+                         subshape=subshape)
         self.niters = niters
 
     @property
@@ -953,12 +997,11 @@ class DAOPhotPSFPhotometry(IterativelySubtractedPSFPhotometry):
         ``y_0`` and ``flux``. `~photutils.psf.prepare_psf_model` can be
         used to prepare any 2D model to match this assumption.
     fitshape : int or length-2 array-like
-        Rectangular shape around the center of a star which will be
-        used to collect the data to do the fitting. Can be an integer
-        to be the same along both axes. For example, 5 is the same as
-        (5, 5), which means to fit only at the following relative pixel
-        positions: [-2, -1, 0, 1, 2]. Each element of ``fitshape`` must
-        be an odd number.
+        Rectangular shape around the center of a star that will be
+        used to define the PSF-fitting region. If ``fitshape`` is a
+        scalar then a square shape of size ``fitshape`` will be used.
+        If ``fitshape`` has two elements, they must be in ``(ny, nx)``
+        order. Each element of ``fitshape`` must be an odd number.
     sigma : float, optional
         Number of standard deviations used to perform sigma clip with a
         `astropy.stats.SigmaClip` object.
@@ -1001,6 +1044,13 @@ class DAOPhotPSFPhotometry(IterativelySubtractedPSFPhotometry):
         List of additional columns for parameters derived by any of the
         intermediate fitting steps (e.g., ``finder``), such as roundness
         or sharpness.
+    subshape : `None`, int, or length-2 array-like
+        Rectangular shape around the center of a star that will be
+        used to define the PSF-subtraction region. If `None`, then
+        ``fitshape`` will be used. If ``subshape`` is a scalar then a
+        square shape of size ``subshape`` will be used. If ``subshape``
+        has two elements, they must be in ``(ny, nx)`` order. Each
+        element of ``subshape`` must be an odd number.
 
     Notes
     -----
@@ -1022,7 +1072,8 @@ class DAOPhotPSFPhotometry(IterativelySubtractedPSFPhotometry):
                  *, sigma=3.0, ratio=1.0, theta=0.0, sigma_radius=1.5,
                  sharplo=0.2, sharphi=1.0, roundlo=-1.0, roundhi=1.0,
                  fitter=LevMarLSQFitter(),
-                 niters=3, aperture_radius=None, extra_output_cols=None):
+                 niters=3, aperture_radius=None, extra_output_cols=None,
+                 subshape=None):
 
         self.crit_separation = crit_separation
         self.threshold = threshold
@@ -1048,4 +1099,5 @@ class DAOPhotPSFPhotometry(IterativelySubtractedPSFPhotometry):
                          psf_model=psf_model, fitshape=fitshape,
                          finder=finder, fitter=fitter, niters=niters,
                          aperture_radius=aperture_radius,
-                         extra_output_cols=extra_output_cols)
+                         extra_output_cols=extra_output_cols,
+                         subshape=subshape)

problem_statement: PSF Photometry: Need different box size for fitting vs. PSF subtraction.
In the PSF photometry classes, currently the fitshape parameter controls both the box size to fit the PSF to a source and the size of the region to use when subtracting that scaled PSF model. Often, we want to use a fit box centered on the core (where SNR is maxed and confusion is less likely) but we want to subtract a much bigger region, especially for bright stars. Currently, this is not possible. 

The code implementation is in photutils.psf.photometry:
```
def star() 
...          
            try:
                image = subtract_psf(image, self.psf_model, param_table,
                                     subshape=self.fitshape)
            except NoOverlapError:
                pass

```
I would recommend allowing subshape to be set as a variable on the object instead of using fitshape by default.
repo: astropy/photutils
test_cmd: pytest --no-header -rA --tb=line --color=no -p no:cacheprovider -W ignore::DeprecationWarning
test_patch: diff --git a/photutils/psf/tests/test_photometry.py b/photutils/psf/tests/test_photometry.py
index 522a55af..ab164bfa 100644
--- a/photutils/psf/tests/test_photometry.py
+++ b/photutils/psf/tests/test_photometry.py
@@ -17,7 +17,7 @@ from numpy.testing import assert_allclose, assert_array_equal, assert_equal
 from photutils.background import MMMBackground, StdBackgroundRMS
 from photutils.datasets import (make_gaussian_prf_sources_image,
                                 make_noise_image)
-from photutils.detection import DAOStarFinder
+from photutils.detection import DAOStarFinder, IRAFStarFinder
 from photutils.psf.groupstars import DAOGroup
 from photutils.psf.models import FittableImageModel, IntegratedGaussianPRF
 from photutils.psf.photometry import (BasicPSFPhotometry, DAOPhotPSFPhotometry,
@@ -830,3 +830,72 @@ def test_photometry_mask_nan():
                       'non-finite values'):
         tbl2 = psfphot(data, init_guesses=init)
         assert tbl == tbl2
+
+
+@pytest.mark.skipif(not HAS_SCIPY, reason='scipy is required')
+def test_photometry_subshape():
+    size = 21
+    cen = (size - 1) // 2
+    sigma = 2.0
+    sources = Table()
+    sources['flux'] = [1]
+    sources['x_0'] = [cen]
+    sources['y_0'] = [cen]
+    sources['sigma'] = [sigma]
+    psf = make_gaussian_prf_sources_image((size, size), sources)
+    psf_model = FittableImageModel(psf)
+
+    sources = Table()
+    sources['flux'] = [2000, 1000]
+    sources['x_0'] = [18, 7]
+    sources['y_0'] = [17, 25]
+    sources['sigma'] = [sigma, sigma]
+    shape = (33, 33)
+    image = make_gaussian_prf_sources_image(shape, sources)
+
+    daogroup = DAOGroup(crit_separation=8)
+    mmm_bkg = MMMBackground()
+    iraffind = IRAFStarFinder(threshold=10, fwhm=5, roundlo=-1, minsep_fwhm=1)
+    fitter = LevMarLSQFitter()
+
+    pobj1 = BasicPSFPhotometry(finder=iraffind, group_maker=daogroup,
+                               bkg_estimator=mmm_bkg, psf_model=psf_model,
+                               fitter=fitter, fitshape=(7, 7),
+                               aperture_radius=5, subshape=(3, 3))
+    pobj2 = BasicPSFPhotometry(finder=iraffind, group_maker=daogroup,
+                               bkg_estimator=mmm_bkg, psf_model=psf_model,
+                               fitter=fitter, fitshape=(7, 7),
+                               aperture_radius=5, subshape=5)
+    pobj3 = BasicPSFPhotometry(finder=iraffind, group_maker=daogroup,
+                               bkg_estimator=mmm_bkg, psf_model=psf_model,
+                               fitter=fitter, fitshape=(7, 7),
+                               aperture_radius=5, subshape=None)
+    pobj4 = BasicPSFPhotometry(finder=iraffind, group_maker=daogroup,
+                               bkg_estimator=mmm_bkg, psf_model=psf_model,
+                               fitter=fitter, fitshape=(7, 7),
+                               aperture_radius=5, subshape=7)
+
+    _ = pobj1(image)
+    _ = pobj2(image)
+    _ = pobj3(image)
+    _ = pobj4(image)
+    resid1 = pobj1.get_residual_image()
+    resid2 = pobj2.get_residual_image()
+    resid3 = pobj3.get_residual_image()
+    resid4 = pobj4.get_residual_image()
+    assert np.sum(resid1) > np.sum(resid2)
+    assert_allclose(np.sum(resid3), np.sum(resid4))
+
+
+@pytest.mark.skipif(not HAS_SCIPY, reason='scipy is required')
+def test_subshape_invalid():
+    basic_phot_obj = make_psf_photometry_objs()[0]
+
+    with pytest.raises(ValueError):
+        basic_phot_obj.subshape = (2, 2)
+    with pytest.raises(ValueError):
+        basic_phot_obj.subshape = 2
+    with pytest.raises(ValueError):
+        basic_phot_obj.subshape = (-1, 0)
+    with pytest.raises(ValueError):
+        basic_phot_obj.subshape = (3, 3, 3)