version 0.0.3

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: BMRMM
 Title: An Implementation of the Bayesian Markov (Renewal) Mixed Models
-Version: 0.0.2
+Version: 0.0.3
 Authors@R: c(
     person(given = "Yutong",family = "Wu",role = c("aut", "cre"),email = "yutong.wu@utexas.edu"),
     person(given = "Abhra",family = "Sarkar",role = c("aut"),email = "abhra.sarkar@utexas.edu"))
@@ -15,6 +15,6 @@ Encoding: UTF-8
 LazyData: true
 RoxygenNote: 7.2.3
 NeedsCompilation: no
-Packaged: 2023-01-06 18:28:26 UTC; yutongwu
+Packaged: 2023-01-21 02:16:23 UTC; yutongwu
 Repository: CRAN
-Date/Publication: 2023-01-06 22:30:42 UTC
+Date/Publication: 2023-01-21 15:10:02 UTC

MD5

@@ -1,24 +1,24 @@
-cc45ca3e95cb08c365a1a105979f0953 *DESCRIPTION
+90accea17ce6cb44a8f1ba3be5040ffc *DESCRIPTION
 8bba2de927cf29bf78c1c70b2787d620 *LICENSE
 ad670136e2ee81545f9e37e70766a927 *NAMESPACE
-aa81672d05125ad19466d0641c7dd5f0 *R/BMRMM.R
+24cfbffd87934e6683245301d13527bb *R/BMRMM.R
 d4df9110d17419d9d90f99237fa27be2 *R/add_isi_as_state.R
 44ff7fe7989c8b106b55acb1d5833afe *R/foxp2.R
 a9eca572c90fb926e8293ee0b7902ffc *R/foxp2_sm.R
-82706e702485cbabddd946f16ab684fe *R/model_cont_isi.R
+ba6c5fa6bb54fbb068ee253779e60278 *R/model_cont_isi.R
 86b0b5963ae673df0b53c249636b4cc4 *R/model_cont_isi_fn.R
-e4c746d0fd91a57707e68c44cf2a1509 *R/model_transition.R
+fc121a7dcc8f6b33e64defe9221b60a4 *R/model_transition.R
 85c58f28917633029b90e6649be93570 *R/model_transition_fn.R
 1ece9a13ceeaa4f1cf6e6916a9d004bc *R/print_and_plot_results.R
 c6ab5e36383884830d30ac97046393ba *R/print_and_plot_results_fn.R
 b215136422855e7e9427f18a2f7d5e4c *data/foxp2.RData
 10837a6621279c6f9fce23e4ab4e801d *data/foxp2_sm.RData
-dbdd5eb1dbf4dd176183e80f85f01d74 *man/BMRMM.Rd
+42b1d0ee2683432379dddced9d6db2b8 *man/BMRMM.Rd
 0a05dc4ec279b6fcb5f2d642f85c97c7 *man/foxp2.Rd
 b7fe6d25d1a3844411b6f405f48b99ad *man/foxp2_sm.Rd
 7e1aaf74710dd80c4c4b66942c7349e7 *man/get_duration_diagnostic_plots.Rd
 4b90a859643c4e774768ffba35a085ab *man/get_estimated_post_mean_and_sd.Rd
-e568a0ee4588270dc1d5956a797e921e *man/get_global_test_results.Rd
+5c2ca8708baf67b9a145ff4a9a8550c4 *man/get_global_test_results.Rd
 eb2e3d23b24147afb9ac81009115359a *man/get_histogram_by_component.Rd
 4bf0938296343d457e633beb7ef2cd82 *man/get_histogram_with_post_mean.Rd
 4c834c1db3c4ba9f6ecdbfd381a749f7 *man/get_lpml_waic_scores.Rd

R/BMRMM.R

@@ -80,6 +80,7 @@ NULL
 #' @param data a data frame containing -- individual ID, covariate values, previous state, current state, duration times (if applicable), in that order
 #' @param num_cov total number of covariates provided in 'data'
 #' @param switch_off_random TRUE if only population-level effects are considered, default is FALSE
+#' @param switch_off_fixed TRUE if only individual-level effects are considered, default is FALSE
 #' @param trans_cov_index indices of covariates that are used for transition probabilities, default is all of the covariates
 #' @param duration_type one of 'None', 'Discrete', 'Continuous', default is 'Continuous'
 #' @param duration_unit the discretization of duration times, only used when 'duration_type' is 'Discrete'
@@ -92,7 +93,8 @@ NULL
 #' @param burnin number of burn-ins for the MCMC iterations, default is simsize/2
 #' 
 #' @export
