Overfitting Affects the Reliability of Radial Velocity Mass Estimates of the V1298 Tau Planets

Planet formation is an uncertain process. Giant planets are thought to form with large radii, inflated due to trapped heat, then cool and contract over the first few hundred million years of their lives (Marley et al. 2007). However, the accretion efficiency of the formation process, which sets the planets' initial entropies and radii, spans orders of magnitude of uncertainty. The processes sculpting the postformation masses and radii of smaller terrestrial exoplanets are also uncertain. Young, terrestrial planets also have uncertain initial entropies, and for highly irradiated planets, the unknown rate of photoevaporation (itself due to uncertainties in a planet's migration history, among other physical unknowns) during and after formation compounds this ambiguity (Lopez et al. 2012; Owen & Wu 2013; Chen & Rogers 2016; Owen 2020).

Measuring the masses, radii, and ages of newly formed planets presents a path forward (Owen 2020). Young moving groups provide rigorous age constraints, and relatively model-independent methods of measuring planetary radii exist for both young directly imaged and transiting planets (for transiting planets in particular, only stellar-radius model dependencies impact the inferred planetary radius). However, in both situations, few model-independent mass measurements exist. For transiting planets, there are two complementary methods for measuring planetary masses: transit timing variations (TTVs), and stellar radial velocity (RV) time series.

Measuring the RV masses of young planets is a difficult task, so some advocate to rely on TTVs alone to measure the masses of young planets. However, not all planets transit, and only planets in multiplanet systems at or near mean motion resonances exhibit TTVs (Fabrycky et al. 2014). Even in systems that do, individual TTV mass posteriors are often covariant, since TTVs to first order constrain the planetary mass ratio (Lithwick et al. 2012; Petigura et al. 2020). In an ideal scenario, both RVs and TTVs would be used to constrain jointly the planetary masses in a given system, reducing the posterior uncertainty and TTV degeneracies.

As the instrumental errors of extremely precise RV instruments approach 10 cm s⁻¹, and as the RV community begins to target more active stars, accurately modeling astrophysical noise is becoming more and more critical. Young stars present a particular challenge. These are highly magnetically active (Johns-Krull 2007), with starspots that occupy significant fractions of the stellar surface and induce RV variations on the order of ∼1 km s⁻¹ (Saar & Donahue 1997). These RV variations are hundreds of times larger than the activity signals of older quiet stars typically targeted by RV surveys and complicate the detection of planet-induced Doppler shifts from even close-in Jupiter-mass planets (e.g., Huerta et al. 2008; Prato et al. 2008).

Other assumptions and/or information can be leveraged to model the activity signal, even if the signal is not easily understandable from the RVs themselves. A widely used practice involves independently constraining the rotation period from a photometric time series, then using an informed prior on the rotation period to model the RVs (e.g., Grunblatt et al. 2015). Other related examples include specifying a quasi-periodic kernel for a Gaussian process regression (GPR) model (i.e., assuming that the stellar activity has a quasi-periodic form), or modeling the RVs jointly with other data sets. The latter approach achieves better model constraints either by explicitly modeling the relationship between the data sets (e.g., Rajpaul et al. 2015) or by sharing hyperparameters between data sets (e.g., Grunblatt et al. 2015; López-Morales et al. 2016).

As is true for every model-fitting process, misspecifying the stellar activity model (i.e., fitting a model that is not representative of the process that generated the data) or allowing too many effective degrees of freedom can lead to overfitting.

Overfitting is a concept ubiquitous in machine learning, and in particular is often used to determine when a model has been optimally trained. One algorithm for determining whether a model is overfitting is as follows: ²³ divide the data into a "training" set and an "evaluation" set (a common split is 80%/20%), and begin optimizing the model using just the training set. At each optimization step, calculate the goodness-of-fit metric for the model on the evaluation set, which is otherwise omitted from the training process altogether. This method of evaluating a model's ability to predict new, or "out-of-sample," data successfully is known as cross-validation (CV).

