Analyzing count variables in individuals and groups: Single level and multilevel models

Abstract

We address two challenges in data analysis of group research. First, data may be clustered (i.e., responses of individual group members are correlated). Second, some dependent variables may consist of integer counts of number of occurrences of an event. Familiar ANOVA and regression models provide nonoptimal analyses in both cases. Standard multilevel (mixed) models yield accurate inference for clustered normally distributed data. Generalized linear models (GLMs), specifically Poisson regression and related models, yield accurate inference for nonclustered count data. New generalized linear mixed models (GLMMs) integrate GLMs with multilevel models, addressing both challenges and yielding accurate inferences for grouped count outcomes. To provide the necessary background for understanding GLMMs, we first introduce GLMs, with detailed coverage in an example of Poisson regression. We then introduce multilevel models. Finally, we develop GLMMs and illustrate in an example their application to clustered count data. Group research may benefit from the flexibility provided by GLMMs.

Keywords

clustered data count data generalized linear mixed models generalized linear models multilevel models

Measures of behavior in psychology and the behavioral sciences often take the form of count variables. Counts sum the number of discrete events in a fixed time period (e.g., 24 hours) and range upward from a minimum of zero. Traditional approaches to the analysis of count variables such as analysis of variance and multiple regression often do not give proper results.

A variety of specialized approaches to count data that provide better answers have been developed, but they are not well known to most psychologists. In addition, group researchers working with data that are clustered face the additional challenge that the responses of those individuals are not independent—the responses of members of the same group are likely to be more similar than those of members of different groups. A new class of statistical models known collectively as the generalized linear mixed model (GLMM) permit researchers to address both issues simultaneously. These new statistical models provide optimal answers both to questions traditionally posed by group researchers and to novel questions that are facilitated by thinking in terms of the new models. The purpose of this article is to initially provide the necessary foundation for understanding statistical models for (a) count data and (b) for clustered (grouped) data. We then present (c) the basics of the new GLMM models. We present two examples of the analysis of count data. The first illustrates the analysis and interpretation of the results of a hypothetical study using individual participants, and the second illustrates the analysis of individual participants who are participating as part of a group. Computer syntax in two common packages (SAS, SPSS) for conducting the analyses of the two examples is provided in the online supplementary material. We also point to more advanced sources on these models for additional information. Our hope is that this article will provide researchers with the basic information necessary to conduct appropriate analyses of count variables in their own research on individuals and groups.

Counts and Their Analysis

Count measures frequently have been used in psychology, ranging from assessments of social network size (Rossotti, Winter, & Watts, 2006) to characterizations of emotional expression (Tov, Ng, Lin, & Qiu, 2013). The methods of analysis of count data for individuals are now well established with chapters (Long, 1997, Chapter 8), didactic articles (Coxe, West, & Aiken, 2009), tutorials in substantive areas (Atkins, Baldwin, Zheng, Gallop, & Neighbors, 2013; Atkins & Gallop, 2007; Berk & McDonald, 2008; Xie, Tao, McHugo, & Drake, 2013) and entire textbooks (Cameron & Trivedi, 2013; Hilbe, 2011, 2014) devoted to the issue. The analysis of clustered data with continuous variables has similarly become well established (Cohen, Cohen, West, & Aiken, 2003, Chapter 13; Hox, 2010; Raudenbush & Bryk, 2002; Snijders & Bosker, 2012). In contrast, the GLMM models that address count data that are clustered (e.g., groups) are a much newer development. Only in the past few years have the statistical models and the software for GLMM become available that can address both issues. Hox (2010) and Snijders and Bosker (2012) have included brief treatments of clustered count outcomes in the most recent versions of their textbooks (see also Hedeker & Gibbons, 2006, Chapter 12). Stroup (2013) has provided an advanced treatment of the analysis of count data using multilevel or mixed models when the data structure contains individual cases that are clustered within groups, or repeated observations of the same individuals (e.g., in daily diary studies).

Given extensive data in the form of count variables, as well as the recent and ongoing statistical developments in the analysis of counts, the first purpose of this article is to provide an overview of the most commonly employed methods for regression analysis of single-level count data, in which participants are independent of one another. The second purpose is to extend the overview to the more recent development of analysis of count data in the mixed model or multilevel framework, in which participants are clustered into groups.

Overview of the Article

We begin our presentation with a consideration of the difficulties in treating count outcomes with data from individuals using familiar ordinary least squares regression (OLS). The regression models for the analysis of count data are examples of a broad collection of regression models known as the generalized linear model (GLM; see Coxe, West, & Aiken, 2013). The GLM subsumes specific regression models for the analysis of many forms of dependent variables other than continuous measures, including binary outcomes, unordered and ordered categories, and counts. We provide a very brief overview of the characteristics of all generalized linear models (GLMs) to set the work on count variables in context; this overview will apply directly to both single-level and multilevel models. We then introduce Poisson regression, the most common analysis applied to count data (Cameron & Trivedi, 2013). We briefly explain three alternatives, quasi-Poisson regression, negative binomial regression, and robust estimation, employed when a fundamental assumption of the Poisson regression model is not met. We consider various aspects of count regression models and their interpretation that diverge from familiar OLS regression. We consider the inherent multiplicative nature of count models, the interpretation of regression coefficients, the nature of interactions, the graphical display of count regression models, the assessment of goodness of fit and gain in prediction. We then provide a brief presentation of key concepts from hierarchical linear (mixed) models with continuous outcomes. With these two foundations, we move to the treatment of count outcomes with clustered data in the multilevel (or mixed model) framework, introducing a new class of models, generalized linear mixed models (GLMMs). We provide an overview of the complexities of GLMMs in both estimation and inference. To make the analyses concrete, we provide illustrations of the analysis of single level and multilevel count data carried out in SPSS and SAS on simulated data sets.

Generalized Linear Models: Analysis of Single Level Count Data

Issues With Traditional Analyses

Count variables do not hold special status as predictors in regression analysis. In contrast, substantial challenges can arise when count data serve as the dependent variable in familiar ordinary least squares (OLS) regression analysis.

OLS regression models the linear relationships of predictors to the criterion, and counts are often nonlinear functions of the predictors.

Count variables frequently have many cases with 0s and other low values; the predicted values for these cases are often negative, outside the bounds of potential values for counts.

Count dependent variables are often skewed. This skewness is often propagated through to the residuals. Nonnormal residuals can lead to biased standard errors, which in turn, bias significance tests and confidence intervals. For statistical inference in OLS regression, normally distributed residuals are assumed.

The variance of a count increases with its mean leading to heteroscedasticity. A distribution with an arithmetic mean count of 1.0 can have little variation in individual counts; with a mean of 5.0, there is much more opportunity for variation among individual counts. Again, OLS estimates of standard errors will be biased (too small), leading to overestimates of predictor significance and to confidence intervals that are too narrow (Kmenta, 1986, p. 278).

Logarithmic transformation of positively skewed dependent variables has long been recommended to reduce or eliminate positive skew (see, e.g., Guilford, 1956, almost 60 years ago). Zero (0) is a frequent count, and the logarithm of zero is minus infinity. One long-standing practice to bypass this problem is to use started logs (Tukey, 1977, p. 396) in which small constant like .1 is added to each score and the logarithm of the resulting value is taken. This practice does not eliminate the pileup of very low scores. Further, treating the logarithm of a count variable as a dependent variable in an OLS regression entirely eliminates neither the nonnormality of residuals nor their troublesome heteroscedasticity. Started logs can still lead to predicted count scores that are below zero (O’Hara & Kotze, 2010). Mathematically, OLS regression coefficients of logged counts have been shown to be biased and inconsistent (i.e., do not approach the corresponding population value with increasing sample size; King, 1988). OLS regression with logged counts is generally not recommended, given the availability of the appropriate and superior statistical models addressed in what follows.

Generalized Linear Models

Statistical models for count data are special cases of the generalized linear model (GLM; Cohen et al., 2003, Chapter 13; Coxe et al., 2013). The GLM is a family of highly flexible regression models that are able to handle a wide range of forms of dependent variables. These dependent variables include continuous, binary, categorical and ordered categorical, time to failure (survival data), and count outcomes, among others. OLS regression (and ANOVA) are also special cases of the GLM.

There are three components to each specific regression model within the GLM. First is the random component (or error structure). This is a theoretical probability distribution that closely resembles the conditional distribution of the dependent variable given values on the predictor(s) (i.e., the distribution of the residuals in the regression analysis). The selected theoretical probability distribution should closely reflect the shape of the conditional distribution and also should reflect the relationship between the mean and variance of the conditional distribution. In OLS regression, the normal distribution is the theoretical distribution; in a normal distribution the mean and variance are independent of one another. The variance of the predicted scores is unrelated to the predicted score (homoscedasticity). In contrast, the mean and variance of count variables are not independent. Counts are typically positively skewed, and the variance of a count measure typically increases with its mean.

The second component is the systematic linear predictor. For all GLMs (e.g., logistic regression for binary outcomes) the regression equation can be expressed in linear form:

f (\hat{Y}) = b_{1} X_{1} + b_{2} X_{2} + \dots + b_{p} X_{p} + b_{0},

