Spectral Fit Residuals as an Indicator to Increase Model Complexity

Spectral fitting of X-ray data has been usually done by minimizing a statistic like χ² or cstat (Avni 1976; Cash 1979). The value of the statistic is useful to evaluate the goodness of the fit, which is also used as a stopping rule in evaluating model complexity. Fits are deemed acceptable if ${\chi }^{2}/\nu \lt 1+3\cdot \sqrt{\tfrac{2}{\nu }}$ (where ν is the degrees of freedom in the fit) or Δcstat<3 · σ _cstat (where ${\sigma }_{{\mathtt{cstat}}}^{2}$ is the estimated variance for the nominal model as given by Kaastra 2017). When the fitting statistics or F-test values drop below the threshold defined by the null distributions, it is commonly acknowledged that further increases in model complexity (or the number of free parameters in the fit) are statistically unsupported. However, these measures are global measures of fit, and leave contiguous ranges in data space with correlated deviations in the residuals. Here we introduce summary statistics that identify the presence of such structures in the residuals and describe how to use them to decide whether the model complexity should be increased.

We use a Chandra/ACIS-S spectrum of the corona of an exoplanet hosting star (HD 179949, ObsID 6122; Acharya et al. 2023) to illustrate the measures (see Figure 1). The data and the spectral fits made with different emission models are shown in the left column of the Figure. The models increase in complexity from top to bottom (see Acharya et al. 2023 for details), with none of the statistical measures (cstat, pct_CuSum, p_area) acceptable for the simplest model, cstat being acceptable for the model in the middle row, and all measures acceptable for the model in the bottom row. Below, we define pct_CuSum and p_area, describe their rationale, and how to compute them and their null distributions.

Zoom In Zoom Out Reset image size

Spectral fit residuals can be expressed as

$\begin{eqnarray*}&&\{{r}_{i}(\hat{\theta })=({m}_{i}(\hat{\theta })-{c}_{i}),\,\,\,i=1,\ldots ,N\},\end{eqnarray*}$

with c_i the observed counts in bin i of N and m_i the predicted model counts in the same bin, generated for an astrophysical model defined by the best-fit parameters $\hat{\theta }$. The cumulative sum (CuSum) at the jth bin is

$\begin{eqnarray}&&\mathrm{CuSum}{(\hat{\theta })}_{j}=\displaystyle \sum _{i=1}^{j}{r}_{i}(\hat{\theta }).\end{eqnarray} \tag{ 1 }$

While the ordering of indices can be reversed, CuSum preserves the sequential order of the bins, thus incorporating extra information ignored in χ² or cstat statistics. The null distribution needed to evaluate the quality of the CuSum of a best-fit can be built in one of two ways:

− generate K mock data sets using the fake_pha() function in Sherpa, and fit using the same model, such that for each mock data set k, we obtain corresponding best-fit parameters θ^k ; or

− obtain K post burn-in draws {θ^k ,  k = 1,...,K} as iterations during a Markov Chain Monte Carlo (MCMC) fit using the get_draws() tool in Sherpa (based on BLoCXS; van Dyk et al. 2001).

Each θ^k yields a CuSum array $\{\mathrm{CuSum}{({\theta }^{k})}_{j},\,j\,=\,1,\ldots ,N\}$, and the resulting sample

$\begin{eqnarray}&&\{\mathrm{CuSum}{({\theta }^{k})}_{j},\,k=1,\ldots ,K,\,j=1,\ldots ,N\}\end{eqnarray} \tag{ 2 }$

defines the null distribution.

This null distribution characterizes the strength of the deviations present in $\mathrm{CuSum}{(\hat{\theta })}_{j}$. We construct two statistics that probe deviation at different scales, pct_CuSum to detect biases in the model continuum and p_area to detect presence/absence of narrow lines.

