The input and output columns are selected by the user in a first phase. Only then can modelization be performed.

Only numeric rows that pass the active search filter are used.

Cross-validation splits the data into \( k = 5\) consecutive blocks (folds) based on the current row order. If your data has a time structure (daily measurements, for instance), the folds follow time: folds 1–4 are trained and fold 5 is validated on the last 20
If the rows are sorted by one of the inputs, the folds will have systematic differences in input range, which can inflate or deflate CV error estimates. For most tabular data with no natural ordering, row order is arbitrary and the results are stable.

Nonlinear regression is implemented by mapping each standardized input through a B-spline basis, converting a nonlinear problem into ordinary linear regression in the transformed space. B-splines have two key advantages over polynomial bases:

• Compact support: each basis function is non-zero only over a short interval, reducing collinearity.
• Numerical stability: the design matrix is much better conditioned than a polynomial Vandermonde matrix.

A 1-D spline space is defined by a knot vector and a polynomial degree \( d\) . The tool uses \( d = 3\) (cubic splines) throughout.

K — the mesh parameter is the number of interior knots per spline term. It is user-configurable (default 4, range 2–20) and controls how finely the input domain is divided. For \( K\) interior knots \( \xi_1 < \xi_2 < \dots < \xi_K\) placed at quantiles of the training data, the full extended knot vector has the boundary values repeated \( d+1\) times:

This gives a spline space of dimension

For the default \( K = 4\) , each input has \( B = 8\) basis functions. For the refined basis (M3), one midpoint per span is inserted, giving \( K' = 2K - 1\) interior knots and \( B' = 2K + 3\) basis functions (e.g. \( K=4 \to K'=7, B'=11\) ; \( K=8 \to K'=15, B'=19\) ).

Basis function \( N_{i,d}(x)\) is defined recursively:

with the convention \( 0/0 = 0\) . The vector \( \phi(x) = (N_{1,d}(x), \dots, N_{B,d}(x))\) is the feature vector for scalar input \( x\) .

An important property is the partition of unity:

This ensures the intercept is not multicollinear with the basis functions.

Interior knots are placed at empirical quantiles of the training data. For \( K\) knots, they are placed at the \( \frac{i}{K+1} \times 100\) percentiles for \( i = 1, \dots, K\) , distributing basis support according to data density. For the default \( K = 4\) this gives the 20th, 40th, 60th, and 80th percentiles. The boundary knots are set at the minimum and maximum training values; inputs outside this range are clamped at prediction time.

Knot placement is recomputed on the training fold only inside each cross-validation fold, preventing data leakage from validation rows into the spline mesh.

4.5.1 What K controls

4.5.2 Effect on the model ladder

4.5.3 The P/N constraint

4.5.4 Reading the diagnostics

4.5.5 When to increase K

4.5.6 When to keep K small

For an additive model with \( m\) inputs, each input \( r \in \{1,\dots,m\}\) gets its own \( B_r\) -dimensional B-spline basis. The design matrix for \( N\) training rows is

where \( \Phi_r \in \mathbb{R}^{N \times B_r}\) is the spline design block for variable \( r\) , and the leading column of 1s is the intercept. The total parameter count per output is

Coefficients \( C \in \mathbb{R}^{P \times n}\) (one column per output) are found by minimizing the penalized residual sum of squares:

where \( \lambda > 0\) is the regularization strength and \( \Omega = \mathrm{diag}(0, 1, \dots, 1)\) leaves the intercept unpenalized and equally penalizes all basis coefficients. The closed-form solution satisfies the normal equations:

This is solved via SVD for numerical robustness.

The regularization parameter is selected from a grid

by 5-fold cross-validation. The \( \lambda\) minimizing mean CV RMSE is used for the final refit on all data.

The \( N\) rows are split into \( k = 5\) folds of equal size (positional: first \( N/5\) rows → fold 1, etc.). For each fold \( f\) :

1. Rows not in fold \( f\) form the training set; rows in fold \( f\) form the validation set.
2. Inputs and outputs are standardized on the training rows; the same mean/std is applied to validation rows.
3. Knot quantiles are computed from training inputs only (no leakage).
4. For PCA models (Branch B), PCA loadings are computed from training inputs only.
5. The model is fitted on the training design matrix at each candidate \( \lambda\) .
6. Validation RMSE (in standardized output space) is recorded.

Mean and standard deviation of CV RMSE across the 5 folds are returned for each \( \lambda\) .

All RMSE values are reported in standardized output space (each output has mean 0, std 1 on training data). This makes scores directly comparable across different outputs:

A CV RMSE near 1.0 means the inputs carry little information about the outputs (or the relationship is dominated by variation not captured by the inputs).

\( P = 1\) . The model predicts the training mean of each output. CV RMSE ≈ 1.0 by construction (predicting the mean on held-out data in standardized space gives RMSE ≈ 1). Serves as the baseline.

\( P = 1 + m\) . Captures linear relationships between inputs and outputs. The design matrix has one column per input plus the intercept.

