Roberto Di MariZsuzsa BakkJennifer OserJouni Kuha Department of Economics and
Business, University of Catania, ItalyDepartment of Methodology and
Statistics, Leiden University, The NetherlandsDepartment of Politics and
Government, Ben-Gurion University, IsraelDepartment of Statistics, London
School of Economics and Political Science, London, UK Corresponding author:
Roberto Di Mari
Address: Corso Italia, 55, 95128, Catania, Italy.
Email<EMAIL_ADDRESS>
# Multilevel latent class analysis with covariates: Analysis of cross-national
citizenship norms with a two-stage approach

###### Abstract

This paper focuses on the substantive application of multilevel LCA to the
evolution of citizenship norms in a diverse array of democratic countries. To
do so, we present a two-stage approach to fit multilevel latent class models:
in the first stage (measurement model construction), unconditional class
enumeration is done separately on both low and high level latent variables,
estimating only a part of the model at a time – hence keeping the remaining
part fixed – and then updating the full measurement model; in the second stage
(structural model construction), individual and/or group covariates are
included in the model. By separating the two parts – first stage and second
stage of model building – the measurement model is stabilized and is allowed
to be determined only by it’s indicators. Moreover, this two-step approach
makes the inclusion/exclusion of a covariate a relatively simple task to
handle. Our proposal amends common practice in applied social science
research, where simple (low–level) LCA is done to obtain a classification of
low-level unit, and this is then related to (low– and high–level) covariates
simply including group fixed effects. Our analysis identifies latent classes
that score either consistently high or consistently low on all measured items,
along with two theoretically important classes that place distinctive emphasis
on items related to engaged citizenship, and duty-based norms.

###### keywords:

multilevel latent class analysis, stepwise estimators, citizenship norms

##  1 Introduction

Latent class analysis (LCA) is an approach used extensively in the social
sciences for creating a grouping based on a set of observed characteristics.
Research using this analytic approach has addressed a variety of substantive
topics, including for example Swiss population opinion about support of
pollution–reduction policies (Hill Kriesi, 2001), party performance in
volatile systems (Mustillo, 2009), electoral incentives on congressional
representation (Grimmer, 2013), and profiles of citizenship norms and
behaviors (Alvarez ., 2017; Hooghe ., 2016; Oser, 2017, 2022). More
methodological applications include modeling why support for extremist parties
is underestimated in surveys (Breen, 2000), and analysis of survey questions
that include a “don’t know” category (Feick, 1989). Examples from other fields
include, for instance, childhood socio–economic circumstances (Zhu ., 2020),
profiles of drinking behavior in US adolescents (Chung ., 2011; Reboussin
Ialongo, 2010), types of body mass index individuals (Brown ., 2020), and
clusters of scientific journals based on bibliometric indicators (Bartolucci,
Dardanoni Peracchi, 2015).

In LCA one of the key assumptions is that sample units are independent after
accounting for their class membership. The independence assumption is
inadequate when data with a hierarchical structure are encountered – as is
often the case in the social sciences. Examples are survey respondents nested
within regions/countries, children nested within families, and pupils within
schools. In such cases, sound data modeling should take the nested data
structure into account.

Multilevel extensions of LCA and mixture models, which are able to take the
nested data structure into account, are available (see, for instance,
Asparouhov Muthén, 2008; J. Vermunt, 2003, 2008). For example multilevel LCA
for nested categorical data, in which a natural grouping is observed, allows
researchers to also classify the groups (level 2 units) based on the
similarity of their members (level–1 units). For nested data such as pupils in
a school, a separate classification of both the pupils and the schools can be
achieved (e.g., Fagginger Auer ., 2016; Mutz Daniel, 2013). While multilevel
LCA is widespread, its use in political science and related fields is very
limited; a recent example can be found in Ruelens Nicaise (2020). See also
Morselli Glaeser (2018); Morselli Passini (2012).

In practice, stepwise estimators are popular in LCA because they separate the
estimation of the measurement and structural models. Substantive researchers
see the estimation of these two parts as distinct steps, which are often even
performed by different researchers at different time points (Bakk Kuha, 2018;
J. Vermunt, 2010). While in single level LCA the recent years have seen a
breakthrough in the development of bias–adjusted three–step estimators
(Asparouhov Muthèn, 2014; Bolck ., 2004; J. Vermunt, 2010), bias adjusted
stepwise estimators are less available in the multilevel settings (Bakk .,
2022). Furthermore, for models with distal outcomes, conceptually, adding the
distal outcome at the same time as the indicators is even counterintuitive:
the outcome defines the LC composition, while the latter should be predicted
by the latent variable (J. Vermunt, 2010). Having a readily available
bias–adjusted stepwise estimator for multilevel LCA is therefore crucial.

Both asymptotic as well as finite–sample properties of such class of stepwise
estimators have been widely studied, and are broadly understood. However,
bias–adjusted three–step estimators have one main drawback. Namely, they
require computing classification error probabilities in their second step,
which entails an explicit classification task. With multilevel models,
grouping occurs at two sequential layers. Therefore, an analogous step 2 in
this context poses substantial problems, whilst classification is often not
the goal of the analysis, but rather a juncture between estimation of
measurement and structural models.

We hereby propose to extend a recently developed two-step estimator for LCA to
the multilevel context (Bakk Kuha, 2018). Simulation results have shown that
two–step estimators (well–known in the econometric literature) have very good
finite sample properties. In addition, without an explicit classification
step, they are more robust in complex scenarios, e.g. under differential item
functioning (Di Mari, R. and Bakk, Z., 2018), or (mild) model
misspecifications (Di Mari ., 2020). In these and other challenging
circumstances, the design of the two–step estimator can be easily adapted (JK.
Vermunt Magidson, 2021).

We focus on the substantive application of multilevel LCA to the important
topic of the evolution of citizenship norms in a diverse array of democratic
countries. The rapid evolution of citizenship norms is considered one of the
driving forces behind political processes in contemporary democracies. An
important theoretical debate is whether more traditional, duty–based
citizenship norms are gradually being replaced by more self–expressive and
civically engaged values, variously referred to as “engaged” (Coffé Van Der
Lippe, 2010; Dalton, 2008), “assertive” (Dalton Welzel, 2014; Welzel ., 2005),
or “self–actualizing” (Bennett, 2008, 2012; Shehata ., 2016). Although
high–quality, cross-national data that is well–suited for latent class
analysis is publicly available to inform this debate, guidance is needed
regarding optimal multilevel modeling procedures to produce robust results.

Common practice in applied social science research is to conduct simple LCA on
the pooled data to obtain a classification of individual units. This is then
related to individual and country level covariates, and the hierarchical data
structure is taken into account by simply including country fixed effects in
the regression. From the stepwise LC modeling literature, we know this
approach to be wrong – in the sense that it produces downward biased
regression estimates. In addition, it can complicate estimation even with a
moderate number of countries, and does not allow for a clustering at the
country level – which is useful from both methodological and
application–oriented perspectives.

Model building in multilevel LCA can be a challenging task for both conceptual
and computational reasons. The first stage of this model building is the
selection of the number of classes at the lower and group level. Related to
this specific aspect, the standard approach is to jointly select the number of
classes at both individual and group levels. This might make the preliminary
stage of model selection extremely time consuming. Lukociene . (2010) suggest
a stepwise approach to class selection, where in step 1 the lower level LCA
model is defined in terms of number of classes. Then keeping the number of
lower level classes fixed, the group level classes are estimated (step 2a).
Next, a last step estimation at the lower level is re–iterated (step 2b) to
re–adjust the number of lower level classes if necessary based on the final
class solution at the higher level.

This model building strategy has the advantage of reducing the computational
burden of the simultaneous selection of the number of classes at both the
lower and group levels which works equally well in practice (see, for example,
J. Vermunt, 2008). Yet, Lukociene . (2010)’s multistep procedure requires
estimating the full multilevel LC model several times at steps 2a and 2b of
their procedure, which can still be time consuming if the choice of the number
of low/high level classes is not informed by some a priori knowledge of the
user, and/or is large. And even worse, with considerable sample size.

In fact, in what follows we argue that step 2a can be carried out in a more
efficient way by fixing the low–level parameters, and only estimating the
high–level ones. This has the advantage of enhancing the computational
efficiency of the whole selection procedure, and stepwise estimators of this
kind are well–known to enhance computation stability even if model complexity
increases (see, for instance, Bakk Kuha, 2018; Di Mari Maruotti, 2021;
Skrondal Kuha, 2012). Then generalized information criteria with well–known
statistical properties, commonly employed for model selection with pseudo-ML
estimators, can be used to select the high–level number of classes (Takeuchi,
1976; see also Bozdogan, 2000; Lv Liu, 2014). Next step 2b, where the full
measurement model is updated given the number of high level classes, becomes
useful to stabilize the measurement model before the inclusion of covariates.

Researchers gain important substantive insights from individual–level
clustering, as well as from multilevel models that identify distinctive
clusters in different group contexts (e.g., countries, schools). While
individual and group level clustering are interesting _per se_ , there is
often substantive interest in trying to predict individual and group clusters
based on covariates. In the multilevel LCA context, when the number of
candidate covariates is large the simultaneous estimator can be impractical:
each time a predictor is added or removed, the full model must be
re–estimated. In addition, algorithms for simultaneous estimation of models
for nested data – like multilevel and Markov models – can be numerically
unstable and slow to converge (Bartolucci, Montanari Pandolfi, 2015).

In the LCA literature it is widely acknowledged that including class
predictors in the stage of selection of the number of classes might alter
class enumeration and meaning (Asparouhov Muthén, 2014; J. Vermunt, 2010). The
general recommendation for individual–level LCA is to add covariates after the
number of classes is fixed (Nylund-Gibson ., 2019). Compared to selection of
the number of classes at the lower and group level, adding predictors of class
membership is seen both conceptually and practically a second stage of model
building. However, general recommendations for including predictors in
multilevel LC models are not available.

This paper introduces a new stepwise technique to fill these gaps: (i) by
amending Lukociene . (2010)’s approach in the phase of class selection.; (ii)
by specifying a two–stage estimator for the inclusion of predictors of class
membership, which extends the two–stage estimation approach for single level
LCA to the multilevel setting (Bakk Kuha, 2018). Operatively, the estimator is
articulated as follows. First, the lower level LC model is estimated for an
increasing number of classes, and the number of classes is chosen such that
some information criterion is optimized alongside substantive considerations
about the number of meaningful classes. Second, while both the number of low
level classes and the low level LC (lower level measurement model) parameters
are kept fixed, the group level number of classes is selected only estimating
the parameters of the equations for the group level measurement model. Third,
parameters of the lower level measurement model are re–estimated, keeping
fixed the group level model, in order to adjust for possible misspecifications
due to higher level grouping. Fourth, the structural model is built by loading
covariates on the (low and high level) latent variables, keeping the
measurement model’s parameters fixed.

Very recently, in a large simulation study, the robustness against violation
of the local independence assumption for the two–stage estimator for
multilevel LC models has been investigated (Bakk ., 2022). When all model
assumptions are met, the two–stage has similar properties to the one-stage
estimator. Where direct effects between the covariate and the items of the
latent class model are present, the two estimators have analogous
performances. Namely, as the severity of the underlying violations increases –
i.e. measurement noninvariance or differential item functioning (DIF) –
ignoring them leads to bias with both approaches. In related contexts, the
general two–step approach was found to be generally more robust to violations
of model assumptions in different parts of the joint model (Di Mari, R. and
Bakk, Z., 2018; Di Mari ., 2020). In the current work, we show how the
two–stage multilevel LCA, thereby specifically proposed for the context of
modeling DIF, can be used in more general settings as well.

This approach has the advantage of making both stages of multilevel LCA
conceptually sounder and practically simpler. We retain Lukociene . (2010)’s
idea of hierarchical class selection, which is known to work well, and we
substitute one of their full information ML steps with a consistent estimator.
Indeed, stepwise estimators of this type are well known to be consistent (Gong
Samaniego, 1981), and estimation of only a part of the model at each stage
greatly reduces the computational burden and enhances stability of the
algorithms that are commonly used to perform ML estimation (Di Mari, R. and
Bakk, Z., 2018; Di Mari ., 2016).

