2.1 Univariate count data

Count data may assume any nonnegative integer value and are typically highly skewed. Count data can be empirically (unconditionally) equidispersed when the empirical mean and variance are equal, empirically underdispersed when the empirical mean exceeds the empirical variance, or empirically overdispersed when the empirical variance exceeds the empirical mean. Realizations of a count variable may be directly observed or unobserved (latent) and may be measured with error in either case (e.g., Cameron & Trivedi, Reference Cameron and Trivedi2013). For example, a pedometer counts steps indirectly as a function of movement (so that step count is latent) and the resulting counts are subject to measurement error arising from the imperfect mapping between detected movement and steps taken. Alternatively, counts of grammatical morphemes in a video-recorded and transcribed oral language sample are directly measured with what is probably a very small amount of error. Count data are also at the ratio level of measurement defined by ordered categories with equal intervals and true zero, where the latter represents an absence of the measured count variable. This is in contrast to the more commonly analyzed ordinal data, whose numeric values, including zero, have no inherent meaning but rather indicate the relative position/ordering of the various levels of the ordinal variable. In a measurement context in which a latent variable gives rise to the observed count variable, a count of zero may represent an absence of the underlying latent variable in addition to an absence of the observed count variable, depending on what the latent variable represents and the probability distribution it is assumed to follow. Consider a latent variable representing symptom severity, for instance. In this scenario, a count of zero may represent absence of the symptom. Alternatively, for a normally distributed latent variable representing an individual’s level of expressive language development, a count of zero may correspond to levels of development falling below a certain threshold along the latent continuum.

In this article, we restrict our attention to count data arising from a single counting process, although count data from multiple response processes can be accommodated. A counting process is a stochastic process describing the nonnegative integer number of events we expect to observe within a given exposure, where the number of observed events cannot decrease with increasing exposure. The exposure quantifies the length of time, space, or number of trials over which events are recorded and must be either a positive real number or a (positive) natural number (e.g., Cameron & Trivedi, Reference Cameron and Trivedi1998, Reference Cameron and Trivedi2013; Hilbe, Reference Hilbe2011). An exposure that may vary across observations (e.g., individuals, measurement occasions) yields an exposure variable. Like counts themselves, an exposure variable (here: the number of utterances sampled from a child at an assessment that defines the length of the oral language sample) may be either directly observed or latent and may be measured with error in either case. When an exposure variable cannot be measured directly, is multi-faceted, and/or is not well-defined, it can sometimes be reconstructed (possibly with error) as a function of a set of measured variables. Measurement error in counts and exposure are considered at length by Cameron and Trivedi (Reference Cameron and Trivedi2013, Chapter 13).

Despite the variability and/or measurement error that are commonly present in exposure, probability distributions for count data implicitly assume an exposure that is fixed and measured consistently across observations. When a consistently measured, fixed exposure is used to collect each observation:

1. (a) All observed responses are on the same scale (a necessary condition for obtaining correct values and interpretations for model parameters, including the event rate and expected count [the mean]);

2. (b) The properties of the maximum likelihood estimators of model parameters are unaffected; and

3. (c) All variability in the observed count response variable is attributable to sampling variability arising from individual differences in model parameters (e.g., Cameron & Trivedi, Reference Cameron and Trivedi2013; Hilbe, Reference Hilbe2011).

In count models, explicitly and properly accounting for an exposure that varies across the units/observations comprising a sample is critical to obtaining correct inferences by achieving (a) and (b) and parsing variability in (c) from sampling variability in the observed responses due to varying exposure (e.g., Cameron & Trivedi, Reference Cameron and Trivedi1998). Failing to properly account for a varying exposure will yield biased parameter estimates, standard errors, and model fit statistics—leading to incorrect inferences—as neither (a), (b), nor (c) will hold. Note, however, that including an exposure variable in a count model will only yield correct inferences if the probability of observing an event per unit exposure is constant (e.g., Cameron & Trivedi, Reference Cameron and Trivedi1998). The opportunity to explicitly incorporate an exposure variable into the probability distribution governing the counting process presents itself through the regression of the distribution mean parameter onto a set of covariates that includes the exposure.

2.2 Probability distributions for a single counting process

Distributions modeling a single counting process lie within the exponential family (summarized in Table 2). For each of these distributions, the mean parameter quantifies the expected response and the variance is expressed as a function of the mean parameter, possibly in addition to a dispersion or “nuisance” parameter. Depending on the distribution, the mean parameter may also represent the expected event rate (e.g., the Poisson and Negative Binomial distributions). Note that different variance functions yield different probability distributions, and different distributions (models) imply different relationships between the mean and variance. For example, the Poisson distribution (de Moivre, Reference de Moivre1711, Reference de Moivre1718; Poisson, Reference Poisson1837) implies a variance that equals the mean (model-implied [conditional] equidispersion). Negative Binomial (NB) distributions (e.g., Cameron & Trivedi, Reference Cameron and Trivedi2013; Hilbe, Reference Hilbe2011) imply a variance that exceeds the mean (model-implied [conditional] overdispersion). Meanwhile, Katz (e.g., Katz, Reference Katz1963), Double Poisson (DP; e.g., Efron, Reference Efron1986), Generalized Poisson (GP; e.g., Consul & Famoye, Reference Consul and Famoye1992; Consul & Jain, Reference Consul and Jain1973; Consul, Reference Consul1989), and Conway–Maxwell–Poisson (CMP; e.g., Conway & Maxwell, Reference Conway and Maxwell1962; Guikema & Coffelt, Reference Guikema and Coffelt2008; Huang, Reference Huang2017; Minka et al., Reference Minka, Shmueli, Kadane, Borle and Boatwright2003; Shmueli et al., Reference Shmueli, Minka, Kadane, Borle and Boatwright2005) distributions can imply a variance that is less than the mean (model-implied [conditional] underdispersion) as well as a variance that equals or exceeds the mean. With that said, distributions other than the Poisson and Negative Binomial with quadratic variance function (denoted as NB2; see Table 2) suffer from notable limitations that have restricted their use in the literature to date.

Table 2 Probability distributions for a single counting process

