## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(CLCM)
library(nnet)

## ----Specify------------------------------------------------------------------
    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_1---------------------------------------------------------------
     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
      
      # 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
 

## ----Generate_2---------------------------------------------------------------
      # Check the LCA 
      prop.table(table(lca))
      # Check LCA across observed Groups:
      prop.table(xtabs( ~ lca + dat$Group), margin = 1)

## ----Generate_3---------------------------------------------------------------
      mod <- multinom(as.factor(lca) ~ dat$Group*dat$Time)
      coef(mod)
      t(Beta)

## ----Generate_4---------------------------------------------------------------
      # 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)


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

## ----Estimate-----------------------------------------------------------------
      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') )    

## ----Results_1----------------------------------------------------------------

      Beta # Generating parameter
      mod$lat.reg.param  # estimate
      

## ----Results_2----------------------------------------------------------------

      transition_matrix_clcm(mod = mod, stratification = FALSE)

      # 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)


## ----Results_3----------------------------------------------------------------

  # Estimated:
  transition_matrix_clcm(mod = mod, stratification = TRUE, covariate = 'Group')

  # 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)
  t(post.true12) %*% post.true22/matrix(colSums(post.true12), nrow = 2^K, ncol = 2^K, byrow = FALSE)
   

## ----Classification-----------------------------------------------------------

      lca.hat <- mod$dat$lca 
      lca.true <- mod$dat$true_lca      
      table(lca.true == lca.hat)
      prop.table(table(lca.true == lca.hat)) 
      xtabs( ~ lca.true + lca.hat)