Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2864505

Formats

Article sections

Authors

Related links

Bayesian Anal. Author manuscript; available in PMC 2010 May 5.

Published in final edited form as:

Bayesian Anal. 2010; 5(1): 189–212.

doi: 10.1214/10-BA508PMCID: PMC2864505

NIHMSID: NIHMS196539

See other articles in PMC that cite the published article.

This work is motivated by a quantitative Magnetic Resonance Imaging study of the relative change in tumor vascular permeability during the course of radiation therapy. The differences in tumor and healthy brain tissue physiology and pathology constitute a notable feature of the image data—spatial heterogeneity with respect to its contrast uptake profile (a surrogate for permeability) and radiation induced changes in this profile. To account for these spatial aspects of the data, we employ a Gaussian hidden Markov random field (MRF) model. The model incorporates a latent set of discrete labels from the MRF governed by a spatial regularization parameter. We estimate the MRF regularization parameter and treat the number of MRF states as a random variable and estimate it via a reversible jump Markov chain Monte Carlo algorithm. We conduct simulation studies to examine the performance of the model and compare it with a recently proposed method using the Expectation-Maximization (EM) algorithm. Simulation results show that the Bayesian algorithm performs as well, if not slightly better than the EM based algorithm. Results on real data suggest that the tumor “core” vascular permeability increases relative to healthy tissue three weeks after starting radiotherapy, which may be an opportune time to initiate chemotherapy and warrants further investigation.

As a non-invasive visualization tool, quantitative Magnetic Resonance Imaging (qMRI) enables researchers to assess pathological and physiological changes of in vivo tissue that cannot be evaluated with conventional anatomic MRI. Recent applications of qMRI include Cao et al. (2005), Moffat et al. (2005), and Hamstra et al. (2005). These applications share a common goal: to use qMRI to predict (local) therapeutic efficacy early during treatment with the aim of individualizing treatment regimens.

Despite recent advances in cancer treatments, the median survival time of patients with high-grade gliomas (a type of brain cancer) has not substantially increased from approximately one year after diagnosis. This is largely attributable to the tight endothelial junctions in the tumor vascular structure (also known as the blood-tumor barrier, BTB) that blocks large molecules, such as traditional chemotherapeutic agents, from entering these tumors. This mechanism also protects healthy tissue in the brain from the cytotoxicity of chemotherapeutic agents (also known as the blood-brain barrier, BBB). Although it is known that radiation can increase vascular permeability via damage to the vascular structure, standard of care sequential radiotherapy followed by chemotherapy has had limited success in treating brain cancer (Medical Research Council, 2001).

Researchers at the University of Michigan, School of Medicine recently conducted a study aimed at determining if the increase in tumor vascular permeability induced by radiation is transient and whether there is a time at which this increase is maximum (Cao et al. 2005). If this hypothesis holds, it would suggest that chemotherapy should begin during the course of radiation therapy instead of after its completion. This was the first study to use quantitative, high-resolution MRI to assess radiation induced increases in tumor vascular permeability among glioma patients (Cao et al. 2005). Eleven subjects with primary, high-grade gliomas participated in the study. T1-weighted MRI^{1} was performed on each patient, with and without contrast enhancement, prior to the beginning of radiotherapy. The same imaging protocol was subsequently performed after the first and third week of radiotherapy, and at one, three and six months after the completion of radiotherapy. The contrast agent used, Gadolinium diethylenetriaminepentaacetic acid (Gd-DTPA), has approximately the same molecule weight as many chemotherapeutic agents used to treat gliomas. Hence, the contrast uptake (quantified as the difference, on the log scale, of the contrast enhanced and pre-enhancement MRI images) was used as a surrogate of vascular permeability to these drugs (Cao et al. 2005). In this paper, we focus on the change in contrast uptake from baseline to week three (the difference between contrast uptake at week 3 and baseline, e.g. Figures 1a and and2a),2a), which is of special interest to the researchers. For simplicity, we hereafter refer to the change in contrast uptake image as the observed image.

Subject 1 results: (a) observed change in contrast uptake (black to white indicating larger increases in contrast uptake); (b) marginal posterior mean of change in contrast uptake _{zi} and (c) standard deviation _{zi}; thresholded **...**

Subject 2 results: (a) observed change in contrast uptake; (b) marginal posterior mean of change in contrast uptake _{zi} and (c) standard deviation _{zi}; thresholded image assuming (d) spatial independence and (e) spatial **...**

It is known that solid tumors are physiologically different from surrounding healthy tissue, and that the contrast uptake, as well as its change, is heterogeneous. Many qMRI analyses ignore the inherent spatial correlation (at the pixel level) in the data, and treat the observed change in contrast uptake at each pixel as independent observations (e.g. Cao et al. 2005; Moffat et al. 2005; Hamstra et al. 2005), which results in biased variance estimates. To model the change in tumor/brain contrast uptake induced by radiation at the pixel-level, we use a model that accounts for the spatial correlation in the data and respects the distinct boundaries between tumor and healthy tissue. We introduce a layer of discrete hidden labels from a Markov random field (MRF, Besag 1974) which accounts for the spatial correlation in the data and avoids over-smoothing. A priori, the MRF encourages spatial continuity but allows for spatial heterogeneity. Like many similar models, we assume the observed data conditional on the hidden labels are independent and normally distributed (e.g. Lei and Udupa, 2002; Liang and Lauterber, 1999, Ch. 8).

The hidden MRF model we employ, known as the Potts model in statistical physics (Potts, 1952), has been applied in diverse areas such as disease mapping (Green and Richardson, 2002), agriculture (Dryden, Scarr, and Taylor, 2003), and landscape genetics (Francois, Ancelet, and Guillot, 2006; Guillot, Estoup, Mortier, and Cosson, 2005). All the above applications share the same feature as ours—spatial heterogeneity. However, our work builds on previous methods by integrating many of the individual ideas that have been previously published.

Firstly, we estimate the spatial regularization parameter, which is often fixed in the MRF models (e.g. Green and Richardson, 2002). Francois, Ancelet, and Guillot (2006) show that results can be sensitive to the choice of this parameter. However, estimation of the spatial regularization parameter requires a corresponding normalizing constant, which is computationally intractable. Some authors use pseudo-likelihood (Besag, 1974), which avoids the need to estimate the normalizing constant. However, this approach tends to overestimate the regularization parameter and over-smooth the data (Melas and Wilson, 2002). Instead, we use thermodynamic integration proposed by Ogata (1989) to estimate the ratio of these constants. This method was subsequently generalized by Gelman and Meng (1998). In 1997, Potamianos and Goutsias conducted a simulation study comparing several methods that approximate the normalizing constant in the Ising model (a two state Potts model) and recommend that Ogata’s thermodynamic integration method be used.

