geo-htseq-paper/scripts/README.Rmd

---
title: "Working with model objects"
output: github_document
---

## Install

- Download and install R [https://www.r-project.org](https://www.r-project.org) (and RStudio [https://www.rstudio.com/products/rstudio/](https://www.rstudio.com/products/rstudio/)).

- Go to R console or open scripts/README.Rmd in RStudio.

- Install these packages (this needs to be done once).

```{r}
if(!require(brms)){
  install.packages(c("brms", "here"))
  library(brms)
  library(here)
}
```

## Run

In R console OR in .Rmd file:

- Load packages to R environment.

```{r, message=FALSE}
library(brms)
library(here)
```

- Load the model object and print model summary.

```{r}
m <- readRDS(here("models/anticons_detool.rds"))
print(m)
```

- Get fitted coefficients with 95% credible intervals.

```{r}
posterior_summary(m)
```

- Get the full posterior.

```{r}
post <- posterior_samples(m)
head(post)
```

- What is the estimated difference in the proportion of anti-conservative p value histograms between DESeq2 and EdgeR?


```{r, eval=FALSE}
posterior_deseq_edger <- inv_logit_scaled(post$b_de_tooldeseq - post$b_de_tooledger)
hist(posterior_deseq_edger, breaks = 40)
```


```{r, echo=FALSE, message=FALSE}
posterior_deseq_edger <- inv_logit_scaled(post$b_de_tooldeseq - post$b_de_tooledger)
png(here("plots/posterior.png"))
hist(posterior_deseq_edger, breaks = 40)
invisible(dev.off())
```

![](plots/posterior.png)


- The posterior summary for the effect size.

```{r}
posterior_summary(posterior_deseq_edger) 
```

```{r, echo=FALSE}
ps <- posterior_summary(posterior_deseq_edger)
```

The estimated effect size is somewhere between `r paste(signif(ps[1,3:4], digits = 3) * 100, collapse = " and ")` percentage points.

- Get data from model object.

```{r}
data <- m$data
head(data)
```

- Extract stan code from model object. This is the fullest model description.

```{r}
stancode(m)
```


## This document

This README.md was generated by running:

```{r, eval=FALSE}
rmarkdown::render("scripts/README.Rmd", output_file = here::here("README.md"))
```