eha: Event History Analysis

Description

Parametric proportional hazards fitting with left truncation and right censoring for common families of distributions, piecewise constant hazards, and discrete models. Parametric accelerated failure time models for left truncated and right censored data. Proportional hazards models for tabular and register data. Sampling of risk sets in Cox regression, selections in the Lexis diagram, bootstrapping. Broström (2022) doi:10.1201/9780429503764.

Details

Eha enhances the recommended survival package in several ways, see the description. The main applications in mind are demography and epidemiology. For standard Cox regression analysis the function coxph in survival is still recommended. The function coxreg in eha in fact calls coxph for the standard kind of analyses.

Author(s)

Maintainer: Göran Broström goran.brostrom@umu.se

Other contributors:

• Jianming Jin [contributor]

References

Broström, G. (2012). Event History Analysis with R, Chapman and Hall/CRC Press, Boca Raton, FL.

See Also

Useful links:

• https://ehar.se/r/eha/

• Report bugs at https://github.com/goranbrostrom/eha/issues

Accelerated Failure Time Regression

Description

The accelerated failure time model with parametric baseline hazard(s). Allows for stratification with different scale and shape in each stratum, and left truncated and right censored data.

Usage

aftreg(
  formula = formula(data),
  data = parent.frame(),
  na.action = getOption("na.action"),
  dist = "weibull",
  init,
  shape = 0,
  id,
  param = c("lifeAcc", "lifeExp"),
  control = list(eps = 1e-08, maxiter = 20, trace = FALSE),
  singular.ok = TRUE,
  model = FALSE,
  x = FALSE,
  y = TRUE
)

Arguments

Details

The parameterization is different from the one used by survreg, when param = "lifeAcc". The result is then true acceleration of time. Then the model is

S(t; a, b, \beta, z) = S_0((t / \exp(b - z\beta))^{\exp(a)})

where S_0 is some standardized survivor function. The baseline parameters a and b are log shape and log scale, respectively. This is for the default parametrization. With the lifeExp parametrization, some signs are changed:

b - z beta

is changed to

b + z beta

. For the Gompertz distribution, the base parametrization is canonical, a necessity for consistency with the shape/scale paradigm (this is new in 2.3).

Value

A list of class "aftreg" with components

Author(s)

Göran Broström

See Also

coxreg, phreg, survreg

Examples

data(mort)
aftreg(Surv(enter, exit, event) ~ ses, param = "lifeExp", data = mort)

Parametric proportional hazards regression

Description

This function is called by aftreg, but it can also be directly called by a user.

Usage

aftreg.fit(X, Y, dist, param, strata, offset, init, shape, id, control, pfixed)

Arguments

Details

See aftreg for more detail.

Value

Author(s)

Göran Broström

See Also

aftreg

Age cut of survival data

Description

For a given age interval, each spell is cut to fit into the given age interval.

Usage

age.window(dat, window, surv = c("enter", "exit", "event"))

Arguments

Details

The window must be in the order (begin, end)

Value

A data frame of the same form as the input data frame, but 'cut' as desired. Intervals exceeding window[2] will be given event = 0. If the selection gives an empty result, NULL is returned, with no warning.

Author(s)

Göran Broström

See Also

cal.window, coxreg, aftreg

Examples

dat <- data.frame(enter = 0, exit = 5.731, event = 1, x = 2)
window <- c(2, 5.3)
dat.trim <- age.window(dat, window)  

Calendar time cut of survival data

Description

For a given time interval, each spell is cut so that it fully lies in the given time interval

Usage

cal.window(dat, window, surv = c("enter", "exit", "event", "birthdate"))

Arguments

Details

The window must be in the order (begin, end)

Value

A data frame of the same form as the input data frame, but 'cut' as desired. Intervals exceeding window[2] will be given event = 0

Author(s)

Göran Broström

See Also

age.window, coxreg, aftreg

Examples

dat <- data.frame(enter = 0, exit = 5.731, event = 1,
birthdate = 1962.505, x = 2)
window <- c(1963, 1965)
dat.trim <- cal.window(dat, window)  

Graphical goodness-of-fit test

Description

Comparison of the cumulative hazards functions for a semi-parametric and a parametric model.

Usage

check.dist(sp, pp, main = NULL, col = 1:2, lty = 1:2, printLegend = TRUE)

Arguments

Details

For the moment only a graphical comparison. The arguments sp and pp may be swapped.

Value

No return value.

Author(s)

Göran Broström

See Also

coxreg and phreg.

Examples

data(mort)
oldpar <- par(mfrow = c(2, 2))
fit.cr <- coxreg(Surv(enter, exit, event) ~ ses, data = mort)
fit.w <- phreg(Surv(enter, exit, event) ~ ses, data = mort)
fit.g <- phreg(Surv(enter, exit, event) ~ ses, data = mort,
dist = "gompertz")
fit.ev <- phreg(Surv(enter, exit, event) ~ ses, data = mort,
dist = "ev")
check.dist(fit.cr, fit.w)
check.dist(fit.cr, fit.g)
check.dist(fit.cr, fit.ev)
par(oldpar)

Check the integrity of survival data.

Description

Check that exit occurs after enter, that spells from an individual do not overlap, and that each individual experiences at most one event.

Usage

check.surv(enter, exit, event, id = NULL, eps = 1e-08)

Arguments

Details

Interval lengths must be strictly positive.

Value

A vector of id's for the insane individuals. Of zero length if no errors.

Author(s)

Göran Broström

See Also

join.spells, coxreg, aftreg

Examples

xx <- data.frame(enter = c(0, 1), exit = c(1.5, 3), event = c(0, 1), id =
c(1,1))
check.surv(xx$enter, xx$exit, xx$event, xx$id)

Child mortality, Skellefteå, Sweeden 1850–1900.

Description

Children born in Skellefteå, Sweden, 1850-1884, are followed fifteen years or until death or out-migration.

Usage

data(child)

Format

A data frame with 26855 children born 1850-1884.

id

An identification number.

m.id

Mother's id.

sex

Sex.

socBranch

Working branch of family (father).

birthdate

Birthdate.

enter

Start age of follow-up, always zero.

exit

Age of departure, either by death or emigration.

event

Type of departure, death = 1, right censoring = 0.

illeg

Born out of marriage ("illegitimate")?

m.age

Mother's age.

Details

The Skellefteå region is a large region in the northern part of Sweden.

Source

Data originate from the Centre for Demographic and Ageing Research, Umeå University, Umeå, Sweden, https://www.umu.se/en/centre-for-demographic-and-ageing-research/.

Examples

fit <- coxreg(Surv(enter, exit, event) ~ sex + socBranch, data = child, coxph = TRUE)
summary(fit)

Graphical comparison of cumulative hazards

Description

Comparison of the estimated baseline cumulative hazards functions for two survival models.

Usage

compHaz(
  fit1,
  fit2,
  main = NULL,
  lty = 1:2,
  col = c("red", "blue"),
  printLegend = TRUE
)

Arguments

Value

No return value.

Author(s)

Göran Broström

See Also

hazards, coxreg, and phreg.

Examples

fit.cr <- coxreg(Surv(enter, exit, event) ~ sex, data = oldmort)
fit.w <- phreg(Surv(enter, exit, event) ~ sex, data = oldmort)
compHaz(fit.cr, fit.w)

Loglihood function (partial likelihood) of a Cox regression

Description

Calculates minus the log likelihood function and its first and second order derivatives for data from a Cox regression model. It is used by coxreg if the argument coxph = FALSE

Usage

coxfunk(beta, X, offset, rs, what = 2)

Arguments

Details

Note that the function returns log likelihood, score vector and minus hessian, i.e. the observed information. The model is

Value

A list with components

Author(s)

Göran Broström

See Also

coxreg

Cox regression

Description

Performs Cox regression with some special attractions, especially sampling of risksets and the weird bootstrap.

Usage

coxreg(formula = formula(data), data = parent.frame(), weights,
subset, t.offset, na.action = getOption("na.action"), init = NULL, method =
c("efron", "breslow", "mppl", "ml"), control = list(eps = 1e-08, maxiter =
25, trace = FALSE), singular.ok = TRUE, model = FALSE, center = NULL, x =
FALSE, y = TRUE, hazards = NULL, boot = FALSE, efrac = 0, geometric = FALSE,
rs = NULL, frailty = NULL, max.survs = NULL, coxph = TRUE)

