A penalized framework for distributed lag non‐linear models

Distributed lag non‐linear models (DLNMs) are a modelling tool for describing potentially non‐linear and delayed dependencies. Here, we illustrate an extension of the DLNM framework through the use of penalized splines within generalized additive models (GAM). This extension offers built‐in model selection procedures and the possibility of accommodating assumptions on the shape of the lag structure through specific penalties. In addition, this framework includes, as special cases, simpler models previously proposed for linear relationships (DLMs). Alternative versions of penalized DLNMs are compared with each other and with the standard unpenalized version in a simulation study. Results show that this penalized extension to the DLNM class provides greater flexibility and improved inferential properties. The framework exploits recent theoretical developments of GAMs and is implemented using efficient routines within freely available software. Real‐data applications are illustrated through two reproducible examples in time series and survival analysis.


Introduction
Distributed lag models (DLMs), originally proposed in econometrics by Almon (1965) and more recently in epidemiology by Schwartz (2000), constitute an elegant analytical framework to describe associations characterized by a delay between an input and a response in time series data. DLMs model the response y t observed at time t in terms of past occurrences x t− of a predictor x. The new quantity , the lag, defines a new space that expresses the temporal structure of the association. In standard DLMs, a parametric function is applied to model the shape of the lag structure, usually polynomials or less often regression splines. However, the estimated shape is dependent on the chosen parametric form, for instance, the degree of the polynomial or the number and locations of the splines knots. More sophisticated smoothing techniques have been proposed to address this issue in DLMs, including penalized splines through generalized additive models (GAMs) (Zanobetti et al., 2000;Muggeo, 2008;Rushworth et al., 2013;Obermeier et al., 2015) or Bayesian approaches with the definition of prior distributions (Welty et al., 2009). While these methods offer greater flexibility and more advanced model selection procedures, they rely on the strong assumption of a linear or linear-threshold doseresponse relationship, and are only applicable to time series data.
Recent work has addressed these limitations. First, Armstrong (2006) and Gasparrini et al. (2010) extended DLMs to distributed lag non-linear models (DLNMs), a framework to describe bi-dimensional dose-lag-response associations potentially varying non-linearly in the dimensions of predictor intensity and lag. Second, Gasparrini (2014) generalized DLMs and DLNMs beyond the time series setting, extending their application to other designs and data structures. However, the current version of DLNMs still requires the user to select the parametric form of the functions expressing the dose-lag-response relationship. Model selection procedures based on information criteria have been proposed, but they lack a solid theoretical basis, and have been shown to partly affect the inferential properties of the estimators (Gasparrini, 2014;Obermeier et al., 2015).
In this contribution, we propose an extended DLNM class developed through penalized splines regression. This development provides greater flexibility for modelling potentially complex bi-dimensional dose-lag-response relationships, and offers built-in model selection and inferential procedures based on recent theoretical work on GAMs. In addition, we extend the methodology further to accommodate a priori assumptions on the shape of the lag structure through the definition of specific penalties. This general framework is applicable to model either linear or non-linear lagged relationships, in various study designs based on either time series or other data structures, and includes most of the models described above as special cases. This extension is fully implemented in freely available software routines.
The article is structured as follows: Section 2 briefly revisits the definition of DLNMs. Section 3 illustrates the extension to a penalized version of DLNMs. In Section 4, we present a simulation study for the assessment of the performance and inferential properties of both standard and extended versions. In Section 5, penalized DLNM are applied in two reproducible illustrative examples. A final discussion is provided in Section 6. Additional information and results are provided in the Web material.

