24.1.1 Data structure

Table 24.1 shows counts of hair color and eye color from a classroom poll of students at the University of Delaware (Snee, 1974). These are the same data as in Table 4.1 (p. 90) but with the original counts instead of proportions, and with levels re-ordered alphabetically. Respondents self-reported their hair color and eye color, using the labels shown in Table 24.1. The cells of the table indicate the frequency with which each combination occurred in the sample. Each respondent fell in one and only one cell of the table. The data to be predicted are the cell counts. The predictors are the nominal variables. This structure is analogous to two-way analysis of variance (ANOVA), which also had two nominal predictors, but had several metric values in each cell instead of a single count.

For data like these, we can ask a number of questions. We could wonder about one predictor at a time, and ask questions such as, “How much more frequent are brown-eyed people than hazel-eyed people?” or “How much more frequent are black-haired people than blond-haired people?” Those questions are analogous to main effects in ANOVA. But usually we display the joint counts specifically because we're interested in the relationship between the variables. We would like to know if the distribution of counts across one predictor is contingent upon the level of the other predictor. For example, does the distribution of hair colors depend on eye color, and, specifically, is the proportion of blond-haired people the same for brown-eyed people as for blue-eyed people? These questions are analogous to interaction contrasts in ANOVA.

24.1.2 Exponential link function

The model can be motivated two different ways. One way is simply to start with the two-way ANOVA model for combining nominal predictors, and find a way to map the predicted value to count data. The predicted μ from the ANOVA model (recall Equations 20.1 and 20.2, p. 584) could be any value from negative to positive infinity, but frequencies are non-negative. Therefore, we must transform the ANOVA predictions to non-negative values, while preserving order. The natural way to do this, mathematically, is with the exponential transformation. But this transformation only gets us to an underlying continuous predicted value, not to the probability of discrete counts. A natural candidate for the needed likelihood distribution is the Poisson (described later), which takes a non-negative value λ and gives a probability for each integer from zero to infinity. But this motivation may seem a bit arbitrary, even if there's nothing wrong with it in principle.

A different motivation starts by treating the cell counts as representative of underlying cell probabilities, and then asking whether the two nominal variables contribute independent influences to the cell probabilities. For example, in Table 24.1, there's a particular marginal probability that hair color is black, and a particular marginal probability that eye color is brown. If hair color and eye color are independent, then the joint probability of black hair and brown eyes is the product of the marginal probabilities. The attributes of hair color and eye color are independent if that relationship holds for every cell in the table. ($\mathrm{R}$ecall the definition of independence in Section 4.4.2, p. 92.) Independence of hair color and eye color means that the proportion of black hair among brown-eyed people is the same as the proportion of black hair among blue-eyed people, and so on for all hair and eye colors.

To check for independence of attributes, we need to estimate the marginal probabilities of the attributes. Denote the marginal (i.e., total) count of the rth row as y_r, and denote the marginal count of the cth column as y_c. Then the marginal proportions are y_r/N and y_c/N, where N is the total of the entire table. If the attributes are independent, then the predicted joint probability, $\widehat{p}\left(r,c\right)$, should equal the product of the marginal probabilities, which means $\widehat{p}\left(r,c\right)=p\left(r\right)\cdot p\left(c\right)$, hence ${\widehat{y}}_{r,c}/N=y_r/N\cdot y_c/N$. Because the models we've been dealing with involve additive combinations, not multiplicative combinations, we convert the multiplicative expression of independence into an additive expression by using the facts that log(a · b) = log(a) + log(b) and exp(log(x)) = x, as follows:

$\begin{array}{ccc}\hfill {\widehat{y}}_{r,c}/N& =& y_r/N\cdot y_c/N\hfill \\ \hfill {\widehat{y}}_{r,c}& =& 1/N\cdot y_r\cdot y_c\hfill \\ \hfill \underset{{\text{λ}}_{r,c}}{\underbrace{{\widehat{y}}_{r,c}}}& =& \exp \left(\underset{{\beta }_0}{\underbrace{\log \left(1/N\right)}}+\underset{{\beta }_r}{\underbrace{\log \left(y_r\right)}}+\underset{{\beta }_c}{\underbrace{\log \left(y_c\right)}}\right)\hfill \end{array}$

If we abstract the form of Equation 24.1 away from the specific counts, we get the equation λ_r,c = exp + (β₀ + β_r + β_c). The idea is that whatever are the values of the β's, the resulting λ's will obey multiplicative independence. An example is shown in Table 24.2. The choice of β's is shown in the margins of the table, and the resulting λ's are shown in the cells of the table. Notice that every row has the same relative probabilities, namely, 10, 100, and 1. In other words, the row and column attributes are independent.