Arguments

Details

The default method, efron, and the alternative, breslow, are both the same as in coxph in package survival. The methods mppl and ml are maximum likelihood, discrete-model, based.

Value

A list of class c("coxreg", "coxph") with components

Warning

The use of rs is dangerous, see note. It can however speed up computing time considerably for huge data sets.

Note

This function starts by creating risksets, if no riskset is supplied via rs, with the aid of risksets. Supplying output from risksets via rs fails if there are any NA's in the data! Note also that it depends on stratification, so rs contains information about stratification. Giving another strata variable in the formula is an error. The same is ok, for instance to supply stratum interactions.

Author(s)

Göran Broström

References

Broström, G. and Lindkvist, M. (2008). Partial partial likelihood. Communications in Statistics: Simulation and Computation 37:4, 679-686.

See Also

coxph, risksets

Examples

 dat <- data.frame(time=  c(4, 3,1,1,2,2,3),
                status=c(1,1,1,0,1,1,0),
                x=     c(0, 2,1,1,1,0,0),
                sex=   c(0, 0,0,0,1,1,1))
 coxreg( Surv(time, status) ~ x + strata(sex), data = dat) #stratified model
 # Same as:
 rs <- risksets(Surv(dat$time, dat$status), strata = dat$sex)
 coxreg( Surv(time, status) ~ x, data = dat, rs = rs) #stratified model
 

Cox regression

Description

Called by coxreg, but a user can call it directly.

Usage

coxreg.fit(
  X,
  Y,
  rs,
  weights,
  t.offset = NULL,
  strats,
  offset,
  init,
  max.survs,
  method = "efron",
  boot = FALSE,
  efrac = 0,
  calc.martres = TRUE,
  control,
  verbose = TRUE,
  calc.hazards = NULL,
  center = NULL
)

Arguments

Details

rs is dangerous to use when NA's are present.

Value

A list with components

Note

It is the user's responsibility to check that indata is sane.

Author(s)

Göran Broström

See Also

coxreg, risksets

Examples

 X <- as.matrix(data.frame(
                x=     c(0, 2,1,4,1,0,3),
                sex=   c(1, 0,0,0,1,1,1)))
 time <- c(1,2,3,4,5,6,7)
 status <- c(1,1,1,0,1,1,0)
 stratum <- rep(1, length(time))

 coxreg.fit(X, Surv(time, status), strats = stratum, max.survs = 6,
     control = list(eps=1.e-4, maxiter = 10, trace = FALSE))

Cox regression

Description

Performs Cox regression with some special attractions, especially sampling of risksets and the weird bootstrap.

Usage

coxreg2(formula = formula(data), data = parent.frame(), weights,
subset, t.offset, na.action = getOption("na.action"), init = NULL, method =
c("efron", "breslow", "mppl", "ml"), control = list(eps = 1e-08, maxiter =
25, trace = FALSE), singular.ok = TRUE, model = FALSE, center = NULL, x =
FALSE, y = TRUE, hazards = NULL, boot = FALSE, efrac = 0, geometric = FALSE,
rs = NULL, frailty = NULL, max.survs = NULL, coxph = TRUE)

Arguments

Details

The default method, efron, and the alternative, breslow, are both the same as in coxph in package survival. The methods mppl and ml are maximum likelihood, discrete-model, based.

Value

A list of class c("coxreg", "coxph") with components

Warning

The use of rs is dangerous, see note. It can however speed up computing time considerably for huge data sets.

Note

This function starts by creating risksets, if no riskset is supplied via rs, with the aid of risksets. Supplying output from risksets via rs fails if there are any NA's in the data! Note also that it depends on stratification, so rs contains information about stratification. Giving another strata variable in the formula is an error. The same is ok, for instance to supply stratum interactions.

Author(s)

Göran Broström

References

Broström, G. and Lindkvist, M. (2008). Partial partial likelihood. Communications in Statistics: Simulation and Computation 37:4, 679-686.

See Also

coxph, risksets

Examples

 dat <- data.frame(time=  c(4, 3,1,1,2,2,3),
                status=c(1,1,1,0,1,1,0),
                x=     c(0, 2,1,1,1,0,0),
                sex=   c(0, 0,0,0,1,1,1))
 coxreg( Surv(time, status) ~ x + strata(sex), data = dat) #stratified model
 # Same as:
 rs <- risksets(Surv(dat$time, dat$status), strata = dat$sex)
 coxreg( Surv(time, status) ~ x, data = dat, rs = rs) #stratified model
 

Creates a minimal representation of a data frame.

Description

Given a data frame with a defined response variable, this function creates a unique representation of the covariates in the data frame, vector (matrix) of responses, and a pointer vector, connecting the responses with the corresponding covariates.

Usage

cro(dat, response = 1)

Arguments

Details

The rows in the data frame are converted to text strings with paste and compared with match.

Value

A list with components

Note

This function is based on suggestions by Anne York and Brian Ripley.

Author(s)

Göran Broström

See Also

match, paste

Examples

dat <- data.frame(y = c(1.1, 2.3, 0.7), x1 = c(1, 0, 1), x2 = c(0, 1, 0))
cro(dat)

The EV Distribution

Description

Density, distribution function, quantile function, hazard function, cumulative hazard function, and random generation for the EV distribution with parameters shape and scale.

Usage

dEV(x, shape = 1, scale = 1, log = FALSE) 
pEV(q, shape = 1, scale = 1, lower.tail = TRUE, log.p = FALSE) 
qEV(p, shape = 1, scale = 1, lower.tail = TRUE, log.p = FALSE) 
hEV(x, shape = 1, scale = 1, log = FALSE) 
HEV(x, shape = 1, scale = 1, log.p = FALSE) 
rEV(n, shape = 1, scale = 1)

Arguments

Details

The EV distribution with scale parameter a and shape parameter \sigma has hazard function given by

h(x) = (b/\sigma)(x/\sigma)^(b-1)\exp((x / \sigma)^b)

for x \ge 0.

Value

dEV gives the density, pEV gives the distribution function, qEV gives the quantile function, hEV gives the hazard function, HEV gives the cumulative hazard function, and rEV generates random deviates.

Invalid arguments will result in return value NaN, with a warning.

Marital fertility nineteenth century

Description

Birth intervals for married women with at least one birth, 19th northern Sweden

Usage

data(fert)

Format

A data frame with 12169 observations the lengths (in years) of birth intervals for 1859 married women with at least one birth. The first interval (parity = 0) is the interval from marriage to first birth.

id

Personal identification number for mother.

parity

Time order of birth interval for the present mother. The interval with parity = 0 is the first, from marriage to first birth.

age

The age of mother at start of interval.

year

The calendar year at start of interval.

next.ivl

The length of the coming time interval.

event

An indicator for whether the next.ivl ends in a new birth (event = 1) or is right censored (event = 0). Censoring occurs when the woman ends her fertility period within her first marriage (marriage dissolution or reaching the age of 48).

prev.ivl

The length of the previous time interval. May be used as explanatory variable in a Cox regression of birth intervals.

ses

Socio-economic status, a factor with levels lower, upper, farmer, and unknown.

parish

The Skelleftea region consists of three parishes, Jorn, Norsjo, and Skelleftea.

Details

The data set contain clusters of dependent observations defined by mother's id.

Source

Data is coming from The Demographic Data Base, Umea University, Umea, Sweden.

References

https://www.umu.se/enheten-for-demografi-och-aldrandeforskning/

Examples

data(fert)
fit <- coxreg(Surv(next.ivl, event) ~ ses + prev.ivl, data = fert, subset =
(parity == 1))
summary(fit)

Frailty experiment

Description

Utilizing GLMM models: Experimental, not exported (yet).

Usage

frail.fit(X, Y, rs, strats, offset, init, max.survs, frailty, control)

Arguments

Constant intensity discrete time proportional hazards

Description

This function is called from coxreg. A user may call it directly.

Usage

geome.fit(X, Y, rs, strats, offset, init, max.survs, method = "ml", control)

Arguments

Value

See the code.

Note

Nothing special

coxreg is a defunct function

Author(s)

Göran Broström

References

See coxreg.

See Also

coxreg

The Gompertz Distribution

