Package {mets}

mets: Analysis of Multivariate Event Times

Description

Implementation of various statistical models for multivariate event history data doi:10.1007/s10985-013-9244-x. Including multivariate cumulative incidence models doi:10.1002/sim.6016, and bivariate random effects probit models (Liability models) doi:10.1016/j.csda.2015.01.014. Modern methods for survival analysis, including regression modelling (Cox, Fine-Gray, Ghosh-Lin, Binomial regression) with fast computation of influence functions.

Implementation of various statistical models for multivariate event history data. Including multivariate cumulative incidence models, and bivariate random effects probit models (Liability models)

Author(s)

Maintainer: Klaus K. Holst klaus@holst.it

Authors:

• Klaus K. Holst klaus@holst.it

• Thomas Scheike

Klaus K. Holst and Thomas Scheike

See Also

Useful links:

• https://kkholst.github.io/mets/

• https://github.com/kkholst/mets

• Report bugs at https://github.com/kkholst/mets/issues

ACTG175, block randomized study from speff2trial package

Description

Data from speff2trial

Format

Randomized study

Source

Hammer et al. 1996, speff2trial package.

Examples

data(ACTG175)

Rates for HPN program for patients of Copenhagen Cohort

Description

Rates for HPN program for patients of Copenhagen Cohort

Format

crbsi: cumulative rate of catheter related bloodstream infection in HPN patients of Copenhagen mechanical: cumulative rate of Mechanical (hole/defect) complication for catheter of HPN patients of Copenhagen trombo: cumulative rate of Occlusion/Thrombosis complication for catheter of HPN patients of Copenhagen terminal: rate of terminal event, patients leaving the HPN program

Source

Estimated data

Clayton-Oakes model with piece-wise constant hazards

Description

Clayton-Oakes frailty model

Usage

ClaytonOakes(
  formula,
  data = parent.frame(),
  cluster,
  var.formula = ~1,
  cuts = NULL,
  type = "piecewise",
  start,
  control = list(),
  var.invlink = exp,
  ...
)

Arguments

Author(s)

Klaus K. Holst

Examples

set.seed(1)
d <- subset(sim_ClaytonOakes(500,4,2,1,stoptime=2,left=2),truncated)
e <- ClaytonOakes(survival::Surv(lefttime,time,status)~x+cluster(~1,cluster),
                  cuts=c(0,0.5,1,2),data=d)
e

d2 <- sim_ClaytonOakes(500,4,2,1,stoptime=2,left=0)
d2$z <- rep(1,nrow(d2)); d2$z[d2$cluster%in%sample(d2$cluster,100)] <- 0
## Marginal=Cox Proportional Hazards model:
## ts <- ClaytonOakes(survival::Surv(time,status)~timereg::prop(x)+cluster(~1,cluster),
##                   data=d2,type="two.stage")
## Marginal=Aalens additive model:
## ts2 <- ClaytonOakes(survival::Surv(time,status)~x+cluster(~1,cluster),
##                    data=d2,type="two.stage")
## Marginal=Piecewise constant:
e2 <- ClaytonOakes(survival::Surv(time,status)~x+cluster(~-1+factor(z),cluster),
                   cuts=c(0,0.5,1,2),data=d2)
e2

e0 <- ClaytonOakes(survival::Surv(time,status)~cluster(~-1+factor(z),cluster),
                   cuts=c(0,0.5,1,2),data=d2)
##ts0 <- ClaytonOakes(survival::Surv(time,status)~cluster(~1,cluster),
##                   data=d2,type="two.stage")
##plot(ts0)
##plot(e0)

e3 <- ClaytonOakes(survival::Surv(time,status)~x+cluster(~1,cluster),cuts=c(0,0.5,1,2),
                   data=d,var.invlink=identity)
e3

Derivatives of the bivariate normal cumulative distribution function

Description

Derivatives of the bivariate normal cumulative distribution function

Usage

Dbvn(p,design=function(p,...) {
     return(list(mu=cbind(p[1],p[1]),
               dmu=cbind(1,1),
               S=matrix(c(p[2],p[3],p[3],p[4]),ncol=2),
               dS=rbind(c(1,0,0,0),c(0,1,1,0),c(0,0,0,1)))  )},                 
     Y=cbind(0,0))

Arguments

Author(s)

Klaus K. Holst

Event history object

Description

Constructur for Event History objects

Usage

Event(time, time2 = TRUE, cause = NULL, cens.code = 0, ...)

Arguments

Details

... content for details

Value

Object of class Event (a matrix)

Author(s)

Klaus K. Holst and Thomas Scheike

Examples

	t1 <- 1:10
	t2 <- t1+runif(10)
	ca <- rbinom(10,2,0.4)
	(x <- Event(t1,t2,ca))

Additive Random effects model for competing risks data for polygenetic modelling

Description

Fits a random effects model describing the dependence in the cumulative incidence curves for subjects within a cluster. Given the gamma distributed random effects it is assumed that the cumulative incidence curves are indpendent, and that the marginal cumulative incidence curves are on additive form

P(T \leq t, cause=1 | x,z) = P_1(t,x,z) = 1- exp( -x^T A(t) - t z^T \beta)

Usage

Grandom_cif(
  cif,
  data,
  cause = NULL,
  cif2 = NULL,
  times = NULL,
  cause1 = 1,
  cause2 = 1,
  cens.code = NULL,
  cens.model = "KM",
  Nit = 40,
  detail = 0,
  clusters = NULL,
  theta = NULL,
  theta.des = NULL,
  weights = NULL,
  step = 1,
  sym = 0,
  same.cens = FALSE,
  censoring.weights = NULL,
  silent = 1,
  var.link = 0,
  score.method = "nr",
  entry = NULL,
  estimator = 1,
  trunkp = 1,
  admin.cens = NULL,
  random.design = NULL,
  ...
)

Arguments

Details

We allow a regression structure for the indenpendent gamma distributed random effects and their variances that may depend on cluster covariates.

random.design specificies the random effects for each subject within a cluster. This is a matrix of 1's and 0's with dimension n x d. With d random effects. For a cluster with two subjects, we let the random.design rows be v_1 and v_2. Such that the random effects for subject 1 is

v_1^T (Z_1,...,Z_d)

, for d random effects. Each random effect has an associated parameter (\lambda_1,...,\lambda_d). By construction subjects 1's random effect are Gamma distributed with mean \lambda_1/v_1^T \lambda and variance \lambda_1/(v_1^T \lambda)^2. Note that the random effect v_1^T (Z_1,...,Z_d) has mean 1 and variance 1/(v_1^T \lambda).

The parameters (\lambda_1,...,\lambda_d) are related to the parameters of the model by a regression construction pard (d x k), that links the d \lambda parameters with the (k) underlying \theta parameters

\lambda = pard \theta

Value

returns an object of type 'random.cif'. With the following arguments:

Author(s)

Thomas Scheike

References

A Semiparametric Random Effects Model for Multivariate Competing Risks Data, Scheike, Zhang, Sun, Jensen (2010), Biometrika.

Cross odds ratio Modelling of dependence for Multivariate Competing Risks Data, Scheike and Sun (2013), Biostatitistics.

Scheike, Holst, Hjelmborg (2014), LIDA, Estimating heritability for cause specific hazards based on twin data

Examples

 ## Reduce Ex.Timings
 d <- sim_nordic_random(1000,delayed=TRUE,
       cordz=1.0,cormz=2,lam0=0.3,country=TRUE)
 times <- seq(50,90,by=10)
 addm <- timereg::comp.risk(Event(time,cause)~-1+factor(country)+cluster(id),data=d,
 times=times,cause=1,max.clust=NULL)

 ### making group indidcator 
 mm <- model.matrix(~-1+factor(zyg),d)

 out1m<-random_cif(addm,data=d,cause1=1,cause2=1,theta=1,
		   theta.des=mm,same.cens=TRUE)
 summary(out1m)
 
 ## this model can also be formulated as a random effects model 
 ## but with different parameters
 out2m<-Grandom_cif(addm,data=d,cause1=1,cause2=1,
		    theta=c(0.5,1),step=1.0,
		    random.design=mm,same.cens=TRUE)
 summary(out2m)
 1/out2m$theta
 out1m$theta
 
 ####################################################################
 ################### ACE modelling of twin data #####################
 ####################################################################
 ### assume that zygbin gives the zygosity of mono and dizygotic twins
 ### 0 for mono and 1 for dizygotic twins. We now formulate and AC model
 zygbin <- d$zyg=="DZ"

 n <- nrow(d)
 ### random effects for each cluster
 des.rv <- cbind(mm,(zygbin==1)*rep(c(1,0)),(zygbin==1)*rep(c(0,1)),1)
 ### design making parameters half the variance for dizygotic components
 pardes <- rbind(c(1,0), c(0.5,0),c(0.5,0), c(0.5,0), c(0,1))

 outacem <-Grandom_cif(addm,data=d,cause1=1,cause2=1,
		same.cens=TRUE,theta=c(0.35,0.15),
            step=1.0,theta.des=pardes,random.design=des.rv)
 summary(outacem)

Influence Functions for phreg objects

Description

Computes the influence functions (IID decomposition) for the regression coefficients and/or the baseline cumulative hazard at a specific time point.

Usage

## S3 method for class 'phreg'
IC(x, type = "robust", all = FALSE, time = NULL, baseline = NULL, ...)

Arguments

Value

A matrix of influence functions. If all=TRUE, columns correspond to regression coefficients and baseline cumulative hazard. Attributes include coef and time.

Author(s)

Thomas Scheike

See Also

phreg, iidBaseline

Simple linear spline

Description

Simple linear spline

Usage

LinSpline(x, knots, num = TRUE, name = "Spline")

Arguments

Author(s)

Thomas Scheike

The TRACE study group of myocardial infarction

Description

The TRACE data frame contains 1877 patients and is a subset of a data set consisting of approximately 6000 patients. It contains data relating survival of patients after myocardial infarction to various risk factors.

Format

This data frame contains the following columns:

id

a numeric vector. Patient code.

status

a numeric vector code. Survival status. 9: dead from myocardial infarction, 0: alive, 7: dead from other causes.

time

a numeric vector. Survival time in years.

chf

a numeric vector code. Clinical heart pump failure, 1: present, 0: absent.

diabetes

a numeric vector code. Diabetes, 1: present, 0: absent.

vf

a numeric vector code. Ventricular fibrillation, 1: present, 0: absent.

wmi

a numeric vector. Measure of heart pumping effect based on ultrasound measurements where 2 is normal and 0 is worst.

sex

a numeric vector code. 1: female, 0: male.

age

a numeric vector code. Age of patient.

Details

sTRACE is a subsample consisting of 300 patients.

tTRACE is a subsample consisting of 1000 patients.

Source

The TRACE study group.

Jensen, G.V., Torp-Pedersen, C., Hildebrandt, P., Kober, L., F. E. Nielsen, Melchior, T., Joen, T. and P. K. Andersen (1997), Does in-hospital ventricular fibrillation affect prognosis after myocardial infarction?, European Heart Journal 18, 919–924.

Examples

data(TRACE)
names(TRACE)

While-Alive Estimands for Recurrent Events

Description

Computes the "While-Alive" estimands for recurrent events in the presence of a terminal event (death). These estimands address the challenge of defining meaningful treatment effects when death prevents further observation of recurrent events.

Usage

WA_recurrent(
  formula,
  data,
  time = NULL,
  cens.code = 0,
  cause = 1,
  death.code = 2,
  trans = NULL,
  cens.formula = NULL,
  augmentR = NULL,
  augmentC = NULL,
  type = NULL,
  marks = NULL,
  ...
)

Arguments

Details

The function estimates two primary quantities:

1. Ratio of Means:

   E(N(\min(D,t))) / E(\min(D,t))

   The expected number of events up to time t (censored by death D) divided by the expected time alive up to t.

2. Mean of Events per Time Unit:

   E(N(\min(D,t)) / \min(D,t))

   The expected rate of events per unit of time alive.

Estimation is based on Inverse Probability of Censoring Weighting (IPCW) to handle administrative censoring and death. The method can be augmented with covariates (double robust estimation) to improve efficiency and robustness.

Value

An object of class "WA" containing:

The object includes influence functions (IID) for all estimators, allowing for further variance calculations or combination with other estimators.

Author(s)

Thomas Scheike

References

Ragni, A., Martinussen, T., & Scheike, T. H. (2023). Nonparametric estimation of the Patient Weighted While-Alive Estimand. arXiv preprint.

Mao, L. (2023). Nonparametric inference of general while-alive estimands for recurrent events. Biometrics, 79(3), 1749–1760.

Schmidli, H., Roger, J. H., & Akacha, M. (2023). Estimands for recurrent event endpoints in the presence of a terminal event. Statistics in Biopharmaceutical Research, 15(2), 238–248.

Examples

data(hfactioncpx12)
dtable(hfactioncpx12,~status)
dd <- WA_recurrent(Event(entry,time,status)~treatment+cluster(id),data=hfactioncpx12,
                   time=2,death.code=2)
summary(dd)

While-Alive Regression for Recurrent Events

Description

Performs regression analysis for the "While-Alive" mean of events per time unit, defined as Z(t) = N(\min(D,t)) / \min(D,t). This function models how covariates affect the rate of recurrent events per unit of time alive.

Usage

WA_reg(
  formula,
  data,
  time = NULL,
  cens.code = 0,
  cause = 1,
  death.code = 2,
  marks = NULL,
  ...,
  trans = 1
)

Arguments

Details

The estimation is based on IPCW (Inverse Probability of Censoring Weighting) and calls binreg after constructing the outcome variable. It supports double robust estimation if covariate augmentation is specified.

Value

An object of class "binreg" containing coefficient estimates, standard errors, confidence intervals, and influence functions for the regression of the event rate per time alive.

Author(s)

Thomas Scheike

References

Ragni, A., Martinussen, T., & Scheike, T. H. (2023). Nonparametric estimation of the Patient Weighted While-Alive Estimand. arXiv preprint.

Examples

data(hfactioncpx12)
hfactioncpx12$age <- rnorm(741)[hfactioncpx12$id] 
dtable(hfactioncpx12,~status)
## exp-link regression 
dd <- WA_reg(Event(entry,time,status)~treatment+age+cluster(id),data=hfactioncpx12,
                    time=2,death.code=2)
summary(dd)

Fast Additive Hazards Model with Robust Standard Errors

Description

Fits a fast Lin-Ying additive hazards model with a possibly stratified baseline. Robust variance is the default variance estimate in the summary.

Usage

aalenMets(formula, data = data, no.baseline = FALSE, ...)

Arguments

Details

Influence functions (IID) follow the numerical order of the given cluster variable. Ordering by $id aligns the IID terms with the dataset order.

Value

An object of class "aalenMets" (extends "phreg") containing:

Author(s)

Thomas Scheike

Examples

data(bmt)
bmt$time <- bmt$time + runif(408) * 0.001
out <- aalenMets(Surv(time, cause == 1) ~ tcell + platelet + age, data = bmt)
summary(out)

## Comparison with timereg::aalen
## out2 <- timereg::aalen(
##   Surv(time, cause == 1) ~ const(tcell) + const(platelet) + const(age),
##   data = bmt
## )
## summary(out2)

Estimation of Concordance in Bivariate Competing Risks Data

Description

Estimates the bivariate cumulative incidence function (concordance) for paired data (e.g., twins, family members) in the presence of competing risks. The function handles both the IPCW (Inverse Probability of Censoring Weighting) estimator and the Aalen-Johansen estimator (via prodlim).

Usage

bicomprisk(
  formula,
  data,
  cause = c(1, 1),
  cens = 0,
  causes,
  indiv,
  strata = NULL,
  id,
  num,
  max.clust = 1000,
  marg = NULL,
  se.clusters = NULL,
  wname = NULL,
  prodlim = FALSE,
  messages = TRUE,
  model,
  return.data = 0,
  uniform = 0,
  conservative = 1,
  resample.iid = 1,
  ...
)

Arguments

Details

The concordance function C(t) is defined as the probability that both members of a pair experience the event of interest by time t:

C(t) = P(T_1 \leq t, T_2 \leq t, \epsilon_1 = k, \epsilon_2 = k)

where T_i is the event time and \epsilon_i is the cause of failure for individual i.

The function supports:

• Stratified analysis (e.g., by zygosity in twin studies).

• Clustering for robust standard errors.

• Both IPCW and Aalen-Johansen estimation methods.

• Resampling for uniform standard errors.

• IID decomposition for further inference (e.g., casewise concordance tests).

Value

An object of class "bicomprisk" (or "bicomprisk.strata" if stratified) containing:

Author(s)

Thomas Scheike, Klaus K. Holst

References

