arXiv CS open access

Implementing neural network mixed-effects models in Template Model Builder (TMB)

2026 · arxiv_cs

arXiv CS · Papers · 2026

Conceptio Open Archive

Analyze & explore related research in Conceptio

Connect this record directly into your research workflow or AI agent. Explore citations across open repositories, query full text, and export structured proof bundles.

Explore in Search →Connect via MCP / API ↗

Implementing neural network mixed-effects models in Template Model Builder (TMB)

arXiv:2608.31133v1 [stat.ML] 31 Aug 2026

Nan Zheng∗1 , Hoi Yiu Cheung1 , Vibhu Sharma1 , James T. Thorson2 , and Noel G. Cadigan3 1

Department of Mathematics and Statistics, Memorial University of Newfoundland, Canada 2 Resource Ecology and Fisheries Management, Alaska Fisheries Science Center, National Marine Fisheries Service, National Oceanic and Atmospheric Administration, USA 3 Centre for Fisheries Ecosystems Research, Fisheries and Marine Institute of Memorial University of Newfoundland, Canada September 1, 2026

Abstract Neural network mixed-effects models (NMMs) have gained traction by combining the strong representation and predictive power of artificial neural networks with the capacity of mixed-effects modeling to capture complex correlation structures. However, existing estimation approaches rely heavily on manual derivations of objective functions and gradients, which inherently forces simplifying approximations and severely constrains the complexity and accuracy of NMMs. In this work, we introduce a general framework for implementing NMMs using Template Model Builder (TMB). By leveraging automatic differentiation and Laplace approximation, TMB requires users to specify only the negative joint log-likelihood and any regularization terms. The framework automatically integrates out random effects and evaluates the marginal objective function alongside its exact gradients, eliminating the need for manual derivations or ad hoc approximations. We demonstrate the efficiency, flexibility, and statistical performance of TMB-based NMMs across two numerical examples, including an application to monotonic NMMs. Reproducible code is provided to facilitate broader adoption.

Keywords: neural network mixed-effects model; TMB; monotonic neural network; feedforward neural network; Laplace approximation; fish maturation model ∗

Email: k33nz@mun.ca

Introduction

The development of technologies for implementing deep neural networks, that is, neural networks with multiple hidden layers, has significantly advanced the fields of artificial intelligence and machine learning (LeCun et al., 2015). Owing to their strong representational and predictive capabilities, neural networks have been increasingly adopted for analyzing conventional statistical datasets. Compared with modern deep learning applications, these datasets are typically much smaller, making relatively shallow network architectures a practical and commonly adopted choice (Hornik et al., 1989; Hastie et al., 2009; Goodfellow et al., 2016). A standard feedforward neural network does not explicitly encode structural priors, such as spatial, temporal, or contextual dependencies, within its architecture. As a result, it must learn these relationships solely from the training data, leading to a hypothesis space that is often substantially larger than necessary for structured learning tasks (Bengio, 2009; Goodfellow et al., 2016). Subsequent developments have addressed this limitation by incorporating domain-specific correlation structures into network architectures. For example, convolutional neural networks exploit translation invariance and spatial hierarchies in image data, recurrent neural networks (RNNs) capture temporal dependencies in sequential data, and transformers model long-range contextual relationships in natural language (see, e.g., Chollet et al., 2022). Although these architectures substantially improve learning efficiency by embedding appropriate structural assumptions, they typically remain data-intensive and often require training data with specific structural characteristics that are not satisfied in many statistical applications. For example, RNNs are generally designed for regularly sampled sequential data and are most naturally applied to collections of sequences with comparable temporal structure, assumptions that are frequently violated by longitudinal datasets with irregular observation times or unequal sequence lengths (Mandel et al., 2023). In contrast, statistical mixed-effects models provide considerably greater flexibility for ac-