Description

Density, distribution function, quantile function, hazard function, cumulative hazard function, and random generation for the Gompertz distribution with parameters shape and scale.

Usage

dgompertz(x, shape = 1, scale = 1, rate, log = FALSE, 
param = c("default", "canonical", "rate")) 
pgompertz(q, shape = 1, scale = 1, rate, lower.tail = TRUE, log.p = FALSE, 
param = c("default", "canonical", "rate")) 
qgompertz(p, shape = 1, scale = 1, rate, lower.tail = TRUE, log.p = FALSE, 
param = c("default", "canonical", "rate")) 
hgompertz(x, shape = 1, scale = 1, rate, log = FALSE, 
param = c("default", "canonical", "rate")) 
Hgompertz(x, shape = 1, scale = 1, rate, log.p = FALSE, 
param = c("default", "canonical", "rate")) 
rgompertz(n, shape = 1, scale = 1, rate, 
param = c("default", "canonical", "rate"))

Arguments

Details

The Gompertz distribution with scale parameter a and shape parameter \sigma has hazard function given by

h(x) = a \exp(x/\sigma)

for x \ge 0. If param = "canonical", then then a –> a/b, so that b is a true scale parameter (for any fixed a), and b is an 'AFT parameter'. If param = "rate", then b –> 1/b.

Value

dgompertz gives the density, pgompertz gives the distribution function, qgompertz gives the quantile function, hgompertz gives the hazard function, Hgompertz gives the cumulative hazard function, and rgompertz generates random deviates.

Invalid arguments will result in return value NaN, with a warning.

Get baseline hazards atoms from fits from

Description

Get baseline hazards atoms from fits from

Usage

hazards(x, cum = TRUE, ...)

Arguments

Value

A list where each component is a two-column matrix representing hazard atoms from one stratum. The first column is event time, and the second column is the corresponding hazard atom.

HISCO to HISCLASS transformation

Description

HISCO to HISCLASS transformation

Usage

HiscoHisclass(hisco, status = NULL, relation = NULL, urban = NULL,
debug = FALSE)

Arguments

Value

A vector of hisclass codes, same length as input hisco.

Author(s)

Göran Broström with translation and modification of a Stata do.

References

Van Leeuwen, M. and Maas, I. (2011). HISCLASS. A Historical International Social Class Scheme. Leuwen University Press.

strata function imported from survival

Description

This function is imported from the survival package. See strata.

Surv function imported from survival

Description

This function is imported from the survival package. See Surv.

Infant mortality and maternal death, Sweeden 1821–1894.

Description

Matched data on infant mortality, from seven parishes in Sweden, 1821–1894.

Usage

data(infants)

Format

A data frame with 80 rows and five variables.

stratum

Triplet No. Each triplet consist of one infant whose mother died (a case), and two controls, i.e, infants whose mother did not die. Matched on covariates below.

enter

Age (in days) of case when its mother died.

exit

Age (in days) at death or right censoring (at age 365 days).

event

Follow-up ends with death (1) or right censoring (0).

mother

dead for cases, alive for controls.

age

Mother's age at infant's birth.

sex

The infant's sex.

parish

Birth parish, either Nedertornea or not Nedertornea.

civst

Civil status of mother, married or unmarried.

ses

Socio-economic status of mother, either farmer or not farmer.

year

Year of birth of the infant.

Details

From 5641 first-born in seven Swedish parishes 1820-1895, from Fleninge in the very south to Nedertorneå in the very north, those whose mother died during their first year of life were selected, in all 35 infants. To each of them, two controls were selected by matching on the given covariates.

Source

Data originate from The Demographic Data Base, Umeå University, Umeå, Sweden, https://www.umu.se/enheten-for-demografi-och-aldrandeforskning/.

References

Broström, G. (1987). The influence of mother's death on infant mortality: A case study in matched data survival analysis. Scandinavian Journal of Statistics 14, 113-123.

Examples

data(infants)
fit <- coxreg(Surv(enter, exit, event) ~ strata(stratum) + mother, data
= infants)
fit
fit.w <- phreg(Surv(enter, exit, event) ~ mother + parish + ses, data =
infants)
summary(fit.w) ## Weibull proportional hazards model.

Straighten up a survival data frame

Description

Unnecessary cut spells are glued together, overlapping spells are "polished", etc.

Usage

join.spells(dat, strict = FALSE, eps = 1e-08)

Arguments

Details

In case of overlapping intervals (i.e., a data error), the appropriate id's are returned if strict is TRUE.

Value

A data frame with the same variables as the input, but individual spells are joined, if possible (identical covariate values, and adjacent time intervals).

Author(s)

Göran Broström

References

Therneau, T.M. and Grambsch, P.M. (2000). Modeling Survival Data: Extending the Cox model. Springer.

See Also

coxreg, aftreg, check.surv

The Loglogistic Distribution

Description

Density, distribution function, quantile function, hazard function, cumulative hazard function, and random generation for the Loglogistic distribution with parameters shape and scale.

Usage

dllogis(x, shape = 1, scale = 1, log = FALSE) 
pllogis(q, shape = 1, scale = 1, lower.tail = TRUE, log.p = FALSE) 
qllogis(p, shape = 1, scale = 1, lower.tail = TRUE, log.p = FALSE) 
hllogis(x, shape = 1, scale = 1, prop = 1, log = FALSE) 
Hllogis(x, shape = 1, scale = 1, prop = 1, log.p = FALSE) 
rllogis(n, shape = 1, scale = 1)

Arguments

Details

The Loglogistic distribution with scale parameter a and shape parameter \sigma has hazard function given by

h(x) = (b/\sigma)(x/\sigma)^(b-1)\exp((x / \sigma)^b)

for x \ge 0.

Value

dllogis gives the density, pllogis gives the distribution function, qllogis gives the quantile function, hllogis gives the hazard function, Hllogis gives the cumulative hazard function, and rllogis generates random deviates.

Invalid arguments will result in return value NaN, with a warning.

The Lognormal Distribution

Description

Density, distribution function, quantile function, hazard function, cumulative hazard function, and random generation for the Lognormal distribution with parameters shape and scale.

Usage

hlnorm(x, meanlog = 0, sdlog = 1, shape = 1 / sdlog, scale = exp(meanlog),
prop = 1, log = FALSE) 
Hlnorm(x, meanlog = 0, sdlog = 1, shape = 1 / sdlog, scale = exp(meanlog), 
prop = 1, log.p = FALSE)

Arguments

Details

The Lognormal distribution with scale parameter a and shape parameter \sigma has hazard function given by

h(x) = (b/\sigma)(x/\sigma)^(b-1)\exp((x / \sigma)^b)

for x \ge 0.

Value

dlnorm gives the density, plnorm gives the distribution function, qlnorm gives the quantile function, hlnorm gives the hazard function, Hlnorm gives the cumulative hazard function, and rlnorm generates random deviates.

Invalid arguments will result in return value NaN, with a warning.

The Log-rank test

Description

Performs the log-rank test on survival data, possibly stratified.

Usage

logrank(Y, group, data = parent.frame())

Arguments

Value

A list of class logrank with components

Note

The test is performed by fitting a Cox regression model and reporting its score test. With tied data, this might be slightly different from the true logrank test, but the difference is unimportant in practice.

Author(s)

Göran Broström

See Also

coxreg, print.logrank.

Examples

fit <- logrank(Y = Surv(enter, exit, event), group = civ, 
data = oldmort[oldmort$region == "town", ])
fit

Rye prices, Scania, southern Sweden, 1801-1894.

Description

The data consists of yearly rye prices from 1801 to 1894. Logged and detrended, so the time series is supposed to measure short term fluctuations in rye prices.

Usage

data(scania)

Format

A data frame with 94 observations in two columns on the following 2 variables.

year

The year the price is recorded.

foodprices

Detrended log rye prices.

Details

The Scanian area in southern Sweden was during the 19th century a mainly rural area.

Source

The Scanian Economic Demographic Database.

References

Jörberg, L. (1972). A History of Prices in Sweden 1732-1914, CWK Gleerup, Lund.

Examples

data(logrye)
summary(logrye)

LaTeX printing of regression results.

Description

This (generic) function prints the LaTeX code of the results of a fit from coxreg, phreg, tpchreg, or aftreg, similar to what xtable does for fits from other functions.

