---
title: "ToxicR_playbook (Scrolling)"
Author: "Hyunsu Ju"
output:
flexdashboard::flex_dashboard:
theme: sandstone
vertical_layout: scroll
social: ["facebook","twitter","menu"]
source_code: embed
navbar:
- {title: "About", href: "https://github.com/NIEHS/ToxicR/wiki/Dichotomous", align: left}
---
```{r setup, include=FALSE}
library(flexdashboard)
library(rstanarm)
options(mc.cores = parallel::detectCores())
library(tidyverse)
library(GGally)
library(bayesplot)
theme_set(bayesplot::theme_default())
SEED=87
```
Page 1
===========
Row1
-----------------
### ToxicR Analysis: Yamamoto 2002 Data
```{r}
library(ToxicR)
library(ggplot2)
Yamamoto<-matrix(0,nrow=4,ncol=3)
colnames(Yamamoto) <- c("Dose","N","Incidence")
Yamamoto[,1] <- c(0,80,497,1530)
Yamamoto[,2] <- c(50,50,50,48)
Yamamoto[,3] <- c(4,2,6,24)
weibull_fit <- single_dichotomous_fit(Yamamoto[,"Dose"],Yamamoto[,"Incidence"],Yamamoto[,"N"],model_type="weibull")
```
```{r}
DT::datatable(Yamamoto,options=list(bPaginate=FALSE))
```
```{r}
weibull_fit$bmd
# Model Parameters
weibull_fit$parameters
# Covariance estimate on those parameters
weibull_fit$covariance
# The prior used
weibull_fit$prior
# Model Parameters
logit_g = weibull_fit$parameters[1]
1/(1+exp(-logit_g))
```
Page 2
===========
Row1
-----------------
### ggplot display
```{r}
library(ggplot2) # This is needed to make changes
# Plot immediately
#plot(weibull_fit)
# Save fit and modify
ggplot_fit <- plot(weibull_fit)
# Transform the x axis
ggplot_fit + scale_x_continuous(trans = "pseudo_log")
# Change labels
ggplot_fit + xlab("(Dose) ppm") + ylab("Tumor Probability") +
scale_x_continuous(trans = "pseudo_log") + ggtitle("Weibull Fit") + theme_classic()
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
#Laplace Explicitly Specified
weibull_fit_MAP <- single_dichotomous_fit(Yamamoto[,"Dose"],Yamamoto[,"Incidence"],Yamamoto[,"N"],model_type="weibull",fit_type="laplace")
# Maximum Likelihood Estimation
weibull_fit_MLE <- single_dichotomous_fit(Yamamoto[,"Dose"],Yamamoto[,"Incidence"],Yamamoto[,"N"],model_type="weibull",fit_type="mle")
# MCMC Estimation
weibull_fit_MCMC <- single_dichotomous_fit(Yamamoto[,"Dose"],Yamamoto[,"Incidence"],Yamamoto[,"N"],model_type="weibull",fit_type="mcmc")
```
Page 3
===========
Row1
-----------------
### Comparison:MCMC,MAP,and MLE
```{r}
library(ggpubr)
figure <- ggarrange(plot(weibull_fit_MAP)+ggtitle(""), plot(weibull_fit_MLE)+ggtitle(""),
plot(weibull_fit_MCMC)+ggtitle(""),
labels = c("MAP", "MLE", "MCMC"),
ncol = 2, nrow = 2)
figure
ggarrange(plot(weibull_fit_MLE)+ggtitle(""), labels=c("Weibull Fit"),
ncol =1, nrow =1)
prior1 <- create_prior_list(normprior(0,10,-100,100), lnormprior(0,1,0,100),
lnormprior(0,1,0,100));
PR1 <- create_dichotomous_prior(prior1,"weibull")
prior2 <- create_prior_list(normprior(0,10,-100,100),
lnormprior(0,0.1,0,100),
lnormprior(0,0.1,0,100));
PR2 <- create_dichotomous_prior(prior2,"weibull")
prior3 <- create_prior_list(normprior(0,10,-100,100),
lnormprior(0,10,0,100),
lnormprior(0,10,0,100));
PR3 <- create_dichotomous_prior(prior3,"weibull")
prior4<- create_prior_list(normprior(0,10,-100,100),
lnormprior(0,10000,0,10000),
lnormprior(0,10000,0,10000));
PR4 <- create_dichotomous_prior(prior4,"weibull")
weibull_fit_prior1 <- single_dichotomous_fit(Yamamoto[,"Dose"],Yamamoto[,"Incidence"],Yamamoto[,"N"],model_type="weibull",fit_type="mcmc",BMR=0.1,samples = 75000,prior=PR1)
weibull_fit_prior2 <- single_dichotomous_fit(Yamamoto[,"Dose"],Yamamoto[,"Incidence"],Yamamoto[,"N"],model_type="weibull",fit_type="mcmc",BMR=0.1,samples = 75000,prior=PR2)
weibull_fit_prior3 <- single_dichotomous_fit(Yamamoto[,"Dose"],Yamamoto[,"Incidence"],Yamamoto[,"N"],model_type="weibull",fit_type="mcmc",BMR=0.1,samples = 75000,prior=PR3)
weibull_fit_prior4 <- single_dichotomous_fit(Yamamoto[,"Dose"],Yamamoto[,"Incidence"],Yamamoto[,"N"],model_type="weibull",fit_type="mcmc", BMR=0.1,samples = 75000,prior=PR4)
## Weibull Model [binomial] Parameter Priors
## ------------------------------------------------------------------------
## Prior [logit(g)]:Normal(mu = 0.00, sd = 2.000) 1[-20.00,20.00]
## Prior [a]:Log-Normal(log-mu = 0.42, log-sd = 0.500) 1[0.00,40.00]
## Prior [b]:Log-Normal(log-mu = 0.00, log-sd = 1.500) 1[0.00,10000.00]
```
Page 4
===========
Row1
-----------------
### Prior Comparison
```{r}
library(ggpubr)
figure <- ggarrange(plot(weibull_fit_prior1)+ggtitle(""),
plot(weibull_fit_prior2)+ggtitle(""),
plot(weibull_fit_prior3)+ggtitle(""),
plot(weibull_fit_MLE)+ggtitle(""),
labels = c("Prior 1", "Prior 2","Prior 3","MLE"),
ncol = 2, nrow = 2)
# compare the curve shape for different priors
figure
```
Row2
-----------------
```{r}
library(tidyverse)
### Compare the BMDs
abmd<-NULL
abmd=rbind(weibull_fit_prior1$bmd,weibull_fit_prior2$bmd,weibull_fit_prior3$bmd,weibull_fit_MCMC$bmd,weibull_fit_MLE$bmd,weibull_fit_prior4$bmd)
abmd<-as.data.frame(abmd)
names(abmd)=c("BMD","BMDL","BMDU")
row.names(abmd)=c("prior1","prior2","prior3","default prior","MLE","prior 4")
abmd
### Extract the BMD CDF for priors
#par(mar=c(1,1,1,1))
#Extract the BMD CDF and remove possible infinities
#that may occur because of asymptotes etc.
#Note ecdf() is the emperical CDF function in R
temp <- ecdf(weibull_fit_prior1$mcmc_result$BMD_samples)
plot(temp,col=1,main="BMD prior comparison",xlab="BMD ppm",xlim=c(0,1550),lwd=2)
# black color for prior1;lnormprior(0,1,0,100)
temp <- ecdf(weibull_fit_prior2$mcmc_result$BMD_samples)
lines(temp,col=2,lwd=2)
# red color for prior2; lnormprior(0,0.1,0,100)
temp <- ecdf(weibull_fit_prior3$mcmc_result$BMD_samples)
lines(temp,col=3,lwd=2)
# green color for prior3: lnormprior(0,10,0,100)
temp <- ecdf(weibull_fit_MCMC$mcmc_result$BMD_samples)
lines(temp,col=4,lwd=2)
# blue color for the default MCMC fit; Log-Normal(log-mu = 0.42, log-sd = 0.500)
## The default MCMC fit
temp <- ecdf(weibull_fit_prior4$mcmc_result$BMD_samples)
lines(temp,col=6,lwd=2)
#For the MLE comparison
#The information is in the MLE/MAP objects object, though it is
# a little different.
#library(tidyverse)
temp <- as.data.frame(weibull_fit_MLE$bmd_dist)
names(temp) <- c("BMD","PROB")
temp <- temp %>% filter(!is.infinite(BMD))
lines(temp$BMD,temp$PROB,col=5,lwd=2)
#> quantile(Yamamoto[,1])
#0% 25% 50% 75% 100%
#0.00 60.00 288.50 755.25 1530.00
abline(v=c(60,288.5,755.25), col=c("black","blue", "red"), lty=c(1,2,3), lwd=c(1))
#abline(v=c(0.1125,0.375,1.225), col=c("black","blue", "red"), lty=c(1,2,3), lwd=c(1,2,3))
# light blue color is MLE
legend("bottomright", inset=.02, legend=c("prior 1", "prior 2","prior 3", "default prior", "MLE","prior 4"), col=1:6,
lty=1, cex=0.3, box.lty=0)
# legend(2.4, 0.2, legend=c("prior 1", "prior 2","prior 3", "default prior", "MLE"), col=1:5,
# lty=1:5, cex=0.8)
#legend("bottomright", inset=.02, legend=c("prior 1", "prior 2","prior 3", "default prior", "MLE"), col=1:5,
# lty=1, cex=0.3, box.lty=0)
#legend(3, 0.4, legend=c("prior 1", "prior 2","prior 3", "defaul prior", "MLE"),col=c("black","red","green","blue","lightblue"),
# lty=1, cex=0.6,box.lty=0)
```
Page 5
===========
Row1
-----------------
```{r}
# Laplace Explicitly Specified
weibull_fit_MCMC_1 <- single_dichotomous_fit(Yamamoto[,"Dose"],Yamamoto[,"Incidence"],Yamamoto[,"N"],model_type="weibull",fit_type="mcmc",samples = 75000)
# Maximum Likelihood Estimation
weibull_fit_MCMC_2<- single_dichotomous_fit(Yamamoto[,"Dose"],Yamamoto[,"Incidence"],Yamamoto[,"N"],model_type="weibull",fit_type="mcmc", BMR=0.05,samples = 7500)
# MCMC Estimation
weibull_fit_MCMC_3 <- single_dichotomous_fit(Yamamoto[,"Dose"],Yamamoto[,"Incidence"],Yamamoto[,"N"],model_type="weibull",fit_type="mcmc", BMR =0.01,samples = 7500)
fig_1 <- plot(weibull_fit_MCMC_1) + ggtitle("") + scale_x_continuous(trans="pseudo_log")
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
fig_2 <- plot(weibull_fit_MCMC_2) + ggtitle("") + scale_x_continuous(trans="pseudo_log")
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
fig_3 <- plot(weibull_fit_MCMC_3) + ggtitle("") + scale_x_continuous(trans="pseudo_log")
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
figure <- ggarrange(fig_1, fig_2, fig_3,
labels = c("BMR=0.1", "BMR=0.05", "BMR=0.01"),
ncol = 2, nrow = 2)
## Warning: Removed 1 rows containing missing values (geom_segment).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_point).
figure
```
