Linear Logistic Scoring Equations for Latent Class and Latent Profile Models: A Simple Method for Classifying New Cases

Abstract

Researchers are often interested in using latent class or latent profile parameter estimates to obtain posterior class membership probabilities for observations other than those of the original sample. In this paper, we demonstrate that these probabilities typically take on the form of linear logistic equations with coefficients which are functions of the original model parameters. In other words, the posterior class membership probabilities can be computed with a prediction formula similar to that of a multinomial logistic regression model. We derive the scoring equations for nominal, ordinal, count, and continuous indicators, as well as investigate models with missing values on class indicators, local dependencies, covariates, or multiple latent variables. In addition to the mathematical derivations of the scoring equations, we describe how either exact or approximate scoring equations can be obtained by estimating a multinomial regression model using a weighted data set.

Keywords:

• Canonical link functions
• classification
• multinomial logistic regression
• out-of-sample prediction
• posterior class membership probabilities

In applications of factor analysis, after selecting and estimating the factor model of interest, one will typically obtain (linear) factor-score equations which can be used to estimate the subjects’ factor scores as a function of the original items included in the model (Bartholomew et al., 2011). An important feature of factor-score equations is that these can be used not only for the subjects in the estimation sample, but also for new subjects, that is, for out-of-sample prediction.

When performing a latent class (LC) analysis, after selecting the final model, one may assign the individuals in the estimation sample to LCs using their posterior class membership probabilities. However, it is not well known that these posterior probabilities can be expressed exactly by a set of linear logistic equations, with “regression” weights which are functions of the original LC model parameters. More specifically, a closed-form expression for the posteriors exists, as a function of the LC model parameters, if the responses are modeled using distributions from the exponential family and with canonical link functions. Availability of a set of scoring equations makes it straightforward to compute the class membership probabilities for subjects which do not belong to the original sample used to estimate the LC model. In this way, one can realize an important goal of many LC analysis applications, namely obtaining out-of-sample class membership predictions. The main advantage of this approach is it allows predicting class memberships of new subjects without the need to use LC analysis software or to program the formula of the estimated LC model. Note that this is similar to what is done in factor analysis, where factor scores are obtained using a linear factor-score formula, without the need to return to the estimates of the factor covariances, factor loadings, and residual covariances.

As far as we know, LatentGOLD (Vermunt & Magidson, 2016; 2021) is currently the only software for LC analysis that allows one to obtain these logistic scoring equations, both in tabular form and in the form of SPSS or R syntax. The aim of this paper is to show how these equations are derived. As will be shown, the slopes of the linear logistic scoring equations are obtained easily, but the expression for the intercept terms (constants) may be somewhat more complex. In many situations, the equations for the posteriors will contain only main effects of the response variables. However, as in quadratic discriminant analysis, when a LC model for continuous responses assumes variances to be class specific, quadratic terms also need to be included, and when the LC model contains covariances/associations which are class specific, interactions are also required. The approach can be extended easily to LC models with covariates and multiple latent variables. More complicated are situations where responses contain missing values (in which case the constants need to be adapted to the missing data pattern), where the model contains direct effects of covariates on the responses (in which case the exact logistic form may no longer hold), and where non-canonical link functions are used (in which case there is no longer any direct relation between the LC model parameters and the scoring equation).

Rather than computing the scoring equations from the LC model parameters, one can also obtain these equations by estimating a multinomial logistic regression model using the posteriors as weights, as done in the LatentGOLD Step3-Scoring option (Vermunt & Magidson, 2016; 2021). This approach has the advantage of increased flexibility in that it is also possible to obtain approximate equations when exact closed-form solutions are not available, or when one prefers a simpler approximate set of scoring equations over more complex exact equations.

Below, we present the scoring equations for models for categorical responses (nominal, ordinal, and counts), models for continuous responses, models with local dependencies, models with covariates, models with missing values on responses, and models with multiple latent variables. We also discuss how the scoring equations can be obtained using a weighted multinomial logistic regression analysis.

1. Latent Class Models for Categorical Responses

Let $Y_j$ denote one of $J$ response variables (or indicators), with $1\le j\le J.$ A particular response for and the number of categories of the $j$th response variable are referred to as $y_j$ and $R_j,$ respectively, with $1\le y_j\le R_j.$ The probability of having a particular set of responses $\mathbf{y}$ is denoted by $P(\mathbf{Y}=\mathbf{y}).$ The discrete latent variable is denoted by $X,$ a particular latent class by $k,$ and the number of classes by $K.$

1.1. Nominal Responses

The standard LC model for nominal responses has the following form (Collins & Lanza, 2010; Goodman, 1974a; 1974b; Hagenaars, 1990; McCutcheon, 1987): $$P(\mathbf{Y}=\mathbf{y})=\sum _{k=1}^{K}P(X=k)\prod _{j=1}^{J}P(Y_j=y_j|X=k),$$ where $P(X=k)$ is the probability of belonging to class $k$ and $P(Y_j=y_j|X=k)$ the conditional probability of giving response $y_j$ on variable $Y_j$ conditional on belonging to class $k.$ These probabilities are often parameterized using logistic equations (Formann, 1992; Heinen, 1996; Magidson & Vermunt, 2004); that is, $$P\left(X=k\right)=\frac{\mathrm{\text{exp}}\hspace{0.17em}\left({\gamma }_k\right)}{D}$$ $$P(Y_j=y_j|X=k)=\frac{\mathrm{\text{exp}}({\alpha }_{y_j}+{\beta }_{y_{j}k})}{E_{\mathit{\text{jk}}}}$$

With $$D=\sum _{k=1}^K\mathrm{e}\mathrm{\text{xp}}\left({\gamma }_k\right)$$ $$E_{\mathit{\text{jk}}}=\sum _{y_j=1}^{R_j}\mathrm{e}\mathrm{\text{xp}}({\alpha }_{y_j}+{\beta }_{y_{j}k})$$

Here, ${\gamma }_k$ are intercept or constant terms in the regression model for $P(X=k),$ and ${\alpha }_{y_j}$ and ${\beta }_{y_{j}k}$ are intercept and slope parameters in the regression model for $P(Y_j=y_j|X=k).$ As always, identifying constraints need to be imposed on the logistic parameters. Typically, they are either restricted to sum to 0 over classes and response categories (referred to as effect coding), or set to 0 for one class and one response category (called dummy coding). The terms $D$ and $E_{\mathit{\text{jk}}}$ are normalizing constants.

