Longitudinal CLCM with Multinomial Latent Regression

The stratified transition matrix is a transition matrix estimated separately for different groups. This is essentially a type of latent regression. Here we regress the latent classes at timepoint 2 onto the observed Group, which allows us to estimate a stratified transition matrix - this is the output of interest when, for example, we look at the differences in treatment arms (observed groups) in the proportion of subjects moving from the “Sick” latent class to the “Not Sick” latent class.

library(CLCM)
library(nnet)

Specify details of the data. Here we have 2,000 subjects belonging to 4 different latent classes responding to 16 items at 2 timepoints.

    N <- 2000
     number.timepoints <- 2
     item.type <- rep( c('Ordinal', 'Nominal', 'Poisson', 'Neg_Binom', 'ZINB', 'ZIP', 'Normal', 'Beta'), 2)
     sim.categories.j <- rep( c(4, 4, 30, 30, 30, 30, NA, NA) , 2)
     J <- length(item.type)
     item.names <- paste0('Item_', 1:J)
     K= 2
     #alpha <- pattern(K)
     Q <- matrix(c(rep(c(1,0), J/2), rep(c(0, 1), J/2)), nrow = J, ncol = K, byrow = TRUE)
     #eta=alpha %*% t(Q)

Generate the subject posteriors via (multinomial) latent regression. The latent class assignments will differ across the two observed groups (Group).

     set.seed(09092026)
     number.groups <- 2
     dat <- data.frame(
                  'USUBJID' = rep(paste0('Subject_', formatC(1:N, width = 4, flag = '0')), length.out= N*number.timepoints),
                  'Group' = rep(paste0('Group_', 1:number.groups), length.out = N*number.timepoints),
                  'Time' = rep(paste0('Time_', 1:number.timepoints), each = N),
                  stringsAsFactors= FALSE)

      # Design Matrix
      XX <- model.matrix( ~ Group*Time, data = dat) 
      # Beta
      Beta <- matrix(0, nrow = ncol(XX), ncol = 2^K - 1, dimnames=list(colnames(XX), rep('LC', 2^K -1)))
      Beta['(Intercept)', ] <- c(0.2, 0.8, 0.4)
      Beta['GroupGroup_2', ] <- 0
      Beta['TimeTime_2', ] <- 0
      Beta['GroupGroup_2:TimeTime_2', ] <- c(0.2, -1.0, -1.6)
      Beta
#>                          LC   LC   LC
#> (Intercept)             0.2  0.8  0.4
#> GroupGroup_2            0.0  0.0  0.0
#> TimeTime_2              0.0  0.0  0.0
#> GroupGroup_2:TimeTime_2 0.2 -1.0 -1.6
      
      # Matrix multiply:
      XB <- XX %*% Beta
      p <- exp(XB)/(1 + apply(exp(XB), 1, sum))
      p <-  cbind(1 - rowSums(p), p)
      lca <-  vector()
      for(i in 1:nrow(p)){
          lca <- c(lca, 
           sample(x = c(1:(2^K)), size = 1, prob = p[i, ])  
                  )
        } #end loop to sample lca

At this point we have used the observed group membership (Group) to generate the (true) posterior distributions/latent class assignments for both timepoints. Let’s check the latent class assignments vs Observed group assignment. Again, we haven’t yet moved into latent class models yet.

      # Check the LCA 
      prop.table(table(lca))
#> lca
#>       1       2       3       4 
#> 0.19725 0.25150 0.33950 0.21175
      # Check LCA across observed Groups:
      prop.table(xtabs( ~ lca + dat$Group), margin = 1)
#>    dat$Group
#> lca   Group_1   Group_2
#>   1 0.4410646 0.5589354
#>   2 0.3916501 0.6083499
#>   3 0.5405007 0.4594993
#>   4 0.6186541 0.3813459