Regarding covariate selection, every time a covariate is added to or removed
from the model, class definitions and the optimal number of classes can
change. By separating the two parts – first stage and second stage of model
building – the measurement model is stabilized and is allowed to be determined
only by indicators. Moreover, this makes the inclusion/exclusion of a
covariate a relatively simpler task to handle, presumably even by different
researchers at different time points on (partially) not overlapping samples.

Our stepwise approach is motivated by the desire to analyze comparative data
on adolescents’ conceptions of good citizenship gathered by surveys conducted
by the International Association for the Evaluation of Educational Achievement
(IEA). Prior studies on citizenship norms have used latent class analysis with
country as a covariate to investigate the evolution of citizenship norms
between 1999 and 2009 (Hooghe Oser, 2015; Oser Hooghe, 2013), and the
determinants of good citizenship for all available countries in the 2009
sample (Hooghe ., 2016). Our analysis in this paper focuses on applying a
multilevel LCA approach to all 24 countries surveyed in 2016, which includes
advanced Western European and Scandinavian democracies often included in
high–quality surveys (e.g., Belgium, Denmark, Finland); a number of South
American countries (e.g., Chile, Colombia, Mexico) and Asian countries (e.g.,
Chinese Taipei, Hong Kong, Korea); as well as a few Eastern European and
lesser studied countries in these sorts of cross–national surveys (e.g.,
Russia, Malta).

The diversity of this high–quality cross–national data is a tremendous
resource that also requires careful modeling to distinguish between predictors
at the individual–level versus the country–level in order to obtain robust
findings for content experts on these topics. Using our stepwise estimator, we
analyze the IEA data, and provide substantive interpretation of the findings
as well as a step–by–step guide for software implementation.

The paper proceeds as follows. In Section 2 we present the modeling framework;
the two–stage estimator is introduced and described in detail in Section 3,
and the IEA data application is presented in Section 4. We conclude with final
remarks regarding the paper’s contributions and topics for future research in
Section 5.

##  2 The multilevel latent class model

###  2.1 Definition of the model

Consider the vector of responses $\mathbf{Y}_{ij}=(Y_{ij1},\dots,Y_{ijK})$,
where $Y_{ijk}$ denotes the response of individual $i$ in group $j$ on the
$k$–th categorical indicator variable, with $1\leq k\leq K$ and $1\leq j\leq
J$, where $K$ denotes the number of categorical indicators and $J$ the number
of level 2 units. In addition, we let $n_{j}$ denote the number of level 1
units within the $j$–th level 2 unit, with $1\leq j\leq J$. For simplicity of
exposition, we focus on dichotomous indicators.

LC analysis assumes that respondents belong to one of the $T$ categories
(“latent classes”) of an underlying categorical latent variable $X$ which
affects the responses (Goodman, 1974; Hagenaars, 1990; McCutcheon, 1987). The
model for $\mathbf{Y}_{ij}$ can then be specified as

$P(\mathbf{Y}_{ij})=\sum\limits_{t=1}^{T}P(X_{ij}=t)P(\mathbf{Y}_{ij}|X_{ij}=t),$
(1)

where the weight $P(X_{ij}=t)$ is the probability of person $i$ in group $j$
to belong to latent class $t$. The term $P(\mathbf{Y}_{ij}|X=t)$ is the
class–specific probability of observing a pattern of responses given that a
person belongs to class $t$. Furthermore we make the “local independence”
assumption that the $K$ indicator variables are independent within the latent
classes, leading to

$P(\mathbf{Y}_{ij})=\sum_{t=1}^{T}P(X_{ij}=t)\prod\limits_{k=1}^{K}P(Y_{ijk}|X_{ij}=t).$
(2)

Note that the general definition in Equation (1) applies to both the standard
and multilevel LC model. To be able to distinguish the simple and multilevel
LC model we must define the model in terms of logit equations. In the simple
LC model

$P(X_{ij}=t)=\frac{\exp({\gamma_{t}})}{1+\sum_{s=2}^{T}\exp(\gamma_{s})},$ (3)

for $1<t\leq T$ – where we have taken the first class as reference – and

$P(Y_{ijk}|X_{ij}=t)=\frac{\exp({\beta_{t}^{k}})}{1+\exp({\beta_{t}^{k}})},$
(4)

for $k=1,\dots,K$. In the simple LC model the parameters $\gamma$ and $\beta$
do not have the subscript $j$, thus assuming the clustering is independent of
the group level classes.

Extending the simple LC model to account for the multilevel data structure is
possible by allowing the parametrizations (3) and (4) to take the grouping
into account as follows

$P(X_{ij}=t)=\frac{\exp({\gamma_{tj}})}{1+\sum_{s=2}^{T}\exp({\gamma_{sj}})}$
(5)

again for $1<t\leq T$ and

$P(Y_{ijk}|X_{ij}=t)=\frac{\exp({\beta_{tj}^{k}})}{1+\exp({\beta_{tj}^{k}})}.$
(6)

for $k=1,\dots,K$.

This is the most general formulation that is equivalent to an unrestricted
multi–group LC model. In most applications however a more restricted version
is used (Lukociene ., 2010; J. Vermunt, 2003) that assumes that
item–conditional probabilities do not depend on the level 2 unit – as in
Equation (4). This is the restriction we will apply also in the data analysis
in the current paper.

###  2.2 Parametrizations of group–specific effects

The multilevel LC model can be seen as a random coefficients logistic
regression model (see, for instance Agresti ., 2000), where the dependent
variable is unobserved and has several observed indicators (J. Vermunt, 2003).
An example of such a model can be a LC model of lower level units of students,
predicted by several indicators of socio–economic status (e.g., parents
education, income level) that are nested in higher level units of schools -–
where clustering is identified at the school level that can be predicted by
additional higher–level variables (e.g., average value per square meter of
apartments in the area). As with random coefficient logistic regression
models, different ways of parametrizing multilevel LC models exist in the
literature (see also Finch French, 2014), mainly falling into the _parametric_
and _nonparametric_ approaches.

In the parametric approach, group–specific effects are assumed to arise from a
certain continuous distribution – typically Gaussian. With $T$ latent classes,
at least two alternatives are available. The first one requires to specify a
($T-1$) multidimensional distribution, say multivariate Gaussian. The second
alternative, as suggested by J. Vermunt (2003), is to work with the following
set of equations

$\gamma_{tj}=\gamma_{t}+\tau_{t}u_{j},$ (7)

with $u_{j}\sim N(0,1)$, for $j=1,\dots,J$ and $t=2,\dots,T$. As before, the
first class is taken as reference and the related parameters ($\gamma_{1}$ and
$\tau_{1}$) are set to zero. Intuitively, Equation (7) specifies random
components sharing the same random effect $u_{j}$ across $t$’s, that is then
scaled within each $t$ by the unknown $\tau_{t}$ (see also J. Vermunt, 2008).

In the nonparametric approach, the continuous distribution is replaced with a
discrete one (Aitkin, 1999; Laird, 1978), that is naturally assumed to be
multinomial with $M$ mass points each with probability $P(W=m)=\pi_{m}$. Thus,
we let $W$ be our multinomial group latent variable; by letting $W_{j}$ be the
value of $W$ for group $j$, in the nonparametric approach the model for the
(individual) latent class probabilities is specified as follows

$P(X_{ij}=t|W_{j}=m)=\frac{\exp({\gamma_{tm}})}{1+\sum_{s=2}^{T}\exp(\gamma_{sm})},$
(8)

where as before the first individual latent class is taken as reference and
its parameters are set to zero for identification. Also the mixing
probabilities $P(W=m)$ can be parametrized by means of logistic regressions as
follows

$P(W=m)=\frac{\exp(\delta_{0m})}{1+\sum_{l=2}^{M}\exp(\delta_{0l})},$ (9)

where parameters for $m=1$ are set to zero for identification and the related
class is set as reference. The path diagram of the resulting multilevel latent
class model is given in Figure 1. Note that $W_{j}$ affects the responses only
through the low level LC variable – and not directly – as in the restricted
multi–group LC model.

$W_{j}$$X_{1j}$$X_{2j}$$X_{n_{j}j}$$\dots$$Y_{1j1}$$Y_{1j2}$$\dots$$Y_{1jK}$$Y_{2j1}$$Y_{2j2}$$\dots$$Y_{2jK}$$Y_{n_{j}j1}$$Y_{n_{j}j2}$$\dots$$Y_{n_{j}jK}$
Figure 1: Multilevel LC model.

Such an approach, which this paper focuses on, allows the user 1) to work with
weaker assumptions, and 2) to reduce the computational burden (J. Vermunt Van
Dijk, 2001) with respect to the parametric approach. In addition, the
nonparametric approach has the byproduct of naturally providing a
classification of the units into a seemingly small number of classes (Aitkin,
1999). Under the nonparametric approach, the total number of free parameters
to be estimated is $M-1+M(T-1)+KT$.

###  2.3 Extending the model with covariates

Level 1 and level 2 covariates can be included to predict individual as well
as group class membership. Denoting one level 2 covariate by $Z_{1j}$ and a
level 1 covariate by $Z_{2ij}$, the model for
$\mathbf{Y}_{ij}|\mathbf{Z}_{j}$, with
$\mathbf{Z}_{j}=(Z_{1j},Z_{2ij})^{\prime}$, can be specified as follows

$P(\mathbf{Y}_{ij}|\mathbf{Z}_{j})=\sum_{m=1}^{M}P(W_{j}=m|Z_{1j})\sum_{t=1}^{T}P(X_{ij}=t|W_{j}=m,Z_{1j},Z_{2ij})\prod_{k=1}^{K}P(Y_{ijk}|X_{ij}=t),$
(10)

where the multinomial logistic regression for $X_{ij}$ with a random intercept
now is as follows

$P(X_{ij}=t|Z_{1j},Z_{2ij},W_{j}=m)=\frac{\exp({\gamma_{0tm}}+\gamma_{1t}Z_{1j}+\gamma_{2t}Z_{2ij})}{\sum_{s=2}^{T}\exp({\gamma_{0sm}}+\gamma_{1s}Z_{1j}+\gamma_{2s}Z_{2ij})}$
(11)

A random slope for the level 1 covariate can be obtained by replacing
$\gamma_{2t}$ by $\gamma_{2mt}$, for $1<t\leq T$. The multinomial logistic
regressions for class membership of group $j$ $P(W_{j}=m)$ is specified now as

$P(W_{j}=m|Z_{1j})=\frac{\exp(\delta_{0m}+\delta_{1m}Z_{1j})}{1+\sum_{l=2}^{L}\exp(\delta_{0l}+\delta_{1l}Z_{1j})}.$
(12)

Assuming a sample of $J$ groups is observed – each with $n_{j}$ individual
units, for $j=1,\dots,J$ – the (log) likelihood function of can be specified
as follows

$\log L(\mathbf{\theta})=\sum_{j=1}^{J}\log P(\mathbf{Y}_{j}|\mathbf{Z}_{j}),$
(13)

where $\mathbf{\theta}$ denotes the vector of all model parameters. The
function (13) is maximized with respect to the unknown $\mathbf{\theta}$ to
find ML estimates. This is typically accomplished by means of iterative
procedures – like the EM or Newton–type algorithms, or a combination of both.

However in practice, the choice of the number of latent classes on level 1 and
level 2 is not so obvious. A generally used recommendation is to use a
stepwise procedure for model selection (Lukociene ., 2010), by first fitting
the LC model at the level 1 – as defined in Equations 3 and 4. Once the number
of classes at the lower level is selected, this number is held fixed and the
number of classes is estimated at the group level. A general recommendation is
that, once the group level classes are selected, these are kept fixed, and
model selection is re–iterated at the lower level one more time before adding
covariates. In the stage of adding covariates the number of classes should be
fixed, also to be in line with general recommendations for LCA with covariates
(Masyn, 2017).

