Package 'popEpi'

Description

This function is only intended to be used within a formula when supplied to e.g. ⁠[survtab_ag]⁠ and should not be used elsewhere.

Usage

adjust(...)

Arguments

Value

Returns a list of promises of the variables supplied which can be evaluated.

Examples

y ~ x + adjust(z)

Description

Aggregates a split Lexis object by given variables and / or expressions into a long-format table of person-years and transitions / end-points. Automatic aggregation over time scales by which data has been split if the respective time scales are mentioned in the aggregation argument to e.g. intervals of calendar time, follow-up time and/or age.

Usage

aggre(
  lex,
  by = NULL,
  type = c("unique", "full"),
  sum.values = NULL,
  subset = NULL,
  verbose = FALSE
)

Arguments

Details

Basics

aggre is intended for aggregation of split Lexis data only. See ⁠[Epi::Lexis]⁠ for forming Lexis objects by hand and e.g. ⁠[Epi::splitLexis]⁠, ⁠[splitLexisDT]⁠, and ⁠[splitMulti]⁠ for splitting the data. ⁠[lexpand]⁠ may be used for simple data sets to do both steps as well as aggregation in the same function call.

Here aggregation refers to computing person-years and the appropriate events (state transitions and end points in status) for the subjects in the data. Hence, it computes e.g. deaths (end-point and state transition) and censorings (end-point) as well as events in a multi-state setting (state transitions).

The result is a long-format data.frame or data.table (depending on options("popEpi.datatable"); see ?popEpi) with the columns pyrs and the appropriate transitions named as fromXtoY, e.g. from0to0 and from0to1 depending on the values of lex.Cst and lex.Xst.

The by argument

The by argument determines the length of the table, i.e. the combinations of variables to which data is aggregated. by is relatively flexible, as it can be supplied as

• a character string vector, e.g. c("sex", "area"), naming variables existing in lex

• an expression, e.g. factor(sex, 0:1, c("m", "f")) using any variable found in lex

