- Applied Bayesian Modeling for the Social Sciences, Dave Armstrong
- Bayesian Methods, Nick Beauchamp
- Statistical & Cognitive Modeling for Formal Semantics, Adrian Brasoveanu
- Introduction to Bayesian Analysis, Brad Carlin
- Doing Bayesian Data Analysis, John Kruschke
- Applied Bayesian Statistics, Marco Steenbergen
- Statistics VI, Tuerlinckx and Vanpaemel (from which this book review arose)
- Bayesian Data Analysis, Uppsala University
- Topics in Quantitative Psychology: Bayesian Data Analysis (no web page), Caren Rotello, U. of Massachusetts
Saturday, June 30, 2012
Courses that use Doing Bayesian Data Analysis?
Instructors seek examples of using the book, Doing Bayesian Data Analysis, as part of a course. A cursory web search yielded these few listed below, but there must be others. For example, my own course web page did not show up in the search, nor did another course from which a published review of the book arose. Please let me know of other courses that use the book, and I will update this list in a subsequent post.
Monday, June 25, 2012
Workshop, July 11-13, U. Wisconsin, Madison
Wednesday, June 20, 2012
Solutions to exercises now available to everyone
Wednesday, June 6, 2012
Mixture of Normal Distributions
In this post I show a simple illustration of a mixture of normal distributions. For the examples, we assume we have metric values that we suppose are generated by a mixture of two different normal distributions, which I'll call clusters. We don't know which datum came from each cluster. Our goal is to estimate the probability that each score came from each of the two clusters, and the means and SD of the normal distributions that describe the clusters.
The model specification (for JAGS): The assumes that the clusters have the same standard deviation, but different means.
model {
# Likelihood:
for( i in 1 : N ) {
y[i] ~ dnorm( mu[i] , tau )
mu[i] <- muOfClust[ clust[i] ]
clust[i] ~ dcat( pClust[1:Nclust] )
}
# Prior:
tau ~ dgamma( 0.01 , 0.01 )
for ( clustIdx in 1: Nclust ) {
muOfClust[clustIdx] ~ dnorm( 0 , 1.0E-10 )
}
pClust[1:Nclust] ~ ddirch( onesRepNclust )
}
The data specification:
# Generate random data from known parameter values:
set.seed(47405)
trueM1 = 100
N1 = 200
trueM2 = 145 # 145 for first example below; 130 for second example
N2 = 200
trueSD = 15
effsz = abs( trueM2 - trueM1 ) / trueSD
y1 = rnorm( N1 )
y1 = (y1-mean(y1))/sd(y1) * trueSD + trueM1
y2 = rnorm( N2 )
y2 = (y2-mean(y2))/sd(y2) * trueSD + trueM2
y = c( y1 , y2 )
N = length(y)
# Must have at least one data point with fixed assignment
# to each cluster, otherwise some clusters will end up empty:
Nclust = 2
clust = rep(NA,N)
clust[which.min(y)]=1 # smallest value assigned to cluster 1
clust[which.max(y)]=2 # highest value assigned to cluster 2
dataList = list(
y = y ,
N = N ,
Nclust = Nclust ,
clust = clust ,
onesRepNclust = rep(1,Nclust)
)
Results when mean of cluster 2 is 3 standard deviations away from mean of cluster 1: The posterior recovers the generating values fairly well.
Results when mean of cluster 2 is 2 standard deviations away from mean of cluster 1: There is lots of uncertainty. See captions for discussion.
The model specification (for JAGS): The assumes that the clusters have the same standard deviation, but different means.
model {
# Likelihood:
for( i in 1 : N ) {
y[i] ~ dnorm( mu[i] , tau )
mu[i] <- muOfClust[ clust[i] ]
clust[i] ~ dcat( pClust[1:Nclust] )
}
# Prior:
tau ~ dgamma( 0.01 , 0.01 )
for ( clustIdx in 1: Nclust ) {
muOfClust[clustIdx] ~ dnorm( 0 , 1.0E-10 )
}
pClust[1:Nclust] ~ ddirch( onesRepNclust )
}
The data specification:
# Generate random data from known parameter values:
set.seed(47405)
trueM1 = 100
N1 = 200
trueM2 = 145 # 145 for first example below; 130 for second example
N2 = 200
trueSD = 15
effsz = abs( trueM2 - trueM1 ) / trueSD
y1 = rnorm( N1 )
y1 = (y1-mean(y1))/sd(y1) * trueSD + trueM1
y2 = rnorm( N2 )
y2 = (y2-mean(y2))/sd(y2) * trueSD + trueM2
y = c( y1 , y2 )
N = length(y)
# Must have at least one data point with fixed assignment
# to each cluster, otherwise some clusters will end up empty:
Nclust = 2
clust = rep(NA,N)
clust[which.min(y)]=1 # smallest value assigned to cluster 1
clust[which.max(y)]=2 # highest value assigned to cluster 2
dataList = list(
y = y ,
N = N ,
Nclust = Nclust ,
clust = clust ,
onesRepNclust = rep(1,Nclust)
)
Results when mean of cluster 2 is 3 standard deviations away from mean of cluster 1: The posterior recovers the generating values fairly well.
 |