The posterior probability of belonging to class $k$ conditional on response vector $\mathbf{y},$ denoted by $P(X=k|\mathbf{Y}=\mathbf{y}),$ can be obtained as follows (Dias & Vermunt, 2008; Goodman, 1974a; 1974b): $$P(X=k|\mathbf{Y}=\mathbf{y})=\frac{P(X=k)\prod _{j=1}^{J}P(Y_j=y_j|X=k)}{\sum _{k=1}^{K}P(X=k)\prod _{j=1}^{J}P(Y_j=y_j|X=k)}.$$

Replacing the model probabilities by their logit equations yields: $$P(X=k|\mathbf{Y}=\mathbf{y})=\frac{\frac{\mathrm{\text{exp}}\left({\gamma }_k\right)}{D}\prod _{j=1}^J\frac{\mathrm{\text{exp}}({\alpha }_{y_j}+{\beta }_{y_{j}k})}{E_{\mathit{\text{jk}}}}}{\sum _{k=1}^K\frac{\mathrm{\text{exp}}\left({\gamma }_k\right)}{D}\prod _{j=1}^J\frac{\mathrm{\text{exp}}({\alpha }_{y_j}+{\beta }_{y_{j}k})}{E_{\mathit{\text{jk}}}}}.$$

This equation can be simplified by removing $D$ and ${\alpha }_{y_j},$ which are redundant because they do not depend on $k.$ Moreover, the product over the $J$ responses can be replaced by a sum over the logs of the terms concerned. This yields $$P(X=k|\mathbf{Y}=\mathbf{y})=\frac{\mathrm{\text{exp}}({\gamma }_k-\sum _{j=1}^J\mathrm{l}\mathrm{\text{og}}(E_{\mathit{\text{jk}}})+\sum _{j=1}^J{\beta }_{y_{j}k})}{\sum _{k=1}^K\mathrm{e}\mathrm{\text{xp}}({\gamma }_k-\sum _{j=1}^J\mathrm{l}\mathrm{\text{og}}(E_{\mathit{\text{jk}}})+\sum _{j=1}^J{\beta }_{y_{j}k})}=\frac{\mathrm{\text{exp}}({\gamma }_k^{*}+\sum _{j=1}^J{\beta }_{y_{j}k})}{\sum _{k=1}^K\mathrm{e}\mathrm{\text{xp}}({\gamma }_k^{*}+\sum _{j=1}^J{\beta }_{y_{j}k})},$$ where ${\gamma }_k^{*}={\gamma }_k-{\sum }_{j=1}^J\mathrm{l}\mathrm{\text{og}}(E_{\mathit{\text{jk}}}).$ Though not really necessary, the ${\gamma }_k^{*}$ parameters may be subjected to the same identifying—effect or dummy coding—constraints as the other parameters.

The above derivation shows that the posterior class membership probabilities can be written as a logistic equation with slopes equal to the LC model logistic regression slopes and with constants equal to the logistic class constants minus the sum of the logs of the normalizing constants. The more difficult part in the computation of these scoring equations is the computation of the constants ${\gamma }_k^{*}.$ But, once we have the scoring equations, we can easily compute the class membership probabilities for any response pattern, including response patterns which were not available in the original data used to estimate the LC model of interest.

It should be noted that the ${\gamma }_k^{*}$ terms are identical to the one-variable parameters for the latent classes in the log-linear formulation of the LC model proposed by Haberman (1979). In this formulation, the joint distribution of the $X$ and $\mathbf{Y},$ $P(X=k,\mathbf{Y}=\mathbf{y}),$ is modelled as follows: $$P(X=k,\mathbf{Y}=\mathbf{y})=\frac{\mathrm{\text{exp}}({\gamma }_k^{*}+\sum _{j=1}^J{\alpha }_{y_j}+\sum _{j=1}^J{\beta }_{y_{j}k})}{F}.$$

The posterior class membership probability is obtained as $P(X=k|\mathbf{Y}=\mathbf{y})=P(X=k,\mathbf{Y}=\mathbf{y})/{\sum }_{k=1}^{K}P(X=k,\mathbf{Y}=\mathbf{y}).$ As can be seen, since ${\sum }_{j=1}^J{\alpha }_{y_j}$ and $F$ cancel, $P(X=k|\mathbf{Y}=\mathbf{y})$ has exactly the form we derived above. Thus, when using Haberman’s log-linear formulation, the constants of the scoring equations are also model parameters. However, an important disadvantage of this formulation is that it is computationally less efficient since parameter estimation involves processing the cell entries in the joint cross-tabulation of $X$ and all $Y_j$ variables. Therefore, this log-linear approach can be used only when the number of response variables is small.

1.2. Ordinal Responses and Counts

In the LatentGOLD program for LC analysis (Vermunt & Magidson, 2016; 2021), ordinal response variables can be modeled using an adjacent-category logit model, that is, using a canonical link function (Agresti, 2002). More specifically, these are multinomial logit models in which the class-indicator association parameters are restricted as follows: ${\beta }_{y_{j}k}={\beta }_{\mathit{\text{jk}}}{y}_j;$ that is, to be nominal-by-linear (Goodman, 1979; Heinen, 1996). This implies that for ordinal variables, $$P(Y_j=y_j|X=k)=\frac{\mathrm{\text{exp}}({\alpha }_{y_j}+{\beta }_{\mathrm{j}k}{y}_j)}{E_{\mathit{\text{jk}}}}.$$

The same restrictions imposed on the ${\beta }_{y_{j}k}$ parameters also apply to the scoring equations; that is, for ordinal variables, we replace the ${\beta }_{{}_{j}k}$ terms in the scoring equations by ${\beta }_{\mathit{\text{jk}}}{y}_j.$ This shows that in the ordinal case, the class-membership logits are linear functions of the item responses. It should be noted that while we assumed that the category scores ranged from 1 to $R_j,$ the adjacent-category logit model allows for any type of scoring. In its more general form, ${\beta }_{y_{j}k}={\beta }_{\mathit{\text{jk}}}{\nu }_{y_j},$ where ${\nu }_{y_j}$ is the score for category $y_j.$

When modeling ordinal variables using other (non-canonical) link functions, such as cumulative logit or cumulative probit link functions, exact expressions for the scoring equations no longer exist. As will be shown below, a possible way out is to estimate the scoring equations treating the response variables as either numeric or nominal predictors of class membership.