-BMRMM <- function(data,num_cov,switch_off_random=FALSE,trans_cov_index=1:num_cov,duration_type='Continuous',
+BMRMM <- function(data,num_cov,switch_off_random=FALSE,switch_off_fixed=FALSE,
+                  trans_cov_index=1:num_cov,duration_type='Continuous',
                   duration_cov_index=1:num_cov,duration_excl_prev_state=FALSE,duration_unit=NULL,
                   duration_num_comp=4,duration_init_shape=rep(1,duration_num_comp),
                   duration_init_rate=rep(1,duration_num_comp),simsize=10000,burnin=simsize/2) {
@@ -125,21 +127,21 @@ BMRMM <- function(data,num_cov,switch_off_random=FALSE,trans_cov_index=1:num_cov
   # get results
   res_duration <- NULL
   if(duration_type=='None') {
-    res_trans <- model_transition(data.frame(cbind(id,trans_covs,prev_state,cur_state)),switch_off_random,simsize,burnin)
+    res_trans <- model_transition(data.frame(cbind(id,trans_covs,prev_state,cur_state)),switch_off_random,switch_off_fixed,simsize,burnin)
   } else if(duration_type=='Discrete') {
     duration <- data[,num_col]
     if(is.null(duration_unit)) {
       stop("Must specify 'duration_unit' if treating duration time as a discrete variable.")
     } else {
       aug_data <- add_isi_as_state(data.frame(cbind(id,trans_covs,prev_state,cur_state,duration)),duration_unit)
-      res_trans <- model_transition(aug_data,switch_off_random,simsize,burnin)
+      res_trans <- model_transition(aug_data,switch_off_random,switch_off_fixed,simsize,burnin)
     }
   } else if(duration_type=='Continuous') {
     duration <- data[,num_col]
     data_for_trans <- data.frame(cbind(id,trans_covs,prev_state,cur_state))
     data_for_duration <- data.frame(cbind(id,duration_covs,duration))
-    res_trans <- model_transition(data_for_trans,switch_off_random,simsize,burnin)
-    res_duration <- model_cont_isi(data_for_duration,switch_off_random,simsize=simsize,burnin=burnin,K=duration_num_comp,duration_init_shape,duration_init_rate)
+    res_trans <- model_transition(data_for_trans,switch_off_random,switch_off_fixed,simsize,burnin)
+    res_duration <- model_cont_isi(data_for_duration,switch_off_random,switch_off_fixed,simsize=simsize,burnin=burnin,K=duration_num_comp,duration_init_shape,duration_init_rate)
   } else {
     stop("'duration_type' must be one of 'None', 'Discrete' and 'Continuous'.")
   }

R/model_cont_isi.R

@@ -2,7 +2,7 @@
 ### This function is used when ISI is modeled as mixture gammas ###
 ###################################################################
 
-model_cont_isi <- function(data,switch_off_random,simsize,burnin,K,init_alpha,init_beta) {
+model_cont_isi <- function(data,switch_off_random,switch_off_fixed,simsize,burnin,K,init_alpha,init_beta) {
 
   ################## process data ##################
   # construct isi data and covs for simulation
@@ -20,7 +20,12 @@ model_cont_isi <- function(data,switch_off_random,simsize,burnin,K,init_alpha,in
   num_combs <- nrow(all_combs)
   z.max <- K # number of mixture components
 
-  M00 <- covs.max
+  if (switch_off_fixed){
+    M00 <- ones(1,p)
+  } else {
+    M00 <- covs.max
+  }
+
   M <- matrix(0,nrow=simsize,ncol=p) # record # of clusters for each iteration
 
   clusts <- kmeans(isi,z.max)$cluster # get clusters via kmeans
@@ -31,11 +36,13 @@ model_cont_isi <- function(data,switch_off_random,simsize,burnin,K,init_alpha,in
   lambdaalpha <- 0.01 # alpha_{isi,0}
   lambda <- array(lambda00,dim=c(z.max,1)) # lambda_{g1,g2,g3}
   # individual probability and lambda
-  if (switch_off_random) {
-    piv <- matrix(1,nrow=max(mouseId),ncol=z.max) # probability for population-level effect
-  } else {
+  if (!switch_off_random && !switch_off_fixed) {
     piv <- matrix(0.8,nrow=max(mouseId),ncol=z.max) # probability for population-level effect
-  }
+  } else if (!switch_off_random) {
+    piv <- matrix(0,nrow=max(mouseId),ncol=z.max) # probability for population-level effect
+  } else { 
+    piv <- matrix(1,nrow=max(mouseId),ncol=z.max) # probability for population-level effect
+  } 
   v <- sample(0:1,num_row,prob=c(piv[1,1],1-piv[1,1]),replace=TRUE) # v=0: exgns, v=1: anmls
   lambda_mice <- matrix(rep(lambda00,each=max(mouseId)),nrow=max(mouseId))
   lambda_mice_alpha <- 0.01 # alpha_{isi}^{0}
@@ -48,11 +55,17 @@ model_cont_isi <- function(data,switch_off_random,simsize,burnin,K,init_alpha,in
   Beta <- matrix(1,nrow=simsize,ncol=z.max) # rate
 
   # cluster assignment for each covariate
-  z.cov <- as.matrix(covs) # initially, each covariate forms its own cluster
-  G <- matrix(0,nrow=p,ncol=max(covs.max))
-  for(pp in 1:p){
-    G[pp,1:covs.max[pp]] <- 1:covs.max[pp]
+  if (switch_off_fixed) {
+    z.cov <- ones(num_row,p) # initially, each covariate forms one cluster
+    G <- matrix(1,nrow=p,ncol=max(covs.max))
+  } else {
+    z.cov <- as.matrix(covs) # initially, each covariate forms its own cluster
+    G <- matrix(0,nrow=p,ncol=max(covs.max))
+    for(pp in 1:p){
+      G[pp,1:covs.max[pp]] <- 1:covs.max[pp]
+    }
   }
+
 
   # pM: prior probabilities of # of clusters for each covariate
   phi <- 0.25
@@ -81,7 +94,7 @@ model_cont_isi <- function(data,switch_off_random,simsize,burnin,K,init_alpha,in
     if(iii%%100==0) print(paste("Duration Times: Iteration ",iii))
 
     # update # cluster for each covariate and cluster mapping
-    if(iii > 1) {
+    if(iii > 1 && !switch_off_fixed) {
       # update cluster indicator z
       M00 <- M[iii,]
       for(pp in 1:p){
@@ -278,7 +291,7 @@ model_cont_isi <- function(data,switch_off_random,simsize,burnin,K,init_alpha,in
       }
     }
 
-    if (!switch_off_random) {
+    if (!switch_off_random || !switch_off_fixed) {
       # update v (v=0: exgns; v=1: anmls)
       mouse_k_pairs <- cbind(mouseId,z.vec)
       prob_exgns <- piv[mouse_k_pairs]*lambda[cbind(z.vec,z00)]
@@ -299,8 +312,10 @@ model_cont_isi <- function(data,switch_off_random,simsize,burnin,K,init_alpha,in
         vMouseMat <- vMat[,,id,]
         piv[id,] <- 1-rbeta(z.max,vMouseMat[1,]+1,vMouseMat[2,]+1)
       }
-
-      # update lambda_mice (lambda^(i))
+    }
+
+    # update lambda_mice (lambda^(i))
+    if (!switch_off_random) {
       for(id in 1:max(mouseId)){
         vMouseCt <- vMat[,,id,][2,]
         lambda_mice[id,] <- rdirichlet(1,lambda_mice_alpha*lambda0+vMouseCt)

R/model_transition.R

@@ -2,7 +2,7 @@
 ### This function is used when ISI is ignored or modeled as a discrete state ###
 ################################################################################
 
-model_transition <- function(data,switch_off_random,simsize,burnin) {
+model_transition <- function(data,switch_off_random,switch_off_fixed,simsize,burnin) {
 
   ################## process data ##################
   # data columns should be: id, selected covariates, prev state, current state
@@ -62,10 +62,12 @@ model_transition <- function(data,switch_off_random,simsize,burnin) {
   log0 <- zeros(N_MCMC,1)
   ind_all <- unique(sortrows(as.matrix(cbind(Xnew,Id))))
 
-  if (switch_off_random) {
-    piv <- ones(max(Id),d0) # probability for population-level effect
-  } else {
+  if (!switch_off_random && !switch_off_fixed) {
     piv <- 0.8*ones(max(Id),d0) # probability for population-level effect
+  } else if (!switch_off_random) {
+    piv <- zeros(max(Id),d0) # probability for population-level effect
+  } else {
+    piv <- ones(max(Id),d0) # probability for population-level effect
   }
   v <- sample(0:1,Ntot,TRUE,prob=c(1-piv[1,1],piv[1,1])) # v=0:anmls; v=1:exgns
 
@@ -123,6 +125,7 @@ model_transition <- function(data,switch_off_random,simsize,burnin) {
   Ts <- c(diff(m00starts),nrow(v0)-m00starts[length(m00starts)]+1)
   C[as.matrix(v00)] <- C[as.matrix(v00)]+Ts
   lambda0 <- C/repmat(colSums(C),d0,1)
+  lambda0[,is.nan(colSums(lambda0))] <- ones(size(lambda0)[1],1)/size(lambda0)[1]
   lambda0_mat_expanded_exgns <- repmat(lambda0,1,prod(dpreds))
   lambda0_mat_expanded_anmls <- repmat(lambda0,1,max(Id))
 
@@ -177,7 +180,7 @@ model_transition <- function(data,switch_off_random,simsize,burnin) {
 
     ### updating z (#clusters)
     M00 <- M[kkk,]
-    if(kkk>1) {
+    if(kkk>1 && !switch_off_fixed) {
       for(j in 1:p) {
         for(k in 1:dpreds[j]) {
           for(l in 1:dpreds[j]) {
@@ -207,7 +210,7 @@ model_transition <- function(data,switch_off_random,simsize,burnin) {
       pi[j,1:dpreds[j]] <- pi[j,1:dpreds[j]]/sum(pi[j,1:dpreds[j]])
     }
 
-    if (!switch_off_random) {
+    if (switch_off_random || switch_off_fixed) {
       ### updating v (v=0: anmls; v=1: exgns)
       prob <- zeros(Ntot,2)
       piv_mat <- matrix(piv,nrow=max(Id))