fix bugs for binary outcome

R/est_lucid.R

@@ -34,9 +34,9 @@
 #' @examples
 #' \dontrun{
 #' set.seed(10)
-#' fit1 <- est.lucid(G = sim1[, 1:10],Z = sim1[, 11:20],Y = as,matrix(sim1[, 21]), 
+#' fit1 <- est.lucid(G = sim1[, 1:10],Z = sim1[, 11:20],Y = as.matrix(sim1[, 21]), 
 #' K = 2, family = "binary")
-#' fit2 <- est.lucid(G = sim1[, 1:10],Z = sim1[, 11:20],Y = as,matrix(sim1[, 21]), 
+#' fit2 <- est.lucid(G = sim1[, 1:10],Z = sim1[, 11:20],Y = as.matrix(sim1[, 21]), 
 #' K = 2, family = "binary", 
 #' tune = def.tune(Select_Z = TRUE, Rho_Z_InvCov = 0.1, Rho_Z_CovMu = 90, 
 #' Select_G = TRUE, Rho_G = 0.02))

R/family.R

@@ -25,17 +25,17 @@ normal <- function(K, ...){
   }
   # update parameters for M step
   f.maxY <- function(Y, r, CoY, K, CoYnames){
-    if(!is.null(CoY)){
+    # if(!is.null(CoY)){
       Set0 <- as.data.frame(cbind(Y, r[, -1], CoY))
       colnames(Set0) <- c("Y", paste0("LC", 2:K), CoYnames)
       Yfit <- glm(as.formula(paste("Y ~", paste(colnames(Set0)[-1], collapse = " + "))), data = Set0, family = gaussian)
       beta <- summary(Yfit)$coefficients[, 1]
       beta[2:K] <- beta[1] + beta[2:K]
       sigma <- rep(sd(residuals(Yfit)), K)
-    } else{
-      beta <- sapply(1:K, function(x){sum(r[, x] * Y) / sum(r[, x]) })
-      sigma <- sqrt(colSums(r * apply(matrix(beta), 1, function(x){(x - Y)^2})) / colSums(r))
-    }
+    # } else{
+    #   beta <- sapply(1:K, function(x){sum(r[, x] * Y) / sum(r[, x]) })
+    #   sigma <- sqrt(colSums(r * apply(matrix(beta), 1, function(x){(x - Y)^2})) / colSums(r))
+    # }
     return(structure(list(beta = beta,
                           sigma = sigma)))
   }
@@ -94,7 +94,9 @@ binary <- function(K, ...){
     Set0 <- as.data.frame(cbind(Y, r[, -1], CoY))
     colnames(Set0) <- c("Y", paste0("LC", 2:K), CoYnames)
     Yfit <- glm(as.formula(paste("Y~", paste(colnames(Set0)[-1], collapse = "+"))), data = Set0, family ="binomial")
-    beta <- c(0, coef(Yfit)[-1])
+    beta <- coef(Yfit)
+    beta[2:K] <- beta[1] + beta[2:K] # log odds for each latent cluster
+    # beta <- c(0, coef(Yfit)[-1])
     return(structure(list(beta = beta,
                           sigma = NULL)))
   }
@@ -107,8 +109,10 @@ binary <- function(K, ...){
     for (i in 1:K) {
       var[[i]] <- sigma[, , i]
     }
-    gamma$beta[1:K] <- (gamma$beta[1:K] - gamma$beta[1:K][index == 1])[index]
-    names(gamma$beta)[1:K] <- paste0("LC", 1:K)
+    ref <- gamma$beta[1:K][index == 1]
+    gamma$beta[1:K] <- (gamma$beta[1:K] - ref)[index]
+    gamma$beta[1] <- ref
+    names(gamma$beta)[1:K] <- c("LC1(reference)", paste0("LC", 2:K))
     return(structure(list(beta = beta, 
                           mu = mu,
                           sigma = var,

R/summary_lucid.R

@@ -109,11 +109,11 @@ f.normal <- function(x, K, se){
 }
 
 f.binary <- function(x, K, se){
-  cat("(1) Y (binary outcome): odds ratio of Y for each latent cluster (covariate) \n")
+  cat("(1) Y (binary outcome): log odds (reference) and log odds ratio of Y for latent cluster (covariate) \n")
   gamma <- as.data.frame(x$beta)
   colnames(gamma) <- "Original"
   if(is.null(se)){
-    gamma$OR <- exp(gamma$Original)
+    gamma$Exp <- exp(gamma$Original)
   } else{
     gamma <- cbind(gamma, se[, 2:4], OR = exp(gamma[, 1]), OR.L = exp(se[, 3]), OR.U = exp(se[, 4]))
   }