Secondly, we treat the number of states of the MRF as a model parameter and employ reversible jump Markov chain Monte Carlo (RJMCMC, Green, 1995; Richardson and Green, 1997) on a large scale dataset. The number of states in the MRF is traditionally assumed known (Khalkighi, Soltanian-Zadeh, and Lucas, 2002; Marroquin, Arce, and Botello, 2003). This is reasonable when there is substantial scientific justification, e.g. segmenting the brain into white matter, gray matter and cerebrospinal fluid. However, the segmentation labels lack a biologically meaningful interpretation in our application. Therefore, it is more appropriate to treat the number of labels as a model parameter. In this manuscript, we implement a complete Bayesian approach using RJMCMC and we propose a novel implementation of the reversible jump proposal inspired by the Swendsen-Wang algorithm (Swendsen and Wang, 1987). We focus on the marginal posterior distribution of change in contrast uptake, rather than prediction of the hidden labels. The labels are used to model spatial correlation and have no intrinsic scientific interpretation of interest.

The rest of this manuscript is organized as follows. In the next section we introduce notation and specify the model. We then discuss its implementation in Section 3. In Section 4, we present results from simulation studies and compare the results with the EM algorithm of Zhang et al. (2008). We also investigate the performance under violations to model assumptions and present results from the motivating example. We conclude by summarizing the strengths and weaknesses of our algorithm, and discuss future work.

Following convention, we use *i* = 1, 2, ···, *N* to index pixels (short for picture element). If pixel *i* and *i*′ share a common edge, we call them neighbors, denoted by *i* ~ *i*′. The set of neighbors of pixel *i* is denoted by *i* = {*i*′: *i*′ ~ *i*}. Let **y** = (*y*_{1}, *y*_{2} ···, *y _{N}*)

We assume the observed data are conditionally independent given the hidden labels,

$$[{y}_{i}\mid {z}_{i}=k,{\mu}_{k},{\sigma}_{k}^{2}]\sim \text{N}({\mu}_{k},{\sigma}_{k}^{2}),\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{where}\phantom{\rule{0.16667em}{0ex}}1\le k\le M.$$

The hidden labels follow a Gibbs distribution with joint probability mass function

$$Pr(\mathbf{Z}=\mathbf{z}\mid \beta ,M)={g}^{-1}(\beta ,M)exp\left\{\beta \sum _{i\sim {i}^{\prime}}\text{I}[{z}_{i}={z}_{{i}^{\prime}}]\right\},$$

(1)

