Evidence For Model Research Articles

Exact gradient evaluation for adaptive quadrature approximate marginal likelihood in mixed models for grouped data

Exact gradient evaluation for adaptive quadrature approximate marginal likelihood in mixed models for grouped data

Psychometric validation of the French version of two scales measuring general (HLS19-Q12) and navigational (HLS19-NAV) health literacy using the Rasch model

BackgroundTo examine the psychometrics properties of the French version of two scales measuring general (HLS19-Q12) and navigational (HLS19-NAV) health literacy (HL) using validation methods based on modern psychometric test theories: a Rasch model analysis.MethodsThe data on representative samples of the French adult population came from the Health Literacy Survey (N = 2 003), conducted in France in two waves (2020 and 2021), and from the third wave of SLAVACO study (N = 2 022), conducted in December 2021. A Rasch analysis was performed using a partial credit model with marginal maximum likelihood estimation adapted to polytomous data. Category probability curves were used to examine whether participants found it consistently difficult to distinguish between response options. Differential item functioning (DIF) was also examined for each item in the two scales as a function of age, gender, and economic status. A Person Separation Index (PSI) of 0.7 has conventionally been considered to be the minimum acceptable PSI level. Chi-square fit statistics, expressed as infit and outfit mean square statistics, were calculated for each item in the two scales. An acceptable fit corresponded to average values between 0.7 and 1.3. The dimensionality of the scales was assessed using Principal Components Analysis (PCA) on the residual correlation matrix to identify whether items are locally independent. Additionally, Smith’s method was employed to detect potential multidimensionality.ResultsThe PSI values for both scales were greater than 0.90. No disordered categories were observed. No evidence of significant DIF was found when associated with gender and economic status. However, two items for both scales had a DIF which was dependent on age but did not have a significant impact on scale scores between younger and older people. Minor local dependence was noted on the HLS19-NAV scale. However, a PCA and Smith’s method demonstrated the unidimensionality of the French version of the scales measuring general and navigational HL.ConclusionsOur results, based on a rigorous statistical analysis, verified the psychometric parameters of the French version of the HLS19-Q12 and HLS19-NAV scales. Despite strong correlation between both scales, each measured a different latent trait. Moreover, despite the fact that presence of DIF was weak, attention must nonetheless be paid when comparing scale response scores between young and older respondents.

Latent Growth Model with Categorical Indicators: Comparison of the Different Estimators

Despite its potential usage in social and behavioral sciences, there has been limited research on Latent Growth Models (LGM) with categorical indicators. The categorical nature of binary or ordinal variables violates the multivariate normality assumption of the traditional maximum likelihood (ML) estimator. Alternative estimators can be employed, such as marginal maximum likelihood (MML), categorical weighted least squares (c-WLS), and categorical diagonally weighted least squares (c-DWLS). A Monte Carlo simulation study was conducted to compare the performance of the three estimators in the context of categorical LGM. The results revealed that the MML estimator performed the best for the correctly specified linear growth model, whereas the robust c-DWLS or c-WLS estimator excelled in correctly specified quadratic growth model. The robust c-DWLS and c-WLS estimators were more sensitive to model misspecifications. The findings of the present study provide important implications and practical recommendations for empirical researchers using categorical LGM.

Handling missing data in variational autoencoder based item response theory.

Recently Variational Autoencoders (VAEs) have been proposed as a method to estimate high dimensional Item Response Theory (IRT) models on large datasets. Although these improve the efficiency of estimation drastically compared to traditional methods, they have no natural way to deal with missing values. In this paper, we adapt three existing methods from the VAE literature to the IRT setting and propose one new method. We compare the performance of the different VAE-based methods to each other and to marginal maximum likelihood estimation for increasing levels of missing data in a simulation study for both three- and ten-dimensional IRT models. Additionally, we demonstrate the use of the VAE-based models on an existing algebra test dataset. Results confirm that VAE-based methods are a time-efficient alternative to marginal maximum likelihood, but that a larger number of importance-weighted samples are needed when the proportion of missing values is large.

Modelling of the time to death of breast cancer patients at Hiwot Fana Specialized University Hospital

