cor.cif: Cross-odds-ratio, OR or RR risk regression for competing...

View source: R/cor.R

cor.cifR Documentation

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(θ,z_1,i,z_2,k,t)=_{default} θ^T z =

with the log cross odds ratio being

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

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

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

and

O(T_1 ≤q t,cause_1=i) = \frac{P_x(T_1 ≤q t,cause_1=i )}{P_x((T_1 ≤q 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 θ^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

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.

Details

The OR dependence measure is given by

OR(i,k) = \log ( \frac{O(T_1 ≤q t,cause_1=i | T_2 ≤q t,cause_2=k)}{ O(T_1 ≤q t,cause_1=i) | T_2 ≤q 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 ≤q t,cause_1=i , T_2 ≤q t,cause_2=k)}{ P(T_1 ≤q t,cause_1=i) P(T_2 ≤q 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 θ^T z.

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.

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

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<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)


mets documentation built on Jan. 17, 2023, 5:12 p.m.