diff --git a/DESCRIPTION b/DESCRIPTION
index 579e88f0..b4de691a 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -4,7 +4,7 @@ Type: Package
 Title: Adjusting for Genetic Confounders in Transcriptome-Wide Association
   Studies Improves Discovery of Risk Genes of Complex Traits
 Date: 2024-06-14
-Version: 0.2.40
+Version: 0.2.41
 Authors@R: c(person("Siming","Zhao",role="aut",email="siming.zhao06@gmail.com"),
              person("Wesley","Crouse",role="aut"),
              person("Sheng","Qian",role="aut",email="shengqian@uchicago.edu"),
diff --git a/NAMESPACE b/NAMESPACE
index 3ba7ab61..ae6067ce 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -18,6 +18,7 @@ export(get_gene_annot_from_ens_db)
 export(get_problematic_genes)
 export(get_region_cor)
 export(load_weights)
+export(make_convergence_plots)
 export(make_locusplot)
 export(map_snp_info_regions)
 export(merge_finemap_regions)
@@ -77,7 +78,6 @@ importFrom(ggplot2,labs)
 importFrom(ggplot2,margin)
 importFrom(ggplot2,scale_alpha_manual)
 importFrom(ggplot2,scale_color_manual)
-importFrom(ggplot2,scale_colour_discrete)
 importFrom(ggplot2,scale_fill_manual)
 importFrom(ggplot2,scale_shape_manual)
 importFrom(ggplot2,scale_size_manual)
diff --git a/R/ctwas_plots.R b/R/ctwas_plots.R
index f9a314b6..7a0d7b65 100644
--- a/R/ctwas_plots.R
+++ b/R/ctwas_plots.R
@@ -238,7 +238,7 @@ make_locusplot <- function(finemap_res,
 
   # PIP panel
   p_pip <- ggplot(loc$data, aes(x=.data$pos/1e6, y=.data$susie_pip, shape=.data$type,
-                       size=.data$object_type, alpha=.data$object_type)) +
+                                size=.data$object_type, alpha=.data$object_type)) +
     geom_point(aes(color=.data$r2_levels, fill=.data$r2_levels)) +
     geom_text_repel(aes(label=.data$label), size=label.text.size, color="black") +
     scale_shape_manual(values = point.shapes) +
@@ -317,3 +317,139 @@ make_locusplot <- function(finemap_res,
   cowplot::plot_grid(p_pvalue, p_pip, p_qtl, p_genes, ncol = 1,
                      rel_heights = panel.heights, align = "v", axis="tblr")
 }