Also for a Poisson and binomial count variable, the scoring equations contain the term ${\beta }_{\mathit{\text{jk}}}{y}_j.$ This can be seen from the fact that the class-specific density of a count variable $Y_j$ takes on the following form: $$P(Y_j=y_j|X=k)\propto \frac{\mathrm{\text{exp}}({\alpha }_j{y}_j+{\beta }_{\mathit{\text{jk}}}{y}_j)}{E_{\mathit{\text{jk}}}},$$ where it should be noted that ${\alpha }_{{}_j}{y}_j$ cancels from the scoring equation. The expression for $\mathrm{\text{log}}{E}_{\mathit{\text{jk}}}$ changes compared to the nominal and ordinal case. For Poisson counts, $\mathrm{\text{log}}{E}_{\mathit{\text{jk}}}=\mathrm{\text{exp}}({\alpha }_j+{\beta }_{\mathit{\text{jk}}})e_j,$ where $e_j$ is the exposure; and for Binomial counts, $\mathrm{\text{log}}{E}_{\mathit{\text{jk}}}=e_j\mathrm{l}\mathrm{\text{og}}(1+\mathrm{\text{exp}}({\alpha }_j+{\beta }_{\mathit{\text{jk}}})),$ where $e_j$ represents the number of trials. This shows that scoring equations should also include terms for the exposure (number of trials) when this number varies across individuals. When the $e_j$ are fixed, as for nominal and ordinal variables, the $\mathrm{\text{log}}{E}_{\mathit{\text{jk}}}$ terms can be included in the constants ${\gamma }_k^{*}.$

1.3. Local Dependencies

Thus far, we assumed that responses are independent within classes. Now we will look at the scoring equations for LC models with local dependencies (Hagenaars, 1988; Magidson & Vermunt, 2004; Oberski et al., 2013). In the most general case, including a local dependency between (nominal) response variables $Y_j$ and $Y_m$ implies that $$P(Y_j=y_j,Y_m=y_m|X=k)=\frac{\mathrm{\text{exp}}({\alpha }_{y_j}+{\alpha }_{y_m}+{\beta }_{y_{j}k}+{\beta }_{y_{m}k}+{\delta }_{y_j{y}_m}+{\lambda }_{y_j{y}_{m}k})}{E_{\mathit{\text{jmk}}}},$$ where $$E_{\mathit{\text{jmk}}}=\sum _{y_j=1}^{R_j}\sum _{y_m=1}^{R_m}\mathrm{e}\mathrm{\text{xp}}({\alpha }_{y_j}+{\alpha }_{y_m}+{\beta }_{y_{j}k}+{\beta }_{y_{m}k}+{\delta }_{y_j{y}_m}+{\lambda }_{y_j{y}_{m}k}).$$

This is a model with an association between $Y_j$ and $Y_m,$ ${\delta }_{y_j{y}_m},$ and an interaction with the latent classes, ${\lambda }_{y_j{y}_{m}k}.$ In other words, it represents a model in which the strength of the local dependency is allowed to vary across classes.

As can be seen, one difference with the local independence model is that the normalizing constants entering in the ${\gamma }_k^{*}$ coefficients of the scoring equations should be computed per set of locally dependent variables. The ${\delta }_{y_j{y}_m}$ term cancels from the scoring equations because it does not depend on the classes. In contrast, the term ${\lambda }_{y_j{y}_{m}k}$ becomes part of the scoring equations, which in the case of class-specific local dependencies not only contains main effects but also interaction terms. Note that when local dependencies are not class-specific, that is, when ${\lambda }_{y_j{y}_{m}k}=0,$ the only remaining difference between local independence and local dependence models concerns the computation of the constants ${\gamma }_k^{*}.$

The scoring equations in local-dependence LC models for ordinal variables are very similar to those for nominal variables. When the ordinal variables are modelled using an adjacent-category logit specification, ${\delta }_{y_j{y}_m}={\delta }_{\mathit{\text{jm}}}{y}_j{y}_m$ and ${\lambda }_{y_j{y}_{m}k}={\lambda }_{\mathit{\text{jmk}}}{y}_j{y}_m.$ The scoring equations will contain the term ${\lambda }_{\mathit{\text{jmk}}}{y}_j{y}_m$ when the interaction parameters ${\lambda }_{\mathit{\text{jmk}}}$ are not fixed to 0.

1.4. Missing Values

When some indicators have missing values, the LC model for the observed values ${\mathbf{Y}}_{\mathit{\text{obs}}}$ can be defined as follows: $$P({\mathbf{Y}}_{\mathit{\text{obs}}}={\mathbf{y}}_{\mathit{\text{obs}}})=\sum _{k=1}^{K}P(X=k)\prod _{j=1}^{J}P(Y_j=y_j|X=k{)}^{r_j},$$ where $r_j=1$ if the response variable concerned is observed and 0 when it has a missing value (Vermunt et al., 2008; Vermunt & Magidson, 2016). Note that this formulation implies that the product is taken over the observed responses only. Therefore, similar to subjects with complete data, the computation of the posteriors for subjects with missing values involves using only their observed responses. This means that the sum ${\sum }_{j=1}^J{\beta }_{y_{j}k}$ should be taken over the observed variables only or, equivalently, that ${\beta }_{y_{j}k}$ should be set to 0 for the missing value category. However, the sum ${\sum }_{j=1}^J\mathrm{l}\mathrm{\text{og}}(E_{\mathit{\text{jk}}})$ which is subtracted from the constants should also be taken over the observed variables only, implying that each pattern of missing data has its own set of constants ${\gamma }_k^{*}.$ A way to account for this is by using the same ${\gamma }_k^{*}$ for all observations but adding a term $\mathrm{\text{log}}(E_{\mathit{\text{jk}}})$ to the scoring equation when variable $j$ has a missing value. In other words, in order to deal with missing data, the scoring equation should be expanded to include the term ${\sum }_{j=1}^J\mathrm{l}\mathrm{\text{og}}(E_{\mathit{\text{jk}}})(1-r_j).$ Note that this approach can be used with any missing data pattern occurring among the new subjects for which one wishes to obtain the posteriors.

A special type of missing data occurs when the LC model is estimated using $J$ variables, but only the first $J_1$ of these are to be used for classification purposes (where $J=J_1+J_2$); for example, this situation may occur if one wishes to ignore the last $J_2$ variables when calculating the classification probabilities because this information will not be available when performing out-of-sample predictions. In this case, the posteriors are obtained as follows: $$P(X=k|{\mathbf{Y}}_1={\mathbf{y}}_1)=\frac{P(X=k)\prod _{j=1}^{J_1}P(Y_j=y_j|X=k)}{\sum _{k=1}^{K}P(X=k)\prod _{j=1}^{J_1}P(Y_j=y_j|X=k)}.$$

