geo-htseq-paper/scripts/README.Rmd

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

## Install

- Download and install R (and 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}
posterior_deseq_edger <- (inv_logit_scaled(post$b_de_tooldeseq - post$b_de_tooledger))
hist(posterior_deseq_edger, breaks = 40)
```


- 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(scales::percent(ps[1,3:4]), collapse = " and ")`.

- 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)
```