Before we move on with the simulation, let’s just do a quick check on whether a simple multinomial logistic regression can recover this. Regress latent class assignment onto observed groups:

      mod <- multinom(as.factor(lca) ~ dat$Group*dat$Time)
#> # weights:  20 (12 variable)
#> initial  value 5545.177444 
#> iter  10 value 5356.080253
#> final  value 5292.792148 
#> converged
      coef(mod)
#>   (Intercept) dat$GroupGroup_2 dat$TimeTime_2 dat$GroupGroup_2:dat$TimeTime_2
#> 2  0.04103973       0.27118184      0.1540957                      -0.1298448
#> 3  0.81158309       0.04191958     -0.1295583                      -0.8533770
#> 4  0.52745497      -0.12405571     -0.2411555                      -1.3903230
      t(Beta)
#>    (Intercept) GroupGroup_2 TimeTime_2 GroupGroup_2:TimeTime_2
#> LC         0.2            0          0                     0.2
#> LC         0.8            0          0                    -1.0
#> LC         0.4            0          0                    -1.6

Recovery seems fine. Let’s create the true posterior distributions and pass to the simulation function.

      # Create true posterior distributions, pass 
      post.true <- matrix(0, nrow = nrow(dat), ncol = 2^K)
      post.true[ cbind(1:nrow(dat), lca) ] <- 1 
      
      # Simulate Data
        set.seed(03102021)
      sim.dat <- simulate_clcm(N = N, number.timepoints = number.timepoints, 
                        Q = Q, 
                        item.names = item.names, 
                        item.type = item.type, 
                        categories.j = sim.categories.j, 
                        post = post.true)

Merge the dataframe with the Observed Group membership with the simulated dataframe. You need the Group variable with the item responses to pass to the model estimation function. If you don’t, the estimation routine won’t have a Group variable to compute the multinomial latent regression with.

      dat.cov <- merge(x = dat, y = sim.dat$dat, by = c('USUBJID', 'Time'))

Specify the latent regression formula - regress latent classes onto Group variable in model estimation by passing lat.reg = list('Time_1' = NULL, 'Time_2' = 'Group'). Note: if we had additional covariates, we would pass the corresponding regression formula, e.g. 'Group + Sex + Age'. This will be passed to as.formula and then passed to a function that fits a regression.

      mod <- clcm(dat = dat.cov, 
                 item.type = item.type, 
                 item.names = item.names, 
                 Q = Q, max.diff = 0.001, 
                 lat.reg = list('Time_1' = NULL, 'Time_2' = 'Group') )    
#> iteration: 1  max diff in item parameter estimates: 8.495015
#> Warning in log(p): NaNs produced
#> iteration: 2  max diff in item parameter estimates: 7.474376
#> iteration: 3  max diff in item parameter estimates: 0.06421
#> iteration: 4  max diff in item parameter estimates: 0.00069


      Beta # Generating parameter
#>                          LC   LC   LC
#> (Intercept)             0.2  0.8  0.4
#> GroupGroup_2            0.0  0.0  0.0
#> TimeTime_2              0.0  0.0  0.0
#> GroupGroup_2:TimeTime_2 0.2 -1.0 -1.6
      mod$lat.reg.param  # estimate
#>                   [,1]       [,2]       [,3]
#> (Intercept)  0.1950155  0.6820203  0.2863678
#> GroupGroup_2 0.1416452 -0.8112243 -1.5142441

Evaluate the transition matrix, first not stratified across the Group variable, and compare that estimated transition matrix to the true transition matrix.


      transition_matrix_clcm(mod = mod, stratification = FALSE)
#>            post_LC_00 post_LC_10 post_LC_01 post_LC_11
#> post_LC_00  0.2616113  0.3093455  0.2800013  0.1490420
#> post_LC_10  0.2259186  0.2995356  0.3222949  0.1522509
#> post_LC_01  0.2336588  0.3083712  0.2853638  0.1726062
#> post_LC_11  0.2098148  0.3053345  0.3244986  0.1603521

      # Compare to Generating parameters
      post.true1 <- post.true[dat$Time == 'Time_1', ]
      post.true2 <- post.true[dat$Time == 'Time_2', ]
      t(post.true1) %*% post.true2/matrix(colSums(post.true1), nrow = 2^K, ncol = 2^K, byrow = FALSE)
