Vig_7.RmdThe stratified transition matrix is a transition matrix estimated separately for different groups. This is essentially for 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.
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 = T)
#eta=alpha %*% t(Q)Generate the subject posteriors via (multinomial) latent regression. The latent class assignments will differ across the two observed groups (Group).
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=F)
# 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.19525 0.27525 0.32350 0.20600
# Check LCA across observed Groups:
prop.table(xtabs( ~ lca + dat$Group), margin = 1)
#> dat$Group
#> lca Group_1 Group_2
#> 1 0.4263764 0.5736236
#> 2 0.3987284 0.6012716
#> 3 0.5625966 0.4374034
#> 4 0.6067961 0.3932039Before 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 <- nnet:::multinom(as.factor(lca) ~ dat$Group*dat$Time)
#> # weights: 20 (12 variable)
#> initial value 5545.177444
#> iter 10 value 5326.348405
#> final value 5271.634859
#> converged
coef(mod)
#> (Intercept) dat$GroupGroup_2 dat$TimeTime_2 dat$GroupGroup_2:dat$TimeTime_2
#> 2 0.3179012 -0.075081667 -0.08304164 0.3120331
#> 3 0.8251576 -0.083503172 -0.08599184 -0.9469774
#> 4 0.3910577 -0.003300142 0.02999722 -1.7651105
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.6Recovery 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 esimtation routine won’t have a Group variable to compute the multinomial latent regression with.
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.698974
#> iteration: 2 max diff in item parameter estimates: 7.626988
#> iteration: 3 max diff in item parameter estimates: 0.058177
#> iteration: 4 max diff in item parameter estimates: 5e-04Evaluate the Latent Regression Parameters:
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.2364976 0.7419877 0.4222014
#> GroupGroup_2 0.2362957 -1.0331095 -1.7749010Evaluate 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 = F)
#> post_LC_00 post_LC_10 post_LC_01 post_LC_11
#> post_LC_00 0.2214193 0.3049894 0.3000771 0.1735143
#> post_LC_10 0.2079710 0.3100054 0.2987354 0.1832881
#> post_LC_01 0.2344934 0.3460567 0.2673978 0.1520522
#> post_LC_11 0.2219060 0.3398452 0.2752494 0.1629994
# 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 = F)
#> [,1] [,2] [,3] [,4]
#> [1,] 0.2215569 0.3053892 0.2994012 0.1736527
#> [2,] 0.2081448 0.3099548 0.2986425 0.1832579
#> [3,] 0.2352941 0.3461012 0.2667579 0.1518468
#> [4,] 0.2210953 0.3387424 0.2758621 0.1643002Appear 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 = T, covariate = 'Group')
#> $Group_1
#> post_LC_00 post_LC_10 post_LC_01 post_LC_11
#> post_LC_00 0.1470123 0.1714995 0.3753373 0.3061509
#> post_LC_10 0.1650507 0.1783732 0.3975230 0.2590531
#> post_LC_01 0.1892237 0.2418851 0.3265645 0.2423268
#> post_LC_11 0.1593103 0.2370078 0.3514829 0.2521990
#>
#> $Group_2
#> post_LC_00 post_LC_10 post_LC_01 post_LC_11
#> post_LC_00 0.2924965 0.4325057 0.2281847 0.04681308
#> post_LC_10 0.2520523 0.4451984 0.1972756 0.10547373
#> post_LC_01 0.2814041 0.4540045 0.2060863 0.05850512
#> post_LC_11 0.2817046 0.4380874 0.2024223 0.07778567
# 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 = F)
#> [,1] [,2] [,3] [,4]
#> [1,] 0.1472393 0.1717791 0.3742331 0.3067485
#> [2,] 0.1651786 0.1785714 0.3973214 0.2589286
#> [3,] 0.1908602 0.2419355 0.3252688 0.2419355
#> [4,] 0.1576763 0.2365145 0.3526971 0.2531120
t(post.true12) %*% post.true22/matrix(colSums(post.true12), nrow = 2^K, ncol = 2^K, byrow = F)
#> [,1] [,2] [,3] [,4]
#> [1,] 0.2923977 0.4327485 0.2280702 0.04678363
#> [2,] 0.2522936 0.4449541 0.1972477 0.10550459
#> [3,] 0.2813370 0.4540390 0.2061281 0.05849582
#> [4,] 0.2817460 0.4365079 0.2023810 0.07936508Again, looks very close.
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)
#>
#> TRUE
#> 4000
prop.table(table(lca.true == lca.hat))
#>
#> TRUE
#> 1
xtabs( ~ lca.true + lca.hat)
#> lca.hat
#> lca.true 1 2 3 4
#> 1 781 0 0 0
#> 2 0 1101 0 0
#> 3 0 0 1294 0
#> 4 0 0 0 824Pretttyy Preeeeeetttty Preeeeettttty good.