+
+#' @title Make convergence plots for the estimated parameters
+#'
+#' @param param a list of cTWAS parameter estimation result from \code{est_param}
+#'
+#' @param gwas_n the sample size of the GWAS summary statistics
+#'
+#' @param title.size font size of the plot title
+#'
+#' @param legend.size font size of the plot legend title
+#'
+#' @param colors colors for different groups, passed as the
+#'   \code{values} input to \code{\link[ggplot2]{scale_color_manual}}.
+#'   If fewer colors than "fits" are given, the colors are recycled.
+#'
+#' @param ncol number of columns in the output plot
+#'
+#' @importFrom ggplot2 ggplot
+#' @importFrom ggplot2 aes
+#' @importFrom ggplot2 geom_point
+#' @importFrom ggplot2 geom_line
+#' @importFrom ggplot2 xlab ylab ggtitle
+#' @importFrom ggplot2 expand_limits
+#' @importFrom ggplot2 theme
+#' @importFrom ggplot2 theme_bw
+#' @importFrom ggplot2 element_text
+#' @importFrom ggplot2 guides
+#' @importFrom ggplot2 guide_legend
+#' @importFrom ggplot2 scale_color_manual
+#' @importFrom cowplot plot_grid theme_cowplot
+#' @importFrom rlang .data
+#'
+#' @export
+make_convergence_plots <- function(param,
+                                   gwas_n,
+                                   title.size = 10,
+                                   legend.size = 8,
+                                   colors = c("#E69F00","#56B4E9","#009E73","#F0E442",
+                                              "#0072B2","#D55E00","#CC79A7", "#999999"),
+                                   ncol = 2){
+
+  # estimated group prior (all iterations)
+  group_prior_iters <- param$group_prior_iters
+  # estimated group prior variance (all iterations)
+  group_prior_var_iters <- param$group_prior_var_iters
+
+  # group size
+  group_size <- param$group_size[rownames(group_prior_iters)]
+
+  # estimated group PVE (all iterations)
+  group_pve_iters <- group_prior_var_iters*group_prior_iters*group_size/gwas_n
+  # estimated enrichment of genes (all iterations)
+  enrichment_iters <- t(sapply(rownames(group_prior_iters)[rownames(group_prior_iters)!="SNP"], function(x){
+    group_prior_iters[rownames(group_prior_iters)==x,]/group_prior_iters[rownames(group_prior_iters)=="SNP"]}))
+
+  # make convergence plots
+
+  # inclusion plot
+  df <- data.frame(niter = rep(1:ncol(group_prior_iters), nrow(group_prior_iters)),
+                   value = unlist(lapply(1:nrow(group_prior_iters), function(x){group_prior_iters[x,]})),
+                   group = rep(rownames(group_prior_iters), each=ncol(group_prior_iters)))
+  factor_levels <- c(setdiff(rownames(group_prior_iters), "SNP"), "SNP")
+  df$group <- factor(df$group, levels = factor_levels)
+
+  p_pi <- ggplot(df, aes(x=.data$niter, y=.data$value, group=.data$group, color=.data$group)) +
+    geom_line() +
+    geom_point() +
+    scale_color_manual(values = colors) +
+    xlab("Iteration") + ylab(bquote(pi)) +
+    ggtitle("Proportion Causal") +
+    theme_cowplot() +
+    theme(plot.title=element_text(size=title.size)) +
+    expand_limits(y=0) +
+    guides(color = guide_legend(title = "Group")) +
+    theme(legend.title = element_text(size=legend.size, face="bold"),
+          legend.text = element_text(size=legend.size))
+
+  # effect size plot
+  df <- data.frame(niter = rep(1:ncol(group_prior_var_iters), nrow(group_prior_var_iters)),
+                   value = unlist(lapply(1:nrow(group_prior_var_iters), function(x){group_prior_var_iters[x,]})),
+                   group = rep(rownames(group_prior_var_iters), each=ncol(group_prior_var_iters)))
+  df$group <- factor(df$group, levels = factor_levels)
+
+  p_sigma2 <- ggplot(df, aes(x=.data$niter, y=.data$value, group=.data$group, color=.data$group)) +
+    geom_line() +
+    geom_point() +
+    scale_color_manual(values = colors) +
+    xlab("Iteration") + ylab(bquote(sigma^2)) +
+    ggtitle("Effect Size") +
+    theme_cowplot() +
+    theme(plot.title=element_text(size=title.size)) +
+    expand_limits(y=0) +
+    guides(color = guide_legend(title = "Group")) +
+    theme(legend.title = element_text(size=legend.size, face="bold"),
+          legend.text = element_text(size=legend.size))
+
+  # PVE plot
+  df <- data.frame(niter = rep(1:ncol(group_pve_iters), nrow(group_pve_iters)),
+                   value = unlist(lapply(1:nrow(group_pve_iters), function(x){group_pve_iters[x,]})),
+                   group = rep(rownames(group_pve_iters), each=ncol(group_pve_iters)))
+  df$group <- factor(df$group, levels = factor_levels)
+
+  p_pve <- ggplot(df, aes(x=.data$niter, y=.data$value, group=.data$group, color=.data$group)) +
+    geom_line() +
+    geom_point() +
+    scale_color_manual(values = colors) +
+    xlab("Iteration") + ylab(bquote(h[G]^2)) +
+    ggtitle("PVE") +
+    theme_cowplot() +
+    theme(plot.title=element_text(size=title.size)) +
+    expand_limits(y=0) +
+    guides(color = guide_legend(title = "Group")) +
+    theme(legend.title = element_text(size=legend.size, face="bold"),
+          legend.text = element_text(size=legend.size))
+
+  # enrichment plot
+  df <- data.frame(niter = rep(1:ncol(enrichment_iters), nrow(enrichment_iters)),
+                   value = unlist(lapply(1:nrow(enrichment_iters), function(x){enrichment_iters[x,]})),
+                   group = rep(rownames(enrichment_iters), each=ncol(enrichment_iters)))
+  df$group <- factor(df$group, levels = factor_levels[factor_levels!="SNP"])
+
+  p_enrich <- ggplot(df, aes(x=.data$niter, y=.data$value, group=.data$group, color=.data$group)) +
+    geom_line() +
+    geom_point() +
+    scale_color_manual(values = colors) +
+    xlab("Iteration") + ylab(bquote(pi[G]/pi[S])) +
+    ggtitle("Enrichment") +
+    theme_cowplot() +
+    theme(plot.title=element_text(size=title.size)) +
+    expand_limits(y=0) +
+    guides(color = guide_legend(title = "Group")) +
+    theme(legend.title = element_text(size=legend.size, face="bold"),
+          legend.text = element_text(size=legend.size))
+
+  cowplot::plot_grid(p_pi, p_sigma2, p_enrich, p_pve, ncol = ncol)
+}
diff --git a/R/ctwas_summarize_parameters.R b/R/ctwas_summarize_parameters.R
index 03016237..e66b88ad 100644
--- a/R/ctwas_summarize_parameters.R
+++ b/R/ctwas_summarize_parameters.R
@@ -1,43 +1,14 @@
-#' @title Summarizes and plots cTWAS parameter estimates
+#' @title Summarizes estimated parameters
 #'