The classic observed behavior is that the goodness-of-fit metrics for both the training and evaluation set improve as the model fits the training data better and better. At a certain point, the model begins to overfit to the training data, and the goodness-of-fit metric for the evaluation data worsens. This is because the model parameters have begun to reproduce the noise in the training set, at the expense of reproducing the signal common to both data sets. A model that is overfitting, then, can be defined as one that predicts the observations in a training set better than those in an evaluation set. An overfitting model fits aspects of the data that are not predictable or common to the entire data set, e.g., noise.

The optimally trained model is selected not by its performance relative to the training data, but by its performance relative to the evaluation data, which was omitted from the training process altogether. Making an analogy to Bayesian model comparison, we could imagine a similar process where the goodness of fit is evaluated for an "evaluation set" left out of the training process (i.e., posterior computation using a Markov Chain Monte Carlo (MCMC) algorithm, nested sampling, etc.) for a series of models. One benefit of this method over, e.g., formal Bayesian model comparison is that it also provides an easily interpretable absolute metric for how well the model fits the data: if the evaluation set goodness of fit is significantly worse than that of the training set, we know the model is misspecified, even if it has (comparatively) the lowest Bayesian evidence.

In this study, we apply the CV technique as defined above to evaluate the predictiveness of one particular model fit to one particular star. This is intended as a case study, aiming to inspire further investigation into the extent of and causes of overfitting in RV modeling of young, active stars.

V1298 Tauri (hereafter V1298 Tau) is a young system of four ≳0.5 R_J planets transiting a K-type pre-main-sequence (PMS) star (David et al. 2019a, 2019b). Very few transiting planets have been discovered around PMS stars (other notable systems being AU Mic, Plavchan et al. 2020; Cale et al. 2021; Klein et al. 2022; Zicher et al. 2022; K2-33, David et al. 2016; DS Tuc, Benatti et al. 2019; Newton et al. 2019; HIP 67522, Rizzuto et al. 2020; Kepler 1627A, Bouma et al. 2022; and TOI 1227, Mann et al. 2022). David et al. (2019b) reported the discovery of V1298 Tau b, a 0.9 R_J object with an orbital period of 24 days. David et al. (2019a) discovered three additional planets in the system: V1298 Tau c at 8 days, V1298 Tau d at 12 days, and a single-transiting object, V1298 Tau e.

It is unclear from their radii (0.5–1 R_J ) alone whether these planets are gas giants that contracted rapidly after forming, or are young terrestrial or mini-Neptune planets, which will lose a large fraction of their atmospheres to photoevaporation (Owen & Wu 2013) and/or core-powered mass loss (Ginzburg et al. 2018) as they age.

Suárez Mascareño et al. (2021), hereafter SM21, presented over 100 RVs of the system with four different instruments, and used a suite of GP models in an effort to isolate the RV signals of the outer two transiting planets (b and e), finding masses of 0.64 and 1.16 M_J , respectively. Combined with radii of 0.91 and 0.78 R_J , this implied high densities of 1.2 and 3.6 g cm⁻³, respectively. Such high densities would require rapid contraction within the 23 Myr lifetime of the system and place the outer V1298 Tau planets at the upper density boundary of even mature field exoplanets, where few theories predict planets to exist. Since V1298 Tau e only transited once while observed with K2, SM21 did not have a period constraint from transits, and derived a period of 40.2 ± 0.9 days from their RV measurements. After the publication of SM21, Feinstein et al. (2022) published updated ephemerides of all four planets using TESS photometry, including a second transit of planet e. They placed a strict lower limit of 42 days on V1298 Tau e's period, in tension with the value from SM21. In addition, Tejada Arevalo et al. (2022) performed a dynamical stability analysis using the mass measurements reported in SM21, finding that 97% of the system configurations consistent with the SM21 posteriors are gravitationally unstable over the lifetime of the system.

SM21 performed a rigorous and state-of-the-art analysis, comparing several complex models with Bayesian methods. However, several independent lines of evidence appear to call the masses they report into question: the tension with formation theory discussed in SM21, the updated planet e orbital period of Feinstein et al. (2022), and the improbability of long-term stability derived by Tejada Arevalo et al. (2022).

The purpose of this paper is two-fold: (1) to demonstrate the use of CV to show that the preferred model of SM21 is overfitting the RVs, and (2) to point out several potential causes of the overfitting, which are not unique to SM21 but common throughout the literature. We do not attempt to update the mass estimates for the V1298 Tau planets in this paper. We also do not attempt to prove that the mass estimates published in SM21 are incorrect, but instead seek to call into question their reliability. Our argument is that future joint models of the stellar activity and planetary signals of V1298 Tau will need to prove their predictiveness in order to be trustworthy.

