IRT–ZIP Modeling for Multivariate Zero-Inflated Count Data

Abstract

This study introduces an item response theory–zero-inflated Poisson (IRT–ZIP) model to investigate psychometric properties of multiple items and predict individuals' latent trait scores for multivariate zero-inflated count data. In the model, two link functions are used to capture two processes of the zero-inflated count data. Item parameters are included to investigate item performance from both propensity and level perspectives. The application of the model was illustrated by analyzing the substance use data from the National Longitudinal Study of Youth. A simulation study based on the empirical data analysis scenario showed that the item parameters can be recovered accurately and precisely with adequate sample sizes. Limitations and future directions are discussed.

Keywords

multivariate zero-inflated count data IRT–ZIP item properties

In psychological and educational research, counts of a behavior or an event in a time interval of specified length often are collected through observations, surveys, or experiments. Cameron and Trivedi (1998) defined an event count as the number of times an event occurs, which is a realization of a nonnegative integer-valued random variable. For example, in a classroom observational study, frequencies of specific behaviors such as effective interactions between students and their teacher are observed and recorded over a period of 10 minutes. Or, in a substance use study, the number of cigarettes smoked by teenagers in a day during the last 30 days can be collected via a survey questionnaire. The Poisson regression model (e.g., Cox & Lewis, 1966; Frome, Kutner, & Beauchamp, 1973; Gart, 1964; McCullagh & Nelder, 1989) is a typical method for analyzing cross-sectional count data. The mixed-effects Poisson factor model (Böckenholt, Kamakura, & Wedel, 2003) can be applied to analyze multivariate count data.

Count data that contain excessive numbers of zeros are called zero-inflated count data. Zero-inflated count data usually come from two distinct processes. One process is related to whether the event has a chance to happen, and the other process is related to the count of the event, given that the event does have a chance to happen. For example, in investigating the number of cigarettes smoked by a participant in a day during the last 30 days, some participants report zeros because they are nonsmokers and some participants report zeros because they happen not to smoke during that period although they are smokers. We define the former zeros (qualitative zeros) as perfect zero states and the latter zeros (quantitative zeros) as zero counts following a Poisson process. Zero-inflated count data can be analyzed using the zero-inflated Poisson model (ZIP model; Lambert, 1992). The ZIP model with random effects (ZIPR model; e.g., Hall, 2000; Lee, Wang, Scott, Yau, & McLachlan, 2006; Liu & Powers, 2007; Min & Agresti, 2005; Rabe-Hesketh & Skrondal, 2007) can be used to analyze multilevel or longitudinal zero-inflated count data.

What to do with measurement errors is always an issue for social and behavioral research. A typical method to reduce the effect of measurement errors is to measure a latent construct by collecting multivariate data using multiple manifest variables. For continuous multivariate data, factor analyses are the typical methods for examining the relationship between manifest variables and latent factors (e.g., Cattell, 1978; Fuller, 1987; McDonald, 1985). For categorical multivariate data, item response theory models (IRT models) are the typical methods to analyze the psychometric properties of items and to obtain latent trait scores for the measured constructs by modeling the relationship between items and trait scores (e.g., Embretson & Reise, 2000; Lord, 1980; Lord & Novick, 1968; Rasch, 1960).

However, how to model multivariate zero-inflated count data to analyze psychometric properties of zero-inflated count items and how to obtain person trait scores from multiple zero-inflated count items were not discussed in previous literature. When analyzing multivariate zero-inflated count data, although one can form and model composite scores, a disadvantage of this approach is that the distribution of the composite scores is extremely skewed right due to excess zeros and the weights of items forming the composite scores are usually arbitrarily set to be 1. Therefore, the purposes of this study are to (a) introduce a model for multivariate zero-inflated count data to analyze psychometric properties of zero-inflated count items and to predict person trait scores from multiple items; (b) explore an empirical data set to illustrate the application of the model; and (c) examine the performance of the proposed model via a simulation study with different sample sizes.

Review of Precursors of the Models

ZIP Model

For handling univariate zero-inflated count data, Lambert (1992) introduced the ZIP model. The ZIP model uses two link functions, a logit link function and a log link function, to capture statistical features of the two processes: the perfect zero state and the Poisson process, for zero-inflated count data. The model is expressed as follows:

Y_{i} ~ {\begin{cases} 0, with probability 1 - p_{i} \\ P o i s s o n (λ_{i}), with probability p_{i} \end{cases} \log (λ_{i}) = x_{i}^{'} β

log (λ_{i}) = x_{i}^{'} β

\log it(p_{i}) = z_{i}^{'} γ,

where Y_i is the observed datum for person i (i = 1, 2, …, n). Therefore, we have

\Pr (Y_{i} = 0) = (1 - p_{i}) + p_{i} e^{- λ_{i}} and P r (Y_{i} = y_{i}) = \frac{p_{i} e^{- λ_{i}} λ_{i}^{y_{i}}}{y_{i}^{!}}, y_{i} = 1, 2, \dots

In this model, the logit link is used to predict the probability of subject i (p_i ) to be in the Poisson process. In addition, the log link is used to predict the expected mean (λ_i) of event counts for subject i, given that the subject is not in the perfect zero state, but rather in the Poisson process. x and z are the vectors of covariates hypothesized to explain interindividual differences in the expected counts and interindividual differences in the types of states. x and z may or may not be the same sets of variables. β and γ are the associated regression parameters.

ZIP Model With Random Effects

For handling nested zero-inflated count data, the ZIP models with random effects (the ZIPR models or multilevel ZIP models) were developed and applied to analyze empirical data (e.g., Hall, 2000; Lee et al., 2006). For example, a modified version of the ZIPR model introduced by Lee, Wang, Scott, Yau, and McLachlan (2006) can be expressed as follows:

Y_{i j} ~ {\begin{cases} 0, with probability 1 - p_{i j} \\ P o i s s o n (λ_{i j}), with probability p_{i j} \end{cases}

\begin{array}{l} \log (λ_{i j}) = x_{i j}^{'} β + w_{i} \\ \log it(p_{i j}) = z_{i j}^{'} γ + u_{i}, \end{array}

where x _ij and z _ij represents the jth observations of the ith individual (i = 1, 2, …, n; j = 1, 2, …, m). β and γ are the corresponding fixed-effects regression parameters. The vectors w _i and u _i are the subject random effects or random intercepts in two equations that are assumed to be independently and normally distributed. Hall’s (2000) model is a reduced model of the above model, which contains a random intercept parameter only in the log equation. By including random effects into the ZIP model, we can account for dependence of observations within clusters. However, a limitation of the ZIPR model is that it only analyzes univariate zero-inflated count variable and item information cannot be obtained because no item parameters are included in the model.

Item Response Models

For analyzing multivariate categorical data, item response modeling (e.g., Embretson & Reise, 2000; Lord, 1980) probably is the most widely used method, especially when individual trait scores and item psychometric features are of interest. There is a series of item response models for different kinds of categorical data. Here, we focus only on the 2-parameter (2PL) item response model for binary data, because it is most pertinent to this study. The 2PL model is used to predict the probability (p_ij ) of a person endorsing an item while considering the item-varying parameters, β_j for item location (or difficulty) parameters and α_j for item slope (or discrimination) parameters, which allow for different weights for different items, and the person-varying latent trait variable θ_i. The 2PL model is expressed as

\log it(p_{i j}) = α_{j} (θ_{i} - β_{j}) .

As α_j increases, the item discriminates more precisely among trait levels. When α_j is fixed to be 1, the 2PL model is reduced to be a Rasch model (Rasch, 1960) or a 1PL model. As β_j increases, the item is more and more difficult to endorse. This 2PL model has been shown to be mathematically equivalent to the confirmatory factor analysis model for binary data (Takane & de Leeuw, 1987). The IRT models can be expressed as generalized mixed/multilevel models (e.g., Adams, Wilson, & Wu, 1997; Rijmen, Tuerlinckx, De Boeck, & Kuppens, 2003). Raudenbush, Johnson, and Sampson (2003) also reexpressed the Rasch model and the 2PL model as a two-level logistic model by including dummy variables indicating item numbers.

Proposed Model

In this section, we propose an IRT zero-inflated Poisson model (IRT–ZIP model) to handle multivariate zero-inflated count data. In the model, two sets of item parameters (item slope parameters and item location parameters) are used to explain item characteristics, and a common latent trait variable is included to obtain information on person trait scores. The model is described as follows.

The observed zero-inflated count data indicate

Y_{i j} ~ {\begin{cases} 0, with probability 1 - p_{i j} \\ P o i s s o n (λ_{i j}), with probability p_{i j} \end{cases}

Therefore,

\begin{matrix} P r (Y_{i j} = 0) = (1 - p_{i j}) + p_{i j} e^{- λ_{i j}} \\ P r (Y_{i j} = y_{i j}) = \frac{p_{i j} e^{- λ_{i j}} λ_{i j}^{y_{i j}}}{y_{i j}!}, y_{i j} = 1, 2, \dots . \end{matrix}

The parameters p_ij and λ_ij can be modeled by

log (λ_{i j}) = η_{i j} = \sum_{m = 1}^{j} α_{1 m} I_{i j m} (θ_{i} - \sum_{m = 1}^{j} β_{1 m} I_{i j m}), and

\log it (p_{i j}) = ξ_{i j} = \sum_{m = 1}^{J} α_{2 m} I_{i j m} (θ_{i} - \sum_{m = 1}^{J} β_{2 m} I_{i j m}),

where

I_{i j m} = {\begin{cases} 1 m = j \\ 0 otherwise, with i = 1, \dots, N; j = 1, \dots, J; m = 1, \dots, J . \end{cases}

In the above model, y_ij is the observed count of item j for person i, p_ij is the probability of being in the Poisson process for person i on item j. λ_ij is the expected frequency for person i on item j, given that person i is in the Poisson process. For both the log link equation and the logit link equation, I_ijm is an indicator (dummy) variable to indicate the item number.

From the IRT perspective, θ_i is the latent trait score for person i and the latent variable θ is assumed to have a normal distribution with mean 0 and variance 1 for identification purpose. α_1m and α_2m are item slope/discrimination parameters. Large item slope values indicate that the items are more discriminating in estimating latent trait scores. If we fix both α_1m and α_2m to be 1, then this model is reduced to a model parallel to the Rasch model (Rasch, 1960). β_1m and β_1m are item location/difficulty parameters. As β_1m (location parameter in the log equation) increases, it is more difficult for a person with a given trait score to reach some amount of expected frequency on the mth item, given that the person is in the Poisson process on the mth item. As β_2m (location parameter in the logit equation) increases, it is less likely for a person with a given trait score to be in the Poisson process on the mth item and thus it is more likely for him or her to be in the perfect zero state on the mth item. From the confirmatory factor analysis perspectives, θ can also be viewed as a latent factor. α_1m and α_2m are factor loadings and α_1mβ_1m and α_2mβ_2m are item thresholds. In addition, from the multilevel or mixed-effects modeling framework, θ can be viewed as a random-effects parameter. Therefore, the IRT–ZIP model can also be viewed as a combination of the ZIPR model and the 2PL IRT model, which can be used to analyze item psychometric properties and predict person trait scores from multivariate zero-inflated count data.

Figure 1 displays the bar graphs of 9 items with different item location parameters for β₁ and β₂ and equal item slope parameter values (α = 1). Clearly, different shapes of empirical distributions and different proportions of zeros were obtained by setting different location parameters. For example, 96% of the data were zeros when both β₁ and β₂ are equal to 2 and only 18% of the data were zeros when both β₁ and β₂ are equal to −2. Figure 2 displays the “Item Characteristic Curves” for both λ and p in two equations for 3 items with different α₁ and α₂ values and equal locations parameter values (β= 0). From Figure 2, we can see that with higher α₁ and α₂ values, the items are more discriminating in estimating latent trait scores.

The advantages of the proposed IRT–ZIP model over the composite scores method, the ZIPR model, and the regular 2PL IRT model in analyzing multivariate zero-inflated data include (a) it accounts for excess zeros in count data by considering two processes and using two link functions, (b) it includes item parameters that allow more psychometric information to be obtained from multivariate items, and (c) it includes a latent variable to predict latent trait scores and the distribution of the trait scores is less skewed by specifying a parametric distribution.

This model bears several assumptions. First, the latent variable is one-dimensional and the manifest items in the model measure this latent construct. Second, the latent trait scores have a normal distribution. Third, the non-perfect zero state is a Poisson process. Finally, the items included in the model are locally independent (independent after controlling person trait-level scores and item parameters) and the individuals in the analysis are independent. Some of these assumptions can be relaxed. For example, the latent trait variable θ can be modified to have a log normal distribution. One can also easily extend the unidimensional model to a multidimensional model by including multiple latent trait variables in the model. For example, one can include one latent variable for the log equation and another different latent variable for the logit equation.

The proposed IRT–ZIP model can be viewed as a generalized multilevel model that can be estimated by the marginal maximum likelihood estimation (MMLE) method. For MMLE, the following marginal likelihood of the data is maximized to obtain parameter estimates for the IRT–ZIP model,

L (α, β | y) = \prod_{i = 1}^{N} L_{i} (α, β | y_{i}) = \prod_{i = 1}^{N} \int p (y_{i} | α, β, θ_{i}) N (θ_{i} | α, β) d θ_{i},

where

p (y_{i} | α, β, θ_{i})

is the conditional probability of observing response patterns of y _i treating the latent variable θ_i as known and N() denotes the standard normal distribution. For the log likelihood function, we have

l = \sum_{i = 1}^{N} \log \int \prod_{j = i}^{J} p (y_{i j} | α, β, θ_{i}) N (θ_{i} | α, β,) d θ_{i},

and from Equation 1a, we have the following conditional probability

p (y_{i j} | α, β, θ_{i}) = {[(1 - p_{i j}) + p_{i j} e^{- λ_{i j}}]}^{1 - u_{i j}} {[\frac{p_{i j} e^{- λ_{i j}} λ_{i j}^{y_{i j}}}{y_{i j}!}]}^{u_{i j}} .

By plugging in Equations 1b and 1c, we can reexpress $p (y_{i j} | α, β, θ_{i})$ by the following equation

{(1 + e^{ξ_{i j}})}^{- 1} \times {(1 - u_{i j}) [1 + e^{ξ_{i j}} \exp (- \exp (η_{i j}))] + u_{i j} \frac{e^{ξ_{i j}} \exp [y_{i j} η_{i j} - \exp (η_{i j})]}{y_{i j}!}} .

Here,

u_{i j} = 1 if y_{i j} \neq 0 and u_{i j} = 0

otherwise.

When estimating the parameters, constraints were imposed on the slope parameters (>0). Because the integral for integrating out random effects in Equation 2 generally does not have a closed form, approximation methods are used. For maximizing the approximated marginal likelihood, the Newton-Raphson algorithm (Rabe-Hesketh & Skrondal, 2007), quasi-Newton algorithm (Pinheiro & Bates, 1995), or Expectation-Maximization algorithm (Hall, 2000) can be used. In this study, the MMLE method with adaptive Gaussian-Hermite quadrature and quasi-Newton algorithm (Pinheiro & Bates, 1995) was used to estimate the proposed model. The estimation methods were implemented in SAS PROC NLMIXED (Wolfinger, 1999) because of its flexibility in specifying general likelihood functions such as the one in Equation 2. Example scripts for estimating an IRT–ZIP model in SAS are contained in the appendix.

Because parameters are estimated in the maximum likelihood framework, the likelihood ratio test can be used to compare nested models. For example, to test whether the slope parameters are equal to 1, two nested models—one with freely estimated slopes and one with slopes = 1—can be fitted to the data. The chi-square difference can be calculated from the log-likelihood values of the two models and compared to the critical value from a chi-square distribution with degrees of freedom equal to the difference in the numbers of parameters from two nested models.

Applications to Substance Use Data

Empirical Data Analysis Design

To illustrate the use of the proposed model, the substance use data from the National Longitudinal Survey of Youth (NLSY97 cohort) were analyzed. The NLSY97 consists of a nationally representative sample of approximately 9,000 youths who were 12 to 17 years old as of December 31, 1996. Round 1 survey took place in 1997, and subsequent surveys were conducted annually. The most recent public data are the Round 8 data collected in 2004.

Three areas of substance use (alcohol use, cigarette use, and marijuana use) were targeted to measure the substance use level of each participant. The frequency of use of each substance was measured by three questions. For example, for measuring cigarette use, three questions were asked: (a) Have you ever smoked a cigarette? (Yes or no answer) (b) During the past 30 days, on how many days did you smoke a cigarette? (c) When you smoked a cigarette during the past 30 days, how many cigarettes did you usually smoke each day? We combined the answers to the three questions into one variable for each substance to count the frequency of substance use.¹ And we will refer to each substance use as an item, given 3 items all together in the following data analysis.

Two empirical research questions are to be answered by analyzing the data using the proposed model. These questions are (a) What are the psychometric properties of those items in measuring the latent substance use variable at each occasion? and (b) How do the item parameters change over time for this sample? To answer these questions, the IRT–ZIP model was fitted to the data at each wave. The psychometric features of each item at each occasion are obtained and compared.

Descriptive Statistics

Adolescents (N = 8,984) ranging in age from 12 to 18 (51.19% male) participated in Round 1 of the study in 1997. From Round 1 to Round 8, the means and standard deviations of the age variable and sample sizes are displayed in Table 1. In Round 1, the ethnic composition for this sample was 49.12% non-Hispanic White, 9.12% Hispanic White, 26.58% African American, and 15.18% reported their ethnicity as American Indian, Asian or Pacific Islander, or Other. For the participants' education levels, in Round 8 (2004), 31.35% of the participants had completed high school. 48.98% of the participants were studying in colleges or had completed college studies. And 19.67% of the participants did not complete their high school studies. Table 1 also gives the descriptive statistics of the substance use data over eight waves of the NLSY97 sample. From Table 1, we can find that for each item across the 8 waves, the percentages of zeros were all substantially large, ranging from 36% (alcohol use at Wave 8) to 91% (marijuana use at Wave 1). Generally, from the percentages of substance users and the average counts in Table 1, the substance use propensities and levels increased over time for this NLSY97 sample. In other words, more and more participants became involved with substance use and the level of substance use was generally increasing with age.

Table 1.

Descriptive Statistics of the Substance Use Data of the NLSY97 Sample

	Wave 1	Wave 2	Wave 3	Wave 4	Wave 5	Wave 6	Wave 7	Wave 8
Frequency of ever used/or not (proportion of ever used)
Cigarettes	3,517/5,436	3,037/5,311	3,065/5,084	3,202/4,825	3,260/4,592	3,380/4,499	3,360/4,370	3,211/4,248
Cigarettes	(39.3%)	(36.4%)	(37.6%)	(39.9%)	(41.5%)	(42.9%)	(43.5%)	(43.0%)
Alcohol	3,837/5,107	4,166/4,173	4,403/3,730	4,684/3,338	4,959/2,885	5,321/2,546	5,404/2,323	5,349/2,054
Alcohol	(42.9%)	(50.0%)	(54.1%)	(58.4%)	(63.2%)	(67.6%)	(69.9%)	(72.2%)
Marijuana	1,795/7,148	1,815/6,500	1,879/6,255	2,015/5,992	2,067/5,774	1,924/5,928	1,744/5,971	1,547/5,811
Marijuana	(20.1%)	(21.8%)	(23.1%)	(25.2%)	(26.4%)	(24.5%)	(22.6%)	(21.0%)
Zero days of substance use in the previous 30 days if ever used
Cigarettes	1,899	864	631	576	499	436	449	344
Alcohol	2,094	1,396	1,023	1,041	874	822	756	609
Marijuana	1,015	778	592	559	531	458	412	368
Mean (standard deviation) of use frequencies
Cigarettes	0.90 (3.19)	1.60 (4.71)	2.15 (5.62)	2.59 (6.13)	2.96 (6.60)	3.61 (6.92)	3.30 (6.75)	3.38 (6.90)
Alcohol	0.85 (3.79)	1.51(4.44)	2.05 (5.21)	2.43 (6.07)	2.76 (6.20)	2.82 (5.78)	3.07 (6.37)	3.03 (6.35)
Marijuana	3.72 (9.76)	3.20 (7.23)	2.85 (6.75)	3.35 (7.70)	3.49 (7.90)	3.40 (8.00)	3.89 (8.42)	3.62 (8.11)
Mean (standard deviation) of age
Age	14.31 (1.48)	15.97 (1.44)	16.92 (1.44)	17.98 (1.44)	18.96 (1.44)	19.97 (1.42)	20.93 (1.44)	21.98 (1.43)

Note. NLSY97 = National Longitudinal Survey of Youth.

Results From the IRT–ZIP Models

To investigate psychometric properties of the 3 items and analyze how the item parameters change over time, the IRT–ZIP models were fitted to data of each wave. Table 2 presents the item parameter estimates with fixed slope parameters (=1) and Table 3 displays the item parameter estimates without any parameter constraints. Key findings from Tables 2 and 3 are as follows.

Table 2.

Item Parameter Estimates From Fitting the IRT–ZIP Model With Fixed Slope Parameters to the NLSY97 Data at Each Wave

	Wave 1	Wave 2	Wave 3	Wave 4	Wave 5	Wave 6	Wave 7	Wave 8
Log equations (location parameters)
Cigarettes	−0.15 (0.02)	−0.57 (0.02)	−0.84 (0.02)	−1.02 (0.02)	−1.19 (0.02)	−1.25 (0.02)	−1.31 (0.02)	−1.32 (0.02)
Alcohol	0.03 (0.02)	−0.36 (0.02)	−0.58 (0.02)	−0.78 (0.02)	−0.86 (0.02)	−0.85 (0.02)	−0.91 (0.02)	−0.88 (0.02)
Marijuana	−0.31 (0.03)	−0.54 (0.03)	−0.67 (0.03)	−1.07 (0.02)	−1.19 (0.02)	−1.35 (0.02)	−1.41 (0.02)	−1.40 (0.03)
Logit equations (location parameters)
Cigarettes	1.69 (0.04)	1.20 (0.03)	0.98 (0.03)	0.85 (0.03)	0.73 (0.03)	0.60 (0.03)	0.57 (0.03)	0.56 (0.03)
Alcohol	1.40 (0.04)	0.59 (0.03)	0.14 (0.03)	0.01 (0.03)	−0.36 (0.03)	−0.67 (0.03)	−0.83 (0.03)	−1.03 (0.04)
Marijuana	3.70 (0.06)	3.47 (0.06)	3.34 (0.05)	3.19 (0.05)	3.15 (0.05)	3.27 (0.05)	3.27 (0.05)	3.44 (0.06)
Goodness of fit
−2 Log likelihood	37,861	51,512	60,536	66,096	70,419	72,890	73,779	70,608

Note. IRT = item response theory; NLSY97 = National Longitudinal Survey of Youth; ZIP = zero-inflated Poisson. The values in the parentheses are the standard errors of parameter estimates. The slope parameters are fixed to be 1.

Table 3.

Item Parameter Estimates From Fitting the IRT–ZIP Model With Free Slope Parameters to the NLSY97 Data at Each Wave

	Wave 1	Wave 2	Wave 3	Wave 4	Wave 5	Wave 6	Wave 7	Wave 8
Log equations (slope parameters)
Cigarettes	1.53 (0.05)	1.45 (0.04)	1.30 (0.03)	1.06 (0.02)	1.05 (0.02)	1.01 (0.02)	0.88(0.02)	0.87 (0.02)
Alcohol	1.62 (0.05)	1.34 (0.03)	1.22 (0.03)	1.20 (0.02)	1.02 (0.02)	0.93 (0.02)	0.97(0.02)	0.95 (0.02)
Marijuana	4.27 (0.29)	2.00 (0.11)	1.17 (0.05)	1.21 (0.05)	0.89 (0.04)	1.45 (0.08)	0.78(0.04)	0.81 (0.04)
Log equations (location parameters)
Cigarettes	0.38 (0.04)	−0.01 (0.04)	−0.38 (0.04)	−0.93 (0.05)	−1.05 (0.05)	−1.21 (0.05)	−1.63 (0.06)	−1.66 (0.06)
Alcohol	0.55 (0.03)	−0.00 (0.03)	−0.31 (0.03)	−0.48 (0.03)	−0.84 (0.04)	−1.01 (0.04)	−0.97 (0.04)	−0.96 (0.04)
Marijuana	1.33 (0.04)	0.49 (0.07)	−0.46 (0.09)	−0.69 (0.08)	−1.54 (0.12)	−0.49 (0.10)	−2.17 (0.04)	−2.01 (0.16)
Logit equations (slope parameters)
Cigarettes	1.54 (0.07)	1.57 (.07)	1.50 (0.06)	1.45 (0.06)	1.32 (0.05)	1.19 (0.05)	1.03 (0.05)	1.09 (0.05)
Alcohol	1.27 (0.07)	1.13 (.06)	1.02 (0.05)	0.84 (0.05)	0.70 (0.05)	0.57 (0.05)	0.42 (0.05)	0.25 (0.06)
Marijuana	0.40 (0.26)	1.26 (.10)	1.38 (0.08)	1.04 (0.07)	1.09 (0.06)	0.96 (0.07)	0.98 (0.06)	0.91 (0.07)
Logit equations (location parameters)
Cigarettes	1.23 (0.03)	0.87 (0.03)	0.75 (0.03)	0.71 (0.03)	0.62 (0.03)	0.55 (0.03)	0.59 (0.03)	0.56 (0.03)
Alcohol	0.98 (0.04)	0.41 (0.04)	0.03 (0.04)	−0.19 (0.05)	−0.58 (0.07)	−1.13 (0.12)	−2.04 (0.29)	−4.09 (0.96)
Marijuana	4.12 (1.48)	2.72 (.13)	2.64 (0.10)	3.04 (0.15)	2.96 (0.14)	3.30 (0.18)	3.36 (0.18)	3.74 (0.23)
Goodness of fit
−2 Log likelihood	36,959	51,114	60,320	65,904	70,302	72,692	73,604	70,371

Note. IRT = item response theory; NLSY97 = National Longitudinal Survey of Youth; ZIP = zero-inflated Poisson. The values in the parentheses are the standard errors of parameter estimates.

First, based on the likelihood ratio test, the models without any parameter constraints significantly fitted the data better than the models with parameter constraints for all eight waves. This indicates that the discrimination ability of different items in estimating trait levels differs by items. However, this conclusion is tempered by the fact that the differences in the fit between the two models are not very large compared to the large sample size in the study.

Second, the order of magnitude of the location parameter estimates was examined within each wave. The results showed that the orders of magnitude of the location parameter estimates were similar across eight waves and between the two nested models. For example, for the location parameter estimates of the log equations, the estimates of the alcohol use items were larger than those of the cigarette use items and the marijuana use items. This indicates that it is less likely for a person with a fixed substance use level to reach the same amount of expected frequency in alcohol use as in cigarette use and marijuana use. This result was consistent with our descriptive statistics in the rows of “use frequencies” in Table 1.² For the location parameters of the logit equations, the estimates of the marijuana use items are all larger than those of the cigarette use items, while the estimates of the cigarette use items are larger than those of the alcohol use items. This indicates that it is less likely for a person with a fixed substance use level to be in the Poisson process on marijuana use than on cigarette use and on alcohol use. In other words, it is less likely on average for a person to be involved in marijuana use than cigarette use or alcohol use. And this result is also consistent with the descriptive statistics in the first six rows in Table 1.

Third, after comparing the item location parameter estimates across eight waves, we find that the three estimates of the log equations are generally getting smaller and smaller, which means that overall the sample had increasing frequencies in substance use. This result is also consistent with the descriptive statistics on “use frequencies” in Table 1. The location parameter estimates of the logit equations for the cigarette use item and alcohol use item are generally decreasing over time, which also means that the participants had increasing likelihoods of becoming involved with substance use or being substance users. The location parameter estimates of the logit equation for marijuana use fluctuated over time in both Tables 2 and 3. The fluctuation in the magnitudes was consistent with the fluctuation in the percentages of participants who ever used marijuana in Table 1.

Figure 3 displays two plots of location parameter estimates from Table 2 across different waves: one for the location parameters of the log equations and the other for the location parameters of the logit equations. From the plots, we can also clearly see the pattern of the order of magnitude within a wave and the change pattern of the item parameter estimates across waves.

In this empirical study, the benefits of having answers from three questions instead of just one question is that we can distinguish whether a participant is a smoker or a nonsmoker by looking at the answer to Question 1 (therefore, we know whether the zero is from a perfect zero state or the Poisson process), whether a participant smoked during those 30 days if he or she is a smoker by looking at the answer to Question 2, and how many cigarettes a participant smoked each day if he or she smoked during the past 30 days by looking at the answer to Question 3. In this way, we can empirically examine the performance of the proposed model on analyzing psychometrics properties of the items by comparing the parameter estimates to the descriptive statistics in Table 1 (see the second and third results above). We can also examine the performance of the proposed model in distinguishing two processes by comparing the predictions from the model to the raw data. The accuracies of predictions from the proposed model were reported in Table 4. For example, using the data from the last wave (year 2004), the predictions from the IRT–ZIP model (with free slope parameters) were compared to the raw data. First, for the predictions on whether the zero is from the perfect zero state or the Poisson process, we found that the proportions of accurate predictions were 72%, 72%, and 79% for the smoking, alcohol use, and marijuana use items, respectively. In terms of predicting the counts of substance use, given that the process is the Poisson process, the correlations between the predictions (predicted λ values) and the real counts from data (excluding those who never used the substance because the process is the Poisson process) are .72, .83, and .54 for the smoking, alcohol use, and marijuana use items, respectively. The prediction results showed that the model can be used to distinguish the two processes and estimate expected counts well.

Table 4.

Accuracies of Predictions From Fitting the IRT–ZIP Model With Free Slope Parameters to the NLSY97 Data at Each Wave

	Wave 1	Wave 2	Wave 3	Wave 4	Wave 5	Wave 6	Wave 7	Wave 8
Accuracies of predictions on states
Cigarettes	70.0%	76.9%	76.0%	75.2%	74.1%	73.3%	70.9%	72.1%
Alcohol	69.7%	73.5%	78.7%	78.4%	85.0%	67.7%	70.0%	72.3%
Marijuana	79.9%	78.2%	76.9%	74.8%	73.6%	75.5%	77.4%	79.5%
Accuracies of predictions on expected counts
Cigarettes	.66	.79	.80	.78	.74	.81	.77	.72
Alcohol	.81	.73	.77	.81	.81	.76	.84	.83
Marijuana	.53	.58	.56	.54	.55	.57	.51	.54

Note. IRT = item response theory; NLSY97 = National Longitudinal Survey of Youth; ZIP = zero-inflated Poisson.

Simulation Study

Simulation Study Design

The purpose of this simulation is to investigate the accuracy and precision of the item parameter estimates under the empirical data analysis scenario and the effects of different sample sizes on the parameter estimates. Data were simulated based on true values adapted from the empirical data analysis parameter estimates of wave 4 (α₁₁ = 1.15, β₁₁ = −0.80, α₂₁ = 1.40, β₂₁ = 0.8; α₁₂ = 1.15, β₁₂ = −0.50, α₂₂ = 0.90, β₂₂ = 0; α₁₃ = 1.30, β₁₃ = −0.5, α₂₃ = 1.15, β₂₃ = 3) with varying sample sizes (N = 200, 500, 1,000, and 2,000). The number of items is three. The proportions of zeros of these three items are 71%, 58%, and 95%. To evaluate the results, both accuracy and precision of parameter estimates were examined. For evaluating accuracy, the means of the empirical distribution of the parameter estimates were compared to the true parameter values and the Root Mean Square (RMS) statistic was calculated for each parameter ( $RMS = \sqrt{\frac{1}{T} \sum_{t = 1}^{T} {(ϑ^{(t)} - ϑ_{0})}^{2}},$ where T is the number of replications, $ϑ^{(t)}$ is the estimate of the tth replication and $ϑ_{0}$ is the true value). For evaluating precision, the means of the estimated standard errors were compared to the standard deviations of the empirical distribution of the estimates. Performance of the proposed model and the composite score method on estimating person trait scores also was compared by investigating the differences between the estimated trait scores and the true trait scores. To quantify the differences, Pearson correlation coefficients between the estimated trait scores and the true values were calculated. In addition, predictions on the probabilities of being in the Possion state (p_i ) and the expected counts (λ_i), given that the person is in the Poisson process are evaluated by examining the Pearson correlations between the predictions from the model and the true values. The number of replications in the simulation study is 200.

Simulation Results

Table 5 displays the item parameter estimates with varying sample sizes. In terms of the accuracy of item parameter estimates, the results showed that the parameters cannot be recovered accurately or precisely when the sample size is relatively small (e.g., N = 200), especially for the parameter β₂₃. However, with adequate sample size (e.g., N = 1,000), the item parameters can be recovered quite well (the maximum relative absolute bias is 3.3% (β₂₃ with N = 1,000). The accuracy of the parameter estimates (e.g., RSM) generally was improved with larger sample sizes. In terms of the precision of item parameter estimates, from Table 5, one can see that (a) the precision of parameter estimates was higher with larger sample sizes and (b) the average parameter standard error estimates (s.e.2) were close to the standard deviations of the parameter estimates (s.e.1) when the sample size is adequate (e.g., N = 1,000). In the empirical study, the sample size is larger than 7,000. Therefore, one can conclude that the item parameters are likely to be estimated accurately and precisely in the empirical study.

Table 5.

Item Parameter Estimates From the MMLE Method by Sample Size for the Simulation

	True Value	N = 200				n = 500				n = 1,000				n = 2,000
	True Value	Mean	RSM	s.e.1	s.e.2	Mean	RSM	s.e.1	s.e.2	Mean	RSM	s.e.1	s.e.2	Mean	RSM	s.e.1	s.e.2
α₁₁	1.15	1.13	0.14	0.14	0.14	1.13	0.10	0.09	0.10	1.13	0.07	0.06	0.06	1.14	0.04	0.04	0.04
α₁₂	1.15	1.13	0.13	0.14	0.13	1.13	0.09	0.08	0.09	1.13	0.07	0.06	0.07	1.14	0.04	0.04	0.04
α₁₃	1.3	1.26	0.33	0.29	0.33	1.27	0.15	0.15	0.15	1.27	0.10	0.10	0.10	1.28	0.07	0.07	0.06
α₂₁	1.4	1.47	0.39	0.39	0.39	1.43	0.22	0.23	0.22	1.42	0.15	0.16	0.15	1.41	0.09	0.11	0.09
α₂₂	0.9	0.97	0.38	0.36	0.37	0.90	0.22	0.21	0.23	0.88	0.16	0.15	0.16	0.89	0.12	0.10	0.12
α₂₃	1.15	1.11	0.46	0.49	0.46	1.13	0.25	0.29	0.25	1.13	0.19	0.20	0.19	1.14	0.15	0.14	0.15
β₁₁	−0.8	−0.85	0.28	0.27	0.27	−0.84	0.18	0.17	0.17	−0.83	0.12	0.11	0.12	−0.81	0.08	0.08	0.08
β₁₂	−0.5	−0.54	0.22	0.21	0.22	−0.53	0.14	0.13	0.14	−0.54	0.12	0.09	0.11	−0.51	0.07	0.06	0.06
β₁₃	−0.5	−0.71	0.68	0.59	0.65	−0.55	0.27	0.27	0.27	−0.54	0.18	0.18	0.18	−0.53	0.12	0.12	0.12
β₂₁	0.8	0.8	0.19	0.19	0.19	0.81	0.11	0.12	0.11	0.8	0.08	0.08	0.08	0.8	0.06	0.06	0.06
β₂₂	0	−0.11	0.76	0.66	0.75	−0.01	0.19	0.18	0.19	0.00	0.11	0.12	0.11	−0.01	0.08	0.08	0.08
β₂₃	3	3.42	1.21	1.51	1.14	3.19	0.67	0.69	0.65	3.10	0.38	0.43	0.37	3.09	0.30	0.30	0.29

Note. MMLE = marginal maximum likelihood estimation RSM = root mean square. Results represent 3 items corresponding to the second subscript on each parameter. The proportions of zeros of these three items are 69%, 58%, and 94%, respectively. When implementing the estimation, starting values were all set to be 1 for the slope parameters and 0 for the location parameters. s.e. 1 is the mean of the standard errors of parameter estimates from different replications, and s.e. 2 is the standard deviation of the parameter estimates from different replications.

Pearson correlations were calculated between the person trait estimates and the true values when the sample size is 2,000. The average correlation was .81 (SD = 0.01). The average correlation between the composite scores and the true values also was calculated and was .54 (SD = 0.05). The correlations indicated that the proposed method should be more appropriate than the composite score method for estimating the substance use of youths.

In terms of accuracies of predictions on the probabilities of being in the Poisson state (p_i ) and the expected counts (λ_i ), Pearson correlations between the predictions from the model and the true values were calculated and averaged across different replications with the sample size of 2,000. The correlations for the probabilities were .85, .79, and .89, respectively, and the correlations for the expected counts are .94, .92, and .93 for the 3 items, respectively. The correlations from the simulation studies are higher than the corresponding correlations from the empirical data mainly because that the simulated data were simulated from the proposed model and thus followed the model assumptions more strictly.

Discussion and Conclusion

This study introduced an IRT–ZIP model to analyze multivariate zero-inflated count data using two link functions: a log link function and a logit link function. By including two link functions, the probability of being in the Poisson process and the expected frequency of an event given that participant is in the Poisson process can be modeled simultaneously. A latent person trait variable is included in the model to measure latent trait scores. Both item slope and location parameters are included in the model to investigate psychometric properties of the multiple items.

The model was applied to analyze the substance use data to illustrate the applications. The empirical results provided rich information regarding both the substance use of adolescence and the model performance. The order of magnitude of the item location parameter estimates in the logit equations indicates that it is less likely for a person with a fixed substance use level to be involved with marijuana use than cigarette use and alcohol use. This result is partly consistent with the “gateway hypothesis” in substance use research (Kandel, 2003; Morral, McCaffrey, & Paddock, 2002). Kandel (2003) described the three propositions in the “gateway hypothesis,” which are sequencing (fixed relationship between two substances such that one substance is regularly initiated before the other), association (initiation of one substance increases the likelihood of initiation of the second substance), and causation (a controversial proposition, use of the first substance causes the use of the second substance). The application of the IRT–ZIP model to the substance use data in this study is based on the association proposition and the results were consistent with the sequencing proposition from this naturalistic population example with a novel method. Overall, the associations and differences among the 3 items measuring latent substance use levels were investigated successfully using the IRT–ZIP model. The consistency between some of the current findings and previous findings supports the reasonableness and reliability of the results from the IRT–ZIP model.

Simulation studies also were conducted to examine the performance of the model under the empirical scenario. The simulation results showed that, similar to regular IRT models (Thissen & Wainer, 1982, 2001), fitting the IRT–ZIP model requires large sample sizes to obtain stable item parameter estimates. For a 2PL IRT model, rules of thumb for the minimum number of examinees required for accurate parameter estimation range from 500 (Hulin, Lissak, & Drasgow, 1982) to 1,000 (Ree & Jensen, 1980). For the proposed IRT–ZIP model, in light of the tested simulation conditions, a sample size of 1,000 seems large enough to recover the parameter estimates. For large-scale psychological and educational observational or survey research, this might not be a big issue. However, in many practical psychological and educational experiments, only a small number of participants' data are collected and thus possibly preclude or at least limit the use of the proposed method. In terms of recovering the person trait scores, the performance of the IRT–ZIP model was better than the composite score method. This is because the IRT–ZIP model considers the natural features of the multivariate zero-inflated count data. By including two link functions, we can model two processes more accurately and validly. The model thus can be useful in large-scale psychological or educational settings where count data on the number of behaviors are available.

Limitations and Future Directions

There are several issues that should be investigated in the future. First, the empirical data characteristics include a large sample size (N > 7,000) and repeated measures making them particularly apt for illustrating the application of the proposed model. However, the data were limited to only 3 items. With 3 items, the correlation between the person trait estimates and the true values was only about .81 under the empirical data scenario, in accord with the simulation results. Further simulations should examine the performance of the model on estimating both item parameters and trait scores with larger numbers of items and different proportions of zeros. Examples of applying the proposed models to problem behaviors or other psychological research with a larger number of items can also be included in future studies. Second, because the empirical data contain repeated measurements, change patterns of substance use levels from adolescence to young adulthood could be investigated here. The proposed IRT–ZIP model can be extended to a three-level mixed-effects model, which involves estimating item parameters, latent person trait scores, and longitudinal change patterns simultaneously in future studies. Third, in the current proposed model, only one common latent variable is included. However, it is possible to include two latent variables to investigate the propensity and level scores individually. Future studies can also be done to investigate the multidimensionality.

In summary, both the simulation results and the empirical results supported the conclusion that the proposed IRT–ZIP model is a valid and feasible method for analyzing multivariate zero-inflated count data more precisely and reasonably. The work so far should encourage researchers to apply the IRT–ZIP model to more substantive research areas and thereby answer more interesting empirical research questions.

Footnotes

Notes

Acknowledgment

The author is grateful for the helpful comments on this study from John R. Nesselroade, John J. McArdle, Steven M. Boker, Karen M. Schmidt, Xiaohui Wang, and the reviewers.

Appendix

References

Adams

R. J.

Wilson

M. L.

(1997). Multilevel item response models: An approach to errors in variables regression. Journal of Educational and Behavioral Statistics, 22, 47–76.

Böckenholt

Kamakura

W. A.

Wedel

(2003). The structure of self-reported emotional experiences: A mixed-effects Poisson factor model. British Journal of Mathematical & Statistical Psychology, 56, 215–229.

Cameron

Trivedi

P. K.

(1998). Regression analysis of count data, Econometric Society Monograph No. 30. Cambridge, England: Cambridge University Press.

Cattell

R. B.

(1978). The scientific use of factor analysis. New York, NY: Plenum.

Cox

D. R.

Lewis

P. A. W.

(1966). The statistical analysis of series of events. (Methuen Monograph). London, England: Methuen.

Embretson

S. E.

Reise

S. P.

(2000). Item response theory for psychologists. Mahwah, NJ: Erlbaum.

Frome

E. L.

Kutner

M. H.

Beauchamp

J. J.

(1973). Regression analysis of Poisson-distributed data. Journal of the American Statistical Association, 68, 935–940.

Fuller

W. A.

(1987). Measurement error models. New York, NY: Wiley.

Gart

J. J.

(1964). The analysis of Poisson regression with an application in virology. Biometrika, 51, 517–521.

10.

Hall

D. B.

(2000). Zero-inflated Poisson and binomial regression with random effects: A case study. Biometrics, 56, 1030–1039.

11.

Hulin

C. L.

Lissak

R. I.

Drasgow

(1982). Recovery of two- and three- parameter logistic item characteristic curves: A Monte Carlo study. Applied Psychological Measurement, 6, 249–260.

12.

Kandel

D. B.

(2003). Does marijuana use cause the use of other drugs?. Journal of the American Medical Association, 289, 482–483.

13.

Lambert

(1992). Zero-inflated Poisson regression, with an application to the defects in manufacturing. Technometrics, 34, 1–14.

14.

Lee

A. H.

Wang

Scott

J. A

Yau

K. K. W.

McLachlan

G. J.

(2006). Multi-level zero-inflated Poisson regression modelling of correlated count data with excess zeros. Statistical Methods in Medical Research, 15, 47–61.

15.

Liu

Powers

D. A.

(2007). Growth curve models for zero-inflated count data: An application to smoking behavior. Structural Equation Modeling, 14, 247–279.

16.

Lord

F. M.

(1980). Applications of item response theory to practical testing problems. Hillsdale, NJ: Erlbaum.

17.

Lord

F. M.

Novick

M. R.

(1968). Statistical theories of mental test scores. Reading, MA: Addison-Wesley.

18.

McCullagh

Nelder

J. A.

(1989). Generalized linear models. 2nd ed. London, England: Chapman and Hall.

19.

McDonald

R. P.

(1985). Factor analysis and related methods. Hillsdale, NJ: Erlbaum.

20.

Min

Y. Y.

Agresti

(2005). Random effect models for repeated measures of zero-inflated count data. Statistical Modeling, 5, 1–19.

21.

Morral

A. R.

McCaffery

D. F.

Paddock

S. M.

(2002). Reassessing the marijuana gateway effect. Addiction, 97, 1493–1504.

22.

Pinheiro

J. C.

Bates

D. M.

(1995). Approximation to the log-likelihood function in the nonlinear mixed-effects model. Journal of Computational and Graphical Statistics, 4, 12–35.

23.

Rabe-Hesketh

Skrondal

(2007). Multilevel and latent variable modeling with composite links and exploded likelihoods. Psychometrika, 72, 123–140.

24.

Rasch

(1960). Probabilistic models for some intelligence and attainment tests. Chicago, IL: University of Chicago Press.

25.

Raudenbush

S. W.

Johnson

Sampson

(2003). A multivariate, multilevel Rasch model with applications to self-reported criminal behavior. Sociological Methodology, 33, 169–211.

26.

Ree

M. J.

Jensen

H. E.

(1980). Effects of sample size on linear equating of item characteristic curve parameters. In Weiss

D. J.

(Ed.), Proceedings of the 1979 computerized adaptive testing conference. Minneapolis: University of Minnesota.

27.

Rijmen

Tuerlinckx

De Boeck

Kuppens

(2003). A nonlinear mixed model framework for item response theory. Psychological Methods, 8, 185–205.

28.

Takane

de Leeuw

(1987). On the relationship between item response theory and factor analysis of discretized variables. Psychometrica, 52, 393–408.

29.

Thissen

Wainer

(1982). Some standard errors in item response theory. Psychometrica, 47, 397–412.

30.

Thissen

Wainer

(2001). Test scoring. Hillsdale, NJ: Erlbaum.

31.

Wolfinger

(1999). Fitting nonlinear mixed models with the new NLMIXED procedure. In SUGI 24 Conference Proceedings, Paper 287. Cary, NC: SAS Institute.