Usage

ltx(
  x,
  caption = NULL,
  label = NULL,
  dr = NULL,
  digits = max(options()$digits - 4, 3),
  ...
)

Arguments

Details

The result is a printout which is (much) nicer than the standard printed output from glm and friends,

Value

LaTeX code version of the results from a run with coxreg, phreg, phreg, or aftreg.

Note

For printing confidence limits, use ltx2.

Author(s)

Göran Broström.

See Also

ltx2, coxreg, phreg, phreg, and aftreg.

Examples

data(oldmort)
fit <- coxreg(Surv(enter, exit, event) ~ civ + sex, data = oldmort)
dr <- drop1(fit, test = "Chisq")
ltx(fit, dr = dr, caption = "A test example.", label = "tab:test1") 

LaTeX alternative printing of regression results.

Description

This (generic) function prints the LaTeX code of the results of a fit from coxreg, phreg, tpchreg, or aftreg.

Usage

ltx2(
  x,
  caption = NULL,
  label = NULL,
  dr = NULL,
  digits = max(options()$digits - 4, 4),
  conf = 0.95,
  keep = NULL,
  ...
)

Arguments

Value

LaTeX code version of the results from a run with coxreg, phreg, phreg, aftreg.

Note

Resulting tables contain estimated hazard ratios and confidence limits instead of regression coefficients and standard errors as in ltx.

Author(s)

Göran Broström.

See Also

xtable, coxreg, phreg, phreg, aftreg, and ltx.

Examples

data(oldmort)
fit <- coxreg(Surv(enter, exit, event) ~ sex, data = oldmort)
ltx2(fit, caption = "A test example.", label = "tab:test1") 

Put communals in "fixed" data frame

Description

Given an ordinary data frame suitable for survival analysis, and a data frame with "communal" time series, this function includes the communal covariates as fixed, by the "cutting spells" method.

Usage

make.communal(
  dat,
  com.dat,
  communal = TRUE,
  start,
  period = 1,
  lag = 0,
  surv = c("enter", "exit", "event", "birthdate"),
  tol = 1e-04,
  fortran = TRUE
)

Arguments

Details

The main purpose of this function is to prepare a data file for use with coxreg, aftreg, and coxph.

Value

The return value is a data frame with the same variables as in the combination of dat and com.dat. Therefore it is an error to have common name(s) in the two data frames.

Note

Not very vigorously tested.

Author(s)

Göran Broström

See Also

coxreg, aftreg, coxph, cal.window

Examples

dat <- data.frame(enter = 0, exit = 5.731, event = 1,
birthdate = 1962.505, x = 2)
## Birth date: July 2, 1962 (approximately).
com.dat <- data.frame(price = c(12, 3, -5, 6, -8, -9, 1, 7))
dat.com <- make.communal(dat, com.dat, start = 1962.000) 

The Gompertz-Makeham Distribution

Description

Density, distribution function, quantile function, hazard function, cumulative hazard function, and random generation for the Gompertz-Makeham distribution with parameters shape and scale.

Usage

dmakeham(x, shape = c(1, 1), scale = 1, log = FALSE)
pmakeham(q, shape = c(1, 1), scale = 1, lower.tail = TRUE, log.p = FALSE)
qmakeham(p, shape = c(1, 1), scale = 1, lower.tail = TRUE, log.p = FALSE)
hmakeham(x, shape = c(1, 1), scale = 1, log = FALSE)
Hmakeham(x, shape = c(1, 1), scale = 1, log.p = FALSE)
rmakeham(n, shape = c(1, 1), scale = 1)

Arguments

Details

The Gompertz-Makeham distribution with shape parameter a and scale parameter \sigma has hazard function given by

h(x) = a[2] + a[1] \exp(x/\sigma)

for x \ge 0.

Value

dmakeham gives the density, pmakeham gives the distribution function, qmakeham gives the quantile function, hmakeham gives the hazard function, Hmakeham gives the cumulative hazard function, and rmakeham generates random deviates.

Invalid arguments will result in return value NaN, with a warning.

Male mortality in ages 40-60, nineteenth century

Description

Males born in the years 1800-1820 and surving at least 40 years in the parish Skellefteå in northern Sweden are followed from their fortieth birthday until death or the sixtieth birthday, whichever comes first.

Usage

data(male.mortality)

Format

A data frame with 2058 observations on the following 6 variables.

id

Personal identification number.

enter

Start of duration. Measured in years since the fortieth birthday.

exit

End of duration. Measured in years since the fortieth birthday.

event

a logical vector indicating death at end of interval.

birthdate

The birthdate in decimal form.

ses

Socio-economic status, a factor with levels lower, upper

Details

The interesting explanatory covariate is ses (socioeconomic status), which is a time-varying covariate. This explains why several individuals are representated by more than one record each. Left trucation and right censoring are introduced this way.

Note

This data set is also known, and accessible, as mort.

Source

Data is coming from The Demographic Data Base, Umea University, Umeå, Sweden.

References

https://www.umu.se/enheten-for-demografi-och-aldrandeforskning/

Examples

data(male.mortality)
fit <- coxreg(Surv(enter, exit, event) ~ ses, data = male.mortality)
summary(fit)

ML proportional hazards regression

Description

Maximum Likelihood estimation of proportional hazards models. Is deprecated, use coxreg instead.

Usage

mlreg(
  formula = formula(data),
  data = parent.frame(),
  na.action = getOption("na.action"),
  init = NULL,
  method = c("ML", "MPPL"),
  control = list(eps = 1e-08, maxiter = 10, n.points = 12, trace = FALSE),
  singular.ok = TRUE,
  model = FALSE,
  center = TRUE,
  x = FALSE,
  y = TRUE,
  boot = FALSE,
  geometric = FALSE,
  rs = NULL,
  frailty = NULL,
  max.survs = NULL
)

Arguments

Details

Method ML performs a true discrete analysis, i.e., one parameter per observed event time. Method MPPL is a compromize between the discrete and continuous time approaches; one parameter per observed event time with multiple events. With no ties in data, an ordinary Cox regression (as with coxreg) is performed.

Value

A list of class c("mlreg", "coxreg", "coxph") with components

Warning

The use of rs is dangerous, see note above. It can however speed up computing time.

Note

This function starts by creating risksets, if no riskset is supplied via rs, with the aid of risksets. This latter mechanism fails if there are any NA's in the data! Note also that it depends on stratification, so rs contains information about stratification. Giving another strata variable in the formula is an error. The same is ok, for instance to supply stratum interactions.

Note futher that mlreg is deprecated. coxreg should be used instead.

Author(s)

Göran Broström

References

Broström, G. (2002). Cox regression; Ties without tears. Communications in Statistics: Theory and Methods 31, 285–297.

See Also

coxreg, risksets

Examples

 dat <- data.frame(time=  c(4, 3,1,1,2,2,3),
                status=c(1,1,1,0,1,1,0),
                x=     c(0, 2,1,1,1,0,0),
                sex=   c(0, 0,0,0,1,1,1))
 mlreg( Surv(time, status) ~ x + strata(sex), data = dat) #stratified model
 # Same as:
 rs <- risksets(Surv(dat$time, dat$status), strata = dat$sex)
 mlreg( Surv(time, status) ~ x, data = dat, rs = rs) #stratified model
 

Male mortality in ages 40-60, nineteenth century

Description

Males born in the years 1800-1820 and surving at least 40 years in the parish Skellefteå in northern Sweden are followed from their fortieth birthday until death or the sixtieth birthday, whichever comes first.

Usage

data(mort)

Format

A data frame with 2058 observations on the following 6 variables.

id

Personal identification number.

enter

Start of duration. Measured in years since the fortieth birthday.

exit

End of duration. Measured in years since the fortieth birthday.

event

a logical vector indicating death at end of interval.

birthdate

The birthdate in decimal form.

ses

Socio-economic status, a factor with levels lower, upper

Details

The interesting explanatory covariate is ses (socioeconomic status), which is a time-varying covariate. This explains why several individuals are representated by more than one record each. Left trucation and right censoring are introduced this way.

Note

This data set is also known, and accessible, as male.mortality

Source

Data is coming from The Demographic Data Base, Umea University, Umeå, Sweden.

References

https://www.umu.se/enheten-for-demografi-och-aldrandeforskning/

