\[ \DeclareMathOperator*{\argmin}{arg\,min} \DeclareMathOperator*{\argmax}{arg\,max} \]
This chapter is currently undergoing conversion and substantial changes will be made without notice. For the sections on practical data analysis we refer, for the time being, to the 2020-edition.
This chapter deals with a range of general problems and questions related to regression modeling and the practical process of data analysis. Some of the general considerations apply to statistics and data analysis in a broad sense, but to keep focus, data analysis is treated in the framework of regression models only.
In the first part of the chapter we position the linear model from Chapter 2 within a more general framework. This makes it possible to formulate some of the central objectives in regression analysis without relying on specific model assumptions.
In the subsequent parts we turn to the process of practical data analysis, which involves a lot of decision making. For an inexperienced data analyst the decision making can be a daunting task, first of all because there is never a set of clearly correct decisions. Instead it appears as if many decisions have to be made on the basis of little knowledge, and whenever a decision is made, it may be open to criticism. It is therefore tempting to hide the decision process, but this can be detrimental to the scientific process.
We will develop methods for understanding and controlling the decision making process involved in practical data analysis. This ranges from grand decisions, like what the overall purpose of the analysis is, to technical decisions, such as how to handle missing values. Central to the decision process is a a solid understanding of how we can interpret an associational regression model and what it can ultimately be used for. For this reason we will also briefly treat some central causal questions and how and when we can use regression models to answer them.
The chapter concludes in Section 7.3.1 by introducing basis expansion techniques to capture nonlinear relations within the framework of the linear model. The techniques are illustrated by a spline basis expansion of the relation between claim size and insurance sum in the fire insurance example.
The chapter primarily relies on the following four R packages.
library(RwR)
library(ggplot2)
library(tibble)
library(broom)
library(splines)
library(dplyr)
Attaching package: 'dplyr'The following objects are masked from 'package:stats':
filter, lagThe following objects are masked from 'package:base':
intersect, setdiff, setequal, unionThere exists a range of generalizations of the linear regression model. It would be misleading to lump them all together under the umbrella term nonlinear regression. Some generalizations focus specifically on nonlinear relations between the response and the predictors, while other generalizations focus on other targets than the conditional expectation, e.g., the median or certain quantiles.
We will consider various generalizations in this book. In Chapter 5 the target is still the conditional expectation – with a particular form of nonlinearity. But just as importantly, the chapter generalizes both the estimation methodology and the statistical theory to account for variance heterogeneity.
In Chapter 11 the target is the conditional hazard function of a survival distribution. One particular class of such survival regression models are models of the conditional expectation of the log-survival time, but most models are not directly targeting a conditional expectation.
In this chapter we treat some general aspects of regression modeling without making any of the linear modeling assumptions A1, A2 or A3, but for clarity, the modeling target will remain the conditional expectation throughout.
Definition 7.1 Let \(X \in \mathbb{R}^p\) and \(Y \in \mathbb{R}\) be random variables and suppose that \(\mathbf{E}(|Y|) < \infty\). The regression function \(m : \mathbb{R}^p \to \mathbb{R}\) is defined by \[ m(x) = \mathbf{E}(Y \mid X = x) \tag{7.1}\]
The regression function \(m\) is in general only uniquely defined almost surely w.r.t. to the distribution of \(X\). Though if it has a continuous version, that version is unique on the support of \(X\). In most applications it is reasonable to assume that \(m\) is continuous.
A general regression model is a set \(\mathcal{M}\) of regression functions. That is, \[ \mathcal{M} \subseteq \{ m \mid m : \mathbb{R}^p \to \mathbb{R} \}. \] We see that the linearity assumption A1 can be expressed in terms of the regression function as \[ m(x) = x^T \beta, \] which means that the corresponding linear regression model can be expressed as \[ \mathcal{M}_{\mathrm{linear}} = \{ x \mapsto x^T \beta \mid \beta \in \mathbb{R}^p \}. \] The linear regression model is parametrized by the parameter \(\beta\) belonging to the parameter set \(\mathbb{R}^p\).
Most regression models are parametrized by a parameter, abstractly denoted \(\theta \in \Theta\), and we then write the model as \[ \mathcal{M} = \{m_\theta \mid \theta \in \Theta \}. \]
Example 7.1 (Dose-response curves) A simple classical example is the logistic function for \(p = 1\): \[ m_{(a,b,c)}(x) = \frac{\beta}{1 + e^{-(b x + a)}} + \alpha \tag{7.2}\] The parameter \(\theta = (a, b, \alpha, \beta)\) is \(4\)-dimensional and the parameter set is \(\Theta = \mathbb{R}^4\). The S-shaped graph of a logistic function can, for instance, model the relation between a dosage (\(x\)) and a continuous response (\(y\)).
Example 7.2 (Neural networks) Let \[ \mu, h: \mathbb{R} \to \mathbb{R} \] be two fixed (nonlinear) functions. With \(\mathbf{B}\) a real \(p \times q\) matrix, \(\beta \in \mathbb{R}^{q}\) and \(\alpha \in \mathbb{R}\) we define the regression function \[ m_{\mathbf{B}, \alpha, \beta} (x) = \mu(\alpha + h(x^T \mathbf{B}) \beta) \] parametrized by \((\mathbf{B}, \alpha, \beta)\). The convention above is that \(h\) is applied coordinatewisely to the \(q\)-dimensional (row)-vector \(x^T \mathbf{B}\).
This is an example of a neural network with a single hidden layer. The matrix \(\mathbf{B}\) parametrizes the input layer mapping the predictor vector \(x \in \mathbb{R}^p\) to \(h(x^T \mathbf{B}) \in \mathbb{R}^{q}\). It is a composition of the parametrized linear map \(x \mapsto x^T \mathbf{B}\) and a fixed coordinatewise nonlinear function \(h\). The number \(\alpha\) and the vector \(\beta\) parametrize the output layer, which is likewise a composition of a linear map and the nonlinear function \(\mu\).
We introduce the notation \(L_{\alpha, \mathbf{B}}^{h} (x) = h(\alpha + x^T \mathbf{B})\) for an activation function \(h\), a \(p \times q\) matrix \(\mathbf{B}\) and a vector \(\alpha \in \mathbb{R}^q\). Then \(L_{\alpha, \mathbf{B}}^h : \mathbb{R}^p \to \mathbb{R}^q\) and it represents a single layer in a neural network. Using this notation we can write the neural network regression function above as a composition of two such maps: \[ m_{\mathbf{B}, \alpha, \beta} = L_{\alpha, \beta}^\mu \circ L_{0, \mathbf{B}}^{h} \]
A general neural network with \(k\) hidden layers can be defined as \[ \begin{align*} m_{\alpha_1, \mathbf{B}_1, \ldots, \alpha_k, \mathbf{B}_k, \alpha, \beta} %& = \mu(h_k( \cdots ( h_2(h_1(x^T \mathbf{B}_1)\mathbf{B}_2 + \beta_2) \cdots \mathbf{B}_k) + \beta_k)\beta + \beta_0) \\ & = L_{\alpha, \beta}^\mu \circ L_{\alpha_k, \mathbf{B}_k}^{h_k} \circ \ldots \circ L_{\alpha_1, \mathbf{B}_1}^{h_1} \end{align*} \] The model is parametrized by \((\alpha_1, \mathbf{B}_1, \ldots, \alpha_k, \mathbf{B}_k, \alpha, \beta)\) where \(\mathbf{B}_{j}\) is a \(p_{j-1} \times p_j\) real matrix (\(p_0 = p\)), \(\alpha_j \in \mathbb{R}^{p_j}\), \(\beta \in \mathbb{R}^{p_k}\) and \(\alpha \in \mathbb{R}\).
It is quite common to use the same activation function for all layers, that is, \(h_1 = h_2 = \ldots = h_k = \mu\). A couple of standard choices are: \[ \begin{align*} \mu(\eta) & = (1 + e^{-\eta})^{-1} \qquad \mathrm{(logistic)} \\ \mu(\eta) & = \log(1 + e^\eta) \qquad \ \mathrm{(softplus)} \\ \mu(\eta) & = \max(0, \eta) \qquad \ \ \ \mathrm{(ReLU)} \end{align*} \]
We will in Chapter 5 treat generalized linear models, which can be regarded as neural networks with no hidden layers. That is, the regression function is of the form \[ m_\beta(x) = \mu(x^T \beta) \] We will in that context see canonical choices of the activation function \(\mu\) dictated by the distribution of the response variable.
Even though the parameter space of a neural network can be very large – in particular for deep neural networks with many hidden layers – the parameter space is formally a finite dimensional space. We might want to choose a model \(\mathcal{M}\) that is even larger, e.g., the set of all continuous functions, the set of all Lipschitz functions, or the set of all smooth functions. These sets cannot be parametrized in a nice way by a finite dimensional parameter1, and we call the corresponding regression models nonparametric.
1 More precisely, they are not finite dimensional differentiable manifolds.
There are substantial theoretical differences between trying to estimate a regression function in a parametric model or in a nonparametric model. From a practical perspective, though, the transition from small parametric models to nonparametric models via large parametric models is much more gradual. We will not deal with estimation methods or theory that are explicitly nonparametric, but it can be conceptually advantageous to define and discuss the purposes and objectives of a regression analysis in a nonparametric way. Doing so, we will be able to define our aims clearly without reference to a particular parametric model – even though we eventually, in this book, develop solutions in terms of specific parametric models.
The most widely used estimator of the regression function is the least squares estimator minimizing the squared error loss. With a model \(\mathcal{M}\) and observations \((X_1, Y_1), \ldots, (X_n, Y_n)\) it is formally defined as \[ \hat{m} = \argmin_{m \in \mathcal{M}} \sum_{i=1}^n (Y_i - m(X_i))^2. \tag{7.3}\] If the model is parametrized by \(\theta \in \Theta\), we can formulate the estimator in terms of estimating \(\theta\) \[ \hat{\theta} = \argmin_{\theta \in \Theta} \sum_{i=1}^n (Y_i - m_{\theta}(X_i))^2. \tag{7.4}\]
Example 7.3 (Birth weight)
There are several practical problems with these formal definitions, though. First, there may be no global minimizer of the sum of squares or there may be many. Second, it may be practically impossible to compute any such global minimizer unless \(\mathcal{M}\) is particularly nice, e.g., \(\mathcal{M} = \mathcal{M}_{\mathrm{linear}}\) as in Theorem 2.1. Third, even if we could compute \(\hat{m}\) it might not be a good estimator if \(\mathcal{M}\) is very large relative to the sample size \(n\).
Many practical estimation algorithms do not attempt to litterally compute \(\hat{m}\) or \(\hat{\theta}\) as defined by (7.3) and (7.4), respectively. Typically, they iteratively find and improve candidate regression functions in the model – using the squared error loss as guide rather than trying to drive it to its minimal value. This can be done in many ways, e.g., by explicitly restricting the model space to an \(n\)-dependent subset \(\mathcal{M}_n \subseteq \mathcal{M}\), or by regularizing the estimator via a penalty function \(\mathrm{pen} : \mathcal{M} \to \mathbb{R}\). In the latter case, the algorithm would then (try to) compute
\[
\hat{m}^{\lambda_n} = \argmin_{m \in \mathcal{M}} \sum_{i=1}^n (Y_i - m(X_i))^2 + \lambda_n \mathrm{pen}(m)
\tag{7.5}\] for an \(n\)-dependent parameter \(\lambda_n > 0\) controlling the tradeoff between the squared error loss and the penalty term. Typically, \(\lambda_n \to 0\) for \(n \to \infty\) so that the regularization has less and less impact on the estimator as the sample size increases. Section 7.3.1 demonstrated by example the use of penalization within the framework of a linear model, where the penalization was used to control the fluctuations of nonlinear interaction terms.
We can interpret the term \(\lambda_n \mathrm{pen}(m)\) as an inductive bias, where we bias the estimator toward small values of \(\mathrm{pen}(m)\). Another way to introduce an inductive bias is by early stopping. If we start the optimization algorithm in a particular subset \(\mathcal{M}_0\) and then stop it before convergence, the result is implicitly biased toward models in \(\mathcal{M}_0\). An algorithm could implement early stopping by monitoring the squared error loss on an independent dataset and stop when that error does not improve anymore. Similarly, the parameter \(\lambda_n\) could be selected by minimizing the squared error loss on an independent dataset over a range of possible values of \(\lambda_n\).
For a notable part of this book we consider parametrized models and loss functions that are nice enough for the estimators we consider to actually be global minimizers, such as \(\hat{\theta}\) in (7.4). Moving beyond those models, regression function estimators used in practice are mostly the result of running a complicated algorithm on data, which does not have a simple analytic representation as, e.g., a minimizer or a solution to an equation. This does not make such estimators bad, but it does make it somewhat more complicated to understand their theoretical properties. The least squares estimator in \(\mathcal{M}_{\mathrm{linear}}\), as treated in detail in Chapter 2, stands out as having a particularly well understood theory, which we can leverage for practical statistical data analysis. We cannot obtain a similarly simple theory for general regression function estimators, but we formulate our modeling objectives in terms of general regression functions in clear and simple way that does not hinge on simplistic model assumptions.
As discussed at some length in the introduction, Chapter 1, one purpose of estimating the regression function \(m\) is to use \(\hat{m}\) to make predictions about \(Y\) from observation of \(X\). This is the primary focus in most of the machine learning literature. Though we will cover this aspect in later chapters – where we discuss quantification and optimization of predictive strength of a model – we also deal at length with other objectives.
One main question we consider is:
Does \(m\) depend on \(x_i\)?
To formulate this precisely, let \(x_{-j} = (x_1, \ldots, x_{j-1}, x_{j+1}, \ldots, x_p)^T\) denote the \((p-1)\)-dimensional vector with the \(j\)-th coordinate removed from \(x\), and let \[ m_{-j}(x_{-j}) = \mathbf{E}(Y \mid X_{-j} = x_{-j}) \] denote the regression function based on the \(X_{-j}\)-predictors only. We can then formulate the hypothesis \[ H_0: \mathbf{E}(Y \mid X) = m_{-j}(X_{-j}). \tag{7.6}\] This hypothesis means that \(m\) does not depend on the \(i\)-th coordinate, and it is an example of a conditional (mean) independence hypothesis. We will develop statistical theory to test \(H_0\) based on various regression modeling frameworks.
Another main objective we consider is to:
Quantify (and estimate) how much \(m\) depends on a particular coordinate \(x_i\).
Obviously, this objective is closely related to the hypothesis above – any sensible quantification should show that \(m\) does not depend on \(x_i\) if \(H_0\) is true. Additionally, such a quantification is naturally tied together with the development of a test statistic of \(H_0\). Testing \(H_0\) is, however, a simpler problem, and the challenge is to develop easily interpretable quantifications that are simultaneously easy to estimate. The regression modeling frameworks we consider will have some natural ways of quantifying such conditional associations.
To aid the general discussion of the two objectives above, we introduce partial dependence functions and partial errors below. To ease notation we will use the convention that \[ m(x_j, x_{-j}) = m((x_1, \ldots, x_{j-1}, x_j, x_{j+1}, \ldots, x_{p})^T), \] where \(m(x_j, x_{-j})\) is simply a shorthand notation for the right hand side above. Recall that \(m(x) = \mathbf{E}(Y \mid X = x)\) denotes the joint regression function of \(Y\) on \(X\).
Definition 7.2 The partial dependence function \(m_j : \mathbb{R} \to \mathbb{R}\) for the \(j\)-th coordinate is \[ m_j(x_j) = \mathbf{E}(m(x_j, X_{-j})) = \mathbf{E}(\mathbf{E}(Y \mid X_j = x_j, X_{-j})). \tag{7.7}\]
The partial dependence function is also known as the adjusted regression function. It is essential for its definition that we integrate out over the marginal distribution of \(X_{-j}\) and not the conditional distribution of \(X_{-j}\) given \(X_{j} = x_j\).
Recall that \(m_{-j}(x_{-j}) = \mathbf{E}(Y \mid X_{-j} = x_{-j})\) is the regression function of \(Y\) on \(X_{-j}\) excluding the \(j\)-th coordinate of \(X\).
Definition 7.3 The \(j\)-th partial error is defined as \[ \varepsilon_{-j} = Y - m_{-j}(X_{-j}). \tag{7.8}\]
Note that the tower property of conditional expectations imply the identity \[ m_{-j}(x_{-j}) = \mathbf{E}( m(X_j, x_{-j}) \mid X_{-j} = x_{-j}), \] which relate \(m_{-j}\) to \(m\). We can now show the following result on implications of the hypothesis \(H_0\) for the partial dependence function and the partial error.
Proposition 7.1 If the regression function \(m(x) = \mathbf{E}(Y \mid X = x)\) does not depend on \(x_j\) then
Proof. It follows directly from (7.7) that if \(m\) does not depend on \(x_j\) then \(m_{j}(x_j) = \mathbf{E}(m(X_{-j}))\) is constant.
To prove the second claim we note that if \(m\) does not depend on the \(j\)-th coordinate, \(m(X_j, x_{-j}) = m(x_j, x_{-j})\) for any fixed \(x_j\), whence \[ m_{-j}(x_{-j}) = \mathbf{E}( m(X_j, x_{-j}) \mid X_{-j} = x_{-j}) = m(x_j, x_{-j}). \]
Finally, if \(m\) does not depend on \(x_j\), \(\varepsilon_{-j} = \varepsilon = Y - \mathbf{E}(Y \mid X)\). Whence by the tower property of conditional expections, \[ \begin{align*} \mathrm{cov}(\varepsilon_{-j}, X_j) & = \mathrm{cov}(\varepsilon, X_j) \\ & = \mathbf{E}( \varepsilon X_j ) \\ & = \mathbf{E}( \mathbf{E}(\varepsilon X_j \mid X)) \\ & = \mathbf{E}( \mathbf{E}(\varepsilon \mid X) X_j) = 0, \end{align*} \] where we have used that \(\mathbf{E}(\varepsilon \mid X) = 0\).
There are various procedures, motivated by Proposition 7.1, that are used in practice to better understand how \(m\) depends on a particular predictor. To describe them, let \(\hat{m}\) and \(\hat{m}_{-j}\) denote estimates of \(m\) and \(m_{-j}\), respectively.
The suggested plots of the partial residuals against either the \(j\)-th predictor or its residual are both attempts to visualize dependencies between \(\varepsilon_{-j}\) and \(X_j\). The plots can be used to graphically detect if \(\mathrm{cov}(\varepsilon_{-j}, X_j) = 0\) or not. We can also use the partial residuals to directly estimate the covariance, either as \[ \tilde{\mathrm{cov}}_j = \frac{1}{n} \sum_{i=1}^n \hat{\varepsilon}_{i, -j} X_{i,j} \tag{7.9}\] or as \[ \hat{\mathrm{cov}}_j = \frac{1}{n} \sum_{i=1}^n \hat{\varepsilon}_{i, -j} (X_{i,j} - \hat{r}(X_{i,-j})). \tag{7.10}\]
Both of these covariance estimates quantify deviations from \(H_0\), and in this way indirectly also how much \(m\) depends on the \(j\)-th coordinate. The estimator (7.9) has the advantage that we do not need to estimate the regression function \(r_j\), but it can be severely biased and \(\hat{\mathrm{cov}}_j\) is preferred in practice.
To better appreciate the general definitions and procedures outlined above, we go through four different examples in increasing complexity, where most of the definitions and quantities above have some explicit expressions in terms of model parameters or components.
Example 7.4 If \(m(x) = x^T \beta\) is linear we find that the partial dependence function is \[ m_j(x_j) = \beta_j x_j + \mathbf{E}(X_{-j})^T \beta_{-j}. \] We see that the partial dependence function is affine with slope \(\beta_j\). It is constant if and only if \(\beta_j = 0\), which is the case if and only if \(m\) does not depend on \(x_j\).
The \(F\)-test statistic defined above coincides for this model with the \(F\)-test statistic defined by Equation 2.11, and it can be used to formally test \(H_0\).
We also find that \[ m_{-j}(x_{-j}) = x_{-j}^T \beta_{-j} + \beta_j \mathbf{E}(X_j \mid X_{-j} = x_{-j}) = x_{-j}^T \beta_{-j} + \beta_j r_j(x_{-j}). \] Therefore \[ m(x) - m_{-j}(x_{-j}) = \beta_j(x_j - r_j(x_{-j})) \] This difference is non-zero only if \(\beta_j \neq 0\). Moreover, \[ \varepsilon_{-j} = \varepsilon + \beta_j(X_j - r_j(X_{-j})) \] and we find that \[ \mathrm{cov}(\varepsilon_{-j}, X_j) = \beta_j \mathbf{V}(X_j - r_j(X_{-j})). \] This covariance is thus proportional to \(\beta_j\), with a proportionality constant \(\mathbf{V}(X_j - r_j(X_{-j}))\) depending on the distribution of \(X\) only.
Example 7.5 The regression function is called partially linear in the \(j\)-th variable if it is of the form \[ m(x) = \beta_j x_j + h_{-j}(x_{-j}). \] In this case, similarly as above, the partial dependence function is an affine function \[ m_j(x_j) = \beta_j x_j + \mathbf{E}(h_{-j}(X_{-j})) \] with slope \(\beta_j\). Moreover, \[ m_{-j}(x_{-j}) = h_{-j}(x_{-j}) + \beta_j \mathbf{E}(X_j \mid X_{-j} = x_{-j}) = h_{-j}(x_{-j}) + \beta_j r_j(x_{-j}), \] and thus \[ m(x) - m_{-j}(x_{-j}) = \beta_j(x_j - r_j(x_{-j})) \] exactly as if \(m\) were linear. By the same argument as above, we get that \[ \mathrm{cov}(\varepsilon_{-j}, X_j) = \beta_j \mathbf{V}(X_j - r_j(X_{-j})). \]
The computations above for a partially linear regression function reveal that the partial dependence plot as well as procedures based on the partial residuals should reveal the same affine dependence of \(m\) on \(x_j\) whether \(m\) is linear or only partially linear. The subtle, but important, practical difference is that this is only true if the estimators used for estimating the partial dependence function or for computing the partial residuals are flexible enough to correctly model the potentially nonlinear functions \(h_{-j}\) and \(m_{-j}\).
Example 7.6 If the regression function has the form \[ m(x) = \beta_j(x_{-j}) x_j + h_{-j}(x_{-j}), \] it is affine in \(x_j\) for any fixed \(x_{-j}\). It is partially linear in the \(j\)-th variable but with a coefficient that depends on \(x_{-j}\). The partial dependence function is the affine function \[ m_j(x_j) = \mathbf{E}(\beta_j(X_{-j})) x_j + \mathbf{E}(h_{-j}(X_{-j})) \] with slope \(\mathbf{E}(\beta_j(X_{-j}))\) being the mean slope of \(x_j\) in \(m\). \[ m_{-j}(x_{-j}) = h_{-j}(x_{-j}) + \beta_j(x_{-j}) \mathbf{E}(X_j \mid X_{-j} = x_{-j}) = h_{-j}(x_{-j}) + \beta_j(x_{-j}) r_j(x_{-j}), \] and thus \[ m(x) - m_{-j}(x_{-j}) = \beta_j(x_{-j})(x_j - r_j(x_{-j})). \] The covariance becomes \[ \begin{align*} \mathrm{cov}(\varepsilon_{-j}, X_j) & = \mathbf{E}(\beta_j(X_{-j})(X_j - r_j(X_{-j}))^2) \\ & = \mathbf{E}(\beta_j(X_{-j})\mathbf{V}(X_j - r_j(X_{-j}) \mid X_{-j})) \\ & = \mathbf{E}(\beta_j(X_{-j})v_j(X_{-j})), \\ \end{align*} \] which we can interpret as a weighted mean of the individual slopes with weights \(v_j(X_{-j}) = \mathbf{V}(X_j - r_j(X_{-j}) \mid X_{-j})\) being conditional variances.
Example 7.7 The regression function is called partially additive in the \(j\)-th variable if it is of the form \[ m(x) = m^0_j(x_j) + h_{-j}(x_{-j}). \] The partial dependence function is then \[ m_j(x_j) = m^0_j(x_j) + \mathbf{E}(h_{-j}(X_{-j})), \] which up to an additive constant coincides with \(m^0_j\). The regression function \(m_{-j}\) becomes \[ m_{-j}(x_{-j}) = h_{-j}(x_{-j}) + \mathbf{E}(m^0(X_j) \mid X_{-j} = x_{-j}) = h_{-j}(x_{-j}) + w_{-j}(x_{-j}) \] where \(w_{-j}(x_{-j}) = \mathbf{E}(m^0_j(X_j) \mid X_{-j} = x_{-j})\). Therefore \[ m(x) - m_{-j}(x_{-j}) = m^0_j(x_j) - w_{-j}(x_{-j}). \] Finally \[ \begin{align*} \mathrm{cov}(\varepsilon_{-j}, X_j) & = \mathrm{cov}(m^0_j(X_j) - w_{-j}(X_{-j}), X_j) \\ & = \mathbf{E}( \mathrm{cov}(m^0_j(X_j), X_j \mid X_{-j})). \end{align*} \]
We found in the earlier analysis of the insurance claims data a lack of model fit. This could be due to a nonlinear relation between log(claims) and log(sum). To handle nonlinear relations between the response and one or more predictors, the predictors (as well as the response) can be nonlinearly transformed before they enter into the linear model. Such pre-modeling transformations extend the scope of the linear model considerably. For the insurance data the original variables were already log-transformed, and it is not obvious which other transformation to use.
x <- seq(11, 21, 0.05)
pol <- as_tibble(cbind(x, poly(x, degree = 5)))
long_pol <- tidyr::pivot_longer(pol, cols = !x, names_to = "basis_fct")
p <- ggplot(NULL, aes(x, value, colour = basis_fct)) +
scale_color_discrete(guide = "none") + xlab("") + ylab("")
p + geom_line(data = long_pol, linewidth = 1)
An alternative to data transformations is to use a small but flexible class of basis functions, which can capture the nonlinearity. One possibility is to use low degree polynomials. This can be done by simply including powers of a predictor as additional predictors. [Due to the parsing of formulas in R, powers have to be wrapped into an I()}-call, e.g., y ~ x + I(x2).] Depending on the scale and range of the predictor, this may work just fine. However, raw powers can result in numerical difficulties. An alternative is to use orthogonal polynomials, which are numerically more well behaved.
Figure 7.1 shows an example of an orthogonal basis of degree 5 polynomials on the range 11–21. This corresponds approximately to the range of log(sum), which is a target for basis expansion in the example. What should be noted in Figure 7.1 is that the behavior near the boundary is quite erratic. This is characteristic for polynomial bases. To achieve flexibility in the central part of the range, the polynomials become erratic close to the boundaries. Extrapolation beyond the boundaries cannot be trusted.
An alternative to polynomials is splines. A spline2 is piecewisely a polynomial, and the pieces are joined together in a sufficiently smooth way. The points where the polynomials are joined together are called knots. The flexibility of a spline is determined by the number and placement of the knots and the degree of the polynomials. A degree \(k\) spline is required to be \(k-1\) times continuously differentiable. A degree 3 spline, also known as a cubic spline, is a popular choice, which thus has a continuous (piecewise linear) second derivative.
2 A spline is also a thin and flexible wood or metal strip used for smooth curve drawing.
b_spline <- as_tibble(cbind(x, splines::bs(x, knots = c(14, 18))))
long_b_spline <- tidyr::pivot_longer(b_spline, cols = !x, names_to = "basis_fct")
p + geom_line(data = long_b_spline, linewidth = 1) +
geom_rug(aes(x = c(11, 14, 18, 21), y = NULL, color = NULL), linewidth = 1)
bs() function with two internal knots at 14 and 18 in addition to two boundary knots at 11 and 21.
Figure 7.2 shows a basis of 5 cubic spline functions. They are so-called \(B\)-splines (basis splines). Note that it is impossible to visually detect the knots where the second derivative is non-differentiable. The degree \(k\) \(B\)-spline basis with \(r\) internal knots has \(k + r\) basis functions. The constant function is not included then, and this is precisely what we want when the basis expansion is used in a regression model including an intercept. As seen in Figure 7.2, the \(B\)-spline basis is also somewhat erratic close to the boundary. For a cubic spline, the behavior close to the boundary can be controlled by requiring that the second and third derivatives are \(0\) at the boundary knots. The result is known as a natural cubic spline. The extrapolation (as a spline) of a natural cubic spline beyond the boundary knots is linear.
Due to the restriction on the derivatives of a natural cubic spline, the basis with \(r\) internal knots has \(r + 1\) basis functions. Thus the basis for the natural cubic splines with \(r + 2\) internal knots has the same number of basis functions as the raw cubic \(B\)-spline basis with \(r\) internal knots. This means in practice that, compared to using raw \(B\)-splines with \(r\) internal knots, we can add two internal knots, and thus increase the central flexibility of the model, while retaining its complexity in terms of \(r + 3\) parameters.
n_spline <- as_tibble(cbind(x, splines::ns(x, knots = c(13, 15, 17, 19))))
long_n_spline <- tidyr::pivot_longer(n_spline, cols = !x, names_to = "basis_fct")
p + geom_line(data = long_n_spline, linewidth = 1) +
geom_rug(aes(x = c(11, 13, 15, 17, 19, 21), y = NULL, color = NULL), linewidth = 1)
ns() function with internal knots at 13, 15, 17 and 19 in addition to the two boundary knots at 11 and 21.
It is possible to use orthogonal polynomial expansions as well as spline expansions in R together with the lm() function via the functions poly(), bs() and ns(). The latter two are in the splines package, which thus has to be loaded.
We will illustrate the usage of basis expansions for claim size modeling using natural cubic splines.
claims_lm_spline_add <- lm(
log(claims) ~ ns(log(sum), knots = c(13, 15, 17, 19)) + grp,
data = claims
)
claims_lm_spline_int <- lm(
log(claims) ~ ns(log(sum), knots = c(13, 15, 17, 19)) * grp,
data = claims
)When we fit a model using basis expansions the coefficients of individual basis functions are rarely interpretable or of interest. Thus reporting the summary including all the estimated coefficients is not particularly informative. Below we only extract the last lines from the summary to report the residual variance and the \(R^2\) quantities.
glance(claims_lm_spline_add) |>
knitr::kable()| r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC | BIC | deviance | df.residual | nobs |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0.06416 | 0.0623 | 1.9185 | 34.509 | 0 | 8 | -8352 | 16724 | 16787 | 14822 | 4027 | 4036 |
glance(claims_lm_spline_int) |>
knitr::kable()| r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC | BIC | deviance | df.residual | nobs |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0.07294 | 0.06762 | 1.9131 | 13.724 | 0 | 23 | -8333 | 16716 | 16874 | 14683 | 4012 | 4036 |
Compared to the additive model without the basis expansion, \(R^2\) (also the adjusted one) is increased a little by the basis expansion. A further increase is observed for the model that includes interactions as well as a basis expansion, and thus gives a unique nonlinear relation between insurance sum and claim size for each trade group. Figure 7.4 shows some diagnostic plots for the additive model, which show that the mean value model is now pretty good, while the residuals are still slightly right skewed. The lack of normality is not a problem for using the sampling distributions to draw inference. However, if we compute prediction intervals, say, then we must account for the skewness in the residual distribution. This can be done by using quantiles for the empirical distribution of the residuals instead of quantiles for the normal distribution.
claims_diag_spline <- augment(claims_lm_spline_add) ## Residuals etc.
gridExtra::grid.arrange(
ggplot(claims_diag_spline, aes(.fitted, .resid), alpha = 0.2) + geom_point() + geom_smooth(),
ggplot(claims_diag_spline, aes(.resid)) + geom_histogram(bins = 40),
ggplot(claims_diag_spline, aes(sample = .std.resid)) + geom_qq() + geom_abline(),
ncol = 3
)Warning in fortify(data, ...): Arguments in `...` must be used.
✖ Problematic argument:
• alpha = 0.2
ℹ Did you misspell an argument name?
Graphical comparisons are typically preferable over parameter comparisons when models involving nonlinear effects and basis expansions are compared.
p0 <- ggplot(
data = claims,
aes(sum, claims),
alpha = I(0.2)
) +
scale_x_log10(
"Insurance sum (DKK)",
breaks = 10^c(6, 7, 8),
labels = c("1M", "10M", "100M")
) +
scale_y_log10(
"Claim size (DKK)",
breaks = 10^c(2, 4, 6, 8),
labels = c("100", "10K", "1M", "100M")
) +
geom_point()Warning in fortify(data, ...): Arguments in `...` must be used.
✖ Problematic argument:
• alpha = I(0.2)
ℹ Did you misspell an argument name?pred_add <- cbind(
claims,
exp(predict(claims_lm_spline_add, interval = "confidence"))
)
pred_int <- cbind(
claims,
exp(predict(claims_lm_spline_int, interval = "confidence"))
)
p0 <- p0 + facet_wrap(~ grp, ncol = 4) + coord_cartesian(ylim = c(100, 1e6))
p0 + geom_line(aes(y = fit), pred_add, color = "blue", linewidth = 2) +
geom_ribbon(aes(ymax = upr, ymin = lwr), pred_add, fill = "blue", alpha = 0.2) +
geom_line(aes(y = fit), pred_int, color = "red", linewidth = 2) +
geom_ribbon(aes(ymax = upr, ymin = lwr), pred_int, fill = "red", alpha = 0.2)
Figure 7.5 shows the fitted models stratified according to trade group. The fitted values were computed using the predict() function in R, which can also give upper and lower confidence bounds on the fitted values. The fitted value for the \(i\)-th observation is \(X_i^T \hat{\beta}\). By Theorem 2.3 a confidence interval for the fitted value can be computed using Equation 2.10 with \(a = X_i\). The confidence intervals reported by predict() for lm objects are computed using this formula, and these are the pointwise confidence bands shown in Figure 7.5. From this figure any gain of the interaction model over the additive model is questionable. The increase of \(R^2\), say, appears to be mostly driven by spurious fluctuations. For trade group 4 the capped claims seem overly influential, and because trade group 3 has relatively few observations the nonlinear fit is poorly determined for this group.
claims_lm_add <- lm(
log(claims) ~ log(sum) + grp,
data = claims
)anova(
claims_lm_add,
claims_lm_spline_add,
claims_lm_spline_int
) |> knitr::kable()| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 4031 | 14902.71 | NA | NA | NA | NA |
| 4027 | 14822.42 | 4 | 80.28605 | 5.484225 | 0.0002115 |
| 4012 | 14683.37 | 15 | 139.04900 | 2.532863 | 0.0009413 |
In the final comparison above using \(F\)-tests we first test the hypothesis that the relation is linear in the additive model, and then we test the additive hypothesis when the nonlinear relation is included. The tests formally reject both hypotheses, but the graphical summary in Figure 7.5 suggests that the nonlinear additive model is preferable anyway.
Introducing a penalty can dampen the spurious fluctuations the arose from the nonlinear basis expansion. We conclude the analysis of the insurance data by illustrating how the penalized least squares fit can be computed, and what effect the penalty has on the least squares solution. First recall that the least squares estimator could have been computed by solving the normal equation directly.
X <- model.matrix(claims_lm_spline_int)
y <- model.response(model.frame(claims_lm_spline_int))
XtX <- crossprod(X) # Faster than but equivalent to t(X) %*% X
Xty <- crossprod(X, y) # Faster than but equivalent to t(X) %*% y
coef_hat <- solve(XtX, Xty)
# Comparison
range(coef_hat - coefficients(claims_lm_spline_int))[1] -8.441248e-11 1.365330e-11The result obtained by solving the normal equation directly is up to numerical errors identical to the solution computed using lm(). The R function solve() calls the Fortran routine DGESV from the LAPACK library for solution of linear equations using LU decomposition with partial pivoting (Gaussian elimination with row permutations).
Likewise, we can compute the penalized estimator by solving the corresponding linear equation. We will do so where we only add a penalty on the coefficients corresponding to the interaction terms in the nonlinear basis expansion. We will use penalty matrices of the form \(\lambda \boldsymbol{\Omega}\) for a fixed matrix \(\boldsymbol{\Omega}\) so that the amount of penalization is controlled by the parameter \(\lambda > 0\).
Omega <- diag(c(rep(0, 9), rep(1, 15)))
coef_hat <- cbind(
solve(XtX + 0.1 * Omega, Xty), # Small penalty
solve(XtX + Omega, Xty), # Medium penalty
solve(XtX + 10 * Omega, Xty) # Large penalty
) Figure 7.6 shows the resulting penalized model compared to the unpenalized interaction model (\(\lambda = 0\)) and the additive model (\(\lambda = \infty\)). A minimal penalization is enough to dampen the most pronounced spurious fluctuation for the fourth trade group considerably.
tmp <- cbind(
coefficients(claims_lm_spline_int),
coef_hat,
c(coefficients(claims_lm_spline_add), rep(0, 15))
)
pred <- cbind(claims[, c("sum", "grp")], exp(X %*% tmp))
pred <- tidyr::pivot_longer(pred, cols = !c(grp, sum), names_to = "lambda")
p0 + geom_line(
aes(y = value, color = lambda),
data = pred,
linewidth = 1) +
theme(legend.position = "top") +
scale_color_discrete(
expression(~lambda),
labels = c(0, 0.1, 1, 10, expression(infinity~~~(additive~model)))
) 
The penalized least squares estimates can also be computed using the QR decomposition via the function lm.fit(). To achieve this we need to augment the model matrix and the response data suitably.
fit_lamb1 <- lm.fit(rbind(X, sqrt(0.1) * Omega), c(y, rep(0, ncol(Omega))))
fit_lamb2 <- lm.fit(rbind(X, sqrt(10) * Omega), c(y, rep(0, ncol(Omega))))
## Comparisons
range(coef_hat[, 1] - coefficients(fit_lamb1)) [1] -5.108247e-12 4.279577e-12range(coef_hat[, 3] - coefficients(fit_lamb2)) [1] -1.565414e-13 2.406964e-13The results are again (numerically) identical to using solve() above. If several penalized estimates using matrices \(\lambda \boldsymbol{\Omega}\) for different choices of \(\lambda\) are to be computed, a third option using a diagonalization of \(\mathbf{X}^T \mathbf{X}\) is beneficial. This is implemented in, e.g., the lm.ridge() function in the package MASS, which unfortunately doesn’t let you choose your own penalty matrix – it only implements \(\boldsymbol{\Omega} = \mathbf{I}\).
In conclusion, we have found evidence in the data for an overall nonlinear relation between insurance sum and claim size plus an additive effect of trade group. For small insurance sums (less than DKK 1M) there is almost no relation, while for larger insurance sums there is roughly a linear relation on a log-log scale. The model fits the data well with a non-normal and slightly right skewed residual distribution, but as a predictive model it is weak, and insurance sum can together with trade group only explain 6-7% of the variation in the (log) claim size distribution.
Basis expansions can be tried if we expect to get something out of it, that is, if we expect that there are some nonlinear relations in the data. On the other hand, if we construct models tailor-made to capture nonlinear relations we have spotted by eye-balling residual plots we run the risk of overfitting random fluctuations. In particular, formal downstream justifications using statistical tests are invalidated. The eye-balling process is a model selection procedure with statistical implications, which are difficult to account for.
We decided to expand gestationalAge using natural cubic splines with three knots in 38, 40, and 42 weeks. The boundary knots were determined by the range of the data set, and were thus 25 and 47. We also expanded age, but we let the ns() function determine the knots automatically for that predictor. The last continuous predictor, alcohol has a very skewed marginal distribution, and it was not judged to be suitable for a standard basis expansion. We also present a test of the nonlinear effect.
data("birth_weight")
birth_weight <- filter(
birth_weight,
weight > 32,
length > 10 & length < 99,
gestationalAge > 18
) |>
na.omit() |>
select(- c(interviewWeek, fetalDeath))
form <- weight ~ gestationalAge +
age +
children +
coffee +
alcohol +
smoking +
abortions +
feverEpisodes
birth_weight_lm <- lm(form, data = birth_weight)
binScale <- scale_fill_continuous(
breaks = c(1, 10, 100, 1000),
low = "gray80",
high = "black",
trans = "log",
guide = "none"
)nsg <- function(x) ns(x, knots = c(38, 40, 42), Boundary.knots = c(25, 47))
form <- weight ~ nsg(gestationalAge) +
ns(age, df = 3) +
children +
coffee +
alcohol +
smoking +
abortions +
feverEpisodes
birth_weight_lm_spline <- lm(form, data = birth_weight)
anova(birth_weight_lm, birth_weight_lm_spline) |> knitr::kable()| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 11139 | 2537626346 | NA | NA | NA | NA |
| 11134 | 2466413420 | 5 | 71212926 | 64.29455 | 0 |
gestationalAge against the main effects model.
The main effects model is a submodel of the nonlinear main effects model. This is not obvious. The nonlinear expansion will result in 4 columns in the model matrix \(\mathbf{X}\), none of which being gestationalAge. However, the linear function is definitely a natural cubic spline, and it is thus in the span of the 4 basis functions. With \(\mathbf{X}'\) the model matrix for the main effects model, it follows that there is a \(C\) such that Equation 2.12 holds. This justifies the use of the \(F\)-test. The conclusion from Table 7.4 is that the nonlinear model is highly significant.
birth_weight_diag <- augment(birth_weight_lm_spline, data = birth_weight)
p1 <- ggplot(birth_weight_diag, aes(.fitted, .std.resid)) +
stat_binhex(bins = 20) + binScale + geom_smooth(linewidth = 1, fill = "blue") +
xlab("fitted values") + ylab("standardized residuals")
p2 <- ggplot(birth_weight_diag, aes(gestationalAge, .std.resid)) +
stat_binhex(bins = 20) + binScale + geom_smooth(linewidth = 1, fill = "blue") +
xlab("gestationalAge") + ylab("")
p3 <- ggplot(birth_weight_diag, aes(sample = .std.resid)) +
geom_abline(intercept = 0, slope = 1, color = "blue", linewidth = 1) +
geom_qq() +
xlab("theoretical quantiles") + ylab("")
gridExtra::grid.arrange(p1, p2, p3, ncol = 3)
gestationalAge expanded using splines.
The positioning of the knots for spline expansion is selected based on the marginal distribution of the predictor3 The knots should be placed reasonably relative to the distribution of the predictor, so that we learn about nonlinearities where there is data to learn from. In this case we placed knots at the median and the 10% and 90% quantiles of the distribution of gestationalAge. The ns() function makes a similar automatic selection of knots based on the marginal distribution of the predictor variables when it is applied in the formula. One just have to specify the number of basis functions using the df argument.
3 Letting the knots be parameters to be estimated is not a statistically viable idea. The resulting estimation problem becomes computationally much more difficult, and it is not straightforward how to adjust subsequent statistical analyses to account for the data adaptive choice of knots.
An eye-ball decision based on Figure 3.8 would be that a single knot around 41 would have done the job, but we refrained from making such a decision. A subsequent test of the nonlinear effect with 1 degrees of freedom would not appropriately take into account how the placement of the knot was made.
Figure 7.7 shows diagnostic plots for the nonlinear main effects model. They show that the inclusion of the nonlinear effect removed the previously observed problem with assumption A1. The error distribution is still not normal, but right skewed with a fatter right tail than the normal distribution. There is a group of extreme residuals for preterm born children, which should be given more attention than they will be given here.
form <- weight ~ nsg(gestationalAge) +
children +
coffee +
smoking
birth_weight_lm_spline_small <- lm(form, data = birth_weight)
anova(birth_weight_lm_spline_small, birth_weight_lm_spline) |>
knitr::kable()| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 11142 | 2469351552 | NA | NA | NA | NA |
| 11134 | 2466413420 | 8 | 2938131 | 1.657931 | 0.1032394 |
Table 7.5 shows that dropping the four weak predictors from the nonlinear main effects model is again borderline and not really significant.
Figure 7.8 shows examples of the fit for the reduced nonlinear main effects model. The figure illustrates the general nonlinear relation between weight and gestationalAge. Differences due to other variables are purely additive in this model, which amounts to translations up or down of the curve. The figure shows a couple of extreme cases; the majority group who have had children before and who do not smoke or drink coffee, and a minority group who have had children before and smoke and drink coffee the most. What we should notice is the wider confidence band on the latter (smoke = 3, coffee = 3) compared to the former, which is explained by the skewness of the predictor distributions. Table 7.6 gives confidence intervals for the remaining parameters based on the nonlinear main effects model.
pred_frame <- expand.grid(
children = factor(1),
smoking = factor(c(1, 3)),
coffee = factor(c(1, 3)),
gestationalAge = seq(25, 47, 0.1),
alcohol = 0,
age = median(birth_weight$age),
feverEpisodes = 0,
abortions = factor(0)
)
pred <- predict(
birth_weight_lm_spline,
newdata = pred_frame,
interval = "confidence"
)
pred_frame <- cbind(pred_frame, pred)
ggplot(pred_frame, aes(gestationalAge, fit, color = coffee)) +
geom_line() + ylab("weight") +
geom_ribbon(
aes(ymin = lwr, ymax = upr, fill = coffee),
alpha = 0.3
) + facet_grid(. ~ smoking, label = label_both)
gestationalAge. Here illustrations of the fitted mean and 95% confidence bands for children=1.
| term | conf.low | conf.high |
|---|---|---|
| children1 | 155.470 | 194.087 |
| coffee2 | -82.608 | -42.431 |
| coffee3 | -193.020 | -87.630 |
| alcohol | -13.167 | 6.509 |
| smoking2 | -125.758 | -75.259 |
| smoking3 | -155.765 | -97.975 |
To conclude the analysis, we will again include interactions. However, the variable gestationalAge is in fact only taking the integer values 25 to 47, and the result of a nonlinear effect coupled with a third order interaction, say, results in an almost saturated model. That is, the interaction model has more or less a separate mean for all observed combinations of the predictors. We choose to consider a less complex model with all second order interactions between gestationalAge and the three factors that we have judged to be strong predictors.
form <- weight ~ (smoking + coffee + children) * nsg(gestationalAge) +
ns(age, df = 3) +
alcohol +
abortions +
feverEpisodes
birth_weight_lm_spline_int <- lm(form, data = birth_weight)pred_frame <- expand.grid(
children = factor(c(0, 1)),
smoking = factor(c(1, 2, 3)),
coffee = factor(c(1, 2, 3)),
gestationalAge = 25:47,
alcohol = 0,
age = median(birth_weight$age),
feverEpisodes = 0,
abortions = factor(0)
)
pred <- predict(birth_weight_lm_spline_int, newdata = pred_frame, interval = "confidence")
pred_frame <- cbind(pred_frame, pred)
pred_frame$fit_spline <- predict(birth_weight_lm_spline, newdata = pred_frame)
ggplot(pred_frame, aes(gestationalAge, fit)) +
geom_point(data = birth_weight, aes(y = weight), size = 0.3, alpha = 0.1) +
facet_grid(coffee ~ children + smoking, label = label_both) +
geom_ribbon(aes(ymin = lwr, ymax = upr), fill = gray(0.85)) +
geom_line(color = "red") + coord_cartesian(ylim = c(0, 5000)) +
geom_line(aes(y = fit_spline), color = "blue") + ylab("weight") +
scale_y_continuous(breaks = c(1000, 2000, 3000, 4000))
Figure 7.9 shows the predicted values for the reduced nonlinear main effects model for all combinations of children, smoking and coffee (blue curves). These curves are all just translations of each other. In addition, the figure shows the predicted values and a 95% confidence band as estimated with the nonlinear interaction model. We observe minor deviations for the reduced main effects model, which can explain the significance of the test, but the deviations appear unsystematic and mostly related to extreme values of gestationalAge. The conclusion is that even though the inclusion of interaction effects is significant, there is little to gain over the reduced nonlinear main effects model.
anova(birth_weight_lm_spline, birth_weight_lm_spline_int) |>
knitr::kable(digits = 10)| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 11134 | 2466413420 | NA | NA | NA | NA |
| 11114 | 2451152208 | 20 | 15261213 | 3.459865 | 2.612e-07 |
Table 7.7 shows that inclusion of interaction terms is still significant by a formal statistical test. It is, however, always a good idea to visualize how the nonlinear interaction model differs from the model with only main effects as above instead of just focusing on the formal test. The test may, on the other hand, support the visualization as a quantification of whether the differences seen in a figure are actually statistically significant.
Exercise 7.1 This is a continuation of Exercise Exercise 2.1.
set.seed(04022016)
factor_labels <- rep(c("A", "B", "C"), each = 50)
formula_data <- tibble(
x = rnorm(150, 0),
z = factor(sample(factor_labels)),
mu = (z == "A") + cos(x) + 0.5 * (z == "B") * sin(x),
y = rnorm(150, mu, 0.5)
)Expansions can be achieved by applying an R function that returns a matrix of basis function evaluations when applied to a vector of predictor values. We consider a spline basis expansion.
library(splines)
model.matrix(y ~ ns(x, df = 4), formula_data) |>
head() (Intercept) ns(x, df = 4)1 ns(x, df = 4)2 ns(x, df = 4)3 ns(x, df = 4)4
1 1 0.75361809 0.169569207 0.08775628 -0.04862859
2 1 0.07603605 -0.191091633 0.44504011 -0.25394848
3 1 0.77817043 0.005536704 0.13653515 -0.07790959
4 1 0.14168625 -0.207697339 0.48371373 -0.27601639
5 1 0.15040919 0.515081509 0.26234462 0.07216468
6 1 0.20011297 0.527404800 0.24317864 0.02930359Before you run the R expression below, answer the following question on the basis of what you have learned.
model.matrix(y ~ ns(x, df = 4) * z, formula_data) |>
head()Figure 7.10 was obtained by fitting the linear model using the model matrix, or rather the formula, from above. The figure illustrates the model fit by comparing the fitted model directly to the true means and via the residuals. Fitted values and residuals were computed using the augment() function from the broom package.
Exercise 7.2 Consider the following matrices
X <- splines::bs(
seq(0, 1, 0.01),
knots = c(0.3, 0.5, 0.7)
)
X_prime <- seq(0, 1, 0.01)of \(B\)-spline basis functions evaluated on a grid in the interval \([0, 1]\). Use Exercise 2.3 to verify that X_prime is in the column space of X. You can use the formula from the exercise to compute \(\mathbf{C}.\) Can you also use lm.fit()?