last bugfix in dco

Author: snaketron

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: ClustIRR
 Type: Package
 Title: Clustering of immune receptor repertoires
-Version: 1.7.5
+Version: 1.7.6
 Authors@R: c(
     person("Simo", "Kitanovski", email = "simokitanovski@gmail.com", 
     role = c("aut", "cre"), comment=c(ORCID="0000-0003-2909-5376")),

R/dco.R

@@ -306,36 +306,21 @@ get_posterior_summaries <- function(cm,
 
     post_beta <- function(f, samples) {
 
-        if(length(samples) == 2) {
-            s <- data.frame(summary(f, par = "beta")$summary)
-            s <- s[, c("mean", "X50.", "X2.5.", "X97.5.", "n_eff", "Rhat")]
-            colnames(s) <- c("mean", "median", "L95", "H95", "n_eff", "Rhat")
-            
-            s$i <- 1:nrow(s)
-            m <- rownames(s)
-            m <- gsub(pattern = paste0("beta\\[|\\]"), replacement = '', x = m)
-            m <- do.call(rbind, strsplit(x = m, split = "\\,"))
-            s$s <- as.numeric(m[,1])
-            s$community <- as.numeric(m[,2])
-            s <- s[order(s$i, decreasing = FALSE),]
-        } 
-        else {
-            s <- data.frame(summary(f, par = "beta")$summary)
-            s <- s[, c("mean", "X50.", "X2.5.", "X97.5.", "n_eff", "Rhat")]
-            colnames(s) <- c("mean", "median", "L95", "H95", "n_eff", "Rhat")
-            
-            # maintain original index order
-            s$i <- 1:nrow(s)
-            m <- rownames(s)
-            m <- gsub(pattern = "beta\\[|\\]", replacement = '', x = m)
-            
-            m <- do.call(rbind, strsplit(x = m, split = "\\,"))
-            s$s <- as.numeric(m[,1])
-            s$community <- as.numeric(m[,2])
-            meta <- data.frame(s = 1:length(samples), sample = samples)
-            s <- merge(x = s, y = meta, by.x = "s", by.y = "s", all.x = TRUE)
-            s <- s[order(s$i, decreasing = FALSE),]
-        }
+        s <- data.frame(summary(f, par = "beta")$summary)
+        s <- s[, c("mean", "X50.", "X2.5.", "X97.5.", "n_eff", "Rhat")]
+        colnames(s) <- c("mean", "median", "L95", "H95", "n_eff", "Rhat")
+        
+        # maintain original index order
+        s$i <- 1:nrow(s)
+        m <- rownames(s)
+        m <- gsub(pattern = "beta\\[|\\]", replacement = '', x = m)
+        
+        m <- do.call(rbind, strsplit(x = m, split = "\\,"))
+        s$s <- as.numeric(m[,1])
+        s$community <- as.numeric(m[,2])
+        meta <- data.frame(s = 1:length(samples), sample = samples)
+        s <- merge(x = s, y = meta, by.x = "s", by.y = "s", all.x = TRUE)
+        s <- s[order(s$i, decreasing = FALSE),]
 
         return(s)
     }

inst/stan/dm.stan

@@ -16,7 +16,6 @@ data {
 
 
 transformed data {
-    int N_beta;
     int N_delta;
     int K_delta;
     K_delta = 0;
@@ -25,34 +24,22 @@ transformed data {
         K_delta = K;
         N_delta = N*(N-1)/2;
     }
-    if(N==2) {
-        N_beta = 1;
-    } else {
-        N_beta = N;
-    }
 }
 
 
 parameters {
     real <lower=0> kappa;
     vector [K] alpha;
-    vector [K] beta [N_beta];
+    vector [K] beta [N];
 }
 
 model {
     target += exponential_lpdf(2+kappa|0.01);
     target += normal_lpdf(alpha|0,3);
 
-    if(N==2) {
-        target += normal_lpdf(beta[1]|0, 1);
-        target += dm_lpmf(y[1]|kappa*softmax(alpha + beta[1]*1));
-        target += dm_lpmf(y[2]|kappa*softmax(alpha + beta[1]*-1));
-    } 
-    else {
-        for(i in 1:N) {
-            target += normal_lpdf(beta[i]|0, 1);
-            target += dm_lpmf(y[i]|kappa*softmax(alpha + beta[i]));
-        }
+    for(i in 1:N) {
+        target += normal_lpdf(beta[i]|0, 1);
+        target += dm_lpmf(y[i]|kappa*softmax(alpha + beta[i]));
     }
 }
 
@@ -64,32 +51,20 @@ generated quantities {
     vector [K_delta] epsilon [N_delta];
     int k;
 
-    if(N==2) {
-        p[1] = dirichlet_rng(kappa*softmax(alpha + beta[1]*1));
-        p[2] = dirichlet_rng(kappa*softmax(alpha + beta[1]*-1));
-        delta[1] = beta[1];
-        epsilon[1] = softmax(alpha+beta[1])-softmax(alpha-beta[1]);
-        log_lik[1] = dm_lpmf(y[1]|kappa*softmax(alpha + beta[1]));
-        log_lik[2] = dm_lpmf(y[2]|kappa*softmax(alpha - beta[1]));
-        y_hat[1] = multinomial_rng(p[1], sum(y[1]));
-        y_hat[2] = multinomial_rng(p[2], sum(y[2]));
-    } 
-    else {
-        for(i in 1:N) {
-            p[i] = dirichlet_rng(kappa*softmax(alpha + beta[i]));
-            y_hat[i] = multinomial_rng(p[i], sum(y[i]));
-            log_lik[i] = dm_lpmf(y[i]|kappa*softmax(alpha + beta[i]));
-        }
-        
-        k = 1;
-        for(i in 1:N) {
-            if(compute_delta==1) {
-                if(i != N) {
-                    for(j in (i+1):N) {
-                        delta[k] = beta[i]-beta[j];
-                        epsilon[k] = softmax(alpha+beta[i])-softmax(alpha+beta[j]);
-                        k = k + 1;
-                    }
+    for(i in 1:N) {
+        p[i] = dirichlet_rng(kappa*softmax(alpha + beta[i]));
+        y_hat[i] = multinomial_rng(p[i], sum(y[i]));
+        log_lik[i] = dm_lpmf(y[i]|kappa*softmax(alpha + beta[i]));
+    }
+    
+    k = 1;
+    for(i in 1:N) {
+        if(compute_delta==1) {
+            if(i != N) {
+                for(j in (i+1):N) {
+                    delta[k] = beta[i]-beta[j];
+                    epsilon[k] = softmax(alpha+beta[i])-softmax(alpha+beta[j]);
+                    k = k + 1;
                 }
             }
         }

man/get_beta_scatterplot.Rd

@@ -80,7 +80,9 @@ gcd <- detect_communities(graph = c$graph,
                           chains = c("CDR3a", "CDR3b"))
 
 # differential community occupancy analysis
-dco <- dco(community_occupancy_matrix = gcd$community_occupancy_matrix)
+dco <- dco(community_occupancy_matrix = gcd$community_occupancy_matrix,
+           mcmc_control = list(mcmc_chains = 2, 
+                               mcmc_cores = 2))
 
 # generate beta violin plots
 beta_scatterplot <- get_beta_scatterplot(beta = dco$posterior_summary$beta,