Nothing
ToxicR: Continuous Analysis
In ToxicR: Analyzing Toxicology Dose-Response Data
knitr::opts_chunk$set(
echo = TRUE,
collapse = TRUE,
comment = "#>"
)
Basics of Toxic R Continuous run.
This file shows ToxicR for dichotomous benchmark dose analyses using both single model and multiple model fits. The first thing we need to do is load the data. From an excel file this can be done using the following R commands:
cont_data <- matrix(0,nrow=5,ncol=4)
colnames(cont_data) <- c("Dose","Mean","N","SD")
cont_data[,1] <- c(0,50,100,200,400)
cont_data[,2] <- c(5.26,5.76,7.13,9.24,9.23)
cont_data[,3] <- c(20,20,20,20,20)
cont_data[,4]<- c(2.23,1.47,2.47,2.24,1.56)
Y <- cont_data[,2:4]
Hill fit
In what follows, I describe fitting the models using summary statistics. One can also fit these models using the original data with no change to the function call except the Y is a $n\times 1$ now one column vector, and the dose vector is also $n\times 1$ and corresponds to each entry in the data vector.
When summary data are used, ToxicR expects $Y$ to be a $n\times 3$ matrix. The first column of $Y$ is the the mean, the second colunn is the number of units on test, and the third column
is the observed standard deviation. In what follows, lets look at the Hill model fit, where the
Hill model is
\begin{align}
f(dose) &= a + \frac{b\times dose}{c^d + [dose]^d}
\end{align}
library(ToxicR)
library(ggplot2)
hill_fit <- single_continuous_fit(cont_data[,"Dose"],Y,
model_type="hill")
Object Extraction
The last line fits the Hill Bayesian Laplace model and puts all of the information one needs in the `hill_fit' object variable. This variable has the same structure as its dichotomous counterpart. By default the procedure assumes normal variance proportional to the mean.
To see what model you fit, simply type the following line:
hill_fit$full_model
hill_fit$prior
Here, we see the Hill model is fit using the Normal-NCV, or normal variance proportional
to the mean. For the Hill model, you can also use 'normal' as a distribution option.
hill_fit <- single_continuous_fit(cont_data[,"Dose"],
cbind(cont_data[,"Mean"],cont_data[,"N"],cont_data[,"SD"]),
model_type="hill",distribution = "normal",
fit_type = "mcmc")
hill_fit$full_model
plot(hill_fit)
Notice how the `distribution' option controls the distribution. For the "hill", "power", and
"polynomial" DR models you can choose either "normal" or "normal-ncv." or the "exp-3" and "exp-5" models you can also choose "lognormal."
exp5_fit <- single_continuous_fit(cont_data[,"Dose"],Y,
model_type="exp-5",distribution = "lognormal",fit_type="laplace")
exp5_fit$full_model
plot(exp5_fit)
Benchmark Dose Types
Absolute Deviation
There are also other types of BMDs you can choose. Here one can specify absolute deviation as
'abs,' which solves
\begin{align}
\mid f(\mbox{BMD}) - f(0) \mid = \mbox{BMR},
\end{align}
where $BMR$ is a specific cutoff value. In the example below, the abosolute difference
of $\mbox{BMR} = 2.$
hill_sd_fit <- single_continuous_fit(cont_data[,"Dose"],Y,
model_type="hill",distribution = "normal-ncv", fit_type="mcmc",
BMD_TYPE="abs",BMR = 2)
hill_sd_fit$full_model
plot(hill_sd_fit)
Standard Deviation
The standard deviation approach is the default approach, and this is the value that solves
\begin{align}
\mid f(\mbox{BMD}) - f(0) \mid = \mbox{BMR} \sigma_0,
\end{align}
and in this definition the $BMR \geq 0$ is the number of standard deviations the
mean changes from no-exposure. In the example below, the $\mbox{BMR} = 1.5.$
hill_sd_fit <- single_continuous_fit(cont_data[,"Dose"],
cbind(cont_data[,"Mean"],cont_data[,"N"],cont_data[,"SD"]),
model_type="hill",distribution = "normal-ncv",fit_type="laplace",
BMD_TYPE="sd",BMR = 1.5)
hill_sd_fit$full_model
plot(hill_sd_fit)
Hybrid Deviation
The hybrid approach is a probabilistic approach that mimics risk for dichotomous
data. Here, the BMD is the value that solves
\begin{align}
\frac{ Pr[y < Y_0 \mid \mbox{BMD}] - Pr[y > Y_0 \mid 0]}{1-Pr[y > Y_0 \mid 0]} = \mbox{BMR} .
\end{align}
Here, $Pr[y > Y_0 \mid 0]$ is the probability that an 'adverse' response is
observed at background dose, and $0 < \mbox{BMR} < 1$ is the increase in probability
of seeing an adverse response at dose $\mbox{BMD}$. For this approach $Pr[y < Y_0 \mid 0]$
must be specified using the "point_p" option. This option is only used when
"hybrid" is specified. The following shows a Hill fit, where the BMD with only a $2.5\%$ chance
of being observed at no dose, but has a $12.25\%$ probability of being observed
at dose $\mbox{BMD}.$ Note: $\frac{0.1225 - 0.025}{0.975} = 0.1.$ For more information
on this approach see Crump [@crump1995calculation].
hill_hybrid_fit <- single_continuous_fit(cont_data[,"Dose"],
cbind(cont_data[,"Mean"],cont_data[,"N"],cont_data[,"SD"]),
model_type="hill",distribution = "normal-ncv",fit_type="mcmc",
BMD_TYPE="hybrid",point_p = 0.025,BMR = 0.1)
hill_hybrid_fit$full_model
plot(hill_hybrid_fit)
Relative Deviation
The relatived deviation approach [@slob2002dose] is the value that solves
\begin{align}
\mid f(\mbox{BMD}) - f(0) \mid = [1 + \mbox{BMR}] \times f(0),
\end{align}
where $\mbox{BMR} > 0$ is the value of increase/decrease from the background.
hill_rd_fit <- single_continuous_fit(cont_data[,"Dose"],
cbind(cont_data[,"Mean"],cont_data[,"N"],cont_data[,"SD"]),
model_type="hill",distribution = "normal-ncv", fit_type="mcmc",
BMD_TYPE="rel",BMR = 0.1,samples = 50000)
hill_rd_fit$full_model
plot(hill_rd_fit)
Note: The above definitions are for for increasing data, `ToxicR' will determine
the direction of adversity and adjust computations accordingly.
Prior Specification
Like the dichotomous case you can modify the priors, to see what prior you want to modify
Again, this is easy we can type the following code to see what our current prior is:
hill_sd_fit$prior
prior <- create_prior_list(normprior(0,1,-100,100),
normprior(0,1,-1e4,1e4),
lnormprior(0,1, 0, 100),
lnormprior(log(1),0.4215,0,18),
lnormprior(0,1,0,100),
normprior(0, 10,-100,100));
p_hill_ncv = create_continuous_prior(prior,"hill","normal-ncv")
prior <- create_prior_list(normprior(0,1,-100,100),
normprior(0,1,-1e4,1e4),
lnormprior(0,1, 0, 100),
lnormprior(log(1),0.4215,0,18),
normprior(0, 10,-100,100));
p_hill_norm = create_continuous_prior(prior,"hill","normal")
When specifying a prior, it is slightly different than the dichotomous case. Here,
we need to specify what distribution we are using too. This adds one or two parameters
depending on the parameters. All models have a parameter called $\log(\sigma^2),$ which
is the variance term, and it is ALWAYS the last parameter in the model, and the
prior is placed over the log of this value. For the 'normal-NCV'
option, the variance is $\sigma^2 f(dose)^\rho,$ i.e., it has an additional parameter $\rho.$
The parameter $\rho$ must have a prior, which is $\mbox{log-Normal}(0,1)$ in the above
code.
hill_sd_a_fit <- single_continuous_fit(cont_data[,"Dose"],
cbind(cont_data[,"Mean"],cont_data[,"N"],cont_data[,"SD"]),
prior=p_hill_ncv,
fit_type="laplace",
BMD_TYPE="sd",BMR = 1.5)
hill_sd_b_fit <- single_continuous_fit(cont_data[,"Dose"],
cbind(cont_data[,"Mean"],cont_data[,"N"],cont_data[,"SD"]),
prior=p_hill_norm,
fit_type="laplace",
BMD_TYPE="sd",BMR = 1.5)
library(ggpubr)
figure <- ggarrange(plot(hill_sd_a_fit)+ggtitle(""),
plot(hill_sd_b_fit)+ggtitle(""),
labels = c("Prior NCV", "Prior Normal"),
ncol = 1, nrow = 2)
figure
In the above code, we didn't need to specify the 'model_type' or 'distribution' because
this is done implicitly when you specify the prior. The above code works for the hill
model, lets try it for a polynomial model.
poly_sd <- single_continuous_fit(cont_data[,"Dose"],Y,
distribution="normal",model_type="polynomial",
degree = 4,
fit_type="laplace",
BMD_TYPE="sd",BMR = 0.5)
plot(poly_sd)
Model Averaging
prior <- create_prior_list(normprior(0,1,-100,100),
normprior(0,1,-1e4,1e4),
lnormprior(0,1, 0, 100),
lnormprior(log(1),0.4215,0,18),
lnormprior(0,1,0,100),
normprior(0, 10,-100,100));
p_hill_ncv = create_continuous_prior(prior,"hill","normal-ncv")
prior <- create_prior_list(normprior(0,1,-100,100),
normprior(0,1,-1e4,1e4),
lnormprior(0,1, 0, 100),
lnormprior(log(1),0.4215,0,18),
normprior(0, 10,-100,100));
p_hill_norm = create_continuous_prior(prior,"hill","normal")
prior <- create_prior_list(normprior(0,1,-100,100),
normprior(0,1,-1e4,1e4),
lnormprior(0,1, 0, 100),
lnormprior(log(1),0.4215,0,18),
normprior(0, 10,-100,100));
p_exp5_norm = create_continuous_prior(prior,"exp-5","normal")
prior <- create_prior_list(normprior(0,1,-100,100),
normprior(0,1,-1e4,1e4),
lnormprior(log(1),0.4215,0,18),
normprior(0, 10,-100,100));
p_power_norm = create_continuous_prior(prior,"power","normal")
prior <- create_prior_list(normprior(0,1,-100,100),
normprior(0,1,-1e4,1e4),
lnormprior(log(1),0.4215,0,18),
normprior(0, 10,-100,100));
p_exp3_norm = create_continuous_prior(prior,"exp-3","normal")
prior_list = list(p_exp3_norm,p_hill_norm,p_exp5_norm,p_power_norm)
With a list of priors one can then run a model average based upon
a user specified model space. Here, all of the same benchmark dose options
of 'single_continuous_fit' apply.
ma_sd_mcmc_2 <- ma_continuous_fit(cont_data[,"Dose"],Y, fit_type= "laplace",
BMD_TYPE="sd",BMR = 0.5,samples = 50000,model_list=prior_list)
plot(ma_sd_mcmc_2)
References
Wout Slob, Dose-Response Modeling of Continuous Endpoints, Toxicological Sciences, Volume 66, Issue 2, April 2002, Pages 298–312
Crump, Kenny S., Calculation of benchmark doses from continuous data, Risk Analysis, Volume 15, Issue 1, February 1995, Pages 79-89
Try the ToxicR package in your browser
Any scripts or data that you put into this service are public.
ToxicR documentation built on Dec. 28, 2022, 3:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( echo = TRUE, collapse = TRUE, comment = "#>" )
Basics of Toxic R Continuous run.
This file shows ToxicR for dichotomous benchmark dose analyses using both single model and multiple model fits. The first thing we need to do is load the data. From an excel file this can be done using the following R commands:
cont_data <- matrix(0,nrow=5,ncol=4) colnames(cont_data) <- c("Dose","Mean","N","SD") cont_data[,1] <- c(0,50,100,200,400) cont_data[,2] <- c(5.26,5.76,7.13,9.24,9.23) cont_data[,3] <- c(20,20,20,20,20) cont_data[,4]<- c(2.23,1.47,2.47,2.24,1.56) Y <- cont_data[,2:4]
Hill fit
In what follows, I describe fitting the models using summary statistics. One can also fit these models using the original data with no change to the function call except the Y is a $n\times 1$ now one column vector, and the dose vector is also $n\times 1$ and corresponds to each entry in the data vector.
When summary data are used, ToxicR expects $Y$ to be a $n\times 3$ matrix. The first column of $Y$ is the the mean, the second colunn is the number of units on test, and the third column is the observed standard deviation. In what follows, lets look at the Hill model fit, where the Hill model is \begin{align} f(dose) &= a + \frac{b\times dose}{c^d + [dose]^d} \end{align}
library(ToxicR) library(ggplot2) hill_fit <- single_continuous_fit(cont_data[,"Dose"],Y, model_type="hill")
Object Extraction
The last line fits the Hill Bayesian Laplace model and puts all of the information one needs in the `hill_fit' object variable. This variable has the same structure as its dichotomous counterpart. By default the procedure assumes normal variance proportional to the mean. To see what model you fit, simply type the following line:
hill_fit$full_model hill_fit$prior
Here, we see the Hill model is fit using the Normal-NCV, or normal variance proportional to the mean. For the Hill model, you can also use 'normal' as a distribution option.
hill_fit <- single_continuous_fit(cont_data[,"Dose"], cbind(cont_data[,"Mean"],cont_data[,"N"],cont_data[,"SD"]), model_type="hill",distribution = "normal", fit_type = "mcmc") hill_fit$full_model plot(hill_fit)
Notice how the `distribution' option controls the distribution. For the "hill", "power", and "polynomial" DR models you can choose either "normal" or "normal-ncv." or the "exp-3" and "exp-5" models you can also choose "lognormal."
exp5_fit <- single_continuous_fit(cont_data[,"Dose"],Y, model_type="exp-5",distribution = "lognormal",fit_type="laplace") exp5_fit$full_model plot(exp5_fit)
Benchmark Dose Types
Absolute Deviation
There are also other types of BMDs you can choose. Here one can specify absolute deviation as 'abs,' which solves \begin{align} \mid f(\mbox{BMD}) - f(0) \mid = \mbox{BMR}, \end{align} where $BMR$ is a specific cutoff value. In the example below, the abosolute difference of $\mbox{BMR} = 2.$
hill_sd_fit <- single_continuous_fit(cont_data[,"Dose"],Y, model_type="hill",distribution = "normal-ncv", fit_type="mcmc", BMD_TYPE="abs",BMR = 2) hill_sd_fit$full_model plot(hill_sd_fit)
Standard Deviation
The standard deviation approach is the default approach, and this is the value that solves \begin{align} \mid f(\mbox{BMD}) - f(0) \mid = \mbox{BMR} \sigma_0, \end{align} and in this definition the $BMR \geq 0$ is the number of standard deviations the mean changes from no-exposure. In the example below, the $\mbox{BMR} = 1.5.$
hill_sd_fit <- single_continuous_fit(cont_data[,"Dose"], cbind(cont_data[,"Mean"],cont_data[,"N"],cont_data[,"SD"]), model_type="hill",distribution = "normal-ncv",fit_type="laplace", BMD_TYPE="sd",BMR = 1.5) hill_sd_fit$full_model plot(hill_sd_fit)
Hybrid Deviation
The hybrid approach is a probabilistic approach that mimics risk for dichotomous data. Here, the BMD is the value that solves \begin{align} \frac{ Pr[y < Y_0 \mid \mbox{BMD}] - Pr[y > Y_0 \mid 0]}{1-Pr[y > Y_0 \mid 0]} = \mbox{BMR} . \end{align} Here, $Pr[y > Y_0 \mid 0]$ is the probability that an 'adverse' response is observed at background dose, and $0 < \mbox{BMR} < 1$ is the increase in probability of seeing an adverse response at dose $\mbox{BMD}$. For this approach $Pr[y < Y_0 \mid 0]$ must be specified using the "point_p" option. This option is only used when "hybrid" is specified. The following shows a Hill fit, where the BMD with only a $2.5\%$ chance of being observed at no dose, but has a $12.25\%$ probability of being observed at dose $\mbox{BMD}.$ Note: $\frac{0.1225 - 0.025}{0.975} = 0.1.$ For more information on this approach see Crump [@crump1995calculation].
hill_hybrid_fit <- single_continuous_fit(cont_data[,"Dose"], cbind(cont_data[,"Mean"],cont_data[,"N"],cont_data[,"SD"]), model_type="hill",distribution = "normal-ncv",fit_type="mcmc", BMD_TYPE="hybrid",point_p = 0.025,BMR = 0.1) hill_hybrid_fit$full_model plot(hill_hybrid_fit)
Relative Deviation
The relatived deviation approach [@slob2002dose] is the value that solves \begin{align} \mid f(\mbox{BMD}) - f(0) \mid = [1 + \mbox{BMR}] \times f(0), \end{align}
where $\mbox{BMR} > 0$ is the value of increase/decrease from the background.
hill_rd_fit <- single_continuous_fit(cont_data[,"Dose"], cbind(cont_data[,"Mean"],cont_data[,"N"],cont_data[,"SD"]), model_type="hill",distribution = "normal-ncv", fit_type="mcmc", BMD_TYPE="rel",BMR = 0.1,samples = 50000) hill_rd_fit$full_model plot(hill_rd_fit)
Note: The above definitions are for for increasing data, `ToxicR' will determine the direction of adversity and adjust computations accordingly.
Prior Specification
Like the dichotomous case you can modify the priors, to see what prior you want to modify Again, this is easy we can type the following code to see what our current prior is:
hill_sd_fit$prior prior <- create_prior_list(normprior(0,1,-100,100), normprior(0,1,-1e4,1e4), lnormprior(0,1, 0, 100), lnormprior(log(1),0.4215,0,18), lnormprior(0,1,0,100), normprior(0, 10,-100,100)); p_hill_ncv = create_continuous_prior(prior,"hill","normal-ncv") prior <- create_prior_list(normprior(0,1,-100,100), normprior(0,1,-1e4,1e4), lnormprior(0,1, 0, 100), lnormprior(log(1),0.4215,0,18), normprior(0, 10,-100,100)); p_hill_norm = create_continuous_prior(prior,"hill","normal")
When specifying a prior, it is slightly different than the dichotomous case. Here,
we need to specify what distribution we are using too. This adds one or two parameters
depending on the parameters. All models have a parameter called $\log(\sigma^2),$ which
is the variance term, and it is ALWAYS the last parameter in the model, and the
prior is placed over the log of this value. For the 'normal-NCV'
option, the variance is $\sigma^2 f(dose)^\rho,$ i.e., it has an additional parameter $\rho.$
The parameter $\rho$ must have a prior, which is $\mbox{log-Normal}(0,1)$ in the above
code.
hill_sd_a_fit <- single_continuous_fit(cont_data[,"Dose"], cbind(cont_data[,"Mean"],cont_data[,"N"],cont_data[,"SD"]), prior=p_hill_ncv, fit_type="laplace", BMD_TYPE="sd",BMR = 1.5) hill_sd_b_fit <- single_continuous_fit(cont_data[,"Dose"], cbind(cont_data[,"Mean"],cont_data[,"N"],cont_data[,"SD"]), prior=p_hill_norm, fit_type="laplace", BMD_TYPE="sd",BMR = 1.5) library(ggpubr) figure <- ggarrange(plot(hill_sd_a_fit)+ggtitle(""), plot(hill_sd_b_fit)+ggtitle(""), labels = c("Prior NCV", "Prior Normal"), ncol = 1, nrow = 2) figure
In the above code, we didn't need to specify the 'model_type' or 'distribution' because this is done implicitly when you specify the prior. The above code works for the hill model, lets try it for a polynomial model.
poly_sd <- single_continuous_fit(cont_data[,"Dose"],Y, distribution="normal",model_type="polynomial", degree = 4, fit_type="laplace", BMD_TYPE="sd",BMR = 0.5) plot(poly_sd)
Model Averaging
prior <- create_prior_list(normprior(0,1,-100,100), normprior(0,1,-1e4,1e4), lnormprior(0,1, 0, 100), lnormprior(log(1),0.4215,0,18), lnormprior(0,1,0,100), normprior(0, 10,-100,100)); p_hill_ncv = create_continuous_prior(prior,"hill","normal-ncv") prior <- create_prior_list(normprior(0,1,-100,100), normprior(0,1,-1e4,1e4), lnormprior(0,1, 0, 100), lnormprior(log(1),0.4215,0,18), normprior(0, 10,-100,100)); p_hill_norm = create_continuous_prior(prior,"hill","normal") prior <- create_prior_list(normprior(0,1,-100,100), normprior(0,1,-1e4,1e4), lnormprior(0,1, 0, 100), lnormprior(log(1),0.4215,0,18), normprior(0, 10,-100,100)); p_exp5_norm = create_continuous_prior(prior,"exp-5","normal") prior <- create_prior_list(normprior(0,1,-100,100), normprior(0,1,-1e4,1e4), lnormprior(log(1),0.4215,0,18), normprior(0, 10,-100,100)); p_power_norm = create_continuous_prior(prior,"power","normal") prior <- create_prior_list(normprior(0,1,-100,100), normprior(0,1,-1e4,1e4), lnormprior(log(1),0.4215,0,18), normprior(0, 10,-100,100)); p_exp3_norm = create_continuous_prior(prior,"exp-3","normal") prior_list = list(p_exp3_norm,p_hill_norm,p_exp5_norm,p_power_norm)
With a list of priors one can then run a model average based upon a user specified model space. Here, all of the same benchmark dose options of 'single_continuous_fit' apply.
ma_sd_mcmc_2 <- ma_continuous_fit(cont_data[,"Dose"],Y, fit_type= "laplace", BMD_TYPE="sd",BMR = 0.5,samples = 50000,model_list=prior_list) plot(ma_sd_mcmc_2)
References
Wout Slob, Dose-Response Modeling of Continuous Endpoints, Toxicological Sciences, Volume 66, Issue 2, April 2002, Pages 298–312
Crump, Kenny S., Calculation of benchmark doses from continuous data, Risk Analysis, Volume 15, Issue 1, February 1995, Pages 79-89
Try the ToxicR package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.