P-Splines with dependent data

Question

I'm working through the generalized additive models book by Simon Wood and I've had a discussion with a friend of mine over how P-Splines estimation would work for dependent data.

For independent data we have the following, staying with the notation from the book,

$y_i=f(x_i)+ \epsilon_i$ where $y_i$ is a response variable, $x_i$ is a covariate, $f$ is a smooth function and $\epsilon_i \sim N(0, \sigma^2)$. Following section 4.2.2 we have that the penalized regression fitting problem is to minimize,

$$||y-X\beta||^2+\lambda \beta^T D \beta$$ with respect to $\beta$ with $\lambda$ being the smoothing parameter, $D$ is a diagonal matrix, $\beta$ is the coefficient for $X$ and where $\beta,X$ can be a vector and matrix respectively depending on how many covariates there are. Now working through the proof to find the penalized least squares estimator of $\beta$,

$$\begin{align} ||y-X\beta||^2+\lambda \beta^T D \beta &=(y-X\beta)^T(y-X\beta)+\lambda \beta^T D \beta \\ &=(y^T- \beta^T X^T)(y-X\beta)+\lambda \beta^T D \beta \\ &=y^Ty-y^TX\beta-\beta^T X^T y+\beta^TX^TX \beta+ \lambda \beta^T D \beta \\ &=y^Ty-2y X^T \beta+\beta^TX^TX \beta+ \lambda \beta^T D \beta \\ \end{align}$$ Now taking the partial derivative with respect to $\beta$ and setting it equal to 0 we have that,

$$\begin{align} \frac{d}{d\beta}||y-X\beta||^2+\lambda \beta^T D \beta &=0-2X^Ty+2X^TX\beta+2 \lambda D \beta \\ 0&=-2X^Ty+2X^TX\beta+2 \lambda D \beta \\ 0&=-X^Ty+X^TX\beta+ \lambda D \beta \\ X^Ty&=(X^TX+ \lambda D) \beta \\ \hat{\beta} &= (X^TX+ \lambda D)^{-1} X^Ty \end{align}$$ Which is the solution for an iid case. The question now becomes how does the derivation work if we know that $\epsilon$ is not iid. So if we redefine the setup as follows,

$y=X\beta+ \epsilon$ with $R$ has the variance-covariance matrix of $\epsilon$ the penalized least squares estimator of $\beta$ now is,

$$\hat{\beta}=(X^T R^{-1}X+ \lambda D)^{-1} X^TR^{-1}y.$$ So my question is how is the solution for $\hat{\beta}$ obtained following a similar derivation above as I can see how $R$ would be used however I am unsure how the $R^{-1}$ is obtained.

Thank you for the assistance.

Answer 1

After some fiddling around with the penalization and referencing "Time Series Analysis and Its Applications: With R Examples, 4th Edition" from Robert H. Shumway & David S. Stoffer I was able to come up with a derivation that follows very similar steps to those that I worked out above. Here's a summary of the modifications:

Since the data is now correlated in some manner, for a time series case we can say that $\epsilon \sim AR(p)$ so now $R$ will be the variance-covariance matrix of $\epsilon$. Thus the penalized regression fitting problem is adjusted such that it is to minimize,

$$||R^{-.5}y-R^{-.5}X \beta||^2 + \lambda \beta^T D \beta$$

The $R^{-.5}$ is introduced into the original problem such that the original model is $R^{-.5}y=R^{-.5} X \beta+ R^{-.5} \epsilon$ if we replace the original $f(x)$ with $X \beta$. Working through the minimization there are several instances when we run into the terms of $\{ R^{-.5} \}^T R^{-.5}$ which turns into $\{ R^{-.5} \}^T R^{-.5}=R^{-1}$. Finally with that property in mind and solving for $\hat{\beta}$ the final result is achieved,

$$\hat{\beta}=(X^TR^{-1}X+ \lambda D)^{-1} X^T R^{-1} y$$