Nothing
R/cor.R
In mets: Analysis of Multivariate Event Times
Defines functions
predictPairPlack
plack.cif2
plack.cif
concordanceCor
concordance.cor
print.cor
coef.cor
summary.cor
print.summary.cor
Grandom.cif
random.cif
or.cif
rr.cif
cor.cif
dep.cif
Documented in
concordance.cor
concordanceCor
cor.cif
Grandom.cif
or.cif
plack.cif
plack.cif2
predictPairPlack
random.cif
rr.cif
summary.cor
dep.cif<-function(cif,data,cause=NULL,model="OR",cif2=NULL,times=NULL,
cause1=1,cause2=1,cens.code=NULL,cens.model="KM",Nit=40,tol=1e-6,detail=0,
clusters=NULL,theta=NULL,theta.des=NULL,step=1,sym=1,weights=NULL,
same.cens=FALSE,censoring.weights=NULL,silent=1,entry=NULL,estimator=1,
trunkp=1,admin.cens=NULL,control=list(),par.func=NULL,dpar.func=NULL,dimpar=NULL,
score.method="nlminb",random.design=NULL,var.link=0,...)
{ ## {{{
## {{{ set up data and design
multi<-0; dscore=1; stab.cens<-FALSE; entry.call<-entry; inverse<-var.link
notaylor<-1; flex.func<-0;
## extract design and time and cause from cif object
time <- cif$response[,"exit"]
if (is.null(cause)) cause <- cif$response[,"cause"] ## attr(cif,"cause");
cause <- as.numeric(cause)
if (is.null(cens.code)) cens.code <- attr(cif,"cens.code")
delta<-(cause!=cens.code)
if (length(cause)!=length(time)) stop("cause and time not of same length\n");
formula <- attr(cif,"Formula")
ldata <- timereg::aalen.des2(formula(delete.response(terms(formula))),data=data,model="aalen")
X<-ldata$X; Z<-ldata$Z;
antpers<-nrow(X);
if (is.null(Z)==TRUE) {npar<-TRUE; semi<-0;} else {Z<-as.matrix(Z); npar<-FALSE; semi<-1;}
if (npar==TRUE) {Z<-matrix(0,antpers,1); pg<-1; fixed<-0;} else {fixed<-1;pg<-ncol(Z);}
ng<-antpers;px<-ncol(X);
if (is.null(times)) times<-cif$cum[,1];
est<-as.matrix(cif$cum);
if (semi==1) gamma<-cif$gamma else gamma<-0;
ntimes<-length(times);
if ((cif$model!="additive") && (cif$model!="fg"))
stop("Marginal Cumulative incidence model, must be either additive or extended Fine-Gray model\n")
cif.model <- switch(cif$model,additive=1,fg=2)
if (cause1!=attr(cif,"causeS")) cat("Cause for marginal model and correlation not the same\n");
if ((cause1[1]!=cause2[1])) {
if (is.null(cif2)==TRUE) stop("Must provide marginal model for both causes");
formula2<-attr(cif2,"Formula");
ldata2 <- timereg::aalen.des2(formula(delete.response(terms(formula2))),data=data,model="aalen");
X2<-ldata2$X; Z2<-ldata$Z;
if (is.null(Z2)==TRUE) {npar2<-TRUE; semi2<-0;} else {Z2<-as.matrix(Z2); npar2<-FALSE; semi2<-1;}
if (npar2==TRUE) {Z2<-matrix(0,antpers,1); pg2<-1; fixed2<-0;} else {fixed2<-1;pg2<-ncol(Z2);}
px2<-ncol(X2);
est2<-cif2$cum;
if (semi2==1) gamma2<-cif2$gamma else gamma2<-0;
est2<-cpred(est2,times);
} else {
X2<-matrix(0,1,1); Z2<-matrix(0,1,1); pg2<-1; px2<-1; semi2<-0; est2<-matrix(0,1,2); gamma2<-0;
npar2<-FALSE;
}
### For truncation
if (is.null(entry)) entry.call <- NULL else entry.call <- 0
if (is.null(entry)) entry <- rep(0,antpers);
cum1<-cpred(rbind(rep(0,px+1),cif$cum),entry)[,-1];
if (cif.model==1) cif1entry <- 1-exp(-apply(X*cum1,1,sum)- (Z %*% gamma )*entry)
else if (cif.model==2) cif1entry <- 1-exp(-apply(X*cum1,1,sum)*exp(Z %*% gamma ))
if ((model!="COR") && (cause1[1]!=cause2[1])) {
cum2<-cpred(rbind(rep(0,px2+1),cif2$cum),entry)[,-1];
if (cif.model==1) cif2entry <- 1-exp(-apply(X2*cum2,1,sum)- (Z2 %*% gamma2 )*entry)
else if (cif.model==2) cif2entry <- 1-exp(-apply(X2*cum2,1,sum)*exp(Z2 %*% gamma2 ))
} else {cum2 <- cum1; cif2entry <- cif1entry;}
## }}}
## {{{ censoring model stuff
cens.weight <- cif$cens.weights ### censoring weights from cif function
Gcxe <- 1;
### when censoring probs given, wants them so cens.model="user.weights"
if (!is.null(censoring.weights)) cens.model <- "user.weights"
if (cens.model!="user.weights") {
if (is.null(cens.weight) ) {
if (cens.model=="KM") { ## {{{
ud.cens<-survival::survfit(Surv(time,cause==cens.code)~+1);
Gfit<-cbind(ud.cens$time,ud.cens$surv)
Gfit<-rbind(c(0,1),Gfit);
Gcx<-cpred(Gfit,time)[,2];
if (is.null(entry.call)==FALSE) Gcxe<-cpred(Gfit,entry)[,2];
Gcx <- Gcx/Gcxe;
Gctimes<- Gcx ## }}}
} else if (cens.model=="cox") { ## {{{
if (npar==TRUE) XZ<-X[,-1] else XZ<-cbind(X,Z)[,-1];
ud.cens<-timereg::cox.aalen(Surv(time,cause==cens.code)~prop(XZ),n.sim=0,robust=0);
Gcx<-cpred(ud.cens$cum,time)[,2];
RR<-exp(XZ %*% ud.cens$gamma)
Gcx<-exp(-Gcx*RR)
Gfit<-rbind(c(0,1),cbind(time,Gcx));
if (is.null(entry.call)==FALSE) {
Gcxe<-cpred(Gfit,entry)[,2];
Gcxe<-exp(-Gcxe*RR)
}
Gcx <- Gcx/Gcxe;
Gctimes<- Gcx ## }}}
} else if (cens.model=="aalen") { ## {{{
if (npar==TRUE) XZ<-X[,-1] else XZ<-cbind(X,Z)[,-1];
ud.cens<-timereg::aalen(Surv(time,cause==cens.code)~XZ,n.sim=0,robust=0);
Gcx<-cpred(ud.cens$cum,time)[,-1];
XZ<-cbind(1,XZ);
Gcx<-exp(-apply(Gcx*XZ,1,sum))
Gcx[Gcx>1]<-1; Gcx[Gcx<0]<-0
Gfit<-rbind(c(0,1),cbind(time,Gcx));
if (is.null(entry.call)==FALSE) {
Gcxe<-cpred(Gfit,entry)[,-1];
Gcxe<-exp(-apply(Gcxe*XZ,1,sum))
Gcxe[Gcxe>1]<-1; Gcxe[Gcxe<0]<-0
}
Gcx <- Gcx/Gcxe;
Gctimes<- Gcx ## }}}
}
} else { Gcx <- cens.weight; Gctimes <- cens.weight;}
} else {
Gcx <- censoring.weights
Gctimes <- cens.weight
}
ntimes<-length(times);
## }}}
## {{{ set up cluster + theta design + define iid variables
if (is.null(clusters)== TRUE) {
## take clusters from cif model
clusters <- attr(cif,"clusters");
antclust<- length(unique(clusters));
max.clust <- attr(cif,"max.clust")
if (attr(cif,"coarse.clust"))
stop("Max.clust should be NULL in marginal model, or clusters should be given at call \n");
} else {
clus<-unique(clusters); antclust<-length(clus);
clusters <- as.integer(factor(clusters, labels = 1:(antclust)))-1;
}
outc <- cluster.index(clusters,index.type=TRUE);
clustsize <- outc$cluster.size
maxclust <- outc$maxclust
clusterindex <- outc$idclust
if (maxclust==1) stop("No clusters given \n");
if (model!="ARANCIF") {
if (is.null(theta.des)==TRUE) { ptheta<-1; theta.des<-matrix(1,antpers,ptheta);
colnames(theta.des) <- "intercept"; } else theta.des<-as.matrix(theta.des);
ptheta<-ncol(theta.des);
if (is.null(dimpar) & is.null(par.func)) dimpar<-ptheta;
}
if (model=="ARANCIF") { ### different parameters for Additive random effects
if (is.null(random.design)) random.design <- matrix(1,antpers,1);
dim.rv <- ncol(random.design);
if (is.null(theta.des)==TRUE) theta.des<-diag(dim.rv);
dimpar <- ncol(theta.des);
if (nrow(theta.des)!=ncol(random.design)) stop("nrow(theta.des)!= ncol(random.design)");
} else random.design <- matrix(0,1,1);
if ( (!is.null(par.func)) && is.null(dimpar) )
stop("Must specify dimension of parameters when specifying R-functions\n")
if ( (is.null(theta)==TRUE) && (model=="ARANCIF")) theta<-rep(0.5,dimpar);
if (is.null(theta)==TRUE) theta<-rep(0.1,dimpar);
if (length(theta)!=dimpar) theta<-rep(theta[1],dimpar);
Biid<-c(); gamma.iid <- 0; B2iid<-c(); gamma2.iid <- 0;
if (notaylor==0) {
for (i in 1:antclust) Biid<-cbind(Biid,cif$B.iid[[i]]);
if (npar==TRUE) gamma.iid<-0 else gamma.iid<-cif$gamma.iid;
if ((model!="COR") && (cause1!=cause2)) {
B2iid<-c()
for (i in 1:antclust) B2iid<-cbind(B2iid,cif2$B.iid[[i]]);
if (npar2==TRUE) gamma2.iid<-0 else gamma2.iid<-cif2$gamma.iid;
} else { B2iid<-Biid; gamma2.iid<-gamma.iid; }
}
var.theta<-hess<-matrix(0,dimpar,dimpar); score<-rep(0,dimpar);
time.pow<-attr(cif,"time.pow");
#if (sum(time.pow)==0) time.pow<-rep(1,pg);
theta.iid<-matrix(0,antclust,dimpar);
if ((nrow(theta.des)!=nrow(X)) && (model!="ARANCIF"))
stop("Dependence design not consistent with data\n");
if (length(trunkp)!=antpers) trunkp <- rep(1,antpers)
if (is.null(weights)==TRUE) weights <- rep(1,antpers);
## }}}
### ## {{{ function and derivative
if (is.null(par.func)==FALSE)
{
flex.func<-1; # use flexible design
if (is.null(dpar.func)) {
cat("Must provide derivative that is used for iid decomposition for SE's\n");
score.method <- "nlminb"
dpar.func <- par.func
}
} ## }}}
Zgamma <- c(Z %*% gamma);
Z2gamma2 <- c(Z2 %*% gamma2);
dep.model <- 0;
dep.model <- switch(model,COR=1,RR=2,OR=3,RANCIF=4,ARANCIF=5)
if (dep.model==0) stop("model must be COR, OR, RR, RANCIF, ARANCIF \n");
if (dep.model<=4) rvdes <- matrix(0,1,1);
obj <- function(par)
{ ## {{{
if (is.null(par.func)==TRUE) {
Xtheta <- theta.des %*% par; DXtheta <- array(0,c(1,1,1));
} else {
Xtheta <- c(); DXtheta <- array(0,c(length(times),antpers,dimpar));
if (is.null(dpar.func))
stop("Must also provide derivative of function wrt to parameters")
s <- 0
for (t in times) {
s <- 1+s
Xttheta <- par.func(par,t,theta.des)
if (length(Xttheta)!=antpers) stop("parfunc(par,t,theta.des) must length n");
Xtheta <- cbind(Xtheta,Xttheta)
Dttheta <- dpar.func(par,t,theta.des);
if (dim(Dttheta)[1]!=antpers || dim(Dttheta)[2]!=dimpar)
stop("dparfunc must return matrix n x dimpar when called on dparfunc(par,t,theta.des)");
DXtheta[s,,] <- Dttheta
}
}
outl<-.Call("cor", ## {{{
itimes=times,iy=time,icause=cause,iCA1=cause1,iKMc=Gcx,
iz=X,iest=matrix(est[,-1],ncol=ncol(est)-1),iZgamma=c(Zgamma),isemi=semi,izsem=Z,
itheta=c(par),iXtheta=Xtheta,iDXtheta=DXtheta,idimDX=dim(DXtheta),
ithetades=theta.des,
icluster=clusters,iclustsize=clustsize,iclusterindex=clusterindex,
iinverse=inverse,iCA2=cause2,
ix2=X2,isemi2=semi2,iest2=as.matrix(est2[,-1]),iZgamma2=c(Z2gamma2),
iflexfunc=flex.func,iiid=iid,isym=sym,iweights=weights,
isamecens=as.numeric(same.cens),istabcens=as.numeric(stab.cens),
iKMtimes=Gctimes,isilent=silent,
icifmodel=cif.model,idepmodel=dep.model,
iestimator=estimator,ientryage=entry,icif1entry=cif1entry,
icif2entry=cif2entry,itrunkp=trunkp,irvdes=random.design,
## Biid,gamma.iid,time.pow,
## B2iid,gamma2.iid,body(htheta),body(dhtheta),new.env(),
PACKAGE="mets"
)
## }}}
attr(outl,"gradient") <-outl$score
if (oout==0) return(outl$ssf) else if (oout==1) return(sum(outl$score^2)) else return(outl)
} ## }}}
p <- theta
iid <- 0; ###no-iid representation for iterations
if (score.method=="nr") { ## {{{
oout <- 2; ### output control for obj
if (Nit>0)
for (i in 1:Nit)
{
oout <- 2
out <- obj(p)
hess <- out$Dscore
if (!is.na(sum(hess))) hessi <- lava::Inverse(out$Dscore) else hessi <- hess
if (detail==1) {## {{{
print(paste("Fisher-Scoring ===================: it=",i));
cat("theta:");print(c(p))
cat("score:");print(c(out$score));
cat("hess:"); print(hess);
}## }}}
delta <- hessi %*% out$score
## do not update last iteration
if (i<Nit) p <- p-delta* step
if (is.nan(sum(out$score))) break;
if (sum(abs(out$score))<tol) break;
if (max(delta)>20) { cat("NR increment > 20, lower step zize, increment= \n"); cat(delta); break; }
}
if (!is.nan(sum(p))) { ## {{{ iid decomposition
oout <- 2
theta <- p
iid <- 1;
out <- obj(p)
score <- out$score
hess <- out$Dscore
} ## }}}
if (detail==1 & Nit==0) {## {{{
print(paste("Fisher-Scoring ===================: final"));
cat("theta:");print(c(p))
cat("score:");print(c(out$score));
cat("hess:"); print(hess);
}## }}}
if (!is.na(sum(hess))) hessi <- lava::Inverse(out$Dscore) else hessi <- diag(nrow(hess))
score1 <- score;
## }}}
} else if (score.method=="nlminb") { ## {{{ nlminb optimizer
iid <- 0; oout <- 0;
tryCatch(opt <- nlminb(theta,obj,control=control),error=function(x) NA)
if (detail==1) print(opt);
iid <- 1;
hess <- numDeriv::hessian(obj,opt$par)
score <- numDeriv::jacobian(obj,opt$par)
hessi <- lava::Inverse(hess);
theta <- opt$par
if (detail==1) cat("iid decomposition\n");
oout <- 2;
out <- obj(opt$par)
score1 <- out$score
## }}}
} else if (score.method=="nlm") { ## {{{ nlm optimizer
iid <- 0; oout <- 0;
tryCatch(opt <- nlm(obj,theta,hessian=TRUE,print.level=detail),error=function(x) NA)
iid <- 1;
hess <- opt$hessian
score <- opt$gradient
if (detail==1) print(opt);
hessi <- lava::Inverse(hess);
theta <- opt$estimate
if (detail==1) cat("iid decomposition\n");
oout <- 2;
out <- obj(opt$estimate)
score1 <- out$score
## }}}
} else stop("score.methods = nlm nlminb nr\n");
theta.iid <- out$theta.iid %*% hessi
if (is.null(par.func)) var.theta <- t(theta.iid) %*% theta.iid else var.theta <- hessi
var.theta <- t(theta.iid) %*% theta.iid
if (is.null(par.func)) thetanames <- colnames(theta.des) else
thetanames <- rep("R-func",dimpar)
ud <- list(theta=theta,score=score,hess=hess,hessi=hessi,var.theta=var.theta,
theta.iid=theta.iid,score1=c(score1),thetanames=thetanames,
brierscore=out$ssf,p11=out$p11);
if (dep.model<=3) class(ud)<-"cor"
else if (dep.model==4) class(ud) <- "randomcif"
else if (dep.model==6) class(ud) <- "randomcif"
else if (dep.model==5) class(ud) <- "randomcifrv"
attr(ud, "Formula") <- formula
attr(ud, "Clusters") <- clusters
attr(ud,"cause1")<-cause1; attr(ud,"cause2")<-cause2
attr(ud,"sym")<-sym;
attr(ud,"inverse")<-inverse;
attr(ud,"antpers")<-antpers;
attr(ud,"antclust")<-antclust;
attr(ud,"var.link")<-var.link;
if (dep.model==4) attr(ud, "Type") <- "randomcif"
if (dep.model==6) attr(ud, "Type") <- "randomcif"
if (model=="COR") attr(ud, "Type") <- "cor"
if (model=="RR") attr(ud, "Type") <- "RR"
if (model=="OR") attr(ud, "Type") <- "OR-cif"
if (dep.model==5) attr(ud, "pardes") <- theta.des
if (dep.model==5) attr(ud, "rv1") <- random.design[1,]
return(ud);
} ## }}}
##' Cross-odds-ratio, OR or RR risk regression for competing risks
##'
##' Fits a parametric model for the log-cross-odds-ratio for the
##' predictive effect of for the cumulative incidence curves for \eqn{T_1}
##' experiencing cause i given that \eqn{T_2} has experienced a cause k :
##' \deqn{
##' \log(COR(i|k)) = h(\theta,z_1,i,z_2,k,t)=_{default} \theta^T z =
##' }
##' with the log cross odds ratio being
##' \deqn{
##' 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
##' \deqn{
##' 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
##' \deqn{
##' 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 \eqn{B^c} is the complement event of \eqn{B},
##' \eqn{P_x} is the distribution given covariates
##' (\eqn{x} are subject specific and \eqn{z} are cluster specific covariates), and
##' \eqn{h()} is a function that is the simple identity
##' \eqn{\theta^T z} by default.
##'
##' The OR dependence measure is given by
##' \deqn{
##' 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
##' \deqn{
##' 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 \eqn{T_1} and \eqn{T_2} can be interchanged and leads to
##' the same cross-odd-ratio (i.e.
##' \eqn{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).
##'
##' \eqn{h()} may be specified as an R-function of the parameters,
##' see example below, but the default is that it is simply \eqn{\theta^T z}.
##'
##' @aliases or.cif rr.cif
##' @param 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
##' @param data a data.frame with the variables.
##' @param cause specifies the causes related to the death
##' times, the value cens.code is the censoring value. When missing it comes from marginal cif.
##' @param times time-vector that specifies the times used for the estimating euqations for the cross-odds-ratio estimation.
##' @param cause1 specificies the cause considered.
##' @param cause2 specificies the cause that is conditioned on.
##' @param cens.code specificies the code for the censoring if NULL then uses the one from the marginal cif model.
##' @param cens.model specified which model to use for the ICPW, KM is Kaplan-Meier alternatively it may be "cox"
##' @param Nit number of iterations for Newton-Raphson algorithm.
##' @param detail if 0 no details are printed during iterations, if 1 details are given.
##' @param clusters specifies the cluster structure.
##' @param theta specifies starting values for the cross-odds-ratio parameters of the model.
##' @param theta.des specifies a regression design for the cross-odds-ratio parameters.
##' @param step specifies the step size for the Newton-Raphson algorithm.
##' @param sym specifies if symmetry is used in the model.
##' @param weights weights for estimating equations.
##' @param par.func parfunc
##' @param dpar.func dparfunc
##' @param dimpar dimpar
##' @param score.method "nlminb", can also use "nr".
##' @param same.cens if true then censoring within clusters are assumed to be the same variable, default is independent censoring.
##' @param censoring.weights these probabilities are used for the bivariate censoring dist.
##' @param silent 1 to suppress output about convergence related issues.
##' @param ... Not used.
##' @return returns an object of type 'cor'. With the following arguments:
##' \item{theta}{estimate of proportional odds parameters of model.}
##' \item{var.theta}{variance for gamma. }
##' \item{hess}{the derivative of the used score.}
##' \item{score}{scores at final stage.}
##' \item{score}{scores at final stage.}
##' \item{theta.iid}{matrix of iid decomposition of parametric effects.}
##' @author 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
##' ###########################################################
##' \donttest{ ## 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)
##' }
##' @export
##' @keywords survival
cor.cif<-function(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,...)
{ ## {{{
fit <- dep.cif(cif=cif,data=data,cause=cause,model="COR",times=times,
cause1=cause1,cause2=cause2,cens.code=cens.code,cens.model=cens.model,Nit=Nit,detail=detail,
clusters=clusters,theta=theta,theta.des=theta.des,par.func=par.func,dpar.func=dpar.func,
dimpar=dimpar,
step=step,sym=sym,weights=weights,
score.method=score.method,same.cens=same.cens,censoring.weights=censoring.weights,silent=silent,...)
fit$call <- match.call()
fit
} ## }}}
##' @export
rr.cif<-function(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, step=1,sym=0,weights=NULL,
same.cens=FALSE,censoring.weights=NULL,silent=1,par.func=NULL,dpar.func=NULL,dimpar=NULL,
score.method="nlminb",entry=NULL,estimator=1,trunkp=1,admin.cens=NULL,...)
{ ## {{{
fit <- dep.cif(cif=cif,data=data,cause=cause,model="RR",cif2=cif2,times=times,
cause1=cause1,cause2=cause2,cens.code=cens.code,cens.model=cens.model,Nit=Nit,detail=detail,
clusters=clusters,theta=theta,theta.des=theta.des, step=step,sym=sym,weights=weights,
same.cens=same.cens,censoring.weights=censoring.weights,silent=silent,
par.func=par.func,dpar.func=dpar.func,dimpar=dimpar, score.method=score.method,
entry=entry,estimator=estimator,trunkp=trunkp,admin.cens=admin.cens,...)
fit$call <- match.call()
fit
} ## }}}
##' @export
or.cif<-function(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, step=1,sym=0, weights=NULL,
same.cens=FALSE,censoring.weights=NULL,silent=1,par.func=NULL,dpar.func=NULL,dimpar=NULL,
score.method="nlminb",entry=NULL,estimator=1,trunkp=1,admin.cens=NULL,...)
{ ## {{{
fit <- dep.cif(cif=cif,data=data,cause=cause,model="OR",cif2=cif2,times=times,
cause1=cause1,cause2=cause2,cens.code=cens.code,cens.model=cens.model,Nit=Nit,detail=detail,
clusters=clusters,theta=theta,theta.des=theta.des, step=step,sym=sym,weights=weights,
same.cens=same.cens,censoring.weights=censoring.weights,silent=silent,
par.func=par.func,dpar.func=dpar.func,dimpar=dimpar,
score.method=score.method,
entry=entry,estimator=estimator,trunkp=trunkp,admin.cens=admin.cens,...)
fit$call <- match.call()
fit
} ## }}}
##' 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 the form
##' \deqn{
##' P(T \leq t, cause=1 | x,z) = P_1(t,x,z) = 1- exp( -x^T A(t) exp(z^T \beta))
##' }
##' We allow a regression structure for the random effects variances that may depend on
##' cluster covariates.
##'
##' @title Random effects model for competing risks data
##' @param cif a model object from the comp.risk function with the
##' marginal cumulative incidence of cause2, i.e., the event that is conditioned on, and whose
##' odds the comparision is made with respect to
##' @param data a data.frame with the variables.
##' @param cause specifies the causes related to the death
##' times, the value cens.code is the censoring value.
##' @param cif2 specificies model for cause2 if different from cause1.
##' @param cause1 cause of first coordinate.
##' @param cause2 cause of second coordinate.
##' @param cens.code specificies the code for the censoring if NULL then uses the one from the marginal cif model.
##' @param cens.model specified which model to use for the ICPW, KM is Kaplan-Meier alternatively it may be "cox"
##' @param Nit number of iterations for Newton-Raphson algorithm.
##' @param detail if 0 no details are printed during iterations, if 1 details are given.
##' @param clusters specifies the cluster structure.
##' @param theta specifies starting values for the cross-odds-ratio parameters of the model.
##' @param theta.des specifies a regression design for the cross-odds-ratio parameters.
##' @param sym 1 for symmetry 0 otherwise
##' @param step specifies the step size for the Newton-Raphson algorith.m
##' @param same.cens if true then censoring within clusters are assumed to be the same variable, default is independent censoring.
##' @param var.link if var.link=1 then var is on log-scale.
##' @param score.method default uses "nlminb" optimzer, alternatively, use the "nr" algorithm.
##' @param entry entry-age in case of delayed entry. Then two causes must be given.
##' @param trunkp gives probability of survival for delayed entry, and related to entry-ages given above.
##' @param ... extra arguments.
##' @return returns an object of type 'cor'. With the following arguments:
##' \item{theta}{estimate of proportional odds parameters of model.}
##' \item{var.theta}{variance for gamma. }
##' \item{hess}{the derivative of the used score.}
##' \item{score}{scores at final stage.}
##' \item{score}{scores at final stage.}
##' \item{theta.iid}{matrix of iid decomposition of parametric effects.}
##' @export
##' @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 (2012), work in progress.
##' @examples
##' \donttest{ ## Reduce Ex.Timings
##' d <- simnordic.random(5000,delayed=TRUE,cordz=0.5,cormz=2,lam0=0.3,country=TRUE)
##' times <- seq(50,90,by=10)
##' add1 <- 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)
##'
##' out1<-random.cif(add1,data=d,cause1=1,cause2=1,theta=1,same.cens=TRUE)
##' summary(out1)
##'
##' out2<-random.cif(add1,data=d,cause1=1,cause2=1,theta=1,
##' theta.des=mm,same.cens=TRUE)
##' summary(out2)
##'
##' #########################################
##' ##### 2 different causes
##' #########################################
##'
##' add2 <- timereg::comp.risk(Event(time,cause)~-1+factor(country)+cluster(id),data=d,
##' times=times,cause=2,max.clust=NULL)
##' out3 <- random.cif(add1,data=d,cause1=1,cause2=2,cif2=add2,sym=1,same.cens=TRUE)
##' summary(out3) ## negative dependence
##'
##' out4 <- random.cif(add1,data=d,cause1=1,cause2=2,cif2=add2,theta.des=mm,sym=1,same.cens=TRUE)
##' summary(out4) ## negative dependence
##' }
##' @keywords survival
##' @author Thomas Scheike
random.cif<-function(cif,data,cause=NULL,cif2=NULL,
cause1=1,cause2=1,cens.code=NULL,cens.model="KM",Nit=40,detail=0,
clusters=NULL,theta=NULL,theta.des=NULL,sym=1,
step=1,same.cens=FALSE,var.link=0,score.method="nr",
entry=NULL,trunkp=1,...)
{ ## {{{
fit <- dep.cif(cif,data=data,cause=cause,model="RANCIF",cif2=cif2,sym=sym,
cause1=cause1,cause2=cause2,cens.code=cens.code,cens.model=cens.model,Nit=Nit,detail=detail,
clusters=clusters,theta=theta,theta.des=theta.des,step=step,same.cens=same.cens,
var.link=var.link,score.method=score.method,entry=entry,trunkp=trunkp,...)
fit$call <- match.call()
fit
} ## }}}
##'Additive Random effects model for competing risks data for polygenetic modelling
##'
##'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
##'\deqn{
##'P(T \leq t, cause=1 | x,z) = P_1(t,x,z) = 1- exp( -x^T A(t) - t z^T \beta)
##'}
##'
##'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
##' \eqn{v_1} and \eqn{v_2}.
##'Such that the random effects for subject
##'1 is \deqn{v_1^T (Z_1,...,Z_d)}, for d random effects. Each random effect
##'has an associated parameter \eqn{(\lambda_1,...,\lambda_d)}. By construction
##'subjects 1's random effect are Gamma distributed with
##'mean \eqn{\lambda_1/v_1^T \lambda}
##'and variance \eqn{\lambda_1/(v_1^T \lambda)^2}. Note that the random effect
##'\eqn{v_1^T (Z_1,...,Z_d)} has mean 1 and variance \eqn{1/(v_1^T \lambda)}.
##'
##'The parameters \eqn{(\lambda_1,...,\lambda_d)}
##'are related to the parameters of the model
##'by a regression construction \eqn{pard} (d x k), that links the \eqn{d}
##'\eqn{\lambda} parameters
##'with the (k) underlying \eqn{\theta} parameters
##'\deqn{
##'\lambda = pard \theta
##'}
##'
##'
##'
##' @param cif a model object from the timereg::comp.risk function with the
##' marginal cumulative incidence of cause2, i.e., the event that is conditioned on, and whose
##' odds the comparision is made with respect to
##' @param data a data.frame with the variables.
##' @param cause specifies the causes related to the death
##' times, the value cens.code is the censoring value.
##' @param cif2 specificies model for cause2 if different from cause1.
##' @param times time points
##' @param cause1 cause of first coordinate.
##' @param cause2 cause of second coordinate.
##' @param cens.code specificies the code for the censoring if NULL then uses the one from the marginal cif model.
##' @param cens.model specified which model to use for the ICPW, KM is Kaplan-Meier alternatively it may be "cox"
##' @param Nit number of iterations for Newton-Raphson algorithm.
##' @param detail if 0 no details are printed during iterations, if 1 details are given.
##' @param clusters specifies the cluster structure.
##' @param theta specifies starting values for the cross-odds-ratio parameters of the model.
##' @param theta.des specifies a regression design for the cross-odds-ratio parameters.
##' @param step specifies the step size for the Newton-Raphson algorith.m
##' @param sym 1 for symmetri and 0 otherwise
##' @param weights weights for score equations.
##' @param same.cens if true then censoring within clusters are assumed to be the same variable, default is independent censoring.
##' @param censoring.weights Censoring probabilities
##' @param silent debug information
##' @param var.link if var.link=1 then var is on log-scale.
##' @param score.method default uses "nlminb" optimzer, alternatively, use the "nr" algorithm.
##' @param entry entry-age in case of delayed entry. Then two causes must be given.
##' @param estimator estimator
##' @param trunkp gives probability of survival for delayed entry, and related to entry-ages given above.
##' @param admin.cens Administrative censoring
##' @param random.design specifies a regression design of 0/1's for the random effects.
##' @param ... extra arguments.
##' @return returns an object of type 'random.cif'. With the following arguments:
##'\item{theta}{estimate of parameters of model.}
##'\item{var.theta}{variance for gamma. }
##'\item{hess}{the derivative of the used score.}
##'\item{score}{scores at final stage.}
##'\item{theta.iid}{matrix of iid decomposition of parametric effects.}
##' @export
##' @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
##' \donttest{ ## Reduce Ex.Timings
##' d <- simnordic.random(5000,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)
##'
##' }
##' @keywords survival
##' @author Thomas Scheike
##' @export
Grandom.cif<-function(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,...)
{ ## {{{
fit <- dep.cif(cif,data=data,cause=cause,model="ARANCIF",cif2=cif2,times=times,
cause1=cause1,cause2=cause2,cens.code=cens.code,cens.model=cens.model,Nit=Nit,detail=detail,
clusters=clusters,theta=theta,theta.des=theta.des,step=step,sym=sym,weights=weights,
same.cens=same.cens,censoring.weights=censoring.weights,silent=silent,var.link=var.link,
score.method=score.method,entry=entry,estimator=estimator,
random.design=random.design,trunkp=trunkp,admin.cens=admin.cens,...)
fit$call <- match.call()
fit
} ## }}}
##' @export
print.summary.cor <- function(x,digits=3,...)
{ ## {{{
if (x$type=="cor") {
cat("Cross odds ratio dependence for competing risks\n\n")
cat("Effect of cause1=",x$cause1," on cause2=",x$cause2,
" under symmetry=",x$sym,fill=TRUE,sep="")
} else if (x$type=="RR") {
cat("Ratio of joint and product of marginals for competing risks\n\n")
cat("Ratio of cumulative incidence for cause1=",x$cause1," and cause2=",x$cause2,sep=" ")
} else if (x$type=="OR-cif") {
cat("OR for dependence for competing risks\n\n")
cat("OR of cumulative incidence for cause1=",x$cause1," and cause2=",x$cause2,sep=" ")
}
cat("\n")
prmatrix(signif(x$estimates,digits))
cat("\n")
if (!is.null(x$marg)) {
cat(paste("Marginal cumulative incidencen",signif(x$marg,digits),"\n"))
prmatrix(signif(x$casewise,digits))
prmatrix(signif(x$concordance,digits))
cat("\n")
}
invisible(x)
} ## }}}
##' Computes concordance and casewise concordance for dependence models for competing risks
##' models of the type cor.cif, rr.cif or or.cif for the given cumulative incidences and the different dependence
##' measures in the object.
##'
##'
##' @title Summary for dependence models for competing risks
##' @param object object from cor.cif rr.cif or or.cif for dependence between competing risks data for two causes.
##' @param marg.cif a number that gives the cumulative incidence in one time point for which concordance and
##' casewise concordance are computed.
##' @param marg.cif2 the cumulative incidence for cause 2 for concordance and
##' casewise concordance are computed. Default is that it is the same as marg.cif.
##' @param digits digits in output.
##' @param ... Additional arguments.
##' @return prints summary for dependence model.
##' \item{casewise}{gives casewise concordance that is, probability of cause 2 (related to cif2) given that cause 1 (related to cif1)
##' has occured.}
##' \item{concordance}{gives concordance that is, probability of cause 2 (related to cif2) and cause 1 (related to cif1).}
##' \item{cif1}{cumulative incidence for cause1.}
##' \item{cif2}{cumulative incidence for cause1.}
##' @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.
##'
##' @author Thomas Scheike
##' @keywords survival
##' @examples
##' ## library("timereg")
##' ## data("multcif",package="mets") # simulated data
##' ## multcif$cause[multcif$cause==0] <- 2
##' ##
##' ## times=seq(0.1,3,by=0.1) # to speed up computations use only these time-points
##' ## add <- timereg::comp.risk(Event(time,cause)~+1+cluster(id),
##' ## data=multcif,n.sim=0,times=times,cause=1)
##' ###
##' ## out1<-cor.cif(add,data=multcif,cause1=1,cause2=1,theta=log(2+1))
##' ## summary(out1)
##' ##
##' ## pad <- predict(add,X=1,se=0,uniform=0)
##' ## summary(out1,marg.cif=pad)
##' @method summary cor
##' @export
summary.cor <- function(object,marg.cif=NULL,marg.cif2=NULL,digits=3,...) { ## {{{
if (!(attr(object,"class") %in% c("cor","randomcif"))) stop("Must be a cor.cif or randomcif object")
if (sum(abs(object$score))>0.001) warning("WARNING: check score for convergence\n")
coefs <- coef(object,...)
ocasewise <- oconcordance <- NULL
if (is.null(marg.cif)==FALSE) { ## {{{
time <- marg.cif$time
marg.cif <- marg.cif$P1
if (is.null(marg.cif2)==FALSE) marg.cif2 <- marg.cif2$P1 else {
if (attr(object,"cause2")==attr(object,"cause1")) marg.cif2 <- marg.cif
else stop("causes not the same and second marginal cif not given\n");
}
pmarg.cif <- marg.cif*marg.cif2
thetav <- coefs[,1]
thetavl <- thetav-1.96*coefs[,2]
thetavu <- thetav+1.96*coefs[,2]
namev <- object$thetanames
if (is.null(namev)) namev <- paste("name",1:length(thetav),sep="")
ocasewise <- oconcordance <- list()
k <- 0
for (theta in thetav) {
k <- k+1
thetal <- thetavl[k];
thetau <- thetavu[k]
if (attr(object,"Type")=="cor") { ## {{{
concordance <- exp(theta)*pmarg.cif/((1-marg.cif)+exp(theta)*marg.cif)
conclower <- exp(thetal)*pmarg.cif/((1-marg.cif)+exp(thetal)*marg.cif)
concup <- exp(thetau)*pmarg.cif/((1-marg.cif)+exp(thetau)*marg.cif)
casewise <- concordance/marg.cif
caselower <- conclower/marg.cif
caseup <- concup/marg.cif
} else if (attr(object,"Type")=="RR") {
casewise<- exp(coefs[,1])*c(marg.cif)
concordance <- exp(coefs[,1])*pmarg.cif
caselower <- marg.cif*exp(thetal)
caseup <- marg.cif*exp(thetau)
conclower <- pmarg.cif* exp(thetal)
concup <- pmarg.cif*exp(thetau)
} else if (attr(object,"Type")=="OR-cif") {
casewise<- plack.cif2(marg.cif,marg.cif,c(theta))/marg.cif
concordance <- plack.cif2(marg.cif,marg.cif,c(theta))
caselower <- plack.cif2(marg.cif,marg.cif,thetal)/marg.cif
caseup <- plack.cif2(marg.cif,marg.cif,thetau)/marg.cif
conclower <- plack.cif2(marg.cif,marg.cif,thetal)
concup <- plack.cif2(marg.cif,marg.cif,thetau)
} else if (attr(object,"Type")=="randomcif") {
theta <- 1/theta
thetal <- 1/thetal
thetau <- 1/thetau
lam <- 1-marg.cif
p11<- 1-lam -lam +lap(theta,2*ilap(theta, lam))
concordance <- p11
conclower <- 1-lam -lam +lap(thetal,2*ilap(thetal, lam))
concup <- 1-lam -lam +lap(thetau,2*ilap(thetau, lam))
casewise <- concordance/marg.cif
caselower <- conclower/marg.cif
caseup <- concup/marg.cif
} ## }}}
outcase <- cbind(time,casewise,caselower,caseup)
outconc <- cbind(time,concordance,conclower,concup)
colnames(outcase) <- c("time","casewise concordance","2.5 %","97.5%")
colnames(outconc) <- c("time","concordance","2.5 %","97.5%")
if (length(thetav)==1) {
ocasewise <- outcase
oconcordance <- outconc
} else {
ocasewise[[k]] <- outcase
oconcordance[[k]] <- outconc
}
}
if (length(thetav)>1) {
if (length(ocasewise)==length(namev)) {
names(ocasewise) <- namev
names(oconcordance) <- namev
}
}
} ## }}}
res <- list(casewise=ocasewise,concordance=oconcordance,estimates=coefs,
marg.cif=marg.cif, marg.cif2=marg.cif2,type=attr(object,"Type"),
sym=attr(object,"sym"),cause1=attr(object,"cause1"),cause2=attr(object,"cause2"))
class(res) <- "summary.cor"
res
} ## }}}
##' @export
coef.cor <- function(object,...)
{ ## {{{
res <- cbind(object$theta, diag(object$var.theta)^0.5)
se<-diag(object$var.theta)^0.5
wald <- object$theta/se
waldp <- (1 - pnorm(abs(wald))) * 2
cor<-exp(object$theta)
res <- as.matrix(cbind(res, wald, waldp,cor,se*cor))
if (attr(object,"Type")=="cor")
colnames(res) <- c("log-Coef.", "SE", "z", "P-val","Cross odds ratio","SE")
else colnames(res) <- c("log-ratio Coef.", "SE", "z", "P-val","Ratio","SE")
if (!is.null(object$thetanames)) rownames(res)<-object$thetanames
return(res)
} ## }}}
##' @export
print.cor<-function(x,digits=3,...)
{ ## {{{
print(x$call);
cat("\n")
print(summary(x));
} ## }}}
##' @export
concordance.cor <- function(object,...) {
concordanceCor(object,...)
}
##' Concordance for Twins
##'
##' The concordance is the probability that both twins have experienced the
##' event of interest and is defined as
##' \deqn{
##' cor(t) = P(T_1 \leq t, \epsilon_1 =1 , T_2 \leq t, \epsilon_2=1)
##' }
##'
##' Similarly, the casewise concordance is
##' \deqn{
##' 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.
##'
##' @title Concordance Computes concordance and casewise concordance
##' @param object Output from the cor.cif, rr.cif or or.cif function
##' @param cif1 Marginal cumulative incidence
##' @param cif2 Marginal cumulative incidence of other cause (cause2) if it is different from cause1
##' @param messages To print messages
##' @param model Specfifies wich model that is considered if object not given.
##' @param coefs Specfifies dependence parameters if object is not given.
##' @param ... Extra arguments, not used.
##' @author Thomas Scheike
##' @aliases concordanceCor concordance.cor
##' @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
##'
##' @export
concordanceCor <- function(object,cif1,cif2=NULL,messages=TRUE,model=NULL,coefs=NULL,...)
{ ## {{{
if (is.null(model)) {
if (!inherits(object, "cor")) stop("Must be a rr.cif, cor.cif or or.cif object")
model <- attr(object,"Type")
}
if (is.null(coefs)) coefs <- coef(object)
if (is.null(coefs)) stop("Must give dependence parameters\n");
if (!is.null(object)) {
cause1 <- attr(object,"cause1");
cause2 <- attr(object,"cause2");
} else cause1 <- cause2 <- 1
if (is.null(cif2)==TRUE) cif2 <- cif1;
if (messages) {
if (model=="cor") { ## {{{
message("Cross odds ratio dependence for competing risks\n\n")
message("Odds of cause1=",cause1," given cause2=",cause2," relative to Odds of cause1=",cause1,"\n",fill=TRUE,sep="")
} else if (model=="RR") {
message("Ratio of joint and product of marginals for competing risks\n\n")
message("Ratio of cumulative incidence for cause1=",cause1," and cause2=",cause2,sep=" ")
} else if (model=="OR-cif") {
message("OR for dependence for competing risks\n\n")
message("OR of cumulative incidence for cause1=",cause1," and cause2=",cause2,sep=" ")
} ## }}}
}
out <- list()
for (k in 1:nrow(coefs)) {
## {{{
if (model=="cor") {
concordance <- exp(coefs[k,1])*cif1*cif2/((1-cif1)+exp(coefs[k,1])*cif1)
conclower <- exp(coefs[k,1]-1.96*coefs[k,2])*cif1*cif2/((1-cif1)+exp(coefs[k,1]-1.96*coefs[k,2])*cif1)
concup <- exp(coefs[k,1]+1.96*coefs[k,2])*cif1*cif2/((1-cif1)+exp(coefs[k,1]+1.96*coefs[k,2])*cif1)
casewise <- concordance/cif1
caselower <- conclower/cif1
caseup <- concup/cif1
} else if (model=="RR") {
casewise<- exp(coefs[k,1])*c(cif2)
concordance <- exp(coefs[k,1])*cif1*cif2
caselower <- cif2*exp(coefs[k,1]-1.96*coefs[k,2])
caseup <- cif2*exp(coefs[k,1]+1.96*coefs[k,2])
conclower <- cif1*cif2* exp(coefs[k,1]-1.96*coefs[k,2])
concup <- cif1*cif2*exp(coefs[k,1]+1.96*coefs[k,2])
} else if (model=="OR-cif") {
thetal <- coefs[k,1]-1.96*coefs[k,2]
thetau <- coefs[k,1]+1.96*coefs[k,2]
casewise<- plack.cif2(cif1,cif2,c(coefs[k,1]))/cif1
concordance <- plack.cif2(cif1,cif2,c(coefs[k,1]))
caselower <- plack.cif2(cif1,cif2,thetal)/cif1
caseup <- plack.cif2(cif1,cif2,thetau)/cif1
conclower <- plack.cif2(cif1,cif2,thetal)
concup <- plack.cif2(cif1,cif2,thetau)
}
outcase <- cbind(c(casewise),c(caselower),c(caseup))
outconc <- cbind(c(concordance),c(conclower),c(concup))
colnames(outcase) <- c("casewise concordance","2.5 %","97.5%")
colnames(outconc) <- c("concordance","2.5 %","97.5%")
## }}}
out[[k]] <- list(concordance=outconc,casewise=outcase)
names(out)[k] <- rownames(coefs)[k]
}
return(out)
} ## }}}
##' .. content for description (no empty lines) ..
##'
##' @title plack Computes concordance for or.cif based model, that is Plackett random effects model
##' @aliases plack.cif2
##' @export
##' @param cif1 Cumulative incidence of first argument.
##' @param cif2 Cumulative incidence of second argument.
##' @param object or.cif object with dependence parameters.
##' @author Thomas Scheike
plack.cif <- function(cif1,cif2,object)
{ ## {{{
coefs <- coef(object)
theta <- exp(object$theta);
cif1 <- c(cif1); cif2 <- c(cif2)
cifs=cif1+cif2;
valn=2*(theta-1);
val1=(1+(theta-1)*(cifs))-( ((1+(theta-1)*cifs))^2-4*cif1*cif2*theta*(theta-1))^0.5;
vali=cif1*cif2;
valr <- vali;
valr[valn!=0] <- val1/valn;
valr <- matrix(valr,length(c(theta)),1)
rownames(valr)=colnames(coefs)
return(valr);
} ## }}}
##' @export
plack.cif2 <- function(cif1,cif2,theta)
{ ## {{{
theta <- exp(c(theta))
cif1 <- c(cif1); cif2 <- c(cif2)
cifs=cif1+cif2;
valn=2*(theta-1);
val1=(1+(theta-1)*(cifs))-( ((1+(theta-1)*cifs))^2-4*cif1*cif2*theta*(theta-1))^0.5;
vali=cif1*cif2;
valr <- vali;
valr[valn!=0] <- val1/valn;
return(valr);
} ## }}}
##' @export
predictPairPlack <- function(cif1,cif2,status1,status2,theta)
{ ## {{{
theta <- exp(c(theta))
cif1 <- c(cif1); cif2 <- c(cif2)
cifs=cif1+cif2;
valn=2*(theta-1);
val1=(1+(theta-1)*(cifs))-( ((1+(theta-1)*cifs))^2-4*cif1*cif2*theta*(theta-1))^0.5;
vali=cif1*cif2;
valr <- vali;
valr[valn!=0] <- val1/valn;
p11 <- valr;
p10 <- cif1-p11
p01 <- cif2-p11
p00 <- 1- p10-p01-p11
#
pred <- (status1==1)*(status2==1)*p11+ (status1==1)*(status2==0)*p10+
(status1==0)*(status2==1)*p01+ (status1==0)*(status2==0)*p00
return(pred);
} ## }}}
##' @export
summary.randomcif<-function (object, ...)
{ ## {{{
if (!inherits(object, "randomcif"))
stop("Must be a random.cif object")
cat("Random effect variance for variation due to clusters\n\n")
cat("Cause", attr(object, "cause1"), "and cause", attr(object,
"cause2"), fill = TRUE)
cat("\n")
if (sum(abs(object$score)) > 0.001) cat("WARNING: check score for convergence")
cat("\n")
coef.randomcif(object, ...)
} ## }}}
##' @export
coef.randomcif<- function (object, digits = 3, ...)
{ ## {{{
res <- cbind(object$theta, diag(object$var.theta)^0.5)
se <- diag(object$var.theta)^0.5
wald <- object$theta/se
waldp <- (1 - pnorm(abs(wald))) * 2
cor <- object$theta + 1
res <- as.matrix(cbind(res, wald, waldp, cor, se))
colnames(res) <- c("Coef.", "SE", "z", "P-val", "Cross odds ratio", "SE")
if (!is.null(object$thetanames)) rownames(res)<-object$thetanames
return(res)
} ## }}}
##' @export
print.randomcif<- function (x , digits = 3, ...)
{ ## {{{
} ## }}}
##' @export
summary.randomcifrv<-function (object, ...)
{ ## {{{
if (!inherits(object, "randomcifrv"))
stop("Must be a random.cifrv object")
cat("Random effect parameters for additive gamma random effects \n\n")
cat("Cause", attr(object, "cause1"), "and cause", attr(object,
"cause2"), fill = TRUE)
cat("\n")
if (sum(abs(object$score)) > 1e-06)
cat("WARNING: check score for convergence")
cat("\n")
res <- coef.randomcifrv(object, ...)
var.link <- attr(object,"var.link");
rv1 <- attr(object,"rv1");
theta.des <- attr(object,"pardes");
if (var.link==1) par <- theta.des %*% exp(object$theta) else par <- theta.des %*% object$theta
if (var.link==1) {
fp <- function(p){ res <- exp(p)/sum(rv1* (theta.des %*% exp(p))); return(res); }
e <- lava::estimate(coef=object$theta,vcov=object$var.theta,f=function(p) fp(p))
pare <- lava::estimate(coef=object$theta,vcov=object$var.theta,f=function(p) exp(p))
res <- list(estimate=res,h=e,exppar=pare)
} else {
fp <- function(p) { p/sum(rv1* (theta.des %*% p)) }
e <- lava::estimate(coef=object$theta,vcov=object$var.theta,f=function(p) fp(p))
res <- list(estimate=res,h=e)
}
res
} ## }}}
##' @export
coef.randomcifrv<- function (object, digits = 3, ...)
{ ## {{{
if (attr(object,"inverse")==1) elog <- 1 else elog <- 0;
if (elog==1) theta <- exp(object$theta) else theta <- object$theta
se <- diag(object$var.theta)^0.5
res <- cbind(object$theta, se)
wald <- object$theta/se
waldp <- (1 - pnorm(abs(wald))) * 2
res <- as.matrix(cbind(res, wald, waldp))
if (elog==0) colnames(res) <- c("Coef.", "SE", "z", "P-val")
if (elog==1) res <- cbind(res,exp(object$theta), exp(object$theta)^2*se)
if (elog==1) colnames(res) <- c("log-parameter","SE","z","P-val","exp(theta)","SE")
if (!is.null(object$thetanames)) rownames(res)<-object$thetanames
prmatrix(signif(res, digits))
} ## }}}
##' @export
print.randomcifrv<- function (x , digits = 3, ...)
{ ## {{{
summary(x, ...)
} ## }}}
## summary.cor<-function(object,digits=3,marg.cif=NULL,...)
## { ## {{{
## if (!inherits(object, "cor")) stop("Must be a cor.cif object")
## if (sum(abs(object$score))>0.001) warning("WARNING: check score for convergence\n")
## coefs <- coef.cor(object,...);
## outcase <- outconc <- NULL
## if (is.null(marg.cif)==FALSE) {
## marg.cif <- max(marg.cif)
## ## {{{
## if (attr(object,"Type")=="cor") {
## concordance <- exp(coefs[,1])*marg.cif^2/((1-marg.cif)+exp(coefs[,1])*marg.cif)
## conclower <- exp(coefs[,1]-1.96*coefs[,2])*marg.cif^2/((1-marg.cif)+exp(coefs[,1]-1.96*coefs[,2])*marg.cif)
## concup <- exp(coefs[,1]+1.96*coefs[,2])*marg.cif^2/((1-marg.cif)+exp(coefs[,1]+1.96*coefs[,2])*marg.cif)
## casewise <- concordance/marg.cif
## caselower <- conclower/marg.cif
## caseup <- concup/marg.cif
## } else if (attr(object,"Type")=="RR") {
## casewise<- exp(coefs[,1])*c(marg.cif)
## concordance <- exp(coefs[,1])*marg.cif^2
## caselower <- marg.cif*exp(coefs[,1]-1.96*coefs[,2])
## caseup <- marg.cif*exp(coefs[,1]+1.96*coefs[,2])
## conclower <- marg.cif^2* exp(coefs[,1]-1.96*coefs[,2])
## concup <- marg.cif^2*exp(coefs[,1]+1.96*coefs[,2])
## } else if (attr(object,"Type")=="OR-cif") {
## thetal <- coefs[,1]-1.96*coefs[,2]
## thetau <- coefs[,1]+1.96*coefs[,2]
## casewise<- plack.cif2(marg.cif,marg.cif,c(coefs[,1]))/marg.cif
## concordance <- plack.cif2(marg.cif,marg.cif,c(coefs[,1]))
## caselower <- plack.cif2(marg.cif,marg.cif,thetal)/marg.cif
## caseup <- plack.cif2(marg.cif,marg.cif,thetau)/marg.cif
## conclower <- plack.cif2(marg.cif,marg.cif,thetal)
## concup <- plack.cif2(marg.cif,marg.cif,thetau)
## }
## outcase <- cbind(casewise,caselower,caseup)
## outconc <- cbind(concordance,conclower,concup)
## rownames(outcase) <- rownames(outconc) <- rownames(coefs)
## colnames(outcase) <- c("casewise concordance","2.5 %","97.5%")
## colnames(outconc) <- c("concordance","2.5 %","97.5%")
## }
## ## }}}
## res <- list(casewise=outcase,concordance=outconc,estimates=coefs,marg=marg.cif,type=attr(object,"Type"),sym=attr(object,"sym"),cause1=attr(object,"cause1"),cause2=attr(object,"cause2"))
## class(res) <- "summary.cor"
## res
## } ## }}}
## coef.cor<-function(object,...)
## { ## {{{
## res <- cbind(object$theta, diag(object$var.theta)^0.5)
## se<-diag(object$var.theta)^0.5
## wald <- object$theta/se
## waldp <- (1 - pnorm(abs(wald))) * 2
## cor<-exp(object$theta)
## res <- as.matrix(cbind(res, wald, waldp,cor,se*cor))
## if (attr(object,"Type")=="cor")
## colnames(res) <- c("log-Coef.", "SE", "z", "P-val","Cross odds ratio","SE")
## else colnames(res) <- c("log-ratio Coef.", "SE", "z", "P-val","Ratio","SE")
## if (is.null((rownames(res)))==TRUE) rownames(res)<-rep(" ",nrow(res))
## return(res)
## } ## }}}
Try the mets package in your browser
Any scripts or data that you put into this service are public.
mets documentation built on Jan. 17, 2023, 5:12 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions predictPairPlack plack.cif2 plack.cif concordanceCor concordance.cor print.cor coef.cor summary.cor print.summary.cor Grandom.cif random.cif or.cif rr.cif cor.cif dep.cif
Documented in concordance.cor concordanceCor cor.cif Grandom.cif or.cif plack.cif plack.cif2 predictPairPlack random.cif rr.cif summary.cor
dep.cif<-function(cif,data,cause=NULL,model="OR",cif2=NULL,times=NULL,
cause1=1,cause2=1,cens.code=NULL,cens.model="KM",Nit=40,tol=1e-6,detail=0,
clusters=NULL,theta=NULL,theta.des=NULL,step=1,sym=1,weights=NULL,
same.cens=FALSE,censoring.weights=NULL,silent=1,entry=NULL,estimator=1,
trunkp=1,admin.cens=NULL,control=list(),par.func=NULL,dpar.func=NULL,dimpar=NULL,
score.method="nlminb",random.design=NULL,var.link=0,...)
{ ## {{{
## {{{ set up data and design
multi<-0; dscore=1; stab.cens<-FALSE; entry.call<-entry; inverse<-var.link
notaylor<-1; flex.func<-0;
## extract design and time and cause from cif object
time <- cif$response[,"exit"]
if (is.null(cause)) cause <- cif$response[,"cause"] ## attr(cif,"cause");
cause <- as.numeric(cause)
if (is.null(cens.code)) cens.code <- attr(cif,"cens.code")
delta<-(cause!=cens.code)
if (length(cause)!=length(time)) stop("cause and time not of same length\n");
formula <- attr(cif,"Formula")
ldata <- timereg::aalen.des2(formula(delete.response(terms(formula))),data=data,model="aalen")
X<-ldata$X; Z<-ldata$Z;
antpers<-nrow(X);
if (is.null(Z)==TRUE) {npar<-TRUE; semi<-0;} else {Z<-as.matrix(Z); npar<-FALSE; semi<-1;}
if (npar==TRUE) {Z<-matrix(0,antpers,1); pg<-1; fixed<-0;} else {fixed<-1;pg<-ncol(Z);}
ng<-antpers;px<-ncol(X);
if (is.null(times)) times<-cif$cum[,1];
est<-as.matrix(cif$cum);
if (semi==1) gamma<-cif$gamma else gamma<-0;
ntimes<-length(times);
if ((cif$model!="additive") && (cif$model!="fg"))
stop("Marginal Cumulative incidence model, must be either additive or extended Fine-Gray model\n")
cif.model <- switch(cif$model,additive=1,fg=2)
if (cause1!=attr(cif,"causeS")) cat("Cause for marginal model and correlation not the same\n");
if ((cause1[1]!=cause2[1])) {
if (is.null(cif2)==TRUE) stop("Must provide marginal model for both causes");
formula2<-attr(cif2,"Formula");
ldata2 <- timereg::aalen.des2(formula(delete.response(terms(formula2))),data=data,model="aalen");
X2<-ldata2$X; Z2<-ldata$Z;
if (is.null(Z2)==TRUE) {npar2<-TRUE; semi2<-0;} else {Z2<-as.matrix(Z2); npar2<-FALSE; semi2<-1;}
if (npar2==TRUE) {Z2<-matrix(0,antpers,1); pg2<-1; fixed2<-0;} else {fixed2<-1;pg2<-ncol(Z2);}
px2<-ncol(X2);
est2<-cif2$cum;
if (semi2==1) gamma2<-cif2$gamma else gamma2<-0;
est2<-cpred(est2,times);
} else {
X2<-matrix(0,1,1); Z2<-matrix(0,1,1); pg2<-1; px2<-1; semi2<-0; est2<-matrix(0,1,2); gamma2<-0;
npar2<-FALSE;
}
### For truncation
if (is.null(entry)) entry.call <- NULL else entry.call <- 0
if (is.null(entry)) entry <- rep(0,antpers);
cum1<-cpred(rbind(rep(0,px+1),cif$cum),entry)[,-1];
if (cif.model==1) cif1entry <- 1-exp(-apply(X*cum1,1,sum)- (Z %*% gamma )*entry)
else if (cif.model==2) cif1entry <- 1-exp(-apply(X*cum1,1,sum)*exp(Z %*% gamma ))
if ((model!="COR") && (cause1[1]!=cause2[1])) {
cum2<-cpred(rbind(rep(0,px2+1),cif2$cum),entry)[,-1];
if (cif.model==1) cif2entry <- 1-exp(-apply(X2*cum2,1,sum)- (Z2 %*% gamma2 )*entry)
else if (cif.model==2) cif2entry <- 1-exp(-apply(X2*cum2,1,sum)*exp(Z2 %*% gamma2 ))
} else {cum2 <- cum1; cif2entry <- cif1entry;}
## }}}
## {{{ censoring model stuff
cens.weight <- cif$cens.weights ### censoring weights from cif function
Gcxe <- 1;
### when censoring probs given, wants them so cens.model="user.weights"
if (!is.null(censoring.weights)) cens.model <- "user.weights"
if (cens.model!="user.weights") {
if (is.null(cens.weight) ) {
if (cens.model=="KM") { ## {{{
ud.cens<-survival::survfit(Surv(time,cause==cens.code)~+1);
Gfit<-cbind(ud.cens$time,ud.cens$surv)
Gfit<-rbind(c(0,1),Gfit);
Gcx<-cpred(Gfit,time)[,2];
if (is.null(entry.call)==FALSE) Gcxe<-cpred(Gfit,entry)[,2];
Gcx <- Gcx/Gcxe;
Gctimes<- Gcx ## }}}
} else if (cens.model=="cox") { ## {{{
if (npar==TRUE) XZ<-X[,-1] else XZ<-cbind(X,Z)[,-1];
ud.cens<-timereg::cox.aalen(Surv(time,cause==cens.code)~prop(XZ),n.sim=0,robust=0);
Gcx<-cpred(ud.cens$cum,time)[,2];
RR<-exp(XZ %*% ud.cens$gamma)
Gcx<-exp(-Gcx*RR)
Gfit<-rbind(c(0,1),cbind(time,Gcx));
if (is.null(entry.call)==FALSE) {
Gcxe<-cpred(Gfit,entry)[,2];
Gcxe<-exp(-Gcxe*RR)
}
Gcx <- Gcx/Gcxe;
Gctimes<- Gcx ## }}}
} else if (cens.model=="aalen") { ## {{{
if (npar==TRUE) XZ<-X[,-1] else XZ<-cbind(X,Z)[,-1];
ud.cens<-timereg::aalen(Surv(time,cause==cens.code)~XZ,n.sim=0,robust=0);
Gcx<-cpred(ud.cens$cum,time)[,-1];
XZ<-cbind(1,XZ);
Gcx<-exp(-apply(Gcx*XZ,1,sum))
Gcx[Gcx>1]<-1; Gcx[Gcx<0]<-0
Gfit<-rbind(c(0,1),cbind(time,Gcx));
if (is.null(entry.call)==FALSE) {
Gcxe<-cpred(Gfit,entry)[,-1];
Gcxe<-exp(-apply(Gcxe*XZ,1,sum))
Gcxe[Gcxe>1]<-1; Gcxe[Gcxe<0]<-0
}
Gcx <- Gcx/Gcxe;
Gctimes<- Gcx ## }}}
}
} else { Gcx <- cens.weight; Gctimes <- cens.weight;}
} else {
Gcx <- censoring.weights
Gctimes <- cens.weight
}
ntimes<-length(times);
## }}}
## {{{ set up cluster + theta design + define iid variables
if (is.null(clusters)== TRUE) {
## take clusters from cif model
clusters <- attr(cif,"clusters");
antclust<- length(unique(clusters));
max.clust <- attr(cif,"max.clust")
if (attr(cif,"coarse.clust"))
stop("Max.clust should be NULL in marginal model, or clusters should be given at call \n");
} else {
clus<-unique(clusters); antclust<-length(clus);
clusters <- as.integer(factor(clusters, labels = 1:(antclust)))-1;
}
outc <- cluster.index(clusters,index.type=TRUE);
clustsize <- outc$cluster.size
maxclust <- outc$maxclust
clusterindex <- outc$idclust
if (maxclust==1) stop("No clusters given \n");
if (model!="ARANCIF") {
if (is.null(theta.des)==TRUE) { ptheta<-1; theta.des<-matrix(1,antpers,ptheta);
colnames(theta.des) <- "intercept"; } else theta.des<-as.matrix(theta.des);
ptheta<-ncol(theta.des);
if (is.null(dimpar) & is.null(par.func)) dimpar<-ptheta;
}
if (model=="ARANCIF") { ### different parameters for Additive random effects
if (is.null(random.design)) random.design <- matrix(1,antpers,1);
dim.rv <- ncol(random.design);
if (is.null(theta.des)==TRUE) theta.des<-diag(dim.rv);
dimpar <- ncol(theta.des);
if (nrow(theta.des)!=ncol(random.design)) stop("nrow(theta.des)!= ncol(random.design)");
} else random.design <- matrix(0,1,1);
if ( (!is.null(par.func)) && is.null(dimpar) )
stop("Must specify dimension of parameters when specifying R-functions\n")
if ( (is.null(theta)==TRUE) && (model=="ARANCIF")) theta<-rep(0.5,dimpar);
if (is.null(theta)==TRUE) theta<-rep(0.1,dimpar);
if (length(theta)!=dimpar) theta<-rep(theta[1],dimpar);
Biid<-c(); gamma.iid <- 0; B2iid<-c(); gamma2.iid <- 0;
if (notaylor==0) {
for (i in 1:antclust) Biid<-cbind(Biid,cif$B.iid[[i]]);
if (npar==TRUE) gamma.iid<-0 else gamma.iid<-cif$gamma.iid;
if ((model!="COR") && (cause1!=cause2)) {
B2iid<-c()
for (i in 1:antclust) B2iid<-cbind(B2iid,cif2$B.iid[[i]]);
if (npar2==TRUE) gamma2.iid<-0 else gamma2.iid<-cif2$gamma.iid;
} else { B2iid<-Biid; gamma2.iid<-gamma.iid; }
}
var.theta<-hess<-matrix(0,dimpar,dimpar); score<-rep(0,dimpar);
time.pow<-attr(cif,"time.pow");
#if (sum(time.pow)==0) time.pow<-rep(1,pg);
theta.iid<-matrix(0,antclust,dimpar);
if ((nrow(theta.des)!=nrow(X)) && (model!="ARANCIF"))
stop("Dependence design not consistent with data\n");
if (length(trunkp)!=antpers) trunkp <- rep(1,antpers)
if (is.null(weights)==TRUE) weights <- rep(1,antpers);
## }}}
### ## {{{ function and derivative
if (is.null(par.func)==FALSE)
{
flex.func<-1; # use flexible design
if (is.null(dpar.func)) {
cat("Must provide derivative that is used for iid decomposition for SE's\n");
score.method <- "nlminb"
dpar.func <- par.func
}
} ## }}}
Zgamma <- c(Z %*% gamma);
Z2gamma2 <- c(Z2 %*% gamma2);
dep.model <- 0;
dep.model <- switch(model,COR=1,RR=2,OR=3,RANCIF=4,ARANCIF=5)
if (dep.model==0) stop("model must be COR, OR, RR, RANCIF, ARANCIF \n");
if (dep.model<=4) rvdes <- matrix(0,1,1);
obj <- function(par)
{ ## {{{
if (is.null(par.func)==TRUE) {
Xtheta <- theta.des %*% par; DXtheta <- array(0,c(1,1,1));
} else {
Xtheta <- c(); DXtheta <- array(0,c(length(times),antpers,dimpar));
if (is.null(dpar.func))
stop("Must also provide derivative of function wrt to parameters")
s <- 0
for (t in times) {
s <- 1+s
Xttheta <- par.func(par,t,theta.des)
if (length(Xttheta)!=antpers) stop("parfunc(par,t,theta.des) must length n");
Xtheta <- cbind(Xtheta,Xttheta)
Dttheta <- dpar.func(par,t,theta.des);
if (dim(Dttheta)[1]!=antpers || dim(Dttheta)[2]!=dimpar)
stop("dparfunc must return matrix n x dimpar when called on dparfunc(par,t,theta.des)");
DXtheta[s,,] <- Dttheta
}
}
outl<-.Call("cor", ## {{{
itimes=times,iy=time,icause=cause,iCA1=cause1,iKMc=Gcx,
iz=X,iest=matrix(est[,-1],ncol=ncol(est)-1),iZgamma=c(Zgamma),isemi=semi,izsem=Z,
itheta=c(par),iXtheta=Xtheta,iDXtheta=DXtheta,idimDX=dim(DXtheta),
ithetades=theta.des,
icluster=clusters,iclustsize=clustsize,iclusterindex=clusterindex,
iinverse=inverse,iCA2=cause2,
ix2=X2,isemi2=semi2,iest2=as.matrix(est2[,-1]),iZgamma2=c(Z2gamma2),
iflexfunc=flex.func,iiid=iid,isym=sym,iweights=weights,
isamecens=as.numeric(same.cens),istabcens=as.numeric(stab.cens),
iKMtimes=Gctimes,isilent=silent,
icifmodel=cif.model,idepmodel=dep.model,
iestimator=estimator,ientryage=entry,icif1entry=cif1entry,
icif2entry=cif2entry,itrunkp=trunkp,irvdes=random.design,
## Biid,gamma.iid,time.pow,
## B2iid,gamma2.iid,body(htheta),body(dhtheta),new.env(),
PACKAGE="mets"
)
## }}}
attr(outl,"gradient") <-outl$score
if (oout==0) return(outl$ssf) else if (oout==1) return(sum(outl$score^2)) else return(outl)
} ## }}}
p <- theta
iid <- 0; ###no-iid representation for iterations
if (score.method=="nr") { ## {{{
oout <- 2; ### output control for obj
if (Nit>0)
for (i in 1:Nit)
{
oout <- 2
out <- obj(p)
hess <- out$Dscore
if (!is.na(sum(hess))) hessi <- lava::Inverse(out$Dscore) else hessi <- hess
if (detail==1) {## {{{
print(paste("Fisher-Scoring ===================: it=",i));
cat("theta:");print(c(p))
cat("score:");print(c(out$score));
cat("hess:"); print(hess);
}## }}}
delta <- hessi %*% out$score
## do not update last iteration
if (i<Nit) p <- p-delta* step
if (is.nan(sum(out$score))) break;
if (sum(abs(out$score))<tol) break;
if (max(delta)>20) { cat("NR increment > 20, lower step zize, increment= \n"); cat(delta); break; }
}
if (!is.nan(sum(p))) { ## {{{ iid decomposition
oout <- 2
theta <- p
iid <- 1;
out <- obj(p)
score <- out$score
hess <- out$Dscore
} ## }}}
if (detail==1 & Nit==0) {## {{{
print(paste("Fisher-Scoring ===================: final"));
cat("theta:");print(c(p))
cat("score:");print(c(out$score));
cat("hess:"); print(hess);
}## }}}
if (!is.na(sum(hess))) hessi <- lava::Inverse(out$Dscore) else hessi <- diag(nrow(hess))
score1 <- score;
## }}}
} else if (score.method=="nlminb") { ## {{{ nlminb optimizer
iid <- 0; oout <- 0;
tryCatch(opt <- nlminb(theta,obj,control=control),error=function(x) NA)
if (detail==1) print(opt);
iid <- 1;
hess <- numDeriv::hessian(obj,opt$par)
score <- numDeriv::jacobian(obj,opt$par)
hessi <- lava::Inverse(hess);
theta <- opt$par
if (detail==1) cat("iid decomposition\n");
oout <- 2;
out <- obj(opt$par)
score1 <- out$score
## }}}
} else if (score.method=="nlm") { ## {{{ nlm optimizer
iid <- 0; oout <- 0;
tryCatch(opt <- nlm(obj,theta,hessian=TRUE,print.level=detail),error=function(x) NA)
iid <- 1;
hess <- opt$hessian
score <- opt$gradient
if (detail==1) print(opt);
hessi <- lava::Inverse(hess);
theta <- opt$estimate
if (detail==1) cat("iid decomposition\n");
oout <- 2;
out <- obj(opt$estimate)
score1 <- out$score
## }}}
} else stop("score.methods = nlm nlminb nr\n");
theta.iid <- out$theta.iid %*% hessi
if (is.null(par.func)) var.theta <- t(theta.iid) %*% theta.iid else var.theta <- hessi
var.theta <- t(theta.iid) %*% theta.iid
if (is.null(par.func)) thetanames <- colnames(theta.des) else
thetanames <- rep("R-func",dimpar)
ud <- list(theta=theta,score=score,hess=hess,hessi=hessi,var.theta=var.theta,
theta.iid=theta.iid,score1=c(score1),thetanames=thetanames,
brierscore=out$ssf,p11=out$p11);
if (dep.model<=3) class(ud)<-"cor"
else if (dep.model==4) class(ud) <- "randomcif"
else if (dep.model==6) class(ud) <- "randomcif"
else if (dep.model==5) class(ud) <- "randomcifrv"
attr(ud, "Formula") <- formula
attr(ud, "Clusters") <- clusters
attr(ud,"cause1")<-cause1; attr(ud,"cause2")<-cause2
attr(ud,"sym")<-sym;
attr(ud,"inverse")<-inverse;
attr(ud,"antpers")<-antpers;
attr(ud,"antclust")<-antclust;
attr(ud,"var.link")<-var.link;
if (dep.model==4) attr(ud, "Type") <- "randomcif"
if (dep.model==6) attr(ud, "Type") <- "randomcif"
if (model=="COR") attr(ud, "Type") <- "cor"
if (model=="RR") attr(ud, "Type") <- "RR"
if (model=="OR") attr(ud, "Type") <- "OR-cif"
if (dep.model==5) attr(ud, "pardes") <- theta.des
if (dep.model==5) attr(ud, "rv1") <- random.design[1,]
return(ud);
} ## }}}
##' Cross-odds-ratio, OR or RR risk regression for competing risks
##'
##' Fits a parametric model for the log-cross-odds-ratio for the
##' predictive effect of for the cumulative incidence curves for \eqn{T_1}
##' experiencing cause i given that \eqn{T_2} has experienced a cause k :
##' \deqn{
##' \log(COR(i|k)) = h(\theta,z_1,i,z_2,k,t)=_{default} \theta^T z =
##' }
##' with the log cross odds ratio being
##' \deqn{
##' 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
##' \deqn{
##' 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
##' \deqn{
##' 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 \eqn{B^c} is the complement event of \eqn{B},
##' \eqn{P_x} is the distribution given covariates
##' (\eqn{x} are subject specific and \eqn{z} are cluster specific covariates), and
##' \eqn{h()} is a function that is the simple identity
##' \eqn{\theta^T z} by default.
##'
##' The OR dependence measure is given by
##' \deqn{
##' 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
##' \deqn{
##' 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 \eqn{T_1} and \eqn{T_2} can be interchanged and leads to
##' the same cross-odd-ratio (i.e.
##' \eqn{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).
##'
##' \eqn{h()} may be specified as an R-function of the parameters,
##' see example below, but the default is that it is simply \eqn{\theta^T z}.
##'
##' @aliases or.cif rr.cif
##' @param 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
##' @param data a data.frame with the variables.
##' @param cause specifies the causes related to the death
##' times, the value cens.code is the censoring value. When missing it comes from marginal cif.
##' @param times time-vector that specifies the times used for the estimating euqations for the cross-odds-ratio estimation.
##' @param cause1 specificies the cause considered.
##' @param cause2 specificies the cause that is conditioned on.
##' @param cens.code specificies the code for the censoring if NULL then uses the one from the marginal cif model.
##' @param cens.model specified which model to use for the ICPW, KM is Kaplan-Meier alternatively it may be "cox"
##' @param Nit number of iterations for Newton-Raphson algorithm.
##' @param detail if 0 no details are printed during iterations, if 1 details are given.
##' @param clusters specifies the cluster structure.
##' @param theta specifies starting values for the cross-odds-ratio parameters of the model.
##' @param theta.des specifies a regression design for the cross-odds-ratio parameters.
##' @param step specifies the step size for the Newton-Raphson algorithm.
##' @param sym specifies if symmetry is used in the model.
##' @param weights weights for estimating equations.
##' @param par.func parfunc
##' @param dpar.func dparfunc
##' @param dimpar dimpar
##' @param score.method "nlminb", can also use "nr".
##' @param same.cens if true then censoring within clusters are assumed to be the same variable, default is independent censoring.
##' @param censoring.weights these probabilities are used for the bivariate censoring dist.
##' @param silent 1 to suppress output about convergence related issues.
##' @param ... Not used.
##' @return returns an object of type 'cor'. With the following arguments:
##' \item{theta}{estimate of proportional odds parameters of model.}
##' \item{var.theta}{variance for gamma. }
##' \item{hess}{the derivative of the used score.}
##' \item{score}{scores at final stage.}
##' \item{score}{scores at final stage.}
##' \item{theta.iid}{matrix of iid decomposition of parametric effects.}
##' @author 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
##' ###########################################################
##' \donttest{ ## 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)
##' }
##' @export
##' @keywords survival
cor.cif<-function(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,...)
{ ## {{{
fit <- dep.cif(cif=cif,data=data,cause=cause,model="COR",times=times,
cause1=cause1,cause2=cause2,cens.code=cens.code,cens.model=cens.model,Nit=Nit,detail=detail,
clusters=clusters,theta=theta,theta.des=theta.des,par.func=par.func,dpar.func=dpar.func,
dimpar=dimpar,
step=step,sym=sym,weights=weights,
score.method=score.method,same.cens=same.cens,censoring.weights=censoring.weights,silent=silent,...)
fit$call <- match.call()
fit
} ## }}}
##' @export
rr.cif<-function(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, step=1,sym=0,weights=NULL,
same.cens=FALSE,censoring.weights=NULL,silent=1,par.func=NULL,dpar.func=NULL,dimpar=NULL,
score.method="nlminb",entry=NULL,estimator=1,trunkp=1,admin.cens=NULL,...)
{ ## {{{
fit <- dep.cif(cif=cif,data=data,cause=cause,model="RR",cif2=cif2,times=times,
cause1=cause1,cause2=cause2,cens.code=cens.code,cens.model=cens.model,Nit=Nit,detail=detail,
clusters=clusters,theta=theta,theta.des=theta.des, step=step,sym=sym,weights=weights,
same.cens=same.cens,censoring.weights=censoring.weights,silent=silent,
par.func=par.func,dpar.func=dpar.func,dimpar=dimpar, score.method=score.method,
entry=entry,estimator=estimator,trunkp=trunkp,admin.cens=admin.cens,...)
fit$call <- match.call()
fit
} ## }}}
##' @export
or.cif<-function(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, step=1,sym=0, weights=NULL,
same.cens=FALSE,censoring.weights=NULL,silent=1,par.func=NULL,dpar.func=NULL,dimpar=NULL,
score.method="nlminb",entry=NULL,estimator=1,trunkp=1,admin.cens=NULL,...)
{ ## {{{
fit <- dep.cif(cif=cif,data=data,cause=cause,model="OR",cif2=cif2,times=times,
cause1=cause1,cause2=cause2,cens.code=cens.code,cens.model=cens.model,Nit=Nit,detail=detail,
clusters=clusters,theta=theta,theta.des=theta.des, step=step,sym=sym,weights=weights,
same.cens=same.cens,censoring.weights=censoring.weights,silent=silent,
par.func=par.func,dpar.func=dpar.func,dimpar=dimpar,
score.method=score.method,
entry=entry,estimator=estimator,trunkp=trunkp,admin.cens=admin.cens,...)
fit$call <- match.call()
fit
} ## }}}
##' 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 the form
##' \deqn{
##' P(T \leq t, cause=1 | x,z) = P_1(t,x,z) = 1- exp( -x^T A(t) exp(z^T \beta))
##' }
##' We allow a regression structure for the random effects variances that may depend on
##' cluster covariates.
##'
##' @title Random effects model for competing risks data
##' @param cif a model object from the comp.risk function with the
##' marginal cumulative incidence of cause2, i.e., the event that is conditioned on, and whose
##' odds the comparision is made with respect to
##' @param data a data.frame with the variables.
##' @param cause specifies the causes related to the death
##' times, the value cens.code is the censoring value.
##' @param cif2 specificies model for cause2 if different from cause1.
##' @param cause1 cause of first coordinate.
##' @param cause2 cause of second coordinate.
##' @param cens.code specificies the code for the censoring if NULL then uses the one from the marginal cif model.
##' @param cens.model specified which model to use for the ICPW, KM is Kaplan-Meier alternatively it may be "cox"
##' @param Nit number of iterations for Newton-Raphson algorithm.
##' @param detail if 0 no details are printed during iterations, if 1 details are given.
##' @param clusters specifies the cluster structure.
##' @param theta specifies starting values for the cross-odds-ratio parameters of the model.
##' @param theta.des specifies a regression design for the cross-odds-ratio parameters.
##' @param sym 1 for symmetry 0 otherwise
##' @param step specifies the step size for the Newton-Raphson algorith.m
##' @param same.cens if true then censoring within clusters are assumed to be the same variable, default is independent censoring.
##' @param var.link if var.link=1 then var is on log-scale.
##' @param score.method default uses "nlminb" optimzer, alternatively, use the "nr" algorithm.
##' @param entry entry-age in case of delayed entry. Then two causes must be given.
##' @param trunkp gives probability of survival for delayed entry, and related to entry-ages given above.
##' @param ... extra arguments.
##' @return returns an object of type 'cor'. With the following arguments:
##' \item{theta}{estimate of proportional odds parameters of model.}
##' \item{var.theta}{variance for gamma. }
##' \item{hess}{the derivative of the used score.}
##' \item{score}{scores at final stage.}
##' \item{score}{scores at final stage.}
##' \item{theta.iid}{matrix of iid decomposition of parametric effects.}
##' @export
##' @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 (2012), work in progress.
##' @examples
##' \donttest{ ## Reduce Ex.Timings
##' d <- simnordic.random(5000,delayed=TRUE,cordz=0.5,cormz=2,lam0=0.3,country=TRUE)
##' times <- seq(50,90,by=10)
##' add1 <- 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)
##'
##' out1<-random.cif(add1,data=d,cause1=1,cause2=1,theta=1,same.cens=TRUE)
##' summary(out1)
##'
##' out2<-random.cif(add1,data=d,cause1=1,cause2=1,theta=1,
##' theta.des=mm,same.cens=TRUE)
##' summary(out2)
##'
##' #########################################
##' ##### 2 different causes
##' #########################################
##'
##' add2 <- timereg::comp.risk(Event(time,cause)~-1+factor(country)+cluster(id),data=d,
##' times=times,cause=2,max.clust=NULL)
##' out3 <- random.cif(add1,data=d,cause1=1,cause2=2,cif2=add2,sym=1,same.cens=TRUE)
##' summary(out3) ## negative dependence
##'
##' out4 <- random.cif(add1,data=d,cause1=1,cause2=2,cif2=add2,theta.des=mm,sym=1,same.cens=TRUE)
##' summary(out4) ## negative dependence
##' }
##' @keywords survival
##' @author Thomas Scheike
random.cif<-function(cif,data,cause=NULL,cif2=NULL,
cause1=1,cause2=1,cens.code=NULL,cens.model="KM",Nit=40,detail=0,
clusters=NULL,theta=NULL,theta.des=NULL,sym=1,
step=1,same.cens=FALSE,var.link=0,score.method="nr",
entry=NULL,trunkp=1,...)
{ ## {{{
fit <- dep.cif(cif,data=data,cause=cause,model="RANCIF",cif2=cif2,sym=sym,
cause1=cause1,cause2=cause2,cens.code=cens.code,cens.model=cens.model,Nit=Nit,detail=detail,
clusters=clusters,theta=theta,theta.des=theta.des,step=step,same.cens=same.cens,
var.link=var.link,score.method=score.method,entry=entry,trunkp=trunkp,...)
fit$call <- match.call()
fit
} ## }}}
##'Additive Random effects model for competing risks data for polygenetic modelling
##'
##'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
##'\deqn{
##'P(T \leq t, cause=1 | x,z) = P_1(t,x,z) = 1- exp( -x^T A(t) - t z^T \beta)
##'}
##'
##'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
##' \eqn{v_1} and \eqn{v_2}.
##'Such that the random effects for subject
##'1 is \deqn{v_1^T (Z_1,...,Z_d)}, for d random effects. Each random effect
##'has an associated parameter \eqn{(\lambda_1,...,\lambda_d)}. By construction
##'subjects 1's random effect are Gamma distributed with
##'mean \eqn{\lambda_1/v_1^T \lambda}
##'and variance \eqn{\lambda_1/(v_1^T \lambda)^2}. Note that the random effect
##'\eqn{v_1^T (Z_1,...,Z_d)} has mean 1 and variance \eqn{1/(v_1^T \lambda)}.
##'
##'The parameters \eqn{(\lambda_1,...,\lambda_d)}
##'are related to the parameters of the model
##'by a regression construction \eqn{pard} (d x k), that links the \eqn{d}
##'\eqn{\lambda} parameters
##'with the (k) underlying \eqn{\theta} parameters
##'\deqn{
##'\lambda = pard \theta
##'}
##'
##'
##'
##' @param cif a model object from the timereg::comp.risk function with the
##' marginal cumulative incidence of cause2, i.e., the event that is conditioned on, and whose
##' odds the comparision is made with respect to
##' @param data a data.frame with the variables.
##' @param cause specifies the causes related to the death
##' times, the value cens.code is the censoring value.
##' @param cif2 specificies model for cause2 if different from cause1.
##' @param times time points
##' @param cause1 cause of first coordinate.
##' @param cause2 cause of second coordinate.
##' @param cens.code specificies the code for the censoring if NULL then uses the one from the marginal cif model.
##' @param cens.model specified which model to use for the ICPW, KM is Kaplan-Meier alternatively it may be "cox"
##' @param Nit number of iterations for Newton-Raphson algorithm.
##' @param detail if 0 no details are printed during iterations, if 1 details are given.
##' @param clusters specifies the cluster structure.
##' @param theta specifies starting values for the cross-odds-ratio parameters of the model.
##' @param theta.des specifies a regression design for the cross-odds-ratio parameters.
##' @param step specifies the step size for the Newton-Raphson algorith.m
##' @param sym 1 for symmetri and 0 otherwise
##' @param weights weights for score equations.
##' @param same.cens if true then censoring within clusters are assumed to be the same variable, default is independent censoring.
##' @param censoring.weights Censoring probabilities
##' @param silent debug information
##' @param var.link if var.link=1 then var is on log-scale.
##' @param score.method default uses "nlminb" optimzer, alternatively, use the "nr" algorithm.
##' @param entry entry-age in case of delayed entry. Then two causes must be given.
##' @param estimator estimator
##' @param trunkp gives probability of survival for delayed entry, and related to entry-ages given above.
##' @param admin.cens Administrative censoring
##' @param random.design specifies a regression design of 0/1's for the random effects.
##' @param ... extra arguments.
##' @return returns an object of type 'random.cif'. With the following arguments:
##'\item{theta}{estimate of parameters of model.}
##'\item{var.theta}{variance for gamma. }
##'\item{hess}{the derivative of the used score.}
##'\item{score}{scores at final stage.}
##'\item{theta.iid}{matrix of iid decomposition of parametric effects.}
##' @export
##' @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
##' \donttest{ ## Reduce Ex.Timings
##' d <- simnordic.random(5000,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)
##'
##' }
##' @keywords survival
##' @author Thomas Scheike
##' @export
Grandom.cif<-function(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,...)
{ ## {{{
fit <- dep.cif(cif,data=data,cause=cause,model="ARANCIF",cif2=cif2,times=times,
cause1=cause1,cause2=cause2,cens.code=cens.code,cens.model=cens.model,Nit=Nit,detail=detail,
clusters=clusters,theta=theta,theta.des=theta.des,step=step,sym=sym,weights=weights,
same.cens=same.cens,censoring.weights=censoring.weights,silent=silent,var.link=var.link,
score.method=score.method,entry=entry,estimator=estimator,
random.design=random.design,trunkp=trunkp,admin.cens=admin.cens,...)
fit$call <- match.call()
fit
} ## }}}
##' @export
print.summary.cor <- function(x,digits=3,...)
{ ## {{{
if (x$type=="cor") {
cat("Cross odds ratio dependence for competing risks\n\n")
cat("Effect of cause1=",x$cause1," on cause2=",x$cause2,
" under symmetry=",x$sym,fill=TRUE,sep="")
} else if (x$type=="RR") {
cat("Ratio of joint and product of marginals for competing risks\n\n")
cat("Ratio of cumulative incidence for cause1=",x$cause1," and cause2=",x$cause2,sep=" ")
} else if (x$type=="OR-cif") {
cat("OR for dependence for competing risks\n\n")
cat("OR of cumulative incidence for cause1=",x$cause1," and cause2=",x$cause2,sep=" ")
}
cat("\n")
prmatrix(signif(x$estimates,digits))
cat("\n")
if (!is.null(x$marg)) {
cat(paste("Marginal cumulative incidencen",signif(x$marg,digits),"\n"))
prmatrix(signif(x$casewise,digits))
prmatrix(signif(x$concordance,digits))
cat("\n")
}
invisible(x)
} ## }}}
##' Computes concordance and casewise concordance for dependence models for competing risks
##' models of the type cor.cif, rr.cif or or.cif for the given cumulative incidences and the different dependence
##' measures in the object.
##'
##'
##' @title Summary for dependence models for competing risks
##' @param object object from cor.cif rr.cif or or.cif for dependence between competing risks data for two causes.
##' @param marg.cif a number that gives the cumulative incidence in one time point for which concordance and
##' casewise concordance are computed.
##' @param marg.cif2 the cumulative incidence for cause 2 for concordance and
##' casewise concordance are computed. Default is that it is the same as marg.cif.
##' @param digits digits in output.
##' @param ... Additional arguments.
##' @return prints summary for dependence model.
##' \item{casewise}{gives casewise concordance that is, probability of cause 2 (related to cif2) given that cause 1 (related to cif1)
##' has occured.}
##' \item{concordance}{gives concordance that is, probability of cause 2 (related to cif2) and cause 1 (related to cif1).}
##' \item{cif1}{cumulative incidence for cause1.}
##' \item{cif2}{cumulative incidence for cause1.}
##' @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.
##'
##' @author Thomas Scheike
##' @keywords survival
##' @examples
##' ## library("timereg")
##' ## data("multcif",package="mets") # simulated data
##' ## multcif$cause[multcif$cause==0] <- 2
##' ##
##' ## times=seq(0.1,3,by=0.1) # to speed up computations use only these time-points
##' ## add <- timereg::comp.risk(Event(time,cause)~+1+cluster(id),
##' ## data=multcif,n.sim=0,times=times,cause=1)
##' ###
##' ## out1<-cor.cif(add,data=multcif,cause1=1,cause2=1,theta=log(2+1))
##' ## summary(out1)
##' ##
##' ## pad <- predict(add,X=1,se=0,uniform=0)
##' ## summary(out1,marg.cif=pad)
##' @method summary cor
##' @export
summary.cor <- function(object,marg.cif=NULL,marg.cif2=NULL,digits=3,...) { ## {{{
if (!(attr(object,"class") %in% c("cor","randomcif"))) stop("Must be a cor.cif or randomcif object")
if (sum(abs(object$score))>0.001) warning("WARNING: check score for convergence\n")
coefs <- coef(object,...)
ocasewise <- oconcordance <- NULL
if (is.null(marg.cif)==FALSE) { ## {{{
time <- marg.cif$time
marg.cif <- marg.cif$P1
if (is.null(marg.cif2)==FALSE) marg.cif2 <- marg.cif2$P1 else {
if (attr(object,"cause2")==attr(object,"cause1")) marg.cif2 <- marg.cif
else stop("causes not the same and second marginal cif not given\n");
}
pmarg.cif <- marg.cif*marg.cif2
thetav <- coefs[,1]
thetavl <- thetav-1.96*coefs[,2]
thetavu <- thetav+1.96*coefs[,2]
namev <- object$thetanames
if (is.null(namev)) namev <- paste("name",1:length(thetav),sep="")
ocasewise <- oconcordance <- list()
k <- 0
for (theta in thetav) {
k <- k+1
thetal <- thetavl[k];
thetau <- thetavu[k]
if (attr(object,"Type")=="cor") { ## {{{
concordance <- exp(theta)*pmarg.cif/((1-marg.cif)+exp(theta)*marg.cif)
conclower <- exp(thetal)*pmarg.cif/((1-marg.cif)+exp(thetal)*marg.cif)
concup <- exp(thetau)*pmarg.cif/((1-marg.cif)+exp(thetau)*marg.cif)
casewise <- concordance/marg.cif
caselower <- conclower/marg.cif
caseup <- concup/marg.cif
} else if (attr(object,"Type")=="RR") {
casewise<- exp(coefs[,1])*c(marg.cif)
concordance <- exp(coefs[,1])*pmarg.cif
caselower <- marg.cif*exp(thetal)
caseup <- marg.cif*exp(thetau)
conclower <- pmarg.cif* exp(thetal)
concup <- pmarg.cif*exp(thetau)
} else if (attr(object,"Type")=="OR-cif") {
casewise<- plack.cif2(marg.cif,marg.cif,c(theta))/marg.cif
concordance <- plack.cif2(marg.cif,marg.cif,c(theta))
caselower <- plack.cif2(marg.cif,marg.cif,thetal)/marg.cif
caseup <- plack.cif2(marg.cif,marg.cif,thetau)/marg.cif
conclower <- plack.cif2(marg.cif,marg.cif,thetal)
concup <- plack.cif2(marg.cif,marg.cif,thetau)
} else if (attr(object,"Type")=="randomcif") {
theta <- 1/theta
thetal <- 1/thetal
thetau <- 1/thetau
lam <- 1-marg.cif
p11<- 1-lam -lam +lap(theta,2*ilap(theta, lam))
concordance <- p11
conclower <- 1-lam -lam +lap(thetal,2*ilap(thetal, lam))
concup <- 1-lam -lam +lap(thetau,2*ilap(thetau, lam))
casewise <- concordance/marg.cif
caselower <- conclower/marg.cif
caseup <- concup/marg.cif
} ## }}}
outcase <- cbind(time,casewise,caselower,caseup)
outconc <- cbind(time,concordance,conclower,concup)
colnames(outcase) <- c("time","casewise concordance","2.5 %","97.5%")
colnames(outconc) <- c("time","concordance","2.5 %","97.5%")
if (length(thetav)==1) {
ocasewise <- outcase
oconcordance <- outconc
} else {
ocasewise[[k]] <- outcase
oconcordance[[k]] <- outconc
}
}
if (length(thetav)>1) {
if (length(ocasewise)==length(namev)) {
names(ocasewise) <- namev
names(oconcordance) <- namev
}
}
} ## }}}
res <- list(casewise=ocasewise,concordance=oconcordance,estimates=coefs,
marg.cif=marg.cif, marg.cif2=marg.cif2,type=attr(object,"Type"),
sym=attr(object,"sym"),cause1=attr(object,"cause1"),cause2=attr(object,"cause2"))
class(res) <- "summary.cor"
res
} ## }}}
##' @export
coef.cor <- function(object,...)
{ ## {{{
res <- cbind(object$theta, diag(object$var.theta)^0.5)
se<-diag(object$var.theta)^0.5
wald <- object$theta/se
waldp <- (1 - pnorm(abs(wald))) * 2
cor<-exp(object$theta)
res <- as.matrix(cbind(res, wald, waldp,cor,se*cor))
if (attr(object,"Type")=="cor")
colnames(res) <- c("log-Coef.", "SE", "z", "P-val","Cross odds ratio","SE")
else colnames(res) <- c("log-ratio Coef.", "SE", "z", "P-val","Ratio","SE")
if (!is.null(object$thetanames)) rownames(res)<-object$thetanames
return(res)
} ## }}}
##' @export
print.cor<-function(x,digits=3,...)
{ ## {{{
print(x$call);
cat("\n")
print(summary(x));
} ## }}}
##' @export
concordance.cor <- function(object,...) {
concordanceCor(object,...)
}
##' Concordance for Twins
##'
##' The concordance is the probability that both twins have experienced the
##' event of interest and is defined as
##' \deqn{
##' cor(t) = P(T_1 \leq t, \epsilon_1 =1 , T_2 \leq t, \epsilon_2=1)
##' }
##'
##' Similarly, the casewise concordance is
##' \deqn{
##' 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.
##'
##' @title Concordance Computes concordance and casewise concordance
##' @param object Output from the cor.cif, rr.cif or or.cif function
##' @param cif1 Marginal cumulative incidence
##' @param cif2 Marginal cumulative incidence of other cause (cause2) if it is different from cause1
##' @param messages To print messages
##' @param model Specfifies wich model that is considered if object not given.
##' @param coefs Specfifies dependence parameters if object is not given.
##' @param ... Extra arguments, not used.
##' @author Thomas Scheike
##' @aliases concordanceCor concordance.cor
##' @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
##'
##' @export
concordanceCor <- function(object,cif1,cif2=NULL,messages=TRUE,model=NULL,coefs=NULL,...)
{ ## {{{
if (is.null(model)) {
if (!inherits(object, "cor")) stop("Must be a rr.cif, cor.cif or or.cif object")
model <- attr(object,"Type")
}
if (is.null(coefs)) coefs <- coef(object)
if (is.null(coefs)) stop("Must give dependence parameters\n");
if (!is.null(object)) {
cause1 <- attr(object,"cause1");
cause2 <- attr(object,"cause2");
} else cause1 <- cause2 <- 1
if (is.null(cif2)==TRUE) cif2 <- cif1;
if (messages) {
if (model=="cor") { ## {{{
message("Cross odds ratio dependence for competing risks\n\n")
message("Odds of cause1=",cause1," given cause2=",cause2," relative to Odds of cause1=",cause1,"\n",fill=TRUE,sep="")
} else if (model=="RR") {
message("Ratio of joint and product of marginals for competing risks\n\n")
message("Ratio of cumulative incidence for cause1=",cause1," and cause2=",cause2,sep=" ")
} else if (model=="OR-cif") {
message("OR for dependence for competing risks\n\n")
message("OR of cumulative incidence for cause1=",cause1," and cause2=",cause2,sep=" ")
} ## }}}
}
out <- list()
for (k in 1:nrow(coefs)) {
## {{{
if (model=="cor") {
concordance <- exp(coefs[k,1])*cif1*cif2/((1-cif1)+exp(coefs[k,1])*cif1)
conclower <- exp(coefs[k,1]-1.96*coefs[k,2])*cif1*cif2/((1-cif1)+exp(coefs[k,1]-1.96*coefs[k,2])*cif1)
concup <- exp(coefs[k,1]+1.96*coefs[k,2])*cif1*cif2/((1-cif1)+exp(coefs[k,1]+1.96*coefs[k,2])*cif1)
casewise <- concordance/cif1
caselower <- conclower/cif1
caseup <- concup/cif1
} else if (model=="RR") {
casewise<- exp(coefs[k,1])*c(cif2)
concordance <- exp(coefs[k,1])*cif1*cif2
caselower <- cif2*exp(coefs[k,1]-1.96*coefs[k,2])
caseup <- cif2*exp(coefs[k,1]+1.96*coefs[k,2])
conclower <- cif1*cif2* exp(coefs[k,1]-1.96*coefs[k,2])
concup <- cif1*cif2*exp(coefs[k,1]+1.96*coefs[k,2])
} else if (model=="OR-cif") {
thetal <- coefs[k,1]-1.96*coefs[k,2]
thetau <- coefs[k,1]+1.96*coefs[k,2]
casewise<- plack.cif2(cif1,cif2,c(coefs[k,1]))/cif1
concordance <- plack.cif2(cif1,cif2,c(coefs[k,1]))
caselower <- plack.cif2(cif1,cif2,thetal)/cif1
caseup <- plack.cif2(cif1,cif2,thetau)/cif1
conclower <- plack.cif2(cif1,cif2,thetal)
concup <- plack.cif2(cif1,cif2,thetau)
}
outcase <- cbind(c(casewise),c(caselower),c(caseup))
outconc <- cbind(c(concordance),c(conclower),c(concup))
colnames(outcase) <- c("casewise concordance","2.5 %","97.5%")
colnames(outconc) <- c("concordance","2.5 %","97.5%")
## }}}
out[[k]] <- list(concordance=outconc,casewise=outcase)
names(out)[k] <- rownames(coefs)[k]
}
return(out)
} ## }}}
##' .. content for description (no empty lines) ..
##'
##' @title plack Computes concordance for or.cif based model, that is Plackett random effects model
##' @aliases plack.cif2
##' @export
##' @param cif1 Cumulative incidence of first argument.
##' @param cif2 Cumulative incidence of second argument.
##' @param object or.cif object with dependence parameters.
##' @author Thomas Scheike
plack.cif <- function(cif1,cif2,object)
{ ## {{{
coefs <- coef(object)
theta <- exp(object$theta);
cif1 <- c(cif1); cif2 <- c(cif2)
cifs=cif1+cif2;
valn=2*(theta-1);
val1=(1+(theta-1)*(cifs))-( ((1+(theta-1)*cifs))^2-4*cif1*cif2*theta*(theta-1))^0.5;
vali=cif1*cif2;
valr <- vali;
valr[valn!=0] <- val1/valn;
valr <- matrix(valr,length(c(theta)),1)
rownames(valr)=colnames(coefs)
return(valr);
} ## }}}
##' @export
plack.cif2 <- function(cif1,cif2,theta)
{ ## {{{
theta <- exp(c(theta))
cif1 <- c(cif1); cif2 <- c(cif2)
cifs=cif1+cif2;
valn=2*(theta-1);
val1=(1+(theta-1)*(cifs))-( ((1+(theta-1)*cifs))^2-4*cif1*cif2*theta*(theta-1))^0.5;
vali=cif1*cif2;
valr <- vali;
valr[valn!=0] <- val1/valn;
return(valr);
} ## }}}
##' @export
predictPairPlack <- function(cif1,cif2,status1,status2,theta)
{ ## {{{
theta <- exp(c(theta))
cif1 <- c(cif1); cif2 <- c(cif2)
cifs=cif1+cif2;
valn=2*(theta-1);
val1=(1+(theta-1)*(cifs))-( ((1+(theta-1)*cifs))^2-4*cif1*cif2*theta*(theta-1))^0.5;
vali=cif1*cif2;
valr <- vali;
valr[valn!=0] <- val1/valn;
p11 <- valr;
p10 <- cif1-p11
p01 <- cif2-p11
p00 <- 1- p10-p01-p11
#
pred <- (status1==1)*(status2==1)*p11+ (status1==1)*(status2==0)*p10+
(status1==0)*(status2==1)*p01+ (status1==0)*(status2==0)*p00
return(pred);
} ## }}}
##' @export
summary.randomcif<-function (object, ...)
{ ## {{{
if (!inherits(object, "randomcif"))
stop("Must be a random.cif object")
cat("Random effect variance for variation due to clusters\n\n")
cat("Cause", attr(object, "cause1"), "and cause", attr(object,
"cause2"), fill = TRUE)
cat("\n")
if (sum(abs(object$score)) > 0.001) cat("WARNING: check score for convergence")
cat("\n")
coef.randomcif(object, ...)
} ## }}}
##' @export
coef.randomcif<- function (object, digits = 3, ...)
{ ## {{{
res <- cbind(object$theta, diag(object$var.theta)^0.5)
se <- diag(object$var.theta)^0.5
wald <- object$theta/se
waldp <- (1 - pnorm(abs(wald))) * 2
cor <- object$theta + 1
res <- as.matrix(cbind(res, wald, waldp, cor, se))
colnames(res) <- c("Coef.", "SE", "z", "P-val", "Cross odds ratio", "SE")
if (!is.null(object$thetanames)) rownames(res)<-object$thetanames
return(res)
} ## }}}
##' @export
print.randomcif<- function (x , digits = 3, ...)
{ ## {{{
} ## }}}
##' @export
summary.randomcifrv<-function (object, ...)
{ ## {{{
if (!inherits(object, "randomcifrv"))
stop("Must be a random.cifrv object")
cat("Random effect parameters for additive gamma random effects \n\n")
cat("Cause", attr(object, "cause1"), "and cause", attr(object,
"cause2"), fill = TRUE)
cat("\n")
if (sum(abs(object$score)) > 1e-06)
cat("WARNING: check score for convergence")
cat("\n")
res <- coef.randomcifrv(object, ...)
var.link <- attr(object,"var.link");
rv1 <- attr(object,"rv1");
theta.des <- attr(object,"pardes");
if (var.link==1) par <- theta.des %*% exp(object$theta) else par <- theta.des %*% object$theta
if (var.link==1) {
fp <- function(p){ res <- exp(p)/sum(rv1* (theta.des %*% exp(p))); return(res); }
e <- lava::estimate(coef=object$theta,vcov=object$var.theta,f=function(p) fp(p))
pare <- lava::estimate(coef=object$theta,vcov=object$var.theta,f=function(p) exp(p))
res <- list(estimate=res,h=e,exppar=pare)
} else {
fp <- function(p) { p/sum(rv1* (theta.des %*% p)) }
e <- lava::estimate(coef=object$theta,vcov=object$var.theta,f=function(p) fp(p))
res <- list(estimate=res,h=e)
}
res
} ## }}}
##' @export
coef.randomcifrv<- function (object, digits = 3, ...)
{ ## {{{
if (attr(object,"inverse")==1) elog <- 1 else elog <- 0;
if (elog==1) theta <- exp(object$theta) else theta <- object$theta
se <- diag(object$var.theta)^0.5
res <- cbind(object$theta, se)
wald <- object$theta/se
waldp <- (1 - pnorm(abs(wald))) * 2
res <- as.matrix(cbind(res, wald, waldp))
if (elog==0) colnames(res) <- c("Coef.", "SE", "z", "P-val")
if (elog==1) res <- cbind(res,exp(object$theta), exp(object$theta)^2*se)
if (elog==1) colnames(res) <- c("log-parameter","SE","z","P-val","exp(theta)","SE")
if (!is.null(object$thetanames)) rownames(res)<-object$thetanames
prmatrix(signif(res, digits))
} ## }}}
##' @export
print.randomcifrv<- function (x , digits = 3, ...)
{ ## {{{
summary(x, ...)
} ## }}}
## summary.cor<-function(object,digits=3,marg.cif=NULL,...)
## { ## {{{
## if (!inherits(object, "cor")) stop("Must be a cor.cif object")
## if (sum(abs(object$score))>0.001) warning("WARNING: check score for convergence\n")
## coefs <- coef.cor(object,...);
## outcase <- outconc <- NULL
## if (is.null(marg.cif)==FALSE) {
## marg.cif <- max(marg.cif)
## ## {{{
## if (attr(object,"Type")=="cor") {
## concordance <- exp(coefs[,1])*marg.cif^2/((1-marg.cif)+exp(coefs[,1])*marg.cif)
## conclower <- exp(coefs[,1]-1.96*coefs[,2])*marg.cif^2/((1-marg.cif)+exp(coefs[,1]-1.96*coefs[,2])*marg.cif)
## concup <- exp(coefs[,1]+1.96*coefs[,2])*marg.cif^2/((1-marg.cif)+exp(coefs[,1]+1.96*coefs[,2])*marg.cif)
## casewise <- concordance/marg.cif
## caselower <- conclower/marg.cif
## caseup <- concup/marg.cif
## } else if (attr(object,"Type")=="RR") {
## casewise<- exp(coefs[,1])*c(marg.cif)
## concordance <- exp(coefs[,1])*marg.cif^2
## caselower <- marg.cif*exp(coefs[,1]-1.96*coefs[,2])
## caseup <- marg.cif*exp(coefs[,1]+1.96*coefs[,2])
## conclower <- marg.cif^2* exp(coefs[,1]-1.96*coefs[,2])
## concup <- marg.cif^2*exp(coefs[,1]+1.96*coefs[,2])
## } else if (attr(object,"Type")=="OR-cif") {
## thetal <- coefs[,1]-1.96*coefs[,2]
## thetau <- coefs[,1]+1.96*coefs[,2]
## casewise<- plack.cif2(marg.cif,marg.cif,c(coefs[,1]))/marg.cif
## concordance <- plack.cif2(marg.cif,marg.cif,c(coefs[,1]))
## caselower <- plack.cif2(marg.cif,marg.cif,thetal)/marg.cif
## caseup <- plack.cif2(marg.cif,marg.cif,thetau)/marg.cif
## conclower <- plack.cif2(marg.cif,marg.cif,thetal)
## concup <- plack.cif2(marg.cif,marg.cif,thetau)
## }
## outcase <- cbind(casewise,caselower,caseup)
## outconc <- cbind(concordance,conclower,concup)
## rownames(outcase) <- rownames(outconc) <- rownames(coefs)
## colnames(outcase) <- c("casewise concordance","2.5 %","97.5%")
## colnames(outconc) <- c("concordance","2.5 %","97.5%")
## }
## ## }}}
## res <- list(casewise=outcase,concordance=outconc,estimates=coefs,marg=marg.cif,type=attr(object,"Type"),sym=attr(object,"sym"),cause1=attr(object,"cause1"),cause2=attr(object,"cause2"))
## class(res) <- "summary.cor"
## res
## } ## }}}
## coef.cor<-function(object,...)
## { ## {{{
## res <- cbind(object$theta, diag(object$var.theta)^0.5)
## se<-diag(object$var.theta)^0.5
## wald <- object$theta/se
## waldp <- (1 - pnorm(abs(wald))) * 2
## cor<-exp(object$theta)
## res <- as.matrix(cbind(res, wald, waldp,cor,se*cor))
## if (attr(object,"Type")=="cor")
## colnames(res) <- c("log-Coef.", "SE", "z", "P-val","Cross odds ratio","SE")
## else colnames(res) <- c("log-ratio Coef.", "SE", "z", "P-val","Ratio","SE")
## if (is.null((rownames(res)))==TRUE) rownames(res)<-rep(" ",nrow(res))
## return(res)
## } ## }}}
Try the mets package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.