-#' @param param a list of cTWAS parameter estimation result from \code{est_param}
+#' @param param a list of parameter estimation result from \code{est_param}
 #'
 #' @param gwas_n the sample size of the GWAS summary statistics
 #'
-#' @param plot if TRUE, return a plot of the estimated parameters
-#'
-#' @param title.size font size of the plot title
-#'
-#' @param legend.size font size of the plot legend title
-#'
-#' @param title.font font of the plot title
-#'
-#' @importFrom logging loginfo
-#' @importFrom ggplot2 ggplot
-#' @importFrom ggplot2 aes
-#' @importFrom ggplot2 geom_point
-#' @importFrom ggplot2 geom_line
-#' @importFrom ggplot2 xlab ylab ggtitle
-#' @importFrom ggplot2 expand_limits
-#' @importFrom ggplot2 theme
-#' @importFrom ggplot2 theme_bw
-#' @importFrom ggplot2 element_text
-#' @importFrom ggplot2 guides
-#' @importFrom ggplot2 guide_legend
-#' @importFrom ggplot2 scale_colour_discrete
-#' @importFrom cowplot plot_grid theme_cowplot
-#' @importFrom rlang .data
+#' @return a list of summarized parameters
 #'
 #' @export
 #'
-summarize_param <- function(param,
-                            gwas_n,
-                            plot = TRUE,
-                            title.size = 10,
-                            legend.size = 8,
-                            title.font = c("bold", "plain", "italic", "bold.italic")){
-
-  title.font <- match.arg(title.font)
+summarize_param <- function(param, gwas_n){
 
   group_prior <- param$group_prior
   # estimated group prior (all iterations)
@@ -60,92 +31,13 @@ summarize_param <- function(param,
     group_prior_iters[rownames(group_prior_iters)==x,]/group_prior_iters[rownames(group_prior_iters)=="SNP"]}))
   enrichment <- enrichment_iters[,ncol(enrichment_iters)]
 