Breast cancer is the most common cause of cancer death and is a frequently diagnosed cancer among women worldwide. It is becoming a challenging health condition in Ethiopia with a high rate of morbidity and mortality. The main aim of this study was to model the time to death in breast cancer patients at Hiwot Fana Specialized University Hospital. A retrospective cohort study was carried out from April 1st, 2020, to April 1st, 2023, and 296 women were included in the study. We used nonparametric methods and Bayesian accelerated failure time models (with Laplace approximation) to identify risk factors and choose a model fitting breast cancer patient data. Model comparison was performed using the marginal likelihood, deviance information criterion and Watanabe Akaike information criterion. From the total of 296 patients in the study, 56 (18.9%) died. The estimated median survival time was 33 months. The log-rank test showed that age group, stage, alcohol consumption, smoking habit, and comorbidity were potential risk factors associated with the time to death in breast cancer patients at the 5% level of significance. The Bayesian Weibull accelerated failure time model was found to be the best fitted model for predicting the survival time of patients with minimum DIC (520.39) and WAIC (521.59) values. The final Bayesian Weibull AFT model with the integrated nested Laplace approximation estimation technique revealed that age group, stage, alcohol consumption, smoking habit, and comorbidity were significantly associated with the time to death in breast cancer patients. Individuals older than 65 years, with stage IV disease, drinking alcohol, smoking cigarettes and having comorbidities had shortened survival times in patients with breast cancer. Hence, Hiwot Fana Specialized University Hospital and related bodies should work on awareness creation to reduce smoking habits and alcohol use as well as give due attention to elderly and stage IV breast cancer patients during intervention.

Earthquake Predictability and Forecast Evaluation Using Likelihood-Based Marginal and Conditional Scores

Abstract Earthquake probability forecasts are typically based on simulations of seismicity generated by statistical (point process) models or direct calculation when feasible. To systematically assess various aspects of such forecasts, the Collaborative Studies on Earthquake Predictability testing center has utilized N- (number), M- (magnitude), S- (space), conditional likelihood-, and T- (Student’s t) tests to evaluate earthquake forecasts in a gridded space–time range. This article demonstrates the correct use of point process likelihood to evaluate forecast performance covering marginal and conditional scores, such as numbers, occurrence times, locations, magnitudes, and correlations among space–time–magnitude cells. The results suggest that for models that only rely on the internal history but not on external observation to do simulation, such as the epidemic-type aftershock sequence model, test and scoring can be rigorously implemented via the likelihood function. Specifically, gridding the space domain unnecessarily complicates testing, and evaluating spatial forecasting directly via marginal likelihood might be more promising.

Iterative Methods for Vecchia-Laplace Approximations for Latent Gaussian Process Models

Latent Gaussian process (GP) models are flexible probabilistic non-parametric function models. Vecchia approximations are accurate approximations for GPs to overcome computational bottlenecks for large data, and the Laplace approximation is a fast method with asymptotic convergence guarantees to approximate marginal likelihoods and posterior predictive distributions for non-Gaussian likelihoods. Unfortunately, the computational complexity of combined Vecchia-Laplace approximations grows faster than linearly in the sample size when used in combination with direct solver methods such as the Cholesky decomposition. Computations with Vecchia-Laplace approximations can thus become prohibitively slow precisely when the approximations are usually the most accurate, i.e., on large data sets. In this article, we present iterative methods to overcome this drawback. Among other things, we introduce and analyze several preconditioners, derive new convergence results, and propose novel methods for accurately approximating predictive variances. We analyze our proposed methods theoretically and in experiments with simulated and real-world data. In particular, we obtain a speed-up of an order of magnitude compared to Cholesky-based calculations and a threefold increase in prediction accuracy in terms of the continuous ranked probability score compared to a state-of-the-art method on a large satellite data set. All methods are implemented in a free C++ software library with high-level Python and R packages. 4

Assessing the Adequacy of Morphological Models using Posterior Predictive Simulations.

