cor.cif | R Documentation |
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.
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, ... )
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. |
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.
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. |
Thomas Scheike
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.
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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.