We have dealt before with additive combinations of row and column influences, in the context of ANOVA. In ANOVA, β₀ is a baseline representing the overall central tendency of the outcomes, and β_r, is a deflection away from baseline due to being in the rth row, and β_c is a deflection away from baseline due to being in the cth column. The deflections are constrained to sum to zero, as required by Equation 20.2 (p. 584), and the example in Table 24.2 respects this constraint.

In ANOVA, when the cell mean is not captured by an additive combination of row and column effects, we include interaction terms, denoted here as β_{r, c}. The interaction terms are constrained so that every row and every column sums to zero: ${\displaystyle {\sum }_r{\beta }_{r,c}=0}$ for all c and ${\displaystyle {\sum }_c{\beta }_{r,c}=0}$ for all r. The key idea is that the interaction terms in the model, which indicate violations of additivity in standard ANOVA, indicate violations of multiplicative independence in exponentiated ANOVA. To summarize, the model of the cell tendencies is Equation 24.1 extended with interaction terms:

${\text{λ}}_{r,c}=\exp \left({\beta }_0+{\beta }_r+{\beta }_c+{\beta }_{r,c}\right)$

with the constraints $\begin{array}{ccc}{\displaystyle \sum _r{\beta }_r=0,}& {\displaystyle \sum _c{\beta }_c=0,}& {\displaystyle \sum _r{\beta }_{r,c}=0}\end{array}$ for all c, and ${\displaystyle \sum _c{\beta }_{r,c}=0}$ for all rIf the researcher is interested in violations of independence, then the interest is on the magnitudes of the β_{r, c} interaction terms, and specifically on meaningful interaction contrasts. The model is especially convenient for this purpose, because specific interaction contrasts can be investigated to determine in more detail where the nonindependence is arising.

24.1.3 Poisson noise distribution

The value of λ_{r, c} in Equation 24.2 is a cell tendency, not a predicted count per se. In particular, the value of λ_{r, c} can be any non-negative real value, but counts can only be integers. What we need is a likelihood function that maps the parameter value λ_{r, c} to the probabilities of possible counts. The Poisson distribution is a natural choice. The Poisson distribution is named after the French mathematician Simon-Denis Poisson (1781–1840) and is defined as

$p\left(y|\text{λ}\right)={\text{λ}}^y\exp \left(-\text{λ}\right)/y!$

where y is a non-negative integer and λ is a non-negative real number. The mean of the Poisson distribution is λ. Importantly, the variance of the Poisson distribution is also λ (i.e., the standard deviation is $\sqrt{\text{λ}}$). Thus, in the Poisson distribution, the variance is completely yoked to the mean. Examples of the Poisson distribution are shown in Figure 24.1. Notice that the distribution is discrete, having masses only at non-negative integer values. Notice that the visual central tendency of the distribution does indeed correspond with λ. And notice that as λ increases, the width of the distribution also increases. The examples in Figure 24.1 use noninteger values of λ to emphasize that λ is not necessarily an integer, even though y is an integer.

The Poisson distribution is often used to model discrete occurrences in time (or across space) when the probability of occurrence is the same at any moment in time (e.g., Sadiku & Tofighi, 1999). For example, suppose that customers arrive at a retail store at an average rate of 35 people per hour. Then the Poisson distribution, with λ = 35, is a model of the probability that any particular number of people will arrive in an hour. As another example, suppose that 11.2% of the students at the University of Delaware in the early 1970s had black hair and brown eyes, and suppose that an average of 600 students per term will fill out a survey. That means an average of 67.2 (= 11.2% · 600) students per term will indicate they have black hair and brown eyes. The Poisson distribution, p(y| λ = 67.2), gives the probability that any particular number of people will give that response in a term.

We will use the Poisson distribution as the likelihood function for modeling the probability of the observed count, y_r,c, given the mean, λ_r,c, from Equation 24.2. The idea is that each particular r, c combination has an underlying average rate of occurrence, λ_r,c. We collect data for a certain period of time, during which we happen to observe particular frequencies, y_r,c of each combination. From the observed frequencies, we infer the underlying average rates.

24.1.4 The complete model and implementation in JAGS

In summary, the preceding discussion says that the predicted tendency of the cell counts is given by Equation 24.2 (including the sum-to-zero constraints), and the probability of the observed cell count is given by the Poisson distribution in Equation 24.3. All we have yet to do is provide prior distributions for the parameters. Fortunately, we have already done this (at least structurally) for two-factor “ANOVA.” The resulting model is displayed in Figure 24.2. The top part shows the prior from two-factor “ANOVA,” reproduced from Figure 20.2 (p. 588). The mathematical expression in the middle of Figure 24.2 is the same as in Figure 20.2 except for the new exponential inverse-link function. The lower part of Figure 24.2 simply illustrates the Poisson noise distribution.

