Skip to content

Commit

Permalink
Convert to meta report to Rmd format intented to use with knitr::knit…
Browse files Browse the repository at this point in the history
…_expand.
  • Loading branch information
AEBilgrau committed Jun 4, 2019
1 parent 77228e0 commit d93fff6
Show file tree
Hide file tree
Showing 3 changed files with 188 additions and 154 deletions.
54 changes: 35 additions & 19 deletions inst/shiny/run-reports.R
Original file line number Diff line number Diff line change
Expand Up @@ -3,33 +3,49 @@ library("rmarkdown")
library("GMCM")

# GENERAL GMCM ----
report_path <- "inst/shiny/www/report_full.R"
report_path <- "inst/shiny/www/report_full.Rmd"

# SPECIAL GMCM ----
report_path <- "inst/shiny/www/report_meta.R"

# Run directly
rmarkdown::render(report_path,
envir = new.env(parent = globalenv()))


# Step-by-step
# Spin document into Rmd
rmd <- spin(hair = report_path, knit = FALSE, format = "Rmd")

# Read yaml
yaml <- rmarkdown::yaml_front_matter(rmd)
params <- yaml$params


report_path <- "inst/shiny/www/report_meta.Rmd"

# Parameters to expand
params <- list(input_file = "../../../data/u133VsExon.csv",
header = TRUE,
sep = ";",
quote = '\"',
model_cols = c("u133", "exon"),
meta_large_vals = FALSE,
init_par = c(pie1 = 0.5, mu = 1, sigma = 1, rho = 0.5),
meta_method = "NM",
meta_max_ite = 50L,
meta_positive_rho = TRUE,
meta_IDR_thres_type = "IDR",
meta_IDR_thres = 1e-04)

# dput object to character string
cput <- function(x) {
f <- tempfile()
dput(x, file = f)
return(readLines(f))
}


expand_args <- c(list(file = report_path), params)
report_expanded <- do.call(knitr::knit_expand, expand_args)
report_expanded_path <- gsub("/report_", "/report_expanded_", report_path)

cat(report_expanded, file = report_expanded_path)

# Render the expanded document
rmarkdown::render(
input = rmd,
input = report_expanded_path,
output_options = list(self_contained = TRUE),
params = params,
envir = new.env(parent = globalenv())
)

# We can also purl
knitr::purl(rmd, output = file.path(dirname(rmd), "report_full_purl.R"), documentation = 0)
knitr::purl(report_expanded_path, output = gsub(".Rmd$", ".R", report_expanded_path),
documentation = 0)


135 changes: 0 additions & 135 deletions inst/shiny/www/report_meta.R

This file was deleted.

153 changes: 153 additions & 0 deletions inst/shiny/www/report_meta.Rmd
Original file line number Diff line number Diff line change
@@ -0,0 +1,153 @@
---
title: "Meta analysis and clustering with special GMCMs"
output: html_document
date: '`r Sys.time()`'
author: "Generated by the GMCM shiny app"
---
<!-- Old YAML params -->
<!-- params: -->
<!-- file: !r file.path(getwd(), "../../../data/freshVsFrozen.csv") -->
<!-- header: TRUE -->
<!-- sep: ";" -->
<!-- quote: '"' -->
<!-- model_cols: !r c("PreVsPost.Fresh.pval", "PreVsPost.Frozen.pval") -->
<!-- meta_large_vals: FALSE -->
<!-- init_par: !r c(pie1 = 0.5, mu = 1, sigma = 1, rho = 0.5) -->
<!-- meta_method: "NM" -->
<!-- meta_max_ite: 50 -->
<!-- meta_positive_rho: TRUE -->
<!-- meta_IDR_thres_type: "IDR" -->
<!-- meta_IDR_thres: 0.05 -->

```{r knit-int, echo=FALSE, include=FALSE}
set.seed(5828993)
```


## Initalisation
The **GMCM**^[1][1]^ package is loaded.