Reconstructing the evolutionary history of different groups of organisms provides insight into how life originated and diversified on Earth. Phylogenetic trees are commonly used to estimate this evolutionary history. Within Bayesian phylogenetics a major step in estimating a tree is in choosing an appropriate model of character evolution. While the most common character data used is molecular sequence data, morphological data remains a vital source of information. The use of morphological characters allows for the incorporation fossil taxa, and despite advances in molecular sequencing, continues to play a significant role in neontology. Moreover, it is the main data source that allows us to unite extinct and extant taxa directly under the same generating process. We therefore require suitable models of morphological character evolution, the most common being the Mk Lewis model. While it is frequently used in both palaeobiology and neontology, it is not known whether the simple Mk substitution model, or any extensions to it, provide a sufficiently good description of the process of morphological evolution. In this study we investigate the impact of different morphological models on empirical tetrapod data sets. Specifically, we compare unpartitioned Mk models with those where characters are partitioned by the number of observed states, both with and without allowing for rate variation across sites and accounting for ascertainment bias. We show that the choice of substitution model has an impact on both topology and branch lengths, highlighting the importance of model choice. Through simulations, we validate the use of the model adequacy approach, posterior predictive simulations, for choosing an appropriate model. Additionally, we compare the performance of model adequacy with Bayesian model selection. We demonstrate how model selection approaches based on marginal likelihoods are not appropriate for choosing between models with partition schemes that vary in character state space (i.e., that vary in Q-matrix state size). Using posterior predictive simulations, we found that current variations of the Mk model are often performing adequately in capturing the evolutionary dynamics that generated our data. We do not find any preference for a particular model extension across multiple data sets, indicating that there is no 'one size fits all' when it comes to morphological data and that careful consideration should be given to choosing models of discrete character evolution. By using suitable models of character evolution, we can increase our confidence in our phylogenetic estimates, which should in turn allow us to gain more accurate insights into the evolutionary history of both extinct and extant taxa.

Generalized Additive Models

Generalized additive models are generalized linear models in which the linear predictor includes a sum of smooth functions of covariates, where the shape of the functions is to be estimated. They have also been generalized beyond the original generalized linear model setting to distributions outside the exponential family and to situations in which multiple parameters of the response distribution may depend on sums of smooth functions of covariates. The widely used computational and inferential framework in which the smooth terms are represented as latent Gaussian processes, splines, or Gaussian random effects is reviewed, paying particular attention to the case in which computational and theoretical tractability is obtained by prior rank reduction of the model terms. An empirical Bayes approach is taken, and its relatively good frequentist performance discussed, along with some more overtly frequentist approaches to model selection. Estimation of the degree of smoothness of component functions via cross validation or marginal likelihood is covered, alongside the computational strategies required in practice, including when data and models are reasonably large. It is briefly shown how the framework extends easily to location-scale modeling, and, with more effort, to techniques such as quantile regression. Also covered are the main classes of smooths of multiple covariates that may be included in models: isotropic splines and tensor product smooth interaction terms.

Analysis of recurrent event data with spatial random effects using a Bayesian approach.

Recurrent event data, which represent the occurrence of repeated incidences, are common in observational studies. Furthermore, collecting possible spatial correlations in health and environmental data is likely to provide more information for risk prediction. This article proposes a comprehensive proportional intensity model considering spatial random effects for recurrent event data using a Bayesian approach. The spatial information for areal data (where the spatial location is known up to a geographic unit such as a county) and georeferenced data (where the location is exactly observed) is examined. A traditional constant baseline intensity function, as well as a flexible piecewise constant baseline intensity function, are both under consideration. To estimate the parameters, a Markov chain Monte Carlo method with the Metropolis-Hastings algorithm and the adaptive Metropolis algorithm are applied. To assess the performance of model fitting, the deviance information criterion and log pseudo marginal likelihood are proposed. Overall, simulation studies demonstrate that the proposed model is significantly better than models that do not consider spatial effects if spatial correlations exist. Finally, our approach is implemented using a dataset related to the recurrence of cardiovascular diseases, which incorporates spatial information.

Patient- and Areal-Level Risk Factors Associated With Lung Cancer Mortality in Victoria, Australia: A Bayesian Spatial Survival Analysis.

In Australia, lung cancer is the leading cause of cancer-related deaths. In Victoria, the mortality risk is assumed to vary across Local Government Areas (LGAs) due to variations in socioeconomic advantage, remoteness, and healthcare accessibility. Thus, we applied Bayesian spatial survival models to examine the geographic variation in lung cancer survival in Victoria. Data on lung cancer cases were extracted from the Victorian Lung Cancer Registry (VLCR). To account for spatial dependence and risk factors of survival in lung cancer patients, we employed a Bayesian spatial survival model. Conditional Autoregressive (CAR) prior was assigned to model the spatial dependence. Deviance Information Criterion (DIC), Watanabe Akaike Information Criterion (WAIC), and Log Pseudo Marginal Likelihood (LPML) were used for model comparison. In the final best-fitted model, the Adjusted Hazard Ratio (AHR) with the 95% Credible Interval (CrI) was reported. The outcome variable was the survival status of lung cancer patients, defined as whether they survived or died during the follow-up period (death was our interest). Our study revealed substantial variations in lung cancer mortality in Victoria. Poor Eastern Cooperative Oncology Group (ECOG) performance status, diagnosed at a regional hospital, Small Cell Lung Cancer (SCLC), advanced age, and advanced clinical stage were associated with a higher risk of mortality, whereas being female, presented at Multidisciplinary Team (MDT) meeting, and diagnosed at a metropolitan private hospital were significantly associated with a lower risk of mortality. Identifying geographical disparities in lung cancer survival may help shape healthcare policy to implement more targeted and effective lung cancer care services.