Using the procedure recommended by Lukociene . (2010), the model is
re–estimated to select the number of group level classes (J. Vermunt, 2003).
Also in this case, to add the covariates, while the number of classes is fixed
the full model will be re–estimated again. Given the complexity of such
multilevel models, 1) estimating the full model multiple times can be time
consuming, and 2) misspecifications in a part of the model may also
destabilize parameters in other parts of the model (J. Vermunt, 2010).

##  3 Stepwise estimation of the multilevel LC model parameters

We propose an efficient stepwise estimator within a two–stage model building
setup. Hereafter we outline the main idea of our two–stage approach, and give
details of each step afterwards.

  1. A)

First stage: unconditional multilevel LC model building (measurement model
construction).

    1. Step 1:

a simple latent class model is fitted on the pooled observations – multilevel
structure of the data is not taken into account (Figure 2a).

    2. Step 2.a:

the group level measurement model parameters of the multilevel LC model are
estimated keeping the lower level measurement model parameters fixed at their
first step values (Figure 2b).

    3. Step 2.b:

the full measurement model is re–estimated to account for possible interaction
effects (Figure LABEL:fig:step3) between low and high level measurement parts
This step is only necessary if the number of low and high level classes is not
known a priori.

  2. B)

Second stage: inclusion of covariates in the model (structural model
construction).

    1. Step 3:

covariates are included in the model (Figure LABEL:fig:step4) and only the
structural model parameters are estimated keeping the measurement model fixed.

$X$$Y_{111}$$Y_{112}$$\dots$$Y_{11K}$$Y_{n_{j}j1}$$Y_{n_{j}j2}$$\dots$$Y_{n_{j}jK}$$Y_{NJ1}$$Y_{NJ2}$$\dots$$Y_{NJK}$
(a) Step 1: simple latent class model on the pooled observations – multilevel
structure of the data is not taken into account.

$W_{j}$$X_{1j}$$X_{2j}$$X_{n_{j}j}$$\dots$$Y_{1j1}$$Y_{1j2}$$\dots$$Y_{1jK}$$Y_{2j1}$$Y_{2j2}$$\dots$$Y_{2jK}$$Y_{n_{j}j1}$$Y_{n_{j}j2}$$\dots$$Y_{n_{j}jK}$fixedfixedfixedfixedfixedfixedfixedfixedfixed
(b) Step 2.a: multilevel latent class model with measurement model kept fixed
at Step 1’s values.

Figure 2: Steps 1 and 2.a of the proposed estimator to fit a multilevel latent
class model. Note that Step 1 is equivalent to simple LCA on the pooled
observations.

$W_{j}$$X_{1j}$$X_{2j}$$X_{n_{j}j}$$\dots$$Y_{1j1}$$Y_{1j2}$$\dots$$Y_{1jK}$$Y_{2j1}$$Y_{2j2}$$\dots$$Y_{2jK}$$Y_{n_{j}j1}$$Y_{n_{j}j2}$$\dots$$Y_{n_{j}jK}$
(a)

$W_{j}$$Z_{1j}$$X_{1j}$$X_{2j}$$X_{n_{j}j}$$\dots$$Z_{2ij}$$Y_{1j1}$$Y_{1j2}$$\dots$$Y_{1jK}$$Y_{2j1}$$Y_{2j2}$$\dots$$Y_{2jK}$$Y_{n_{j}j1}$$Y_{n_{j}j2}$$\dots$$Y_{n_{j}jK}$fixedfixedfixedfixedfixedfixedfixedfixedfixed
(b)

Figure 3: Steps 2.b and 3 of the proposed estimator to fit a multilevel latent
class model.

###  3.1 Stage 1. Measurement model building

####  3.1.1 Step 1: Simple LC model

In this step (Figure 2a) a simple LC model is estimated on the pooled data
$T_{\text{max}}$ times, where $T_{\text{max}}$ is pre–specified by the user.
We let
$\mathbf{\theta}_{1}^{(T)}=(\beta_{2}^{1},\dots,\beta_{T}^{1},\dots,\beta_{2}^{K},\dots,\beta_{T}^{K})^{\prime}$
for each choice of $T=1,\dots,T_{\max}$. Under the parametrizations (3) and
(4), and a sample of $N$ observations – where $N=\sum_{j=1}^{J}n_{j}$ – the
(log) likelihood function of the first step model can be specified as follows

$\log L(\mathbf{\theta}_{1})=\sum_{i=1}^{N}\log P(\mathbf{Y}_{ij}),$ (14)

which we maximize in order to find the ML estimates of the LC model
parameters, that we call
$\widehat{\mathbf{\theta}}_{1}^{(T)}=(\widehat{\beta}_{2}^{1},\dots,\widehat{\beta}_{T}^{1},\dots,\widehat{\beta}_{2}^{K},\dots,\widehat{\beta}_{T}^{K})^{\prime}$.
Then, the optimal number of classes $T^{*}$ is selected such that
$T^{*}=\min\limits_{T=1,\dots,T_{\text{max}}}I(T)$ – where $I(T)$ is some
information criterion like AIC or BIC – along with the corresponding ML
estimates $\widehat{\boldsymbol{\theta}}_{1}$ from which we have suppressed
the superscript (T) for simplicity of notation.

####  3.1.2 Step 2.a: Multilevel LC model

In this step (Figure 2b), the group level measurement model parameters of the
multilevel LC model are estimated keeping measurement model parameters at the
lower level fixed at $\widehat{\boldsymbol{\theta}}_{1}$. Similarly as for
Step 1, this step has to be carried out $M_{\text{max}}$ times, where
$M_{\text{max}}$ is a pre–specified by the user number of latent classes for
$W$. We let
$\mathbf{\boldsymbol{\theta}}_{2}^{(M)}=(\delta_{2},\dots,\delta_{M},\dots,\gamma_{021},\dots,\gamma_{0T1},\dots,\gamma_{0TM})^{\prime}$
for each choice of $M=1,\dots,M_{\max}$. Under the parametrizations (4), (8)
and (9), and a sample of $J$ groups – each with $n_{j}$ individual units, for
$j=1,\dots,J$ – the (log) likelihood function of the step 2.a model can be
written as follows

$\log
L(\mathbf{\boldsymbol{\theta}}_{2}|\boldsymbol{\theta}_{1}=\widehat{\boldsymbol{\theta}}_{1}^{T^{*}})=\sum_{j=1}^{J}\log
P(\mathbf{Y}_{j}),$ (15)

where $|\boldsymbol{\theta}_{1}=\widehat{\boldsymbol{\theta}}_{1}$ indicates
that the measurement model parameters are kept fixed at
$\widehat{\boldsymbol{\theta}}_{1}$. The function (15) is maximized with
respect to the unknown $\mathbf{\boldsymbol{\theta}}_{2}$ to find ML
estimates. Note that only $M-1$ parameters are estimated in this step, out of
$M-1+M(T-1)+KT$ of the full multilevel LC model.

The fact that some of the model parameters are kept fixed in this pseudo–ML
step brings in a certain degree of misspecification. Therefore we propose
finding the optimal number of high–level groups $M$ by minimizing the Takeuchi
Information Criterion (TIC, Takeuchi, 1976). The TIC is defined as follows

$\text{TIC}=-2\log
L(\widehat{\boldsymbol{\theta}}_{2}|\boldsymbol{\theta}_{1}=\widehat{\boldsymbol{\theta}}_{1})+2tr(\widehat{\mathbf{H}}),$
(16)

where $\widehat{\mathbf{H}}=\widehat{\mathbf{F}}^{-1}\widehat{\mathbf{R}}$ is
an estimate of the contrast matrix and $\widehat{\mathbf{F}}^{-1}$ and
$\widehat{\mathbf{R}}$ are respectively the estimates of inverse Fisher
information matrix in inner productor (Hessian form) and of the outer product
form of the Fisher information matrix. Note that $\widehat{\mathbf{H}}$ is the
well–known Lagrange multiplier test statistic and $tr(\widehat{\mathbf{H}})$
incorporates the effects of model complexity and model misspecification in a
single term.

On a general note, the use of TIC is indeed very broad, and covers diverse
areas of statistics, including regression and generalized linear models (Lv
Liu, 2014), as well as latent variable models (Vrieze, 2012; see also Ranalli
Rocci, 2016).

####  3.1.3 Step 2.b: Re–update the measurement model

In this step (Figure LABEL:fig:step3) the measurement model is re–estimated
for a given number of high–level classes. We argue that such a step is
necessary especially if the number of low–level classes is unknown, to avoid
overestimating $T$ due to unmodeled interactions between low and high level
variability or poor lower–level class separation (Lukociene ., 2010). In those
cases, with $M$ selected from the previous step, this step is carried out
$T_{\text{max}}$ times.

Given a sample of $J$ groups – each with $n_{j}$ individual units, for
$j=1,\dots,J$, under the parametrizations (4), (8) and (9), the (log)
likelihood function of the step 2.b model can be written as follows

$\log L(\mathbf{\theta}_{1},\theta_{2})=\sum_{j=1}^{J}\log P(\mathbf{Y}_{j}).$
(17)

This specification, with respect to that of Step 1, now takes the multilevel
structure of the data into account.

###  3.2 Stage 2. The structural model with covariates

####  3.2.1 Step 3: including predictors for class memberships

As a next step the covariates can be added to the model (Figure
LABEL:fig:step4). At this stage also a decision needs to be taken whether a
stepwise approach is preferred (adding first lower level covariates, and after
fixing those adding at the group level) or all covariates can be added in a
single step. The benefit of the first option can be robustness, however no
simulation or theoretical results are available. For sake of conciseness, we
will present the simultaneous step 3 – its split counterpart can be derived
analogously.

Let us define
$\mathbf{\theta}_{3}=(\gamma_{12},\dots,\gamma_{1T},\gamma_{22},\dots,\gamma_{2T})^{\prime}$.
With the parametrizations specified in Equations (4), 11 and (12), the model
log–likelihood can be written as follows

$\log
L(\mathbf{\theta}_{2},\mathbf{\theta}_{3}|\theta_{1}=\widehat{\theta}_{1})=\sum_{j=1}^{J}\log
P(\mathbf{Y}_{j}|\mathbf{Z}_{j}),$ (18)

which we maximize with respect to $\mathbf{\theta}_{2}$ and
$\mathbf{\theta}_{3}$ keeping $\mathbf{\theta}_{1}$ fixed at its step 2.b
values $\widehat{\theta}_{1}$.

###  3.3 Statistical properties of the two–stage estimator

Our two–stage estimator is an instance of pseudo maximum likelihood estimation
(Gong Samaniego, 1981). Such estimators are consistent under very general
regularity conditions (see, for instance, Gourieroux Monfort, 1995). Using our
notation, let the true parameter vector be
$\operatorname{\boldsymbol{\theta}}^{\star}=(\operatorname{\boldsymbol{\theta}}_{1}^{\star},\operatorname{\boldsymbol{\theta}}_{2}^{\star},\operatorname{\boldsymbol{\theta}}_{3}^{\star})^{\prime}$.
Then our two–stage estimate $\widehat{\operatorname{\boldsymbol{\theta}}}$ is
consistent for $\operatorname{\boldsymbol{\theta}}^{\star}$ if the ML of
$\operatorname{\boldsymbol{\theta}}$ is itself consistent for
$\operatorname{\boldsymbol{\theta}}^{\star}$ and, most importantly, if 1)
$\operatorname{\boldsymbol{\theta}}_{1}$ and
$\operatorname{\boldsymbol{\theta}}_{2}$ and
$\operatorname{\boldsymbol{\theta}}_{3}$ can vary independently of each other,
and 2) $\widehat{\operatorname{\boldsymbol{\theta}}}_{1}$ is consistent for
$\operatorname{\boldsymbol{\theta}}_{1}^{\star}$. These conditions are
satisfied in our case: 1) is satisfied since
$\operatorname{\boldsymbol{\theta}}_{2}$ and
$\operatorname{\boldsymbol{\theta}}_{3}$ are independent from
$\operatorname{\boldsymbol{\theta}}_{1}$; 2) is satisfied as well, since the
parameter estimates at step 1 is a ML estimate of a simple LC model.