with \( B = K + 4\) basis functions per input (\( K\) interior knots, degree 3). \( P = 1 + (K+4)m\) . For the default \( K = 4\) this gives \( B = 8\) and \( P = 1 + 8m\) . Captures smooth nonlinear main effects but no interaction between inputs. Each input is modeled independently by a cubic spline.

Same structure as M2 but the knot vector of each term is refined by inserting one midpoint per span, giving \( K' = 2K - 1\) interior knots and \( B' = 2K + 3\) basis functions per input.

\( P = 1 + (2K+3)m\) . For the default \( K = 4\) this gives \( K' = 7\) , \( B' = 11\) , and \( P = 1 + 11m\) .

The refined knots are a strict superset of the coarse knots (nested spaces: \( \mathcal{F}_{M2} \subseteq \mathcal{F}_{M3}\) ).

Extends M3 by adding a tensor-product spline block for one pair of inputs \( (r, s)\) :

where \( \otimes\) denotes the outer product of the two \( K\) -dimensional coarse feature vectors, giving \( K^2\) additional columns. For the default \( K = 4\) this is 16 columns.

\( P_{M4} = P_{M3} + K^2\) .

The pair \( (r, s)\) is selected greedily: all candidate pairs are evaluated on M3 residuals with a fixed \( \lambda = 1\) , and the pair giving the largest residual reduction is chosen.

Adds a second tensor-product block for a second greedily-selected pair \( (r', s') \ne (r, s)\) .

\( P_{M5} = P_{M4} + K^2\) .

For \( m = 1\) , the ladder stops at M3 (no pairs exist). M5 requires \( m \geq 3\) .

Parameter count formulas (m inputs, K interior knots):

Increasing K from 4 to 8 raises M2 from 41 to 61 parameters (for m=5), and causes interaction blocks to grow from 16 to 64 columns — always check that \( P \lesssim N/5\) before running with a large K.

Given standardized inputs \( \tilde X \in \mathbb{R}^{N \times m}\) , the thin SVD gives

where the columns of \( V \in \mathbb{R}^{m \times m}\) are the PC directions (sorted by descending singular value). The first \( q\) PC scores are

The fraction of variance explained by the first \( q\) PCs is \( \sum_{i=1}^q \sigma_i^2 / \sum_i \sigma_i^2\) .

PCA is refit inside each cross-validation fold on the training rows only, then applied to validation rows. This avoids leakage of global variance structure into the validation estimates.

For each number of PCs \( q \in \{1, 2, 3\}\) (up to \( \min(m,3)\) ), the following models are evaluated:

Interaction entries appear only for \( q \geq 2\) . Interaction pairs are fixed as \( (\text{PC}_0, \text{PC}_1)\) , \( (\text{PC}_0, \text{PC}_2)\) , … — always anchored to the highest-variance PCs.

Total Branch B entries: \( q_{\max}(q_{\max}+3)/2\) (5 for \( q_{\max}=2\) , 9 for \( q_{\max}=3\) ).

Branch B is useful when:

• the inputs are correlated and a few principal components capture most of the variance,
• you want to know whether the outputs are primarily driven by a global contrast (PC1) vs. finer distinctions (PC2, PC3),
• you suspect interactions but only in the dominant variance directions.

It will generally perform similarly to Branch A when the inputs are already near-independent, but may outperform Branch A when there is strong multicollinearity.

After running rungs, the recommended model is the simplest rung satisfying both:

