Nothing
tests/testthat/test-survival.R
In cvasi: Calibration, Validation, and Simulation of TKTD Models
#
# Fixed version of morse function, original doesn't work in morse package
# version <= 3.2.5
#
Surv.SD_Cext <- function(Cw, time, kk, kd, z, hb) {
time.prec = dplyr::lag(time, 1) ; time.prec[1] = time[1] # time[1] = tprec[1]
# Using log transfomration: log(a+b) = log(a) + log(1+b/a) may prevent the numerical issues raised by exponential
diff.int = (exp(time %*% t(kd)) + exp(time.prec %*% t(kd)) )*Cw/2 * (time-time.prec) #OK time[1]-tprec[1] = 0
D = kd * exp(-kd %*% t(time)) * t(apply(diff.int, 2, cumsum))
lambda = kk * pmax(D-z,0) + hb # the pmax function is important here for elementwise maximum with 0 and D[i,j]-z ATTENTION: pmax(0,D) != pmax(D,0)
## fix for new behaviour of dplyr::lag()
lambda.prec = dplyr::lag(lambda[1,], 1 ) ; lambda.prec[1] = lambda[1]
#lambda.prec = lag(lambda, 1 ) ; lambda.prec[1] = lambda[1]
## fix end
int.lambda = t(t((lambda + lambda.prec)/2) * (time-time.prec))
S <- exp(-t(apply(int.lambda,1,cumsum)))
return(S)
}
test_that("GUTS-RED constant exposure", {
# rjags not readily available on macOS
skip_on_os("mac")
skip_if_not_installed("morse")
# GUTS model verification as conducted in
# EFSA Scientific Opinion on TKTD models, pp. 36
# doi:10.2903/j.efsa.2018.5377
#
# Begin EFSA code: 01.GUTS-implementation-check.R
#
exp <- data.frame(seq(0,7,0.01),
c(rep(5,401), rep(0,300)))
colnames(exp) <- c("time","conc")
MF <- 1 # Multiplication factor of the exposure profile
kD <- 0.3 # unit: [time^-1]
hb <- 0 # background mortality rate [time^-1]
SOT_SD <- Surv.SD_Cext ## FIXED, morse function contains a bug preventing calculations atm
#SOT_SD <- morse:::Surv.SD_Cext
bw <- 0.5 # unit: 1/[D]
zw <- 2.5 # unit: [D]
SOT_SD_cst_simu <- SOT_SD(Cw=exp$conc, time=exp$time, kk=bw, kd=kD, z=zw, hb=hb)
SOT_IT <- morse:::Surv.IT_Cext
mw <- 2.5 # unit: [D]
beta <- 2 # dimensionless
SOT_IT_cst_simu <- SOT_IT(Cw=exp$conc, time=exp$time, kd=kD, hb=hb, alpha=mw, beta=beta)
#
# EFSA code end
#
# Compare with morse's standard prediction
GUTS_RED_IT() %>%
set_param(c(kd=kD,hb=hb,alpha=mw,beta=beta)) %>%
set_exposure(exp) %>%
survival(approx="constant") -> git
GUTS_RED_SD() %>%
set_param(c(kd=kD,hb=hb,kk=bw,z=zw)) %>%
set_exposure(exp) %>%
survival(approx="constant") -> gsd
expect_lt(mae(SOT_IT_cst_simu[1,],git$survival), 1e-5, "IT:S approx")
expect_lt(rmse(SOT_IT_cst_simu[1,],git$survival), 1e-5, "IT:S approx")
expect_lt(mae(SOT_SD_cst_simu[1,],gsd$survival), 1e-5, "SD:S approx")
expect_lt(rmse(SOT_SD_cst_simu[1,],gsd$survival), 1e-5, "SD:S approx")
# Compare with morse's ODE prediction
# Fit of our approach should be better as morse will also use the ODE solver
GUTS_RED_IT() %>%
set_param(c(kd=kD,hb=hb,alpha=mw,beta=beta)) %>%
set_exposure(exp) %>%
survival(approx="linear", hmax=NULL) -> git
GUTS_RED_SD() %>%
set_param(c(kd=kD,hb=hb,kk=bw,z=zw)) %>%
set_exposure(exp) %>%
survival(approx="linear", hmax=NULL) -> gsd
morse:::SurvSD_ode(exp$conc,exp$time,replicate=1,kk=bw,kd=kD,z=zw,hb=hb,mcmc_size=1) -> ode.sd
morse:::SurvIT_ode(exp$conc,exp$time,replicate=1,kd=kD,hb=hb,alpha=mw,beta=beta,mcmc_size=1) -> ode.it
expect_lt(mae(ode.it[,1],git$survival), 1e-10, "IT:S ode")
expect_lt(rmse(ode.it[,1],git$survival), 1e-10, "IT:S ode")
expect_lt(mae(ode.sd[,1],gsd$survival), 1e-10, "SD:D ode")
expect_lt(rmse(ode.sd[,1],gsd$survival), 1e-10, "SD:D ode")
})
test_that("GUTS-RED time-variable exposure", {
# rjags not readily available on macOS
skip_on_os("mac")
skip_if_not_installed("morse")
# GUTS model verification as conducted in
# EFSA Scientific Opinion on TKTD models, pp. 36
# doi:10.2903/j.efsa.2018.5377
#
# Begin EFSA code: 01.GUTS-implementation-check.R
#
exp <- data.frame(seq(0,14,0.01),
c(rep(0,301), rep(5,200), rep(0,900)))
colnames(exp) <- c("time","conc")
kD <- 0.3 # unit: [time^-1]
hb <- 0 # background mortality rate [time^-1]
SOT_IT <- morse:::Surv.IT_Cext
SOT_SD <- Surv.SD_Cext ## FIXED, morse function contains a bug preventing calculations atm
#SOT_SD <- morse:::Surv.SD_Cext
bw <- 0.5 # unit: 1/[D]
zw <- 2.5 # unit: [D]
mw <- 2.5 # unit: [D]
beta <- 2 # dimensionless
#
# EFSA code end
#
# Simulate GUTS-RED-IT
for (i in 2:50){
MF <- i
morse:::SurvIT_ode(exp$conc*MF,exp$time,replicate=1,kd=kD,hb=hb,alpha=mw,beta=beta,mcmc_size=1) -> ode.it
GUTS_RED_IT() %>%
set_param(c(kd=kD,hb=hb,alpha=mw,beta=beta)) %>%
set_exposure(data.frame(t=exp$time,p=exp$conc*MF)) %>%
survival(hmax=NULL) -> git
expect_lt(nmae(ode.it[,1],git$survival), 1e-5, paste0("IT:S, MF=",MF))
expect_lt(nrmse(ode.it[,1],git$survival), 1e-5, paste0("IT:S, MF=",MF))
}
# Simulate GUTS-RED-SD
for (i in 2:50){
MF <- i
morse:::SurvSD_ode(exp$conc*MF,exp$time,replicate=1,kk=bw,kd=kD,z=zw,hb=hb,mcmc_size=1) -> ode.sd
GUTS_RED_SD() %>%
set_param(c(kd=kD,hb=hb,kk=bw,z=zw)) %>%
set_exposure(data.frame(t=exp$time,p=exp$conc*MF)) %>%
survival(hmax=NULL) -> gsd
expect_lt(nmae(ode.sd[,1],gsd$survival), 1e-5, paste0("SD:S, MF=",MF))
expect_lt(nrmse(ode.sd[,1],gsd$survival), 1e-5, paste0("SD:S, MF=",MF))
}
})
test_that("GUTS-RED extreme cases", {
# rjags not readily available on macOS
skip_on_os("mac")
skip_if_not_installed("morse")
# GUTS model verification as conducted in
# EFSA Scientific Opinion on TKTD models, pp. 36
# doi:10.2903/j.efsa.2018.5377
#
# Begin EFSA code: 01.GUTS-implementation-check.R
#
MF <- 1
kD <- 0.3 # unit: [time^-1]
hb <- 0 # background mortality rate
bw <- 0.5 # unit: 1/[D]
zw <- 2.5 # unit: [D]
mw <- 2.5 # unit: [D]
beta <- 2 # dimensionless
SOT_SD <- Surv.SD_Cext ## FIXED, morse function contains a bug preventing calculations atm
#SOT_SD <- morse:::Surv.SD_Cext
SOT_IT <- morse:::Surv.IT_Cext
# Zero exposure (conc=0)
exp <- data.frame(seq(0,7,0.01), rep(0,701))
colnames(exp) <- c("time","conc")
SOT_SD_simu_0 <- SOT_SD(Cw=exp$conc, time=exp$time, kk=bw, kd=kD, z=zw, hb=hb)
SOT_IT_simu_0 <- SOT_IT(Cw=exp$conc, time=exp$time, kd=kD,hb=hb, alpha=mw, beta=beta)
#
# EFSA code end
#
morse:::SurvSD_ode(exp$conc,exp$time,replicate=1,kk=bw,kd=kD,z=zw,hb=hb,mcmc_size=1) -> ode.sd
morse:::SurvIT_ode(exp$conc,exp$time,replicate=1,kd=kD,hb=hb,alpha=mw,beta=beta,mcmc_size=1) -> ode.it
GUTS_RED_IT() %>%
set_param(c(kd=kD,hb=hb,alpha=mw,beta=beta)) %>%
set_exposure(exp) %>%
survival(approx="constant", hmax=NULL) -> git
GUTS_RED_SD() %>%
set_param(c(kd=kD,hb=hb,kk=bw,z=zw)) %>%
set_exposure(exp) %>%
survival(approx="constant", hmax=NULL) -> gsd
expect_lt(nmae(SOT_IT_simu_0[1,],git$survival), 1e-10, "IT:S approx,zero exp")
expect_lt(nrmse(SOT_IT_simu_0[1,],git$survival), 1e-10, "IT:S approx,zero exp")
expect_lt(nmae(SOT_SD_simu_0[1,],gsd$survival), 1e-10, "SD:S approx,zero exp")
expect_lt(nrmse(SOT_SD_simu_0[1,],gsd$survival), 1e-10, "SD:S approx,zero exp")
expect_lt(nmae(ode.it[,1],git$survival), 1e-10, "IT:S ODE,zero exp")
expect_lt(nrmse(ode.it[,1],git$survival), 1e-10, "IT:S ODE,zero exp")
expect_lt(nmae(ode.sd[,1],gsd$survival), 1e-10, "SD:S ODE,zero exp")
expect_lt(nrmse(ode.sd[,1],gsd$survival), 1e-10, "SD:S ODE,zero exp")
rm(git,gsd,ode.it,ode.sd,exp,SOT_IT_simu_0,SOT_SD_simu_0)
#
# Begin EFSA code: 01.GUTS-implementation-check.R
#
# High exposure over 4 days (conc=100)
exp <- data.frame(seq(0,7,0.01),
c(rep(100,401), rep(0,300)))
colnames(exp) <- c("time","conc")
SOT_SD_simu_100 <- SOT_SD(Cw=exp$conc, time=exp$time, kk=bw, kd=kD, z=zw, hb=hb)
SOT_IT_simu_100 <- SOT_IT(Cw=exp$conc, time=exp$time, kd=kD,hb=hb, alpha=mw, beta=beta)
#
# EFSA code end
#
morse:::SurvSD_ode(exp$conc,exp$time,replicate=1,kk=bw,kd=kD,z=zw,hb=hb,mcmc_size=1) -> ode.sd
morse:::SurvIT_ode(exp$conc,exp$time,replicate=1,kd=kD,hb=hb,alpha=mw,beta=beta,mcmc_size=1) -> ode.it
GUTS_RED_IT() %>%
set_param(c(kd=kD,hb=hb,alpha=mw,beta=beta)) %>%
set_exposure(exp) %>%
survival(approx="constant", hmax=NULL) -> git
GUTS_RED_SD() %>%
set_param(c(kd=kD,hb=hb,kk=bw,z=zw)) %>%
set_exposure(exp) %>%
survival(approx="constant", hmax=NULL) -> gsd
## morse's approximations have low precision
expect_lt(nmae(SOT_IT_simu_100[1,],git$survival), 1e-5, "IT:S approx,zero exp")
expect_lt(nrmse(SOT_IT_simu_100[1,],git$survival), 1e-5, "IT:S approx,zero exp")
expect_lt(nmae(SOT_SD_simu_100[1,],gsd$survival), 1e-3, "SD:S approx,zero exp")
expect_lt(nrmse(SOT_SD_simu_100[1,],gsd$survival), 1e-3, "SD:S approx,zero exp")
## morse's ODE models fare much better
expect_lt(nmae(ode.it[,1],git$survival), 1e-6, "IT:S ODE,zero exp")
expect_lt(nrmse(ode.it[,1],git$survival), 1e-6, "IT:S ODE,zero exp")
expect_lt(nmae(ode.sd[,1],gsd$survival), 1e-10, "SD:S ODE,zero exp")
expect_lt(nrmse(ode.sd[,1],gsd$survival), 1e-10, "SD:S ODE,zero exp")
rm(git,gsd,ode.it,ode.sd,exp,SOT_IT_simu_100,SOT_SD_simu_100)
#
# Begin EFSA code: 01.GUTS-implementation-check.R
#
# Only background mortality (conc=0, hb=0.05)
hb <- 0.05
exp <- data.frame(seq(0,7,0.01), rep(0,701))
colnames(exp) <- c("time","conc")
SOT_SD_simu_BG <- SOT_SD(Cw=exp$conc, time=exp$time, kk=bw,kd=kD, z=zw, hb=hb)
SOT_IT_simu_BG <- SOT_IT(Cw=exp$conc, time=exp$time, kd=kD,hb=hb, alpha=mw, beta=beta)
#
# EFSA code end
#
morse:::SurvSD_ode(exp$conc,exp$time,replicate=1,kk=bw,kd=kD,z=zw,hb=hb,mcmc_size=1) -> ode.sd
morse:::SurvIT_ode(exp$conc,exp$time,replicate=1,kd=kD,hb=hb,alpha=mw,beta=beta,mcmc_size=1) -> ode.it
GUTS_RED_IT() %>%
set_param(c(kd=kD,hb=hb,alpha=mw,beta=beta)) %>%
set_exposure(exp) %>%
survival(approx="constant") -> git
GUTS_RED_SD() %>%
set_param(c(kd=kD,hb=hb,kk=bw,z=zw)) %>%
set_exposure(exp) %>%
survival(approx="constant") -> gsd
expect_lt(nmae(SOT_IT_simu_BG[1,],git$survival), 1e-10, "IT:S approx,zero exp")
expect_lt(nrmse(SOT_IT_simu_BG[1,],git$survival), 1e-10, "IT:S approx,zero exp")
expect_lt(nmae(SOT_SD_simu_BG[1,],gsd$survival), 1e-10, "SD:S approx,zero exp")
expect_lt(nrmse(SOT_SD_simu_BG[1,],gsd$survival), 1e-10, "SD:S approx,zero exp")
expect_lt(nmae(ode.it[,1],git$survival), 1e-10, "IT:S ODE,zero exp")
expect_lt(nrmse(ode.it[,1],git$survival), 1e-10, "IT:S ODE,zero exp")
expect_lt(nmae(ode.sd[,1],gsd$survival), 1e-10, "SD:S ODE,zero exp")
expect_lt(nrmse(ode.sd[,1],gsd$survival), 1e-10, "SD:S ODE,zero exp")
rm(git,gsd,ode.it,ode.sd,exp,SOT_IT_simu_BG,SOT_SD_simu_BG)
})
test_that("GUTS-RED sensitivity analysis", {
# rjags not readily available on macOS
skip_on_os("mac")
skip_if_not_installed("morse")
# GUTS model verification as conducted in
# EFSA Scientific Opinion on TKTD models, pp. 36
# doi:10.2903/j.efsa.2018.5377
#
# Begin EFSA code: 01.GUTS-implementation-check.R
#
kD <- 0.3 # unit: [time^-1]
hb <- 0 # background mortality rate
bw <- 0.5 # unit: 1/[D]
zw <- 2.5 # unit: [D]
mw <- 2.5 # unit: [D]
beta <- 2 # dimensionless
exp <- data.frame(seq(0,5,0.01),
c(rep(0,301), rep(15,200)))
colnames(exp) <- c("time","conc")
percent <- -75:100
#
# EFSA code end
#
GUTS_RED_IT() %>%
set_param(c(kd=kD,hb=hb,alpha=mw,beta=beta)) %>%
set_exposure(exp) -> git
GUTS_RED_SD() %>%
set_param(c(kd=kD,hb=hb,kk=bw,z=zw)) %>%
set_exposure(exp) -> gsd
# Variation of kD
kD <- 0.3*(1+percent/100)
for (i in 1:length(percent)){
morse:::SurvIT_ode(exp$conc,exp$time,replicate=1,kd=kD[i],hb=hb,alpha=mw,beta=beta,mcmc_size=1) %>% tail(1) -> ode.it
git %>% set_param(c(kd=kD[i])) %>% survival(hmax=NULL) %>% tail(1) -> e.it
expect_lt(nmae(ode.it[,1],e.it$survival), 1e-5, "IT:S ODE,kd varying")
morse:::SurvSD_ode(exp$conc,exp$time,replicate=1,kk=bw,kd=kD[i],z=zw,hb=hb,mcmc_size=1) %>% tail(1) -> ode.sd
gsd %>% set_param(c(kd=kD[i])) %>% survival(hmax=NULL) %>% tail(1) -> e.sd
expect_lt(nmae(ode.sd[,1],e.sd$survival), 1e-5, "SD:S ODE,kd varying")
}
rm(i,ode.sd,ode.it,e.it,e.sd)
# Variation of bw (model SD)
kD <- 0.3
bw <- 0.5*(1+percent/100)
for (i in 1:length(percent)){
morse:::SurvSD_ode(exp$conc,exp$time,replicate=1,kk=bw[i],kd=kD,z=zw,hb=hb,mcmc_size=1) %>% tail(1) -> ode.sd
gsd %>% set_param(c(kk=bw[i])) %>% survival(hmax=NULL) %>% tail(1) -> e.sd
expect_lt(nmae(ode.sd[,1],e.sd$survival), 1e-5, "SD:S ODE,kk varying")
}
rm(i,ode.sd,e.sd)
# Variation of zw (model SD)
bw <- 0.5
zw <- 2.5*(1+percent/100)
for (i in 1:length(percent)){
morse:::SurvSD_ode(exp$conc,exp$time,replicate=1,kk=bw,kd=kD,z=zw[i],hb=hb,mcmc_size=1) %>% tail(1) -> ode.sd
gsd %>% set_param(c(z=zw[i])) %>% survival(hmax=NULL) %>% tail(1) -> e.sd
## !! reduced precision !!
expect_lt(nmae(ode.sd[,1],e.sd$survival), 5e-5, "SD:S ODE,z varying")
}
rm(i,ode.sd,e.sd)
# Variation of mw (model IT)
zw <- 2.5
mw <- 2.5*(1+percent/100)
for (i in 1:length(percent)){
morse:::SurvIT_ode(exp$conc,exp$time,replicate=1,kd=kD,hb=hb,alpha=mw[i],beta=beta,mcmc_size=1) %>% tail(1) -> ode.it
git %>% set_param(c(alpha=mw[i])) %>% survival(hmax=NULL) %>% tail(1) -> e.it
expect_lt(nmae(ode.it[,1],e.it$survival), 1e-5, "IT:S ODE,alpha varying")
}
rm(i,ode.it,e.it)
# Variation of beta (model IT)
mw <- 2.5
beta <- 2*(1+percent/100)
for (i in 1:length(percent)){
morse:::SurvIT_ode(exp$conc,exp$time,replicate=1,kd=kD,hb=hb,alpha=mw,beta=beta[i],mcmc_size=1) %>% tail(1) -> ode.it
git %>% set_param(c(beta=beta[i])) %>% survival(hmax=NULL) %>% tail(1) -> e.it
expect_lt(nmae(ode.it[,1],e.it$survival), 1e-5, "IT:S ODE,beta varying")
}
rm(i,ode.it,e.it)
})
Try the cvasi package in your browser
Any scripts or data that you put into this service are public.
cvasi documentation built on Sept. 11, 2024, 5:21 p.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.
#
# Fixed version of morse function, original doesn't work in morse package
# version <= 3.2.5
#
Surv.SD_Cext <- function(Cw, time, kk, kd, z, hb) {
time.prec = dplyr::lag(time, 1) ; time.prec[1] = time[1] # time[1] = tprec[1]
# Using log transfomration: log(a+b) = log(a) + log(1+b/a) may prevent the numerical issues raised by exponential
diff.int = (exp(time %*% t(kd)) + exp(time.prec %*% t(kd)) )*Cw/2 * (time-time.prec) #OK time[1]-tprec[1] = 0
D = kd * exp(-kd %*% t(time)) * t(apply(diff.int, 2, cumsum))
lambda = kk * pmax(D-z,0) + hb # the pmax function is important here for elementwise maximum with 0 and D[i,j]-z ATTENTION: pmax(0,D) != pmax(D,0)
## fix for new behaviour of dplyr::lag()
lambda.prec = dplyr::lag(lambda[1,], 1 ) ; lambda.prec[1] = lambda[1]
#lambda.prec = lag(lambda, 1 ) ; lambda.prec[1] = lambda[1]
## fix end
int.lambda = t(t((lambda + lambda.prec)/2) * (time-time.prec))
S <- exp(-t(apply(int.lambda,1,cumsum)))
return(S)
}
test_that("GUTS-RED constant exposure", {
# rjags not readily available on macOS
skip_on_os("mac")
skip_if_not_installed("morse")
# GUTS model verification as conducted in
# EFSA Scientific Opinion on TKTD models, pp. 36
# doi:10.2903/j.efsa.2018.5377
#
# Begin EFSA code: 01.GUTS-implementation-check.R
#
exp <- data.frame(seq(0,7,0.01),
c(rep(5,401), rep(0,300)))
colnames(exp) <- c("time","conc")
MF <- 1 # Multiplication factor of the exposure profile
kD <- 0.3 # unit: [time^-1]
hb <- 0 # background mortality rate [time^-1]
SOT_SD <- Surv.SD_Cext ## FIXED, morse function contains a bug preventing calculations atm
#SOT_SD <- morse:::Surv.SD_Cext
bw <- 0.5 # unit: 1/[D]
zw <- 2.5 # unit: [D]
SOT_SD_cst_simu <- SOT_SD(Cw=exp$conc, time=exp$time, kk=bw, kd=kD, z=zw, hb=hb)
SOT_IT <- morse:::Surv.IT_Cext
mw <- 2.5 # unit: [D]
beta <- 2 # dimensionless
SOT_IT_cst_simu <- SOT_IT(Cw=exp$conc, time=exp$time, kd=kD, hb=hb, alpha=mw, beta=beta)
#
# EFSA code end
#
# Compare with morse's standard prediction
GUTS_RED_IT() %>%
set_param(c(kd=kD,hb=hb,alpha=mw,beta=beta)) %>%
set_exposure(exp) %>%
survival(approx="constant") -> git
GUTS_RED_SD() %>%
set_param(c(kd=kD,hb=hb,kk=bw,z=zw)) %>%
set_exposure(exp) %>%
survival(approx="constant") -> gsd
expect_lt(mae(SOT_IT_cst_simu[1,],git$survival), 1e-5, "IT:S approx")
expect_lt(rmse(SOT_IT_cst_simu[1,],git$survival), 1e-5, "IT:S approx")
expect_lt(mae(SOT_SD_cst_simu[1,],gsd$survival), 1e-5, "SD:S approx")
expect_lt(rmse(SOT_SD_cst_simu[1,],gsd$survival), 1e-5, "SD:S approx")
# Compare with morse's ODE prediction
# Fit of our approach should be better as morse will also use the ODE solver
GUTS_RED_IT() %>%
set_param(c(kd=kD,hb=hb,alpha=mw,beta=beta)) %>%
set_exposure(exp) %>%
survival(approx="linear", hmax=NULL) -> git
GUTS_RED_SD() %>%
set_param(c(kd=kD,hb=hb,kk=bw,z=zw)) %>%
set_exposure(exp) %>%
survival(approx="linear", hmax=NULL) -> gsd
morse:::SurvSD_ode(exp$conc,exp$time,replicate=1,kk=bw,kd=kD,z=zw,hb=hb,mcmc_size=1) -> ode.sd
morse:::SurvIT_ode(exp$conc,exp$time,replicate=1,kd=kD,hb=hb,alpha=mw,beta=beta,mcmc_size=1) -> ode.it
expect_lt(mae(ode.it[,1],git$survival), 1e-10, "IT:S ode")
expect_lt(rmse(ode.it[,1],git$survival), 1e-10, "IT:S ode")
expect_lt(mae(ode.sd[,1],gsd$survival), 1e-10, "SD:D ode")
expect_lt(rmse(ode.sd[,1],gsd$survival), 1e-10, "SD:D ode")
})
test_that("GUTS-RED time-variable exposure", {
# rjags not readily available on macOS
skip_on_os("mac")
skip_if_not_installed("morse")
# GUTS model verification as conducted in
# EFSA Scientific Opinion on TKTD models, pp. 36
# doi:10.2903/j.efsa.2018.5377
#
# Begin EFSA code: 01.GUTS-implementation-check.R
#
exp <- data.frame(seq(0,14,0.01),
c(rep(0,301), rep(5,200), rep(0,900)))
colnames(exp) <- c("time","conc")
kD <- 0.3 # unit: [time^-1]
hb <- 0 # background mortality rate [time^-1]
SOT_IT <- morse:::Surv.IT_Cext
SOT_SD <- Surv.SD_Cext ## FIXED, morse function contains a bug preventing calculations atm
#SOT_SD <- morse:::Surv.SD_Cext
bw <- 0.5 # unit: 1/[D]
zw <- 2.5 # unit: [D]
mw <- 2.5 # unit: [D]
beta <- 2 # dimensionless
#
# EFSA code end
#
# Simulate GUTS-RED-IT
for (i in 2:50){
MF <- i
morse:::SurvIT_ode(exp$conc*MF,exp$time,replicate=1,kd=kD,hb=hb,alpha=mw,beta=beta,mcmc_size=1) -> ode.it
GUTS_RED_IT() %>%
set_param(c(kd=kD,hb=hb,alpha=mw,beta=beta)) %>%
set_exposure(data.frame(t=exp$time,p=exp$conc*MF)) %>%
survival(hmax=NULL) -> git
expect_lt(nmae(ode.it[,1],git$survival), 1e-5, paste0("IT:S, MF=",MF))
expect_lt(nrmse(ode.it[,1],git$survival), 1e-5, paste0("IT:S, MF=",MF))
}
# Simulate GUTS-RED-SD
for (i in 2:50){
MF <- i
morse:::SurvSD_ode(exp$conc*MF,exp$time,replicate=1,kk=bw,kd=kD,z=zw,hb=hb,mcmc_size=1) -> ode.sd
GUTS_RED_SD() %>%
set_param(c(kd=kD,hb=hb,kk=bw,z=zw)) %>%
set_exposure(data.frame(t=exp$time,p=exp$conc*MF)) %>%
survival(hmax=NULL) -> gsd
expect_lt(nmae(ode.sd[,1],gsd$survival), 1e-5, paste0("SD:S, MF=",MF))
expect_lt(nrmse(ode.sd[,1],gsd$survival), 1e-5, paste0("SD:S, MF=",MF))
}
})
test_that("GUTS-RED extreme cases", {
# rjags not readily available on macOS
skip_on_os("mac")
skip_if_not_installed("morse")
# GUTS model verification as conducted in
# EFSA Scientific Opinion on TKTD models, pp. 36
# doi:10.2903/j.efsa.2018.5377
#
# Begin EFSA code: 01.GUTS-implementation-check.R
#
MF <- 1
kD <- 0.3 # unit: [time^-1]
hb <- 0 # background mortality rate
bw <- 0.5 # unit: 1/[D]
zw <- 2.5 # unit: [D]
mw <- 2.5 # unit: [D]
beta <- 2 # dimensionless
SOT_SD <- Surv.SD_Cext ## FIXED, morse function contains a bug preventing calculations atm
#SOT_SD <- morse:::Surv.SD_Cext
SOT_IT <- morse:::Surv.IT_Cext
# Zero exposure (conc=0)
exp <- data.frame(seq(0,7,0.01), rep(0,701))
colnames(exp) <- c("time","conc")
SOT_SD_simu_0 <- SOT_SD(Cw=exp$conc, time=exp$time, kk=bw, kd=kD, z=zw, hb=hb)
SOT_IT_simu_0 <- SOT_IT(Cw=exp$conc, time=exp$time, kd=kD,hb=hb, alpha=mw, beta=beta)
#
# EFSA code end
#
morse:::SurvSD_ode(exp$conc,exp$time,replicate=1,kk=bw,kd=kD,z=zw,hb=hb,mcmc_size=1) -> ode.sd
morse:::SurvIT_ode(exp$conc,exp$time,replicate=1,kd=kD,hb=hb,alpha=mw,beta=beta,mcmc_size=1) -> ode.it
GUTS_RED_IT() %>%
set_param(c(kd=kD,hb=hb,alpha=mw,beta=beta)) %>%
set_exposure(exp) %>%
survival(approx="constant", hmax=NULL) -> git
GUTS_RED_SD() %>%
set_param(c(kd=kD,hb=hb,kk=bw,z=zw)) %>%
set_exposure(exp) %>%
survival(approx="constant", hmax=NULL) -> gsd
expect_lt(nmae(SOT_IT_simu_0[1,],git$survival), 1e-10, "IT:S approx,zero exp")
expect_lt(nrmse(SOT_IT_simu_0[1,],git$survival), 1e-10, "IT:S approx,zero exp")
expect_lt(nmae(SOT_SD_simu_0[1,],gsd$survival), 1e-10, "SD:S approx,zero exp")
expect_lt(nrmse(SOT_SD_simu_0[1,],gsd$survival), 1e-10, "SD:S approx,zero exp")
expect_lt(nmae(ode.it[,1],git$survival), 1e-10, "IT:S ODE,zero exp")
expect_lt(nrmse(ode.it[,1],git$survival), 1e-10, "IT:S ODE,zero exp")
expect_lt(nmae(ode.sd[,1],gsd$survival), 1e-10, "SD:S ODE,zero exp")
expect_lt(nrmse(ode.sd[,1],gsd$survival), 1e-10, "SD:S ODE,zero exp")
rm(git,gsd,ode.it,ode.sd,exp,SOT_IT_simu_0,SOT_SD_simu_0)
#
# Begin EFSA code: 01.GUTS-implementation-check.R
#
# High exposure over 4 days (conc=100)
exp <- data.frame(seq(0,7,0.01),
c(rep(100,401), rep(0,300)))
colnames(exp) <- c("time","conc")
SOT_SD_simu_100 <- SOT_SD(Cw=exp$conc, time=exp$time, kk=bw, kd=kD, z=zw, hb=hb)
SOT_IT_simu_100 <- SOT_IT(Cw=exp$conc, time=exp$time, kd=kD,hb=hb, alpha=mw, beta=beta)
#
# EFSA code end
#
morse:::SurvSD_ode(exp$conc,exp$time,replicate=1,kk=bw,kd=kD,z=zw,hb=hb,mcmc_size=1) -> ode.sd
morse:::SurvIT_ode(exp$conc,exp$time,replicate=1,kd=kD,hb=hb,alpha=mw,beta=beta,mcmc_size=1) -> ode.it
GUTS_RED_IT() %>%
set_param(c(kd=kD,hb=hb,alpha=mw,beta=beta)) %>%
set_exposure(exp) %>%
survival(approx="constant", hmax=NULL) -> git
GUTS_RED_SD() %>%
set_param(c(kd=kD,hb=hb,kk=bw,z=zw)) %>%
set_exposure(exp) %>%
survival(approx="constant", hmax=NULL) -> gsd
## morse's approximations have low precision
expect_lt(nmae(SOT_IT_simu_100[1,],git$survival), 1e-5, "IT:S approx,zero exp")
expect_lt(nrmse(SOT_IT_simu_100[1,],git$survival), 1e-5, "IT:S approx,zero exp")
expect_lt(nmae(SOT_SD_simu_100[1,],gsd$survival), 1e-3, "SD:S approx,zero exp")
expect_lt(nrmse(SOT_SD_simu_100[1,],gsd$survival), 1e-3, "SD:S approx,zero exp")
## morse's ODE models fare much better
expect_lt(nmae(ode.it[,1],git$survival), 1e-6, "IT:S ODE,zero exp")
expect_lt(nrmse(ode.it[,1],git$survival), 1e-6, "IT:S ODE,zero exp")
expect_lt(nmae(ode.sd[,1],gsd$survival), 1e-10, "SD:S ODE,zero exp")
expect_lt(nrmse(ode.sd[,1],gsd$survival), 1e-10, "SD:S ODE,zero exp")
rm(git,gsd,ode.it,ode.sd,exp,SOT_IT_simu_100,SOT_SD_simu_100)
#
# Begin EFSA code: 01.GUTS-implementation-check.R
#
# Only background mortality (conc=0, hb=0.05)
hb <- 0.05
exp <- data.frame(seq(0,7,0.01), rep(0,701))
colnames(exp) <- c("time","conc")
SOT_SD_simu_BG <- SOT_SD(Cw=exp$conc, time=exp$time, kk=bw,kd=kD, z=zw, hb=hb)
SOT_IT_simu_BG <- SOT_IT(Cw=exp$conc, time=exp$time, kd=kD,hb=hb, alpha=mw, beta=beta)
#
# EFSA code end
#
morse:::SurvSD_ode(exp$conc,exp$time,replicate=1,kk=bw,kd=kD,z=zw,hb=hb,mcmc_size=1) -> ode.sd
morse:::SurvIT_ode(exp$conc,exp$time,replicate=1,kd=kD,hb=hb,alpha=mw,beta=beta,mcmc_size=1) -> ode.it
GUTS_RED_IT() %>%
set_param(c(kd=kD,hb=hb,alpha=mw,beta=beta)) %>%
set_exposure(exp) %>%
survival(approx="constant") -> git
GUTS_RED_SD() %>%
set_param(c(kd=kD,hb=hb,kk=bw,z=zw)) %>%
set_exposure(exp) %>%
survival(approx="constant") -> gsd
expect_lt(nmae(SOT_IT_simu_BG[1,],git$survival), 1e-10, "IT:S approx,zero exp")
expect_lt(nrmse(SOT_IT_simu_BG[1,],git$survival), 1e-10, "IT:S approx,zero exp")
expect_lt(nmae(SOT_SD_simu_BG[1,],gsd$survival), 1e-10, "SD:S approx,zero exp")
expect_lt(nrmse(SOT_SD_simu_BG[1,],gsd$survival), 1e-10, "SD:S approx,zero exp")
expect_lt(nmae(ode.it[,1],git$survival), 1e-10, "IT:S ODE,zero exp")
expect_lt(nrmse(ode.it[,1],git$survival), 1e-10, "IT:S ODE,zero exp")
expect_lt(nmae(ode.sd[,1],gsd$survival), 1e-10, "SD:S ODE,zero exp")
expect_lt(nrmse(ode.sd[,1],gsd$survival), 1e-10, "SD:S ODE,zero exp")
rm(git,gsd,ode.it,ode.sd,exp,SOT_IT_simu_BG,SOT_SD_simu_BG)
})
test_that("GUTS-RED sensitivity analysis", {
# rjags not readily available on macOS
skip_on_os("mac")
skip_if_not_installed("morse")
# GUTS model verification as conducted in
# EFSA Scientific Opinion on TKTD models, pp. 36
# doi:10.2903/j.efsa.2018.5377
#
# Begin EFSA code: 01.GUTS-implementation-check.R
#
kD <- 0.3 # unit: [time^-1]
hb <- 0 # background mortality rate
bw <- 0.5 # unit: 1/[D]
zw <- 2.5 # unit: [D]
mw <- 2.5 # unit: [D]
beta <- 2 # dimensionless
exp <- data.frame(seq(0,5,0.01),
c(rep(0,301), rep(15,200)))
colnames(exp) <- c("time","conc")
percent <- -75:100
#
# EFSA code end
#
GUTS_RED_IT() %>%
set_param(c(kd=kD,hb=hb,alpha=mw,beta=beta)) %>%
set_exposure(exp) -> git
GUTS_RED_SD() %>%
set_param(c(kd=kD,hb=hb,kk=bw,z=zw)) %>%
set_exposure(exp) -> gsd
# Variation of kD
kD <- 0.3*(1+percent/100)
for (i in 1:length(percent)){
morse:::SurvIT_ode(exp$conc,exp$time,replicate=1,kd=kD[i],hb=hb,alpha=mw,beta=beta,mcmc_size=1) %>% tail(1) -> ode.it
git %>% set_param(c(kd=kD[i])) %>% survival(hmax=NULL) %>% tail(1) -> e.it
expect_lt(nmae(ode.it[,1],e.it$survival), 1e-5, "IT:S ODE,kd varying")
morse:::SurvSD_ode(exp$conc,exp$time,replicate=1,kk=bw,kd=kD[i],z=zw,hb=hb,mcmc_size=1) %>% tail(1) -> ode.sd
gsd %>% set_param(c(kd=kD[i])) %>% survival(hmax=NULL) %>% tail(1) -> e.sd
expect_lt(nmae(ode.sd[,1],e.sd$survival), 1e-5, "SD:S ODE,kd varying")
}
rm(i,ode.sd,ode.it,e.it,e.sd)
# Variation of bw (model SD)
kD <- 0.3
bw <- 0.5*(1+percent/100)
for (i in 1:length(percent)){
morse:::SurvSD_ode(exp$conc,exp$time,replicate=1,kk=bw[i],kd=kD,z=zw,hb=hb,mcmc_size=1) %>% tail(1) -> ode.sd
gsd %>% set_param(c(kk=bw[i])) %>% survival(hmax=NULL) %>% tail(1) -> e.sd
expect_lt(nmae(ode.sd[,1],e.sd$survival), 1e-5, "SD:S ODE,kk varying")
}
rm(i,ode.sd,e.sd)
# Variation of zw (model SD)
bw <- 0.5
zw <- 2.5*(1+percent/100)
for (i in 1:length(percent)){
morse:::SurvSD_ode(exp$conc,exp$time,replicate=1,kk=bw,kd=kD,z=zw[i],hb=hb,mcmc_size=1) %>% tail(1) -> ode.sd
gsd %>% set_param(c(z=zw[i])) %>% survival(hmax=NULL) %>% tail(1) -> e.sd
## !! reduced precision !!
expect_lt(nmae(ode.sd[,1],e.sd$survival), 5e-5, "SD:S ODE,z varying")
}
rm(i,ode.sd,e.sd)
# Variation of mw (model IT)
zw <- 2.5
mw <- 2.5*(1+percent/100)
for (i in 1:length(percent)){
morse:::SurvIT_ode(exp$conc,exp$time,replicate=1,kd=kD,hb=hb,alpha=mw[i],beta=beta,mcmc_size=1) %>% tail(1) -> ode.it
git %>% set_param(c(alpha=mw[i])) %>% survival(hmax=NULL) %>% tail(1) -> e.it
expect_lt(nmae(ode.it[,1],e.it$survival), 1e-5, "IT:S ODE,alpha varying")
}
rm(i,ode.it,e.it)
# Variation of beta (model IT)
mw <- 2.5
beta <- 2*(1+percent/100)
for (i in 1:length(percent)){
morse:::SurvIT_ode(exp$conc,exp$time,replicate=1,kd=kD,hb=hb,alpha=mw,beta=beta[i],mcmc_size=1) %>% tail(1) -> ode.it
git %>% set_param(c(beta=beta[i])) %>% survival(hmax=NULL) %>% tail(1) -> e.it
expect_lt(nmae(ode.it[,1],e.it$survival), 1e-5, "IT:S ODE,beta varying")
}
rm(i,ode.it,e.it)
})
Try the cvasi 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.