• a list (fully or partially named) of expressions, e.g. ⁠list(gender = factor(sex, 0:1, c("m", "f"), area)⁠

Note that expressions effectively allow a variable to be supplied simply as e.g. by = sex (as a symbol/name in R lingo).

The data is then aggregated to the levels of the given variables or expression(s). Variables defined to be time scales in the supplied Lexis are processed in a special way: If any are mentioned in the by argument, intervals of them are formed based on the breaks used to split the data: e.g. if age was split using the breaks c(0, 50, Inf), mentioning age in by leads to creating the age intervals ⁠[0, 50)⁠ and ⁠[50, Inf)⁠ and aggregating to them. The intervals are identified in the output as the lower bounds of the appropriate intervals.

The order of multiple time scales mentioned in by matters, as the last mentioned time scale is assumed to be a survival time scale for when computing event counts. E.g. when the data is split by the breaks list(FUT = 0:5, CAL = c(2008,2010)), time lines cut short at CAL = 2010 are considered to be censored, but time lines cut short at FUT = 5 are not. See Return.

Aggregation types (styles)

It is almost always enough to aggregate the data to variable levels that are actually represented in the data (default aggre = "unique"; alias "non-empty"). For certain uses it may be useful to have also "empty" levels represented (resulting in some rows in output with zero person-years and events); in these cases supplying aggre = "full" (alias "cartesian") causes aggre to determine the Cartesian product of all the levels of the supplied by variables or expressions and aggregate to them. As an example of a Cartesian product, try

merge(1:2, 1:5).

Value

A long data.frame or data.table of aggregated person-years (pyrs), numbers of subjects at risk (at.risk), and events formatted fromXtoY, where X and X are states transitioning from and to or states at the end of each lex.id's follow-up (implying X = Y). Subjects at risk are computed in the beginning of an interval defined by any Lexis time scales and mentioned in by, but events occur at any point within an interval.

When the data has been split along multiple time scales, the last time scale mentioned in by is considered to be the survival time scale with regard to computing events. Time lines cut short by the extrema of non-survival-time-scales are considered to be censored ("transitions" from the current state to the current state).

Author(s)

Joonas Miettinen

See Also

⁠[aggregate]⁠ for a similar base R solution, and ⁠[ltable]⁠ for a data.table based aggregator. Neither are directly applicable to split Lexis data.

Other aggregation functions: as.aggre(), lexpand(), setaggre(), summary.aggre()

Examples

## form a Lexis object
library(Epi)
data(sibr)
x <- sibr[1:10,]
x[1:5,]$sex <- 0 ## pretend some are male
x <- Lexis(data = x,
           entry = list(AGE = dg_age, CAL = get.yrs(dg_date)),
           exit = list(CAL = get.yrs(ex_date)),
           entry.status=0, exit.status = status)
x <- splitMulti(x, breaks = list(CAL = seq(1993, 2013, 5),
                                 AGE = seq(0, 100, 50)))

## these produce the same results (with differing ways of determining aggre)
a1 <- aggre(x, by = list(gender = factor(sex, 0:1, c("m", "f")),
             agegroup = AGE, period = CAL))

a2 <- aggre(x, by = c("sex", "AGE", "CAL"))

a3 <- aggre(x, by = list(sex, agegroup = AGE, CAL))

## returning also empty levels
a4 <- aggre(x, by = c("sex", "AGE", "CAL"), type = "full")

## computing also expected numbers of cases
x <- lexpand(sibr[1:10,], birth = bi_date, entry = dg_date,
             exit = ex_date, status = status %in% 1:2,
             pophaz = popmort, fot = 0:5, age = c(0, 50, 100))
x$d.exp <- with(x, lex.dur*pop.haz)
## these produce the same result
a5 <- aggre(x, by = c("sex", "age", "fot"), sum.values = list(d.exp))
a5 <- aggre(x, by = c("sex", "age", "fot"), sum.values = "d.exp")
a5 <- aggre(x, by = c("sex", "age", "fot"), sum.values = d.exp)
## same result here with custom name
a5 <- aggre(x, by = c("sex", "age", "fot"),
             sum.values = list(expCases = d.exp))

## computing pohar-perme weighted figures
x$d.exp.pp <- with(x, lex.dur*pop.haz*pp)
a6 <- aggre(x, by = c("sex", "age", "fot"),
             sum.values = c("d.exp", "d.exp.pp"))
## or equivalently e.g. sum.values = list(expCases = d.exp, expCases.p = d.exp.pp).

Description

Given a character vector, checks if all names are present in names(data). Throws error if stops=TRUE, else returns FALSE if some variable name is not present.

Usage

all_names_present(data, var.names, stops = TRUE, msg = NULL)

Arguments

Value

TRUE if all var.names are in data, else FALSE,

Author(s)

Joonas Miettinen

See Also

⁠[robust_values]⁠

Description

Utilities to transform objects between array, data.frame, and survival::ratetable.

Usage

long_df_to_array(x, stratum.col.nms, value.col.nm)

long_df_to_ratetable(
  x,
  stratum.col.nms,
  value.col.nm,
  dim.types,
  cut.points = NULL
)

long_dt_to_array(x, stratum.col.nms, value.col.nm)

long_dt_to_ratetable(
  x,
  stratum.col.nms,
  value.col.nm,
  dim.types,
  cut.points = NULL
)

array_to_long_df(x)

array_to_long_dt(x)

array_to_ratetable(x, dim.types, cut.points = NULL)

ratetable_to_array(x)

ratetable_to_long_df(x)

ratetable_to_long_dt(x)

Arguments

Details

• long_df_to_array: converts a long-format data.frame to an array with one or more dimensions

• long_df_to_ratetable: calls long_df_to_array and then array_to_ratetable

• long_dt_to_array: simply asserts that x is a data.table and calls long_df_to_array

• long_dt_to_ratetable: calls long_dt_to_array and then array_to_ratetable

• array_to_long_df: converts an array with one or more dimensions into a long-format data.frame; any dimnames are used to name and fill the stratifying columns; for dimensions without a name, ".dX" is used for stratifying column number X; for each k, if there are no contents in dimnames(x)[[k]], the elements of seq(dim(x)[k]) are used to fill the corresponding stratifying column; the value column always has the name "value"

• array_to_long_dt: calls array_to_long_df and converts result to a data.table for convenience

• array_to_ratetable: converts an array to a survival::ratetable

• ratetable_to_array: converts a survival::ratetable to an array

• ratetable_to_long_df: calls ratetable_to_array and then array_to_long_df

• ratetable_to_long_dt: calls ratetable_to_array and then array_to_long_dt

Value

• long_df_to_array: an array

• long_df_to_ratetable: a survival::ratetable

• long_dt_to_array: an array

• long_dt_to_ratetable: a survival::ratetable

• array_to_long_df: an data.frame

• array_to_long_dt: an data.table

• array_to_ratetable: a survival::ratetable

• ratetable_to_array: an array

• ratetable_to_long_df: a data.frame

• ratetable_to_long_dt: a data.table

Examples

long_dt <- popEpi::popmort
arr <- long_df_to_array(long_dt, c("agegroup", "year", "sex"), "haz")
rt <- array_to_ratetable(arr, dim.types = c(2L, 4L, 1L))

arr2 <- ratetable_to_array(rt)
long_df2 <- array_to_long_df(arr2)

identical(sort(long_dt[["haz"]]), sort(long_df2[["value"]]))

Description

Coerces an R object to an aggre object, identifying the object as one containing aggregated counts, person-years and other information.

Usage

as.aggre(x, values = NULL, by = NULL, breaks = NULL, ...)

## S3 method for class 'data.frame'
as.aggre(x, values = NULL, by = NULL, breaks = NULL, ...)

## S3 method for class 'data.table'
as.aggre(x, values = NULL, by = NULL, breaks = NULL, ...)

## Default S3 method:
as.aggre(x, ...)

Arguments

Value

Returns a copy of x with attributes set to those of an object of class "aggre".

Methods (by class)

• as.aggre(data.frame): Coerces a data.frame to an aggre object

• as.aggre(data.table): Coerces a data.table to an aggre object

• as.aggre(default): Default method for as.aggre (stops computations if no class-specific method found)

Author(s)

Joonas Miettinen

See Also

Other aggregation functions: aggre(), lexpand(), setaggre(), summary.aggre()

Examples

library("data.table")
df <- data.frame(sex = rep(c("male", "female"), each = 5),
                 obs = rpois(10, rep(7,5, each=5)),
                 pyrs = rpois(10, lambda = 10000))
dt <- as.data.table(df)

df <- as.aggre(df, values = c("pyrs", "obs"), by = "sex")
dt <- as.aggre(dt, values = c("pyrs", "obs"), by = "sex")

class(df)
class(dt)

BL <- list(fot = 0:5)
df <- data.frame(df)
df <- as.aggre(df, values = c("pyrs", "obs"), by = "sex", breaks = BL)

Description

Coerces an yrs object to a Date object. Some loss of information comes if year.length = "approx" was set when using ⁠[get.yrs]⁠, so the transformation back to Date will not be perfect there. With year.length = "actual" the original values are perfectly retrieved.

Usage

## S3 method for class 'yrs'
as.Date(x, ...)

Arguments

Value

A vector of Date values based on the input fractional years.

Author(s)

Joonas Miettinen

See Also

⁠[get.yrs]⁠

Examples

data("sire", package = "popEpi")

## approximate year lengths: here 20 % have an extra day added
sire$dg_yrs <- get.yrs(sire$dg_date)
summary(sire$dg_yrs)
dg_date2 <- as.Date(sire$dg_yrs)
summary(as.numeric(dg_date2 - as.Date(sire$dg_date)))

## using actual year lengths
sire$dg_yrs <- get.yrs(sire$dg_date, year.length = "actual")
summary(sire$dg_yrs)
dg_date2 <- as.Date(sire$dg_yrs)
summary(as.numeric(dg_date2 - as.Date(sire$dg_date)))

Description

Convenience function for using ⁠[data.table::dcast.data.table]⁠; inputs are character strings (names of variables) instead of a formula.

Usage

cast_simple(data = NULL, columns = NULL, rows = NULL, values = NULL)

Arguments

Details

This function is just a small interface for dcast / dcast.data.table and less flexible than the originals.

Note that all data.table objects are also data.frame objects, but that each have their own dcast method. ⁠[data.table::dcast.data.table]⁠ is faster.

If any values in value.vars need to be aggregated, they are aggregated using sum. See ?dcast.

Value

A data.table just like ⁠[data.table::dcast]⁠.

Author(s)

Matti Rantanen, Joonas Miettinen

Examples

library("data.table")
## e.g. silly counts from a long-format table to a wide format
test <- data.table::copy(popEpi::sire)
test$dg_y <- year(test$dg_date)
test$ex_y <- year(test$ex_date)
tab <- ltable(test, c("dg_y","ex_y"))
cast_simple(tab, columns='dg_y', rows="ex_y", values="obs")

Description

Selects lowest values of each factor after cut() based on the assumption that the value starts from index 2 and end in comma ",".

Usage

cut_bound(t, factor = TRUE)

Arguments

Details

type = 'factor': "[50,52)" -> "50-51" OR "[50,51)" -> "50"

type = 'numeric': lowest bound in numeric.

Value

If factor = TRUE, returns a character vector; else returns a numeric vector.

Author(s)

Matti Rantanen

Examples

cut_bound("[1900, 1910)") ## "1900-1909"

Description

Several functions in popEpi have support for direct standardization of estimates. This document explains the usage of weighting with those functions.

Details

Direct standardization is performed by computing estimates of E by the set of adjusting variables A, to which a set of weights W is applicable. The weighted average over A is then the direct-adjusted estimate of E (⁠E*⁠).

To enable both quick and easy as well as more rigorous usage of direct standardization with weights, the weights arguments in popEpi can be supplied in several ways. Ability to use the different ways depends on the number of adjusting variables.

The weights are always handled internally to sum to 1, so they do not need to be scaled in this manner when they are supplied. E.g. counts of subjects in strata may be passed.

Basic usage - one adjusting variable

In the simple case where we are adjusting by only one variable (e.g. by age group), one can simply supply a vector of weights:

FUN(weights = c(0.1, 0.25, 0.25, 0.2, 0.2))

which may be stored in advance:

w <- c(0.1, 0.25, 0.25, 0.2, 0.2)

FUN(weights = w)

The order of the weights matters. popEpi functions with direct adjusting enabled match the supplied weights to the adjusting variables as follows: If the adjusting variable is a factor, the order of the levels is used. Otherwise, the alphabetic order of the unique values is used (try sort to see how it works). For clarity and certainty we recommend using factor or numeric variables when possible. character variables should be avoided: to see why, try sort(15:9) and sort(as.character(15:9)).

It is also possible to supply a character string corresponding to one of the age group standardization schemes integrated into popEpi:

• 'europe_1976_18of5' - European std. population (1976), 18 age groups

• 'nordic_2000_18of5' - Nordic std. population (2000), 18 age groups

• 'world_1966_18of5' - world standard (1966), 18 age groups

• 'world_2000_18of5' - world standard (2000), 18 age groups

• 'world_2000_20of5' - world standard (2000), 20 age groups

• 'world_2000_101of1' - world standard (2000), 101 age groups

Additionally, ⁠[ICSS]⁠ contains international weights used in cancer survival analysis, but they are not currently usable by passing a string to weights and must be supplied by hand.

You may also supply weights = "internal" to use internally computed weights, i.e. usually simply the counts of subjects / person-time experienced in each stratum. E.g.

FUN(weights = "world_2000_18of5")

will use the world standard population from 2000 as weights for 18 age groups, that your adjusting variable is assumed to contain. The adjusting variable must be coded in this case as a numeric variable containing 1:18 or as a factor with 18 levels (coded from the youngest to the oldest age group).

More than one adjusting variable

In the case that you employ more than one adjusting variable, separate weights should be passed to match to the levels of the different adjusting variables. When supplied correctly, "grand" weights are formed based on the variable-specific weights by multiplying over the variable-specific weights (e.g. if men have w = 0.5 and the age group 0-4 has w = 0.1, the "grand" weight for men aged 0-4 is 0.5*0.1). The "grand" weights are then used for adjusting after ensuring they sum to one.

When using multiple adjusting variables, you are allowed to pass either a named list of weights or a data.frame of weights. E.g.

WL <- list(agegroup = age_w, sex = sex_w)

FUN(weights = WL)

where age_w and sex_w are numeric vectors. Given the conditions explained in the previous section are satisfied, you may also do e.g.

WL <- list(agegroup = "world_2000_18of", sex = sex_w)

FUN(weights = WL)

and the world standard pop is used as weights for the age groups as outlined in the previous section.

Sometimes using a data.frame can be clearer (and it is fool-proof as well). To do this, form a data.frame that repeats the levels of your adjusting variables by each level of every other adjusting variable, and assign the weights as a column named "weights". E.g.

wdf <- data.frame(sex = rep(0:1, each = 18), agegroup = rep(1:18, 2))

wdf$weights <- rbinom(36, size = 100, prob = 0.25)

FUN(weights = wdf)

If you want to use the counts of subjects in strata as the weights, one way to do this is by e.g.

wdf <- as.data.frame(x$V1, x$V2, x$V3) names(wdf) <- c("V1", "V2", "V3", "weights")

Author(s)

Joonas Miettinen

References

Source of the Nordic standard population in 5-year age groups (also contains European & 1966 world standards): https://www-dep.iarc.fr/NORDCAN/english/glossary.htm

Source of the 1976 European standard population:

Waterhouse, J.,Muir, C.S.,Correa, P.,Powell, J., eds (1976). Cancer Incidence in Five Continents, Vol. III. IARC Scientific Publications, No. 15, Lyon, IARC. ISBN: 9789283211150

Source of 2000 world standard population in 1-year age groups: https://seer.cancer.gov/stdpopulations/stdpop.singleages.html

See Also

Other weights: ICSS, stdpop101, stdpop18

Other popEpi argument evaluation docs: flexible_argument

Description

Convert factor variable with numbers as levels into a numeric variable

Usage

fac2num(x)

Arguments

Details

For example, a factor with levels c("5","7") is converted into a numeric variable with values c(5,7).

Value

A numeric vector based on the levels of x.

Source

Stackoverflow thread

See Also

⁠[robust_values]⁠

Examples

## this is often not intended
as.numeric(factor(c(5,7))) ## result: c(1,2)
## but this
fac2num(factor(c(5,7))) ## result: c(5,7)

## however
as.numeric(factor(c("5","7","a"))) ## 1:3

suppressWarnings(
  fac2num(factor(c("5","7","a"))) ## c(5,7,NA)
)

Description

Certain arguments in popEpi can be passed in multiple ways. This document shows the usage and a pitfall in the usage of such flexible arguments.

Details

Flexible arguments in popEpi are used to pass variables existing in your data or in the environment where the function is used (for everyday users this is the global environment - in simple terms, where your data is / your work space). The flexible arguments are modelled after the by argument in data.tables - see ?data.table. There are many ways to supply the same information to certain functions in popEpi, but the possible ways listed below may be limited in some of them to only allow for using only a part of them.

Everyday usage

Most commonly you may pass variable names as character strings, e.g.

FUN(arg = c("V1", "V2"), data = x)

which may be stored in advance:

vars <- c("V1", "V2")

FUN(arg = vars, data = x)

where x contains those variables. You may also supply variable names as symbols:

FUN(arg = V1, data = x)

Or as a list of symbols (similarly to as in ⁠[aggregate]⁠):

FUN(arg = list(V1, V2), data = x)

Or as a list of expressions:

FUN(arg = list(V1 + 1, factor(V2)), data = x)

A formula without a left-hand-side specified is sometimes allowed as well:

FUN(arg = ~ I(V1 + 1) + factor(V2), data = x)

Using a symbol or a list of symbols/expressions typically causes the function to look for the variable(s) first in the supplied data (if any) and then where the function was called. For everyday users this means you might define e.g.

V3 <- factor(letters)

and do e.g.

FUN(arg = list(V1 + 1, factor(V2), V3), data = x)

provided V1 and V2 exist in x or in the function calling environment.

A pitfall

There is one way to use flexible arguments incorrectly: By supplying the name of a variable which exists both in the supplied data and the calling environment, and intending the latter to be used. E.g.

vars <- c("V2")

FUN(arg = V3, data = x)

where x has a column named vars. This causes the function to use x$vars and NOT x$V2.

Advanced

Function programmers are advised to pass character strings whenever possible. To fool-proof against conflicts as described in the section above, refer to the calling environment explicitly when passing the variable containing the character strings:

TF <- environment() ## current env to refer to

vars <- c("V1", "V2")

FUN(arg = TF$vars, data = x)

Even if x has columns named vars and TF, using TF$vars does not use those columns but only evaluates TF$vars in the calling environment. This is made possible by the fact that data is always passed as a data.frame, within which evaluation of expressions using the dollar operator is not possible. Therefore it is safe to assume the data should not be used. However, lists of expressions will not be checked for dollar use and will fail in conflict situations:

TF <- environment() ## current env to refer to

vars <- letters[1:5]

x <- data.frame(vars = 1:5, TF = 5:1, V1 = 10:6)

FUN(arg = list(TF$vars, V1), data = x)

On the other hand you may typically also pass quoted (⁠[quote]⁠) or substituted ⁠[substitute]⁠ expressions etc., where the env$object trick will work as well:

q <- quote(list(vars, V1))

FUN(arg = TF$q, data = x)

This works even with

a <- 1:5

V1 <- quote(TF$a)

FUN(arg = TF$V1, data = x)

So no conflicts should occur.

Author(s)

Joonas Miettinen

See Also

Other popEpi argument evaluation docs: direct_standardization

Examples

data(sire)
## prepare data for e.g. 5-year "period analysis" for 2008-2012
## note: sire is a simulated cohort integrated into popEpi.
BL <- list(fot=seq(0, 5, by = 1/12))
x <- lexpand(sire, birth = bi_date, entry = dg_date, exit = ex_date,
             status = status %in% 1:2,
             breaks = BL)

x <- aggre(x, by = fot)

## silly example of referring to pyrs data by fixed character string;
## its possible that the real name wont be fixed in a real-life application.
pyrs <- "actual_pyrs"
TF <- environment()
x$actual_pyrs <- as.numeric(x$pyrs)
x$pyrs <- 1

## this works (uses actual_pyrs eventually)
st <- survtab_ag(fot ~ 1, data = x, surv.type = "surv.obs",
                 pyrs = TF$pyrs, d = from0to1,
                 surv.method = "hazard")
## this would be wrong (sees expression 'pyrs' and uses that column,
## which is not what is intended here)
st <- survtab_ag(fot ~ 1, data = x, surv.type = "surv.obs",
                 pyrs = pyrs, d = from0to1,
                 surv.method = "hazard")

Description

Using Date objects, calculates given dates as fractional years.

Usage

get.yrs(x, year.length = "approx", ...)

Arguments

Details

x should preferably be a Date or IDate object, although it can also be a character string variable which is coerced internally to Date format using ⁠[as.Date.character]⁠.

When year.length = 'actual', fractional years are calculated as year + (day_in_year-1)/365 for non-leap-years and as year + (day_in_year-1)/366 for leap years. If year.length = 'approx', fractional years are always calculated as in year + (day_in_year-1)/365.242199.

There is a slight difference, then, between the two methods when calculating durations between fractional years. For meticulous accuracy one might instead want to calculate durations using dates (days) and convert the results to fractional years.

Note that dates are effectively converted to fractional years at 00:00:01 o'clock:

get.yrs("2000-01-01") = 2000, and get.yrs("2000-01-02") = 2000 + 1/365.242199.

Value

A numeric vector of fractional years.

Author(s)

Joonas Miettinen

See Also

⁠[Epi::cal.yr]⁠, ⁠[as.Date.yrs]⁠, ⁠[as.Date]⁠

Examples

data("sire")
sire$dg_yrs <- get.yrs(sire$dg_date)
summary(sire$dg_yrs)

## see: ?as.Date.yrs
dg_date2 <- as.Date(sire$dg_yrs)
summary(as.numeric(dg_date2 - as.Date(sire$dg_date)))

## Epi's cal.yr versus get.yrs
d <- as.Date("2000-01-01")
Epi::cal.yr(d) ## 1999.999
get.yrs(d) ## 2000

## "..." passed on to as.Date, so character / numeric also accepted as input
## (and whatever else as.Date accepts)
get.yrs("2000-06-01")
get.yrs("20000601", format = "%Y%m%d")
get.yrs("1/6/00", format = "%d/%m/%y")

get.yrs(100, origin = "1970-01-01")

Description

Contains three sets age-standardisation weights for age-standardized survival (net, relative or observed).

Format

data.table with columns

• age - lower bound of the age group

• ICSS1 - first set of weights, sums to 100 000

• ICSS2 - second set of weights, sums to 100 000

• ICSS3 - third set of weights, sums to 100 000

Source

ICSS weights (US National Cancer Institute website)

Corazziari, Isabella, Mike Quinn, and Riccardo Capocaccia. "Standard cancer patient population for age standardising survival ratios." European Journal of Cancer 40.15 (2004): 2307-2316.

See Also

Other popEpi data: meanpop_fi, popmort, sibr, sire, stdpop101, stdpop18

Other weights: direct_standardization, stdpop101, stdpop18

Examples

## aggregate weights to a subset of age groups
data(ICSS)
cut <- c(0, 30, 50, 70, Inf)
agegr <- cut(ICSS$age, cut, right = FALSE)
aggregate(ICSS1~agegr, data = ICSS, FUN = sum)

Description

Given a vector or column of year values (numeric or integer), ⁠[is_leap_year]⁠ returns a vector of equal length of logical indicators, i.e. a vector where corresponding leap years have value TRUE, and FALSE otherwise.

Usage

is_leap_year(years)

Arguments

Value

A logical vector where TRUE indicates a leap year.

Author(s)

Joonas Miettinen

Examples

## can be used to assign new columns easily, e.g. a dummy indicator column
df <- data.frame(yrs=c(1900,1904,2005,1995))
df$lyd <- as.integer(is_leap_year(df$yrs))

## mostly it is useful as a condition or to indicate which rows have leap years
which(is_leap_year(df$yrs)) # 2
df[is_leap_year(df$yrs),] # 2nd row

Description

Test whether obj inherits one of Date or IDate.

Usage

is.Date(obj)

Arguments

Value

TRUE if obj is of class "Date" or "IDate".

Author(s)

Joonas Miettinen

See Also

⁠[get.yrs]⁠, ⁠[is_leap_year]⁠, ⁠[as.Date]⁠

Description

Make ⁠[Epi::Lexis]⁠ objects.

Usage

Lexis_fpa(
  data,
  birth = NULL,
  entry = NULL,
  exit = NULL,
  entry.status = NULL,
  exit.status = NULL,
  subset = NULL,
  ...
)

Lexis_dt(...)

Arguments

Value

popEpi::Lexis_fpa

Returns a Lexis object with the additional class data.table. It has the usual columns that Lexis objects have, and with time scale columns fot, per, and age. They are calculated as

fot = entry - entry (to ensure correct format, e.g. difftime)

per = entry

and

age = entry - birth.

popEpi::Lexis_dt

Returns a Lexis object that is also a data.table — therefore output has the classes c("Lexis", "data.table", "data.frame").

Functions

popEpi::Lexis_fpa

popEpi::Lexis_fpa collects data from its inputs to to call ⁠[Epi::Lexis]⁠. This is a convenience function for making a Lexis object with the time scales fot, per, and age.

popEpi::Lexis_dt

popEpi::Lexis_dt performs the following steps:

• Calls ⁠[Epi::Lexis]⁠.

• Calls ⁠[data.table::setDT]⁠ to make output into a data.table, ⁠[data.table::setattr]⁠ to set its class to c("Lexis", "data.table", "data.frame"), and ⁠[data.table::setkeyv]⁠ to sort the output by lex.id and the time scales.

• Returns the Lexis / data.table object.

News for version 0.5.0

Lexis_fpa output is now also a data.table, so the complete class vector is c("Lexis", "data.table", "data.frame").

New function Lexis_dt, a wrapper of Epi::Lexis which also sets class to c("Lexis", "data.table", "data.frame").

See Also

Other Lexis_functions: lexis_merge(), lexis_split_merge_aggregate_by_stratum()

Examples

# popEpi::Lexis_fpa
data("sire", package = "popEpi")
lex <- Lexis_fpa(sire,
                 birth = "bi_date",
                 entry = dg_date,
                 exit = ex_date + 1L,
                 exit.status = "status")

## some special cases
myVar <- "bi_date"
l <- list(myVar = "bi_date")
sire$l <- sire$myVar <- 1

## conflict: myVar taken from data when "bi_date" was intended
lex <- Lexis_fpa(sire,
                 birth = myVar,
                 entry = dg_date,
                 exit = ex_date + 1L,
                 exit.status = "status")

## no conflict with names in data
lex <- Lexis_fpa(sire,
                 birth = l$myVar,
                 entry = dg_date,
                 exit = ex_date + 1L,
                 exit.status = "status")
stopifnot(
  identical(class(lex), c("Lexis", "data.table", "data.frame")),
  Epi::timeScales(lex) %in% c("fot", "per", "age")
)


# popEpi::Lexis_dt
lex_1 <- popEpi::Lexis_dt(
  data = popEpi::sire,
  entry = list(
    ts_fut = 0L,
    ts_age = as.integer(dg_date - bi_date),
    ts_cal = as.integer(dg_date)
  ),
  exit = list(ts_cal = as.integer(ex_date)),
  entry.status = 0L,
  exit.status = status
)
stopifnot(
  class(lex_1) == c("Lexis", "data.table", "data.frame"),
  Epi::timeScales(lex_1) %in% c("ts_fut", "ts_age", "ts_cal")
)

lex_2 <- popEpi::Lexis_dt(
  data = popEpi::sire,
  entry = list(
    ts_fut = 0L,
    ts_age = as.integer(dg_date - bi_date),
    ts_cal = as.integer(dg_date)
  ),
  duration = as.integer(ex_date - dg_date),
  entry.status = 0L,
  exit.status = status
)
stopifnot(
  class(lex_2) == c("Lexis", "data.table", "data.frame"),
  Epi::timeScales(lex_2) %in% c("ts_fut", "ts_age", "ts_cal")
)

Description

Function(s) to merge data into Lexis objects intelligently when merging is (partially) based on time scales.

Usage

lexis_merge(
  lexis,
  merge_dt,
  merge_dt_by,
  merge_dt_harmonisers = NULL,
  subset = NULL,
  optional_steps = NULL,
  lex_dur_multiplier = NULL
)

Arguments

Value

popEpi::lexis_merge

Returns lexis after adding value columns from merge_dt into lexis.

Functions

popEpi::lexis_merge

popEpi::lexis_merge can be used to merge additional information into Lexis data, allowing the use of the Lexis time scales in the merge. The typical use-case is to split Lexis data and then merge population (expected) hazards to the subject-intervals. popEpi::lexis_merge performs the following steps:

• Run optional_steps[["on_entry"]](call_env = call_env, eval_env = eval_env) if that optional_steps element exists. call_env is the environment where this function was called and eval_env is the temporary environment in which the commands that this function consists of are evaluated.

• Run on.exit(optional_steps[["on_exit"]](call_env = call_env, eval_env = eval_env)) if that optional_steps element exists.

• If is.null(merge_dt_harmonisers), attempt to automatically determine the harmonisers making use of cut for each time scale column to merge by as follows:

  • If merge_dt[[col_nm]] is a factor column, look at the levels() of the column. If they all look like ⁠[x, y[⁠, ⁠]x, y]⁠, or x - y, extract breaks from them, infer whether right = FALSE or TRUE, and use the levels() as labels. E.g. factor("[60, 61[") produces list(breaks = 60:61, right = FALSE, labels = "[60, 61[").

  • If merge_dt[[col_nm]] contains numbers we define cut breaks as the unique values of merge_dt[[col_nm]] and as the ceiling max(merge_dt[[col_nm]]) + last_diff. Here last_diff is the difference between the highest and second highest values. E.g. c(1950, 1960, 1970:2020, 2021) for merge_dt[[col_nm]] containing unique values c(1950, 1960, 1970:2020).

  • If merge_dt[[col_nm]] is not of class integer, numeric, or factor, an error is raised because we don't know how to automatically form a harmoniser.

  • Run optional_steps[["pre_default_harmoniser_creation"]](call_env = call_env, eval_env = eval_env, lapply_eval_env = lapply_eval_env) if that optional_steps element exists. Here make_harmoniser_eval_env is similar to eval_env but it is the evaluation environment of the function passed to lapply which attempts to make each default harmoniser.

  • With the cut breaks defined, the automatically created harmoniser becomes a cut call with arguments inferred above and with

    • x = col + lex.dur * lex_dur_multiplier, where col is the current column.

    • dig.lab = 10.

  • Depending on what we are harmonising with, the harmoniser will return either a factor with labels based on the cut call, e.g. ts_cal = 2021.524 turns into ⁠[2021, 2022[⁠ with both 2021 and 2022 in breaks, or the identified interval's lower bound such as 2021.

  • Run optional_steps[["post_default_harmoniser_creation"]](call_env = call_env, eval_env = eval_env, lapply_eval_env = lapply_eval_env) if that optional_steps element exists.

• Run optional_steps[["post_merge_dt_harmonisers"]](call_env = call_env, eval_env = eval_env) if that optional_steps element exists.

• Armed with either user-defined or automatically created merge_dt_harmonisers, they are each evaluated to create a temporary data.table with harmonised data from lexis. This is performed via eval with envir = lexis and enclos = harmoniser_eval_env where harmoniser_eval_env is a temporary environment which contains every argument passed to lexis_merge as well as eval_env and call_env which you can read about in the detailed explanation of optional_steps. lexis is a subset if !is.null(subset). However, if a column has no harmoniser at this point then it is used as-is. For instance there is no need to harmonise stratifying columns such as area because they are not changed by splitting the Lexis data.

• Run optional_steps[["post_harmonisation"]](call_env = call_env, eval_env = eval_env) if that optional_steps element exists.

• Then we perform the actual merge between merge_dt and the harmonised data. The value columns in merge_dt are added into lexis either using data.table::set for a data.table or a simple lexis[, merge_value_col_nms] <- join_result assignment. Note that data.table::set modifies the lexis object in-place.

• Run optional_steps[["post_merge"]](call_env = call_env, eval_env = eval_env) if that optional_steps element exists.

• Each merged-in column from merge_dt (all columns not in merge_dt_by) are inspected for missing values. If there are any, an error is raised. This usually occurs if merge_dt does not contain data for all data in lexis. For instance it only covers years 1950-2020 but lexis contains also data for 2021. This error helps you to spot those problems early isntead of producing nonsense results downstream.

• Returns lexis after adding value columns from merge_dt into lexis.

See Also

Other Lexis_functions: lexis_funs, lexis_split_merge_aggregate_by_stratum()

Examples

# popEpi::lexis_merge
lexis <- Epi::Lexis(
  entry = list(ts_fut = 0.0, ts_cal = 2010.3, ts_age = 56.8),
  exit = list(ts_cal = 2024.9999),
  entry.status = 0L,
  exit.status = 0L
)
lexis$sex <- 0L
lexis <- popEpi::splitMulti(
  data = lexis,
  breaks = list(ts_fut = seq(0, 3, 1 / 12))
)
my_merge_dt <- data.table::CJ(sex = 0:1, ts_age = 0:100, ts_cal = 2000:2025)
data.table::set(
  x = my_merge_dt,
  j = "merge_value",
  value = runif(nrow(my_merge_dt))
)
# lexis is also a data.table and is modified in-place. no need to keep
# output of popEpi::lexis_merge call.
popEpi::lexis_merge(
  lexis = lexis,
  merge_dt = my_merge_dt,
  merge_dt_by = c("sex", "ts_age", "ts_cal")
)
stopifnot(
  "merge_value" %in% names(lexis),
  !is.na(lexis[["merge_value"]])
)
data.table::set(
  x = lexis,
  j = "merge_value",
  value = NULL
)
popEpi::lexis_merge(
  lexis = lexis,
  merge_dt = my_merge_dt,
  merge_dt_by = c("sex", "ts_age", "ts_cal"),
  merge_dt_harmonisers = list(
    ts_cal = quote(as.integer(ts_cal)),
    ts_age = quote(as.integer(ts_age))
  )
)
stopifnot(
  "merge_value" %in% names(lexis),
  !is.na(lexis[["merge_value"]])
)

# special factor columns in merge_dt
data.table::set(
  x = lexis,
  j = "merge_value",
  value = NULL
)
data.table::set(
  x = my_merge_dt,
  j = c("ts_age", "ts_cal"),
  value = list(
    cut(my_merge_dt[["ts_age"]], breaks = 0:101, right = FALSE),
    cut(my_merge_dt[["ts_cal"]], breaks = 2000:2026, right = FALSE,
        dig.lab = 4)
  )
)
popEpi::lexis_merge(
  lexis = lexis,
  merge_dt = my_merge_dt,
  merge_dt_by = c("sex", "ts_age", "ts_cal")
)
stopifnot(
  "merge_value" %in% names(lexis),
  !is.na(lexis[["merge_value"]])
)

Description

Function(s) which split, merge, and aggregate a Lexis object in one go.

Usage

lexis_split_merge_aggregate_by_stratum(
  lexis,
  breaks,
  aggre_exprs,
  aggre_by = NULL,
  merge_dt = NULL,
  merge_dt_by = NULL,
  merge_optional_args = NULL,
  subset = NULL,
  weight_col_nm = NULL,
  split_lexis_column_exprs = NULL,
  breaks_collapse_args = NULL,
  optional_steps = NULL
)

Arguments

Details

popEpi::lexis_split_merge_aggregate_by_stratum can be used to split Lexis (⁠[Epi::Lexis]⁠) data, merge something to it after the merge, and then perform an aggregation step. The following steps are performed:

• Handle aggre_exprs as follows:

  • An aggre_exprs element that is string such as "n_events" is replaced with the appropriate expression retrieved from a table of pre-defined expressions (see below).

  • If a string element is transition-specific, e.g. n_events_[0, 1], we first turn it into its general form (e.g. n_events_[x, y]). Additionally, any uses of x and y in the fetched expression is replaced with the correct states, e.g. sum(lex.Cst == x & lex.Xst == y) is turned into sum(lex.Cst == 0 & lex.Xst == 1).

  • An aggre_exprs element that is an expression already (e.g. my_n_events = quote(sum(lex.Cst == 0 & lex.Xst != 0))) does not receive the same treatment regarding the states and when you write your own expressions you must specify them yourself. Therefore e.g. my_n_events = quote(sum(lex.Cst == x & lex.Xst != y)) will not work.

  • Regardless of the type of the aggre_exprs element, we remove every use of ⁠ * iw⁠ in the expression if individual weights are not used.

  • Table of general aggregation expressions known to popEpi:

• If weight_col_nm was supplied, include in subset only records where lexis[[weight_col_nm]] > 0.

• Call optional_steps[["on_entry"]](eval_env = eval_env, call_env = call_env) if that optional_steps element exists. eval_env is the temporary evaluation environment of popEpi::lexis_split_merge_aggregate_by_stratum which contains all contains all the arguments of popEpi::lexis_split_merge_aggregate_by_stratum and call_env is the environment where it was called.

• Call on.exit(optional_steps[["on_exit"]](eval_env = eval_env, call_env = call_env)) if that optional_steps element exists.

• For each stratum in aggre_by:

  • Run optional_steps[["stratum_on_entry"]](stratum_eval_env = stratum_eval_env, eval_env = eval_env, call_env = call_env) if that optional_steps element exists. stratum_eval_env is the environment where the stratum-specific steps are performed.

  • Loop over every box defined by all other time scales except the last one. E.g. if breaks = list(ts_cal = c(2001, 2004, 2007), ts_fut = 0:5) then we perform what follows below separately for the boxes ⁠]2001, 2004]⁠, ⁠]2004, 2007]⁠.

    • Run optional_steps[["stratum_pre_split"]](stratum_eval_env = stratum_eval_env, eval_env = eval_env, call_env = call_env) if that optional_steps element exists.

    • Using a subset lexis containing data only for the current stratum (e.g. sex = 0), drop data outside the current stratum box (e.g. ts_cal outside of ⁠]2001, 2004]⁠). Also crop follow-up to the end of the box. This does not split the data yet.

    • If !is.null(breaks_collapse_args), run popEpi::lexis_breaks_collapse_1d on the last element of breaks. Arguments lexis and breaks are set automatically. This makes it possible to "collapse" the breaks for each stratum and stratum box combination separately.

    • Run ⁠[splitMulti]⁠ for the second time, this time splitting the data by the last time scale in breaks. Remember that we are inside the box defined by the other time scales already.

    • Run optional_steps[["stratum_post_split"]](stratum_eval_env = stratum_eval_env, eval_env = eval_env, call_env = call_env) if that optional_steps element exists.

    • Run ⁠[lexis_merge]⁠ with merge_dt, merge_dt_by, and merged_optional_args, if merge_dt has been supplied.

    • Run optional_steps[["stratum_post_merge"]](stratum_eval_env = stratum_eval_env, eval_env = eval_env, call_env = call_env) if that optional_steps element exists.

    • Evaluate aggre_exprs as follows:

      • First collect variables mentioned in aggre_exprs using ⁠[all.vars]⁠.

      • There is a pre-defined table of expressions that are evaluated and added as new columns in the split data if they appear in any of the aggre_exprs. See below for the table.

      • Table of expressions which create new columns into split data:

   + It can happen that e.g. a survival interval has absolutely no data in
     it. Especially in sparse data and with delayed entry. For the
     pre-specified aggregation expressions such as `n_events` we ensure
     that empty time scale boxes have value zero. Your custom aggregation
     expressions will result in `NA` values in empty time scale boxes.
