Hi !
I am new to bayesian statistics so I might misunderstood completely things here.
I am trying to fit my dose escalation preclinical data to a logistic regression and I would like to estimate the dose and associated uncertainty for a target protection probability of 70%
Data:
Dose the continuous predictor
Protected the binary outcome of survival
Model:
Protected ~ Binomial(1,\pi)
Logit(\pi) = \alpha + \beta * Dose
Am I allowed to compute a posterior for the Dose fusing the posterior distribution of the model ? Something like :
Posterior(Dose|\pi=0.7) = \frac{Logit(0.7)-Posterior(\alpha)}{Posterior(\beta)}
Or should it be specified in the brm call to be drawn during sampling ? Or am I misunderstanding something in the bayesian approach ?
The code with simulated data, and the posterior computation I described:
library(brms)
library(dplyr)
library(tidyr)
library(magrittr)
library(ggplot2)
#simulated data
n_mice=10 #mice number per Dose
Dose=c(6,36,66,96,126,156,186,216) #dose value
proba = c(0.37,0.51,0.65,0.77,0.85,0.91,0.95,0.97) #probability of protection for each dose
Protected=rbinom(n=n_mice*length(Dose),
size=1,
prob=rep(proba,each=n_mice))
data=tibble(Dose=rep(Dose,each=n_mice),
Protected)
#visualize data
data%>%
group_by(Dose)%>%
summarise(Alive=sum(Protected),
Dead=n()-Alive)%>%
ungroup()%>%
pivot_longer(cols=c('Alive','Dead'),
names_to='Status',
values_to = 'Count')%>%
mutate(Status=factor(Status,levels=c('Dead','Alive')))%>%
ggplot(aes(x=Dose,y=Status,alpha=Count,size=Count,color=Status))+
geom_point(shape=15)+
#scale_color_manual()+
theme_classic()
#Bayesian fitting
fit_bayes = brm ( Protected ~ Dose ,
data = data ,
family = bernoulli ("logit") ,
warmup= 500,
iter =4000 ,
chains=1,
cores=2)
summary(fit_bayes)
plot(fit_bayes)
#Posterior
target=0.7 # define a protection level of 70%
post=as_draws_df(fit_bayes) # extract posterior draws of model parameters
post_dose=(logit_scaled(target)-post$b_Intercept)/post$b_Dose #combination of the model parameter posteriors
hist(post_dose,
1000) # posterior distribution of the dose
mean(post_dose) # estimation of the dose for a target protection
quantile(post_dose,c(0.025,0.975)) #credible interval for the dose
Many thanks for your help!