commodating unbalanced longitudinal data, irregular observation schedules, and complex sampling designs while remaining effective with relatively small sample sizes. Motivated by the complementary strengths of these two paradigms, neural network mixed-effects models (NMMs) have emerged as a natural extension that combines the representational power of neural networks with the flexibility and interpretability of statistical mixed-effects models (Tandon et al., 2006; Tran et al., 2017; Mandel et al., 2023). A NMM can be generally formulated as µ = g{fN (X) + Z⊤ b}, where fN (X) represents the output of a neural network with input X, Z is the design matrix associated with the random effects (REs) b, µ denotes the vector of parameters governing the observational model (e.g., the mean vector), and g(·) is an elementwise transformation analogous to the inverse link function in generalized linear mixed-effects models (GLMMs). Throughout this paper, boldface notation is used for vectors and matrices. The REs b are typically excluded from the neural network fN (X), as their inclusion may lead to prohibitive computational costs. Together with the observational model f (D|µ) for responses D, the RE prior f (b), and the weight regularization term freg (w), the joint objective function is f (D, b) = f (D|µ)f (b)freg (w).

(1)

The weight regularization term penalizes excessively large neural network weights, thereby reducing overfitting and improving the generalization ability of the model (see, e.g., Krogh and Hertz, 1991). Before applying optimization algorithms, such as stochastic gradient descent, in deep learning, the REs b must be appropriately accommodated to obtain the objective function and its gradients. To this end, Tran et al. (2017) perform estimation in a Bayesian framework using variational approximation, whereas Mandel et al. (2023) replace the conditional likelihood function f (D|µ) with a quasi-likelihood function and integrate out b using the Laplace approximation to obtain the marginal objective function. In this work, we adopt the latter approach because of its computational simplicity and the availability of efficient modeling packages, such as the Template Model Builder (TMB; 3

Kristensen et al., 2016) within R (R Core Team, 2024), which automatically integrates out the REs from the joint likelihood and provides the marginal likelihood together with its gradients. By leveraging automatic differentiation alongside sparse matrix techniques and parallel computing, TMB enables rapid likelihood-based inference even in high-dimensional settings, and offers a highly competitive alternative to Bayesian methods in terms of speed, scalability, and model complexity (Thygesen et al., 2017). Without leveraging modern computational frameworks such as TMB, Mandel et al. (2023) relied on manual derivations of the marginal objective function and its corresponding gradient to enable numerical optimization. This approach introduces several key limitations. First, as discussed previously, the quasi-likelihood function, rather than the full likelihood, must be used to facilitate manual derivation. Second, the approximation requires omission of a term involving the determinant of the covariance matrix of the REs and the Hessian matrix of the joint objective function with respect to the REs. Neglecting this term may reduce the accuracy and efficiency of inference for certain model parameters. Third, the analytical complexity increases substantially as the number of neural network layers grows, limiting the scalability of the approach (to two layers in Mandel et al., 2023) and making the derivation and implementation increasingly difficult and error-prone. Notably, Mandel et al. (2023) omitted the derivative of the RE estimators with respect to the model parameters when differentiating the objective function (see, e.g., Eq. 7 in Kristensen et al., 2016), an omission that may lead to suboptimal optimization or less precise parameter estimates. To address these limitations, we propose and demonstrate the use of TMB for implementing NMMs. Rather than using a quasi-likelihood approximation, we directly provide the joint objective function (1) to TMB. TMB then automatically and efficiently computes the Laplace approximation to the marginal objective function, together with its corresponding gradients, without omitting any terms. Consequently, the proposed approach improves

computational efficiency, inferential accuracy, ease of implementation, and scalability relative to the manual derivation approach. The remainder of this paper is organized as follows. Section 2 details the proposed NMM framework and its implementation using TMB. Section 3 replicates the simulation studies from Mandel et al. (2023), demonstrating the computational efficiency and superior accuracy of our approach. In Section 4, we apply both the NMM and GLMM to an American plaice maturity dataset, illustrating the ease of implementing monotonic neural networks within this framework. Finally, Section 5 offers concluding remarks and directions for future research.

Methods

Figure 1: The architecture of a multi-layer NMM.

Our model architecture is displayed in Fig. 1 and closely follows the feed-forward network framework of Mandel et al. (2023). Let xi ∈ Rp denote the input vector for the ith observation and µi ∈ R be its corresponding univariate output. For a network with L hidden layers, where layer l ∈ {1, . . . , L} consists of dl units, each layer sequentially applies