• 1.  
  pct_CuSum: We compute the equal-tail 90% point-wise range $[{C}_{J}^{05},{C}_{J}^{95}]$ of $\mathrm{CuSum}{({\theta }^{k})}_{j\,=\,J}$ for each detector channel J in the specified passband. Then we evaluate the percentage of bins, pct_CuSum, where the observed CuSum exceeds the 90% bounds, i.e., if n channels have $\mathrm{CuSum}{(\hat{\theta })}_{J}\lt {C}_{J}^{05}$ or $\mathrm{CuSum}{(\hat{\theta })}_{J}\gt {C}_{J}^{95}$, then pct${}_{\mathrm{CuSum}}=100\cdot \tfrac{n}{N}$. If pct_CuSum ≳ 10%, the model used is considered inadequate and we suggest that more complex models with more free parameters should be explored. If, on the other hand, pct_CuSum ≪ 10% this is a sign of overfitting, and thus less complex models are favored (see the middle column of Figure 1).
• 2.  
  p_area: For each of the k = 1, ..., K CuSum arrays in the null, we calculate the total area that falls outside the 5%–95% bounds across all the channels,

  $\begin{eqnarray}&&{\mathrm{area}}^{k}=\displaystyle \sum _{J}{Z}_{J}^{+}{(\mathrm{CuSum}{({\theta }^{k})}_{J}-{C}_{J}^{95})+{Z}_{J}^{-}({C}_{J}^{05}-\mathrm{CuSum}({\theta }^{k})}_{J}),\end{eqnarray} \tag{ 3 }$

  where the indicator indices Z⁺ = Z⁻ = 1 if the corresponding terms in the summation are +ve, and 0 otherwise. We then compute a similar area measure for the best-fit spectrum,

  $\begin{eqnarray}&&\mathrm{area}{(\hat{\theta })=\displaystyle \sum _{J}{Z}_{J}^{+}{(\mathrm{CuSum}(\hat{\theta })}_{J}-{C}_{J}^{95})+{Z}_{J}^{-}({C}_{J}^{05}-\mathrm{CuSum}(\hat{\theta })}_{J}).\end{eqnarray} \tag{ 4 }$

  The set {area^k ,  k = 1,...,K} now forms a null distribution of excess areas to which the value of $\mathrm{area}(\hat{\theta })$ can be compared to obtain a p-value. The statistic p_area is the probability that $\mathrm{area}(\hat{\theta })\gt {\mathrm{area}}^{k}$. We consider the best-fit model an inadequate fit if p_area ≪ 0.05, as an indication that a large deviation is present in $\mathrm{CuSum}{(\hat{\theta })}_{j}$ (see the right column of Figure 1). This process is similar to the posterior predictive p-value calibration procedure described by Protassov et al. (2002).

In Figure 1, we show the results of the residuals analysis for three best-fit models for an exemplar dataset of HD 179949 (Acharya et al. 2023). The top row represents a single temperature APEC model with variable metallicity (1m). The Δcstat = + 3.3 · σ _cstat , indicates a poor fit. This is supported by the pct_CuSum = 40.6% and a p-value = 0.0. The middle row represents a 2-temperature model with metallicities scaling together (2m). The Δcstat = + 2.1 · σ _cstat , indicates an acceptable model. The argument for improvements is solidified by pct_CuSum = 34.2% and p-value = 0.0. The bottom row represents a 2-temperature model with abundances grouped by First Ionization Potential and varying in several groups (2v). This gives Δcstat = + 1.3 · σ _cstat , pct_CuSum = 13.5% and p-value = 0.09. This model is therefore accepted as the best representation of the data.

pct_CuSum is sensitive to broad differences such as would be present in a biased fit of the continuum, and p_area to small scale changes such as appear when spectral lines or abundance differences are mis-estimated. Our method guards against overfitting to statistical fluctuations in the data by specifying criteria for acceptability that are tied to the null distributions. We note that our arguments and method are general, and can be extended to any nonlinear fitting scenario to supplement fitting procedures.

We thank Doug Burke (CXC/SDS) for suggestions with the code development.

Code Availability: The code as described here is publicly available on Zenodo (Acharya & Kashyap 2023). See https://github.com/anshuman1998/csresid for updates.