From standard theory we know that our two–stage estimator is less efficient
than the full information (one–step) ML estimator. In order to correctly
estimate the true variance of our estimator, we suggest to use non–parametric
bootstrap.

Regarding model selection, simultaneous model selection on a LC multilevel
model without covariates should always be considered as the statistical first
best. However, although it is the gold standard, the simultaneous selection
strategy becomes quickly unfeasible even for moderate $T$ and $M$. This is why
hierarchical model selection has been proposed in the literature by Lukociene
. (2010). All their steps are full ML steps, whereas we suggest a shortcut in
step 2a. In our step 2a only the high–level parameters are estimated, and the
low level parameters from the previous step are kept fixed. This is an
instance of pseudo ML estimation: such an approach has the advantage of
enhancing convergence speed, and also enhances computation stability (Bakk
Kuha, 2018), especially when a priori information on the number of high level
classes is not available. Then, available information criteria for
misspecified models (Takeuchi, 1976; see also Bozdogan, 2000; Lv Liu, 2014)
can be used instead of standard AIC or BIC to correctly assess the number of
high–level classes.

###  3.4 Model selection: an extended approach

In the previous subsections, we have advocated the use of information criteria
to select the number of classes. There can be real–data problems where one
specification undoubtedly dominates the others: in such cases, the final
choice is unquestionable. However, in certain occasions “there is no
particular reason to choose a single best model according to some criterion.
Rather it makes more sense to deselect models that are obviously poor,
maintaining a subset for further consideration. Sometimes this subset might
consist of a single model, but sometimes perhaps not” (Kadane Lazar, 2004).

Consequently, a wise use of information criteria should be complemented with
(i) the evaluation of class separation, and, perhaps most importantly, (ii)
the inspection of comparable solutions (in terms of information criteria and
class separation measures) to locate the most substantively meaningful model
configuration(s) (see, among others, Magidson Vermunt, 2004).

Several options are available to assess class separation in latent variable
models. An extensive review can be found in J. Vermunt Magidson (2016); see
also Farcomeni . (2021) for an example with nested data. The most extensively
used measure in the social sciences is Magidson (1981)’s entropy–based R2. In
the multilevel context, the latter can be meaningfully defined at both lower
and higher levels - see, e.g., (Lukociene ., 2010).

More specifically, let

o

$P(X_{ij}=t,W_{j}=m|\mathbf{Y}_{ij},Z_{1j},Z_{2ij})$ be the joint (bivariate)
posterior probability that individual $i$, within country $j$, belongs to
low–level class $t$, and high–level class $m$, given the observed responses,
and possibly covariates;

o

$P(W_{j}=m|\mathbf{Y}_{ij},Z_{1j})$ be the posterior probability that country
$j$ belongs to high–level class $m$, given the observed responses, and
possibly covariates;

o

$P(X_{ij}=t|\mathbf{Y}_{ij},Z_{1j},Z_{2ij})=\sum_{m=1}^{M}P(X_{ij}=t,W_{j}=m|\mathbf{Y}_{ij},Z_{1j},Z_{2ij})$
be the marginal posterior probability that individual $i$, within country $j$,
belongs to low–level class $t$, given the observed responses, and possibly
covariates.

The first two quantities are already available within the Expectation step of
the EM algorithm used for ML estimation of the second stage model (further
technical details can be found in J. Vermunt, 2003). The marginal posteriors
$P(X_{ij}=t|\mathbf{Y}_{ij},Z_{1j},Z_{2ij})$, as noted above, can be obtained
without effort from the joint posteriors, by summing over $W_{j}$.

The entropy–based R2 measures the proportional reduction of entropy when both
responses and covariates are available, compared to the case where these are
unknown. Its general formulation is as follows

$\text{R}^{2}_{\text{entr,low/high}}=\frac{E_{\text{TOT,low/high}}-\bar{E}_{\text{low/high}}}{E_{\text{TOT,low/high}}},$
(19)

where the subscript “low/high” indicates that the general definition can be
applied to both hidden layers, and “$E$” stands for entropy.

At lower level,
$E_{\text{TOT,low}}=\sum_{t=1}^{T}-\widehat{p}_{t}\log\widehat{p}_{t}$, where
$\widehat{p}_{t}=\frac{1}{N}\sum_{j=1}^{J}\sum_{i=1}^{n_{j}}P(X_{ij}=t|\mathbf{Y}_{ij},Z_{1j},Z_{2ij})$,
and $N$ is the total sample size;
$\bar{E}_{\text{low}}=\frac{1}{N}\sum_{j=1}^{J}\sum_{i=1}^{n_{j}}-P(X_{ij}=t|\mathbf{Y}_{ij},Z_{1j},Z_{2ij})\log\left\\{P(X_{ij}=t|\mathbf{Y}_{ij},Z_{1j},Z_{2ij})\right\\}.$

At higher level,
$E_{\text{TOT,high}}=\sum_{m=1}^{M}-\widehat{\omega}\log\widehat{\omega}$,
where $\widehat{\omega}=\frac{1}{J}P(W_{j}=m|\mathbf{Y}_{ij},Z_{1j})$, and
$\bar{E}_{\text{high}}=\frac{1}{J}-P(W_{j}=m|\mathbf{Y}_{ij},Z_{1j})\log\left\\{P(W_{j}=m|\mathbf{Y}_{ij},Z_{1j})\right\\}$.

In order to compare and/or rank similar solutions, any decision process, based
on information criteria, and on class separation measures, substantively
benefits from a comparative evaluation of the key model features. These
include characteristics of the measurement, e.g. the lower and higher level
different item–class profiles, as well as the distinct structural model
estimates that are obtained for different choices of the number of groups.
Such an extended approach is an important part of our empirical strategy, and,
more in general, what we recommend using in multilevel latent class analysis
with external predictors.

##  4 Multilevel LCA of cross–national citizenship norms

###  4.1 Data description

As part of a comprehensive evaluation of education systems, the IEA conducted
surveys in 1999, 2009 and 2016 in school classes of 14–year olds to
investigate civic education with the same scientific rigor as the evaluation
of more traditional educational skills of language and mathematics.

The target population of the survey is defined as all students in Grade 8
(approximately 14–year old), provided that the average age of students in this
grade was 13.5 years or above at the time of the assessment. If, on average,
students in Grade 8 were below 13.5 year old, Grade 9 became the target
population. The school samples were designed as stratified two–stage cluster
samples; schools were randomly selected at the first stage with probability
proportional to size, and intact classrooms were sampled at the second stage.

Typically, each country aimed for a sample size of 150 schools, the total
number of sampled students ranging between 3000 and 4500. The sample
participation requirement for the survey was 85 percent participation of the
selected schools and 85 percent of the selected students within the
participating schools, or a weighted overall participation rate of 75 percent.
Countries that did not meet the required response rate are excluded from the
sample. Additional details can be found in Schulz, Carstens . (2018).

More in general, the IEA’s survey design is comparable over time, and the
involvement of educational authorities in the participating countries’
school–based fieldwork ensures a high level of representativeness for a
respondent–based survey. The resulting data are therefore not compromised by
sampling challenges that researchers increasingly face in efforts to obtain
representative survey data of adults through common approaches of random
sampling through telephone, face–to–face interviews, or online panel samples.
Due to this internationally comparable and representative sampling design,
data from these surveys has been used to inform research on a variety of
topics related to civic attitudes over time, ranging from attitudes towards
migrants (Amnå Zetterberg, 2010; Munck ., 2018) to attitudes on citizenship
learning and education studies (Barber Torney-Purta, 2012; Knowles ., 2018),
and tolerance and civic engagement (Sandoval-Hernández ., 2018; Torney-Purta,
2002).

The survey includes a comprehensive battery of questions on citizenship norms
that can be validly compared across time and place: in 1999 in 28 countries
(Schulz ., 2004), in 2009 in 38 countries (Schulz ., 2011), and in 2016 in 24
countries (Schulz, Ainley ., 2018). The present study focuses on the third
wave of the survey that was conducted in 2016. The full sample comprises 94603
units, but we focus on complete cases only - 87123 respondents, as about 8% of
the full sample has missing values on at least one indicator. More in detail,
each country is measured from 1300 up to more than 7000 times by low–level
units over 12 indicators (see Figure 1). This, according to recent
recommendations (Park Yu, 2018), is a safe environment for multilevel LCA with
the nonparametric approach, even with low class distinctness and high
complexity of the latent structure – where this is defined by the number of
latent clusters and classes at the higher and lower levels.

Country | sample size
---|---
Belgium | 2750
Bulgaria | 2682
Chile | 4753
Colombia | 4992
Denmark | 5692
Germany | 1313
Dominican Republic | 2779
Estonia | 2770
Finland | 3037
Hong Kong | 2553
Croatia | 3655
Italy | 3274
Republic of Korea | 2557
Lithuania | 3422
Latvia | 3000
Mexico | 4987
Malta | 3317
Netherlands | 2692
Norway | 5740
Peru | 4713
Russia | 7049
Slovenia | 2664
Sweden | 2828
Taiwan | 3904
Table 1: Number of respondents per country of the third wave (2016) of the IEA
survey used for the analysis.

Questions regarding citizenship norms in all three waves asked respondents to
explain their understanding of what a good adult citizen is or does. The
survey then lists a variety of activities for respondents to rate in terms of
how important these activities are in order to be considered a good adult
citizen. The twelve items range from obeying the law and voting in elections,
to protecting the environment and defending human rights. The items are
presented in random order, so that respondents can indicate their own
individual set of specific preferences. The survey question regarding
citizenship norms is worded as follows: (Köhler ., 2018) Q23 How important are
the following behaviors for being a good adult citizens? Respondent options
include “very important”; “quite important”; “not very important”; and “not
important at all”. The options “very important” and “quite important” are
coded as 1, and “not very important” and “not important at all” are coded as
0.

Citizenship norm items included in the analysis are as follows: Obeying the
law (obey), taking part in activities promoting human rights (rights),
participating in activities to benefit people in the community (local),
working hard (work), taking part in activities to protect the environment
(envir), voting in every election (vote), learning about the country’s history
(history), showing respect for government representatives (respect), following
the political issues in the newspaper on the radio or on tv (news),
participating in a peaceful protest against a law believed to be unjust
(protest), engaging in political discussions (discuss), and joining a
political party (party). See appendix for full question wording for the
citizenship norm indicators.

Covariates included in the analysis are customary determinants of citizenship
norms from the literature at the individual–level of socio–economic measures
and country–level measure of gross domestic product (GDP) per capita. For the
individual–level covariates, prior studies have found gender and
socio–economic characteristics to be important covariates with citizenship
norms, with engaged adolescents more likely to be female and have greater
socio–economic advantages. We follow standard practice of research on
adolescents to operationalize socio–economic status through a proxy measure of
the number of books at home, along with measures of respondents’ educational
expectations, and parental education (Andrew Hauser, 2011; Hooghe ., 2016).
See the appendix for additional information on covariate question wording and
variable categories.

Data preparation and model fit have been carried out with Latent GOLD 5.1 (J.
Vermunt Magidson, 2016) in combination with R111Data preparation, as well as
estimation code and syntax files for replication are available from the
corresponding author upon request..

###  4.2 Results

In this section we first report on the measurement model resulting from the
first stage of the two–stage estimator. Then, we relate the low and high level
latent classes to covariates and – for sake of comparison – we present the
regression outputs for the simultaneous (one–stage) estimator and the
computationally more efficient two–stage estimator. Subsequently, we compare
both estimators for the multilevel latent class models to the results obtained
with the naïve estimator – which relates external variables to class
membership through a simple multinomial regression of the estimated vector of
class memberships on covariates, including country fixed effects.