The structure of this paper is as follows: in Section 2, we review the literature data scrutinized in this paper and describe one additional contemporaneous RV data set taken with Keck/HIRES. In Section 3, we demonstrate that the preferred model of SM21 is overfit. Section 4 discusses several potential causes of this overfitting, and advises modelers on how to detect and/or avoid these subtle pitfalls. In particular, Section 4.4 argues that differential rotation is an important effect for V1298 Tau, and must be modeled carefully. We conclude in Section 5. We also provide an Appendix that provides a geometric interpretation of how GPR penalizes complexity. All of the code used to create the plots in this work are publicly available on GitHub, ²⁴ and a copy of these materials is deposited on Zenodo. ²⁵

We downloaded EVEREST-processed (Luger et al. 2016, 2018) K2 light curves for V1298 Tau using the lightkurve package (Lightkurve Collaboration et al. 2018). We used built-in lightkurve functions to remove NaNs, remove outliers, and normalize the data.

We obtained ground-based LCO photometry verbatim from SM21.

We obtained TESS light curves from Feinstein et al. (2022), who combined time series photometry of V1298 Tau from TESS Sectors 43 and 44. Feinstein et al. (2022) used the 2 minute light curve created by the Science Processing Operations Center pipeline (SPOC; Jenkins et al. 2016), and binned those observations to 10 minutes. We normalized the data for each TESS orbit separately, following Feinstein et al. (2022).

We obtained CARMENES and HARPS-N RVs directly from SM21. We note that SM21 excluded infrared-arm CARMENES RVs in their analysis, and we do the same here.

HARPS-N RVs are wavelength calibrated using a ThAr lamp, and the HARPS-N spectrograph covers 360–690 nm.

The visible arm of the CARMENES instrument covers the spectral range 520–960 nm, and spectra from this instrument are wavelength calibrated using a Fabry–Perot etalon, anchored using hollow cathode lamps.

Between 2018 November 16 and 2020 February 6, we obtained 36 RVs using the HIRES spectrograph on the Keck I telescope (Vogt et al. 1994). Wavelength calibration was performed by passing starlight through a warm iodine cell, and data reduction was performed using the California Planet Search pipeline described in Howard et al. (2010), which is adapted from Butler et al. (1996). All HIRES RVs used in this study are given in Table 1. Some of these RVs were previously published in Johnson et al. (2022), and the processing is identical in that paper and this. The same stellar template, constructed from two stellar spectra taken on 2019 October 24 UT without the iodine in the light path, was used to derive RVs in both studies. In-transit RVs from that study have been excluded here. Spectra were typically taken using the C2 decker (14'' × 0861), which enables sky subtraction, and is the CPS HIRES observer "decker of choice" for stars fainter than V ∼ 10. However, a CPS HIRES observer "rule of thumb" is to use the shorter B5 decker (35 × 0861) in poor seeing conditions, as the Doppler pipeline sky subtraction algorithm is unreliable when the stellar point-spread function fills the slit. Sky subtraction is not performed under such conditions. Accordingly, seven RVs published here were calculated from spectra using the B5 decker. In both modes, HIRES has a resolving power of ∼60,000, and the iodine cell spectral grasp translates to contributions to the RVs from wavelengths between 500 and 620 nm (Butler et al. 1996).

The mathematical formalism in this section is essentially identical to that of Cale et al. (2021; see their Section 3.2), but was developed independently. We encourage readers to compare our explanations, and we ask readers to also cite Cale et al. (2021) whenever referencing Section 4.1 of this paper.