an affine transformation followed by an element-wise activation function fa (·) (see, e.g., Chollet et al., 2022). Specifically, for the ith observation, the output of the first hidden (1)

layer, αi ∈ Rd1 , is expressed as (1)

αi = fa {w(1) xi + δ (1) },

(2)

where w(1) ∈ Rd1 ×p and δ (1) ∈ Rd1 represent the weight matrix and bias vector, respectively. For each hidden layer l ∈ {2, . . . , L}, the activation output is given by (l)

(l−1)

αi = fa {w(l) αi

+ δ (l) },

(3)

with parameters w(l) ∈ Rdl ×dl−1 and δ (l) ∈ Rdl . Finally, the univariate output of the neural network is related to the observational parameter µi via (L)

µi (bi ) = g{w αi

+ δ + Z⊤ i bi },

(4)

where w ∈ R1×dL , δ ∈ R, and Zi and bi denote the design matrix and RE vector associated with the ith observation, respectively. The REs bi are assumed to follow a multivariate normal (MVN) distribution, inducing dependence among the components of the ith observation. Dependence across observations can be introduced either by sharing one or more RE components among observations or by specifying a correlated MVN distribution for the REs bi across observations. In Mandel et al. (2023), µi denotes the conditional mean of the response under an exponential-family distribution given bi , which enables the use of a quasi-likelihood approximation to the true conditional likelihood. In contrast, our framework is not restricted to exponential-family distributions. The observational distribution can be specified more generally, and µi need not represent the conditional mean; rather, it may denote any governing parameter of the observational distribution. A comparison can also be made to the GLMM framework (Breslow and Clayton, 1993), in which the conditional mean is linked to the linear predictor via ⊤ µi = g{x⊤ i β + β0 + Zi bi },

(5)

where β represents the vector of fixed effects, µi denotes the conditional mean of the ith observation given the REs, and g(·) is the inverse link function. The NMM defined in Eqs. (2)–(4) replaces the fixed-effects portion of the linear predictor with the output of a feed-forward neural network. Because neural networks offer significantly greater expressiveness than linear predictors, the NMM is expected to yield better representational and predictive capabilities compared to the standard GLMM. ⊤ ⊤ Let there be n observations, and let B = (b⊤ 1 , . . . , bn ) denote the collection of all REs.

The marginal likelihood l for NMM defined by (2)–(4) is expressed as # Z "Y n f {Di |µi (bi )} f (B)dB, l(θ|D1 , . . . , Dn ) =

(6)

i=1

where Di denotes the ith observation, f {Di | µi (bi )} is its conditional distribution given the REs, and f (B) is the prior (marginal) distribution of B. The vector θ collects all model parameters, including the neural network weights and biases, observational model parameters, and prior distribution parameters. To evaluate the integral over the REs in Eq. (6), we employ TMB. TMB efficiently approximates the marginal likelihood via the Laplace approximation, scaling to tens of thousands of REs by exploiting the sparsity in their joint distribution with the observations. The model parameters are estimated by minimizing the objective function − log{l(θ|D1 , . . . , Dn )} + λ(W⊤ W + ∆⊤ ∆),

(7)

where W and ∆ denote the vectors of all neural network weights and biases, respectively, and λ > 0 is a regularization hyperparameter. The second term introduces L2 regularization, which penalizes larger parameter magnitudes to mitigate overfitting and enhance generalization. Following Mandel et al. (2023), we employ an L2 penalty to prioritize predictive accuracy over feature selection (Tibshirani, 1996; Goodfellow et al., 2016) while ensuring differentiability.

Simulation study

We assess the proposed TMB implementation of the NMM by reproducing the nonlinear simulation experiment of Mandel et al. (2023), considering both a continuous and a binary response. In each case the network is provided only with the raw predictors and must recover a nonlinear signal from them, and predictive performance is assessed by leaveone-out cross-validation on the final observation. To facilitate direct comparison with the results of Mandel et al. (2023), we adopt the notation used in their simulation studies throughout this section.

3.1

Simulation design