-  outlist <- list(group_size = group_size,
-                  group_prior = group_prior,
-                  group_prior_var = group_prior_var,
-                  enrichment = enrichment,
-                  group_pve = group_pve,
-                  total_pve = sum(group_pve),
-                  attributable_pve = group_pve/sum(group_pve))
-
-  if (plot){
-
-    # inclusion plot
-    df <- data.frame(niter = rep(1:ncol(group_prior_iters), nrow(group_prior_iters)),
-                     value = unlist(lapply(1:nrow(group_prior_iters), function(x){group_prior_iters[x,]})),
-                     group = rep(rownames(group_prior_iters), each=ncol(group_prior_iters)))
-    factor_levels <- c(setdiff(rownames(group_prior_iters), "SNP"), "SNP")
-    df$group <- factor(df$group, levels = factor_levels)
-
-    p_pi <- ggplot(df, aes(x=.data$niter, y=.data$value, group=.data$group, color=.data$group)) +
-      geom_line() +
-      geom_point() +
-      xlab("Iteration") + ylab(bquote(pi)) +
-      ggtitle("Proportion Causal") +
-      theme_cowplot() +
-      theme(plot.title=element_text(size=title.size)) +
-      expand_limits(y=0) +
-      guides(color = guide_legend(title = "Group")) +
-      theme(legend.title = element_text(size=legend.size, face=title.font),
-            legend.text = element_text(size=legend.size))
-
-    # effect size plot
-    df <- data.frame(niter = rep(1:ncol(group_prior_var_iters), nrow(group_prior_var_iters)),
-                     value = unlist(lapply(1:nrow(group_prior_var_iters), function(x){group_prior_var_iters[x,]})),
-                     group = rep(rownames(group_prior_var_iters), each=ncol(group_prior_var_iters)))
-    df$group <- factor(df$group, levels = factor_levels)
-
-    p_sigma2 <- ggplot(df, aes(x=.data$niter, y=.data$value, group=.data$group, color=.data$group)) +
-      geom_line() +
-      geom_point() +
-      xlab("Iteration") + ylab(bquote(sigma^2)) +
-      ggtitle("Effect Size") +
-      theme_cowplot() +
-      theme(plot.title=element_text(size=title.size)) +
-      expand_limits(y=0) +
-      guides(color = guide_legend(title = "Group")) +
-      theme(legend.title = element_text(size=legend.size, face=title.font),
-            legend.text = element_text(size=legend.size))
-
-    # PVE plot
-    df <- data.frame(niter = rep(1:ncol(group_pve_iters), nrow(group_pve_iters)),
-                     value = unlist(lapply(1:nrow(group_pve_iters), function(x){group_pve_iters[x,]})),
-                     group = rep(rownames(group_pve_iters), each=ncol(group_pve_iters)))
-    df$group <- factor(df$group, levels = factor_levels)
-
-    p_pve <- ggplot(df, aes(x=.data$niter, y=.data$value, group=.data$group, color=.data$group)) +
-      geom_line() +
-      geom_point() +
-      xlab("Iteration") + ylab(bquote(h[G]^2)) +
-      ggtitle("PVE") +
-      theme_cowplot() +
-      theme(plot.title=element_text(size=title.size)) +
-      expand_limits(y=0) +
-      guides(color = guide_legend(title = "Group")) +
-      theme(legend.title = element_text(size=legend.size, face=title.font),
-            legend.text = element_text(size=legend.size))
-
-    # enrichment plot
-    df <- data.frame(niter = rep(1:ncol(enrichment_iters), nrow(enrichment_iters)),
-                     value = unlist(lapply(1:nrow(enrichment_iters), function(x){enrichment_iters[x,]})),
-                     group = rep(rownames(enrichment_iters), each=ncol(enrichment_iters)))
-    df$group <- factor(df$group, levels = factor_levels[factor_levels!="SNP"])
-
-    p_enrich <- ggplot(df, aes(x=.data$niter, y=.data$value, group=.data$group, color=.data$group)) +
-      geom_line() +
-      geom_point() +
-      xlab("Iteration") + ylab(bquote(pi[G]/pi[S])) +
-      ggtitle("Enrichment") +
-      theme_cowplot() +
-      theme(plot.title=element_text(size=title.size)) +
-      expand_limits(y=0) +
-      guides(color = guide_legend(title = "Group")) +
-      theme(legend.title = element_text(size=legend.size, face=title.font),
-            legend.text = element_text(size=legend.size)) +
-      scale_colour_discrete(drop = FALSE)
-
-    outlist$convergence_plot <- plot_grid(p_pi, p_sigma2, p_enrich, p_pve)
-  }
+  res <- list(group_size = group_size,
+              group_prior = group_prior,
+              group_prior_var = group_prior_var,
+              enrichment = enrichment,
+              group_pve = group_pve,
+              total_pve = sum(group_pve),
+              attributable_pve = group_pve/sum(group_pve))
 