There is a difference between correlated measurements that are allowed to have different GP amplitudes and data sets that share GP hyperparameters but are themselves uncorrelated. We are motivated to stress this distinction by the need in RV time series fitting to write down the joint likelihood of a model applied to data sets taken from several different instruments. As a concrete example, let us consider three fictional RV data points, the first two from HIRES and the next one from CARMENES, to which we would like to fit a GP model. Because of the different bandpasses of HIRES and CARMENES, we might expect the same stellar activity signal to have a different amplitude when observed by these two instruments. However, we might expect the time characteristics of the signals to be identical. In other words, we expect the CARMENES activity signal to be a scalar multiple of the HIRES activity signal. ²⁸ As discussed in Section 3, this assumption is borne out, at least to first order, in observations of other active stars at different wavelengths (see, e.g., Mahmud et al. 2011, who investigated the RV activity of the T-Tauri object Hubble I 4 with contemporaneous infrared and optical spectra taken with different instruments), but this point warrants further scrutiny. Comparing the variability of active stars with different instruments, as well as the variability of the Sun with different solar instruments, is an important endeavor.

Another important caveat is the use of different techniques for computing RVs from stellar spectra (e.g., the iodine/forward-modeling technique of HIRES versus the simultaneous reference/CCF technique of CARMENES and HARPS-N). Switching from one of these techniques to another is not expected to affect an astronomer's ability to recover common Keplerian signals, but spot activity is not a simple Doppler shift. More work is needed to understand and model spot activity at the spectral level. We proceed by assuming that modeling the same spectrum using an iodine/forward-model and with a simultaneous ThAr lamp reference (as an example) will only change the effective wavelength range of the spectrum that is used to compute the RV, and therefore affect only the amplitude of spot-induced variations.

Assuming linearly related GPs for different instruments, we can write down the joint covariance matrix for our three fictional data points, allowing unique amplitude terms a_C and a_H for each data set, and assuming an arbitrary kernel function k_i,j describing the covariance between RVs at times t_i and t_j :

$$\begin{eqnarray}{C}_{\mathrm{joint}}=\left(\begin{array}{ccc}{a}_{{\rm{H}}}^{2}{k}_{\mathrm{0,0}} & {a}_{{\rm{H}}}^{2}{k}_{\mathrm{0,1}} & {a}_{{\rm{H}}}{a}_{{\rm{C}}}{k}_{\mathrm{0,2}}\\ {a}_{{\rm{H}}}^{2}{k}_{\mathrm{1,0}} & {a}_{{\rm{H}}}^{2}{k}_{\mathrm{1,1}} & {a}_{{\rm{H}}}{a}_{{\rm{C}}}{k}_{\mathrm{1,2}}\\ {a}_{{\rm{C}}}{a}_{{\rm{H}}}{k}_{\mathrm{2,0}} & {a}_{{\rm{C}}}{a}_{{\rm{H}}}{k}_{\mathrm{2,1}} & {a}_{{\rm{C}}}^{2}{k}_{\mathrm{2,2}}\end{array}\right).\end{eqnarray} \tag{ 1 }$$

Optimizing the hyperparameters of a fit that uses this covariance matrix to define the GP likelihood will give the desired result.

SM21, following many other fits in the literature, constructed an independent covariance matrix for each RV instrument in their data set and summed the log(likelihoods) given by these together. This allows each RV data set to be independent; i.e., a data point taken by HIRES is not correlated with a data point taken at exactly the same time by CARMENES. Figures 11 and 12 illustrate the difference between these two likelihood definitions for data for a different object (chosen because it is easier to see the effect using this data set).

This assumption of independent data for each instrument effectively adds additional free parameters to a model, and makes it more susceptible to overfitting. This is also why, in Figures 2 and 3, we could demonstrate that the preferred SM21 model was overfitting by comparing the model prediction conditioned on HARPS-N data to the CARMENES data; the CARMENES data influenced the final values of the hyperparameters, since they were shared between the two GPs, but otherwise the data sets were treated as independent.

 To illustrate the effects discussed in this paper, we used a modified version of radvel (Fulton et al. 2016), built on tinygp (Foreman-Mackey et al. 2022), that treats the models for different instruments as correlated, but allows each instrument its own GP amplitude, white noise jitter term, and RV zero-point offset term. ²⁹ The difference between the previous version of radvel and this version is also illustrated in Figures 11 and 12 in the Appendix. ³⁰

Future work should continue to test this assumption by obtaining simultaneous (or near-simultaneous) RVs for a variety of stellar types with different instruments, across a wide range of bandpasses.