Cross-fade sampling: extremely efficient Bayesian inversion for a variety of geophysical problems

SUMMARY This paper introduces cross-fade sampling, a computationally efficient Markov Chain Monte Carlo simulation method that uses a semi-analytical approach to quickly solve Bayesian inverse problems that do not themselves have an analytical solution. Cross-fading is efficient in two ways. First, it requires fewer samples to obtain the same quality simulation of the target probability density function (PDF). Secondly, it is much faster to evaluate the posterior probability of each sample than conventional sampling methods for simulating Bayesian posterior PDFs. Conventional methods require evaluating the prior probability (which describes your a priori constraints) and data likelihood (which describes the fit between the observations and the predictions of the model) for each sample model. However, cross-fading does not require evaluating the data likelihood, meaning that ‘big data’ can be fit with zero additional computational cost. Further, the cross-fading approach can be used to calculate the marginal likelihood associated with a model design, facilitating model comparison and Bayesian model averaging. Topics covered in this paper include derivation of the cross-fade approach and how it can be used to simulate Bayesian posterior PDFs and compute the marginal likelihood, discussion of the class of problems to which cross-fading can be applied (with examples from earthquake statistics, earthquake ground motion modelling, volcanic eruption forecasting, and finite fault slip modelling), demonstration of efficiency relative to existing sampling methods and discussion of how cross-fading can be used to account for prediction errors (i.e. epistemic errors) as part of the geophysical inverse problem.

Spike Count Analysis for MultiPlexing Inference (SCAMPI).

Understanding how neurons encode multiple simultaneous stimuli is a fundamental question in neuroscience. We have previously introduced a novel theory of stochastic encoding patterns wherein a neuron's spiking activity dynamically switches among its constituent single-stimulus activity patterns when presented with multiple stimuli (Groh et al., 2024). Here, we present an enhanced, comprehensive statistical testing framework for such " multiplexing " or " code juggling ". Our new approach evaluates whether dual-stimulus responses can be accounted for as mixtures of Poissons either anchored to or bounded by single-stimulus benchmarks. Our enhanced framework improves upon previous methods in two key ways. First, it introduces a stronger set of foils for multiplexing, including an "overreaching" category that captures overdispersed activity patterns unrelated to the single-stimulus benchmarks, reducing false detection of multiplexing/code-juggling. Second, it detects faster fluctuations - i.e. at sub-trial timescales - that would have been overlooked before. We utilize a Bayesian inference framework, considering the hypothesis with the highest posterior probability as the winner, and employ predictive recursion marginal likelihood method for the involving nonparametric density estimation. Reanalysis of previous findings confirms the general observation of "code juggling" and indicates that such juggling may well occur on faster timescales than previously suggested. We further confirm that juggling is more prevalent in (a) the inferotemporal face patch system for combinations of face stimuli than for faces and non-face objects; and (b) the primary visual cortex for distinct vs fused objects.

The irtQ R package: a user-friendly tool for item response theorybased test data analysis and calibration.

Computerized adaptive testing (CAT) has become a widely adopted test design for high-stakes licensing and certification exams, particularly in the health professions in the United States, due to its ability to tailor test difficulty in real time, reducing testing time while providing precise ability estimates. A key component of CAT is item response theory (IRT), which facilitates the dynamic selection of items based on examinees' ability levels during a test. Accurate estimation of item and ability parameters is essential for successful CAT implementation, necessitating convenient and reliable software to ensure precise parameter estimation. This paper introduces the irtQ R package (http://CRAN.R-project.org/), which simplifies IRTbased analysis and item calibration under unidimensional IRT models. While it does not directly simulate CAT, it provides essential tools to support CAT development, including parameter estimation using marginal maximum likelihood estimation via the expectation-maximization algorithm, pretest item calibration through fixed item parameter calibration and fixed ability parameter calibration methods, and examinee ability estimation. The package also enables users to compute item and test characteristic curves and information functions necessary for evaluating the psychometric properties of a test. This paper illustrates the key features of the irtQ package through examples using simulated datasets, demonstrating its utility in IRT applications such as test data analysis and ability scoring. By providing a user-friendly environment for IRT analysis, irtQ significantly enhances the capacity for efficient adaptive testing research and operations. Finally, the paper highlights additional core functionalities of irtQ, emphasizing its broader applicability to the development and operation of IRT-based assessments.