As can be seen, only the slope parameters and the normalizing constants of the first $J_1$ response variables will enter into the scoring equations.

1.5. An Example with Five Dichotomous Indicators

Table 1 provides an example illustrating the computation of the scoring equations for an application with five dichotomous response variables. It concerns the model with three latent classes estimated for the LatentGOLD “political.sav” demo data set. The upper part of Table 1 gives the estimates of the model parameters ${\gamma }_k,$ ${\alpha }_{y_j},$ and ${\beta }_{y_{j}k}$ using dummy coding with the parameters for the first class and the first item category fixed to 0. The lower part gives the values of ${\gamma }_k^{*}$ and $\mathrm{\text{log}}(E_{\mathit{\text{jk}}}),$ where for consistency we also use dummy coding for the $\mathrm{\text{log}}\left(E_{\mathit{\text{jk}}}\right)$ terms. To obtain the $\mathrm{\text{log}}(E_{1k})$ values for the first item, we first compute ${\mathrm{\text{exp}}(\alpha }_{y_1}+{\beta }_{y_{1}k}),$ which for all three classes equals 1.0000 for $y_1=1,$ and 2.6533, 0.4452, and 1.4312 for $y_1=2.$ Next, we sum the obtained values across the two item categories and take the log, yielding 1.29562, 0.3682, and 0.8884 for the three classes. Because of the dummy coding, we subtract the value of the first class, yielding the reported $\mathrm{\text{log}}(E_{1k})$ values 0.0000, −0.9275, and −0.4073. The $\mathrm{\text{log}}\left(E_{\mathit{\text{jk}}}\right)$ values of all items are subtracted from the ${\gamma }_k$ values to get the intercepts ${\gamma }_k^{*}$ of the scoring equations, and the slopes ${\beta }_{y_{j}k}$ can be used without any modification in the scoring equations. In the case of a missing value, the slope parameters for the item concerned equal $\mathrm{\text{log}}\left(E_{\mathit{\text{jk}}}\right).$ Appendix A shows R code generated for this application by LatentGOLD, which can be used for classifying new observations.

Table 1. Latent class model parameters and scoring equation parameters for the political.sav data example.

LC model parameters
		Intercept X = 1	Intercept X = 2	Intercept X = 3
X		0.0000	−0.0723	−0.5173
	Intercept	Slopes X = 1	Slopes X = 2	Slopes X = 3
Y1 = 1	0.0000	0.0000	0.0000	0.0000
Y1 = 2	0.9758	0.0000	−1.7853	−0.6173
Y2 = 1	0.0000	0.0000	0.0000	0.0000
Y2 = 2	−0.3534	0.0000	−3.0502	−0.2328
Y3 = 1	0.0000	0.0000	0.0000	0.0000
Y3 = 2	−1.7062	0.0000	0.5660	3.6819
Y4 = 1	0.0000	0.0000	0.0000	0.0000
Y4 = 2	0.2028	0.0000	−0.7463	3.0609
Y5 = 1	0.0000	0.0000	0.0000	0.0000
Y5 = 2	2.0140	0.0000	−3.0398	−1.0034
Scoring equation parameters
		Intercept X = 1	Intercept X = 2	Intercept X = 3
X		0.0000	3.4186	−3.6425
		E values for X = 1	E values for X = 2	E values for X = 3
Y1		0.0000	−0.9275	−0.4073
Y2		0.0000	−0.4993	−0.0896
Y3		0.0000	0.1106	1.9387
Y4		0.0000	−0.3418	2.5015
Y5		0.0000	−1.8329	−0.8183

2. Other Types of Latent Class and Mixture Models

2.1. Continuous Responses

Now, let us turn to LC or mixture models for continuous response variables (McLachlan & Peel, 2000), also referred as latent profile models. In a local independence model with normal within-class distributions with possibly unequal variances, the response distributions have the following form: $$P(Y_j=y_j|X=k)\propto \mathrm{\text{exp}}(-\frac{1}{2}\mathrm{l}\mathrm{\text{og}}{\sigma }_{\mathit{\text{jk}}}^2-\frac{1}{2}\frac{{\mu }_{\mathit{\text{jk}}}{\mu }_{\mathit{\text{jk}}}}{{\sigma }_{\mathit{\text{jk}}}^2}+\frac{{\mu }_{\mathit{\text{jk}}}{y}_j}{{\sigma }_{\mathit{\text{jk}}}^2}-\frac{1}{2}\frac{y_j{y}_j}{{\sigma }_{\mathit{\text{jk}}}^2}),$$ where ${\mu }_{\mathit{\text{jk}}}$ and ${\sigma }_{\mathit{\text{jk}}}^2$ denote the mean and variance of $Y_j$ in latent class $k.$ It can be seen that in the construction of the logistic scoring equations, the terms $-\frac{1}{2}\mathrm{l}\mathrm{\text{og}}{\sigma }_{\mathit{\text{jk}}}^2$ and $-\frac{1}{2}\frac{{\mu }_{\mathit{\text{jk}}}{\mu }_{\mathit{\text{jk}}}}{{\sigma }_{\mathit{\text{jk}}}^2},$ which do not contain the response, become part of the constants. Moreover, the equations will contain the linear and quadratic terms $\frac{{\mu }_{\mathit{\text{jk}}}}{{\sigma }_{\mathit{\text{jk}}}^2}{y}_j$ and $-\frac{1}{2{\sigma }_{\mathit{\text{jk}}}^2}{y}_j{y}_j.$

When variances are assumed to be equal across classes, the first and last term of the above univariate normal distribution become $-\frac{1}{2}\mathrm{l}\mathrm{\text{og}}{\sigma }_j^2$ and $-\frac{1}{2}\frac{y_j{y}_j}{{\sigma }_j^2},$ respectively, implying that these cancel from the scoring equations because they do not depend on the class. This yields a set of scoring equations similar to those obtained in linear discriminant analysis (Hastie et al., 2008).