For each of m = 100 subjects, the number of observations ni is drawn from a Poisson(6) + 2 distribution, ensuring at least one training and one validation observation per subject. For each observation we generate six independent binary predictors xij = (xij1 , . . . , xij6 )⊤ , each from a Bernoulli(0.5) distribution, and form three nonlinear features Zij1 = 1(xij1 + xij2 = 1),

Zij2 = 1(xij3 + xij4 ̸= 1),

Zij3 = 1(xij5 + xij6 = 1),

(8)

which encode exclusive-or (XOR) type interactions that a linear predictor cannot represent. Writing Zij = (Zij1 , Zij2 , Zij3 )⊤ , the continuous response is generated as yij = Z⊤ ij β + bi + ϵij ,

β ∼ MVN(2 · 1, I),

(9)

and the binary response as logit{P(yij = 1)} = Z⊤ ij β + bi + ϵij ,

β ∼ MVN(3 · 1, I),

(10)

where the subject-level random intercept is bi ∼ N (0, τ ) and ϵij ∼ N (0, 0.052 ), and β is drawn once per replicate. The variance component τ governs the strength of withinsubject correlation and is swept over {0, 0.1, 0.2, 0.3, 0.4, 0.5} for the continuous response and {0, 4, 8, 12, 16, 20} for the binary response, following Mandel et al. (2023). Crucially, 8

the model is fit using only the raw predictors xij ; the features Zij are never supplied, so the network must recover the XOR structure of (8) on its own. Predictive accuracy is evaluated by leave-one-out cross-validation: the final observation of each subject is withheld, the model is fit to the remaining observations, and the withheld observation is predicted as µ̂ij = g{fˆN (xij ) + b̂i }, combining the fitted network output with the subject’s estimated random intercept. Performance is summarized by the mean squared prediction error (MSPE) for the continuous response and by the area under the receiver operating characteristic curve (AUC) for the binary response, averaged across the 100 withheld observations. Each configuration is replicated 500 times.

3.2

Model and estimation

We implement the two-hidden-layer network of Mandel et al. (2023), denoted GNMM-2, with d1 = 10 units in the first hidden layer and d2 = 5 in the second, a six-dimensional input, and a scalar output, as in Eqs. (2)–(4). The activation function fa is the logistic sigmoid, and the output transformation g is the identity for the continuous response and the inverse logit for the binary response. A single scalar random intercept enters the output layer for each subject (Zi ≡ 1 in Eq. (4)), with prior bi ∼ N (0, τ ). The observational distribution is Gaussian with an estimated residual variance for the continuous response and Bernoulli for the binary response, as specified in Eqs. (9) and (10), respectively. Estimation follows the objective (7) with regularization parameter λ = 0.005, matching the value used for the two-hidden-layer network in Mandel et al. (2023). The random intercepts are integrated out by TMB via the Laplace approximation, and the resulting marginal objective is optimized using nlminb function in R; standard deviations are parameterized on the log scale to enforce positivity, and λ is held fixed. In contrast to the manual quasi-likelihood derivation of Mandel et al. (2023), the TMB implementation requires only a specification of the joint objective (1), automatically deriving the marginal

objective and its gradients.

3.3

Results

Figures 2 and 3 present the predictive performance of the TMB implementation of GNMM-2 as a function of τ , overlaid on the corresponding GNMM-2 results reported by Mandel et al. (2023). For the continuous response, the TMB MSPE is essentially flat at approximately 0.003 across all values of τ , closely tracking the irreducible noise floor σϵ2 = 0.0025: the fitted network recovers the XOR signal almost exactly and the random intercept absorbs the within-subject correlation, leaving only observation noise. For the binary response, the TMB AUC rises from approximately 0.89 at τ = 0 to 0.91 at τ = 20, increasing with τ because a larger random-intercept variance renders the subject-specific estimate b̂i more informative for the withheld observation. In both settings the TMB implementation attains markedly better predictive performance than the GNMM-2 results reported by Mandel et al. (2023), whose MSPE is approximately 0.5 and whose AUC ranges from approximately 0.78 to 0.88.

