context("ODE")
test_that("We correctly convert week and year into a date", {
expect_that(as.character(getTimeFromWeekYear( 35, 1974 )), equals("1974-08-26 12:00:00"))
expect_that(as.character(getTimeFromWeekYear( 3, 1973 )), equals("1973-01-15 12:00:00"))
expect_that(as.character(getTimeFromWeekYear( 53, 1973 )), equals("1973-12-31 12:00:00"))
expect_that(as.character(getTimeFromWeekYear( 53, 1973 )), equals(as.character(getTimeFromWeekYear( 1, 1974 ))))
expect_that(as.character(getTimeFromWeekYear( 3, 1972 )), equals("1972-01-17 12:00:00"))
})
test_that("Runge Kutta ode solver works correctly",
{
if (!exists(".runRKF"))
{
# In certain situation this function is hidden
# (i.e. for devtools::check(), but not for devtools::test()
skip(".runRKF not available (hidden)")
}
#library("deSolve")
#predprey <- function(t,y,prs)
#{
# r <- 1.5
# a <- 1
# g <- 1
# m <- 3
# list(c(r*y[1]-a*y[1]*y[2],g*y[1]*y[2]-m*y[2]))
#}
#yini <- c(10,5)
#sol <- deSolve::ode(y=yini, times=seq(0,20,0.1),func=predprey,parms=c())
sol <- .runStep()
rkf <- .runRKF(step_size=0.1)
expect_lt( sum(abs(sol[,2:3]-rkf[,2:3]))/nrow(sol),0.05 )
})
test_that("Runge Kutta ode solver works with a large initial integration dt",
{
if (!exists(".runRKF"))
{
# In certain situation this function is hidden
# (i.e. for devtools::check(), but not for devtools::test()
skip(".runRKF not available (hidden)")
}
sol <- .runStep(step_size=0.1)
rkf <- .runRKF(step_size=0.1,h_step=1.0)
expect_lt( sum(abs(sol[,2:3]-rkf[,2:3]))/nrow(sol),0.05 )
})
test_that("Return the correct ODE results", {
if (!exists(".infection.model"))
{
# In certain situation this function is hidden
# (i.e. for devtools::check(), but not for devtools::test()
skip(".infection.model not available (hidden)")
}
data("age_sizes")
data("vaccine_calendar")
data("polymod_uk")
data("mcmcsample")
odes <- .infection.model( age_sizes[,1],
vaccine_calendar,
as.matrix(polymod_uk[mcmcsample$contact_ids+1,]),
mcmcsample$parameters$susceptibility,
mcmcsample$parameters$transmissibility,
mcmcsample$parameters$init_pop,
c(0.8,1.8) )
# Check all sums
sums <- c(119002.481089555, 633877.330662713, 1463334.37622813, 2870258.66403954, 4391396.3601774, 2786119.98841863, 496487.103476314, 2531.54735113481, 36283.6707667877, 155096.824666694, 263690.71146717, 415426.823720658, 528526.506836539, 358207.414064678)
odeSums <- colSums(odes[,2:15])
for( i in 1:length(sums) )
{
ratio <- odeSums[i]/sums[i]
expect_lt( ratio, 1.03 )
expect_gt( ratio, 0.97 )
}
expect_equal( nrow(odes), 364 )
# Check max value
maxs <- c(3132.20022124067, 16720.2726361168, 38489.4726884184, 77065.0627636948, 115970.502880153, 73496.0794524567, 13004.8663674856, 66.6225150332303, 956.806903893294, 4077.66057862465, 7073.81188424645, 10955.7871500801, 13891.6328086064, 9360.94718096195)
odeMaxs <- sapply( seq(2,15), function(x) max(odes[,x]) )
for( i in 1:length(maxs) )
{
ratio <- odeMaxs[i]/maxs[i]
expect_lt( ratio, 1.03 )
expect_gt( ratio, 0.97 )
}
# Check time of max
maxts <- c(118, 117, 118, 114, 117, 117, 118, 118, 117, 118, 114, 117, 117, 118)
odeMaxts <- sapply( seq(2,15), function(x) which.max(odes[,x]) )
for( i in 1:length(maxts) )
{
dt <- abs(odeMaxts[i]-maxts[i])
expect_lt( dt, 5 )
}
})
test_that("Return the correct ODE results, when using a different interval", {
if (!exists(".infection.model"))
{
# In certain situation this function is hidden
# (i.e. for devtools::check(), but not for devtools::test()
skip(".infection.model not available (hidden)")
}
data("age_sizes")
data("vaccine_calendar")
data("polymod_uk")
data("mcmcsample")
odes <- .infection.model( age_sizes[,1],
vaccine_calendar,
as.matrix(polymod_uk[mcmcsample$contact_ids+1,]),
mcmcsample$parameters$susceptibility,
mcmcsample$parameters$transmissibility,
mcmcsample$parameters$init_pop,
c(0.8,1.8),
7 )
# Check all sums
sums <- c(119002.481089555, 633877.330662713, 1463334.37622813, 2870258.66403954, 4391396.3601774, 2786119.98841863, 496487.103476314, 2531.54735113481, 36283.6707667877, 155096.824666694, 263690.71146717, 415426.823720658, 528526.506836539, 358207.414064678)
odeSums <- colSums(odes[,2:15])
for( i in 1:length(sums) )
{
ratio <- odeSums[i]/sums[i]
expect_lt( ratio, 1.05 )
expect_gt( ratio, 0.97 )
}
expect_equal( nrow(odes), 52 )
# Check max value
maxs <- c(3132.20022124067, 16720.2726361168, 38489.4726884184, 77065.0627636948, 115970.502880153, 73496.0794524567, 13004.8663674856, 66.6225150332303, 956.806903893294, 4077.66057862465, 7073.81188424645, 10955.7871500801, 13891.6328086064, 9360.94718096195)
odeMaxs <- sapply( seq(2,15), function(x) max(odes[,x]) )
for( i in 1:length(maxs) )
{
ratio <- odeMaxs[i]/(7*maxs[i])
expect_lt( ratio, 1.03 )
expect_gt( ratio, 0.97 )
}
# Check time of max
maxts <- c(118, 117, 118, 114, 117, 117, 118, 118, 117, 118, 114, 117, 117, 118)
odeMaxts <- sapply( seq(2,15), function(x) which.max(odes[,x]) )
for( i in 1:length(maxts) )
{
dt <- abs(odeMaxts[i]*7-maxts[i]) # Take into account different time
expect_lt( dt, 6 )
}
})
test_that("ODE works correctly with the new vaccine date vector", {
if (!exists(".infection.model"))
{
# In certain situation this function is hidden
# (i.e. for devtools::check(), but not for devtools::test()
skip(".infection.model not available (hidden)")
}
data("age_sizes")
data("polymod_uk")
data("mcmcsample")
data("vaccine_calendar")
test.vac <- vaccine_calendar
# We would take first of the month, but because we want to directly
# compare with the old way we use their dates
#test.vac[["dates"]] <- c(as.Date("1999-10-01"), as.Date("1999-11-01"),
# as.Date("1999-12-01"), as.Date("2000-01-01"),
# as.Date("2000-02-01"))
test.vac[["dates"]] <- c(as.Date("1970-10-07"), as.Date("1970-11-07"),
as.Date("1970-12-07"), as.Date("1971-01-07"),
as.Date("1971-02-07"))
test.vac[["calendar"]] <- matrix(c(test.vac[["calendar"]][1,],
test.vac[["calendar"]][32,],
test.vac[["calendar"]][62,],
test.vac[["calendar"]][93,]),ncol=21,byrow=TRUE)
odes <- .infection.model( age_sizes[,1],
test.vac,
as.matrix(polymod_uk[mcmcsample$contact_ids+1,]),
mcmcsample$parameters$susceptibility,
mcmcsample$parameters$transmissibility,
mcmcsample$parameters$init_pop,
c(0.8,1.8) )
# Check all sums
sums <- c(119002.481089555, 633877.330662713, 1463334.37622813, 2870258.66403954, 4391396.3601774, 2786119.98841863, 496487.103476314, 2531.54735113481, 36283.6707667877, 155096.824666694, 263690.71146717, 415426.823720658, 528526.506836539, 358207.414064678)
odeSums <- colSums(odes[,2:15])
for( i in 1:length(sums) )
{
ratio <- odeSums[i]/sums[i]
expect_lt( ratio, 1.03 )
expect_gt( ratio, 0.97 )
}
# Check max value
maxs <- c(3132.20022124067, 16720.2726361168, 38489.4726884184, 77065.0627636948, 115970.502880153, 73496.0794524567, 13004.8663674856, 66.6225150332303, 956.806903893294, 4077.66057862465, 7073.81188424645, 10955.7871500801, 13891.6328086064, 9360.94718096195)
odeMaxs <- sapply( seq(2,15), function(x) max(odes[,x]) )
for( i in 1:length(maxs) )
{
ratio <- odeMaxs[i]/maxs[i]
expect_lt( ratio, 1.03 )
expect_gt( ratio, 0.97 )
}
# Check time of max
maxts <- c(118, 117, 118, 114, 117, 117, 118, 118, 117, 118, 114, 117, 117, 118)
odeMaxts <- sapply( seq(2,15), function(x) which.max(odes[,x]) )
for( i in 1:length(maxts) )
{
dt <- abs(odeMaxts[i]-maxts[i])
expect_lt( dt, 5 )
}
})
test_that("We can set population sizes etc with infectionODEs", {
data("age_sizes")
data("polymod_uk")
data("mcmcsample")
data("vaccine_calendar")
test.vac <- vaccine_calendar
# We would take first of the month, but because we want to directly
# compare with the old way we use their dates
#test.vac[["dates"]] <- c(as.Date("1999-10-01"), as.Date("1999-11-01"),
# as.Date("1999-12-01"), as.Date("2000-01-01"),
# as.Date("2000-02-01"))
test.vac[["dates"]] <- c(as.Date("1970-10-07"), as.Date("1970-11-07"),
as.Date("1970-12-07"), as.Date("1971-01-07"),
as.Date("1971-02-07"))
test.vac[["calendar"]] <- matrix(c(test.vac[["calendar"]][1,],
test.vac[["calendar"]][32,],
test.vac[["calendar"]][62,],
test.vac[["calendar"]][93,]),ncol=21,byrow=TRUE)
age.groups <- stratify_by_age( age_sizes[,1],
c(1,5,15,25,45,65) )
risk.ratios <- matrix( c(
0.021, 0.055, 0.098, 0.087, 0.092, 0.183, 0.45,
0, 0, 0, 0, 0, 0, 0
), ncol=7, byrow=T )
popv <- stratify_by_risk(
age.groups, risk.ratios );
initial.infected <- rep( 10^mcmcsample$parameters$init_pop, 7 )
initial.infected <- stratify_by_risk(
initial.infected, risk.ratios );
odes <- infectionODEs( popv, initial.infected,
test.vac,
contact_matrix( as.matrix(polymod_uk[mcmcsample$contact_ids+1,]), age_sizes[,1], c(1,5,15,25,45,65) ),
mcmcsample$parameters$susceptibility,
mcmcsample$parameters$transmissibility,
c(0.8,1.8), 1 )
# Check all sums
sums <- c(119002.481089555, 633877.330662713, 1463334.37622813, 2870258.66403954, 4391396.3601774, 2786119.98841863, 496487.103476314, 2531.54735113481, 36283.6707667877, 155096.824666694, 263690.71146717, 415426.823720658, 528526.506836539, 358207.414064678)
odeSums <- colSums(odes[,2:15])
for( i in 1:length(sums) )
{
ratio <- odeSums[i]/sums[i]
expect_lt( ratio, 1.03 )
expect_gt( ratio, 0.97 )
}
# Check max value
maxs <- c(3132.20022124067, 16720.2726361168, 38489.4726884184, 77065.0627636948, 115970.502880153, 73496.0794524567, 13004.8663674856, 66.6225150332303, 956.806903893294, 4077.66057862465, 7073.81188424645, 10955.7871500801, 13891.6328086064, 9360.94718096195)
odeMaxs <- sapply( seq(2,15), function(x) max(odes[,x]) )
for( i in 1:length(maxs) )
{
ratio <- odeMaxs[i]/maxs[i]
expect_lt( ratio, 1.03 )
expect_gt( ratio, 0.97 )
}
# Check time of max
maxts <- c(118, 117, 118, 114, 117, 117, 118, 118, 117, 118, 114, 117, 117, 118)
odeMaxts <- sapply( seq(2,15), function(x) which.max(odes[,x]) )
for( i in 1:length(maxts) )
{
dt <- abs(odeMaxts[i]-maxts[i])
expect_lt( dt, 5 )
}
})
test_that("infectionODEs works correctly with different interval", {
data("age_sizes")
data("polymod_uk")
data("mcmcsample")
data("vaccine_calendar")
test.vac <- vaccine_calendar
# We would take first of the month, but because we want to directly
# compare with the old way we use their dates
#test.vac[["dates"]] <- c(as.Date("1999-10-01"), as.Date("1999-11-01"),
# as.Date("1999-12-01"), as.Date("2000-01-01"),
# as.Date("2000-02-01"))
test.vac[["dates"]] <- c(as.Date("1970-10-07"), as.Date("1970-11-07"),
as.Date("1970-12-07"), as.Date("1971-01-07"),
as.Date("1971-02-07"))
test.vac[["calendar"]] <- matrix(c(test.vac[["calendar"]][1,],
test.vac[["calendar"]][32,],
test.vac[["calendar"]][62,],
test.vac[["calendar"]][93,]),ncol=21,byrow=TRUE)
age.groups <- stratify_by_age( age_sizes[,1],
c(1,5,15,25,45,65) )
risk.ratios <- matrix( c(
0.021, 0.055, 0.098, 0.087, 0.092, 0.183, 0.45,
0, 0, 0, 0, 0, 0, 0
), ncol=7, byrow=T )
popv <- stratify_by_risk(
age.groups, risk.ratios );
initial.infected <- rep( 10^mcmcsample$parameters$init_pop, 7 )
initial.infected <- stratify_by_risk(
initial.infected, risk.ratios );
# Need to separate into age groups...
odes <- infectionODEs( popv, initial.infected,
test.vac,
contact_matrix( as.matrix(polymod_uk[mcmcsample$contact_ids+1,]), age_sizes[,1], c(1,5,15,25,45,65) ),
mcmcsample$parameters$susceptibility,
mcmcsample$parameters$transmissibility,
c(0.8,1.8), 7 )
# Check all sums
sums <- c(119002.481089555, 633877.330662713, 1463334.37622813, 2870258.66403954, 4391396.3601774, 2786119.98841863, 496487.103476314, 2531.54735113481, 36283.6707667877, 155096.824666694, 263690.71146717, 415426.823720658, 528526.506836539, 358207.414064678)
odeSums <- colSums(odes[,2:15])
for( i in 1:length(sums) )
{
ratio <- odeSums[i]/sums[i]
expect_lt( ratio, 1.03 )
expect_gt( ratio, 0.97 )
}
expect_equal( nrow(odes), 52 )
# Check max value
maxs <- c(3132.20022124067, 16720.2726361168, 38489.4726884184, 77065.0627636948, 115970.502880153, 73496.0794524567, 13004.8663674856, 66.6225150332303, 956.806903893294, 4077.66057862465, 7073.81188424645, 10955.7871500801, 13891.6328086064, 9360.94718096195)
odeMaxs <- sapply( seq(2,15), function(x) max(odes[,x]) )
for( i in 1:length(maxs) )
{
ratio <- odeMaxs[i]/(7*maxs[i])
expect_lt( ratio, 1.03 )
expect_gt( ratio, 0.96 )
}
# Check time of max
maxts <- c(118, 117, 118, 114, 117, 117, 118, 118, 117, 118, 114, 117, 117, 118)
odeMaxts <- sapply( seq(2,15), function(x) which.max(odes[,x]) )
for( i in 1:length(maxts) )
{
dt <- abs(odeMaxts[i]*7-maxts[i]) # Take into account different time
expect_lt( dt, 9 )
}
})
test_that("Second risk group works as expected", {
data("age_sizes")
data("polymod_uk")
data("mcmcsample")
data("vaccine_calendar")
test.vac <- vaccine_calendar
# We would take first of the month, but because we want to directly
# compare with the old way we use their dates
#test.vac[["dates"]] <- c(as.Date("1999-10-01"), as.Date("1999-11-01"),
# as.Date("1999-12-01"), as.Date("2000-01-01"),
# as.Date("2000-02-01"))
test.vac[["dates"]] <- c(as.Date("1970-10-07"), as.Date("1970-11-07"),
as.Date("1970-12-07"), as.Date("1971-01-07"),
as.Date("1971-02-07"))
test.vac[["calendar"]] <- matrix(c(test.vac[["calendar"]][1,],
test.vac[["calendar"]][32,],
test.vac[["calendar"]][62,],
test.vac[["calendar"]][93,]),ncol=21,byrow=TRUE)
test.vac[["calendar"]][,15:21] <- test.vac[["calendar"]][,8:14] # Use same calendar for Preg as for high risk
age.groups <- stratify_by_age( age_sizes[,1],
c(1,5,15,25,45,65) )
risk.ratios <- matrix( c(
0, 0, 0, 0, 0, 0, 0,
0.021, 0.055, 0.098, 0.087, 0.092, 0.183, 0.45
), ncol=7, byrow=T )
popv <- stratify_by_risk(
age.groups, risk.ratios );
initial.infected <- rep( 10^mcmcsample$parameters$init_pop, 7 )
initial.infected <- stratify_by_risk(
initial.infected, risk.ratios );
# Need to separate into age groups...
odes <- infectionODEs( popv, initial.infected,
test.vac,
contact_matrix( as.matrix(polymod_uk[mcmcsample$contact_ids+1,]), age_sizes[,1], c(1,5,15,25,45,65) ),
mcmcsample$parameters$susceptibility,
mcmcsample$parameters$transmissibility,
c(0.8,1.8), 1 )
# Check all sums
sums <- c(119002.481089555, 633877.330662713, 1463334.37622813, 2870258.66403954, 4391396.3601774, 2786119.98841863, 496487.103476314, 0,0,0,0,0,0,0, 2531.54735113481, 36283.6707667877, 155096.824666694, 263690.71146717, 415426.823720658, 528526.506836539, 358207.414064678)
odeSums <- colSums(odes[,2:22])
for( i in 1:length(sums) )
{
if( odeSums[i] != 0 & sums[i] != 0 )
{
ratio <- odeSums[i]/sums[i]
expect_lt( ratio, 1.03 )
expect_gt( ratio, 0.97 )
}
}
# Check max value
maxs <- c(3132.20022124067, 16720.2726361168, 38489.4726884184, 77065.0627636948, 115970.502880153, 73496.0794524567, 13004.8663674856, 0,0,0,0,0,0,0, 66.6225150332303, 956.806903893294, 4077.66057862465, 7073.81188424645, 10955.7871500801, 13891.6328086064, 9360.94718096195)
odeMaxs <- sapply( seq(2,22), function(x) max(odes[,x]) )
for( i in 1:length(maxs) )
{
if( odeMaxs[i] != 0 & maxs[i] != 0 )
{
ratio <- odeMaxs[i]/maxs[i]
expect_lt( ratio, 1.03 )
expect_gt( ratio, 0.97 )
}
}
# Check time of max
maxts <- c(118, 117, 118, 114, 117, 117, 118,
0,0,0,0,0,0,0,
118, 117, 118, 114, 117, 117, 118)
odeMaxts <- sapply( seq(2,22), function(x) which.max(odes[,x]) )
for( i in 1:length(maxts) )
{
dt <- abs(odeMaxts[i]-maxts[i])
expect_lt( dt, 5 )
}
})
test_that("infectionODEs works with less than 3 risk groups", {
data(age_sizes)
data(polymod_uk)
ag <- stratify_by_age(age_sizes$V1, limits=c(65)) # c( 43670500, 8262600 )
population <- stratify_by_risk( ag, matrix(c(0.01,0.4),nrow=1) ) # c( 43233795, 4957560, 436705, 3305040 )
ag <- c(1000,1000)
initial.infected <- stratify_by_risk( ag, matrix(c(0.01,0.4),nrow=1) )
vaccine_calendar <- list(
"efficacy" = c(0.7,0.3),
"calendar" = matrix(c(0,0.007,0.001,0.007),nrow=1),
"dates" = c(as.Date("2010-10-01"), as.Date("2011-02-01")) # begin and end date
)
# Polymod data is subdivided in seven age groups
poly <- polymod_uk[,c(1,2,3,9)]
poly[,3] <- rowSums(polymod_uk[,3:8])
contacts <- contact_matrix(as.matrix(poly), age_sizes$V1, c(65))
susceptibility <- c( 0.7, 0.3 ) # Different for different ages
transmissibility <- 0.17 # Same for all ages
infection_delays <- c( 0.8, 1.8 ) # 0.8 and 1.8 day.
odes <- infectionODEs( population, initial.infected,
vaccine_calendar, contacts,
susceptibility, transmissibility, infection_delays, 7 )
expect_equal(ncol(odes), 5)
expect_equal(format(odes$Time[1],format="%Y"),"2010" );
comp <- mapply( function(x,y) difftime(y,x)-7, odes$Time[1:(nrow(odes)-1)], odes$Time[2:nrow(odes)] )
# Up to two time diffs will be slightly more/less than 7 days due to summer time
expect_lt(sum(comp!=0),3)
expect_lt( sum(abs(comp)), 2.1/24 )
})
test_that("parameter_mapping works", {
pm <- parameter_mapping(epsilon = c(1,1,2,2,3), psi = 4, transmissibility = 5,
susceptibility = c(6,6,6,7,7,7,8), initial_infected = 9)
expect_equal(pm$epsilon, c(1,1,2,2,3))
pm <- parameter_mapping(psi = 4, transmissibility = 5,
susceptibility = c(6,6,6,7,7,7,8), initial_infected = 9)
expect_equal(pm$susceptibility, c(6,6,6,7,7,7,8))
# Test with named parameters
pars <- c(0.1,0.2,0.3,0.4,0.5,0.6, 0.7)
names(pars) <- c("epsilon", "epsilon", "psi", "transmissibiltiy", "susceptibility1", "susceptibility2", "initial")
pm <- parameter_mapping(parameters = pars)
expect_equal(pm$epsilon, c(1,2))
expect_equal(pm$susceptibility, c(5,6))
# Test with list of parameter names
pm <- parameter_mapping(parameters = c("epsilon", "epsilon", "psi", "transmissibiltiy", "susceptibility1", "susceptibility2", "initial"))
expect_equal(pm$epsilon, c(1,2))
expect_equal(pm$susceptibility, c(5,6))
# Test with vector of parameters (assume simplest model in this case)
names(pars) <- NULL
pm <- parameter_mapping(parameters = pars)
expect_equal(pm$epsilon, c(1,2))
expect_equal(pm$susceptibility, c(5,6))
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.