In the more general case of multivariate normal responses with unrestricted covariances ${\mathrm{\Sigma }}_k,$ the LC model becomes $$P\left(\mathbf{Y}=\mathbf{y}\right)=\sum _{k=1}^{K}P\left(X=k\right)P\left(\mathbf{Y}=\mathbf{y}\right|X=k),$$ with $$P(\mathbf{Y}=\mathbf{y}|X=k)\propto \mathrm{\text{exp}}(-\frac{1}{2}\mathrm{l}\mathrm{\text{og}}\left|{\mathbf{\Sigma }}_k\right|-\frac{1}{2}{\mathbf{\mu }}_k^{\prime }{\mathbf{\Sigma }}_k^{-1}{\mathbf{\mu }}_k+{\mathbf{\mu }}_k^{\prime }{\mathbf{\Sigma }}_k^{-1}\mathbf{y}-\frac{1}{2}{\mathbf{y}}^{\prime }{\mathbf{\Sigma }}_k^{-1}\mathbf{y}).$$

As can be seen, the scoring equations now not only contain linear and quadratic terms, but also interaction terms are needed. More specifically, denoting an entry of ${\mathbf{\Sigma }}_k^{-1}$ by $a_{\mathit{\text{mjk}}},$ the weights for $y_j,$ $y_j^2,$ and $y_j{y}_m\mathrm{ }$ are ${\sum }_{m=1}^J\mathrm{ }{\mu }_{\mathit{\text{mk}}}{a}_{\mathit{\text{jmk}}},$ $-\frac{1}{2}{a}_{\mathit{\text{jjk}}},$ and $-a_{\mathit{\text{jmk}}},$ respectively. The first two terms of the multivariate normal density, which do not depend on the responses, become part of the constants. When variances and covariances are equal across classes, we again have equations with main effects only, as in linear discriminant analysis.

Various kinds of restricted mixtures of multivariate normal distributions have been proposed in which constraints are imposed on the class-specific means and/or covariances. Examples include mixture factor models (McLachlan & Peel, 2000; Yung, 1997), mixture structural equation models (Dolan & Van der Maas, 1997), mixture models with constrained eigenvalue decompositions of ${\mathrm{\Sigma }}_k$ (Banfield & Raftery, 1993), and mixture growth models (Muthén, 2004). For these models, the same scoring equations can be used as when means and covariances are unrestricted.

2.2. An Example with Three Continuous Indicators

Table 2 reports the model parameters and the scoring equations for a LC model with three continuous indicators (Glucose, Insulin, and SSPG) from the LatentGOLD “diabetes.dat” demo data set. It is a three-class model with a free residual covariance between the first two class indicators and with class-specific residual (co)variances. The scoring equation for this model contains not only linear terms, but also quadratic terms as well as the interaction terms between $y_1$ and $y_2.$

Table 2. Latent class model parameters and scoring equation parameters for the diabetes.dat data example.

LC model parameters
	Intercept X = 1	Intercept X = 2	Intercept X = 3
X	0	−0.6927	−1.036
	Parameters X = 1	Parameters X = 2	Parameters X = 3
Mean Y1	91.2315	104.0049	234.7598
Mean Y2	359.2211	495.0568	1121.0893
Mean Y3	163.1271	309.4323	76.9772
Variance Y1	76.4770	230.0891	5005.9106
Variance Y2	2669.7454	14844.5520	73551.0945
Variance Y3	2421.4506	22966.5152	2224.5020
Covariance Y1-Y2	96.4624	1279.9240	17910.7089
Scoring equation parameters
	Intercept X = 1	Intercept X = 2	Intercept X = 3
X	0.0000	42.6430	56.8566
	Slopes X = 1	Slopes X = 2	Slopes X = 3
Y1	0.0000	−0.5599	−1.1314
Y2	0.0000	−0.1066	−0.0661
Y3	0.0000	−0.0539	−0.0328
Y1*Y1	0.0000	0.0027	0.0061
Y2*Y2	0.0000	0.0002	−0.0001
Y3*Y3	0.0000	0.0001	0.0001
Y1*Y2	0.0000	0.0002	0.0000

2.3. Covariates

When covariates are included in the model, the latent class probabilities are typically modeled as a logistic function of the covariates (Bandeen-Roche et al, 1997; Dayton & Macready, 1988; Yamaguchi, 2000). That is, $$P(X=k|\mathbf{z})=\frac{\mathrm{\text{exp}}({\gamma }_{0k}+\sum _{p=1}^P{\gamma }_{\mathit{\text{pk}}}{z}_p)}{\sum _{k=1}^K\mathrm{e}\mathrm{\text{xp}}({\gamma }_{0k}+\sum _{p=1}^P{\gamma }_{\mathit{\text{pk}}}{z}_p)}.$$

Here, $\mathbf{z}$ denotes the vector of covariates, and ${\gamma }_{0k}$ and ${\gamma }_{\mathit{\text{pk}}}$ represent the constants and the regression parameters for covariate $z_p.$

Since the denominator does not depend on the class, it cancels from the formula for the posterior class membership probability and thus also from the scoring equations. Assuming the response variables are nominal, the posterior probability of class membership given responses and covariates becomes: $$P(X=k|\mathbf{z})=\frac{\mathrm{\text{exp}}({\gamma }_{0k}^{*}+\sum _{j=1}^J{\beta }_{y_{j}k}+\sum _{p=1}^P{\gamma }_{\mathit{\text{pk}}}{z}_p)}{\sum _{k=1}^K\mathrm{e}\mathrm{\text{xp}}({\gamma }_{0k}^{*}+\sum _{j=1}^J{\beta }_{y_{j}k}+\sum _{p=1}^P{\gamma }_{\mathit{\text{pk}}}{z}_p)}.$$

That is, the covariate terms can simply be added to the scoring equations.

Covariates may also have direct effects on the indicators. Let us assume we have a single covariate $z$ which has a direct effect on the categorical response variable $Y_j;$ that is, $$P(Y_j=y_j|X=k,z)=\frac{\mathrm{\text{exp}}({\alpha }_{y_j}+{\beta }_{y_{j}k}+{\delta }_{y_j}z)}{E_{\mathit{\text{jk}}|z}},$$ where $$E_{\mathit{\text{jk}}|z}=\sum _{y_j=1}^{R_j}\mathrm{e}\mathrm{\text{xp}}({\alpha }_{y_j}+{\beta }_{y_{j}k}+{\delta }_{y_j}z).$$

