diff --git a/multivariate-fine-mapping/mixture_prior.ipynb b/multivariate-fine-mapping/mixture_prior.ipynb index 75d3ac5..50b23d9 100644 --- a/multivariate-fine-mapping/mixture_prior.ipynb +++ b/multivariate-fine-mapping/mixture_prior.ipynb @@ -71,6 +71,12 @@ " --datadir $m --name `basename $m` &> bovy_$m.log\n", "```\n", "\n", + "### Plot results\n", + "\n", + "```\n", + "sos run mixture_prior.ipynb plot_U --model-data /path/to/mixture_model.rds --cwd output\n", + "```\n", + "\n", "Notice that for production use, each `sos run` command should be submitted to the cluster as a job." ] }, @@ -452,6 +458,8 @@ "# A list of models where we only update the scales and not the matrices\n", "# A typical choice is to estimate scales only for canonical components\n", "parameter: scale_only = []\n", + "# Tolerance for change in likelihood\n", + "parameter: ud_tol_lik = 1e-3\n", "input: [f\"{cwd}/{name}.rds\"] + [f\"{cwd}/{name}.{m}.rds\" for m in mixture_components]\n", "output: f'{cwd}/{name}.{ud_method}{\"_unconstrained\" if len(scale_only) == 0 else \"\"}{(\".\" + os.path.basename(resid_cor)[:-4]) if resid_cor.is_file() else \"\"}.rds'\n", "task: trunk_workers = 1, walltime = '36h', trunk_size = 1, mem = '10G', cores = 4, tags = f'{_output:bn}'\n", @@ -481,7 +489,7 @@ " library(udr)\n", " message(paste(\"Running ${ud_method.upper()} via udr package for\", length(U), \"mixture components\"))\n", " f0 = ud_init(X = as.matrix(dat$strong.z), V = V, U_scaled = U_scaled, U_unconstrained = U, n_rank1=0)\n", - " res = ud_fit(f0, X = na.omit(f0$X), control = list(unconstrained.update = \"${ud_method}\", resid.update = 'none', scaled.update = \"fa\", maxiter=5000, tol = 1e-06), verbose=TRUE)\n", + " res = ud_fit(f0, X = na.omit(f0$X), control = list(unconstrained.update = \"${ud_method}\", resid.update = 'none', scaled.update = \"fa\", maxiter=5000, tol.lik = ${ud_tol_lik}), verbose=TRUE)\n", " # add back col and row names to U\n", " # https://github.com/stephenslab/udr/issues/9\n", " for (i in 1:length(res$U)) {\n", @@ -499,6 +507,7 @@ "outputs": [], "source": [ "[ed]\n", + "parameter: ed_tol = 1e-6\n", "input: [f\"{cwd}/{name}.rds\"] + [f\"{cwd}/{name}.{m}.rds\" for m in mixture_components]\n", "output: f'{cwd}/{name}.ed_bovy{(\".\" + os.path.basename(resid_cor)[:-4]) if resid_cor.is_file() else \"\"}.rds'\n", "task: trunk_workers = 1, walltime = '36h', trunk_size = 1, mem = '10G', cores = 4, tags = f'{_output:bn}'\n", @@ -515,7 +524,7 @@ " # Fit mixture model using ED code by J. Bovy\n", " mash_data = mashr::mash_set_data(dat$strong.z, V=V)\n", " message(paste(\"Running ED via J. Bovy's code for\", length(U), \"mixture components\"))\n", - " res = mashr:::bovy_wrapper(mash_data, U, logfile=${_output:nr}, tol = 1e-06)\n", + " res = mashr:::bovy_wrapper(mash_data, U, logfile=${_output:nr}, tol = ${ed_tol})\n", " saveRDS(list(U=res$Ulist, w=res$pi, loglik=scan(\"${_output:n}_loglike.log\")), ${_output:r})" ] },