* Run
  `optional_steps[["stratum_post_aggregation"]](stratum_eval_env = stratum_eval_env, eval_env = eval_env, call_env = call_env)`
  if that `optional_steps` element exists.

• Replace NA values in value columns starting with t_ or with n_ with zeroes. E.g. t_at_risk will then be zero in intervals with no subjects in follow-up. Even a custom expression such as t_at_risk_squared will be handled in this manner.

• Set proper data.table attributes on the resulting big table and sort it by the stratifying variables.

• Store as metadata a list with names stratum_col_nms, ts_col_nms, and value_col_nms into the attribute named lexis_split_merge_aggregate_by_stratum_meta, where stratum_col_nms = names(aggre_by), ts_col_nms = names(breaks), and value_col_nms are the names of the columns resulting from aggre_exprs.

• Run optional_steps[["post_aggregation"]](eval_env = eval_env, call_env = call_env) if that optional_steps element exists.

• Return a data.table with stratum columns as specified via aggre_by, time scale columns specified via breaks and value columns as specified via aggre_exprs.

Value

Returns a data.table with stratum columns as specified via aggre_by, time scale columns specified via breaks and value columns as specified via aggre_exprs.

See Also

Other Lexis_functions: lexis_funs, lexis_merge()

Examples

# popEpi::split_merge_aggregate_by_stratum
make_pm <- function() {
  pm <- data.table::copy(popEpi::popmort)
  data.table::setnames(
    pm,
    c("year", "agegroup", "haz"),
    c("ts_cal", "ts_age", "h_exp")
  )
  data.table::setkeyv(pm, c("sex", "ts_cal", "ts_age"))
  data.table::setkeyv(pm, c("sex", "ts_cal", "ts_age"))
  return(pm[])
}