Examples

data(mort)
fit <- coxreg(Surv(enter, exit, event) ~ ses, data = mort)
summary(fit)

Create an oe object

Description

Create an oe ("occurrence/exposure") object, used as a response variable in a model formula specifically in tpchreg.

Usage

oe(count, exposure)

Arguments

See Also

tpchreg.

Old age mortality, Sundsvall, Sweden, 1860-1880.

Description

The data consists of old age life histories from 1 January 1860 to 31 december 1880, 21 years. Only (parts of) life histories above age 60 is considered.

Usage

data(oldmort)

Format

A data frame with 6508 observations from 4603 persons on the following 13 variables.

id

Identification number.

enter

Start age for the interval.

exit

Stop age for the interval.

event

Indicator of death; equals TRUE if the person died at the end of the interval, FALSE otherwise.

birthdate

Birthdate as a real number (i.e., "1765-06-27" is 1765.490).

m.id

Mother's identification number.

f.id

Father's identification number.

sex

Gender, a factor with levels male female

civ

Civil status, a factor with levels unmarried married widow

ses.50

Socio-economic status at age 50, a factor with levels middle unknown upper farmer lower

birthplace

a factor with levels parish region remote

imr.birth

Infant mortality rate at birth in the region of birth

region

Subregion of Sundsvall, a factor with levels town industry rural

Details

The Sundsvall area in mid-Sweden was during the 19th century a fast growing forest industry. At the end of the century, it was one of the largest sawmill area in Europe. The town Sundsvall is fast growing part of the region and center for the commerse.

Source

The Demographic Data Base, Umeå University, Sweden.

References

Edvinsson, S. (2000). The Demographic Data Base at Umeå University: A resource for historical studies. In Hall, McKaa, and Thorvaldsen (eds), "Handbook of International Historical Microdata for Population Research", Minnesota Population Center, Minneapolis.

Examples

data(oldmort)
summary(oldmort)
## maybe str(oldmort) ; plot(oldmort) ...

The Piecewise Constant Hazards distribution.

Description

Density, distribution function, quantile function, hazard function, cumulative hazard function, mean, and random generation for the Piecewice Constant Hazards (pch) distribution.

Usage

ppch(q, cuts, levels, lower.tail = TRUE, log.p = FALSE)
dpch(x, cuts, levels, log = FALSE)
hpch(x, cuts, levels, log = FALSE)
Hpch(x, cuts, levels, log.p = FALSE)
qpch(p, cuts, levels, lower.tail = TRUE, log.p = FALSE)
mpch(cuts, levels)
rpch(n, cuts, levels)

Arguments

Details

The pch distribution has a hazard function that is piecewise constant on intervals defined by cutpoints

0 < c_1 < \cdots < c_n < \infty, n \ge 0

If n = 0, this reduces to an exponential distribution.

Value

dpch gives the density, ppch gives the distribution function, qpch gives the quantile function, hpch gives the hazard function, Hpch gives the cumulative hazard function, mpch gives the mean, and rpch generates random deviates.

Note

the parameter levels must have length at least 1, and the number of cut points must be one less than the number of levels.

Piecewise Constant Proportional Hazards Regression

Description

Proportional hazards model with piecewise constant baseline hazard(s). Allows for stratification and left truncated and right censored data.

Usage

pchreg(
  formula = formula(data),
  data = parent.frame(),
  na.action = getOption("na.action"),
  cuts = NULL,
  init,
  control = list(eps = 1e-08, maxiter = 20, trace = FALSE),
  singular.ok = TRUE,
  model = FALSE,
  x = FALSE,
  y = TRUE
)

Arguments

Value

A list of class "pchreg" with components

Author(s)

Göran Broström

See Also

phreg, coxreg, link{aftreg}.

Examples

## Not run: 
dat <- age.window(oldmort, c(60, 80))
fit <- pchreg(Surv(enter, exit, event) ~ ses.50 + sex, 
data = dat, cuts = seq(60, 80, by = 4))
summary(fit)

fit.cr <- coxreg(Surv(enter, exit, event) ~ ses.50 + sex, data = dat)
check.dist(fit.cr, fit, main = "Cumulative hazards")

## End(Not run)

Period statistics

Description

Calculates occurrence / exposure rates for time periods given by period and for ages given by age.

Usage

perstat(surv, period, age = c(0, 200))

Arguments

Value

A list with components

Author(s)

Göran Broström

See Also

piecewise

Loglihood function of a proportional hazards regression

Description

Calculates minus the log likelihood function and its first and second order derivatives for data from a Weibull regression model.

Usage

phfunc(
  beta = NULL,
  lambda,
  p,
  X = NULL,
  Y,
  offset = rep(0, length(Y)),
  ord = 2,
  pfixed = FALSE,
  dist = "weibull"
)

Arguments

Details

Note that the function returns log likelihood, score vector and minus hessian, i.e. the observed information. The model is

S(t; p, \lambda, \beta, z) = S_0((t / \lambda)^p)^{e^(z \beta)}

Value

A list with components

Author(s)

Göran Broström

See Also

phreg

Parametric Proportional Hazards Regression

Description

Proportional hazards model with parametric baseline hazard(s). Allows for stratification with different scale and shape in each stratum, and left truncated and right censored data.

Usage

phreg(
  formula = formula(data),
  data = parent.frame(),
  na.action = getOption("na.action"),
  dist = "weibull",
  cuts = NULL,
  init,
  shape = 0,
  param = c("canonical", "rate"),
  control = list(eps = 1e-08, maxiter = 20, trace = FALSE),
  singular.ok = TRUE,
  model = FALSE,
  x = FALSE,
  y = TRUE
)

Arguments

Details

The parameterization is the same as in coxreg and coxph, but different from the one used by survreg (which is not a proportional hazards modelling function). The model is

S(t; a, b, \beta, z) = S_0((t/b)^a)^{\exp((z-mean(z))\beta)}

where S0 is some standardized survivor function.

Value

A list of class c("phreg", "coxreg") with components

Warning

The lognormal and loglogistic distributions are included on an experimental basis for the moment. Use with care, results may be unreliable!

The gompertz distribution has an exponentially increasing hazard function under the canonical parametrization. This may cause instability in the convergence of the fitting algorithm in the case of near-exponential data. It may be resolved by using param = "rate".

Note

The lognormal and loglogistic baseline distributions are extended to a three-parameter family by adding a "proportionality" parameter (multiplying the baseline hazard function). The log of the estimated parameter turns up as '(Intercept)' in the printed output. The reason for this extension is that the standard lognormal and loglogistic distributions are not closed under proportional hazards.

Author(s)

Göran Broström

See Also

coxreg, check.dist, link{aftreg}.

Examples

data(mort)
fit <- phreg(Surv(enter, exit, event) ~ ses, data = mort)
fit
plot(fit)
fit.cr <- coxreg(Surv(enter, exit, event) ~ ses, data = mort)
check.dist(fit.cr, fit)

Parametric proportional hazards regression

Description

This function is called by phreg, but it can also be directly called by a user.

Usage

phreg.fit(X, Y, dist, strata, offset, init, shape, control)

Arguments

Details

See phreg for more detail.

Value

Author(s)

Göran Broström

See Also

phreg

Piecewise hazards

Description

Calculate piecewise hazards, no. of events, and exposure times in each interval indicated by cutpoints.

Usage

piecewise(enter, exit, event, cutpoints)

Arguments

Details

Exact calculation.

Value

A list with components

Author(s)

Göran Broström

See Also

perstat

Plots output from an AFT regression

Description

Just a simple plot of the hazard (cumulative hazard, density, survival) functions for each stratum.

Usage

## S3 method for class 'aftreg'
plot(
  x,
  fn = c("haz", "cum", "den", "sur"),
  main = NULL,
  xlim = NULL,
  ylim = NULL,
  xlab = "Duration",
  ylab = "",
  col,
  lty,
  printLegend = TRUE,
  new.data = x$means,
  ...
)

Arguments

Details

The plot is drawn at the mean values of the covariates, by default.

Value

No return value.

Author(s)

Göran Broström

See Also

aftreg

Examples

y <- rllogis(40, shape = 1, scale = 1)
x <- rep(c(1,1,2,2), 10)
fit <- aftreg(Surv(y, rep(1, 40)) ~ x, dist = "loglogistic")
plot(fit)