-  return(outlist)
+  return(res)
 }
diff --git a/man/summarize_param.Rd b/man/summarize_param.Rd
index a21a9f9b..d52c8e26 100644
--- a/man/summarize_param.Rd
+++ b/man/summarize_param.Rd
@@ -2,30 +2,18 @@
 % Please edit documentation in R/ctwas_summarize_parameters.R
 \name{summarize_param}
 \alias{summarize_param}
-\title{Summarizes and plots cTWAS parameter estimates}
+\title{Summarizes estimated parameters}
 \usage{
-summarize_param(
-  param,
-  gwas_n,
-  plot = TRUE,
-  title.size = 10,
-  legend.size = 8,
-  title.font = c("bold", "plain", "italic", "bold.italic")
-)
+summarize_param(param, gwas_n)
 }
 \arguments{
-\item{param}{a list of cTWAS parameter estimation result from \code{est_param}}
+\item{param}{a list of parameter estimation result from \code{est_param}}
 
 \item{gwas_n}{the sample size of the GWAS summary statistics}
-
-\item{plot}{if TRUE, return a plot of the estimated parameters}
-
-\item{title.size}{font size of the plot title}
-
-\item{legend.size}{font size of the plot legend title}
-
-\item{title.font}{font of the plot title}
+}
+\value{
+a list of summarized parameters
 }
 \description{
-Summarizes and plots cTWAS parameter estimates
+Summarizes estimated parameters
 }
diff --git a/tests/testthat/test-ctwas_summarize_parameters.R b/tests/testthat/test-ctwas_summarize_parameters.R
index dd1b4249..d07bc372 100644
--- a/tests/testthat/test-ctwas_summarize_parameters.R
+++ b/tests/testthat/test-ctwas_summarize_parameters.R
@@ -5,7 +5,7 @@ test_that("summarize_param works", {
   param <- ctwas_res$param
   precomputed_ctwas_parameters <- readRDS(system.file("extdata/sample_data", "LDL_example.ctwas_parameters_no_plots.RDS", package = "ctwas"))
 
-  ctwas_parameters <- summarize_param(param, gwas_n, plot = FALSE)
+  ctwas_parameters <- summarize_param(param, gwas_n)
   expect_equal(ctwas_parameters, precomputed_ctwas_parameters)
 
 })
diff --git a/vignettes/simple_ctwas_tutorial.Rmd b/vignettes/simple_ctwas_tutorial.Rmd
index fb3b25b7..eaea456f 100644
--- a/vignettes/simple_ctwas_tutorial.Rmd
+++ b/vignettes/simple_ctwas_tutorial.Rmd
@@ -245,10 +245,9 @@ ctwas_parameters$total_pve
 ctwas_parameters$attributable_pve
 ```
 
-
-Plot convergence
+Make convergence plots of estimated parameters
 ```{r convergence_plot, fig.width = 7, fig.height=5}
-ctwas_parameters$convergence_plot
+make_convergence_plots(param, gwas_n)
 ```
 
 
diff --git a/vignettes/summarizing_ctwas_results.Rmd b/vignettes/summarizing_ctwas_results.Rmd
index 54487369..852e8618 100644
--- a/vignettes/summarizing_ctwas_results.Rmd
+++ b/vignettes/summarizing_ctwas_results.Rmd
@@ -92,13 +92,12 @@ ctwas_parameters$attributable_pve
 ```
 
 
-Plot convergence
+Make convergence plots of estimated parameters
 ```{r convergence_plot, fig.width = 7, fig.height=5}
-ctwas_parameters$convergence_plot
+make_convergence_plots(param, gwas_n)
 ```
 
 
-
 ## Summarizing cTWAS results