make_column_icss_ag <- function(age) {
  cut(
    age,
    breaks = c(0, 60, 70, 80, Inf),
    right = FALSE,
    labels = c("0-59", "60-69", "70-79", "80+")
  )
}

make_standard_weight_dt <- function() {
  return(popEpi::ICSS[
    j = list(
      weight = as.double(sum(.SD[["ICSS1"]]))
    ),
    keyby = list(
      icss_ag = make_column_icss_ag(popEpi::ICSS[["age"]])
    )
  ][])
}

make_sire <- function() {
  sire <- popEpi::sire
  sire <- sire[
    sire[["dg_date"]] < sire[["ex_date"]] &
      data.table::between(
        sire[["ex_date"]],
        as.Date("1999-01-01"),
        as.Date("2003-12-31"),
        incbounds = TRUE
      ) &
      (get.yrs(sire[["ex_date"]]) - get.yrs(sire[["bi_date"]])) < 100
  ]
  sire[j = "my_stratum" := sample(2L, size = nrow(sire), replace = TRUE)]
  sire <- sire[
    j = .SD[as.integer(seq(1L, .N, length.out = 50L))],
    keyby = "my_stratum"
  ]
  # you can also use popEpi::Lexis_dt
  sire <- Epi::Lexis(
    entry = list(
      ts_cal = popEpi::get.yrs(dg_date),
      ts_age = popEpi::get.yrs(dg_date) - popEpi::get.yrs(bi_date),
      ts_fut = 0.0
    ),
    duration = popEpi::get.yrs(ex_date) - popEpi::get.yrs(dg_date),
    entry.status = 0L,
    exit.status = status,
    data = sire
  )
  sire[["icss_ag"]] <- make_column_icss_ag(sire[["dg_age"]])
  sire[["individual_weight"]] <- popEpi::surv_individual_weights(
    df = sire,
    standard_weight_dt = wdt
  )
  sire
}