Figure 2: Mean squared prediction error (MSPE) of the TMB implementation of GNMM-2 for the continuous response (solid, blue), as a function of the random-intercept variance τ , overlaid on the GNMM-2 results reported by Mandel et al. (2023) (grey dashed, read from their Figure 2). Both panels show the same two curves for the nonlinear data-generating mechanism and a first hidden layer of 10 units. Panel (a) uses the full y-axis range [0, 4] of Mandel et al. (2023): on this scale the TMB curve is pinned to the horizontal axis, clearly far below the reported GNMM-2 curve. Panel (b) shows the same data on a restricted yaxis [0, 0.6], which makes the difference explicit: the TMB MSPE (≈ 0.003) lies just above the axis, well below the reported GNMM-2 MSPE (≈ 0.45–0.55).

Figure 3: Area under the ROC curve (AUC) of the TMB implementation of GNMM-2 for the binary response (solid), as a function of the random-intercept variance τ , overlaid on the GNMM-2 results reported by Mandel et al. (2023) (grey dashed, read from their Figure 3). Both correspond to the nonlinear data-generating mechanism and a first hidden layer of 10 units; the axes match the scale of Mandel et al. (2023).

Real data example: maturity data for 3LNO American plaice

To illustrate the implementation of the NMM using real data, we fitted the model to maturity data for female American plaice (Hippoglossoides platessoides) collected during annual surveys conducted by Fisheries and Oceans Canada (DFO) in the Northwest Atlantic Fisheries Organization (NAFO) Divisions 3L, 3N, and 3O (3LNO) from 1978 to 2015. Fish born in the same year (i.e., belonging to the same cohort) are generally exposed to similar environmental conditions and other shared influences throughout their life history. Therefore, cohort (birth year), rather than survey year, was used as the temporal index. In total, the dataset comprises 58 cohorts spanning birth years from 1958 to 2015. We further partition the 58 cohorts into cohort segments (cohort seg) of eight consecutive cohorts each. Specifically, the first eight cohorts are assigned cohort seg 1, the next eight cohorts are assigned cohort seg 2, and so forth. Under this partitioning, the final cohort segment contains only three cohorts. To avoid such a small segment, we merge it with the penultimate cohort segment (cohort seg 7). Consequently, the data are divided into seven cohort segments, with the first six segments each containing eight consecutive cohorts and the seventh segment containing 11 consecutive cohorts. The NMM consisted of two hidden layers, each containing twenty nodes. Input for the ith fish was an eight-dimensional vector comprising age as the first component, followed by a seven-dimensional one-hot encoded representation of the cohort segment (cohort seg number), where the active category was assigned a value of 1 and all remaining categories were set to 0. The activation function fa in Fig. 1 is the logistic sigmoid. The output yields the conditional maturation probability for the ith fish given the REs, specifying Eq. (4) as (2)

Pi (ai |θ, bci ) = inv logit{w αi + δ + bci },

(11)

where ai and ci denote respectively the age and cohort index for the ith individual, 13

inv logit(·) represents the inverse logit transformation function, and cohort-specific REs (b1 , . . . , b58 ) are modeled using a first-order autoregressive [AR(1)] process across the 58 cohorts. Maturity was recorded as a binary response variable (di ), where 1 indicates a mature fish and 0 indicates an immature fish. Accordingly, the observational model f {Di |µi (bi )} in Eq. (6) takes the Bernoulli form: f {di |Pi } = {Pi }di {1 − Pi }1−di . The objective function is defined by Eqs. (6)-(7). The maturation probability in Eq. (11) is constrained to be a non-decreasing function of age. This monotonicity constraint is imposed by generalizing the Lagrange multiplier method using the Karush-Kuhn-Tucker (KKT) conditions (see, e.g., Karush, 1939; Chen and Ye, 2022). Specifically, let ∂P(a|θ, bc ) h(θ) = sup − , ∂a a,c

(12)

denote the minimum derivative of the maturity probability with respect to age a across all the ages and cohort segments (c). The generalized Lagrangian function is then formulated as L(θ, η) = − log{l(θ|d1 , . . . , dn )} + λ(W⊤ W + ∆⊤ ∆) + ηh(θ),

(13)

