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,6.13,8.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)
#>
#>
#> _______ _____
#> |__ __| 🤓 | __ \
#> | | _____ ___ ___| |__) |
#> | |/ _ \ \/ / |/ __| _ /
#> | | (_) > <| | (__| | \ \
#> |_|\___/_/\_\_|\___|_| \_\
#> 22.8.1.0.2
#> ___
#> | |
#> / \ ____()()
#> /☠☠☠\ / 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.
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
#> [1] "Model: Hill Distribution: Normal-NCV"
hill_fit$prior
#> $priors
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 1 1.0000000 1.0000000 -100 100
#> [2,] 1 0.0000000 1.0000000 -100 100
#> [3,] 2 0.0000000 2.0000000 0 100
#> [4,] 2 0.4700036 0.4214036 0 18
#> [5,] 2 0.0000000 1.0000000 0 100
#> [6,] 1 0.8349247 1.0000000 -18 18
#>
#> $model
#> [1] "Hill Model [normal-ncv]"
#>
#> $parameters
#> [1] "a" "b" "c" "d" "rho"
#> [6] "log(sigma^2)"
#>
#> $mean
#> [1] "hill"
#>
#> attr(,"class")
#> [1] "BMD_Bayes_continuous_model"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
#> [1] "Model: Hill Distribution: Normal"
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
#> [1] "Model: Exponential-5 Distribution: Log-Normal"
plot(exp5_fit)
Note: The above definitions are for for increasing data, `ToxicR’ will determine the direction of adversity and adjust computations accordingly.