pm <- make_pm()
wdt <- make_standard_weight_dt()
sire <- make_sire()

# note that we use here 1-year survival intervals which are too wide for
# practical survival estimation. we also aggregate for period analysis
# for demonstration purposes.
bl <- list(ts_cal = c(1999, 2004), ts_fut = 0:5)

# using some pre-defined aggregation expressions
agdt <- popEpi::lexis_split_merge_aggregate_by_stratum(
  lexis = sire,
  breaks = bl,
  aggre_exprs = c("n_events", "t_at_risk", "n_events_[0, 1]"),
  aggre_by = "my_stratum"
)
stopifnot(
  c("n_events", "t_at_risk", "n_events_[0, 1]") %in% names(agdt)
)

# using some pre-defined aggregation expressions + using expected hazard.
agdt <- popEpi::lexis_split_merge_aggregate_by_stratum(
  lexis = sire,
  breaks = bl,
  aggre_exprs = c(
    "n_events", "t_at_risk", "n_events_[0, 1]",
    "n_events_exp_e2"
  ),
  aggre_by = "my_stratum",
  merge_dt = pm,
  merge_dt_by = c("sex", "ts_cal", "ts_age")
)

# using custom expressions --- `h_exp` is merged from `merge_dt`
agdt <- popEpi::lexis_split_merge_aggregate_by_stratum(
  lexis = sire,
  breaks = bl,
  aggre_exprs = list(
    "n_events", "t_at_risk", "n_events_[0, 1]",
    "n_events_exp_e2",
    "my_n_events" = quote(sum(lex.Cst == 0 & lex.Xst != 0)),
    "my_n_events_exp_e2" = quote(sum(h_exp * lex.dur))
  ),
  aggre_by = "my_stratum",
  merge_dt = pm,
  merge_dt_by = c("sex", "ts_cal", "ts_age")
)
stopifnot(
  agdt[["my_n_events"]] == agdt[["n_events"]],
  all.equal(agdt[["my_n_events_exp_e2"]], agdt[["n_events_exp_e2"]])
)

# using custom expressions --- some pre-defined variables which are added
# automatically
agdt <- popEpi::lexis_split_merge_aggregate_by_stratum(
  lexis = sire,
  breaks = bl,
  aggre_exprs = list(
    "n_events_pp",
    "my_n_events_pp" = quote(sum((lex.Cst == 0 & lex.Xst != 0) * pp)),
    "n_at_risk_eff",
    "my_n_at_risk_eff" = quote(sum(
      in_follow_up_at_interval_start +
        0.5 * entered_late_during_interval -
        0.5 * left_early_during_interval
    )),
    "my_n_at_risk_eff_2" = quote({
      # custom definition for detecting late entry. in this example we add
      # a larger tolerance than is used normally. you can find the normal
      # tolerance used in the table of expressions which create columns into
      # the split lexis dataset.
      ts_fut_floor <- eval_env[["box_dt"]][["ts_fut_start"]][box_id]
      distance_from_floor <- (ts_fut - ts_fut_floor)
      entered_late_during_interval <- distance_from_floor > 1e-3
      sum(
        in_follow_up_at_interval_start +
          0.5 * entered_late_during_interval -
          0.5 * left_early_during_interval
      )
    })
  ),
  aggre_by = "my_stratum",
  merge_dt = pm,
  merge_dt_by = c("sex", "ts_cal", "ts_age")
)
stopifnot(
  all.equal(agdt[["my_n_events_pp"]], agdt[["n_events_pp"]]),
  all.equal(agdt[["n_at_risk_eff"]], agdt[["my_n_at_risk_eff"]]),
  all.equal(agdt[["n_at_risk_eff"]], agdt[["my_n_at_risk_eff_2"]])
)

# using custom expressions --- also a custom column at the split lexis level.
# this is for demonstration of what is possible and probably no-one would
# need such an h_exp_mean in a real application.
agdt <- popEpi::lexis_split_merge_aggregate_by_stratum(
  lexis = sire,
  breaks = bl,
  aggre_exprs = list(
    "my_stat" = quote(sum(lex.dur * h_exp_mean))
  ),
  aggre_by = "my_stratum",
  merge_dt = pm,
  merge_dt_by = c("sex", "ts_cal", "ts_age"),
  split_lexis_column_exprs = list(
    h_exp_mean = quote({
      dt <- data.table::setDT(list(lex.id = lex.id, h_exp = h_exp))
      h_exp_lag1 <- dt[
        j = list(lag1 = data.table::shift(h_exp, n = 1L, type = "lag")),
        by = "lex.id"
      ][["lag1"]]
      h_exp_lag1[is.na(h_exp_lag1)] <- h_exp[is.na(h_exp_lag1)]
      h_exp_lead1 <- dt[
        j = list(lead1 = data.table::shift(h_exp, n = 1L, type = "lead")),
        by = "lex.id"
      ][["lead1"]]
      h_exp_lead1[is.na(h_exp_lead1)] <- h_exp[is.na(h_exp_lead1)]
      (h_exp + h_exp_lag1 + h_exp_lead1) / 3
    })
  )
)
stopifnot(
  "my_stat" %in% names(agdt),
  !is.na(agdt[["my_stat"]])
)

# this function was written for survival estimation but is is in fact more
# general. for instance we can look at a calendar year / age grid.
agdt <- popEpi::lexis_split_merge_aggregate_by_stratum(
  lexis = sire,
  breaks = list(ts_cal = 1999:2004, ts_age = seq(0, 100, 10)),
  aggre_exprs = list(
    "n_events_exp_e2"
  ),
  merge_dt = pm,
  merge_dt_by = c("sex", "ts_cal", "ts_age")
)
stopifnot(
  c("ts_cal_start", "ts_age_start") %in% names(agdt)
)
agdt_wide <- data.table::dcast(
  data = agdt,
  formula = ts_age_start ~ ts_cal_start,
  value.var = "n_events_exp_e2"
)
# print(round(agdt_wide, 2))
# Key: <ts_age_start>
#     ts_age_start  1999  2000  2001  2002  2003
#            <num> <num> <num> <num> <num> <num>
#  1:            0  0.00  0.00  0.00  0.00  0.00
#  2:           10  0.00  0.00  0.00  0.00  0.00
#  3:           20  0.00  0.00  0.00  0.00  0.00
#  4:           30  0.00  0.00  0.00  0.00  0.00
#  5:           40  0.00  0.00  0.00  0.00  0.00
#  6:           50  0.01  0.01  0.01  0.00  0.00
#  7:           60  0.04  0.04  0.03  0.01  0.02
#  8:           70  0.33  0.28  0.31  0.16  0.05
#  9:           80  1.36  1.04  0.89  0.48  0.20
# 10:           90  0.39  0.29  0.47  0.37  0.56

# this is what happens when there is an empty interval
agdt <- popEpi::lexis_split_merge_aggregate_by_stratum(
  lexis = popEpi::Lexis_dt(
    entry = list(ts_fut = c(0.0, 3.5)),
    exit = list(ts_fut = c(1.5, 4.5)),
    entry.status = 0L,
    exit.status = 1L
  ),
  breaks = list(ts_fut = 0:5),
  aggre_exprs = list(
    "t_at_risk",
    "n_events",
    my_n_events = quote(sum(lex.Cst == 0 & lex.Xst != 0))
  )
)
stopifnot(
  identical(agdt[["t_at_risk"]], c(1.0, 0.5, 0.0, 0.5, 0.5)),
  identical(agdt[["n_events"]], c(0L, 1L, 0L, 0L, 1L)),
  identical(agdt[["my_n_events"]], c(0L, 1L, NA, 0L, 1L))
)

# this is what happens when an empty stratum is requested
agdt <- popEpi::lexis_split_merge_aggregate_by_stratum(
  lexis = sire,
  aggre_by = data.table::data.table(sex = 0:2),
  breaks = list(ts_fut = 0:5),
  aggre_exprs = list(
    "t_at_risk",
    "my_stat" = quote(sum(lex.dur ^ 2)),
    "t_at_risk_squaredd" = quote(sum(lex.dur ^ 2))
  )
)
stopifnot(
  nrow(agdt) == length(0:2) * (length(0:5) - 1),
  agdt[["t_at_risk"]][!agdt[["sex"]] %in% sire[["sex"]]] == 0,
  is.na(agdt[["my_stat"]][!agdt[["sex"]] %in% sire[["sex"]]]),
  agdt[["t_at_risk_squared"]][!agdt[["sex"]] %in% sire[["sex"]]] == 0
)

# and this is what happens when no records are in the dataset
agdt <- popEpi::lexis_split_merge_aggregate_by_stratum(
  lexis = sire,
  subset = rep(FALSE, nrow(sire)),
  aggre_by = data.table::data.table(sex = 0:2),
  breaks = list(ts_fut = 0:5),
  aggre_exprs = list(
    "t_at_risk",
    "my_stat" = quote(sum(lex.dur ^ 2))
  )
)
stopifnot(
  nrow(agdt) == length(0:2) * (length(0:5) - 1),
  agdt[["t_at_risk"]] == 0,
  is.na(agdt[["my_stat"]])
)

# and this is what happens when no records are in the dataset
agdt <- popEpi::lexis_split_merge_aggregate_by_stratum(
  lexis = sire,
  subset = rep(FALSE, nrow(sire)),
  aggre_by = data.table::data.table(sex = 0:2),
  breaks = list(ts_fut = 0:5),
  aggre_exprs = list(
    "t_at_risk",
    "my_stat" = quote(sum(lex.dur ^ 2))
  )
)
stopifnot(
  nrow(agdt) == length(0:2) * (length(0:5) - 1),
  agdt[["t_at_risk"]] == 0,
  is.na(agdt[["my_stat"]])
)

Description

Given subject-level data, data is split by calendar time (per), age, and follow-up time (fot, from 0 to the end of follow-up) into subject-time-interval rows according to given breaks and additionally processed if requested.

Usage