Scheike, T. H.; Holst, K. K. & Hjelmborg, J. B. (2014). Estimating twin concordance for bivariate competing risks twin data. Statistics in Medicine, 33, 1193-1204.

See Also

test_casewise, casewise

Examples

library("timereg")

## Simulated data example
prt <- sim_nordic_random(2000,delayed=TRUE,ptrunc=0.7,
	      cordz=0.5,cormz=2,lam0=0.3)
## Bivariate competing risk, concordance estimates
p11 <- bicomprisk(Event(time,cause)~strata(zyg)+id(id),data=prt,cause=c(1,1))

p11mz <- p11$model$"MZ"
p11dz <- p11$model$"DZ"
par(mfrow=c(1,2))
## Concordance
plot(p11mz,ylim=c(0,0.1));
plot(p11dz,ylim=c(0,0.1));

## Entry time, truncation weighting
### Other weighting procedure
prtl <-  prt[!prt$truncated,]
prt2 <- ipw2(prtl,cluster="id",same.cens=TRUE,
     time="time",cause="cause",entrytime="entry",
     pairs=TRUE,strata="zyg",obs.only=TRUE)

prt22 <- fast.reshape(prt2,id="id")

prt22$event <- (prt22$cause1==1)*(prt22$cause2==1)*1
prt22$timel <- pmax(prt22$time1,prt22$time2)
ipwc <- timereg::comp.risk(Event(timel,event)~-1+factor(zyg1),
  data=prt22,cause=1,n.sim=0,model="rcif2",times=50:90,
  weights=prt22$weights1,cens.weights=rep(1,nrow(prt22)))

p11wmz <- ipwc$cum[,2]
p11wdz <- ipwc$cum[,3]
lines(ipwc$cum[,1],p11wmz,col=3)
lines(ipwc$cum[,1],p11wdz,col=3)

Fits Clayton-Oakes or bivariate Plackett (OR) models for binary data using marginals that are on logistic form. If clusters contain more than two times, the algoritm uses a compososite likelihood based on all pairwise bivariate models.

Description

The pairwise pairwise odds ratio model provides an alternative to the alternating logistic regression (ALR).

Usage

binomial_twostage(
  margbin,
  data = parent.frame(),
  method = "nr",
  detail = 0,
  clusters = NULL,
  silent = 1,
  weights = NULL,
  theta = NULL,
  theta.des = NULL,
  var.link = 0,
  var.par = 1,
  var.func = NULL,
  iid = 1,
  notaylor = 1,
  model = "plackett",
  marginal.p = NULL,
  beta.iid = NULL,
  Dbeta.iid = NULL,
  strata = NULL,
  max.clust = NULL,
  se.clusters = NULL,
  numDeriv = 0,
  random.design = NULL,
  pairs = NULL,
  dim.theta = NULL,
  additive.gamma.sum = NULL,
  pair.ascertained = 0,
  case.control = 0,
  no.opt = FALSE,
  twostage = 1,
  beta = NULL,
  ...
)

Arguments

Details

The reported standard errors are based on a cluster corrected score equations from the pairwise likelihoods assuming that the marginals are known. This gives correct standard errors in the case of the Odds-Ratio model (Plackett distribution) for dependence, but incorrect standard errors for the Clayton-Oakes types model (that is also called "gamma"-frailty). For the additive gamma version of the standard errors are adjusted for the uncertainty in the marginal models via an iid deomposition using the iid() function of lava. For the clayton oakes model that is not speicifed via the random effects these can be fixed subsequently using the iid influence functions for the marginal model, but typically this does not change much.

For the Clayton-Oakes version of the model, given the gamma distributed random effects it is assumed that the probabilities are indpendent, and that the marginal survival functions are on logistic form

logit(P(Y=1|X)) = \alpha + x^T \beta

therefore conditional on the random effect the probability of the event is

logit(P(Y=1|X,Z)) = exp( -Z \cdot Laplace^{-1}(lamtot,lamtot,P(Y=1|x)) )

Can also fit a structured additive gamma random effects model, such the ACE, ADE model for survival data:

Now random.design specificies the random effects for each subject within a cluster. This is a matrix of 1's and 0's with dimension n x d. With d random effects. For a cluster with two subjects, we let the random.design rows be v_1 and v_2. Such that the random effects for subject 1 is

v_1^T (Z_1,...,Z_d)

, for d random effects. Each random effect has an associated parameter (\lambda_1,...,\lambda_d). By construction subjects 1's random effect are Gamma distributed with mean \lambda_j/v_1^T \lambda and variance \lambda_j/(v_1^T \lambda)^2. Note that the random effect v_1^T (Z_1,...,Z_d) has mean 1 and variance 1/(v_1^T \lambda). It is here asssumed that lamtot=v_1^T \lambda is fixed over all clusters as it would be for the ACE model below.

The DEFAULT parametrization uses the variances of the random effecs (var.par=1)

\theta_j = \lambda_j/(v_1^T \lambda)^2

For alternative parametrizations (var.par=0) one can specify how the parameters relate to \lambda_j with the function

Based on these parameters the relative contribution (the heritability, h) is equivalent to the expected values of the random effects \lambda_j/v_1^T \lambda

Given the random effects the probabilities are independent and on the form

logit(P(Y=1|X)) = exp( - Laplace^{-1}(lamtot,lamtot,P(Y=1|x)) )

with the inverse laplace of the gamma distribution with mean 1 and variance lamtot.

The parameters (\lambda_1,...,\lambda_d) are related to the parameters of the model by a regression construction pard (d x k), that links the d \lambda parameters with the (k) underlying \theta parameters

\lambda = theta.des \theta

here using theta.des to specify these low-dimension association. Default is a diagonal matrix.

Author(s)

Thomas Scheike

References

Two-stage binomial modelling

Examples

data(twinstut)
twinstut0 <- subset(twinstut, tvparnr<4000)
twinstut <- twinstut0
twinstut$binstut <- (twinstut$stutter=="yes")*1
theta.des <- model.matrix( ~-1+factor(zyg),data=twinstut)
margbin <- glm(binstut~factor(sex)+age,data=twinstut,family=binomial())
bin <- binomial_twostage(margbin,data=twinstut,var.link=1,
         clusters=twinstut$tvparnr,theta.des=theta.des,detail=0)
summary(bin)

twinstut$cage <- scale(twinstut$age)
theta.des <- model.matrix( ~-1+factor(zyg)+cage,data=twinstut)
bina <- binomial_twostage(margbin,data=twinstut,var.link=1,
		         clusters=twinstut$tvparnr,theta.des=theta.des)
summary(bina)

theta.des <- model.matrix( ~-1+factor(zyg)+factor(zyg)*cage,data=twinstut)
bina <- binomial_twostage(margbin,data=twinstut,var.link=1,
		         clusters=twinstut$tvparnr,theta.des=theta.des)
summary(bina)

### use of clayton oakes binomial additive gamma model
###########################################################
 ## Reduce Ex.Timings
data <- sim_binClaytonOakes_family_ace(1000,2,1,beta=NULL,alpha=NULL)
margbin <- glm(ybin~x,data=data,family=binomial())
margbin

head(data)
data$number <- c(1,2,3,4)
data$child <- 1*(data$number==3)

### make ace random effects design
out <- ace_family_design(data,member="type",id="cluster")
out$pardes
head(out$des.rv)

bints <- binomial_twostage(margbin,data=data,
     clusters=data$cluster,detail=0,var.par=1,
     theta=c(2,1),var.link=0,
     random.design=out$des.rv,theta.des=out$pardes)
summary(bints)

data <- sim_binClaytonOakes_twin_ace(1000,2,1,beta=NULL,alpha=NULL)
out  <- twin_polygen_design(data,id="cluster",zygname="zygosity")
out$pardes
head(out$des.rv)
margbin <- glm(ybin~x,data=data,family=binomial())

bintwin <- binomial_twostage(margbin,data=data,
     clusters=data$cluster,var.par=1,
     theta=c(2,1),random.design=out$des.rv,theta.des=out$pardes)
summary(bintwin)
concordanceTwinACE(bintwin)

Binomial Regression for Censored Competing Risks Data

Description

Fits a binomial regression model for a specific time point in the presence of right-censored data and competing risks. This function implements the Inverse Probability of Censoring Weighted (IPCW) estimating equation approach.

Usage

binreg(
  formula,
  data,
  cause = 1,
  time = NULL,
  beta = NULL,
  type = c("II", "I"),
  offset = NULL,
  weights = NULL,
  cens.weights = NULL,
  cens.model = ~+1,
  se = TRUE,
  kaplan.meier = TRUE,
  cens.code = 0,
  no.opt = FALSE,
  method = "nr",
  augmentation = NULL,
  outcome = c("cif", "rmst", "rmtl"),
  model = c("default", "logit", "exp", "lin"),
  Ydirect = NULL,
  ...
)

Arguments

Details

The model estimates the probability:

P(T \leq t, \epsilon=1 | X ) = \text{expit}( X^T \beta)

Based on binomial regresion IPCW response estimating equation:

X ( \Delta^{ipcw}(t) I(T \leq t, \epsilon=1 ) - expit( X^T beta)) = 0

where

