testthat::context("SIR")
testthat::test_that("SIR w/ coh=ref=popEpi::sire", {
## don't skip on CRAN
sire2 <- copy(popEpi::sire)
sire2[, agegroup := cut(dg_age, breaks = c(0:17*5, Inf))]
levels(sire2$agegroup) <- 1:18
sire2[, ex_y := year(ex_date)]
sire2[, dur := as.integer(ex_date-dg_date)/365.242199]
ltre <- ltable(sire2[status != 0], c("ex_y", "agegroup"),
expr = list(obs=.N, pyrs=sum(dur)))
setDT(ltre)
ltre[is.na(pyrs), pyrs := 0]
sibr2 <- copy(popEpi::sibr)
sibr2[, agegroup := cut(dg_age, breaks = c(0:17*5, Inf))]
levels(sibr2$agegroup) <- 1:18
sibr2[, ex_y := year(ex_date)]
sibr2[, dur := as.integer(ex_date-dg_date)/365.242199]
ltbr <- ltable(sibr2[status != 0], c("ex_y", "agegroup"),
expr = list(obs=.N, pyrs=sum(dur)))
setDT(ltbr)
ltbr[is.na(pyrs), pyrs := 0]
sl <- sir(coh.data=ltre, coh.obs="obs", coh.pyrs="pyrs",
ref.data=ltre, ref.obs="obs", ref.pyrs="pyrs",
adjust= c("agegroup","ex_y"))
plot(sl)
## SIR w/ coh=ref=popEpi::sire
## don't skip on CRAN
testthat::expect_equal(sl$sir, 1)
testthat::expect_equal(sl$pyrs, 13783.81, tolerance=0.01)
testthat::expect_equal(sl$expected, 4595)
testthat::expect_equal(sl$observed, 4595)
sl <- sir(coh.data=ltre, coh.obs="obs", coh.pyrs="pyrs",
ref.data=ltbr, ref.obs="obs", ref.pyrs="pyrs",
adjust= c("agegroup","ex_y"))
## SIR w/ coh=ref=popEpi::sire"
testthat::expect_equal(sl$sir, 1.39, tolerance=0.01)
testthat::expect_equal(sl$pyrs, 13783.81, tolerance=0.01)
testthat::expect_equal(sl$expected, 3305.04, tolerance=0.01)
testthat::expect_equal(sl$observed, 4595)
})
# SIR mstate, subset + lexpand aggre ---------------------------------------
# same model
c <- lexpand( popEpi::sire[dg_date<ex_date,], status = status, birth = bi_date, exit = ex_date, entry = dg_date,
breaks = list(per = 1990:2013, age = 0:100, fot = c(0,10,20,Inf)),
aggre = list(fot, agegroup = age, year = per, sex) )
# different models
c2 <- lexpand( popEpi::sire[dg_date<ex_date,], status = status, birth = bi_date, exit = ex_date, entry = dg_date,
breaks = list(per = 1990:2010, age = 0:100, fot = c(0,10,20,Inf)),
aggre = list(fot, agegroup = age, year = per, sex) )
testthat::test_that("misc things work", {
ms1 <- sir( coh.data = c, coh.obs = c('from0to1'), coh.pyrs = 'pyrs',
ref.data = data.table(popEpi::popmort), ref.rate = 'haz',
adjust = c('agegroup','year','sex'),
print = list(age.cat = cut(agegroup, c(0,75,120)), per = cut(year, c(1990,2010,2014), dig.lab = 8) ))
plot(ms1)
# only one print level with no variation
testthat::expect_message( ms2 <- sir( coh.data = c, coh.obs = c('from0to1'), coh.pyrs = 'pyrs', conf.type= 'profile',
ref.data = data.table(popEpi::popmort), ref.rate = 'haz',
adjust = c('agegroup','year','sex'),
print = list(sex)) )
testthat::expect_equal(attributes(ms2)$sir.meta$conf.type, 'profile')
# coef and confint ---
testthat::expect_equal(as.numeric(coef(ms1)), ms1$sir)
testthat::expect_equal(data.table(confint(ms1)), data.table('2.5 %' = ms1$sir.lo, '97.5 %' = ms1$sir.hi), tolerance = 0.001)
})
testthat::test_that("confidence intervals an other options works", {
# CI's
ci1 <- sir( coh.data = c, coh.obs = c('from0to1'), coh.pyrs = 'pyrs', conf.type = 'profile',
ref.data = data.table(popEpi::popmort), ref.rate = 'haz',
adjust = c('agegroup','year','sex'),
print = list(FOT = fot, age.cat = cut(agegroup, c(0,75,120)), per = cut(year, c(1990,2010,2014), dig.lab = 8) ))
testthat::expect_message( ci2 <- sir( coh.data = c, coh.obs = c('from0to1'), coh.pyrs = 'pyrs', conf.type = 'wald',
ref.data = data.table(popEpi::popmort), ref.rate = 'haz',
adjust = c('agegroup','year','sex'),
print = list(FOT = fot, sex, age.cat = cut(agegroup, c(0,75,120)), per = cut(year, c(1990,2010,2014), dig.lab = 8) )))
ci3 <- sir( coh.data = c, coh.obs = c('from0to1'), coh.pyrs = 'pyrs', conf.type = 'univariate',
ref.data = data.table(popEpi::popmort), ref.rate = 'haz',
adjust = c('agegroup','year','sex'),
print = list(FOT = fot, age.cat = cut(agegroup, c(0,75,120)), per = cut(year, c(1990,2010,2014), dig.lab = 8) ))
testthat::expect_equal(attributes(ci1)$sir.meta$conf.type,'profile')
testthat::expect_equal(attributes(ci2)$sir.meta$conf.type,'wald')
testthat::expect_equal(attributes(ci3)$sir.meta$conf.type,'univariate')
testthat::expect_message( sir( coh.data = c, coh.obs = c('from0to1'), coh.pyrs = 'pyrs', conf.type = 'wald',
ref.data = data.table(popEpi::popmort), ref.rate = 'haz',
adjust = c('agegroup','year','sex'),
print = list(sex )))
testthat::expect_message( ci4 <- sir( coh.data = c, coh.obs = c('from0to2'), coh.pyrs = 'pyrs', conf.type = 'profile',
ref.data = data.table(popEpi::popmort), ref.rate = 'haz',
adjust = c('agegroup','year','sex'),
print = list(dg = cut(agegroup, c(0,5,10,20,40,75,120)) )) )
testthat::expect_equal(attributes(ci4)$sir.meta$conf.type,'wald')
testthat::expect_message( ci5 <- sir( coh.data = c, coh.obs = c('from0to2'), coh.pyrs = 'pyrs', conf.type = 'profile',
ref.data = data.table(popEpi::popmort), ref.rate = 'haz',
adjust = c('agegroup','year','sex'),
print = list(sex=sex,FOT = fot,
age.cat = cut(agegroup, c(0,50,75,120)),
per = cut(year, c(1990,2010,2014), dig.lab = 8) )) )
# ear
ea <- sir( coh.data = c, coh.obs = c('from0to2'), coh.pyrs = 'pyrs', conf.type = 'wald', EAR = TRUE,
ref.data = data.table(popEpi::popmort), ref.rate = 'haz',
adjust = c('agegroup','year','sex'), conf.level = 0.95)
testthat::expect_equal(ea[,EAR], round((1490-1482.13)/39905.92*1000,3))
# Update
ud1 <- sir( coh.data = c, coh.obs = c('from0to1'), coh.pyrs = 'pyrs', conf.type = 'profile',
ref.data = data.table(popEpi::popmort), ref.rate = 'haz',
adjust = c('agegroup','year','sex'),
print = list(fot))
ud2 <- update(ud1, print = list(fot, age.cat = cut(agegroup, c(0,75,120))))
ud3 <- sir( coh.data = c, coh.obs = c('from0to1'), coh.pyrs = 'pyrs', conf.type = 'profile',
ref.data = data.table(popEpi::popmort), ref.rate = 'haz',
adjust = c('agegroup','year','sex'),
print = list(fot, age.cat = cut(agegroup, c(0,75,120))))
testthat::expect_equal(ud2, ud3)
})
testthat::test_that("SIR works with multistate aggregated lexpand data", {
## don't skip on CRAN
testthat::expect_warning(
se <- sir( coh.data = c, coh.obs = c('from0to1','from0to2'), coh.pyrs = 'pyrs',
subset = year %in% 1990:2009,
ref.data = data.table(popEpi::popmort), ref.rate = 'haz',
adjust = c('agegroup','year','sex'), print =c('cause','fot'), mstate = 'cause')
)
dummy.coh <- data.table(agegroup = 1:18, sex = 1, from0to1 = round(runif(18)),
from0to2 = round(runif(18)+0.4),
pyrs = (rnorm(18)+100)*200)
dummy.ref <- data.table(agegroup = rep(1:18,times=2), sex = 1, obs = floor(runif(36)*100),
categ2 = rep(c(1,2),each=18),pyrs = (rnorm(36)+100)*2000)
sm <- sir( coh.data = dummy.coh, coh.obs = c('from0to1','from0to2'), coh.pyrs = 'pyrs',
ref.data = dummy.ref, ref.obs = 'obs', ref.pyrs = 'pyrs', mstate = 'categ2',
adjust = c('agegroup','sex','categ2'), print =c('categ2'))
setkey(12345)
dummy.coh2 <- data.table(agegroup = rep(1:18,times=6), sex = rep(1:2, each= 18), obs = floor(runif(18*6)*100), level = rep(1:6,each=18),pyrs = (rnorm(36)+100)*2000)
sg <- sir( coh.data = dummy.coh2, coh.obs = c('obs'), coh.pyrs = 'pyrs',
adjust = c('agegroup'), print =c('level'))
testthat::expect_equal(attributes(sg)$sir.meta$pooled[,sir], 1)
testthat::expect_equal(sg[,.N], 6)
s1 <- sir( coh.data = c, coh.obs = c('from0to1'), coh.pyrs = 'pyrs', subset = year %in% 1990:2009,
ref.data = data.table(popEpi::popmort), ref.rate = 'haz',
adjust = c('agegroup','year','sex'), print =c('fot'))
s2 <- sir( coh.data = c2, coh.obs = c('from0to2'), coh.pyrs = 'pyrs',
ref.data = data.table(popEpi::popmort), ref.rate = 'haz',
adjust = c('agegroup','year','sex'), print =c('fot'))
s12 <- rbind( cbind(cause=1L, data.table(s1)),
cbind(cause=2L, data.table(s2)))
# compare by hand calculated sir
r <- merge(c2, data.table(popEpi::popmort), c('sex','agegroup','year'), all.x=TRUE)
setDT(r)
r[, exp := haz*pyrs]
est <- r[, list(observed=sum(from0to1, na.rm=TRUE),expected=sum(exp, na.rm=TRUE)), by=.(fot)]
est <- round(est,4)
testthat::expect_is(object = se, class = 'sir')
testthat::expect_equivalent(se, s12)
setDT(s1)
testthat::expect_equivalent(s1[,1:3, with=FALSE], est)
})
# SIR utils ---------------------------------------------------------------
testthat::test_that('Util functions work', {
testthat::expect_equal( poisson.ci(5,5)$rate, 1)
testthat::expect_equal(p.round(0.0000001), "< 0.001")
})
testthat::test_that("data.table subsettin works", {
ss <- sir( coh.data = c2, coh.obs = c('from0to2'), coh.pyrs = 'pyrs',
ref.data = data.table(popEpi::popmort), ref.rate = 'haz',
adjust = c('agegroup','year','sex'), print =c('fot'))
testthat::expect_equal(ss[1,.N], 1)
testthat::expect_equal(c(ss[,.(fot)]), list(fot = c(0,10)))
})
# sir splines -------------------------------------------------------------
if (requireNamespace("splines")) {
testthat::test_that("SIR spline throws errors correctly", {
popEpi:::skip_normally()
sp0 <- suppressWarnings(try(sirspline( coh.data = c, coh.obs = 'from0to2', coh.pyrs = 'pyrs',
subset = year %in% 1990:2008,
ref.data = data.table(popEpi::popmort), ref.rate = 'haz',
adjust = c('agegroup','year','sex'), print = NULL,
spline=c('agegroup','year','fot') )))
sp1 <- suppressWarnings(try(sirspline( coh.data = c, coh.obs = c('from0to1','from0to2'), coh.pyrs = 'pyrs',
subset = year %in% 1990:2008,
ref.data = data.table(popEpi::popmort), ref.rate = 'haz',
adjust = c('agegroup','year','sex'), print =c('cause'),
mstate = 'cause', spline=c('agegroup','year','fot') )))
sp2 <- suppressWarnings(try(sirspline( coh.data = c, coh.obs = c('from0to1','from0to2'), coh.pyrs = 'pyrs',
subset = year %in% 1990:2008,
ref.data = data.table(popEpi::popmort), ref.rate = 'haz', dependent.spline=FALSE,
adjust = c('agegroup','year','sex'), print =c('cause'),
mstate = 'cause', spline=c('agegroup','year','fot') )))
sp3 <- suppressWarnings(try(sirspline( coh.data = c, coh.obs = c('from0to1','from0to2'), coh.pyrs = 'pyrs',
subset = year %in% 1990:2008,
ref.data = data.table(popEpi::popmort), ref.rate = 'haz',
adjust = c('agegroup','year','sex'), print =c('cause'),
mstate = 'cause', spline='agegroup') ))
sp4 <- suppressWarnings(try(sirspline( coh.data = c, coh.obs = c('from0to1','from0to2'), coh.pyrs = 'pyrs',
subset = year %in% 1990:2008,
ref.data = data.table(popEpi::popmort), ref.rate = 'haz', reference.points = c(2000,4),
adjust = c('agegroup','year','sex'), print =c('cause'),
mstate = 'cause', spline=c('agegroup','year','fot') )))
testthat::expect_is( object = sp1, class = 'sirspline')
testthat::expect_is( object = sp2, class = 'sirspline')
testthat::expect_is( object = sp3, class = 'sirspline')
testthat::expect_is( object = sp4, class = 'sirspline')
})
}
testthat::test_that("print accepts a function and subset works", {
popEpi:::skip_normally()
testthat::expect_warning(
pl1 <- sir( coh.data = c, coh.obs = c('from0to1','from0to2'), coh.pyrs = 'pyrs',
subset = year %in% 1990:2008,
ref.data = data.table(popEpi::popmort), ref.rate = 'haz',
adjust = c('agegroup','year','sex'), print = list(year.int = findInterval(year,c(1989,2000,2010))),
mstate = 'cause')
)
testthat::expect_warning(
pl2 <- sir( coh.data = c, coh.obs = c('from0to1','from0to2'), coh.pyrs = 'pyrs',
subset = year %in% 1990:2008,
ref.data = data.table(popEpi::popmort), ref.rate = 'haz',
adjust = c('agegroup','year','sex'), print = list(year.int = findInterval(year,c(1989,2000,2010))),
mstate = 'cause')
)
setDT(pl1)
testthat::expect_equal( pl1[,year.int], 1:2)
testthat::expect_is(object=pl2, 'sir')
})
testthat::test_that("sir_ratio", {
## don't skip on CRAN
dt1 <- data.table(obs = rep(c(5,7), 10),
pyrs = rep(c(250,300,350,400), 5),
var = 1:20)
dt2 <- data.table(obs = rep(c(10,13), 10),
pyrs = rep(c(250,300,350,400), 5),
var = 1:20)
Ref <- data.table(obs = rep(c(50,70,80,100), 5),
pyrs = rep(c(2500,3000,3500,4000), 5),
var = 1:20)
s1 <- sir(coh.data = dt1, coh.obs = obs, coh.pyrs = pyrs,
ref.data = Ref, ref.obs = obs, ref.pyrs = pyrs,
adjust = var)
s2 <- sir(coh.data = dt2, coh.obs = obs, coh.pyrs = pyrs,
ref.data = Ref, ref.obs = obs, ref.pyrs = pyrs,
adjust = var)
## SIR w/ coh=ref=popEpi::sire
## don't skip on CRAN
# Test the example from the book (statistics with confidence):
testthat::expect_equal( sir_ratio(x = c(64, 45.6), y = c(25,23.7), digits=2)[2:3], c(lower=0.83, upper=2.21))
# other tests
testthat::expect_message(sir_ratio(s1, s2, digits = 2, alternative = 'less', conf.level = 0.95, type = 'asymptotic'))
testthat::expect_equal(sir_ratio(s1, s2), sir_ratio(c(120,150), c(230,150)), tolerance=0.001)
testthat::expect_equal(sir_ratio(s1, c(230,150)), sir_ratio(c(120,150), s2), tolerance=0.001)
testthat::expect_equal(sir_ratio(s1, s2)[1], c(sir_ratio =round((120/150)/(230/150),3)), tolerance=0.001)
})
testthat::test_that("sir_exp and sir_ag and sir_lexis are working", {
## don't skip on CRAN
# from0to1 ei oo vakio
BL <- list(fot = 0:5, per = c("2003-01-01","2008-01-01", "2013-01-01"))
suppressMessages(
# aggregated with breaks
x1 <- lexpand(sire, breaks = BL, status = status != 0,
birth = bi_date, entry = dg_date, exit = ex_date,
pophaz=popmort,
aggre=list(sex, period = per, surv.int = fot))
)
suppressMessages(
# no aggreate but breaks
x2 <- lexpand(sire, breaks = BL, status = status != 0,
birth = bi_date, entry = dg_date, exit = ex_date,
pophaz=popmort)
)
suppressMessages(
x3 <- lexpand(sire, status = status != 0,
birth = bi_date, entry = dg_date, exit = ex_date,
pophaz=popmort)
)
suppressMessages(
x4 <- lexpand(sire, status = status != 0,
birth = bi_date, entry = bi_date, exit = ex_date,
pophaz=popmort)
)
# library(Epi)
# # no aggreate but breaks
# suppressMessages(
# x5 <- Lexis(entry = list(FUT = 0, AGE = 0, CAL = get.yrs(bi_date)),
# exit = list(CAL = get.yrs(ex_date)),
# data = sire[sire$dg_date < sire$ex_date, ],
# exit.status = factor(status, levels = 0:2,
# labels = c("alive", "canD", "othD")),
# merge = TRUE)
# )
# AGGRE
testthat::expect_message( a1 <- sir_ag(x1, obs = 'from0to1', conf.type = 'profile') )
testthat::expect_message( a0 <- sir_ag(x1, obs = from0to1, conf.type = 'profile') )
a0 <- sir_ag(x1, obs = 'from0to1', print = list(period))
a0 <- sir_ag(x1, obs = 'from0to1', print = list(per=period))
testthat::expect_message( a0 <- sir_ag(x1, obs = from0to1, exp = d.exp, pyrs = pyrs) )
a0 <- sir_ag(x1, obs = from0to1, exp = d.exp, pyrs = pyrs, print = period)
testthat::expect_message(a2 <- sir_ag(x1, print = attr(x, "aggre.meta")$by))
testthat::expect_message(a3 <- sir_ag(x1))
testthat::expect_equal(data.table(a2),data.table(a3))
a4 <- sir_ag(x1, obs = 'from0to1', print = surv.int)
a4b <- sir_ag(x1, obs = 'from0to1', print = c('surv.int','period'))
a5 <- sir_ag(x1, obs = 'from0to1', print = list(surv.int,period))
a6 <- sir_ag(x1, obs = 'from0to1', print = c('period','surv.int'))
setcolorder(a5, c(2,1,3:9))
setkey(a5, period)
setkey(a6, period)
testthat::expect_equal(data.table(a5)[order(period, surv.int)], data.table(a6)[order(period, surv.int)])
testthat::expect_error(sir_ag(x2, print = NULL))
# UPDATE
a4c <- update(a4, print = c('surv.int','period')) # update to l2
testthat::expect_equal(a4c,a4b)
# LEXIS
testthat::expect_error(sir_lex(x1))
l1 <- sir_lex(x2, print = attr(x, "aggre.meta")$by)
l2 <- sir_lex(x2, print = 'per')
l2 <- sir_lex(x2, print = per)
l3 <- sir_lex(x2, print = NULL)
testthat::expect_equal(l3[,.N], 1)
l1 <- sir_lex(x2, print = NULL)
l2 <- sir_lex(x2, print = 'per')
l2 <- sir_lex(x2, print = per)
l3 <- sir_lex(x2, print = NULL)
# breaks...
BL1 <- list(fot = 0:5, per = c(2003,2008,2013))
b1 <- sir_lex(x3, breaks = BL1, print = c('fot','per'))
BL <- list(fot = 0:5, per = c('2003-01-01','2008-01-01','2013-01-01'), age = c(0,75,100))
b3 <- sir_lex(x3, breaks = BL, print = c('fot','per','age'))
testthat::expect_equal(b3[,.N], 24)
BL2 <- list(fot = 0:5, age = c(0,50,70,100))
b3 <- sir_lex(x2, breaks = BL2, print = c('fot','age')) # already split data
testthat::expect_equal(b3[,.N], 15)
# Character date
BL <- list(fot = 0:5, per = c("2003-01-01","2008-01-01", "2013-01-01"))
b5 <- sir_lex(x3, breaks = BL, print = c('fot','per'))
BL <- list(fot = 0:5, per = c(2002.999, 2007.999, 2013.001))
b6 <- sir_lex(x3, breaks = BL, print = c('fot','per'))
testthat::expect_equal(data.table(b5), data.table(b6), tolerance = 0.001)
# EXPECTED
# without pyrs;
e0 <- sir_exp(x1, obs = "from0to1", exp = "d.exp", pyrs = , print = period)
e1 <- sir_exp(x1, obs = "from0to1", exp = "d.exp", pyrs = "pyrs", print = period)
e1 <- sir_exp(x1, obs = from0to1, exp = d.exp, pyrs = pyrs, print = period)
e2 <- sir_exp(x2, obs = "lex.Xst", exp = list(exp = pop.haz*lex.dur), pyrs = "lex.dur",
print = list(period = cut(per, c(2003,2008, 2013), right = FALSE)))
e1 <- data.table(e1)
e2 <- data.table(e2)
testthat::expect_equal(e1[,.(pyrs, sir, sir.lo, sir.hi, p_value)], e2[,.(pyrs, sir, sir.lo, sir.hi, p_value)])
# WARNINGS (when missing pop.haz)
testthat::expect_error( f1 <- sir_lex(x4) )
testthat::expect_error( f1 <- sir_exp(x4, obs ='lex.Xst', exp = list(exp = lex.dur*pop.haz), pyrs ='lex.dur') )
# SIR_RATIO
testthat::expect_equal(sir_ratio(a2, a3)[1], c(sir_ratio = 1))
plot(a2)
plot(a4)
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.