lexpand(
  data,
  birth = NULL,
  entry = NULL,
  exit = NULL,
  event = NULL,
  status = status != 0,
  entry.status = NULL,
  breaks = list(fot = c(0, Inf)),
  id = NULL,
  overlapping = TRUE,
  aggre = NULL,
  aggre.type = c("unique", "cartesian"),
  drop = TRUE,
  pophaz = NULL,
  pp = TRUE,
  subset = NULL,
  merge = TRUE,
  verbose = FALSE,
  ...
)

Arguments

Details

Basics

⁠[lexpand]⁠ splits a given data set (with e.g. cancer diagnoses as rows) to subintervals of time over calendar time, age, and follow-up time with given time breaks using ⁠[splitMulti]⁠.

The dataset must contain appropriate Date / IDate format or other numeric variables that can be used as the time variables.

You may take a look at a simulated cohort ⁠[sire]⁠ as an example of the minimum required information for processing data with lexpand.

Many arguments can be supplied as a character string naming the appropriate variable (e.g. "sex"), as a symbol (e.g. sex) or as an expression (e.g. factor(sex, 0:1, c("m", "f"))) for flexibility.

Breaks

You should define all breaks as left inclusive and right exclusive time points (e.g.⁠[a,b)⁠ ) for 1-3 time dimensions so that the last member of a breaks vector is a meaningful "final upper limit", e.g. per = c(2002,2007,2012) to create a last subinterval of the form ⁠[2007,2012)⁠.

All breaks are explicit, i.e. if drop = TRUE, any data beyond the outermost breaks points are dropped. If one wants to have unspecified upper / lower limits on one time scale, use Inf: e.g. breaks = list(fot = 0:5, age = c(0,45,Inf)). Breaks for per can also be given in Date / IDateformat, whereupon they are converted to fractional years before used in splitting.

The age time scale can additionally be automatically split into common age grouping schemes by naming the scheme with an appropriate character string:

• "18of5": age groups 0-4, 5-9, 10-14, ..., 75-79, 80-84, 85+

• "20of5": age groups 0-4, 5-9, 10-14, ..., 85-89, 90-94, 95+

• "101of1": age groups 0, 1, 2, ..., 98, 99, 100+

Time variables

If any of the given time variables (birth, entry, exit, event) is in any kind of date format, they are first coerced to fractional years before splitting using ⁠[get.yrs]⁠ (with year.length = "actual").

Sometimes in e.g. SIR/SMR calculation one may want the event time to differ from the time of exit from follow-up, if the subject is still considered to be at risk of the event. If event is specified, the transition to status is moved to event from exit using ⁠[Epi::cutLexis]⁠. See Examples.

The status variable

The statuses in the expanded output (lex.Cst and lex.Xst) are determined by using either only status or both status and entry.status. If entry.status = NULL, the status at entry is guessed according to the type of variable supplied via status: For numeric variables it will be zero, for factors the first level (levels(status)[1]) and otherwise the first unique value in alphabetical order (sort(unique(status))[1]).

Using numeric or factor status variables is strongly recommended. Logical expressions are also allowed (e.g. status = my_status != 0L) and are converted to integer internally.

Merging population hazard information

To enable computing relative/net survivals with ⁠[survtab]⁠ and ⁠[relpois]⁠, lexpand merges an appropriate population hazard data (pophaz) to the expanded data before dropping rows outside the specified time window (if drop = TRUE). pophaz must, for this reason, contain at a minimum the variables named agegroup, year, and haz. pophaz may contain additional variables to specify different population hazard levels in different strata; e.g. popmort includes sex. All the strata-defining variables must be present in the supplied data. lexpand will automatically detect variables with common names in the two datasets and merge using them.

Currently year must be an integer variable specifying the appropriate year. agegroup must currently also specify one-year age groups, e.g. popmort specifies 101 age groups of length 1 year. In both year and agegroup variables the values are interpreted as the lower bounds of intervals (and passed on to a cut call). The mandatory variable haz must specify the appropriate average rate at the person-year level; e.g. haz = -log(survProb) where survProb is a one-year conditional survival probability will be the correct hazard specification.

The corresponding pophaz population hazard value is merged by using the mid points of the records after splitting as reference values. E.g. if age=89.9 at the start of a 1-year interval, then the reference age value is 90.4 for merging. This way we get a "typical" population hazard level for each record.

Computing Pohar-Perme weights

If pp = TRUE, Pohar-Perme weights (the inverse of cumulative population survival) are computed. This will create the new pp variable in the expanded data. pp is a reserved name and lexpand throws exception if a variable with that name exists in data.

When a survival interval contains one or several rows per subject (e.g. due to splitting by the per scale), pp is cumulated from the beginning of the first record in a survival interval for each subject to the mid-point of the remaining time within that survival interval, and that value is given for every other record that a given person has within the same survival interval.

E.g. with 5 rows of duration 1/5 within a survival interval ⁠[0,1)]⁠, pp is determined for all records by a cumulative population survival from 0 to 0.5. The existing accuracy is used, so that the weight is cumulated first up to the end of the second row and then over the remaining distance to the mid-point (first to 0.4, then to 0.5). This ensures that more accurately merged population hazards are fully used.

Event not at end of follow-up & overlapping time lines

event may be used if the event indicated by status should occur at a time differing from exit. If event is defined, cutLexis is used on the data set after coercing it to the Lexis format and before splitting. Note that some values of event are allowed to be NA as with cutLexis to accommodate observations without an event occurring.

Additionally, setting overlapping = FALSE ensures that (irrespective of using event) the each subject defined by id only has one continuous time line instead of possibly overlapping time lines if there are multiple rows in data by id.

Aggregating

Certain analyses such as SIR/SMR calculations require tables of events and person-years by the unique combinations (interactions) of several variables. For this, aggre can be specified as a list of such variables (preferably factor variables but not mandatory) and any arbitrary functions of the variables at one's disposal. E.g.

aggre = list(sex, agegr = cut(dg_age, 0:100))

would tabulate events and person-years by sex and an ad-hoc age group variable. Every ad-hoc-created variable should be named.

fot, per, and age are special reserved variables which, when present in the aggre list, are output as categories of the corresponding time scale variables by using e.g.

cut(fot, breaks$fot, right=FALSE).

This only works if the corresponding breaks are defined in breaks or via "...". E.g.

aggre = list(sex, fot.int = fot) with

breaks = list(fot=0:5).

The output variable fot.int in the above example will have the lower limits of the appropriate intervals as values.

aggre as a named list will output numbers of events and person-years with the given new names as categorizing variable names, e.g. aggre = list(follow_up = fot, gender = sex, agegroup = age).

The output table has person-years (pyrs) and event counts (e.g. from0to1) as columns. Event counts are the numbers of transitions (lex.Cst != lex.Xst) or the lex.Xst value at a subject's last record (subject possibly defined by id).

If aggre.type = "unique" (alias "non-empty"), the above results are computed for existing combinations of expressions given in aggre, but also for non-existing combinations if aggre.type = "cartesian" (alias "full"). E.g. if a factor variable has levels ⁠"a", "b", "c"⁠ but the data is limited to only have levels ⁠"a", "b"⁠ present (more than zero rows have these level values), the former setting only computes results for ⁠"a", "b"⁠, and the latter also for "c" and any combination with other variables or expression given in aggre. In essence, "cartesian" forces also combinations of variables used in aggre that have no match in data to be shown in the result.

If aggre is not NULL and pophaz has been supplied, lexpand also aggregates the expected counts of events, which appears in the output data by the reserved name d.exp. Additionally, having pp = TRUE causes lexpand to also compute various Pohar-Perme weighted figures necessary for computing Pohar-Perme net survivals with ⁠[survtab_ag]⁠. This can be slow, so consider what is really needed. The Pohar-Perme weighted figures have the suffix .pp.

[0,1)]: R:0,1) [survtab_ag]: R:survtab_ag

Value

If aggre = NULL, returns a data.table or data.frame (depending on options("popEpi.datatable"); see ?popEpi) object expanded to accommodate split observations with time scales as fractional years and pophaz merged in if given. Population hazard levels in new variable pop.haz, and Pohar-Perme weights as new variable pp if requested.

If aggre is defined, returns a long-format data.table/data.frame with the variable pyrs (person-years), and variables for the counts of transitions in state or state at end of follow-up formatted fromXtoY, where X and Y are the states transitioned from and to, respectively. The data may also have the columns d.exp for expected numbers of cases and various Pohar-Perme weighted figures as identified by the suffix .pp; see Details.

Author(s)

Joonas Miettinen

See Also

⁠[Epi::Lexis]⁠, ⁠[popmort]⁠

Other splitting functions: splitLexisDT(), splitMulti()

Other aggregation functions: aggre(), as.aggre(), setaggre(), summary.aggre()

Examples

## prepare data for e.g. 5-year cohort survival calculation
x <- lexpand(sire, breaks=list(fot=seq(0, 5, by = 1/12)),
             birth = bi_date, entry = dg_date, exit = ex_date,
             status =  status != 0, pophaz=popmort)

## prepare data for e.g. 5-year "period analysis" for 2008-2012
BL <- list(fot = seq(0, 5, by = 1/12), per = c("2008-01-01", "2013-01-01"))
x <- lexpand(sire, breaks = BL,
             birth = bi_date, entry = dg_date, exit = ex_date,
             pophaz=popmort, status =  status != 0)

## aggregating
BL <- list(fot = 0:5, per = c("2003-01-01","2008-01-01", "2013-01-01"))
ag <- lexpand(sire, breaks = BL, status = status != 0,
             birth = bi_date, entry = dg_date, exit = ex_date,
              aggre=list(sex, period = per, surv.int = fot))

## aggregating even more
ag <- lexpand(sire, breaks = BL, status = status != 0,
              birth = bi_date, entry = dg_date, exit = ex_date,
              aggre=list(sex, period = per, surv.int = fot),
              pophaz = popmort, pp = TRUE)

## using "..."
x <- lexpand(sire, fot=0:5, status =  status != 0,
             birth = bi_date, entry = dg_date, exit = ex_date,
             pophaz=popmort)

x <- lexpand(sire, fot=0:5, status =  status != 0,
             birth = bi_date, entry = dg_date, exit = ex_date,
             aggre=list(sex, surv.int = fot))

## using the "event" argument: it just places the transition to given "status"
## at the "event" time instead of at the end, if possible using cutLexis
x <- lexpand(sire, status = status, event = dg_date,
             birth = bi_date, entry = dg_date, exit = ex_date,)

## aggregating with custom "event" time
## (the transition to status is moved to the "event" time)
x <- lexpand(sire, status = status, event = dg_date,
             birth = bi_date, entry = dg_date, exit = ex_date,
             per = 1970:2014, age = c(0:100,Inf),
             aggre = list(sex, year = per, agegroup = age))

Description

Plot SIR spline lines with R base graphics

Usage

## S3 method for class 'sirspline'
lines(x, conf.int = TRUE, print.levels = NA, select.spline, ...)

Arguments

Details

In lines.sirspline most of graphical parameters is user adjustable. Desired spline variable can be selected with select.spline and only one can be plotted at a time. The spline variable can include several levels, e.g. gender (these are the levels of print from sirspline). All levels are printed by default, but a specific level can be selected using argument print.levels. Printing the levels separately enables e.g. to give different colours for each level.