As can be seen, in this model, the normalizing constants depend on the covariate value, meaning that we no longer have a single $\mathrm{\text{log}}{E}_{\mathit{\text{jk}}}$ which can be subtracted from ${\gamma }_k.$ Because the $\mathrm{\text{log}}{E}_{\mathit{\text{jk}}|z}$ are neither linear functions of the covariate values, they cannot be added to the linear term for the covariate concerned. In other words, the exact linear logistic representation of the posterior probabilities collapses in this situation, though, as discussed below and in Appendix B, it may still be used as an approximation. An exception is the situation in which the covariate is a nominal or dichotomous variable, in which case exact scoring equations can still be obtained by subtracting $\mathrm{\text{log}}{E}_{\mathit{\text{jk}}|z}$ from the ${\gamma }_{\mathit{\text{pk}}}$ terms of the covariate concerned.

Note that when the direct effect of a covariate is allowed to be class specific, or equivalently, when an interaction term is included, the indicator-covariate interaction should also be added to the scoring equation. Again, the scoring equations will be exact only when the covariate concerned is nominal or dichotomous. For example, this is the specification used in multiple-group LC models in which response probabilities may be allowed to differ across subgroups for one or more indicators (Clogg & Goodman, 1984; Eid et al., 2003; Kankaras et al., 2010).

2.4. Multiple Latent Variables

Suppose the LC model contains two latent variables $X_1$ and $X_2$ instead of one, so that we have a LC Factor or Discrete Factor model. Such a model, has the following form (Goodman, 1974b; Hagenaars, 1990; Magidson & Vermunt, 2001; Vermunt & Magidson, 2005): $$P(\mathbf{Y}=\mathbf{y})=\sum _{k_1=1}^{K_1}\sum _{k_2=1}^{K_2}P(X_1=k_1,X_2=k_2)\prod _{j=1}^{J}P(Y_j=y_j|X_1=k_1,X_2=k_2).$$

As in the single latent variable case, the $P(X_1=k_1,X_2=k_2)$ and $P(Y_j=y_j|X_1=k_1,X_2=k_2)$ can be modelled using logistic regression models (Magidson & Vermunt, 2001). For example, $$P\left(X_1=k_1,X_2=k_2\right)=\frac{\mathrm{\text{exp}}\hspace{0.17em}\left({\gamma }_{k_1}^1+{\gamma }_{k_2}^2+{\gamma }_{k_1{k}_2}^{12}\right)}{D}$$ $$P(Y_j=y_j|X_1=k_1,X_2=k_2)=\frac{\mathrm{\text{exp}}({\alpha }_{y_j}+{\beta }_{y_j{k}_1}^1+{\beta }_{y_j{k}_2}^2)}{E_{j{k}_1{k}_2}}.$$

Also in this case, the posterior probabilities can be written as functions of the LC model parameters; i.e., $$P(X_1=k_1,X_2=k_2|\mathbf{Y}=\mathbf{y})=\frac{\mathrm{\text{exp}}({\gamma }_{k_1{k}_2}^{*}+\sum _{j=1}^J{\beta }_{y_j{k}_1}^1+\sum _{j=1}^J{\beta }_{y_j{k}_2}^2)}{\sum _{k_1=1}^{K_1}\sum _{k_2=1}^{K_2}\mathrm{e}\mathrm{\text{xp}}({\gamma }_{k_1{k}_2}^{*}+\sum _{j=1}^J{\beta }_{y_j{k}_1}^1+\sum _{j=1}^J{\beta }_{y_j{k}_2}^2)},$$ where the ${\gamma }_{k_1{k}_2}^{*}$ contain the $\gamma$ terms and the normalizing constants $E_{j{k}_1{k}_2}.$

Note that the above equation for two latent variables can easily be generalized to an arbitrary number of $S$ latent variables. The logistic scoring equation then becomes: $$P(X_1=k_1,\dots ,X_S=k_S|\mathbf{Y}=\mathbf{y})=\frac{\mathrm{\text{exp}}({\gamma }_{k_1\dots k_S}^{*}+\sum _{s=1}^S\sum _{j=1}^J{\beta }_{y_j{k}_s}^s)}{\sum _{k_1=1}^{K_1}\dots \sum _{k_S=1}^{K_S}\mathrm{e}\mathrm{\text{xp}}({\gamma }_{k_1\dots k_S}^{*}+\sum _{s=1}^S\sum _{j=1}^J{\beta }_{y_j{k}_s}^s)}.$$

It should be noted that the marginal posterior probabilities $P(X_s=k_s|\mathbf{Y}=\mathbf{y}),$ which are obtained by collapsing over the other latent variables, cannot be written as logistic functions. However, logistic approximations of the marginal posteriors may be precise enough in most applications. Below, we discuss how such approximations can be obtained.

3. Estimating the Scoring Equations Using Logistic Regression Analysis

Rather than computing the scoring equations from the parameters of the LC model, it is also possible to obtain these equations posthoc using a standard routine for multinomial logistic regression analysis. This involves the following three steps:

1. After selection of the final LC model, save the posterior class membership probabilities to an output file. This is a feature available in all software packages for LC analysis.

2. Create an expanded data set with $K$ records per subject, which contains a column with the class number taking on values from $1$ to $K,$ a column with the posterior probability for the person and class concerned, and columns for the response variables and covariates used in the LC model, the latter columns containing the same values repeated in each of the K records for each subject.

3. Estimate a logistic regression model in which the posteriors are used as weights. The class number is the dependent variable, and the responses and covariates are the predictors.

Depending on the situation, in the third step the responses and covariates are modeled as either nominal or numeric predictors, quadratic and/or interaction effects are added, and missing value dummies are included. For count variables, one should include the exposures (or total number of trials) as additional numeric predictors when these differ across individuals. Steps 2 and 3 are automated in the LatentGOLD program (Vermunt & Magidson, 2016; 2021), and are called Step3-Scoring.

This approach can be used not only for the posthoc computation of the exact scoring equations, but also for obtaining approximate scoring equations. This is useful when an exact form does not exist, such as when direct effects of numeric covariates on indicators were included in the LC model or when non-canonical link functions were used for ordinal variables, as well as the situation where one prefers a set of simplified equations, say without quadratic or interaction terms, that are almost as good as the exact ones. An example of the latter can be seen in Table 3, which reports the approximate scoring equations for the diabetes.dat example presented above, but leaving out the quadratic terms of $y_1$ and $y_2$ and their interaction terms. The approximate equations predict the class memberships almost as well as the exact equations; that is, the entropy R² equals 0.817, while its original value equals 0.833.

Table 3. Approximate scoring equation parameters for the diabetes.dat data example.

	Intercept X = 1	Intercept X = 2	Intercept X = 3