Note: All models are members of the exponential family and have support $\mathbb {N}_{0}$ . 1PEF denotes the one-parameter exponential family. 2PEF denotes the two-parameter exponential family. Y denotes the count response variable. $\mathbb {E}\left (Y\right )$ denotes the model-implied (conditional) mean. $\mathbb {V}\left (Y\right )$ denotes the model-implied (conditional) variance, which is also called the variance function as the variance is a function of the mean. For all models, $\phi $ denotes the dispersion or “nuisance” parameter, where the interpretation and parameter space of $\phi $ depends on the distribution. In the Poisson, NB, and $\text {CMP}_{\mu }$ distributions, $\mu $ is the mean parameter. In the DP distribution, the parameter $\widetilde {\mu }$ is “similar to the mean parameter” (Cameron & Trivedi, Reference Cameron and Trivedi2013). In the CMP and $\text {CMP}_{\mu }$ distributions, $\lambda $ is the rate parameter. For the NB family of distributions, the variance is a function of the mean and dispersion parameters (i.e., $\mathbb {V}\left (Y\right )=\omega (\mu ,\phi )$ ). NB1 denotes the Negative Binomial distribution with linear variance function $\mu + \phi \mu $ , NB2 denotes the Negative Binomial distribution with quadratic variance function $\mu + \phi \mu ^{2}$ , and NBp denotes the Negative Binomial distribution with variance function $\mu + \phi \mu ^{p}$ .

2.3 Regression models

Incorporating a regression model for the mean parameter of a counting process distribution permits the expected response to be a function of exposure and other predictors. For counting processes in the two-parameter exponential family (2PEF), a regression model for the dispersion parameter may also be specified, where the total set of regression equations specified may or may not have predictors or coefficients in common. Regression models for a count response arising from a single counting process distribution for which an exact expression for the mean response is available, such as the Poisson and NB distributions (see Table 2), can be formulated within a generalized linear model (GLM), generalized nonlinear model (GNLM), generalized linear mixed model (GLMM), or generalized nonlinear mixed model (GNLMM) framework, depending on whether the relations among predictors are linear or nonlinear and whether the coefficients of the predictors are fixed, random, or some combination thereof (e.g., Agresti, Reference Agresti2013; Cameron & Trivedi, Reference Cameron and Trivedi2013; Fitzmaurice et al., Reference Fitzmaurice, Laird and Ware2011; Hilbe, Reference Hilbe2011; McCullagh & Nelder, Reference McCullagh and Nelder1989; Nelder & Wedderburn, Reference Nelder and Wedderburn1972; Vonesh, Reference Vonesh2012). With that said, few examples exist in the literature of a nonlinear growth function of linear and nonlinear random effects connected to observed counts via a nonlinear link function due to the computational difficulties that may arise in model estimation.

For a GLM, GNLM, GLMM, or GNLMM regressed on the mean parameter of a counting process distribution: (1) the systematic component may include an offset—the natural log of an exposure variable with corresponding regression coefficient fixed at one; (2) the random component is parameterized in terms of the mean and dispersion of the counting process (see Table 2); and (3) a natural log link function is used to connect the random and systematic components. Although other link functions that ensure the mean is always positive may be used with count responses (e.g., Wedel et al., Reference Wedel, Böckenholt and Kamakura2003), the natural log link function is generally preferred as it is the canonical link function for the Poisson mean parameter. Note that including an offset in the systematic component of the regression on the mean response allows the mean to be expressed as the product of exposure and the hazard rate quantifying the expected count per unit exposure, where either or both may be functions of observed and/or latent variables (e.g., Cameron & Trivedi, Reference Cameron and Trivedi2013).

Note that one may observe equidispersion, overdispersion, or underdispersion under a fitted model when the empirical (observed) variance equals, exceeds, or is less than the model-implied (conditional) variance, respectively. (The discerning reader may observe the distinction between under-, equi-, and overdispersion that is empirical, model-implied, or arising when a model is fit to data is rarely made explicit in the literature.) The presence of under- or overdispersion is understood to yield incorrect standard errors, incorrect model fit statistics, and, as a result, fallacious inferences with real, practical implications (e.g., Cameron & Trivedi, Reference Cameron and Trivedi2013; Hilbe, Reference Hilbe2011). Under- or overdispersion arising from model misspecification may additionally yield incorrect inferences about the response process but can be addressed through careful reconsideration of the systematic component of the regression model for the mean response and subsequent selection of a more appropriate conditional probability distribution (see Table 2).

Given the computational challenges involved in implementing other counting process distributions enumerated in Table 2, Poisson and NB2 distributions remain popular choices across a variety of scientific applications. Furthermore, given the frequency with which overdispersion is observed in practice, NB2 models present an attractive choice for use in modeling change processes measured by the repeated measurement of a single count variable over time. Using a conditional response distribution in the 2PEF additionally permits investigation of the joint behavior of the mean response and unexplained variability among responses (manifesting as overdispersion) over time by fitting a first-order LGM to the mean responses and a separate trajectory (that need not be a GLM or GNLM) to the time-varying dispersion parameters. As such, although the proposed first-order LGM may be specified using any counting process distribution in the exponential family, we illustrate our approach using the NB2 distribution.

2.4 First-order multiphase SLCM for count data

Typically, first-order LGMs are those that model a single observed indicator repeatedly measured at a set of time points or occasions for a sample of individuals. We begin this section by explicating the notation used to define the response data, the measurement occasions, and individuals. For additional clarity, we couch explication of the notation in the empirical example of modeling counts of Brown’s (Reference Brown1973) second grammatical morpheme (BGM2) “in” (see Table 1) collected in the longitudinal assessment of GAE morphosyntactic development over the course of early childhood (Figure 1).

Let $y_{iw_{i}}$ denote the number of “in” morphemes produced by child i out of $x_{iw_{i}}$ utterances sampled at chronological age $t_{w_{i}}$ months. As implied by the notation, the chronological ages (in months) at which children were assessed varied across children both within and among corpora (see Table 3 and Figure 2). This is because, in each corpus, oral language was sampled at certain, planned chronological ages that differed across corpora, and children within a corpus were sometimes also assessed at slightly different ages than planned and/or were missing planned assessments.Footnote ¹ As implied by the notation $x_{iw_{i}}$ , the number of sampled utterances (exposure) varied also across children and assessments (Figure 3). A spaghetti plot of the rates at which the morpheme “in” is produced within an oral language sample in the overall sample of children is given in Figure 2 with a lowess smooth of the mean function superimposed.

Table 3 Sample characteristics

Note: CHILDES denotes the Child Language Data Exchange System (MacWhinney, Reference MacWhinney2000); DOI denotes Digital Object Identifier; n denotes the number of children sampled from the indicated corpus; NR denotes Not Reported; and $W_{i}$ denotes the number of measurements taken over time from individual (child) i. In the Activity column, TP denotes oral language data were collected through engagement in a toy play activity, G denotes a group activity, M denotes a meal activity, B denotes a book activity, and N denotes a narrative activity. Samples drawn from listed corpora contained typically developing (TD) monolinguistic native speakers of General American English (GAE) aged [18, 72) months (i.e., aged [1.5, 6) years) with at least 25 Mean Length of Utterance (MLU)-eligible sampled utterances.

