Skip to contents

This vignette describes the data generation process for a two-profile CULTA model with a covariate. The CULTA model incorporates covariates, latent categorical variables, trait components, state components, and profile-specific means to simulate longitudinal data with latent profile transitions. Data generation is demonstrated at the end of the vignette using the GenCULTA2Profiles function, which integrates all components into a single function call.

Let i{1,,n}i \in \left\{ 1, \ldots, n \right\} denote the index for individuals, let t{0,,m1}t \in \left\{ 0, \ldots, m - 1 \right\} denote the index measurement occasions, let k{1,,p}k \in \left\{ 1, \ldots, p \right\} denote the index items, and let c{0,1}c \in \left\{ 0, 1 \right\} be the index of the two latent profiles (profile 0 and profile 1). Let qq be the trait dimension, q=1q = 1 in this context.

# dimensions
n # number of individuals
#> [1] 10000
m # measurement occasions
#> [1] 6
p # number of items
#> [1] 4
q # common trait dimension
#> [1] 1

Covariate

The covariate is generated from a normal distribution with mean μX\mu_X and σX\sigma_X variance.

# covariate parameters
mu_x 
#> [1] 11.4009
sigma_x
#> [1] 24.67566

Latent Categorical Variables

Latent categorical variables represent profile membership for each individual at each measurement occasion. In a two-profile model, profile membership is influenced by covariates and previous profile status, following a logistic formulation. We distinguish between:

  • Initial profile membership (baseline time point)
  • Profile transitions across subsequent time points

We describe both components below.

Initial Profile Membership

For the first measurement occasion (t=0t = 0), profile membership is determined by the following log-odds for belonging to profile 0 (with profile 1 as the reference category): (ν0+κ0×Covariate0).\begin{equation} \left( \begin{array}{cc} \nu_{0} + \kappa_{0} \times \mathrm{Covariate} & 0 \\ \end{array} \right) . \end{equation} The corresponding probability of belonging to each profile is given by: (exp(ν0+κ0×Covariate)exp(ν0+κ0×Covariate)+11exp(ν0+κ0×Covariate)+1).\begin{equation} \left( \begin{array}{cc} \frac{ \exp \left( \nu_{0} + \kappa_{0} \times \mathrm{Covariate} \right) }{ \exp \left( \nu_{0} + \kappa_{0} \times \mathrm{Covariate} \right) + 1 } & \frac{1}{ \exp \left( \nu_{0} + \kappa_{0} \times \mathrm{Covariate} \right) + 1 } \\ \end{array} \right) . \end{equation} Profile membership at the first occasion is sampled based on these probabilities.

Profile Transitions

For subsequent occasions (t=1,,m1t = 1, \ldots, m - 1), profile transitions depend on the profile at the previous occasion and the covariate. The log-odds for transitioning to profile 0 at time tt are given by: (α0+β00+γ00×Covariate0α0+γ01×Covariate0).\begin{equation} \left( \begin{array}{cc} \alpha_{0} + \beta_{00} + \gamma_{00} \times \mathrm{Covariate} & 0 \\ \alpha_{0} + \gamma_{01} \times \mathrm{Covariate} & 0 \\ \end{array} \right) . \end{equation} The probability of transitioning to each profile is computed as: (exp(α0+β00+γ00×Covariate)exp(α0+β00+γ00×Covariate)+11exp(α0+β00+γ00×Covariate)+1exp(α0+γ01×Covariate)exp(α0+γ01×Covariate)+11exp(α0+γ01×Covariate)+1).\begin{equation} \left( \begin{array}{cc} \frac{ \exp \left( \alpha_{0} + \beta_{00} + \gamma_{00} \times \mathrm{Covariate} \right) }{ \exp \left( \alpha_{0} + \beta_{00} + \gamma_{00} \times \mathrm{Covariate} \right) + 1 } & \frac{1}{ \exp \left( \alpha_{0} + \beta_{00} + \gamma_{00} \times \mathrm{Covariate} \right) + 1 } \\ \frac{ \exp \left( \alpha_{0} + \gamma_{01} \times \mathrm{Covariate} \right) }{ \exp \left( \alpha_{0} + \gamma_{01} \times \mathrm{Covariate} \right) + 1 } & \frac{1}{ \exp \left( \alpha_{0} + \gamma_{01} \times \mathrm{Covariate} \right) + 1 } \\ \end{array} \right) . \end{equation} Profile membership for each subsequent time point is sampled using these transition probabilities, based on the individual’s covariate value and previous profile.