X	0.0000	−9.8317	−25.7983
	Slopes X = 1	Slopes X = 2	Slopes X = 3
Y1	0.0000	0.0435	0.0626
Y2	0.0000	0.0236	0.0416
Y3	0.0000	−0.0816	−0.0173
Y3*Y3	0.0000	0.0003	−0.0001

When ordinal variables are modeled using non-canonical link functions, such as cumulative logit or probit models, we have two options. Option 1 is to compute the exact scoring equations by treating the response variables as nominal predictors in the posthoc logistic regression analysis; that is, by making use of the fact that the estimated $P(Y_j=y_j|X=k)$ based on an ordinal model can be reproduced perfectly by an unrestricted multinomial model. Option 2 is to estimate the scoring equations using the response variables as numeric predictors, which in fact implies that the estimated $P(Y_j=y_j|X=k)$ from the original LC model are approximated by an adjacent-category logit model.

As shown above, in LC models with multiple latent variables, an exact set of logistic scoring equations exists for the joint class membership probabilities, but not for the marginal class membership probabilities. The posthoc estimation method can also be used to obtain approximate scoring equations based on the marginal posteriors. Applying these equations will be simpler than first computing the joint and subsequently collapsing over the other latent variables, especially with models containing more than two latent variables. The quality of the resulting approximation can be assessed by a goodness-of-fit measure.

4. Discussion

As in continuous latent variables models, in LC models it is important to have a simple scoring rule for predicting a person’s value on the latent variable. In this paper, we showed that for LC models this scoring rule has the form of a linear logistic equation, with weights which are simple functions of the original LC model parameters. We derived the exact scoring equations for nominal, ordinal, count, and continuous response variables, for local independence and local dependence models, for models with covariates, for models with multiple latent variables, and for models with missing values on some of the indicators. Moreover, we discussed several situations in which exact scoring equations may not exist, such as LC models with direct effects of covariates on the indicators and LC models in which the conditional response distributions are restricted using regression models based on non-canonical link functions.

We also explained how to compute exact or approximate scoring equations with the saved posterior probabilities from any LC analysis program. This can be achieved with standard routines for logistic regression analysis. In practice, this may be much easier than computing the scoring equations from the LC model parameters, where the constants from the scoring equations may be somewhat more tedious to obtain.

While not discussed explicitly, the computation of the scoring equations proceeds in exactly the same manner in LC models for mixed responses; that is, in LC models for combinations of nominal, ordinal, count, and continuous indicators (Hennig & Liao, 2013; Hunt & Jorgensen, 1999; Vermunt & Magidson, 2002). The only thing that needs to be done in the computation of the scoring equations is to collect the terms for the different indicators, irrespective of their scale types. When using the posthoc method based on a logistic regression analysis, things are even easier. Nominal indicators are used as nominal predictors, and ordinal, count and continuous indicators as numeric predictors. Depending on the situation, quadratic and/or interaction terms may also need to be included.

The scoring equations discussed in this article can be used to obtain point estimates for the posterior probabilities, not only for subjects in the original sample, but also for new subjects. However, an issue not dealt with in this paper is the uncertainty about those estimates. Since the “regression” weights of the scoring equation are sample estimates, when deriving a prediction, it would be better to take into account this sampling variability. Note that the weights are functions of the original model parameters, for which we have the estimated asymptotic variance-covariance matrix. A possible approach to obtain the covariance matrix of the weights involves sampling say 100 parameter sets from their estimated multivariate normal distribution and computing 100 sets of weights. Other options to explore are the delta method and bootstrapping (Dias & Vermunt, 2008). Our future research will focus on this important topic.

Disclosure Statement

Jeroen K. Vermunt and Jay Magidson are co-developers of the LatentGOLD software.

Appendix A:

R Code Generated by LatentGOLD for the First Example Application

With the LatentGOLD output options “ScoringEquations” and “WriteRsyntax= <filename>“, one can request an R syntax file that can be used to classify new observations. The variable names in the “political.sav” data file are sys_resp, ideo_lev, rep_pot, prot_app, and conv_par. The lg_scoring function consists of three parts:

1. First, it creates the variables to be used as “predictors” in the scoring equations. For ordinal and continuous indicators, these are copies of the variables in the data set, and for categorical variables, these are dummies for the response categories. In addition, dummies are created for missing values.

2. Then, it computes the class-specific linear terms using the variables created in part 1 and the scoring equations’ parameters.

3. Subsequently, the linear terms are exponentiated and transformed to posterior probabilities.

The function’s output are the posterior probabilities and the modal class. Below, you find the R code, which ends with example code calling the lg_scoring function to add classification information to a data set.

## Scoring function to be called per record