Another practice that may have made the SM21 preferred fit susceptible to overfitting involves constructing a GP kernel with one term at the rotation period and another term at its first harmonic. In other words, the SM21 preferred model kernel has the following form:

$$\begin{eqnarray}&&{C}_{{ij}}={f}_{1}({P}_{{\rm{rot}}})+{f}_{2}({P}_{{\rm{rot}}}/2).\end{eqnarray} \tag{ 2 }$$

To understand the motivation for this, we first need to scrutinize the RV signal in Fourier space. Figure 6 shows a Lomb–Scargle periodogram of all RV data presented in SM21, zooming in on two important parts of the period space. There are four extremely significant peaks in the RVs, which can all be explained with a single periodic signal at 2.91 days, the rotation period identified by SM21. Along with a strong peak at 2.91 days (hereafter P_rot), there is a signal at P_rot/2, which is often observed in RVs of stars showing starspot-induced variability (Nava et al. 2020). The other two strongly significant peaks can be explained as 1 day aliases of P_rot and P_rot/2. In other words, the dominant RV signal is periodic, but requires a two-component sinusoidal fit (i.e., it needs more terms in its Fourier expansion) in order for the fit to reproduce the shape of the curve. This is visualized in Figure 7, which shows the RVs phase folded to P_rot. In summary, the RV curve comprises a single periodic pattern, but that pattern is not a simple sinusoid.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The preferred SM21 model kernel sums two approximately quasi-periodic terms, one at P_rot and one at P_rot/2, because the approximate quasi-periodic kernel used in SM21 (Equation (1), derived in Foreman-Mackey et al. 2017) is less flexible than the standard quasi-periodic kernel (SM21, Equation (3)). In other words, the approximate kernel is less capable of fitting nonsinusoidal shapes. However, each term was modeled with its own independent exponential decay timescale. This adds an additional free parameter to the fit, which exacerbates the potential for overfitting.

The most straightforward way to address this is to construct a model with fewer unnecessary free parameters, for example by equating the parameters L₁ and L₂ in SM21's Equation (1). A more complicated suggestion, which would be an excellent avenue for further study, is to leverage the correlation between the photometry and RVs, following, for example, Rajpaul et al. (2015). This requires assuming (or fitting for) a relationship between a photometric data point and an RV data point at the same time. Our preliminary investigations along these lines indicate that the FF' formalism, which models an RV signal as a function of a simultaneous photometric (F) data set and the time derivative of the photometric data set (F'; Aigrain et al. 2012), does not allow for a good phenomenological match between the LCO photometry and the contemporaneous RVs, but the derivative of the LCO photometry appears to fit better (i.e., it appears possible to model the RV curve as a linear combination of the F' component only). ³¹ Future work could write down a joint GP formalism that models RVs as the time derivative of the photometry (such a formalism would be very similar to that of Rajpaul et al. 2015).

Regardless, in order to be confident in the relationship between the photometry and the RVs, as well as to pick out the components of the RVs that do not occur at P_rot, we suggest a very-high-cadence (several observations per night) RV follow-up campaign with contemporaneous photometry ³² in order to develop a high-fidelity model of the stellar variability. ³³ It is important to note that this campaign need not be performed by an RV instrument with 30 cm s⁻¹ precision; Johnson et al. (2022) demonstrated 6–7 m s⁻¹ rms precision with HIRES over several hours, even though the stars moves by hundreds of meters per second over even a single night. This level of instrumental RV error should be sufficient to understand the stellar activity, so long as the cadence is as high as possible.

A Keplerian signal has five free parameters (semiamplitude, eccentricity, argument of periastron, time of periastron, and period). A model with two Keplerian signals therefore has 10 additional free parameters than a model without. To first order, more free parameters means more model flexibility. This problem can be addressed using model comparison, which penalizes complexity. However, if there is unmodeled noise in the data, including additional Keplerian signals in the model can lead to overfitting; for example, high-eccentricity Keplerian models have similar properties to delta functions, which have relatively "flat" RV curves, except for a spike in RV near periastron. With insufficient sampling, outlier data points can be overfit with eccentric Keplerian signals.