We selected the number of classes based on the procedure described in Section
3. A preliminary assessment of $T$ (Step 1) was done by fitting a simple LC
model on the pooled data for $T=1,\dots,10$. As displayed by the model
summaries in Table 4, the minimum BIC model is achieved at $T=10$. This is
unfortunately only a ceiling effect, as by increasing $T$ up to 15 lowers the
BIC further. Balancing BIC values with a class separation measure (the
entropy–based $R^{2}$), a fairly good choice of $T$ is 5.

We thus proceed to Step 2a, fixing $T=5$ and selecting $M$ between 1 and 4 (we
do not go beyond 4 to avoid class uninterpretability and parameter
proliferation). Model summaries are displayed in Table 5. Together with TIC
values, we also report the BIC with penalty term computed based on the
high–level sample size following Lukociene . (2010)’s recommendation, and the
standard BIC based on low–level sample size. The best model according to TIC –
and all other criteria – is with $M=3$.

Table 6 reports the summary information for model fits of Step 2b. In this
case, the minimum BIC model is achieved at $T=10$, although the best solutions
balancing model fit with parsimony and class separation are with T ranging
between 4 and 6. Arguably, in terms of class interpretation, $T=4$ offers the
most substantively meaningful model configuration (see Figure 4). This is also
consistent with previous research on prior waves of the same survey conducted
in 1999 and 2009 (Hooghe ., 2016; Hooghe Oser, 2015; Oser Hooghe, 2013). The 2
new classes identified in the 5- and 6-class models offer more nuanced
variations on the 4-class solution but are less theoretically coherent in
relation to the literature, and tend more to mean scores on all measures.

For the 5-class solution (see 5), based on class sizes, the new class seems to
emerge almost completely from taking the “duty” class in the 4-class solution
that was 25.58% of the population, and splitting it in the 5-class solution
into LC2 (duty, 16%), and the new LC5, 12%. This can be identified as somewhat
between engaged and duty - we could also call it “ mainstream”, since its
values are around the mean prevalence of the indicators in the overall sample.
The 6-class solution (see 6) includes a new class - cluster 4 - that is akin
to the 5-class solution’s new “mainstream” group but with some elements of the
“maximalist” group; and the additional “mainstream” class (cluster 6) has
similar emphases, with particularly high scores on news and discuss.

Overall, we find that results from the simultaneous and two–stage estimators
are very similar, which is to be expected in general when model assumptions
are not violated. However we report that the simultaneous estimator is about
10 times slower than that of the two–stage estimator in terms of computing
time. In addition, in line with the cited literature on stepwise estimators,
we observed a general enhanced algorithmic stability and speed of convergence
of the two–stage estimator (in terms of number of algorithm iterations) with
respect to the simultaneous estimator.

The LC response probabilities – as after step 2b – are reported in Table 2
(right part)222Figure 4 gives a visual representation of the class profiles..
The findings point to four main distinctive citizenship norms. Two of these
latent classes are distinctively characterized by conditional probabilities
that are consistently high (the “maximal” group, 44%) or consistently low
(“subject” group, 4%). The remaining two latent classes are characterized by
distinct emphases of the type of indicators of citizenship they view as most
important for good citizenship, and relate to two theoretically important
concepts in the literature on citizenship norms. The “engaged” group (27%)
emphasizes the types of norms most associated with more recently prevalent
engaged attitudes (e.g., protecting the environment) while placing less
emphasis on more traditional citizenship activities (e.g., party membership).
Alternatively, the “duty” class (26%) places a high emphasis on the importance
of voting, discussing politics, and party activity, but expresses relatively
low interest in engaged norms of protecting human rights, and engaging in
activities to benefit the local community.

| HL class 1 | HL class 2 | HL class 3 | LL class 1 | LL class 2 | LL class 3 | LL class 4
---|---|---|---|---|---|---|---
| (25%) | (50%) | (25%) | Maximal (44%) | Engaged (27%) | Subject (4%) | Duty (26%)
obey | 0.89 | 0.94 | 0.91 | 0.98 | 0.91 | 0.47 | 0.89
rights | 0.75 | 0.91 | 0.92 | 0.99 | 0.97 | 0.15 | 0.68
local | 0.72 | 0.88 | 0.88 | 0.97 | 0.92 | 0.12 | 0.64
work | 0.80 | 0.87 | 0.84 | 0.93 | 0.83 | 0.46 | 0.78
envir | 0.77 | 0.92 | 0.94 | 0.99 | 1.00 | 0.27 | 0.69
vote | 0.80 | 0.89 | 0.84 | 0.98 | 0.80 | 0.29 | 0.79
history | 0.73 | 0.88 | 0.88 | 0.96 | 0.91 | 0.28 | 0.64
respect | 0.77 | 0.84 | 0.81 | 0.89 | 0.81 | 0.24 | 0.77
news | 0.72 | 0.83 | 0.69 | 0.96 | 0.58 | 0.15 | 0.72
protest | 0.53 | 0.74 | 0.68 | 0.89 | 0.65 | 0.12 | 0.41
discuss | 0.38 | 0.55 | 0.34 | 0.76 | 0.12 | 0.03 | 0.34
party | 0.30 | 0.44 | 0.30 | 0.60 | 0.16 | 0.04 | 0.25
Maximal (LL class 1) | 0.22 | 0.62 | 0.31 |  |  |  |
Engaged (LL class 2) | 0.14 | 0.17 | 0.56 |  |  |  |
Subject (LL class 3) | 0.06 | 0.03 | 0.04 |  |  |  |
Duty (LL class 4) | 0.58 | 0.18 | 0.09 |  |  |  |

Table 2: The lower and higher level class latent class model averaged across
all countries.

At the group class level, item probabilities are reported in Table 2 (left
part). The main distinction is between a group class with relatively low
scores on all indicators (HL class 1, which is characterized by a predominance
of more advanced democracies, e.g., Belgium, Denmark, Sweden); a group class
with higher scores on all indicators (HL class 2, which is characterized by a
predominance of least developed democracies, e.g., Croatia, Russia, Taiwan);
and a group with mid–range scores that is characterized by mid–range levels of
democratic development (HL class 3, Bulgaria, Chile, Colombia). These findings
are consistent with prior related research (e.g., Hooghe ., 2016). Detailed
country profiles are available in the appendix.

In Table 3 we report results of the structural model parameters for the
one–stage, two–stage and naïve estimators. Overall, we observe that parameter
values are similar for all estimators. Also inference on the parameter values
leads to the same conclusions in most cases, although the naïve estimator
seems to show the well-documented tendency to have an inflated type–I error
(Di Mari ., 2016; J. Vermunt, 2010); see, for instance, parameter and SE
estimates for Education goal and Non–native language level.

The regression findings are consistent with prior research on this topic. The
engaged and maximal groups are more female, while the duty and subject groups
are more male. When socio–economic distinctions are identified, the engaged
norm is generally characterized by higher educational and socio–economic
status measures than the other identified norms.

| Higher level structural model
---|---
| Class 2 |  |  | Class 3 |  |  |
| one–stage | two–stage |  |  | one–stage | two–stage |  |  |
log GDP | -7.33∗∗∗ | -7.33∗∗ |  |  | -4.67∗ | -4.67∗∗ |  |  |
| (3.34) | (3.31) |  |  | (2.73) | (2.36) |  |  |
| Lower level structural model
| Class 2 (Engaged) | Class 3 (Subject) | Class 4 (Duty)
| one–stage | two–stage | naïve | one–stage | two–stage | naïve | one–stage | two–stage | naïve
SES proxy | -0.26∗∗ | -0.25∗∗ | 0.05 | -0.25∗∗ | -0.24∗ | -0.42∗∗∗ | -0.12∗∗ | -0.12∗∗ | -0.15∗∗∗
| (0.09) | (0.09) | (0.03) | (0.10) | (0.12) | (0.09) | (0.04) | (0.04) | (0.04)
Female | 0.37∗∗∗ | 0.38∗∗∗ | 0.39∗∗∗ | -0.46∗∗ | -0.47∗∗ | -0.52∗∗∗ | -0.22 | -0.2∗ | -0.19∗∗∗
| (0.11) | (0.10) | (0.05) | (0.15) | (0.15) | (0.12) | (0.14) | (0.10) | (0.06)
Education goal | 0.14 | 0.16 | 0.12∗ | -0.65∗∗∗ | -0.64∗∗∗ | -0.53∗∗∗ | -0.35∗∗∗ | -0.35∗∗∗ | -0.34∗∗∗
| (0.13) | (0.12) | (0.06) | (0.05) | (0.08) | (0.10) | (0.05) | (0.08) | (0.05)
Mother education | 0.00 | 0.01 | 0.02 | 0.03 | 0.04 | 0.12 | 0.01 | 0.01 | 0.00
| (0.07) | (0.08) | (0.04) | (0.04) | (0.06) | (0.10) | (0.03) | (0.04) | (0.04)
Father education | -0.10 | -0.09 | 0.09∗ | -0.07 | -0.07 | -0.15 | -0.09∗∗ | -0.09∗ | -0.06
| (0.08) | (0.09) | (0.04) | (0.08) | (0.08) | (0.09) | (0.03) | (0.04) | (0.04)
Non-native language level | -0.19 | -0.22 | 0.01 | 0.14 | 0.15 | 0.63∗ | 0.26 | 0.27 | 0.20
| (0.17) | (0.18) | (0.11) | (0.35) | (0.25) | (0.25) | (0.25) | (0.30) | (0.10)

Table 3: Structural model parameter estimates for the simultaneous (one–stage)
estimator, the two–stage estimator and the naïve estimator, with 4 lower level
and 3 higher level classes. The number of books available at home is taken as
socio–economic status (ses) proxy. First lower and upper level classes are
taken as reference for the multinomial logistic regressions. Country log GDP
is included also as low–level control for the naïve approach (not reported).
Standard errors in parentheses (standard errors for the two–stage estimator
are obtained with non parametric bootstrap on the 4th step). ∗∗∗
$p$–value$<$0.01, ∗∗ $p$–value$<$0.05, ∗ $p$–value$<$0.1

##  5 Discussion

In this paper we have proposed a two–stage estimator for multilevel LCA to
model citizenship norms with the IEA data We have shown that the proposed
estimator greatly simplifies model construction compared to the simultaneous
estimator, although they both produce almost equal results in terms of
structural model parameter estimates. This is to be expected if model
assumptions are not violated. On a practical note, we report that the
simultaneous estimator was about 10 times slower than that of the two–stage
estimator in terms of computing time to produce parameter estimates More in
general, we observed a general enhanced algorithmic stability in terms of
convergence to stationary points with no boundary values, and speed of
convergence, in terms of number of algorithm iterations, of the two–stage
estimator with respect to the simultaneous estimator, in line with the cited
literature on stepwise estimators.

In our empirical analysis we have also included results from the approach
(naïve approach) that is current common practice in applied social science
research. We report that we observed non-convergence (the maximum number of
iterations is reached but the difference of two subsequent log-likelihood
values was never smaller than or equal to 10-8) and boundary parameter
estimates due to the large number of country fixed effects that are included
in the model. In addition, although the overall conclusions are the same, the
naïve estimator seems to show the well-documented tendency to have an inflated
type–I error (Di Mari ., 2016; J. Vermunt, 2010). Our recommendation is to use
the naïve approach only if the number of group–level units is small, and to
generally rely on multilevel modeling.

Consistent with prior research on citizenship norms, the findings identify
several distinctive normative types. In addition to the theoretically relevant
classes that emphasize engaged versus duty–based types of norms, the findings
also identify classes that consistently score either high or low on all
normative items. The socio–demographic covariates are consistent with prior
research, confirming that the engaged norm adherents among adolescents tends
to be more female, and characterized by higher socio–economic status.