# profile membership and transition parameters
nu_0
#> [1] -3.563
kappa_0
#> [1] 0.122
alpha_0
#> [1] -3.586
beta_00
#> [1] 2.25
gamma_00
#> [1] 0.063
gamma_10
#> [1] 0.094

Trait Components

The trait variate captures between-person differences and is composed of a shared (common) component and item-specific (unique) components. The full decomposition is given by: Traiti=CommonTraitLoading×CommonTraiti+UniqueTraiti.\begin{equation} \mathrm{Trait}_{i} = \mathrm{Common\ Trait\ Loading} \times \mathrm{Common\ Trait}_{i} + \mathrm{Unique\ Trait}_{i} . \end{equation} We describe each component below.

Common Trait

The common trait CommonTraiti\mathrm{Common\ Trait}_{i} represents shared individual differences that influence all items uniformly. It is drawn from a normal distribution with mean μT\mu_T and variance ψT\psi_T: CommonTraiti𝒩(μT,ψT)\begin{equation} \mathrm{Common\ Trait}_{i} \sim \mathcal{N} \left( \mu_T, \psi_T \right) \end{equation}

The influence of the common trait on each item is determined by the p×qp \times q common trait loading,

Unique Traits

The unique trait component UniqueTraitk,i\mathrm{Unique\ Trait}_{k, i} captures item-specific stable differences and is drawn from a multivariate normal distribution: UniqueTraiti𝒩(𝛍p,𝚿p×p)\begin{equation} \mathrm{Unique\ Trait}_{i} \sim \mathcal{N} \left( \boldsymbol{\mu}_p, \boldsymbol{\Psi}_{p \times p} \right) \end{equation}

Combined Trait Variate

The trait variate for item kk and individual ii is obtained by combining the common and unique trait components: Traitk,i=CommonTraitLoadingk×CommonTraiti+UniqueTraitk,i.\begin{equation} \mathrm{Trait}_{k, i} = \mathrm{Common\ Trait\ Loading}_{k} \times \mathrm{Common\ Trait}_{i} + \mathrm{Unique\ Trait}_{k, i} . \end{equation} The common trait component introduces shared variance across items, while the unique trait component allows for item-specific differences not explained by the common trait.

# trait parameters
psi_t
#>      [,1]
#> [1,]  0.1
mu_t
#> [1] 0
psi_p
#>      [,1] [,2] [,3] [,4]
#> [1,]  0.1  0.0  0.0  0.0
#> [2,]  0.0  0.1  0.0  0.0
#> [3,]  0.0  0.0  0.5  0.0
#> [4,]  0.0  0.0  0.0  0.5
mu_p
#> [1] 0 0 0 0
common_trait_loading
#>      [,1]
#> [1,]    1
#> [2,]    1
#> [3,]    1
#> [4,]    1

State Components

The state variate is composed of two parts: a common state shared across items, and unique states specific to each item. The full decomposition is given by: Statek,i,t=CommonStateLoadingk×CommonStatei,t+UniqueStatek,i,t.\begin{equation} \mathrm{State}_{k, i, t} = \mathrm{Common\ State\ Loading}_{k} \times \mathrm{Common\ State}_{i, t} + \mathrm{Unique\ State}_{k, i, t} . \end{equation} We describe each component below.

Common State