A common worry in the RV modeling community is that using GPR to model stellar activity will "soak up" Keplerian signals, leading to underestimates of Keplerian RV semiamplitudes (discussed in Aigrain & Foreman-Mackey 2022), even when modeled jointly. However, we find evidence for the opposite effect in the SM21 preferred fit: that the Keplerian signals function as extra parameters that make the model susceptible to overfitting, and the GP is forced to compensate. Examining Figure 8, which shows the contributions to the mean model prediction from the Keplerians and the activity-only portion of the mean GP model, ³⁴ we find that the activity model interferes with the Keplerian model where RV data exist. This is seen most readily when smoothing the activity model over several rotation periods (effectively low-pass filtering the activity model).

Zoom In Zoom Out Reset image size

We can explain this behavior by imagining that there is some unmodeled noise source in the data that is inconsistent with Keplerian motion or quasi-periodic variability (see the next section). If some unphysical combination of parameters fits the data better at an epoch with many data points that is affected by this noise source, this may outweigh the negative Bayesian evidence contributions from (1) the added complexity and (2) the worse fit at epochs with fewer data points. We would then expect the Keplerian model to oversubtract at epochs with fewer data points (e.g., around JD = 1725 in Figure 8).

This effect suggests that the Keplerian signals in the SM21 preferred fit are not a viable description of the RV variability at timescales greater than the rotation period. More effort certainly needs to be spent understanding this phenomenon, but in the meantime we suggest performing CV tests in order to detect overfitting of this nature.

The previous subsections all argue that the preferred SM21 fit had too many free parameters (or effective free parameters) that allowed the model to overfit. In other words, we have argued that a simpler model (one for which the GP predictions for each instrument are scalar multiples of each other, a single period is present in the kernel, and no Keplerian signals are present in the model) would be more predictive, albeit perhaps with larger uncertainties. In this section, we suggest that this much simpler proposed model is still insufficient, because the host star has multiple, differentially rotating, active regions.

Differential rotation may not be the unmodeled noise source that we propose is affecting the SM21 preferred fit. The conclusions of this paper do not change if this is true. We discuss it here because it is potentially widely relevant, especially for young stars. We call for more work on modeling and understanding differential rotation in RVs.

4.4.1. Evidence for a Strong Differential Rotation Signal from Photometry

In the K2 and TESS photometry of V1298 Tau (Figure 9), two periodic signals of different amplitudes are visible by eye. These peaks are coherent in phase toward the end of both baselines, producing a larger overall photometric variability amplitude. Although each baseline covers only a portion of the beat periods implied by these different periods coming into and out of phase, the beating "envelope" is still easily distinguished. To guide the eye, we overplotted the shape of the beating envelope formed by the three dominant periods in the Lomb–Scargle periodogram of the K2 data.

Zoom In Zoom Out Reset image size

Multiple closely related periodicities are also apparent in the periodograms of the K2 and TESS data (and the LCO data, albeit at lower significance, potentially due to the lower cadence of that data set; Figure 10). In particular, over both the K2 and TESS baselines, dominant periodicities at 2.85 and 2.92 days, respectively, and two less prominent periodicities (one at a larger period, and one at a smaller period) are present. The multiple periodicities in the light curve, visible both in the shape of the beating envelope and in Fourier space, have often been interpreted as smoking gun evidence of differential rotation (see, e.g., Lanza et al. 1994; Frasca et al. 2011). It is important to note, however, that short spot lifetimes may also produce the observed photometric pattern, and have been shown in simulations to be easily confused with differential rotation (see, e.g., Basri & Shah 2020). Longer photometric temporal baselines than are available in the photometric data presented in this paper are needed to distinguish between the two. The conclusion of this section (that there is a noise source visible in photometry that is unmodeled in the SM21 preferred model) would remain unchanged in this case, but this interpretation has important implications for future modeling efforts. That the signals arise from a close binary is ruled out by the multiple nearby periods in the light curve (rotation of two tidally extended binary stars can produce a similar pattern, but with a single period), while astroseismic pulsations are ruled out by the amplitude and period of the variability; V1298 Tau is a 1.2 M_⊙ PMS star with $\mathrm{log}g$ = 4.48 (SM21), which we would expect to be oscillating on scales of minutes and ≲1 ppt, not days and 20 ppt (Chaplin & Miglio 2013; see their Figure 3).

Zoom In Zoom Out Reset image size

4.4.2. Effect on RVs