| Upper panel: Data with underlying normal generators. Lower panel: For each datum, the posterior probability that it is assigned to cluster 2. |
 |
| Marginal posterior on cluster means and SD. |
 |
| Pairs plot of cluster means and SD. |
Results when mean of cluster 2 is 2 standard deviations away from mean of cluster 1: There is lots of uncertainty. See captions for discussion.
 |
| Lower panel: Notice that the lowest and highest data values have fixed cluster assignments, but all the other data values have posterior probabilities of cluster assignment noticeably far from 0 or 1. |
 |
| Notice the bimodal distribution of sigma (SD). |
 |
| Notice in the in the right column that when sigma is small, around 15, then the cluster means are near their true generating values. But when sigma is large, then the cluster means get close together. Essentially, there is a bimodal posterior: Either there are two clusters, with smaller sigma and distinct means, or there is one cluster, with larger sigma and both cluster means set near the mean of the one cluster. |
Friday, June 1, 2012
Beta distribution parameterized by mode instead of mean
In this post, I describe how it is easier to intuit the beta distribution in terms of its mode than its mean. This is especially handy when specifying a prior beta distribution.
(In a previous post, I explained how it is easier to intuit the gamma distribution in terms of its mode instead of its mean.)
A problem with using the mean to describe a distribution is that for skewed distributions, the mean may be far from the mode, but the mode may be what we intuitively want as the "descriptive handle" on the distribution, and therefore the mean is not a good surrogate for the description of central tendency. Especially when we are specifying a prior distribution, we may want to express our intuition in terms of the mode of the prior instead of the mean.
For a beta distribution with shape parameters a and b, the mode is (a-1)/(a+b-2). Suppose we have a desired mode, and we want to determine the corresponding shape parameters. Here's the solution. First, we express the "certainty" of the estimate in terms of the equivalent prior sample size,
Here are a few examples:
The book expressed beta distributions in terms of mean and certainty instead of mode and certainty; cf. Eqn. 5.5, p. 83, where m denoted the mean and n denoted the certainty instead of k used here.
(In a previous post, I explained how it is easier to intuit the gamma distribution in terms of its mode instead of its mean.)
A problem with using the mean to describe a distribution is that for skewed distributions, the mean may be far from the mode, but the mode may be what we intuitively want as the "descriptive handle" on the distribution, and therefore the mean is not a good surrogate for the description of central tendency. Especially when we are specifying a prior distribution, we may want to express our intuition in terms of the mode of the prior instead of the mean.
For a beta distribution with shape parameters a and b, the mode is (a-1)/(a+b-2). Suppose we have a desired mode, and we want to determine the corresponding shape parameters. Here's the solution. First, we express the "certainty" of the estimate in terms of the equivalent prior sample size,
k=a+b, with k≥2.
The certainty must be at least 2 because it essentially assumes that the prior contains at least one "head" and one "tail," which is to say that we know each outcome is at least possible. Then a little algebra reveals:
a = mode * (k-2) + 1
b = (1-mode) * (k-2) + 1
Here are a few examples:
The book expressed beta distributions in terms of mean and certainty instead of mode and certainty; cf. Eqn. 5.5, p. 83, where m denoted the mean and n denoted the certainty instead of k used here.
Sunday, May 27, 2012
Bayesian estimation supersedes the t test
Bayesian estimation for two groups provides complete distributions of credible values for the effect size, group means and their difference, standard deviations and their difference, and the normality of the data. The method handles outliers. The decision rule can accept the null value (unlike traditional t tests) when certainty in the estimate is high (unlike Bayesian model comparison using Bayes factors). The method also yields precise estimates of statistical power for various research goals. The software and programs are free, and run on Macintosh, Linux, and Windows platforms. See this linked page for the latest paper and the software.
Search tag: BEST (for Bayesian estimation)
(This is the second revision of a the first version announced at a previous post.)
Saturday, May 19, 2012
Split-Plot Design in JAGS: Revised version
A previous post reported an analysis of a "split plot" design, in which one factor is between subjects and a second factor is within subjects. That program has now been revised, and the advantage of Bayesian analysis over NHST has been confirmed.
My heartfelt thanks to Wolf Vanpaemel for pointing out an error in one of the contrast analyses in the originally posted program. I have now corrected the error and expanded the program to handle the intended contrast. Subsequent discussion with Wolf also prompted me to think harder about whether the model I specified really matched the model used by NHST analyses, and whether the apparent advantage of the Bayesian analysis (reported in the previous post) really made sense. They do, as I explain below. I also revised other aspects of the program to incorporate the improvements in the most recent forms of the other ANOVA-type programs, as reported at this previous post. (By the way, you may recall that Wolf and co-author Francis Tuerlinckx wrote a review of the book.)
A reminder of the split plot design:
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.
(More info appears at the original post.)
An advantage of Bayesian analysis:
The original post showed that the main effect of age, that is, the between-subject factor, was credibly non-zero, which agreed with the result from NHST, which had a p value a little under .05. The estimated magnitude of difference between Old and Young closely matched the marginal means in the data. But the 95% HDI on the difference was very precise compared with the NHST confidence interval. In other words, the Bayesian analysis was much more powerful. The diference was dramatic, so I expressed some doubt in the original post. But I am now convinced that the difference is real and correct. Essentially, the Bayesian analysis estimates all the parameters simultaneously, and creates a single, fixed posterior distribution over multi-dimensional parameter space. We then examine the posterior from difference perspectives, such as, for example, collapsing across all dimensions except Old minus Young, or collapsing across all dimensions except Four degrees minus Zero degrees. We do not "re-scale" the posterior for different perspectives on the parameters. In NHST, on the other hand, different tests can have different error terms in their F ratios, which means that different comparisons are scaled differently. In particular, the F ratio for Old versus Young uses a different error term (i.e., denominator) than the F ratio for Four degrees versus Zero degrees, and it turns out that the former error term is larger than the latter error term, hence the former test is not as powerful, relatively speaking. Bayes is better!
The models in the Bayesian analysis and in classical ANOVA:
In classical ANOVA, the model analyzes the data variance into five components plus residual noise. (See Maxwell & Delaney, Ch. 12, cited above.) The five components are: (1) the (between-subect) main effect, A, (2) the (within-subject) main effect, B, (3) the main effect of subject within levels of A, S/A, (4) the interaction AxB, and (5) the interaction of B with subjects within levels of A, BxS/A. It turns out that the five components actually use up all the variance, so there is zero residual noise remaining. Another way of thinking about it is that the BxS/A interaction term cannot be identified separately from the noise, because there is only one measurement per cell of the design. (That is also why the BxS/A term is used as the error term for some of the F ratios.) Therefore the model I've used in the Bayesian analysis does not include a BxS/A term:
The corrected contrast analysis:
The original post included a few different contrast analyses. One was programmed incorrectly, because I mistakenly put it in the section for interaction contrasts, instead of creating a new section for so-called "simple" contrasts. While correcting the code, I also re-expressed the contrasts in the more meaningful new format reported at this previous post. Thus, instead of just a matrix of numerical contrast coefficients without clear meaning, the new code uses logical referencing with the level names, like this:
Other revisions in the program:
I also changed the prior specification in the model, as explained at this previous post. For example, here is the prior on the main effect of A:
Thanks again to Wolf for finding the error in the contrast analysis and spurring the revisions!
Get the revised program here and the data file here.
My heartfelt thanks to Wolf Vanpaemel for pointing out an error in one of the contrast analyses in the originally posted program. I have now corrected the error and expanded the program to handle the intended contrast. Subsequent discussion with Wolf also prompted me to think harder about whether the model I specified really matched the model used by NHST analyses, and whether the apparent advantage of the Bayesian analysis (reported in the previous post) really made sense. They do, as I explain below. I also revised other aspects of the program to incorporate the improvements in the most recent forms of the other ANOVA-type programs, as reported at this previous post. (By the way, you may recall that Wolf and co-author Francis Tuerlinckx wrote a review of the book.)
A reminder of the split plot design:
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) |
An advantage of Bayesian analysis:
The original post showed that the main effect of age, that is, the between-subject factor, was credibly non-zero, which agreed with the result from NHST, which had a p value a little under .05. The estimated magnitude of difference between Old and Young closely matched the marginal means in the data. But the 95% HDI on the difference was very precise compared with the NHST confidence interval. In other words, the Bayesian analysis was much more powerful. The diference was dramatic, so I expressed some doubt in the original post. But I am now convinced that the difference is real and correct. Essentially, the Bayesian analysis estimates all the parameters simultaneously, and creates a single, fixed posterior distribution over multi-dimensional parameter space. We then examine the posterior from difference perspectives, such as, for example, collapsing across all dimensions except Old minus Young, or collapsing across all dimensions except Four degrees minus Zero degrees. We do not "re-scale" the posterior for different perspectives on the parameters. In NHST, on the other hand, different tests can have different error terms in their F ratios, which means that different comparisons are scaled differently. In particular, the F ratio for Old versus Young uses a different error term (i.e., denominator) than the F ratio for Four degrees versus Zero degrees, and it turns out that the former error term is larger than the latter error term, hence the former test is not as powerful, relatively speaking. Bayes is better!
The models in the Bayesian analysis and in classical ANOVA:
In classical ANOVA, the model analyzes the data variance into five components plus residual noise. (See Maxwell & Delaney, Ch. 12, cited above.) The five components are: (1) the (between-subect) main effect, A, (2) the (within-subject) main effect, B, (3) the main effect of subject within levels of A, S/A, (4) the interaction AxB, and (5) the interaction of B with subjects within levels of A, BxS/A. It turns out that the five components actually use up all the variance, so there is zero residual noise remaining. Another way of thinking about it is that the BxS/A interaction term cannot be identified separately from the noise, because there is only one measurement per cell of the design. (That is also why the BxS/A term is used as the error term for some of the F ratios.) Therefore the model I've used in the Bayesian analysis does not include a BxS/A term:
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 model has no BxS term because that would leave zero noise variance.
}
The model remains unchanged from the original program; the news here is explicitly explaining it.
The original post included a few different contrast analyses. One was programmed incorrectly, because I mistakenly put it in the section for interaction contrasts, instead of creating a new section for so-called "simple" contrasts. While correcting the code, I also re-expressed the contrasts in the more meaningful new format reported at this previous post. Thus, instead of just a matrix of numerical contrast coefficients without clear meaning, the new code uses logical referencing with the level names, like this:
simpleContrastList = list(
A1B1vA2B1 = ( outer(aLvlNames=="A2(Old)",bLvlNames=="B1(Zero)")
- outer(aLvlNames=="A1(Young)",bLvlNames=="B1(Zero)") )
)
That code simply means to compute the difference of "A2(Old),B1(Zero)" and "A1(Young),B1(Zero)", that is, the difference of Old and Young at Zero degrees angle. If you haven't previously tried programming a contrast, that code might still look mysterious, but once you get the hang of it, it's more meaningful and resistant to inadvertent re-ordering of levels than the old format: simpleContrastList = list(
A1B1vA2B1 = matrix( c(1,0,0,-1,0,0) , nrow=2 , byrow=TRUE )
)
While that old form is shorter, it is less meaningful because it does not show the names of the levels being referred to, and it is easily broken if the levels get re-ordered by different operating system conventions for alphabetizing. One more change: The computation of the simple contrast, at the end of the program, relies on a new section of code.Other revisions in the program:
I also changed the prior specification in the model, as explained at this previous post. For example, here is the prior on the main effect of A:
for ( j in 1:NaLvl ) { a[j] ~ dnorm( 0.0 , aTau ) }
aTau <- 1 / pow( aSD , 2 )
aSD ~ dgamma(1.221,0.003683) # mode=60,sd=300. Change for scale of data!
The estimate of the standard deviation of the A effect, aSD, has a gamma prior that has a mode at 60 with declining density toward zero, instead of a "generic" uniform or folded-t or gamma that has notable mass near zero.Thanks again to Wolf for finding the error in the contrast analysis and spurring the revisions!
Get the revised program here and the data file here.
Subscribe to:
Posts (Atom)