Figure 2 Rate at which the morpheme “in” is produced within an oral language sample by chronological age and corpus. Note: Brown’s (Reference Brown1973) grammatical morphemes are summarized in Table 1. Sample characteristics are provided in Table 3.

Figure 3 Relationship between chronological age and number of sampled utterances. Note: Sample characteristics are provided in Table 3.

Let $\mathbf {y}_{i}=\{y_{iw_{i}}: w_{i} = 1,\dots ,W_{i}\}$ , $\mathbf {x}_{i}=\{x_{iw_{i}}: w_{i} = 1,\dots ,W_{i}\}$ , and $\mathbf {t}_{i}=\{t_{iw_{i}}: w_{i} = 1,\dots ,W_{i}\}$ denote the sets of observed responses (BGM2 counts), exposures (numbers of sampled utterances), and measurement times (chronological ages of assessment in months), respectively, collected from the repeated measurement of a single count item/indicator (i.e., the number of “in” morphemes produced) over the course of $W_{i}$ occasions for individual i, where the observed counts are nonnegative integers and the exposures and measurement times are positive real numbers. The number of assessments per child ( $W_{i}$ ) ranged from one to seven, inclusive, with corresponding ages of assessment ( $t_{iw_{i}}$ ) ranging from 18 to 71 months, inclusive (Table 3). Chronological age of assessment was binned into 3-month assessment windows (intervals), yielding a total of $W=18$ unique measurement occasions in the overall sample of $N=1,084$ children: [18,21), [21,24), [24,27), [27,30), [30,33), [33,36), [36,39), [39,42), [42,45), [45,48), [48,51), [51,54), [54,57), [57,60), [60,63), [63,66), [66,69), and [69,72). In the analysis dataset, these measurement windows were coded using the lower bound of each 3-month interval, centered at age 18 months, yielding

$$ \begin{align*} \mathbf{t} = \{t_{w}: w=1,\dotsc,W\} = \{0, 3, 6, 9, 12, 15, 18, 21, 24, 27, 30, 33, 36, 39, 42, 45, 48, 51\}, \end{align*} $$

so that w indexes unique measurement occasion (unique chronological age of assessment) within the overall sample of N individuals. As such, data collected within a given interval are treated as though collected at the chronological age indicated by the lower bound of that interval in the analysis and interpretation of model parameters. For example, BGM2 counts collected at ages [18, 21) months are treated as though collected at age 18 months. Since not all children were assessed within each 3-month interval, the ages of assessment for child i ( $\mathbf {t}_{i}$ ) are a subset of the unique ages of assessment in the overall sample ( $\mathbf {t}$ ), such that

$$ \begin{align*} \mathbf{t}_{i} = \{t_{iw_{i}}: w_{i}=1,\dotsc,W_{i} \leq W=18\} \in \mathbf{t}, \quad \text{where} \hspace{0.2cm} i=1,\dotsc,N. \end{align*} $$

If a child had more than one assessment within a 3-month interval, one of those assessments was randomly selected for analysis so that each child contributed at most one assessment per interval (i.e., data are cross-sectional within each 3-month assessment window). Note that the use of age intervals spanning a few months is a fairly common practice in the clinical assessment of expressive language development because of dynamic and rapid growth in grammar acquisition during the earliest stages of child language development (e.g., Carrow-Woolfolk, Reference Carrow-Woolfolk2011; Leadholm & Miller, Reference Leadholm and Miller1992; Miller & Chapman, Reference Miller and Chapman1981; Pavelko & Owens, Reference Pavelko and Owens2017; Rice et al., Reference Rice, Smolik, Perpich, Thompson, Rytting and Blossom2010; Scarborough et al., Reference Scarborough, Rescorla, Tager-Flusberg, Fowler and Sudhalter1991; Sparrow et al., Reference Sparrow, Cicchetti and Saulnier2016). Here, 3-month intervals were used in order to balance (a) choosing intervals small enough to maximize the amount of longitudinal data taken from each child, the number of unique measurement occasions in the overall sample, and granularity with respect to chronological age with (b) choosing intervals large enough to contain enough children to permit the estimation of dispersion across children within each age interval.

2.4.1 Measurement model

Due to the myriad of components unique to the modeling of count response data as well as the complexity of the notation, we present an example of the proposed first-order multiphase LGM as a path diagram in Figure 4. Following the typical structural equation modeling convention for path diagrams, squares represent observed variable indicators, circles represent latent variables, single-headed arrows denote directed relations among observed and latent variables, and double-headed arrows denote variances of and covariances between variables. Other notation germane to the explication of the measurement and structural components of the LGM in Figure 4 is detailed and explained shortly.

Figure 4 Path diagram for a first-order linear–linear latent growth model fit to repeated measurements of a single count response variable that conditionally follows a distribution in the two-parameter exponential family at each measurement occasion. Note: A solid, single-line, black arrow indicates a structural relationship with an identity link function. A solid, single-line, red arrow indicates a structural relationship with a non-identity (e.g., natural log) link function. A dashed black arrow from A to B indicates A gives rise to B directly and/or indirectly. A solid, double-line, black arrow from A to B indicates A generates B.

Suppose $Y_{iw_{i}}$ measures individual i’s level of target construct $\theta $ (e.g., level of morphosyntactic development) at time $t_{iw_{i}}=t_{w}$ , where corresponding observed response $y_{iw_{i}}$ is generated from a counting process distribution in the 2PEF with density

(1) $$ \begin{align} f_{2\text{PEF}}\left(y_{iw_{i}} \mid \eta_{iw_{i}}, \phi_{w} \right) = \text{exp}\left\{ a(\eta_{iw_{i}},\phi_{w}) + b(y_{iw_{i}},\phi_{w}) + c(\eta_{iw_{i}},\phi_{w})t(y_{iw_{i}}) \right\} \end{align} $$

with time-specific dispersion parameter $\phi _{w}$ and natural/canonical parameter $\eta _{iw_{i}} = \ln \mu _{iw_{i}}$ , where $\mu _{iw_{i}}$ represents individual i’s expected response at occasion $w_{i}$ . For the empirical example, we use the $\text {NB2}(\mu _{iw_{i}},\phi _{w})$ density in Equation (2).