Value

Always returns NULL invisibly. This function is called for its side effects.

Author(s)

Matti Rantanen

See Also

Other sir functions: plot.sirspline(), sir(), sir_exp(), sir_ratio(), sirspline()

Description

Plots the observed (with extrapolation) and expected survival curves for all strata in an object created by ⁠[survmean]⁠

Usage

## S3 method for class 'survmean'
lines(x, ...)

Arguments

Details

This function is intended to be a workhorse for ⁠[plot.survmean]⁠. If you want finer control over the plotted curves, extract the curves from the survmean output using

attr(x, "curves")

where x is a survmean object.

Value

Always returns NULL invisibly. This function is called for its side effects.

Author(s)

Joonas Miettinen

See Also

Other survmean functions: Surv(), plot.survmean(), survmean()

Description

Plot lines from a survtab object

Usage

## S3 method for class 'survtab'
lines(x, y = NULL, subset = NULL, conf.int = TRUE, col = NULL, lty = NULL, ...)

Arguments

Value

Always returns NULL invisibly. This function is called for its side effects.

Author(s)

Joonas Miettinen

See Also

Other survtab functions: Surv(), plot.survtab(), print.survtab(), summary.survtab(), survtab(), survtab_ag()

Examples

data(sire)
data(sibr)
si <- rbind(sire, sibr)
si$period <- cut(si$dg_date, as.Date(c("1993-01-01", "2004-01-01", "2013-01-01")), right = FALSE)
si$cancer <- c(rep("rectal", nrow(sire)), rep("breast", nrow(sibr)))
x <- lexpand(si, birth = bi_date, entry = dg_date, exit = ex_date,
             status = status %in% 1:2,
             fot = 0:5, aggre = list(cancer, period, fot))
st <- survtab_ag(fot ~ cancer + period, data = x,
                 surv.method = "lifetable", surv.type = "surv.obs")

plot(st, "surv.obs", subset = cancer == "breast", ylim = c(0.5, 1), col = "blue")
lines(st, "surv.obs", subset = cancer == "rectal", col = "red")

## or
plot(st, "surv.obs", col = c(2,2,4,4), lty = c(1, 2, 1, 2))

Description

selects lowest values of each factor after cut() based on that the value starts from index 2 and end in comma ",".

Usage

lower_bound(cut)

Arguments

Value

A numeric vector.

Author(s)

Matti Rantanen

Description

ltable makes use of data.table capabilities to tabulate frequencies or arbitrary functions of given variables into a long format data.table/data.frame. expr.by.cj is the equivalent for more advanced users.

Usage

ltable(
  data,
  by.vars = NULL,
  expr = list(obs = .N),
  subset = NULL,
  use.levels = TRUE,
  na.rm = FALSE,
  robust = TRUE
)

expr.by.cj(
  data,
  by.vars = NULL,
  expr = list(obs = .N),
  subset = NULL,
  use.levels = FALSE,
  na.rm = FALSE,
  robust = FALSE,
  .SDcols = NULL,
  enclos = parent.frame(1L),
  ...
)

Arguments

Details

Returns expr for each unique combination of given by.vars.

By default makes use of any and all ⁠[levels]⁠ present for each variable in by.vars. This is useful, because even if a subset of the data does not contain observations for e.g. a specific age group, those age groups are nevertheless presented in the resulting table; e.g. with the default expr = list(obs = .N) all age group levels are represented by a row and can have obs = 0.

The function differs from the vanilla ⁠[table]⁠ by giving a long format table of values regardless of the number of by.vars given. Make use of e.g. ⁠[cast_simple]⁠ if data needs to be presented in a wide format (e.g. a two-way table).

The rows of the long-format table are effectively Cartesian products of the levels of each variable in by.vars, e.g. with by.vars = c("sex", "area") all levels of area are repeated for both levels of sex in the table.

The expr allows the user to apply any function(s) on all levels defined by by.vars. Here are some examples:

• .N or list(.N) is a function used inside a data.table to calculate counts in each group

• list(obs = .N), same as above but user assigned variable name

• list(sum(obs), sum(pyrs), mean(dg_age)), multiple objects in a list

• list(obs = sum(obs), pyrs = sum(pyrs)), same as above with user defined variable names

If use.levels = FALSE, no levels information will be used. This means that if e.g. the agegroup variable is a factor and has 18 levels defined, but only 15 levels are present in the data, no rows for the missing levels will be shown in the table.

na.rm simply drops any rows from the resulting table where any of the by.vars values was NA.

Value

A data.table of statistics (e.g. counts) stratified by the columns defined in by.vars.

Functions

• expr.by.cj(): Somewhat more streamlined ltable with defaults for speed. Explicit determination of enclosing environment of data.

Author(s)

Joonas Miettinen, Matti Rantanen

See Also

⁠[table]⁠, ⁠[cast_simple]⁠, ⁠[data.table::melt]⁠

Examples

data("sire", package = "popEpi")
sr <- sire
sr$agegroup <- cut(sr$dg_age, breaks=c(0,45,60,75,85,Inf))
## counts by default
ltable(sr, "agegroup")

## any expression can be given
ltable(sr, "agegroup", list(mage = mean(dg_age)))
ltable(sr, "agegroup", list(mage = mean(dg_age), vage = var(dg_age)))

## also returns levels where there are zero rows (expressions as NA)
ltable(sr, "agegroup", list(obs = .N,
                            minage = min(dg_age),
                            maxage = max(dg_age)),
       subset = dg_age < 85)

#### expr.by.cj
expr.by.cj(sr, "agegroup")

## any arbitrary expression can be given
expr.by.cj(sr, "agegroup", list(mage = mean(dg_age)))
expr.by.cj(sr, "agegroup", list(mage = mean(dg_age), vage = var(dg_age)))

## only uses levels of by.vars present in data
expr.by.cj(sr, "agegroup", list(mage = mean(dg_age), vage = var(dg_age)),
           subset = dg_age < 70)

## .SDcols trick
expr.by.cj(sr, "agegroup", lapply(.SD, mean),
           subset = dg_age < 70, .SDcols = c("dg_age", "status"))

Description

Mean population counts in Finland year, sex, and age group.

Format

data.table with columns

• sex gender coded as male, female (0, 1)

• year calendar year 1981-2016

• agegroup - coded 0 to 100; one-year age groups

• meanpop the mean population count; that is, the mean of the annual population counts of two consecutive years; e.g. for 1990 meanpop is the mean of population counts for 1990 and 1991 (counted at 1990-01-01 and 1991-01-01, respectively)

Source

Statistics Finland

See Also

Other popEpi data: ICSS, popmort, sibr, sire, stdpop101, stdpop18

Description

Given a ⁠data.table DT⁠, replaces any NA values in the variables given in vars in DT. Takes a copy of the original data and returns the modified copy.

Usage

na2zero(DT, vars = NULL)

Arguments

Details

Given a data.table object, converts NA values to numeric (double) zeros for all variables named in vars or all variables if vars = NULL.

Value

A copy of DT where NA values have been replaced with zero.

Author(s)

Joonas Miettinen

Description

Plot rate estimates with confidence intervals lines using R base graphics

Usage

## S3 method for class 'rate'
plot(x, conf.int = TRUE, eps = 0.2, left.margin, xlim, ...)

Arguments

Details

This is limited explanatory tool but most graphical parameters are user adjustable.

Value

Always returns NULL invisibly. This function is called for its side effects.

Author(s)

Matti Rantanen

Description

Plot SIR estimates with error bars

Usage

## S3 method for class 'sir'
plot(
  x,
  conf.int = TRUE,
  ylab,
  xlab,
  xlim,
  main,
  eps = 0.2,
  abline = TRUE,
  log = FALSE,
  left.margin,
  ...
)

Arguments

Details

Plot SIR estimates and confidence intervals

• univariate - plots SIR with univariate confidence intervals

• model - plots SIR with Poisson modelled confidence intervals

Customize Normal plot parameters can be passed to plot. These can be a vector when plotting error bars:

• pch - point type

• lty - line type

• col - line/point colour

• lwd - point/line size

Tips for plotting splines It's possible to use plot to first draw the confidence intervals using specific line type or colour and then plotting again the estimate using lines(... , conf.int = FALSE) with different settings. This works only when plot.type is 'splines'.

Value

Always returns NULL invisibly. This function is called for its side effects.

Author(s)

Matti Rantanen

See Also

⁠[sir]⁠, ⁠[sirspline]⁠

Examples

# Plot SIR estimates
# plot(sir.by.gender, col = c(4,2), log=FALSE, eps=0.2, lty=1, lwd=2, pch=19,
#      main = 'SIR by gender', abline=TRUE)

Description

Plot SIR splines using R base graphics.

Usage

## S3 method for class 'sirspline'
plot(x, conf.int = TRUE, abline = TRUE, log = FALSE, type, ylab, xlab, ...)

Arguments

Details

In plot.sirspline almost every graphical parameter are user adjustable, such as ylim, xlim. plot.sirsplines calls lines.splines to add lines.

The plot axis without lines can be plotted using option type = 'n'. On top of the frame it's then possible to add a grid, abline or text before plotting the lines (see: sirspline).

Value

Always returns NULL invisibly. This function is called for its side effects.

Author(s)

Matti Rantanen

See Also

Other sir functions: lines.sirspline(), sir(), sir_exp(), sir_ratio(), sirspline()

Description

Plots the observed (with extrapolation) and expected survival curves for all strata in an object created by ⁠[survmean]⁠

Usage

## S3 method for class 'survmean'
plot(x, ...)

Arguments

Details

For examples see ⁠[survmean]⁠. This function is intended only for graphically inspecting that the observed survival curves with extrapolation and the expected survival curves have been sensibly computed in survmean.

If you want finer control over the plotted curves, extract the curves from the survmean output using

attr(x, "curves")

where x is a survmean object.

Value

Always returns NULL invisibly. This function is called for its side effects.

Author(s)

Joonas Miettinen

See Also

Other survmean functions: Surv(), lines.survmean(), survmean()

Description

Plotting for survtab objects

Usage

## S3 method for class 'survtab'
plot(
  x,
  y = NULL,
  subset = NULL,
  conf.int = TRUE,
  col = NULL,
  lty = NULL,
  ylab = NULL,
  xlab = NULL,
  ...
)

Arguments

Value

Always returns NULL invisibly. This function is called for its side effects.

Author(s)

Joonas Miettinen

See Also

Other survtab functions: Surv(), lines.survtab(), print.survtab(), summary.survtab(), survtab(), survtab_ag()

Examples

data(sire)
data(sibr)
si <- rbind(sire, sibr)
si$period <- cut(si$dg_date, as.Date(c("1993-01-01", "2004-01-01", "2013-01-01")), right = FALSE)
si$cancer <- c(rep("rectal", nrow(sire)), rep("breast", nrow(sibr)))
x <- lexpand(si, birth = bi_date, entry = dg_date, exit = ex_date,
             status = status %in% 1:2,
             fot = 0:5, aggre = list(cancer, period, fot))
st <- survtab_ag(fot ~ cancer + period, data = x,
                 surv.method = "lifetable", surv.type = "surv.obs")

