<!DOCTYPE html> <html xmlns="http://www.w3.org/1999/xhtml"> <head> <meta charset="utf-8" /> <meta http-equiv="Content-Type" content="text/html; charset=utf-8" /> <meta name="generator" content="pandoc" /> <meta name="author" content="Joyce Hsiao" /> <title>Predicting cell times for individuals</title> <script src="site_libs/jquery-1.11.3/jquery.min.js"></script> <meta name="viewport" content="width=device-width, initial-scale=1" /> <link href="site_libs/bootstrap-3.3.5/css/cosmo.min.css" rel="stylesheet" /> <script src="site_libs/bootstrap-3.3.5/js/bootstrap.min.js"></script> <script src="site_libs/bootstrap-3.3.5/shim/html5shiv.min.js"></script> <script src="site_libs/bootstrap-3.3.5/shim/respond.min.js"></script> <script src="site_libs/navigation-1.1/tabsets.js"></script> <link href="site_libs/highlightjs-9.12.0/textmate.css" rel="stylesheet" /> <script src="site_libs/highlightjs-9.12.0/highlight.js"></script> <link href="site_libs/font-awesome-5.0.13/css/fa-svg-with-js.css" rel="stylesheet" /> <script src="site_libs/font-awesome-5.0.13/js/fontawesome-all.min.js"></script> <script src="site_libs/font-awesome-5.0.13/js/fa-v4-shims.min.js"></script> <style type="text/css">code{white-space: pre;}</style> <style type="text/css"> pre:not([class]) { background-color: white; } </style> <script type="text/javascript"> if (window.hljs) { hljs.configure({languages: []}); hljs.initHighlightingOnLoad(); if (document.readyState && document.readyState === "complete") { window.setTimeout(function() { hljs.initHighlighting(); }, 0); } } </script> <style type="text/css"> h1 { font-size: 34px; } h1.title { font-size: 38px; } h2 { font-size: 30px; } h3 { font-size: 24px; } h4 { font-size: 18px; } h5 { font-size: 16px; } h6 { font-size: 12px; } .table th:not([align]) { text-align: left; } </style> </head> <body> <style type = "text/css"> .main-container { max-width: 940px; margin-left: auto; margin-right: auto; } code { color: inherit; background-color: rgba(0, 0, 0, 0.04); } img { max-width:100%; height: auto; } .tabbed-pane { padding-top: 12px; } .html-widget { margin-bottom: 20px; } button.code-folding-btn:focus { outline: none; } </style> <style type="text/css"> /* padding for bootstrap navbar */ body { padding-top: 51px; padding-bottom: 40px; } /* offset scroll position for anchor links (for fixed navbar) */ .section h1 { padding-top: 56px; margin-top: -56px; } .section h2 { padding-top: 56px; margin-top: -56px; } .section h3 { padding-top: 56px; margin-top: -56px; } .section h4 { padding-top: 56px; margin-top: -56px; } .section h5 { padding-top: 56px; margin-top: -56px; } .section h6 { padding-top: 56px; margin-top: -56px; } </style> <script> // manage active state of menu based on current page $(document).ready(function () { // active menu anchor href = window.location.pathname href = href.substr(href.lastIndexOf('/') + 1) if (href === "") href = "index.html"; var menuAnchor = $('a[href="' + href + '"]'); // mark it active menuAnchor.parent().addClass('active'); // if it's got a parent navbar menu mark it active as well menuAnchor.closest('li.dropdown').addClass('active'); }); </script> <div class="container-fluid main-container"> <!-- tabsets --> <script> $(document).ready(function () { window.buildTabsets("TOC"); }); </script> <!-- code folding --> <div class="navbar navbar-default navbar-fixed-top" role="navigation"> <div class="container"> <div class="navbar-header"> <button type="button" class="navbar-toggle collapsed" data-toggle="collapse" data-target="#navbar"> <span class="icon-bar"></span> <span class="icon-bar"></span> <span class="icon-bar"></span> </button> <a class="navbar-brand" href="index.html">fucci-seq</a> </div> <div id="navbar" class="navbar-collapse collapse"> <ul class="nav navbar-nav"> <li> <a href="index.html">Home</a> </li> <li> <a href="about.html">About</a> </li> <li> <a href="license.html">License</a> </li> </ul> <ul class="nav navbar-nav navbar-right"> <li> <a href="https://github.com/jdblischak/fucci-seq"> <span class="fa fa-github"></span> </a> </li> </ul> </div><!--/.nav-collapse --> </div><!--/.container --> </div><!--/.navbar --> <!-- Add a small amount of space between sections. --> <style type="text/css"> div.section { padding-top: 12px; } </style> <div class="fluid-row" id="header"> <h1 class="title toc-ignore">Predicting cell times for individuals</h1> <h4 class="author"><em>Joyce Hsiao</em></h4> </div> <div id="TOC"> <ul> <li><a href="#pipeline">Pipeline</a></li> <li><a href="#top-x-cyclical-genes">Top X cyclical genes</a></li> <li><a href="#prediction-error">Prediction error</a></li> <li><a href="#session-information">Session information</a></li> </ul> </div> <!-- The file analysis/chunks.R contains chunks that define default settings shared across the workflowr files. --> <!-- Update knitr chunk options --> <!-- Insert the date the file was last updated --> <p><strong>Last updated:</strong> 2018-07-18</p> <!-- Insert the code version (Git commit SHA1) if Git repository exists and R package git2r is installed --> <p><strong>Code version:</strong> 91067ec</p> <hr /> <div id="pipeline" class="section level2"> <h2>Pipeline</h2> <ol style="list-style-type: decimal"> <li><p>Split data to validation samples and training samples.</p> <ul> <li>The split samples are saved as <code>data/ind_XX_withheld.rds</code> and <code>data/ind_XX_training.rds</code>.<br /> </li> <li>The code used to generated teh split samples is saved in <code>code/working/makedata_supervised.R</code>.</li> </ul></li> <li><p>In the training sample, repeat the following 5 times. Each time, take samples from 4 individuals, fit trendfilter and identified the top X cyclical genes, and the compute the error in predicting cell time of the held-out individual, as a function of the top X cyclical genes.</p> <ul> <li><code>code/run_methods.R</code>: runs all methods, including supervised, unsupervised trendfilter, unsupervised bspline, unsupervised lowess, seurat) and compute prediction error margin.<br /> </li> <li><code>code/working/job_run_methods.train.ind.R</code> calls <code>code/run_methods.R</code></li> <li><code>code/working/job_run_methods.train.ind.sbatch</code> submits jobs one fold at a time.</li> </ul></li> <li><p>Based on the training sample results, select the top best X genes and fit to the validation sample. Compute and compare the methods’ performance on the validation sample.</p></li> </ol> <hr /> </div> <div id="top-x-cyclical-genes" class="section level2"> <h2>Top X cyclical genes</h2> <p>time based on fucci only</p> <pre class="r"><code># observed data library(gplots) inds <- c("NA19098", "NA18511", "NA18870", "NA19101", "NA18855", "NA19160") ngenes <- c(5, seq(10, 11040, 10)) for (j in 1:length(inds)) { ind <- inds[j] gene_names <- rownames(readRDS(paste0("../data/results/ind_",ind,"_data_training_cyclical_genes.fold.", 1,".rds"))) genes_list <- lapply(1:length(ngenes), function(i) { ngene <- ngenes[i] tmp <- do.call(cbind, lapply(1:5, function(fold) { fl_name <- paste0("../data/results/ind_",ind,"_data_training_cyclical_genes.fold.", fold,".rds") df <- readRDS(fl_name) top_list <- rownames(df)[order(df$pve,decreasing = T)[1:ngene]] rownames(df) %in% top_list }) ) rownames(tmp) <- gene_names return(rownames(tmp)[rowSums(tmp)>=4]) }) names(genes_list) <- ngenes saveRDS(genes_list, file = paste0("../data/results/ind_",ind,"_results_topgenes.rds")) } ind_NA19098_topgenes <- readRDS("../data/results/ind_NA19098_results_topgenes.rds") ind_NA18511_topgenes <- readRDS("../data/results/ind_NA18511_results_topgenes.rds") ind_NA18870_topgenes <- readRDS("../data/results/ind_NA18870_results_topgenes.rds") ind_NA19101_topgenes <- readRDS("../data/results/ind_NA19101_results_topgenes.rds") ind_NA18855_topgenes <- readRDS("../data/results/ind_NA18855_results_topgenes.rds") ind_NA19160_topgenes <- readRDS("../data/results/ind_NA19160_results_topgenes.rds") topgenes_list <- list(ind_NA19098_topgenes, ind_NA18511_topgenes, ind_NA18870_topgenes, ind_NA19101_topgenes, ind_NA18855_topgenes, ind_NA19160_topgenes) library(dplyr) topgenes_combined <- lapply(1:length(ind_NA19160_topgenes), function(i) { ll <- lapply(topgenes_list, "[[", i) ll <- unlist(ll) tab <- table(ll) names(tab)[which(tab>=4)] }) names(topgenes_combined) <- names(ind_NA19160_topgenes) saveRDS(topgenes_combined, file="../data/results/ind_results_topgenes.rds")</code></pre> <p>time based on fucci and dapi only</p> <pre class="r"><code># observed data library(gplots) inds <- c("NA19098", "NA18511", "NA18870", "NA19101", "NA18855", "NA19160") ngenes <- c(5, seq(10, 11040, 10)) for (j in 1:length(inds)) { ind <- inds[j] gene_names <- rownames(readRDS(paste0("../data/results/triple_ind_",ind, "_data_training_cyclical_genes.fold.", 1,".rds"))) genes_list <- lapply(1:length(ngenes), function(i) { ngene <- ngenes[i] tmp <- do.call(cbind, lapply(1:5, function(fold) { fl_name <- paste0("../data/results/triple_ind_",ind, "_data_training_cyclical_genes.fold.", fold,".rds") df <- readRDS(fl_name) top_list <- rownames(df)[order(df$pve,decreasing = T)[1:ngene]] rownames(df) %in% top_list }) ) rownames(tmp) <- gene_names return(rownames(tmp)[rowSums(tmp)>=4]) }) names(genes_list) <- ngenes saveRDS(genes_list, file = paste0("../data/results/ind_",ind,"_results_topgenes_triple.rds")) } ind_NA19098_topgenes <- readRDS("../data/results/ind_NA19098_results_topgenes_triple.rds") ind_NA18511_topgenes <- readRDS("../data/results/ind_NA18511_results_topgenes_triple.rds") ind_NA18870_topgenes <- readRDS("../data/results/ind_NA18870_results_topgenes_triple.rds") ind_NA19101_topgenes <- readRDS("../data/results/ind_NA19101_results_topgenes_triple.rds") ind_NA18855_topgenes <- readRDS("../data/results/ind_NA18855_results_topgenes_triple.rds") ind_NA19160_topgenes <- readRDS("../data/results/ind_NA19160_results_topgenes_triple.rds") topgenes_list <- list(ind_NA19098_topgenes, ind_NA18511_topgenes, ind_NA18870_topgenes, ind_NA19101_topgenes, ind_NA18855_topgenes, ind_NA19160_topgenes) library(dplyr) topgenes_combined <- lapply(1:length(ind_NA19160_topgenes), function(i) { ll <- lapply(topgenes_list, "[[", i) ll <- unlist(ll) tab <- table(ll) names(tab)[which(tab>=4)] }) names(topgenes_combined) <- names(ind_NA19160_topgenes) saveRDS(topgenes_combined, file="../data/results/ind_results_topgenes_triple.rds")</code></pre> <p>Compare number of genes selected at Top X cyclical genes.</p> <pre class="r"><code>topgenes_double <- readRDS(file="../data/results/ind_results_topgenes.rds") topgenes_triple <- readRDS(file="../data/results/ind_results_topgenes_triple.rds") topgenes_overlap <- sapply(1:length(topgenes_double), function(i) { intersect(topgenes_double[[i]], topgenes_triple[[i]]) }) names(topgenes_overlap) <- names(topgenes_double) plot(names(topgenes_double)[1:21], sapply(topgenes_double, length)[1:21], xlab="Top X cyclical genes", pch=16, col=1, cex=.7, ylim=c(0,60), ylab="Selected genes", main = "Number of selected genes") points(names(topgenes_double)[1:21], sapply(topgenes_triple, length)[1:21], col="blue", pch=16, cex=.7) points(names(topgenes_double)[1:21], sapply(topgenes_overlap, length)[1:21], col="1", pch=1, cex=.7) legend("topleft", legend=c("double", "triple", "both"), col=c("black", "blue", "black"), pch=c(16, 16, 1))</code></pre> <p><img src="figure/method-train-ind.Rmd/unnamed-chunk-3-1.png" width="672" style="display: block; margin: auto;" /></p> <hr /> </div> <div id="prediction-error" class="section level2"> <h2>Prediction error</h2> <p>summarize results</p> <pre class="r"><code>diff_time_wrapper <- function(results_list) { methods_list <- sapply(names(results_list), function(x) strsplit(x, split=".", fixed=TRUE)[[1]][2]) diff_time_list <- do.call(rbind, lapply(1:length(results_list), function(i) { diff_time <- results_list[[i]]$diff_time diff_mean <- mean(diff_time/2/pi) # diff_se <- sd(diff_time/2/pi)/sqrt(ncol(results_list[[1]]$Y)) return(data.frame(diff_mean=diff_mean, # diff_se=diff_se, methods=methods_list[i])) }) ) return(diff_time_list) } library(gplots) ngenes <- c(5, seq(10,200, by=10)) inds <- c("NA19098", "NA18511", "NA18870", "NA19101", "NA18855", "NA19160") train_top <- lapply(1:length(inds), function(j) { ind <- inds[j] out <- do.call(rbind, lapply(1:length(ngenes), function(i) { ngene <- ngenes[i] train_topX <- do.call(rbind, lapply(1:5, function(fold) { # print(ind) # print(ngene) # print(fold) fl_name <- paste0("../data/results/ind_",ind,"_results_train.fold.",fold, ".top",ngene,".rds") df <- readRDS(fl_name) out <- diff_time_wrapper(df$fit.test) out$fold <- fold return(out) }) ) train_topX$ngenes <- ngene #return(train_topX) agg_mn <- aggregate(diff_mean ~ methods, data=train_topX, FUN=mean) agg_sd <- aggregate(diff_mean ~ methods, data=train_topX, FUN=sd) obj <- data.frame(methods=agg_mn$methods, diff_mean=agg_mn$diff_mean, diff_se=agg_sd$diff_mean/sqrt(5)) obj$ngenes <- ngene return(obj) }) ) out$ind <- ind return(out) }) names(train_top) <- inds saveRDS(train_top, file = "../output/method-train-ind.Rmd/diff_time_all.rds") ngenes <- c(5, seq(10,200, by=10)) inds <- c("NA19098", "NA18511", "NA18870", "NA19101", "NA18855", "NA19160") train_triple <- lapply(1:length(inds), function(j) { ind <- inds[j] out <- do.call(rbind, lapply(1:length(ngenes), function(i) { ngene <- ngenes[i] train_topX <- do.call(rbind, lapply(1:5, function(fold) { # print(ind) # print(ngene) # print(fold) fl_name <- paste0("../data/results/triple_ind_",ind,"_results_train.fold.",fold, ".top",ngene,".rds") df <- readRDS(fl_name) out <- diff_time_wrapper(df$fit.test) out$fold <- fold return(out) }) ) train_topX$ngenes <- ngene #return(train_topX) agg_mn <- aggregate(diff_mean ~ methods, data=train_topX, FUN=mean) agg_sd <- aggregate(diff_mean ~ methods, data=train_topX, FUN=sd) obj <- data.frame(methods=agg_mn$methods, diff_mean=agg_mn$diff_mean, diff_se=agg_sd$diff_mean/sqrt(5)) obj$ngenes <- ngene return(obj) }) ) out$ind <- ind return(out) }) names(train_triple) <- inds saveRDS(train_triple, file = "../output/method-train-ind.Rmd/diff_time_all_triple.rds") # ngenes <- c(5, seq(10,200, by=10)) # train_top_permute <- do.call(rbind, lapply(1:length(ngenes), function(i) { # ngene <- ngenes[i] # train_topX <- do.call(rbind, lapply(1:5, function(fold) { # fl_name <- paste0("../data/results/results_train_permute_oncyclical.fold.", # fold,".top",ngene,".rds") # df <- readRDS(fl_name) # out <- diff_time_wrapper(df$fit.test) # out$fold <- fold # return(out) # }) ) # train_topX$ngenes <- ngene # #return(train_topX) # agg_mn <- aggregate(diff_mean ~ methods, # data=train_topX, FUN=mean) # agg_sd <- aggregate(diff_mean ~ methods, # data=train_topX, FUN=sd) # # obj <- data.frame(methods=agg_mn$methods, # diff_mean=agg_mn$diff_mean, # diff_se=agg_sd$diff_mean/sqrt(5)) # obj$ngenes <- ngene # return(obj) # }) ) # # saveRDS(train_top_permute, # file = "../output/method-train-classifiers-all.Rmd/diff_time_fucci_on_rando.rds")</code></pre> <p>Results</p> <pre class="r"><code>train_double <- readRDS("../output/method-train-ind.Rmd/diff_time_all.rds") train_double <- do.call(rbind, train_double) train_double <- subset(train_double, methods=="supervised") train_double$ind <- as.factor(train_double$ind) train_triple <- readRDS("../output/method-train-ind.Rmd/diff_time_all_triple.rds") train_triple <- do.call(rbind, train_triple) train_triple <- subset(train_triple, methods=="supervised") train_triple$ind <- as.factor(train_triple$ind)</code></pre> <p>time based on fucci</p> <pre class="r"><code>library(ggplot2) ggplot(train_double, aes(x=ngenes, y=diff_mean, col=ind)) + # geom_errorbar(train_top, # mapping=aes(ymin=diff_mean-diff_se, ymax=diff_mean+diff_se, # col=ind)) + geom_point() + geom_line(lty=3) + ylab("Error margin (% arc length)") + xlab("Top X cyclical genes") + ggtitle("Prediction error in the training set") + ylim(0,.25)</code></pre> <p><img src="figure/method-train-ind.Rmd/unnamed-chunk-5-1.png" width="672" style="display: block; margin: auto;" /></p> <pre class="r"><code>ggplot(train_double, aes(x=ngenes, y=diff_mean, col=ind)) + # geom_errorbar(train_top, # mapping=aes(ymin=diff_mean-diff_se, ymax=diff_mean+diff_se, # col=ind)) + geom_point() + geom_line(lty=3) + ylab("Error margin (% arc length)") + xlab("Top X cyclical genes") + ggtitle("Prediction error in the training set") </code></pre> <p><img src="figure/method-train-ind.Rmd/unnamed-chunk-5-2.png" width="672" style="display: block; margin: auto;" /></p> <pre class="r"><code>ggplot(train_double, aes(x=ngenes, y=diff_mean, col=ind)) + geom_errorbar(train_double, mapping=aes(ymin=diff_mean-diff_se, ymax=diff_mean+diff_se, col=ind)) + geom_point() + geom_line(lty=3) + ylab("Error margin (% arc length)") + xlab("Top X cyclical genes") + ggtitle("Prediction error in the training set") </code></pre> <p><img src="figure/method-train-ind.Rmd/unnamed-chunk-5-3.png" width="672" style="display: block; margin: auto;" /></p> <p>time based on fucci + dapi</p> <pre class="r"><code>library(ggplot2) ggplot(train_triple, aes(x=ngenes, y=diff_mean, col=ind)) + # geom_errorbar(train_top, # mapping=aes(ymin=diff_mean-diff_se, ymax=diff_mean+diff_se, # col=ind)) + geom_point() + geom_line(lty=3) + ylab("Error margin (% arc length)") + xlab("Top X cyclical genes") + ggtitle("Prediction error in the training set") + ylim(0,.25)</code></pre> <p><img src="figure/method-train-ind.Rmd/unnamed-chunk-6-1.png" width="672" style="display: block; margin: auto;" /></p> <pre class="r"><code>library(ggplot2) ggplot(train_triple, aes(x=ngenes, y=diff_mean, col=ind)) + # geom_errorbar(train_top, # mapping=aes(ymin=diff_mean-diff_se, ymax=diff_mean+diff_se, # col=ind)) + geom_point() + geom_line(lty=3) + ylab("Error margin (% arc length)") + xlab("Top X cyclical genes") + ggtitle("Prediction error in the training set") </code></pre> <p><img src="figure/method-train-ind.Rmd/unnamed-chunk-6-2.png" width="672" style="display: block; margin: auto;" /></p> <pre class="r"><code>ggplot(train_triple, aes(x=ngenes, y=diff_mean, col=ind)) + geom_errorbar(train_double, mapping=aes(ymin=diff_mean-diff_se, ymax=diff_mean+diff_se, col=ind)) + geom_point() + geom_line(lty=3) + ylab("Error margin (% arc length)") + xlab("Top X cyclical genes") + ggtitle("Prediction error in the training set") </code></pre> <p><img src="figure/method-train-ind.Rmd/unnamed-chunk-6-3.png" width="672" style="display: block; margin: auto;" /></p> <hr /> </div> <div id="session-information" class="section level2"> <h2>Session information</h2> <pre class="r"><code>sessionInfo()</code></pre> <pre><code>R version 3.4.3 (2017-11-30) Platform: x86_64-pc-linux-gnu (64-bit) Running under: Scientific Linux 7.4 (Nitrogen) Matrix products: default BLAS/LAPACK: /software/openblas-0.2.19-el7-x86_64/lib/libopenblas_haswellp-r0.2.19.so locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=en_US.UTF-8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] ggplot2_2.2.1 gplots_3.0.1 loaded via a namespace (and not attached): [1] Rcpp_0.12.17 knitr_1.20 magrittr_1.5 [4] munsell_0.4.3 colorspace_1.3-2 rlang_0.1.6 [7] stringr_1.2.0 plyr_1.8.4 caTools_1.17.1 [10] tools_3.4.3 grid_3.4.3 gtable_0.2.0 [13] KernSmooth_2.23-15 git2r_0.21.0 htmltools_0.3.6 [16] gtools_3.5.0 lazyeval_0.2.1 yaml_2.1.16 [19] rprojroot_1.3-2 digest_0.6.15 tibble_1.4.2 [22] bitops_1.0-6 evaluate_0.10.1 rmarkdown_1.10 [25] labeling_0.3 gdata_2.18.0 stringi_1.1.6 [28] pillar_1.1.0 compiler_3.4.3 scales_0.5.0 [31] backports_1.1.2 </code></pre> </div> <!-- Adjust MathJax settings so that all math formulae are shown using TeX fonts only; see http://docs.mathjax.org/en/latest/configuration.html. This will make the presentation more consistent at the cost of the webpage sometimes taking slightly longer to load. Note that this only works because the footer is added to webpages before the MathJax javascript. --> <script type="text/x-mathjax-config"> MathJax.Hub.Config({ "HTML-CSS": { availableFonts: ["TeX"] } }); </script> <hr> <p> This <a href="http://rmarkdown.rstudio.com">R Markdown</a> site was created with <a href="https://github.com/jdblischak/workflowr">workflowr</a> </p> <hr> <!-- To enable disqus, uncomment the section below and provide your disqus_shortname --> <!-- disqus <div id="disqus_thread"></div> <script type="text/javascript"> /* * * CONFIGURATION VARIABLES: EDIT BEFORE PASTING INTO YOUR WEBPAGE * * */ var disqus_shortname = 'rmarkdown'; // required: replace example with your forum shortname /* * * DON'T EDIT BELOW THIS LINE * * */ (function() { var dsq = document.createElement('script'); dsq.type = 'text/javascript'; dsq.async = true; dsq.src = '//' + disqus_shortname + '.disqus.com/embed.js'; (document.getElementsByTagName('head')[0] || document.getElementsByTagName('body')[0]).appendChild(dsq); })(); </script> <noscript>Please enable JavaScript to view the <a href="http://disqus.com/?ref_noscript">comments powered by Disqus.</a></noscript> <a href="http://disqus.com" class="dsq-brlink">comments powered by <span class="logo-disqus">Disqus</span></a> --> </div> <script> // add bootstrap table styles to pandoc tables function bootstrapStylePandocTables() { $('tr.header').parent('thead').parent('table').addClass('table table-condensed'); } $(document).ready(function () { bootstrapStylePandocTables(); }); </script> <!-- dynamically load mathjax for compatibility with self-contained --> <script> (function () { var script = document.createElement("script"); script.type = "text/javascript"; script.src = "https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML"; document.getElementsByTagName("head")[0].appendChild(script); })(); </script> </body> </html>