where η ≥ 0 is the KKT multiplier corresponding to the inequality constraint. According to the KKT conditions, the optimization procedure begins by setting η = 0 to find the initial minimizer θ̂ of the generalized Lagrangian function (13) with respect to θ. The algorithm then evaluates whether the monotonicity constraint h(θ̂) ≤ 0 is satisfied. If the condition is violated, η is increased by a small increment and θ̂ is re-estimated; this process is repeated iteratively until the condition h(θ̂) ≤ 0 holds. TMB efficiently computes the generalized Lagrangian function (13) and its gradient with respect to θ for optimization via nlminb function in R. As a baseline for comparison, we consider a GLMM, in which the conditional maturation

probability for the ith fish, given the REs, is modeled as: Pi (ai |β0 , βsci , bci ) = inv logit(β0 + βsci ai + bci ),

(14)

where sci is the cohort seg number for cohort ci , and β0 and βsci ’s are model parameters to estimate. The principal distinction between the NMM (11) and the GLMM (14) lies in the age effect: while the latter incorporates a strictly linear predictor with respect to age, the former accommodates flexible, nonlinear functional relationships. We used the 3L data (31,427 observations) as the training set, the 3N data (23,392 observations) as the validation set, and the 3O data (24,104 observations) as the test set. We evaluated model performance on the validation and test sets using the average negative loglikelihood, referred to as the binary cross-entropy (BCE) in the machine learning literature. For instance, the BCE for the validation set, denoted by BCEvalid , is defined as Nvalid h i 1 X BCEvalid = − di log{Pi(ai | φ̂, b̂ci )} + (1 − di ) log{1 − Pi(ai | φ̂, b̂ci )} , (15) Nvalid i=1 where i = 1, . . . , Nvalid indexes the validation observations, and φ̂ and b̂ci represent the estimated model parameters and REs, respectively. We implemented both the NMM and GLMM using TMB, which efficiently evaluates the marginal objective function and its gradient. These quantities were then supplied to nlminb() for optimization of the marginal objective function. During optimization, the models were fitted to the training data, while BCEvalid was recorded at each optimization step. We identified the point at which BCEvalid ceased to improve as the point of best generalization. The corresponding parameter and RE estimates were subsequently taken as the final estimates. Prior to this final fit, we performed a grid search over the regularization hyperparameter λ, and at each grid point, the KKT multiplier η was determined using the algorithm described previously. This procedure yielded λ = 0.001 and η = 1.7. Based on the optimized parameter and RE estimates, we evaluated the BCE on the test set, BCEtest , to compare the predictive performance of the NMM and GLMM. The NMM achieved a slightly lower BCEtest of 0.249, compared with 0.250 for the GLMM. 15

TMB automatically predicts the REs using their posterior modes and predicts functions of the REs and fixed parameters using plug-in predictors. Figure 4 compares the predicted age-specific maturation probabilities for the 50th cohort between the NMM and the GLMM. While the GLMM yields predictions that are symmetric around the inflection point, the NMM accommodates mild asymmetry to better capture the underlying data.

1.00

Maturation Probability

0.75

0.50

0.25

Model GLMM NMM

0.00 2

Age

Figure 4: Predicted probabilities of maturation by age for the 50th cohort, using the NMM and GLMM.

Conclusion

Due to its enormous implementation and computation complexity, efficient software and hardware tools are always one of the central problem in artificial intelligence. In this work, we propose to implement neural network mixed-effects models (NMMs) in statistical modeling using Template Model Builder (TMB). In this approach, users need to specify only the negative joint log-likelihood and any regularization terms. The framework automatically integrates out random effects and evaluates the marginal objective function alongside its exact gradients, eliminating the need for manual derivations or ad hoc approximations. Replicating the simulation studies of Mandel et al. (2023), the results in Section 3 highlight several advantages of our proposed methodology. First, adopting a full likelihood formulation instead of a quasi-likelihood approach increases statistical efficiency. Second, by fully automating the derivation process, our method eliminates human error and avoids the heuristic omissions or approximations typically needed for manual solutions. As a result, our approach offers enhanced estimation accuracy, improved implementation efficiency, and superior scalability toward complex NMMs. The empirical results in Section 4 demonstrate the simplicity of implementing constrained NMMs within the TMB framework. Specifically, TMB enables the efficient computation of the monotonicity-constrained objective function and its gradients for the downstream optimization algorithm. Furthermore, even without extensive architecture tuning, this application underscores the representational capacity of NMMs in modeling intricate, nonlinear input–output dynamics. In the current implementation, feedforward neural networks are coded manually; however, extending this manual approach to convolutional, recurrent, or transformer architectures is impractical. Immediate future work will focus on integrating Keras (see, e.g., Chapter 3 of Chollet et al., 2022), an R interface to TensorFlow, with RTMB (Kristensen, 2024), which provides native R automatic differentiation for TMB without C++ com-