(2) $$\begin{align} f_{\text{NB2}}\left(y_{iw_{i}} \mid \mu_{iw_{i}}, \phi_{w} \right) &= \frac{\Gamma\left(y_{iw_{i}} + \phi_{w}^{-1}\right)}{\Gamma\left(y_{iw_{i}} + 1\right)\Gamma\left(\phi_{w}^{-1}\right)} \left[\frac{\phi_{w}^{-1}}{\phi_{w}^{-1} + \mu_{iw_{i}}}\right]^{\phi_{w}^{-1}} \left[\frac{\mu_{iw_{i}}}{\phi_{w}^{-1} + \mu_{iw_{i}}}\right]^{y_{iw_{i}}} \end{align}$$

Individual i’s expected response at chronological age $t_{iw_{i}}=t_{iw}$ months is expressed as a function (Equation (3)) of exposure $x_{iw_{i}}$ , individual i’s levels of the target construct $\theta _{iw_{i}}$ , item intercept $\xi _{0}$ , item slope $\xi _{1}$ , and error $\epsilon _{iw_{i}}$ .

(3) $$\begin{align} \eta_{iw_{i}} = \ln \mu_{iw_{i}} &= \ln x_{iw_{i}} + \xi_{0} + \xi_{1}\theta_{iw_{i}} + \epsilon_{iw_{i}} \end{align}$$

Note the subscript $w_{i}$ in $\theta _{iw_{i}}$ indicates Equation (3) applies to those chronological ages at which child i is actually assessed, though the child’s latent construct exists continuously across all chronological ages (including $\mathbf {t}$ ), measured or not, so that $\boldsymbol {\theta }_{i} = (\theta _{i1},\dotsc ,\theta _{iW})^{\top }$ . The item intercept quantifies the log expected response (expected count) per unit exposure at the population average level of the target construct (i.e., when $\theta _{iw_{i}}=0$ ), while the item slope captures the strength of the (positive) linear relation between $\theta _{iw_{i}}$ and log expected response, $\ln \mu _{iw_{i}j}$ . Note that the expected count per unit exposure increases multiplicatively by an order of $\text {exp}\left (\xi _{0}\right )$ as $\xi _{0}$ increases, holding $\theta _{iw_{i}}$ constant. Also note that item parameters $\boldsymbol {\xi }=\{\xi _{k} : k=0,1\}$ do not vary across individuals or measurement occasions to achieve measurement invariance.

The error, $\epsilon _{iw_{i}}$ , in the linear predictor of the mean response in Equation (3) may be non-zero due to the omission of important predictors of the mean response, misspecification of the structural model, and/or measurement error in the offset, $\ln x_{iw_{i}}$ . This error propagates down to the observed data level manifesting as overdispersion (unexplained variability) in the observed responses (e.g., Cameron & Trivedi, Reference Cameron and Trivedi2013). The set of time-specific errors for individual i, $\boldsymbol {\epsilon }_{i} $ , is normally distributed with zero mean vector and covariance matrix $\boldsymbol {\Omega }_{i}$ .

$$\begin{align*} \boldsymbol{\epsilon}_{i} = (\epsilon_{i1}, \dotsc, \epsilon_{iW_{i}})^{\top} \overset{iid}{\sim} \mathcal{N}\left(\mathbf{0},\boldsymbol{\Omega}_{i}\right) \end{align*}$$

The set of time-specific latent constructs for individual i follows a different multivariate normal distribution with mean vector $\boldsymbol {\theta }$ and covariance matrix $\boldsymbol {\Theta }$ .

$$\begin{align*} \boldsymbol{\theta}_{i} = (\theta_{i1},\dotsc,\theta_{iW})^{\top} \overset{iid}{\sim} \mathcal{N}\left(\boldsymbol{\theta}, \boldsymbol{\Theta} \right) \end{align*}$$

2.4.2 Structural model

Change in an individual’s level of the target construct over time follows a theoretically defensible growth model, f, expressed as a function of measurement times, $\mathbf {t}$ , and individual latent growth factors, $\boldsymbol {\beta }_{i}$ .

(4) $$\begin{align} \boldsymbol{\theta}_{i} &= f\left(\mathbf{t}, \boldsymbol{\beta}_{i}\right) + \boldsymbol{\delta}_{i} \end{align}$$

The set of disturbances, $\boldsymbol {\delta }_{i} = \{\delta _{iw} : w=1,\dotsc ,W\}$ , in Equation (4) represents the set of time-specific regression errors induced by misspecifying the true (data-generating) trajectory of $\boldsymbol {\theta }_{i}$ over $\mathbf {t}$ . This could occur by omitting predictors of an individual’s level of the target construct and/or misspecifying the functional form of f. The disturbances are assumed to jointly follow a multivariate normal distribution with zero mean vector and covariance matrix $\boldsymbol {\Sigma }$ .

$$\begin{align*} \boldsymbol{\delta}_{i} =\left(\delta_{i1},\dotsc,\delta_{iW}\right)^{\top} \overset{iid}{\sim} \mathcal{N}\left(\mathbf{0},\boldsymbol{\Sigma}\right) \end{align*}$$

In Equation (5), we assume that f is a piecewise growth model having the general form,

(5) $$ \begin{align} f\left(\mathbf{t}, \boldsymbol{\beta}_{i}\right) &= \begin{cases} f_{1}\left(\mathbf{t}, \boldsymbol{\beta}_{i}\right), & \quad \mathbf{t} \le \gamma_{1i} \\ f_{2}\left(\mathbf{t}, \boldsymbol{\beta}_{i}\right), & \quad \gamma_{1i} < \mathbf{t} \le \gamma_{2i} \\ \quad \vdots & \quad \vdots \\ f_{D}\left(\mathbf{t}, \boldsymbol{\beta}_{i}\right), & \quad \mathbf{t}> \gamma_{(D-1)i} \end{cases}. \end{align} $$

The trajectory in each of the $D> 1$ phases may have a distinct functional form that need not be a polynomial. Moreover, the trajectory need not be monotonic within a phase nor over the entire measurement period. The changepoints (a.k.a., join points or knots, denoted by $\boldsymbol {\gamma }_{i}=\{\gamma _{id} : d=1,\dotsc ,D-1\}$ ) indicate the times of transition from phase to phase. These transition times may be unknown (i.e., parameters to be estimated from the data) and may vary across individuals (i.e., have random effects). The transition from one phase to the next may be discontinuous (e.g., a jump up or a drop down), continuous but abrupt (zero-order continuity), or continuous and gradual/smooth (first-order continuity or greater, where higher orders of continuity correspond to greater degrees of smoothness).