The DLNM Framework
In time series data, DLMs and DLNMs model a response y t measured at time t = 1, . . . , m in terms of lagged occurrences of a predictor x t , represented by the vector q t = [x t− 0 , . . . , x t−L ] T , with 0 and L as minimum and maximum lags, respectively, (Gasparrini et al., 2010). The framework can be extended beyond the time series setting by including an additional indexing structure, allowing each response y i,t , with i = 1, . . . , n, to depend on equally spaced lagged values Here, each observation identified by the index i, for instance, a specific subject followed in time, refers to a different exposure profile that determines the lagged dose pattern at time t. Time series data represent a special case where i = t, while extensions to more complex designs such as survival or repeated-measures longitudinal data are straightforward. See Gasparrini (2014) for a more detailed explanation of the extension beyond time series data and related algebraic definitions.
The association is represented through a function s, defined as: (1) Here, the bi-dimensional dose-lag-response function f · w(x, ) is composed of two marginal functions: the standard dose-response function f (x), and the additional lag-response function w( ) that models the lag structure in the space = [ 0 , . . . , L] T . Parameterization of f and w is obtained by applying known basis transformations to the vectors q i,t and , producing marginal basis matrices R i,t and C with dimensions (L − 0 + 1) × v x and (L − 0 + 1) × v , respectively. Identifiability constraints require a reparameterization of R (see Section 3.2). The function s, here termed cross-basis function and parameterized by coefficients η, is constructed by: with w i,t as a set of known transformations derived from A i,t , which in turn is computed by a row-wise Kronecker product between the two basis matrices (Eilers et al., 2006), as: with 1 j as a vector of 1's with length j. The n × v x · v crossbasis matrix W, obtained by applying (1)-(3) to the full set of n observations, can be included in the design matrix of standard regression models, such as generalized linear models (GLMs) or Cox proportional hazard models, to estimate the parameters η.
The dose-lag-response surface can be recovered by predicting effectsβ x, on a grid of predictor values x and lag . For ease of interpretation,β x, are defined as specific contrasts of f ·w(x, ) by centering the dose-response function f (x) to a reference value of the predictor x. These effectsβ x, are interpreted in the usual scale of risk ratio or difference. In particular, the analysis commonly focuses on specific summaries, such as estimated lag-response associations at a given predictor value, or the overall dose-response association obtained by cumulating the risk across the lag period. Algebraic and interpretational details are given elsewhere (Armstrong, 2006;Gasparrini et al., 2010;Gasparrini, 2014).

