Fix a bug related to dummify

R/clm_MI_CR.R

@@ -17,8 +17,8 @@
 #'   \code{thin}.
 #'
 
-clm_MI_CR <- function(object, treatment, start=NULL, end=NULL, thin=NULL,
-                       exclude_chains=NULL, subset=FALSE, seed=NULL, mess=FALSE,...)
+clm_MI_CR <- function(object, treatment, start=NULL, end=NULL, thin=NULL,exclude_chains=NULL,
+                      subset=FALSE, ord_cov_dummy=TRUE, seed=NULL, mess=FALSE,...)
   {
   if (!inherits(object, "remiod")) stop("Use only with 'remiod' objects.")
 
@@ -77,7 +77,7 @@ clm_MI_CR <- function(object, treatment, start=NULL, end=NULL, thin=NULL,
     varname = vimp[j]
     coefs <- coef_list[[vimp[j]]]
     coefs <- merge(coefs, vdv, by="varname", all.x=TRUE)
-    #coefs$var[is.na(coefs$var)] = coefs$varname
+    coefs$var[is.na(coefs$var)] = coefs$varname
 
     misind_j = subset(misind, mvar==vimp[j])
 
@@ -94,27 +94,32 @@ clm_MI_CR <- function(object, treatment, start=NULL, end=NULL, thin=NULL,
       cvar = unique(coefs$var)[which(unique(coefs$var) %in% misind_h$mvar)]
 
       if (length(cvar)>0) {
-
-        dsmat0 = lapply(cvar, function(vk){
-          ## in CR, treatment effect need to removed for subjects in treated arm
-          contr_vk = attr(refs[[vk]],"contr_matrix")
-          contr_vk = data.frame(lev=rownames(contr_vk), contr_vk)
-          mcc = MCMC[, misind_h$mc_imp_col_nam[misind_h$mvar==vk],drop=F]
-          colnames(mcc) = "lev"
-          left_join(mcc,contr_vk, by="lev")[,-1]
-        })
-        dsmat0 = as.matrix(do.call(cbind.data.frame, dsmat0))
-        if (any(is.na(dsmat0))) dsmat0[is.na(dsmat0)] <- 0
-
-        dsmat1 = lapply(cvar, function(vk){
-          contr_vk = attr(refs[[vk]],"contr_matrix")
-          contr_vk = data.frame(lev=rownames(contr_vk), contr_vk)
-          mcc = MCMC_CR[, misind_h$mc_imp_col_nam[misind_h$mvar==vk],drop=F]
-          colnames(mcc) = "lev"
-          left_join(mcc,contr_vk, by="lev")[,-1]
-        })
-        dsmat1 = as.matrix(do.call(cbind.data.frame, dsmat1))
-        if (any(is.na(dsmat1))) dsmat1[is.na(dsmat1)] <- 0
+        if (ord_cov_dummy){
+          dsmat0 = lapply(cvar, function(vk){
+            ## in CR, treatment effect need to removed for subjects in treated arm
+            contr_vk = attr(refs[[vk]],"contr_matrix")
+            contr_vk = data.frame(lev=rownames(contr_vk), contr_vk)
+            mcc = MCMC[, misind_h$mc_imp_col_nam[misind_h$mvar==vk],drop=F]
+            colnames(mcc) = "lev"
+            left_join(mcc,contr_vk, by="lev")[,-1]
+          })
+          dsmat0 = as.matrix(do.call(cbind.data.frame, dsmat0))
+          if (any(is.na(dsmat0))) dsmat0[is.na(dsmat0)] <- 0
+
+          dsmat1 = lapply(cvar, function(vk){
+            contr_vk = attr(refs[[vk]],"contr_matrix")
+            contr_vk = data.frame(lev=rownames(contr_vk), contr_vk)
+            mcc = MCMC_CR[, misind_h$mc_imp_col_nam[misind_h$mvar==vk],drop=F]
+            colnames(mcc) = "lev"
+            left_join(mcc,contr_vk, by="lev")[,-1]
+          })
+          dsmat1 = as.matrix(do.call(cbind.data.frame, dsmat1))
+          if (any(is.na(dsmat1))) dsmat1[is.na(dsmat1)] <- 0
+        } else {
+          misind_ho = misind_h[match(cvar, misind_h$mvar), , drop = FALSE] # ordering monitored imp cols
+          dsmat0 = MCMC[, misind_ho$mc_imp_col_nam,drop=F]
+          dsmat1 = MCMC_CR[, misind_ho$mc_imp_col_nam,drop=F]
+        }
 
         coef_t = coefs$coef[coefs$var %in% cvar]
         eta = eta - apply(MCMC[,coef_t,drop=FALSE] * dsmat0, 1, sum) +

R/clm_MI_J2R.R

@@ -1,6 +1,7 @@
 #' Apply Jump-to-Reference(J2R) Method to Update JAGS MCMC outputs under MAR for Cumulative Logistic Model
 #'
 #' Internal function to obtain Jump-to-Reference(J2R) MCMC from an MAR object.
+#' @inheritParams commParams
 #' @param object an object of class remiod
 #' @param treatment the variable name of treatment. Reference level of treatment should be coded as 0.
 #' @param start first iteration to be used.
@@ -17,8 +18,8 @@
 #'   \code{thin}.
 #'
 
-clm_MI_J2R <- function(object, treatment, start=NULL, end=NULL, thin=NULL,
-                      exclude_chains=NULL,subset=FALSE, seed=NULL, mess=FALSE,...)
+clm_MI_J2R <- function(object, treatment, start=NULL, end=NULL, thin=NULL, exclude_chains=NULL,
+                       subset=FALSE, ord_cov_dummy=TRUE, seed=NULL, mess=FALSE,...)
   {
   if (!inherits(object, "remiod")) stop("Use only with 'remiod' objects.")
 
@@ -85,7 +86,7 @@ clm_MI_J2R <- function(object, treatment, start=NULL, end=NULL, thin=NULL,
     varname = vimp[j]
     coefs <- coef_list[[vimp[j]]]
     coefs <- merge(coefs, vdv, by="varname", all.x=TRUE)
-    #coefs$var[is.na(coefs$var)] = coefs$varname
+    coefs$var[is.na(coefs$var)] = coefs$varname
 
     misind_j = subset(misind, mvar==vimp[j])
 
@@ -120,27 +121,32 @@ clm_MI_J2R <- function(object, treatment, start=NULL, end=NULL, thin=NULL,
       cvar = unique(coefs$var)[which(unique(coefs$var) %in% misind_h$mvar)]
 
       if (length(cvar)>0){
-        dsmat0 = lapply(cvar, function(vk){
-          ## in J2R, treatment effect need to removed for subjects in treated arm
-          contr_vk = attr(refs[[vk]],"contr_matrix")
-          contr_vk = data.frame(lev=rownames(contr_vk), contr_vk)
-          mcc = MCMC[, misind_h$mc_imp_col_nam[misind_h$mvar==vk],drop=F]
-          colnames(mcc) = "lev"
-          left_join(mcc,contr_vk, by="lev")[,-1]
-        })
-        dsmat0 = as.matrix(do.call(cbind.data.frame, dsmat0))
-        if (any(is.na(dsmat0))) dsmat0[is.na(dsmat0)] <- 0
-
-        dsmat1 = lapply(cvar, function(vk){
-          contr_vk = attr(refs[[vk]],"contr_matrix")
-          contr_vk = data.frame(lev=rownames(contr_vk), contr_vk)
-          mcc = MCMC_J2R[, misind_h$mc_imp_col_nam[misind_h$mvar==vk],drop=F]
-          colnames(mcc) = "lev"
-          left_join(mcc,contr_vk, by="lev")[,-1]
-        })
-        dsmat1 = as.matrix(do.call(cbind.data.frame, dsmat1))
-        if (any(is.na(dsmat1))) dsmat1[is.na(dsmat1)] <- 0
-
+        if (ord_cov_dummy){
+          dsmat0 = lapply(cvar, function(vk){
+            ## in J2R, treatment effect need to removed for subjects in treated arm
+            contr_vk = attr(refs[[vk]],"contr_matrix")
+            contr_vk = data.frame(lev=rownames(contr_vk), contr_vk)
+            mcc = MCMC[, misind_h$mc_imp_col_nam[misind_h$mvar==vk],drop=F]
+            colnames(mcc) = "lev"
+            left_join(mcc,contr_vk, by="lev")[,-1]
+          })
+          dsmat0 = as.matrix(do.call(cbind.data.frame, dsmat0))
+          if (any(is.na(dsmat0))) dsmat0[is.na(dsmat0)] <- 0
+
+          dsmat1 = lapply(cvar, function(vk){
+            contr_vk = attr(refs[[vk]],"contr_matrix")
+            contr_vk = data.frame(lev=rownames(contr_vk), contr_vk)
+            mcc = MCMC_J2R[, misind_h$mc_imp_col_nam[misind_h$mvar==vk],drop=F]
+            colnames(mcc) = "lev"
+            left_join(mcc,contr_vk, by="lev")[,-1]
+          })
+          dsmat1 = as.matrix(do.call(cbind.data.frame, dsmat1))
+          if (any(is.na(dsmat1))) dsmat1[is.na(dsmat1)] <- 0
+        } else {
+          misind_ho = misind_h[match(cvar, misind_h$mvar), , drop = FALSE] # ordering monitored imp cols
+          dsmat0 = MCMC[, misind_ho$mc_imp_col_nam,drop=F]
+          dsmat1 = MCMC_J2R[, misind_ho$mc_imp_col_nam,drop=F]
+        }
         coef_t = coefs$coef[coefs$var %in% cvar]
         eta = eta - apply(MCMC[,coef_t,drop=FALSE] * dsmat0, 1, sum) +
           apply(MCMC[,coef_t,drop=FALSE] * dsmat1, 1, sum)

R/clm_MI_delta.R

@@ -21,8 +21,8 @@
 #'   \code{thin}.
 #'
 
-clm_MI_delta <- function(object, treatment, delta=0, start=NULL, end=NULL, thin=NULL,
-                       exclude_chains=NULL, subset=FALSE, seed=NULL, mess=FALSE,...)
+clm_MI_delta <- function(object, treatment, delta=0, start=NULL, end=NULL, exclude_chains=NULL,
+                         thin=NULL, subset=FALSE, ord_cov_dummy=TRUE, seed=NULL, mess=FALSE,...)
   {
   if (!inherits(object, "remiod")) stop("Use only with 'remiod' objects.")
 
@@ -78,7 +78,7 @@ clm_MI_delta <- function(object, treatment, delta=0, start=NULL, end=NULL, thin=
     varname = vimp[j]
     coefs <- coef_list[[vimp[j]]]
     coefs <- merge(coefs, vdv, by="varname", all.x=TRUE)
-    #coefs$var[is.na(coefs$var)] = coefs$varname
+    coefs$var[is.na(coefs$var)] = coefs$varname
 
     misind_j = subset(misind, mvar==vimp[j])
 
@@ -95,29 +95,35 @@ clm_MI_delta <- function(object, treatment, delta=0, start=NULL, end=NULL, thin=
       cvar = unique(coefs$var)[which(unique(coefs$var) %in% misind_h$mvar)]
 
       if (length(cvar)>0) {
-
-        dsmat0 = lapply(cvar, function(vk){
-          ## in J2R, treatment effect need to removed for subjects in treated arm
-          contr_vk = attr(refs[[vk]],"contr_matrix")
-          contr_vk = data.frame(lev=rownames(contr_vk), contr_vk)
-          mcc = MCMC[, misind_h$mc_imp_col_nam[misind_h$mvar==vk],drop=F]
-          colnames(mcc) = "lev"
-          left_join(mcc,contr_vk, by="lev")[,-1]
-        })
-        dsmat0 = as.matrix(do.call(cbind.data.frame, dsmat0))
-        if (any(is.na(dsmat0))) dsmat0[is.na(dsmat0)] <- 0
-
-        dsmat1 = lapply(cvar, function(vk){
-          contr_vk = attr(refs[[vk]],"contr_matrix")
-          contr_vk = data.frame(lev=rownames(contr_vk), contr_vk)
-          mcc = MCMC_delta[, misind_h$mc_imp_col_nam[misind_h$mvar==vk],drop=F]
-          colnames(mcc) = "lev"
-          left_join(mcc,contr_vk, by="lev")[,-1]
-        })
-        dsmat1 = as.matrix(do.call(cbind.data.frame, dsmat1))
-        if (any(is.na(dsmat1))) dsmat1[is.na(dsmat1)] <- 0
+        if (ord_cov_dummy){
+          dsmat0 = lapply(cvar, function(vk){
+            ## in J2R, treatment effect need to removed for subjects in treated arm
+            contr_vk = attr(refs[[vk]],"contr_matrix")
+            contr_vk = data.frame(lev=rownames(contr_vk), contr_vk)
+            mcc = MCMC[, misind_h$mc_imp_col_nam[misind_h$mvar==vk],drop=F]
+            colnames(mcc) = "lev"
+            left_join(mcc,contr_vk, by="lev")[,-1]
+          })
+          dsmat0 = as.matrix(do.call(cbind.data.frame, dsmat0))
+          if (any(is.na(dsmat0))) dsmat0[is.na(dsmat0)] <- 0
+
+          dsmat1 = lapply(cvar, function(vk){
+            contr_vk = attr(refs[[vk]],"contr_matrix")
+            contr_vk = data.frame(lev=rownames(contr_vk), contr_vk)
+            mcc = MCMC_delta[, misind_h$mc_imp_col_nam[misind_h$mvar==vk],drop=F]
+            colnames(mcc) = "lev"
+            left_join(mcc,contr_vk, by="lev")[,-1]
+          })
+          dsmat1 = as.matrix(do.call(cbind.data.frame, dsmat1))
+          if (any(is.na(dsmat1))) dsmat1[is.na(dsmat1)] <- 0
+        } else {
+          misind_ho = misind_h[match(cvar, misind_h$mvar), , drop = FALSE] # ordering monitored imp cols
+          dsmat0 = MCMC[, misind_ho$mc_imp_col_nam,drop=F]
+          dsmat1 = MCMC_delta[, misind_ho$mc_imp_col_nam,drop=F]
+        }
 
         coef_t = coefs$coef[coefs$var %in% cvar]
+
         eta = eta - apply(MCMC[,coef_t,drop=FALSE] * dsmat0, 1, sum) +
           apply(MCMC[,coef_t,drop=FALSE] * dsmat1, 1, sum)
       }

R/get_MI_RB.R

@@ -32,8 +32,8 @@
 #'
 #'
 get_MI_RB <- function(object, treatment, method=c("MAR","J2R","CR","delta"), delta=0,
-                      exclude_chains=NULL, start=NULL, end=NULL, seed=NULL,
-                      thin=NULL, subset=FALSE, include=TRUE, mess=TRUE, ...)
+                      exclude_chains=NULL, start=NULL, end=NULL, seed=NULL, thin=NULL,
+                      subset=FALSE, include=TRUE, ord_cov_dummy=TRUE, mess=TRUE, ...)
   {
   if(!missing(method) & length(method)>1) stop("Only one 'method' allowed.")
   method <- match.arg(method)
@@ -91,7 +91,7 @@ get_MI_RB <- function(object, treatment, method=c("MAR","J2R","CR","delta"), del
   MCMC = mcUpdateFun(object, treatment=treatment, delta=delta, seed = seed,
                      start = start, end = end, thin = thin,
                      subset = subset, exclude_chains = exclude_chains,
-                     mess = mess)
+                     ord_cov_dummy=ord_cov_dummy, mess = mess)
 
   # prepare a list of copies of the original data
   df_list <- list()

R/tang_MI_RB.R

@@ -16,7 +16,7 @@
 #' @param treatment treatment variable.
 
 tang_MI_RB = function(object, dtimp, treatment, method="MAR", delta=0, ord_cov_dummy=FALSE,
-                      exclude_chains=NULL, include=FALSE){
+                      exclude_chains=NULL, include=FALSE, thin=1){
   Mlist = get_Mlist(object)
   infolist = object$info_list
   coefs = object$coef_list
@@ -76,13 +76,29 @@ tang_MI_RB = function(object, dtimp, treatment, method="MAR", delta=0, ord_cov_d
   }
   else {
     bmcmc = object$MCMC
+    thin.old = attr(bmcmc[[1]], "mcpar")[3]
+    if (thin > thin.old) thin.new = thin %/% thin.old
+    else thin.new = thin.old
+
     chains = seq_along(bmcmc)
     if (!is.null(exclude_chains)) {
       chains <- chains[-exclude_chains]
       bmcmc = bmcmc[chains]
+      dtimc = dtimp[chains]
+    }
+
+    if (thin > 1) {
+      bmcmt = window(bmcmc, thin=thin)
+
+      dtimp =  lapply(dtimc, function(dx) {
+        dk = seq(1, max(dx$Imputation_), by= thin.new)
+        dxs = subset(dx, Imputation_ %in% dk)
+        dxs$Imputation_ = 1 + (dxs$Imputation_ %/% thin.new)
+        dxs
+      })
     }
 
-    bmcmc = lapply(bmcmc, function(x) as.matrix(x))
+    bmcmc = lapply(bmcmt, function(x) as.matrix(x))
     dtims = list()
 
     for (k in 1:length(bmcmc)){

R/tang_seq_helper.R

@@ -308,6 +308,7 @@ beta_mat_to_list = function(beta, coefs, infolist){
     names(cexpv) = v
     clist = c(clist, cexpv)
   }
+
   return(list(b=blist, cexps=clist))
 }
 

R/tang_seq_imp.R

@@ -150,7 +150,12 @@ tang_seq_imp = function(object, beta.init=NULL, ord_cov_dummy, seed = 1234, rinv
     if (t>burnin){
       tt = t - burnin
 
-      if (tt %% thin == 0){
+      if (thin==1){
+        dtimp = cbind(Imputation_ = tt, datat)
+        if (tt==1) dtimps = dtimp
+        else dtimps = rbind(dtimps, dtimp)
+      }
+      else if (tt %% thin == 0){
         ts = ts + 1
         dtimp = cbind(Imputation_ = ts, datat)
         if (ts==1) dtimps = dtimp

R/updateMI.R

@@ -47,24 +47,26 @@ updateMI <- function(object, method=c("MAR","J2R","CR","delta"), delta=0,
     t0 <- Sys.time()
     if (algorithm == "tang_seq"){
       dimp = tang_MI_RB(object=mcsamp, dtimp=dtimp, treatment=trtvar, method=method,
-                        ord_cov_dummy=ord_cov_dummy,
-                        exclude_chains=exclude_chains, include=include)
+                        ord_cov_dummy=ord_cov_dummy, exclude_chains=exclude_chains,
+                        include=include, thin=thin)
     } else{
       dimp = get_MI_RB(object=mcsamp, delta=delta, treatment=trtvar, method=method,
                        exclude_chains=exclude_chains, thin=thin, include=include,
-                       start=start, end=end, seed = seed)
+                       ord_cov_dummy=ord_cov_dummy,start=start, end=end, seed = seed)
     }
     t1 <- Sys.time()
   } else {
-    future_info <- get_future_info()
-
-    run_delta <- ifelse(future_info$parallel, delta_parallel, delta_seq)
+    if (length(delta)==1) run_delta <- delta_seq
+    else {
+      future_info <- get_future_info()
+      run_delta <- ifelse(future_info$parallel, delta_parallel, delta_seq)
+    }
 
     t0 <- Sys.time()
     dimp <- run_delta(object=mcsamp, dtimp=dtimp, algorithm=algorithm, delta=delta,
                       treatment=trtvar, method=method, n_workers = future_info$workers,
                       exclude_chains=exclude_chains, thin=thin, include=include,
-                      start=start, end=end, seed = seed)
+                      ord_cov_dummy=ord_cov_dummy,start=start, end=end, seed = seed)
     t1 <- Sys.time()
   }