In Chapter 13, program minNforHDIpower.R is used to compute the minimum sample size needed for desired power. The program that was previously posted had a few glitches that have now been corrected. (Specifically, the source()-ing HDIofICDF.R needed quotes, and the call to HDIofICDF needed to include the argument credMass=HDImass. I also modified the condition at the beginning to use xor().) The results reported in the book are correct. You can download the updated program from the link in the book website.
If there are still any problems, please let me know.
Many thanks to Nico Corea for bringing this issue to my attention.
Friday, January 13, 2012
Tuesday, January 10, 2012
New procedure required for installing BRugs and OpenBUGS
Before R 2.14, OpenBUGS came packaged with BRugs, so all you had to do was install BRugs, and OpenBUGS magically came along for the ride. Now, you must install OpenBUGS as a separate step in the process. Of course, if you are using JAGS, you do not need to install OpenBUGS or BRugs. But if you want to use BRugs and OpenBUGS, here's what to do:
1. Install R from this site. Use the 32-bit Windows version.
2. Install OpenBUGS from this site.
3. Invoke R (32-bit Windows), and type install.packages("BRugs")
I recommend using JAGS instead, as explained at this previous blog post.
Warning message produced by R 2.14.1 in multiple linear regression programs
R version 2.14.1 produces a warning message that did not occur in earlier versions of R. Specifically, the sd(y) function now produces a warning message if y is a column matrix and not a vector. This situation arises in the various programs for multiple linear regression. I have therefore replaced the offending uses of sd(y) with apply(y,2,sd). You can download the updated programs from the link in the book website.
If I missed any occurrences in these or other programs, please let me know.
Many thanks to Agustin Alonso Rodriguez for bringing this issue to my attention.
If I missed any occurrences in these or other programs, please let me know.
Many thanks to Agustin Alonso Rodriguez for bringing this issue to my attention.
Friday, January 6, 2012
Complete example of right censoring in JAGS (with rjags)
In this post I provide a complete and simple example of how to do estimation of right-censored data in JAGS (with rjags). I could not find a complete, simple example anywhere (on the web or in print), and it took me lots of trial and error to figure out the tacit assumptions, which I try to make explicit here. (Of course, there still might be errors or infelicities. Let me know.)
Before launching into programming details, let me illustrate why this is important. Right-censored data can occur in many situations. For example, when measuring response times, the subject is given a limited temporal window in which to respond, and if no response has occurred by the end of the window, then the trial is aborted and the response is recorded as censored. If we assume that the subject was "on task" but merely didn't have enough time to respond before the experimenter got bored, then it would be inappropriate to omit this trial from the modeled data. In other words, the trial is providing genuine information about the response time, namely, that it takes longer than the experimenter's censoring limit. If we omit this trial from the modeled data, then our estimates will be too small!
Here is an example. The data are 500 values generated from a normal distribution, rescaled so that the sample mean is 100.0 and the sample standard deviation is 15.0. Then I marked any value above +1SD as NA. In other words, I censored about 15% of the data in the right tail. Suppose we model the data with a normal likelihood and estimate its μ and σ parameters. We would like the estimates to recover the true generating values of μ=100.0 and σ=15.0. But if we simply omit the censored values (as if they weren't part of the data), the estimates look like this:
Notice that the estimated parameter values are too low. This makes sense, because the model is trying merely to describe the uncensored data shown in grey, without any attempt to take into account the censored data.
On the other hand, if we tell JAGS that there are censored values that it should take into account, the estimates look like this:
Notice that the estimated parameter values are very accurate. This improvement in accuracy (and meaningfulness!) is why we want to be able to model censored data.
The basic idea for implementing in JAGS is presented briefly (one paragraph) in the JAGS user manual. The idea starts with defining an indicator vector that encodes, for every data value, whether the data value was censored or not. The indicator value, here named isCensored, is 1 for censored values and 0 for non-censored values. Then that indicator value is modeled thus: isCensored ~ dinterval(y,censorLimit) It took me a while to infer from that syntax the key idea: When isCensored is 1 and y is missing, then all JAGS knows is that y is somewhere above censorLimit, so JAGS effectively imputes a random value for y from the likelihood function (defined separately, in the usual way). Here are the tacit assumptions that I eventually figured out:
Before launching into programming details, let me illustrate why this is important. Right-censored data can occur in many situations. For example, when measuring response times, the subject is given a limited temporal window in which to respond, and if no response has occurred by the end of the window, then the trial is aborted and the response is recorded as censored. If we assume that the subject was "on task" but merely didn't have enough time to respond before the experimenter got bored, then it would be inappropriate to omit this trial from the modeled data. In other words, the trial is providing genuine information about the response time, namely, that it takes longer than the experimenter's censoring limit. If we omit this trial from the modeled data, then our estimates will be too small!
Here is an example. The data are 500 values generated from a normal distribution, rescaled so that the sample mean is 100.0 and the sample standard deviation is 15.0. Then I marked any value above +1SD as NA. In other words, I censored about 15% of the data in the right tail. Suppose we model the data with a normal likelihood and estimate its μ and σ parameters. We would like the estimates to recover the true generating values of μ=100.0 and σ=15.0. But if we simply omit the censored values (as if they weren't part of the data), the estimates look like this:
 |
| Posterior estimates when censored data are omitted. Estimates are too low. (ROPE is shown only for visual landmarks; it's not particularly meaningful in this case.) |
On the other hand, if we tell JAGS that there are censored values that it should take into account, the estimates look like this:
 |
| Posterior estimates when censored data are modeled. Estimates are accurate. (ROPE is shown only for visual landmarks; it's not particularly meaningful in this case.) |
The basic idea for implementing in JAGS is presented briefly (one paragraph) in the JAGS user manual. The idea starts with defining an indicator vector that encodes, for every data value, whether the data value was censored or not. The indicator value, here named isCensored, is 1 for censored values and 0 for non-censored values. Then that indicator value is modeled thus: isCensored ~ dinterval(y,censorLimit) It took me a while to infer from that syntax the key idea: When isCensored is 1 and y is missing, then all JAGS knows is that y is somewhere above censorLimit, so JAGS effectively imputes a random value for y from the likelihood function (defined separately, in the usual way). Here are the tacit assumptions that I eventually figured out:
- Censored data must be recorded as NA, not as the value of censoring limit.
- When explicitly initializing the chains, the censored values of the data must be explicitly initialized (to values above the censoring limits)!
Sunday, January 1, 2012
Now in JAGS! Now in JAGS!
I have created JAGS versions of all the BUGS programs in Doing Bayesian Data Analysis. Unlike BUGS, JAGS runs on MacOS, Linux, and Windows. JAGS has other features that make it more robust and user-friendly than BUGS. I recommend that you use the JAGS versions of the programs. Please let me know if you encounter any errors or inaccuracies in the programs.
JAGS is an MCMC sampler much like BUGS. To communicate with JAGS from R, a library called rjags is loaded into R. The diagram below illustrates the analogous roles of the sampling programs and the interface libraries:
Even if you use Windows, JAGS is nicer than BUGS for many reasons. One nicety is that JAGS has a dynamic progress bar that tells you how much of the requested chain has been sampled. JAGS also seems more robust and actually works in many situations where BUGS mysteriously fails; e.g., the equals(,) function works in JAGS where it would not work in BUGS.
It is easy to install JAGS and rjags. First, go to the JAGS web site and follow the instructions for downloading and installing JAGS for your operating system. Then, at the command line in R, type install.packages("rjags"). Done!
It is easy to get the JAGS versions of the programs for Doing Bayesian Data Analysis. JAGS versions of the programs use the same name as the BUGS versions, but with the string "Bugs" or "BRugs" replaced with "Jags". (All of the original programs are still available.) A zip file with all the programs and data files is available here. A list of individual programs is available here; click on the column header "last modified" twice to get the most recent files to appear at the top of the list. [Update January 12, 2012: A few programs (BernBeta...) were inadvertently missing from the upload of January 1. I've now posted them with the other programs. Thanks to a reader for pointing out the missing files.]
As explained in other blog posts, there are other revisions and additions to the programs.
* There is a new program for split-plot designs, called SplitPlotJags.R. See this blog post.
* The new versions of the programs do no thinning but use longer chains (defaulting to chain lengths of 50000, but you can use even longer chains for real research reports). See this blog post.
* The ANOVA-like programs consider only the sum-to-zero (STZ) versions of the parameters, because autocorrelation in the pre-STZ parameters is irrelevant. Use the versions of the programs with filenames that end with "STZ". See this blog post.
* Instead of saving plots using dev.copy2eps, the programs use savePlot. It is easy to change the savePlot command to save in whatever format you prefer, such as jpg/jpeg. Each program also begins with a check of the user's operating system to redefine the windows() command if necessary. Both of these changes were recommended by comments from readers.
JAGS is an MCMC sampler much like BUGS. To communicate with JAGS from R, a library called rjags is loaded into R. The diagram below illustrates the analogous roles of the sampling programs and the interface libraries:
 |
| JAGS, OpenBUGS, and WinBUGS are MCMC samplers. Each has a corresponding package for communicating with R. The book used OpenBUGS with BRugs. But the new programs feature JAGS with rjags. |
It is easy to install JAGS and rjags. First, go to the JAGS web site and follow the instructions for downloading and installing JAGS for your operating system. Then, at the command line in R, type install.packages("rjags"). Done!
It is easy to get the JAGS versions of the programs for Doing Bayesian Data Analysis. JAGS versions of the programs use the same name as the BUGS versions, but with the string "Bugs" or "BRugs" replaced with "Jags". (All of the original programs are still available.) A zip file with all the programs and data files is available here. A list of individual programs is available here; click on the column header "last modified" twice to get the most recent files to appear at the top of the list. [Update January 12, 2012: A few programs (BernBeta...) were inadvertently missing from the upload of January 1. I've now posted them with the other programs. Thanks to a reader for pointing out the missing files.]
As explained in other blog posts, there are other revisions and additions to the programs.
* There is a new program for split-plot designs, called SplitPlotJags.R. See this blog post.
* The new versions of the programs do no thinning but use longer chains (defaulting to chain lengths of 50000, but you can use even longer chains for real research reports). See this blog post.
* The ANOVA-like programs consider only the sum-to-zero (STZ) versions of the parameters, because autocorrelation in the pre-STZ parameters is irrelevant. Use the versions of the programs with filenames that end with "STZ". See this blog post.
* Instead of saving plots using dev.copy2eps, the programs use savePlot. It is easy to change the savePlot command to save in whatever format you prefer, such as jpg/jpeg. Each program also begins with a check of the user's operating system to redefine the windows() command if necessary. Both of these changes were recommended by comments from readers.
[Added 28 January: Complete installation instructions are listed here.]
Wednesday, December 28, 2011
Split-Plot Design in JAGS (preliminary version)
For many months I've been meaning to create the code for a split-plot design. The basic split-plot design has each subject participate in every level of factor B, but only one level of factor A. In this blog post, I report an example of a hierarchical Bayesian approach to a split-plot design, coded in JAGS (not BUGS). As this is my first go at this sort of design, please consider this a preliminary version of the analysis, and please provide feedback regarding any errors or infelicities. I'm also looking for people to do analogous NHST analyses for comparison, and I'm looking for other data sets for additional examples.
To make things concrete, consider an example provided by Maxwell & Delaney (2004, Designing Experiments and Analyzing Data: A Model Comparison Perspective, 2nd Edition, Erlbaum; Ch. 12). (As I've mentioned elsewhere, if you must learn NHST, their book is a great resource.) A perceptual psychologist is studying response times for identifying letters flashed on a screen. The letters can be rotated off of vertical by zero degrees, four degrees, or eight degrees. Every subject responds many times to letters at each of the three angles. The experimenter (for unknown reasons) analyzes only the median response time of each subject at each angle. Thus, each subject contributes only one datum at each level of angle (factor B). There are two types of subjects: "young" and "old." Age of subject is factor A. The structure of the data is shown below. Y is the median response time, in milliseconds. There are 10 subjects per Age group.
Split-plot designs are difficult to analyze in NHST because it is challenging to select an appropriate "error term" for constructing F ratios of different effects. Ch. 12 of Maxwell & Delaney is devoted largely to explaining the various options for error terms and the different F ratios (and different p values) that result. Split-plot designs get even more difficult when they are unbalanced, that is, when there are different numbers of subjects in the different levels of factor A. And, of course, when we do multiple comparisons we have to worry over which corrections to use.
In a Bayesian approach, on the other hand, there is no need to worry over the selection of error terms because we never construct F ratios. And there is no difficulty with unbalanced designs because the parameters are estimated based on whatever data are provided. The only challenge, in the hierarchical Bayesian approach that I prefer, is the construction of sum-to-zero parameterized effects. In all the ANOVA designs described in the book, the parameter values are converted so the effects, as deflections from baseline, sum to zero. (In case you haven't seen it yet, there is an important update regarding the sum-to-zero computations in the book's ANOVA programs, linked here.) With some effort, I figured out a way to convert the parameter estimates to sum-to-zero versions. It is done in R, outside of JAGS.
The model involves main effects for factor A and factor B, an interaction effect for AxB, and a main effect for subject within level of A. The basic model specification looks like this:
for ( i in 1:Ntotal ) {
y[i] ~ dnorm( mu[i] , tau )
mu[i] <- base + a[aLvl[i]] + s[sLvl[i]] + b[bLvl[i]] + axb[aLvl[i],bLvl[i]]
}
The index i goes over rows of the data table, and Ntotal is the total number of data points. The data values, y[i], are assumed to be normally distributed around the predicted mean mu[i]. The predicted mean is the overall baseline, base, plus a deflection a[aLvl[i]] due to being in level aLvl of factor A, plus a deflection for each subject, plus a deflection for factor B, plus a deflection due to interaction of factors A and B. There is a hierarchical prior on each type of deflection, so that the variance of the deflections is estimated.
The tricky part comes in converting the parameter values to sum-to-zero values. Essentially, at each step in the chain the predicted mu is computed, then the cell means for the AxB table are computed. From the cell means, the baseline, A, B, and AxB deflections are computed. And, from the cell means, the within-level deflections of the subjects are computed. (The difficult part for me was figuring out that I should compute the cell means of the AxB table first, and compute everything else, including the subject effects, relative to those cell means.)
The first figure (below) shows marginals of the posterior standard deviations on the various effects and the noise standard deviation. The estimates of the standard deviations are not of primary interest, but here they are:
The next plot (below) shows the baseline, factor A deflections, factor B deflections, and AxB interaction deflections. Notice that the deflections do indeed sum to zero within each type:
The deflections are informative primarily when they are entered into contrasts. For example, the main effect of factor A is quite large:
This agrees with the conclusion from NHST analysis from Maxwell and Delaney, who reported F(1,18)=7.28, p=.0147. But the Bayesian posterior seems much more powerful than the NHST p value. Did I do something wrong, or is the Bayesian analysis simply much more powerful in this design?
We can also do contrasts, such as comparing young versus old (A1 vs A2) at the angle of zero (B1) only. The result is shown in the left panel below, where it can be seen that a difference of zero is within the 95% HDI:
This result agrees with the conclusion from NHST reported by Maxwell & Delaney (pp. 602-605), with F=3.16 and p=.092, or F=3.08 and p still >.05, depending on the choice of error term.
The right panel above shows a contrast that investigates the difference of quadratic trends across the age groups. A quadratic trend for the three levels of angle has contrast coefficients of 1,-2,1. The question is whether the contrast is different across the two levels of age. The Bayesian analysis shows that a difference of zero is within the 95% HDI. This result agrees with conclusion from NHST reported by Maxwell & Delaney (pp. 605-607), with F(1,54)=3.101 and p>.05, or F(1,18)=4.038 with p>.05, depending again on the choice of error term.
As usual, the above contrasts would have to be corrected for multiple comparisons in an NHST approach.
It is seamless to analyze an unbalanced design in the Bayesian approach. For example, the analysis runs completely unchanged (with different results) if we simply delete the last subject from the data table, so that there are 10 subjects in level A1 but only 9 subjects in level A2.
The data file is linked here, and the program is linked here. To run the program, just put it in the same folder as the data file and make sure that R has that folder as its working directory. The program uses JAGS and rjags, not OpenBUGS and BRugs. You must first install JAGS. JAGS runs on all platforms, including Mac and Linux! And then in R you must install rjags; just type install.packages("rjags").
Please do comment!
To make things concrete, consider an example provided by Maxwell & Delaney (2004, Designing Experiments and Analyzing Data: A Model Comparison Perspective, 2nd Edition, Erlbaum; Ch. 12). (As I've mentioned elsewhere, if you must learn NHST, their book is a great resource.) A perceptual psychologist is studying response times for identifying letters flashed on a screen. The letters can be rotated off of vertical by zero degrees, four degrees, or eight degrees. Every subject responds many times to letters at each of the three angles. The experimenter (for unknown reasons) analyzes only the median response time of each subject at each angle. Thus, each subject contributes only one datum at each level of angle (factor B). There are two types of subjects: "young" and "old." Age of subject is factor A. The structure of the data is shown below. Y is the median response time, in milliseconds. There are 10 subjects per Age group.
| Y | Subj | Angle | Age |
| 450 | 1 | B1(Zero) | A1(Young) |
| 510 | 1 | B2(Four) | A1(Young) |
| 630 | 1 | B3(Eight) | A1(Young) |
| 390 | 2 | B1(Zero) | A1(Young) |
| 480 | 2 | B2(Four) | A1(Young) |
| 540 | 2 | B3(Eight) | A1(Young) |
| ... | ... | ... | ... |
| 510 | 10 | B1(Zero) | A1(Young) |
| 540 | 10 | B2(Four) | A1(Young) |
| 660 | 10 | B3(Eight) | A1(Young) |
| 420 | 11 | B1(Zero) | A2(Old) |
| 570 | 11 | B2(Four) | A2(Old) |
| 690 | 11 | B3(Eight) | A2(Old) |
| 600 | 12 | B1(Zero) | A2(Old) |
| 720 | 12 | B2(Four) | A2(Old) |
| 810 | 12 | B3(Eight) | A2(Old) |
| ... | ... | ... | ... |
| 510 | 20 | B1(Zero) | A2(Old) |
| 690 | 20 | B2(Four) | A2(Old) |
| 810 | 20 | B3(Eight) | A2(Old) |
Split-plot designs are difficult to analyze in NHST because it is challenging to select an appropriate "error term" for constructing F ratios of different effects. Ch. 12 of Maxwell & Delaney is devoted largely to explaining the various options for error terms and the different F ratios (and different p values) that result. Split-plot designs get even more difficult when they are unbalanced, that is, when there are different numbers of subjects in the different levels of factor A. And, of course, when we do multiple comparisons we have to worry over which corrections to use.
In a Bayesian approach, on the other hand, there is no need to worry over the selection of error terms because we never construct F ratios. And there is no difficulty with unbalanced designs because the parameters are estimated based on whatever data are provided. The only challenge, in the hierarchical Bayesian approach that I prefer, is the construction of sum-to-zero parameterized effects. In all the ANOVA designs described in the book, the parameter values are converted so the effects, as deflections from baseline, sum to zero. (In case you haven't seen it yet, there is an important update regarding the sum-to-zero computations in the book's ANOVA programs, linked here.) With some effort, I figured out a way to convert the parameter estimates to sum-to-zero versions. It is done in R, outside of JAGS.
The model involves main effects for factor A and factor B, an interaction effect for AxB, and a main effect for subject within level of A. The basic model specification looks like this:
for ( i in 1:Ntotal ) {
y[i] ~ dnorm( mu[i] , tau )
mu[i] <- base + a[aLvl[i]] + s[sLvl[i]] + b[bLvl[i]] + axb[aLvl[i],bLvl[i]]
}
The index i goes over rows of the data table, and Ntotal is the total number of data points. The data values, y[i], are assumed to be normally distributed around the predicted mean mu[i]. The predicted mean is the overall baseline, base, plus a deflection a[aLvl[i]] due to being in level aLvl of factor A, plus a deflection for each subject, plus a deflection for factor B, plus a deflection due to interaction of factors A and B. There is a hierarchical prior on each type of deflection, so that the variance of the deflections is estimated.
The tricky part comes in converting the parameter values to sum-to-zero values. Essentially, at each step in the chain the predicted mu is computed, then the cell means for the AxB table are computed. From the cell means, the baseline, A, B, and AxB deflections are computed. And, from the cell means, the within-level deflections of the subjects are computed. (The difficult part for me was figuring out that I should compute the cell means of the AxB table first, and compute everything else, including the subject effects, relative to those cell means.)
The first figure (below) shows marginals of the posterior standard deviations on the various effects and the noise standard deviation. The estimates of the standard deviations are not of primary interest, but here they are:
 |
| January 8, 2012: It was pointed out to me that the SD estimates for aSD (age), bSD (Angle), and axbSD are not very stable, even for chains of length 100,000. The reason is that only a few values inform the estimates. Specifically, for aSD, only two values inform it, namely a[1] and a[2], so it is very uncertain. And for bSD, only three values inform it, namely b[1], b[2], and b[3], so it too is very uncertain. I therefore included a warning message in the display. . Estimates of the effects (graphs below) are stable. |
The next plot (below) shows the baseline, factor A deflections, factor B deflections, and AxB interaction deflections. Notice that the deflections do indeed sum to zero within each type:
The deflections are informative primarily when they are entered into contrasts. For example, the main effect of factor A is quite large:
We can also do contrasts, such as comparing young versus old (A1 vs A2) at the angle of zero (B1) only. The result is shown in the left panel below, where it can be seen that a difference of zero is within the 95% HDI:
The right panel above shows a contrast that investigates the difference of quadratic trends across the age groups. A quadratic trend for the three levels of angle has contrast coefficients of 1,-2,1. The question is whether the contrast is different across the two levels of age. The Bayesian analysis shows that a difference of zero is within the 95% HDI. This result agrees with conclusion from NHST reported by Maxwell & Delaney (pp. 605-607), with F(1,54)=3.101 and p>.05, or F(1,18)=4.038 with p>.05, depending again on the choice of error term.
As usual, the above contrasts would have to be corrected for multiple comparisons in an NHST approach.
It is seamless to analyze an unbalanced design in the Bayesian approach. For example, the analysis runs completely unchanged (with different results) if we simply delete the last subject from the data table, so that there are 10 subjects in level A1 but only 9 subjects in level A2.
See revised version here.
The data file is linked here, and the program is linked here. To run the program, just put it in the same folder as the data file and make sure that R has that folder as its working directory. The program uses JAGS and rjags, not OpenBUGS and BRugs. You must first install JAGS. JAGS runs on all platforms, including Mac and Linux! And then in R you must install rjags; just type install.packages("rjags").
Please do comment!
Thursday, December 8, 2011
Some more nice reviews on Amazon.com
It's greatly appreciated when people go to the effort to write a nice review on Amazon. It's appreciated not only by the author :-) but crucially also by prospective readers who are trying to decide whether the book is worth getting. Here are some excerpts from some recent reviews on Amazon.com
This is one of the best written and accessible statistics books I've come across. Obviously, a lot of thinking went into coming up with examples and intuitive explanation of various ideas. I was consistently amazed at author's ability to not just say how something is done but why it is done that way using simple examples. I've read far more mathematically sophisticated explanations of statiscal modeling but, in this book,I felt I was allowed to peek into the mind of previous authors as to what they were really thinking when writing down their math formulas. (Posted November 11, 2011 by Davar314, San Francisco, CA)
As far as I am concerned, if you write a book this good, you get to put whatever you like on the cover - puppies, Angelina Jolie, even members of the metal band "Das Kruschke". While reading "DBDA" - reading *and* stepping through the code examples - will not make you a "Bayesian black-belt", it's impressive how much information it *will* give you - the book is almost 700 pages, after all - and you don't need (but it helps) to have tried to get the hang of the "Bayesian stuff" with other books to appreciate how friendly and effective this one is. (The author's explanation of the Metropolis algorithm is a good example). At the risk of sounding grandiose, the book just might do for Bayesian methods what Apple's original Mac did for the personal computer; here's hoping. (Posted December 7, 2011 by Dimitri Shvorob)Click here for the full text of all the reviews on Amazon. Thanks again, reviewers, for the nice comments and for helping prospective readers.
Subscribe to:
Posts (Atom)