#>           [,1]      [,2]      [,3]      [,4]
#> [1,] 0.2621951 0.3079268 0.2804878 0.1493902
#> [2,] 0.2258883 0.2994924 0.3223350 0.1522843
#> [3,] 0.2334218 0.3090186 0.2851459 0.1724138
#> [4,] 0.2099237 0.3053435 0.3244275 0.1603053

Appear to be very similar, very good recovery of the true generating parameters. Now let’s check the estimated transition matrix, stratified on the Group variable:


  # Estimated:
  transition_matrix_clcm(mod = mod, stratification = TRUE, covariate = 'Group')
#> $Group_1
#>            post_LC_00 post_LC_10 post_LC_01 post_LC_11
#> post_LC_00  0.2037142  0.2213342  0.3594356  0.2155160
#> post_LC_10  0.1781415  0.1954175  0.3678050  0.2586361
#> post_LC_01  0.1755545  0.2287569  0.3323286  0.2633600
#> post_LC_11  0.1766021  0.2226087  0.3852151  0.2155741
#> 
#> $Group_2
#>            post_LC_00 post_LC_10 post_LC_01 post_LC_11
#> post_LC_00  0.3214143  0.4002540  0.1979520 0.08037972
#> post_LC_10  0.2637161  0.3819060  0.2862907 0.06808716
#> post_LC_01  0.2915595  0.3877066  0.2385636 0.08217030
#> post_LC_11  0.2488048  0.4024503  0.2532206 0.09552431

  # Compare to Generating parameters:
  post.true11 <- post.true[dat$Time == 'Time_1' & dat$Group == 'Group_1', ]
  post.true12 <- post.true[dat$Time == 'Time_1' & dat$Group == 'Group_2', ]
  post.true21 <- post.true[dat$Time == 'Time_2' & dat$Group == 'Group_1', ]
  post.true22 <- post.true[dat$Time == 'Time_2' & dat$Group == 'Group_2', ]
  t(post.true11) %*% post.true21/matrix(colSums(post.true11), nrow = 2^K, ncol = 2^K, byrow = FALSE)
#>           [,1]      [,2]      [,3]      [,4]
#> [1,] 0.2035928 0.2215569 0.3592814 0.2155689
#> [2,] 0.1781609 0.1954023 0.3678161 0.2586207
#> [3,] 0.1755319 0.2287234 0.3324468 0.2632979
#> [4,] 0.1766784 0.2226148 0.3851590 0.2155477
  t(post.true12) %*% post.true22/matrix(colSums(post.true12), nrow = 2^K, ncol = 2^K, byrow = FALSE)
#>           [,1]      [,2]      [,3]       [,4]
#> [1,] 0.3229814 0.3975155 0.1987578 0.08074534
#> [2,] 0.2636364 0.3818182 0.2863636 0.06818182
#> [3,] 0.2910053 0.3888889 0.2380952 0.08201058
#> [4,] 0.2489627 0.4024896 0.2531120 0.09543568

Next, compare true classification (latent class assignment, lca) with estimates:


      lca.hat <- mod$dat$lca 
      lca.true <- mod$dat$true_lca      
      table(lca.true == lca.hat)
#> 
#> FALSE  TRUE 
#>     1  3999
      prop.table(table(lca.true == lca.hat)) 
#> 
#>   FALSE    TRUE 
#> 0.00025 0.99975
      xtabs( ~ lca.true + lca.hat)
#>         lca.hat
#> lca.true    1    2    3    4
#>        1  789    0    0    0
#>        2    0 1006    0    0
#>        3    1    0 1357    0
#>        4    0    0    0  847