Penalized DLNMs
A penalized extension of DLNM can be described within the family of GAMs (Hastie and Tibshirani, 1990;Wood, 2006a). These models extend the strong parametric form of GLMs by allowing the linear predictor to include flexible smooth functions of the covariates. A recent development of GAMs defines smooth components through penalized regression splines, using low-rank basis terms and a simple form of penalized likelihood (Wood, 2006b). This definition provides theoretically well-grounded estimators implemented using efficient and numerically stable routines (Wood, 2008(Wood, , 2011. We refer to the two versions of the method as unpenalized and penalized DLNMs, sometimes using the shortcuts GLMs and GAMs, respectively.

Penalized Likelihood
In unpenalized models, the dose-lag-response association defined by a DLNM can be estimated by maximizing the model likelihood l(η, γ) in terms of the model parameters [η T , γ T ] T , with η corresponding to coefficients of the crossbasis and γ to coefficients of additional covariate terms in the model, respectively. The idea underlying the extension to penalized DLNM is to form a richly parameterized crossbasis, and then to apply penalties through its parameters η to smooth the dose-lag-response surface. Following similar developments for tensor product bi-dimensional smoothing (Currie et al., 2004), a penalized version l p (η, γ, λ) of the model likelihood is obtained by: Here, the penalization of η is obtained through penalty matrices S x and S and penalty (or smoothing) parameters λ = [λ x , λ ] T that control the degree of smoothness of the surface. The definition in (4) offers several advantages. First, it allows different degrees of penalization along the two dimensions of the dose-lag-response function f ·w(x t− , ), by independently calibrating the smoothness in the two marginal spaces x and through λ x and λ , respectively. In addition, a mix of penalized and unpenalized functions can be defined, for example, when strong parametric assumptions can be made for either f (x) or w( ) in (1)-(3), with the exclusion of the related smoothing parameter and penalty matrix from (4). The framework proposed above, therefore, includes previously proposed penalized DLM (Zanobetti et al., 2000;Obermeier et al., 2015) by specifying a linear unpenalized f (x). Extension to models with multiple cross-basis terms or additional penalized terms are straightforward.

Choice of the Smoother
Smooth terms in a GAM can be specified by different smoothers, characterized by alternative basis functions and penalties. In penalized DLNMs, the choice of the smoother determines the basis transformations used to generate R i,t and C in (2)-(3) and the penalties that form S x and S in (4). Here, we describe two options, although others are available (Wood, 2006a).
The first smoother, labeled ps, is based on P-splines (Eilers and Marx, 1996), which offer good performance in multidimensional smoothing and both simplicity and flexibility in the penalty definition. The basis matrix of this smoother is composed of v B-splines of degree p, defined by v + p + 1 equally spaced knots. The smoothing is obtained by penalizing the difference between coefficients corresponding to adjacent splines using a difference order d. The penalty matrix is derived as: Examples of the first two orders are: A second smoother, labeled cr, is based on cubic regression splines with penalties on second derivatives. As described in Wood (2006a), for computational convenience the basis matrix of this smoother is derived using a special parameterization of v natural cubic splines, where only one spline is non-zero at each of the v knots. These knots can be placed everywhere along the range of the predictor, and by default at equally spaced quantiles. The smoothing is induced by penalizing the second derivative of the function, with more complex computation required to derive the penalty matrix S (Green and Silverman, 1994).
The tensor product form of the cross-basis in (1)-(3) requires constraints to ensure the identifiability of the regression parameters. Specifically, the identifiability constraints are absorbed into R and S x through a reparameterization. This step coincides with the exclusion of the intercept from f (x) in unpenalized DLNMs, and ensures that the n × v x · v crossbasis matrix W can be full rank. Additional details are given in Wood (2006a) and the references cited above.

Alternative Penalties on the Lag Structure
Specific assumptions can be made about the shape of the relationship in the lag dimension. These assumptions can be incorporated through additional penalties, which fall into two main categories. First, varying ridge penalties can be imposed to shrink different parts of the lag-response curve towards the null value. These type of penalties can be used with either ps or cr smoothers, and takes two alternative forms: Here, P is a pre-specified diagonal matrix of weights p, which in (7a) are applied directly to the v coefficients η, while in (7b) are chosen for the L − 0 + 1 lags and mapped into η through the basis matrix C defined in (3). These were previously discussed in Muggeo (2008) and Obermeier et al. (2015).
A second type are varying difference penalties that can be applied to enforce a different degree of smoothness along the lag-response curve. These penalties are naturally defined for ps smoothers, and while technically applicable with cr smoothers as well, they are less theoretically grounded for the latter. They take the forms: where P defines weights p for v − d and L − 0 + 1 − d differences between coefficients in (8a) and lags in (8b), respectively, while D d is a matrix defined in (6) of consistent dimensions. Single or multiple penalties for the lag structure as in (5) or (7)-(8) can be imposed in the same model by defining for each of them the smoothing parameters λ and penalty matrices S in (4). See Sections 4-5 for specific examples.

Estimation
After the model has been defined by the choice of basis terms and penalty matrices, maximization of the penalized log-likelihood l p (η, γ, λ) in (4) is solved through standard estimation methods for GAMs (Wood, 2006a). Briefly, a penalized iterative reweighted least square (P-IRLS) method is integrated with multiple smoothing parameter selection to estimate the degree of smoothness. Alternative methods are available for the estimation of the smoothing parameters λ within P-IRLS, such as generalized cross validation (GCV), unbiased risk estimator (UBRE, essentially scaled AIC) and (restricted) maximum likelihood (REML and ML), all implemented using reliable and computationally efficient routines (Wood, 2008(Wood, , 2011. Simulations indicate that REML and ML are superior in terms of mean-square error performance and smoothing properties (Wood, 2011).
Approximate point-wise confidence intervals of the doselag-response surface and its summaries are computed from the estimated posterior (co)variance matrix of the coefficientŝ η, derived using empirical Bayesian estimators (Marra and Wood, 2012). These account for the inherent bias affecting smooth terms and have been shown to provide confidence intervals with across-the-function frequentist coverage close to nominal. Although the estimators applied here neglect the uncertainty in the estimation of the smoothing parameters λ, this has little effect on interval performance in real-data settings (Marra and Wood, 2012). The smoothness of the dose-lag-response surface can be quantified in terms of effective degrees of freedom (edf), with the high boundary usually represented by the product of the dimensions of the two marginal basis matrices, v x × v (when λ x = λ = 0), and the lower boundary determined by the product of their null space dimensions (when λ x → ∞ and λ → ∞). The null space dimension of each marginal basis is equal to the order of the penalty, that is, the difference order d in ps or the order of derivative (usually 2) in cr smoothers, respectively, minus any constraint (Wood, 2006a).

Simulation Study
To assess the performance and inferential properties of different versions of penalized DLNMs and to compare them with the standard unpenalized approach, we performed a simulation study based on scenarios of dose-lag-response surfaces with varying degree of complexity.

Simulation Setting
The predictor x t was represented by the daily temperature series in Chicago within the period 1987-2000 (Samet et al., 2000), standardized over the range 0-10. For each replicate, we simulated an outcome series y t of daily mortality counts, with t = 1, . . . , 5114, from a Poisson distribution with mean μ t , using: We repeated the simulations in three scenarios j = 1, 2, 3, with the dose-lag-response function f j ·w j (x, ) over lag 0-40 described by: r Scenario 1: a simple plane; r Scenario 2: a shape resembling previously estimated temperature-mortality associations; r Scenario 3: a complex wiggly surface.
A graphical representation of these three scenarios is offered in Figure 1, with algebraic details provided in Web Appendix A. The intercept α j was used as a signal-to-noise parameter to achieve a Pearson correlation coefficient between μ t and y t of approximately 0.5 in each scenario.
For each simulated series, we fitted alternative models where the second term in (9) is replaced by a cross-basis s(x t , . . . , x t−40 ). The primary model, simply labeled gam, used a penalized DLNM with ps smoothers of rank v = 10 (minus constraints), degree 3 (cubic B-splines) and second-order difference (d = 2) penalties for each marginal dimension, estimated by REML. Previous research (Wood, 2006a) indicates that basis dimension is not critical, if large enough to fit the underlying marginal shape, while the smoother and estimator were chosen for their flexibility and inferential performance, respectively. This model was compared with: • Alternative estimators: glm-aic, defined by (unpenalized) quadratic B-spline functions with the optimal number of equally spaced knots selected by minimizing AIC among combinations producing 1-10 df (minus constraints) in each dimension; gam-aic, fitted by replacing REML with a UBRE-AIC estimator.
r Alternative smoothers: gam-cr, defined by replacing the ps with cr smoothers; gam-ps 2,1 , with difference penalties of order 2 and 1 for f (x) and w( ), respectively.
These additional/alternative penalties follow the assumption of a lag-response that approaches the null value at the end of the lag period, or that is smoother at longer lags. See Muggeo (2008) and Obermeier et al. (2015) for details.
We assessed the performance of the eight models above using 1000 simulation replicates in each of the three scenarios depicted in Figure 1, by comparing the across-the-surface cov-  (2012)

Results of the Simulation Study
Results are illustrated in Table 1 and Figure 2. Table 1 reports the average computing time and edf, the empirical coverage of 95% confidence intervals and the empirical RMSE relative to the gam model across the surfaces. Figure 2 displays the estimated dose-response and lag-response curves corresponding to the bold black lines across the surfaces in Figure 1 for the models gam, glm-aic, and gam-add last . The same graphical representation for the other models is provided in Web Figures S2-S3 in the online Supporting Information. In all the scenarios, penalized DLNMs appear superior to the unpenalized counterpart. In particular, glm-aic shows higher RMSE (as also suggested by the wigglier curves in Figure 2), and a substantial under-coverage due to unaccounted additional variability of the model selection procedure, consistently with what previously reported (Sylvestre and Abrahamowicz, 2009;Gasparrini, 2014). The REML estimator exhibits a slightly better performance if compared to UBRE-AIC in gam-aic, with the latter showing higher RMSE and some evidence of undersmoothing, especially in the simplest scenario. Alternative smoothers in gam-cr and gam-ps 2,1 provide similar outputs, with the latter performing better in the plane scenario, which is consistent with its null space of 1 edf.
The model gam-add last shows an improved performance in the second scenario, where the extended flat region (see Figure 2) is well fitted through the addition of a varying ridge penalty, which also helps identify the correct lag period even when the interval is extended well beyond it, as previously reported (Obermeier et al., 2015). This doubly penalized model performs well also in the other scenarios that do not match the assumption of the penalty, with only a minor bias produced in the plane scenario, as noticeable in the last part of the estimated lag-response curves in Figures 2 and S3. This good performance is due to the possibility of virtually removing the additional penalty by estimating a very low smoothing parameter λ . Models gam-alt quad and gam-alt exp , where the standard penalty was removed, perform well in the second and third scenarios, but the latter fails to fit the plane dose-lag-response surface, which is not compatible with its strong assumptions about form of the lag-response shape (see Web Figure S3).
Generally, penalized models show across-the-surface coverage close to the nominal value, although some undercoverage is evident for some models in the second scenario (see also Web Figures S4-S5 in the online Supporting Information). In addition, gam-aic fails to converge in 1.4% of replicates of the plane scenario, where the simulated surface represents the null space dimension of the tensor product smoother. However, the analysis of non-convergent models does not identify problems with point estimates and coverage.

Two Examples
As an illustration of the application of penalized DLNMs in different study designs, we replicate two published analyses. The reader can refer to the original publications for details on the analytical methods and data (Gasparrini and Leone, 2014;Gasparrini, 2014).

Outdoor Temperature and All-Cause Mortality
The first example illustrates the application of penalized DLNMs in time series data, using daily series from London in the period 1993-2006. Specifically, the relationship between counts of all-cause mortality y t at day t and outdoor temperature x t− , accounting for up to 25 days of lag, was estimated with a quasi-Poisson GLM of form: log[E(y t )] = α + s(x t , . . . , x t−25 ; η) + g(t; γ) + 6 j=1 δ j w j,t , (10) with g as natural cubic splined defined by 10 df/year accounting for seasonal and long term trends, and w j as an indicator of day of the week. In the original analysis (Gasparrini and Table 1) in 1000 replicates. The panels represent the dose-response (rows 1-3) and lag-response curves (rows 4-6) corresponding to the bold black lines in the three simulated surfaces (by column) in Figure 1. Continuous gray, and dashed red and continuous black lines represent the fit from 25 random replicates, the average across all replicates, and the true simulated curves, respectively.
Leone, 2014), the dependency was modelled with an unpenalized DLNM within a GLM, using a cross-basis function s with 4 × 5 = 20 df composed of quadratic B-splines defined by 2 equally spaced internal knots for the dose-response function f (x) and natural cubic splines by three equally spaced internal knots in the log scale plus intercept for the lag-response function w( ). Boundary knots were placed by default at the ranges. We replicated the analysis using a penalized DLNM within a GAM with a REML estimator, specifying marginal ps smoothers with dimension 10 (minus constraints) for both spaces. Penalization of f (x) was enforced through a default second-order difference penalty as in (5). Extending previous models (Muggeo, 2008;Obermeier et al., 2015), we applied a double varying penalty to w( ) using a second-order difference form (8b) with p k = k 2 for k = 0, . . . , 23, and a ridge penalty of form (7a) with p = [0, 0, 0, 0, 0, 0, 1, 1, 1, 1] T . These choices are motivated by the assumption of a shape that is smoother at longer lags and approaches the null value at the end of the lag period.
The GAM used 35.45 edf to model the dose-lag-response surface, and suggests a strong and short-term association with heat and a more delayed association with cold temperatures, consistently with previous results . The estimates, reported in the first row of Figure 3, are very similar to those from the original analysis, replicated in the second row. However, it is interesting to note the effect of the double varying penalty in the estimated lag-response at 29 • C, with the curve shrinking toward the null at lags higher than 15. In addition, while the cross-basis specification of the unpenalized DLNM was originally defined a priori, an AIC-based selection suggests a very complex and implausible model with 10 × 9 = 90 df, with estimates illustrated in the third row of Figure 3.
As previously mentioned (Section 3.1), the flexibility of this modelling framework allows a mix of penalized and unpenalized functions. As an example, we replaced the ps smoother for f (x) with an unpenalized double-threshold function, that is, linear splines which model a straight relationship below 17 • C and above 21 • C, and a flat region in between. Results are displayed in the last row of Figure 3. This model uses only 10.64 edf to define the dose-lag-response surface, although this comes at the price of making strong parametric assumptions for one of the two spaces. The same approach can be used to specify simpler penalized DLMs, by selecting a linear function as f (x).

Occupational Radon and Lung Cancer Mortality
The second example describes the extension of penalized DLNMs to individual time-to-event data, using a cohort of 3347 miners working in the Colorado Plateau mines, with follow-up at December 31, 1982. Specifically, the association between an indicator of occurrence of lung cancer death y i,t for subject i at age t, and yearly occupational radon exposure x i,t− , measured in working-level months (WLM), with lags of 2-40 years, was estimated with a Poisson GLM of form: This GLM approximates the Cox proportional hazard model applied in the original analysis (Gasparrini, 2014) by splitting the follow-up time of each individual into 1-year periods, and modelling the baseline risk with a cubic B-spline function g(t) with 5 df. This allows the use of penalized splines implemented within GAMs with survival data. Other terms in the model are a cross-basis function s z to control for the lagged effect of smoking z, and a linear term for calendar year c. In the original analysis, the association with radon was modelled with a cross-basis function s x composed of quadratic B-splines functions with a single internal knot at 59.4 WLM/year and 13.3 years of lag for f (x) and w( ), respectively, and boundary knots at the respective ranges. The intercept was excluded from the latter, assuming no effect for exposures experienced within the first two years. This model, using a total of 9 df to define the association, was selected by minimizing AIC.
The analysis was replicated with a penalized DLNM using a GAM with a REML estimator, using marginal cr smoothers with dimension 11 (minus constraints) and 10 for exposure and lag spaces, respectively. The use of the cr smoother allows placing the knots of the dose-response function f (x) at equally spaced intervals in the log scale, accounting for the highly skewed distribution of radon exposure, and allows excluding the intercept in s( ) following previous assumptions. In addition to the default penalty on the second derivative, enforced in both spaces, we added a varying ridge penalty of form (7b) to w( ) with p = 1 if > 30 and 0 otherwise, thus assuming no additional risk 30 years after the exposure to radon.
Results are displayed in Figure 4. The penalized DLNM (first row) indicates a peak in lung cancer mortality risk approximately 11 years after the exposure to radon. The nonlinear dose-response shows how the risk flattens out above 50 WLM/year. The model used a total of 8.03 edf to describe the association. These findings are consistent with the unpenalized DLNM fitted with a GLM (second row of Figure 4), which closely approximates the original estimates from the Cox model illustrated by Gasparrini (2014, Figure 2). However, the addition of the ridge penalty in the GAM produces more precise estimates at the end of the lag period, suggesting that the risk completely disappears 30 years after the exposure.

Discussion
In this contribution, we describe a penalized framework for DLNMs that provides significant developments to this modelling class, through built-in smoothness selection of potentially complex marginal functions and the flexible definition of penalties to accommodate assumptions on the lag structure. This method includes previous smoothing approaches for simpler DLMs (Zanobetti et al., 2000;Obermeier et al., 2015) as special cases, and fully extends the penalized approaches to bi-dimensional dose-lag-response surfaces. The DLNM framework unifies methods proposed to investigate lagged associations in different research fields, beyond time series analysis in environmental research. For instance, these include case-control studies in cancer epidemiology (Thomas, 1988;Hauptmann et al., 2000;Berhane et al., 2008;Richardson, 2009)   summarizing the association between occupational radon exposure and lung cancer mortality, estimated by a GAM with an additional varying ridge penalty in the lag space and a GLM with AIC-based selection (as in Gasparrini (2014)) (by row). Colorado Plateau Uranium miners cohort, follow-up at December 31, 1982.

and survival analysis in
pharmaco-epidemiology (Sylvestre and Abrahamowicz, 2009;Abrahamowicz et al., 2012) This penalized version addresses the problem of choosing the appropriate degree of complexity of the DLNM. This is a critical limitation of traditional unpenalized DLNMs, for which current selection methods are not effective (as demonstrated with the first example in Section 5.1), and produce less efficient estimators (as illustrated in the simulation study in Section 4). This penalized extension is based on well-grounded theoretical results and estimation methods, recently discussed (Wood, 2006a(Wood, , 2008(Wood, , 2011, it can be performed with stable and efficient routines implemented in freely available software (Wood, 2006a), and it shows improved inferential properties if compared to the standard unpenalized version.
The results confirm the good inferential properties of REML and UBRE-AIC estimators, with the former appearing slightly superior (Wood, 2011), and the similar performance of alternative types of smoothers (Wood, 2006a). The latter can be selected due to convenient characteristics, such as the possibility of including varying difference penalties with the ps smoother (see Section 5.1) or the flexibility in the knots placement and exclusion of intercept with the cr smoother (see Section 5.2). In particular, the inclusion of additional penalties on the lag dimension provides a way to accommodate realistic assumptions on the underlying shape. These addi-tional penalties can be selected based on prior knowledge, and do not represent strong constraints on the lag-response shape, as their influence can be calibrated through the estimate of smoothing parameters. As previously suggested (Muggeo, 2008;Obermeier et al., 2015) and shown in the second scenario of the simulation study, additional penalties can improve the model fit and make the model less sensitive to the choice of the lag period.
Some limitations must be acknowledged. The issue of penalizing complex bi-dimensional functions has been investigated in a limited set of simulated scenarios and two real-data examples. Also, simulations show some issue with nonconvergence in the simplest scenario, where the selected edf tend to be close the null space of the cross-basis function, although this problem does not seem to critically affect inference. The penalized approach substantially improves the coverage properties of the confidence intervals, even though in some scenarios and for some models the empirical coverage falls short of the nominal value. In addition, the method presented here shares a known limitation of GAMs, which tend to select simpler (i.e., smoother) models when the statistical power decreases. Finally, smoothing methods for dose-lag-response relationships are difficult to validate, as the lag dimension is not directly observable in the data, thus preventing the use of standard techniques such as residual analysis. These issues will be hopefully addressed in future research.
Penalized DLNMs can be further extended to varyingcoefficients models, describing dose-lag-response relationships that change in time or along the space of other predictors, as previously described for simpler penalized DLMs (Rushworth et al., 2013) or unpenalized DLNMs . In addition, the DLNM class has interesting links with penalized functional regression, where a functional outcome (say the shape of the dose-response) is allowed to vary depending on a functional predictor (say the lag dimension) (McLean et al., 2014;Scheipl et al., 2015), providing the possibility of further extensions through this established modelling framework.
Lagged associations occur almost universally in biomedical research, and well beyond. Penalized DLNMs offer a flexible modelling class to describe these phenomena, avoiding biases due to incorrect assumptions about the lag structure, sometimes made when using simpler approaches, and potentially extending the knowledge of the association under study. The recent extension of DLNMs beyond time series data (Gasparrini, 2014) unifies and extends methods proposed in different study designs and paves the way for original and promising applications of this modelling framework.

Supplementary Materials
Web Appendices, Web Figures, and R code are available at the Biometrics website on Wiley OnlineLibrary. In addition to Web Appendices A-B, referenced in Section 4, Web Appendix C briefly describes the software implementation in the R package dlnm. The R code fully reproduces the simulation studies and the two examples, with an updated version available from GitHub and the personal website of the first author (see Web Appendix C).