The common state CommonStatei,t\mathrm{Common\ State}_{i, t} evolves over time following a first-order autoregressive process: CommonStatei,t=ϕc×CommonStatei,t1+ζi,t.\begin{equation} \mathrm{Common\ State}_{i, t} = \phi_{c} \times \mathrm{Common\ State}_{i, t - 1} + \zeta_{i, t} . \end{equation} The initial common state is drawn from a normal distribution: CommonStatei,0𝒩(0,ψs0).\begin{equation} \mathrm{Common\ State}_{i, 0} \sim \mathcal{N} \left( 0, \psi_{s_{0}} \right) . \end{equation} The innovation term ζi,t\zeta_{i, t} is normally distributed: ζi,t𝒩(0,ψs).\begin{equation} \zeta_{i, t} \sim \mathcal{N} \left( 0, \psi_{s} \right) . \end{equation} The autoregressive parameter ϕc\phi_{c} depends on latent profile membership cc: ϕc=ϕ0+(ϕ1ϕ0)c.\begin{equation} \phi_{c} = \phi_{0} + \left( \phi_{1} - \phi_{0} \right) c . \end{equation} Here, ϕ0\phi_{0} and ϕ1\phi_{1} represent the autoregressive coefficients for profiles coded as 0 and 1, respectively.

Unique State

The unique state UniqueStatek,i,t\mathrm{Unique\ State}_{k, i, t} captures item-specific deviations and is drawn from a multivariate normal distribution: UniqueStatei,t𝒩(0,𝛉)\begin{equation} \mathrm{Unique\ State}_{i, t} \sim \mathcal{N} \left( 0, \boldsymbol{\theta} \right) \end{equation} where 𝛉\boldsymbol{\theta} is the item-level covariance matrix for the unique state component.

Combined State Variate

The state variate for item kk, individual ii, and time tt combines the common and unique state components: Statek,i,t=CommonStateLoadingk×CommonStatei,t+UniqueStatek,i,t\begin{equation} \mathrm{State}_{k, i, t} = \mathrm{Common\ State\ Loading}_{k} \times \mathrm{Common\ State}_{i, t} + \mathrm{Unique\ State}_{k, i, t} \end{equation} The common state loading parameter CommonStateLoadingk\mathrm{Common\ State\ Loading}_{k} controls the influence of the shared state on each item.

# state parameters
common_state_loading
#>      [,1]
#> [1,]    1
#> [2,]    1
#> [3,]    1
#> [4,]    1
phi_0
#> [1] 0
phi_1
#> [1] 0.311
psi_s0
#> [1] 1
psi_s
#> [1] 0.25
theta
#>      [,1] [,2] [,3] [,4]
#> [1,] 0.15 0.00 0.00 0.00
#> [2,] 0.00 0.15 0.00 0.00
#> [3,] 0.00 0.00 0.15 0.00
#> [4,] 0.00 0.00 0.00 0.15

Trait-State Plus Profile Specific Means

The observed variable is given by Yk,i,t=μk,c+Traitk,i+Statek,i,t\begin{equation} Y_{k, i, t} = \mu_{k, c} + \mathrm{Trait}_{k, i} + \mathrm{State}_{k, i, t} \end{equation} where μk,c\mu_{k, c} is the profile specific mean, while Traitk,i\mathrm{Trait}_{k, i} and Statek,i,t\mathrm{State}_{k, i, t} correspond to the trait and state components of the model.

# profile-specific means
mu_profile
#>       [,1]   [,2]
#> [1,] 2.253 -0.278
#> [2,] 1.493 -0.165
#> [3,] 1.574 -0.199
#> [4,] 1.117 -0.148

Data Generation

The GenCULTA2Profiles function puts together all the components described above to generate data using a single function call.

# complete list of R function arguments

# random seed for reproducibility
set.seed(42)

# dimensions
n # number of individuals
#> [1] 10000
m # measurement occasions
#> [1] 6
p # number of items
#> [1] 4
q # common trait dimension
#> [1] 1

# covariate parameters
mu_x 
#> [1] 11.4009
sigma_x
#> [1] 24.67566

# profile membership and transition parameters
nu_0
#> [1] -3.563
kappa_0
#> [1] 0.122
alpha_0
#> [1] -3.586
beta_00
#> [1] 2.25
gamma_00
#> [1] 0.063
gamma_10
#> [1] 0.094