pilation. Synthesizing these tools will allow us to construct sophisticated deep learning mixed-effects models. A secondary direction involves extending these models to incorporate complex spatiotemporal structures via random effects. Additionally, migrating our implementation of NMMs from a strictly CPU-bound environment to a hybrid CPU/GPU workflow will optimize computational performance and support the increased overhead of these advanced architectures.

Codes availability The TMB code supporting the findings of this study is publicly available at https://github.com/ZhengNan18/Implementing-neural-network-mixed-effects-models-in-TemplateModel-Builder-TMB-

Acknowledgments Research funding to NZ was provided by the Natural Sciences and Engineering Research Council of Canada [RGPIN-2024-04746]. Research funding to NC was provided by 1) the Ocean Choice International Industry Research Chair program at the Marine Institute of Memorial University of Newfoundland, 2) the Ocean Frontier Institute, through an award from the Canada First Research Excellence Fund, and 3) the Natural Sciences and Engineering Research Council of Canada [RGPIN-2016-04307].

References Bengio, Y. (2009). Learning deep architectures for ai. Foundations and Trends® in Finance 2 (1), 1–127.

Breslow, N. and D. Clayton (1993). Approximate inference in generalized linear mixed models. Journal of the American statistical Association 88 (421), 9–25. Chen, D. and W. Ye (2022). Monotonic neural additive models: Pursuing regulated machine learning models for credit scoring. In Proceedings of the third ACM international conference on AI in finance, pp. 70–78. Chollet, F., T. Kalinowski, and J. J. Allaire (2022). Deep learning with R. Simon and Schuster. Goodfellow, I., Y. Bengio, and A. Courville (2016). Deep Learning. MIT Press. Hastie, T., R. Tibshirani, J. Friedman, et al. (2009). The elements of statistical learning. Hornik, K., M. Stinchcombe, and H. White (1989). Multilayer feedforward networks are universal approximators. Neural networks 2 (5), 359–366. Karush, W. (1939). Minima of functions of several variables with inequalities as side constraints. M. Sc. Dissertation. Dept. of Mathematics, Univ. of Chicago. Kristensen, K. (2024). RTMB: R Bindings for ’TMB’. R package version 1.5. Kristensen, K., A. Nielsen, C. W. Berg, H. Skaug, and B. Bell (2016). Tmb: automatic differentiation and laplace approximation. Journal of Statistical Software 70 (5), 1–21. Krogh, A. and J. Hertz (1991). A simple weight decay can improve generalization. Advances in neural information processing systems 4. LeCun, Y., Y. Bengio, and G. Hinton (2015). Deep learning. nature 521 (7553), 436–444. Mandel, F., R. P. Ghosh, and I. Barnett (2023). Neural networks for clustered and longitudinal data using mixed effects models. Biometrics 79 (2), 711–721.

R Core Team (2024). R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. Tandon, R., S. Adak, and J. A. Kaye (2006). Neural networks for longitudinal studies in alzheimer’s disease. Artificial intelligence in medicine 36 (3), 245–255. Thygesen, U. H., C. M. Albertsen, C. W. Berg, K. Kristensen, and A. Nielsen (2017). Validation of ecological state space models using the laplace approximation. Environmental and Ecological Statistics 24 (2), 317–339. Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society Series B: Statistical Methodology 58 (1), 267–288. Tran, M.-N., N. Nguyen, D. Nott, and R. Kohn (2017). Random effects models with deep neural network basis functions: Methodology and computation. Technical report.

Implementing neural network mixed-effects models in Template Model Builder (TMB)

Analyze & explore related research in Conceptio

Related documents