(1.1)

where $f (\hat{Y})$ represents the predicted value. In the linear form the predicted score is some function of the predicted score that would arise in the units of the observed score. For counts, the predicted score in the linear equation is the natural logarithm of the predicted count,

[ [ $f (\hat{Y}) = \ln (\hat{Y})]$ . In the special case of OLS regression, the linear predictor is in the same units as the actual dependent variable [ [ $f (\hat{Y}) = \hat{Y}$ ] because the predicted scores in OLS regression are in the same units as the original dependent variable scores.

The third component of each GLM is the link function. The link function relates the units of the linear predictor (i.e., predicted score) to the units of the observed criterion score. In OLS regression, if we predict children’s weight in pounds from their food intake (X₁) and height (X₂) using the linear regression equation $\hat{Y} = b_{1} X_{1} + b_{2} X_{2} + b_{0}$ , both the observed value Y and the predicted value $\hat{Y}$ are in the same units (pounds). Thus, for OLS regression the link function is the identity function—the predicted and observed scores are in the same units. For count outcome variables, the linear predictor form predicts the natural logarithm of the predicted count, $l n (\hat{Y}) = b_{1} X_{1} + b_{2} X_{2} + b_{0}$ , whereas the outcome is measured as a raw count (0, 1, 2, …). For count variables the units in which the predicted score is expressed are a transformation of the units of the observed criterion score, the natural logarithm of the count dependent variable. Thus, the link function is the natural log.

Poisson Regression for Count Dependent Variables

The most basic model for the prediction of count dependent variables, Poisson regression, includes each of these three GLM components. The first component is the random component (or error structure). The most fundamental probability distribution (or “benchmark parametric model”; Cameron & Trivedi, 2013, p. 3) for characterizing counts is the Poisson probability distribution (named for French mathematician Siméon Denis Poisson). The Poisson distribution describes the distribution of the number of independent events that are generated in a particular time period t, for example, the number of aggressive acts exhibited by each of a number of children during a 30-minute school recess. The number of events generated depends on a rate (mean) parameter μ, which characterizes the mean number of events, here, mean aggressive acts averaged across all the children. Assuming that the time period is constant for all cases, the expression for the Poisson probability (Pr) of the occurrence of a particular number of events y is given as follows (Cameron & Trivedi, 2013, p. 3):

\Pr (Y = y) = \frac{e^{- μ} μ^{y}}{y!} for y = 0, 1, 2 \dots .

(1.2)

where Y is the count variable and y is a specific numerical count (0, 1, 2, …). The Poisson probability distribution has the special property that the mean and the variance of the distribution are equal. The mean and the variance both equal μ, the rate parameter, that is $σ^{2} = μ$ . This condition is referred to as equidispersion. A figure illustrating the Poisson distribution with varying rate parameters is given in Coxe et al. (2009, p. 123). There are two noteworthy characteristics of Poisson distributions: (a) the range of expected counts in the distribution increases with an increasing mean μ, and (b) the skew decreases with increasing mean μ.

The second component is the linear predictor. It is the linear form of the Poisson regression equation in which the predicted score is a linear combination of the predictors, as shown in Equation (1.3):

\ln (\hat{C o u n t)} = b_{1} X_{1} + b_{2} X_{2} + \dots + b_{p} X_{p} + b_{0}

(1.3)

This linear regression equation has the same familiar form as OLS regression. However, there are two fundamental differences: (a) The Poisson distribution and not the normal distribution is used for statistical inference; (b) the predicted score is the natural logarithm of the predicted count, not the count itself.

The third component of Poisson GLM is the link function, the transformation that relates the units of the predicted score in the linear predictor, $\ln (\hat{C o u n t)}$ to the original count dependent variable units. The link function in Poisson regression is the natural logarithm.

Two Forms of the Poisson Regression Equation and Interpretation of Regression Coefficients

Every regression model in the GLM family can be written in two forms. One form is the linear form, as in Equation (1.3), in which the predicted score (linear predictor) is a linear combination of the predictors and in which the predicted score is a nonlinear function of the observed dependent variable (for Poisson, the natural logarithm). The second form is a form in which the predicted score is in the units of the measured outcome. The Poisson regression equation may be written in exponentiated form, such that the predicted score is in actual count units:

\hat{C o u n t} = e^{b_{1} X_{1} + b_{2} X_{2} + \dots + b_{p} X_{p} + b_{0}}, or equivalently,

(1.4a)

\hat{C o u n t} = e^{b_{1} X_{1}} e^{b_{2} X_{2}} \dots e^{b_{p} X_{p}} e^{b_{0}} .

(1.4b)

The form of Equation (1.4b) is useful because it allows researchers to make predictions in the metric of the original scores, here counts, which can sometimes be easier to interpret. Equation (1.4b) also reveals a fundamental characteristic of Poisson regression. The predicted counts are a multiplicative function of the predictors, not an additive function as in OLS regression. Poisson regression is thus a nonlinear form of regression analysis. We will see the regression coefficients in Poisson regression given in both forms, as $b_{j}$ as in Equation (1.3) and as $e^{b_{j}}$ for predictor $X_{j}$ as in Equations (1.4a) and (1.4b). The form $b_{j}$ is the amount added to $\ln (\hat{C o u n t)}$ the logarithm of the predicated count, for each one point increase in predictor $X_{j}$ , as in Equation (1.3). The form $e^{b_{j}}$ is the amount by which $\hat{C o u n t}$ , the predicted count, is multiplied for each one point increase in predictor $X_{j}$ , as is clear in Equations (1.4a) and (1.4b).

Hypothetical Example: In-Group Out-Group Bias

To illustrate Poisson regression, we present the data of a fictitious experiment addressing in-group out-group bias (see e.g., Brewer, 1979). Participants are members of an in-group, students attending an excellent public university (Purdue University, the first author’s PhD alma mater). In all, 500 Purdue students provide ratings of strength of their identification with their university on a 1 to 7 scale, where 7 = strongest identification. They then review an advertising campaign video for a second public university. The campaign stresses the virtues and excellence of the second university—saying that it is the flagship university of the state, with internationally known faculty, outstanding research centers, and exceptionally qualified students who are successful following graduation. Competition between the in-group (Purdue University students) and the other university is manipulated. In a high competition condition (n = 250), the second university is the rival flagship public university (Indiana University) in the same state. Alternatively, in the low competition condition (n = 250), the comparable university is one with which there is no rivalry, located in another region of the country (University of Virginia). Students are randomly assigned to condition (high vs. low competition) after having provided their ratings of identification with their own university. Once they view the advertising campaign video, they write a 100-word narrative of their overall impressions of the other university. Students are run individually with no opportunity to interact with other participants throughout the experiment. The dependent variable is the count of the number of derogatory remarks students include in their narratives about the other university. The two independent variables are identification with one’s own university and the competition condition.

Figure 1 shows the observed distribution of counts of disparaging remarks per participant over all 500 participants, and within each competition condition. The distribution of counts in a Poisson distribution with the same arithmetic mean number of counts as in the observed distribution is also displayed. The Poisson distribution is a good, though not perfect approximation to the data. As seen in Figure 1a, there are a greater number of smaller counts (zero or one derogatory remark) and a greater number of higher counts (six to eight derogatory remarks) compared to the Poisson distribution. A comparison of Figures 1b versus 1c shows that this deviation in number of counts from the Poisson distribution is attributable to the manipulation of competition.

Figure 1.

Distribution of actual counts of derogatory remarks versus counts expected from a Poisson distribution with the same arithmetic mean counts as the observed distribution. Solid lines represent the data; dashed lines represent the values predicted by the model.

Table 1 provides the means and variances of number of disparaging remarks overall and in each condition as a function of identification and condition. There are four important observations of note. First, the mean number of disparaging remarks increases as identification increases. Second, the within cell variance increases as the mean increases, characteristic of the heteroscedasticity typically observed in count data and in violation of a fundamental assumption of OLS regression. Third, the arithmetic mean count in the high competition condition is over twice that of the low competition condition ( $\bar{X}$ = 3.02 vs. 1.31, respectively), indicating that the high competition manipulation increased the number of derogatory remarks. Fourth, both overall and in each competition condition, the variance of the number of counts exceeds the mean number of counts (e.g., overall mean of 2.17 vs. overall variance of 3.91). This situation in which the conditional variance exceeds the conditional mean contrasts with the theoretical Poisson probability distribution in which the mean and variance are equal; this situation is known as overdispersion. For now, we focus on the analysis of the simulated data with basic Poisson regression, ignoring the overdispersion, understanding that test statistics will be positively biased. We do this to focus on the interpretation of key features of the basic analysis—the structure of overall tests of significance of prediction of the full model, tests of the contribution of individual predictors, assessment of goodness of fit, interpretation of individual coefficients, and interactions between predictors. Following our consideration of key features of basic Poisson regression, we return to the issue of overdispersion.

Table 1.

Mean and variance of number of disparaging remarks in narratives as a function of identification with one’s university and competition with the second university (fictitious data). Means with standard deviations in parentheses and number of cases.

Identification level	Full sample		Low competition		High competition
	(n = 500 participants)		(n = 250 participants)		(n = 250 participants)
	Mean	Variance	Mean	Variance	Mean	Variance
	(Number of cases)		(Number of cases)		(Number of cases)
1 = lowest	1.08	1.16	0.95	0.87	1.24	1.50
	(73)		(40)		(33)
2	1.18	1.27	1.05	1.28	1.32	1.26
	(74)		(40)		(34)
3	1.60	1.44	1.18	0.75	1.90	1.73
	68		(28)		(40)
4	2.53	4.00	1.39	1.58	3.37	4.19
	(66)		(28)		(38)
5	2.64	2.96	1.71	1.42	3.85	2.37
	(78)		(44)		(34)
6	2.90	5.10	1.33	2.51	4.22	3.49
	(66)		(30)		(36)
7 = highest	3.24	7.16	1.53	1.84	5.20	6.05
	(75)		(40)		(35)
Overall	2.17	3.91	1.31	1.48	3.02	4.89
	(500)		(250)		(250)

Estimation, Tests of Model Fit and of Individual Predictors, and Fit Indices GLMs

Estimation

The parameters of GLMs (i.e., regression coefficients and regression intercept) including Poisson regression are estimated using maximum likelihood estimation (MLE) rather than ordinary least squares estimation used in OLSregression. In MLE parameter estimates are selected that maximize the likelihood of the observed relationships in the data. Each regression model has a maximum likelihood function, a mathematical expression that describes the likelihood (similar to probability) that the observed data could be produced by specific values of the parameters in the model. The estimates of the regression coefficients and intercept are selected that maximize the likelihood function (see Enders, 2005, for a detailed overview). In the case of OLS regression, maximum likelihood estimates and OLS estimates of parameters are identical, again bringing OLS regression into the GLM family.

Overall Model Fit and Contribution of Individual Predictors

The measure of overall model fit in GLMs differs from the familiar squared multiple correlation $(R^{2})$ from OLS regression. In GLMs including Poisson regression the measure of overall model fit is termed deviance. In fact, in OLS regression, $S S_{r e s i d u a l}$ is identical to deviance. Deviance is a measure of the lack of fit (or “badness” of fit) of a GLM in comparison to a perfectly fitting model, which is termed a fully saturated model. Evaluation focuses on the reduction in deviance from a less complete model to a more complete model (more predictors). The models being compared must be nested models—that is, the less complete model must contain only a subset of predictors from the more complete model. Evaluation begins with the assessment of the deviance (D₀) of a model containing only an intercept (null model) which has the maximum deviance. We then estimate the deviance of a more complete model (D_p) that includes the intercept plus one or more predictors. The decrease in deviance from the null model to the model containing p predictors ( $D_{0} - D_{p})$ provides a test statistic for the significance of gain in prediction attributable to the p predictors. This measure ( $D_{0} - D_{p})$ is distributed as chi square with degrees of freedom equal to the number of predictors added to the model, here $χ$ ² (p). This is not the familiar Pearson chi square but is rather a likelihood ratio chi square, a component of maximum likelihood estimation.¹ This test statistic can also be used for testing the gain in prediction by the addition of m predictors to an equation already containing p predictors, ( $D_{p} - D_{p + m})$ which is distributed as χ²(m). Adding only a single predictor yields a test of significance of the addition of that single predictor. An alternate test statistic for the contribution of each individual predictor is also available, the Wald χ² with one degree of freedom for each predictor. The Wald χ² is the square of a z-statistic, Wald χ² = ${(\frac{b_{j}}{s_{b_{j}}})}^{2}$ with 1 degree of freedom for the predictor. The Wald χ² and the decrease in deviance will be equal in large (asymptotic) samples. The Wald χ² is the test statistic reported in both SPSS and SAS for individual predictors. The likelihood ratio χ²(1) = $(D_{p - 1} - D_{p})$ for individual predictors is recommended because it may have higher power than the Wald χ ²(1) for sample sizes realized in group research.

Fit Indices

There is no single measure of goodness of fit for nonlinear GLMs such as Poisson regression. In OLS regression we have an orthogonal partition of SS_Y, the total criterion variation, into predictable and residual variation, $S S_{\hat{Y}}$ and $S S_{Y - \hat{Y}}$ , respectively. Then the squared multiple correlation $R^{2} = \frac{S S_{\hat{Y}}}{S S_{Y}} = 1 - \frac{S S_{Y - \hat{Y}}}{S S_{Y}}$ . However, $S S_{\hat{Y}}$ and $S S_{Y - \hat{Y}}$ are not orthogonal in nonlinear GLMs. Effect size measures patterned on $R^{2}$ are employed; however, they do not have the desirable property of $R^{2}$ in OLS regression that $R^{2}$ must either remain the same or increase with addition of predictors and must always be between 0 and 1. In Poisson regression, effect size measures can even decrease as more predictors are added. Cameron and Trivedi (2013) and others refer to the fit indices employed with GLMs as pseudo-indices. The proportional reduction in deviance by the addition of p predictors to a model containing only the intercept serves as an overall model index of model fit in GLMs of a model with p predictors, ( $D_{0} - D_{p}) / D_{0}$ . This index can be calculated for the contribution of a set of m predictors to p predictors already in the regression equation, ( $D_{p} - D_{p + m}) / D_{p}$ or to the contribution of an individual predictor in an equation containing p predictors in all ( $D_{p - 1} - D_{p}) / D_{p - 1}$ . This is analogous to calculating $R^{2}$ change from linear regression.

Numerical Example: Poisson Regression Analyses

To illustrate Poisson regression we report four Poisson analyses in Table 2. All results (i.e., regression coefficients) are summarized in the linear form of Poisson regression, predicting the natural logarithm of the predicted count, for example:

\ln (\hat{c o u n t}) = b_{1} I d e n t_{c} + b_{2} C o n d + b_{3} I d e n t_{c} \times C o n d + b_{0} .

Table 2.

A series of Poisson regression models in linear form predicting $\ln (\hat{Y})$ , the natural logarithm of the predicted number of derogatory remarks based on identification with the home university, competition condition, and the Condition x Identification interaction. Coefficients with standard error in parentheses. Each row represents a different model (see text).

Analysis	Intercept	Condition^a	Identification^a	Interaction	Model deviance	Likelihood ratio chi-square (df)	Proportional reduction in deviance
Poisson regression models in linear form predicting $l n \hat{Y}$
A. Poisson, intercept only	.773 (.030)	—	—	—	913.504	—
B. Poisson, intercept plus identification	.703 (.033)	—	.186* (.016)	—	766.777	χ²(1) = 146.727^*	B vs. A^b .161
C. Poisson, intercept plus condition, identification	.617 (.035)	.835* (.066)	.188* (.016)	—	593.333	χ²(2) = 320.171^*	C vs. B .226
D. Poisson, intercept, condition, identification, plus interaction (full model)	.624 (.035)	.731* (.069)	.162* (.017)	.163* (.033)	569.757	χ²(3) = 343.747^*	D vs. C .040

Condition and identification are centered at their grand means (+.5, −.5 for high and low competition conditions, respectively; mean identification = 2.17). The interaction is formed as the cross-product of the centered predictors.

Proportional reduction in deviance = (913.504 – 766.777)/913.504 = .161

p < .001 by Wald chi-square test for individual coefficients and for likelihood ratio chi square for full model.

Ident_c is centered Identification ( $I d e n t_{c} = I d e n t - \bar{I d e n t)}$ , centered at the grand mean = 2.17. Cond is condition (coded +.5, −.5, for high vs. low competition, respectively). Use of +.5, −.5 coding yields a value of the unstandardized b₂ coefficient for condition that equals the difference between the means of the high minus low competition conditions.² Finally, Ident $\times C o n d$ is the interaction formed by multiplying the two terms (as in OLS regression with interactions; Aiken & West, 1991). In Table 2, the standard error of each regression coefficient is reported in parentheses just beneath the coefficient.

Model A of Table 2 is the null model, containing only the intercept; this model yields the null deviance $D_{0}$ = 913.504, against which Model B is compared to assess gain in prediction by the addition of the single predictor Ident_c (Model B). Testing the addition of the single predictor Ident_c in Model B yields a substantial reduction in deviance $(D_{0} - D_{p})$ = 913.504 – 766.777 = likelihood ratio χ²(1) = 146.727, p < .001. The Wald χ² is not reported in Table 2, but can be directly calculated from the information given, specifically the regression coefficient in the linear form of the regression equation $b_{j}$ and the standard error of this coefficient. $b_{j}$ for Ident_c in Model B and its standard error yield Wald χ²(1) = (.186/.016)² =135.141, a slightly smaller value than that of the likelihood ratio χ²(1) = 146.727. The proportional reduction in deviance from this single predictor $(D_{0} - D_{p}) / D_{0}$ = (913.504 – 766.777)/913.504 = .161, or a 16% reduction in deviance. Model C gives the prediction from both condition and identification, with model deviance = 593.333. Comparison of Model C to Model B in Table 2 shows the marked contribution of the Cond manipulation over the Ident_c variable. The final Model D in Table 2 shows the significant contribution of the Ident_c by Cond interaction over the two individual predictors, χ² (1) = 593.333 – 569.757 = 23.58, p < .01, a significant reduction in deviance of 23.58/593.333 = .040, or 4% gain in prediction from addition of the interaction to additive Model C.

The sequential testing we have shown before, adding each predictor and the interaction in turn was employed to illustrate the concept of reduction in deviance or gain in prediction in the likelihood ratio context. An alternative approach is to estimate the full regression equation containing all predictors. Then likelihood ratio testing is accomplished by dropping each predictor in turn and examining the contribution of that individual predictor to the equation containing all p predictors, that is, testing significance of the individual predictor in the full equation with a 1 df χ² likelihood ratio test. This is the familiar approach of Type III sums of squares in which the contribution of each predictor is considered with all other predictors held constant.

Two Forms of the Poisson Regression Equation

Often it is more convenient to interpret the results in terms of the original metric of counts rather than the transformed ln(counts). Next, we show the computation of predicted counts from the exponential form of the Poisson regression equation.

Additive Poisson Model

Consider Model C reported in linear predictor form in Table 2:

\ln (\hat{C o u n t)} = . 188 I d e n t_{c} + . 835 C o n d + . 617

We now reexpress this equation in the exponential form predicting the count itself, employing Equation (3b)

\hat{C o u n t} = e^{. 188 I d e n t_{c}} e^{. 835 C o n d} e^{. 617}

The values $e^{b_{j}}$ in the exponentiated regression equation are all given in standard statistical software along with the linear regression coefficients $b_{j}$ . The intercept $e^{. 617} = 1.853$ is the predicted arithmetic mean count for an individual with scores of 0 on both predictors. The Ident_c score = 0 at the mean Ident of the entire sample since Cond is coded +.5 and −.5 for the high versus low competition conditions, respectively. There is slight underestimation of the empirical mean count = 2.17, given in Table 1. For a case with a centered identification score of 2 (or, equivalently a raw identification score of 6), in the high competition condition (coded + .5), the predicted count is $\hat{C o u n t}$ = $e^{. 188 (2)} e^{. 835 (. 5)} e^{. 617}$ = 1.207 × 2.305 × 1.853 = 5.156, an overestimation of the mean observed count of 4.22 in the high competition condition for individuals with an Ident score of 6 (or Ident_c = 2).

Interactive Poisson Model

The addition of the interaction term yielding Model D in Table 2 increases the predicted scores by multiplying the predicted count by a value greater than 1. In the linear form of Model D with grand mean centered predictors we have:

\ln (\hat{C o u n t)} = . 162 I d e n t_{c} + . 731 C o n d + . 163 I d e n t_{c} \times C o n d + . 624

In the exponentiated form of Equation (1.4b) we have the following equation, which shows the multiplicative nature of the interaction.

\hat{C o u n t} = e^{. 162 I d e n t_{c}} e^{. 731 C o n d} e^{. 163 I d e n t_{c} \times C o n d} e^{. 624}

The intercept $e^{. 624} = 1.866$ is the predicted arithmetic mean count for an individual with scores of 0 on both predictors, at zero (0) on Ident_c or at the mean of the 500 cases on Ident, still a slight underestimate of the mean count of 2.17. For a case with a centered identification (Ident_c) score of 2 (or, equivalently Ident = 6), in the high competition condition (coded + .5), we have the predicted count as

\hat{C o u n t} = e^{. 162 (2)} e^{. 731 (. 5)} e^{. 163 (2) (. 5)} e^{. 624} = 1.383 \times 1.441 \times 1.177 \times 1.866 = 4.38

For an individual in high competition and centered identification of 3 (or, equivalently, raw identification = 7)

\hat{C o u n t} = e^{. 162 (3)} e^{. 731 (. 5)} e^{. 163 (3) (. 5)} e^{. 624} = 1.626 \times 1.441 \times 1.277 \times 1.866 = 5.58 .

Portraying Poisson Regression Results Graphically

Figure 2 provides illustrations of the predicted scores in the linear and exponential forms of the Poisson equations for the numerical example, (i.e., $\ln (\hat{c o u n t})$ and $\hat{c o u n t}$ , respectively, both without and with the Cond x Ident_c interaction). In addition, we portray the results of traditional, but nonoptimal OLS regression analyses showing predicted counts ( $\hat{c o u n t}$ ) in models without and with the interaction. Models without and with the interaction are shown in the first and second rows, respectively, of Figure 2. Columns A, B, and C of Figure 2 represent OLS regression, Poisson regression in the exponential form (Equation 1.4b), and Poisson regression in linear form, respectively (Equation 1.3). The high and low competition conditions appear separately in each graph. The Poisson regressions in exponential form (column 2) illustrate the nonlinear nature of the Poisson regression. Even without the interaction, the Poisson regression predicts a curvilinear (accelerating) increase in counts as a function of identification, with a predicted count of approximately 5 for individuals with high identification in the high competition condition. The interaction between identification and condition is positive, indicating the increasing acceleration of the predicted counts in the high competition condition and raising the predicted count for high identification individuals in the high competition condition to 5.58. In contrast, the corresponding OLS regression does not reflect this positive acceleration in counts with increasing identification. At first glance, the linear forms of the two Poisson regression models (column 3) appear to be very similar to the corresponding OLS regression models (column 1). The Y-axes are roughly of the same length; ln(6) is just under 2.0, where 6 is the highest value on the Y axis for predicted counts, and 2.0 is the highest value on the Y axis for the natural log of the predicted counts. However, recall that the dependent variables are different in the two models: ln(counts) in the Poisson models and counts in the OLS regression model. Note also that OLS regression with the interaction yields essentially the same predicted count in the high and low competition conditions for individuals with the lowest identification score of 1, whereas the Poisson regression model shows that the mean counts for these individuals differ across competition conditions, as is actually the case in the data (see Table 1).

Figure 2.

OLS regression versus Poisson regression both without interaction (row 1) and with interaction (row 2) predicting counts (columns A and B) and the natural logarithm of counts, ln(count), (column C).

Addressing Overdispersion: Alternative Regression Models for Counts

There are multiple sources of overdispersion that may produce datasets in which $σ^{2} > μ$ . First is nonindependence of cases. Participants are assumed to be independent of one another, that is, there is no clustering of cases into groups that differ in rate parameters. Group data can produce overdispersion. Second, the occurrence of an event is assumed not to modify the probability that the event will occur again. For many types of behavior, this assumption is violated. One episode of binge drinking may beget another; symptoms associated with a particular disease may be correlated; one derogatory remark hurled at another may mark the onset of a barrage of such remarks. Contagion (i.e., that the occurrence of one event increases the probability of additional events) can result in overdispersion. Third, missing predictors that account for differences in rate parameters across individuals can cause overdispersion. Our illustrative experimental example contained a manipulation of competition that produced different rates μ (arithmetic mean counts) in the two groups, as can be seen in Figure 1. Omitting the competition condition as a predictor would result in unexplained heterogeneity in the rates across the two groups.

When count data are overdispersed, regression coefficients are still unbiased. However, estimates of standard errors are too small, resulting in positively biased tests of significance. Three alternative analysis options exist that correct for overdispersion and produce improved standard error estimates: (a) quasi-Poisson regression which includes a scale parameter (> 1) that represents the excess variance as a multiplicative function of the rate parameter; (b) negative binomial regression, which generates estimates of standard errors as acombination of a Poisson probability distribution plus a second probability distribution, the gamma distribution, and (c) Poisson regression with robust standard errors. The first two approaches are explained and illustrated in Coxe et al. (2009) with SPSS and SAS syntax. SPSS and SAS syntax to accomplish these analyses is also given in the online supplementary material.

Additional Matters

Space limitations preclude a full discussion of other important topics. These include the analysis of data with excess zeros, meaning count data in which there are many more zeros than expected from a Poisson distribution, likely due to individuals in a sample who never exhibit a behavior (e.g., a substantial proportion of teetotalers in a sample estimating number of drinks consumed in an evening). Zero-inflated Poisson and hurdle models address the analysis of such data (see Cameron & Trivedi, 2013, Chapter 4; Coxe et al., 2009; Long, 1997, Chapter 8). Regression diagnostics for examination of individual cases have been developed for GLMs. Fox (2008, Chapter 15) and Coxe et al. (2009) provide accounts. On these and other topics Cameron and Trivedi (2013) and Hilbe (2011, 2014) provide more extended, technical treatments.

Generalized Linear Mixed Models: The Analysis of Count Data in the Mixed Model or Multilevel Framework

The regression models within the GLM family all assume that the residuals are independent of one another. In research on individuals this is assured by having each individual participate alone. In contrast, research on group processes often involves individuals interacting in groups, frequently with an outcome measured on each individual. The data gathered from such research typically have a clustered structure. By clustering is meant that individuals within a particular group have scores on the dependent variable, independent variable, or both that are more similar to one another than we would expect from the scores of randomly constituted groups of individuals.

Models for the treatment of clustered data in psychology date back to the publication of Bryk and Raudenbush (1992). However, the history of the analysis of clustered data extends far back into the history of the analysis of variance with contributions from many disciplines. Consequently, models for the analysis of clustered data, either individuals nested within groups, or repeated observations nested within individuals, go by many names—random coefficient regression models, mixed models, multilevel models, and hierarchical linear models. The multiple traditions and use of multiple nomenclatures arising from different disciplinary origins are fully evident in the analysis of GLMs with clustered data. The treatment of clustered data in the generalized linear model framework is a relatively recent development, still evolving at the present time. As a class, these statistical models are referred to by some authors as generalized linear mixed models (GLMM; Stroup, 2013); by others as hierarchical generalized linear models (Snidjers & Bosker, 2012). Here we briefly develop the central concepts underlying the analysis of clustered data required to understand GLMMs—the concepts of fixed versus random coefficients in regression analysis, the concept of mixed models, and the conceptual translation between mixed and multilevel models. Our rationale is that familiarity with these constructs across statistical traditions is necessary in order to navigate the literature on GLMMs, as it has evolved mainly outside psychology. The analysis of count outcomes in data sets with inherent clustering is one example of GLMMs.

A Modified Hypothetical Example With Clustered Data

We continue the fictitious in-group out-group example of counts of derogatory remarks in essays of college students as a function of their identification with their own university and their assignment to one of two experimental competition conditions. One aspect of the study is modified. After students complete measures of identification with their home university (Purdue), the students are randomly assigned to groups of five students each. From 500 participants in all, we have 100 groups. The 100 groups are then randomly assigned to one of the two competition conditions, with 50 groups per condition. The members of each group view the videos about the other university and then discuss the materials with one another. Each student then writes an individual 100-word essay evaluating the materials. The interaction among group members will typically produce clustering of the counts of individual participants’ derogatory remarks within the groups, so that participants within a single group will have counts that are more similar to one another than would be observed among randomly selected individuals. We caution that the simulated data for this second example are completely distinct from the simulated data for the first sample. No comparison should be made between the first and second example, since effect sizes are not simulated to be the same across data sets.

Fixed, Random, Mixed, and Multilevel Regression Models

Fixed effects regression

The regression models we have considered thus far in the GLM framework, including OLS regression, are fixed effects regression models. It is assumed that in the population, there is one single value of each regression parameter, that is, the regression coefficients and the regression intercept, $β_{1}, β_{2}, \dots β_{p,} β_{0}$ . In GLM models the dependent variable score $Y_{i}$ for each individual case $i$ from the population is assumed to be produced by these population parameters. Returning to the simulated example, consider the linear predictor form of the population Poisson regression model with the one predictor identification:

\ln (\hat{c o u n t_{i}}) = β_{1} I d e n t i f i c a t i o n_{i} + β_{0}

(2.1)

Individual members of the population vary in level of identification, an individual difference variable $X_{i}$ , but one fixed population regression coefficient $β_{1}$ captures the relationship of identification to the production of derogatory remarks for the whole population. Moreover, the single regression intercept $β_{0}$ applies to all members of the population. In sum, the fixed effects regression model has one set of population regression parameters that applies to all individuals within the population, and individual predicted scores are produced by that one set of parameters.

Random coefficient regression

Random coefficient regression is the core regression model applied to the analysis of clustered data. In our example the clusters are created in the experimental setting through assignment of individuals to small groups in which the individuals interact. In other research, the clustering might stem from the use of preexisting groups such as campus or community groups. Both laboratory groups and preexisting groups might differ in their mean number of derogatory remarks (the rate parameter in Poisson regression), for example, if the group contained a vocal member who was particularly derogatory.³ Consider again the relationship of the single independent variable identification to the production of derogatory remarks. This relationship might also differ in the groups. Nonetheless, we still consider the whole population of students pooled over all clusters as being represented by fixed parameters $β_{0}$ and $β_{1}$ for the population intercept and slope, respectively, in the prediction of number of derogatory remarks from identification. At the same time, each group $j$ can be characterized more specifically as having its own intercept $b_{0 j}$ and regression coefficient $b_{1 j}$ . These group intercepts and regression coefficients are assumed to vary randomly around the fixed population parameters $β_{0}$ and $β_{1}$ (random slopes and intercepts are illustrated in Cohen et al., 2003, p. 540). We assume that these distributions are normal in the population so the distribution of intercepts is characterized⁴ by a mean of $β_{0}$ with variance $τ_{00}$ and the distribution of slopes is characterized by a mean of $β_{1}$ and variance $τ_{11}$ . If we consider the observed groups as a random sample of all possible groups in the population, then the variances of the observed $b_{0 j}$ and $b_{1 j}$ coefficients across (five-person laboratory) groups in the experiment yield an estimate of the variance of these coefficients in the population of all groups. The variances of the observed coefficients are referred to as the variance components of the model. A regression model that captures random variation in regression coefficients across groups is termed a random coefficient regression model.

Mixed models

The term mixed models is used to describe regression models that contain both fixed and random effects. The overall estimates of the population regression coefficients $β_{0}, β_{1}, \dots β_{p}$ are estimates of fixed effects. Note that the regression coefficients $β_{1}$ … $β_{p}$ plus the regression constant $β_{0}$ contain no subscript for group; they are fixed coefficients that pertain to the population of all groups. The estimates of the variances of the intercepts and slopes across groups, (i.e., the variance components) are estimates of the random effects. Mixed models are also called random coefficient models.

A model may be specified in which the intercepts vary across groups (e.g., mean number of derogatory remarks differs across groups), whereas the relationship of a predictor to the outcome is constant across groups (e.g., the same relationship of identification to number of derogatory comments). Such a model is termed a random intercept model. Alternatively, both the intercepts and slopes may be specified as random, yielding a random slope and intercept model. In practice, we examine the statistical significance of these variance components. Often a model specified initially to have both random slopes and intercepts may be simplified to random intercepts alone if there is no evidence for random slopes. Finally, if it is found that neither the slopes nor the intercepts exhibit significant random variation, the model reduces from a mixed model to a fixed model, that is, GLMM simplifies to GLM.

The variance components in a mixed model may represent important theoretical aspects of the model. If group membership is driven by important forces (e.g., communities, see e.g., Rios, Aiken, & Zautra, 2012), then the variance components may tell us whether theoretically important relationships vary as a function of group membership (random slopes), and whether the groups differ in levels on important outcomes (random intercepts).

Multilevel models

In the multilevel framework, familiar to psychologists, we conceptualize clustered data as existing at multiple hierarchical levels with measurements specifically tied to each of the multiple levels. In group data, the lowest level of aggregation is the individual case (Level 1); the data for the individual cases are aggregated into separate groups (Level 2). In our example, identification is a Level 1 predictor, measured separately on each individual. Competition (high, low) is a Level 2 predictor that applies equally to all groups (conceptually all members of a group have the same score on competition, which is the code for level of condition to which the group was assigned). Such a model will contain a regression coefficient for each predictor, $β_{1}$ and $β_{2}$ , respectively. The model may also contain a cross-level interaction between competition at Level 2 and identification at Level 1. Conceptually, this interaction indicates that the relationship of identification to the count outcomes of individuals at Level 1 is moderated by the competition condition at Level 2. The regression coefficients for identification, competition condition, and their interaction are the fixed effects of the multilevel model. The variance components represent the random effects in the multilevel model; the same terminology is used in the multilevel modeling literature as in the mixed model literature. In the multilevel model we add “Level 1, Level 2” thinking to the model structure—some variables (here, identification) are measured on individuals; other variables (here, competition) are measured (manipulated) on groups. The variance components (variances of the regression coefficients and regression intercept across groups) are associated with the Level 1 predictors.⁵ The mathematics of the mixed and multilevel models are identical; only the conceptualization changes.

A final note: The intraclass correlation

In multilevel models with continuous outcome variables, researchers often calculate the intraclass correlation coefficient (ICC) as a useful index of the degree of clustering. The ICC is based on the extent of variance of the intercepts (i.e., between class variance in the ANOVA, often noted as $τ_{00}$ in the multilevel literature) relative to the total variance of the individual outcomes. The total variance is the sum of the between class variance and the pooled within class variance in ANOVA, the latter often noted as $σ_{ϵ}^{2} .$ Specifically, the ICC is the proportion of total variance among the individual cases that is attributable to the clustering, where the ICC = $\frac{τ_{00}}{τ_{00} + σ_{ε}^{2}}$ . We will see later that this measure is not defined for clustering of count outcomes. An alternative assessment and index of clustering is provided next.

Analysis of Generalized Linear Mixed Models (GLMMs)

All the multilevel (or mixed) model structures introduced here—fixed population regression coefficients, variance components that reflect the variance of the group coefficients around corresponding fixed population coefficients, random intercept models, random slope and intercept models, cross-level interactions—apply to GLMMs. However, a number of differences and complexities are introduced into the analysis of GLMMs that do not exist for GLMs. These complexities are primarily in the realms of parameter estimation and statistical inference. Of importance, the Poisson distribution is a one-parameter distribution in which the mean and variance are equal and represented by a single parameter μ. Once the mean or rate parameter μ is specified, there is no independent pooled within class variance (at Level 1) in a count GLMM, that is, $σ_{ϵ}^{2}$ does not exist. Therefore ICC is not defined for count GLMMs.

Estimation

As explained previously, GLMs employ maximum likelihood estimation (MLE) of model parameters. MLE requires that a likelihood function can be stated mathematically that gives the likelihood of the observed data given specific values of the parameters of the model. Once the function is specified, it must be solved for the parameter estimates, that is, the estimates must be found that maximize the likelihood function. In some cases, there is an analytic solution for the estimates, that is, a set of equations can be derived that give the solution to the estimates, as in the equations for the regression coefficients and intercept in OLS regression. In general for GLMs, there is no analytic solution, and the solution to the coefficient estimates is found iteratively (i.e., through systematic trial and search algorithms implemented by computer). GLMs are fixed models; there are no random components. GLMMs introduce a major new layer of complexity in finding parameter estimates. By definition, GLMMs all contain random effects due to the clustering of the data. Theoretically all possible values of the random effects must be incorporated into the estimation procedure. The estimation process becomes slow, even with a fast computer, or even computationally infeasible (Bolker et al., 2009). It is necessary to use an approximation to obtain a solution. A number of approaches to estimation have been developed for GLMMs, which can be organized into two broad classes—linearization methods and integral approximation (Stroup, 2013, p. 140). Bolker et al. (2009) provide a very accessible, nonmathematical overview of these methods. Stroup (2013) provides extensive mathematical detail. Here we provide a brief introduction.

Linearization methods

We have seen that the GLMs are nonlinear models; GLMMs are as well. The linearization approach to estimation for GLMMs uses a linear transformation of the nonlinear outcome variable for parameter estimation, essentially transforming the nonlinear estimation problem into a simpler linear estimation problem (Stroup, 2013, Chapter 4). This linearization approach is referred to as penalized quasilikelihood (PQL) or pseudolikelihood (PL). PL is widely implemented in current software (it is currently the only estimation method available in SPSS for GLMMs and the default in SAS for GLMMs). The method yields what is referred to as a quasilikelihood, rather than a true likelihood, which raises issues of its utility for inference. Stroup (2013) cautions that PL may not yield usable estimates of the likelihood; Bolker et al. (2009) indicates that PL often gives biased estimates when the random effects (variance components) are large, and more specifically, that PL does not work well in Poisson models with very low arithmetic mean counts. Further, PL may not yield parameter estimates in count GLMMs when the data contain many zeroes (see also Silva & Tenreyro, 2010).

Integral approximation

In general, maximum likelihood estimation involves solving for the area under the likelihood function. Integral approximation employs numerical approaches to estimate the area under the likelihood function. The integral approximation methods estimate the actual likelihood by approximating the area under complex curves with a series of rectangles. That the integral approximation approaches estimate the true likelihood rather than a quasilikelihood gives these methods an advantage in inference. Two methods of integral approximation are Laplace approximation and adaptive Gauss–Hermite quadrature; the term adaptive refers to a method that controls approximation error by improving the steps to estimation in the course of the estimation process. Both are more accurate than PL but slower than PL. Gauss–Hermite is more accurate than Laplace approximation, but it is even slower and is thus restricted to use with models than have only a few random effects.

Inference for fixed effects in GLMMs

With mixed models, including GLMMs, we have two distinct classes of parameters to be tested—the fixed effects (regression coefficients) and the random effects (the variance components). We present a brief overview of inference here. Our overview explains what researchers will encounter in statistical testing with GLMMs and what results researchers will see and not see in GLMM computer output—for example, whether likelihood ratio (LR) tests even exist, depending on the method of estimation employed, and why we see F and t tests rather than the test statistics we encountered for GLMs.

We focus first on fixed effects. There are two approaches to hypothesis testing and confidence interval estimation in linear models, which we have encountered before. First is the likelihood ratio (LR) approach which computes model deviances and differences in deviances, presented earlier in the GLM section of this paper. This approach yields LR χ2 tests of significance of prediction from full models, from subsets of predictors, and from individual predictors. Second is the Wald based approach, which yields the tests of significance of prediction for the full model as well as z-tests or, equivalently, one degree of freedom Wald χ² tests for individual predictors. Wald tests of individual predictors are standardly reported in software for the analysis of GLMs (and in familiar structural equation modeling software).

For GLMMs estimated by linearization methods such as penalized likelihood, LR tests are not available—PL does not estimate the true likelihood. Therefore, only Wald tests are reported for PL. In contrast, for GLMMs estimated by integral approximation, both the LR and the Wald approach to testing are available; testing can proceed with the same likelihood ratio χ² tests we encountered for GLMs. This result occurs because integral approximation estimates true likelihoods. The deviance measures on which statistical inference is based are functions of these estimates of the true likelihoods. As pointed out earlier in the discussion of GLMs, the difference between the deviance of a model with a set of parameters and the deviance of a reduced model with a subset of the same parameters is a LR χ² test. The LR χ² tests we encountered for GLMs can be used in the same way in GLMMs—for overall prediction, for prediction by a subset of predictors, for prediction by an individual predictor. There is one caveat to the use of LR tests in GLMMs estimated with integral approximation: For highly complex models, integral approximation may be computationally very intense and slow, in which case Wald tests, rather than LR tests, are often used (Stroup, 2013, pp. 159–160).

The transition from GLMs to GLMMs also yields a change in the statistics that are reported. For GLMs standard software reports one degree of freedom Wald χ² s for tests of significance of individual predictors. In contrast, for GLMMs standard software reports t tests of the significance of individual fixed effects regression coefficients. In addition, in GLMMs F tests are reported for Type III sums of squares estimates of fixed effects. There are relationships between the tests. For example, these F tests for Type III sums of squares are the squares of the t tests reported for the individual coefficients, since both test the significance of each predictor all other predictors held constant. Cameron and Trivedi (2013, p. 51) suggest that the change in labeling tests from Wald z to t in some software is for consistency with other familiar software (e.g., for linear regression). In fact, regardless of labeling as z versus t, Wald tests are z tests only asymptotically. The Type III F tests for individual predictor significance are asymptotically equal to one degree of freedom likelihood ratio chi square tests. In the numerical example with clustered data, the likelihood ratio test for reduction in deviance due to the addition of condition (from Model C to Model D in Table 3) is χ²(1) = 6.19; the Type III F test for condition is F(1, 98) = 6.40. These values would be theoretically expected to be equal at asymptote.

Table 3.

Series of multilevel (mixed model analyses) to examine fixed and random effects and final model estimates. Each row represents a different model (see text).

3a. Series of multilevel models.

Analysis	Fixed effects	Variance components	Model deviance	Likelihood ratio chi square (df)	Proportional reduction in deviance
A. Unclustered, intercept only	Intercept	—	4775.77	—	—
B. Clustered, intercept only	Intercept	Intercept	2763.64	χ²(1)= 2012.13**	B vs. A = .42^a
C. Clustered, intercept, identification	Intercept, identification	Intercept	2149.29	χ²(1) = 614.35**	C vs. B = .29
D. Clustered, intercept, identification, condition	Intercept, identification, condition	Intercept	2143.10	χ²(1) = 6.19*	D vs. C = .002
E. Clustered, intercept, identification, condition, identification by condition	Intercept, identification, condition, Condition x Identification	Intercept	2142.77	χ²(1) = .33 ns	E vs. D = .00015

Note. Condition and identification are centered at their grand means (+.5, −.5 for high and low competition conditions, respectively; mean identification = 2.17). The interaction is formed as the cross-product of the centered predictors.

Proportional reduction in deviance = (4775.77 – 2763.64)/4775.77 = .42.

p < .001; *p < .05.

3b. Final model estimates of fixed and random effects: Variance component: Intercept variance $\hat{τ_{o o}}$ = .4998.

Fixed effects	Estimate	Standard error	df	t value	p
Intercept	1.1029	.0768	98	14.35	< .0001
Identification	0.2948	.0132	398	22.32	< .0001
Condition	0.3701	.1527	98	2.42	.0172
Ident x Condition	0.0151	.0265	398	0.57	.5689

There is complexity in the computation of the appropriate degrees of freedom for statistical tests of fixed effects in GLMMs. Thinking back to simple fixed effects ANOVA and OLS regression, there is only one error term for the whole analysis, $M S_{w i t h i n c e l l}$ or $M S_{r e s i d u a l}$ , respectively, that is, the estimate of the one variance component $σ_{ϵ}^{2}$ . In models with random effects, the error terms for some effects are combinations of variance components, referred to as pooled error terms. Although training on classic mixed model ANOVA in psychology was common in earlier times and represented in the textbooks of that earlier era, most training now focuses on fixed models, and thus this combining of error terms may be unfamiliar. Pooled error terms require specialized degrees of freedom reflecting the pooling of heterogeneous sources of variance. Readers may be familiar with Satterthwaite approximate degrees of freedom for pooled error terms in ANOVA. The Satterthwaite approximate degrees of freedom can be used with GLMMs, but only with pseudolikelihood (PL) and not integral approximation estimation methods. The Kenward–Roger approximate df is preferred to the Satterthwaite df because it yields more accurate Type 1 error rates in mixed models. However, like Satterthwaite, Kenward–Roger approximate df can only be used with PL estimation. With integral approximation methods, one approach is to separate the degrees of freedom into between class and within class degrees of freedom, as in ANOVA (Schluchter & Elashoff, 1990).

Inference for variance components in GLMMs

As we pointed out, variance components may play an important role in theory, for example, in uncovering variability across groups in important aspects of relationships among constructs. The testing of significance of variance components in GLMMs has its own complexity. Stroup (2013, Chapter 5) states forceful conclusions about strategies for testing variance components; we reiterate them here. First of all, inference for variance components does not exist with PL estimation. Second, Wald statistics for testing the significance of variance components are grossly conservative and should not be used. Third, likelihood ratio testing of significance of variance components can be carried out only when integral approximation is used to estimate the model likelihood.

Numerical Example, Clustered Data

We analyzed the clustered data set, examining both fixed and random effects to arrive at a final model characterizing experimental outcomes. (Again we caution that the data sets from the two examples are distinct, and comparisons should not be made across data sets with regard to effect sizes of fixed effects.)

We employed integral approximation with adaptive quadrature, so that estimates of the likelihood function and therefore model deviances would be available. We addressed the significance of both fixed and random effects in a sequence of analyses summarized in Table 3.

While there is no estimate of the ICC in Poisson regression, one can assess the magnitude of clustering by estimating two models. The first is Model A of Table 3, an “empty” fixed model containing only the fixed intercept, with no clustering of the data. The second is Model B, a random coefficient model in which the overall intercept is estimated, clustering is specified, intercepts are considered random, and the variance component $τ_{00}$ of the intercept is included. In Table 3, there was 42% reduction in deviance from model A to model B, indicating substantial clustering on the dependent variable⁶ (i.e., a large variance component of the group intercepts, $τ_{00}$ ). In model C, identification was added as a fixed predictor to model B; model C thus contains a fixed slope and intercept along with the random intercept variance component $τ_{00}$ . The addition of identification yielded a further 29% decrease in deviance. In a subsequent model, not shown in Table 3, the variance $τ_{11}$ of the slopes for the regression of derogatory comments on identification was added as a second variance component. While this model converged and produced parameter estimates, the variance of the slopes was estimated to be zero; there was no change in the model deviance from the model without the random slope variance (i.e., the deviance for the model with both random slope and intercept was 2149.29, the same value as for model C in Table 3 with only the random intercept). Model C was retained as base model for further comparison. In Model D, the condition main effect was added to the model containing identification, the fixed intercept, and the random intercept variance component $τ_{00}$ , yielding a small but significant reduction in deviance. Finally, in Model E, the addition of the Identification x Condition interaction to the model yielded no reduction in deviance (i.e., a nonsignificant interaction).

As mentioned before, an alternate strategy to this step-by-step introduction of fixed and random effects is to begin with a complete model containing all fixed and random effects and to step down, eliminating nonsignificant effects. One would examine the variance components first with all the fixed effects in the model, asking whether there is evidence of random slopes. If not, then the random slope parameter for random slopes would be eliminated from the model. Likelihood ratio chi square tests would be employed. The challenge with this strategy is the models with all variance components included might well not be able to be estimated, as would have been the case in this example, in which the slope variance was clearly zero.

Conditional versus marginal coefficient estimates

The parameter estimates of the regression coefficients (fixed effects) add one final layer of complexity. The parameter estimates presented here are conditional estimates as opposed to marginal estimates. This subtle distinction becomes prominent in clustered models with nonnormally distributed outcomes (i.e., GLMMs) due to the nonlinear transformation of predicted values that occurs in GLMMs. (This distinction is related to the distinction between conditional and marginal probability in contingency tables, hence the similar terminology; see Agresti, 2007, Chapter 2.)

Conditional estimation, as in the models for clustered data that we have presented here, provides estimates of the parameter coefficients, conditional upon the random effects (i.e., the clusters). These models allow us to estimate the intercept and regression coefficient for each cluster in the data. The intercept and the regression coefficient estimate for the Level 1 predictor that are reported can be thought of as modeling the effect for the average cluster. Conditional estimates are reported in software (e.g., SAS GLIMMIX) when we include variance components in our model, such as the intercept variance included in models B through E of Table 3.

Marginal estimation provides an estimate of the slope for each Level 1 predictor and the intercept in the total population, ignoring clustering (i.e., ignoring the random effects). The most prominent approach to estimating marginal effects is generalized estimating equations (GEEs).

While conditional and marginal models give numerically different results, have different interpretations, and may lead to different statistical conclusions, it is important to note that neither model is the “correct” model in all situations. Conditional and marginal models estimate different effects and so produce different values; the choice of model depends on the study design and research questions. For example, Hubbard et al. (2010) recommend using conditional models when the cluster is theoretically relevant and the variability of the effect across clusters (reflected in the variance of the random intercept and random slope) is of interest. This is often the case in psychological research, where the cluster (i.e., group, neighborhood, or family) is of particular interest to the researcher. Marginal models are recommended when the larger population from which clusters are drawn is of particular interest and the cluster and variability of the effect across clusters is considered a nuisance. Additionally, conditional models provide information about the cluster-level variability, whereas marginal models do not. We refer interested readers to more detailed technical articles comparing conditional versus marginal models, such as Hu, Goldberg, Hedeker, Flay, and Pentz (1998) and Hubbard et al. (2010). We note that marginal models are often treated as synonymous with GEE in current literature.

In the numerical example presented in this article, individuals are clustered into groups, resulting in nonindependence. In this hypothetical study, subjects are randomly assigned to each group. This results in groups that are largely exchangeable; the differences between groups (on both outcome and predictor) are attributable to random error and are probably not of particular research interest. It may not be surprising that conditional models presented in Table 3 and marginal models (not presented here) for this example lead to similar numerical results and identical statistical conclusions.

Contrast this situation with a community study, in which individuals have been sampled or selected from different neighborhoods. In this situation, neighborhoods are not exchangeable. Neighborhoods are different from one another, and the effect of a predictor on an outcome may vary across neighborhoods. Researchers are interested in both the differences among neighborhoods and what predicts these differences. In such a study conditional models would be preferred and would likely produce different results from marginal models. In research comparing existing groups, the characteristics of individual groups are often important (for example, in a study in which different religious groups are selected). Here conditional models that capture the unique characteristics of groups as part of the modeling appear to be the appropriate choice.

We have summarized the final fixed effects part of the regression model with conditional coefficients in Table 3b. Readers will note the partitioning of degrees of freedom into what are termed between degrees of freedom—those for the Level 2 fixed effects (here condition and the overall intercept, with 98 df per effect) and within degrees of freedom—those for the Level 1 effects (here, identification and the Identification x Condition interaction, with 398 df per effect).

The t tests of the fixed effects mirror the results of the likelihood ratio tests of significance of the fixed effects. In Table 3, we show the likelihood ratio chi square tests χ²(1) = 614.35, 6.19, and .33 for the sequential addition of identification (Model C), condition (Model D), and the Identification x Condition interaction (Model E), respectively. Not shown in Table 3 is the sequence of Wald t tests that correspond to this sequence of χ²(1) tests. The Wald t(399) = 13.83 (p < .01) for the addition of identification (Model C), Wald t(98) = 2.53, (p < .05), for the addition of condition (Model D), and Wald t(398) = 0.57, ns, for the Identification x Condition interaction. The results reported in Table 3b provide the Wald t values for the fixed effects in final Model E.

Summary and Conclusion

Research with groups (clustered individuals) poses challenges for analysis, whether the groups under investigation are intact preexisting groups or are groups created in the experimental setting. Research that involves count dependent variables poses a second set of challenges for analysis. The implication of either grouped (clustered) data or count outcomes for analysis is that standard fixed effects analysis of variance and OLS regression will not produce proper results; different statistical models are required. The challenges produced by clustering of data into groups are addressed with multilevel (or mixed) statistical models. The challenges produced by count data are addressed in the generalized linear model (GLM), a broad class of regression models designed to handle a variety of (nonclustered) dependent variables that are not continuous and normally distributed. In group research with count dependent variables these two sets of challenges intersect. The challenges are simultaneously addressed in a new class of statistical models, the generalized linear mixed model (GLMM), which integrates GLM models with multilevel (mixed) models. In order to employ GLMMs effectively for data analysis the researcher must have familiarity with both the structure of GLMs and the structure of multilevel (mixed models).

We began this article with an explanation of the problems of employing OLS regression for count outcomes. We then provided a broad overview of GLMs, introducing the structure of all GLMs. We considered in detail the analysis of count outcomes for independent (nonclustered, single-level data) with the most common analysis for count data, Poisson regression, a member of the GLM family. We showed the two forms in which Poisson regression equations are expressed: (a) a linear equation in which the predicted score is an additive function of the predictors and in which the form of the predicted scores is the natural logarithm of the predicted count, and (b) an exponential equation in which the predicted count is a multiplicative function of the predictors and the form of the predicted score is a predicted count. We provided an illustration of an analysis of a count outcome in a simulated group experiment, including showing the impact of interaction between predictors in comparison with predictor interactions in OLS regression. We then moved to the analysis of clustered data. We provided a basic understanding of the effects of clustering, and the features of the random coefficient (multilevel) models that address clustering by capturing the variation that may exist across groups. After providing the foundation in the two necessary areas, we introduced the new generation of generalized linear mixed models that simultaneously address clustering and count outcomes.

These GLMM models represent the current state of the art, and as such they are in an active state of development. The GLMM introduces new complexities of which users must be aware, among them the existence of multiple approaches to parameter estimation (linearization, integral approximation), multiple approaches to inference (likelihood ratios, Wald statistics), multiple estimates of degrees of freedom for statistical inference, multiple estimates of regression coefficients (conditional, marginal) by which we judge the outcomes of our research. We provided overviews of these multiple approaches, and explained how they are intertwined—the choice of one approach for estimation determines the approaches available for inference and the estimates of regression coefficients. Finally, we presented an example of the analysis of multilevel count data using generalized linear mixed models. We illustrated choices we made in analysis, among them, the use of integral approximation methods for estimation and the reporting of conditional estimates of parameters. We illustrated the interpretation of the coefficients, tests of significance, and measures of fit. In the online supplemental material, we have provided SPSS and SAS computer code for single level Poisson regression, and SAS computer code for GLMM. Data sets for the examples also are provided in the supplemental materials. Generalized linear mixed models offer a set of methods for analyzing multiple count data that provide proper answers to key research questions. These models also hold the promise of stimulating researchers to pose new types of questions and providing the statistical tools to answer those questions.

Footnotes

Funding

This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors.

Notes

References

Agresti

A. A.

(2007). An introduction to categorical data analysis (2nd ed.). Hoboken, NJ: Wiley.

Aiken

L. S.

West

S. G.

(1991). Multiple regression: Testing and interpreting interactions. Newbury Park, CA: Sage.

Atkins

D. C.

Baldwin

S. A.

Zheng

Gallop

R. J.

Neighbors

(2013). A tutorial on count regression and zero-altered count models for longitudinal substance use data. Psychology of Addictive Behaviors, 27, 166–177. doi:10.1037/a0029508

Atkins

D. C.

Gallop

R. J.

(2007). Rethinking how family researchers model infrequent outcomes: A tutorial on count regression and zero-inflated models. Journal of Family Psychology, 21, 726–735. doi:10.1037/0893–3200.21.4.726

Berk

McDonald

(2008). Overdispersion and Poisson regression. Journal of Quantitative Criminology, 24, 269–284. doi:10.1007/s10940–008–9048–4

Bolker

B. M.

Brooks

M. E.

Clark

C. J.

Geange

S. W.

Poulsen

J. R.

Stevens

M. H. H.

White

J. S.

(2009). Generalized linear mixed models: A practical guide for ecology and evolution. Trends in Ecology and Evolution, 24, 127–135. doi:10.1016

Brewer

M. B.

(1979). In-group bias in the minimal intergroup situation: A cognitive-motivational analysis. Psychological Bulletin, 86, 307–324. doi:10.1037/0033–2909.86.2.307

Bryk

A. S.

Raudenbush

S. W.

(1992). Hierarchical linear models in the social and behavioral sciences. Newbury Park, CA: Sage.

Cameron

A. C.

Trivedi

P. K.

(2013). Regression analysis of count data (2nd ed.). New York, NY: Cambridge University Press.

10.

Cohen

West

S. G.

Aiken

L. S.

(2003). Applied multiple regression/correlation analysis for the behavioral sciences (3rd ed.). Mahwah, NJ: Erlbaum.

11.

Coxe

West

S. G.

Aiken

L. S.

(2009). The analysis of count data: A gentle introduction to Poisson regression and its alternatives. Journal of Personality Assessment, 91, 121–136. doi:10.1080/00223890802634175

12.

Coxe

West

S. G.

Aiken

L. S.

(2013). Generalized linear models. In Little

(Ed), Oxford handbook of quantitative methods (pp. 26–51). New York, NY: Oxford University Press.

13.

Enders

C. K.

(2005). Maximum likelihood estimation. In Everitt

B. S.

Howell

D. C.

(Eds.), Encyclopedia of statistics in the behavioral sciences (pp. 1164–1170). Chichester, UK: Wiley.

14.

Fox

(2008). Applied regression analysis and generalized linear models (2nd ed.). Thousand Oaks, CA: Sage.

15.

Guilford

J. P.

(1956). Fundamental statistics in psychology and education. New York, NY: McGraw Hill.

16.

Hedeker

Gibbons

R. D.

(2006). Longitudinal data analysis. Hoboken, NJ: Wiley.

17.

Hilbe

J. M.

(2011). Negative binomial regression (2nd ed.). New York, NY: Cambridge University Press.

18.

Hilbe

J. M.

(2014). Modeling count data. New York, NY: Cambridge University Press.

19.

Hox

J. J.

(2010). Multilevel analysis: Techniques and applications (2nd ed.). New York, NY: Routledge.

20.

F. B.

Goldberg

Hedeker

Flay

B. R.

Pentz

M. A.

(1998). Comparison of population-averaged and subject-specific approaches for analyzing repeated binary outcomes. American Journal of Epidemiology, 147, 694–703.

21.

Hubbard

A. E.

Ahern

Fleischer

N. L.

van der Laan

Lippman

S. A.

Jewell

… Satariano

W. A.

(2010). To GEE or not to GEE: Comparing population average and mixed models for estimating the associations between neighborhood risk factors and health. Epidemiology, 21, 467–474. doi:10.1097/EDE.0b013e3181caeb90

22.

King

(1988). Statistical models for political science event counts: Bias in conventional procedures and evidence for the exponential Poisson regression model. American Journal of Political Science, 32, 838–863. Retrieved from http://www.jstor.org/stable/2111248

23.

Kmenta

(1986). Elements of econometrics (2nd ed.). New York, NY: McMillan.

24.

Long

J. S.

(1997). Regression models for categorical and limited dependent variables. Thousand Oaks, CA: Sage.

25.

O’Hara

R. B.

Kotze

D. J.

(2010). Do not log–transform count data. Methods in Ecology and Evolution, 1, 118–122. doi:10.1111/j.2041-210X.2010.00021.x

26.

Raudenbush

S. W.

Bryk

A. S.

(2002). Hierarchical linear models: Applications and data analysis methods (2nd ed.). Thousand Oaks, CA: Sage.

27.

Rios

Aiken

L. S.

Zautra

(2012). Neighborhood contexts and the mediating role of neighborhood social cohesion on health and psychological distress among Hispanic and non-Hispanic residents. Annals of Behavioral Medicine, 43, 50–61. doi:10.1007/s12160-011-9306-9

28.

Rossotti

N. G.

Winter

D. A.

Watts

M. H.

(2006). Trust and dependency in younger an older people. In Caputi

Foster

Viney

(Eds.), Personal construct psychology: New ideas (pp. 163–175). Hoboken, NJ: Wiley.

29.

Schluchter

M. D.

Elashoff

J. T.

(1990). Small sample adjustments to tests with unbalanced repeated measures assuming several covariance structures. Journal of Statistical Computation and Simulation, 37, 69–87. doi:10.1080/00949659008811295

30.

Silva

J. J. C. S.

Tenreyro

(2010). On the existence of the maximum likelihood estimates in Poisson regression. Economics Letters, 107, 310–312. doi:10.1016/j.econlet.2010.02.020

31.

Snijders

T. A. B.

Bosker

(2012). Multilevel analysis: An introduction to basic and advanced multilevel modeling (2nd ed.). Thousand Oaks, CA: Sage.

32.

Stoel

R. D.

Garre

F. G.

Dolan

Wittenboer

(2006). On the likelihood ratio test in structural equation modeling when parameters are subject to boundary constraints. Psychological Methods, 11, 439–455. doi:10.1037/1082-989X.11.4.439

33.

Stroup

W. W.

(2013). Generalized linear mixed models: Modern concepts, methods, and applications. Boca Raton, FL: CRC Press.

34.

Tov

K. L.

Lin

Qiu

(2013). Detecting well-being via computerized content of brief diary entries. Psychological Assessment, 25, 1069–1078. doi:10.1037/a0033007

35.

Tukey

(1977). Exploratory data analysis. Reading, MA: Addison-Wesley.

36.

Verbeke

Molenberghs

(2000). Linear mixed models for longitudinal research. New York, NY: Springer.

37.

West

B. T.

Welch

K. B.

Galecki

A. T.

(2007). Linear mixed models: A practical guide using statistical software. Boca Raton, FL: Chapman & Hall/CRC.

38.

Xie

Tao

McHugo

G. J.

Drake

R. E.

(2013). Comparing statistical methods for analyzing skewed longitudinal count data with many zeros: An example of smoking cessation. Journal of Substance Abuse Treatment, 45, 99–108. doi:10.1016/j.jsat.2013.01.005

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.10 MB