Plot method for coxreg objects

Description

A plot of a baseline function of a coxreg fit is produced, one curve for each stratum. A wrapper for plot.survfit.

Usage

## S3 method for class 'coxreg'
plot(
  x,
  fn = c("cum", "surv", "log", "loglog"),
  conf.int = FALSE,
  fig = TRUE,
  xlim = NULL,
  ylim = NULL,
  main = NULL,
  xlab = "Duration",
  ylab = "",
  col = 1,
  lty = 1,
  printLegend = TRUE,
  ...
)

Arguments

Value

An object of class hazdata containing the coordinates of the curve(s).

Plots of hazdata objects.

Description

Baseline hazards estimates.

Usage

## S3 method for class 'hazdata'
plot(
  x,
  strata = NULL,
  fn = c("cum", "surv", "log", "loglog"),
  fig = TRUE,
  xlim = NULL,
  ylim = NULL,
  main = NULL,
  xlab = "",
  ylab = "",
  col = "black",
  lty = 1,
  printLegend = TRUE,
  ...
)

Arguments

Details

It is also possible to have as first argument an object of type "coxreg", given that it contains a component of type "hazdata".

Value

A list where the elements are two-column matrices, one for each stratum in the model. The first column contains risktimes, and the second the y coordinates for the requested curve(s).

Note

x is a list where each element is a two-column matrix. The first column contains failure times, and the second column contains the corresponding 'hazard atoms'.

Author(s)

Göran Broström

Examples

time0 <- numeric(50)
group <- c(rep(0, 25), rep(1, 25))
x <- runif(50, -0.5, 0.5)
time1 <- rexp( 50, exp(group) )
event <- rep(1, 50)
fit <- coxreg(Surv(time0, time1, event) ~ x + strata(group), method = "ml")
plot(fit$hazards, col = 1:2, fn = "surv", xlab = "Duration")
## Same result as:
## plot(fit, col = 1:2, fn = "sur", xlab = "Duration")

Plots of hazdata objects.

Description

Baseline hazards estimates.

Usage

## S3 method for class 'logrank'
plot(
  x,
  fn = c("cum", "surv", "log", "loglog"),
  xlim = NULL,
  ylim = NULL,
  main = NULL,
  xlab = "",
  ylab = "",
  col = "black",
  lty = 1,
  printLegend = TRUE,
  ...
)

Arguments

Details

It is also possible to have as first argument an object of type "coxreg", given that it contains a component of type "hazdata".

Value

A list where the elements are two-column matrices, one for each stratum in the model. The first column contains risktimes, and the second the y coordinates for the requested curve(s).

Note

x is a list where each element is a two-column matrix. The first column contains failure times, and the second column contains the corresponding 'hazard atoms'.

Author(s)

Göran Broström

Examples

fit <- logrank(Surv(enter, exit, event), group = civ, data = oldmort[oldmort$region == "town", ])
plot(fit)

Plots output from a phreg regression

Description

Plot(s) of the hazard, density, cumulative hazards, and/or the survivor function(s) for each stratum.

Usage

## S3 method for class 'phreg'
plot(
  x,
  fn = c("haz", "cum", "den", "sur"),
  main = NULL,
  xlim = NULL,
  ylim = NULL,
  xlab = "Duration",
  ylab = "",
  col,
  lty,
  printLegend = TRUE,
  score = 1,
  fig = TRUE,
  ...
)

Arguments

Value

No return value if fig = TRUE, otherwise the cumulative hazards function (coordinates), given fn = "cum".

Note

Reference hazard is given by the fit; zero for all covariates, and the reference category for factors.

Author(s)

Göran Broström

See Also

phreg

Examples

y <- rllogis(40, shape = 1, scale = 1)
x <- rep(c(1,1,2,2), 10)
fit <- phreg(Surv(y, rep(1, 40)) ~ x, dist = "loglogistic")
plot(fit)

Plots output from a tpchreg regression

Description

Plot(s) of the hazard, cumulative hazards, and/or the survivor function(s) for each stratum.

Usage

## S3 method for class 'tpchreg'
plot(
  x,
  fn = c("haz", "cum", "sur"),
  log = "",
  main = NULL,
  xlim = NULL,
  ylim = NULL,
  xlab = "Duration",
  ylab = "",
  col,
  lty,
  printLegend = TRUE,
  ...
)

Arguments

Value

No return value.

Author(s)

Göran Broström

See Also

tpchreg

Plots output from a Weibull regression

Description

Plot(s) of the hazard, density, cumulative hazards, and/or the survivor function(s) for each stratum.

Usage

## S3 method for class 'weibreg'
plot(
  x,
  fn = c("haz", "cum", "den", "sur"),
  main = NULL,
  xlim = NULL,
  ylim = NULL,
  xlab = NULL,
  ylab = NULL,
  new.data = x$means,
  ...
)

Arguments

Details

The plot is drawn at the mean values of the covariates.

Value

No return value

Author(s)

Göran Broström

See Also

phreg, weibreg

Examples

y <- rweibull(4, shape = 1, scale = 1)
x <- c(1,1,2,2)
fit <- weibreg(Surv(y, c(1,1,1,1)) ~ x)
plot(fit)

Graphical comparing of cumulative hazards

Description

Comparison of the cumulative hazards functions for a semi-parametric and parametric models.

Usage

plotHaz(
  sp,
  pp,
  interval,
  main = NULL,
  xlab = "Time",
  ylab = "Cum. hazards",
  leglab,
  col = c("blue", "red"),
  lty = 1:2,
  ylim,
  log = "",
  printLegend = TRUE
)

Arguments

Details

For the moment only a graphical comparison. The arguments sp and pp may be swapped.

Value

No return value.

Author(s)

Göran Broström

See Also

check.dist, coxreg and phreg.

Examples

data(mort)
op <- par(mfrow = c(1, 2))
fit.cr <- coxreg(Surv(enter, exit, event) ~ ses, data = mort)
fit.w <- phreg(Surv(enter, exit, event) ~ ses, data = mort)
fit.g <- phreg(Surv(enter, exit, event) ~ ses, data = mort,
dist = "gompertz")
plotHaz(fit.cr, fit.w, interval = c(0, 20), main = "Weibull")
plotHaz(fit.cr, fit.g, main = "Gompertz")
par(op)

Prints aftreg objects

Description

The hazard, the cumulative hazard, the density, and the survivor baseline functions are plotted.

Usage

## S3 method for class 'aftreg'
print(x, digits = max(options()$digits - 4, 3), ...)

Arguments

Value

No value is returned.

Note

Doesn't work for threeway or higher order interactions. Use print.coxph in that case.

Author(s)

Göran Broström

See Also

phreg, print.coxph

Prints coxreg objects

Description

More "pretty-printing" than print.coxph, which is a fall-back for 'difficult' objects.

Usage

## S3 method for class 'coxreg'
print(x, digits = max(options()$digits - 4, 3), ...)

Arguments

Details

Doesn't work with three-way and higher interactions, in which case print.coxph is used.

Value

No value is returned.

Author(s)

Göran Broström

See Also

coxreg, print.coxph

Prints logrank objects

Description

The result of logrank is printed

Usage

## S3 method for class 'logrank'
print(x, digits = max(options()$digits - 4, 6), ...)

Arguments

Value

The input is returned invisibly.

Author(s)

Göran Broström

See Also

logrank, coxreg

Prints phreg objects

Description

The hazard, the cumulative hazard, the density, and the survivor baseline functions are plotted.

Usage

## S3 method for class 'phreg'
print(x, digits = max(options()$digits - 4, 3), ...)

Arguments

Value

No value is returned.

Note

Doesn't work for threeway or higher order interactions. Use print.coxph in that case.

Author(s)

Göran Broström

See Also

phreg, print.coxph

Prints a summary of the content of a set of risk sets.

Description

Given the output from risksets, summary statistics are given for it.

Usage

## S3 method for class 'risksets'
print(x, ...)

Arguments

Value

No value is returned; the function prints summary statistics of risk sets.

Note

There is no summary.risksets yet. On the TODO list.

Author(s)

Göran Broström

See Also

risksets

Examples

rs <- with(mort, risksets(Surv(enter, exit, event)))
print(rs)