# trait parameters
psi_t
#>      [,1]
#> [1,]  0.1
mu_t
#> [1] 0
psi_p
#>      [,1] [,2] [,3] [,4]
#> [1,]  0.1  0.0  0.0  0.0
#> [2,]  0.0  0.1  0.0  0.0
#> [3,]  0.0  0.0  0.5  0.0
#> [4,]  0.0  0.0  0.0  0.5
mu_p
#> [1] 0 0 0 0
common_trait_loading
#>      [,1]
#> [1,]    1
#> [2,]    1
#> [3,]    1
#> [4,]    1

# state parameters
common_state_loading
#>      [,1]
#> [1,]    1
#> [2,]    1
#> [3,]    1
#> [4,]    1
phi_0
#> [1] 0
phi_1
#> [1] 0.311
psi_s0
#> [1] 1
psi_s
#> [1] 0.25
theta
#>      [,1] [,2] [,3] [,4]
#> [1,] 0.15 0.00 0.00 0.00
#> [2,] 0.00 0.15 0.00 0.00
#> [3,] 0.00 0.00 0.15 0.00
#> [4,] 0.00 0.00 0.00 0.15

# profile-specific means
mu_profile
#>       [,1]   [,2]
#> [1,] 2.253 -0.278
#> [2,] 1.493 -0.165
#> [3,] 1.574 -0.199
#> [4,] 1.117 -0.148
data <- GenCULTA2Profiles(
  n = n,
  m = m,
  mu_x = mu_x,
  sigma_x = sigma_x,
  nu_0 = nu_0,
  kappa_0 = kappa_0,
  alpha_0 = alpha_0,
  beta_00 = beta_00,
  gamma_00 = gamma_00,
  gamma_10 = gamma_10,
  mu_t = mu_t,
  psi_t = psi_t,
  mu_p = mu_p,
  psi_p = psi_p,
  common_trait_loading = common_trait_loading,
  common_state_loading = common_state_loading,
  phi_0 = phi_0,
  phi_1 = phi_1,
  psi_s0 = psi_s0,
  psi_s = psi_s,
  theta = theta, 
  mu_profile = mu_profile
)

Model Fitting

The FitCULTA2Profiles function fits the model using Mplus. Note: This function requires that Mplus is already installed on the system. To speed up model fitting, consider using the ncores argument to leverage multiple cores.

fit <- FitCULTA2Profiles(data = data)
summary(fit)
#>               est     se        z      p    2.5%   97.5%
#> mu_10      2.2690 0.0117 193.2129 0.0000  2.2460  2.2920
#> mu_20      1.5140 0.0112 134.9006 0.0000  1.4920  1.5360
#> mu_30      1.6058 0.0128 125.1805 0.0000  1.5806  1.6309
#> mu_40      1.1275 0.0126  89.7203 0.0000  1.1029  1.1522
#> lambda_t2  1.0505 0.0381  27.5961 0.0000  0.9759  1.1251
#> lambda_s2  1.0063 0.0048 209.3670 0.0000  0.9969  1.0157
#> lambda_t3  1.0209 0.0481  21.2367 0.0000  0.9267  1.1151
#> lambda_s3  1.0041 0.0047 214.8586 0.0000  0.9949  1.0133
#> lambda_t4  1.0261 0.0474  21.6425 0.0000  0.9332  1.1190
#> lambda_s4  1.0121 0.0048 209.9015 0.0000  1.0027  1.0216
#> theta_11   0.1508 0.0014 110.9646 0.0000  0.1481  0.1535
#> theta_22   0.1496 0.0013 117.9673 0.0000  0.1471  0.1521
#> theta_33   0.1484 0.0012 118.8785 0.0000  0.1460  0.1509
#> theta_44   0.1493 0.0013 114.0088 0.0000  0.1468  0.1519
#> phi_0     -0.0140 0.0170  -0.8248 0.4095 -0.0473  0.0193
#> psi_t      0.0896 0.0052  17.0664 0.0000  0.0793  0.0999
#> psi_p_11   0.1004 0.0032  30.9448 0.0000  0.0940  0.1067
#> psi_p_22   0.1004 0.0035  28.7053 0.0000  0.0935  0.1073
#> psi_p_33   0.4963 0.0081  61.0379 0.0000  0.4804  0.5122
#> psi_p_44   0.4938 0.0082  60.5071 0.0000  0.4778  0.5098
#> psi_s0     0.9992 0.0189  52.9458 0.0000  0.9622  1.0362
#> psi_s      0.2481 0.0028  87.6672 0.0000  0.2425  0.2536
#> mu_11     -0.2779 0.0058 -47.6439 0.0000 -0.2893 -0.2665
#> mu_21     -0.1652 0.0059 -28.2113 0.0000 -0.1767 -0.1538
#> mu_31     -0.1830 0.0086 -21.3750 0.0000 -0.1998 -0.1662
#> mu_41     -0.1516 0.0085 -17.7637 0.0000 -0.1683 -0.1349
#> phi_1      0.3169 0.0059  53.6573 0.0000  0.3053  0.3285
#> nu_0      -3.6911 0.1177 -31.3505 0.0000 -3.9218 -3.4603
#> alpha_0   -3.5634 0.0551 -64.7026 0.0000 -3.6713 -3.4555
#> kappa_0    0.1274 0.0080  15.8995 0.0000  0.1117  0.1431
#> beta_00    2.1226 0.1132  18.7586 0.0000  1.9008  2.3443
#> gamma_00   0.0677 0.0064  10.5087 0.0000  0.0550  0.0803
#> gamma_10   0.0942 0.0040  23.5322 0.0000  0.0864  0.1021

