REML for Two-Dimensional P-Splines

We propose a new method based on residual (or restricted) maximum likelihood (REML) for P-splines. Existing methods use a transformation of P-splines to a mixed model; in the new model it is shown that a more direct method can be used keeping the sparse structure of P-splines. The method is illustrated with a two-dimensional example using the R-package LMMsolver on CRAN. We will show that for this example LMMsolver is several orders of magnitude faster than other methods, which use a transformation of P-splines to mixed models where the sparse structure of P-splines is lost.