```{r load-packages, include=TRUE}
#install.packages("GMCM") # Uncomment to install the GMCM package
library("GMCM")
```


If **GMCM** is *not* installed, please uncomment the above line and rerun the script.

## Load data
The data is loaded and the first rows are printed

```{r load-data, include=TRUE, echo=TRUE}
ds <- read.table(file = "{{input_file}}",
header = {{header}},
sep = "{{sep}}",
quote = "\{{quote}}")
head(ds, n = 4)
```


Next, the data is subset to the columns of interest.

```{r select-data, include=TRUE, echo=TRUE}
x <- ds[, {{cput(model_cols)}}]
head(x, n = 2)
```

The rescaling is modified if the 'large values indicate strong evidence' checkbox was checked or not.
The checkbox was `{{ ifelse(meta_large_vals, "checked.", "unchecked, so we reverse ranking.") }}`

```{r preprocess-data}
{{ ifelse(meta_large_vals, "u <- Uhat(x)", "u <- Uhat(-x)") }}
```


## Initial parameters
The initial parameters, as chosen in the application, are set to

```{r show-initial-params, include=TRUE, echo=TRUE}
init_par <- {{cput(init_par)}}
```


## Model fitting
With the data loaded and defined initial parameters, the model is can be fitted.

```{r fit_model, error=TRUE}
par <- fit.meta.GMCM(u = u,
init.par = init_par,
method = "{{meta_method}}",
max.ite = {{meta_max_ite}},
positive.rho = {{meta_positive_rho}},
verbose = FALSE)
print(par)
```


I.e. we set the fitting method to `"{{meta_method}}"` with a maximum number of iterations of `{{meta_max_ite}}`.

## Meta analysis with unsupervised clustering
The estimated parameters are used to calculate the local and adjusted irreproducibility discovery rates:

```{r compute_probs}
meta_IDR_thres <- {{meta_IDR_thres}}
out <- get.IDR(x, par = par,
threshold = meta_IDR_thres) # Compute IDR
str(out)
# Use idr instead of IDR to threshold if specified
out <- get.IDR(u, par = par, threshold = meta_IDR_thres)
below <- (out[["{{meta_IDR_thres_type}}"]] < meta_IDR_thres)
out$l <- sum(below)
out$Khat <- ifelse(below, 2, 1)
```

By default, `get.IDR` computes thresholds on the adjusted IDR.
The local irreproducibility discovery rate correspond to the posterior probability of the point originating from the irreproducible component.

## Results
The classes are counted by

```{r classes_table}
table(out$Khat)
```

The results are also displayed by plotting

```{r plot_results}
plot(x, col = out$Khat, asp = 1) # Plot of raw values
plot(u, col = out$Khat, asp = 1) # Plot of copula values
z <- GMCM:::qgmm.marginal(u, theta = meta2full(par, d = ncol(u))) # Estimate latent process
plot(z, col = out$Khat, asp = 1) # Plot of estimated latent process
```

The fitted `par` object can be converted to a `theta` object which can be plotted directly:

```{r plot_theta}
plot(meta2full(par, d = ncol(u)))
```


### Session information
This report was generated using **rmarkdown**^[2][2]^ and **knitr**^[3][3]^ under the session
given below. The report utilizes [parameterized reports][2] and [`knitr::spin`][3].


```{r session-info}
sessionInfo()
```


### References
Please cite the **GMCM** paper^[1][1]^ if you use the package or shiny app.

```{r citation, echo=FALSE, results='asis'}
cites <- lapply(c("GMCM", "knitr", "rmarkdown"), citation)
fmt_cites <- unlist(lapply(cites, format, style = "text"))
cat(paste0(" ", seq_along(fmt_cites), ". ", fmt_cites, "\n"))
```

[1]: http://doi.org/10.18637/jss.v070.i02
[2]: https://bookdown.org/yihui/rmarkdown/parameterized-reports.html
[3]: https://yihui.name/knitr/demo/stitch/

0 comments on commit d93fff6

Please sign in to comment.