Individual growth factors $\boldsymbol {\beta }_{i} = (\boldsymbol {\psi }_{i}^{\top }, \boldsymbol {\varphi }_{i}^{\top })^{\top }$ include p parameters, $\boldsymbol {\psi }_{i}$ , that enter the function in Equation (5) linearly and q parameters, $\boldsymbol {\varphi }_{i}$ , that enter the function in a nonlinear manner. Here, we define a growth parameter as being linear if the first partial derivative of f taken with respect to the growth parameter does not include the parameter. Alternatively, a growth parameter is considered to be nonlinear if the first partial derivative of f taken with respect to the growth parameter includes the parameter. Nonlinear growth factors include—but need not be limited to—unknown, individual-specific changepoints.

Individual linear and nonlinear growth factors may each be expressed as a linear combination of population growth parameters governing the population average trajectory—linear and nonlinear fixed effects $\boldsymbol {\psi }$ and $\boldsymbol {\varphi }$ —and individual i’s linear and nonlinear random effects, $\mathbf {u}_{i}$ and $\mathbf {g}_{i}$ , respectively. Random effects $\mathbf {b}_{i} = (\mathbf {u}_{i}^{\top }, \mathbf {g}_{i}^{\top })^{\top }$ are assumed to jointly follow a multivariate normal distribution with zero mean vector and symmetric covariance matrix T.

(6) $$\begin{align} \boldsymbol{\beta}_{i} = \begin{pmatrix} \boldsymbol{\psi}_{i} \\ \boldsymbol{\varphi}_{i} \end{pmatrix} = \begin{pmatrix} \boldsymbol{\psi} \\ \boldsymbol{\varphi} \end{pmatrix} + \begin{pmatrix} \mathbf{u}_{i} \\ \mathbf{g}_{i} \end{pmatrix}, \quad \text{where } \begin{pmatrix} \mathbf{u}_{i} \\ \mathbf{g}_{i} \end{pmatrix} \overset{iid}{\sim} \mathcal{N} \begin{pmatrix} \begin{bmatrix} \textbf{0} \\ \textbf{0} \end{bmatrix}, \; \begin{bmatrix} \textbf{T}_{uu} & \textbf{T}_{ug} \\ \textbf{T}_{gu} & \textbf{T}_{gg} \end{bmatrix} \end{pmatrix} \end{align}$$

Note that elements of T may be constrained at zero for theoretical but also potentially computational (pragmatic) reasons.

Following S. A. Blozis & Harring (Reference Blozis and Harring2016) and S. A. Blozis & Harring (Reference Blozis and Harring2017), the growth function for individual i in Equation (5) can be reformulated as a SLCM defined by a first-order Taylor series expansion taken with respect to the parameters of the mean growth function and linearly weighted by a set of individual-specific weights (i.e., random effects, $\textbf {b}_{i}$ ).

(7) $$\begin{align} f\left(\mathbf{t}, \boldsymbol{\beta}_{i}\right) &\approx f(\mathbf{t}, \boldsymbol{\beta}) + \boldsymbol{\Lambda}(\mathbf{t}, \boldsymbol{\beta}) \mathbf{b}_{i} \end{align}$$

The columns (i.e., basis functions) of $\boldsymbol {\Lambda }(\mathbf {t}, \boldsymbol {\beta })$ are the the first-order partial derivatives of the mean (i.e., target) function, $f(\mathbf {t}, \boldsymbol {\beta })$ .

(8) $$\begin{align} \boldsymbol{\Lambda}= \boldsymbol{\Lambda}(\mathbf{t}, \boldsymbol{\beta}) = \frac{\partial f(\mathbf{t}, \boldsymbol{\beta})}{\partial\boldsymbol{\beta}} =\left( \frac{\partial f(\mathbf{t}, \boldsymbol{\beta})}{\partial{\psi}_{1}} \cdots \frac{\partial f(\mathbf{t}, \boldsymbol{\beta})}{\partial{\psi}_{p}} \hspace{0.2cm} \frac{\partial f(\mathbf{t}, \boldsymbol{\beta})}{\partial{\varphi}_{1}} \cdots \frac{\partial f(\mathbf{t}, \boldsymbol{\beta})}{\partial{\varphi}_{q}} \right) \end{align}$$

As a SLCM, $f(\mathbf {t}, \boldsymbol {\beta })$ is assumed to be invariant to a constant scaling factor (see Shapiro & Browne, Reference Shapiro and Browne1987, Condition 2). Consequentially, there is a set of parameters, denoted here by $\boldsymbol {\alpha }$ , such that

(9) $$ \begin{align} f(\mathbf{t}, \boldsymbol{\beta}) = \boldsymbol{\Lambda} \boldsymbol{\alpha}. \end{align} $$

Because $\boldsymbol {\Lambda }$ is the set of first-order partial derivatives of $f(\mathbf {t}, \boldsymbol {\beta })$ taken with respect to $\boldsymbol {\beta }$ (see Equation (8)), elements of parameter vector $\boldsymbol {\alpha }$ can be obtained by solving the linear equations in Equation (9). It turns out that solving these linear equations results in setting all parameters in $\boldsymbol {\alpha }$ that enter nonlinearly to 0 (i.e., $\boldsymbol {\varphi }=\mathbf {0}$ ). This permits the recovery of the target function (Preacher & Hancock, Reference Preacher and Hancock2015). Then in the individual-level model in Equation (7), $\boldsymbol {\Lambda }\boldsymbol {\alpha }$ can be substituted for the mean function, $f(\mathbf {t}, \boldsymbol {\beta })$ . Thus, the individual-level model can be re-expressed as

$$ \begin{align*} f\left(\mathbf{t}, \boldsymbol{\beta}_{i}\right) &= \boldsymbol{\Lambda} \boldsymbol{\alpha} + \boldsymbol{\Lambda} \textbf{b}_{i} = \boldsymbol{\Lambda} \boldsymbol{\eta}_{i}, \end{align*} $$

where $\boldsymbol {\eta }_{i} = \boldsymbol {\alpha } + \textbf {b}_{i}$ .

Note that when a single count indicator is measured at each occasion, the item intercept is constrained at zero $\left (\xi _{0}=0\right )$ and item slope at one $\left (\xi _{1}=1\right )$ to achieve model identification while preserving meaningful interpretation of the growth parameters.