Statistical physical view of statistical inference in Bayesian linear regression model.

This paper considers similarities between statistical physics and Bayes inference through the Bayesian linear regression model. Some similarities have been discussed previously, such as the analogy between the marginal likelihood in Bayes inference and the partition function in statistical mechanics. In particular, this paper considers the proposal to associate discrete sample size with inverse temperature [C. H. LaMont and P. A. Wiggins, Phys. Rev. E 99, 052140 (2019)2470-004510.1103/PhysRevE.99.052140]. The previous study suggested that incorporating this similarity motivates the derivation of analogs of thermodynamic functions such as energy and entropy. The study also anticipated that those analogous functions have potential to describe Bayes estimation from physical points of view and to provide physical insights into mechanisms of estimation. This paper incorporates a macroscopic perspective as an asymptotics similar to the thermodynamic limit into the previous suggestion. Its motivation stems from the statistical mechanical concept of deriving thermodynamic functions that characterize macroscopic properties of macroscopic systems. This incorporation not only allows analogs of macroscopic thermodynamic functions to be considered but also suggests a candidate for an analog of inverse temperature with continuity, which is partly consistent with the previous proposal to associate the discrete sample size with inverse temperature. On the basis of this suggestion, we analyze analogs of macroscopic thermodynamic functions for a Bayesian linear regression model which is the basis of various machine learning models. We further investigate, through the behavior of these functions, how Bayes estimation is described from the perspective of physics and what kind of physical insight is obtained. As a result, the estimation of regression coefficients, which is the primary task of regression, appears to be described by the physical picture of balance between decreasing energy and increasing entropy as in equilibrium states of thermodynamic systems. More specifically we observe the physical view of Bayes inference as follows: the estimation succeeds where the effect of decreasing energy is dominant at low temperature. On the other hand, the estimation fails where the effect of increasing entropy is dominant at high temperature.

Handling Out‐of‐Sample Areas to Estimate the Unemployment Rate at Local Labour Market Areas in Italy

SummaryUnemployment rate estimates for small areas are used to efficiently support the distribution of services and the allocation of resources, grants and funding. A Fay–Herriot type model is the most used tool to obtain these estimates. Under this approach out‐of‐sample areas require some synthetic estimates. As the geographical context is extremely important for analysing local economies, in this paper, we allow for area random effects to be spatially correlated. The spatial model parameters are estimated by a marginal likelihood method and are used to predict in‐sample as well as out‐of‐sample areas. Extensive simulation experiments are used to assess the impact of the auto‐regression parameter and of the rate of out‐of‐sample areas on the performance of this approach. The paper concludes with an illustrative application on real data from the Italian Labour Force Survey in which the estimation of the unemployment rate in each Local Labour Market Area is addressed.

An independent base composition of each rate class for improved likelihood-based phylogeny estimation; the 5rf model.

The combination of a time reversible Markov process with a "hidden" mixture of gamma distributed relative site rates plus invariant sites have become the most favoured options for likelihood and other probabilistic models of nucleotide evolution (e.g., tr4gi which approximates a gamma with four rate classes). However, these models assume a homogeneous and stationary distribution of nucleotide (character or base) frequencies. Here, we explore the potential benefits and pitfalls of allowing each rate category (rate class) of a 4gi mixture model to have its own base frequencies. This is achieved by starting each of the five rate classes, at the tree's root, with its own free choice of nucleotide frequencies to create a 4gi5rf model or a 5rf model in shorthand. We assess the practical identifiability of this approach with a BEAST 2 implementation, aiming to determine if it can accurately estimate credibility intervals and expected values for a wide range of plausible parameter values. Practical identifiability, as distinguished from mathematical identifiability, gauges the model's ability to identify parameters in real-world scenarios, as opposed to theoretically with infinite data. One of the most common types of phylogenetic data is mitochondrial DNA (mtDNA) protein coding sequence. It is often assumed current models analyse robustly such data and that higher likelihood/posterior probability models do better. However, this abstract shows that vertebrate mtDNA remains a very difficult type of data to fully model, and that dramatically higher likelihoods do not mean a model is measurably more accurate with respect to recovering key parameters of biological interest (e.g., monophyletic groups, their support and their ages). The 4gi5rf model considerably improves marginal likelihoods and seems to reverse some apparent errors exacerbated by the 4gi model, while introducing others. Problems appear to be linked to non-stationary DNA repair processes that alter the mutation/substitution spectra across lineages and time. We also show such problems are not unique to mtDNA and are encountered in analysing nuclear sequences. Non-stationarity of DNA repair processes mutation/substitution spectra thus pose an active challenge to obtaining reliable inferences of relationships and divergence times near the root of placental mammals, for example. An open source implementation is available under the LGPL 3.0 license in the beastbooster package for BEAST 2, available from https://github.com/rbouckaert/beastbooster.