where I[·] is the indicator function (I[*z _{i}* =

$$\stackrel{\sim}{Pr}(\mathbf{Z}=\mathbf{z}\mid \beta ,M)=\prod _{i=1}^{N}\frac{exp\{\beta {\sum}_{{i}^{\prime}\in {\partial}_{i}}\text{I}[{z}_{{i}^{\prime}}={z}_{i}]\}}{{\sum}_{k=1}^{M}exp\{\beta {\sum}_{{i}^{\prime}\in {\partial}_{i}}\text{I}[{z}_{{i}^{\prime}}=k]\}},$$

as an approximation to Pr(**Z** = **z** |*β*, *M*). However, Barker and Rayner (1997) show that, under certain circumstances, this pseudo-likelihood may lead to an improper posterior distribution. Others completely avoid estimating it by assuming fixed values for both *β* and *M*, the results of which may depend on the choice of *β* (Francois, Ancelet, and Guillot, 2006). We use thermodynamic integration^{2} (see Appendix 1 for details) as proposed by Ogata (1989), brought to the attention of the statistical community by Gelman and Meng (1998), and subsequently used by Green and Richardson (2002) and Higdon (1998), among others.

It follows that the joint distribution of the data and hidden labels is

$$f(\mathbf{y},\mathbf{Z}=\mathbf{z}\mid \mathit{\mu},{\mathit{\sigma}}^{2},\beta ,M)=Pr(\mathbf{Z}=\mathbf{z}\mid \beta ,M)\prod _{i=1}^{N}{(2\pi {\sigma}_{{z}_{i}}^{2})}^{-1/2}exp[-0.5{({y}_{i}-{\mu}_{{z}_{i}})}^{2}/{\sigma}_{{z}_{i}}^{2}].$$

Although we assume conditional independence of the observed data given the hidden labels, one can show that the observed data are marginally correlated because of the correlation introduced by the MRF labels.

A priori, regularization parameter, *β*, follows a uniform distribution on the interval [0, *β*_{max}] with *β*_{max} = 3. The number of components, *M*, is given a uniform distribution on the set of integers {*M*_{min}, …, *M*_{max}}, where *M*_{min} = 1 and *M*_{max} = 20. The component means, *μ _{k}* are independently and identically distributed (i.i.d.) Unif[

We illustrate the model structure in a Directed Acyclic Graph (DAG) in Figure 3 and write in vector form ** μ** = (

By Bayes’ Theorem, the joint posterior distribution is proportional to

$$\prod _{k=1}^{M}\{{\beta}_{\sigma}^{2.1}{({\sigma}_{k}^{-2})}^{2.1-1}exp(-{\beta}_{\sigma}{\sigma}_{k}^{-2})\}\times \prod _{i=1}^{N}{({\sigma}_{{z}_{i}}^{2})}^{-1/2}exp[-0.5{({y}_{i}-{\mu}_{{z}_{i}})}^{2}/{\sigma}_{{z}_{i}}^{2}]\times {g}^{-1}(\beta ,M)exp\left\{\beta \sum _{i\sim {i}^{\prime}}\text{I}[{z}_{i}={z}_{{i}^{\prime}}]\right\}\times {\beta}_{\sigma}^{2.5-1}exp(-5{\beta}_{\sigma})\times M!.$$

(2)

Note that the above joint distribution is invariant to permutations of the component labels conditional on *M*, and therefore the component means and variances are not identifiable. In many fixed-dimension problems with segmentation as the main goal, a constraint such as *μ*_{1} < *μ*_{2} < ··· < *μ _{k}* is imposed to resolve the identifiability issue. However, when

The conditional posterior distributions of the model parameters are

$$[{\mu}_{k}\mid \xb7]\sim \text{N}\left({N}_{k}^{-1}\sum _{i\in {D}_{k}}{y}_{i},{\sigma}_{k}^{2}/{N}_{k}\right)$$

(3)

$$[{\sigma}_{k}^{2}\mid \xb7]\sim \text{Inv}-\text{Gamma}\left(0.5{N}_{k}+{\alpha}_{\sigma},0.5\sum _{i\in {D}_{k}}{({y}_{i}-{\mu}_{k})}^{2}+{\beta}_{\sigma}\right)$$

(4)

$$[{\beta}_{\sigma}\mid \xb7]\sim \text{Gamma}\left(a+M{\alpha}_{\sigma},b+\sum _{k=1}^{M}{\sigma}_{k}^{-2}\right)$$

(5)

$$\pi (\beta \mid \xb7)\propto {g}^{-1}(\beta ,M)exp\left\{\beta \sum _{i\sim {i}^{\prime}}\text{I}[{z}_{i}={z}_{{i}^{\prime}}]\right\}$$

(6)

$$Pr[\mathbf{Z}=\mathbf{z}\mid \xb7]\propto exp\left\{\beta \sum _{i\sim {i}^{\prime}}\text{I}[{z}_{i}={z}_{{i}^{\prime}}]\right\}\prod _{i=1}^{N}{\sigma}_{{z}_{i}}^{-1}exp\{-0.5{({y}_{i}-{\mu}_{{z}_{i}})}^{2}/{\sigma}_{{z}_{i}}^{2}\}$$

(7)

for *k* = 1, 2, …, *M*, where *D _{k}* = {

Our goal is to establish the underlying change in contrast uptake, *μ _{zi}*, and characterize it with its posterior mean,
${\eta}_{i}={\sum}_{k=1}^{M}{\mu}_{k}Pr({Z}_{i}=k\mid \mathbf{y})$, and variance,
${\psi}_{i}^{2}={\sum}_{k=1}^{M}{({\mu}_{k}-{\eta}_{i})}^{2}Pr({Z}_{i}=k\mid \mathbf{y})$. Both of these quantities are estimated via Markov chain Monte Carlo (MCMC). Let

We initialize the Monte Carlo chain with a relatively large number of components, *M* = 15. The initial values of ** μ** are evenly spaced on the range of observed data ([

The parameters ** μ**,

When updating the spatial regularization parameter at the *t*th iteration, *β*^{(}^{t}^{)}, we use a Gaussian proposal distribution centered at *β*^{(}^{t}^{−1)},
${\beta}^{\ast}\sim \text{N}({\beta}^{(t-1)},{\sigma}_{\text{pro}}^{2})$, where the variance of the proposal,
${\sigma}_{\text{pro}}^{2}$, is adaptively tuned during burn-in to give an acceptance rate of approximately 35%. Since we use a symmetric proposal distribution, the acceptance probability only depends on the ratio of the conditional posterior distributions in Equation (6). The ratios of normalizing constants are estimated before-hand on a grid of values for *β* and *M* (*β* = 0, 0.05, ···, 3.00 and *M* = 1, 2, ···, 20, see Appendix 1 for details), interpolating between values when necessary.

The kernel of the conditional posterior distribution of **Z**, Equation (7), consists of two parts—an interaction term (between neighboring pixels) and a likelihood term (sometimes referred to as the external field, Higdon, 1998). A pixel-wise updating scheme may not mix well in the presence of the interaction term (Higdon, 1998). Hence, we use the Swendsen-Wang algorithm, an efficient sampling scheme designed to speed up the mixing of the Potts model (see Appendix 2 for details).

Applications of RJMCMC in Potts models include Green and Richardson (2002), who analyze disease mapping data using a trans-dimensional proposal that emulates the spatial dependence structure of the Potts model; and Dryden, Scarr, and Taylor (2003), who apply a similar scheme to analyze weed and crop images. However, we found that these proposals do not work well for our application: with tens of thousands of pixels per image, the acceptance rate using these proposals is practically zero.

Therefore, we propose a novel implementation of the trans-dimensional move inspired by the Swendsen-Wang algorithm. We first randomly choose between a split and a merge proposal with

$${\text{P}}_{\text{split}}(M)=\{\begin{array}{ll}0\hfill & \text{if}\phantom{\rule{0.16667em}{0ex}}M={M}_{max}\hfill \\ 1\hfill & \text{if}\phantom{\rule{0.16667em}{0ex}}M={M}_{min}\hfill \\ 0.5\hfill & \text{otherwise}\hfill \end{array}\phantom{\rule{0.38889em}{0ex}}\text{and}\phantom{\rule{0.16667em}{0ex}}{\text{P}}_{\text{merge}}(M)=\{\begin{array}{ll}0\hfill & \text{if}\phantom{\rule{0.16667em}{0ex}}M={M}_{min}\hfill \\ 1\hfill & \text{if}\phantom{\rule{0.16667em}{0ex}}M={M}_{max}\hfill \\ 0.5\hfill & \text{otherwise}\hfill \end{array}.$$

If a split proposal is chosen, we randomly pick a component *k* (1 ≤ *k* ≤ *M*) to split, i.e.
${\text{P}}_{\text{select}}^{\text{split}}(k)=1/M$. The dimension of the parameter space increases by two accordingly. To match the increase in dimension, we introduce two independent random variables, *u*_{1}, *u*_{2} ~ Beta(2, 2), and define a bijective transformation
$({\mu}_{{k}_{1}}^{\ast},{\mu}_{{k}_{2}}^{\ast},{{\sigma}_{{k}_{1}}^{\ast}}^{2},{{\sigma}_{{k}_{2}}^{\ast}}^{2})=\psi ({\mu}_{k},{\sigma}_{k}^{2},{u}_{1},{u}_{2})$ that matches the first two moments,

$$\begin{array}{ll}{\mu}_{{k}_{1}}^{\ast}={\mu}_{k}-{u}_{1}{\sigma}_{k},\hfill & {\mu}_{{k}_{2}}^{\ast}={\mu}_{k}+{u}_{1}{\sigma}_{k},\hfill \\ {{\sigma}_{{k}_{1}}^{\ast}}^{2}=2{u}_{2}(1-{u}_{1}^{2}){\sigma}_{k}^{2},\hfill & {{\sigma}_{{k}_{2}}^{\ast}}^{2}=2(1-{u}_{2})(1-{u}_{1}^{2}){\sigma}_{k}^{2}.\hfill \end{array}$$

We denote the proposed set of parameters with a superscript^{*}:
${\mathit{\mu}}^{\ast}={({\mu}_{1}^{\ast},{\mu}_{2}^{\ast},\cdots ,{\mu}_{M+1}^{\ast})}^{\text{T}}$ and
${\mathit{\sigma}}^{\ast 2}={({\sigma}_{1}^{\ast 2},{\sigma}_{2}^{\ast 2},\cdots ,{\sigma}_{M+1}^{\ast 2})}^{\text{T}}$. When there exists a component *h* (1 ≤ *h* ≤ *M* + 1) satisfying
${\mu}_{k1}^{\ast}<{\mu}_{h}^{\ast}<{\mu}_{k2}^{\ast}$, the proposal is not likely to be accepted. We therefore reject such proposals. The common hidden labels define contiguous equivalence classes on the lattice. We base the proposed labeling on the equivalence classes. In contrast to the Swendsen-Wang update, we do not further break contiguous same-labeled regions with stochastic “bonds” (Appendix 2). Rather, we partition the configuration deterministically. This is equivalent to a Swendsen-Wang update with an infinite regularization parameter (*β* = ∞) where all same-labeled neighbors are “bonded”. Suppose there are *L* such equivalence classes *ξ*_{1}, *ξ*_{2}, ···, *ξ _{L}* with current labels

$${p}_{k}^{l}\propto \prod _{i\in {\xi}_{l}}{\sigma}_{k}^{\ast -1}exp\{-0.5{\sigma}_{k}^{\ast -2}{({y}_{i}-{\mu}_{k}^{\ast})}^{2}\},\phantom{\rule{0.38889em}{0ex}}\text{for}\phantom{\rule{0.16667em}{0ex}}k=1,2,\cdots ,M+1,$$

(8)

where the
${p}_{k}^{l}$ are normalized so that
${\sum}_{k=1}^{M+1}{p}_{k}^{l}=1$ for all *l*. Let
${E}_{l}^{\ast}$ denote the proposed label for class *ξ _{l}* from Equation (8), i.e.
${p}_{{E}_{l}^{\ast}}^{l}=Pr({z}_{i}^{\ast}={E}_{l}^{\ast},\forall i\in {\xi}_{l})$. The allocation probability in the current state conditional on the equivalence classes is
${\text{P}}_{\text{alloc}}(\mathbf{z})={\mathrm{\Pi}}_{l=1}^{L}{p}_{{E}_{l}^{\ast}}^{l}$.

We also compute the allocation probability for the proposed configuration given the same equivalence classes,
${\text{P}}_{\text{alloc}}({\mathbf{z}}^{\ast})={\mathrm{\Pi}}_{l=1}^{L}{p}_{{E}_{l}}^{l\ast}$, where
${p}_{k}^{l\ast}\propto {\mathrm{\Pi}}_{i\in {\xi}_{l}}{\sigma}_{k}^{-1}exp\{-0.5{\sigma}_{k}^{-2}{({y}_{i}-{\mu}_{k})}^{2}\}$, for *k* = 1, 2, ···, *M*.

It follows that the acceptance probability for the proposed split move is min(1, *A*) with

$$A=\frac{Pr({\mathbf{z}}^{\ast},{\mathit{\mu}}^{\ast},{\mathit{\sigma}}^{2\ast},{\beta}_{\sigma},\beta ,M+1\mid \mathbf{y})}{Pr(\mathbf{z},\mathit{\mu},{\mathit{\sigma}}^{2},{\beta}_{\sigma},\beta ,M\mid \mathbf{y})}\times \frac{{\text{P}}_{\text{merge}}(M+1){\text{P}}_{\text{select}}^{\text{merge}}({k}_{1},{k}_{2}){\text{P}}_{\text{alloc}}({\mathbf{z}}^{\ast})}{{\text{P}}_{\text{split}}(M){\text{P}}_{\text{select}}^{\text{split}}(k){\text{P}}_{\text{alloc}}(\mathbf{z})b({u}_{1})b({u}_{2})}\times \left|\frac{\partial ({\mu}_{{k}_{1}}^{\ast},{\mu}_{{k}_{2}}^{\ast},{{\sigma}_{{k}_{1}}^{\ast}}^{2},{{\sigma}_{{k}_{2}}^{\ast}}^{2})}{\partial ({\mu}_{k},{\sigma}_{k}^{2},{u}_{1},{u}_{2})}\right|.$$

(9)

Here *b*(·) is the density function of the Beta(2, 2) distribution. The first line is the product of the posterior ratio and the proposal ratio, while the second line is the Jacobian of the bijective transformation.

We first randomly pick a pair of components *k*_{1} and *k*_{2} with adjacent means,
${\text{P}}_{\text{select}}^{\text{merge}}({k}_{1},{k}_{2})=1/(M-1)$. The parameters for the new component are determined as the inverse of the split proposal, i.e.

$${\mu}_{k}^{\ast}=0.5({\mu}_{{k}_{1}}+{\mu}_{{k}_{2}})\phantom{\rule{0.38889em}{0ex}}\text{and}\phantom{\rule{0.38889em}{0ex}}{\mu}_{k}^{\ast 2}+{\sigma}_{k}^{\ast 2}=0.5({\mu}_{{k}_{1}}^{2}+{\sigma}_{{k}_{1}}^{2}+{\mu}_{{k}_{2}}^{2}+{\sigma}_{{k}_{2}}^{2}).$$

We denote the proposed set of parameters by ${\mathit{\mu}}^{\ast}={({\mu}_{1}^{\ast},{\mu}_{2}^{\ast},\cdots ,{\mu}_{M-1}^{\ast})}^{\text{T}}$ and ${\mathit{\sigma}}^{\ast 2}={({\sigma}_{1}^{\ast 2},{\sigma}_{2}^{\ast 2},\cdots ,{\sigma}_{M-1}^{\ast 2})}^{\text{T}}$.

We propose a new configuration **z**^{*} with *M* − 1 components based on the same-labeled contiguous equivalence classes. The new label
${E}_{l}^{\ast}$ of class *ξ _{l}* is drawn from
${p}_{k}^{l}\propto {\mathrm{\Pi}}_{i\in {\xi}_{l}}{\sigma}_{k}^{\ast -1}\times exp\{-0.5{\sigma}_{k}^{\ast -2}{({y}_{i}-{\mu}_{k}^{\ast})}^{2}\}$, for

Hence, the acceptance probability of the merge proposal is min(1, *B*), where

$$B=\frac{Pr({\mathbf{z}}^{\ast},{\mathit{\mu}}^{\ast},{\mathit{\sigma}}^{2\ast},{\beta}_{\sigma},\beta ,M-1\mid \mathbf{y})}{Pr(\mathbf{z},\mathit{\mu},{\mathit{\sigma}}^{2},{\beta}_{\sigma},\beta ,M\mid \mathbf{y})}\times \frac{{\text{P}}_{\text{split}}(M-1){\text{P}}_{\text{select}}^{\text{split}}(k){\text{P}}_{\text{alloc}}({\mathbf{z}}^{\ast})b({u}_{1})b({u}_{2})}{{\text{P}}_{\text{merge}}(M){\text{P}}_{\text{select}}^{\text{merge}}({k}_{1},{k}_{2}){\text{P}}_{\text{alloc}}(\mathbf{z})}\times \left|\frac{\partial ({\mu}_{k}^{\ast},{\sigma}_{k}^{\ast 2},{u}_{1},{u}_{2})}{\partial ({\mu}_{{k}_{1}},{\mu}_{{k}_{2}},{\sigma}_{{k}_{1}}^{2},{\sigma}_{{k}_{2}}^{2})}\right|.$$

We note that identifiability issues may still arise under the constraints *μ*_{1} < *μ*_{2} < ··· < *μ _{M}*, when

To evaluate the performance of the proposed method, we start with simulation studies under the model assumption of conditionally independent noise, and then investigate model performance under violations of this assumption. We define the image Mean Squared Error (MSE) of pixel intensities by
$\text{MSE}({\mu}_{z})={N}^{-1}{\sum}_{i=1}^{N}{({\widehat{\eta}}_{i}-{\mu}_{{z}_{i}}^{\text{true}})}^{2}$, where
${\mu}_{{z}_{i}}^{\text{true}}$ is the simulated true value. We generate multiple datasets and compute the average MSE across the datasets, i.e.
$\overline{\text{MSE}}({\mu}_{z})={L}^{-1}{\sum}_{l=1}^{L}\text{MSE}({\mu}_{z}^{l})$, where
$\text{MSE}({\mu}_{z}^{l})$ is the image MSE of the *l*th simulated dataset (1 ≤ *l* ≤ *L*). The results are compared with an EM algorithm (Zhang et al. 2008) as well as a fixed dimension MCMC. We also present results from the motivating example.

We divide a 128 × 128 lattice into 16 regions of various shapes (Figure 4). We assume the observed pixel values within the same region follow the same distribution, drawn randomly from eight candidate Gaussian distributions with means ** μ** = (−3.5, −2.5, −1.5, −0.5, 0.5, 1.5, 2.5, 3.5)

One set of simulations: the skeleton (top left) the “true” scene (
${\mathit{\mu}}_{z}^{\text{true}}$, top middle), the observed images (**y**, left), the posterior mean (_{i}, middle) and standard deviation (_{i}, right) of pixel **...**

Zhang et al. (2008) report biased parameter estimates when the pixels are treated as independent observations, which is equivalent to a Gaussian mixture model (GMM) with equal component weights. We find similar results with *M* = 8 and *β* = 0: the average MSEs are
${\overline{\text{MSE}}}_{\text{GMM}}^{M=8}({\mu}_{z})=0.333\pm 0.012$ (mean ± standard deviation).

The EM algorithm (Zhang et al. 2008) yields smaller MSEs than the GMM by modeling the spatial correlation, i.e.
${\overline{\text{MSE}}}_{\text{EM}}^{M=8}({\mu}_{z})=0.032\pm 0.024$ (Table 1). However, our results also show that the number of components has a sizable impact on the average MSE. When the number of components is under-specified, e.g. *M* = 6, the performance of the EM algorithm is not satisfactory with
${\overline{\text{MSE}}}_{\text{EM}}^{M=6}({\mu}_{z})=0.107\pm 0.010$. On the other hand, when *M* is over-specified, e.g. *M* = 12, the model fit improves with extra degrees of freedom, i.e.
${\overline{\text{MSE}}}_{\text{EM}}^{M=12}({\mu}_{z})=0.014\pm 0.002$. These results present a potential problem with the EM algorithm—if both small MSE and correct estimation of *M* are important, the EM algorithm does not yield a “best” model under both criteria.

$\overline{\mathit{MSE}}({\mu}_{z})$ of the proposed method vs. fixed dimension MCMC and the EM algorithm.

As an alternative, parallel to the EM algorithm, we also implement a fixed dimension MCMC algorithm (see results in Table 1). As with the EM algorithm, the choice of *M* impacts the average MSE. When *M* = 8, the average MSE is smaller than that of the EM algorithm, i.e.
${\overline{\text{MSE}}}_{\text{MCMC}}^{M=8}({\mu}_{z})=0.020\pm 0.017$. The same holds true when *M* = 6 with
${\overline{\text{MSE}}}_{\text{MCMC}}^{M=6}({\mu}_{z})=0.103\pm 0.005$. When the data are over-fitted with *M* = 12, the average MSE decreases to
${\overline{\text{MSE}}}_{\text{MCMC}}^{M=12}({\mu}_{z})=0.013\pm 0.002$. In summary, the MCMC algorithm has smaller average MSE than the EM under all three choices for *M*. Furthermore, the MCMC algorithm appears to be less variable as evidenced by the smaller standard deviations of the
$\overline{\text{MSE}}$.

We then run our proposed RJMCMC algorithm for 50000 iterations and discard the first half as burn-ins. Figure 4 displays the posterior mean (* _{i}*) and standard deviation (

We reason that the Bayesian algorithms explore the posterior distribution more efficiently and are less prone to getting stuck in local modes, whereas the EM algorithm only guarantees convergence to a local mode. Furthermore, we reason that the additional flexibility of changing between models (*M*) afforded by the RJMCMC algorithm, can further aid in exploring the posterior. We argue that one significant benefit of the proposed algorithm is that it produces satisfactory results (and in many cases better results) without the risk of over- or under-specifying the number of components.

In Section 2.2, we have chosen diffused prior distributions that cover reasonable ranges for most parameters. We now investigate alternative values for *α _{σ}*:

We also examine our proposed algorithm under various signal to noise ratios (SNR). The simulation set in the above section has signal level Δ* _{μ}* =

We evaluate the robustness of the method to violations of the conditional independence assumption by applying Gaussian smoothing kernels on the simulated datasets. We follow the same Gaussian kernel parameters as Zhang et al. (2008), i.e. *σ* = 0.42, 0.85, 1.70, and 3.40. Larger values of *σ* indicate smoother images.

The “edge preservation” of the proposed algorithm is more evident with lighter smoothing (the bottom two rows of Figure 4). The average MSEs for all degrees of smoothing investigated are displayed in Table 1. Within each algorithm, the average MSE increases as smoothing gets heavier. When smoothing is light (*σ* = 0.42), the fixed dimension MCMC achieves smaller MSEs compared to the EM algorithm using the same value for *M*. Our algorithm results in the smallest average MSE regardless of the choice of *M* in the other two algorithms. However, when the smoothing is heavy, e.g. *σ* ≥ 0.85, none of the algorithms has satisfactory performance. Thus we suggest that no smoothing of the data be performed prior to running any of these algorithms.

In the motivating study, eleven subjects received standard of care radiotherapy and MRI images were obtained before, during, and after treatment. All images were registered to anatomical Computed Tomography (CT) images obtained for treatment planning purposes. Due to space limitations, we only display the images on the same two subjects as Zhang et al. (2008).

In the original analysis, Cao et al. (2005) divided the tumor into two regions, one with relatively high pre-treatment contrast uptake and the other with relatively low initial contrast uptake. This more-or-less divides the tumor into a “core” (low initial contrast uptake) surrounded by an “annulus” (with high initial contrast uptake). Tumor cores are typically hypoxic (with low oxygen content) due to a lack of blood supply, which is known to have a protective effect against damage due to both radiation and chemotherapy and may be a source of tumor re-growth. Hence, the focus was on demonstrating that radiation therapy has a transient effect on the tumor “core” with respect to increasing the contrast uptake. The time when this increase is greatest may suggest an optimal time for initiating chemotherapy, allowing for more effective delivery of chemotherapeutic agents. Following Zhang et al. (2008), we use the 95th percentile of healthy tissue contrast uptake at baseline as the threshold to divide the tumor into core and annulus.

We then run the proposed algorithm for 50000 iterations after 25000 burn-in iterations. The acceptance rate of the trans-dimensional proposals is 12.5% for subject 1. The posterior mean of change in contrast uptake and associated standard deviation (Figure 1b and 1c) delineates the heterogeneous response of the tumor. The results from subject 2 are similar (Figure 2b and 2c), as are the results from all subjects (not shown).

As discussed in the introduction, large increases in contrast uptake may suggest more effective delivery of chemotherapy drugs. Hence, an optimal time to deliver chemotherapy may exist when the contrast uptake increase in the tumor is large relative to healthy tissue. One way to quantify this is to define a threshold and compare the proportions of healthy and diseased tissue that exceed this threshold. Although, there is no established threshold, we apply an illustrative threshold, 0.06, used by Zhang et al. (2008). For subject 1, 8.2% of the normal tissue exceeds the threshold, while 59.5% of the tumor “core” has an increase in contrast uptake above the threshold, compared to 15.6% in the tumor “annulus” (Figure 1e). Similarly, for subject 2, 4.0% of healthy tissue exceeds the threshold. While 54.3% of the tumor “core” has a change in contrast uptake that exceeds the threshold, compared to 12.7% in the “annulus” (Figure 2e). The percentages are similar to those reported by Zhang et al. (2008).

Since the choice of the above threshold is arbitrary and is not supported by biological justifications, we propose to use the Mann Whitney U statistic to summarize the differential change in contrast uptake of the tumor and healthy tissue. It can be interpreted as the probability that a randomly selected pixel in the tumor has a larger mean change in contrast uptake than a randomly selected pixel in healthy tissue. A larger U statistic indicates a larger increase in contrast uptake of the tumor relative to healthy tissue which is desirable (i.e. more contrast is entering the tumor). There is a connection between thresholds and the U statistic. The U statistic is equivalent to the empirical area under the receiver operating characteristics curve, or AUC (Section 4.3.4 of Pepe, 2003).

We compute the U statistic at each iteration for the tumor “core” versus healthy tissue and the tumor “annulus” versus healthy tissue. The 95% posterior credible interval of the U statistic between the tumor “core” and healthy tissue is (0.74, 0.76) for subject 1, which suggests the “core” region has a substantial increase, on average, in contrast uptake compared to healthy tissue. The tumor “core” of subject 2 is also well separated from healthy tissue with 95% posterior credible interval (0.81, 0.82). The U statistic between the tumor “annulus” and healthy tissue is under 0.5 for both subjects: (0.23, 0.25) for subject 1 and (0.25, 0.26) for subject 2, suggesting the “annulus” on average has decreased contrast uptake compared to healthy tissue. The U statistic for all eleven subjects are listed in Table 2. Most subjects demonstrate significant increase in contrast uptake in the tumor “core”. We also investigate the sensitivity of the results to alternative cutoff of tumor core vs. annulus, i.e. using 90% and 97.5% percentile of the baseline healthy tissue contrast uptake instead of 95% percentile. The results are listed in Table 2.

In this article, we implement a statistical imaging model that respects the key feature of spatial heterogeneity in the qMRI data. Compared with previous work, we integrate many ideas that have been previously discussed individually. First, we estimate the normalizing constant via thermal integration, instead of using the pseudo-likelihood. Second, we estimate the spatial regularization parameter from the data, rather than holding it fixed. Third, we acknowledge that there is no clear substantive knowledge regarding the number of components, given the lack of biological interpretation of the hidden labels. We therefore marginalize over both *M* and the labels, and focus on the marginal posterior distribution of the pixel change in contrast uptake, which is of primary scientific interest. Furthermore, the split scheme proposed by Richardson and Green (1997) is inadequate for the scale of data, we therefore propose a non-trivial implementation inspired by the Swendsen-Wang algorithm. Finally, we propose using the U statistic to summarize the differential change in contrast uptake of the tumor “core/annulus” versus healthy tissue.

We find the performance of our RJMCMC algorithm improves on the EM algorithm and the fixed dimension MCMC algorithm regardless of the number of components specified in the latter two algorithms. Furthermore, our algorithm is preferable when there is no clear basis for the choice of *M*, since the loss associated with a miss-specified *M* can be severe with both the EM and the fixed dimension MCMC algorithms.

There are alternatives to the proposed algorithm. For example, birth-and-death MCMC (BDMCMC, Stephens 2000) is an alternative trans-dimensional MCMC algorithm. Another alternative could be a spatial infinite mixture model (Guillot, Estoup, Mortier, and Cosson, 2005, Francois, Ancelet, and Taylor, 2006), which extends infinite mixture models using the Dirichlet process prior to the spatial setting. However, implementation of the these alternatives may be challenging giving the large volume of data.

Finally, we would like to address that it is a complicated issue to determine the optimal timing of administering chemotherapy with respect to radiation. The image model discussed in this manuscript provides crucial information, but its role is not to be over-stated. It is more appropriately considered as a piece of a large puzzle. Future research in the image model itself as well as the underlying biological science are much needed. One future direction of our research is in extending the image model by incorporating the repeated images that captures the longitudinal profile of tumor/brain contrast uptake.

This work was partially funded by NIH grants PO1 CA087684-5 and PO1 CA59827-11A1.

The outline of the estimating procedure follows the general method found in Ogata (1989) and Gelman and Meng (1998) with appropriate notational changes. For a given *M* (number of components), we wish to estimate the log ratio of normalizing constants *λ _{M}*(

$$\begin{array}{l}\frac{dln[g(\beta ,M)]}{d\beta}=\frac{dln\left({\sum}_{\mathbf{z}\in {\mathcal{Z}}_{M}}exp\left\{\beta {\sum}_{i\sim j}\text{I}({z}_{i}={z}_{j})\right\}\right)}{d\beta}\\ =\sum _{\mathbf{z}\in {\mathcal{Z}}_{M}}\frac{{\scriptstyle \frac{d}{d\beta}}exp\left\{\beta {\sum}_{i\sim j}\text{I}({z}_{i}={z}_{j})\right\}}{{\sum}_{\mathbf{z}\in {\mathcal{Z}}_{M}}exp\left\{\beta {\sum}_{i\sim j}\text{I}({z}_{i}={z}_{j})\right\}}\times \frac{exp\left\{\beta {\sum}_{i\sim j}\text{I}({z}_{i}={z}_{j})\right\}}{exp\left\{\beta {\sum}_{i\sim j}\text{I}({z}_{i}={z}_{j})\right\}}\\ =\sum _{\mathbf{z}\in {\mathcal{Z}}_{M}}\left\{Pr(\mathbf{Z}=\mathbf{z}\mid \beta ,M)\sum _{i\sim j}\text{I}({z}_{i}={z}_{j})\right\}\\ ={\text{E}}_{\mathbf{Z}}\left\{\sum _{i\sim j}\text{I}({z}_{i}={z}_{j})\right\}.\end{array}$$

Therefore,

$${\lambda}_{M}(a,b)={\int}_{a}^{b}\frac{dln[g(\beta ,M)]}{d\beta}d\beta ={\int}_{a}^{b}{\text{E}}_{\mathbf{Z}}\left\{\sum _{i\sim j}\text{I}({z}_{i}={z}_{j})\right\}d\beta .$$

(10)

In order to estimate *λ _{M}*(

For the RJMCMC proposal, the ratio ln[*g*(*β*, *M*)/*g*(*β*, *M* + 1)] is required. This can be estimated using the fact that ln[*g*(0, *M*)] = ln Σ_{z}_{} 1 = *N* ln(*M*), where *N* is the number of pixels, and the following identity:

$$ln[g(\beta ,M)/g(\beta ,M+1)]={\lambda}_{M}(0,\beta )-{\lambda}_{M+1}(0,\beta )+ln[g(0,M)]+ln[g(0,M+1)].$$

The Swendsen-Wang algorithm stochastically partitions the configuration into same-labeled contiguous regions such that the label for these regions can be updated independently. This clever idea is implemented in the following three steps.

- We first generate independently random variables,
*u*_{ii}_{′}~ Uniform(0, exp{*β*I[*z*=_{i}*z*_{i}_{′}]}), for each pair of neighbors*i*~*i*′. Clearly, the joint distribution of all auxiliary variables**u**= {*u*_{ii}_{′}}_{i}_{~}_{i}_{′}is*f*(**u**) = Π_{i}_{~}_{i}_{′}exp {−*β*I[*z*=_{i}*z*_{i}_{′}]} I[0 ≤*u*_{ii}_{′}≤ exp{*β*I[*z*=_{i}*z*_{i}_{′}]}]. - By Bayes theorem, the posterior distribution of the labels conditional on the auxiliary variables is$$Pr(\mathbf{Z}=\mathbf{z}\mid \mathbf{u},\mathit{\theta})\propto \prod _{i=1}^{N}{\sigma}_{{z}_{i}}^{-1}exp\{-0.5{\sigma}_{{z}_{i}}^{-2}{({y}_{i}-{\mu}_{{z}_{i}})}^{2}\}\prod _{i\sim {i}^{\prime}}\text{I}[0\le {u}_{i{i}^{\prime}}\le exp\{\beta \text{I}[{z}_{i}={z}_{{i}^{\prime}}]\}].$$(11)The second term in Equation (11) defines the range over which the posterior distribution is non-zero. More concretely, when 1 ≤
*u*_{ii}_{′}≤ exp{*β*I[*z*=_{i}*z*_{i}_{′}]} holds, it builds a virtual stochastic “bond” between*i*and*i*′ that requires both pixels to assume the same label. That is if pixel*i*and*i*′ assume the same label in the current configuration, they are “bonded” to assume a common label (could be different from the current one) with probability 1 −*e*^{−}. The bonds define an equivalence relation, and partition a configuration into contiguous regions of bonded pixels. We denote these equivalence classes by^{β}*C*_{1},*C*_{2}, …,*C*. According to the definition, a component consists of one or more such equivalence classes._{J} - According to Equation (11), the labels of the equivalence classes can be updated independently. The new label ${C}_{j}^{(t+1)}$ of class
*j*satisfies$$Pr({C}_{j}^{(t+1)}=k\mid {\mathbf{u}}^{(t)},{\mathit{\theta}}^{(t)})\propto \prod _{i\in {C}_{j}^{(t)}}{\sigma}_{k}^{-1(t)}exp\left\{-.5{\sigma}_{k}^{-2(t)}{\left({y}_{i}-{\mu}_{k}^{(t)}\right)}^{2}\right\}.$$

^{1}A T1-weighted MRI uses short repetition time and echo time in which paramagnetic agents (e.g. Gadolinium) appear brighter in the image. Note that the qMRI images derived from T1-weighted MRI reflect the change in the intensities of the MRI images, and therefore reflect the change in the density of the paramagnetic agents, e.g. contrast uptake of Gd-DTPA.

^{2}Technically, we only need and compute the ratio of normalizing constants in this manuscript, as detailed in Appendix 1. However, the normalizing constant itself can be readily approximated from the ratios (Gelmen and Meng, 1998).

- Medical Research Council. Randomized Trial of Procarbazine, Lomustine, and Vincristine in the Adjuvant Treatment of High-Grade Astrocytoma: A Medical Research Council Trial. J Clin Oncol. 2001;19(2):509–518. [PubMed]
- Barker SA, Rayner PJW. Unsupervised Image Segmentation Using Markov Random Field Models. 1997.
- Besag J. Spatial Interaction and the Statistical Analysis of Lattice Systems (with discussions) Journal of the Royal Statistical Society Series B. 1974;36:192–236.
- Cao Y, Tsien CI, Shen Z, Tatro DS, Ten Haken R, Kessler ML, Chenevert TL, Lawrence TS. Use of Magnetic Resonance Imaging to Assess Blood-Brain/Blood-Glioma Barrier Opening During Conformal Radiotherapy. Journal of Clinical Oncology. 2005;23(18):4127–4136. [PubMed]
- Cappé O, Robert CP, Rydén T. Reversible Jump, Birth-and-Death and More General Continuous Time Markov Chain Monte Carlo Samplers. Journal of the Royal Statistical Society B. 2003;65:679–700.
- Chalmond B. An Iterative Gibbsian Technique for Reconstruction of Mary Images. Pattern Recognition. 1989;22:747–761.
- Dryden IL, Scarr MR, Taylor CC. Bayesian Texture Segmentation of Weed and Crop Images Using Reversible Jump Markov Chain Monte Carlo Methods. Journal of Royal Statistical Society, Series C. 2003;52:31–50.
- Francois O, Ancelet S, Guillot G. Bayesian Clustering Using Hidden Markov Random Fields in Spatial Population Genetics. Genetics. 2006;174:805–816. [PubMed]
- Gelman A, Meng X. Simulating Normalizing Constants: from IMportance Sampling to Bridge Sampling to Path Sampling. Statistical Science. 1998;13(2):163–185.
- Green PJ. Reversible Jump Markov Chain Monte Carlo Computation and Bayesian Model Determination. Biometrika. 1995;82(4):711–732.
- Green PJ, Richardson S. Hidden Markov Models and Disease Mapping. Journal of the American Statistical Association. 2002;97:1055–1070.
- Guillot G, Estoup A, Mortier F, Cosson JF. A Spatial Statistical Model for Landscape Genetics. Genetics. 2005;170(3):1261–80. [PubMed]
- Hamstra DA, Chenevert TL, Moffat BA, Johnson TD, Meyer CR, Mukherji SK, Quint DJ, Gebarski SS, Fan X, Tsien CI, Lawrence TS, Junck L, Rehemtulla A, Ross BD. Evaluation of the functional diffusion map as an early biomarker of time-to-progression and overall survival in high-grade glioma. PNSA. 2005;102(46):16759–16764. [PubMed]
- Higdon D. Auxiliary Variable Methods for Markov Chain Monte Carlo with Applications. Journal of the American Statistical Association. 1998:93.
- Khalighi MM, Soltanian-Zadeh HS, Lucas C. Unsupervised MRI Segmentation with Spatial Connectivity. Proceedings of SPIE, Medical Imaging; 2002. pp. 1742–1750.
- Liang ZP, Lauterbur PC. Principles of Magnetic Resonance Imaging: A Signal Processing Perspective. Section 2 Chapter 8. SPIE press book; 1999.
- Lei T, Udupa JK. Statistical Properties of X-ray CT and MRI from Imaging Physics to Image Statistics. Proceedings of the SPIE, Medical Imaging; 2002. pp. 82–93.
- Marroquin JL, Arce E, Botello S. Hidden Markov Measure Field Models for Image Segmentation. IEEE Transactions on Pattern Analysis and Machine Intelligence. 2003;25:1380–1387.
- Melas DE, Wilson SP. Double Markov Random Fields and Bayesian Image Segmentation. Institute of Electrical and Electronics Engineers Transactions on Signal Processing. 2002;50:357–365.
- Moffat BA, Chenevert TL, Lawrence TS, Meyer CR, Johnson TD, Dong Q, Tsien C, Mukherji S, Quint DJ, Gebarski SS, Robertson PL, Junck LR, Rehemtulla A, Ross BD. Functional diffusion map: A noninvasive MRI biomarker for early stratification of clinical brain tumor response. PNSA. 2005;102(15):5524–5529. [PubMed]
- Ogata Y. A Monte Carlo Method for High Dimensional Integration. Numerische Mathematik. 1989;55:137–157.
- Pepe MS. The Statistical Evaluation of Medical Tests for Classification and Prediction. Chapter 4.3 Oxford Press; 2003.
- Potts RB. Some generalized order-disorder transformations. Proceedings of the Cambridge Philosophical Society. 1952;48:106–109.
- Richardson S, Green PJ. On Bayesian Analysis of Mixtures with an Unknown Number of Components. Journal of the Royal Statistical Society Series B. 1997;59(4):731–792.
- Stephens M. Bayesian Analysis of Mixture Models with an Unknown Number of Components - An Alternative to Reversible Jump Methods. The Annals of Statistics. 2000;28:40–74.
- Swendsen RH, Wang JS. Nonuniversal Critical Dynamics in Monte Carlo Simulations. Physical Review Letters. 1987;58:86–88. [PubMed]
- Won CS, Derin H. Unsupervised Segmentation of Noisy and Textured Images Using Markov Random Fields. Graphical Models and Image Processing. 1992;54:308–328.
- Zhang X, Johnson TD, Little RJA, Cao Y. Quantitative Magnetic Resonance Image Analysis via the EM Algorithm with Stochastic Variation. Annals of Applied Statistics. 2008;2(2):736–755. [PMC free article] [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |