inst/example/compare_flat.R

In zhenkewu/slamR: Structured Latent Attribute Models in R

# compare DMR with flattened analysis
# Mar 04, 2020

\dontrun{
  library(slamR)
  library(gplots) # heatmap.2
  rm(list=ls())

  production_dir <- "/Users/zhenkewu/Dropbox/OptumInsight\ Data\ and\ Restricted\
Latent\ Class\ Model/DMR_R_code/"

  N <- 15000 # sample size.
  err_prob <- 0.2 # noise level.

  # shrinkage algorithm parameters:
  thres_c <- 0.01  # used in the E step (modified EM for log-type penalized likelihood).
  thres   <- 0.5/N # for thresholding at the end of the modified EM algoritm or
  # general fractional power (equivalent formulation of pen-likelihood)
  # variational (E step for attributes use Dirichlet variational family) EM.

  C1 <- 3 # level 1 # of latent attribute patterns.


  # Simulation settings:
  A_set1   <- rbind(c(0,0,1),c(1,1,0),c(1,1,1))
  Q1_small <- rbind(diag(1,3),c(1,1,0),c(0,0,1))

  t(get_ideal_resp(Q1_small,A_set1))

  Q1 <- do.call("rbind", rep(list(Q1_small), 20))

  J1 <- nrow(Q1) # level 1 dimension, J1 (here = 200)
  K1 <- ncol(Q1)

  # level 2 structural matrix:
  Q2 <- vector("list",3)
  Q2[[1]] <- do.call("rbind", rep(list(rbind(diag(1,2),c(1,0),c(0,1),c(1,1))), 200))
  Q2[[2]] <- do.call("rbind", rep(list(rbind(diag(1,2),c(1,0),c(1,1),c(1,1))), 200))
  Q2[[3]] <- do.call("rbind", rep(list(rbind(diag(1,2),c(0,1),c(1,1),c(1,1))), 200))

  J2 <- nrow(Q2[[1]])
  K2 <- ncol(Q2[[1]]) # here we use identical K2 in the simulation, though need not be.


  # specify the taxonomy via D_mat
  D_mat <- matrix(0,J1,J2)
  J_ratio <- J2/J1
  for (j1 in 1:J1){
    D_mat[j1,((j1-1)*J_ratio+1):(j1*J_ratio)] <- 1
  }

  heatmap.2(D_mat,dendrogram='none', Rowv=FALSE, Colv=FALSE,trace='none')

  # model parameters:
  p1 <- c(0.5,0.3,0.2) # in paper c(0.3,0.2,0.5)
  p2 <- vector("list",C1)

  # design 1:
  A_set2 <- vector("list",C1)
  A_set2[[1]] <- rbind(c(0,0),c(0,1),c(1,0),c(1,1))
  A_set2[[2]] <- rbind(c(0,0),c(0,1),c(1,1))
  A_set2[[3]] <- rbind(c(0,0),c(1,0),c(1,1))

  p2[[1]] <- c(0.3,0.2,0.3,0.2)
  p2[[2]] <- c(0.4,0.2,0.4)
  p2[[3]] <- c(0.3,0.3,0.4)

  c1 <- rep((1-err_prob),J1)
  g1 <- rep(err_prob,J1)

  # let the probabilities of different classes in the first
  # level to have different item parameters at level 2:
  c2 <- vector("list",C1)
  for (cc in 1:C1){
    c2[[cc]] <- rep(1-err_prob,J2)
  }
  g2 <- vector("list",C1)
  for (cc in 1:C1){
    g2[[cc]] <- rep(err_prob,J2)
  }

  make_list(N,J1,J2,K1,K2)

  #
  # generate data
  #

  # level 1: coarser level
  set.seed(0513)
  res <- generate_X_fromA(N,A_set1,p1,Q1,c1,g1)
  X1 <- res$X
  Z1 <- res$Z
  X_true1 <- res$X_true
  rm("res")

  # level 2: finer level
  X2 <- matrix(0,nrow=N,ncol=J2)
  Z2 <- matrix(0,nrow=N,ncol=K2)

  for (cc in 1:C1){
    ind_cc <- which(bin2ind(Z1)==bin2ind(A_set1[cc,]))
    set.seed(0513)
    res <- generate_X_tax2(X1[ind_cc,,drop=FALSE],
                           A_set2[[cc]],
                           D_mat,
                           p2[[cc]], Q2[[cc]],c2[[cc]],g2[[cc]])

    X2[ind_cc,] <- res$X2
    Z2[ind_cc,] <- res$Z
  }

  rm("res")
  # model fitting:

  # stage 1 fitting:
  set.seed(0513)
  res <- get_initial_s(N,Q1,Z1)
  Q_ini1 <- res$Q_ini
  Z_ini1 <- res$Z_ini

  max_iter <- 50


  # This is only for tunning:
  #X=X1;Z_ini=Z_ini1;Q_ini=Q_ini1;max_iter=50
  time1 <- Sys.time()
  res <- adg_em(X1,Z_ini1,Q_ini1,max_iter,err_prob)
  Sys.time()-time1

  Q_est <- res$Q_arr[[length(res$Q_arr)]] # still may not identically recover Q?
  Z_est <- res$Z_est
  Z_candi <- res$Z_candi
  rm(res)

  check_complete(Q_est)
  # if not complete force to be complete.
  if (check_complete(Q_est)$is_complete==0){
    Q_est[1:K1,] <- diag(1,K1)
  }
  # ZW: does thia matter for estimating Q by forcing completeness?


  ## checking:
  cat("==look at estimated Q\n==")
  sum(sum(abs(get_ideal_resp(Q1,A_set1)-get_ideal_resp(Q_est,A_set1))>0))

  cat("==look at initial Q\n==")
  sum(sum(abs(get_ideal_resp(Q1,A_set1)-get_ideal_resp(Q_ini1,A_set1))>0))

  #estimated unique latent attribute profiles:
  unique(Z_est)
  table(bin2ind(Z_est)) # this is the candidates after screening; input for shrinkage.
  ## end of checking

  # shrinkage estimation at coarser level (level 1):
  c_ini1 <- c1
  g_ini1 <- g1

  lambda_vec <- seq(-0.2,-4.2,by=-0.4) # grid of lambda values in pen-likelihood formulation.

  res <- perform_shrink(X1, Q_est, Z_candi,
                        lambda_vec, c_ini1, g_ini1, 0)

  A_final1 <- res$A_final

  rm(res)

  res <- get_em_classify(X1,Q_est,A_final1,err_prob)
  Z_shrink1 <- res$Z_shrink


  pattern1 <- A_final1
  profile1 <- res$Z_shrink
  Q1_est   <- Q_est

  # end: 1st coarser fitting <---------------------------- end of first level fitting.

  #
  # begin 2nd finer resolution fitting:
  #
  table(bin2ind(Z1))        # true profiles.
  table(bin2ind(Z_shrink1)) # estimated. identical? this is excellent!

  C1_hat <- nrow(A_final1)
  pattern2_est <- vector("list",C1_hat)
  profile2_est <- vector("list",C1_hat)
  Q2_est       <- vector("list",C1_hat)

  Z_ini2 <- matrix(0,nrow=N,ncol=K2)
  Z_shrink2 <- matrix(0,nrow=N,ncol=K2) # currently we are assuming the same K2.


  must_maxiter <- 0
  set.seed(0513)

  for (cc in 1:C1_hat){
    ind_cc_est <- apply(Z_shrink1,1,function(v) all(v==A_final1[cc,]))
    X2_cc <- X2[ind_cc_est,,drop=FALSE]
    X1_cc <- X1[ind_cc_est,,drop=FALSE]

    #res <- get_initial_n(sum(ind_cc_est),Q2[[cc]],Z2[ind_cc_est,,drop=FALSE])
    # function not programmed.
    # caveat: the cc here may not match with the cc in the orginal simulation
    # A_final1 rows may be ordered differently than A_set1. Need to match!
    # in matlab - unique automatically order by row; A_set1 is ordered by row.



    #initialization matters: if using wrong Q2, Z2, might have problems:
    # the final Z_est might not be in Z_ini, for example.
    res <- get_initial_n(sum(ind_cc_est),Q2[[cc]],Z2[ind_cc_est,,drop=FALSE])
    Q_ini_t <- res$Q_ini
    Z_ini_t <- res$Z_ini

    max_iter <- 50

    res <- adg_em(X2_cc,Z_ini_t,Q_ini_t,
                  max_iter,err_prob,must_maxiter,D_mat,X1_cc)

    Z_est <- res$Z_est
    Z_candi <- res$Z_candi
    Q_arr   <- res$Q_arr

    Q_est <- Q_arr[[length(Q_arr)]]
    check_complete(Q_est)

    if (check_complete(Q_est)$is_complete==0){
      Q_est[1:K2,] <- diag(1,K2)
    }

    Q2_est[[cc]] <- Q_est


    ## some checks:
    table(bin2ind(Z_ini_t))

    unique(Z_est)
    table(bin2ind(Z_est))

    ind_cc <- apply(Z1,1,function(v) all(v==A_set1[cc,]))
    table(bin2ind(Z2[ind_cc,]))
    ## check end.

    # no shrinkage: after screening the latent attributes are estimated well.
    A_final <- Z_candi
    # but in level 1, we used PEM to choose A, plain EM for EBIC.

    table(bin2ind(Z2[ind_cc,]))
    table(bin2ind(Z_est))

    pattern2_est[[cc]] <- A_final
    profile2_est[[cc]] <- Z_est

    Z_ini2[ind_cc_est,] <- Z_ini_t
    Z_shrink2[ind_cc_est,] <- Z_est
  }


  #
  # NEXT: look at estimation results.
  #


  sum(abs(profile1-Z1),1)

  image(f(profile1-Z1))

  par(mfrow=c(1,2))
  image(f(Z_ini2-Z2))
  image(f(Z_shrink2-Z2))

  # combine
  par(mfrow=c(1,2))
  image(f(cbind(Z_ini1,Z_ini2)-cbind(Z1,Z2)))
  image(f(cbind(profile1,Z_shrink2)-cbind(Z1,Z2)))

  #combine results from doubly-multireolution clustering
  sum(abs(cbind(profile1,Z_shrink2)-cbind(Z1,Z2)))
  Z_multi_res <- cbind(profile1,Z_shrink2)

  for (cc in 1:C1){
    ind_cc <- apply(Z1,1,function(v) all(v==A_set1[cc,]))

    print(table(bin2ind(Z2[ind_cc,])))
    print(table(bin2ind(profile2_est[[cc]])))
  }


  ind_cc_arr <- vector("list",C1)
  png(file.path(production_dir,"data_level2.png"))
  par(mfrow=c(1,3))
  for (cc in 1:C1){
    ind_cc_arr[[cc]] <- apply(Z1,1,function(v) all(v==A_set1[cc,]))
    image(f(X2[ind_cc_arr[[cc]],]))
  }
  dev.off()

  # check tree constraint:
  all(all(X2 <= X1 %*% D_mat))

  #look at clustering at the 1st level:

  png(file.path(production_dir,"clustering_level1.png"))
  par(mfrow=c(1,2))
  image(f(Z_ini1-Z1),main="pattern_ini-pattern_true")
  image(f(Z_shrink1-Z1),main="pattern_est-pattern_true")
  dev.off()

  # second level clustering:
  png(file.path(production_dir,"clustering_level2.png"))
  par(mfrow=c(C1,2))
  for (cc in 1:C1){
    image(f(Z_ini2[ind_cc_arr[[cc]],]-Z2[ind_cc_arr[[cc]],]),
          main=paste0(paste(A_set1[cc,],collapse = ""),"; FINER 2: pattern_ini-pattern_true"))
    image(f(Z_shrink2[ind_cc_arr[[cc]], ]-Z2[ind_cc_arr[[cc]], ]),
          main=paste0(paste(A_set1[cc,],collapse = ""),"; FINER 2: pattern_est-pattern_true"))
  }
  dev.off()

  #
  # for flattened pattersn:
  #
  Z_flat <- cbind(Z1,Z2)
  A_set_flat <- unique_sort_binmat(Z_flat)


  table(bin2ind(Z_flat))
  A_set_all <- unique_sort_binmat(Z_flat)

  #### fitting begins:
  ## fitting flattened model
  set.seed(0513)

  # get initial values:
  res <- get_initial_n(N,Q1,Z1)
  Q_ini1 <- res$Q_ini
  Z_ini1 <- res$Z_ini

  res <- get_initial_n(N,Q2[[1]],Z2)
  Q_ini2 <- res$Q_ini
  Z_ini2 <- res$Z_ini

  K_flat <- K1+K2

  Z_flat_ini <- cbind(Z_ini1,Z_ini2)

  # try different initialization for Q_flat:
  Q_flat_ini <- rbind(cbind(Q_ini1,matrix(runif(J1*K2)<0.5,nrow=J1,ncol=K2)),
                      cbind(matrix(runif(J1*K2)<0.5,nrow=J2,ncol=K1),Q_ini2))

  ## begin real fitting:
  set.seed(0513)
  # one-stage shrinkage of flattened data
  X <- cbind(X1,X2)

  mat_iter <- 50
  time1 <- Sys.time()
  res <- adg_em(X, Z_flat_ini, Q_flat_ini, max_iter, err_prob) # currently slow.
  Sys.time()-time1

  Z_est <- res$Z_est
  Z_candi <- res$Z_candi
  Q_arr   <- res$Q_arr

  Q_est <- Q_arr[[length(Q_arr)]]

  check_complete(Q_est)

  unique_sort_binmat(Z_est)
  table(bin2ind(Z_est))

  # check if screened patterns include the true flattened patterns.
  # check_pattern(A_set_flat,Z_candi)

  # look at clustering
  png(file.path(production_dir,"clustering_level2_flattened.png"))
  par(mfrow=c(1,3))
  image(f(Z_flat_ini),main="pattern_ini")
  image(f(Z_est),main="pattern_est")
  image(f(cbind(Z1,Z2)),main="pattern_true")
  dev.off()

  # difference
  png(file.path(production_dir,"difference_from_truth.png"))
  par(mfrow=c(1,3))
  image(f(Z_flat_ini-cbind(Z1,Z2)),main="initial-true")
  image(f(Z_est-cbind(Z1,Z2)),main="FLAT: est-true")
  image(f(Z_multi_res-cbind(Z1,Z2)),main="MULTI: est-true")
  dev.off()

  # # look at Q
  # png(file.path(production_dir,"Q_flat.png"))
  # par(mfrow=c(C1,3))
  # image(f(Q_flat_ini),main="Q flat ini")
  # image(f(Q_est),main="FLAT: est")
  # image(f(rbind(cbind(Q1,matrix(0,J1,K2)),
  #               cbind(matrix(0,J2,K1),Q2[[1]]))),main="MULTI: truth") # <-- change to same Q20
  # dev.off()

  save.image(file.path(production_dir,"example.RDATA"))

}