Whereby the two–stage estimator has been used with the IEA data, more
generally this estimator is well suited for all data sets with a similar
nested structure – and possibly several candidate covariates. As a general
recommendation based on our findings, although simultaneous estimators are
theoretically to be preferred over stepwise estimators for their statistical
properties, the two–stage estimator for LC multilevel models is a more
efficient and stable alternative. An issue that is worth still further
investigation using multilevel LCA (both one and two–step) is model selection
and model selection criteria. While AIC and BIC are readily available in
commercial software for multilevel models or, like the TIC, can easily be
derived from available output, other more sophisticated model selection
approaches such as the bootstrap likelihood–ratio test Nylund . (2007) have
not yet been extended to these family of models.

Similar two–stage (or two–step) approaches have been already proposed in other
contexts, like simple LC or latent Markov modeling. A multi-stage
(hierarchical) class selection approach has been explored by Lukociene .
(2010) in multilevel models, but the two–stage estimator is new in this
context: multilevel data have their own specific features, which require
appropriate tools. Note that, although longitudinal data can be viewed as
nested data (individuals are high–level units, time points are low–level
units), the modeling approaches are conceptually and practically different.

In this paper we implemented one of the two main ways we described as
approaches for conducting step 2.b, namely we fix all parameters in $W$ and
$X\mid W$ multinomial logistic regressions, and both LV number of classes and
re–estimate the level 1 measurement parameters. The second approach we
describe is to conduct a “full” 2.b by fixing only parameters in the
multinomial logistic regression for $W$, and re–estimate $X\mid W$ and $Y\mid
X$ multinomial logistic regression parameters for, say,
$T=1,\dots,T_{\text{max}}$. Although the latter approach is more
computationally intensive, it could be useful in specific situations as
described by Lukociene . (2010): “when the lower–level classes are poorly
separated, information on the higher–level classes may help in finding the
correct number of lower–level classes, provided that the higher–level classes
are well separated themselves” (Lukociene ., 2010).

One of the main working assumption of LC model parameter estimators is local
independence of the indicators given the latent variable(s). We have assumed
that this assumption holds in the IEA application. Local fit measures are
available for multilevel latent class models (Nagelkerke ., 2016, 2017),
although the issue of detecting measurement non–invariance is not yet resolved
as these are not asymptotically chi squared. Future research could focus on
the two–stage estimator’s behavior and adaptations to cases where the local
independence assumption is violated (e.g., due to direct effects of covariates
on the indicators).

Another issue with two–step estimation for multilevel LC models is obtaining
SE estimates that take into account uncertainty about the fixed parameters.
Standard errors based on the Hessian of the last step are known to
underestimate the true underlying variability of stepwise estimators (Bakk .,
2014). To obtain correct standard errors without increasing much the
computational burden, we bootstrapped (non–parametrically) only stage 2 – also
in line with guidelines for stepwise LCA (Bakk ., 2016). In general, non
parametric bootstrap can be a “safe” option using available softwares.
However, quicker options can be thought of, based on sandwich formulas (White,
1982) and correcting for underestimation of the true variability with additive
corrections (Oberski Satorra, 2013). Investigating theoretically correct SE
estimation – possibly with closed–form formulas – is an important topic for
future work.

##  Appendix A Additional tables and figures

Figure 4: The basic latent class model (measurement model) at the lower level

##  Appendix B Model summaries

$T$ | log–likelihood | BIC | Npar | entr. $R^{2}$
---|---|---|---|---
1 | -501351.0962 | 1002838.693 | 12 | 1
2 | -457351.5799 | 914987.5366 | 25 | 0.7213
3 | -448212.3953 | 896857.0434 | 38 | 0.6561
4 | -444754.1647 | 890088.4583 | 51 | 0.6239
5 | -442942.6846 | 886613.3741 | 64 | 0.6039
6 | -441147.5251 | 883170.9312 | 77 | 0.6065
7 | -440338.5716 | 881700.9001 | 90 | 0.5798
8 | -439962.9764 | 881097.5857 | 103 | 0.5732
9 | -439698.4559 | 880716.4206 | 116 | 0.5486
10 | -439499.6612 | 880466.7072 | 129 | 0.5694
Table 4: Summary statistics for model selection at low level (Step 1). Selected $T$ in bold.  $M$ | log–likelihood | BIC | BIC ($n_{j}$) | TIC | Npar | entr. $R^{2}$ (low–lev)
---|---|---|---|---|---|---
1 | -1702481.198 | 3405007.867 | 3404975.108 | 3404962.415 | 4 | 1
2 | -1701985.349 | 3404073.007 | 3403999.300 | 3403969.852 | 9 | 0.792
3 | -1701799.232 | 3403757.612 | 3403642.957 | 3403597.627 | 14 | 0.814
4 | -1701985.579 | 3404187.145 | 3404031.541 | 3403971.185 | 19 | 0.833
Table 5: Summary statistics for model selection at high level (Step 2a), with $nj=5$ as selected from Step 1. Selected $M$ in bold.  $T$ | log–likelihood | BIC | Npar | entr. $R^{2}$
---|---|---|---|---
1 | -484857.7962 | 969874.7406 | 14 | 1.0000
2 | -446868.0265 | 894065.7169 | 29 | 0.6893
3 | -436544.3193 | 873588.8183 | 44 | 0.6508
4 | -432088.8293 | 864848.3541 | 59 | 0.6401
5 | -429566.1149 | 859973.4411 | 74 | 0.7136
6 | -427043.1994 | 855098.1259 | 89 | 0.7534
7 | -424477.0087 | 850136.2604 | 104 | 0.7415
8 | -422633.8249 | 846620.4086 | 119 | 0.7387
9 | -422010.3229 | 845543.9203 | 134 | 0.6991
10 | -421339.4679 | 844372.7262 | 149 | 0.6917
Table 6: Summary statistics for model selection at low level (Step 2b), $M=3$.
Selected $T$ in bold.

##  Appendix C Output from 5 and 6 low–level class solutions

Figure 5: The basic latent class model (measurement model) at the lower level, with $T=5$, and $M=3$. Figure 6: The basic latent class model (measurement model) at the lower level, with $T=6$, and $M=3$. | Class 2 | Class 3 | Class 4 | Class 5
---|---|---|---|---
| (Engaged) | (Mainstream) | (Duty) | (Subject)
SES proxy | -0.19∗∗∗ | -0.10 | -0.23∗∗∗ | -0.43∗∗∗
| (0.05) | (0.08) | (0.05) | (0.11)
Female | 0.40∗∗∗ | 0.30∗∗∗ | -0.24∗∗∗ | -0.43∗∗
| (0.10) | (0.05) | (0.06) | (0.14)
Education goal | 0.26 | 0.01 | -0.34∗∗∗ | -0.59∗∗∗
| (0.15) | (0.14) | (0.06) | (0.10)
Mother education | 0.05 | -0.04 | 0.01 | 0.08
| (0.04) | (0.08) | (0.04) | (0.10)
Father education | -0.2∗ | -0.03 | -0.13∗∗ | -0.05
| (0.08) | (0.05) | (0.04) | (0.10)
Non-native language level | -0.52∗∗ | 0.07 | 0.16 | -0.28
| (0.19) | (0.18) | (0.26) | (0.43)
Table 7: Lower–level structural model parameter estimates for the two–stage
estimator, with 5 lower level and 3 higher level classes. First class taken as
reference. Standard errors, obtained with non parametric bootstrap on the 4th
step, in parentheses. ∗∗∗ $p$–value$<$0.01, ∗∗ $p$–value$<$0.05, ∗
$p$–value$<$0.1

| Class 2 | Class 3 | Class 4 | Class 5 | Class 6
---|---|---|---|---|---
| (Mainstream-discussion) | (Mainstream-maximal) | (Duty) | (Engaged) | (Subject)
SES proxy | -0.07 | -0.16∗∗∗ | -0.2∗∗∗ | -0.33∗ | -0.43∗∗∗
| (0.08) | (0.03) | (0.04) | (0.13) | (0.1)
Female | 0.37∗∗∗ | 0.5∗∗∗ | -0.17∗ | 0.16 | -0.38∗∗
| (0.08) | (0.11) | (0.09) | (0.3) | (0.14)
Education goal | 0.02 | 0.17 | -0.37∗∗ | 0.15 | -0.58∗∗∗
| (0.15) | (0.23) | (0.11) | (0.31) | (0.14)
Mother education | -0.04 | 0.01 | 0.03 | 0.01 | 0.08
| (0.11) | (0.07) | (0.06) | (0.11) | (0.11)
Father education | -0.02 | -0.13 | -0.09∗∗∗ | -0.31 | -0.07
| (0.09) | (0.09) | (0.02) | (0.17) | (0.11)
Non-native language level | 0.15 | -0.59∗ | 0.25 | -0.32 | -0.28
| (0.23) | (0.26) | (0.25) | (0.42) | (0.37)

Table 8: Lower–level structural model parameter estimates for the two–stage
estimator, with 6 lower level and 3 higher level classes. First class taken as
reference. Standard errors, obtained with non parametric bootstrap on the 4th
step, in parentheses. ∗∗∗ $p$–value$<$0.01, ∗∗ $p$–value$<$0.05, ∗
$p$–value$<$0.1

##  Appendix D Survey questions

The survey question regarding citizenship norms is worded as follows: (Köhler
., 2018)

Q23 How important are the following behaviors for being a good adult citizen?
Respondent options: “very important”; “quite important”; “not very important”;
“not important at all”.

IEA variable | Indicator | Questionnaire text
---|---|---
IS3G23L | obey | Always obeying the law
IS3G23I | rights | Taking part in activities
|  | promoting human rights
IS3G23H | local | Participating in activities to
|  | benefit people in
|  | the local community
IS3G23K | work | Working hard
IS3G23J | envir | Taking part in activities to
|  | protect the environment
IS3G23A | vote | Voting in every
|  | national election
IS3G23C | history | Learning about the
|  | country’s history
IS3G23E | respect | Showing respect for
|  | government representatives
IS3G23D | news | Following political issues
|  | in the newspaper, on the
|  | radio, on TV, or
|  | on the Internet
IS3G23G | protest | Participating in peaceful
|  | protests against laws
|  | believed to be unjust
IS3G23F | discuss | Engaging in political discussions
IS3G23B | party | Joining a political party