(10) $$\begin{align} \eta_{iw_{i}} = \ln \mu_{iw_{i}} &= \ln x_{iw_{i}} + \xi_{0} + \xi_{1}\left[ f\left(t_{iw_{i}}, \boldsymbol{\beta}_{i}\right) + \delta_{iw_{i}}\right] + \epsilon_{iw_{i}} \nonumber \\ &= \ln x_{iw_{i}} + f\left(t_{iw_{i}}, \boldsymbol{\beta}_{i}\right) + \varepsilon_{iw_{i}}, \quad \text{where} \hspace{0.2cm} t_{iw_{i}}=t_{w} \end{align}$$

Additionally, regression disturbance $\delta _{iw_{i}}$ is absorbed into measurement error $\epsilon _{iw_{i}}$ , impacting the interpretation of $\varepsilon _{iw_{i}}$ and any resulting overdispersion.

$$\begin{align*} \boldsymbol{\delta}_{i} + \boldsymbol{\epsilon}_{i} = \boldsymbol{\varepsilon}_{i} =(\varepsilon_{i1},\dotsc,\varepsilon_{iW_{i}})^{\prime} \overset{iid}{\sim} \mathcal{N}\left(\mathbf{0},\boldsymbol{\Xi}_{i}\right) \end{align*}$$

2.4.3 Dispersion trajectory

If a conditional distribution in the 2PEF, such as the NB2 distribution, is utilized for the count response measured repeatedly over time, the magnitude of conditional (model-implied) dispersion may notably vary over time and possibly even follow it’s own trajectory. Let $\boldsymbol {\phi }=\{\phi _{w} : w=1,\dotsc , W \}$ denote the complete set of freely estimated dispersion parameters corresponding to the W unique measurement times $\mathbf {t}=\{t_{w} : w=1,\dotsc , W \}$ in the overall sample of N individuals. Where appropriate, one may fit a trajectory with coefficients $\boldsymbol {\omega }$ to time-varying dispersion parameters $\boldsymbol {\phi }$ to describe change in dispersion over time.

(11) $$\begin{align} \boldsymbol{\phi} = g\left(\mathbf{t},\boldsymbol{\omega}\right) \end{align}$$

The trajectory fit to the time-specific dispersion parameters need not be a GLM nor utilize polynomial growth, but whatever function is ultimately used, it must capture the essential characteristics of the freely estimated dispersion parameters across time. For example, in modeling counts of Brown’s (Reference Brown1973) second grammatical morpheme (BGM2) “in” (see Table 1) collected in the longitudinal assessment of GAE morphosyntactic development over the course of early childhood (Figure 1), change in model-implied overdispersion over time might be described by a linear function fit to the natural log of the NB2 dispersion parameters

(12) $$ \begin{align} \ln \phi_{w} = \omega_{1} + \omega_{2} \hspace{0.05cm} t_{w}, \end{align} $$

where $w=1,\dotsc ,18$ and $t_{w}=0,\dotsc ,51$ . Alternatively, a more precipitous decline over the first several months of early childhood might be achieved through an exponential decay function with a nonzero asymptote (see Equation (13)) to describe change in dispersion over time, where $\omega _{1}+\omega _{3}$ quantifies the dispersion parameter at age 18 months when $t_{w}=0$ , $\omega _{2}$ is the decay factor such that the NB2 dispersion parameter decreases by $100(1-\omega _{2})\%$ with every 1 month increase in chronological age after age 18 months, and $\omega _{3}$ is the nonzero asymptote quantifying the NB2 dispersion parameter as children age beyond early childhood,

(13) $$ \begin{align} \phi_{w} = \omega_{1}\left(\omega_{2}^{t_{w}}\right) + \omega_{3}, \quad \text{where} \hspace{0.2cm} w=1,\dotsc,18 \hspace{0.2cm} \text{and} \hspace{0.2cm} t_{w}=0,\dotsc,51. \end{align} $$

Note that the dispersion parameters are time-specific (as indicated by the subscript w) but do not vary across individuals (as indicated by the omission of i from the subscript). As such, any trajectory fit to the dispersion parameters exists only at the population level and is specified by imposing constraints on the time-varying dispersion parameters, reducing the number of freely estimated dispersion-specific parameters from W to the length of $\boldsymbol {\omega }$ . As such, fitting a trajectory to the dispersion parameters may convey notable parsimony in addition to the ability to describe a change process of substantive interest, where gains in parsimony for a given function g increase as the number of unique measurement occasions (W) in the overall sample of N individuals increases.

3.1 Model assumptions

As with any fully parametric LGM in which parametric probability distributions are assumed for both the latent variables (random effects) and observed variables (indicators/items), the following assumptions are implied when fitting the proposed first-order multiphase SLCM to longitudinal count data. First, it is assumed that the model is correctly specified (e.g., Agresti, Reference Agresti2013; Cameron & Trivedi, Reference Cameron and Trivedi1998, Reference Cameron and Trivedi2013; Hilbe, Reference Hilbe2011; McCullagh & Nelder, Reference McCullagh and Nelder1989; McNeish & Kelley, Reference McNeish and Kelley2019; Vonesh, Reference Vonesh2012; Woods & Thissen, Reference Woods and Thissen2006), including correct specification of:

1. 1. The joint distribution of the random effects;

2. 2. The fixed and random effects included in the linear predictor of the mean response and the relations among them;

3. 3. The conditional distribution assumed for the response; and

4. 4. The link function connecting the mean response to its linear predictor.

Additionally, several assumptions are made about the target population and sample of individuals from whom the count responses are collected. First, the sample from which model parameters are to be estimated is assumed to be both homogeneous (i.e., all sampled individuals come from the same population; OECD, 2004) and representative (i.e., the sample is selected probabilistically and the composition of the sample is “typical” of the population with respect to certain, specified characteristics of interest upon which inferences will be based; OECD, 2002). Second, sampled individuals are assumed to be independent and sampled measurements (count responses) are assumed to be conditionally independent within an individual (e.g., McCulloch, Reference McCulloch2003; Vonesh, Reference Vonesh2012; Woods & Thissen, Reference Woods and Thissen2006).

Several measurement-specific assumptions are also made. First, measurement invariance is assumed across individuals and occasions. Second, the IRT/IFA assumption of monotonicity applies here, which posits that the probability of endorsing a given response category or higher increases as the level of the latent construct measured by the item increases. For a count response, the assumption of monotonicity implies that the expected response (expected count) increases as the level of the latent construct increases. Third, as measurement error is not the focus of this research, it is assumed that count responses and exposures are directly measured with, at most, minimal error that is uncorrelated with predictors of the mean response (e.g., Cameron & Trivedi, Reference Cameron and Trivedi2013).