Parameter Recovery

Parameter recovery was assessed by calculating the differences between the population values and the estimated parameters.

Parameter Recovery
Parameter Estimate Difference
mu_10 2.253 2.2689973 -0.0159973
mu_20 1.493 1.5140045 -0.0210045
mu_30 1.574 1.6057853 -0.0317853
mu_40 1.117 1.1275251 -0.0105251
lambda_t2 1.000 1.0504634 -0.0504634
lambda_s2 1.000 1.0062885 -0.0062885
lambda_t3 1.000 1.0209131 -0.0209131
lambda_s3 1.000 1.0040943 -0.0040943
lambda_t4 1.000 1.0260833 -0.0260833
lambda_s4 1.000 1.0121351 -0.0121351
theta_11 0.150 0.1508128 -0.0008128
theta_22 0.150 0.1496218 0.0003782
theta_33 0.150 0.1484347 0.0015653
theta_44 0.150 0.1493491 0.0006509
phi_0 0.000 -0.0140050 0.0140050
psi_t 0.100 0.0895956 0.0104044
psi_p_11 0.100 0.1003680 -0.0003680
psi_p_22 0.100 0.1004050 -0.0004050
psi_p_33 0.500 0.4962964 0.0037036
psi_p_44 0.500 0.4937720 0.0062280
psi_s0 1.000 0.9991861 0.0008139
psi_s 0.250 0.2480841 0.0019159
mu_11 -0.278 -0.2779016 -0.0000984
mu_21 -0.165 -0.1652498 0.0002498
mu_31 -0.199 -0.1829942 -0.0160058
mu_41 -0.148 -0.1515997 0.0035997
phi_1 0.311 0.3169217 -0.0059217
nu_0 -3.563 -3.6910557 0.1280557
alpha_0 -3.586 -3.5634025 -0.0225975
kappa_0 0.122 0.1273813 -0.0053813
beta_00 2.250 2.1225611 0.1274389
gamma_00 0.063 0.0676685 -0.0046685
gamma_10 0.094 0.0942463 -0.0002463

Methods

The FitCULTA2Profiles function has a number of methods including the following:

  • print
  • summary
  • converged
  • confint
  • coef
  • vcov
  • logLik
  • AIC
  • BIC
  • entropy
  • anova

Mplus Script Used

TITLE:
  2-Profile CULTA with Covariate;

DATA:
  FILE = culta_data.dat;