The model is implemented in JAGS as a variation on the model for two-factor “ANOVA.” As in that model, the sum-to-zero constraint is imposed by re-centering the deflection parameters that are found by MCMC without the constraint. The deflection parameters that are directly sampled by the MCMC process are denoted a1, a2, and a1a2, and the sum-to-zero versions are denoted b1, b2, and b1b2. To understand the following JAGS code, it is helpful to understand that the data file has one row per cell, such that the ith line of the data file contains the count y[i] from table row x1[i] and column x2[i]. The first part of the JAGS model specification directly implements the Poisson noise distribution and exponential inverse-link function:

model {
 for ( i in 1:Ncell ) {
 y[i] ~ dpois( lambda[i] )
 lambda[i] <- exp( a0 + a1[x1[i]] + a2[x2[i]] + a1a2[x1[i],x2[i]] )
}

The next part of the JAGS model specifies the prior distribution on the baseline and deflection parameters. The prior is supposed to be broad on the scale of the data, but we must be careful about what scale is being modeled by the baseline and deflections. The counts are being directly described by λ, but it is log(λ) being described by the baseline and deflections. Thus, the prior on the baseline and deflections should be broad on the scale of the logarithm of the data. To establish a generic baseline, consider that if the data points were distributed equally among the cells, the mean count would be the total count divided by the number of cells. The biggest possible standard deviation across cells would occur when all the counts were loaded into a single cell and all the other cells were zero. We take the logarithms of those values to define the following constants:

yLogMean = log( sum(y) / (Nx1Lvl*Nx2Lvl) )
yLogSD = log( sd( c( rep(0,Ncell-1), sum(y) ) ) )
agammaShRa = unlist( gammaShRaFromModeSD( mode=yLogSD, sd=2*yLogSD ) )

Those constants are used to set broad priors, as follows:

a0 ~ dnorm( yLogMean , 1/(yLogSD*2)^2 )
for ( j1 in 1:Nx1Lvl ) { a1[j1] ~ dnorm( 0.0 , 1/a1SD^2 ) }
a1SD ~ dgamma (agammaShRa[1],agammaShRa[2])
for ( j2 in 1:Nx2Lvl ) { a2[j2] ~ dnorm(0.0 , 1/a2SD^2 ) }
a2SD ~ dgamma(agammaShRa[1],agammaShRa[2])
for ( j1 in 1:Nx1Lvl ) { for ( j2 in 1:Nx2Lvl ) {
 a1a2[j1,j2] ~ dnorm(0.0 , 1/a1a2SD^2 )
} }
a1a2SD ~ dgamma(agammaShRa[1],agammaShRa[2])

Near the end of the JAGS model specification, the baseline and deflections are recentered to respect the sum-to-zero constraints, exactly as was done for two-factor “ANOVA”:

# Convert a0,a1[],a2[],a1a2[,] to sum-to-zero b0,b1[],b2[],b1b2[,] :
for ( j1 in 1:Nx1Lvl ) { for ( j2 in 1:Nx2Lvl ) {
. m[j1,j2] <- a0 + a1[j1] + a2[j2] + a1a2[j1,j2] # cell means
} }
b0 <- mean( m[1:Nx1Lvl,1:Nx2Lvl] )
for ( j1 in 1:Nx1Lvl ) { b1[j1] <- mean( m[j1,1:Nx2Lvl] ) - b0 }
for ( j2 in 1:Nx2Lvl ) { b2[j2] <- mean( m[1:Nx1Lvl,j2] ) - b0 }
for ( j1 in 1:Nx1Lvl ) { for ( j2 in 1:Nx2Lvl ) {
 b1b2[j1,j2] <- m[j1,j2] - ( b0 + b1[j1] + b2[j2] )
} }

Finally, a novel section of the model specification computes posterior predicted cell probabilities. Each cell mean (which was computed when re-centering for sum to zero) is exponentiated to compute the predicted count for each cell, and then the counts are normalized to compute the predicted proportion in each cell, denoted ppx1x2p. The predicted marginal proportions are also computed, as follows:

 # Compute predicted proportions:
 for ( j1 in 1:Nx1Lvl ) { for ( j2 in 1:Nx2Lvl ) {
 expm[j1,j2] <- exp(m[j1,j2])
 ppx1x2p[j1,j2] <- expm[j1,j2]/sum(expm[1:Nx1Lvl,1:Nx2Lvl])
 } }
 for ( j1 in 1:Nx1Lvl ) { ppx1p[j1] <- sum(ppx1x2p[j1,1:Nx2Lvl]) }
 for ( j2 in 1:Nx2Lvl ) { ppx2p[j2] <- sum(ppx1x2p[1:Nx1Lvl,j2]) }
}

The model specification and other functions are defined in the file Jags-Ycount-Xnom2fac-MpoissonExp.R, and a high-level script that calls the functions is the file Jags-Ycount-Xnom2fac-MpoissonExp-Example.R.