- Email: [email protected]
Bayesian analysis of the stimulus–response curve
Motor Unit Number Estimation (MUNE) and Quantitative EMG (Supplements to Clinical Neurophysiology, Vol. 60) Editors: M.B. Bromberg # 2009 Elsevier B.V. All rights reserved
Chapter 4
Bayesian analysis of the stimulus–response curve P.G. Ridalla,1,*, A.N. Pettitta, P.A. McCombeb and R.D. Hendersonc a
Faculty of Science, School of Mathematical Sciences, Queensland University of Technology, GPO Box 2434, Brisbane, QLD 4001, Australia b University of Queensland, Centre for Clinical Research, Royal Brisbane Hospital, Herston Queensland 4029, Australia
c
Department of Neurology, Royal Brisbane and Women’s Hospital, Herston, QLD 4029, and University of Queensland, Brisbane, Australia 1 Present address: Department of Mathematics and Statistics, Fylde College, Lancaster University, Lancaster LA1 4YF, UK
1. Introduction In this paper we present our Bayesian method for carrying out motor unit number estimation (MUNE) using stimulus–response data collected from surface electrophysiological recordings. We formulate and justify our assumptions in Ridall et al. (2006) and base these on available scientific evidence, as outlined in the previous paper. The object of our methodology is to express the uncertainty about the number of motor units in a muscle in terms of a posterior distribution. From studies taken over time, these posterior distributions can be used to estimate the rate of loss of motor units. A by-product of this method of MUNE is that it provides a means of tracking a given population of motor units over time. Examples of the para*
Correspondence to: Dr. P.G. Ridall (at present address), Department of Mathematics and Statistics, Fylde College, Lancaster University, Lancaster LA1 4YF, UK. Tel.: +44 (0) 1524-592302; Fax: +44 (0) 1524-593429; E-mail: [email protected]
meters that can be tracked over time include the distribution of the excitability parameters, the single motor unit action potentials and the between and within motor unit variability. In Section 2, we describe the fundamentals of Bayesian methodology and terminology and discuss why it is particularly appropriate for the MUNE problem. In Section 3, we outline the three essential components of the model. In Section 4, we explain how our computational method is able to explore the variety of possible models. In Section 5 we show how the output of the Bayesian analysis can be used to measure the rate of loss of motor units in patients with amyotrophic lateral sclerosis (ALS). Finally, in Section 6 we present a brief summary of our conclusions. 2. Bayesian statistics and directed acyclic graphs (DAGs) The attraction of the Bayesian paradigm is that it is a natural and logical way of carrying out statistical reasoning free from the ad hoc thinking that currently permeates traditional statistics.
48 The Bayesian paradigm recognises that for any statistical problem there are many types of uncertainty and that each of these types of uncertainty can be expressed by an appropriate probability distribution. Uncertainty in a parameter, before the data are collected, is expressed using a prior distribution. Once the data have been collected and their variability taken into account through the likelihood, this prior knowledge is updated to obtain what is known as the posterior distribution, the focus of our interest. The primary aim of the Bayesian statistical method of MUNE is to find the posterior probability of the number of motor units in a muscle. However, in a biologically realistic model there are many other unknowns that also need to be estimated. Each unknown must be found by assuming an appropriate prior distribution and then finding its posterior distribution. Sometimes we refer to these parameters as nuisance parameters and we describe what they are and how to deal with them below. A convenient and non-mathematical way of describing the assumptions in our model is using a directed acyclic graph (DAG). It shows the relationship between the unknown stochastic parameters and the data that are assumed to be known. In the graph unknowns are drawn as circular nodes and square nodes are used for known data. We now show how a DAG can be used to illustrate the following statistical concepts: (a) the likelihood, (b) conditional independence, (c) nuisance parameters, and finally (d) Bayesian inference and the derivation of the posterior distribution. These concepts are needed for a satisfactory understanding of how our method of MUNE works. Figure 4.1A shows an example of likelihood. The likelihood, represented by a square node, is the probability of the data (Y) given an unknown parameter (l) and is denoted by P(Y|l). Here l could be a Poisson mean and Y an observation of count data. Likelihood is not solely a Bayesian idea and many non-Bayesians use it, for example to find maximum likelihood estimates. However, Bayesians apply a different type of reasoning. We use panel 4.1D to explain this. Fig. 4.1B represents a plate (the outer rectangle), contain-
ing a group of observations (Yi) that (depicted in this way) are assumed to be conditionally independent, given l. In a DAG, the arrows initially point from the priors (the unknowns) to the data shown in the likelihood. A prior distribution P(l) is able to express the degree of knowledge of a given parameter before the data are collected or taken into account. This degree of knowledge can range from a total lack of knowledge to strong prior knowledge. Figure 4.1C shows the added complication of assuming the existence of a nuisance parameter, S. This nuisance parameter may be of no interest in itself, but its calculation is needed because it provides the link between the unknown and the observed data. The DAG which we use for MUNE has many nuisance parameters. Figure 4.1D shows a re-expression of the model this time with the unknown (the quantity of interest) as the subject. This reversal of arrows (from the direction of the likelihood) is a simple way of showing how Bayesians carry out their statistical reasoning and is a simple and natural idea. It gives the various probability distributions (prior, likelihood) the same mathematical status. The posterior distribution is the probability of the unknown parameter given the data and is shown as P(l|Y). This expresses knowledge about an unknown obtained from the data and is obtained as a result of applying Bayes theorem. This cannot be done without a prior distribution on the parameter of interest. In very simple examples the posterior distribution can be expressed in algebraic terms or the posterior distribution can be assumed to have a known distribution — for example, a normal distribution with mean and variance known in terms of the data. However, in the presence of many nuisance parameters, this often cannot be done and in these cases we must resort to a standard computational technique called Markov chain Monte Carlo (MCMC). MCMC has, over the last 20 years, become an increasingly useful computational tool because of greatly increased computer speeds. MCMC involves an iterative process with the chain taking sequentially changed values. The values of the chain are initialized and then
49
Fig. 4.1 Panel (A) shows a simple likelihood. Panel (B) shows a likelihood consisting of conditionally independent variables. Panel (C) shows the likelihood in the presence of a nuisance parameter (in this case S). Panel (D) shows the posterior distribution with the direction of the arrows in the opposite direction to that of the likelihood.
changed until convergence is attained. Convergence implies that the initial values no longer affect the values produced from the chain. Once convergence is obtained the output of the MCMC is from the posterior distribution. One way of representing the values produced from the posterior is by a histogram of values generated from the posterior (Fig. 4.2). For MUNE, the mode of the posterior gives us the best estimate of motor unit numbers with the spread of the histogram expressing the uncertainty in that estimate. The histogram shown in Fig. 4.2 is taken from 100,000 values generated from the chain. The mode, or the most likely number of motor units, is found to be N = 13 units. The 95% credible interval in this case is
[12, 15]. This means that this interval contains 95% of the posterior distribution and has a different meaning to a classical confidence interval. 3. The elements of our statistical model In the specific problem of MUNE, uncertainty exists at distinct levels and each is related to motor unit physiology. Firstly, there is uncertainty at the level of the axonal threshold which describes whether a motor unit will fire when a stimulus is applied. Secondly, there are two types of variation that contribute to the uncertainty in the compound muscle action potential (CMAP) at the observation level. These are the between motor unit variation (uncertainty at the level of the single motor
50
Fig. 4.2
The posterior probability of the number of motor units, N, expressed as a histogram.
unit response) and the within motor unit variation. Our distribution of CMAPs must be robust enough to describe outlying observations that occur, particularly in the latter stages of ALS. 3.1. The depolarisation of the axon membrane by a stimulus We use the DAG in Fig. 4.3 to illustrate this part of the model. Unit k only fires at time t if, and only if, the stimulus St exceeds the threshold of that motor unit, tk,t, at that instance of time. The variable skt indicates whether a motor unit is on or off at that time. The threshold is not assumed fixed in time but to vary independently as a Gaussian variable with its own mean threshold, mk, and variance parameter, uk. These two parameters are called the excitability parameters. We illustrate the meaning of these two parameters in Fig. 4.4, using the modal output of the Markov chain used to analyse the data shown in the top panel of Fig. 4.4. The first of these, the mean threshold of a motor unit, mk, is defined as the stimulus for which a motor unit has a 50% probability of firing.
The first two mean thresholds of the data shown in the top panel are shown in the second panel of Fig. 4.4. The second set of parameters, the variances of the thresholds, uk, determine the steepness of the corresponding excitability curves and in doing so define the range over which the unit displays probabilistic firing. The region of probabilistic firing (ROPF) is shown in the bottom panel of Fig. 4.4B as a double arrow for the fourth unit. The first three excitability curves on the bottom panel of Fig. 4.4 are not as steep as the rest indicating that the first three variance threshold parameters, u1, u2 and u3, are much larger than those for the rest of the motor units. 3.2. The response at the neuromuscular junction We let the average area of a single motor unit action potential (MUAP) be mk and refer to this as the size of a particular unit. The variance of the distribution of mk is referred to as the between motor unit variation. We assume that each single MUAP is normally distributed about its mean, mk and with a variance s2 common to all units.
51
Fig. 4.3
The DAG for showing the relationship of the relevant variables for describing the depolarisation of an axon.
CMAP area (μVms)
8000 6000 4000 2000 0 45
50
55
60 Stimulus (mA)
65
70
75
1
p(St>τkt)
0.8
m2 m1 R.O.P.F.
0.6
p(St>τkt)=0.5
0.4 0.2 0 45
50
55
60 Stimulus (St)
65
70
75
Fig. 4.4 The top panel shows the stimulus–response curves plotted for the data. The bottom panel shows the excitability curves obtained from the MCMC output.
52 (2) 40
30
30 Frequency
Frequency
(1) 40
20 10 0
20 10
500
1000
1500
0
2000
Single MUAp (μVms) 40
30
30
20 10 0
2000
(4)
40
Frequency
Frequency
(3)
500 1000 1500 Single MUAp (μVms)
20 10
500
1000
1500
2000
Single MUAp (μVms)
0
500 1000 1500 Single MUAp (μVms)
2000
Fig. 4.5 The diagram shows how the between unit variation or the distribution of the single MUAP size (taken from MCMC output) for a particular patient taken on four occasions over an approximate 16-month period. Note how the proportion of large sized units increases as the disease progresses.
In ALS the within unit variation is of interest because it appears to evolve with the disease. Fig. 4.5 shows an example of the evolution of the between unit variability of an ALS patient over a 16-month period. Note that both the mean and variance of the single motor unit potentials increase with time. 3.3. The additivity of the mean and variance of CMAP area Fig. 4.6 is a DAG describing the variability of motor units at the observation level. The expected CMAP area (given knowledge of the motor units that fire, skt, their sizes, mk, together with the mean value of the baseline noise, mb) is expressed as: mTt ¼ mb þ
N X k¼1
mk;t S k;t :
The variability around the expected value was more problematical. In our initial models we assumed that the observations of CMAP area, yt, were conditionally independent and distributed normally about the expected CMAP, mtT, with variance Vt. The variance of the CMAP area about its expected value (Vt), was explained in terms of the within unit variability, s2 , the number of units that fire, nt, and the baseline variability using the expression: V t ¼ nt s2 þ s2b : The expressions for mt and Vt show how the mean and variance of the CMAP increase as the number of units firing increase. This can also be also be seen by examining Fig. 4.7. A closer look at the data has shown us that a normal distribution has insufficient probability mass in its tails (or extremes) to explain the CMAP area. To create a model more robust to outlying observations, we therefore allow the variances of each observation to be determined by a
14000
14000
12000
12000
10000
10000
Expected CMAP
CMAP area
Fig. 4.6 The DAG shows the relevant variables at the observation level. The CMAP is expressed as a normal distribution with its own mean and variance. Both the mean and the variance depend on the motor units that are firing at a given time. The node qt is introduced to explain the heavy tails in the distribution of the CMAP which gives the CMAP a marginal t-distribution with n degrees of freedom.
8000
6000
8000
6000
4000
4000
2000
2000
0 25
30
35
40 45 Stimulus
50
55
0
0
20
40
60 80 100 Every third observation
120
140
Fig. 4.7 The stimulus–response curve is shown on the left. The diagram on the right shows the expected CMAP, mtT, of every third observation. The contributing single MUAPs are also shown. The actual CMAPs are denoted by a + symbol. Note that the variance of the observations about their expected values increases as the number of units firing increases.
54
Fig. 4.8 Our assumptions about the structure of the model are displayed in a DAG. In RJMCMC the stochastic nodes are initialised and each is updated sequentially so that the chain of values from each node converges to the posterior distribution.
divisor, qt, from a gamma(a,a) distribution with unit mean. Thus, in this improved model, the variance at time t can be expressed as: Vt ¼
nt s2 þ s2 : qt
The parameter of the gamma distribution, a, controls the size of the tails of the distribution. This gives the CMAP marginally a t-distribution with 2a degrees of freedom (df ). In our model we set a = 2 giving a t-distribution with 4 df which is commonly used by statisticians for robust models. All of the assumptions of our
model can be expressed in a DAG and are shown in Fig. 4.9. The CMAP is shown as a square node and is denoted by yt where t is an instance in time. The stimulus is also fixed and also denoted by a square node denoted by St. 4. The reversible jump Markov Chain Monte Carlo (RJMCMC) MCMC is a computational algorithm from which we are able to simulate from the posterior distributions of each parameter. After initialization, each of the stochastic nodes of Fig. 4.9 is updated so that the
55 100 90 80 70
N
60 50 40 30 20 10 0
0
100
200 300 time in days (t)
400
500
Fig. 4.9 The posterior distributions of the number of motor units taken from stimulus–response data from a patient over almost a 2-year period are shown. Each posterior distribution is summarised by a mode shown as a short horizontal bar and a 95% credible interval shown as a vertical line. We have fitted an exponential curve to calculate the exponential rate of decline and from that we have calculated the expected half-life for the motor units from a particular muscle on a particular patient.
resultant Markov chain converges to the correct posterior distribution. A feature of our Markov chain is that the number of parameters or the dimension of the model (N) changes (Ridall et al., 2007). This is a particular kind of MCMC known as reversible jump MCMC (RJMCMC) (Green, 1995). A well designed RJMCMC should be able to stochastically visit all possible models and needs to take into account the peculiarities of a particular model. In the case of MUNE, we have specifically designed an RJMCMC that is able to switch between models that use alternation and those that do not use it. 5. Using the output of the RJMCMC to measure decline Figure 4.9 shows how RJMCMC can be performed on serial studies to obtain the rate of exponential decline. The posterior distributions of the number of motor units at various instances
of time are summarised by 95% credible intervals (using vertical bars) and the modal estimates (using short horizontal lines). A weighted exponential fit gives us an estimate of the half-life of a motor unit. In the example shown in Fig. 4.9 the half-life is estimated as 268 days with a 95% credible interval of between 246 and 295 days. 6. Conclusions The aim of this paper has been to show how the stimulus response–curve can be used both to conduct a new type of MUNE and also to estimate quantities that could be of interest to the neurologist. To do this we have constructed a model that is biologically realistic and follows from motor unit physiology (Ridall et al., 2006). The computational algorithm for MUNE is based on the RJMCMC methodology of Green (1995) and described in Ridall et al. (2006, 2007). This results in simulation from the posterior distribution of
56 the number of remaining motor units supplying a given muscle. A useful measure of decline of a patient can be obtained from these posterior distributions by applying an exponential fit and calculating the expected half-life. In addition to carrying out MUNE, our RJMCMC is also able to find the posterior distributions of parameters of interest describing the population of motor units being studied. These include estimates of the between motor unit variation, the within unit variation, the mean and variance of the size of the units and the distribution of the excitability parameters.
References Green, P.J. (1995) Reversible jump Markov chain Monte Carlo computation and Bayesian model determination. Biometrika, 82: 711–732. Ridall, P., Pettitt, A.N., Henderson, R.D. and McCombe, P.A. (2006) Motor unit number estimation — a Bayesian approach. Biometrics, doi: 10.1111/j.1541–0420.2006.00577.x. Ridall, P.G., Pettitt, A.N., Friel, N., McCombe, P.A. and Henderson, R.D. (2007) Motor unit number estimation using reversible jump Markov chain Monte Carlo. Appl. Statist., 56: 1–26.
Chapter 4
Bayesian analysis of the stimulus–response curve P.G. Ridalla,1,*, A.N. Pettitta, P.A. McCombeb and R.D. Hendersonc a
Faculty of Science, School of Mathematical Sciences, Queensland University of Technology, GPO Box 2434, Brisbane, QLD 4001, Australia b University of Queensland, Centre for Clinical Research, Royal Brisbane Hospital, Herston Queensland 4029, Australia
c
Department of Neurology, Royal Brisbane and Women’s Hospital, Herston, QLD 4029, and University of Queensland, Brisbane, Australia 1 Present address: Department of Mathematics and Statistics, Fylde College, Lancaster University, Lancaster LA1 4YF, UK
1. Introduction In this paper we present our Bayesian method for carrying out motor unit number estimation (MUNE) using stimulus–response data collected from surface electrophysiological recordings. We formulate and justify our assumptions in Ridall et al. (2006) and base these on available scientific evidence, as outlined in the previous paper. The object of our methodology is to express the uncertainty about the number of motor units in a muscle in terms of a posterior distribution. From studies taken over time, these posterior distributions can be used to estimate the rate of loss of motor units. A by-product of this method of MUNE is that it provides a means of tracking a given population of motor units over time. Examples of the para*
Correspondence to: Dr. P.G. Ridall (at present address), Department of Mathematics and Statistics, Fylde College, Lancaster University, Lancaster LA1 4YF, UK. Tel.: +44 (0) 1524-592302; Fax: +44 (0) 1524-593429; E-mail: [email protected]
meters that can be tracked over time include the distribution of the excitability parameters, the single motor unit action potentials and the between and within motor unit variability. In Section 2, we describe the fundamentals of Bayesian methodology and terminology and discuss why it is particularly appropriate for the MUNE problem. In Section 3, we outline the three essential components of the model. In Section 4, we explain how our computational method is able to explore the variety of possible models. In Section 5 we show how the output of the Bayesian analysis can be used to measure the rate of loss of motor units in patients with amyotrophic lateral sclerosis (ALS). Finally, in Section 6 we present a brief summary of our conclusions. 2. Bayesian statistics and directed acyclic graphs (DAGs) The attraction of the Bayesian paradigm is that it is a natural and logical way of carrying out statistical reasoning free from the ad hoc thinking that currently permeates traditional statistics.
48 The Bayesian paradigm recognises that for any statistical problem there are many types of uncertainty and that each of these types of uncertainty can be expressed by an appropriate probability distribution. Uncertainty in a parameter, before the data are collected, is expressed using a prior distribution. Once the data have been collected and their variability taken into account through the likelihood, this prior knowledge is updated to obtain what is known as the posterior distribution, the focus of our interest. The primary aim of the Bayesian statistical method of MUNE is to find the posterior probability of the number of motor units in a muscle. However, in a biologically realistic model there are many other unknowns that also need to be estimated. Each unknown must be found by assuming an appropriate prior distribution and then finding its posterior distribution. Sometimes we refer to these parameters as nuisance parameters and we describe what they are and how to deal with them below. A convenient and non-mathematical way of describing the assumptions in our model is using a directed acyclic graph (DAG). It shows the relationship between the unknown stochastic parameters and the data that are assumed to be known. In the graph unknowns are drawn as circular nodes and square nodes are used for known data. We now show how a DAG can be used to illustrate the following statistical concepts: (a) the likelihood, (b) conditional independence, (c) nuisance parameters, and finally (d) Bayesian inference and the derivation of the posterior distribution. These concepts are needed for a satisfactory understanding of how our method of MUNE works. Figure 4.1A shows an example of likelihood. The likelihood, represented by a square node, is the probability of the data (Y) given an unknown parameter (l) and is denoted by P(Y|l). Here l could be a Poisson mean and Y an observation of count data. Likelihood is not solely a Bayesian idea and many non-Bayesians use it, for example to find maximum likelihood estimates. However, Bayesians apply a different type of reasoning. We use panel 4.1D to explain this. Fig. 4.1B represents a plate (the outer rectangle), contain-
ing a group of observations (Yi) that (depicted in this way) are assumed to be conditionally independent, given l. In a DAG, the arrows initially point from the priors (the unknowns) to the data shown in the likelihood. A prior distribution P(l) is able to express the degree of knowledge of a given parameter before the data are collected or taken into account. This degree of knowledge can range from a total lack of knowledge to strong prior knowledge. Figure 4.1C shows the added complication of assuming the existence of a nuisance parameter, S. This nuisance parameter may be of no interest in itself, but its calculation is needed because it provides the link between the unknown and the observed data. The DAG which we use for MUNE has many nuisance parameters. Figure 4.1D shows a re-expression of the model this time with the unknown (the quantity of interest) as the subject. This reversal of arrows (from the direction of the likelihood) is a simple way of showing how Bayesians carry out their statistical reasoning and is a simple and natural idea. It gives the various probability distributions (prior, likelihood) the same mathematical status. The posterior distribution is the probability of the unknown parameter given the data and is shown as P(l|Y). This expresses knowledge about an unknown obtained from the data and is obtained as a result of applying Bayes theorem. This cannot be done without a prior distribution on the parameter of interest. In very simple examples the posterior distribution can be expressed in algebraic terms or the posterior distribution can be assumed to have a known distribution — for example, a normal distribution with mean and variance known in terms of the data. However, in the presence of many nuisance parameters, this often cannot be done and in these cases we must resort to a standard computational technique called Markov chain Monte Carlo (MCMC). MCMC has, over the last 20 years, become an increasingly useful computational tool because of greatly increased computer speeds. MCMC involves an iterative process with the chain taking sequentially changed values. The values of the chain are initialized and then
49
Fig. 4.1 Panel (A) shows a simple likelihood. Panel (B) shows a likelihood consisting of conditionally independent variables. Panel (C) shows the likelihood in the presence of a nuisance parameter (in this case S). Panel (D) shows the posterior distribution with the direction of the arrows in the opposite direction to that of the likelihood.
changed until convergence is attained. Convergence implies that the initial values no longer affect the values produced from the chain. Once convergence is obtained the output of the MCMC is from the posterior distribution. One way of representing the values produced from the posterior is by a histogram of values generated from the posterior (Fig. 4.2). For MUNE, the mode of the posterior gives us the best estimate of motor unit numbers with the spread of the histogram expressing the uncertainty in that estimate. The histogram shown in Fig. 4.2 is taken from 100,000 values generated from the chain. The mode, or the most likely number of motor units, is found to be N = 13 units. The 95% credible interval in this case is
[12, 15]. This means that this interval contains 95% of the posterior distribution and has a different meaning to a classical confidence interval. 3. The elements of our statistical model In the specific problem of MUNE, uncertainty exists at distinct levels and each is related to motor unit physiology. Firstly, there is uncertainty at the level of the axonal threshold which describes whether a motor unit will fire when a stimulus is applied. Secondly, there are two types of variation that contribute to the uncertainty in the compound muscle action potential (CMAP) at the observation level. These are the between motor unit variation (uncertainty at the level of the single motor
50
Fig. 4.2
The posterior probability of the number of motor units, N, expressed as a histogram.
unit response) and the within motor unit variation. Our distribution of CMAPs must be robust enough to describe outlying observations that occur, particularly in the latter stages of ALS. 3.1. The depolarisation of the axon membrane by a stimulus We use the DAG in Fig. 4.3 to illustrate this part of the model. Unit k only fires at time t if, and only if, the stimulus St exceeds the threshold of that motor unit, tk,t, at that instance of time. The variable skt indicates whether a motor unit is on or off at that time. The threshold is not assumed fixed in time but to vary independently as a Gaussian variable with its own mean threshold, mk, and variance parameter, uk. These two parameters are called the excitability parameters. We illustrate the meaning of these two parameters in Fig. 4.4, using the modal output of the Markov chain used to analyse the data shown in the top panel of Fig. 4.4. The first of these, the mean threshold of a motor unit, mk, is defined as the stimulus for which a motor unit has a 50% probability of firing.
The first two mean thresholds of the data shown in the top panel are shown in the second panel of Fig. 4.4. The second set of parameters, the variances of the thresholds, uk, determine the steepness of the corresponding excitability curves and in doing so define the range over which the unit displays probabilistic firing. The region of probabilistic firing (ROPF) is shown in the bottom panel of Fig. 4.4B as a double arrow for the fourth unit. The first three excitability curves on the bottom panel of Fig. 4.4 are not as steep as the rest indicating that the first three variance threshold parameters, u1, u2 and u3, are much larger than those for the rest of the motor units. 3.2. The response at the neuromuscular junction We let the average area of a single motor unit action potential (MUAP) be mk and refer to this as the size of a particular unit. The variance of the distribution of mk is referred to as the between motor unit variation. We assume that each single MUAP is normally distributed about its mean, mk and with a variance s2 common to all units.
51
Fig. 4.3
The DAG for showing the relationship of the relevant variables for describing the depolarisation of an axon.
CMAP area (μVms)
8000 6000 4000 2000 0 45
50
55
60 Stimulus (mA)
65
70
75
1
p(St>τkt)
0.8
m2 m1 R.O.P.F.
0.6
p(St>τkt)=0.5
0.4 0.2 0 45
50
55
60 Stimulus (St)
65
70
75
Fig. 4.4 The top panel shows the stimulus–response curves plotted for the data. The bottom panel shows the excitability curves obtained from the MCMC output.
52 (2) 40
30
30 Frequency
Frequency
(1) 40
20 10 0
20 10
500
1000
1500
0
2000
Single MUAp (μVms) 40
30
30
20 10 0
2000
(4)
40
Frequency
Frequency
(3)
500 1000 1500 Single MUAp (μVms)
20 10
500
1000
1500
2000
Single MUAp (μVms)
0
500 1000 1500 Single MUAp (μVms)
2000
Fig. 4.5 The diagram shows how the between unit variation or the distribution of the single MUAP size (taken from MCMC output) for a particular patient taken on four occasions over an approximate 16-month period. Note how the proportion of large sized units increases as the disease progresses.
In ALS the within unit variation is of interest because it appears to evolve with the disease. Fig. 4.5 shows an example of the evolution of the between unit variability of an ALS patient over a 16-month period. Note that both the mean and variance of the single motor unit potentials increase with time. 3.3. The additivity of the mean and variance of CMAP area Fig. 4.6 is a DAG describing the variability of motor units at the observation level. The expected CMAP area (given knowledge of the motor units that fire, skt, their sizes, mk, together with the mean value of the baseline noise, mb) is expressed as: mTt ¼ mb þ
N X k¼1
mk;t S k;t :
The variability around the expected value was more problematical. In our initial models we assumed that the observations of CMAP area, yt, were conditionally independent and distributed normally about the expected CMAP, mtT, with variance Vt. The variance of the CMAP area about its expected value (Vt), was explained in terms of the within unit variability, s2 , the number of units that fire, nt, and the baseline variability using the expression: V t ¼ nt s2 þ s2b : The expressions for mt and Vt show how the mean and variance of the CMAP increase as the number of units firing increase. This can also be also be seen by examining Fig. 4.7. A closer look at the data has shown us that a normal distribution has insufficient probability mass in its tails (or extremes) to explain the CMAP area. To create a model more robust to outlying observations, we therefore allow the variances of each observation to be determined by a
14000
14000
12000
12000
10000
10000
Expected CMAP
CMAP area
Fig. 4.6 The DAG shows the relevant variables at the observation level. The CMAP is expressed as a normal distribution with its own mean and variance. Both the mean and the variance depend on the motor units that are firing at a given time. The node qt is introduced to explain the heavy tails in the distribution of the CMAP which gives the CMAP a marginal t-distribution with n degrees of freedom.
8000
6000
8000
6000
4000
4000
2000
2000
0 25
30
35
40 45 Stimulus
50
55
0
0
20
40
60 80 100 Every third observation
120
140
Fig. 4.7 The stimulus–response curve is shown on the left. The diagram on the right shows the expected CMAP, mtT, of every third observation. The contributing single MUAPs are also shown. The actual CMAPs are denoted by a + symbol. Note that the variance of the observations about their expected values increases as the number of units firing increases.
54
Fig. 4.8 Our assumptions about the structure of the model are displayed in a DAG. In RJMCMC the stochastic nodes are initialised and each is updated sequentially so that the chain of values from each node converges to the posterior distribution.
divisor, qt, from a gamma(a,a) distribution with unit mean. Thus, in this improved model, the variance at time t can be expressed as: Vt ¼
nt s2 þ s2 : qt
The parameter of the gamma distribution, a, controls the size of the tails of the distribution. This gives the CMAP marginally a t-distribution with 2a degrees of freedom (df ). In our model we set a = 2 giving a t-distribution with 4 df which is commonly used by statisticians for robust models. All of the assumptions of our
model can be expressed in a DAG and are shown in Fig. 4.9. The CMAP is shown as a square node and is denoted by yt where t is an instance in time. The stimulus is also fixed and also denoted by a square node denoted by St. 4. The reversible jump Markov Chain Monte Carlo (RJMCMC) MCMC is a computational algorithm from which we are able to simulate from the posterior distributions of each parameter. After initialization, each of the stochastic nodes of Fig. 4.9 is updated so that the
55 100 90 80 70
N
60 50 40 30 20 10 0
0
100
200 300 time in days (t)
400
500
Fig. 4.9 The posterior distributions of the number of motor units taken from stimulus–response data from a patient over almost a 2-year period are shown. Each posterior distribution is summarised by a mode shown as a short horizontal bar and a 95% credible interval shown as a vertical line. We have fitted an exponential curve to calculate the exponential rate of decline and from that we have calculated the expected half-life for the motor units from a particular muscle on a particular patient.
resultant Markov chain converges to the correct posterior distribution. A feature of our Markov chain is that the number of parameters or the dimension of the model (N) changes (Ridall et al., 2007). This is a particular kind of MCMC known as reversible jump MCMC (RJMCMC) (Green, 1995). A well designed RJMCMC should be able to stochastically visit all possible models and needs to take into account the peculiarities of a particular model. In the case of MUNE, we have specifically designed an RJMCMC that is able to switch between models that use alternation and those that do not use it. 5. Using the output of the RJMCMC to measure decline Figure 4.9 shows how RJMCMC can be performed on serial studies to obtain the rate of exponential decline. The posterior distributions of the number of motor units at various instances
of time are summarised by 95% credible intervals (using vertical bars) and the modal estimates (using short horizontal lines). A weighted exponential fit gives us an estimate of the half-life of a motor unit. In the example shown in Fig. 4.9 the half-life is estimated as 268 days with a 95% credible interval of between 246 and 295 days. 6. Conclusions The aim of this paper has been to show how the stimulus response–curve can be used both to conduct a new type of MUNE and also to estimate quantities that could be of interest to the neurologist. To do this we have constructed a model that is biologically realistic and follows from motor unit physiology (Ridall et al., 2006). The computational algorithm for MUNE is based on the RJMCMC methodology of Green (1995) and described in Ridall et al. (2006, 2007). This results in simulation from the posterior distribution of
56 the number of remaining motor units supplying a given muscle. A useful measure of decline of a patient can be obtained from these posterior distributions by applying an exponential fit and calculating the expected half-life. In addition to carrying out MUNE, our RJMCMC is also able to find the posterior distributions of parameters of interest describing the population of motor units being studied. These include estimates of the between motor unit variation, the within unit variation, the mean and variance of the size of the units and the distribution of the excitability parameters.
References Green, P.J. (1995) Reversible jump Markov chain Monte Carlo computation and Bayesian model determination. Biometrika, 82: 711–732. Ridall, P., Pettitt, A.N., Henderson, R.D. and McCombe, P.A. (2006) Motor unit number estimation — a Bayesian approach. Biometrics, doi: 10.1111/j.1541–0420.2006.00577.x. Ridall, P.G., Pettitt, A.N., Friel, N., McCombe, P.A. and Henderson, R.D. (2007) Motor unit number estimation using reversible jump Markov chain Monte Carlo. Appl. Statist., 56: 1–26.