VARIABLE:
  NAMES =
    id
    x
    y1t0
    y2t0
    y3t0
    y4t0
    y1t1
    y2t1
    y3t1
    y4t1
    y1t2
    y2t2
    y3t2
    y4t2
    y1t3
    y2t3
    y3t3
    y4t3
    y1t4
    y2t4
    y3t4
    y4t4
    y1t5
    y2t5
    y3t5
    y4t5
  ;
  USEVARIABLES =
    x
    y1t0
    y2t0
    y3t0
    y4t0
    y1t1
    y2t1
    y3t1
    y4t1
    y1t2
    y2t2
    y3t2
    y4t2
    y1t3
    y2t3
    y3t3
    y4t3
    y1t4
    y2t4
    y3t4
    y4t4
    y1t5
    y2t5
    y3t5
    y4t5
  ;
  IDVARIABLE = id;
  CLASSES =
    c0(2)
    c1(2)
    c2(2)
    c3(2)
    c4(2)
    c5(2)
  ;

ANALYSIS:
  TYPE = MIXTURE;
  STARTS = 20 4;
  STSCALE = 5;
  STITERATIONS = 10;
  PROCESS = 1;
  MODEL = NOCOV;

MODEL:
  %OVERALL%
    ! common trait ------------------------------------------------------

    !! factor loadings
    !!! t = 0
    t BY y1t0@1;
    t BY y2t0 (lambdat2);
    t BY y3t0 (lambdat3);
    t BY y4t0 (lambdat4);
    !!! t = 1
    t BY y1t1@1;
    t BY y2t1 (lambdat2);
    t BY y3t1 (lambdat3);
    t BY y4t1 (lambdat4);
    !!! t = 2
    t BY y1t2@1;
    t BY y2t2 (lambdat2);
    t BY y3t2 (lambdat3);
    t BY y4t2 (lambdat4);
    !!! t = 3
    t BY y1t3@1;
    t BY y2t3 (lambdat2);
    t BY y3t3 (lambdat3);
    t BY y4t3 (lambdat4);
    !!! t = 4
    t BY y1t4@1;
    t BY y2t4 (lambdat2);
    t BY y3t4 (lambdat3);
    t BY y4t4 (lambdat4);
    !!! t = 5
    t BY y1t5@1;
    t BY y2t5 (lambdat2);
    t BY y3t5 (lambdat3);
    t BY y4t5 (lambdat4);

    !! latent mean
    [ t@0 ];

    !! latent variance
    t (psit);

    ! unique traits -----------------------------------------------------

    !! factor loadings
    !!! k = 1
    u1 BY y1t0@1;
    u1 BY y1t1@1;
    u1 BY y1t2@1;
    u1 BY y1t3@1;
    u1 BY y1t4@1;
    u1 BY y1t5@1;
    !!! k = 2
    u2 BY y2t0@1;
    u2 BY y2t1@1;
    u2 BY y2t2@1;
    u2 BY y2t3@1;
    u2 BY y2t4@1;
    u2 BY y2t5@1;
    !!! k = 3
    u3 BY y3t0@1;
    u3 BY y3t1@1;
    u3 BY y3t2@1;
    u3 BY y3t3@1;
    u3 BY y3t4@1;
    u3 BY y3t5@1;
    !!! k = 4
    u4 BY y4t0@1;
    u4 BY y4t1@1;
    u4 BY y4t2@1;
    u4 BY y4t3@1;
    u4 BY y4t4@1;
    u4 BY y4t5@1;

    !! latent means
    [ u1@0 ];
    [ u2@0 ];
    [ u3@0 ];
    [ u4@0 ];

    !! latent variances
    u1 (psip1);
    u2 (psip2);
    u3 (psip3);
    u4 (psip4);

    ! common states -----------------------------------------------------

    !! factor loadings
    !!! t = 0
    s0 BY y1t0@1;
    s0 BY y2t0 (lambdas2);
    s0 BY y3t0 (lambdas3);
    s0 BY y4t0 (lambdas4);
    !!! t = 1
    s1 BY y1t1@1;
    s1 BY y2t1 (lambdas2);
    s1 BY y3t1 (lambdas3);
    s1 BY y4t1 (lambdas4);
    !!! t = 2
    s2 BY y1t2@1;
    s2 BY y2t2 (lambdas2);
    s2 BY y3t2 (lambdas3);
    s2 BY y4t2 (lambdas4);
    !!! t = 3
    s3 BY y1t3@1;
    s3 BY y2t3 (lambdas2);
    s3 BY y3t3 (lambdas3);
    s3 BY y4t3 (lambdas4);
    !!! t = 4
    s4 BY y1t4@1;
    s4 BY y2t4 (lambdas2);
    s4 BY y3t4 (lambdas3);
    s4 BY y4t4 (lambdas4);
    !!! t = 5
    s5 BY y1t5@1;
    s5 BY y2t5 (lambdas2);
    s5 BY y3t5 (lambdas3);
    s5 BY y4t5 (lambdas4);

    !! latent means
    [ s0@0 ];
    [ s1@0 ];
    [ s2@0 ];
    [ s3@0 ];
    [ s4@0 ];
    [ s5@0 ];

    !! latent variance of s0
    s0 (psis0);

    !! variance of the process noise
    s1 (psis);
    s2 (psis);
    s3 (psis);
    s4 (psis);
    s5 (psis);

    ! unique states -----------------------------------------------------

    !! variances
    !!! t = 0
    y1t0 (theta11);
    y2t0 (theta22);
    y3t0 (theta33);
    y4t0 (theta44);
    !!! t = 1
    y1t1 (theta11);
    y2t1 (theta22);
    y3t1 (theta33);
    y4t1 (theta44);
    !!! t = 2
    y1t2 (theta11);
    y2t2 (theta22);
    y3t2 (theta33);
    y4t2 (theta44);
    !!! t = 3
    y1t3 (theta11);
    y2t3 (theta22);
    y3t3 (theta33);
    y4t3 (theta44);
    !!! t = 4
    y1t4 (theta11);
    y2t4 (theta22);
    y3t4 (theta33);
    y4t4 (theta44);
    !!! t = 5
    y1t5 (theta11);
    y2t5 (theta22);
    y3t5 (theta33);
    y4t5 (theta44);

    ! constrained intercepts --------------------------------------------

    !! t = 0
    [ y1t0@0 ];
    [ y2t0@0 ];
    [ y3t0@0 ];
    [ y4t0@0 ];
    !! t = 1
    [ y1t1@0 ];
    [ y2t1@0 ];
    [ y3t1@0 ];
    [ y4t1@0 ];
    !! t = 2
    [ y1t2@0 ];
    [ y2t2@0 ];
    [ y3t2@0 ];
    [ y4t2@0 ];
    !! t = 3
    [ y1t3@0 ];
    [ y2t3@0 ];
    [ y3t3@0 ];
    [ y4t3@0 ];
    !! t = 4
    [ y1t4@0 ];
    [ y2t4@0 ];
    [ y3t4@0 ];
    [ y4t4@0 ];
    !! t = 5
    [ y1t5@0 ];
    [ y2t5@0 ];
    [ y3t5@0 ];
    [ y4t5@0 ];

    ! lta ---------------------------------------------------------------

    !! initial profile membership
    [ c0#1 ] (nu0);
    c0#1 ON x (kappa0);

    !! profile transitions
    [ c1#1 ] (alpha0);
    [ c2#1 ] (alpha0);
    [ c3#1 ] (alpha0);
    [ c4#1 ] (alpha0);
    [ c5#1 ] (alpha0);
    c1#1 ON c0#1 (beta00);
    c2#1 ON c1#1 (beta00);
    c3#1 ON c2#1 (beta00);
    c4#1 ON c3#1 (beta00);
    c5#1 ON c4#1 (beta00);

MODEL c0:
  %c0#1%
    ! profile specific means
    [ y1t0 ] (mu10);
    [ y2t0 ] (mu20);
    [ y3t0 ] (mu30);
    [ y4t0 ] (mu40);

    ! covariate
    c1 ON x (gamma00);

  %c0#2%
    ! profile specific means
    [ y1t0 ] (mu11);
    [ y2t0 ] (mu21);
    [ y3t0 ] (mu31);
    [ y4t0 ] (mu41);

    ! covariate
    c1 ON x (gamma10);

MODEL c1:
  %c1#1%
    ! profile specific means
    [ y1t1 ] (mu10);
    [ y2t1 ] (mu20);
    [ y3t1 ] (mu30);
    [ y4t1 ] (mu40);

    ! covariate
    c2 ON x (gamma00);

    ! inertia
    s1 ON s0 (phi0);

  %c1#2%
    ! profile specific means
    [ y1t1 ] (mu11);
    [ y2t1 ] (mu21);
    [ y3t1 ] (mu31);
    [ y4t1 ] (mu41);

    ! covariate
    c2 ON x (gamma10);

    ! inertia
    s1 ON s0 (phi1);

MODEL c2:
  %c2#1%
    ! profile specific means
    [ y1t2 ] (mu10);
    [ y2t2 ] (mu20);
    [ y3t2 ] (mu30);
    [ y4t2 ] (mu40);

    ! covariate
    c3 ON x (gamma00);

    ! inertia
    s2 ON s1 (phi0);

  %c2#2%
    ! profile specific means
    [ y1t2 ] (mu11);
    [ y2t2 ] (mu21);
    [ y3t2 ] (mu31);
    [ y4t2 ] (mu41);

    ! covariate
    c3 ON x (gamma10);

    ! inertia
    s2 ON s1 (phi1);

MODEL c3:
  %c3#1%
    ! profile specific means
    [ y1t3 ] (mu10);
    [ y2t3 ] (mu20);
    [ y3t3 ] (mu30);
    [ y4t3 ] (mu40);

    ! covariate
    c4 ON x (gamma00);

    ! inertia
    s3 ON s2 (phi0);

  %c3#2%
    ! profile specific means
    [ y1t3 ] (mu11);
    [ y2t3 ] (mu21);
    [ y3t3 ] (mu31);
    [ y4t3 ] (mu41);

    ! covariate
    c4 ON x (gamma10);

    ! inertia
    s3 ON s2 (phi1);

MODEL c4:
  %c4#1%
    ! profile specific means
    [ y1t4 ] (mu10);
    [ y2t4 ] (mu20);
    [ y3t4 ] (mu30);
    [ y4t4 ] (mu40);

    ! covariate
    c5 ON x (gamma00);

    ! inertia
    s4 ON s3 (phi0);

  %c4#2%
    ! profile specific means
    [ y1t4 ] (mu11);
    [ y2t4 ] (mu21);
    [ y3t4 ] (mu31);
    [ y4t4 ] (mu41);

    ! covariate
    c5 ON x (gamma10);

    ! inertia
    s4 ON s3 (phi1);

MODEL c5:
  %c5#1%
    ! profile specific means
    [ y1t5 ] (mu10);
    [ y2t5 ] (mu20);
    [ y3t5 ] (mu30);
    [ y4t5 ] (mu40);

    ! inertia
    s5 ON s4 (phi0);

  %c5#2%
    ! profile specific means
    [ y1t5 ] (mu11);
    [ y2t5 ] (mu21);
    [ y3t5 ] (mu31);
    [ y4t5 ] (mu41);

    ! inertia
    s5 ON s4 (phi1);


MODEL CONSTRAINT:
  ! means for the first profile are higher than the second
  mu10 > mu11;
  mu20 > mu21;
  mu30 > mu31;
  mu40 > mu41;

  ! make sure variances are greater than zero
  psit > 0;
  psip1 > 0;
  psip2 > 0;
  psip3 > 0;
  psip4 > 0;
  psis0 > 0;
  psis > 0;
  theta11 > 0;
  theta22 > 0;
  theta33 > 0;
  theta44 > 0;

OUTPUT:
  TECH1
  TECH3
  TECH4
  TECH7
  TECH8
  TECH12
  TECH13
  TECH15
  ENTROPY;

SAVEDATA:
  ESTIMATES = culta_34rXapIX7srOsFvdz2Ab_estimates.dat;
  RESULTS = culta_34rXapIX7srOsFvdz2Ab_results.dat;
  TECH3 = culta_34rXapIX7srOsFvdz2Ab_tech3.dat;
  TECH4 = culta_34rXapIX7srOsFvdz2Ab_tech4.dat;
  FILE = culta_34rXapIX7srOsFvdz2Ab_cprobs.dat;
  SAVE = CPROBABILITIES;