Prints summary.aftreg objects

Description

Prints summary.aftreg objects

Usage

## S3 method for class 'summary.aftreg'
print(x, digits = max(getOption("digits") - 3, 3), short = FALSE, ...)

Arguments

Value

No value is returned.

Author(s)

Göran Broström

See Also

aftreg, summary.aftreg

Prints summary.coxreg objects

Description

Prints summary.coxreg objects

Usage

## S3 method for class 'summary.coxreg'
print(x, digits = 3, short = FALSE, ...)

Arguments

Value

No value is returned.

Author(s)

Göran Broström

See Also

coxreg, summary.coxreg

Prints summary.phreg objects

Description

Prints summary.phreg objects

Usage

## S3 method for class 'summary.phreg'
print(x, digits = max(getOption("digits") - 3, 3), ...)

Arguments

Value

No value is returned.

Author(s)

Göran Broström

See Also

phreg, summary.phreg

Prints summary.tpchreg objects

Description

Prints summary.tpchreg objects

Usage

## S3 method for class 'summary.tpchreg'
print(x, digits = max(getOption("digits") - 3, 3), ...)

Arguments

Value

No value is returned.

Author(s)

Göran Broström

See Also

tpchreg, summary.tpchreg

Prints tpchreg objects

Description

More "pretty-printing"

Usage

## S3 method for class 'tpchreg'
print(x, digits = max(options()$digits - 4, 3), ...)

Arguments

Details

Doesn't work with three-way or higher interactions.

Value

No value is returned.

Author(s)

Göran Broström

See Also

tpchreg, print.coxreg

Prints weibreg objects

Description

The hazard, the cumulative hazard, the density, and the survivor baseline functions are plotted.

Usage

## S3 method for class 'weibreg'
print(x, digits = max(options()$digits - 4, 3), ...)

Arguments

Value

No value is returned.

Note

Doesn't work for threeway or higher order interactions. Use print.coxph in that case.

Author(s)

Göran Broström

See Also

weibreg, print.coxph

Retrieves regression tables

Description

Retrieves regression tables

Usage

regtable(x, digits = 3, short = TRUE, ...)

Arguments

Value

A character data frame, ready to print in various formats.

Note

Should not be used if interactions present.

Author(s)

Göran Broström

See Also

coxreg, summary.coxreg

Finds the compositions and sizes of risk sets

Description

Focus is on the risk set composition just prior to a failure.

Usage

risksets(
  x,
  strata = NULL,
  max.survs = NULL,
  members = TRUE,
  collate_sets = FALSE
)

Arguments

Details

If the input argument max.survs is left alone, all survivors are accounted for in all risk sets.

Value

A list with components (if collate_sets = FALSE)

Note

Can be used to "sample the risk sets".

Author(s)

Göran Broström

See Also

table.events, coxreg.

Examples

 enter <- c(0, 1, 0, 0)
 exit <- c(1, 2, 3, 4)
 event <- c(1, 1, 1, 0)
 risksets(Surv(enter, exit, event))

Old age mortality, Scania, southern Sweden, 1813-1894.

Description

The data consists of old age life histories from 1 January 1813 to 31 december 1894. Only (parts of) life histories above age 50 is considered.

Usage

data(scania)

Format

A data frame with 1931 observations from 1931 persons on the following 9 variables.

id

Identification number (enumeration).

enter

Start age for the interval.

exit

Stop age for the interval.

event

Indicator of death; equals TRUE if the person died at the end of the interval, FALSE otherwise.

birthdate

Birthdate as a real number (i.e., "1765-06-27" is 1765.490).

sex

Gender, a factor with levels male female.

parish

One of five parishes in Scania, coded 1, 2, 3, 4, 5. Factor.

ses

Socio-economic status at age 50, a factor with levels upper and lower.

immigrant

a factor with levels no region and yes.

Details

The Scanian area in southern Sweden was during the 19th century a mainly rural area.

Source

The Scanian Economic Demographic Database, Lund University, Sweden.

References

https://www.lusem.lu.se/organisation/research-centres/centre-economic-demography/cedpop-databases-ced

Examples

data(scania)
summary(scania)

Prints aftreg objects

Description

Prints aftreg objects

Usage

## S3 method for class 'aftreg'
summary(object, ...)

Arguments

Author(s)

Göran Broström

See Also

print.coxreg

Examples

## The function is currently defined as
function (object, ...) 
print(object)

A summary of coxreg objects.

Description

A summary of coxreg objects.

Usage

## S3 method for class 'coxreg'
summary(object, ...)

Arguments

Author(s)

Göran Broström

See Also

print.coxreg

Examples

fit <- coxreg(Surv(enter, exit, event) ~ sex + civ, data = oldmort)
summary(fit)
 

A summary of phreg objects.

Description

A summary of phreg objects.

Usage

## S3 method for class 'phreg'
summary(object, ...)

Arguments

Author(s)

Göran Broström

See Also

print.phreg

Examples

fit <- phreg(Surv(enter, exit, event) ~ sex + civ, 
data = oldmort[oldmort$region == "town", ])
summary(fit)
 

Summary for tpchreg objects

Description

Summary for tpchreg objects

Usage

## S3 method for class 'tpchreg'
summary(object, ci = FALSE, level = 0.95, ...)

Arguments

Author(s)

Göran Broström

See Also

tpchreg

Examples

## The function is currently defined as
## function (object, ...) 

Prints a weibreg object

Description

This is the same as print.weibreg

Usage

## S3 method for class 'weibreg'
summary(object, ...)

Arguments

Author(s)

Göran Broström

See Also

print.weibreg

Examples

## The function is currently defined as
function (object, ...) 
print(object)

Split a survival object at specified durations.

Description

Given a survival object, (a matrix with two or three columns) and a set of specified cut times, split each record into multiple subrecords at each cut time. The new survival object will be in ‘counting process’ format, with an enter time, exit time, and event status for each record.

Usage

SurvSplit(Y, cuts)

Arguments

Value

A list with components

Note

This function is used in phreg for the piecewise constant hazards model. It uses age.window for each interval.

Author(s)

Göran Broström

See Also

survSplit, age.window.

Examples

##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.

## The function is currently defined as
function(Y, cuts){
    if (NCOL(Y) == 2) Y <- cbind(rep(0, NROW(Y)), Y)
    indat <- cbind(Y, 1:NROW(Y), rep(-1, NROW(Y)))
    colnames(indat) <- c("enter", "exit", "event", "idx", "ivl")
    n <- length(cuts)
    cuts <- sort(cuts)
    if ((cuts[1] <= 0) || (cuts[n] == Inf))
        stop("'cuts' must be positive and finite.")
    cuts <- c(0, cuts, Inf)
    n <- n + 1
    out <- list()
    indat <- as.data.frame(indat)
    for (i in 1:n){
        out[[i]] <- age.window(indat, cuts[i:(i+1)])
        out[[i]]$ivl <- i
        out[[i]] <- t(out[[i]])
    }
    Y <- matrix(unlist(out), ncol = 5, byrow = TRUE)
    colnames(Y) <- colnames(indat)
    list(Y = Y[, 1:3],
         ivl = Y[, 5],
         idx = Y[, 4]
         )
  }

Swedish death data, 1969-2020.

Description

A data frame containing data on the number of deaths by sex, age and year, Sweden 1969-2020.

Usage

data(swedeaths)

Format

A data frame with 5 variables and 10504 observations.

age

Numerical with integer values 0-100, representing achieved age in years during the actual calendar year. The highest value, 100, represents ages 100 and above.

sex

A factor with two levels, "women" and "men".

year

Calendar year.

deaths

Number of deaths by age, sex, and year.

id

Created by the reshape procedure, see Details.

Details

Data are downloaded from Statistics Sweden in the form of a csv file and and in that process converted to a data frame. Variable names are translated from Swedish, and some of them are coverted to factors. Each numeric column contains the number of deaths by sex and age. The original data set is in wide form and then converted to long format.

Source

Statistics Sweden, https://www.scb.se.

See Also

swepop, tpchreg

Examples

summary(swedeaths)
## maybe str(swedeaths) ...

Swedish population data, 1969-2020.

Description

A data frame containing data on the population size by sex, age and year, Sweden 1969-2020.

Usage

data(swepop)

Format

A data frame with 5 variables and 10504 observations.

age

Numerical with integer values 0-100, representing achieved age in years during the actual calendar year.The highest value, 100, represents ages 100 and above.

sex

A factor with two levels, "women" and "men".

year

Calendar year.

pop

Average population by age, sex, and year.

id

Created by the reshape procedure, see Details.

Details

Data are downloaded from Statistics Sweden in the form of a csv file and converted to a data frame. Variable names are translated from Swedish, and some of them are coverted to factors. The variable pop contains the average population by sex and age, calculated by taking the mean value of the population size at December 31 the previous year and December 31 the current year. The original data contain the sizes at the end of each year. The original data set is in wide format and converted to long format by reshape.

Source

Statistics Sweden, https://www.scb.se.

See Also

swedeaths

Examples

summary(swepop)
## maybe str(swepop) ...

Calculating failure times, risk set sizes and No. of events in each risk set

Description

From input data of the 'interval' type, with an event indicator, summary statistics for each risk set (at an event time point) are calculated.

Usage

table.events(enter = rep(0, length(exit)), exit, event, strict = TRUE)

Arguments

Value

A list with components

Author(s)

Göran Broström

See Also

risksets

Examples

exit = c(1,2,3,4,5)
event = c(1,1,0,1,1)
table.events(exit = exit, event = event)

Transforms a "survival" data frame into a data frame suitable for binary (logistic) regression

Description

The result of the transformation can be used to do survival analysis via logistic regression. If the cloglog link is used, this corresponds to a discrete time analogue to Cox's proportional hazards model.

Usage

toBinary(
  dat,
  surv = c("enter", "exit", "event"),
  strats,
  max.survs = NROW(dat)
)

Arguments

Details

toBinary calls risksets in the eha package.

Value

Returns a data frame expanded risk set by risk set. The three "survival variables" are replaced by a variable named event (which overwrites an eventual variable by that name in the input). Two more variables are created, riskset and orig.row.

Note

The survival variables must be three. If you only have exit and event, create a third containing all zeros.

Author(s)

Göran Broström

See Also

coxreg, glm.

Examples

enter <- rep(0, 4)
exit <- 1:4
event <- rep(1, 4)
z <- rep(c(-1, 1), 2)
dat <- data.frame(enter, exit, event, z)
binDat <- toBinary(dat)
dat
binDat
coxreg(Surv(enter, exit, event) ~ z, method = "ml", data = dat)
## Same as:
summary(glm(event ~ z + riskset, data = binDat, family = binomial(link = cloglog)))

Convert time in years since "0000-01-01" to a date.

Description

This function uses as.Date and a simple linear transformation.

Usage

toDate(times)

Arguments

Value

A vector of dates as character strings of the type "1897-05-21".

Author(s)

Göran Broström

See Also

toTime

Examples

##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.
toDate(1897.357)

Calculate duration in years from "0000-01-01" to a given date

Description

Given a vector of dates, the output is a vector of durations in years since "0000-01-01".

Usage

toTime(dates)

Arguments

Value

A vector of durations, as decribed above.

Author(s)

Göran Broström

See Also

toDate

Examples

##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.


toTime(c("1897-05-16", "1901-11-21"))

Transform survival data to tabular form

Description

Transform a "survival data frame" to tabular form aggregating number of events and exposure time by time intervals and covariates.

Usage

toTpch(formula, data, cuts, enter = "enter", exit = "exit",
event = "event", episode = "age")

Arguments

Details

If cuts is missing, nothing is done. Internally, this function first calls survival::survSplit and then stats::aggregate.

Value

A data frame with exposure time and number of events aggregated by time intervals and covariates. If all covariates are factors,this usually results in a huge reduction of the size of thedata frame, but otherwise the size of the output may be larger than the size of the input data frame

Note

Episodes, or parts of episodes, outside min(cuts), max(cuts) are cut off. With continuous covariates, consider rounding them so that the number of distinct oberved values is not too large.

Author(s)

Göran Broström

Proportional hazards regression with piecewise constant hazards and tabular data.

Description

Proportional hazards regression with piecewise constant hazards and tabular data.

Usage

tpchreg(formula, data, time, weights, last, subset, na.action, 
contrasts = NULL, start.coef = NULL, 
control = list(epsilon = 1.e-8, maxit = 200, trace = FALSE))

Arguments

Note

The interpretation of cuts is different from that in hpch. This is intentional.

See Also

oe.

Examples

sw <- swepop
sw$deaths <- swedeaths$deaths
fit <- tpchreg(oe(deaths, pop) ~ strata(sex) + I(year - 2000), 
time = age, last = 101, data = sw[sw$year >= 2000, ])
summary(fit)

Weibull Regression

Description

Proportional hazards model with baseline hazard(s) from the Weibull family of distributions. Allows for stratification with different scale and shape in each stratum, and left truncated and right censored data.

Usage

weibreg(
  formula = formula(data),
  data = parent.frame(),
  na.action = getOption("na.action"),
  init,
  shape = 0,
  control = list(eps = 1e-04, maxiter = 10, trace = FALSE),
  singular.ok = TRUE,
  model = FALSE,
  x = FALSE,
  y = TRUE,
  center = TRUE
)

Arguments

Details

The parameterization is the same as in coxreg and coxph, but different from the one used by survreg. The model is

h(t; a, b, \beta, z) = (a/b) (t/b)^{a-1} exp(z\beta)

This is in correspondence with Weibull. To compare regression coefficients with those from survreg you need to divide by estimated shape (\hat{a}) and change sign. The p-values and test statistics are however the same, with one exception; the score test is done at maximized scale and shape in weibreg.

This model is a Weibull distribution with shape parameter a and scale parameter b \exp(-z\beta / a)

Value

A list of class c("weibreg", "coxreg") with components

Warning

The print method print.weibreg doesn't work if threeway or higher order interactions are present.

Note further that covariates are internally centered, if center = TRUE, by this function, and this is not corrected for in the output. This affects the estimate of \log(scale), but nothing else. If you don't like this, set center = FALSE.

Note

This function is not maintained, and may behave in unpredictable ways. Use phreg with dist = "weibull" (the default) instead! Will soon be declared deprecated.

Author(s)

Göran Broström

See Also

phreg, coxreg, print.weibreg

Examples

 dat <- data.frame(time = c(4, 3, 1, 1, 2, 2, 3),
                status = c(1, 1, 1, 0, 1, 1, 0),
                x = c(0, 2, 1, 1, 1, 0, 0),
                sex = c(0, 0, 0, 0, 1, 1, 1))
 weibreg( Surv(time, status) ~ x + strata(sex), data = dat) #stratified model

Weibull regression

Description

This function is called by weibreg, but it can also be directly called by a user.

Usage

weibreg.fit(X, Y, strata, offset, init, shape, control, center = TRUE)

Arguments

Details

See weibreg for more detail.

Value

Author(s)

Göran Broström

See Also

weibreg

The (Cumulative) Hazard Function of a Weibull Distribution

Description

hweibull calculates the hazard function of a Weibull distribution, and Hweibull calculates the corresponding cumulative hazard function.

Usage

hweibull(x, shape, scale = 1, log = FALSE)

Arguments

Details

See dweibull.

Value

The (cumulative) hazard function, evaluated at x.

Author(s)

Göran Broström

See Also

pweibull

Examples

hweibull(3, 2, 1)
dweibull(3, 2, 1) / pweibull(3, 2, 1, lower.tail = FALSE)
Hweibull(3, 2, 1)
-pweibull(3, 2, 1, log.p = TRUE, lower.tail = FALSE)

Loglihood function of a Weibull regression

Description

Calculates minus the log likelihood function and its first and second order derivatives for data from a Weibull regression model. Is called by weibreg.

Usage

wfunk(
  beta = NULL,
  lambda,
  p,
  X = NULL,
  Y,
  offset = rep(0, length(Y)),
  ord = 2,
  pfixed = FALSE
)

Arguments

Details

Note that the function returns log likelihood, score vector and minus hessian, i.e. the observed information. The model is

h(t; p, \lambda,\beta, z) = p / \lambda (t / \lambda)^{(p-1)}\exp{(-( t / \lambda)^p})\exp(z\beta)

This is in correspondence with dweibull.

Value

A list with components

Author(s)

Göran Broström

See Also

weibreg