cor.cif: 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,
  ...
)

Value

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

theta: estimate of proportional odds parameters of model.
var.theta: variance for gamma.
hess: the derivative of the used score.
score: scores at final stage.
score: scores at final stage.
theta.iid: matrix of iid decomposition of parametric effects.

Arguments

cif: a model object from the timereg::comp.risk function with the marginal cumulative incidence of cause1, i.e., the event of interest, and whose odds the comparision is compared to the conditional odds given cause2
data: a data.frame with the variables.
cause: specifies the causes related to the death times, the value cens.code is the censoring value. When missing it comes from marginal cif.
times: time-vector that specifies the times used for the estimating euqations for the cross-odds-ratio estimation.
cause1: specificies the cause considered.
cause2: specificies the cause that is conditioned on.
cens.code: specificies the code for the censoring if NULL then uses the one from the marginal cif model.
cens.model: specified which model to use for the ICPW, KM is Kaplan-Meier alternatively it may be "cox"
Nit: number of iterations for Newton-Raphson algorithm.
detail: if 0 no details are printed during iterations, if 1 details are given.
clusters: specifies the cluster structure.
theta: specifies starting values for the cross-odds-ratio parameters of the model.
theta.des: specifies a regression design for the cross-odds-ratio parameters.
step: specifies the step size for the Newton-Raphson algorithm.
sym: specifies if symmetry is used in the model.
weights: weights for estimating equations.
par.func: parfunc
dpar.func: dparfunc
dimpar: dimpar
score.method: "nlminb", can also use "nr".
same.cens: if true then censoring within clusters are assumed to be the same variable, default is independent censoring.
censoring.weights: these probabilities are used for the bivariate censoring dist.
silent: 1 to suppress output about convergence related issues.
...: Not used.

Author

Thomas Scheike

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$.

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

Run this code

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
###########################################################
 ## Reduce Ex.Timings
set.seed(100)
prt<-simnordic.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

Run the code above in your browser using DataLab