plot(st, "surv.obs", subset = cancer == "breast", ylim = c(0.5, 1), col = "blue")
lines(st, "surv.obs", subset = cancer == "rectal", col = "red")

## or
plot(st, "surv.obs", col = c(2,2,4,4), lty = c(1, 2, 1, 2))

Description

Computes confidence intervals for Poisson rates

Usage

poisson.ci(x, pt = 1, conf.level = 0.95)

Arguments

Value

A data.frame with columns

• x: arg x

• pt: arg pt

• rate: result of x / pt

• lower: lower bound of CI

• upper: upper bound of CI

• conf.level: arg conf.level

Author(s)

epitools

Examples

poisson.ci(x = 4, pt = 5, conf.level = 0.95)

Description

Several functions in popEpi make use of population or expected hazards in computing the intended estimates (e.g. ⁠[survtab]⁠). This document explains using such data sets in this package.

Details

Population hazard data sets (pophaz for short) in popEpi should be data.frames in the "long" format where one of the columns must be named haz (for hazard), and other columns define the values or levels in variables relating to subjects in your data. For example, ⁠[popmort]⁠ contains Finnish population mortality hazards by sex, calendar year, and 1-year age group.

The names of the columns should match to the names of the variables that you have in your subject-level data. Time variables in your pophaz may also correspond to Lexis time scales; see ⁠[survtab]⁠.

Any time variables (as they usually have) should be coded consistently: When using fractional years in your data, the time variables in your pophaz must also be coded in fractional years. When using e.g. Dates in your data, ensure that the pophaz time variables are coded at the level of days (or Dates for calendar time).

The haz variable in your pophaz should also be coded consistently with the used time variables. E.g. haz values in life-tables reported as deaths per person-year should be multiplied by 365.25 when using day-level time variables. Typically you'll have calendar time and age expressed in years, which means haz should be expressed as the number of deaths per person-year.

If you have your population hazards in a ratetable object usable by functions in survival and relsurv, you may transform them to long-format data.frames using ⁠[ratetable_to_long_dt]⁠. Ensure, however, that the created haz column is coded at the right level (events per days or years typically).

National statistical institutions, the WHO, and e.g. the Human Life-Table Database supply life-table data.

Author(s)

Joonas Miettinen

Description

Population mortality rates in Finland 1951 - 2013 in 101 age groups and by gender. This is an example of a population hazard table as used in popEpi; for the general help page, see ⁠[pophaz]⁠.

Format

data.table with columns

• sex gender coded as male, female (0, 1)

• year calendar year

• agegroup - coded 0 to 100; one-year age groups

• haz the average population mortality rate per person-year (d/(pyrs), where d is the number of deaths and pyrs is the person-years)

Source

Statistics Finland

See Also

⁠[pophaz]⁠

Other popEpi data: ICSS, meanpop_fi, sibr, sire, stdpop101, stdpop18

Description

prepExpo uses a Lexis object of periods of exposure to fill gaps between the periods and overall entry and exit times without accumulating exposure time in periods of no exposure, and splits the result if requested.

Usage

prepExpo(
  lex,
  freezeScales = "work",
  cutScale = "per",
  entry = min(get(cutScale)),
  exit = max(get(cutScale)),
  by = "lex.id",
  breaks = NULL,
  freezeDummy = NULL,
  subset = NULL,
  verbose = FALSE,
  ...
)

Arguments

Details

prepExpo is a convenience function for the purpose of eventually aggregating person-time and events in categories of not only normally progressing ⁠[Epi::Lexis]⁠ time scales but also some time scales which should not progress sometimes. For example a person may work at a production facility only intermittently, meaning exposure time (to work-related substances for example) should not progress outside of periods of work. This allows for e.g. a correct aggregation of person-time and events by categories of cumulative time of exposure.

Given an ⁠[Epi::Lexis]⁠ object containing rows (time lines) where a subject is exposed to something (and NO periods without exposure), fills any gaps between exposure periods for each unique combination of by and the subject-specific "ultimate" entry and exit times, "freezes" the cumulative exposure times in periods of no exposure, and splits data using breaks passed to ⁠[splitMulti]⁠ if requested. Results in a (split) Lexis object where freezeScales do not progress in time periods where no exposure was recorded in lex.

This function assumes that entry and exit arguments are the same for each row within a unique combination of variables named in by. E.g. with by = "lex.id" only each lex.id has a unique value for entry and exit at most.

The supplied breaks split the data using splitMulti, with the exception that breaks supplied concerning any frozen time scales ONLY split the rows where the time scales are not frozen. E.g. with freezeScales = "work", breaks = list(work = 0:10, cal = 1995:2010) splits all rows over "cal" but only non-frozen rows over "work".

Only supports frozen time scales that advance and freeze contemporaneously: e.g. it would not currently be possible to take into account the cumulative time working at a facility and the cumulative time doing a single task at the facility, if the two are not exactly the same. On the other hand one might use the same time scale for different exposure types, supply them as separate rows, and identify the different exposures using a dummy variable.

Value

Returns a Lexis object that has been split if breaks is specified. The resulting time is also a data.table if options("popEpi.datatable") == TRUE (see: ?popEpi)

Description

Function(s) to compute prevalence statistics.

Usage

prev_lexis(
  lexis,
  observation_time_points,
  stratum_breaks,
  aggre_by = NULL,
  subset = NULL,
  merge_dt = NULL,
  merge_optional_args = NULL
)

Arguments

Value

popEpi::prev_lexis

Returns a data.table with

• Stratifying columns defined via aggre_by (if any),

• Stratifying time scale columns defined via stratum_breaks (if any),

• The observation time scale column defined via observation_time_points,

• n_prev, the number of subjects in follow-up, and

• n_prev_eff, the above plus the "extrapolated" number in follow-up, if this was requested.

Functions

popEpi::prev_lexis

popEpi::prev_lexis can be used to compute numbers of (potentially effective numbers of) subjects remaining in follow-up at arbitrary time points. It performs the following steps:

• For each observation time point in observation_time_points[[1]]:

  • Calls ⁠[lexis_split_merge_aggregate_by_stratum]⁠ with the lexis, subset, aggre_by supplied to popEpi::prev_lexis and with aggre_exprs = list(n_prev = quote(.N)), and breaks = stratum_breaks. This produces a table stratified by aggre_by and all time scales used in stratum_breaks. The only value column at this point is n_prev, the number of subjects in follow-up at the given observation time point.

  • If !is.null(merge_dt), collect subjects in lexis who were censored before the current observation time point. This is defined as those in lexis who have lexis[["lex.Cst"]] == lexis[["lex.Xst"]] and lexis[[obs_ts_col_nm]] + lexis[["lex.dur"]] < obs_tp, where obs_ts_col_nm = names(observation_time_points) and obs_tp is the current observation time point. These collected subjects are the ones we need to extrapolate to the current observation time point.

  • If there are no observations that need to be extrapolated then we end the extrapolation step early and simply use n_prev_eff = n_prev. Otherwise proceed as described below.

  • If inherits(merge_dt, "list"), collect arguments from merge_dt to call ⁠[surv_lexis]⁠. Arguments lexis, subset, and aggre_by are assigned internally to the respective arguments supplied to prev_lexis except subset additionally is used to detect only those cases who entered follow-up (any time) before the current observation time point. It may be worth emphasizing that lexis really is the input argument lexis and not only those whose follow-up we need to "extrapolate". We also set estimators = "S_ch". Any arguments that are passed via merge_dt override even the internally set arguments. The arguments not set internally or supplied by the user via merge_dt make use of the defaults of ⁠[surv_lexis]⁠. E.g. merge_dt = list(aggre_by = NULL, breaks = list(ts_fut = futs)) are both passed to ⁠[surv_lexis]⁠ despite what aggre_by was for prev_lexis.

  • If inherits(merge_dt, "list"), we also interpolate survival estimates in the table we have just produced so that there are always 1000 intervals starting from zero and ending where the last survival estimate is available. This is peformed to ensure that the survival estimates we merge will be available at a sufficient "resolution". For instance, even if the survival estimates are for one-year intervals due to sparsity, we get a reasonable survival estimate for a subject who was censored at 1.01 — from somewhere close to 1.01 instead from 2.0.

  • Merge (for the first time) merge_dt with the collected subjects at the original exit time of each subject. This yields the survival probability for each subject at exit, in math S(t_e) (starting from zero — delayed entry is not supported).

  • Merge merge_dt for the second time, this time at the current prevalence observation time point such as at ts_cal = 2009.999. In math this is S(t_p) where t_p is the prevalence observation time point.

  • With both S(t_e) and S(t_p) available, our "extrapolated" or "effective" number of being in follow-up is between zero and one for each subject and defined simply as the conditional survival up to t_p starting from t_e, S(t_p|t_e) = S(t_p) / S(t_e). E.g. S(t_p) / S(t_e) = 0.8 / 0.9 ~ 0.8888889.

  • Call ⁠[lexis_split_merge_aggregate_by_stratum]⁠ for the second time, this time with the subjects collected for extrapolation, and sum the number of extrapolated subjects in follow-up into a table with the identical stratification as the one created before.

  • Add column n_prev_eff as the sum of the number of extrapolated subjects and the original n_prev into the first table we created.

• Collect the observation time point-specific results into one big table and return it.

Description

Print method function for aggre objects; see ⁠[as.aggre]⁠ and ⁠[aggre]⁠.

Usage

## S3 method for class 'aggre'
print(x, subset = NULL, ...)

Arguments

Value

Always returns NULL invisibly. This function is called for its side effects.

Author(s)

Joonas Miettinen

Description

Print method function for rate objects; see ⁠[rate]⁠.

Usage

## S3 method for class 'rate'
print(x, subset = NULL, ...)

Arguments

Value

Always returns NULL invisibly. This function is called for its side effects.

Author(s)

Matti Rantanen

Description

Print method function for survtab objects; see ⁠[survtab_ag]⁠.

Usage

## S3 method for class 'survtab'
print(x, subset = NULL, ...)

Arguments

Value

Always returns NULL invisibly. This function is called for its side effects.

Author(s)

Joonas Miettinen

See Also

Other survtab functions: Surv(), lines.survtab(), plot.survtab(), summary.survtab(), survtab(), survtab_ag()

Description

rate calculates adjusted rates using preloaded weights data or user specified weights.

Usage

rate(
  data,
  obs = NULL,
  pyrs = NULL,
  print = NULL,
  adjust = NULL,
  weights = NULL,
  subset = NULL
)

Arguments

Details

Input data needs to be in aggregated format with observations and person-years. For individual data use ⁠[lexpand]⁠, or ⁠[ltable]⁠ and merge person-years manually.

The confidence intervals are based on the normal approximation of the logarithm of the rate. The variance of the log rate that is used to derive the confidence intervals is derived using the delta method.

Value

Returns a data.table with observations, person-years, rates and adjusted rates, if available. Results are stratified by print. Adjusted rates are identified with suffix .adj and .lo and .hi are for confidence intervals lower and upper 95% bounds, respectively. The prefix SE. stands for standard error.

Author(s)

Matti Rantanen, Joonas Miettinen

See Also

⁠[lexpand]⁠, ⁠[ltable]⁠