Lastly, combining the model-, sample-, and measurement-specific assumptions enumerated above, we assume errors are uncorrelated (mutually independent) among individuals at a given measurement time as well as across occasions within each individual after conditioning on the growth trajectory, so that $\varepsilon _{iw_{i}} \overset {iid}{\sim } \mathcal {N}(0,\sigma _{w}^{2})$ in Equation (10) at time $t_{iw_{i}}=t_{w}$ (such as is shown in Figure 4) and

$$ \begin{align*} \boldsymbol{\Xi} &= \begin{pmatrix} \sigma_{1}^{2} & & & \\ & \sigma_{2}^{2} & & \\ & & \ddots & \\ & & & \sigma_{W}^{2} \end{pmatrix}. \end{align*} $$

3.2 Model fit assessment

The goal of model fit assessment is to identify model assumptions that appear to be notably violated based on available empirical data (that are, hopefully, homogeneous and representative of the target population). Identifying sources of model-data misfit can, in conjunction with theoretical considerations, inform re-specification of the model such that more valid—and therefore more useful—inferences about the data-generation process may be drawn.

3.3 Model identification

The mean structure of the proposed first-order multiphase SLCM fit to count responses is comprised of all freely estimated population growth parameters in $\boldsymbol {\beta }$ . For conditional response (counting process) distributions in the one-parameter exponential family (e.g., the Poisson distribution), the covariance structure is comprised of all freely estimated, unique growth factor variances and covariances (freely estimated unique elements of $\textbf {T}$ ). For conditional response distributions in the 2PEF (e.g., the NB2 distribution), the covariance structure additionally includes occasion-specific dispersion parameters $\boldsymbol {\phi } = \{\phi _{w}: w=1,\dotsc ,W\}$ or, if a growth trajectory is imposed on $\boldsymbol {\phi }$ , the parameters of said trajectory (i.e., $\boldsymbol{\omega}$ ). Note that the elements of factor loading matrix $\boldsymbol {\Lambda }$ are not freely estimated but rather are functions of measurement times $\mathbf {t}$ and freely estimated population growth parameters in $\boldsymbol {\beta }$ . Likewise, expected counts $\boldsymbol {\mu }= \{\mu _{iw_{i}}: w_{i}=1,\dotsc , W_{i}; \hspace {0.1cm} i=1,\dotsc , N\}$ , linear predictors of the mean response $\boldsymbol {\eta }= \{\eta _{iw_{i}}: w_{i}=1,\dotsc , W_{i}; \hspace {0.1cm} i=1,\dotsc , N\}$ , and error variances $\boldsymbol {\sigma }^{2}=\{\sigma _{w}^{2}: w=1,\dotsc , W\}$ are part of the hypothesized model but are not freely estimated model parameters.

Since the proposed first-order SLCM for count responses is a CFA model with a mean structure (e.g., Browne, Reference Browne1993; Kline, Reference Kline2016), the number of observations is $W(W+3)/2$ , where, as previously noted, W is the number of observed (count) variables or, equivalently for this model, the number of unique measurement occasions in the overall sample of N individuals (e.g., Kline (Reference Kline2016, Rule 15.5)). Additionally, since “the identification status of a mean structure must be considered separately from that of the covariance structure” (Kline, Reference Kline2016), the mean and covariance structures must each be overidentified in order for the model as a whole to be overidentified. To ensure the mean structure of the proposed model is overidentified, the total number of unique measurement occasions W in the sample of N individuals must exceed the number of freely estimated population growth parameters in $\boldsymbol {\beta }$ (e.g., Kline, Reference Kline2016). Likewise, to ensure the covariance structure is overidentified, $W(W+1)/2$ must exceed the number of freely estimated unique elements of $\mathbf {T}$ and—for conditional response distributions in the 2PEF—the number of occasion-specific dispersion parameters (W) or dispersion-specific regression coefficients.

3.4 Model estimation

The estimation of latent variable models (and GLMMs) for count responses can be challenging because each variable in the systematic component for the mean response has a nonlinear relationship with the conditional mean due to the use of a non-identity link function (e.g., Olsen & Schafer, Reference Olsen and Schafer2001; Vonesh, Reference Vonesh2012). As a result, there is generally no closed form solution to either the marginal log-likelihood or marginal moments, so that contemporary model estimation approaches aim to either:

1. 1. Approximate the marginal moments of an approximate quasi-likelihood function or approximate the marginal quasi-likelihood function corresponding to specified first- and possibly also second-order conditional moments through linearization (via Taylor series expansion);

2. 2. Maximize the marginal log-likelihood via numerical integration or simulation; or

3. 3. Approximate the posterior distribution for the model parameters given the observed data via fully Bayesian approaches, where—like the log-likelihood—the posterior distribution typically does not have a closed form solution.

These different estimation strategies have different analytic objectives and require different assumptions. A thorough treatment of these various approaches to model estimation may be found in, for example, Vonesh (Reference Vonesh2012), Bolker et al. (Reference Bolker, Brooks, Clark, Geange, Poulsen, Stevens and White2009), and Hoff (Reference Hoff2009).

Of these different approaches, the most popular by far for latent variable models for count responses has been marginal maximum likelihood (MML) estimation implemented via numerical integration (e.g., Beisemann, Reference Beisemann2022; Beisemann et al., Reference Beisemann, Forthmann and Doebler2024; Forthmann & Doebler, Reference Forthmann and Doebler2021; Forthmann et al., Reference Forthmann, Gühne and Doebler2020; Hung, Reference Hung2012; Jansen, Reference Jansen1994, Reference Jansen1995; Jansen & van Duijn, Reference Jansen and van Duijn1992; H. Liu, Reference Liu2007; Magnus & Thissen, Reference Magnus and Thissen2017; Rabe-Hesketh et al., Reference Rabe-Hesketh, Skrondal and Pickles2004; Shiyko et al., Reference Shiyko, Li and Rindskopf2012; Wang, Reference Wang2010). Fully parametric MML via numerical integration can yield consistent and asymptotically unbiased parameter estimates, even when data are missing at random (MAR; Rubin, Reference Rubin1976) or otherwise unbalanced (e.g., Asparouhov & Muthén, Reference Asparouhov and Muthén2012; De Boeck & Wilson, Reference De Boeck and Wilson2004; Gunes & Chen, Reference Gunes and Chen2014; H. Liu, Reference Liu2007; Rabe-Hesketh et al., Reference Rabe-Hesketh, Skrondal and Pickles2004; Vonesh, Reference Vonesh2012) or when counts are small or underdispersed ( and therefore manifestly more discrete; Vonesh, Reference Vonesh2012). In addition, MML permits the computation of likelihood-based information criteria (e.g., the LRT, AIC, and BIC), which remain the most powerful diagnostic tools available for use with multivariate count data, greatly facilitating the detection of model-data misfit (e.g., Magnus & Thissen, Reference Magnus and Thissen2017; Vonesh, Reference Vonesh2012). Lastly, although MML estimation via numerical integration can be computationally intensive, it need not be prohibitively so. Smith & Blozis (Reference Smith and Blozis2014), for example, demonstrate that very few quadrature nodes may be needed, especially when adaptive Gauss–Hermite (AGH) quadrature is used.

