# Model fitting Whittaker fits GAMs by **penalized iteratively reweighted least squares (P-IRLS)**, a generalization of weighted least squares that handles non-Gaussian responses by iterating over a sequence of penalized weighted regression problems. This page describes the fitting algorithm, smoothness selection criteria, and the practical controls you have over the fitting process. # The P-IRLS algorithm For a given set of smoothing parameters \boldsymbol{\lambda} = (\lambda_1, \ldots, \lambda_J), fitting proceeds in five steps: 1. **Initialize**: set \$= \$ `family.initialize(y)`, \eta = g(\mu). The initialization is family-specific: for Gaussian data \mu = y, for Poisson \mu = y + 0.1, for binomial \mu = (y + 0.5)/2. 2. **Pseudo-data**: compute the working response z = \eta + (y - \mu)\\ g'(\mu) This is a first-order Taylor expansion of the link function around the current \mu, converting the non-Gaussian problem into a weighted least-squares problem. 3. **Working weights**: compute W = \frac{1}{\bigl\[g'(\mu)\bigr\]^2 \\ V(\mu)} where V(\mu) is the variance function of the family. These weights account for both the link function and the heteroscedasticity of the response. 4. **Penalized WLS**: solve the penalized weighted least-squares system \bigl(X^\top W X + \textstyle\sum\_{j=1}^{J} \lambda_j S_j\bigr)\\\hat\beta = X^\top W z where S_j is the penalty matrix for the j-th smooth term. This is the core computational step. The penalty matrices S_j are positive semi-definite and penalize the roughness of each smooth, pulling f_j toward the penalty null space (typically low-order polynomials). 5. **Update**: compute \eta = X\hat\beta and \mu = g^{-1}(\eta). Repeat steps 2-5 until the relative change in penalized deviance falls below 10^{-7}. ## The Gaussian identity-link special case For the Gaussian family with an identity link, P-IRLS reduces to a **single** penalized least-squares solve. The link function is the identity (g(\mu) = \mu), so g'(\mu) = 1, and the variance function is constant (V(\mu) = 1). The working weights become W = I and the pseudo-data become z = y, so the system collapses to: \bigl(X^\top X + \textstyle\sum_j \lambda_j S_j\bigr)\\\hat\beta = X^\top y No iteration is required. This makes Gaussian GAMs substantially faster to fit than their non-Gaussian counterparts. # Smoothness selection The roughness of each smooth is controlled by a non-negative smoothing parameter \lambda_j. Larger \lambda_j produces a smoother f_j while \lambda_j = 0 removes the penalty entirely (interpolating spline). Rather than setting \lambda_j by hand, Whittaker selects it automatically by optimizing one of three criteria. ## GCV (generalized cross-validation) GCV selects \boldsymbol\lambda by minimizing: \text{GCV}(\boldsymbol\lambda) = \frac{n\\ D(y,\hat\mu)}{\[n - \text{tr}(\mathbf{A})\]^2} where D is the deviance and \mathbf{A} = X(X^\top W X + \sum_j \lambda_j S_j)^{-1}X^\top W is the influence (hat) matrix. The trace \text{tr}(\mathbf{A}) is the effective degrees of freedom of the model. GCV is an approximation to leave-one-out cross-validation that avoids refitting the model n times. It is minimized over \log\boldsymbol\lambda using a Newton method with analytical gradient and Hessian. **Properties**: - GCV is well-understood theoretically and widely used. - It tends to undersmooth in finite samples -- occasional fits with too many effective degrees of freedom. - It can be sensitive to influential observations because it weights all residuals equally. ``` python model = wk.GAM("y ~ s(x)") model.fit(data, method="GCV") ``` ## REML (restricted maximum likelihood) REML treats the smooth coefficients as random effects and maximizes the restricted log-likelihood: \ell\_{\text{REML}}(\boldsymbol\lambda, \phi) = -\frac{1}{2}\Bigl\[ n \log(2\pi\phi) + \frac{D}{\phi} + \log\bigl\|X^\top W X + \textstyle\sum_j \lambda_j S_j\bigr\| - \log\bigl\|\textstyle\sum_j \lambda_j S_j^{+}\bigr\| \Bigr\] where S_j^{+} denotes the pseudoinverse restricted to the range space of S_j, and \phi is the scale parameter. **Advantages over GCV**: - REML has better finite-sample properties: it is less prone to undersmoothing. - It is more stable when the true function is close to the penalty null space. - REML accounts for the uncertainty in \hat\beta when estimating \phi, producing better variance estimates. - It is the recommended default in both Whittaker and `mgcv`. ``` python model = wk.GAM("y ~ s(x)") model.fit(data, method="REML") # recommended default ``` ## ML (maximum likelihood) ML maximizes the full marginal log-likelihood rather than the restricted version: \ell\_{\text{ML}}(\boldsymbol\lambda, \phi) = -\frac{1}{2}\Bigl\[ n \log(2\pi\phi) + \frac{D}{\phi} + \log\bigl\|X^\top W X + \textstyle\sum_j \lambda_j S_j\bigr\| \Bigr\] The key difference from REML is the absence of the \log\|\sum_j \lambda_j S_j^{+}\| term. This means ML does not adjust for the degrees of freedom consumed by the fixed effects, which can lead to slight undersmoothing relative to REML -- the same bias that makes ML variance estimates biased downward in classical linear models. **When to use ML**: ML is appropriate when you need to compare models with different fixed-effect structures using likelihood ratio tests, because REML likelihoods are not comparable across models with different fixed effects. ``` python model = wk.GAM("y ~ s(x)") model.fit(data, method="ML") ``` ## Comparing GCV, REML, and ML ``` python import numpy as np import whittaker as wk # Generate noisy data with a smooth underlying function rng = np.random.default_rng(23) n = 150 x = np.linspace(0, 2 * np.pi, n) y_true = np.sin(x) + 0.3 * np.cos(3 * x) y = y_true + rng.normal(0, 0.4, n) data = {"x": x, "y": y} # Fit with each method model_gcv = wk.GAM("y ~ s(x, k=20)") model_gcv.fit(data, method="GCV") model_reml = wk.GAM("y ~ s(x, k=20)") model_reml.fit(data, method="REML") model_ml = wk.GAM("y ~ s(x, k=20)") model_ml.fit(data, method="ML") # Compare EDF and smoothing parameters print("Method | EDF total | Smoothing param") print("--------|-----------|----------------") print(f"GCV | {model_gcv.edf_total:9.3f} | {model_gcv.smoothing_params[0]:.4f}") print(f"REML | {model_reml.edf_total:9.3f} | {model_reml.smoothing_params[0]:.4f}") print(f"ML | {model_ml.edf_total:9.3f} | {model_ml.smoothing_params[0]:.4f}") ``` Method | EDF total | Smoothing param --------|-----------|---------------- GCV | 11.054 | 0.6030 REML | 9.853 | 1.0530 ML | 9.853 | 1.0530 GCV typically selects the smallest smoothing parameter (largest EDF), REML the largest (fewest EDF), and ML falls between the two. ``` python import altair as alt # Build predictions for each method x_plot = np.linspace(0, 2 * np.pi, 200) preds_gcv = model_gcv.predict({"x": x_plot}) preds_reml = model_reml.predict({"x": x_plot}) preds_ml = model_ml.predict({"x": x_plot}) # Assemble plot data plot_records = [] for i in range(len(x_plot)): plot_records.append({"x": float(x_plot[i]), "fit": float(preds_gcv.values[i]), "Method": "GCV"}) plot_records.append({"x": float(x_plot[i]), "fit": float(preds_reml.values[i]), "Method": "REML"}) plot_records.append({"x": float(x_plot[i]), "fit": float(preds_ml.values[i]), "Method": "ML"}) # Observed data obs_records = [{"x": float(x[i]), "y": float(y[i])} for i in range(n)] # True function (recomputed on the plot grid) y_true_plot = np.sin(x_plot) + 0.3 * np.cos(3 * x_plot) true_records = [{"x": float(x_plot[i]), "y": float(y_true_plot[i])} for i in range(len(x_plot))] points = alt.Chart({"values": obs_records}).mark_circle( size=12, opacity=0.25, color="gray" ).encode( x=alt.X("x:Q", title="x"), y=alt.Y("y:Q", title="y"), ) lines = alt.Chart({"values": plot_records}).mark_line(strokeWidth=2).encode( x="x:Q", y="fit:Q", color=alt.Color("Method:N", scale=alt.Scale( domain=["GCV", "REML", "ML"], range=["#e45756", "#4c78a8", "#72b7b2"] )), strokeDash=alt.StrokeDash("Method:N", scale=alt.Scale( domain=["GCV", "REML", "ML"], range=[[0], [0], [4, 4]] )), ) true_line = alt.Chart({"values": true_records}).mark_line( color="black", strokeDash=[2, 2], strokeWidth=1, opacity=0.5 ).encode(x="x:Q", y="y:Q") (points + true_line + lines).properties( width="container", height=300, title="Smoothness selection: GCV vs REML vs ML" ) ``` The black dashed line is the true function. REML (blue) tends to produce the smoothest fit, GCV (red) the most flexible, and ML (teal, dashed) sits in between. > **Tip: When in doubt, use REML** > > REML is the recommended default for nearly all applications. It is less prone to overfitting than GCV and produces more reliable confidence intervals. Switch to GCV only if you have a specific reason (e.g., reproducing a legacy analysis) or to ML when you need comparable likelihoods across models with different fixed effects. # Fixed smoothing parameters Sometimes you want to bypass automatic smoothness selection and fix \lambda_j at known values. This is useful for: - **Reproducing results** from another analysis where \lambda was determined externally. - **Sensitivity analysis**: checking how the fit changes as you vary \lambda. - **Simulation studies**: fitting at a known truth. Pass a list of smoothing parameters to `fit()` in the same order as the smooth terms appear in the formula: ``` python # Fix the smoothing parameter for s(x) model_fixed = wk.GAM("y ~ s(x, k=20)") model_fixed.fit(data, smoothing_params=[1.0]) print(f"Fixed lambda: {model_fixed.smoothing_params}") print(f"EDF: {model_fixed.edf_total:.3f}") ``` Fixed lambda: [1.0] EDF: 9.960 When `smoothing_params` is provided, the [method](../reference/CausalGAM.md#whittaker.CausalGAM.method) argument is ignored -- no optimization over \lambda takes place. ``` python # Multiple smooths: one lambda per smooth rng = np.random.default_rng(0) n = 200 x1 = np.linspace(0, 1, n) x2 = rng.uniform(0, 1, n) y2 = np.sin(4 * x1) + 0.5 * x2**2 + rng.normal(0, 0.3, n) data2 = {"x1": x1, "x2": x2, "y": y2} model_two = wk.GAM("y ~ s(x1) + s(x2)") model_two.fit(data2, smoothing_params=[1.0, 0.5]) print(f"Lambda for s(x1): {model_two.smoothing_params[0]:.2f}") print(f"Lambda for s(x2): {model_two.smoothing_params[1]:.2f}") ``` Lambda for s(x1): 1.00 Lambda for s(x2): 0.50 > **Warning: Warning** > > Fixing smoothing parameters at inappropriate values can produce severely under- or overfit models. Use automatic selection (REML) unless you have a compelling reason to fix \lambda. # Variable selection with `select=True` Standard penalized regression cannot shrink a smooth term to zero -- the penalty null space (the space of functions not penalized, usually polynomials up to a given order) always survives. The **double penalty** approach adds an extra penalty on the null space, allowing the entire smooth to be penalized toward zero: \text{penalty}\_j = \lambda_j \boldsymbol\beta^\top S_j \boldsymbol\beta + \lambda_j^{\*} \boldsymbol\beta^\top S_j^{\*} \boldsymbol\beta where S_j^{\*} penalizes the null space of S_j and \lambda_j^{\*} is an additional smoothing parameter estimated alongside \lambda_j. When `select=True`, Whittaker adds this extra penalty for every smooth term. If a predictor has no effect, both \lambda_j and \lambda_j^{\*} can grow large enough to drive the smooth's EDF to approximately zero, effectively removing it from the model. ``` python # Variable selection example: x3 is pure noise rng = np.random.default_rng(23) n = 300 x1 = np.linspace(0, 2 * np.pi, n) x2 = rng.uniform(0, 1, n) x3 = rng.normal(0, 1, n) # noise predictor y_sel = np.sin(x1) + x2**2 + rng.normal(0, 0.3, n) data_sel = {"x1": x1, "x2": x2, "x3": x3, "y": y_sel} # Without variable selection model_nosel = wk.GAM("y ~ s(x1) + s(x2) + s(x3)") model_nosel.fit(data_sel, method="REML") # With variable selection model_sel = wk.GAM("y ~ s(x1) + s(x2) + s(x3)") model_sel.fit(data_sel, method="REML", select=True) print("Without select=True:") print(model_nosel.summary()) print("\nWith select=True:") print(model_sel.summary()) ``` Without select=True: GAM fit summary ============================================================ Formula: y ~ s(x1) + s(x2) + s(x3) Family: Gaussian(link='identity') Inference: REML Observations: 300 Coefficients: 28 Parametric coefficients: Term Estimate Std.Err t value p-value ------------------------ ---------- ---------- ---------- ---------- (Intercept) 0.3331 0.0165 20.128 < 1e-16 Approximate significance of smooth terms: Term EDF Ref.df Chi.sq p-value ------------------------ ------ ------ ---------- ---------- s(x1) 7.57 8 1924.348 < 1e-16 s(x2) 2.52 3 253.698 < 1e-16 s(x3) 1.90 2 2.021 0.3641 Total EDF: 12.99 Scale est: 0.082157 Deviance: 23.5797 Null dev: 201.2325 Dev. expl: 88.3% GCV score: 0.085877 AIC: 114.62 BIC: 162.74 With select=True: GAM fit summary ============================================================ Formula: y ~ s(x1) + s(x2) + s(x3) Family: Gaussian(link='identity') Inference: REML Observations: 300 Coefficients: 28 Parametric coefficients: Term Estimate Std.Err t value p-value ------------------------ ---------- ---------- ---------- ---------- (Intercept) 0.3331 0.0166 20.086 < 1e-16 Approximate significance of smooth terms: Term EDF Ref.df Chi.sq p-value ------------------------ ------ ------ ---------- ---------- s(x1) 7.55 8 1914.498 < 1e-16 s(x2) 2.46 3 250.370 < 1e-16 s(x3) 0.00 1 0.000 0.9969 Total EDF: 11.01 Scale est: 0.082505 Deviance: 23.8434 Null dev: 201.2325 Dev. expl: 88.2% GCV score: 0.085647 AIC: 113.90 BIC: 154.67 ``` python # Grouped bar chart comparing per-term EDF with and without select=True term_labels = ["s(x1)", "s(x2)", "s(x3)"] edf_nosel = model_nosel.edf edf_sel = model_sel.edf bar_records = [] for label, edf_std, edf_shrink in zip(term_labels, edf_nosel, edf_sel): bar_records.append({"Term": label, "Model": "Standard", "EDF": float(edf_std)}) bar_records.append({"Term": label, "Model": "select=True", "EDF": float(edf_shrink)}) alt.Chart({"values": bar_records}).mark_bar().encode( x=alt.X("Term:N", title="Smooth term", axis=alt.Axis(labelAngle=0)), y=alt.Y("EDF:Q", title="Effective degrees of freedom"), color=alt.Color("Model:N", scale=alt.Scale( domain=["Standard", "select=True"], range=["#4c78a8", "#e45756"] )), xOffset="Model:N", ).properties(width="container", height=300, title="Variable selection: per-term EDF") ``` With `select=True`, the EDF for `s(x3)` should be driven close to zero, confirming that Whittaker identified the noise predictor. The EDF values for `s(x1)` and `s(x2)` remain largely unchanged. > **Tip: When to use variable selection** > > Use `select=True` when you have many candidate predictors and want the model to automatically determine which ones have a nonlinear effect. It is especially useful in exploratory analysis where you are unsure which predictors are relevant. For confirmatory analysis where the model structure is known, standard fitting (without `select=True`) is appropriate. # Observation weights Observation weights enter the penalized WLS system by modifying the weight matrix. If you supply a weight vector w_i, the working weight matrix becomes \tilde{W} = \text{diag}(w_i) \cdot W, and the penalized WLS system is: \bigl(X^\top \tilde{W} X + \textstyle\sum_j \lambda_j S_j\bigr)\\\hat\beta = X^\top \tilde{W} z This allows you to: - **Downweight outliers**: set w_i \< 1 for observations you suspect are contaminated. - **Account for known precision**: when observations have different known variances, set w_i = 1 / \sigma_i^2. - **Handle aggregated data**: when each row represents n_i observations, set w_i = n_i. ``` python # Weights example: downweight outlier observations rng = np.random.default_rng(23) n = 100 x = np.linspace(0, 2 * np.pi, n) y_w = np.sin(x) + rng.normal(0, 0.3, n) # Add some outliers outlier_idx = [10, 30, 50, 70, 90] y_w[outlier_idx] += rng.choice([-3, 3], size=len(outlier_idx)) # Create weights: 1.0 for normal observations, 0.1 for outliers w = np.ones(n) w[outlier_idx] = 0.1 data_w = {"x": x, "y": y_w} # Fit without weights model_unweighted = wk.GAM("y ~ s(x)") model_unweighted.fit(data_w, method="REML") # Fit with weights model_weighted = wk.GAM("y ~ s(x)") model_weighted.fit(data_w, method="REML", weights=w) print(f"Unweighted scale: {model_unweighted.scale:.4f}") print(f"Weighted scale: {model_weighted.scale:.4f}") print(f"Unweighted EDF: {model_unweighted.edf_total:.3f}") print(f"Weighted EDF: {model_weighted.edf_total:.3f}") ``` Unweighted scale: 0.5376 Weighted scale: 0.1516 Unweighted EDF: 5.587 Weighted EDF: 7.000 ``` python # Scatter with outliers highlighted plus weighted vs unweighted fitted curves x_fit = np.linspace(0, 2 * np.pi, 200) preds_unw = model_unweighted.predict({"x": x_fit}) preds_w = model_weighted.predict({"x": x_fit}) is_outlier = np.zeros(n, dtype=bool) is_outlier[outlier_idx] = True scatter_records = [ {"x": float(x[i]), "y": float(y_w[i]), "Type": "Outlier" if is_outlier[i] else "Normal"} for i in range(n) ] fit_records = [] for i in range(len(x_fit)): fit_records.append({"x": float(x_fit[i]), "fit": float(preds_unw.values[i]), "Model": "Unweighted"}) fit_records.append({"x": float(x_fit[i]), "fit": float(preds_w.values[i]), "Model": "Weighted"}) points = alt.Chart({"values": scatter_records}).mark_circle().encode( x=alt.X("x:Q", title="x"), y=alt.Y("y:Q", title="y"), color=alt.Color("Type:N", scale=alt.Scale( domain=["Normal", "Outlier"], range=["gray", "#e45756"] )), size=alt.Size("Type:N", scale=alt.Scale( domain=["Normal", "Outlier"], range=[30, 80] ), legend=None), opacity=alt.condition( alt.datum.Type == "Outlier", alt.value(0.8), alt.value(0.3), ), ) lines = alt.Chart({"values": fit_records}).mark_line(strokeWidth=2).encode( x="x:Q", y="fit:Q", color=alt.Color("Model:N", scale=alt.Scale( domain=["Unweighted", "Weighted"], range=["#e45756", "#4c78a8"] )), strokeDash=alt.StrokeDash("Model:N", scale=alt.Scale( domain=["Unweighted", "Weighted"], range=[[4, 4], [0]] )), ) (points + lines).properties( width="container", height=300, title="Effect of observation weights on outlier robustness" ) ``` The weighted fit should produce a lower scale estimate and fewer effective degrees of freedom because it is not distorted by the outliers. # Convergence ## What convergence means The P-IRLS algorithm converges when the relative change in penalized deviance between successive iterations falls below the tolerance threshold (10^{-7} by default). Specifically, the algorithm stops when: \frac{\|D^{(k)} - D^{(k-1)}\|}{\|D^{(k-1)}\| + 0.1} \< \text{tol} The 0.1 in the denominator prevents numerical issues when the deviance is near zero. ## Checking the fit After fitting, inspect the model to verify the fit is reasonable: ``` python model_check = wk.GAM("y ~ s(x)") model_check.fit({"x": x, "y": np.sin(x) + rng.normal(0, 0.3, n)}, method="REML") print(f"EDF total: {model_check.edf_total:.1f}") print(f"Scale: {model_check.scale:.4f}") print(f"Deviance: {model_check.deviance:.2f}") ``` EDF total: 7.6 Scale: 0.0930 Deviance: 8.59 For Gaussian models with an identity link, the P-IRLS algorithm converges in a single step. For non-Gaussian models, convergence typically occurs within 5-15 iterations. > **Warning: What to do if the fit looks wrong** > > 1. **Check the data**: fitting problems are often caused by complete separation (in logistic regression), extreme outliers, or a poor choice of family. > 2. **Reduce `k`**: an overly flexible model can be unstable. Reducing the basis dimension can help. > 3. **Try a different method**: GCV and REML can behave differently on ill-conditioned problems. > 4. **Run `model.check()`**: the basis dimension adequacy test can reveal underfitting (see [Diagnostics](diagnostics.md)). # Effective degrees of freedom (EDF) The effective degrees of freedom measures the complexity of each smooth term. It is defined as: \text{EDF}\_j = \text{tr}(\mathbf{A}\_j) where \mathbf{A}\_j is the block of the hat matrix corresponding to the j-th smooth. The total model EDF is \text{EDF}\_{\text{total}} = \sum_j \text{EDF}\_j + p, where p is the number of parametric (unpenalized) coefficients including the intercept. **Interpreting EDF**: | EDF value | Interpretation | |----|----| | \approx 1 | The smooth is approximately linear | | \approx 2 | The smooth is approximately quadratic | | 3-5 | Moderate curvature | | \> 5 | High complexity: consider checking the data for artifacts | | \approx k - 1 | The smooth is using nearly all available basis functions: increase `k` | ``` python # EDF interpretation model_edf = wk.GAM("y ~ s(x, k=20)") model_edf.fit({"x": x, "y": np.sin(x) + rng.normal(0, 0.3, n)}, method="REML") print(f"Total EDF: {model_edf.edf_total:.3f}") print(f"Scale estimate: {model_edf.scale:.4f}") ``` Total EDF: 10.293 Scale estimate: 0.0745 > **Note: Note** > > EDF is not the same as the basis dimension `k`. The parameter `k` sets an upper bound on complexity, while the smoothing parameter \lambda (selected by REML or GCV) determines how much of that capacity is actually used. Setting `k=20` does not mean the smooth uses 20 degrees of freedom -- it may use only 4 or 5 if the data support a simpler shape. See the [smooth terms](smooths.md) page for guidance on choosing `k`. # Scale estimation The scale parameter \hat\phi measures the residual variability after accounting for the smooth effects. For the Gaussian family, \hat\phi = \hat\sigma^2 is the residual variance. It is estimated by: \hat\phi = \frac{D(y, \hat\mu)}{n - \text{EDF}\_{\text{total}}} where D is the deviance and \text{EDF}\_{\text{total}} is the total effective degrees of freedom. This is the analogue of s^2 = \text{RSS} / (n - p) in ordinary least squares, replacing p with the (fractional) EDF. For families with a known scale (Poisson, binomial), \phi is fixed at 1 and is not estimated. ``` python # Scale estimation print(f"Scale (phi-hat): {model_edf.scale:.4f}") print(f"Deviance: {model_edf.deviance:.4f}") print(f"GCV score: {model_edf.gcv_score:.4f}") ``` Scale (phi-hat): 0.0745 Deviance: 6.6825 GCV score: 0.0830 The scale estimate directly affects confidence intervals and p-values: standard errors are proportional to \sqrt{\hat\phi}, so an overestimated scale produces wider confidence bands. # Practical guidance ## When to use REML vs GCV | Situation | Recommendation | |----|----| | Default, general-purpose fitting | **REML** | | Reproducing a legacy analysis that used GCV | GCV | | Comparing models with different fixed effects via LRT | ML | | Very large n (\> 50,000) | GCV may be faster (both give similar results) | | Sparse data or small n | **REML** (GCV tends to undersmooth) | ## When to increase `k` The basis dimension `k` should be large enough that the smooth can capture the true function shape. After fitting, run `model.check()` (see [Model diagnostics](diagnostics.md)) to test whether `k` is adequate: ``` python model.check() ``` If the k-index is below 1 with a significant p-value, double `k` and re-fit: ``` python model = wk.GAM("y ~ s(x, k=20)") model.fit(data, method="REML") model.check() ``` > **Tip: A safe workflow for production models** > > 1. Start with the defaults: `method="REML"`, `k=10`. > 2. Fit the model and run `model.check()`. > 3. If any smooth fails the k-index test, increase `k` for that term and re-fit. > 4. Review the summary: are the EDF values plausible for the data? > 5. If you suspect irrelevant predictors, re-fit with `select=True`. > 6. Generate [predictions](prediction.md) and [diagnostic plots](diagnostics.md). # Summary of fit attributes After calling `model.fit()`, the following attributes are available: | Attribute | Type | Description | |----|----|----| | `model.smoothing_params` | `list[float]` | Estimated (or fixed) \lambda_j values | | `model.edf_total` | `float` | Total effective degrees of freedom | | `model.scale` | `float` | Estimated scale parameter \hat\phi | | `model.deviance` | `float` | Model deviance | | `model.gcv_score` | `float` | GCV score (computed regardless of method) | # Where to go next - **[Prediction and inference](prediction.md)**: standard errors, confidence intervals, and term-level predictions from a fitted model. - **[Model diagnostics](diagnostics.md)**: residual plots, k-index tests, and `model.check()`. - **[Smoothing parameter sensitivity](sensitivity.md)**: verify that predictions are robust to the specific smoothing parameters chosen by REML, GCV, or ML. - **[Variational inference](variational-inference.md)**: a principled Bayesian posterior for non-Gaussian families, using `method="VI"`. - **[Smooth terms](smooths.md)**: choosing the right basis and basis dimension. - **[Response families](families.md)**: how the family and link function affect the P-IRLS algorithm.