DichotomousBad
library(ToxicR)
library(knitr)
dose <- c(0,62.5,125,250)
obs <- c(1,9,8,14)
n <- c(50,50,50,50)
a <- single_dichotomous_fit(dose,obs,n,model_type = "log-probit",fit_type = "mcmc",samples = 300000)This dataset is similar to the continuous, but it is a data set used in regulatory toxicology (NIOSH and EPA TOSCA have used it), so we can't just dismiss it as an edge case. NIOSH used model averaging because the data are 'poor.' You must either pick a model that fits the data 'best' or use a simplified model. Let's look at three models.
a <- single_dichotomous_fit(dose,obs,n,model_type ="log-probit",fit_type = "mle")
b <- single_dichotomous_fit(dose,obs,n,model_type ="qlinear",fit_type = "mle")
c <- single_dichotomous_fit(dose,obs,n,model_type ="weibull",fit_type = "mle")The three fits look as follows:
summary(a)## Summary of single model fit (MLE) using ToxicR
## Model: Log-Probit
##
## BMD: 29.36 (0.00, 101.33) 90.0% CI
##
## Model GOF
## --------------------------------------------------
## X^2 P-Value
## Test: X^2 GOF 0.717 0.397
summary(b)## Summary of single model fit (MLE) using ToxicR
## Model: Quantal-Linear
##
## BMD: 78.60 (54.07, 145.06) 90.0% CI
##
## Model GOF
## --------------------------------------------------
## X^2 P-Value
## Test: X^2 GOF 2.676 0.262
summary(c)## Summary of single model fit (MLE) using ToxicR
## Model: Weibull
##
## BMD: 27.99 (0.00, 105.29) 90.0% CI
##
## Model GOF
## --------------------------------------------------
## X^2 P-Value
## Test: X^2 GOF 0.684 0.408
AIC1 <- -a$maximum + 2*3
AIC2 <- -b$maximum + 2*2
AIC3 <- -c$maximum + 2*3
AICs <- list(name=c("log-probit","quantal linear","weibull"),
AIC = c(AIC1,AIC2,AIC3))
kable(AICs,digits=3)
|
|
From a ''best model' point of view, one would argue that the quantal linear model is suitable. Now that is fine, but something is disconcerting with this choice. First, the quantal linear model's BMD is much greater than the other two, but the log-probit and Weibull model's BMD is zero. Why is this?
a <- single_dichotomous_fit(dose,obs,n,model_type ="log-probit",fit_type = "mle")
plot(a)
dose <- c(0,62.5,125,250)/250For this example, the issue is twofold. First,
alphas <- seq(0.0001,1,0.0001)
likes <- rep(0,length(alphas))
########################################
for (ii in 1:length(likes)){
alph = alphas[ii]
##############################
mle_bin_likelihood <- function(p){
prob <- logprobit_dr(p)
lprobs <- obs*log(prob)+(n-obs)*log(1-prob)
return(-sum(lprobs))
}
###############
logprobit_dr <- function(p){
returnV <- p[1] + (1-p[1])*pnorm(p[2]+alph*log(dose))
return(returnV)
}
likes[ii] <- -(optim(c(0.02,-2.2016318),mle_bin_likelihood,lower = c(0.0001,-3),
upper = c(1-1e-8,3),method="L-BFGS-B")$value)
}
#alphas[2:length(alphas)] = alphas[1:(length(alphas)-1)]
likes[1] = likes[length(alphas)]
#likes[1] = likes[length(likes)]
#alphas[1]=0.001
plot(alphas,likes,type='l',axes=FALSE,ylab="Log-Likelihood",
xlab = "Values of `b`",col="#33638DFF")
polygon(alphas,likes,col="#33638DFF")
axis(1,seq(0,5,0.33))
axis(2,seq(-87.5,-80.5,0.5))
a <- single_dichotomous_fit(dose*250,obs,n,model_type ="log-probit",fit_type = "laplace")
summary(a)## Summary of single model fit (Bayesian-MAP) using ToxicR
##
##
## BMD: 97.84 (45.04, 232.34) 90.0% CI
##
## Model GOF
## --------------------------------------------------
## X^2 P-Value
## Test: X^2 GOF 3.995 0.06
This likelihood has an actual maximum, but no significant difference exists between values of 'b' at zero vs. the maximum. As
First, the issue is there is very little knowledge about this slope parameter for a wide range of values, and a little bit of information goes a long way.
## Summary of single model fit (Bayesian-MAP) using ToxicR
##
##
## BMD: 97.84 (45.04, 232.34) 90.0% CI
##
## Model GOF
## --------------------------------------------------
## X^2 P-Value
## Test: X^2 GOF 3.995 0.06
Now our estimate of
prior <- create_prior_list(normprior(0,2,-20,20),
normprior(0,1,-40,40),
lnormprior(log(1.1),0.5, 0, 100));
p_lprior = create_dichotomous_prior(prior,"log-probit")
a <- single_dichotomous_fit(dose*250,obs,n,prior=p_lprior,fit_type = "laplace")
summary(a)## Summary of single model fit (Bayesian-MAP) using ToxicR
##
##
## BMD: 73.37 (29.90, 161.44) 90.0% CI
##
## Model GOF
## --------------------------------------------------
## X^2 P-Value
## Test: X^2 GOF 2.261 0.183
Here, we can see that the model now fits the data. So, how do we appropriately model this data? That is where model averaging comes in.
ma_fit <- ma_dichotomous_fit(dose*250,obs,n)
summary(ma_fit)## Summary of single MA BMD
##
## Individual Model BMDS
## Model BMD (BMDL, BMDU) Pr(M|Data)
## ___________________________________________________________________________________________
## Quantal-Linear 85.69 (57.73 ,155.91) 0.473
## Hill 64.34 (14.53 ,165.52) 0.156
## Probit 140.61 (103.33 ,313.86) 0.100
## Weibull 75.66 (26.64 ,181.41) 0.092
## Log-Logistic 73.78 (29.87 ,150.81) 0.079
## Gamma 98.65 (50.08 ,206.65) 0.056
## Logistic 160.96 (111.72 ,NaN) 0.031
## Log-Probit 97.84 (45.04 ,232.34) 0.012
## Multistage 74.68 (54.66 ,105.29) 0.001
## ___________________________________________________________________________________________
## Model Average BMD: 87.74 (31.26, 185.46) 90.0% CI
The lower bound is closer to this 'arbitrary' prior that was chosen just for the 'log-probit' model.
Notice how the BMDL is now lower than the quantal linear model estimate