1. \( \text{cvRmse} \leq \text{cvRmse}_{\text{best}} \times 1.02\) (within 2
   
2. \( \kappa < 10^8\) (conditioning below the disqualification limit).

This deliberately favors simpler models: if M1 and M3 give nearly the same CV RMSE, M1 is chosen because it has fewer parameters and is more likely to generalize.

The scatter plot displays every run rung as a dot with coordinates:

• x-axis: CV RMSE (accuracy, lower is better),
• y-axis: \( \log_{10} \kappa\) (ill-conditioning, lower is better).

The ideal corner is top-left: perfect accuracy (\( \text{cvRmse} \to 0\) ) with perfect conditioning (\( \log_{10}\kappa \to 0\) ).

Two reference lines mark practical limits:

• Dashed red — the no-skill baseline at \( \text{cvRmse} = 1.0\) . Models to the right of this line predict worse than the training mean.
• Dashed orange — the conditioning limit at \( \log_{10}\kappa = 8\) (\( \kappa = 10^8\) ). Models above this line are numerically unreliable.

The axes are anchored at these absolute thresholds (not auto-scaled), so the chart is visually stable as you add more rungs.

To rank rungs, each point is mapped to a normalized distance from the ideal origin:

This combines accuracy and robustness on a common scale. Both thresholds are physically meaningful:

• dividing by 1.0 puts no-skill performance at distance 1 on the x-axis,
• dividing by 8.0 puts the conditioning limit at distance 1 on the y-axis.

The rung with the smallest \( d\) is ranked first.

Climbing the ladder increases the nominal parameter count \( P\) , but the ridge regularizer prevents all \( P\) parameters from being active. The effective degrees of freedom (EDF, shown in the ranking table) is the actual complexity the fit exercises:

At \( \lambda = 0\) the model uses all \( P\) parameters (EDF = P); as \( \lambda \to \infty\) the model collapses to the intercept (EDF = 1).

What low EDF means in practice:

• EDF = 1 — the model predicts a constant; all basis coefficients have been regularized to zero. The fit is maximally robust but carries no information.
• EDF ≈ 2–3 — roughly one or two input directions carry a non-trivial signal; the rest is suppressed.
• EDF close to P — little regularization was needed; the data supports all the basis functions.

The tradeoff in a nutshell. Adding more basis functions (going from M1 to M2, M3, …) gives the model capacity to represent complex shapes, but:

• if the signal is weak, the regularizer compensates by shrinking most new coefficients to near-zero, and EDF barely increases over the simpler model;
• the design matrix gains columns, raising \( \kappa\) and moving the scatter-plot dot upward;
• CV RMSE may not improve, or may even worsen, because the extra flexibility is fitting noise.

The scatter plot makes this visible: climbing the ladder on data with a weak signal moves points right and up — worse accuracy and worse conditioning — without reducing EDF meaningfully. The right response is not to add more model complexity, but to reconsider the inputs (see §12.6).

A CV RMSE near 1.0 across all models means the inputs have little predictive power over the outputs given the current data. This can happen because:

• Confounding seasonal cycles: if both inputs and outputs follow an annual cycle (e.g., temperature and evaporation), the direct correlation between them at the daily scale is masked. The model explains variance due to seasonality, not the physical relationship you want to study.

Solution: deseasonalize both inputs and outputs by subtracting a fitted seasonal mean (e.g., monthly or weekly averages) before running the modelization. The residuals will then expose the within-cycle correlation.

• Missing drivers: the true cause of output variation is not in the input columns. Adding more relevant inputs can help.

• Too short a series: with few data points (\( N < 5P\) ) the cross-validation estimate is noisy and the model cannot generalize.

A high \( \kappa\) for a given model class usually means either:

• the inputs are highly collinear (correlations near ±1), or
• the model has too many parameters relative to the data range.

Increasing \( \lambda\) (which the λ grid search handles automatically) will reduce \( \kappa\) by adding a stabilizing diagonal to \( G_\lambda\) . If even the maximum \( \lambda = 100\) leaves \( \kappa > 10^8\) , the model class is too complex for the available data — try a simpler rung.

For time-series data (e.g., daily measurements), keep the rows in chronological order. The 5-fold split will then validate the last 20

• If \( m \leq 2\) , both branches fit roughly the same models; Branch A is simpler to interpret.
• If \( m \geq 3\) and the inputs are correlated, try Branch B first to see whether 1–2 PCs already capture most of the relevant signal.
• Compare the best Branch A and Branch B rungs on the scatter plot: the rung closest to the top-left corner is preferred.

A model can be accurate without being meaningful, and meaningful without being accurate. The CV RMSE measures only predictive accuracy; it says nothing about whether the learned relationship reflects a genuine causal or physical link between inputs and outputs.

Confounding is the most common source of misleading accuracy. If both an input and an output are driven by a third variable that is not in the model, the fit will capture the shared variation due to that third variable rather than the direct relationship. The model will appear to predict well in-sample but will fail whenever the confounding variable behaves differently.

Conversely, a genuine physical relationship between two variables may be weak in the data if both are dominated by a shared confounding cycle. In that case CV RMSE stays near 1.0 across all model classes — the scatter plot shows all dots near the right edge — and EDF collapses to 1–2 regardless of which rung you run. This is not a failure of the fitting algorithm; it is a correct signal that the inputs, as presented, carry little information about the outputs beyond what the confounder already explains.

Recognizing confounding:

• All model classes give CV RMSE ≈ 1.0 and EDF ≈ 1–2 despite having many input columns.
• The inputs and outputs share a visible trend or periodic structure (yearly cycle, business cycle, growth trend, etc.).
• Removing or controlling for that trend substantially changes the apparent correlation between inputs and outputs.

Remedies:

1. Add the confounder as an explicit input. If the shared trend can be represented by a column already in the table (e.g., a date index, a calendar variable, a known covariate), include it. The spline model will then explain the shared trend through that column, and the remaining inputs will be estimated on the residuals.

1. Detrend or deseasonalize before modeling. Subtract a fitted trend or periodic mean from both inputs and outputs, then run the modelization on the residuals. This forces the model to explain within-cycle variation only.

1. Restrict the data range. Use the DataTable filter to select a sub-period or sub-group in which the confounding trend is approximately constant. The model fitted on the filtered rows estimates the relationship conditional on that range.

In all cases, the goal is to present the model with inputs and outputs whose shared variation is primarily due to the relationship you wish to study, not to a third variable driving both.

If the DataTable search filter is active, only the visible rows are used for fitting. This can be used deliberately (e.g., restrict fitting to a season or experimental condition), but the resulting model is only valid for inputs in that filtered range.