Bayes factors for longitudinal model assessment via power posteriors

Bayes factor, defined as the ratio of the marginal likelihood functions of two competing models, is the natural Bayesian procedure for model selection. Marginal likelihoods are usually computationally demanding and complex. This scenario is particularly cumbersome in linear mixed models (LMMs) because the marginal likelihood functions involve integrals of large dimensions determined by the number of parameters and the number of random effects, which in turn increase with the number of individuals in the sample. Power posterior is an attractive proposal in the context of the Markov chain Monte Carlo algorithms that allows expressing marginal likelihoods as one-dimensional integrals over the unit range. This paper explores the use of power posteriors in LMMs and discusses its behavior through three simulation studies and a real data set on European sardine landings in the Mediterranean Sea.

Ab initio uncertainty quantification in scattering analysis of microscopy.

Estimating parameters from data is a fundamental problem in physics, customarily done by minimizing a loss function between a model and observed statistics. In scattering-based analysis, it is common to work in the reciprocal space. Researchers often employ their domain expertise to select a specific range of wave vectors for analysis, a choice that can vary depending on the specific case. We introduce another paradigm that defines a probabilistic generative model from the beginning of data processing and propagates the uncertainty for parameter estimation, termed the ab initio uncertainty quantification (AIUQ). As an illustrative example, we demonstrate this approach with differential dynamic microscopy (DDM) that extracts dynamical information through minimizing a loss function for the squared differences of the Fourier-transformed intensities, at a selected range of wave vectors. We first show that the conventional way of estimation in DDM is equivalent to fitting a temporal variogram in the reciprocal space using a latent factor model as the generative model. Then we derive the maximum marginal likelihood estimator, which optimally weighs the information at all wave vectors, therefore eliminating the need to select the range of wave vectors. Furthermore, we substantially reduce the computational cost of computing the likelihood function without approximation, by utilizing the generalized Schur algorithm for Toeplitz covariances. Simulated studies of a wide range of dynamical systems validate that the AIUQ method improves estimation accuracy and enables model selection with automated analysis. The utility of AIUQ is also demonstrated by three distinct sets of experiments: first in an isotropic Newtonian fluid, pushing limits of optically dense systems compared to multiple particle tracking; next in a system undergoing a sol-gel transition, automating the determination of gelling points and critical exponent; and lastly, in discerning anisotropic diffusive behavior of colloids in a liquid crystal. These studies demonstrate that the new approach does not require manually selecting the wave vector range and enables automated analysis.

Sparse Bayesian learning using TMB (Template Model Builder)

Sparse Bayesian Learning, and more specifically the Relevance Vector Machine (RVM), can be used in supervised learning for both classification and regression problems. Such methods are particularly useful when applied to big data in order to find a sparse (in weight space) representation of the model. This paper demonstrates that the Template Model Builder (TMB) is an accurate and flexible computational framework for implementation of sparse Bayesian learning methods.The user of TMB is only required to specify the joint likelihood of the weights and the data, while the Laplace approximation of the marginal likelihood is automatically evaluated to numerical precision. This approximation is in turn used to estimate hyperparameters by maximum marginal likelihood. In order to reduce the computational cost of the Laplace approximation we introduce the notion of an “active set” of weights, and we devise an algorithm for dynamically updating this set until convergence, similar to what is done in other RVM type methods. We implement two different methods using TMB; the RVM and the Probabilistic Feature Selection and Classification Vector Machine method, where the latter also performs feature selection. Experiments based on benchmark data show that our TMB implementation performs comparable to that of the original implementation, but at a lower implementation cost. TMB can also calculate model and prediction uncertainty, by including estimation uncertainty from both latent variables and the hyperparameters. In conclusion, we find that TMB is a flexible tool that facilitates implementation and prototyping of sparse Bayesian methods.