Streamlined Variational Inference for Higher Level Group-Specific Curve Models
M. Menictas, T.H. Nolan, D.G. Simpson, M.P. Wand
TL;DR
This paper develops a systematic approach to streamline variational inference for complex multi-level group-specific curve models, enabling analysis of hierarchical data such as ultrasound technology measurements.
Contribution
It introduces a structured method for variational inference in multi-level models, extending previous sparse matrix techniques to higher levels.
Findings
Effective inference for two- and three-level models demonstrated.
Framework facilitates analysis of hierarchical data structures.
Pattern generalizes to models with more than three levels.
Abstract
A two-level group-specific curve model is such that the mean response of each member of a group is a separate smooth function of a predictor of interest. The three-level extension is such that one grouping variable is nested within another one, and higher level extensions are analogous. Streamlined variational inference for higher level group-specific curve models is a challenging problem. We confront it by systematically working through two-level and then three-level cases and making use of the higher level sparse matrix infrastructure laid down in Nolan and Wand (2018). A motivation is analysis of data from ultrasound technology for which three-level group-specific curve models are appropriate. Whilst extension to the number of levels exceeding three is not covered explicitly, the pattern established by our systematic approach sheds light on what is required for even higher level…
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7| naïve | streamlined | naïve/streamlined | |
|---|---|---|---|
| 100 | 75 (1.21) | 0.748 (0.0334) | 100 |
| 200 | 660 (7.72) | 1.490 (0.0491) | 442 |
| 300 | 2210 (22.00) | 2.260 (0.0567) | 974 |
| 400 | 5180 (92.20) | 3.040 (0.0718) | 1700 |
| 500 | NA | 3.780 (0.0593) | NA |
| 0.635 | 2.900 | 16.90 | 95.00 |
| (0.183) | (0.391) | (1.92) | (4.92) |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsStatistical Methods and Inference · Soil Geostatistics and Mapping · Statistical Methods and Bayesian Inference
Streamlined Variational Inference for
Higher Level Group-Specific Curve Models
By M. Menictas1, T.H. Nolan1, D.G. Simpson2 and M.P. Wand1
University of Technology Sydney1 and University of Illinois2
11th March, 2019
Abstract
A two-level group-specific curve model is such that the mean response of each member of a group is a separate smooth function of a predictor of interest. The three-level extension is such that one grouping variable is nested within another one, and higher level extensions are analogous. Streamlined variational inference for higher level group-specific curve models is a challenging problem. We confront it by systematically working through two-level and then three-level cases and making use of the higher level sparse matrix infrastructure laid down in Nolan & Wand (2018). A motivation is analysis of data from ultrasound technology for which three-level group-specific curve models are appropriate. Whilst extension to the number of levels exceeding three is not covered explicitly, the pattern established by our systematic approach sheds light on what is required for even higher level group-specific curve models.
Keywords: longitudinal data analysis, multilevel models, panel data, mean field variational Bayes.
1 Introduction
We provide explicit algorithms for fitting and approximate Bayesian inference for multilevel models involving, potentially, thousands of noisy curves. The algorithms include covariance parameter estimation and allow for pointwise credible intervals around the fitted curves. Contrast function fitting and inference is also supported by our approach. Both two-level and three-level situations are covered, and a template for even higher level situations is laid down.
Models and methodology for statistical analyses of grouped data for which the basic unit is a noisy curve continues to be an important area of research. A driving force is rapid technological change which is resulting in the generation of curve-type data at fine resolution levels. Examples of such technology include accelerometers (e.g. Goldsmith et al., 2015) personal digital assistants (e.g. Trail et al., 2014) and quantitative ultrasound (e.g. Wirtzfeld et al., 2015). In some applications curve-type data have higher levels of grouping, with groups at one level nested inside other groups. Our focus here is streamlined variational inference for such circumstances.
Some motivating data is shown in Figure 1 from an experiment involving quantitative ultrasound technology. Each curve corresponds to a logarithmically transformed backscatter coefficient over a fine grid of frequency values for tumors in laboratory mice, with exactly one tumor per mouse. The backscatter/frequency curves are grouped according to one of 5 slices of the same tumor, corresponding to probe locations. The slices are grouped according to being from one of 10 tumors. We refer to such data as three-level data with frequency measurements at level 1, slices being the level 2 groups and tumors constituting the level 3 groups. The gist of this article is efficient and flexible variational fitting and inference for such data, that scales well to much larger multilevel data sets. Indeed, our algorithms are linear in the number of groups at both level 2 and level 3. Simulation study results given later in this article show that curve-type data with thousands of groups can be analyzed quickly using our new methodology. Depending on sample sizes and implementation language, fitting times range from a few seconds to a few minutes. In contrast, naïve implementations become infeasible when the number of groups are in the several hundreds due to storage and computational demands.
We work with a variant of group-specific curve models that at least go back to Donnelly, Laird & Ware (1995). Other contributions of this type include Brumback & Rice (1998), Verbyla et al. (1999), Wang (1998) and Zhang et al. (1998). The specific formulation that we use is that given by Durban et al. (2005) which involves an embedding within the class of linear mixed models (e.g. Robinson, 1991) with low-rank smoothing splines used for flexible function modelling and fitting.
Even though approximate Bayesian variational inference is our overarching goal, we also provide an important parallelism involving classical frequentist inference. Contemporary mixed model software such as nlme() (Pinheiro et al., 2018) and lme4() (Bates et al., 2015) in the R language provide streamlined algorithms for obtaining the best linear unbiased predictions of fixed and random effects in multilevel mixed models with details given in, for example, Pinheiro & Bates (2000). However, the sub-blocks of the covariance matrices required for construction of pointwise confidence interval bands around the estimated curves are not provided by such software. In the variational Bayesian analog, these sub-blocks are required for covariance parameter fitting and inference which, in turn, are needed for curve estimation. A significant contribution of this article is streamlined computation for both the best linear unbiased predictors and its corresponding covariance computation. Similar mathematical results lead to the mean field variational Bayesian inference equivalent. We present explicit ready-to-code algorithms for both two-level and three-level group-specific curve models. Extensions to higher level models could be derived using the blueprint that we establish here. Nevertheless, the algebraic overhead is increasingly burdensome with each increment in the number of levels. It is prudent to treat each multilevel case separately and here we already require several pages to cover two-level and three-level group-specific curve models. To our knowledge, this is the first article to provide streamlined algorithms for fitting three-level group-specific curve models.
Another important aspect of our group-specific curve fitting algorithms is the fact that they make use of the SolveTwoLevelSparseLeastSquares and SolveThreeLevelSparseLeastSquares algorithms developed for ordinary linear mixed models in Nolan et al. (2018). This realization means that the algorithms listed in Sections 2 and 3 are more concise and code-efficient: there is no need to repeat the implementation of these two fundamental algorithms for stable QR-based solving of higher level sparse linear systems. Sections S.11–S.12 of the web-supplement provide details on the SolveThreeLevelSparseLeastSquares and SolveThreeLevelSparseLeastSquares algorithms.
Section 2 deals with the two-level case and the three-level case is covered in Section 3. In Section 4 we provide some assessments concerning the accuracy and speed of the new variational inference algorithms.
2 Two-Level Models
The simplest version of group-specific curve models involves the pairs where is the th value of the predictor variable within the th group and is the corresponding value of the response variable. We let denote the number of groups and denote the number of predictor/response pairs within the th group. The Gaussian response two-level group specific curve model is
[TABLE]
where the smooth function is the global regression mean function and the smooth functions , , allow for flexible group-specific deviations from . As in Durban et al. (2005), we use mixed model-based penalized basis functions to model and the . Specifically,
[TABLE]
where \{z_{\mbox{\rm\tiny gbl,k}}(\cdot):1\leq k\leq K_{\mbox{\rm\tiny gbl}}\} and are suitable sets of basis functions. Splines and wavelet families are the most common choices for the z_{\mbox{\rm\tiny gbl,k}}(\cdot) and . In our illustrations and simulation studies we use the canonical cubic O’Sullivan spline basis as described in Section 4 of Wand & Ormerod (2008), which corresponds to a low-rank version of classical smoothing splines (e.g. Wahba, 1990). The variance parameters and control the effective degrees of freedom used for the global mean and group-specific deviation functions respectively. Lastly, is a unstructured covariance matrix for the coefficients of the group-specific linear deviations.
We also use the notation:
[TABLE]
for the vectors of predictors and responses corresponding to the th group. Notation such as z_{\mbox{\rm\tiny gbl,1}}(\boldsymbol{x}_{i}) denotes the vector containing z_{\mbox{\rm\tiny gbl,1}}(x_{ij}) values, .
2.1 Best Linear Unbiased Prediction
Model (1) is expressible as a Gaussian response linear mixed model as follows:
[TABLE]
where
[TABLE]
are the fixed effects design matrix and coefficients, corresponding to the linear component of . The random effects design matrix and corresponding random effects vector are partitioned according to
[TABLE]
where \boldsymbol{u}_{\mbox{\rm\tiny gbl}}=[u_{\mbox{\rm\tiny gbl,1}}\ \ldots\ u_{\mbox{\rm\tiny gbl,K_{\mbox{\rm\tiny gbl}}}}]^{T} are the coefficients corresponding to the non-linear component of , \boldsymbol{u}_{\mbox{\rm\tiny lin,i}}=[u_{\mbox{\rm\tiny lin,i0}}\ u_{\mbox{\rm\tiny lin,i1}}]^{T} are the coefficients corresponding to the linear component of and \boldsymbol{u}_{\mbox{\rm\tiny grp,i}}=[u_{\mbox{\rm\tiny grp,i1}}\ \ldots\ u_{\mbox{\rm\tiny grp,iK_{\mbox{\rm\tiny grp}}}}]^{T} are the coefficients corresponding to the non-linear component of , . In (4), \boldsymbol{Z}_{\mbox{\rm\tiny gbl}}\equiv\mathop{\mbox{\rm stack}}_{1\leq i\leq m}(\boldsymbol{Z}_{\mbox{\rm\tiny gbl,i}}) and the matrices \boldsymbol{Z}_{\mbox{\rm\tiny gbl,i}} and \boldsymbol{Z}_{\mbox{\rm\tiny grp,i}}, , contain, respectively, spline basis functions for the global mean function and the th group deviation functions . Specifically,
[TABLE]
for . The corresponding fixed and random effects vectors are
[TABLE]
Hence, the full random effects covariance matrix is
[TABLE]
Next define the matrices
[TABLE]
The best linear unbiased predictor of and corresponding covariance matrix are
[TABLE]
This covariance matrix grows quadratically in , so its storage becomes infeasible for large numbers of groups. However, only the following sub-blocks are required for adding pointwise confidence intervals to curve estimates:
[TABLE]
As in Nolan, Menictas & Wand (2019), we define the generic two-level sparse matrix to be determination of the vector which minimizes the least squares criterion
[TABLE]
with having the two-level sparse form
[TABLE]
In (10), for any , the matrices , \overset{\ \mbox{{\bullet}}}{\boldsymbol{B}}_{i} and each have the same number of rows. The numbers of columns in and \overset{\ \mbox{{\bullet}}}{\boldsymbol{B}}_{i} are arbitrary whereas the are column vectors. In addition to solving for , the sub-blocks of corresponding to the non-sparse regions of are included in our definition of a two-level sparse matrix least squares problem. Algorithm 2 of Nolan et al. (2018) provides a stable and efficient solution to this problem and labels it the SolveTwoLevelSparseLeastSquares algorithm. Section S.11 of the web-supplement contains details regarding this algorithm. In Nolan et al. (2018) we used SolveTwoLevelSparseLeastSquares for fitting two-level linear mixed models. However, precisely the same algorithm can be used for fitting two-level group-specific curve models because of:
Result 1**.**
Computation of and each of the sub-blocks of listed in (8) are expressible as solutions to the two-level sparse matrix least squares problem:
[TABLE]
where the non-zero sub-blocks and , according to the notation in (10), are for :
[TABLE]
with each of these matrices having rows and with having columns and \overset{\ \mbox{{\bullet}}}{\boldsymbol{B}}_{i} having columns. The solutions are
[TABLE]
and
[TABLE]
[TABLE]
A derivation of Result 1 is given in Section S.1 of the web-supplement. Algorithm 1 encapsulates streamlined best linear unbiased prediction computation together with coefficient covariance matrix sub-blocks of interest.
2.2 Mean Field Variational Bayes
We now consider the following Bayesian extension of (3) and (5):
[TABLE]
Here the vector and symmetric positive definite matrix are hyperparameters corresponding to the prior distribution on and
[TABLE]
are hyperparameters for the variance and covariance matrix parameters. Details on the Inverse G-Wishart distribution, and the Inverse- special case, are given in Section S.3 of the web-supplement. The auxiliary variable is defined so that has a Half- distribution with degrees of freedom parameter and scale parameter , with larger values of corresponding to greater noninformativity. Analogous comments apply to the other standard deviation parameters. Setting \nu_{\mbox{\rm\tiny{\boldsymbol{\Sigma}}}}=2 leads to the correlation parameter in having a Uniform distribution on (Huang & Wand, 2013).
Throughout this article we use generically to denote a density function corresponding to random quantities in Bayesian models such as (11). For example, denotes the prior density function of and denotes the density function of conditional on . Now consider the following mean field restriction on the joint posterior density function of all parameters in (11):
[TABLE]
Here, generically, each denotes an approximate posterior density function of the random vector indicated by its argument according to the mean field restriction (12). Then application of the minimum Kullback-Leibler divergence equations (e.g. equation (10.9) of Bishop, 2006) leads to the optimal -density functions for the parameters of interest being as follows:
[TABLE]
The optimal -density parameters are determined via an iterative coordinate ascent algorithm, with details given in Section S.5 of this article’s web-supplement. The stopping criterion is based on the variational lower bound on the marginal likelihood (e.g. Bishop, 2006; Section 10.2.2) and denoted . Its logarithmic form and derivation are given in Section S.6 of the web-supplement.
Note that updates for and may be written
[TABLE]
where
[TABLE]
For increasingly large numbers of groups the matrix approaches a size that is untenable for random access memory storage on standard 2020s workplace computers. However, only the following relatively small sub-blocks of are required for variational inference concerning the variance and covariance matrix parameters:
[TABLE]
For a streamlined mean field variational Bayes algorithm, we appeal to:
Result 2**.**
The mean field variational Bayes updates of and each of the sub-blocks of in (15) are expressible as a two-level sparse matrix least squares problem of the form:
[TABLE]
where the non-zero sub-blocks and , according to the notation in (10), are, for ,
[TABLE]
and
[TABLE]
with each of these matrices having rows and with having columns and \overset{\ \mbox{{\bullet}}}{\boldsymbol{B}}_{i} having columns. The solutions are
[TABLE]
[TABLE]
and
[TABLE]
Algorithm 2 utilizes Result 2 to facilitate streamlined computation of the variational parameters.
Lastly, we note that Algorithm 2 is loosely related to Algorithm 2 of Lee & Wand (2016). One difference is that we are treating the Gaussian, rather than Bernoulli, response situation here. In addition, we are using the recent sparse multilevel matrix results of Nolan & Wand (2018) which are amenable to higher level extensions, such as the three-level group specific curve model treated in Section 3.
2.3 Contrast Function Extension
In many curve-type data applications the data can be categorized as being from two or more types. Of particular interest in such circumstances are contrast function estimates and accompanying standard errors. The streamlined approaches used in Algorithms 1 and 2 still apply for the contrast function extension regardless of the number of categories. The two category situation, where there is a single contrast function, is described here. The extension to higher numbers of categories is straightforward.
Suppose that the pairs are from one of two categories, labeled and , and introduce the indicator variable data:
[TABLE]
Then penalized spline models for the global mean and deviation functions for each category are
[TABLE]
and
[TABLE]
This allows us to estimate the global contrast function
[TABLE]
The distributions on the random coefficients are
[TABLE]
and
[TABLE]
independently of each other. In this two-category extension, the matrix is an unstructured covariance matrix.
Algorithms 1 and 2 can be used to achieve streamlined fitting and inference for the contrast curve extension, but with key matrices having new definitions. Firstly, the , \boldsymbol{Z}_{\mbox{\rm\tiny gbl,i}} and \boldsymbol{Z}_{\mbox{\rm\tiny grp,i}} matrices need to become:
[TABLE]
[TABLE]
and
[TABLE]
where is the vector of values. In the case of best linear unbiased prediction the updates for the and \overset{\ \mbox{{\bullet}}}{\boldsymbol{B}}_{i} matrices in Algorithm 1 need to be replaced by:
[TABLE]
and the output coefficient vectors change to
[TABLE]
In the case of mean field variational Bayes the updates of the and \overset{\ \mbox{{\bullet}}}{\boldsymbol{B}}_{i} matrices in Algorithm 2 need to be replaced by:
[TABLE]
and
[TABLE]
A contrast curves adjustment to the mean field variational Bayes updates is also required for some of the covariance matrix parameters. However, these calculations are comparatively simple and analogous to those given in Section S.5.
We demonstrate the use of Algorithm 2 in this setting for data from a longitudinal study on adolescent somatic growth. More detail on this data can be found in Pratt et al. (1989). The variables of interest are
[TABLE]
for and . The subjects are categorized into black ethnicity and white ethnicity and comparison of mean height between the two populations is of interest. Algorithm 2 is seen to have good agreement with the data in each sub-panel of the top two plots in Figure 2. The bottom panels of Figure 2 show the estimated height gap between black and white adolescents as a function of age. For the females, there is a significant height difference only at 16-17 years old. Between 5 and 15 years, there is no obvious height difference. For the males, it is highest and (marginally) statistically significant up to about 14 years of age, peaking at 13 years of age. Between 17 and 20 years old there is no discernible height difference between the two populations.
3 Three-Level Models
The three-level version of group-specific curve models corresponds to curve-type data having two nested groupings. For example, the data in each panel of Figure 1 are first grouped according to slice, which is the level 2 group, and the slices are grouped according to tumor which is the level 3 group. We denote predictor/response pairs as where is the th value of the predictor variable in the th level 3 group and th level 2 group and is the corresponding value of the response variable. We let denote the number of level 3 groups, denote that number of level 2 groups in the th level 3 group and denote the the number of units within the th level 2 group. The Figure 1 data, which happen to be balanced, are such that
[TABLE]
The Gaussian response three-level group specific curve model for such data is
[TABLE]
where the smooth function is the global mean function, the functions, , allow for group-specific deviations according to membership of the th level 3 group and the , and allow for an additional level of group-specific deviations according to membership of the th level 2 group within the th level 3 group. The mixed model-based penalized spline models for these functions are
[TABLE]
with all random effect distributions independent of each other. For this three-level case we have three bases:
[TABLE]
The variance and covariance matrix parameters are analogous to the two-level model. For example, \boldsymbol{\Sigma}_{\mbox{\rm\tinyg}} and \boldsymbol{\Sigma}_{\mbox{\rm\tinyh}} are both unstructured matrices corresponding to the linear components of the and respectively.
The following notation is useful for setting up the required design matrices: if is a set of matrices each having the same number of columns then
[TABLE]
We then define, for and ,
[TABLE]
3.1 Best Linear Unbiased Prediction
Model (17) is expressible as a Gaussian response linear mixed model as follows:
[TABLE]
where the design matrices are
[TABLE]
and
[TABLE]
where
[TABLE]
and the matrices \boldsymbol{Z}_{\mbox{\rm\tiny gbl,ij}}, \boldsymbol{Z}^{\mbox{\rm\tinyg}}_{\mbox{\rm\tiny grp,ij}} and \boldsymbol{Z}^{\mbox{\rm\tinyh}}_{\mbox{\rm\tiny grp,ij}}, , , contain, respectively, spline basis functions for the global mean function , the th level one group deviation functions and th level two group deviation functions . Specifically,
[TABLE]
The fixed and random effects vectors are
[TABLE]
with \boldsymbol{u}^{\mbox{\rm\tinyg}}_{\mbox{\rm\tiny grp,i}}, \boldsymbol{u}^{\mbox{\rm\tinyh}}_{\mbox{\rm\tiny lin,ij}} and \boldsymbol{u}^{\mbox{\rm\tinyh}}_{\mbox{\rm\tiny grp,ij}} defined similarly and the covariance matrix of is
[TABLE]
We define matrices in a similar way to what is given in (6). The best linear unbiased predictor of and corresponding covariance matrix are as shown in (7), but, with entries as described in this section. This covariance matrix grows quadratically in both and the s, and so, storage becomes impractical for large numbers of level 2 and level 3 groups. However, only certain sub-blocks are required for the addition of pointwise confidence intervals to curve estimates. In particular, we only require the non-zero sub-blocks of the general three-level sparse matrix given in Section 3 of Nolan & Wand (2018) that correspond to . In the case of the three-level Gaussian response linear model, Nolan & Wand’s
[TABLE]
As described in Nolan, Menictas & Wand (2019), the SolveThreeLevelSparseLeastSquares algorithm arises in the special case where is the minimizer of the least squares problem given in equation (9), where has the three-level sparse form and is partitioned according to that shown in equation (7) of Nolan & Wand (2018). This algorithm can be used for fitting three-level group-specific curve models by making use of Result 3.
Result 3**.**
Computation of and each of the sub-blocks of listed in (8) are expressible as the three-level sparse matrix least squares form:
[TABLE]
where the non-zero sub-blocks and , according to the notation in Section 3.1 of Nolan & Wand (2018), are for and :
[TABLE]
[TABLE]
with each of these matrices having {\tilde{o}}_{ij}=o_{ij}+K_{\mbox{\rm\tiny gbl}}+2+K_{\mbox{\rm\tiny grp}}^{\mbox{\rm\tinyg}}+2+K_{\mbox{\rm\tiny grp}}^{\mbox{\rm\tinyh}} rows and with having columns, \overset{\ \mbox{{\bullet}}}{\boldsymbol{B}}_{i} having q_{1}=2+K_{\mbox{\rm\tiny grp}}^{\mbox{\rm\tinyg}} columns and \overset{\mbox{{\bullet}}\mbox{{\bullet}}}{\boldsymbol{B}}_{ij} having q_{2}=2+K_{\mbox{\rm\tiny grp}}^{\mbox{\rm\tinyh}} columns. The solutions are
[TABLE]
[TABLE]
[TABLE]
[TABLE]
[TABLE]
and
[TABLE]
A derivation of Result 3 is given in Section S.7 of the web-supplement. Result 3 combined with Theorem 4 of Nolan & Wand (2018) leads to Algorithm 3. The SolveThreeLevelSparseLeastSquares algorithm is given in Section S.12.
3.2 Mean Field Variational Bayes
A Bayesian extension of (20) and (22) is:
[TABLE]
The following mean field restriction is imposed on the joint posterior density function of all parameters in (23):
[TABLE]
The optimal -density functions for the parameters of interest are
[TABLE]
The optimal -density parameters are determined through an iterative coordinate ascent algorithm, details of which are given in Section S.10 of the web-supplement. As in the two-level case, the updates for and may be written in the same form as (13) but with a three-level version of the matrix and
[TABLE]
For large numbers of level 2 and level 3 groups, ’s size becomes infeasible to deal with. However, only relatively small sub-blocks of are needed for variational inference regarding the variance and covariance parameters. These sub-block positions correspond to the non-zero sub-block positions of a general three-level sparse matrix defined in Section 3 of Nolan & Wand (2018). Here, Nolan & Wand’s
[TABLE]
We appeal to Result 4 for a streamlined mean field variational Bayes algorithm.
Result 4**.**
The mean field variational Bayes updates of and each of the sub-blocks of in (26) are expressible as a three-level sparse matrix least squares problem of the form:
[TABLE]
where the non-zero sub-blocks and , according to the notation in Section 3.1 of Nolan & Wand (2018), are for and .
[TABLE]
[TABLE]
with each of these matrices having {\tilde{o}}_{ij}=o_{ij}+2+K_{\mbox{\rm\tiny gbl}}+2+K_{\mbox{\rm\tiny grp}}^{\mbox{\rm\tinyg}}+2+K_{\mbox{\rm\tiny grp}}^{\mbox{\rm\tinyh}} rows and with having columns, \overset{\ \mbox{{\bullet}}}{\boldsymbol{B}}_{i} having q_{1}=2+K_{\mbox{\rm\tiny grp}}^{\mbox{\rm\tinyg}} columns and \overset{\mbox{{\bullet}}\mbox{{\bullet}}}{\boldsymbol{B}}_{ij} having q_{2}=2+K_{\mbox{\rm\tiny grp}}^{\mbox{\rm\tinyh}} columns. The solutions are
[TABLE]
[TABLE]
[TABLE]
[TABLE]
[TABLE]
and
[TABLE]
Algorithm 4 makes use of Result 4 to facilitate streamlined computation of all variational parameters in the three-level group specific curves model.
Figure 3 provides illustration of Algorithm 4 by showing the fits to the Figure 1 ultrasound data. Posterior mean curves and (narrow) 99% pointwise credible intervals are shown. As discussed in the next section, such fits can be obtained rapidly and accurately and Algorithm 4 is scalable to much larger data sets of the type illustrated by Figures 1 and 3.
4 Accuracy and Speed Assessment
In this section we provide some assessment of the accuracy and speed of the inference delivered by streamlined variational inference for group-specific curves models.
4.1 Accuracy Assessment
Mean field restrictions such as (12) and (24) imply that there is some loss of accuracy in inference produced by Algorithms 2 and 4. However, at least for the Gaussian response case treated here, approximate parameter orthogonality between the coefficient parameters and covariance parameters from likelihood theory implies that such restrictions are mild and mean field accuracy is high. Figure 4 corroborates this claim by assessing accuracy of the mean function estimates and 95% credible intervals at the median values of frequency for each panel in Figure 3. As a benchmark we use Markov chain Monte Carlo-based inference via the rstan package (Guo et al., 2018). After a warmup of size we retained Markov chain Monte Carlo samples from the mean function and median frequency posterior distributions and used kernel density estimation to approximate the corresponding posterior density function. For a generic univariate parameter , the accuracy of an approximation to is defined to be
[TABLE]
The percentages in the top right-hand panel of Figure 4 correspond to (27) with replacement of by the above-mentioned kernel density estimate. In this case accuracy is seen to be excellent, with accuracy percentages between 97% and 99% for all 40 curves.
4.2 Speed Assessment
We also conducted some simulation studies to assess the speed of streamlined variational higher level group-specific curve models, in terms of both comparative advantage over naïve implementation and absolute performance. The focus of these studies was variational inference in the two-level case and to probe maximal speed potential Algorithm 2 was implemented in the low-level computer language Fortran 77. An implementation of the naïve counterpart of Algorithm 2, involving storage and direct calculations concerning the full matrix, was also carried out. We then simulated data according to model (1) with ,
[TABLE]
where, for each , , and are, respectively, random draws from the distribution and the sets and . The level-2 sample sizes generated randomly from the set and the level-1 sample sizes ranging over the set . All data were generated from a Uniform distribution over the unit interval. Table 1 summarizes the timings based on 100 replications with the number of mean field variational Bayes iterations fixed at 50. The study was run on a MacBook Air laptop with a 2.2 gigahertz processor and 8 gigabytes of random access memory.
For ranging from to we see that the naïve to streamlined ratios increase from about to . When the naïve implementation fails to run due to its excessive storage demands. In contrast, the streamlined fits are produced in about seconds. It is clear that streamlined variational inference is to be preferred and is the only option for large numbers of groups.
We then obtained timings for the streamlined algorithm for becoming much larger, taking on values , , and . The iterations in Algorithm 2 were stopped when the relative increase in the marginal log-likelihood fell below . The average and standard deviation times in seconds over 100 replications are shown in Table 2. We see that the computational times are approximately linear in . Even with twelve and a half thousand groups, Algorithm 2 is able to deliver fitting and inference on a contemporary laptop computer in about one and a half minutes.
Acknowledgments
This research was partially supported by the Australian Research Council Discovery Project DP140100441. The ultrasound data was provided by the Bioacoustics Research Laboratory, Department of Electrical and Computer Engineering, University of Illinois at Urbana-Champaign, Illinois, U.S.A.
References
Atay-Kayis, A. & Massam, H. (2005). A Monte Carlo method for computing marginal likelihood in nondecomposable Gaussian graphical models. Biometrika, 92, 317–335.
Bates, D., Mächler, M., Bolker, B. and Walker, S. (2015). Fitting linear mixed-effects models using lme4. Journal of Statistical Software, 67(1), 1–48.
Bishop, C.M. (2006). Pattern Recognition and Machine Learning. New York: Springer.
Brumback, B.A. and Rice, J.A. (1998). Smoothing spline models for the analysis of nested and crossed samples of curves (with discussion). Journal of the American Statistical Association, 93, 961–994.
Donnelly, C.A., Laird, N.M. and Ware, J.H. (1995). Prediction and creation of smooth curves for temporally correlated longitudinal data. Journal of the American Statistical Association, 90, 984–989.
Durban, M., Harezlak, J., Wand, M.P. and Carroll, R.J. (2005). Simple fitting of subject-specific curves for longitudinal data. Statistics in Medicine, 24, 1153–1167.
Goldsmith, J., Zipunnikov, V. and Schrack, J. (2015). Generalized multilevel function-on-scalar regression and principal component analysis. Biometrics, 71, 344–353.
Guo, J., Gabry, J. and Goodrich, B. (2018). rstan: R interface to Stan. R package version 2.18.2.
Huang, A. and Wand, M.P. (2013). Simple marginally noninformative prior distributions for covariance matrices. Bayesian Analysis, 8, 439–452.
Lee, C.Y.Y. and Wand, M.P. (2016). Variational inference for fitting complex Bayesian mixed effects models to health data. Statistics in Medicine, 35, 165–188.
Nolan, T.H., Menictas, M. and Wand, M.P. (2019). Streamlined computing for variational inference with higher level random effects. Unpublished manuscript submitted to the arXiv.org e-Print archive; on hold as of 11th March 2019. Soon to be posted also on http://matt-wand.utsacademics.info/statsPapers.html
Nolan, T.H. and Wand, M.P. (2018). Solutions to sparse multilevel matrix problems. Unpublished manuscript available at https://arxiv.org/abs/1903.03089.
Pinheiro, J.C. and Bates, D.M. (2000). Mixed-Effects Models in S and S-PLUS. New York: Springer.
Pinheiro, J., Bates, D., DebRoy, S., Sarkar, D. and R Core Team. (2018). nlme: Linear and nonlinear mixed effects models. R package version 3.1.
http://cran.r-project.org/package=nlme.
Pratt, J.H., Jones, J.J., Miller, J.Z., Wagner, M.A. and Fineberg, N.S. (1989). Racial differences in aldosterone excretion and plasma aldosterone concentrations in children. New England Journal of Medicine, 321, 1152–1157.
Robinson, G.K. (1991). That BLUP is a good thing: the estimation of random effects. Statistical Science, 6, 15–51.
Trail, J.B., Collins, L.M., Rivera, D.E., Li, R., Piper, M.E. and Baker, T.B. (2014). Functional data analysis for dynamical system identification of behavioral processes. Psychological Methods, 19(2), 175–187.
Verbyla, A.P., Cullis, B.R., Kenward, M.G. and Welham, S.J. (1999). The analysis of designed experiments and longitudinal data by using smoothing splines (with discussion). Applied Statistics, 48, 269–312.
Wahba, G. (1990). Spline Models for Observational Data. Philadelphia: Society for Industrial and Applied Mathematics.
Wand, M.P. and Ormerod, J.T. (2008). On semiparametric regression with O’Sullivan penalized splines. Australian and New Zealand Journal of Statistics, 50, 179–198.
Wand, M.P. and Ormerod, J.T. (2011). Penalized wavelets: embedding wavelets into semiparametric regression. Electronic Journal of Statistics, 5, 1654–1717.
Wang, Y. (1998). Mixed effects smoothing spline analysis of variance. Journal of the Royal Statistical Society, Series B, 60, 159–174.
Wirtzfeld, L.A., Ghoshal, G., Rosado-Mendez, I.M., Nam, K., Park, Y., Pawlicki, A.D., Miller, R.J., Simpson, D.G., Zagzebski, J.A., Oelze, M.I., Hall, T.J. and O’Brien, W.D. (2015). Quantitative ultrasound comparison of MAT and 4T1 mammary tumors in mice and rates across multiple imaging systems. Journal of Ultrasound Medicine, 34, 1373–1383.
Zhang, D., Lin, X., Raz, J. and Sowers, M. (1998). Semi-parametric stochastic mixed models for longitudinal data. Journal of the American Statistical Association, 93, 710–719.
Web-Supplement for:
Streamlined Variational Inference for
Higher Level Group-Specific Curve Models
By M. Menictas1, T.H. Nolan1, D.G. Simpson2 and M.P. Wand1
University of Technology Sydney1 and University of Illinois2
S.1 Derivation of Result 1
Straightforward algebra can be used to verify that
[TABLE]
where and have sparse forms (10) with non-zero sub-blocks equal to
[TABLE]
Therefore, in view of (7) and (8),
[TABLE]
S.2 Derivation of Algorithm 1
Algorithm 1 is simply a proceduralization of Result 1.
S.3 The Inverse G-Wishart and Inverse Distributions
The Inverse G-Wishart corresponds to the matrix inverses of random matrices that have a G-Wishart distribution (e.g. Atay-Kayis & Massam, 2005). For any positive integer , let be an undirected graph with nodes labeled and set consisting of sets of pairs of nodes that are connected by an edge. We say that the symmetric matrix respects if
[TABLE]
A random matrix has an Inverse G-Wishart distribution with graph and parameters and symmetric matrix , written
[TABLE]
if and only if the density function of satisfies
[TABLE]
over arguments such that is symmetric and positive definite and respects . Two important special cases are
[TABLE]
for which the Inverse G-Wishart distribution coincides with the ordinary Inverse Wishart distribution, and
[TABLE]
for which the Inverse G-Wishart distribution coincides with a product of independent Inverse Chi-Squared random variables. The subscripts of and reflect the fact that is a full matrix and is a diagonal matrix in each special case.
The case corresponds to the ordinary Inverse Wishart distribution. However, with message passing in mind, we will work with the more general Inverse G-Wishart family throughout this article.
In the special case the graph and the Inverse G-Wishart distribution reduces to the Inverse Chi-Squared distributions. We write
[TABLE]
for this special case with and scalar.
S.4 Derivation of Result 2
It is straightforward to verify that the and updates, given at (13), may be written as
[TABLE]
where and have the forms (10) with
[TABLE]
and
[TABLE]
Result 2 immediately follows from Theorem 2 of Nolan & Wand (2018).
S.5 Derivation of Algorithm 2
We provide expressions for the -densities for mean field variational Bayesian inference for the parameters in (11), with product density restriction (12). Arguments analogous to those given in, for example, Appendix C of Wand & Ormerod (2011) lead to:
[TABLE]
where
[TABLE]
with , and defined via (14),
[TABLE]
where and
[TABLE]
where \boldsymbol{C}_{\mbox{\rm\tiny gbl,i}}\equiv[\boldsymbol{X}_{i}\,\boldsymbol{Z}_{\mbox{\rm\tiny gbl,i}}], \boldsymbol{C}_{\mbox{\rm\tiny grp,i}}\equiv[\boldsymbol{X}_{i}\,\boldsymbol{Z}_{\mbox{\rm\tiny grp,i}}], and with reciprocal moment
[TABLE]
where and
[TABLE]
with reciprocal moment
[TABLE]
where and
[TABLE]
with reciprocal moment
[TABLE]
where \xi_{\mathfrak{q}(\boldsymbol{\Sigma})}=\nu_{\mbox{\rm\tiny{\boldsymbol{\Sigma}}}}+2+m
[TABLE]
with inverse moment ,
[TABLE]
where ,
[TABLE]
with reciprocal moment ,
[TABLE]
where ,
[TABLE]
with reciprocal moment ,
[TABLE]
where ,
[TABLE]
with reciprocal moment and
[TABLE]
where \xi_{\mathfrak{q}(\boldsymbol{A}_{\mbox{\rm\tiny\boldsymbol{\Sigma}}})}=\nu_{\mbox{\rm\tiny{\boldsymbol{\Sigma}}}}+2,
[TABLE]
with inverse moment \boldsymbol{M}_{\mathfrak{q}(\boldsymbol{A}_{\mbox{\rm\tiny\boldsymbol{\Sigma}}}^{-1})}=\xi_{\mathfrak{q}(\boldsymbol{A}_{\mbox{\rm\tiny\boldsymbol{\Sigma}}})}\boldsymbol{\Lambda}_{\mathfrak{q}(\boldsymbol{A}_{\mbox{\rm\tiny\boldsymbol{\Sigma}}})}^{-1}.
S.6 Marginal Log-Likelihood Lower Bound and Derivation
The expression for the lower bound on the marginal log-likelihood for Algorithm 2 is
[TABLE]
Derivation: The lower-bound on the marginal log-likelihood is achieved through the following expression:
[TABLE]
[TABLE]
First we note that
[TABLE]
where
[TABLE]
and
[TABLE]
Therefore,
[TABLE]
The remainder of the expectations in (S.8) are expressed as:
[TABLE]
[TABLE]
[TABLE]
[TABLE]
[TABLE]
[TABLE]
[TABLE]
[TABLE]
[TABLE]
[TABLE]
[TABLE]
[TABLE]
[TABLE]
[TABLE]
[TABLE]
[TABLE]
[TABLE]
[TABLE]
In the summation of each of these terms, note that the coefficient of is
[TABLE]
The coefficient of is
[TABLE]
The coefficient of is
[TABLE]
The coefficient of is
[TABLE]
The coefficient of is
[TABLE]
The coefficient of is
[TABLE]
The coefficient of is
[TABLE]
The coefficient of E_{\mathfrak{q}}\{\log|\boldsymbol{A}_{\mbox{\rm\tiny\boldsymbol{\Sigma}}}|\} is
[TABLE]
Therefore these terms can be dropped and then the cancellations led by the above expectations leads to the lower bound expression in (S.7).
S.7 Derivation of Result 3
If and have the same forms given by equation (7) in Nolan & Wand (2018) with
[TABLE]
[TABLE]
then straightforward algebra leads to
[TABLE]
where
[TABLE]
and as defined in (22). The remainder of the derivation of Result 3 is analogous to that of Result 1.
S.8 Derivation of Algorithm 3
Algorithm 3 is simply a proceduralization of Result 3.
S.9 Derivation of Result 4
It is straightforward to verify that the and updates, given at (13) but with as given in (25), may be written as
[TABLE]
where and have the forms given by equation (7) in Nolan & Wand (2018) with
[TABLE]
[TABLE]
Result 4 immediately follows from Theorem 4 of Nolan & Wand (2018).
S.10 Derivation of Algorithm 4
We provide expressions for the -densities for mean field variational Bayesian inference for the parameters in (23) with product density restriction (24).
[TABLE]
where
[TABLE]
with , \boldsymbol{o}_{\mbox{{\bf MFVB}}}\equiv\left[\begin{array}[]{c}\boldsymbol{\Sigma}_{\boldsymbol{\beta}}^{-1}\boldsymbol{\mu}_{\boldsymbol{\beta}}\\[4.30554pt] \boldsymbol{0}\end{array}\right] and as given in (25).
[TABLE]
where and
[TABLE]
[TABLE]
[TABLE]
where \boldsymbol{C}_{\mbox{\rm\tiny gbl,ij}}\equiv[\boldsymbol{X}_{ij}\,\boldsymbol{Z}_{\mbox{\rm\tiny gbl,ij}}], \boldsymbol{C}^{\mbox{\rm\tinyg}}_{\mbox{\rm\tiny grp,ij}}\equiv[\boldsymbol{X}_{ij}\,\boldsymbol{Z}^{\mbox{\rm\tinyg}}_{\mbox{\rm\tiny grp,ij}}], \boldsymbol{C}^{\mbox{\rm\tinyh}}_{\mbox{\rm\tiny grp,ij}}\equiv[\boldsymbol{X}_{ij}\,\boldsymbol{Z}^{\mbox{\rm\tinyh}}_{\mbox{\rm\tiny grp,ij}}] and with reciprocal moment
[TABLE]
where and
[TABLE]
with reciprocal moment
[TABLE]
where \xi_{\mathfrak{q}(\sigma_{\mbox{\rm\tiny grp,g}}^{2})}=\nu_{\mbox{\rm\tiny grp, g}}+mK_{\mbox{\rm\tiny grp}}^{\mbox{\rm\tinyg}} and
[TABLE]
with reciprocal moment \mu_{\mathfrak{q}(1/{\sigma_{\mbox{\rm\tiny grp,g}}^{2}})}=\xi_{\mathfrak{q}(\sigma_{\mbox{\rm\tiny grp,g}}^{2})}/\lambda_{\mathfrak{q}(\sigma_{\mbox{\rm\tiny grp,g}}^{2})},
[TABLE]
where \xi_{\mathfrak{q}(\boldsymbol{\Sigma}_{\mbox{\rm\tinyg}})}=\nu_{\mbox{\rm\tiny{\boldsymbol{\Sigma}_{\mbox{\rm\tiny}}}}}+2+m and
[TABLE]
with inverse moment \boldsymbol{M}_{\mathfrak{q}(\boldsymbol{\Sigma}_{\mbox{\rm\tinyg}}^{-1})}=(\xi_{\mathfrak{q}(\boldsymbol{\Sigma}_{\mbox{\rm\tinyg}})}-1)\boldsymbol{\Lambda}_{\mathfrak{q}(\boldsymbol{\Sigma}_{\mbox{\rm\tinyg}})}^{-1},
[TABLE]
where \xi_{\mathfrak{q}(\sigma_{\mbox{\rm\tiny grp,h}}^{2})}=\nu_{\mbox{\rm\tiny grp, h}}+K_{\mbox{\rm\tiny grp}}^{\mbox{\rm\tinyh}}\sum_{i=1}^{m}n_{i} and
[TABLE]
with reciprocal moment \mu_{\mathfrak{q}(1/{\sigma_{\mbox{\rm\tiny grp,h}}^{2}})}=\xi_{\mathfrak{q}(\sigma_{\mbox{\rm\tiny grp,h}}^{2})}/\lambda_{\mathfrak{q}(\sigma_{\mbox{\rm\tiny grp,h}}^{2})},
[TABLE]
where \xi_{\mathfrak{q}(\boldsymbol{\Sigma}_{\mbox{\rm\tinyh}})}=\nu_{\mbox{\rm\tiny{\boldsymbol{\Sigma}_{\mbox{\rm\tiny}}}}}+2+\sum_{i=1}^{m}n_{i} and
[TABLE]
with inverse moment \boldsymbol{M}_{\mathfrak{q}(\boldsymbol{\Sigma}_{\mbox{\rm\tinyh}}^{-1})}=(\xi_{\mathfrak{q}(\boldsymbol{\Sigma}_{\mbox{\rm\tinyh}})}-1)\boldsymbol{\Lambda}_{\mathfrak{q}(\boldsymbol{\Sigma}_{\mbox{\rm\tinyh}})}^{-1},
[TABLE]
where ,
[TABLE]
with reciprocal moment ,
[TABLE]
where ,
[TABLE]
with reciprocal moment ,
[TABLE]
where \xi_{\mathfrak{q}(a_{\mbox{\rm\tiny grp, g}})}=\nu_{\mbox{\rm\tiny grp, g}}+1,
[TABLE]
with reciprocal moment \mu_{\mathfrak{q}(1/a_{\mbox{\rm\tiny grp, g}})}=\xi_{\mathfrak{q}(a_{\mbox{\rm\tiny grp, g}})}/\lambda_{\mathfrak{q}(a_{\mbox{\rm\tiny grp, g}})} and
[TABLE]
where \xi_{\mathfrak{q}(\boldsymbol{A}_{\mbox{\rm\tiny\boldsymbol{\Sigma}{\mbox{\rm\tiny}}}})}=\nu_{\mbox{\rm\tiny{\boldsymbol{\Sigma}{\mbox{\rm\tiny}}}}}+2,
[TABLE]
with inverse moment \boldsymbol{M}_{\mathfrak{q}(\boldsymbol{A}_{\mbox{\rm\tiny\boldsymbol{\Sigma}{\mbox{\rm\tiny}}}}^{-1})}=\xi_{\mathfrak{q}(\boldsymbol{A}_{\mbox{\rm\tiny\boldsymbol{\Sigma}{\mbox{\rm\tiny}}}})}\boldsymbol{\Lambda}_{\mathfrak{q}(\boldsymbol{A}_{\mbox{\rm\tiny\boldsymbol{\Sigma}_{\mbox{\rm\tiny}}}})}^{-1},
[TABLE]
where \xi_{\mathfrak{q}(a_{\mbox{\rm\tiny grp, h}})}=\nu_{\mbox{\rm\tiny grp, h}}+1,
[TABLE]
with reciprocal moment \mu_{\mathfrak{q}(1/a_{\mbox{\rm\tiny grp, h}})}=\xi_{\mathfrak{q}(a_{\mbox{\rm\tiny grp, h}})}/\lambda_{\mathfrak{q}(a_{\mbox{\rm\tiny grp, h}})} and
[TABLE]
where \xi_{\mathfrak{q}(\boldsymbol{A}_{\mbox{\rm\tiny\boldsymbol{\Sigma}{\mbox{\rm\tiny}}}})}=\nu_{\mbox{\rm\tiny{\boldsymbol{\Sigma}{\mbox{\rm\tiny}}}}}+2
[TABLE]
with inverse moment \boldsymbol{M}_{\mathfrak{q}(\boldsymbol{A}_{\mbox{\rm\tiny\boldsymbol{\Sigma}{\mbox{\rm\tiny}}}}^{-1})}=\xi_{\mathfrak{q}(\boldsymbol{A}_{\mbox{\rm\tiny\boldsymbol{\Sigma}{\mbox{\rm\tiny}}}})}\boldsymbol{\Lambda}_{\mathfrak{q}(\boldsymbol{A}_{\mbox{\rm\tiny\boldsymbol{\Sigma}_{\mbox{\rm\tiny}}}})}^{-1}.
S.11 The SolveTwoLevelSparseLeastSquares Algorithm
The SolveTwoLevelSparseLeastSquares is listed in Nolan et al. (2018) and based on Theorem 2 of Nolan & Wand (2018). Given its centrality to Algorithms 1 and 2 we list it again here. The algorithm solves a sparse version of the the least squares problem:
[TABLE]
which has solution where where and have the following structure:
[TABLE]
The sub-matrices corresponding to the non-zero blocks of are labelled according to:
[TABLE]
with denoting sub-blocks that are not of interest. The SolveTwoLevelSparseLeastSquares algorithm is given in Algorithm S.1.
S.12 The SolveThreeLevelSparseLeastSquares Algorithm
The SolveThreeLevelSparseLeastSquares, listed in Nolan et al. (2018) is a proceduralization of Theorem 4 of Nolan & Wand (2018). Since it is central to Algorithms 3 and 4 we list it here. The SolveThreeLevelSparseLeastSquares algorithm is concerned with solving the sparse three-level version of
[TABLE]
with the solution where where and have the following structure:
[TABLE]
and
[TABLE]
The three-level sparse matrix inverse problem involves determination of the sub-blocks of corresponding to the non-zero sub-blocks of . Our notation for these sub-blocks is illustrated by
[TABLE]
for the , and case. The symbol denotes sub-blocks that are not of interest. The SolveThreeLevelSparseLeastSquares algorithm is given in Algorithm S.2.