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