Paper deep dive
Distributed Lag Neural Additive Models
Calle Helmersson, Shivang Pandey, Leonardo Olivetti, Elena Raffetti
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 93%
Last extracted: 9/13/2026, 5:18:11 AM
Summary
The paper introduces Distributed Lag Neural Additive Models (DLNAMs), a neural-network-based extension of Distributed Lag Non-linear Models (DLNMs) designed to learn nonlinear exposure-lag response surfaces. DLNAMs replace prespecified spline cross-bases with learned neural components, utilizing exp-centered input layers, smooth activations, and subnetwork mixtures to achieve local adaptivity and smoothness while maintaining additive interpretability. The study demonstrates through simulations and real-world applications (Chicago temperature-mortality and multi-country malaria data) that DLNAMs generally outperform traditional DLNM variants (QAIC, QBIC, P-DLNM, T-DLNM) in recovering response functions, particularly in complex, nonseparable scenarios, with better calibration and lower bias.
Entities (9)
Relation Signals (7)
DLNAM → appliedto → Chicago NMMAPS
confidence 95% · Applications to Chicago temperature–mortality data... test whether the same component architecture can be reused
DLNAM → appliedto → DHS/MIS Malaria
confidence 95% · DLNAM was applied to 348,565 childhood malaria outcomes from Demographic and Health Surveys (DHS) and Malaria Indicator Surveys (MIS)
DLNAM → isananalogueof → DLNM
confidence 95% · We introduce Distributed Lag Neural Additive Models (DLNAMs), neural-additive analogues of Distributed Lag Non-linear Models (DLNMs)
Shivang Pandey → iscorrespondingauthorof → DLNAM
confidence 90% · Corresponding author: Shivang Pandey... We introduce Distributed Lag Neural Additive Models (DLNAMs)
DLNAM → outperforms → T-DLNM
confidence 90% · In simulations, DLNAMs generally outperformed DLNM comparators, including penalized and treed variants
DLNAM → outperforms → DLNM
confidence 90% · In simulations, DLNAMs generally outperformed DLNM comparators... in recovering known response functions
DLNAM → usescomponent → ExU
confidence 90% · The DLNAM architecture combines exp-centered unit (ExU) input layers [12] with smooth activations
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:We introduce Distributed Lag Neural Additive Models (DLNAMs), neural-additive analogues of Distributed Lag Non-linear Models (DLNMs) for learning nonlinear effects distributed over lags. DLNAMs replace a prespecified spline cross-basis with neural components that learn exposure--lag response surfaces, avoiding choices of basis family, dimension, and knot placement while preserving additive interpretability and familiar distributed-lag summaries. Exp-centered input layers, smooth activations, and learned subnetwork mixtures produce smooth, locally adaptive representations; pointwise uncertainty combines a conditional last-layer Laplace approximation with between-member ensemble variation. In simulations, DLNAMs generally outperformed DLNM comparators, including penalized and treed variants, in recovering known response functions, with lower bias, stronger boundary recovery, and better-calibrated cumulative intervals; gains were largest for more demanding functions. The architecture performed consistently across sample sizes, outcome families, lag horizons, and jointly fitted multi-exposure settings, retaining recovery performance as exposures were added; fit-specific changes were largely confined to optimization, and applications recovered established empirical patterns.
Tags
Links
- Source: https://arxiv.org/abs/2609.07381v1
- Canonical: https://arxiv.org/abs/2609.07381v1
Trouble viewing inline? Open PDF directly →
Full Text
163,325 characters extracted from source content.
Expand or collapse full text
Distributed Lag Neural Additive Models Calle Helmersson Affiliation: Department of Mathematics, KTH Royal Institute of Technology, Stockholm, Sweden Affiliation: Department of Global Public Health, Karolinska Institutet, Stockholm, Sweden Shivang Pandey Affiliation: Department of Global Public Health, Karolinska Institutet, Stockholm, Sweden Affiliation: Swedish Centre for Impacts of Climate Extremes (climes), Uppsala University, Uppsala, Sweden Affiliation: Corresponding author: Shivang Pandey, shivang.pandey@ki.se Leonardo Olivetti Affiliation: Swedish Centre for Impacts of Climate Extremes (climes), Uppsala University, Uppsala, Sweden Affiliation: Department of Earth Sciences, Uppsala University, Uppsala, Sweden Affiliation: Centre of Natural Hazards and Disaster Science (CNDS), Uppsala University, Uppsala, Sweden Elena Raffetti Affiliation: Department of Global Public Health, Karolinska Institutet, Stockholm, Sweden Affiliation: Swedish Centre for Impacts of Climate Extremes (climes), Uppsala University, Uppsala, Sweden Affiliation: Department of Earth Sciences, Uppsala University, Uppsala, Sweden Affiliation: British Heart Foundation Cardiovascular Epidemiology Unit, Department of Public Health and Primary Care, University of Cambridge, Cambridge, UK Affiliation: Victor Phillip Dahdaleh Heart and Lung Research Institute, University of Cambridge, Cambridge, UK Abstract We introduce Distributed Lag Neural Additive Models (DLNAMs), neural-additive analogues of Distributed Lag Non-linear Models (DLNMs) for learning nonlinear effects distributed over lags. DLNAMs replace a prespecified spline cross-basis with neural components that learn exposure–lag response surfaces, avoiding choices of basis family, dimension, and knot placement while preserving additive interpretability and familiar distributed-lag summaries. Exp-centered input layers, smooth activations, and learned subnetwork mixtures produce smooth, locally adaptive representations; pointwise uncertainty combines a conditional last-layer Laplace approximation with between-member ensemble variation. In simulations, DLNAMs generally outperformed DLNM comparators, including penalized and treed variants, in recovering known response functions, with lower bias, stronger boundary recovery, and better-calibrated cumulative intervals; gains were largest for more demanding functions. The architecture performed consistently across sample sizes, outcome families, lag horizons, and jointly fitted multi-exposure settings, retaining recovery performance as exposures were added; fit-specific changes were largely confined to optimization, and applications recovered established empirical patterns. 1 Introduction When an outcome depends on a history of past exposure rather than on its current value alone, the scientific object is a response function over both exposure level and elapsed time. Such distributed-lag relationships arise whenever a system responds gradually, transiently, or with delay. Formalized early in econometrics [1], they are widely studied in environmental epidemiology, where the timing of temperature, pollution, and other exposures can matter alongside their magnitude [2, 3]. Distributed Lag Non-linear Models (DLNMs) provide a widely used framework for modeling nonlinear effects distributed over past exposures within an additive, interpretable structure [4, 5, 6]. The DLNM framework represents exposure–lag relationships through smooth additive functions, most commonly implemented in environmental epidemiology using spline cross-bases [4, 5, 6]. These combine marginal basis functions over exposure and lag to construct a smooth response surface. Their flexibility is therefore shaped by choices such as basis family, degrees of freedom, and knot placement. These choices may be guided by substantive expertise or selected from candidate specifications using quasi-likelihood information criteria (QAIC or QBIC). P-DLNM reduces sensitivity to smoothness specification through penalization, but local variation is still expressed within a prescribed spline representation [7]. T-DLNM instead learns local partitions of the exposure–lag domain, gaining substantially greater local adaptivity but at the cost of a piecewise-constant response [8]. Each approach addresses a different facet of the representation problem, yet smoothness, local adaptivity, and analysis-specific representation remain in tension. Neither combines a smooth, locally adaptive response with a representation whose internal features are learned during estimation. We introduce Distributed Lag Neural Additive Models (DLNAMs), in which each lagged predictor contributes a learned neural component within an additive predictor. DLNAMs replace prespecified spline representations with learned smooth functions and retain the principal scientific objects through which distributed-lag models are interpreted: exposure–lag surfaces, cumulative response functions, and covariate effects. Additional lagged predictors enter as further learned components, with the core architecture generally reusable as the model expands rather than requiring the response representation to be redesigned for each added predictor. Prediction-oriented work has used distributed-lag exposure histories as inputs to generic artificial neural networks [9]; DLNAMs instead estimate a separate exposure–lag response surface for each exposure, from which the usual distributed-lag summaries are obtained. The organizing principle is that additivity and functional representation are separable: the additive predictor makes each component separately reportable, irrespective of the function class used to represent it. This separation has an additive-model precedent. Generalized Linear Models (GLMs) use linear components [10], Generalized Additive Models (GAMs) replace them with smooth basis-expanded functions [11], and Neural Additive Models (NAMs) learn the component functions with neural networks [12]. In the cross-sectional setting, NAMs achieve accuracy comparable to state-of-the-art GAMs while retaining additive interpretability and the architectural flexibility of neural networks [12]. The distributed-lag lineage runs in parallel, from Distributed Lag Models (DLMs) to DLNMs [1, 4]; DLNAMs complete it with the neural-component analogue. Cross-sectional:GLM⟶GAM⟶NAM,Distributed lag:DLM⟶DLNM⟶DLNAM. array[]l@ cCross-sectional:&GLM& &GAM& &NAM,\\ Distributed lag:&DLM& &DLNM& &DLNAM. array Each step replaces the functional representation of the components while leaving the additive structure intact. Constructing that analogue is not straightforward. A generic neural substitution is insufficient: a useful component must reconcile local adaptivity with smoothness, stable estimation, additive reporting, and uncertainty quantification. The DLNAM architecture combines exp-centered unit (ExU) input layers [12] with smooth activations [13, 14] and a learned mixture of subnetworks. Because exposure and lag enter a surface component jointly, the scalar ExU construction of NAMs is extended to multiple input coordinates; alternative multivariate constructions are compared directly. Complete fitted models are ensembled across independent initializations [15]. Pointwise uncertainty combines conditional last-layer Laplace approximations with between-member ensemble variation [16, 17]. Together, these design choices make DLNAMs structured neural estimators rather than generic feed-forward substitutes for splines. We evaluate DLNAMs by how well they recover the exposure–lag response functions of scientific interest. Simulation studies assess recovery of cumulative responses and full exposure–lag surfaces, interval calibration, architectural ablations, and additive component recovery under correlated concurrent exposures. Applications to Chicago temperature–mortality data and a large multi-country malaria analysis test whether the same component architecture can be reused across markedly different sample sizes, outcome families, and lag structures. We additionally quantify recurring representation-selection burden and fixed-budget computational scaling as concurrent lagged components are added. Taken together, these experiments test whether representation learning can augment distributed-lag modeling without surrendering the scientific structure that makes its effects interpretable. 2 Results The principal simulation comprised one separable data-generating process (DGP) and three structurally demanding, nonseparable DGPs with localized or heterogeneous exposure–lag structure (Fig. S1). We compare DLNAM with four DLNM-family estimators: DLNM (QAIC), DLNM (QBIC), P-DLNM, and T-DLNM (Section 4.7). Across Monte Carlo replicates, we evaluated recovery of the cumulative response and full exposure–lag surface, interval calibration, and performance in sparsely supported exposure tails (Section 4.6). Here, recovery denotes agreement between an estimated target and its known data-generating counterpart. 2.1 Recovery Across DGPs 2–4, DLNAM attains the lowest error of the five estimators over the whole exposure grid, the interior, and the boundary region (Fig. 1B; Table S4). On the separable DGP 1, T-DLNM attains the lowest whole-grid and boundary RMSE, while DLNAM is lowest in the interior; T-DLNM’s advantage is driven by lower variance rather than lower bias. DLNAM also has the smallest across-DGP range in whole-grid RMSE. Its gains concentrate on DGPs 2–4, consistent with representation learning being most useful for demanding local or heterogeneous structure, although the design does not isolate nonseparability, localization, or any single dimension of complexity as the cause. Figure 1: Simulation Study: Model Comparison. Monte Carlo simulation of four DGPs, with R=200R=200 replicates per DGP. (A) Estimated cumulative relative-risk functions and data-generating functions, shown as pointwise Monte Carlo means across fitted replicates. (B) Cumulative log-relative-risk RMSE by method and exposure region, partitioned into squared bias and variance. (C) Empirical coverage of pointwise 95%95\% intervals. Error bars denote 95% Monte Carlo confidence intervals. Boundary denotes exposure values below the 5th or above the 95th percentile of the reference exposure distribution; interior denotes its complement. Associated numerical results are reported in Tables S4 and S5. The broad pattern persists when recovery is scored on the full exposure–lag surface—a target emphasized in both comparator studies [7, 8]—with the treed fit remaining competitive at the boundary (Fig. S4; Table S10). The cumulative gains therefore do not arise from lag-specific errors cancelling in the sum. 2.2 Bias and Variance The bias–variance decomposition exposes distinct failure modes (Fig. 1B; Table S5). DLNAM has the lowest whole-grid squared bias in DGPs 2–4, while T-DLNM has the lowest whole-grid variance in DGPs 1–2. On DGPs 2–4, rising comparator error is driven chiefly by squared bias; DLNAM’s residual cumulative error is predominantly variance on every DGP. The same contrast appears in the sparsely supported exposure tails: fitted cross-bases can oscillate there, whereas the learned components tend to flatten. Oscillation can displace the fitted response; flattening instead tends to attenuate a true tail effect. 2.3 Calibration No estimator attains nominal coverage of pointwise 95%95\% intervals throughout the design (Fig. 1C). On the cumulative curve, the DLNAM is closest to nominal coverage in 11 of the 12 DGP-by-region cells; the exception is the DGP 1 boundary, where T-DLNM is closer to nominal (Table S4). Coverage and interval width must be read jointly because coverage can be raised mechanically by widening intervals (Table S6). On the full exposure–lag surface, every estimator with higher coverage than DLNAM also has wider intervals; on the cumulative curve the same holds apart from two DGP-by-region cells. T-DLNM makes the trade-off most visible on the surface target, where it combines the highest coverage with the widest intervals and larger error. DLNAM’s calibration advantage therefore does not come from uniformly wider intervals. Adding between-member ensemble variation to the last-layer Laplace variance raises DLNAM coverage in every DGP for both cumulative and full-surface targets, with the Laplace-only result shown in Table S7. Fit-to-fit instability thus carries material information for interval calibration. 2.4 Architecture Ablation The ablation isolates which architectural elements matter for recovery (Fig. 2; Table S9). The reference specification denotes the complete DLNAM architecture used in the principal analysis. Replacing the ExU input layer with a linear input layer produces the largest increase in error across all four DGPs and sharply reduces coverage, with the largest degradation on DGPs 2–4. The learned input scales therefore play a role analogous to marginal resolution in a cross-basis, except that they are fitted parameters rather than a dimension declared in advance. Figure 2: Simulation Study: Architecture Ablation. Monte Carlo simulation of four DGPs, with R=50R=50 replicates per DGP. (A) Estimated cumulative relative-risk functions and data-generating functions, shown as pointwise Monte Carlo means across fitted replicates. (B) Cumulative log-relative-risk RMSE by method and exposure region, partitioned into squared bias and variance. (C) Empirical coverage of pointwise 95%95\% intervals. Error bars denote 95% Monte Carlo confidence intervals. Boundary denotes exposure values below the 5th or above the 95th percentile of the reference exposure distribution; interior denotes its complement. Ablations are applied one at a time and are not factorial. Associated numerical results are reported in Table S9. The remaining ablations expose different failure modes. Replacing the mixture with the standard single-network formulation increases error and depresses coverage across all four DGPs; replacing the smooth activations increases error but leaves coverage at or above the full model’s. The input-layer and mixture ablations are most costly on DGPs 2–4, whereas the smooth-activation ablation costs most when the generating surface itself varies smoothly. These patterns indicate that recovery is not explained by generic neural-network capacity alone. ExU supplies local sensitivity, the mixture composes learned functions, and the smooth hidden activations govern how that flexibility is expressed. Because the ablations are one-at-a-time, they do not support factorial attribution of independent effects. 2.5 Joint-Exposure Recovery With four active exposure–lag components fitted concurrently alongside a fifth correlated null exposure, DLNAM achieved the lowest cumulative error of the five estimators for every active component over the whole exposure grid, the interior, and the boundary region (Fig. 3A, B). Whole-grid degradation ratios remained near 1 across all four active components relative to the corresponding single-exposure experiments (Fig. 3D). Region-specific ratios were somewhat higher in the interior and lower at the boundary, indicating redistribution of error across the exposure domain rather than a uniform joint-fitting penalty. Because T-DLNM is target-specific in this experiment, its larger degradation pertains to the practical multi-exposure formulation used here, not to a hypothetical fully joint treed model. Coverage of the active components remained below nominal, but the DLNAM was closest to nominal in all 12 active component-by-region cells (Fig. 3C; Table S15). For the correlated null exposure, DLNAM and DLNM (QBIC) showed the least leakage by whole-grid RMSE (0.0124 and 0.0123, respectively), while DLNAM had coverage closest to nominal (0.955). Figure 3: Simulation Study: Model Comparison (Joint). Monte Carlo simulation jointly fitting the four DGPs as concurrent exposure–lag components, with R=50R=50 replicates, alongside one correlated null exposure. (A) Estimated cumulative relative-risk functions and data-generating functions for the four active components, shown as pointwise Monte Carlo means across fitted replicates. (B) Cumulative log-relative-risk RMSE by method and exposure region, partitioned into squared bias and variance. (C) Empirical coverage of pointwise 95%95\% intervals. (D) RMSE degradation relative to the corresponding single-exposure fits. Error bars denote 95% Monte Carlo confidence intervals. Boundary denotes exposure values below the 5th or above the 95th percentile of the reference exposure distribution; interior denotes its complement. Associated numerical results, including the correlated null exposure, are reported in Table S15. Together, these results support stable recovery of additive cumulative components under the specified correlation structure without establishing robustness to arbitrary collinearity or exposure–exposure interactions. 2.6 Chicago NMMAPS Temperature and Mortality Using 4638 complete days from the Chicago component of the National Morbidity, Mortality, and Air Pollution Study (NMMAPS) and a 30-day lag window, the DLNAM recovers the familiar temperature–mortality pattern [6] (Fig. 4): a minimum-mortality temperature with risk rising toward both extremes, heat acting within days, and cold accumulating over one to two weeks. All five estimators agree on the broad cumulative shape but place the minimum across a range of several degrees, showing the fitted optimum’s sensitivity to specification. They nonetheless differ in how the heat association is distributed across lags, with the DLNAM surface spreading a lower elevation across more of the early lags. The Chicago NMMAPS data are the worked example distributed with the dlnm package [6], and the published analysis supplies an empirical reference under a different representation and adjustment set. DLNAM produced a similar overall shape, minimum location, and heat-response magnitude without explicit cross-basis selection. In the sparse cold tail, its curve remains within the envelope of the likelihood-matched fits while its intervals are the narrowest among them. This observational, in-sample application tests whether the architecture can recover a familiar empirical pattern; it is not intended to establish causal, predictive, or benchmark superiority. Figure 4: Chicago NMMAPS: Model Comparison. Cumulative relative-risk functions and corresponding temperature–lag response surfaces are shown for the fitted methods. Shaded bands around the cumulative functions denote pointwise 95%95\% intervals. All methods use the same adjustment variables. Analysis specifications are reported in Table S3. 2.7 Multi-Country DHS/MIS Malaria DLNAM was applied to 348,565 childhood malaria outcomes from Demographic and Health Surveys (DHS) and Malaria Indicator Surveys (MIS) across 26 sub-Saharan African countries, under a logit link over 6 monthly lags (Fig. 5). One model was fitted per reported exposure using the adjustment set specified for that target in the source analysis [18]. Within each fit, the adjusting exposures retain full exposure–lag surfaces; the specific-humidity model therefore contains 5 concurrently lagged surfaces. A DLNM is fitted under the same adjustment sets, with adjusting exposures reduced to scalar lag means. The two analyses broadly agree on the response shapes and temperature optimum, with the largest divergence for specific humidity. We do not interpret that divergence as evidence in favor of either model because the analyses also differ in adjustment representation and treatment of the survey hierarchy: weight-decayed level-specific effects in the DLNAM and variance-component random intercepts in the DLNM fit. The application demonstrates that several fully lagged exposures can be represented together at this scale. At the lag-specific level, the fitted DLNAM temperature surface also concentrates the elevated association around the shortest lags, consistent with the source analysis. Figure 5: DHS/MIS Malaria: Model Comparison. Cumulative odds-ratio functions and corresponding DLNAM and DLNM exposure–lag response surfaces are shown for five environmental exposures. Shaded bands around the cumulative functions denote pointwise 95%95\% intervals. Each column reports the target-specific fit for that exposure; the analyses differ in adjustment representation and hierarchy. Analysis specifications are reported in Table S3. 2.8 Computational Scaling A fixed-budget benchmark measured absolute fitting time and scaling as concurrently lagged exposures were added (Fig. 6; Table S3). On the reported workstation, DLNAM was slower in absolute wall-clock time than the comparators for one exposure. At a fixed 2500-epoch DLNAM budget, runtime increased 4.3-fold from one to four exposures, compared with 6.2-fold for T-DLNM, 12.5-fold for DLNM (QBIC), 23.7-fold for DLNM (QAIC), and 33.5-fold for P-DLNM. The production joint DLNAM used 5000 epochs; under approximately linear scaling with training budget, this would imply an approximately 8.6-fold increase over the one-exposure benchmark. This extrapolation exceeds the observed T-DLNM increase, although its multi-exposure implementation uses target-specific treed fits and represents remaining exposures by spline adjustment terms rather than fitting a single joint treed model. The benchmark therefore compares implemented procedures rather than establishing an unconditional efficiency ranking; absolute timings depend on hardware and implementation. Figure 6: Computational Scaling: Model Comparison. Absolute end-to-end fit time and runtime relative to each method’s one-exposure fit are shown across one to four exposures, with R=3R=3 repeats per method and exposure count. Markers denote the median across repeats; bars denote the observed minimum–maximum range. Benchmark specifications are reported in Table S3. 3 Discussion The central methodological result is the model class itself: DLNAMs learn smooth exposure–lag response surfaces as neural-additive components while preserving the additive structure and familiar scientific summaries of DLNMs. On the more structurally demanding surfaces, comparator error rises chiefly through squared bias, whereas DLNAM’s remaining cumulative error is predominantly variance. Representation bias systematically displaces the reported curve when the fitted function class cannot express the target; sampling variability, by contrast, can recede with additional information. The same contrast is especially visible in sparsely supported exposure tails. On the more structurally demanding DGPs, DLNAM retained lower boundary error and its learned components tended to flatten rather than display the oscillatory edge behavior seen in some spline fits. This suggests a practical advantage where support is weakest, without implying that spline estimators are generically unstable at boundaries. Some remaining variability reflects random initialization and can be reduced by enlarging the ensemble, at additional computational cost. Recovery alone is insufficient: nominal interval formulas do not guarantee nominal coverage, so we treat calibration as an empirical property, as do the comparator method studies [7, 8]. In the principal simulation, DLNAM is closest to nominal cumulative coverage in 11 of 12 DGP-by-region cells; in the joint experiment it is closest in all 12 active component-by-region cells. On the full-surface target, every estimator with higher coverage than DLNAM also reports wider intervals, with only two exceptions to the analogous pattern on the cumulative target. The advantage is therefore one of calibration relative to interval width, not coverage in isolation. The interval construction is deliberately pragmatic: combining a conditional last-layer Laplace approximation with between-member ensemble variation avoids repeated bootstrap refitting and the posterior-sampling burden of Markov chain Monte Carlo. Conditional on the learned representation, final-layer inference reduces to a low-dimensional generalized linear model calculation, giving a neural-additive analogue of pointwise DLNM inference [6]. These are not simultaneous bands or prediction intervals and do not propagate all uncertainty in the learned representation; the ensemble term captures the part expressed as disagreement across independently initialized fits. A neural estimator exposes a large nominal hyperparameter space, but nominal choice is not recurring choice. Operationally, the relevant question is which decisions must be reopened from one analysis to the next. DLNMs require the analyst to prescribe or select a cross-basis representation, whether through substantive judgment or a selection procedure. P-DLNM learns effective smoothness within a penalized representation, but the marginal smoother and nominal basis dimension remain modeling choices. T-DLNM instead learns local partition structure under a specified tree-ensemble, prior, and sampler configuration [7, 8]. A spline cross-basis fixes the coordinates of the function space before estimation, after which fitting estimates their coefficients; DLNAM learns its internal features jointly with the response. Across the analyses here, the same architecture remained effective over two outcome families, three lag horizons, and one to five simultaneous exposures; fit-specific changes were concentrated primarily in optimization—training budget, learning rate, or batching—when the fitting problem changed materially (Section 4.7). This does not establish universal tuning insensitivity; it does show that a large nominal neural design space need not translate into a comparably large recurring specification burden. The ablation shows that learning the representation is nontrivial: recovery depends on how the neural component is structured, not merely on replacing a cross-basis with a larger feed-forward network. Exposure and lag enter jointly through a multivariate ExU input layer, several jointly trained subnetworks compose the component, and smooth hidden activations regulate how local sensitivity propagates through it. Removing any of these elements degrades recovery with a distinct error or coverage signature. The contribution is therefore architectural as well as representational: local sensitivity, compositional flexibility, smoothness, and stable fitting must coexist inside an additive component. Independent complete fits are then ensembled across initializations to temper fit-to-fit instability. The ExU comparison and subnetwork decomposition are reported in Supplementary Texts S6 and S5. The additive form is especially important when several lagged inputs are modeled together. A cross-basis analysis can also accommodate multiple lagged predictors, but every additional surface introduces another representation to specify. In DLNAMs, additional lagged predictors enter as further learned components. The joint simulation tests component recovery under correlated exposures; the malaria analysis shows that the construction remains feasible at application scale. The substantive scientific decisions—which exposures to include, what to adjust for, and which contrasts to report—remain. What need not be redesigned for every added exposure is the learned response representation. The malaria analysis still fits one model per reported target exposure, but within each fit every adjusting exposure retains its full distributed-lag surface rather than being reduced to the scalar lag-mean adjustment used in the comparison. These advantages do not imply that DLNAMs are preferable in every setting. When modeling a single exposure with a comparatively simple smooth surface adequately represented by a low-dimensional cross-basis, a DLNM is simpler, competitive, and cheaper to fit in this benchmark. Abrupt critical-window structure can likewise favor a treed representation over a smooth estimator [8]. Learning the representation entails greater computational demand in the reported setup. Under a fixed training budget, however, DLNAM runtime increased less steeply as concurrent exposures were added than the comparator implementations. The larger training budget used for the production joint fit cautions against treating that benchmark as an unconditional efficiency ranking. Comparison with T-DLNM is further qualified by its target-specific rather than fully joint multi-exposure implementation. Absolute runtime also depends on hardware and software implementation, so the reported timings should not be interpreted as a fixed computational property of the model class. The evidence is necessarily finite in scope: the principal comparison uses R=200R=200 replicates and the architecture and joint-exposure studies R=50R=50. Monte Carlo uncertainty remains relevant for closely spaced quantities, although in the joint study Monte Carlo standard errors are small relative to the RMSE separations among active components, so the principal recovery ranking remains well resolved. The joint experiment evaluates a single correlation structure and therefore does not establish how recovery changes across a broader range of between-exposure dependence. The simulations span several distinct surface geometries but do not define a universal threshold at which a learned neural representation becomes preferable to a simpler smoother. The applications target recovery of response structure rather than prediction; out-of-sample evaluation is secondary to that aim and can be nontrivial with dependent lag histories. The model remains additive across exposures. Selected pairwise terms extend readily to additive models, including NAMs [19, 12], but admit non-equivalent distributed-lag formulations; [20] provide one structured alternative, whose integration with DLNAMs remains future work. Two further compositional extensions follow directly. Outcome-specific mixing weights over shared subnetworks define a Multitask DLNAM in the manner of multitask NAMs [12]. For one outcome it reduces to the Mixture DLNAM, of which the standard DLNAM is a special case. Shared subnetworks induce cross-outcome regularization, potentially improving recovery for related outcomes, but confine outcome-specific flexibility to the mixing weights. A Hierarchical DLNAM adds partially pooled group-specific deviations to global components, reducing to the Mixture DLNAM when these vanish and paralleling hierarchical DLNMs [21]. Hierarchical DLNMs provide a single-stage alternative to the two-stage pooling common in multicountry analyses [2]; the malaria analysis motivates the corresponding DLNAM because country and survey cluster currently enter only through level-specific intercepts. Neither extension was implemented or evaluated here; formal definitions are given in Supplementary Text S7. Beyond these extensions, the reported specification was not exhaustively optimized; the ablations show that changes to individual modules can materially alter recovery, while its modular construction allows modules to evolve without redesigning the model. Component-specific shrinkage [17] and fuller propagation of uncertainty in the learned representation [16] are particularly immediate directions suggested by the present results. The reported architecture is therefore an effective first realization rather than a performance ceiling. At the application level, fitted DLNAMs could provide the temperature–mortality response component in operational frameworks that translate weather forecasts into real-time health forecasts [22], using trajectories from AI weather models such as Pangu-Weather and FourCastNet [23, 24], reflecting the growing role of deep learning in weather forecasting [25]. The present results therefore identify a specific role for DLNAMs: settings in which the response structure is sufficiently demanding that a prescribed representation may become restrictive, several lagged predictors must be modeled together, or repeated cross-basis specification is itself a substantial part of the analysis. Simpler established estimators remain attractive when a modest smooth representation is adequate, particularly where computation is limiting. More broadly, DLNAMs provide a route for translating advances in deep-learning methodology into distributed-lag models while preserving additive interpretability and familiar distributed-lag summaries. 4 Methods 4.1 Distributed Lag Neural Additive Model Let Ytt=1n\Y_t\_t=1^n, xtt=1n\x_t\_t=1^n, and utt=1n\u_t\_t=1^n denote outcomes, exposures, and additional covariates, respectively, with Yt∈Y_t , xt∈ℝMx_t ^M and ut∈ℝKu_t ^K. Writing ℒ=(ℒ1,…,ℒM)L=(L_1,…,L_M) for the collection of finite integer-lag sets, define xtℒ=((xt−ℓ,m)ℓ∈ℒm)m=1M,ℋt=(xtℒ,ut),x_t^L= ((x_t- ,m)_ _m )_m=1^M, _t=(x_t^L,u_t), and let μt=(Yt∣ℋt) _t=E(Y_t _t) denote the conditional mean outcome given the exposure history and additional covariates. For an invertible link function g, the DLNAM predictor is g(μt)=α+∑m=1M∑ℓ∈ℒmfmθmx(xt−ℓ,m,ℓ)+∑k=1Khkθku(utk).g( _t)=α+ _m=1^M _ _mf_m _m^x(x_t- ,m, )+ _k=1^Kh_k _k^u(u_tk). (1) Recurring notation is summarized in Supplementary Table S1. The exposure–lag response surfaces fmθmxf_m _m^x and covariate shape functions hkθkuh_k _k^u are, for continuous inputs, represented by feed-forward neural networks, with θ=(α,θx,θu)θ=(α,θ^x,θ^u) denoting the complete parameter vector. In DLNMs, the exposure–lag function is instead defined through a prescribed cross-basis combining marginal basis expansions over exposure and lag, with the resulting contribution evaluated and summed over the lagged exposure history [4, 5, 6]. Written in terms of the response surface implied by that cross-basis, DLNMs share the predictor form of Equation (1), with the defining distinction lying in the representation of fmθmxf_m _m^x: DLNAMs learn this function directly as a neural component, while terms corresponding to hkθkuh_k _k^u may enter parametrically or through separate smooths. For identifiability and interpretation, x⋆∈ℝMx ^M and u⋆∈ℝKu ^K define the reference point, and fitted components are expressed as contrasts satisfying fmθmx(xm⋆,ℓ)=0,ℓ∈ℒm,hkθku(uk⋆)=0,f_m _m^x(x_m , )=0, _m, h_k _k^u(u_k )=0, with the resulting constants absorbed into α, so that α is the linear predictor at the reference values and each component is interpreted as a contrast on the link scale. Motivated by the compositional subnetwork architecture of multitask NAMs [12], our primary specification, the Mixture DLNAM, generalizes the standard DLNAM by replacing the direct parameterization of each component with a learned linear combination of jointly trained subnetworks, fmθmx(⋅,⋅)=∑s=1Smxωmsxf~msϑmsx(⋅,⋅),hkθku(⋅)=∑s=1Skuωksuh~ksϑksu(⋅).f_m _m^x(·,·)= _s=1^S_m^x _ms^x f_ms _ms^x(·,·), h_k _k^u(·)= _s=1^S_k^u _ks^u h_ks _ks^u(·). (2) Under this parameterization, ϑ=(ϑx,ϑu) =( ^x, ^u) collects the subnetwork parameters and ω=(ωx,ωu)ω=(ω^x,ω^u) their mixing weights, so that θ=(α,ϑ,ω)θ=(α, ,ω). The mixing weights are unconstrained rather than restricted to form a convex combination, allowing subnetworks to reinforce or partially offset one another. The standard DLNAM is nested within Equation (2) as the special case in which each component contains a single subnetwork and the mixing weights form a convex combination. Allowing the sole mixing weight to remain unconstrained leaves the resulting function class unchanged, since that scalar can be absorbed into the output layer, but changes the optimization parameterization: its independent initialization rescales the component and can temper its early contribution and gradients. The distinction between the standard DLNAM and a one-subnetwork Mixture DLNAM is therefore computational rather than functional. To reduce variability attributable to random initialization, we ensemble B independently initialized DLNAMs. Mixture composition and ensembling operate at different levels: subnetworks are optimized jointly within each component of a fitted model, whereas the ensemble averages complete fitted models. Supplementary Text S5 illustrates how weighted subnetworks combine through reinforcement and partial cancellation within one ensemble member. The ensemble is combined on the link scale, ηt=1B∑b=1Bηt(θ(b)),ηt(θ(b))=g(μt(θ(b))). _t= 1B _b=1^B _t\! (θ^(b) ), _t\! (θ^(b) )=g\! ( _t\! (θ^(b) ) ). (3) By additivity, the same averaging applies componentwise to centered exposure–lag functions, covariate effects, and their derived cumulative contrasts. Reported summaries are obtained only after this averaging by applying the appropriate transformation from the link scale; Section 4.4 extends the ensemble construction to pointwise uncertainty. 4.2 Neural Component Architecture For continuous components that must resolve local nonlinear structure, we use exp-centered hidden units (ExU) [12] in the input layer; simpler covariate components may instead use linear inputs. Let z∈ℝpz ^p denote the input to a network component. A standard input layer maps z affinely and then applies an elementwise activation ϕφ, h:ℝp→ℝdh:R^p ^d, h(z)=ϕ(Wz+b),W∈ℝd×p,b∈ℝd,h(z)=φ(Wz+b), W ^d× p, b ^d, where W is the weight matrix and b the bias vector. ExU preserves this map-then-activation structure but reparameterizes the input map through log-scale weights and unit-specific locations. The scalar ExU used in NAMs acts on a univariate input zjz_j. For hj:ℝ→ℝdjh_j:R ^d_j and wj,bj∈ℝdjw_j,b_j ^d_j, hj(zj)=ϕ(ewj⊙(zjdj−bj)),h_j(z_j)=φ (e^w_j (z_j1_d_j-b_j) ), (4) where ⊙ denotes the Hadamard product, wjw_j are log-scale weights, and bjb_j are location biases. Exponential scaling gives each unit locally adjustable sensitivity around bjb_j, allowing the input map to resolve thresholds and rapid transitions. For multivariate inputs, we extend the scalar ExU coordinate-wise by concatenation, with d=∑j=1pdjd= _j=1^pd_j: h(z)=(h1(z1)⊤,…,hp(zp)⊤)⊤∈ℝd.h(z)= (h_1(z_1) ,…,h_p(z_p) ) ^d. (5) Alternative multivariate generalizations are compared in Supplementary Text S6. For an exposure–lag surface, h(x,ℓ)=(hx(x)⊤,hℓ(ℓ)⊤)⊤h(x, )=(h_x(x) ,h_ ( ) ) , with the layer width split as evenly as possible between the two coordinates, dx=⌈d/2⌉d_x= d/2 and dℓ=⌊d/2⌋d_ = d/2 . The architecture pairs ExU with smooth activations rather than the capped ReLU-n used in the original NAM construction [12], where capping localizes individual ExU units to a restricted input range. We consider two closely related smooth activation functions, Mish [13] and SiLU [14], both of which provide smooth ReLU-like nonlinearities. Mish is particularly attractive for its smooth profile and because its original study reported performance matching or improving upon SiLU (Swish-1) across several deep-network benchmarks. Their activation functions are defined as ϕMish(z)=ztanh(softplus(z)),ϕSiLU(z)=zσ(z), _Mish(z)=z\, \! (softplus(z) ), _SiLU(z)=z\,σ(z), where softplus(z)=log(1+ez)softplus(z)= (1+e^z) is the Softplus function and σ(z)=(1+e−z)−1σ(z)=(1+e^-z)^-1 is the logistic sigmoid. Their continuous derivatives preserve a smooth reported response. The reported fits use Mish; SiLU is a close alternative from the same family and was not separately evaluated. The division of labor is deliberate: ExU controls the location and scale of sensitivity, while the hidden activation controls how smoothly that sensitivity propagates through the component. This architectural distinction is evaluated in the ablation study (Fig. 2). 4.3 Optimization Parameters are estimated by minimizing the average negative log-likelihood, ℒ(θ)=−1n∑t=1nlogpYt∣μt(θ),L(θ)=- 1n _t=1^n p\Y_t _t(θ)\, (6) where p denotes the outcome likelihood or working likelihood as appropriate, ηt(θ) _t(θ) denotes the predictor in Equation (1), and μt(θ)=g−1(ηt(θ)) _t(θ)=g^-1( _t(θ)). Optimization uses AdamW [26] under cosine annealing [27], with gradient-norm clipping and full-batch updates unless otherwise stated. In the quasi-Poisson application, optimization proceeds under a Poisson working likelihood [28], with overdispersion estimated post hoc from the Pearson dispersion and carried into the quasi-likelihood uncertainty calculation of Section 4.4. The reported fits use weight decay, with coefficient λw _w, as their principal explicit regularizer, together with a surface-roughness penalty in the Chicago analysis. Standard dropout, subnetwork dropout, and an output penalty were also implemented and examined in preliminary simulation work, but did not yield consistent improvements in recovery and were therefore inactive in the reported analyses. The Chicago roughness penalty is applied because the unpenalized fit is sensitive to the training budget beyond roughly 3000 epochs, whereas under the penalty the fitted surface is stable across budgets from 2500 to 10000. Beyond explicit regularization, the architecture and fitting procedure shape the learned representation in several distinct ways: smooth activations favor smooth component functions, mixture composition represents each component through jointly trained subnetworks rather than a single network, and ensemble averaging reduces variation attributable to random initialization. Optimizer, scheduler, batching, gradient-clipping, and regularization settings are reported in Table S2. Input scaling is matched to the parameterization of the first layer and is treated as part of the fitted model. Inputs to ExU layers are min–max scaled to [0,1][0,1], aligning the observed domain with the initialization range of the location biases, whereas linear inputs are standardized by their analyzed-sample mean and standard deviation. Exposure scaling pools the raw values supplied across a term’s observed lag columns, giving each exposure a common scale across its lag history; lag coordinates are scaled using the fixed lag range used for fitting. Data-dependent scaling constants are estimated once within each fitted data set, and all transformations are then held fixed for training and effect extraction, so evaluation-grid values are transformed using the fitted scaling constants rather than rescaled to the evaluation grid. Initialization likewise reflects the role of each parameter in the architecture. ExU log-weights are drawn from (μw,σw2)N( _w, _w^2) and location biases from (0,1)U(0,1), so that μw _w controls the initial sensitivity scale through ewe^w while σw _w disperses that sensitivity across units. Ordinary linear layers use Kaiming uniform initialization [29], and mixing weights are initialized independently from (0,σω2)N(0, _ω^2) across components, subnetworks, and ensemble members. 4.4 Last-Layer Laplace Confidence Intervals Uncertainty is quantified with a target-specific last-layer Laplace approximation applied separately to each ensemble member, related to the Laplace treatment of NAMs developed by [17]. The approximation conditions on the learned representation and, in the mixture formulation, the mixing weights; the resulting member-level conditional uncertainty is combined with between-member variation. Supplementary Text S2 gives the full derivation and regularity conditions, and frequentist calibration is evaluated directly in simulation. For a given inferential target, partition the parameter vector as θ=(ψ,ξ)θ=(ψ,ξ), where ξ∈ℝqξ ^q contains the global intercept and the final-layer parameters of the corresponding component, while ψ contains all remaining parameters. Conditional on the fitted representation ψ ψ, the predictor is affine in ξ, reducing the calculation to a generalized linear model in a low-dimensional last-layer parameterization. The mixing weights remain conditioned because allowing them to vary jointly with the corresponding output layers would introduce a bilinear scale non-identifiability. Let η(ξ)=(η1(ξ),…,ηn(ξ))⊤η(ξ)=( _1(ξ),…, _n(ξ)) denote the vector of linear predictors across the n observations, with each ηt _t given by Equation (1), and define the last-layer Jacobian Φ=∂η(ξ)∂ξ⊤∈ℝn×q, = ∂η(ξ)∂ξ ^n× q, (7) which plays the same design-matrix role as the cross-basis in a DLNM. Let W=diag(w1,…,wn)W=diag(w_1,…,w_n), where wt=(dμt/dηt)2/V(μt)w_t=(d _t/d _t)^2/V( _t) are the generalized-linear-model working weights evaluated at the fitted means and V is the variance function. Let φ denote the dispersion or working-dispersion factor appropriate to the fitted outcome model, taking its family-specified value when fixed and being estimated where required. The dispersion-adjusted last-layer information is H=Φ⊤WΦ/φ^.H= W / . (8) Because η(ξ)η(ξ) is affine in ξ, Φ is exact and constant and second derivatives of the predictor with respect to ξ vanish, so no predictor-curvature term enters the Hessian. For canonical likelihood fits, H coincides with the observed likelihood information; for non-canonical links, it is the corresponding Fisher information, while for working-likelihood fits it is interpreted as dispersion-scaled working information. For the post-hoc Laplace approximation, a zero-mean Gaussian prior with precision λ is placed on the final-layer parameters but not on the intercept, a split encoded by the diagonal selector PξP_ξ. The precision λ λ is selected by a MacKay evidence fixed-point approximation with the unpenalized intercept profiled out [30, 31]. Because λ λ is selected after training, the fitted ξ ξ is not re-optimized under the Gaussian prior; it is retained as the expansion center and treated as an approximate posterior mode. The conditional last-layer covariance is Σ^ξ=(H+λ^Pξ+ϵI)−1, _ξ= (H+ λP_ξ+ε I )^-1, (9) where ϵ>0ε>0 is a numerical ridge. The inferential targets are the centered components of Equation (1): the exposure–lag surface fmψ^,ξf_m ψ,ξ, its cumulative Fmψ^,ξ(⋅)=∑ℓ∈ℒmfmψ^,ξ(⋅,ℓ),F_m ψ,ξ(·)= _ _mf_m ψ,ξ(·, ), (10) and the covariate shape function hkψ^,ξh_k ψ,ξ. For τψ^,ξ∈fmψ^,ξ,Fmψ^,ξ,hkψ^,ξτ ψ,ξ∈ \f_m ψ,ξ,F_m ψ,ξ,h_k ψ,ξ \ evaluated at point ζ, let Jζ=∂τψ^,ξ(ζ)∂ξ⊤J_ζ= ∂τ ψ,ξ(ζ)∂ξ denote its gradient. Conditional on ψ ψ, these targets are affine in ξ, so JζJ_ζ is constant and covariance propagation is exact, with no delta-method linearization error. The pointwise conditional standard deviation on the link scale is sdτψ^,ξ(ζ)=JζΣ^ξJζ⊤1/2.sd\τ ψ,ξ(ζ)\= \J_ζ _ξJ_ζ \^1/2. (11) Pointwise 1−α1-α intervals are first formed on the link-contrast scale and, where the reported estimand is a monotone transformation of that contrast, their endpoints are transformed accordingly. When an absolute linear predictor is reported on the mean scale, the transformation is g−1g^-1. For the ensemble, write ψ^=(ψ^1,…,ψ^B),ξ^=(ξ^1,…,ξ^B) ψ=( ψ_1,…, ψ_B), ξ=( ξ_1,…, ξ_B) for the collections of fitted member-specific conditioned parameters and last-layer expansion centers, respectively. Applying the ensemble construction of Equation (3) to a reported target, the interval is centered at the ensemble mean and its variance combines average member-level conditional uncertainty with between-member variation, τ^ψ^,ξ^(ζ)=1B∑b=1Bτ^bψ^b,ξ^b(ζ),^(ζ)=1B∑b=1BJbζΣ^ξbJbζ⊤⏟within-member+1B∑b=1Bτ^bψ^b,ξ^b(ζ)−τ^ψ^,ξ^(ζ)2⏟between-member. τ ψ, ξ(ζ)= 1B _b=1^B τ_b ψ_b, ξ_b(ζ), V(ζ)= 1B _b=1^BJ_bζ _ _bJ_bζ _within-member+ 1B _b=1^B \ τ_b ψ_b, ξ_b(ζ)- τ ψ, ξ(ζ) \^2_between-member. (12) Combining the conditional last-layer Laplace uncertainty with between-member variation, the reported pointwise 1−α1-α confidence interval is τ^ψ^,ξ^(ζ)±z1−α/2^(ζ)1/2, τ ψ, ξ(ζ)± z_1-α/2 V(ζ)^1/2, (13) with the same transformation to the reported scale as above. The first term in Equation (12) averages member-specific conditional last-layer variances without further division by B, while the second records between-member variation. Consequently, τ^ψ^,ξ^(ζ) τ ψ, ξ(ζ) and ^(ζ) V(ζ) are the mean and variance of the equally weighted mixture of member-level Gaussian Laplace approximations. The reported interval is a Gaussian moment-matched summary of this generally non-Gaussian mixture; its variance is therefore a constructed uncertainty measure, not the sampling variance of the ensemble mean or full posterior uncertainty in the learned representation. 4.5 Simulation Design For each DGP, we simulated a single-exposure distributed lag Poisson model with lag horizon L=14L=14, n=5000n=5000 daily observations, and reference exposure x⋆=20x =20: Yt∣ℋt∼Poisson(μt),logμt=α+∑ℓ=0Lf(xt−ℓ,ℓ),Y_t _t ( _t), _t=α+ _ =0^Lf(x_t- , ), with α=log50α= 50. Exposure followed a seasonal autoregressive process, xt=20+12sin(2πt/365)+at,at=0.8at−1+εt,εt∼(0,22),x_t=20+12 (2π t/365)+a_t, a_t=0.8\,a_t-1+ _t, _t (0,2^2), with a1a_1 drawn from its stationary distribution. Each replicate included an additional L-observation exposure prehistory, leaving n=5000n=5000 analyzed observations with complete lag histories. Effects were evaluated on the equally spaced grid =0,0.2,…,40,G=\0,0.2,…,40\, which includes the reference x⋆x . The four generating surfaces f1,…,f4f_1,…,f_4 were constructed from sigmoidal exposure responses, Gaussian bumps, and normalized lag kernels; Supplementary Text S3 gives their exact definitions. DGP 1 is separable, with a smooth exposure response modulated by a single exponentially decaying lag kernel across the exposure domain. DGPs 2–4 are nonseparable and differ in how lag structure varies across the exposure domain: DGP 2 combines early and late effects with opposing sign at intermediate lags, DGP 3 combines a narrow high-exposure response with a broader low-exposure response on a slower kernel, and DGP 4 places sharp responses at both extremes of the exposure range with distinct dominant lag timing. All four data-generating surfaces are smooth and continuous, so the spline-based comparators are evaluated on targets compatible with their intended function class rather than on discontinuous or otherwise deliberately unfavorable alternatives. 4.6 Evaluation Metrics The likelihood depends on each exposure history through the aggregate contribution ∑ℓ=0Lf(xt−ℓ,ℓ) _ =0^Lf(x_t- , ). Under the log link, we take the cumulative effect F(x)=∑ℓ=0Lf(x,ℓ)F(x)= _ =0^Lf(x, ) as the primary estimand, corresponding to the cumulative log-relative-risk summary conventionally reported in distributed lag nonlinear model analyses [6, 7, 8]. This contrast holds exposure fixed across lags and sums the resulting surface values. We additionally evaluate f as a complementary diagnostic of lag-specific allocation. Because individual surface locations are not observed by the likelihood in isolation, surface scoring evaluates a finer decomposition of the fitted effect rather than a more directly fitted or intrinsically superior target. Let τ∈F,fτ∈\F,f\ denote the evaluation target, with target-specific grids F=,f=×0,…,L.G_F=G, _f=G×\0,…,L\. Across R Monte Carlo replicates, define drτ(ζ)=τ^r(ζ)−τ(ζ),d¯τ(ζ)=1R∑r=1Rdrτ(ζ),ζ∈τ.d_r^τ(ζ)= τ_r(ζ)-τ(ζ), d^τ(ζ)= 1R _r=1^Rd_r^τ(ζ), ζ _τ. For each target, summaries are reported over the full grid, the boundary region below the 55th or above the 9595th empirical percentile of a fixed reference exposure series, and its interior complement. The reference series is generated once from the simulation exposure process at the nominal sample size and held fixed across replicates and estimators, ensuring identical evaluation regions between fits. Surface regions are defined by the exposure coordinate and include every lag. Writing τ⊆τA_τ _τ for any such region, RMSEτ(τ)=1R∑r=1R1|τ|∑ζ∈τdrτ(ζ)2,Biasτ2(τ)=1|τ|∑ζ∈τd¯τ(ζ)2,RMSE_τ(A_τ)= 1R _r=1^R 1|A_τ| _ζ _τ\d_r^τ(ζ)\^2, _τ^2(A_τ)= 1|A_τ| _ζ _τ\ d^τ(ζ)\^2, and Varτ(τ)=1R∑r=1R1|τ|∑ζ∈τdrτ(ζ)−d¯τ(ζ)2.Var_τ(A_τ)= 1R _r=1^R 1|A_τ| _ζ _τ\d_r^τ(ζ)- d^τ(ζ)\^2. Under these definitions, RMSEτ2(τ)=Biasτ2(τ)+Varτ(τ),RMSE_τ^2(A_τ)=Bias_τ^2(A_τ)+Var_τ(A_τ), so squared bias and variance partition the mean squared error exactly; the reported bias and variance shares are their corresponding proportions of RMSEτ2(τ)RMSE_τ^2(A_τ). If Irτ(ζ)I_r^τ(ζ) denotes the pointwise 95%95\% interval for target τ in replicate r, empirical coverage is Covτ(τ)=1R∑r=1R1|τ|∑ζ∈τ(ζ)∈Irτ(ζ).Cov_τ(A_τ)= 1R _r=1^R 1|A_τ| _ζ _τ1 \τ(ζ)∈ I_r^τ(ζ) \. For DLNAM, Irτ(ζ)I_r^τ(ζ) is the combined last-layer Laplace and ensemble interval of Equation (13). This evaluates pointwise coverage of the fixed generating function, not simultaneous coverage over τG_τ. Monte Carlo uncertainty is calculated from replicate-level summaries, thereby accounting for dependence among evaluation-grid points within each replicate rather than treating them as independent Monte Carlo replicates. For any reported summary S S with Monte Carlo standard error MCSE(S^)MCSE( S), plotted 95%95\% Monte Carlo confidence intervals use the normal approximation S^±1.96MCSE(S^). S± 1.96\,MCSE( S). Coverage MCSEs use replicate-level variability directly; RMSE uncertainty is obtained by applying the delta method to replicate-wise mean squared error, while regional squared-bias and variance uncertainties use delete-one-replicate jackknife estimates. The complete replicate-level, delta-method, and jackknife calculations are given in Supplementary Text S4. 4.7 Model Specifications and Fitting All estimators were fitted to the same replicate data sets, evaluated on the same grids, and centered at the same reference value. DLNAM and the spline-based comparators used the log link. DLNAM minimized the negative average Poisson log-likelihood in Equation (6); DLNM (QAIC), DLNM (QBIC), and P-DLNM were fitted as quasi-Poisson regressions, with the Poisson log-likelihood at the fitted means entering the QAIC and QBIC calculations. The likelihood difference for T-DLNM is described below. Comparator configurations and, where applicable, tuning ranges were anchored in the originating methodological papers and their simulation studies rather than optimized for the present DGPs; basis- or outcome-specific adaptations and departures required by the application or multi-exposure setting are stated explicitly below. The comparisons should therefore be read as benchmarks against established, literature-grounded implementations in the settings they were designed to address, rather than against deliberately weak baselines. DLNAM. The reported DLNAM used 3 subnetworks per component and 3 independently initialized complete fits, with uncertainty constructed as in Section 4.4. Each surface subnetwork comprised a 128-unit ExU input layer, one 128-unit dense hidden layer with a Mish activation, and a linear output. ExU log-weight initialization used μx=1.5 _x=1.5 for exposure and μℓ=2.5 _ =2.5 for lag, with σw=0.5 _w=0.5 for both. The larger μℓ _ gives greater initial sensitivity along the lag coordinate, reflecting the expectation of more rapid variation at short lags that also motivates logarithmic lag-knot placement in conventional DLNMs [6]; the ExU scales themselves remain learned parameters. This surface-component architecture was fixed before the comparison and carried across the reported analyses, while optimization budgets were allowed to vary where required by problem size. Developed in preliminary simulation work, the configuration is treated as an empirically supported specification rather than as a globally optimal design. The principal simulation and Chicago analysis used 5000 configured full-batch epochs, with no early stopping. Study-specific optimization departures for the joint-exposure simulation and malaria analysis are described in their respective sections. Exact parameter counts and remaining fitting settings are given in Table S2. DLNM. The two criterion-selected variants, DLNM (QAIC) and DLNM (QBIC), follow the conventional DLNM framework of [4] and the dlnm implementation [6]: quasi-Poisson log-linear fitting with an exposure–lag cross-basis, natural cubic spline marginal bases, equally spaced exposure knots, logarithmically spaced lag knots, and information-criterion selection of basis complexity. Criterion-based selection of cross-basis complexity has also been evaluated in simulation studies [5]. The candidate range was anchored in prior methodological simulation studies: [7] selected the complexity of an unpenalized DLNM over combinations spanning 1–10 degrees of freedom in each marginal dimension, a range subsequently used for the corresponding spline comparator in the T-DLNM simulation study of [8]. We use the corresponding range 2,…,10\2,…,10\ in each margin, excluding one degree of freedom because under the present natural-spline construction it contains no interior knot and reduces the margin to a linear term. Thus, in simulation, both criteria are evaluated over 2,…,102\2,…,10\^2, giving 81 candidate cross-bases per replicate. Chicago uses the narrower range 2,…,52\2,…,5\^2. Sensitivity analysis of this model structure reports that increasing the number of exposure knots produces a visibly less smooth exposure–response curve attributable to overfitting, while the lag dimension is comparatively insensitive to basis choice [4]; the applied range is restricted accordingly. On this series both criteria select the largest available exposure dimension, so the upper limit acts as a specification choice rather than as a search boundary, and widening it reproduces that instability, with QAIC placing minimum mortality at −19.7∘C-19.7\, C. Pointwise Wald intervals are formed from the fitted cross-basis contrasts and coefficient covariance. QAIC and QBIC differ in their complexity penalty, with QBIC penalizing larger cross-bases more strongly at the present sample sizes; exact criterion formulas are reported in Supplementary Text S9. For these quasi-Poisson fits, the estimated dispersion enters both the coefficient covariance and the criterion calculations. Penalized DLNM. P-DLNM follows the penalized framework of [7], implemented using the dlnm cross-basis construction and mgcv penalized-GAM machinery [6, 32]. Marginal cubic P-spline bases are coupled through the framework’s cross-basis penalty construction, with second-order difference penalties and smoothing parameters estimated by restricted maximum likelihood. In simulation, we use the rank-10 marginal specification used in the P-DLNM simulation study, giving 10×1010× 10; the same rank-10, second-order-penalty, REML specification was subsequently used as the penalized spline comparator in the T-DLNM simulation study [8]. In simulation, realized effective degrees of freedom reached at most 57% of the available rank. On the Chicago series the penalty does not play the same role: realized effective degrees of freedom track the imposed rank closely, reaching 78% of the basis at the specification used, and both fitted complexity and the roughness of the cumulative curve increase with rank rather than stabilizing. Marginal rank is therefore a specification choice on this series rather than an upper bound resolved by the penalty. Chicago accordingly uses a fixed 5×55× 5 specification, matching the dimension available to the criterion-selected variants, so that the spline comparators differ in how complexity is controlled rather than in the dimension available to them. Pointwise Wald intervals use the fitted quasi-Poisson coefficient covariance, including its estimated dispersion. Treed DLNM. T-DLNM uses the nonlinear Treed Distributed Lag Model of [8], implemented in the authors’ dlmtree package [33]. We retain the tree priors and sampler configuration used in the authors’ simulation study rather than retuning T-DLNM to the present DGPs; exact tree and Markov chain settings are reported in Table S3. Because the dlmtree version used here does not accept a count likelihood in its nonlinear implementation, T-DLNM was fitted with a Gaussian likelihood for log(1+Yt) (1+Y_t). It therefore serves as a benchmark for the flexible treed method class rather than as a likelihood-matched quasi-Poisson comparator. Its exponentiated contrasts in [log(1+Y)]E[ (1+Y)] are used as a comparison-scale proxy for relative risk; they are not algebraically identical to contrasts in log(Y) \E(Y)\. Posterior credible intervals were extracted on the common evaluation grid and reference. Implementation-specific fixed-effect and retry details are given in Supplementary Text S9. 4.8 Joint-Exposure Simulation The multi-exposure recovery study reuses the four data-generating surfaces of Section 4.5 as concurrent additive effects. Five exposure series share the marginal seasonal AR(1) structure of the single-exposure simulation; the four active series have seasonal phase shifts 0,π/6,π/3,π/20,π/6,π/3,π/2, while the null series has phase 00. All five innovation processes are equicorrelated at 0.50.5, and the null exposure has an identically zero true surface. The Poisson log mean is the intercept plus the four active distributed-lag contributions. The final study comprises R=50R=50 independent data sets of n=5000n=5000 observations with the same lag horizon, reference value, and evaluation grid as the principal simulation. The joint DLNAM contains one exposure–lag component per exposure and otherwise retains the specification of the principal simulation. For the larger five-component model, the initial learning rate was reduced from 8⋅10−48· 10^-4 to 7⋅10−47· 10^-4, cosine annealing still reached 1⋅10−41· 10^-4, and training was extended to 50005000 full-batch epochs. Preliminary stability checks at the baseline optimization setting showed occasional high-loss convergence and material under-training after 25002500 epochs. The comparator implementations reflect the feasible multi-exposure form of each method. DLNM (QAIC) and DLNM (QBIC) use one sequential coordinate-wise sweep, initialized at 4×44× 4 marginal dimensions for all five components and updated in the order DGP 1, DGP 2, DGP 3, DGP 4, then the null exposure. For each exposure, its marginal dimensions are selected over 2,…,102\2,…,10\^2 while the current dimensions of the other exposures are held fixed, followed by one joint quasi-Poisson fit. This practical strategy is potentially order-dependent and is not an exhaustive combinatorial search: with 81 candidate cross-bases per exposure, the coordinate-wise sweep requires 5×815× 81 fits, whereas an exhaustive joint search over the five components would require 815≈3.5×10981^5≈ 3.5× 10^9. P-DLNM fits all penalized cross-bases jointly with REML. T-DLNM is fitted separately for each target exposure, with the other concurrent exposures represented by fixed 4×44× 4 natural-spline cross-bases; its target component retains the Gaussian likelihood on log(1+Y) (1+Y) used in the single-exposure comparison. Recovery is scored for each active component on the cumulative target using the full-grid, interior, and boundary regions and pointwise coverage definition of Section 4.6. The null component is scored against zero as a leakage diagnostic under the specified correlation structure. For each active exposure and evaluation region, degradation is defined as the ratio of joint-fit RMSE to the corresponding single-exposure RMSE, with values above one denoting worse recovery in the joint setting. When uncertainty is attached to this ratio, Monte Carlo uncertainty from both RMSE estimates is propagated using the delta-method approximation defined in Supplementary Text S4. 4.9 Chicago Temperature and Mortality Analysis We analyzed the Chicago mortality series distributed with the R dlnm package using the adjustment specification of [4], with PM10PM_10 in place of carbon monoxide, which is not included in the distributed Chicago series. The analysis is not intended to reproduce a published model: all estimators use the same adjustment variables, while their functional representations differ. Daily all-cause mortality is modeled as a function of temperature over lags 00–3030, centered at the observed median of 10.3∘C10.3\, C. Shared adjustment variables are date, day of week, lag 00–11 mean dew point temperature, ozone, and PM10PM_10. Dropping days with any missing model variable leaves 4638 complete days. DLNM (QAIC), DLNM (QBIC), and P-DLNM use quasi-Poisson regression. DLNAM uses a Poisson working likelihood with post-hoc Pearson dispersion φ^=1.03 =1.03 for the reported Laplace intervals, while T-DLNM retains its Gaussian model for log(1+Yt) (1+Y_t). DLNAM retains the exposure–lag surface architecture used in simulation, while the adjustment components use architectures tailored to their roles. Date is represented by a deeper ExU trend network, with an additional 128-unit hidden layer and a larger ExU log-weight mean, μw=4.5 _w=4.5, to accommodate variation across multiple temporal scales; the remaining continuous covariates use separate 32-unit Mish networks with linear inputs after standardization. This allocation follows the relative flexibility of the adjustment specification itself, which assigns 98 degrees of freedom to the long-term trend and at most three to the remaining continuous covariates. The spline and treed comparators follow Section 4.7, including its Chicago-specific dimensions. Full architecture, adjustment, and comparator specifications are reported in Tables S2 and S3. The analysis is observational and in-sample. It assesses whether the estimators yield a familiar empirical temperature–mortality pattern on a canonical DLNM data set and therefore serves as an empirical plausibility analysis rather than an evaluation of causal effects or forecasting performance. 4.10 Multi-Exposure Malaria Analysis The second application uses 348,565 individual childhood malaria outcomes from a large multi-country survey [18]. Five environmental exposures are modeled: mean temperature, precipitation, soil moisture, actual evapotranspiration, and specific humidity. DLNAM uses a Bernoulli likelihood with a logit link and 6 monthly lags spanning months 1–6 before testing, with effects centered at each exposure’s median. One DLNAM is fitted per reported exposure using the target-specific adjustment set of the source study. Exposures treated there as confounders enter alongside the target, whereas mediators are omitted. The adjustment set therefore varies by target, ranging from one additional exposure for temperature and precipitation to four for specific humidity. Table S3 gives the exact target-to-adjustment mapping; each reported curve comes from the fit in which that exposure is the target. In the DLNAM fits, every adjusting exposure retains a full exposure–lag surface, so a target-specific fit contains up to 5 concurrently lagged surfaces. Month and survey year enter through 32-unit Mish covariate networks with standardized linear inputs, while survey cluster and country enter as weight-decayed level-specific effects, permitting joint minibatch optimization; variance-component random intercepts would require a hierarchical extension not implemented here. Fits were trained for 50 configured epochs using minibatches containing 1% of the data rather than the full-batch updates used in the smaller analyses. For the reported malaria intervals, the Laplace approximation includes the target exposure’s final-layer parameters together with the global intercept; all remaining fitted parameters, including the adjusting-exposure networks, month/year covariate networks, and cluster/country effects, are conditioned upon. The DLNM uses the same target-specific environmental adjustment sets and is likewise fitted one target at a time. The target exposure uses a natural-spline cross-basis with 3 degrees of freedom in each exposure and lag margin over months 1–6, while adjusting exposures enter as scalar lag means. Month is represented by a 4-degree-of-freedom spline, survey year enters linearly, and cluster and country enter as random intercepts estimated by glmmTMB [34]. Thus, the two analyses use the same environmental control sets but differ in the functional representation of those controls and in their treatment of the survey hierarchy: cluster and country effects are shrunk by a fixed weight-decay penalty in the DLNAM and by an estimated variance component in the DLNM, so both are pooled but only the latter estimates the degree of pooling. The comparison therefore contrasts two complete analysis specifications rather than isolating the effect of estimator choice. The application assesses feasibility and interpretation in a large, multi-exposure observational setting rather than comparative recovery. Full architecture, optimization, adjustment, and sampling specifications are reported in Table S2. Declarations Acknowledgments. We are grateful to Manuel Martellini O Nocentini, whose suggestions helped shape the experimental design, comparator selection, and evaluation criteria of the simulation study. Funding. This research was supported by Formas, the Swedish Research Council for Sustainable Development (grant nos. 2023-01774 and 2022-01845); Forte, the Swedish Research Council for Health, Working Life and Welfare (grant no. 2024-00833); and the Swedish Research Council (grant nos. 2025-02804 and 2022-06599). Computational resources. The computational work presented in this study was supported by high-performance computing resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS), funded by the Swedish Research Council. These resources were used for model development, training, simulation experiments, and computational testing. Author contributions. Conceptualization: C.H., S.P., and L.O.; Data curation: C.H. and S.P.; Formal analysis: C.H.; Investigation: C.H. and S.P.; Methodology: C.H., S.P., and L.O.; Project administration: C.H., S.P., L.O., and E.R.; Resources: S.P. and L.O.; Software: C.H. and S.P.; Supervision: S.P., L.O., and E.R.; Validation: C.H., S.P., L.O., and E.R.; Visualization: C.H.; Writing – original draft: C.H.; Writing – review and editing: C.H., S.P., L.O., and E.R. All authors reviewed and approved the final version of the manuscript. Competing interests. The authors declare no competing interests. Data availability. The simulated data sets generated for this study can be reproduced using the simulation code (see Code availability) and are also available from the corresponding author upon reasonable request. The Chicago NMMAPS data set is publicly available through the dlnm R package. The malaria data are subject to the data-access conditions of the original study and may be available upon request from the authors of that study. Code availability. The source code implementing DLNAM, together with scripts for the simulation studies and analyses presented in this manuscript, is publicly available on GitHub at https://github.com/shivangpandey31/DLNAM, and is provided without the underlying restricted-access data for the malaria analysis. The version corresponding to this manuscript is archived on Zenodo (DOI: https://doi.org/10.5281/zenodo.22288964). Ethics statement. This methodological study involved no new collection of individual-level participant data. The Chicago NMMAPS data set is publicly available. The malaria analysis used previously collected data, for which the relevant ethical approvals and informed-consent procedures are described in the original study and corresponding data sources. AI use. Generative artificial intelligence tools were used during the development of this work to assist with coding and to improve sentence structure, grammar, and language clarity. All AI-assisted code and text were reviewed and verified by the authors, who take full responsibility for the final manuscript, analyses, and interpretation of the results. Prior dissemination. This article is based on and extends work first presented in C.H.’s master’s thesis at KTH Royal Institute of Technology, carried out at Karolinska Institutet and supervised by S.P., L.O., and E.R. References [1] Shirley Almon “The Distributed Lag Between Capital Appropriations and Expenditures” In Econometrica 33.1, 1965, p. 178–196 DOI: 10.2307/1911894 [2] Antonio Gasparrini et al. “Mortality risk attributable to high and low ambient temperature: a multicountry observational study” In The Lancet 386.9991, 2015, p. 369–375 DOI: 10.1016/S0140-6736(14)62114-0 [3] Krishnan Bhaskaran et al. “Time series regression studies in environmental epidemiology” In International Journal of Epidemiology 42.4, 2013, p. 1187–1195 DOI: 10.1093/ije/dyt092 [4] Antonio Gasparrini, Ben Armstrong and Michael. Kenward “Distributed lag non-linear models” In Statistics in Medicine 29.21, 2010, p. 2224–2234 DOI: 10.1002/sim.3940 [5] Antonio Gasparrini “Modeling exposure-lag-response associations with distributed lag non-linear models” In Statistics in Medicine 33.5, 2014, p. 881–899 DOI: 10.1002/sim.5963 [6] Antonio Gasparrini “Distributed lag linear and non-linear models in R: the package dlnm” In Journal of Statistical Software 43.8, 2011, p. 1–20 DOI: 10.18637/jss.v043.i08 [7] Antonio Gasparrini, Fabian Scheipl, Ben Armstrong and Michael. Kenward “A penalized framework for distributed lag non-linear models” In Biometrics 73.3, 2017, p. 938–948 DOI: 10.1111/biom.12645 [8] Daniel Mork and Ander Wilson “Treed distributed lag nonlinear models” In Biostatistics 23.3, 2022, p. 754–771 DOI: 10.1093/biostatistics/kxaa051 [9] Chao-Yu Guo, Tse-Wei Liu and Yi-Hau Chen “A novel cross-validation strategy for artificial neural networks using distributed-lag environmental factors” In PLOS ONE 16.1, 2021, p. e0244094 DOI: 10.1371/journal.pone.0244094 [10] Peter McCullagh and John. Nelder “Generalized Linear Models” London: Chapman & Hall/CRC, 1989 [11] Trevor. Hastie and Robert. Tibshirani “Generalized Additive Models” London: Chapman & Hall/CRC, 1990 [12] Rishabh Agarwal et al. “Neural Additive Models: Interpretable Machine Learning with Neural Nets” In Advances in Neural Information Processing Systems 34, 2021, p. 4699–4711 [13] Diganta Misra “Mish: A Self Regularized Non-Monotonic Activation Function” In Proceedings of the 31st British Machine Vision Conference (BMVC), 2020 DOI: 10.5244/C.34.191 [14] Stefan Elfwing, Eiji Uchibe and Kenji Doya “Sigmoid-weighted linear units for neural network function approximation in reinforcement learning” In Neural Networks 107, 2018, p. 3–11 DOI: 10.1016/j.neunet.2017.12.012 [15] Balaji Lakshminarayanan, Alexander Pritzel and Charles Blundell “Simple and Scalable Predictive Uncertainty Estimation using Deep Ensembles” In Advances in Neural Information Processing Systems 30, 2017 [16] Erik Daxberger et al. “Laplace Redux – Effortless Bayesian Deep Learning” In Advances in Neural Information Processing Systems 34, 2021, p. 20089–20103 [17] Kouroche Bouchiat et al. “Improving Neural Additive Models with Bayesian Principles” In Proceedings of the 41st International Conference on Machine Learning (ICML) 235, Proceedings of Machine Learning Research PMLR, 2024, p. 4416–4443 [18] Manuel Martellini et al. “Climatic Drivers of Malaria Risk in Children Under Five: A Large-Scale Analysis of Individual-Level Data for 350,000 Children in 26 Sub-Saharan African Countries” Preprint, posted 22 June 2026 In medRxiv, 2026 DOI: 10.64898/2026.06.17.26355855 [19] Yin Lou, Rich Caruana, Johannes Gehrke and Giles Hooker “Accurate Intelligible Models with Pairwise Interactions” In Proceedings of the 19th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD), 2013, p. 623–631 DOI: 10.1145/2487575.2487579 [20] Yin-Hsiu Chen, Bhramar Mukherjee and Veronica. Berrocal “Distributed lag interaction models with two pollutants” In Journal of the Royal Statistical Society Series C: Applied Statistics 68.1, 2019, p. 79–97 DOI: 10.1111/rssc.12297 [21] Theo Economou et al. “A unifying modelling approach for hierarchical distributed lag models” Preprint, posted 18 July 2024 In arXiv, 2024 DOI: 10.48550/arXiv.2407.13374 [22] Malcolm. Mistry and Antonio Gasparrini “Real-time forecast of temperature-related excess mortality at small-area level: towards an operational framework” In Environmental Research: Health 2.3, 2024, p. 035011 DOI: 10.1088/2752-5309/ad5f51 [23] Kaifeng Bi et al. “Accurate medium-range global weather forecasting with 3D neural networks” In Nature 619, 2023, p. 533–538 DOI: 10.1038/s41586-023-06185-3 [24] Jaideep Pathak et al. “FourCastNet: A Global Data-driven High-resolution Weather Model using Adaptive Fourier Neural Operators” Preprint, posted 22 February 2022 In arXiv, 2022 DOI: 10.48550/arXiv.2202.11214 [25] Michael Eisenstein “From quantum computing to mRNA therapeutics: seven technologies to watch in 2026” Technology Feature In Nature 649.8098, 2026, p. 1065–1069 DOI: 10.1038/d41586-026-00188-6 [26] Ilya Loshchilov and Frank Hutter “Decoupled Weight Decay Regularization” In 7th International Conference on Learning Representations (ICLR), 2019 [27] Ilya Loshchilov and Frank Hutter “SGDR: Stochastic Gradient Descent with Warm Restarts” In 5th International Conference on Learning Representations (ICLR), 2017 [28] Robert.. Wedderburn “Quasi-likelihood functions, generalized linear models, and the Gauss–Newton method” In Biometrika 61.3, 1974, p. 439–447 DOI: 10.1093/biomet/61.3.439 [29] Kaiming He, Xiangyu Zhang, Shaoqing Ren and Jian Sun “Delving Deep into Rectifiers: Surpassing Human-Level Performance on ImageNet Classification” In Proceedings of the IEEE International Conference on Computer Vision (ICCV), 2015, p. 1026–1034 DOI: 10.1109/ICCV.2015.123 [30] David.. MacKay “Bayesian Interpolation” In Neural Computation 4.3, 1992, p. 415–447 DOI: 10.1162/neco.1992.4.3.415 [31] David.. MacKay “A Practical Bayesian Framework for Backpropagation Networks” In Neural Computation 4.3, 1992, p. 448–472 DOI: 10.1162/neco.1992.4.3.448 [32] Simon. Wood “Generalized Additive Models: An Introduction with R” Boca Raton: Chapman & Hall/CRC, 2017 [33] Seongwon Im, Ander Wilson and Daniel Mork “Structured Bayesian Regression Tree Models for Estimating Distributed Lag Effects: The R Package dlmtree” In The R Journal 17.1, 2025, p. 136–159 DOI: 10.32614/RJ-2025-007 [34] Mollie. Brooks et al. “glmmTMB Balances Speed and Flexibility Among Packages for Zero-inflated Generalized Linear Mixed Modeling” In The R Journal 9.2, 2017, p. 378–400 DOI: 10.32614/RJ-2017-066 [35] Agustinus Kristiadi, Matthias Hein and Philipp Hennig “Being Bayesian, Even Just a Bit, Fixes Overconfidence in ReLU Networks” In Proceedings of the 37th International Conference on Machine Learning (ICML) 119, Proceedings of Machine Learning Research PMLR, 2020, p. 5436–5446 [36] Alexander Immer, Maciej Korzepa and Matthias Bauer “Improving predictions of Bayesian neural nets via local linearization” In Proceedings of the 24th International Conference on Artificial Intelligence and Statistics (AISTATS) 130, Proceedings of Machine Learning Research PMLR, 2021, p. 703–711 Supplementary Materials S1 Notation Symbol Meaning A. Data and Model YtY_t Outcome at Time t xtx_t, utu_t Exposure Vector; Covariate Vector x⋆x , u⋆u Exposure and Covariate Reference Values ℒmL_m, ℒL Exposure-Specific Lag Set; Collection of Lag Sets xtℒx_t^L Lagged Exposure History ℋtH_t Conditioning History μt _t, ηt _t Conditional Mean; Linear Predictor g Invertible Link Function n, M, K Number of Observations; Number of Exposures; Number of Covariates B. Additive Components α Model Intercept; Reference-Value Linear Predictor fmf_m, hkh_k Exposure–Lag Response Surface; Covariate Shape Function FmF_m Cumulative Exposure–Response Component f~ms f_ms, h~ks h_ks Exposure–Lag Subnetwork; Covariate Subnetwork SmxS_m^x, SkuS_k^u Numbers of Exposure and Covariate Subnetworks C. Parameters, Mixture, and Ensemble θ, ϑ Complete Parameter Vector; Subnetwork Parameters ω Mixing Weights λw _w Weight-Decay Coefficient B Number of Ensemble Members ψ Conditioned Parameters ξ Target-Specific Last-Layer Parameters D. Optimization E Number of Training Epochs η, ηmin _ Learning Rate; Minimum Learning Rate E. Last-Layer Inference ζ Target Evaluation Point τ Inferential Target Φ Last-Layer Predictor Jacobian W GLM Working-Weight Matrix φ Dispersion or Working-Dispersion Factor H Dispersion-Adjusted Last-Layer Information λ, PξP_ξ Gaussian Prior Precision; Penalization Selector Σ^ξ _ξ Conditional Last-Layer Laplace Covariance JζJ_ζ Target Jacobian V Moment-Matched Ensemble Variance F. Simulation Evaluation G Evaluation Grid R Number of Monte Carlo Replicates Table S1: Principal Notation. Recurring notation in the main and supplementary methods; symbols local to individual derivations are defined where they occur. S2 Last-Layer Laplace Confidence Intervals Notation follows the main text. This section gives the full derivation of the target-specific last-layer Laplace construction, its propagation to the reported inferential targets, and its combination across ensemble members. S2.1 Setup, Partition, and Scope For a given inferential target, partition the parameter vector as θ=(ψ,ξ),ξ∈ℝq,θ=(ψ,ξ), ξ ^q, where ξ contains the global intercept and the final-layer parameters of the corresponding component, while ψ contains all remaining parameters: the ExU input layers, hidden layers, other model components and, in the mixture formulation, the mixing weights ω. The conditional calculation fixes ψ at its fitted value ψ ψ, leaving ξ free; ξ ξ denotes its fitted value and serves as the expansion center below. The construction follows the principle of Laplace-approximated Neural Additive Models (LA-NAMs) [17]: form a Gaussian approximation in parameter space and propagate it through the additive components to obtain function-space uncertainty. LA-NAM develops both linearized and last-layer variants for cross-sectional NAM feature networks; DLNAM uses the last-layer construction for distributed-lag components, conditioning on the learned representation and, through ψ ψ, the fitted mixing weights. The partition is chosen so that the conditional model becomes a generalized linear model in the last-layer parameters. Given ψ ψ, each subnetwork contributing to the selected component maps its input to a fixed hidden-feature vector, while all components outside the selected target are fixed. The selected component is therefore a linear combination of fixed features, with the conditioned mixing weights acting only as fixed multipliers, and the predictor is affine in ξ. Conditional curvature and function-space propagation consequently reduce to generalized linear model calculations with a Gaussian prior, avoiding high-dimensional full-network curvature and restricting the local Gaussian approximation to those last-layer parameters [16, 35]. Conditioning on the mixing weights is also required for an affine, identified last-layer parameterization. In the mixture formulation, fmθmx(⋅,⋅)=∑s=1Smxωmsxf~msϑmsx(⋅,⋅),f_m _m^x(·,·)= _s=1^S_m^x _ms^x f_ms _ms^x(·,·), so a mixing weight and the corresponding subnetwork’s final linear layer enter the predictor only through their product. Allowing both to vary within ξ would make the predictor bilinear and scale-non-identifiable: rescaling ωmsx↦cωmsx _ms^x c\, _ms^x while multiplying the corresponding final-layer parameters by c−1c^-1 leaves the fitted function, and hence the likelihood or working likelihood, unchanged. The curvature is therefore singular along this direction. Including ω in ψ and conditioning on ψ ψ removes the indeterminacy at the cost of omitting within-member mixing-weight uncertainty; Section S2.7 returns to this limitation. Write the vector of linear predictors across the n observations as η(ξ)=(η1(ξ),…,ηn(ξ))⊤,η(ξ)= ( _1(ξ),…, _n(ξ) ) , where ηt(ξ)=α+∑m=1M∑ℓ∈ℒmfmψ^,ξ(xt−ℓ,m,ℓ)+∑k=1Khkψ^,ξ(utk), _t(ξ)=α+ _m=1^M _ _mf_m ψ,ξ(x_t- ,m, )+ _k=1^Kh_k ψ,ξ(u_tk), with ηt(ξ)=gμt(ξ) _t(ξ)=g\ _t(ξ)\ and fitted means μ^t=g−1ηt(ξ^). μ_t=g^-1\ _t( ξ)\. Define the last-layer Jacobian Φ=∂η(ξ)∂ξ⊤∈ℝn×q, = ∂η(ξ)∂ξ ^n× q, (S1) which maps perturbations of the last-layer parameters to changes in the linear predictors and plays the same design-matrix role as the cross-basis in a DLNM. Because η(ξ)η(ξ) is affine in ξ, Φ is exact, constant, and independent of the point at which it is evaluated. Throughout, = c= denotes equality, and ≈ c≈ approximation, up to an additive constant that does not depend on the free variable of the expression: ξ in Section S2.2 and λ in Section S2.4. S2.2 Prior and Laplace Approximation Partition the last-layer parameters as ξ=(ξ0,ξP),ξ=( _0, _P), where ξ0 _0 is the global intercept and ξP _P contains the remaining final-layer parameters. Because the intercept is left unpenalized, define Pξ=diag(0,1,…,1)∈ℝq×q.P_ξ=diag(0,1,…,1) ^q× q. The number of penalized parameters is therefore qP=tr(Pξ)=q−1,ξP∈ℝqP,q_P=tr(P_ξ)=q-1, _P ^q_P, and ξ⊤Pξξ=‖ξP‖22.ξ P_ξ=\| _P\|_2^2. For the post-hoc Laplace construction, place an isotropic zero-mean Gaussian prior with precision λ on ξP _P and a flat, improper prior on ξ0 _0, ξP∣λ∼qP(0,λ−1IqP),p(ξ0)∝1. _P λ _q_P\! (0,λ^-1I_q_P ), p( _0) 1. (S2) Thus −logp(ξ∣λ)=λ2ξ⊤Pξξ−qP2logλ.- p(ξ λ)\; c=\; λ2ξ P_ξ- q_P2 λ. Leaving the intercept unpenalized lets the data determine the baseline response. Its improper normalization is independent of λ and therefore does not affect its selection in Section S2.4. The term −qPlog(λ)/2-q_P (λ)/2 affects neither the mode nor the curvature but must be retained in the evidence because it depends on λ. For a likelihood fit, write ℓ(ξ)=logp(Y∣ξ,ψ^) (ξ)= p(Y ξ, ψ) for the conditional log-likelihood. By Bayes’ rule, p(ξ∣Y,ψ^,λ)∝p(Y∣ξ,ψ^)p(ξ∣λ),p(ξ Y, ψ,λ) p(Y ξ, ψ)\,p(ξ λ), so the negative log-posterior is Ψ(ξ)=−ℓ(ξ)−logp(ξ∣λ). (ξ)\; c=\;- (ξ)- p(ξ λ). (S3) For a working-likelihood fit, the same expression is used with ℓ interpreted as the corresponding conditional working criterion; the resulting construction is therefore a Gaussian working analogue rather than a literal posterior from a normalized likelihood. The precision λ λ is selected after network training by the evidence-style procedure of Section S2.4. The fitted ξ ξ is therefore retained rather than re-optimized under the Gaussian prior and is treated as an approximate posterior mode for likelihood fits, or as the corresponding working-mode analogue for working-likelihood fits. Thus, ∇ξΨ(ξ^)≈0, _ξ ( ξ)≈ 0, with accuracy determined by how closely the trained solution satisfies the post-hoc penalized score condition at λ λ. Unlike the alternating parameter–hyperparameter optimization of LA-NAM [17], the present construction holds ξ ξ and H fixed while selecting λ, leaving the reported point estimate unchanged. Expanding about ξ ξ, Ψ(ξ)=Ψ(ξ^)+∇ξΨ(ξ^)⊤(ξ−ξ^)+12(ξ−ξ^)⊤ℐξ(ξ−ξ^)+(‖ξ−ξ^‖3),ℐξ=∇ξ2Ψ(ξ^). (ξ)= ( ξ)+ _ξ ( ξ) (ξ- ξ)+ 12(ξ- ξ) I_ξ(ξ- ξ)+O (\|ξ- ξ\|^3 ), _ξ= _ξ^2 ( ξ). (S4) The remainder is of third order under the usual smoothness of the conditional likelihood or working criterion as a function of the linear predictor. Once ψ ψ is fixed and η is affine in ξ, this condition places no further smoothness requirement on the network representation itself. Under the approximate-mode assumption, the linear term is neglected and the quadratic expansion defines the member-level Gaussian Laplace approximation ξ|Y,ψ^≈(ξ^,Σ^ξ),Σ^ξ=ℐξ−1,ξ Y, ψ\ ≈\ N\! ( ξ, _ξ ), _ξ=I_ξ^-1, for likelihood fits [31, 16]; for working-likelihood fits, the same expression is interpreted as its Gaussian working analogue. S2.3 Last-Layer Curvature Because the likelihood or working criterion depends on ξ only through η(ξ)η(ξ), the chain rule gives ∇ξ(−ℓ)=Φ⊤∇η(−ℓ), _ξ(- )= _η(- ), and ∇ξ2(−ℓ)=Φ⊤[∇η2(−ℓ)]Φ+∑t=1n∂(−ℓ)∂ηt∇ξ2ηt. _ξ^2(- )= [ _η^2(- ) ] + _t=1^n ∂(- )∂ _t _ξ^2 _t. Because η(ξ)η(ξ) is affine in ξ, ∇ξ2ηt=0 _ξ^2 _t=0 for every t, so no predictor-curvature term enters the Hessian. The second term above therefore vanishes identically, and no Gauss–Newton approximation is required for the parameterization used here. Had ξ entered nonlinearly, as it would if the mixing weights were allowed to vary jointly with their corresponding output layers, this term would generally be nonzero; discarding it on the basis of its vanishing expectation would then constitute a Gauss–Newton or Fisher approximation [36]. For a generalized linear model [10], let W=diag(w1,…,wn),wt=(dμt/dηt)2V(μt),W=diag(w_1,…,w_n), w_t= (d _t/d _t)^2V( _t), with the weights evaluated at the fitted means and V denoting the variance function. To see how these weights enter the curvature, for an individual observation the likelihood or quasi-likelihood score with respect to the linear predictor has the form ∂ℓt∂ηt=Yt−μtφV(μt)dμtdηt. ∂ _t∂ _t= Y_t- _t V( _t) d _td _t. Differentiating once more gives a term −1φ(dμt/dηt)2V(μt)- 1 (d _t/d _t)^2V( _t) together with terms proportional to Yt−μtY_t- _t. Because (Yt−μt∣μt)=0,E(Y_t- _t _t)=0, the latter vanish in expectation, yielding −[∂2ℓt∂ηt2|μt]=wtφ,-E [ ∂^2 _t∂ _t^2\, |\, _t ]= w_t , or, across observations, −[∇η2ℓ∣μ]=W/φ.-E [ _η^2 μ ]=W/ . Thus W/φW/ is the expected information on the linear-predictor scale, and Φ pulls this information back to the last-layer parameterization. For likelihood fits with canonical links, the residual-dependent terms in the observed second derivative vanish and the observed and expected information coincide. For non-canonical likelihood fits, W/φW/ gives the corresponding Fisher information; for working-likelihood fits it gives the analogous working information. Let φ denote the dispersion or working-dispersion factor appropriate to the fitted outcome model, taking its family-specified value when fixed and being estimated where required. The dispersion-adjusted last-layer information is therefore H=Φ⊤WΦ/φ^.H= W / . (S5) When dispersion is fixed by the fitted family, its family-specified value is used. When it is estimated in the reported analyses, we use a Pearson-residual estimate. For each ensemble member b=1,…,Bb=1,…,B, with fitted means μ^bt μ_bt, φ^b=1ν∑t=1n(Yt−μ^bt)2V(μ^bt),φ^=max1,1B∑b=1Bφ^b. _b= 1ν _t=1^n (Y_t- μ_bt)^2V( μ_bt), = \1,\, 1B _b=1^B _b \. (S6) The lower bound φ^≥1 ≥ 1 prevents apparent underdispersion from understating uncertainty. We use the pragmatic large-sample approximation ν=n−1ν=n-1. This does not treat the conditioned representation as parameter-free: the effective last-layer dimension γ derived below quantifies the penalized last-layer directions supported by the data, not the residual degrees of freedom of the complete trained network, and is therefore not substituted for ν. An exact full-network correction is unavailable, and the distinction is negligible at the sample sizes considered here. Differentiating the Gaussian prior twice gives ∇ξ[−logp(ξ∣λ)]=λPξξ,∇ξ2[−logp(ξ∣λ)]=λPξ. _ξ[- p(ξ λ)]=λ P_ξ, _ξ^2[- p(ξ λ)]=λ P_ξ. A Gaussian prior therefore contributes its precision additively to the conditional curvature. At the selected precision λ λ, and with a small numerical ridge ϵ>0ε>0, ℐξ=H+λ^Pξ+ϵI,Σ^ξ=ℐξ−1.I_ξ=H+ λP_ξ+ε I, _ξ=I_ξ^-1. (S7) This is the conditional last-layer covariance given ψ ψ, with Φ playing the role of the design matrix. In implementation, the computation is organized on the unscaled information as Σ^ξ=φ^(Φ⊤WΦ+φ^λ^Pξ+ϵI)−1, _ξ= ( W + λP_ξ+ε I )^-1, which is equivalent to Equation (S7) up to the scaling convention for the numerical ridge. Factoring out the dispersion rescales the prior term to φ^λ^Pξ λP_ξ in this implementation form, while ϵε serves only as a numerical safeguard. Although Φ⊤WΦ W may be rank-deficient because constant columns induced by final-layer biases are collinear with the global intercept, the prior supplies curvature in penalized directions and the ridge additionally protects against numerical ill-conditioning, including in the unpenalized intercept direction. S2.4 Evidence-Style Prior Precision For likelihood fits, λ is selected by a post-hoc approximation to type-I maximum likelihood, or empirical Bayes. The target criterion is the evidence for λ, with the last-layer parameters integrated out, p(Y∣λ)∝∫λp(Y∣ξ,ψ^)p(ξ∣λ)ξ,p(Y λ)\; _λ\; p(Y ξ, ψ)p(ξ λ)\,dξ, rather than the prior density evaluated at ξ ξ, whose maximization over λ would make no reference to the data. Here ∝λ _λ denotes equality up to a multiplicative constant independent of λ, including the arbitrary normalization of the flat intercept prior. In the implemented post-hoc calculation, the trained expansion center ξ ξ and working information H are held fixed as λ is updated. The procedure is therefore a MacKay evidence fixed-point approximation rather than literal maximization of a fully re-optimized marginal likelihood [30, 31]. For working-likelihood fits, where no normalized likelihood and hence no literal marginal likelihood is available, the same fixed-point calculation is applied to the working information and interpreted as the corresponding evidence-style analogue. Using ℐξ(λ)=H+λPξ,I_ξ(λ)=H+λ P_ξ, with the numerical ridge omitted for selection of λ, the λ-dependent terms of the Laplace log evidence under the fixed-ξ ξ, fixed-H approximation are logp(Y∣λ)≈ℓ(ξ^)−λ2‖ξ^P‖22+qP2logλ−12log|ℐξ(λ)|. p(Y λ)\; c≈\; ( ξ)- λ2\| ξ_P\|_2^2+ q_P2 λ- 12 |I_ξ(λ)|. The three λ-dependent terms balance the Gaussian-prior normalization, the compatibility of the prior scale with the fitted final-layer parameters, and the log-determinant complexity term arising from the Gaussian integral. Differentiating with respect to λ and using ∂λlog|ℐξ(λ)|=tr[ℐξ(λ)−1Pξ] ∂λ |I_ξ(λ)|=tr [I_ξ(λ)^-1P_ξ ] gives −12‖ξ^P‖22+qP2λ−12tr[ℐξ(λ)−1Pξ]=0.- 12\| ξ_P\|_2^2+ q_P2λ- 12tr [I_ξ(λ)^-1P_ξ ]=0. The dependence of the expansion center and working information on λ is omitted by construction. Multiplying by 2λ2λ gives λ‖ξ^P‖22=γ(λ),λ\| ξ_P\|_2^2=γ(λ), where γ(λ)=qP−λtr[(H+λPξ)−1Pξ]=∑i=1qPeiei+λ∈[0,qP),γ(λ)=q_P- [(H+λ P_ξ)^-1P_ξ ]= _i=1^q_P e_ie_i+λ∈[0,q_P), (S8) and hence the MacKay fixed point is λ^=γ(λ^)‖ξ^P‖22. λ= γ( λ)\| ξ_P\|_2^2. (S9) Here e1,…,eqPe_1,…,e_q_P are the eigenvalues of the penalized information after profiling out the unpenalized intercept. Partitioning H into intercept (00) and penalized (P) components gives H+λPξ=(H00H0PHP0HPP+λI).H+λ P_ξ= pmatrixH_00&H_0P\\ H_P0&H_P+λ I pmatrix. The raw block HPPH_P does not by itself represent the information available for the penalized parameters when the intercept is allowed to adjust, because changes in the penalized coefficients may be partially offset by changes in the intercept. Profiling out the intercept removes this shared direction, leaving the Schur complement HP=HPP−HP0H00−1H0P.H_P=H_P-H_P0H_00^-1H_0P. (S10) Thus HPH_P is the information for the penalized parameters after accounting for their dependence on the unpenalized intercept. By the block-inverse identity, the penalized block of (H+λPξ)−1(H+λ P_ξ)^-1 is (HP+λI)−1,(H_P+λ I)^-1, and consequently tr[(H+λPξ)−1Pξ]=tr[(HP+λI)−1].tr [(H+λ P_ξ)^-1P_ξ ]=tr [(H_P+λ I)^-1 ]. The raw submatrix HPPH_P would suffice only if the intercept were orthogonal to the penalized directions, which reference centering does not ensure. Equation (S8) identifies γ as the effective number of penalized directions constrained by the data: each eigendirection contributes ei/(ei+λ)e_i/(e_i+λ), approaching one when the data dominate the prior and zero when the prior dominates. If the log-determinant term were omitted, the fixed point would reduce to λ=qP/‖ξ^P‖22,λ=q_P/\| ξ_P\|_2^2, treating all qPq_P penalized directions as constrained. Retaining the determinant replaces qPq_P by the effective count γ. Equivalently, λ^−1=‖ξ^P‖22γ λ^-1= \| ξ_P\|_2^2γ estimates the prior variance as the mean squared fitted weight per effective penalized direction, structurally analogous to degrees-of-freedom corrections in linear regression and to marginal-likelihood selection of smoothing parameters in penalized spline models. Under the fixed expansion-center and fixed-information approximation, the eie_i are computed once from HPH_P and do not vary with λ. To characterize the solution, define F(λ)=λ‖ξ^P‖22−γ(λ).F(λ)=λ\| ξ_P\|_2^2-γ(λ). In the nondegenerate case with ‖ξ^P‖2>0\| ξ_P\|_2>0 and information in at least one penalized direction, F′(λ)=‖ξ^P‖22+∑i=1qPei(ei+λ)2>0,F (λ)=\| ξ_P\|_2^2+ _i=1^q_P e_i(e_i+λ)^2>0, while F(λ)F(λ) is negative as λ↓0λ 0 and diverges to +∞+∞ as λ→∞λ→∞. The stationarity equation therefore has a unique positive solution. Equation (S9) is solved by the iteration λ(r+1)=γ(λ(r))‖ξ^P‖22,r=0,1,…,λ^(r+1)= γ(λ^(r))\| ξ_P\|_2^2, r=0,1,…, from an initial λ(0)>0λ^(0)>0 until numerical convergence. For an ensemble, γ and ‖ξ^P‖22\| ξ_P\|_2^2 are summed across members before the update, yielding one shared pooled fixed-point rule across representations fitted to the same data. This pooling is not a product of independent marginal likelihood contributions. S2.5 Inferential Targets and Propagation The inferential targets are the centered components of Equation (1): the exposure–lag surface fmψ^,ξf_m ψ,ξ, its cumulative Fmψ^,ξ(⋅)=∑ℓ∈ℒmfmψ^,ξ(⋅,ℓ),F_m ψ,ξ(·)= _ _mf_m ψ,ξ(·, ), and the covariate shape function hkψ^,ξh_k ψ,ξ. Under the reference-centering convention, fmψ^,ξ(xm⋆,ℓ)=0for all ℓ∈ℒm,hkψ^,ξ(uk⋆)=0,f_m ψ,ξ(x_m , )=0 all _m, h_k ψ,ξ(u_k )=0, so these functions are contrasts on the link scale. For τψ^,ξ∈fmψ^,ξ,Fmψ^,ξ,hkψ^,ξ,τ ψ,ξ∈ \f_m ψ,ξ,F_m ψ,ξ,h_k ψ,ξ \, let τψ^,ξ _G ψ,ξ denote its evaluations on a finite grid G, and define J=∂τψ^,ξ∂ξ⊤,Jζ=∂τψ^,ξ(ζ)∂ξ⊤,ζ∈.J_G= ∂ _G ψ,ξ∂ξ , J_ζ= ∂τ ψ,ξ(ζ)∂ξ , ζ . Conditional on ψ ψ, every target is affine in ξ, so τψ^,ξ=τψ^,ξ^+J(ξ−ξ^) _G ψ,ξ= _G ψ, ξ+J_G(ξ- ξ) holds exactly and J_G is constant. An affine transformation of a Gaussian is Gaussian; hence Σ^=JΣ^ξJ⊤,sdτψ^,ξ(ζ)=JζΣ^ξJζ⊤1/2,ζ∈. _G=J_G _ξJ_G , \τ ψ,ξ(ζ)\= \J_ζ _ξJ_ζ \^1/2, ζ . (S11) Thus covariance propagation is exact conditional on ψ ψ, with no delta-method linearization error. Because the inferential targets are reference-centered contrasts, the intercept column of J_G is zero, although the intercept remains in ℐξI_ξ. Evaluating the conditional target at the fitted expansion center gives τ^ψ^,ξ^(ζ)≡τψ^,ξ(ζ)|ξ=ξ τ ψ, ξ(ζ)≡ .τ ψ,ξ(ζ) |_ξ= ξ. The member-level pointwise 1−α1-α interval on the link-contrast scale is therefore τ^ψ^,ξ^(ζ)±z1−α/2sdτψ^,ξ(ζ). τ ψ, ξ(ζ)± z_1-α/2sd\τ ψ,ξ(ζ)\. Where the reported estimand is a monotone transformation of the link-scale contrast, the interval endpoints are transformed accordingly. When an absolute linear predictor is reported on the mean scale, the transformation is g−1g^-1. In the reported analyses, exponentiation therefore maps log-link contrasts to relative risks and logit-link contrasts to odds ratios, while the identity link requires no transformation. S2.6 Ensembles For an ensemble of B independently initialized members [15], write ψ^=(ψ^1,…,ψ^B),ξ^=(ξ^1,…,ξ^B) ψ=( ψ_1,…, ψ_B), ξ=( ξ_1,…, ξ_B) for the collections of fitted member-specific conditioned parameters and last-layer expansion centers, respectively. Applying the ensemble construction of Equation (3) to a reported inferential target gives τ^ψ^,ξ^(ζ)=1B∑b=1Bτ^bψ^b,ξ^b(ζ),^(ζ)=1B∑b=1BJbζΣ^ξbJbζ⊤⏟within-member+1B∑b=1Bτ^bψ^b,ξ^b(ζ)−τ^ψ^,ξ^(ζ)2⏟between-member. τ ψ, ξ(ζ)= 1B _b=1^B τ_b ψ_b, ξ_b(ζ), V(ζ)= 1B _b=1^BJ_bζ _ _bJ_bζ _within-member+ 1B _b=1^B \ τ_b ψ_b, ξ_b(ζ)- τ ψ, ξ(ζ) \^2_between-member. (S12) Combining conditional last-layer Laplace uncertainty with between-member variation, the reported pointwise 1−α1-α confidence interval is τ^ψ^,ξ^(ζ)±z1−α/2^(ζ)1/2, τ ψ, ξ(ζ)± z_1-α/2 V(ζ)^1/2, with the same transformation to the reported scale as above. The first term in Equation (S12) averages member-specific conditional last-layer variances without further division by B. Dividing it by B would instead estimate the variance of an arithmetic mean under an independence interpretation and would mechanically narrow the interval as more initializations were fitted, whereas the present construction represents an equally weighted mixture of member-level conditional approximations. The second term records between-member variation; its 1/B1/B denominator is the variance convention of that finite mixture, not an unbiased estimate of a hypothetical population variance over random initializations. Consequently, τ^ψ^,ξ^(ζ) τ ψ, ξ(ζ) and ^(ζ) V(ζ) are the mean and variance of the equally weighted mixture of member-level Gaussian Laplace approximations. The reported interval is a Gaussian moment-matched summary of this generally non-Gaussian mixture. Because the members are fitted to the same data, ^(ζ) V(ζ) is a constructed uncertainty measure, not the sampling variance of the ensemble mean or full posterior uncertainty in the learned representation. Its between-member term captures only fit instability expressed through disagreement among independently initialized representations and last-layer centers. S2.7 Regularity Conditions and Scope The procedure as a whole is approximate, but conditional on ψ ψ several parts of the calculation are exact. Because the predictor is affine in ξ, the Jacobian Φ is exact and constant and the predictor-curvature term in the Hessian vanishes identically, so no Gauss–Newton approximation is required. The inferential targets are likewise affine in ξ, making the propagation in Equation (S11) exact with no delta-method linearization error. For likelihood fits with canonical links, the likelihood contribution to the last-layer curvature is also exact. The remaining assumptions and approximations are: (i) the trained ξ ξ, which is not re-optimized after selection of λ λ, is treated as an approximate posterior mode for likelihood fits, or as the corresponding working-mode analogue for working-likelihood fits, so that ∇ξΨ(ξ^)≈0 _ξ ( ξ)≈ 0; (i) ℐξI_ξ is positive definite, secured numerically by the prior together with the ridge; (i) the conditional likelihood or working criterion is three times continuously differentiable in the linear predictor over the relevant domain, ensuring the third-order expansion in Equation (S4); (iv) for non-canonical likelihood fits, the Fisher information adequately represents the observed information; (v) for working-likelihood fits, the dispersion-scaled working information provides an adequate curvature surrogate; and (vi) at the ensemble level, the equally weighted member mixture is adequately summarized by the Gaussian interval obtained by matching its mean and variance. The member-level Laplace approximations are conditional on ψ^b ψ_b, including the fitted mixing weights in the mixture formulation, and therefore do not constitute a posterior over the complete parameter vector. Relaxing this conditioning would require the linearized Laplace variant, which linearizes the complete network about the trained parameters and so propagates uncertainty in ψ ψ as well [36, 17]. For distributed-lag components this would require Jacobians with respect to all network parameters at every evaluated exposure–lag point, and a correspondingly larger and denser information matrix; the last-layer construction is retained here for that reason. The reported ensemble confidence intervals additionally incorporate between-member variation, which captures fit instability not represented by the member-level conditional variance, including that arising from differences in learned representations and fitted last-layer centers, but only insofar as it appears as disagreement among independently initialized fits of the same data. Their frequentist adequacy is consequently an empirical question and is evaluated by coverage against the generating functions in the simulation studies. The reported intervals are pointwise confidence intervals for the inferential targets fmf_m, FmF_m, and hkh_k, not simultaneous confidence bands; statements about an entire curve or surface would require an additional multiplicity adjustment. They quantify uncertainty in the fitted response functions rather than prediction of future outcomes. For a fully specified outcome distribution, prediction intervals could instead be obtained from quantiles of the predictive distribution, additionally integrating over uncertainty in the fitted linear predictor. A mean–variance specification alone, as in quasi-likelihood models, does not determine predictive quantiles and would require an additional distributional assumption. S3 Data-Generating Surfaces The outcome model, exposure process, reference value, and evaluation grid are defined in Section 4.5. Let σ(z)=1+exp(−z)−1,ρc,δ(z)=exp−12(z−cδ)2,σ(z)=\1+ (-z)\^-1, _c,δ(z)= \- 12 ( z-cδ )^2 \, denote the logistic sigmoid and Gaussian bump, respectively, and define the normalized lag kernels Dκ(ℓ)=exp(−ℓ/κ)∑ℓ′=0Lexp(−ℓ′/κ),Qc,δ(ℓ)=ρc,δ(ℓ)∑ℓ′=0Lρc,δ(ℓ′).D_κ( )= (- /κ) _ =0^L (- /κ), Q_c,δ( )= _c,δ( ) _ =0^L _c,δ( ). The four data-generating surfaces were f1(x,ℓ) f_1(x, ) =0.50σ((15−x)/2.5)+0.60σ((x−25)/1.5)D3(ℓ), = \0.50\,σ((15-x)/2.5)+0.60\,σ((x-25)/1.5) \D_3( ), f2(x,ℓ) f_2(x, ) =0.70σ(x−30)Q1,1.5(ℓ)+σ((12−x)/3)0.70Q9,2(ℓ)−0.35Q5,1.5(ℓ), =0.70\,σ(x-30)\,Q_1,1.5( )+σ((12-x)/3) \0.70\,Q_9,2( )-0.35\,Q_5,1.5( ) \, f3(x,ℓ) f_3(x, ) =0.70ρ32,2(x)Q1,1.5(ℓ)+0.40σ((10−x)/3)D4(ℓ), =0.70\, _32,2(x)\,Q_1,1.5( )+0.40\,σ((10-x)/3)\,D_4( ), f4(x,ℓ) f_4(x, ) =0.70σ(x−28)Q1,1(ℓ)+0.60σ(12−x)Q9,1.5(ℓ). =0.70\,σ(x-28)\,Q_1,1( )+0.60\,σ(12-x)\,Q_9,1.5( ). As throughout, the generating surfaces are reported and scored after centering at x⋆=20x =20. DGP 1 is separable: a smooth exposure response is multiplied by a single exponentially decaying lag kernel. DGP 2 combines a sharp high-exposure response acting at early lags with a low-exposure response acting later and an opposing-sign contribution at intermediate lags. DGP 3 combines a narrow high-exposure response with a broader low-exposure arm on a slower lag kernel. DGP 4 places sharp responses at both extremes of the exposure range, with different dominant lag timing. Figure S1 shows the resulting exposure–lag surfaces and cumulative curves, the latter being the primary estimand of the simulation. Figure S1: Simulation Study: Data-Generating Processes. Cumulative relative-risk functions and corresponding exposure–lag relative-risk surfaces for the four DGPs. S4 Monte Carlo Uncertainty Monte Carlo uncertainty is calculated from replicate-level summaries, thereby accounting for dependence among evaluation-grid points within each replicate rather than treating them as independent Monte Carlo replicates. For any reported summary that is an average of replicate-wise quantities, S^=S¯=1R∑r=1RSr, S= S= 1R _r=1^RS_r, the Monte Carlo standard error is MCSE(S^)=1R1R−1∑r=1R(Sr−S¯)2.MCSE( S)= 1 R 1R-1 _r=1^R(S_r- S)^2. The plotted 95%95\% Monte Carlo confidence interval is then S^±1.96MCSE(S^). S± 1.96\,MCSE( S). For pointwise coverage of target τ over evaluation region τA_τ, define the replicate-level coverage summary Crτ(τ)=1|τ|∑ζ∈τ(ζ)∈Irτ(ζ).C_r^τ(A_τ)= 1|A_τ| _ζ _τ1 \τ(ζ)∈ I_r^τ(ζ) \. Then Covτ(τ)=1R∑r=1RCrτ(τ),Cov_τ(A_τ)= 1R _r=1^RC_r^τ(A_τ), and its MCSE follows directly from the replicate-level formula above. This construction preserves arbitrary dependence among coverage indicators across evaluation points within the same fitted replicate. For RMSE, let Mrτ(τ)=1|τ|∑ζ∈τdrτ(ζ)2M_r^τ(A_τ)= 1|A_τ| _ζ _τ\d_r^τ(ζ)\^2 denote the replicate-wise mean squared error and let M¯τ(τ)=1R∑r=1RMrτ(τ), M^τ(A_τ)= 1R _r=1^RM_r^τ(A_τ), so that RMSE^τ(τ)=M¯τ(τ). RMSE_τ(A_τ)= M^τ(A_τ). The replicate-level formula gives MCSEM¯τ(τ)=1R1R−1∑r=1RMrτ(τ)−M¯τ(τ)2,MCSE \ M^τ(A_τ) \= 1 R 1R-1 _r=1^R \M_r^τ(A_τ)- M^τ(A_τ) \^2, and first-order delta-method propagation through the square-root transformation gives MCSERMSE^τ(τ)≈MCSEM¯τ(τ)2M¯τ(τ).MCSE \ RMSE_τ(A_τ) \≈ MCSE \ M^τ(A_τ) \2 M^τ(A_τ). For the joint-to-single degradation ratio D^=RMSE^jointRMSE^single, D= RMSE_joint RMSE_single, the implemented first-order propagation sets the covariance between the two Monte Carlo RMSE estimates to zero, giving MCSE(D^)≈D^MCSE(RMSE^joint)RMSE^joint2+MCSE(RMSE^single)RMSE^single2.MCSE( D)≈ D \ MCSE( RMSE_joint) RMSE_joint \^2+ \ MCSE( RMSE_single) RMSE_single \^2. No covariance term is included in this ratio approximation. Regional squared bias and variance depend on the complete collection of replicates through d¯τ(ζ) d^τ(ζ), rather than being simple averages of fixed replicate-wise summaries. Their Monte Carlo uncertainty is therefore estimated by delete-one-replicate jackknifing. For T∈Biasτ2(τ),Varτ(τ),T∈ \Bias_τ^2(A_τ),Var_τ(A_τ) \, let T(−r)T_(-r) denote the estimate obtained after omitting replicate r, and define T¯(−⋅)=1R∑r=1RT(−r). T_(-·)= 1R _r=1^RT_(-r). The jackknife Monte Carlo standard error is MCSEJK(T)=R−1R∑r=1RT(−r)−T¯(−⋅)2.MCSE_JK(T)= R-1R _r=1^R \T_(-r)- T_(-·) \^2. These calculations are applied separately to the cumulative and full-surface targets, τ∈F,fτ∈\F,f\, and to each reported evaluation region τA_τ. S5 Mixture Subnetwork Decomposition Each Mixture DLNAM component is a learned linear combination of S separately initialized, jointly trained subnetworks, fm(⋅,⋅)=∑s=1Sωmsf~ms(⋅,⋅).f_m(·,·)= _s=1^S _ms\, f_ms(·,·). The ablation of Section 2 removes this mixture construction and returns to the standard single-network formulation of Equation (1); although recovery degrades, the ablation does not show how the subnetworks combine within a fitted component. Figure S2 provides that decomposition for one independently initialized ensemble member fitted to one replicate of each DGP, without ensemble averaging. The decomposition indicates that the mixture distributes functional complexity across weighted subnetwork contributions rather than requiring any single subnetwork to reproduce the complete exposure–response function. Individual contributions can emphasize different features of the response and combine through reinforcement or partial cancellation on the log-relative-risk scale. Together with the ablation results, this supports a representational role for the mixture construction beyond simple parameter expansion. The decomposition is illustrative rather than representative and does not imply that particular subnetwork roles are identifiable or stable across fits. Figure S2: Mixture DLNAM: Subnetwork Decomposition. Cumulative relative-risk functions for the fitted Mixture DLNAM and its three weighted subnetwork contributions. S6 Multivariate ExU Input Layers The scalar ExU of Equation (4) maps one input coordinate to djd_j units, so components with p>1p>1 inputs require a multivariate extension; for a distributed-lag surface, z=(x,ℓ)⊤z=(x, ) with p=2p=2. We compare the reported coordinate-wise concatenation of Equation (5) with two unified alternatives that allow coordinate mixing within the input layer. Throughout, ϕφ acts elementwise, ⊙ denotes the Hadamard product, k1_k is the all-ones vector of length k, and exponentiation of vector- or matrix-valued quantities is elementwise; eWe^W does not denote a matrix exponential. Components specified with linear inputs instead use the standard input layer, which mixes all coordinates in a single affine map, hLin(z)=ϕ(Wz+b),W∈ℝd×p,b∈ℝd,h^Lin(z)=φ(Wz+b), W ^d× p, b ^d, with dp+ddp+d input-layer parameters and no unit-specific localization. It is not among the layers compared below, which concern the multivariate ExU parameterization used for exposure–lag surfaces. The concatenation layer applies an independent scalar ExU to each coordinate and stacks the resulting outputs, hC(z)=(h1(z1)⊤,…,hp(zp)⊤)⊤∈ℝd,hj(zj)=ϕ(ewj⊙(zjdj−bj)),h^C(z)= (h_1(z_1) ,…,h_p(z_p) ) ^d, h_j(z_j)=φ (e^w_j (z_j1_d_j-b_j) ), where wj,bj∈ℝdjw_j,b_j ^d_j and d=∑j=1pdjd= _j=1^pd_j. For an exposure–lag surface, the width is split as dx=⌈d/2⌉d_x= d/2 and dℓ=⌊d/2⌋d_ = d/2 . Each first-layer unit therefore depends on only one coordinate, with exposure–lag interaction formed in subsequent hidden layers. The two unified layers instead allow every first-layer unit to depend on all coordinates. With log-weights W∈ℝp×dW ^p× d, the unified local-bias layer assigns a separate location to each coordinate–unit pair through B∈ℝp×dB ^p× d, hL(z)=ϕ([eW⊙(zd⊤−B)]⊤p)∈ℝd.h^L(z)=φ ( [e^W (z1_d -B) ] 1_p ) ^d. This retains the unit-specific localization of the scalar ExU and reduces exactly to Equation (4) when p=1p=1. The unified shared-bias layer instead uses one location per coordinate, b∈ℝpb ^p, shared across units, hS(z)=ϕ((eW)⊤(z−b))∈ℝd.h^S(z)=φ ((e^W) (z-b) ) ^d. When p=1p=1, all units share a single location and the scalar ExU is not recovered; the construction is included as a simpler unified multivariate parameterization with one shared location per coordinate. The three layers are compared at matched output width d, not matched parameter count. Their input-layer weights and biases number 2d2d, 2pd2pd, and pd+ppd+p for concatenation, unified local-bias, and unified shared-bias, respectively. Input scaling and initialization settings otherwise match the reported configuration of Section 4.7. Figure S3: Simulation Study: ExU Comparison. Monte Carlo simulation of four DGPs, with R=50R=50 replicates per DGP. (A) Estimated cumulative relative-risk functions and data-generating functions, shown as pointwise Monte Carlo means across fitted replicates. (B) Cumulative log-relative-risk RMSE by method and exposure region, partitioned into squared bias and variance. (C) Empirical coverage of pointwise 95%95\% intervals. Error bars denote 95% Monte Carlo confidence intervals. Boundary denotes exposure values below the 5th or above the 95th percentile of the reference exposure distribution; interior denotes its complement. Associated RMSE and coverage results are reported in Table S8. Across DGPs 1–4 at R=50R=50 replicates, the reported concatenation layer gives the lowest error in every DGP for both the cumulative and full-surface targets (Figs. S3 and S6; Tables S8 and S13). The unified shared-bias layer, which removes unit-specific localization by sharing one location per coordinate across units, performs worst, giving the highest error on both the cumulative and full-surface targets and the lowest cumulative coverage. The unified local-bias layer retains unit-specific localization while allowing coordinate mixing in the input layer and gives intermediate cumulative error. These results indicate that the subsequent hidden layers are sufficient to represent the exposure–lag interactions present in the tested DGPs, while introducing multivariate coordinate mixing within the ExU input layer provides no recovery advantage. The comparison therefore supports the coordinate-wise concatenation used in the reported DLNAM against these two unified alternatives, without implying that it is universally optimal. S7 Prospective Extensions S7.1 Multitask DLNAM Let Yt∈OY_t ^O denote the outcomes and define μt=(Yt∣ℋt) _t=E(Y_t _t). For an invertible componentwise link function g, the Multitask DLNAM predictor is go(μto)=αo+∑m=1M∑ℓ∈ℒmfmθmox(xt−ℓ,m,ℓ)+∑k=1Khkθkou(utk),g_o( _to)= _o+ _m=1^M _ _mf_m _mo^x(x_t- ,m, )+ _k=1^Kh_k _ko^u(u_tk), (S13) where the outcome-specific components are fmθmox(⋅,⋅)=∑s=1Smxωmsoxf~msϑmsx(⋅,⋅),hkθkou(⋅)=∑s=1Skuωksouh~ksϑksu(⋅).f_m _mo^x(·,·)= _s=1^S_m^x _mso^x f_ms _ms^x(·,·), h_k _ko^u(·)= _s=1^S_k^u _kso^u h_ks _ks^u(·). (S14) The subnetwork parameters ϑ are shared, while the intercepts αo _o and mixing weights ωo=(ωox,ωou) _o=( _o^x, _o^u) are outcome-specific. This is the distributed-lag analogue of the multitask NAM construction [12]. For one outcome, Equations (S13) and (S14) reduce to the Mixture DLNAM. The shared representation induces cross-outcome regularization while limiting outcome-specific component flexibility to the mixing weights; it does not itself specify residual dependence. S7.2 Hierarchical DLNAM For grouped outcomes Ytc∈Y_tc , define μtc=(Ytc∣ℋt) _tc=E(Y_tc _t). For an invertible link function g, the Hierarchical DLNAM predictor is g(μtc)=α+αc+∑m=1M∑ℓ∈ℒmfmθmx(xt−ℓ,m,ℓ)+fmcθmcx(xt−ℓ,m,ℓ)+∑k=1Khkθku(utk)+hkcθkcu(utk),g( _tc)=α+ _c+ _m=1^M _ _m \f_m _m^x(x_t- ,m, )+f_mc _mc^x(x_t- ,m, ) \+ _k=1^K \h_k _k^u(u_tk)+h_kc _kc^u(u_tk) \, (S15) where the first component in each brace is global and follows Equation (2), while the group-specific deviations are fmcθmcx(⋅,⋅)=∑s=1Smxωmscxf~mscϑmscx(⋅,⋅),hkcθkcu(⋅)=∑s=1Skuωkscuh~kscϑkscu(⋅).f_mc _mc^x(·,·)= _s=1^S_m^x _msc^x f_msc _msc^x(·,·), h_kc _kc^u(·)= _s=1^S_k^u _ksc^u h_ksc _ksc^u(·). (S16) Partial pooling is induced by ωmscx∼(0,(σmsx)2),ωkscu∼(0,(σksu)2),αc∼(0,σα2). _msc^x (0,( _ms^x)^2 ), _ksc^u (0,( _ks^u)^2 ), _c (0, _α^2). (S17) This is the global-plus-deviation construction of hierarchical DLNMs [21], expressed through global and group-specific neural subnetworks. Unlike the Multitask construction, the deviation subnetworks are not shared across groups. The Gaussian hierarchy shrinks their mixing weights toward zero and thereby pools each group toward the global components. Because a deviation subnetwork and its mixing weight enter only through their product, the pooling standard deviations are defined relative to the scale of the deviation subnetworks, which the weight-decay penalty fixes in fitting. As the group-level variances approach zero, Equation (S15) reduces to a common Mixture DLNAM. Nested or crossed groupings follow by adding the corresponding deviation terms and pooling distributions. Full uncertainty propagation for the group-specific deviations and pooling parameters would require extending the last-layer construction of Section S2, which conditions on fitted mixing weights and learned representations. S8 Surface-Target Comparisons The main text emphasizes the cumulative curve as the primary simulation estimand and applied summary. Here we repeat the model comparison and architecture ablation for the full exposure–lag surface, scoring every exposure-grid point at every lag as defined in Section 4.6. Figure S4 and Tables S10, S11, and S12 give the estimator comparison; Fig. S5 and Table S14 give the ablation. As discussed in Section 4.6, surface scoring provides a complementary, finer-grained assessment of how the fitted effect is allocated across exposure and lag; it is not treated as a superior recovery target to the cumulative summary. Figure S4: Simulation Study: Model Comparison (Surface). Monte Carlo simulation of four DGPs, with R=200R=200 replicates per DGP. (A) Log-relative-risk RMSE by method and exposure region, partitioned into squared bias and variance. (B) Empirical coverage of pointwise 95%95\% intervals. Error bars denote 95% Monte Carlo confidence intervals. Boundary denotes exposure values below the 5th or above the 95th percentile of the reference exposure distribution; interior denotes its complement. All lags are included within each exposure region. Associated numerical results are reported in Tables S10 and S11. Figure S5: Simulation Study: Architecture Ablation (Surface). Monte Carlo simulation of four DGPs, with R=50R=50 replicates per DGP. (A) Log-relative-risk RMSE by method and exposure region, partitioned into squared bias and variance. (B) Empirical coverage of pointwise 95%95\% intervals. Error bars denote 95% Monte Carlo confidence intervals. Boundary denotes exposure values below the 5th or above the 95th percentile of the reference exposure distribution; interior denotes its complement. All lags are included within each exposure region. Ablations are applied one at a time and are not factorial. Associated numerical results are reported in Table S14. Figure S6: Simulation Study: ExU Comparison (Surface). Monte Carlo simulation of four DGPs, with R=50R=50 replicates per DGP. (A) Log-relative-risk RMSE by method and exposure region, partitioned into squared bias and variance. (B) Empirical coverage of pointwise 95%95\% intervals. Error bars denote 95% Monte Carlo confidence intervals. Boundary denotes exposure values below the 5th or above the 95th percentile of the reference exposure distribution; interior denotes its complement. All lags are included within each exposure region. Associated RMSE and coverage results are reported in Table S13. S9 Analysis Specifications This section collects supplementary implementation details and the study-specific specifications used in the principal simulation and applications. Tables S2 and S3 summarize the study design, DLNAM fitting, comparator, adjustment, and computing specifications; settings shared across analyses are listed once, with study-specific departures shown separately. For the criterion-selected cross-basis DLNMs, a candidate with k cross-basis coefficients and fitted dispersion φ is scored by QAIC=−2ℓ(μ^)+2φ^k,QBIC=−2ℓ(μ^)+log(n)φ^k.QAIC=-2 ( μ)+2 k, =-2 ( μ)+ (n)\, k. The two criteria differ only in the penalty applied to cross-basis complexity. The nonlinear T-DLNM implementation requires a two-dimensional fixed-effect design, so a single centered linear index was included as a fixed-effect regressor; it forms no part of the DGP. Fits that failed to initialize were re-attempted with a new seed, up to ten times. Of 800800 treed fits in the principal model comparison, 800800 completed, 8383 required more than one attempt (at most 77), and 00 were excluded. Setting Simulation Chicago Malaria A. Analysis Design Family / Link Poisson; Log Quasi-Poisson; Log Bernoulli; Logit Dispersion Fixed: 1 Pearson: φ^=1.03 =1.03 Fixed: 1 Observations n=5000n=5000 n=4638n=4638 n=348,565n=348,565 Exposures / Lags 1; 0–14 1; 0–30 5; 1–6 Reference x⋆=20x =20 Median Median Setting Applicability Specification B. DLNAM Architecture and Fitting Ensemble Members, Subnetworks All B=3B=3; S=3S=3 Surface Component All ExU Input: 128; Hidden Layer: 128; Mish; Linear Output Exposure Input All Min–Max Scaling; Concatenation ExU; (μw,σw)=(1.5,0.5)( _w, _w)=(1.5,0.5) Lag Input All Min–Max Scaling; Concatenation ExU; (μw,σw)=(2.5,0.5)( _w, _w)=(2.5,0.5) Trend Network Chicago ExU; Widths: 128–128–128; Mish; (μw,σw)=(4.5,0.5)( _w, _w)=(4.5,0.5) Covariate Networks Chicago Width: 32; Mish; Linear Input; Z-Score Malaria Month + Year; Width: 32; Mish; Linear Input; Z-Score Grouping Encoding Malaria Cluster + Country; Learned Scalar per Level; Indexed Lookup Last-Layer Parameters All Global Intercept + Target Final Layer Mixing Weights All σω=0.1 _ω=0.1; Unconstrained Optimizer Simulation, Chicago, Malaria AdamW; η=8⋅10−4η=8· 10^-4; ηmin=1⋅10−4 _ =1· 10^-4; Gradient Clip =10=10 Joint Simulation AdamW; η=7⋅10−4η=7· 10^-4; ηmin=1⋅10−4 _ =1· 10^-4; Gradient Clip =10=10 Scheduler All Cosine Annealing; Step: Per Update Epochs, Batching Simulation E=2500E=2500; Full Batch Chicago E=5000E=5000; Full Batch Joint Simulation E=5000E=5000; Full Batch Malaria E=50E=50; Minibatch 1%1\% Regularization All Weight Decay: λw=1⋅10−4 _w=1· 10^-4; Dropout: 0; Subnetwork Dropout: 0; Output Penalty: 0 Simulation, Malaria Surface Roughness Penalty: 0 Chicago Surface Roughness Penalty: 0.003 Uncertainty All Last-Layer Laplace + Between-Member Variation Table S2: Study Design and DLNAM Specifications. Analysis design and DLNAM fitting specifications for the principal simulation and applications; joint-exposure departures are shown where applicable. Setting Applicability Specification A. DLNM Exposure Margin Simulation, Chicago Natural Cubic Spline; Equally Spaced Knots Malaria Natural Cubic Spline; 33 df; Default Knots Lag Margin Simulation, Chicago Natural Cubic Spline; Log-Spaced Knots Malaria Natural Cubic Spline; 33 df Dimension Grid Simulation (vx,vℓ)∈2,…,102(v_x,v_ )∈\2,…,10\^2 Chicago (vx,vℓ)∈2,…,52(v_x,v_ )∈\2,…,5\^2 Selection Criterion Simulation, Chicago QAIC; QBIC Dispersion / Intervals Simulation, Chicago Pearson Dispersion; Pointwise Wald Intervals B. P-DLNM Marginal Bases Simulation P-Spline; 10×1010× 10 Chicago P-Spline; 5×55× 5 Smoothing Simulation, Chicago Second-Order Difference Penalties; REML Intervals Simulation, Chicago Pointwise Wald Intervals C. T-DLNM Likelihood Simulation, Chicago Gaussian; Response: log(1+Yt) (1+Y_t) Trees, Exposure Splits Simulation, Chicago 2020 Trees; Package Priors; 3030 Exposure Splits Sampler Simulation, Chicago Burn-In: 5,0005,000; Post-Burn: 15,00015,000; Thinning: 1010; Retry Cap: 1010 Intervals Simulation, Chicago Posterior Credible Intervals D. Adjustment Trend, Seasonality Chicago Date Spline; 9898 df Malaria Month Spline; 4 df Long-Term Trend Malaria Survey Year; Linear Other Covariates Chicago Day of Week: Categorical (Monday Reference); Dew Point: Spline, 3 df; Ozone: Linear; PM10PM_10: Linear; Exposure Means: Lags 0–1 Target / Controls Malaria Temperature: Precipitation; Precipitation: Temperature; Soil Moisture: Temperature + Precipitation; Evapotranspiration: Soil Moisture + Temperature + Precipitation; Specific Humidity: Soil Moisture + Temperature + Precipitation + Evapotranspiration Multi-Exposure Representation Malaria DLNM: Target-Specific Fit; Controls: Scalar Lag Means; DLNAM: Target-Specific Fit; Controls: Full Lag Surfaces Hierarchy Malaria DLNM: Random Intercepts, Estimated Variance (Cluster, Country); DLNAM: Level-Specific Effects, Fixed Ridge Penalty E. Software and Computing Software Simulation Python 3.10.12; PyTorch 2.5.0a0+e000cf0ad9.nv24.10; NumPy 1.24.4; pandas 2.2.2; R 4.6.1 Applications, Scaling Python 3.11.9; PyTorch 2.13.0+cu126; NumPy 2.4.6; pandas 3.0.3; R 4.5.0 All dlnm 2.4.10; mgcv 1.9.4; dlmtree 1.1.1; glmmTMB 1.1.14 DLNAM Device All CUDA Main MC Compute Simulation Alvis, Chalmers University of Technology Workstation Applications, Scaling Intel Core i7-9700K (8 Logical Cores); 16 GB RAM; NVIDIA GeForce GTX 1070 (8 GB) Scaling Benchmark Computational Scaling n=5000n=5000; M=1,…,4M=1,…,4; R=3R=3; 25002500 Epochs Timing Definition Computational Scaling Construction/Preprocessing + Fitting + Uncertainty + Cumulative Extraction; Excludes Data Generation + Startup; T-DLNM: Sum of Target-Specific Fits Table S3: Comparator and Analysis Specifications. Comparator, adjustment, software, and computing specifications for the principal simulation, applications, and computational-scaling benchmark. RMSE (MCSE) Coverage (MCSE) DGP Model Total Interior Boundary Total Interior Boundary DGP 1 DLNAM 0.0167 (0.0005) 0.0130 (0.0004) 0.0232 (0.0010) 0.850 (0.012) 0.815 (0.014) 0.931 (0.009) DLNM (QAIC) 0.0225 (0.0007) 0.0147 (0.0005) 0.0346 (0.0015) 0.795 (0.012) 0.756 (0.016) 0.888 (0.011) DLNM (QBIC) 0.0238 (0.0005) 0.0191 (0.0005) 0.0323 (0.0010) 0.633 (0.012) 0.563 (0.016) 0.798 (0.012) P-DLNM 0.0256 (0.0005) 0.0220 (0.0004) 0.0326 (0.0011) 0.654 (0.009) 0.545 (0.012) 0.909 (0.009) T-DLNM 0.0155 (0.0002) 0.0153 (0.0002) 0.0161 (0.0005) 0.752 (0.007) 0.661 (0.009) 0.966 (0.007) DGP 2 DLNAM 0.0172 (0.0004) 0.0120 (0.0003) 0.0256 (0.0009) 0.922 (0.007) 0.913 (0.008) 0.943 (0.008) DLNM (QAIC) 0.0287 (0.0007) 0.0192 (0.0005) 0.0435 (0.0013) 0.801 (0.011) 0.777 (0.013) 0.858 (0.010) DLNM (QBIC) 0.0370 (0.0008) 0.0269 (0.0005) 0.0537 (0.0017) 0.560 (0.012) 0.510 (0.014) 0.676 (0.013) P-DLNM 0.0440 (0.0008) 0.0364 (0.0005) 0.0579 (0.0020) 0.490 (0.008) 0.402 (0.010) 0.696 (0.011) T-DLNM 0.0345 (0.0006) 0.0250 (0.0003) 0.0502 (0.0013) 0.646 (0.008) 0.615 (0.010) 0.719 (0.009) DGP 3 DLNAM 0.0240 (0.0014) 0.0114 (0.0003) 0.0404 (0.0028) 0.935 (0.006) 0.935 (0.008) 0.938 (0.008) DLNM (QAIC) 0.0514 (0.0016) 0.0293 (0.0008) 0.0826 (0.0032) 0.642 (0.013) 0.635 (0.016) 0.656 (0.011) DLNM (QBIC) 0.0538 (0.0024) 0.0306 (0.0010) 0.0867 (0.0048) 0.624 (0.014) 0.614 (0.017) 0.646 (0.012) P-DLNM 0.0981 (0.0016) 0.0640 (0.0009) 0.1504 (0.0036) 0.294 (0.009) 0.239 (0.010) 0.422 (0.009) T-DLNM 0.0504 (0.0015) 0.0337 (0.0006) 0.0764 (0.0030) 0.476 (0.009) 0.445 (0.010) 0.549 (0.010) DGP 4 DLNAM 0.0223 (0.0005) 0.0157 (0.0004) 0.0329 (0.0010) 0.895 (0.007) 0.894 (0.008) 0.898 (0.010) DLNM (QAIC) 0.0363 (0.0006) 0.0301 (0.0007) 0.0478 (0.0011) 0.533 (0.008) 0.468 (0.008) 0.685 (0.013) DLNM (QBIC) 0.0351 (0.0005) 0.0285 (0.0004) 0.0471 (0.0010) 0.531 (0.007) 0.463 (0.008) 0.692 (0.011) P-DLNM 0.0399 (0.0012) 0.0369 (0.0014) 0.0463 (0.0012) 0.485 (0.009) 0.391 (0.010) 0.704 (0.012) T-DLNM 0.0359 (0.0006) 0.0290 (0.0004) 0.0483 (0.0013) 0.645 (0.008) 0.581 (0.009) 0.795 (0.012) Table S4: Simulation Study: Model Comparison. Cumulative log-relative-risk RMSE and pointwise 95% interval coverage by DGP, method, and exposure region, with Monte Carlo standard errors in parentheses; R=200R=200. 103×Bias210^3×Bias^2 (MCSE) 103×Variance10^3×Variance (MCSE) DGP Model Total Interior Boundary Total Interior Boundary DGP 1 DLNAM 0.0769 (0.0082) 0.0837 (0.0081) 0.0609 (0.0163) 0.2023 (0.0134) 0.0845 (0.0055) 0.4790 (0.0421) DLNM (QAIC) 0.0637 (0.0100) 0.0658 (0.0084) 0.0588 (0.0230) 0.4444 (0.0321) 0.1497 (0.0091) 1.1368 (0.1008) DLNM (QBIC) 0.1657 (0.0156) 0.1836 (0.0174) 0.1238 (0.0277) 0.4012 (0.0203) 0.1802 (0.0087) 0.9204 (0.0592) P-DLNM 0.2916 (0.0173) 0.3271 (0.0186) 0.2082 (0.0347) 0.3648 (0.0196) 0.1551 (0.0100) 0.8576 (0.0585) T-DLNM 0.1403 (0.0043) 0.1730 (0.0051) 0.0634 (0.0105) 0.1008 (0.0063) 0.0607 (0.0045) 0.1950 (0.0162) DGP 2 DLNAM 0.0291 (0.0040) 0.0365 (0.0042) 0.0117 (0.0060) 0.2664 (0.0149) 0.1063 (0.0055) 0.6427 (0.0447) DLNM (QAIC) 0.1467 (0.0143) 0.1261 (0.0095) 0.1951 (0.0360) 0.6766 (0.0400) 0.2407 (0.0175) 1.7009 (0.1188) DLNM (QBIC) 0.6632 (0.0485) 0.4472 (0.0249) 1.1707 (0.1331) 0.7073 (0.0375) 0.2774 (0.0140) 1.7178 (0.1199) P-DLNM 1.2549 (0.0489) 1.0649 (0.0374) 1.7013 (0.1476) 0.6778 (0.0381) 0.2632 (0.0168) 1.6521 (0.1172) T-DLNM 0.9377 (0.0302) 0.5337 (0.0109) 1.8869 (0.0996) 0.2537 (0.0174) 0.0933 (0.0063) 0.6306 (0.0532) DGP 3 DLNAM 0.1029 (0.0208) 0.0097 (0.0014) 0.3219 (0.0698) 0.4742 (0.0554) 0.1198 (0.0067) 1.3073 (0.1812) DLNM (QAIC) 1.1575 (0.0888) 0.4120 (0.0302) 2.9095 (0.2583) 1.4815 (0.1202) 0.4461 (0.0284) 3.9148 (0.3881) DLNM (QBIC) 1.2833 (0.1040) 0.4436 (0.0366) 3.2567 (0.3071) 1.6143 (0.1927) 0.4917 (0.0338) 4.2523 (0.6377) P-DLNM 6.5289 (0.2732) 3.3630 (0.1293) 13.9687 (0.8815) 3.0937 (0.2279) 0.7324 (0.0465) 8.6428 (0.7251) T-DLNM 1.8178 (0.0928) 0.9466 (0.0351) 3.8650 (0.2569) 0.7208 (0.0837) 0.1916 (0.0164) 1.9644 (0.2673) DGP 4 DLNAM 0.1859 (0.0164) 0.0850 (0.0073) 0.4228 (0.0485) 0.3092 (0.0170) 0.1600 (0.0109) 0.6599 (0.0430) DLNM (QAIC) 0.5428 (0.0227) 0.4653 (0.0215) 0.7251 (0.0408) 0.7744 (0.0374) 0.4406 (0.0285) 1.5588 (0.0972) DLNM (QBIC) 0.4869 (0.0141) 0.4079 (0.0072) 0.6726 (0.0401) 0.7443 (0.0337) 0.4031 (0.0236) 1.5463 (0.0914) P-DLNM 0.5275 (0.0298) 0.5991 (0.0439) 0.3592 (0.0312) 1.0675 (0.0776) 0.7618 (0.0676) 1.7860 (0.1207) T-DLNM 0.8138 (0.0327) 0.6520 (0.0182) 1.1938 (0.0949) 0.4721 (0.0244) 0.1889 (0.0148) 1.1378 (0.0698) Table S5: Simulation Study: Bias and Variance. Squared bias and variance of cumulative log-relative-risk by DGP, method, and exposure region, with Monte Carlo standard errors in parentheses; R=200R=200. Relative RMSE Mean Interval Width DGP Model Total Interior Boundary Total Interior Boundary DGP 1 DLNAM 1.00 0.78 1.39 0.046 0.034 0.076 DLNM (QAIC) 1.35 0.88 2.07 0.053 0.034 0.098 DLNM (QBIC) 1.42 1.14 1.93 0.043 0.027 0.080 P-DLNM 1.53 1.31 1.95 0.054 0.035 0.100 T-DLNM 0.93 0.91 0.96 0.042 0.030 0.069 DGP 2 DLNAM 1.00 0.70 1.49 0.053 0.038 0.088 DLNM (QAIC) 1.67 1.11 2.53 0.068 0.046 0.120 DLNM (QBIC) 2.15 1.57 3.13 0.051 0.033 0.093 P-DLNM 2.56 2.12 3.37 0.059 0.038 0.107 T-DLNM 2.01 1.46 2.92 0.054 0.036 0.094 DGP 3 DLNAM 1.00 0.47 1.68 0.059 0.039 0.105 DLNM (QAIC) 2.14 1.22 3.44 0.071 0.047 0.129 DLNM (QBIC) 2.24 1.27 3.61 0.071 0.047 0.128 P-DLNM 4.08 2.66 6.26 0.073 0.042 0.145 T-DLNM 2.10 1.40 3.18 0.056 0.038 0.099 DGP 4 DLNAM 1.00 0.70 1.48 0.062 0.048 0.095 DLNM (QAIC) 1.63 1.35 2.15 0.062 0.042 0.107 DLNM (QBIC) 1.58 1.28 2.12 0.060 0.041 0.104 P-DLNM 1.79 1.66 2.08 0.058 0.038 0.105 T-DLNM 1.61 1.30 2.17 0.066 0.045 0.114 Table S6: Simulation Study: Relative Error and Interval Width. Relative RMSE and mean pointwise 95% interval width by DGP, method, and exposure region; R=200R=200. Relative RMSE uses the DLNAM total RMSE within each DGP as the common denominator. Interval width is on the log-relative-risk scale. Coverage – Reported (MCSE) Coverage – Laplace (MCSE) DGP Model Total Interior Boundary Total Interior Boundary DGP 1 DLNAM 0.850 (0.012) 0.815 (0.014) 0.931 (0.009) 0.816 (0.012) 0.774 (0.015) 0.917 (0.010) DLNM (QAIC) 0.795 (0.012) 0.756 (0.016) 0.888 (0.011) – – – DLNM (QBIC) 0.633 (0.012) 0.563 (0.016) 0.798 (0.012) – – – P-DLNM 0.654 (0.009) 0.545 (0.012) 0.909 (0.009) – – – T-DLNM 0.752 (0.007) 0.661 (0.009) 0.966 (0.007) – – – DGP 2 DLNAM 0.922 (0.007) 0.913 (0.008) 0.943 (0.008) 0.892 (0.008) 0.877 (0.010) 0.925 (0.009) DLNM (QAIC) 0.801 (0.011) 0.777 (0.013) 0.858 (0.010) – – – DLNM (QBIC) 0.560 (0.012) 0.510 (0.014) 0.676 (0.013) – – – P-DLNM 0.490 (0.008) 0.402 (0.010) 0.696 (0.011) – – – T-DLNM 0.646 (0.008) 0.615 (0.010) 0.719 (0.009) – – – DGP 3 DLNAM 0.935 (0.006) 0.935 (0.008) 0.938 (0.008) 0.903 (0.008) 0.899 (0.009) 0.913 (0.009) DLNM (QAIC) 0.642 (0.013) 0.635 (0.016) 0.656 (0.011) – – – DLNM (QBIC) 0.624 (0.014) 0.614 (0.017) 0.646 (0.012) – – – P-DLNM 0.294 (0.009) 0.239 (0.010) 0.422 (0.009) – – – T-DLNM 0.476 (0.009) 0.445 (0.010) 0.549 (0.010) – – – DGP 4 DLNAM 0.895 (0.007) 0.894 (0.008) 0.898 (0.010) 0.806 (0.010) 0.792 (0.011) 0.837 (0.012) DLNM (QAIC) 0.533 (0.008) 0.468 (0.008) 0.685 (0.013) – – – DLNM (QBIC) 0.531 (0.007) 0.463 (0.008) 0.692 (0.011) – – – P-DLNM 0.485 (0.009) 0.391 (0.010) 0.704 (0.012) – – – T-DLNM 0.645 (0.008) 0.581 (0.009) 0.795 (0.012) – – – Table S7: Simulation Study: DLNAM Interval Decomposition. Empirical pointwise 95% interval coverage for the reported interval and, for DLNAM, the conditional last-layer Laplace interval, by DGP, method, and exposure region, with Monte Carlo standard errors in parentheses; R=200R=200. RMSE (MCSE) Coverage (MCSE) DGP Encoder Total Interior Boundary Total Interior Boundary DGP 1 Concatenation 0.0169 (0.0013) 0.0126 (0.0008) 0.0243 (0.0029) 0.863 (0.023) 0.834 (0.029) 0.931 (0.020) Unified Shared Bias 0.0185 (0.0011) 0.0158 (0.0010) 0.0237 (0.0017) 0.877 (0.022) 0.861 (0.026) 0.912 (0.021) Unified Local Bias 0.0172 (0.0008) 0.0150 (0.0008) 0.0215 (0.0014) 0.760 (0.028) 0.713 (0.033) 0.871 (0.024) DGP 2 Concatenation 0.0172 (0.0009) 0.0111 (0.0005) 0.0264 (0.0018) 0.940 (0.009) 0.935 (0.012) 0.951 (0.013) Unified Shared Bias 0.0325 (0.0018) 0.0257 (0.0012) 0.0445 (0.0031) 0.778 (0.024) 0.767 (0.025) 0.806 (0.030) Unified Local Bias 0.0214 (0.0012) 0.0175 (0.0010) 0.0286 (0.0021) 0.831 (0.021) 0.808 (0.026) 0.885 (0.023) DGP 3 Concatenation 0.0237 (0.0032) 0.0113 (0.0007) 0.0397 (0.0061) 0.937 (0.013) 0.939 (0.016) 0.934 (0.019) Unified Shared Bias 0.0442 (0.0015) 0.0329 (0.0015) 0.0632 (0.0024) 0.788 (0.021) 0.787 (0.026) 0.790 (0.019) Unified Local Bias 0.0346 (0.0024) 0.0167 (0.0008) 0.0580 (0.0046) 0.813 (0.019) 0.845 (0.021) 0.738 (0.030) DGP 4 Concatenation 0.0215 (0.0011) 0.0158 (0.0012) 0.0311 (0.0017) 0.903 (0.013) 0.907 (0.015) 0.892 (0.020) Unified Shared Bias 0.0651 (0.0022) 0.0620 (0.0019) 0.0720 (0.0033) 0.511 (0.030) 0.515 (0.030) 0.502 (0.039) Unified Local Bias 0.0392 (0.0024) 0.0387 (0.0024) 0.0405 (0.0029) 0.749 (0.032) 0.726 (0.033) 0.804 (0.035) Table S8: Simulation Study: ExU Comparison. Cumulative log-relative-risk RMSE and pointwise 95% interval coverage by DGP, method, and exposure region, with Monte Carlo standard errors in parentheses; R=50R=50. RMSE (MCSE) Coverage (MCSE) DGP Configuration Total Interior Boundary Total Interior Boundary DGP 1 Reference 0.0169 (0.0013) 0.0126 (0.0008) 0.0243 (0.0029) 0.863 (0.023) 0.834 (0.029) 0.931 (0.020) ExU Ablation 0.0570 (0.0034) 0.0584 (0.0045) 0.0535 (0.0027) 0.729 (0.028) 0.697 (0.030) 0.805 (0.032) Mixture Ablation 0.0252 (0.0011) 0.0181 (0.0006) 0.0368 (0.0020) 0.779 (0.021) 0.738 (0.021) 0.875 (0.028) Smoothness Ablation 0.0502 (0.0028) 0.0231 (0.0011) 0.0848 (0.0055) 0.944 (0.008) 0.967 (0.009) 0.891 (0.013) DGP 2 Reference 0.0172 (0.0009) 0.0111 (0.0005) 0.0264 (0.0018) 0.940 (0.009) 0.935 (0.012) 0.951 (0.013) ExU Ablation 0.1402 (0.0021) 0.0726 (0.0008) 0.2313 (0.0047) 0.114 (0.018) 0.129 (0.021) 0.080 (0.016) Mixture Ablation 0.0342 (0.0009) 0.0280 (0.0005) 0.0454 (0.0018) 0.595 (0.019) 0.543 (0.017) 0.717 (0.031) Smoothness Ablation 0.0429 (0.0022) 0.0228 (0.0008) 0.0703 (0.0042) 0.963 (0.005) 0.981 (0.004) 0.922 (0.012) DGP 3 Reference 0.0237 (0.0032) 0.0113 (0.0007) 0.0397 (0.0061) 0.937 (0.013) 0.939 (0.016) 0.934 (0.019) ExU Ablation 0.1869 (0.0051) 0.1031 (0.0012) 0.3033 (0.0104) 0.468 (0.022) 0.496 (0.025) 0.403 (0.024) Mixture Ablation 0.0911 (0.0042) 0.0453 (0.0009) 0.1515 (0.0085) 0.512 (0.020) 0.508 (0.021) 0.522 (0.032) Smoothness Ablation 0.0458 (0.0023) 0.0233 (0.0008) 0.0759 (0.0046) 0.966 (0.004) 0.981 (0.003) 0.930 (0.010) DGP 4 Reference 0.0215 (0.0011) 0.0158 (0.0012) 0.0311 (0.0017) 0.903 (0.013) 0.907 (0.015) 0.892 (0.020) ExU Ablation 0.1540 (0.0046) 0.0954 (0.0021) 0.2409 (0.0095) 0.414 (0.031) 0.525 (0.041) 0.152 (0.019) Mixture Ablation 0.0918 (0.0014) 0.0856 (0.0013) 0.1050 (0.0023) 0.265 (0.012) 0.264 (0.011) 0.268 (0.020) Smoothness Ablation 0.0364 (0.0014) 0.0234 (0.0014) 0.0561 (0.0026) 0.965 (0.007) 0.974 (0.009) 0.943 (0.011) Table S9: Simulation Study: Architecture Ablation. Cumulative log-relative-risk RMSE and pointwise 95% interval coverage by DGP, method, and exposure region, with Monte Carlo standard errors in parentheses; R=50R=50. RMSE (MCSE) Coverage (MCSE) DGP Model Total Interior Boundary Total Interior Boundary DGP 1 DLNAM 0.0060 (0.0001) 0.0052 (0.0001) 0.0076 (0.0001) 0.903 (0.004) 0.900 (0.004) 0.912 (0.004) DLNM (QAIC) 0.0053 (0.0001) 0.0042 (0.0001) 0.0072 (0.0001) 0.784 (0.007) 0.776 (0.007) 0.802 (0.007) DLNM (QBIC) 0.0050 (0.0000) 0.0042 (0.0000) 0.0065 (0.0001) 0.686 (0.005) 0.670 (0.006) 0.722 (0.005) P-DLNM 0.0051 (0.0000) 0.0044 (0.0000) 0.0064 (0.0001) 0.837 (0.004) 0.818 (0.004) 0.880 (0.005) T-DLNM 0.0081 (0.0001) 0.0072 (0.0001) 0.0099 (0.0002) 0.928 (0.005) 0.921 (0.005) 0.945 (0.005) DGP 2 DLNAM 0.0076 (0.0001) 0.0057 (0.0001) 0.0108 (0.0001) 0.892 (0.003) 0.908 (0.004) 0.857 (0.004) DLNM (QAIC) 0.0129 (0.0001) 0.0097 (0.0001) 0.0184 (0.0002) 0.790 (0.004) 0.811 (0.004) 0.740 (0.004) DLNM (QBIC) 0.0173 (0.0001) 0.0129 (0.0001) 0.0248 (0.0002) 0.639 (0.004) 0.656 (0.005) 0.599 (0.005) P-DLNM 0.0099 (0.0001) 0.0080 (0.0001) 0.0135 (0.0001) 0.880 (0.003) 0.882 (0.003) 0.875 (0.004) T-DLNM 0.0111 (0.0001) 0.0085 (0.0001) 0.0157 (0.0002) 0.929 (0.003) 0.927 (0.003) 0.933 (0.003) DGP 3 DLNAM 0.0066 (0.0001) 0.0052 (0.0000) 0.0092 (0.0001) 0.927 (0.003) 0.930 (0.003) 0.920 (0.003) DLNM (QAIC) 0.0105 (0.0002) 0.0074 (0.0002) 0.0155 (0.0004) 0.914 (0.003) 0.916 (0.004) 0.911 (0.004) DLNM (QBIC) 0.0101 (0.0002) 0.0093 (0.0000) 0.0118 (0.0006) 0.770 (0.005) 0.754 (0.006) 0.807 (0.005) P-DLNM 0.0138 (0.0002) 0.0105 (0.0001) 0.0195 (0.0004) 0.830 (0.003) 0.831 (0.003) 0.827 (0.004) T-DLNM 0.0105 (0.0002) 0.0080 (0.0001) 0.0148 (0.0004) 0.912 (0.003) 0.907 (0.003) 0.925 (0.003) DGP 4 DLNAM 0.0137 (0.0001) 0.0111 (0.0001) 0.0186 (0.0001) 0.847 (0.003) 0.863 (0.003) 0.809 (0.005) DLNM (QAIC) 0.0157 (0.0002) 0.0127 (0.0001) 0.0212 (0.0002) 0.730 (0.004) 0.740 (0.005) 0.706 (0.005) DLNM (QBIC) 0.0177 (0.0001) 0.0143 (0.0001) 0.0239 (0.0002) 0.607 (0.004) 0.631 (0.004) 0.549 (0.004) P-DLNM 0.0158 (0.0001) 0.0129 (0.0001) 0.0212 (0.0001) 0.809 (0.004) 0.810 (0.004) 0.807 (0.004) T-DLNM 0.0140 (0.0002) 0.0119 (0.0001) 0.0179 (0.0002) 0.864 (0.005) 0.840 (0.006) 0.921 (0.005) Table S10: Simulation Study: Model Comparison (Surface). Log-relative-risk RMSE and pointwise 95% interval coverage by DGP, method, and exposure region, with Monte Carlo standard errors in parentheses; R=200R=200. 103×Bias210^3×Bias^2 (MCSE) 103×Variance10^3×Variance (MCSE) DGP Model Total Interior Boundary Total Interior Boundary DGP 1 DLNAM 0.0131 (0.0008) 0.0110 (0.0006) 0.0179 (0.0011) 0.0234 (0.0008) 0.0161 (0.0006) 0.0405 (0.0014) DLNM (QAIC) 0.0057 (0.0008) 0.0045 (0.0006) 0.0086 (0.0013) 0.0219 (0.0014) 0.0131 (0.0010) 0.0426 (0.0027) DLNM (QBIC) 0.0155 (0.0001) 0.0123 (0.0001) 0.0230 (0.0002) 0.0095 (0.0003) 0.0053 (0.0002) 0.0194 (0.0007) P-DLNM 0.0156 (0.0005) 0.0127 (0.0004) 0.0224 (0.0008) 0.0105 (0.0003) 0.0070 (0.0003) 0.0190 (0.0007) T-DLNM 0.0124 (0.0014) 0.0107 (0.0011) 0.0164 (0.0021) 0.0536 (0.0024) 0.0418 (0.0019) 0.0812 (0.0037) DGP 2 DLNAM 0.0241 (0.0006) 0.0116 (0.0003) 0.0536 (0.0014) 0.0336 (0.0010) 0.0209 (0.0006) 0.0633 (0.0019) DLNM (QAIC) 0.0553 (0.0026) 0.0320 (0.0015) 0.1103 (0.0055) 0.1118 (0.0047) 0.0618 (0.0029) 0.2293 (0.0104) DLNM (QBIC) 0.2139 (0.0092) 0.1199 (0.0051) 0.4348 (0.0193) 0.0864 (0.0068) 0.0474 (0.0036) 0.1780 (0.0149) P-DLNM 0.0506 (0.0013) 0.0332 (0.0008) 0.0917 (0.0029) 0.0482 (0.0012) 0.0303 (0.0008) 0.0902 (0.0027) T-DLNM 0.0419 (0.0017) 0.0232 (0.0009) 0.0858 (0.0038) 0.0824 (0.0028) 0.0494 (0.0020) 0.1601 (0.0053) DGP 3 DLNAM 0.0138 (0.0004) 0.0114 (0.0002) 0.0195 (0.0010) 0.0301 (0.0008) 0.0155 (0.0005) 0.0646 (0.0021) DLNM (QAIC) 0.0218 (0.0015) 0.0095 (0.0007) 0.0508 (0.0044) 0.0876 (0.0044) 0.0445 (0.0024) 0.1889 (0.0103) DLNM (QBIC) 0.0696 (0.0020) 0.0731 (0.0017) 0.0616 (0.0045) 0.0334 (0.0035) 0.0141 (0.0016) 0.0785 (0.0104) P-DLNM 0.1188 (0.0038) 0.0796 (0.0021) 0.2111 (0.0116) 0.0724 (0.0034) 0.0317 (0.0010) 0.1681 (0.0104) T-DLNM 0.0399 (0.0017) 0.0235 (0.0007) 0.0785 (0.0048) 0.0706 (0.0030) 0.0406 (0.0017) 0.1414 (0.0075) DGP 4 DLNAM 0.1530 (0.0020) 0.0995 (0.0013) 0.2789 (0.0038) 0.0356 (0.0010) 0.0227 (0.0006) 0.0661 (0.0020) DLNM (QAIC) 0.1137 (0.0016) 0.0777 (0.0010) 0.1983 (0.0033) 0.1332 (0.0046) 0.0837 (0.0031) 0.2496 (0.0097) DLNM (QBIC) 0.2171 (0.0067) 0.1452 (0.0042) 0.3861 (0.0126) 0.0971 (0.0065) 0.0595 (0.0039) 0.1854 (0.0134) P-DLNM 0.1832 (0.0032) 0.1196 (0.0020) 0.3327 (0.0061) 0.0670 (0.0018) 0.0467 (0.0014) 0.1147 (0.0032) T-DLNM 0.0623 (0.0045) 0.0513 (0.0031) 0.0883 (0.0081) 0.1325 (0.0045) 0.0903 (0.0032) 0.2318 (0.0081) Table S11: Simulation Study: Bias and Variance (Surface). Squared bias and variance of log-relative-risk by DGP, method, and exposure region, with Monte Carlo standard errors in parentheses; R=200R=200. Relative RMSE Mean Interval Width DGP Model Total Interior Boundary Total Interior Boundary DGP 1 DLNAM 1.00 0.86 1.27 0.017 0.014 0.025 DLNM (QAIC) 0.87 0.69 1.18 0.011 0.008 0.018 DLNM (QBIC) 0.83 0.70 1.08 0.009 0.007 0.014 P-DLNM 0.85 0.73 1.07 0.013 0.010 0.019 T-DLNM 1.34 1.20 1.64 0.023 0.019 0.033 DGP 2 DLNAM 1.00 0.75 1.42 0.020 0.016 0.031 DLNM (QAIC) 1.70 1.27 2.43 0.028 0.021 0.042 DLNM (QBIC) 2.28 1.70 3.26 0.020 0.016 0.031 P-DLNM 1.31 1.05 1.78 0.026 0.021 0.038 T-DLNM 1.47 1.12 2.06 0.033 0.025 0.052 DGP 3 DLNAM 1.00 0.78 1.38 0.019 0.014 0.030 DLNM (QAIC) 1.58 1.11 2.34 0.024 0.018 0.037 DLNM (QBIC) 1.53 1.41 1.79 0.012 0.009 0.020 P-DLNM 2.09 1.59 2.94 0.023 0.018 0.034 T-DLNM 1.59 1.21 2.24 0.021 0.016 0.033 DGP 4 DLNAM 1.00 0.80 1.35 0.022 0.017 0.032 DLNM (QAIC) 1.14 0.93 1.54 0.030 0.024 0.045 DLNM (QBIC) 1.29 1.04 1.74 0.021 0.017 0.032 P-DLNM 1.15 0.94 1.54 0.029 0.024 0.042 T-DLNM 1.02 0.87 1.30 0.035 0.027 0.053 Table S12: Simulation Study: Relative Error and Interval Width (Surface). Relative RMSE and mean pointwise 95% interval width by DGP, method, and exposure region; R=200R=200. Relative RMSE uses the DLNAM total RMSE within each DGP as the common denominator. Interval width is on the log-relative-risk scale. RMSE (MCSE) Coverage (MCSE) DGP Encoder Total Interior Boundary Total Interior Boundary DGP 1 Concatenation 0.0063 (0.0002) 0.0055 (0.0001) 0.0079 (0.0002) 0.888 (0.010) 0.883 (0.011) 0.901 (0.010) Unified Shared Bias 0.0111 (0.0003) 0.0100 (0.0003) 0.0133 (0.0004) 0.949 (0.005) 0.946 (0.006) 0.958 (0.006) Unified Local Bias 0.0079 (0.0002) 0.0070 (0.0002) 0.0098 (0.0003) 0.929 (0.005) 0.920 (0.005) 0.950 (0.006) DGP 2 Concatenation 0.0075 (0.0001) 0.0056 (0.0001) 0.0107 (0.0002) 0.898 (0.006) 0.912 (0.007) 0.864 (0.008) Unified Shared Bias 0.0140 (0.0002) 0.0124 (0.0002) 0.0172 (0.0003) 0.891 (0.007) 0.890 (0.008) 0.893 (0.007) Unified Local Bias 0.0100 (0.0002) 0.0085 (0.0002) 0.0129 (0.0003) 0.900 (0.007) 0.894 (0.009) 0.912 (0.009) DGP 3 Concatenation 0.0066 (0.0001) 0.0051 (0.0001) 0.0091 (0.0003) 0.926 (0.006) 0.927 (0.007) 0.921 (0.007) Unified Shared Bias 0.0157 (0.0004) 0.0129 (0.0003) 0.0209 (0.0005) 0.938 (0.006) 0.938 (0.007) 0.940 (0.007) Unified Local Bias 0.0091 (0.0002) 0.0066 (0.0001) 0.0132 (0.0004) 0.962 (0.004) 0.966 (0.004) 0.952 (0.006) DGP 4 Concatenation 0.0138 (0.0001) 0.0111 (0.0001) 0.0185 (0.0002) 0.849 (0.007) 0.865 (0.007) 0.810 (0.009) Unified Shared Bias 0.0221 (0.0003) 0.0195 (0.0003) 0.0273 (0.0004) 0.838 (0.007) 0.838 (0.008) 0.839 (0.008) Unified Local Bias 0.0159 (0.0004) 0.0137 (0.0003) 0.0203 (0.0005) 0.872 (0.007) 0.872 (0.008) 0.873 (0.009) Table S13: Simulation Study: ExU Comparison (Surface). Log-relative-risk RMSE and pointwise 95% interval coverage by DGP, method, and exposure region, with Monte Carlo standard errors in parentheses; R=50R=50. RMSE (MCSE) Coverage (MCSE) DGP Configuration Total Interior Boundary Total Interior Boundary DGP 1 Reference 0.0063 (0.0002) 0.0055 (0.0001) 0.0079 (0.0002) 0.888 (0.010) 0.883 (0.011) 0.901 (0.010) ExU Ablation 0.0177 (0.0003) 0.0172 (0.0003) 0.0190 (0.0003) 0.511 (0.017) 0.511 (0.020) 0.512 (0.016) Mixture Ablation 0.0095 (0.0002) 0.0085 (0.0001) 0.0115 (0.0002) 0.884 (0.004) 0.872 (0.005) 0.910 (0.004) Smoothness Ablation 0.0110 (0.0002) 0.0077 (0.0002) 0.0163 (0.0004) 0.955 (0.004) 0.967 (0.003) 0.927 (0.007) DGP 2 Reference 0.0075 (0.0001) 0.0056 (0.0001) 0.0107 (0.0002) 0.898 (0.006) 0.912 (0.007) 0.864 (0.008) ExU Ablation 0.0321 (0.0001) 0.0237 (0.0000) 0.0463 (0.0002) 0.206 (0.010) 0.229 (0.010) 0.153 (0.009) Mixture Ablation 0.0109 (0.0002) 0.0084 (0.0001) 0.0151 (0.0002) 0.841 (0.007) 0.852 (0.007) 0.816 (0.009) Smoothness Ablation 0.0108 (0.0002) 0.0079 (0.0002) 0.0155 (0.0003) 0.950 (0.003) 0.963 (0.003) 0.922 (0.005) DGP 3 Reference 0.0066 (0.0001) 0.0051 (0.0001) 0.0091 (0.0003) 0.926 (0.006) 0.927 (0.007) 0.921 (0.007) ExU Ablation 0.0245 (0.0003) 0.0192 (0.0002) 0.0337 (0.0008) 0.586 (0.017) 0.619 (0.019) 0.507 (0.019) Mixture Ablation 0.0130 (0.0003) 0.0091 (0.0001) 0.0193 (0.0006) 0.839 (0.006) 0.830 (0.007) 0.860 (0.006) Smoothness Ablation 0.0113 (0.0002) 0.0076 (0.0001) 0.0170 (0.0004) 0.964 (0.003) 0.976 (0.002) 0.935 (0.005) DGP 4 Reference 0.0138 (0.0001) 0.0111 (0.0001) 0.0185 (0.0002) 0.849 (0.007) 0.865 (0.007) 0.810 (0.009) ExU Ablation 0.0372 (0.0002) 0.0303 (0.0002) 0.0496 (0.0004) 0.365 (0.017) 0.382 (0.019) 0.322 (0.016) Mixture Ablation 0.0200 (0.0002) 0.0173 (0.0002) 0.0252 (0.0002) 0.676 (0.009) 0.661 (0.012) 0.712 (0.006) Smoothness Ablation 0.0130 (0.0002) 0.0102 (0.0002) 0.0179 (0.0003) 0.921 (0.004) 0.935 (0.005) 0.889 (0.006) Table S14: Simulation Study: Architecture Ablation (Surface). Log-relative-risk RMSE and pointwise 95% interval coverage by DGP, method, and exposure region, with Monte Carlo standard errors in parentheses; R=50R=50. RMSE (MCSE) Coverage (MCSE) Component Estimator Total Interior Boundary Total Interior Boundary DGP 1 DLNAM 0.0167 (0.0009) 0.0137 (0.0006) 0.0215 (0.0017) 0.854 (0.020) 0.813 (0.021) 0.938 (0.023) DLNM (QAIC) 0.0264 (0.0015) 0.0182 (0.0009) 0.0383 (0.0028) 0.793 (0.023) 0.748 (0.028) 0.887 (0.024) DLNM (QBIC) 0.0320 (0.0021) 0.0227 (0.0014) 0.0457 (0.0038) 0.639 (0.026) 0.588 (0.030) 0.746 (0.032) P-DLNM 0.0330 (0.0018) 0.0265 (0.0012) 0.0436 (0.0035) 0.617 (0.023) 0.525 (0.025) 0.809 (0.029) T-DLNM 0.0512 (0.0029) 0.0404 (0.0021) 0.0684 (0.0052) 0.607 (0.024) 0.551 (0.024) 0.725 (0.037) DGP 2 DLNAM 0.0173 (0.0007) 0.0143 (0.0005) 0.0223 (0.0014) 0.890 (0.013) 0.858 (0.016) 0.956 (0.014) DLNM (QAIC) 0.0304 (0.0014) 0.0211 (0.0011) 0.0439 (0.0025) 0.778 (0.019) 0.760 (0.025) 0.815 (0.021) DLNM (QBIC) 0.0405 (0.0015) 0.0291 (0.0013) 0.0575 (0.0028) 0.582 (0.024) 0.541 (0.029) 0.668 (0.027) P-DLNM 0.0474 (0.0015) 0.0380 (0.0013) 0.0626 (0.0031) 0.507 (0.021) 0.431 (0.023) 0.666 (0.028) T-DLNM 0.0596 (0.0029) 0.0403 (0.0014) 0.0871 (0.0059) 0.629 (0.020) 0.600 (0.023) 0.690 (0.033) DGP 3 DLNAM 0.0241 (0.0014) 0.0136 (0.0005) 0.0382 (0.0029) 0.888 (0.015) 0.891 (0.018) 0.880 (0.020) DLNM (QAIC) 0.0459 (0.0024) 0.0288 (0.0021) 0.0700 (0.0039) 0.646 (0.025) 0.644 (0.031) 0.650 (0.024) DLNM (QBIC) 0.0471 (0.0025) 0.0293 (0.0021) 0.0721 (0.0042) 0.639 (0.025) 0.635 (0.029) 0.648 (0.025) P-DLNM 0.0971 (0.0024) 0.0693 (0.0019) 0.1398 (0.0054) 0.273 (0.015) 0.213 (0.016) 0.403 (0.020) T-DLNM 0.0716 (0.0027) 0.0502 (0.0013) 0.1040 (0.0056) 0.483 (0.020) 0.434 (0.022) 0.591 (0.027) DGP 4 DLNAM 0.0224 (0.0008) 0.0170 (0.0008) 0.0308 (0.0016) 0.847 (0.017) 0.824 (0.021) 0.894 (0.015) DLNM (QAIC) 0.0368 (0.0013) 0.0297 (0.0008) 0.0483 (0.0024) 0.506 (0.018) 0.439 (0.016) 0.645 (0.028) DLNM (QBIC) 0.0377 (0.0013) 0.0301 (0.0010) 0.0500 (0.0022) 0.494 (0.018) 0.430 (0.016) 0.626 (0.029) P-DLNM 0.0463 (0.0019) 0.0394 (0.0017) 0.0582 (0.0029) 0.440 (0.017) 0.361 (0.016) 0.606 (0.027) T-DLNM 0.0722 (0.0026) 0.0643 (0.0027) 0.0864 (0.0040) 0.428 (0.020) 0.365 (0.023) 0.560 (0.027) Null Exposure DLNAM 0.0124 (0.0008) 0.0092 (0.0006) 0.0177 (0.0014) 0.955 (0.012) 0.964 (0.013) 0.935 (0.022) DLNM (QAIC) 0.0167 (0.0016) 0.0114 (0.0011) 0.0249 (0.0032) 0.897 (0.023) 0.879 (0.031) 0.937 (0.019) DLNM (QBIC) 0.0123 (0.0008) 0.0079 (0.0006) 0.0188 (0.0012) 0.928 (0.029) 0.930 (0.031) 0.924 (0.028) P-DLNM 0.0171 (0.0013) 0.0119 (0.0010) 0.0253 (0.0019) 0.841 (0.042) 0.850 (0.043) 0.821 (0.042) T-DLNM 0.0603 (0.0044) 0.0459 (0.0024) 0.0845 (0.0089) 0.549 (0.035) 0.507 (0.036) 0.644 (0.045) Table S15: Simulation Study: Model Comparison (Joint). Cumulative log-relative-risk RMSE and pointwise 95% interval coverage by component, method, and exposure region, with Monte Carlo standard errors in parentheses; R=50R=50.