4.1 Data

Expressive language data were taken from oral language sampled from 1,084 typically developing monolinguistic native speakers of GAE aged 1.5 to 5 years, inclusive (i.e., aged $[18,71]$ months), drawn from across 23 corpora in the Child Language Data Exchange System (CHILDES; MacWhinney, Reference MacWhinney2000), accessible at https://childes.talkbank.org/. Sample characteristics are provided in Table 3.

The number of children sampled from each corpus varied notably across corpora. The percentages of male and female children also varied across corpora but were comparable in the overall sample (46.86% male, 45.20% female, and 7.93% not reported). Oral language was sampled through engagement in either a toy play, narrative, group, book, or meal activity, where the range of activities varied across and sometimes also within corpora. The number of sampled utterances (exposure) varied across corpora, children within corpora, and assessments within child. Sampled assessments were required to contain a minimum of 25 utterances of sufficient quality (i.e., at least 25 “mean length of utterance (MLU)-eligible” utterances) to ensure a child had sufficient opportunity to demonstrate his/her level of morphosyntactic development. No notable relationship was discerned between a child’s chronological age in months (strongly related to level of morphosyntactic development in typically developing young children) and the number of sampled utterances in either the overall sample or within corpora (Figure 4), facilitating stable model estimation, the selection of an appropriate model (through stable parameter estimates, accurate standard errors, and the resulting apparent significance of individual parameter estimates), and the interpretation of estimated growth parameters as intended.

4.2 Analysis

Expressive language data were retrieved from North American English corpora in CHILDES (see Table 3) at TalkBank.org and extracted using Computerized Language ANalysis (CLAN) software (MacWhinney, Reference MacWhinney2000). To inform the selection of an appropriate conditional response distribution (and thus measurement model), several analyses were conducted, where candidate distributions included the Poisson and NB2.Footnote ² First, the empirical mean and variance of the number of “in” morphemes produced were computed within each 3-month age interval and plotted over time (i.e., over the course of early childhood) to check for empirical under-, equi-, or overdispersion. Second, Poisson and NB2 models were fit to the cross-sectional data within each 3-month age bracket in R (R Core Team, Reference Team2023) using the glm function in the stats package and glm.nb function in the MASS package (Venables & Ripley, Reference Venables and Ripley2002), respectively, with the number of sampled utterances (the exposure) included in the linear predictor of the mean response (Equation (14)). The relative fit of the Poisson and NB2 models was evaluated by visually comparing corresponding KLD, Pearson, AIC, and BIC statistics within each 3-month age interval. Note that likelihood-based information criteria are being used to compare the overall fit of the Poisson and NB2 models instead of the LRT as the Poisson distribution is not a special case (i.e., a constrained version) of the NB2 distribution. Rather, it is a limiting distribution as the strictly positive NB2 dispersion parameter tends to zero (e.g., Casella & Berger, Reference Casella and Berger2002).

(14) $$\begin{align} \hat{\eta}_{iw_{i}} &= \widehat{\ln \mu}_{iw_{i}} = \ln x_{iw_{i}} + \hat{\xi}_{0w}, \quad \text{where} \hspace{0.2cm} t_{iw_{i}}=t_{w} \end{align}$$

After selecting a measurement model, scatterplots and cubic smoothing splines of the cross-sectionally estimated model-implied mean log expected BGM2 counts and, if applicable, dispersion parameters were plotted over time to inform the selection of the functional form of the population average trajectory and, if applicable, the dispersion trajectory, respectively. A multiphase SLCM was subsequently fit to the longitudinal data in Mplus ^© Version 8.5Footnote ³ using MML estimation with parameter estimates and corresponding standard errors robust to both nonnormality of the count responses and dependence among observations by specifying TYPE=COMPLEX and ESTIMATOR=MLR in the ANALYSIS command in conjunction with use of the CLUSTER option in the VARIABLE command (Muthén & Muthén, Reference Muthén and Muthén2017). Robust standard errors were computed using a sandwich estimator and—by default—the observed information matrix evaluated at the maximum likelihood estimates of model parameters (i.e., the Hessian matrix; Muthén & Muthén, Reference Muthén and Muthén2017). MML estimation was implemented via the Expectation-Maximization (EM) algorithm (e.g., Bock & Aitkin, Reference Bock and Aitkin1981) with adaptive Gauss-Hermite quadrature (e.g., Cai, Reference Cai2010; Rabe-Hesketh et al., Reference Rabe-Hesketh, Skrondal and Pickles2004; Schilling & Bock, Reference Schilling and Bock2005; Vonesh, Reference Vonesh2012; Wang, Reference Wang2010) with the default of 15 integration points per dimension of integration.

Because the number of sampled utterances (the exposure) varied across children and assessments within child, each measurement model fit to BGM2 counts throughout this analysis included an offset in the linear predictor of the mean response (the expected BGM2 count) constructed as the natural log of the number of utterances sampled from a child at an assessment with regression coefficient fixed at one. An alpha level of 0.05 was used for all hypothesis tests. All tables and figures were produced in R version 4.3.1 (R Core Team, Reference Team2023). Likewise, all data processing/management and descriptive and exploratory analyses were conducted in version 4.3.1 (R Core Team, Reference Team2023). Mplus ^© model results were harvested and read into R using the MplusAutomation package (Hallquist & Wiley, Reference Hallquist and Wiley2018) and analyzed using the MASS package (Venables & Ripley, Reference Venables and Ripley2002). All datasets, computer code, and Supplementary Material are provided on the OSF website for this project at: https://osf.io/j6rp7/?view_only=c6a9f74add0d476dbeb1e7abd6a76cb2.