## References

  * Agresti . (2000) agresti2000Agresti, A., Booth, J., Hobert, J. Caffo, B.  2000\.  Random-effects modeling of categorical response data Random-effects modeling of categorical response data.  Sociological Methodology30127–80.
  * Aitkin (1999) aitkin1999Aitkin, M.  1999\.  A general maximum likelihood analysis of variance components in generalized linear models A general maximum likelihood analysis of variance components in generalized linear models.  Biometrics551117–128.
  * Alvarez . (2017) alvarez2017Alvarez, R., Levin, I. Núñez, L.  2017\.  The four faces of political participation in Argentina: Using latent class analysis to study political behavior The four faces of political participation in Argentina: Using latent class analysis to study political behavior.  The Journal of Politics7941386–1402.
  * Amnå Zetterberg (2010) amnaa2010Amnå, E. Zetterberg, P.  2010\.  A political science perspective on socialization research: Young Nordic citizens in a comparative light A political science perspective on socialization research: Young nordic citizens in a comparative light.  L. Sherrod, J. Torney-Purta C. Flanagan (), Handbook of Youth’s Civic Engagement Handbook of youth’s civic engagement ( 69-100).  Wiley, New York.
  * Andrew Hauser (2011) andrew2011Andrew, M. Hauser, R.  2011\.  Adoption? Adaptation? Evaluating the formation of educational expectations Adoption? adaptation? evaluating the formation of educational expectations.  Social Forces902497–520.
  * Asparouhov Muthén (2008) asparouhov2008Asparouhov, T. Muthén, B.  2008\.  Multilevel Mixture Models Multilevel mixture models.  G. Hancock K. Samuelsen (), Advances in Latent Variable Mixture Models Advances in latent variable mixture models ( 27-75).  Charlotte, NC: Information Age Publishing.
  * Asparouhov Muthèn (2014) asparouhov14bAsparouhov, T. Muthèn, B.  2014\.  Auxiliary Variables in Mixture Modeling: Three-Step Approaches Using Mplus Auxiliary variables in mixture modeling: Three-step approaches using mplus.  Structural Equation Modeling: A Multidisciplinary Journal213329-341.
  * Asparouhov Muthén (2014) asparouhov2012auxiliaryAsparouhov, T. Muthén, B.  2014\.  Auxiliary variables in mixture modeling: Three-step approaches using Mplus Auxiliary variables in mixture modeling: Three-step approaches using Mplus.  Structural Equation Modeling: A Multidisciplinary Journal21.
  * Bakk . (2022) bakk2022Bakk, Z., Di Mari, R., Oser, J. Kuha, J.  2022\.  Two-stage multilevel latent class analysis with covariates in the presence of direct effects Two-stage multilevel latent class analysis with covariates in the presence of direct effects.  Structural Equation Modeling: A Multidisciplinary Journal292267–277.
  * Bakk Kuha (2018) bakk+kuha18Bakk, Z. Kuha, J.  2018\.  Two-Step Estimation of Models Between Latent Classes and External Variables Two-step estimation of models between latent classes and external variables.  Psychometrika83871–892.
  * Bakk . (2014) zsuzsa2Bakk, Z., Oberski, D. Vermunt, J.  2014\.  Relating latent class assignments to external variables: standard errors for correct inference Relating latent class assignments to external variables: standard errors for correct inference.  Political Analysis224520–540.
  * Bakk . (2016) bakk2016Bakk, Z., Oberski, D. Vermunt, J.  2016\.  Relating latent class membership to continuous distal outcomes: Improving the LTB approach and a modified three-step implementation Relating latent class membership to continuous distal outcomes: Improving the ltb approach and a modified three-step implementation.  Structural Equation Modeling: A Multidisciplinary Journal232278–289.
  * Barber Torney-Purta (2012) barber2012Barber, C. Torney-Purta, J.  2012\.  Comparing attitudes in the 1999 and 2009 IEA civic and citizenship education studies: Opportunities and limitations illustrated in five countries Comparing attitudes in the 1999 and 2009 IEA civic and citizenship education studies: Opportunities and limitations illustrated in five countries.  JSSE-Journal of Social Science Education11147-63.
  * Bartolucci, Dardanoni Peracchi (2015) bartolucci2015Bartolucci, F., Dardanoni, V. Peracchi, F.  2015\.  Ranking scientific journals via latent class models for polytomous item response data Ranking scientific journals via latent class models for polytomous item response data.  Journal of the Royal Statistical Society: Series A (Statistics in Society)17841025-1049.
  * Bartolucci, Montanari Pandolfi (2015) bartolucci2015_3Bartolucci, F., Montanari, GE. Pandolfi, S.  2015\.  Three-step estimation of latent Markov models with covariates Three-step estimation of latent markov models with covariates.  Computational Statistics & Data Analysis83287–301.
  * Bennett (2008) bennett2008Bennett, W.  2008\.  Civic Life Online: Learning How Digital Media Can Engage Youth Civic life online: Learning how digital media can engage youth.  The John D. and Catherine T. MacArthur Foundation Series on Digital Media and Learning. Cambridge, MA: MIT Press.
  * Bennett (2012) bennett2012Bennett, W.  2012\.  The personalization of politics: Political identity, social media, and changing patterns of participation The personalization of politics: Political identity, social media, and changing patterns of participation.  The Annals of the American Academy of Political and Social Science644120–39.
  * Bolck . (2004) Bolck:04Bolck, A., Croon, M. Hagenaars, J.  2004\.  Estimating Latent Structure Models with Categorical Variables: One-Step versus Three-Step Estimators Estimating latent structure models with categorical variables: One-step versus three-step estimators.  Political Analysis123-27.
  * Bozdogan (2000) bozdogan2000Bozdogan, H.  2000\.  Akaike’s information criterion and recent developments in information complexity Akaike’s information criterion and recent developments in information complexity.  Journal of Mathematical Psychology44162–91.
  * Breen (2000) breenBreen, R.  2000\.  Why Is Support for Extreme Parties Underestimated by Surveys? A Latent Class Analysis Why Is Support for Extreme Parties Underestimated by Surveys? A Latent Class Analysis.  British Journal of Political Science30375-82.
  * Brown . (2020) brown2020Brown, S., Greene, W. Harris, M.  2020\.  A novel approach to latent class modelling: identifying the various types of body mass index individuals A novel approach to latent class modelling: identifying the various types of body mass index individuals.  Journal of the Royal Statistical Society: Series A.
  * Chung . (2011) chung2011Chung, H., Anthony, J. Schafer, J.  2011\.  Latent class profile analysis: an application to stage sequential processes in early onset drinking behaviours Latent class profile analysis: an application to stage sequential processes in early onset drinking behaviours.  Journal of the Royal Statistical Society: Series A (Statistics in Society)1743689–712.
  * Coffé Van Der Lippe (2010) coffe2010Coffé, H. Van Der Lippe, T.  2010\.  Citizenship Norms in Eastern Europe Citizenship Norms in Eastern Europe.  Social Indicators Research963479–496.
  * Dalton (2008) dalton2008Dalton, R.  2008\.  Citizenship norms and the expansion of political participation Citizenship norms and the expansion of political participation.  Political Studies56176–98.
  * Dalton Welzel (2014) dalton2014Dalton, R. Welzel, C.  2014\.  The civic culture transformed: From allegiant to assertive citizens The civic culture transformed: From allegiant to assertive citizens.  Cambridge University Press.
  * Di Mari, R. and Bakk, Z. (2018) dimaribakkDi Mari, R. and Bakk, Z.  2018\.  Mostly Harmless Direct Effects: A Comparison of Different Latent Markov Modeling Approaches Mostly harmless direct effects: A comparison of different latent markov modeling approaches.  Structural Equation Modeling: A Multidisciplinary Journal253467-483.
  * Di Mari . (2020) dimaribakkpunzo2020Di Mari, R., Bakk, Z. Punzo, A.  2020\.  A random-covariate approach for distal outcome prediction with latent class analysis A random-covariate approach for distal outcome prediction with latent class analysis.  Structural Equation Modeling: A Multidisciplinary Journal273351–368.
  * Di Mari Maruotti (2021) dimari2021Di Mari, R. Maruotti, A.  2021\.  A two–step estimator for generalized linear models for longitudinal data with time-varying measurement error A two–step estimator for generalized linear models for longitudinal data with time-varying measurement error.  Advances in Data Analysis and Classification1–28.
  * Di Mari . (2016) dimari2016Di Mari, R., Oberski, D. Vermunt, J.  2016\.  Bias-adjusted three-step latent Markov modeling with covariates Bias-adjusted three-step latent markov modeling with covariates.  Structural Equation Modeling: A Multidisciplinary Journal235649–660.
  * Fagginger Auer . (2016) faggingerFagginger Auer, MF., Hickendorff, M., Van Putten, CM., Bèguin, AA. Heiser, WJ.  2016\.  Multilevel latent class analysis for large-scale educational assessment data: Exploring the relation between the curriculum and Students’ mathematical strategies. Multilevel latent class analysis for large-scale educational assessment data: Exploring the relation between the curriculum and students’ mathematical strategies.  Applied Measurement in Education29144-159.
  * Farcomeni . (2021) farcomeni2021Farcomeni, A., Ranalli, M. Viviani, S.  2021\.  Dimension reduction for longitudinal multivariate data by optimizing class separation of projected latent Markov models Dimension reduction for longitudinal multivariate data by optimizing class separation of projected latent markov models.  Test302462–480.
  * Feick (1989) feickFeick, L.  1989\.  Latent class analysis of survey questions that include don’t know responses Latent class analysis of survey questions that include don’t know responses.  Public Opinion Quarterly53525-47.
  * Finch French (2014) finch2014Finch, W. French, B.  2014\.  Multilevel latent class analysis: Parametric and nonparametric models Multilevel latent class analysis: Parametric and nonparametric models.  The Journal of Experimental Education823307–333.
  * Gong Samaniego (1981) gong1981pseudoGong, G. Samaniego, F.  1981\.  Pseudo maximum likelihood estimation: Theory and applications Pseudo maximum likelihood estimation: Theory and applications.  The Annals of Statistics861–869.
  * Goodman (1974) goodman:74Goodman, L.  1974\.  The Analysis of Systems of Qualitative Variables When Some of the Variables Are Unobservable. Part I: A Modified Latent Structure Approach The analysis of systems of qualitative variables when some of the variables are unobservable. Part I: A modified latent structure approach.  American Journal of Sociology79-259.
  * Gourieroux Monfort (1995) gourieroux1995Gourieroux, C. Monfort, A.  1995\.  Statistics and Econometric Models Statistics and econometric models ( 1).  Cambridge University Press.
  * Grimmer (2013) grimmer2013appropriatorsGrimmer, J.  2013\.  Appropriators not position takers: The distorting effects of electoral incentives on congressional representation Appropriators not position takers: The distorting effects of electoral incentives on congressional representation.  American Journal of Political Science.
  * Hagenaars (1990) hagenaars:90Hagenaars, J.  1990\.  Categorical longitudinal data- Loglinear analysis of panel, trend and cohort data Categorical longitudinal data- loglinear analysis of panel, trend and cohort data.  Newbury Park, CA:Sage.
  * Hill Kriesi (2001) hill2001classificationHill, J. Kriesi, H.  2001\.  Classification by opinion-changing behavior: A mixture model approach Classification by opinion-changing behavior: A mixture model approach.  Political Analysis94301–324.
  * Hooghe Oser (2015) hooghe2015Hooghe, M. Oser, J.  2015\.  The rise of engaged citizenship: The evolution of citizenship norms among adolescents in 21 countries between 1999 and 2009 The rise of engaged citizenship: The evolution of citizenship norms among adolescents in 21 countries between 1999 and 2009.  International Journal of Comparative Sociology56129–52.
  * Hooghe . (2016) hooghe2016Hooghe, M., Oser, J. Marien, S.  2016\.  A comparative analysis of ‘good citizenship’: A latent class analysis of adolescents’ citizenship norms in 38 countries A comparative analysis of ‘good citizenship’: A latent class analysis of adolescents’ citizenship norms in 38 countries.  International Political Science Review371115–129.
  * Kadane Lazar (2004) kadane2004Kadane, JB. Lazar, NA.  2004\.  Methods and criteria for model selection Methods and criteria for model selection.  Journal of the American Statistical Association99465279–290.
  * Knowles . (2018) knowles2018Knowles, R., Torney-Purta, J. Barber, C.  2018\.  Enhancing citizenship learning with international comparative research: Analyses of IEA civic education datasets Enhancing citizenship learning with international comparative research: Analyses of iea civic education datasets.  Citizenship Teaching & Learning1317-30.
  * Köhler . (2018) kohler2018Köhler, H., Weber, S., Brese, F., Schulz, W. Carstens, R.  2018\.  ICCS 2016 user guide for the international database: IEA International Civic and Citizenship Education Study 2016 ICCS 2016 user guide for the international database: IEA International Civic and Citizenship Education Study 2016.  Amsterdam: The International Association for the Evaluation of Educational Achievement (IEA).
  * Laird (1978) laird1978Laird, N.  1978\.  Nonparametric maximum likelihood estimation of a mixing distribution Nonparametric maximum likelihood estimation of a mixing distribution.  Journal of the American Statistical Association73364805–811.
  * Lukociene . (2010) Lukociene.10Lukociene, O., Varriale, R. Vermunt, J.  2010\.  The simultaneous decision(s) about the number of lower- and higher-level classes in multilevel latent class analysis The simultaneous decision(s) about the number of lower- and higher-level classes in multilevel latent class analysis.  Sociological Methodology401247–283.
  * Lv Liu (2014) lv2014Lv, J. Liu, JS.  2014\.  Model selection principles in misspecified models Model selection principles in misspecified models.  Journal of the Royal Statistical Society: Series B: Statistical Methodology141–167.
  * Magidson (1981) magidson1981Magidson, J.  1981\.  Qualitative variance, entropy, and correlation ratios for nominal dependent variables Qualitative variance, entropy, and correlation ratios for nominal dependent variables.  Social Science Research102177–194.
  * Magidson Vermunt (2004) magidsonvermunt2004Magidson, J. Vermunt, J.  2004\.  Latent class models Latent class models.  D. Kaplan (), The Sage handbook of quantitative methodology for the social sciences The sage handbook of quantitative methodology for the social sciences ( 175–198).  Sage.
  * Masyn (2017) masyn:difMasyn, K.  2017\.  Measurement Invariance and Differential Item Functioning in Latent Class Analysis With Stepwise Multiple Indicator Multiple Cause Modeling Measurement invariance and differential item functioning in latent class analysis with stepwise multiple indicator multiple cause modeling.  Structural Equation Modeling: A Multidisciplinary Journal242180-197.
  * McCutcheon (1987) mccutcheon:87McCutcheon, A.  1987\.  Latent Class Analysis Latent class analysis.  Newbury Park, CA:Sage.
  * Morselli Glaeser (2018) morselli2018Morselli, D. Glaeser, S.  2018\.  Economic conditions and social trust climates in Europe over ten years: An ecological analysis of change Economic conditions and social trust climates in Europe over ten years: An ecological analysis of change.  Journal of Trust Research8168–86.
  * Morselli Passini (2012) morselli2012Morselli, D. Passini, S.  2012\.  Disobedience and support for democracy: Evidences from the World Values Survey Disobedience and support for democracy: Evidences from the world values survey.  The Social Science Journal493284–294.
  * Munck . (2018) munck2018Munck, I., Barber, C. Torney-Purta, J.  2018\.  Measurement invariance in comparing attitudes toward immigrants among youth across Europe in 1999 and 2009: The alignment method applied to IEA CIVED and ICCS Measurement invariance in comparing attitudes toward immigrants among youth across europe in 1999 and 2009: The alignment method applied to iea cived and iccs.  Sociological Methods & Research474687–728.
  * Mustillo (2009) mustillo2009modelingMustillo, T.  2009\.  Modeling new party performance: A conceptual and methodological approach for volatile party systems Modeling new party performance: A conceptual and methodological approach for volatile party systems.  Political Analysis173311–332.
  * Mutz Daniel (2013) mutz2013Mutz, R. Daniel, H.  2013\.  University and student segmentation: Multilevel latent-class analysis of students? Attitudes towards research methods and statistics University and student segmentation: Multilevel latent-class analysis of students? Attitudes towards research methods and statistics.  British Journal of Educational Psychology832280–304.
  * Nagelkerke . (2016) nagelkerke2016Nagelkerke, E., Oberski, DL. Vermunt, JK.  2016\.  Goodness-of-fit of multilevel latent class models for categorical data Goodness-of-fit of multilevel latent class models for categorical data.  Sociological Methodology461252–282.
  * Nagelkerke . (2017) nagelkerke2017Nagelkerke, E., Oberski, DL. Vermunt, JK.  2017\.  Power and type I error of local fit statistics in multilevel latent class analysis Power and type I error of local fit statistics in multilevel latent class analysis.  Structural Equation Modeling: A Multidisciplinary Journal242216–229.
  * Nylund . (2007) nylund2007Nylund, KL., Asparouhov, T. Muthén, BO.  2007\.  Deciding on the number of classes in latent class analysis and growth mixture modeling: A Monte Carlo simulation study Deciding on the number of classes in latent class analysis and growth mixture modeling: A Monte Carlo simulation study.  Structural equation modeling: A multidisciplinary Journal144535–569.
  * Nylund-Gibson . (2019) nylund_masynNylund-Gibson, K., Grimm, R. Masyn, K.  2019\.  Prediction from Latent Classes: A Demonstration of Different Approaches to Include Distal Outcomes in Mixture Models Prediction from latent classes: A demonstration of different approaches to include distal outcomes in mixture models.  Structural Equation Modeling: A Multidisciplinary Journal266967-985.
  * Oberski Satorra (2013) oberskiOberski, D. Satorra, A.  2013\.  Measurement Error Models with Uncertainty about the Error Variance Measurement error models with uncertainty about the error variance.  Structural Equation Modeling: A Multidisciplinary Journal20.
  * Oser (2017) oser2017Oser, J.  2017\.  Assessing how participators combine acts in their “political tool kits”: A person-centered measurement approach for analyzing citizen participation Assessing how participators combine acts in their “political tool kits”: A person-centered measurement approach for analyzing citizen participation.  Social Indicators Research1331235–258.
  * Oser (2022) oser2022Oser, J.  2022\.  Protest as one political act in individuals’ participation repertoires: Latent class analysis and political participant types Protest as one political act in individuals’ participation repertoires: Latent class analysis and political participant types.  American Behavioral Scientist664510–532.
  * Oser Hooghe (2013) oser2013Oser, J. Hooghe, M.  2013\.  The Evolution of Citizenship Norms among Scandinavian Adolescents, 1999–2009 The evolution of citizenship norms among scandinavian adolescents, 1999–2009.  Scandinavian Political Studies364320–346.
  * Park Yu (2018) park2018Park, J. Yu, HT.  2018\.  Recommendations on the sample sizes for multilevel latent class models Recommendations on the sample sizes for multilevel latent class models.  Educational and Psychological Measurement785737–761.
  * Ranalli Rocci (2016) ranalli2016Ranalli, M. Rocci, R.  2016\.  Standard and Novel Model Selection Model Selection Criteria in the Pairwise Likelihood PEstimation of a Mixture Model for Ordinal Data Standard and novel model selection model selection criteria in the pairwise likelihood pestimation of a mixture model for ordinal data.  AF. Wilhelm HA. Kestler (), Analysis of Large and Complex Data Analysis of large and complex data ( 53–68).  Springer International Publishing.
  * Reboussin Ialongo (2010) reboussin2010Reboussin, B. Ialongo, N.  2010\.  Latent transition models with latent class predictors: attention deficit hyperactivity disorder subtypes and high school marijuana use Latent transition models with latent class predictors: attention deficit hyperactivity disorder subtypes and high school marijuana use.  Journal of the Royal Statistical Society: Series A (Statistics in Society)1731145–164.
  * Ruelens Nicaise (2020) ruelens2020Ruelens, A. Nicaise, I.  2020\.  Investigating a typology of trust orientations towards national and European institutions: A person-centered approach Investigating a typology of trust orientations towards national and european institutions: A person-centered approach.  Social Science Research102414.
  * Sandoval-Hernández . (2018) sandoval2018Sandoval-Hernández, A., Magdalena Isac, M. Miranda, D.  2018\.  Teaching tolerance in a globalized world Teaching tolerance in a globalized world.  Springer Nature.
  * Schulz . (2011) schulz2011Schulz, W., Ainley, J. Fraillon, J.  2011\.  ICCS 2009 technical report ICCS 2009 technical report.  Amsterdam, The Netherlands: International Association for the Evaluation of Educational Achievement (IEA).
  * Schulz, Ainley . (2018) schulz2018Schulz, W., Ainley, J., Fraillon, J., Losito, B., Agrusti, G. Friedman, T.  2018\.  Becoming citizens in a changing world: IEA International Civic and Citizenship Education Study 2016 international report Becoming citizens in a changing world: IEA International Civic and Citizenship Education Study 2016 international report.  Springer.
  * Schulz, Carstens . (2018) schulz2018IEASchulz, W., Carstens, R., Losito, B. Fraillon, J.  2018\.  ICCS 2016 technical report Iccs 2016 technical report.  Amsterdam: International Association for the Evaluation of Educational Achievement (IEA).
  * Schulz . (2004) schulz2004Schulz, W., Sibberns, H., Foy, P., Husfeldt, V., Lehmann, R., Nikolova, R. Torney-Purta, J.  2004\.  IEA civic education study technical report Iea civic education study technical report.  International Association for the Evaluation of Educaional Achivement.
  * Shehata . (2016) shehata2016Shehata, A., Ekström, M. Olsson, T.  2016\.  Developing self-actualizing and dutiful citizens: Testing the AC-DC model using panel data among adolescents Developing self-actualizing and dutiful citizens: Testing the AC-DC model using panel data among adolescents.  Communication Research4381141–1169.
  * Skrondal Kuha (2012) skrondal2012Skrondal, A. Kuha, J.  2012\.  Improved regression calibration Improved regression calibration.  Psychometrika774649–669.
  * Takeuchi (1976) takeuchi1976Takeuchi, K.  1976\.  Distribution of informational statistics and a criterion of model fitting. 153 12-18 Distribution of informational statistics and a criterion of model fitting. 153 12-18.  Suri-Kagaku (Mathematical Sciences)15312–18 (Japanese).
  * Torney-Purta (2002) torney2002Torney-Purta, J.  2002\.  The school’s role in developing civic engagement: A study of adolescents in twenty-eight countries The school’s role in developing civic engagement: A study of adolescents in twenty-eight countries.  Applied Developmental Science64203–212.
  * J. Vermunt (2003) vermunt.03Vermunt, J.  2003\.  Multilevel Latent Class Models Multilevel latent class models.  Sociological Methodology331213-239.
  * J. Vermunt (2008) vermunt2008Vermunt, J.  2008\.  Latent class and finite mixture models for multilevel data sets Latent class and finite mixture models for multilevel data sets.  Statistical Methods in Medical Research17133–51.
  * J. Vermunt (2010) vermunt:10Vermunt, J.  2010\.  Latent Class Modeling with Covariates: Two Improved Three-Step Approaches Latent class modeling with covariates: Two improved three-step approaches.  Political Analysis18450-469.
  * J. Vermunt Magidson (2016) LG51technicalVermunt, J. Magidson, J.  2016\.  Technical Guide for Latent GOLD 5.1: Basic, Advanced and Syntax Technical guide for Latent GOLD 5.1: Basic, Advanced and Syntax.  Belmont, MAStatistical Innovations.
  * J. Vermunt Van Dijk (2001) vermunt2001nonparametricVermunt, J. Van Dijk, L.  2001\.  A nonparametric random-coefficients approach: The latent class regression model A nonparametric random-coefficients approach: The latent class regression model.  Multilevel Modelling Newsletter1326–13.
  * JK. Vermunt Magidson (2021) vermunt2021Vermunt, JK. Magidson, J.  2021\.  How to perform three-step latent class analysis in the presence of measurement non-invariance or differential item functioning How to perform three-step latent class analysis in the presence of measurement non-invariance or differential item functioning.  Structural Equation Modeling: A Multidisciplinary Journal283356–364.
  * Vrieze (2012) vrieze2012Vrieze, SI.  2012\.  Model selection and psychological theory: a discussion of the differences between the Akaike information criterion (AIC) and the Bayesian information criterion (BIC). Model selection and psychological theory: a discussion of the differences between the akaike information criterion (aic) and the bayesian information criterion (bic).  Psychological Methods172228.
  * Welzel . (2005) welzel2005Welzel, C., Inglehart, R. Deutsch, F.  2005\.  Social capital, voluntary associations and collective action: which aspects of social capital have the greatest civic payoff? Social capital, voluntary associations and collective action: which aspects of social capital have the greatest civic payoff?  Journal of Civil Society12121–146.
  * White (1982) white1982maximumWhite, H.  1982\.  Maximum likelihood estimation of misspecified models Maximum likelihood estimation of misspecified models.  Econometrica: Journal of the Econometric Society5011–25.
  * Zhu . (2020) zhu2020Zhu, Y., Steele, F. Moustaki, I.  2020\.  A multilevel structural equation model for the interrelationships between multiple latent dimensions of childhood socio-economic circumstances, partnership transitions and mid-life health A multilevel structural equation model for the interrelationships between multiple latent dimensions of childhood socio-economic circumstances, partnership transitions and mid-life health.  Journal of the Royal Statistical Society: Series A (Statistics in Society)18331029–1050.