lg_scoring<-function(dat) {

# Part 1: Create variables to be used as predictors in scoring

# equations

if(is.na(dat$sys_resp)) {

 sys_resp_lg_1<-0;sys_resp_lg_2<-0;sys_resp_lg_m<-1

}

else {

 if(dat$sys_resp==1) {

  sys_resp_lg_1<-1;sys_resp_lg_2<-0;sys_resp_lg_m<-0

 }

 else if(dat$sys_resp==2) {

  sys_resp_lg_1<-0;sys_resp_lg_2<-1;sys_resp_lg_m<-0

 }

 else {

  sys_resp_lg_1<-0;sys_resp_lg_2<-0;sys_resp_lg_m<-1

 }

}

# The same is done for the other 4 indicators 

…

# Part 2: Compute the class-specific linear terms 

Cluster_lg_1<-(0)+

  (0)*sys_resp_lg_1+(0)*sys_resp_lg_2+

  (0)*ideo_lev_lg_1+(0)*ideo_lev_lg_2+

  (0)*rep_pot_lg_1+(0)*rep_pot_lg_2+

  (0)*prot_app_lg_1+(0)*prot_app_lg_2+

  (0)*conv_par_lg_1+(0)*conv_par_lg_2+

  (0)*sys_resp_lg_m+(0)*ideo_lev_lg_m+

  (0)*rep_pot_lg_m+(0)*prot_app_lg_m+

  (0)*conv_par_lg_m

Cluster_lg_2<-(3.4185551)+

  (0)*sys_resp_lg_1+(-1.7853117)*sys_resp_lg_2+

  (0)*ideo_lev_lg_1+(-3.0502076)*ideo_lev_lg_2+

  (0)*rep_pot_lg_1+(0.56595846)*rep_pot_lg_2+

  (0)*prot_app_lg_1+ (-0.74630356)*prot_app_lg_2+

  (0)*conv_par_lg_1+(-3.039846)*conv_par_lg_2+

  (-0.9274766)*sys_resp_lg_m+(-0.49927555)*ideo_lev_lg_m+

  (0.11060958)*rep_pot_lg_m+ (-0.34180273)*prot_app_lg_m+

  (-1.8328984)*conv_par_lg_m

Cluster_lg_3<-(-3.6424675)+

  (0)*sys_resp_lg_1+ (-0.61732246)*sys_resp_lg_2+

  (0)*ideo_lev_lg_1+ (-0.23276128)*ideo_lev_lg_2+

  (0)*rep_pot_lg_1+ (3.6818864)*rep_pot_lg_2+

  (0)*prot_app_lg_1+ (3.0608866)*prot_app_lg_2+

  (0)*conv_par_lg_1+ (-1.0034281)*conv_par_lg_2+

  (-0.40726727)*sys_resp_lg_m+(-0.089565904)*ideo_lev_lg_m+

  (1.9387485)*rep_pot_lg_m+(2.5015355)*prot_app_lg_m+

  (-0.81826931)*conv_par_lg_m

# Part 3: Compute odds from logits, as well as modal class and

# probabilities from odds

max_lg<-Cluster_lg_1

if(Cluster_lg_2 > max_lg) {

 max_lg<-Cluster_lg_2

}

if(Cluster_lg_3 > max_lg) {

 max_lg<-Cluster_lg_3

}

Cluster_lg_1<-exp(Cluster_lg_1-max_lg)

Cluster_lg_2<-exp(Cluster_lg_2-max_lg)

Cluster_lg_3<-exp(Cluster_lg_3-max_lg)

max_lg<-Cluster_lg_1

Cluster_lg_modal<-1

if(Cluster_lg_2 > max_lg) {

 max_lg<-Cluster_lg_2; Cluster_lg_modal<-2

}

if(Cluster_lg_3 > max_lg) {

 max_lg<-Cluster_lg_3; Cluster_lg_modal<-3

}

sum_lg<-Cluster_lg_1 + Cluster_lg_2 + Cluster_lg_3

Cluster_lg_1<-Cluster_lg_1/sum_lg

Cluster_lg_2<-Cluster_lg_2/sum_lg

Cluster_lg_3<-Cluster_lg_3/sum_lg

return(list(

 ”Cluster_modal”=Cluster_lg_modal,

 ”Cluster_1”=Cluster_lg_1,

 ”Cluster_2”=Cluster_lg_2,

 ”Cluster_3”=Cluster_lg_3

))

}

## Example of call of scoring function in a loop over records

outdata<-inpdata

for(i in 1:nrow(outdata)){

scoring<-lg_scoring(outdata[i,])

outdata[i,”Cluster_modal”]<-scoring$Cluster_modal

outdata[i,”Cluster_1”]<-scoring$Cluster_1

outdata[i,”Cluster_2”]<-scoring$Cluster_2

outdata[i,”Cluster_3”]<-scoring$Cluster_3

}

As an example, let us take a subject with sys_resp = NA, ideo_lev = 1, rep_pot = 2, prot_app = 2, and conv_par = NA. For this data pattern, in part 1, the dummies sys_resp_lg_m, ideo_lev_lg_1, rep_pot_lg_2, prot_app_lg_2, and conv_par_lg_m are set to 1, and the remaining ones to 0. By summing the weights for which the dummies equal 1, part 2 yields the values 0, −1.8329, and −.8183 for the three linear equations. Transforming these to probabilities in part 3, gives posteriors .1095, .1766, and .7139 for the three classes.

Appendix B: Taylor Approximation of the Normalizing Constants with a Covariate Having a Direct Effect on a Categorical Indicator

As shown in the main text, when a numeric covariate z has a direct effect on a categorical indicator, the normalizing constant of the indicator concerned becomes: $$E_{\mathit{\text{jk}}|z}=\sum _{y_j=1}^{R_j}\mathrm{e}\mathrm{\text{xp}}\left({\alpha }_{y_j}+{\beta }_{y_{j}k}+{\delta }_{y_j}z\right),$$ and will thus depend on the value of z. As a result, the scoring equations will no longer be linear logistic. However, a possible way out is to approximate this term using a Taylor expansion.

For simplicity, assume covariate z is centered, and thus has a mean of 0. The second-order Taylor approximation of $\mathrm{\text{log}}\hspace{0.17em}{E}_{\mathit{\text{jk}}|z}$ at $z=0$ equals: $$\mathrm{\text{log}}\hspace{0.17em}{}_{\mathit{\text{jk}}|z}=\mathrm{\text{log}}\hspace{0.17em}{E}_{\mathit{\text{jk}}|z=0}+\frac{\mathrm{d }\mathrm{\text{log}}{E}_{\mathit{\text{jk}}|z=0}}{\mathrm{d}z}z+\frac{{\mathrm{d}}^2\mathrm{ }\mathrm{\text{log}}{E}_{\mathit{\text{jk}}|z=0}}{{\mathrm{d}}^{2}z}/2\mathrm{ }{z}^2,$$ with $$\frac{\mathrm{d }\mathrm{\text{log}}{E}_{\mathit{\text{jk}}|z=0}}{\mathrm{d}z}=\sum _{y_j=1}^{R_j}P\left(Y_j=y_j\right|X=k,z=0){\delta }_{y_j}$$

and $$\frac{{\mathrm{d}}^2\mathrm{ }\mathrm{\text{log}}{E}_{\mathit{\text{jk}}|z=0}}{{\mathrm{d}}^{2}z}=\sum _{y_j=1}^{R_j}P\left(Y_j=y_j\right|X=k,z=0\left)\right[{\delta }_{y_j}-\frac{\mathrm{d }\mathrm{\text{log}}{E}_{\mathit{\text{jk}}|z=0}}{\mathrm{d}z}]{\delta }_{y_j}.$$

The term $\mathrm{\text{log}}\hspace{0.17em}{E}_{\mathit{\text{jk}}|z=0}$ is subtracted from the intercept of class k, the $\frac{\mathrm{d }\mathrm{\text{log}}{E}_{\mathit{\text{jk}}|z=0}}{\mathrm{d}z}$ are subtracted from the linear terms for z, and the $-\frac{{\mathrm{d}}^2\mathrm{ }\mathrm{\text{log}}{E}_{\mathit{\text{jk}}|z=0}}{{\mathrm{d}}^{2}z}/2$ appear as quadratic terms in the scoring equations.