\Delta^{ipcw}(t) = I((min(t,T)< C)/G_c(min(t,T)-)

is IPCW adjustment of the response

Y(t)= I(T \leq t, \epsilon=1 )

. Two types of estimators are available:

• type="I": Solves the standard IPCW estimating equation.

• type="II": Adds a censoring augmentation term for efficiency gains, solving

  X \int E(Y(t)| T>s)/G_c(s) d \hat M_c

  .

Alternatively, logitIPCW performs a standard logistic regression with IPCW weights applied directly to the likelihood. Thus solving

X I(min(T_i,t) < G_i)/G_c(min(T_i ,t)) ( I(T \leq t, \epsilon=1 ) - expit( X^T beta)) = 0

a standard logistic regression with IPCW weights.

The variance estimation is based on squared influence functions, with options for naive variance (assuming known censoring) and robust variance (accounting for censoring model estimation).

Censoring model may depend on strata (cens.model=~strata(gX)).

Value

An object of class "binreg" containing coefficients, variance-covariance matrix, influence functions (iid), and model details.

References

• Blanche PF, Holt A, Scheike T (2022). "On logistic regression with right censored data, with or without competing risks, and its use for estimating treatment effects." Lifetime data analysis, 29, 441–482.

• Scheike TH, Zhang MJ, Gerds TA (2008). "Predicting cumulative incidence probability by direct binomial regression." Biometrika, 95(1), 205–220.

Author(s)

Thomas Scheike

Examples

data(bmt); bmt$time <- bmt$time+runif(408)*0.001
# logistic regresion with IPCW binomial regression 
out <- binreg(Event(time,cause)~tcell+platelet,bmt,time=50)
summary(out)
head(iid(out))

predict(out,data.frame(tcell=c(0,1),platelet=c(1,1)),se=TRUE)

outs <- binreg(Event(time,cause)~tcell+platelet,bmt,time=50,cens.model=~strata(tcell,platelet))
summary(outs)

## glm with IPCW weights 
outl <- logitIPCW(Event(time,cause)~tcell+platelet,bmt,time=50)
summary(outl)

##########################################
### risk-ratio of different causes #######
##########################################
data(bmt)
bmt$id <- 1:nrow(bmt)
bmt$status <- bmt$cause
bmt$strata <- 1
bmtdob <- bmt
bmtdob$strata <-2
bmtdob <- dtransform(bmtdob,status=1,cause==2)
bmtdob <- dtransform(bmtdob,status=2,cause==1)
###
bmtdob <- rbind(bmt,bmtdob)
dtable(bmtdob,cause+status~strata)

cif1 <- cif(Event(time,cause)~+1,bmt,cause=1)
cif2 <- cif(Event(time,cause)~+1,bmt,cause=2)
plot(cif1)
plot(cif2,add=TRUE,col=2)

cifs1 <- binreg(Event(time,cause)~tcell+platelet+age,bmt,cause=1,time=50)
cifs2 <- binreg(Event(time,cause)~tcell+platelet+age,bmt,cause=2,time=50)
summary(cifs1)
summary(cifs2)

cifdob <- binreg(Event(time,status)~-1+factor(strata)+
	 tcell*factor(strata)+platelet*factor(strata)+age*factor(strata)
	 +cluster(id),bmtdob,cause=1,time=50,cens.model=~strata(strata))
summary(cifdob)
head(iid(cifdob)) 

newdata <- data.frame(tcell=1,platelet=1,age=0,strata=1:2,id=1)
riskratio <- function(p) {
  cifdob$coef <- p
  p <- predict(cifdob,newdata,se=0)
  return(p[1]/p[2])
}
lava::estimate(cifdob,f=riskratio)

predict(cifdob,newdata)
(p1 <- predict(cifs1,newdata))
(p2 <- predict(cifs2,newdata))
p1[1,1]/p2[1,1]

Average Treatment Effect for Censored Competing Risks Data using Binomial Regression

Description

Estimates the average treatment effect (ATE) E(Y(1) - Y(0)) for censored competing risks data using binomial regression with Inverse Probability of Censoring Weighting (IPCW).

Usage

binregATE(
  formula,
  data,
  cause = 1,
  time = NULL,
  beta = NULL,
  treat.model = ~+1,
  cens.model = ~+1,
  offset = NULL,
  weights = NULL,
  cens.weights = NULL,
  se = TRUE,
  type = c("II", "I"),
  kaplan.meier = TRUE,
  cens.code = 0,
  no.opt = FALSE,
  method = "nr",
  augmentation = NULL,
  outcome = c("cif", "rmst", "rmtl"),
  model = c("default", "logit", "exp", "lin"),
  Ydirect = NULL,
  typeATE = "II",
  ...
)

Arguments

Details

Under standard causal assumptions, the ATE can be estimated. These assumptions include:

• Consistency: The observed outcome equals the potential outcome under the observed treatment.

• Ignorability: (Y(1), Y(0)) \perp A | X (treatment assignment is independent of potential outcomes given covariates).

• Positivity: All treatment levels have non-zero probability given covariates.

The first covariate in the competing risks regression model must be the treatment variable, which should be coded as a factor. If the factor has more than two levels, multinomial logistic regression (mlogit) is used for propensity score modeling. In the absence of censoring, this reduces to ordinary logistic regression.

The ATE is estimated using standard doubly robust estimating equations that are IPCW-censoring adjusted. As an alternative to binomial regression, logitIPCWATE provides an IPCW-weighted version of standard logistic regression.

When typeATE = "II", the estimating equation is augmented with:

(A/\pi(X)) \int E( O(t) | T \geq t, S(X))/ G_c(t,S(X)) d \hat M_c(s)

when estimating the mean outcome for the treated group.

Value

An object of class c("binreg", "ATE") containing:

References

• Blanche PF, Holt A, Scheike T (2022). "On logistic regression with right censored data, with or without competing risks, and its use for estimating treatment effects." Lifetime Data Analysis, 29, 441–482.

Author(s)

Thomas Scheike

See Also

binreg, logitIPCWATE, logitATE, binregG

[kumarsim()] [kumarsimRCT()]

Examples

data(bmt)
dfactor(bmt)  <-  ~.

brs <- binregATE(Event(time,cause)~tcell.f+platelet+age,bmt,time=50,cause=1,
  treat.model=tcell.f~platelet+age)
summary(brs)
head(brs$riskDR.iid)
head(brs$riskG.iid)

brsi <- binregATE(Event(time,cause)~tcell.f+tcell.f*platelet+tcell.f*age,bmt,time=50,cause=1,
  treat.model=tcell.f~platelet+age)
summary(brsi)
head(brs$riskDR.iid)
head(brs$riskG.iid)

Estimate Casewise Concordance Using Binomial Regression

Description

Estimates the casewise concordance based on concordance and marginal estimates obtained from binreg objects. Uses cluster-based IID for standard errors, which are often better than those from casewise (which can be conservative).

Usage

binregCasewise(concbreg, margbreg, zygs = c("DZ", "MZ"), newdata = NULL, ...)

Arguments

Value

A list containing:

Author(s)

Thomas Scheike

Examples

data(prt)
prt <- force_same_cens(prt,cause="status")

dd <- bicompriskData(Event(time, status)~strata(zyg)+id(id), data=prt, cause=c(2, 2))
newdata <- data.frame(zyg=c("DZ","MZ"),id=1)

## concordance 
bcif1 <- binreg(Event(time,status)~-1+factor(zyg)+cluster(id), data=dd,
                time=80, cause=1, cens.model=~strata(zyg))
pconc <- predict(bcif1,newdata)

## marginal estimates 
mbcif1 <- binreg(Event(time,status)~cluster(id), data=prt, time=80, cause=2)
mc <- predict(mbcif1,newdata)

cse <- binregCasewise(bcif1,mbcif1)
cse

G-Estimator for Binomial Regression Model (Standardized Estimates)

Description

Computes the G-estimator (G-formula) for standardized risk estimates based on a fitted binreg object. The G-estimator standardizes predictions over the covariate distribution in the data:

\hat F(t, A=a) = n^{-1} \sum_{i=1}^n \hat F(t, A=a, Z_i)

Usage

binregG(x, data, Avalues = NULL, varname = NULL)

Arguments

Details

This function assumes that the first covariate in the original model formula represents the treatment or exposure variable (A). It calculates the marginal risk for specified values of A by averaging the conditional predictions over the observed covariate distribution Z.

The function returns influence functions for these risk estimates, allowing for the computation of standard errors and confidence intervals.

If the first covariate is a factor, contrasts between all levels are computed automatically. If it is continuous, specific values must be provided via Avalues.

Value

An object of class "survivalG" containing:

References

• Blanche PF, Holt A, Scheike T (2022). "On logistic regression with right censored data, with or without competing risks, and its use for estimating treatment effects." Lifetime Data Analysis, 29, 441–482.

Author(s)

Thomas Scheike

See Also

binreg, binregATE

Examples

data(bmt); bmt$time <- bmt$time+runif(408)*0.001
bmt$event <- (bmt$cause!=0)*1

b1 <- binreg(Event(time,cause)~age+tcell+platelet,bmt,cause=1,time=50)
sb1 <- binregG(b1,bmt,Avalues=c(0,1,2))
summary(sb1)

Percentage of Years Lost Due to a Cause Regression

Description

Estimates the percentage of the restricted mean time lost (RMTL) that is attributable to a specific cause and models how covariates affect this percentage using IPCW regression.

Usage

binregRatio(
  formula,
  data,
  cause = 1,
  time = NULL,
  beta = NULL,
  type = c("III", "II", "I"),
  offset = NULL,
  weights = NULL,
  cens.weights = NULL,
  cens.model = ~+1,
  se = TRUE,
  relative.to.causes = NULL,
  kaplan.meier = TRUE,
  cens.code = 0,
  no.opt = FALSE,
  method = "nr",
  augmentation = NULL,
  outcome = c("rmtl", "cif"),
  model = c("logit", "exp", "lin"),
  Ydirect = NULL,
  ...
)

Arguments

Details

Let the total years lost be Y = t - \min(T, t) and the years lost due to cause 1 be Y_1 = I(\epsilon=1) (t - \min(T, t)). The function models the ratio:

\text{logit}\left( \frac{E(Y_1 | X)}{E(Y | X)} \right) = X^T \beta

Estimation is based on a binomial regression IPCW response estimating equation:

X \left( \Delta^{\text{ipcw}}(t) \left( Y \cdot \text{expit}(X^T \beta) - Y_1 \right) \right) = 0

where \Delta^{\text{ipcw}}(t) = I(\min(t,T) < C) / G_c(\min(t,T)) is the IPCW adjustment.

The function supports three types of estimators:

• "I": Classical outcome IPCW regression (no augmentation).

• "II": Adds a censoring augmentation term X \int E(Z(t)| T>s)/G_c(s) d \hat M_c to improve efficiency (requires an initial estimate of \beta).

• "III": Adds a more complex augmentation term separating the expectations of Y and Y_1 for further efficiency gains.

The variance is based on the squared influence functions (IID). A "naive" variance (assuming known censoring) is also provided for comparison.

Value

An object of class "binreg" and "ratio" containing:

Author(s)

Thomas Scheike

References

Scheike, T. & Tanaka, S. (2025). Restricted mean time lost ratio regression: Percentage of restricted mean time lost due to specific cause. WIP.

See Also

resmeanIPCW, binreg

Examples

data(bmt); bmt$time <- bmt$time+runif(408)*0.001

rmtl30 <- rmstIPCW(Event(time,cause!=0)~platelet+tcell+age, bmt, time=30, cause=1, outcome="rmtl")
rmtl301 <- rmstIPCW(Event(time,cause)~platelet+tcell+age, bmt, time=30, cause=1)
rmtl302 <- rmstIPCW(Event(time,cause)~platelet+tcell+age, bmt, time=30, cause=2)

estimate(rmtl30)
estimate(rmtl301)
estimate(rmtl302)

## Percentage of total RMTL due to cause 1
rmtlratioI <- rmtlRatio(Event(time,cause)~platelet+tcell+age, bmt, time=30, cause=1)
summary(rmtlratioI)

newdata <- data.frame(platelet=1, tcell=1, age=1)
pp <- predict(rmtlratioI, newdata)
pp

## Percentage of total cumulative incidence due to cause 1
cifratio <- binregRatio(Event(time,cause)~platelet+tcell+age, bmt, time=30, cause=1, model="cif")
summary(cifratio)
pp <- predict(cifratio, newdata)
pp

Two-Stage Randomization for Survival or Competing Risks Data

Description

Estimates the average treatment effect E(Y(i,j)) of treatment regime (i,j) under two-stage randomization. The estimator can be augmented using information from both randomizations and dynamic censoring augmentation to improve efficiency.

Usage

binregTSR(
  formula,
  data,
  cause = 1,
  time = NULL,
  cens.code = 0,
  response.code = NULL,
  augmentR0 = NULL,
  treat.model0 = ~+1,
  augmentR1 = NULL,
  treat.model1 = ~+1,
  augmentC = NULL,
  cens.model = ~+1,
  estpr = c(1, 1),
  response.name = NULL,
  offset = NULL,
  weights = NULL,
  cens.weights = NULL,
  beta = NULL,
  kaplan.meier = TRUE,
  no.opt = FALSE,
  method = "nr",
  augmentation = NULL,
  outcome = c("cif", "rmst", "rmst-cause"),
  model = "exp",
  Ydirect = NULL,
  return.dataw = 0,
  pi0 = 0.5,
  pi1 = 0.5,
  cens.time.fixed = 1,
  outcome.iid = 1,
  meanCs = 0,
  ...
)

Arguments

Details

The method solves the estimating equation:

\frac{I(\min(T_i,t) < G_i)}{G_c(\min(T_i,t))} I(T \leq t, \epsilon=1) - AUG_0 - AUG_1 + AUG_C - p(i,j) = 0

where:

• AUG_0 = \frac{A_0(i) - \pi_0(i)}{\pi_0(i)} X_0 \gamma_0 uses covariates from augmentR0

• AUG_1 = \frac{A_0(i)}{\pi_0(i)} \frac{A_1(j) - \pi_1(j)}{\pi_1(j)} X_1 \gamma_1 uses covariates from augmentR1

• AUG_C = \int_0^t \gamma_c(s)^T (e(s) - \bar e(s)) \frac{1}{G_c(s)} dM_c(s) is the censoring augmentation

Standard errors are estimated using the influence function of all estimators, enabling tests of differences to be computed subsequently. The method handles both survival data and competing risks data, and supports multiple treatment levels.

Value

An object of class "binregTSR" containing:

Author(s)

Thomas Scheike

References

Scheike, T. H. (2024). Two-stage randomization analysis for survival data. mets package documentation.

See Also

binreg, phreg_rct

Examples

ddf <- mets:::gsim(200,covs=1,null=0,cens=1,ce=2)

bb <- binregTSR(Event(entry,time,status)~+1+cluster(id),ddf$datat,time=2,cause=c(1),
        cens.code=0,treat.model0=A0.f~+1,treat.model1=A1.f~A0.f,
        augmentR1=~X11+X12+TR,augmentR0=~X01+X02,
        augmentC=~A01+A02+X01+X02+A11t+A12t+X11+X12+TR,
        response.code=2)
summary(bb) 

Bivariate Probit model

Description

Bivariate Probit model

Usage

biprobit(
  x,
  data,
  id,
  rho = ~1,
  num = NULL,
  strata = NULL,
  eqmarg = TRUE,
  indep = FALSE,
  weights = NULL,
  weights.fun = function(x) ifelse(any(x <= 0), 0, max(x)),
  randomeffect = FALSE,
  vcov = "robust",
  pairs.only = FALSE,
  allmarg = !is.null(weights),
  control = list(trace = 0),
  messages = 1,
  constrain = NULL,
  table = pairs.only,
  p = NULL,
  ...
)

Arguments

Examples

data(prt)
prt0 <- subset(prt,country=="Denmark")
a <- biprobit(cancer~1+zyg, ~1+zyg, data=prt0, id="id")
predict(a, newdata=lava::Expand(prt, zyg=c("MZ")))
## b <- biprobit(cancer~1+zyg, ~1+zyg, data=prt0, id="id",pairs.only=TRUE)
## predict(b,newdata=lava::Expand(prt,zyg=c("MZ","DZ")))

 ## Reduce Ex.Timings
n <- 2e3
x <- sort(runif(n, -1, 1))
y <- rmvn(n, c(0,0), rho=cbind(tanh(x)))>0
d <- data.frame(y1=y[,1], y2=y[,2], x=x)
dd <- fast.reshape(d)

a <- biprobit(y~1+x,rho=~1+x,data=dd,id="id")
summary(a, mean.contrast=c(1,.5), cor.contrast=c(1,.5))
with(predict(a,data.frame(x=seq(-1,1,by=.1))), plot(p00~x,type="l"))

pp <- predict(a,data.frame(x=seq(-1,1,by=.1)),which=c(1))
plot(pp[,1]~pp$x, type="l", xlab="x", ylab="Concordance", lwd=2, xaxs="i")
lava::confband(pp$x,pp[,2],pp[,3],polygon=TRUE,lty=0,col=lava::Col(1))

pp <- predict(a,data.frame(x=seq(-1,1,by=.1)),which=c(9)) ## rho
plot(pp[,1]~pp$x, type="l", xlab="x", ylab="Correlation", lwd=2, xaxs="i")
lava::confband(pp$x,pp[,2],pp[,3],polygon=TRUE,lty=0,col=lava::Col(1))
with(pp, lines(x,tanh(-x),lwd=2,lty=2))

xp <- seq(-1,1,length.out=6); delta <- mean(diff(xp))
a2 <- biprobit(y~1+x,rho=~1+I(cut(x,breaks=xp)),data=dd,id="id")
pp2 <- predict(a2,data.frame(x=xp[-1]-delta/2),which=c(9)) ## rho
lava::confband(pp2$x,pp2[,2],pp2[,3],center=pp2[,1])




## Time
## Not run: 
    a <- biprobit.time(cancer~1, rho=~1+zyg, id="id", data=prt, eqmarg=TRUE,
                       cens.formula=Surv(time,status==0)~1,
                       breaks=seq(75,100,by=3),fix.censweights=TRUE)

    a <- biprobit.time2(cancer~1+zyg, rho=~1+zyg, id="id", data=prt0, eqmarg=TRUE,
                       cens.formula=Surv(time,status==0)~zyg,
                       breaks=100)

    #a1 <- biprobit.time2(cancer~1, rho=~1, id="id", data=subset(prt0,zyg=="MZ"), eqmarg=TRUE,
    #                   cens.formula=Surv(time,status==0)~1,
    #                   breaks=100,pairs.only=TRUE)

    #a2 <- biprobit.time2(cancer~1, rho=~1, id="id", data=subset(prt0,zyg=="DZ"), eqmarg=TRUE,
    #                    cens.formula=Surv(time,status==0)~1,
    #                    breaks=100,pairs.only=TRUE)

## End(Not run)

Block sampling

Description

Sample blockwise from clustered data

Usage

blocksample(data, size, idvar = NULL, replace = TRUE, ...)

Arguments

Details

Original id is stored in the attribute 'id'

Value

data.frame

Author(s)

Klaus K. Holst

Examples

d <- data.frame(x=rnorm(5), z=rnorm(5), id=c(4,10,10,5,5), v=rnorm(5))
(dd <- blocksample(d,size=20,~id))
attributes(dd)$id

## Not run: 
blocksample(data.table::data.table(d),1e6,~id)

## End(Not run)


d <- data.frame(x=c(1,rnorm(9)),
               z=rnorm(10),
               id=c(4,10,10,5,5,4,4,5,10,5),
               id2=c(1,1,2,1,2,1,1,1,1,2),
               v=rnorm(10))
dsample(d,~id, size=2)
dsample(d,.~id+id2)
dsample(d,x+z~id|x>0,size=5)

The Bone Marrow Transplant Data

Description

Bone marrow transplant data with 408 rows and 5 columns.

Format

The data has 408 rows and 5 columns.

cause

a numeric vector code. Survival status. 1: dead from treatment related causes, 2: relapse , 0: censored.

time

a numeric vector. Survival time.

platelet

a numeric vector code. Plalelet 1: more than 100 x 10^9 per L, 0: less.

tcell

a numeric vector. T-cell depleted BMT 1:yes, 0:no.

age

a numeric vector code. Age of patient, scaled and centered ((age-35)/15).

Source

Simulated data

References

NN

Examples

data(bmt)
names(bmt)

Liability model for twin data

Description

Liability-threshold model for twin data

Usage

bptwin(
  x,
  data,
  id,
  zyg,
  DZ,
  group = NULL,
  num = NULL,
  weights = NULL,
  weights.fun = function(x) ifelse(any(x <= 0), 0, max(x)),
  strata = NULL,
  messages = 1,
  control = list(trace = 0),
  type = "ace",
  eqmean = TRUE,
  pairs.only = FALSE,
  samecens = TRUE,
  allmarg = samecens & !is.null(weights),
  stderr = TRUE,
  robustvar = TRUE,
  p,
  indiv = FALSE,
  constrain,
  varlink,
  ...
)

Arguments

Author(s)

Klaus K. Holst

See Also

twinlm, twinlm.time, twinlm.strata, twinsim

Examples

data(twinstut)
b0 <- bptwin(stutter~sex,
        data=droplevels(
          subset(twinstut, zyg%in%c("mz","dz") & tvparnr<5e3)
        ),
        id="tvparnr",zyg="zyg",DZ="dz",type="ae")
summary(b0)

CALGB 8923, twostage randomization SMART design

Description

Data from CALGB 8923

Format

Data from smart design id: id of subject status : 1-death, 2-response for second randomization, 0-censoring A0 : treament at first randomization A1 : treament at second randomization At.f : treament given at record (A0 or A1) TR : time of response sex : 0-males, 1-females consent: 1 if agrees to 2nd randomization, censored if not R: 1 if response trt1: A0 trt2: A1

Source

https://github.com/ycchao/code_Joint_model_SMART

Examples

data(calgb8923)

Estimate Casewise Concordance from prodlim Objects

Description

Estimates the casewise concordance based on concordance and marginal estimates derived from prodlim objects. Unlike test_casewise, this function does not perform hypothesis testing but focuses on estimation and plotting.

Usage

casewise(conc, marg, cause.marg)

Arguments

Value

An object of class "casewise" containing:

Author(s)

Thomas Scheike

See Also

test_casewise, bicomprisk

Examples

 ## Reduce Ex.Timings
library(prodlim)
data(prt);
prt <- force_same_cens(prt,cause="status")

### marginal cumulative incidence of prostate cancer
outm <- prodlim(Hist(time,status)~+1,data=prt)

times <- 60:100
cifmz <- predict(outm,cause=2,time=times,newdata=data.frame(zyg="MZ"))
cifdz <- predict(outm,cause=2,time=times,newdata=data.frame(zyg="DZ"))

### concordance for MZ and DZ twins
cc <- bicomprisk(Event(time,status)~strata(zyg)+id(id),data=prt,cause=c(2,2),prodlim=TRUE)
cdz <- cc$model$"DZ"
cmz <- cc$model$"MZ"

cdz <- casewise(cdz,outm,cause.marg=2) 
cmz <- casewise(cmz,outm,cause.marg=2)

plot(cmz,ci=NULL,ylim=c(0,0.5),xlim=c(60,100),legend=TRUE,col=c(3,2,1))
par(new=TRUE)
plot(cdz,ci=NULL,ylim=c(0,0.5),xlim=c(60,100),legend=TRUE)
summary(cdz)
summary(cmz)

Casewise Concordance from Concordant/Discordant Counts

Description

Computes casewise concordance probability and confidence interval from counts of concordant and discordant pairs using a binomial GLM.

Usage

casewise_bin(nc, nd)

Arguments

Value

A list with p.casewise (estimated probability) and ci.casewise (confidence interval).

Cumulative Incidence with Robust Standard Errors

Description

Computes cumulative incidence functions with robust standard errors.

Usage

cif(formula, data = data, cause = 1, cens.code = 0, death.code = NULL, ...)

Arguments

Value

An object of class "cif" (extends "phreg") containing:

Author(s)

Thomas Scheike

Examples

data(bmt)
bmt$cluster <- sample(1:100, 408, replace = TRUE)
out1 <- cif(Event(time, cause) ~ +1, data = bmt, cause = 1)
out2 <- cif(Event(time, cause) ~ +1 + cluster(cluster), data = bmt, cause = 1)

par(mfrow = c(1, 2))
plot(out1, se = TRUE)
plot(out2, se = TRUE)

Restricted Mean Time Lost for Competing Risks

Description

Computes the Restricted Mean Time Lost (RMTL) for competing risks based on the integrated Aalen-Johansen estimator.

Usage

cif_yearslost(formula, data = data, cens.code = 0, times = NULL, ...)

Arguments

Details

A set of time points can be given to be returned in the summary, but the function computes years-lost for all event times and can be plotted/viewed. The RMTL for a specific time-point can also be computed using the rmstIPCW function.

Value

An object of class "resmean_phreg" containing:

Author(s)

Thomas Scheike

Examples

data(bmt)
bmt$time <- bmt$time + runif(408) * 0.001

## Years lost decomposed into causes
drm1 <- cif_yearslost(Event(time, cause) ~ strata(tcell, platelet), data = bmt, times = c(40, 50))
par(mfrow = c(1, 2))
plot(drm1, cause = 1, se = 1)
plot(drm1, cause = 2, se = 1)
summary(drm1)
estimate(drm1, cause = 1)
estimate(drm1, cause = 2)

## Comparing populations
drm1 <- cif_yearslost(Event(time, cause) ~ strata(tcell, platelet), data = bmt, times = 40)
summary(drm1, contrast = list(1:4))
e1 <- estimate(drm1)
estimate(e1, rbind(c(1, -1, 0, 0)))

Cumulative Incidence Function (CIF) Regression

Description

Fits a regression model for the cumulative incidence function (CIF) in the presence of competing risks. Supports two link functions:

• propodds=1 (default): Logistic link model (logit of CIF), providing Odds Ratio (OR) interpretations.

• propodds=NULL: Fine-Gray (cloglog) regression model, providing subdistribution hazard ratio interpretations.

Usage

cifreg(
  formula,
  data,
  propodds = 1,
  cause = 1,
  cens.code = 0,
  no.codes = NULL,
  death.code = NULL,
  ...
)

Arguments

Details

For the Fine-Gray model, the score equations are:

\int (X - E(t)) Y_1(t) w(t) dM_1

summed over clusters and returned as iid.naive (multiplied by the inverse of the second derivative). Here,

w(t) = G(t) (I(T_i \wedge t < C_i)/G_c(T_i \wedge t))

,

E(t) = S_1(t)/S_0(t)

, and

S_j(t) = \sum X_i^j Y_{i1}(t) w_i(t) \exp(X_i^T \beta)

.

The full influence function (IID decomposition) for the beta coefficients includes a censoring term:

\int (X - E(t)) Y_1(t) w(t) dM_1 + \int q(s)/p(s) dM_c

which is returned as the iid component.

For the logistic link model, standard errors may be slightly underestimated because uncertainty from the recursive baseline estimation is not fully accounted for. For smaller datasets, it is recommended to use the prop.odds.subdist function from the timereg package (which uses more efficient weights) or to bootstrap the standard errors.

Value

An object of class "cifreg" (extending "phreg") containing:

Author(s)

Thomas Scheike

See Also

cifregFG, recreg, gofFG

Examples

## data with no ties
data(bmt,package="mets")
bmt$time <- bmt$time+runif(nrow(bmt))*0.01
bmt$id <- 1:nrow(bmt)

## logistic link  OR interpretation
or=cifreg(Event(time,cause)~tcell+platelet+age,data=bmt,cause=1)
summary(or)
par(mfrow=c(1,2))
plot(or)
nd <- data.frame(tcell=c(1,0),platelet=0,age=0)
por <- predict(or,nd)
plot(por)

## approximate standard errors 
por <-mets:::predict.phreg(or,nd)
plot(por,se=1)

## Fine-Gray model
fg=cifregFG(Event(time,cause)~tcell+platelet+age,data=bmt,cause=1)
summary(fg)
##fg=recreg(Event(time,cause)~tcell+platelet+age,data=bmt,cause=1,death.code=2)
##summary(fg)
plot(fg)
nd <- data.frame(tcell=c(1,0),platelet=0,age=0)
pfg <- predict(fg,nd,se=1)
plot(pfg,se=1)

## bt <- iidBaseline(fg,time=30)
## bt <- IIDrecreg(fg$cox.prep,fg,time=30)

## not run to avoid timing issues
## gofFG(Event(time,cause)~tcell+platelet+age,data=bmt,cause=1)

sfg <- cifregFG(Event(time,cause)~strata(tcell)+platelet+age,data=bmt,cause=1)
summary(sfg)
plot(sfg)

### predictions with CI based on iid decomposition of baseline and beta
### these are used in the predict function above
fg <- cifregFG(Event(time,cause)~tcell+platelet+age,data=bmt,cause=1)
Biid <- iidBaseline(fg,time=20)
pfg1 <- FGprediid(Biid,nd)
pfg1

Fine-Gray Cumulative Incidence Function Regression

Description

Convenience wrapper for cifreg that specifically fits the Fine-Gray model by setting propodds=NULL. This provides subdistribution hazard ratio interpretations.

Usage

cifregFG(formula, data, propodds = NULL, ...)

Arguments

Value

An object of class "cifreg" (extending "phreg") with propodds=NULL.

Author(s)

Thomas Scheike

See Also

cifreg

Finds subjects related to same cluster

Description

Finds subjects related to same cluster

Usage

cluster_index(
  clusters,
  index.type = FALSE,
  num = NULL,
  Rindex = 0,
  mat = NULL,
  return.all = FALSE,
  code.na = NA
)

Arguments

Author(s)

Klaus Holst, Thomas Scheike

References

Cluster indeces

See Also

familycluster_index familyclusterWithProbands.index

Examples

i<-c(1,1,2,2,1,3)
d<- cluster_index(i)
print(d)

type<-c("m","f","m","c","c","c")
d<- cluster_index(i,num=type,Rindex=1)
print(d)

Coarsen Cluster Identifiers

Description

Reduces the number of unique clusters by binning into quantile groups.

Usage

coarse_clust(clusters, max.clust = 100)

Arguments

Value

Integer vector of coarsened cluster indices (0-based).

Concordance Computes concordance and casewise concordance

Description

Concordance for Twins

Usage

concordanceCor(
  object,
  cif1,
  cif2 = NULL,
  messages = TRUE,
  model = NULL,
  coefs = NULL,
  ...
)

Arguments

Details

The concordance is the probability that both twins have experienced the event of interest and is defined as

cor(t) = P(T_1 \leq t, \epsilon_1 =1 , T_2 \leq t, \epsilon_2=1)

Similarly, the casewise concordance is

casewise(t) = \frac{cor(t)}{P(T_1 \leq t, \epsilon_1=1) }

that is the probability that twin "2" has the event given that twins "1" has.

Author(s)

Thomas Scheike

References

Estimating twin concordance for bivariate competing risks twin data Thomas H. Scheike, Klaus K. Holst and Jacob B. Hjelmborg, Statistics in Medicine 2014, 1193-1204

Estimating Twin Pair Concordance for Age of Onset. Thomas H. Scheike, Jacob V B Hjelmborg, Klaus K. Holst, 2015 in Behavior genetics DOI:10.1007/s10519-015-9729-3

Cross-odds-ratio, OR or RR risk regression for competing risks

Description

Fits a parametric model for the log-cross-odds-ratio for the predictive effect of for the cumulative incidence curves for T_1 experiencing cause i given that T_2 has experienced a cause k :

\log(COR(i|k)) = h(\theta,z_1,i,z_2,k,t)=_{default} \theta^T z =

with the log cross odds ratio being

COR(i|k) = \frac{O(T_1 \leq t,cause_1=i | T_2 \leq t,cause_2=k)}{ O(T_1 \leq t,cause_1=i)}

the conditional odds divided by the unconditional odds, with the odds being, respectively

O(T_1 \leq t,cause_1=i | T_2 \leq t,cause_1=k) = \frac{ P_x(T_1 \leq t,cause_1=i | T_2 \leq t,cause_2=k)}{ P_x((T_1 \leq t,cause_1=i)^c | T_2 \leq t,cause_2=k)}

and

O(T_1 \leq t,cause_1=i) = \frac{P_x(T_1 \leq t,cause_1=i )}{P_x((T_1 \leq t,cause_1=i)^c )}.

Here B^c is the complement event of B, P_x is the distribution given covariates (x are subject specific and z are cluster specific covariates), and h() is a function that is the simple identity \theta^T z by default.

Usage

cor_cif(
  cif,
  data,
  cause = NULL,
  times = NULL,
  cause1 = 1,
  cause2 = 1,
  cens.code = NULL,
  cens.model = "KM",
  Nit = 40,
  detail = 0,
  clusters = NULL,
  theta = NULL,
  theta.des = NULL,
  step = 1,
  sym = 0,
  weights = NULL,
  par.func = NULL,
  dpar.func = NULL,
  dimpar = NULL,
  score.method = "nlminb",
  same.cens = FALSE,
  censoring.weights = NULL,
  silent = 1,
  ...
)

Arguments

Details

The OR dependence measure is given by

OR(i,k) = \log ( \frac{O(T_1 \leq t,cause_1=i | T_2 \leq t,cause_2=k)}{ O(T_1 \leq t,cause_1=i) | T_2 \leq t,cause_2=k)}

This measure is numerically more stabile than the COR measure, and is symetric in i,k.

The RR dependence measure is given by

RR(i,k) = \log ( \frac{P(T_1 \leq t,cause_1=i , T_2 \leq t,cause_2=k)}{ P(T_1 \leq t,cause_1=i) P(T_2 \leq t,cause_2=k)}

This measure is numerically more stabile than the COR measure, and is symetric in i,k.

The model is fitted under symmetry (sym=1), i.e., such that it is assumed that T_1 and T_2 can be interchanged and leads to the same cross-odd-ratio (i.e. COR(i|k) = COR(k|i)), as would be expected for twins or without symmetry as might be the case with mothers and daughters (sym=0).

h() may be specified as an R-function of the parameters, see example below, but the default is that it is simply \theta^T z.

Value

returns an object of type 'cor'. With the following arguments:

Author(s)

Thomas Scheike

References

Cross odds ratio Modelling of dependence for Multivariate Competing Risks Data, Scheike and Sun (2012), Biostatistics.

A Semiparametric Random Effects Model for Multivariate Competing Risks Data, Scheike, Zhang, Sun, Jensen (2010), Biometrika.

Examples

 ## Reduce Ex.Timings
library("timereg")
data(multcif);
multcif$cause[multcif$cause==0] <- 2
zyg <- rep(rbinom(200,1,0.5),each=2)
theta.des <- model.matrix(~-1+factor(zyg))

times=seq(0.05,1,by=0.05) # to speed up computations use only these time-points
add <- timereg::comp.risk(Event(time,cause)~+1+cluster(id),data=multcif,cause=1,
               n.sim=0,times=times,model="fg",max.clust=NULL)
add2 <- timereg::comp.risk(Event(time,cause)~+1+cluster(id),data=multcif,cause=2,
               n.sim=0,times=times,model="fg",max.clust=NULL)

out1 <- cor_cif(add,data=multcif,cause1=1,cause2=1)
summary(out1)

out2 <- cor_cif(add,data=multcif,cause1=1,cause2=1,theta.des=theta.des)
summary(out2)

##out3 <- cor_cif(add,data=multcif,cause1=1,cause2=2,cif2=add2)
##summary(out3)
###########################################################
# investigating further models using parfunc and dparfunc
###########################################################
set.seed(100)
prt<-sim_nordic_random(2000,cordz=2,cormz=5)
prt$status <-prt$cause
table(prt$status)

times <- seq(40,100,by=10)
cifmod <- timereg::comp.risk(Event(time,cause)~+1+cluster(id),data=prt,
                    cause=1,n.sim=0,
                    times=times,conservative=1,max.clust=NULL,model="fg")
theta.des <- model.matrix(~-1+factor(zyg),data=prt)

parfunc <- function(par,t,pardes)
{
par <- pardes %*% c(par[1],par[2]) +
       pardes %*% c( par[3]*(t-60)/12,par[4]*(t-60)/12)
par
}
head(parfunc(c(0.1,1,0.1,1),50,theta.des))

dparfunc <- function(par,t,pardes)
{
dpar <- cbind(pardes, t(t(pardes) * c( (t-60)/12,(t-60)/12)) )
dpar
}
head(dparfunc(c(0.1,1,0.1,1),50,theta.des))

names(prt)
or1 <- or_cif(cifmod,data=prt,cause1=1,cause2=1,theta.des=theta.des,
              same.cens=TRUE,theta=c(0.6,1.1,0.1,0.1),
              par.func=parfunc,dpar.func=dparfunc,dimpar=4,
              score.method="nr",detail=1)
summary(or1)

## cor1 <- cor_cif(cifmod,data=prt,cause1=1,cause2=1,theta.des=theta.des,
##                 same.cens=TRUE,theta=c(0.5,1.0,0.1,0.1),
##                 par.func=parfunc,dpar.func=dparfunc,dimpar=4,
##                 control=list(trace=TRUE),detail=1)
## summary(cor1)

### piecewise contant OR model
gparfunc <- function(par,t,pardes)
{
	cuts <- c(0,80,90,120)
	grop <- diff(t<cuts)
paru  <- (pardes[,1]==1) * sum(grop*par[1:3]) +
    (pardes[,2]==1) * sum(grop*par[4:6])
paru
}

dgparfunc <- function(par,t,pardes)
{
	cuts <- c(0,80,90,120)
	grop <- diff(t<cuts)
par1 <- matrix(c(grop),nrow(pardes),length(grop),byrow=TRUE)
parmz <- par1* (pardes[,1]==1)
pardz <- (pardes[,2]==1) * par1
dpar <- cbind( parmz,pardz)
dpar
}
head(dgparfunc(rep(0.1,6),50,theta.des))
head(gparfunc(rep(0.1,6),50,theta.des))

or1g <- or_cif(cifmod,data=prt,cause1=1,cause2=1,
               theta.des=theta.des, same.cens=TRUE,
               par.func=gparfunc,dpar.func=dgparfunc,
               dimpar=6,score.method="nr",detail=1)
summary(or1g)
names(or1g)
head(or1g$theta.iid)

Compute cumulative event counts as time-dependent covariates

Description

For each subject and each row of data, counts the number of prior events of specified types in the recurrent event history. The resulting count columns can be used as time-dependent covariates in subsequent models, e.g. to capture event-history dependence in the recurrent event rate.

Usage

count_history(
  data,
  status = "status",
  id = "id",
  types = 1,
  names.count = "Count",
  lag = TRUE,
  multitype = FALSE,
  marks = NULL
)

Arguments

Details

When lag = TRUE (default), the count at each row reflects events that occurred strictly before the current time point (i.e. N(t-)), making it suitable as a left-continuous covariate in counting-process models.

Value

The input data frame data with one or more new integer columns appended. With multitype = FALSE, columns are named paste0(names.count, k) for each k in types; with multitype = TRUE, a single column named paste0(names.count, types[1]) is added. An internal bookkeeping column lbnr__id is also added.

Author(s)

Thomas Scheike

Examples

data(hfactioncpx12)
hf <- hfactioncpx12
dtable(hf, ~status)

## Separate counts for event types 1 and 2
rr <- count_history(hf, types = 1:2, id = "id", status = "status")
dtable(rr, ~"Count*" + status, level = 1)

Cumulative Odds Regression for Discrete Time Data

Description

A wrapper function for interval_logitsurv_discrete that simplifies the interface for discrete time-to-event data where the event time is observed exactly or as a factor level. It converts a factor response into an interval-censored format internally.

Usage

cumoddsreg(formula, data, ...)

Arguments

Value

An object of class "cumoddsreg" with the same structure as interval_logitsurv_discrete.

Author(s)

Thomas Scheike

See Also

interval_logitsurv_discrete

aggregating for for data frames

Description

aggregating for for data frames

Usage

daggregate(
  data,
  y = NULL,
  x = NULL,
  subset,
  ...,
  fun = "summary",
  regex = mets.options()$regex,
  missing = FALSE,
  remove.empty = FALSE,
  matrix = FALSE,
  silent = FALSE,
  na.action = na.pass,
  convert = NULL
)

Arguments

Examples

data("sTRACE")
daggregate(iris, "^.e.al", x="Species", fun=cor, regex=TRUE)
daggregate(iris, Sepal.Length+Petal.Length ~Species, fun=summary)
daggregate(iris, log(Sepal.Length)+I(Petal.Length>1.5) ~ Species,
                 fun=summary)
daggregate(iris, "*Length*", x="Species", fun=head)
daggregate(iris, "^.e.al", x="Species", fun=tail, regex=TRUE)
daggregate(sTRACE, status~ diabetes, fun=table)
daggregate(sTRACE, status~ diabetes+sex, fun=table)
daggregate(sTRACE, status + diabetes+sex ~ vf+I(wmi>1.4), fun=table)
daggregate(iris, "^.e.al", x="Species",regex=TRUE)
dlist(iris,Petal.Length+Sepal.Length ~ Species |Petal.Length>1.3 & Sepal.Length>5,
            n=list(1:3,-(3:1)))
daggregate(iris, I(Sepal.Length>7)~Species | I(Petal.Length>1.5))
daggregate(iris, I(Sepal.Length>7)~Species | I(Petal.Length>1.5),
                 fun=table)

dsum(iris, .~Species, matrix=TRUE, missing=TRUE)

par(mfrow=c(1,2))
data(iris)
drename(iris) <- ~.
daggregate(iris,'sepal*'~species|species!="virginica",fun=plot)
daggregate(iris,'sepal*'~I(as.numeric(species))|I(as.numeric(species))!=1,fun=summary)

dnumeric(iris) <- ~species
daggregate(iris,'sepal*'~species.n|species.n!=1,fun=summary)

Calculate summary statistics grouped by

Description

Calculate summary statistics grouped by variable

Usage

dby(
  data,
  INPUT,
  ...,
  ID = NULL,
  ORDER = NULL,
  SUBSET = NULL,
  SORT = 0,
  COMBINE = !REDUCE,
  NOCHECK = FALSE,
  ARGS = NULL,
  NAMES,
  COLUMN = FALSE,
  REDUCE = FALSE,
  REGEX = mets.options()$regex,
  ALL = TRUE
)

Arguments

Details

Calculate summary statistics grouped by

dby2 for column-wise calculations

Author(s)

Klaus K. Holst and Thomas Scheike

Examples

n <- 4
k <- c(3,rbinom(n-1,3,0.5)+1)
N <- sum(k)
d <- data.frame(y=rnorm(N),x=rnorm(N),
   id=rep(seq(n),k),num=unlist(sapply(k,seq))
)
d2 <- d[sample(nrow(d)),]

dby(d2, y~id, mean)
dby(d2, y~id + order(num), cumsum)

dby(d,y ~ id + order(num), dlag)
dby(d,y ~ id + order(num), dlag, ARGS=list(k=1:2))
dby(d,y ~ id + order(num), dlag, ARGS=list(k=1:2), NAMES=c("l1","l2"))

dby(d, y~id + order(num), mean=mean, csum=cumsum, n=length)
dby(d2, y~id + order(num), a=cumsum, b=mean, N=length,
  l1=function(x) c(NA,x)[-length(x)]
)

dby(d, y~id + order(num), nn=seq_along, n=length)
dby(d, y~id + order(num), nn=seq_along, n=length)

d <- d[,1:4]
dby(d, x<0) <- list(z=mean)
d <- dby(d, is.na(z), z=1)

f <- function(x) apply(x,1,min)
dby(d, y+x~id, min=f)

dby(d,y+x~id+order(num), function(x) x)

f <- function(x) { cbind(cumsum(x[,1]),cumsum(x[,2]))/sum(x)}
dby(d, y+x~id, f)

## column-wise
a <- d
dby2(a, mean, median, REGEX=TRUE) <- '^[y|x]'~id
a
## wildcards 
dby2(a,'y*'+'x*'~id,mean) 


## subset
dby(d, x<0) <- list(z=NA)
d
dby(d, y~id|x>-1, v=mean,z=1)
dby(d, y+x~id|x>-1, mean, median, COLUMN=TRUE)

dby2(d, y+x~id|x>0, mean, REDUCE=TRUE)

dby(d,y~id|x<0,mean,ALL=FALSE)

a <- iris
a <- dby(a,y=1)
dby(a,Species=="versicolor") <- list(y=2)

summary, tables, and correlations for data frames

Description

summary, tables, and correlations for data frames

Usage

dcor(data, y = NULL, x = NULL, use = "pairwise.complete.obs", ...)

Arguments

Author(s)

Klaus K. Holst and Thomas Scheike

Examples

data("sTRACE",package="timereg")
dt<- sTRACE
dt$time2 <- dt$time^2
dt$wmi2 <- dt$wmi^2
head(dt)

dcor(dt)

dcor(dt,~time+wmi)
dcor(dt,~time+wmi,~vf+chf)
dcor(dt,time+wmi~vf+chf)

dcor(dt,c("time*","wmi*"),~vf+chf)

Cutting, sorting, rm (removing), rename for data frames

Description

Cut variables, if breaks are given these are used, otherwise cuts into using group size given by probs, or equispace groups on range. Default is equally sized groups if possible

Usage

dcut(
  data,
  y = NULL,
  x = NULL,
  breaks = 4,
  probs = NULL,
  equi = FALSE,
  regex = mets.options()$regex,
  sep = NULL,
  na.rm = TRUE,
  labels = NULL,
  all = FALSE,
  ...
)

Arguments

Author(s)

Klaus K. Holst and Thomas Scheike

Examples

data("sTRACE",package="timereg")
sTRACE$age2 <- sTRACE$age^2
sTRACE$age3 <- sTRACE$age^3

mm <- dcut(sTRACE,~age+wmi)
head(mm)

mm <- dcut(sTRACE,catage4+wmi4~age+wmi)
head(mm)

mm <- dcut(sTRACE,~age+wmi,breaks=c(2,4))
head(mm)

mm <- dcut(sTRACE,c("age","wmi"))
head(mm)

mm <- dcut(sTRACE,~.)
head(mm)

mm <- dcut(sTRACE,c("age","wmi"),breaks=c(2,4))
head(mm)

gx <- dcut(sTRACE$age)
head(gx)


## Removes all cuts variables with these names wildcards
mm1 <- drm(mm,c("*.2","*.4"))
head(mm1)

## wildcards, for age, age2, age4 and wmi
head(dcut(mm,c("a*","?m*")))

## with direct asignment
drm(mm) <- c("*.2","*.4")
head(mm)

dcut(mm) <- c("age","*m*")
dcut(mm) <- ageg1+wmig1~age+wmi
head(mm)

############################
## renaming
############################

head(mm)
drename(mm, ~Age+Wmi) <- c("wmi","age")
head(mm)
mm1 <- mm

## all names to lower
drename(mm1) <- ~.
head(mm1)

## A* to lower
mm2 <-  drename(mm,c("A*","W*"))
head(mm2)
drename(mm) <- "A*"
head(mm)

dd <- data.frame(A_1=1:2,B_1=1:2)
funn <- function(x) gsub("_",".",x)
drename(dd) <- ~.
drename(dd,fun=funn) <- ~.
names(dd)

Dermal ridges data (families)

Description

Data on dermal ridge counts in left and right hand in (nuclear) families

Format

Data on 50 families with ridge counts in left and right hand for moter, father and each child. Family id in 'family' and gender and child number in 'sex' and 'child'.

Source

Sarah B. Holt (1952). Genetics of dermal ridges: bilateral asymmetry in finger ridge-counts. Annals of Eugenics 17 (1), pp.211–231. DOI: 10.1111/j.1469-1809.1952.tb02513.x

Examples

data(dermalridges)
fast.reshape(dermalridges,id="family",varying=c("child.left","child.right","sex"))

Dermal ridges data (monozygotic twins)

Description

Data on dermal ridge counts in left and right hand in (nuclear) families

Format

Data on dermal ridge counts (left and right hand) in 18 monozygotic twin pairs.

Source

Sarah B. Holt (1952). Genetics of dermal ridges: bilateral asymmetry in finger ridge-counts. Annals of Eugenics 17 (1), pp.211–231. DOI: 10.1111/j.1469-1809.1952.tb02513.x

Examples

data(dermalridgesMZ)
fast.reshape(dermalridgesMZ,id="id",varying=c("left","right"))

The Diabetic Retinopathy Data

Description

The data was colleceted to test a laser treatment for delaying blindness in patients with dibetic retinopathy. The subset of 197 patiens given in Huster et al. (1989) is used.

Format

This data frame contains the following columns:

id

a numeric vector. Patient code.

agedx

a numeric vector. Age of patient at diagnosis.

time

a numeric vector. Survival time: time to blindness or censoring.

status

a numeric vector code. Survival status. 1: blindness, 0: censored.

trteye

a numeric vector code. Random eye selected for treatment. 1: left eye 2: right eye.

treat

a numeric vector. 1: treatment 0: untreated.

adult

a numeric vector code. 1: younger than 20, 2: older than 20.

Source

Huster W.J. and Brookmeyer, R. and Self. S. (1989) MOdelling paired survival data with covariates, Biometrics 45, 145-56.

Examples

data(diabetes)
names(diabetes)

Split a data set and run function

Description

Split a data set and run function

Usage

divide_conquer(func = NULL, data, size, splits, id = NULL, ...)

Arguments

Author(s)

Thomas Scheike, Klaus K. Holst

Examples

## avoid dependency on timereg
## library(timereg)
## data(TRACE)
## res <- divide_conquer(prop.odds,TRACE,
## 	     formula=Event(time,status==9)~chf+vf+age,n.sim=0,size=200)

Lag operator

Description

Lag operator

Usage

dlag(data, x, k = 1, combine = TRUE, simplify = TRUE, names, ...)

Arguments

Examples

d <- data.frame(y=1:10,x=c(10:1))
dlag(d,k=1:2)
dlag(d,~x,k=0:1)
dlag(d$x,k=1)
dlag(d$x,k=-1:2, names=letters[1:4])

list, head, print, tail

Description

listing for data frames

Usage

dprint(data, y = NULL, n = 0, ..., x = NULL)

Arguments

Author(s)

Klaus K. Holst and Thomas Scheike

Examples

m <- lava::lvm(letters)
d <- lava::sim(m, 20)

dlist(d,~a+b+c)
dlist(d,~a+b+c|a<0 & b>0)
## listing all : 
dlist(d,~a+b+c|a<0 & b>0,n=0)
dlist(d,a+b+c~I(d>0)|a<0 & b>0)
dlist(d,.~I(d>0)|a<0 & b>0)
dlist(d,~a+b+c|a<0 & b>0, nlast=0)
dlist(d,~a+b+c|a<0 & b>0, nfirst=3, nlast=3)
dlist(d,~a+b+c|a<0 & b>0, 1:5)
dlist(d,~a+b+c|a<0 & b>0, -(5:1))
dlist(d,~a+b+c|a<0 & b>0, list(1:5,50:55,-(5:1)))
dprint(d,a+b+c ~ I(d>0) |a<0 & b>0, list(1:5,50:55,-(5:1)))

Regression for data frames with dutility call

Description

Regression for data frames with dutility call

Usage

dreg(
  data,
  y,
  x = NULL,
  z = NULL,
  x.oneatatime = TRUE,
  x.base.names = NULL,
  z.arg = c("clever", "base", "group", "condition"),
  fun. = lm,
  summary. = summary,
  regex = FALSE,
  convert = NULL,
  doSummary = TRUE,
  special = NULL,
  equal = TRUE,
  test = 1,
  ...
)

Arguments

Author(s)

Klaus K. Holst, Thomas Scheike

Examples

##'
data(iris)
dat <- iris
drename(dat) <- ~.
names(dat)
set.seed(1)
dat$time <- runif(nrow(dat))
dat$time1 <- runif(nrow(dat))
dat$status <- rbinom(nrow(dat),1,0.5)
dat$S1 <- with(dat, Surv(time,status))
dat$S2 <- with(dat, Surv(time1,status))
dat$id <- 1:nrow(dat)

mm <- dreg(dat, "*.length"~"*.width"|I(species=="setosa" & status==1))
mm <- dreg(dat, "*.length"~"*.width"|species+status)
mm <- dreg(dat, "*.length"~"*.width"|species)
mm <- dreg(dat, "*.length"~"*.width"|species+status,z.arg="group")

 ## Reduce Ex.Timings
y <- "S*"~"*.width"
xs <- dreg(dat, y, fun.=phreg)
## xs <- dreg(dat, y, fun.=survdiff)

y <- "S*"~"*.width"
xs <- dreg(dat, y, x.oneatatime=FALSE, fun.=phreg)

## under condition
y <- S1~"*.width"|I(species=="setosa" & sepal.width>3)
xs <- dreg(dat, y, z.arg="condition", fun.=phreg)
xs <- dreg(dat, y, fun.=phreg)

## under condition
y <- S1~"*.width"|species=="setosa"
xs <- dreg(dat, y, z.arg="condition", fun.=phreg)
xs <- dreg(dat, y, fun.=phreg)

## with baseline  after |
y <- S1~"*.width"|sepal.length
xs <- dreg(dat, y, fun.=phreg)

## by group by species, not working
y <- S1~"*.width"|species
ss <- split(dat, paste(dat$species, dat$status))

xs <- dreg(dat, y, fun.=phreg)

## species as base, species is factor so assumes that this is grouping
y <- S1~"*.width"|species
xs <- dreg(dat, y, z.arg="base", fun.=phreg)

##  background var after | and then one of x's at at time
y <- S1~"*.width"|status+"sepal*"
xs <- dreg(dat, y, fun.=phreg)

##  background var after | and then one of x's at at time
##y <- S1~"*.width"|status+"sepal*"
##xs <- dreg(dat, y, x.oneatatime=FALSE, fun.=phreg)
##xs <- dreg(dat, y, fun.=phreg)

##  background var after | and then one of x's at at time
##y <- S1~"*.width"+factor(species)
##xs <- dreg(dat, y, fun.=phreg)
##xs <- dreg(dat, y, fun.=phreg, x.oneatatime=FALSE)

y <- S1~"*.width"|factor(species)
xs <- dreg(dat, y, z.arg="base", fun.=phreg)

y <- S1~"*.width"|cluster(id)+factor(species)
xs <- dreg(dat, y, z.arg="base", fun.=phreg)
xs <- dreg(dat, y, z.arg="base", fun.=survival::coxph)

## under condition with groups
y <- S1~"*.width"|I(sepal.length>4)
xs <- dreg(subset(dat, species=="setosa"), y,z.arg="group",fun.=phreg)

## under condition with groups
y <- S1~"*.width"+I(log(sepal.length))|I(sepal.length>4)
xs <- dreg(subset(dat, species=="setosa"), y,z.arg="group",fun.=phreg)

y <- S1~"*.width"+I(dcut(sepal.length))|I(sepal.length>4)
xs <- dreg(subset(dat,species=="setosa"), y,z.arg="group",fun.=phreg)

ff <- function(formula,data,...) {
 ss <- survfit(formula,data,...)
 kmplot(ss,...)
 return(ss)
}

if (interactive()) {
dcut(dat) <- ~"*.width"
y <- S1~"*.4"|I(sepal.length>4)
par(mfrow=c(1, 2))
xs <- dreg(dat, y, fun.=ff)
}

relev levels for data frames

Description

levels shows levels for variables in data frame, relevel relevels a factor in data.frame

Usage

drelevel(
  data,
  y = NULL,
  x = NULL,
  ref = NULL,
  newlevels = NULL,
  regex = mets.options()$regex,
  sep = NULL,
  overwrite = FALSE,
  ...
)

Arguments

Author(s)

Klaus K. Holst and Thomas Scheike

Examples

data(mena)
dstr(mena)
dfactor(mena)  <- ~twinnum
dnumeric(mena) <- ~twinnum.f

dstr(mena)

mena2 <- drelevel(mena,"cohort",ref="(1980,1982]")
mena2 <- drelevel(mena,~cohort,ref="(1980,1982]")
mena2 <- drelevel(mena,cohortII~cohort,ref="(1980,1982]")
dlevels(mena)
dlevels(mena2)
drelevel(mena,ref="(1975,1977]")  <-  ~cohort
drelevel(mena,ref="(1980,1982]")  <-  ~cohort
dlevels(mena,"coh*")
dtable(mena,"coh*",level=1)

### level 1 of zyg as baseline for new variable
drelevel(mena,ref=1) <- ~zyg
drelevel(mena,ref=c("DZ","[1973,1975]")) <- ~ zyg+cohort
drelevel(mena,ref=c("DZ","[1973,1975]")) <- zygdz+cohort.early~ zyg+cohort
### level 2 of zyg and cohort as baseline for new variables
drelevel(mena,ref=2) <- ~ zyg+cohort
dlevels(mena)

##################### combining factor levels with newlevels argument

dcut(mena,labels=c("I","II","III","IV")) <- cat4~agemena
dlevels(drelevel(mena,~cat4,newlevels=1:3))
dlevels(drelevel(mena,ncat4~cat4,newlevels=3:2))
drelevel(mena,newlevels=3:2) <- ncat4~cat4
dlevels(mena)

dlevels(drelevel(mena,nca4~cat4,newlevels=list(c(1,4),2:3)))

drelevel(mena,newlevels=list(c(1,4),2:3)) <- nca4..2 ~ cat4
dlevels(mena)

drelevel(mena,newlevels=list("I-III"=c("I","II","III"),"IV"="IV")) <- nca4..3 ~ cat4
dlevels(mena)

drelevel(mena,newlevels=list("I-III"=c("I","II","III"))) <- nca4..4 ~ cat4
dlevels(mena)

drelevel(mena,newlevels=list(group1=c("I","II","III"))) <- nca4..5 ~ cat4
dlevels(mena)

drelevel(mena,newlevels=list(g1=c("I","II","III"),g2="IV")) <- nca4..6 ~ cat4
dlevels(mena)

Remove Special Terms from a Formula

Description

Removes terms such as strata() or cluster() from a formula/terms object.

Usage

drop.specials(x, components, ...)

Arguments

Value

The updated formula with an attribute "variables" containing the extracted variable names from removed specials.

Sort data frame

Description

Sort data according to columns in data frame

Usage

dsort(data, x, ..., decreasing = FALSE, return.order = FALSE)

Arguments

Value

data.frame

Examples

data(data="hubble",package="lava")
dsort(hubble, "sigma")
dsort(hubble, hubble$sigma,"v")
dsort(hubble,~sigma+v)
dsort(hubble,~sigma-v)

## with direct asignment
dsort(hubble) <- ~sigma-v

Simple linear spline

Description

Constructs simple linear spline on a data frame using the formula syntax of dutils that is adds (x-cuti)* (x>cuti) to the data-set for each knot of the spline. The full spline is thus given by x and spline variables added to the data-set.

Usage

dspline(
  data,
  y = NULL,
  x = NULL,
  breaks = 4,
  probs = NULL,
  equi = FALSE,
  regex = mets.options()$regex,
  sep = NULL,
  na.rm = TRUE,
  labels = NULL,
  all = FALSE,
  ...
)

Arguments

Author(s)

Thomas Scheike

Examples

data(TRACE)
TRACE <- dspline(TRACE,~wmi,breaks=c(1,1.3,1.7))
cca <- survival::coxph(Surv(time,status==9)~age+vf+chf+wmi,data=TRACE)
cca2 <- survival::coxph(Surv(time,status==9)~age+wmi+vf+chf+
                           wmi.spline1+wmi.spline2+wmi.spline3,data=TRACE)
anova(cca,cca2)

nd <- data.frame(age=50,vf=0,chf=0,wmi=seq(0.4,3,by=0.01))
nd <- dspline(nd,~wmi,breaks=c(1,1.3,1.7))
pl <- predict(cca2,newdata=nd)
plot(nd$wmi,pl,type="l")

tables for data frames

Description

tables for data frames

Usage

dtable(
  data,
  y = NULL,
  x = NULL,
  ...,
  level = -1,
  response = NULL,
  flat = TRUE,
  total = FALSE,
  prop = FALSE,
  summary = NULL
)

Arguments

Author(s)

Klaus K. Holst and Thomas Scheike

Examples

data("sTRACE",package="timereg")

dtable(sTRACE,~status)
dtable(sTRACE,~status+vf)
dtable(sTRACE,~status+vf,level=1)
dtable(sTRACE,~status+vf,~chf+diabetes)

dtable(sTRACE,c("*f*","status"),~diabetes)
dtable(sTRACE,c("*f*","status"),~diabetes, level=2)
dtable(sTRACE,c("*f*","status"),level=1)

dtable(sTRACE,~"*f*"+status,level=1)
dtable(sTRACE,~"*f*"+status+I(wmi>1.4)|age>60,level=2)
dtable(sTRACE,"*f*"+status~I(wmi>0.5)|age>60,level=1)
dtable(sTRACE,status~dcut(age))

dtable(sTRACE,~status+vf+sex|age>60)
dtable(sTRACE,status+vf+sex~+1|age>60, level=2)
dtable(sTRACE,.~status+vf+sex|age>60,level=1)
dtable(sTRACE,status+vf+sex~diabetes|age>60)
dtable(sTRACE,status+vf+sex~diabetes|age>60, flat=FALSE)

dtable(sTRACE,status+vf+sex~diabetes|age>60, level=1)
dtable(sTRACE,status+vf+sex~diabetes|age>60, level=2)

dtable(sTRACE,status+vf+sex~diabetes|age>60, level=2, prop=1, total=TRUE)
dtable(sTRACE,status+vf+sex~diabetes|age>60, level=2, prop=2, total=TRUE)
dtable(sTRACE,status+vf+sex~diabetes|age>60, level=2, prop=1:2, summary=summary)

Transform that allows condition

Description

Defines new variables under condition for data frame

Usage

dtransform(data, ...)

Arguments

Examples

data(mena)

xx <- dtransform(mena,ll=log(agemena)+twinnum)

xx <- dtransform(mena,ll=log(agemena)+twinnum,agemena<15)
xx <- dtransform(xx  ,ll=100+agemena,ll2=1000,agemena>15)
dsummary(xx,ll+ll2~I(agemena>15))

event_split (SurvSplit).

Description

Constructs start stop formulation of event time data after a variable in the data.set. Similar to SurvSplit of the survival package but can also split after random time given in data frame.

Usage

event_split(
  data,
  time = "time",
  status = "status",
  cuts = "cuts",
  name.id = "id",
  name.start = "start",
  cens.code = 0,
  order.id = TRUE,
  time.group = FALSE
)

Arguments

Author(s)

Thomas Scheike

Examples

set.seed(1)
d <- data.frame(event=round(5*runif(5),2),start=1:5,time=2*1:5,
		status=rbinom(5,1,0.5),x=1:5)
d

d0 <- event_split(d,cuts="event",name.start=0)
d0

dd <- event_split(d,cuts="event")
dd
ddd <- event_split(dd,cuts=3.5)
ddd
event_split(ddd,cuts=5.5)

### successive cutting for many values 
dd <- d
for  (cuts in seq(2,3,by=0.3)) dd <- event_split(dd,cuts=cuts)
dd

###########################################################################
### same but for situation with multiple events along the time-axis
###########################################################################
d <- data.frame(event1=1:5+runif(5)*0.5,start=1:5,time=2*1:5,
		status=rbinom(5,1,0.5),x=1:5,start0=0)
d$event2 <- d$event1+0.2
d$event2[4:5] <- NA 
d

d0 <- event_split(d,cuts="event1",name.start="start",time="time",status="status")
d0
###
d00 <- event_split(d0,cuts="event2",name.start="start",time="time",status="status")
d00

Event split with two time-scales, time and gaptime

Description

Cuts time for two time-scales, as event.split

Usage

event_split2(
  data,
  time = "time",
  status = "status",
  entry = "start",
  cuts = "cuts",
  name.id = "id",
  gaptime = NULL,
  gaptime.entry = NULL,
  cuttime = c("time", "gaptime"),
  cens.code = 0,
  order.id = TRUE
)

Arguments

Author(s)

Thomas Scheike

Examples

rr  <- data.frame(time=c(500,1000),start=c(0,500),status=c(1,1),id=c(1,1))
rr$gaptime <-  rr$time-rr$start
rr$gapstart <- 0

rr1 <- event_split2(rr,cuts=600,cuttime="time",   gaptime="gaptime",gaptime.entry="gapstart")
rr2 <- event_split2(rr1,cuts=100,cuttime="gaptime",gaptime="gaptime",gaptime.entry="gapstart")

dlist(rr1,start-time+status+gapstart+gaptime~id)
dlist(rr2,start-time+status+gapstart+gaptime~id)

Extract survival estimates from lifetable analysis

Description

Summary for survival analyses via the 'lifetable' function

Usage

eventpois(
  object,
  ...,
  timevar,
  time,
  int.len,
  confint = FALSE,
  level = 0.95,
  individual = FALSE,
  length.out = 25
)

Arguments

Details

Summary for survival analyses via the 'lifetable' function

Author(s)

Klaus K. Holst

Extend Cumulative Hazard Functions to Common Time Range

Description

Extends a collection of cumulative hazard functions so that they all cover the same time range. Cumulative hazards that end before the maximum observed time are extrapolated linearly using either a user-supplied rate or the average hazard rate estimated from the existing cumulative hazard.

Usage

extendCums(cumA, cumB, extend = NULL)

Arguments

Value

A list of cumulative hazard matrices extended to the common maximum time.

See Also

rchaz, rcrisk, mets-simulation

Finds all pairs within a cluster (famly) with the proband (case/control)

Description

second column of pairs are the probands and the first column the related subjects

Usage

familyclusterWithProbands_index(
  clusters,
  probands,
  index.type = FALSE,
  num = NULL,
  Rindex = 1
)

Arguments

Author(s)

Klaus Holst, Thomas Scheike

References

Cluster indeces

See Also

familycluster_index cluster_index

Examples

i<-c(1,1,2,2,1,3)
p<-c(1,0,0,1,0,1)
d<- familyclusterWithProbands_index(i,p)
print(d)

Finds all pairs within a cluster (family)

Description

Finds all pairs within a cluster (family)

Usage

familycluster_index(clusters, index.type = FALSE, num = NULL, Rindex = 1)

Arguments

Author(s)

Klaus Holst, Thomas Scheike

References

Cluster indeces

See Also

cluster_index familyclusterWithProbands_index

Examples

i<-c(1,1,2,2,1,3)
d<- familycluster_index(i)
print(d)

Fast approximation

Description

Fast approximation

Usage

fast.approx(
  time,
  new.time,
  equal = FALSE,
  type = c("nearest", "right", "left"),
  sorted = FALSE,
  ...
)

Arguments

Author(s)

Klaus K. Holst

Examples

id <- c(1,1,2,2,7,7,10,10)
fast.approx(unique(id),id)

t <- 0:6
n <- c(-1,0,0.1,0.9,1,1.1,1.2,6,6.5)
fast.approx(t,n,type="left")

Fast Cluster Index Conversion

Description

Converts a cluster variable to consecutive numeric indices using compiled code.

Usage

fast.cluster(x, ...)

Arguments

Value

Integer vector of 0-based cluster indices.

Fast pattern

Description

Fast pattern

Usage

fast.pattern(x, y, categories = 2, ...)

Arguments

Author(s)

Klaus K. Holst

Examples

X <- matrix(rbinom(100,1,0.5),ncol=4)
fast.pattern(X)

X <- matrix(rbinom(100,3,0.5),ncol=4)
fast.pattern(X,categories=4)

Fast reshape

Description

Fast reshape/tranpose of data

Usage

fast.reshape(
  data,
  varying,
  id,
  num,
  sep = "",
  keep,
  idname = "id",
  numname = "num",
  factor = FALSE,
  idcombine = TRUE,
  labelnum = FALSE,
  labels,
  regex = mets.options()$regex,
  dropid = FALSE,
  ...
)

Arguments

Author(s)

Thomas Scheike, Klaus K. Holst

Examples

m <- lava::lvm(c(y1,y2,y3,y4)~x)
d <- lava::sim(m,5)
d
fast.reshape(d,"y")
fast.reshape(fast.reshape(d,"y"),id="id")

##### From wide-format
(dd <- fast.reshape(d,"y"))
## Same with explicit setting new id and number variable/column names
## and seperator "" (default) and dropping x
fast.reshape(d,"y",idname="a",timevar="b",sep="",keep=c())
## Same with 'reshape' list-syntax
fast.reshape(d,list(c("y1","y2","y3","y4")),labelnum=TRUE)

##### From long-format
fast.reshape(dd,id="id")
## Restrict set up within-cluster varying variables
fast.reshape(dd,"y",id="id")
fast.reshape(dd,"y",id="id",keep="x",sep=".")

#####
x <- data.frame(id=c(5,5,6,6,7),y=1:5,x=1:5,tv=c(1,2,2,1,2))
x
(xw <- fast.reshape(x,id="id"))
(xl <- fast.reshape(xw,c("y","x"),idname="id2",keep=c()))
(xl <- fast.reshape(xw,c("y","x","tv")))
(xw2 <- fast.reshape(xl,id="id",num="num"))
fast.reshape(xw2,c("y","x"),idname="id")

### more generally:
### varying=list(c("ym","yf","yb1","yb2"), c("zm","zf","zb1","zb2"))
### varying=list(c("ym","yf","yb1","yb2")))

##### Family cluster example
d <- mets:::sim_BinFam(3)
d
fast.reshape(d,var="y")
fast.reshape(d,varying=list(c("ym","yf","yb1","yb2")))

d <- lava::sim(lava::lvm(~y1+y2+ya),10)
d
(dd <- fast.reshape(d,"y"))
fast.reshape(d,"y",labelnum=TRUE)
fast.reshape(dd,id="id",num="num")
fast.reshape(dd,id="id",num="num",labelnum=TRUE)
fast.reshape(d,c(a="y"),labelnum=TRUE) ## New column name


##### Unbalanced data
m <- lava::lvm(c(y1,y2,y3,y4)~ x+z1+z3+z5)
d <- lava::sim(m,3)
d
fast.reshape(d,c("y","z"))

##### not-varying syntax:
fast.reshape(d,-c("x"))

##### Automatically define varying variables from trailing digits
fast.reshape(d)

##### Prostate cancer example
data(prt)
head(prtw <- fast.reshape(prt,"cancer",id="id"))
ftable(cancer1~cancer2,data=prtw)
rm(prtw)

Fast Reshape from Long to Wide Format

Description

Reshapes clustered long-format data to wide format efficiently using compiled code.

Usage

faster.reshape(data, clusters, index.type = FALSE, num = NULL, Rindex = 1)

Arguments

Value

A wide-format matrix.

Generate Random Fold Indices for Cross-Validation

Description

Splits n observations into random folds for cross-validation.

Usage

folds(n, folds = 10)

Arguments

Value

A list of integer vectors, each containing indices for one fold.

Force Same Censoring Within Clusters

Description

Enforces the same censoring time within clusters (pairs) by censoring both subjects at the minimum censoring time when censoring occurs first.

Usage

force.same.cens(
  data,
  id = "id",
  time = "time",
  cause = "cause",
  entrytime = NULL,
  cens.code = 0
)

force_same_cens(
  data,
  id = "id",
  time = "time",
  cause = "cause",
  entrytime = NULL,
  cens.code = 0
)

Arguments

Value

A data.frame with enforced same censoring.

IPTW GLM, Inverse Probabibilty of Treatment Weighted GLM

Description

Fits GLM model with treatment weights

w(A)= \sum_a I(A=a)/P(A=a|X)

, computes standard errors via influence functions that are returned as the IID argument. Propensity scores are fitted using either logistic regression (glm) or the multinomial model (mlogit) when more than two categories for treatment. The treatment needs to be a factor and is identified on the rhs of the "treat.model".

Usage

glm_IPTW(
  formula,
  data,
  treat.model = NULL,
  family = binomial(),
  id = NULL,
  weights = NULL,
  estpr = 1,
  pi0 = 0.5,
  ...
)

Arguments

Details

Also works with cluster argument.

Author(s)

Thomas Scheike

Goodness-of-Fit for Cox PH Regression (Proportionality)

Description

Performs cumulative score process residual tests for the proportional hazards (PH) assumption in Cox regression. The test statistics are based on the cumulative score process:

U(t) = \int_0^t (X_i - E(t) ) d \hat M_i(s)

where \hat M_i(s) are the martingale residuals.

Usage

## S3 method for class 'phreg'
gof(object, n.sim = 1000, silent = 1, robust = NULL, ...)

Arguments

Details

P-values are computed using the Lin, Wei, and Ying (1993) resampling method, which simulates the asymptotic distribution of the supremum of the score process under the null hypothesis of proportional hazards.

The function supports two types of simulation:

• Standard: Uses dN_i (counting process increments) for simulation.

• Robust: Uses \hat M_i(t) (martingale residuals) adjusted for clustering if a cluster() term is present in the model.

Value

An object of class "gof.phreg" containing:

Author(s)

Thomas Scheike and Klaus K. Holst

References

Lin, D. Y., Wei, L. J., & Ying, Z. (1993). Checking the Cox model with cumulative sums of martingale-based residuals. Biometrika, 80(3), 557-572.

See Also

gofM_phreg, gofZ_phreg

Examples

data(sTRACE)

m1 <- phreg(Surv(time,status==9)~vf+chf+diabetes, data=sTRACE) 
gg <- gof(m1)
gg
par(mfrow=c(1,3))
plot(gg)

m1 <- phreg(Surv(time,status==9)~strata(vf)+chf+diabetes, data=sTRACE) 
gg <- gof(m1)

## Robust simulations with cluster
sTRACE$id <- 1:500
m1 <- phreg(Surv(time,status==9)~vf+chf+diabetes+cluster(id), data=sTRACE) 
gg <- gof(m1)
gg

Goodness-of-Fit for Cox Covariates (Model Matrix)

Description

Tests the functional form of covariates in a Cox PH model by computing cumulative residuals against a user-specified model matrix. This helps detect non-linear effects or time-varying coefficients (interaction with time).

Usage

gofM_phreg(
  formula,
  data,
  offset = NULL,
  weights = NULL,
  modelmatrix = NULL,
  n.sim = 1000,
  silent = 1,
  ...
)

Arguments

Details

The test statistic is:

U(t) = \int_0^t K^T d \hat M

where K is the model matrix (e.g., a set of basis functions for a continuous covariate) and \hat M are the martingale residuals.

P-values are based on the Lin, Wei, and Ying (1993) resampling method. The plot shows whether the residuals are consistent with the model across the range of the covariate.

Value

An object of class "gof.phreg" containing:

Author(s)

Thomas Scheike and Klaus K. Holst

References

Lin, D. Y., Wei, L. J., & Ying, Z. (1993). Checking the Cox model with cumulative sums of martingale-based residuals. Biometrika, 80(3), 557-572.

See Also

gof.phreg, gofZ_phreg, cumContr

Examples

data(TRACE)
set.seed(1)
TRACEsam <- blocksample(TRACE, idvar="id", replace=FALSE, 100)
dcut(TRACEsam) <- ~. 
mm <- model.matrix(~-1+factor(wmicat.4), data=TRACEsam)
m1 <- gofM_phreg(Surv(time,status==9)~vf+chf+wmi, data=TRACEsam, modelmatrix=mm)
summary(m1)
if (interactive()) {
par(mfrow=c(2,2))
plot(m1)
}

## Cumulative sums in covariates via design matrix
mm <- mets:::cumContr(TRACEsam$wmi, breaks=10, equi=TRUE)
m1 <- gofM_phreg(Surv(time,status==9)~strata(vf)+chf+wmi, data=TRACEsam,
         modelmatrix=mm, silent=0)
summary(m1)

Goodness-of-Fit for Cox Covariates (Linearity)

Description

Tests the functional form of continuous covariates in a Cox PH model to check for linearity. It computes cumulative residuals evaluated at a grid of covariate values z.

Usage

gofZ_phreg(
  formula,
  data,
  vars = NULL,
  offset = NULL,
  weights = NULL,
  breaks = 50,
  equi = FALSE,
  n.sim = 1000,
  silent = 1,
  ...
)

Arguments

Details

The test statistic is:

U(z, \tau) = \int_0^\tau K(z)^T d \hat M

where K(z) is a design matrix based on indicator functions I(Z_i \leq z_l) for a grid of points z_l.

The p-value is valid but depends on the chosen grid. As the number of break points increases, this test converges to the original Lin, Wei, and Ying test for linearity.

Value

An object of class "gof.phreg" with type "Zmodelmatrix" containing:

Author(s)

Thomas Scheike and Klaus K. Holst

See Also

gofM_phreg, cumContr

Examples

data(TRACE)
set.seed(1)
TRACEsam <- blocksample(TRACE, idvar="id", replace=FALSE, 100)

## Test linearity of continuous covariates
 ## Reduce Ex.Timings
m1 <- gofZ_phreg(Surv(time,status==9)~strata(vf)+chf+wmi+age, data=TRACEsam)
summary(m1) 
plot(m1, type="z")

Create Group Contingency Table from Clustered Data

Description

Creates a contingency table by group from paired/clustered data, optionally combining lower and upper triangles.

Usage

grouptable(
  data,
  id,
  group,
  var,
  lower = TRUE,
  labels,
  order,
  group.labels,
  group.order,
  combine = " & ",
  ...
)

Arguments

Value

A table or list of tables.

haplo fun data

Description

haplo fun data

Format

hapfreqs : haplo frequencies haploX: covariates and response for haplo survival discrete survival ghaplos: haplo-types for subjects of haploX data

Source

Estimated data

Discrete Time-to-Event Haplotype Analysis

Description

Performs cycle-specific logistic regression to estimate haplotype effects on discrete time-to-event data, accounting for phase ambiguity. Given observed genotypes G and unobserved haplotypes H, the method integrates (mixes out) over the possible haplotype configurations using the conditional probabilities P(H|G).

Usage

haplo_surv_discrete(
  X = NULL,
  y = "y",
  time.name = "time",
  Haplos = NULL,
  id = "id",
  desnames = NULL,
  designfunc = NULL,
  beta = NULL,
  no.opt = FALSE,
  method = "NR",
  stderr = TRUE,
  designMatrix = NULL,
  response = NULL,
  idhap = NULL,
  design.only = FALSE,
  covnames = NULL,
  fam = binomial,
  weights = NULL,
  offsets = NULL,
  idhapweights = NULL,
  ...
)

Arguments

Details

The survival function is computed by averaging over the possible haplotypes:

S(t|x,G) = E[ S(t|x,H) | G ] = \sum_{h \in G} P(h|G) S(t|x,h)

The discrete hazard function is modeled using logistic regression:

\text{logit}(P(T=t | T \geq t, x, h)) = \alpha_t + x(h) \beta

where \alpha_t are time-specific intercepts (baseline hazards), x(h) is the regression design constructed from covariates and haplotypes h=(h_1, h_2), and \beta are the regression coefficients.

The likelihood is maximized numerically. Standard errors are computed assuming that P(H|G) is known (i.e., ignoring the uncertainty in haplotype estimation).

The design matrix over possible haplotypes is constructed by merging the covariate data X with the haplotype probabilities Haplos and applying a user-defined designfunc.

Value

An object of class "haplosurvd" containing:

Author(s)

Thomas Scheike

References

Scheike, T. H. (2024). Discrete time survival analysis with haplotype effects. mets package documentation.

Examples

## Some haplotypes of interest
types <- c("DCGCGCTCACG","DTCCGCTGACG","ITCAGTTGACG","ITCCGCTGAGG")

## Some haplotype frequencies for simulations 
data(haplo)
hapfreqs <- haplo$hapfreqs 

www <- which(hapfreqs$haplotype %in% types)
hapfreqs$freq[www]

baseline <- hapfreqs$haplotype[9]
baseline

## Design function: indicator for presence of any 'types' haplotype
designftypes <- function(x, sm=0) {
  hap1 <- x[1]
  hap2 <- x[2]
  if (sm == 0) y <- 1 * ((hap1 == types) | (hap2 == types))
  if (sm == 1) y <- 1 * (hap1 == types) + 1 * (hap2 == types)
  return(y)
}

tcoef <- c(-1.93110204, -0.47531630, -0.04118204, -1.57872602, -0.22176426, -0.13836416,
           0.88830288, 0.60756224, 0.39802821, 0.32706859)

ghaplos <- haplo$ghaplos
haploX  <- haplo$haploX

haploX$time <- haploX$times
Xdes <- model.matrix(~ factor(time), haploX)
colnames(Xdes) <- paste("X", 1:ncol(Xdes), sep="")
X <- dkeep(haploX, ~ id + y + time)
X <- cbind(X, Xdes)
Haplos <- dkeep(ghaplos, ~ id + "haplo*" + p)
desnames <- paste("X", 1:6, sep="")   # Six X's related to 6 cycles 
out <- haplo_surv_discrete(X=X, y="y", time.name="time",
         Haplos=Haplos, desnames=desnames, designfunc=designftypes) 
names(out$coef) <- c(desnames, types)
out$coef
summary(out)

hfaction, subset of block randomized study HF-ACtion from WA package

Description

Data from HF-action trial slightly modified from WA package, consisting of 741 nonischemic patients with baseline cardiopulmonary test duration less than or equal to 12 minutes.

Format

Randomized study status : 1-event, 2-death, 0-censoring treatment : 1/0

Source

WA package, Connor et al. 2009

Examples

data(hfactioncpx12)

Influence Functions or IID Decomposition of Baseline

Description

Computes the influence functions for the baseline cumulative hazard (and optionally regression coefficients) for phreg, recreg, or cifregFG objects.

Usage

iidBaseline(
  object,
  time = NULL,
  ft = NULL,
  fixbeta = NULL,
  beta.iid = NULL,
  tminus = FALSE,
  ...
)

Arguments

Details

The decomposition is based on the formula:

\sum_i \int_0^t \frac{f(s)}{S_0(s)} dM_{ki}(s) - P(t) \beta_k

where k denotes the stratum and i the cluster.

Value

An object of class "iidBaseline" containing:

Author(s)

Thomas Scheike

See Also

phreg, recreg, cifreg

Inverse Laplace Transform Helper

Description

Computes the inverse Laplace transform for gamma frailty simulation.

Usage

ilap(theta, t)

Arguments

Value

Numeric value of the inverse Laplace transform.

Discrete Time-to-Event Analysis with Interval Censoring

Description

Fits a cumulative odds model for discrete time-to-event data, handling interval censoring where the event time is known only to lie within an interval (t_l, t_r]. The model assumes:

\text{logit}(P(T \leq t | x)) = \log(G(t)) + x \beta

where G(t) is the baseline cumulative odds function and \beta are the regression coefficients. This is equivalent to:

P(T \leq t | x) = \frac{G(t) \exp(x \beta)}{1 + G(t) \exp(x \beta)}

Usage

interval_logitsurv_discrete(
  formula,
  data,
  beta = NULL,
  no.opt = FALSE,
  method = "NR",
  stderr = TRUE,
  weights = NULL,
  offsets = NULL,
  exp.link = 1,
  increment = 1,
  ...
)

Arguments

Details

The baseline G(t) is parameterized as the cumulative sum of exponentials (G(t) = \sum \exp(\alpha)), ensuring positivity. The regression coefficients describe the log-odds of the event occurring by time t.

The likelihood is maximized over the observed intervals:

L = \prod_i [ P(T_i > t_{il} | x_i) - P(T_i > t_{ir} | x_i) ]

where t_{il} and t_{ir} are the left and right endpoints of the interval for subject i. Right-censored intervals have t_{ir} = \infty.

Value

An object of class "cumoddsreg" containing:

Author(s)

Thomas Scheike

References

Scheike, T. H. (2024). Discrete time survival analysis with interval censoring. mets package documentation.

See Also

cumoddsreg, predictlogitSurvd, simlogitSurvd

Examples

data(ttpd) 
dtable(ttpd,~entry+time2)

out <- interval_logitsurv_discrete(Interval(entry,time2)~X1+X2+X3+X4,ttpd)
summary(out)
head(iid(out)) 

pred <- predictlogitSurvd(out,se=FALSE)
plotSurvd(pred)

ttpd <- dfactor(ttpd,fentry~entry)
out <- cumoddsreg(fentry~X1+X2+X3+X4,ttpd)
summary(out)

Inverse Probability of Censoring Weights

Description

Internal function. Calculates Inverse Probability of Censoring Weights (IPCW) and adds them to a data.frame

Usage

ipw(
  formula,
  data,
  cluster,
  same.cens = FALSE,
  obs.only = FALSE,
  weight.name = "w",
  trunc.prob = FALSE,
  weight.name2 = "wt",
  indi.weight = "pr",
  cens.model = "aalen",
  pairs = FALSE,
  theta.formula = ~1,
  ...
)

Arguments

Author(s)

Klaus K. Holst

Examples

## Not run: 
data("prt",package="mets")
prtw <- ipw(Surv(time,status==0)~country, data=prt[sample(nrow(prt),5000),],
            cluster="id",weight.name="w")
plot(0,type="n",xlim=range(prtw$time),ylim=c(0,1),xlab="Age",ylab="Probability")
count <- 0
for (l in unique(prtw$country)) {
    count <- count+1
    prtw <- prtw[order(prtw$time),]
    with(subset(prtw,country==l),
         lines(time,w,col=count,lwd=2))
}
legend("topright",legend=unique(prtw$country),col=1:4,pch=-1,lty=1)

## End(Not run)

Inverse Probability of Censoring Weights

Description

Internal function. Calculates Inverse Probability of Censoring and Truncation Weights and adds them to a data.frame

Usage

ipw2(
  data,
  times = NULL,
  entrytime = NULL,
  time = "time",
  cause = "cause",
  same.cens = FALSE,
  cluster = NULL,
  pairs = FALSE,
  strata = NULL,
  obs.only = TRUE,
  cens.formula = NULL,
  cens.code = 0,
  pair.cweight = "pcw",
  pair.tweight = "ptw",
  pair.weight = "weights",
  cname = "cweights",
  tname = "tweights",
  weight.name = "indi.weights",
  prec.factor = 100,
  ...
)

Arguments

Author(s)

Thomas Scheike

Examples

library("timereg")
set.seed(1)
d <- sim_nordic_random(5000,delayed=TRUE,ptrunc=0.7,
      cordz=0.5,cormz=2,lam0=0.3,country=FALSE)
d$strata <- as.numeric(d$country)+(d$zyg=="MZ")*4
times <- seq(60,100,by=10)
## c1 <- timereg::comp.risk(Event(time,cause)~1+cluster(id),data=d,cause=1,
## 	model="fg",times=times,max.clust=NULL,n.sim=0)
## mm=model.matrix(~-1+zyg,data=d)
## out1<-random_cif(c1,data=d,cause1=1,cause2=1,same.cens=TRUE,theta.des=mm)
## summary(out1)
## pc1 <- predict(c1,X=1,se=0)
## plot(pc1)
## 
## dl <- d[!d$truncated,]
## dl <- ipw2(dl,cluster="id",same.cens=TRUE,time="time",entrytime="entry",cause="cause",
##            strata="strata",prec.factor=100)
## cl <- timereg::comp.risk(Event(time,cause)~+1+
## 		cluster(id),
##  		data=dl,cause=1,model="fg",
## 		weights=dl$indi.weights,cens.weights=rep(1,nrow(dl)),
##           times=times,max.clust=NULL,n.sim=0)
## pcl <- predict(cl,X=1,se=0)
## lines(pcl$time,pcl$P1,col=2)
## mm=model.matrix(~-1+factor(zyg),data=dl)
## out2<-random_cif(cl,data=dl,cause1=1,cause2=1,theta.des=mm,
##                  weights=dl$weights,censoring.weights=rep(1,nrow(dl)))
## summary(out2)

Extract Event (Jump) Times

Description

Extracts unique event times from survival data, optionally restricting to concordant pairs or specific causes.

Usage

jumptimes(
  time,
  status = TRUE,
  id,
  cause,
  sample,
  sample.all = TRUE,
  strata = NULL,
  num = NULL,
  ...
)

Arguments

Value

Sorted vector of event times.

Kaplan-Meier with Robust Standard Errors

Description

Computes Kaplan-Meier estimates with robust standard errors. Robust variance is the default and is obtained from the predict call.

Usage

km(formula, data = data, km = TRUE, ...)

Arguments

Value

An object of class "km" (extends "predictphreg") containing:

Author(s)

Thomas Scheike

Examples

data(sTRACE)
sTRACE$cluster <- sample(1:100, 500, replace = TRUE)
out1 <- km(Surv(time, status == 9) ~ strata(vf, chf), data = sTRACE)
out2 <- km(Surv(time, status == 9) ~ strata(vf, chf) + cluster(cluster), data = sTRACE)

summary(out1, times = 1:3)
summary(out2, times = 1:3)

par(mfrow = c(1, 2))
plot(out1, se = TRUE)
plot(out2, se = TRUE)

Life-course plot

Description

Life-course plot for event life data with recurrent events

Usage

lifecourse(
  formula,
  data,
  id = "id",
  group = NULL,
  type = "l",
  lty = 1,
  col = 1:10,
  alpha = 0.3,
  lwd = 1,
  recurrent.col = NULL,
  recurrent.lty = NULL,
  legend = NULL,
  pchlegend = NULL,
  by = NULL,
  status.legend = NULL,
  place.sl = "bottomright",
  xlab = "Time",
  ylab = "",
  add = FALSE,
  ...
)

Arguments

Author(s)

Thomas Scheike, Klaus K. Holst

Examples

data = data.frame(id=c(1,1,1,2,2),start=c(0,1,2,3,4),slut=c(1,2,4,4,7),
                  type=c(1,2,3,2,3),status=c(0,1,2,1,2),group=c(1,1,1,2,2))
ll = lifecourse(Event(start,slut,status)~id,data,id="id")
ll = lifecourse(Event(start,slut,status)~id,data,id="id",recurrent.col="type")

ll = lifecourse(Event(start,slut,status)~id,data,id="id",group=~group,col=1:2)
op <- par(mfrow=c(1,2))
ll = lifecourse(Event(start,slut,status)~id,data,id="id",by=~group)
par(op)
legends=c("censored","pregnant","married")
ll = lifecourse(Event(start,slut,status)~id,data,id="id",group=~group,col=1:2,status.legend=legends)

Life table

Description

Create simple life table

Usage

## S3 method for class 'matrix'
lifetable(x, strata = list(), breaks = c(),
   weights=NULL, confint = FALSE, ...)

 ## S3 method for class 'formula'
lifetable(x, data=parent.frame(), breaks = c(),
   weights=NULL, confint = FALSE, ...)

Arguments

Author(s)

Klaus K. Holst

Examples

library(timereg)
data(TRACE)

d <- with(TRACE,lifetable(Surv(time,status==9)~sex+vf,breaks=c(0,0.2,0.5,8.5)))
lava::estimate(glm(events ~ offset(log(atrisk))+factor(int.end)*vf + sex*vf,
            data=d,poisson))

Proportional Odds Survival Model

Description

Fits a semiparametric proportional odds model where:

\mbox{logit}(S(t|x)) = \log(\Lambda(t)) + x \beta

Thus, covariate effects represent the odds ratio (OR) of survival.

Usage

logitSurv(formula, data, offset = NULL, weights = NULL, ...)

Arguments

Details

This is equivalent to using a hazards model:

Z \lambda(t) \exp(x \beta)

where Z is gamma distributed with mean and variance 1.

Value

An object of class "phreg" with propodds=1.

Author(s)

Thomas Scheike

References

Eriksson, Frank, Li, Jianing, Scheike, Thomas, and Zhang, Mei-Jie (2015). "The proportional odds cumulative incidence model for competing risks." Biometrics, 71(3), 687–695.

Examples

data(TRACE)
dcut(TRACE) <- ~.
out1 <- logitSurv(Surv(time, status == 9) ~ vf + chf + strata(wmicat.4), data = TRACE)
summary(out1)
gof(out1)
plot(out1)

Mediation analysis in survival context

Description

Mediation analysis in survival context with robust standard errors taking the weights into account via influence function computations. Mediator and exposure must be factors. This is based on numerical derivative wrt parameters for weighting. See vignette for more examples.

Usage

mediatorSurv(
  survmodel,
  weightmodel,
  data = data,
  wdata = wdata,
  id = "id",
  silent = TRUE,
  ...
)

Arguments

Author(s)

Thomas Scheike

Examples

n <- 400
dat <- kumarsimRCT(n,rho1=0.5,rho2=0.5,rct=2,censpar=c(0,0,0,0),
          beta = c(-0.67, 0.59, 0.55, 0.25, 0.98, 0.18, 0.45, 0.31),
    treatmodel = c(-0.18, 0.56, 0.56, 0.54),restrict=1)
dfactor(dat) <- dnr.f~dnr
dfactor(dat) <- gp.f~gp
drename(dat) <- ttt24~"ttt24*"
dat$id <- 1:n
dat$ftime <- 1

weightmodel <- fit <- glm(gp.f~dnr.f+preauto+ttt24,data=dat,family=binomial)
wdata <- medweight(fit,data=dat)

### fitting models with and without mediator
aaMss2 <- binreg(Event(time,status)~gp+dnr+preauto+ttt24+cluster(id),data=dat,time=50,cause=2)
aaMss22 <- binreg(Event(time,status)~dnr+preauto+ttt24+cluster(id),data=dat,time=50,cause=2)

### estimating direct and indirect effects (under strong strong assumptions) 
aaMss <- binreg(Event(time,status)~dnr.f0+dnr.f1+preauto+ttt24+cluster(id),
                data=wdata,time=50,weights=wdata$weights,cause=2)
## to compute standard errors , requires numDeriv
ll <- mediatorSurv(aaMss,fit,data=dat,wdata=wdata)
summary(ll)
## not run bootstrap (to save time)
## bll <- BootmediatorSurv(aaMss,fit,data=dat,k.boot=500)

Computes mediation weights

Description

Computes mediation weights for either binary or multinomial mediators. The important part is that the influence functions can be obtained to compute standard errors.

Usage

medweight(
  fit,
  data = data,
  var = NULL,
  name.weight = "weights",
  id.name = "id",
  ...
)

Arguments

Author(s)

Thomas Scheike

The Melanoma Survival Data

Description

The melanoma data frame has 205 rows and 7 columns. It contains data relating to survival of patients after operation for malignant melanoma collected at Odense University Hospital by K.T. Drzewiecki.

Format

This data frame contains the following columns:

no

a numeric vector. Patient code.

status

a numeric vector code. Survival status. 1: dead from melanoma, 2: alive, 3: dead from other cause.

days

a numeric vector. Survival time.

ulc

a numeric vector code. Ulceration, 1: present, 0: absent.

thick

a numeric vector. Tumour thickness (1/100 mm).

sex

a numeric vector code. 0: female, 1: male.

Source

Andersen, P.K., Borgan O, Gill R.D., Keiding N. (1993), Statistical Models Based on Counting Processes, Springer-Verlag.

Drzewiecki, K.T., Ladefoged, C., and Christensen, H.E. (1980), Biopsy and prognosis for cutaneous malignant melanoma in clinical stage I. Scand. J. Plast. Reconstru. Surg. 14, 141-144.

Examples

data(melanoma)
names(melanoma)

Menarche data set

Description

Menarche data set

Source

Simulated data

Simulation Helper Functions

Description

Internal simulation functions used for generating data from various survival, competing risks, frailty, and twin/family models. These functions are primarily intended for use in examples and testing.

Survival and Competing Risks

simrchaz

Simulate from cumulative hazard via inverse CDF.

simul_cifs

Simulate competing risks data from cumulative incidence functions.

simlogitSurvd

Simulate survival data using logistic model.

kumarsim

Simulate competing risks data (Kumaraswamy-type).

kumarsimRCT

Simulate RCT competing risks data (Kumaraswamy-type).

Clayton-Oakes and Frailty Models

sim_ClaytonOakes_family_ace

Simulate family data from Clayton-Oakes ACE model.

sim_ClaytonOakes_twin_ace

Simulate twin data from Clayton-Oakes ACE model.

sim_Compete_twin_ace

Simulate twin competing risks data with ACE frailty.

sim_Compete_simple

Simulate competing risks data with shared frailty.

sim_Frailty_simple

Simulate survival data with shared frailty.

sim_SurvFam

Simulate family survival data with shared frailty.

Binomial Twin/Family Models

sim_BinPlack

Simulate paired binary data (Plackett model).

sim_BinFam

Simulate binary family data.

sim_BinFam2

Simulate binary family data (alternative parameterization).

sim_binClaytonOakes_family_ace

Simulate binary family data (Clayton-Oakes ACE).

sim_binClaytonOakes_twin_ace

Simulate binary twin data (Clayton-Oakes ACE).

sim_binClaytonOakes_pairs

Simulate binary paired data (Clayton-Oakes).

sim_bptwin

Simulate from a biprobit twin model.

Nordic Twin Studies

sim_nordictwin

Simulate Nordic twin registry data.

sim_nordic_random

Simulate Nordic twin data with random effects.

corsim_prostate_random

Simulate correlated prostate cancer data with random effects.

Set global options for mets

Description

Extract and set global parameters of mets.

Usage

mets.options(...)

Arguments

Details

• regex: If TRUE character vectors will be interpreted as regular expressions (dby, dcut, ...)

• silent: Set to FALSE to disable various output messages

Value

list of parameters

Examples

## Not run: 
mets.options(regex=TRUE)

## End(Not run)

Migraine data

Description

Migraine data

Multinomial Regression Based on phreg

Description

Fits a multinomial regression model for a categorical outcome with K levels:

P_i = \frac{ \exp( X^\top \beta_i ) }{ \sum_{j=1}^K \exp( X^\top \beta_j ) }

for i=1, \dots, K, where \beta_1 = 0 (baseline category).

Usage

mlogit(formula, data, offset = NULL, weights = NULL, fix.X = FALSE, ...)

Arguments

Details

This ensures that \sum_j P_j = 1. The model is fitted using the phreg function by expanding the data into a long format with strata for each category.

The coefficients represent the log-Relative-Risk (log-RR) relative to the baseline group (the first level of the factor, which can be reset using relevel).

Standard errors are computed based on the sandwich form:

D U^{-1} \left( \sum U_i^2 \right) D U^{-1}

where U is the score vector and D is the derivative matrix.

Influence functions (possibly robust) can be obtained via the iid() function. The response variable should be a factor.

Can also fit a cumulative odds model as a special case of interval_logitsurv_discrete.

Value

An object of class "mlogit" (extending "phreg") containing:

Author(s)

Thomas Scheike

Examples

data(bmt)
bmt$id <- sample(200,408,replace=TRUE)
dfactor(bmt) <- cause1f~cause
drelevel(bmt,ref=3) <- cause3f~cause
dlevels(bmt)

mreg <- mlogit(cause1f~+1+cluster(id),bmt)
summary(mreg)
head(iid(mreg))
dim(iid(mreg))

mreg <- mlogit(cause1f~tcell+platelet,bmt)
summary(mreg)
head(iid(mreg))
dim(iid(mreg))

mreg3 <- mlogit(cause3f~tcell+platelet,bmt)
summary(mreg3)

## inverse information standard errors 
lava::estimate(coef=mreg3$coef,vcov=mreg3$II)

## predictions based on seen response or not 
## all probabilities
head(predict(mreg,response=FALSE))
head(predict(mreg))
## using newdata 
newdata <- data.frame(tcell=c(1,1,1),platelet=c(0,1,1),cause1f=c("2","2","0"))
## only probability of seen response 
predict(mreg,newdata)
## without response
predict(mreg,newdata,response=FALSE)
## given indexx of P(Y=j)
predict(mreg,newdata,Y=c(1,2,3))
##  reponse not given 
newdata <- data.frame(tcell=c(1,1,1),platelet=c(0,1,1))
predict(mreg,newdata)

Multivariate Cumulative Incidence Function example data set

Description

Multivariate Cumulative Incidence Function example data set

Source

Simulated data

np data set

Description

np data set

Source

Simulated data

Concordance Probability from Twostage Model

Description

Computes concordance probability (joint probability of both subjects experiencing the event) given dependence parameters and random-effect variance structures from a twostage model.

Computes concordance probabilities for twin ACE/ADE models from a binomial twostage object.

Functions for constructing random-effects design matrices for twin and family models. These designs specify the genetic (A), dominance (D), common environment (C), and unique environment (E) variance components.

Usage

p11_binomial_twostage_RV(
  theta,
  rv1,
  rv2,
  p1,
  p2,
  pardes,
  ags = NULL,
  link = 0,
  i = 1,
  j = 1
)

concordanceTwostage(
  theta,
  p,
  rv1,
  rv2,
  theta.des,
  additive.gamma.sum = NULL,
  link = 0,
  var.par = 0,
  ...
)

concordanceTwinACE(
  object,
  rv1 = NULL,
  rv2 = NULL,
  xmarg = NULL,
  type = "ace",
  ...
)

kendall_ClaytonOakes_twin_ace(parg, parc, K = 10000, test = 0)

kendall.ClaytonOakes.twin.ace(x, y, ...)

kendall_normal_twin_ace(parg, parc, K = 10000)

ascertained_pairs(pairs, data, cr.models, bin = FALSE)

twin.polygen.design(x, ...)

ace_family_design(
  data,
  id = "id",
  member = "type",
  mother = "mother",
  father = "father",
  child = "child",
  child1 = "child",
  type = "ace",
  ...
)

make_pairwise_design(pairs, kinship, type = "ace")

Arguments

Details

twin_polygen_design creates a polygenic random-effects design for twin pairs, distinguishing MZ and DZ twins.

twin.polygen.design is an alias for twin_polygen_design.

ace_family_design creates designs for nuclear families (mother, father, children).

make_pairwise_design creates pairwise random-effects designs for arbitrary kinship structures.

concordanceTwostage computes concordance probabilities from a twostage model.

concordanceTwinACE computes concordance from a twin ACE model.

kendall_ClaytonOakes_twin_ace and kendall_normal_twin_ace compute Kendall's tau for Clayton-Oakes and normal-frailty twin ACE models respectively.

ascertained_pairs identifies ascertained (affected) pairs in clustered survival data.

p11_binomial_twostage_RV computes the joint probability P(T1<=t, T2<=t) for the additive gamma binary random effects model.

Value

A list of concordance tables, one per pair, each containing pmat (2x2 probability matrix), casewise (casewise concordance), and marg (marginal probabilities).

A list of concordance tables per zygosity group.

A list with components:

Author(s)

Klaus K. Holst, Thomas Scheike

Fast Cox Proportional Hazards Regression

Description

Fits a Cox proportional hazards model using a fast C++ backend. Robust variance (sandwich estimator) is the default variance estimate in the summary.

Usage

phreg(formula, data, offset = NULL, weights = NULL, ...)

Arguments

Details

The influence functions (IID decomposition) follow the numerical order of the given cluster variable. Ordering the results by '$id' will align the IID terms with the original dataset order.

Value

An object of class "phreg" containing:

Author(s)

Klaus K. Holst, Thomas Scheike

See Also

plot.phreg, predict.phreg, robust_phreg

Examples

data(TRACE)
dcut(TRACE) <- ~.
# Fit model with clustering
out1 <- phreg(Surv(time, status == 9) ~ wmi + age + strata(vf, chf) + cluster(id), data = TRACE)
summary(out1)

# Plotting baselines
par(mfrow = c(1, 2))
plot(out1)

# Computing robust variance for baseline
rob1 <- robust_phreg(out1)
plot(rob1, se = TRUE, robust = TRUE)

# IID decomposition for regression parameters
head(iid(out1))

# IID decomposition for baseline at a specific time-point
Aiiid <- iid(out1, time = 30)
head(Aiiid)

# Combined IID decomposition (beta and baseline)
dd <- iidBaseline(out1, time = 30)
head(dd$beta.iid)
head(dd$base.iid)

# Stratified model
outs <- phreg(Surv(time, status == 9) ~ strata(vf, wmicat.4) + cluster(id), data = TRACE)
summary(outs)
par(mfrow = c(1, 2))
plot(outs)

IPTW Cox Regression (Inverse Probability of Treatment Weighted)

Description

Fits a Cox model with treatment weights

w(A) = \sum_a I(A=a)/\pi(a|X)

, where

\pi(a|X) = P(A=a|X)

.

Usage

phreg_IPTW(
  formula,
  data,
  treat.model = NULL,
  treat.var = NULL,
  weights = NULL,
  estpr = 1,
  pi0 = 0.5,
  se.cluster = NULL,
  ...
)

Arguments

Details

Standard errors are computed via influence functions that are returned as the IID argument. Propensity scores are fitted using either logistic regression (glm) or the multinomial model (mlogit) when there are more than two treatment categories.

The treatment variable must be a factor and is identified on the RHS of the treat.model. Recurrent events can be considered with a start-stop structure, requiring cluster(id). Robust standard errors are computed in all cases.

Time-dependent propensity score weights can be computed when treat.var is used. This weight be 1 at the time of first (A_0) and 2nd treatment (A_1), then uses weights

w_0(A_0) * w_1(A_1)^{t>T_r}

where

T_r

is time of 2nd randomization. The weights are constructed using a glm or mlogit model based on the data where treat.var=1. The propensity score can be constructed for any number of treatments in a similar manner.

Value

An object of class "phreg" with additional IPTW components:

Author(s)

Thomas Scheike

Examples

data <- mets:::simLT(0.7, 100, beta = 0.3, betac = 0, ce = 1, betao = 0.3)
dfactor(data) <- Z.f ~ Z
out <- phreg_IPTW(Surv(time, status) ~ Z.f, data = data, treat.model = Z.f ~ X)
summary(out)

Lu-Tsiatis More Efficient Log-Rank for Randomized Studies with Baseline Covariates

Description

Efficient implementation of the Lu-Tsiatis improvement using baseline covariates, extended to competing risks and recurrent events. The results are almost equivalent to the speffSurv function of the speff2trial package in the survival case. A dynamic censoring augmentation regression is also computed to gain additional efficiency from the censoring augmentation. The function handles two-stage randomizations and recurrent events (start,stop) with cluster structure.

Usage

phreg_rct(
  formula,
  data,
  cause = 1,
  cens.code = 0,
  typesR = c("R0", "R1", "R01"),
  typesC = c("C", "dynC"),
  weights = NULL,
  augmentR0 = NULL,
  augmentR1 = NULL,
  augmentC = NULL,
  treat.model = ~+1,
  RCT = TRUE,
  treat.var = NULL,
  km = TRUE,
  level = 0.95,
  cens.model = NULL,
  estpr = 1,
  pi0 = 0.5,
  base.augment = FALSE,
  return.augmentR0 = FALSE,
  mlogit = FALSE,
  ...
)

Arguments

Details

The method improves the efficiency of the log-rank test by utilizing auxiliary baseline covariates to reduce variance, particularly useful in randomized clinical trials (RCTs) where covariate adjustment can increase power.

Value

An object of class "phreg_rct" containing: