Quarto enables you to weave together content and executable code into a finished document. To learn more about Quarto see https://quarto.org.
Running Code
When you click the Render button a document will be generated that includes both content and the output of embedded code. You can embed code like this:
Registered S3 method overwritten by 'GGally':
method from
+.gg ggplot2
library(bayesplot)
This is bayesplot version 1.9.0
- Online documentation and vignettes at mc-stan.org/bayesplot
- bayesplot theme set to bayesplot::theme_default()
* Does _not_ affect other ggplot2 plots
* See ?bayesplot_theme_set for details on theme setting
theme_set(bayesplot::theme_default())SEED=87
You can add options to executable code like this
_______ _____
|__ __| 🤓 | __ \
| | _____ ___ ___| |__) |
| |/ _ \ \/ / |/ __| _ /
| | (_) > <| | (__| | \ \
|_|\___/_/\_\_|\___|_| \_\
___
| |
/ \ ____()()
/☠☠☠\ / xx
(_____) `~~~~~\_;m__m.__>o
THE SOFTWARE IS PROVIDED AS IS, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,
INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE
OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
The echo: false option disables the printing of code (only output is displayed).
library(ggplot2) # This is needed to make changes# Plot immediately#plot(weibull_fit)# Save fit and modifyggplot_fit <-plot(weibull_fit)# Transform the x axisggplot_fit +scale_x_continuous(trans ="pseudo_log")
Scale for 'x' is already present. Adding another scale for 'x', which will
replace the existing scale.
Scale for 'x' is already present. Adding another scale for 'x', which will
replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will## replace the existing scale.#Laplace Explicitly Specifiedweibull_fit_MAP <-single_dichotomous_fit(Yamamoto[,"Dose"],Yamamoto[,"Incidence"],Yamamoto[,"N"],model_type="weibull",fit_type="laplace")# Maximum Likelihood Estimationweibull_fit_MLE <-single_dichotomous_fit(Yamamoto[,"Dose"],Yamamoto[,"Incidence"],Yamamoto[,"N"],model_type="weibull",fit_type="mle")# MCMC Estimationweibull_fit_MCMC <-single_dichotomous_fit(Yamamoto[,"Dose"],Yamamoto[,"Incidence"],Yamamoto[,"N"],model_type="weibull",fit_type="mcmc")
### 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 Rtemp <-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 fittemp <-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 MLElegend("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)
