R/ipcw-residualmean.r
In scheike/timereg: Flexible Regression Models for Survival Data
Defines functions
vcov.resmean
coef.resmean
res.mean
Documented in
coef.resmean
res.mean
#' Residual mean life (restricted)
#'
#' Fits a semiparametric model for the residual life (estimator=1): \deqn{ E(
#' \min(Y,\tau) -t | Y>=t) = h_1( g(t,x,z) ) } or cause specific years lost of
#' Andersen (2012) (estimator=3) \deqn{ E( \tau- \min(Y_j,\tau) | Y>=0) =
#' \int_0^t (1-F_j(s)) ds = h_2( g(t,x,z) ) } where \eqn{Y_j = \sum_j Y
#' I(\epsilon=j) + \infty * I(\epsilon=0)} or (estimator=2) \deqn{ E( \tau-
#' \min(Y_j,\tau) | Y<\tau, \epsilon=j) = h_3( g(t,x,z) ) = h_2(g(t,x,z))
#' F_j(\tau,x,z) }{} where \eqn{F_j(s,x,z) = P(Y<\tau, \epsilon=j | x,z )} for a
#' known link-function \eqn{h()} and known prediction-function \eqn{g(t,x,z)}
#'
#' Uses the IPCW for the score equations based on \deqn{ w(t)
#' \Delta(\tau)/P(\Delta(\tau)=1| T,\epsilon,X,Z) ( Y(t) - h_1(t,X,Z)) } and
#' where \eqn{\Delta(\tau)}{} is the at-risk indicator given data and requires a
#' IPCW model.
#'
#' Since timereg version 1.8.4. the response must be specified with the
#' \code{\link{Event}} function instead of the \code{\link{Surv}} function and
#' the arguments.
#'
#' @param formula a formula object, with the response on the left of a '~'
#' operator, and the terms on the right. The response must be a survival object
#' as returned by the `Event' function. The status indicator is not important
#' here. Time-invariant regressors are specified by the wrapper const(), and
#' cluster variables (for computing robust variances) by the wrapper cluster().
#' @param data a data.frame with the variables.
#' @param cause For competing risk models specificies which cause we consider.
#' @param restricted gives a possible restriction times for means.
#' @param times specifies the times at which the estimator is considered.
#' Defaults to all the times where an event of interest occurs, with the first
#' 10 percent or max 20 jump points removed for numerical stability in
#' simulations.
#' @param Nit number of iterations for Newton-Raphson algorithm.
#' @param clusters specifies cluster structure, for backwards compability.
#' @param gamma starting value for constant effects.
#' @param n.sim number of simulations in resampling.
#' @param weighted Not implemented. To compute a variance weighted version of
#' the test-processes used for testing time-varying effects.
#' @param model "additive", "prop"ortional.
#' @param detail if 0 no details are printed during iterations, if 1 details
#' are given.
#' @param interval specifies that we only consider timepoints where the
#' Kaplan-Meier of the censoring distribution is larger than this value.
#' @param resample.iid to return the iid decomposition, that can be used to
#' construct confidence bands for predictions
#' @param cens.model specified which model to use for the ICPW, KM is
#' Kaplan-Meier alternatively it may be "cox" or "aalen" model for further
#' flexibility.
#' @param cens.formula specifies the regression terms used for the regression
#' model for chosen regression model. When cens.model is specified, the
#' default is to use the same design as specified for the competing risks
#' model. "KM","cox","aalen","weights". "weights" are user specified weights
#' given is cens.weight argument.
#' @param time.pow specifies that the power at which the time-arguments is
#' transformed, for each of the arguments of the const() terms, default is 1
#' for the additive model and 0 for the proportional model.
#' @param time.pow.test specifies that the power the time-arguments is
#' transformed for each of the arguments of the non-const() terms. This is
#' relevant for testing if a coefficient function is consistent with the
#' specified form A_l(t)=beta_l t^time.pow.test(l). Default is 1 for the
#' additive model and 0 for the proportional model.
#' @param silent if 0 information on convergence problems due to non-invertible
#' derviates of scores are printed.
#' @param conv gives convergence criterie in terms of sum of absolute change of
#' parameters of model
#' @param estimator specifies what that is estimated.
#' @param cens.weights censoring weights for estimating equations.
#' @param conservative for slightly conservative standard errors.
#' @param weights weights for estimating equations.
#' @return returns an object of type 'comprisk'. With the following arguments:
#' \item{cum}{cumulative timevarying regression coefficient estimates are
#' computed within the estimation interval.} \item{var.cum}{pointwise variances
#' estimates. } \item{gamma}{estimate of proportional odds parameters of
#' model.} \item{var.gamma}{variance for gamma. } \item{score}{sum of absolute
#' value of scores.} \item{gamma2}{estimate of constant effects based on the
#' non-parametric estimate. Used for testing of constant effects.}
#' \item{obs.testBeq0}{observed absolute value of supremum of cumulative
#' components scaled with the variance.} \item{pval.testBeq0}{p-value for
#' covariate effects based on supremum test.} \item{obs.testBeqC}{observed
#' absolute value of supremum of difference between observed cumulative process
#' and estimate under null of constant effect.} \item{pval.testBeqC}{p-value
#' based on resampling.} \item{obs.testBeqC.is}{observed integrated squared
#' differences between observed cumulative and estimate under null of constant
#' effect.} \item{pval.testBeqC.is}{p-value based on resampling.}
#' \item{conf.band}{resampling based constant to construct 95\% uniform
#' confidence bands.} \item{B.iid}{list of iid decomposition of non-parametric
#' effects.} \item{gamma.iid}{matrix of iid decomposition of parametric
#' effects.} \item{test.procBeqC}{observed test process for testing of
#' time-varying effects} \item{sim.test.procBeqC}{50 resample processes for for
#' testing of time-varying effects} \item{conv}{information on convergence for
#' time points used for estimation.}
#' @author Thomas Scheike
#' @references
#' Andersen (2013), Decomposition of number of years lost according
#' to causes of death, Statistics in Medicine, 5278-5285.
#'
#' Scheike, and Cortese (2015), Regression Modelling of Cause Specific Years Lost,
#'
#' Scheike, Cortese and Holmboe (2015), Regression Modelling of Restricted
#' Residual Mean with Delayed Entry,
#' @keywords survival
#' @examples
#'
#' data(bmt);
#' tau <- 100
#'
#' ### residual restricted mean life
#' out<-res.mean(Event(time,cause>=1)~factor(tcell)+factor(platelet),data=bmt,cause=1,
#' times=0,restricted=tau,n.sim=0,model="additive",estimator=1);
#' summary(out)
#'
#' out<-res.mean(Event(time,cause>=1)~factor(tcell)+factor(platelet),data=bmt,cause=1,
#' times=seq(0,90,5),restricted=tau,n.sim=0,model="additive",estimator=1);
#' par(mfrow=c(1,3))
#' plot(out)
#'
#' ### restricted years lost given death
#' out21<-res.mean(Event(time,cause)~factor(tcell)+factor(platelet),data=bmt,cause=1,
#' times=0,restricted=tau,n.sim=0,model="additive",estimator=2);
#' summary(out21)
#' out22<-res.mean(Event(time,cause)~factor(tcell)+factor(platelet),data=bmt,cause=2,
#' times=0,restricted=tau,n.sim=0,model="additive",estimator=2);
#' summary(out22)
#'
#'
#' ### total restricted years lost
#' out31<-res.mean(Event(time,cause)~factor(tcell)+factor(platelet),data=bmt,cause=1,
#' times=0,restricted=tau,n.sim=0,model="additive",estimator=3);
#' summary(out31)
#' out32<-res.mean(Event(time,cause)~factor(tcell)+factor(platelet),data=bmt,cause=2,
#' times=0,restricted=tau,n.sim=0,model="additive",estimator=3);
#' summary(out32)
#'
#'
#' ### delayed entry
#' nn <- nrow(bmt)
#' entrytime <- rbinom(nn,1,0.5)*(bmt$time*runif(nn))
#' bmt$entrytime <- entrytime
#'
#' bmtw <- prep.comp.risk(bmt,times=tau,time="time",entrytime="entrytime",cause="cause")
#'
#' out<-res.mean(Event(time,cause>=1)~factor(tcell)+factor(platelet),data=bmtw,cause=1,
#' times=0,restricted=tau,n.sim=0,model="additive",estimator=1,
#' cens.model="weights",weights=bmtw$cw,cens.weights=1/bmtw$weights);
#' summary(out)
#'
#' @export
res.mean<-function(formula,data=parent.frame(),cause=1,restricted=NULL,times=NULL,Nit=50,
clusters=NULL,gamma=0,n.sim=0,weighted=0,model="additive",detail=0,interval=0.01,resample.iid=1,
cens.model="KM",cens.formula=NULL,time.pow=NULL,time.pow.test=NULL,silent=1,conv=1e-6,estimator=1,cens.weights=NULL,
conservative=1,weights=NULL){
## {{{
# restricted residual mean life models, using IPCW
# two models, additive and proportional
# trans=1 E(min(Y,tau) - t | Y>= t) = ( x' b(b)+ z' gam (tau-t))
# trans=2 E(min(Y,tau) - t | Y>= t) = exp( x' b(b)+ z' gam (tau-t))
# unrestricted residual mean life models, using IPCW
# trans=1 E(Y - t | Y>= t) = ( x' b(b)+ z' gam )
# trans=2 E(Y - t | Y>= t) = exp( x' b(b)+ z' gam )
if (model=="additive") trans<-1;
if (model=="prop") trans<-2;
if (model=="additive") Nit=2; ### one for estimation and one for residuals
# estimator =1 E(min(Y,tau)-t | Y>=t) or E(Y - t | Y>=t)
# estimator =2 Year lost due to cause 1 up to time tau given event
# estimator =2 E( tau - min(T,tau) | T <= tau, epsilon=j)
# estimator =3 PKA Years lost due to cause 1 up to time tau
# estimator =3 E( tau - min(T_j,tau) | T>t ) = \int_t^tau (1-F_j(s)) ds / S(t)
# estimator =3 E( tau - min(T,tau) | T <= tau, epsilon=j, T>t) * F_j(tau)
cens.tau <- 1
line <- 0
m<-match.call(expand.dots = FALSE);
m$gamma<-m$times<-m$cause<-m$Nit<-m$weighted<-m$n.sim<-
### m$cens.tau <-
m$model<- m$detail<- m$cens.model<-m$time.pow<-m$silent<-
m$interval<- m$clusters<-m$resample.iid<-m$restricted <- m$weights <-
m$time.pow.test<-m$conv<-m$estimator <- m$cens.weights <-
m$conservative <- m$cens.formula <- NULL
special <- c("const","cluster")
if (missing(data)) {
Terms <- terms(formula, special)
} else {
Terms <- terms(formula, special, data = data)
}
m$formula <- Terms
m[[1]] <- as.name("model.frame")
m <- eval(m, parent.frame())
mt <- attr(m, "terms")
intercept <- attr(mt, "intercept")
Y <- model.extract(m, "response")
## {{{ Event stuff
cens.code <- attr(Y,"cens.code")
if (ncol(Y)==3) stop("Left-truncation, through weights \n");
time2 <- eventtime <- Y[,1]
status <- delta <- Y[,2]
event <- (status==cause)
entrytime <- rep(0,length(time2))
if (sum(event)==0) stop("No events of interest in data\n");
## }}}
if (n.sim==0) sim<-0 else sim<-1; antsim<-n.sim;
des<-read.design(m,Terms)
X<-des$X; Z<-des$Z; npar<-des$npar; px<-des$px; pz<-des$pz;
covnamesX<-des$covnamesX; covnamesZ<-des$covnamesZ;
if(is.null(clusters)){ clusters <- des$clusters}
if(is.null(clusters)){
clusters <- 0:(nrow(X) - 1)
antclust <- nrow(X)
} else {
clusters <- as.integer(factor(clusters))-1
antclust <- length(unique(clusters))
}
pxz <-px+pz;
### always survival case
status; ### cause==1 and cause==0 censoring
if (is.null(restricted)) tau <- max(time2)+0.1 else tau <- restricted
if (is.null(times)) {
times<-sort(unique(time2[status==cause]));
if (tau>0) times <- times[times<tau];
###times<-times[-c(1:5)];
} else times <- sort(times);
n<-nrow(X); ntimes<-length(times);
if (npar==TRUE) {Z<-matrix(0,n,1); pg<-1; fixed<-0;} else {fixed<-1;pg<-pz;}
### non-cens at restricted time
delta<-(status!=cens.code)
if (!is.null(restricted)) {
if (cens.tau==1) {
time2tau <- pmin(time2,restricted)
deltatau <- rep(1,n)
deltatau[status==cens.code & time2<restricted] <- 0
} else { time2tau <- time2; deltatau <- delta*1;}
} else { time2tau <- time2; deltatau <- delta*1;}
if (is.null(weights)==TRUE) weights <- rep(1,n);
### cens model ## {{{
if (cens.model=="KM") { ## {{{
ud.cens<-survfit(Surv(time2,status==cens.code)~+1);
Gfit<-cbind(ud.cens$time,ud.cens$surv)
Gfit<-rbind(c(0,1),Gfit);
Gcx<-Cpred(Gfit,time2tau)[,2];
Gtimes<-Cpred(Gfit,times)[,2];
Gctimes<-Cpred(Gfit,times)[,2];
Gctimes[Gctimes<=0] <- 1 ## }}}
} else if (cens.model=="cox") { ## {{{
if (npar==TRUE) XZ<-X[,-1] else XZ<-cbind(X,Z)[,-1];
ud.cens<-cox.aalen(Surv(time2,status==cens.code)~prop(XZ),n.sim=0,robust=0);
Gcx<-Cpred(ud.cens$cum,time2tau)[,2];
RR<-exp(XZ %*% ud.cens$gamma)
Gcx<-exp(-Gcx*RR)
Gfit<-rbind(c(0,1),cbind(time2,Gcx));
Gctimes<-Cpred(Gfit,times,strict=TRUE)[,2];
Gctimes[Gctimes<=0] <- 1 ## }}}
} else if (cens.model=="aalen") { ## {{{
if (npar==TRUE) XZ <- X[,-1] else XZ <- cbind(X,Z)[,-1];
ud.cens <- aalen(Surv(time2,status==cens.code)~XZ,n.sim=0,robust=0);
Gcx <- Cpred(ud.cens$cum,time2tau)[,-1];
Gcx1 <- Gcx
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(time2,Gcx));
Gctimes<-Cpred(Gfit,times)[,2];
Gctimes[Gctimes<=0] <- 1 ## }}}
} else if (cens.model=="weights") {
if (length(weights)!=n) stop(paste("Weights length=",length(weights),"do not have length \n",length(time2),"problems with missing values?\n"));
Gcx <- cens.weights
ord2 <- order(time2)
Gctimes <- Cpred(cbind(time2[ord2],Gcx[ord2]),times)
Gctimes[Gctimes<=0] <- 1
} else { stop('Unknown censoring model') }
## }}}
if (resample.iid == 1) {
biid <- double(ntimes* antclust * px);
gamiid<- double(antclust *pg);
} else {
gamiid <- biid <- NULL;
}
ps<-px;
hess<-matrix(0,ps,ps); var<-score<-matrix(0,ntimes,ps+1);
if (sum(gamma)==0) gamma<-rep(0,pg); gamma2<-rep(0,ps);
test<-matrix(0,antsim,3*ps); testOBS<-rep(0,3*ps); unifCI<-c();
testval<-c(); rani<--round(runif(1)*10000);
Ut<-matrix(0,ntimes,ps+1); simUt<-matrix(0,ntimes,50*ps);
var.gamma<-matrix(0,pg,pg);
pred.covs.sem<-0
if (is.null(time.pow)==TRUE & !is.null(restricted)) time.pow<-rep(1,pg);
if (is.null(time.pow.test)==TRUE & !is.null(restricted)) time.pow<-rep(1,pg);
### if (is.null(time.pow)==TRUE & model=="additive") time.pow<-rep(1,pg);
if (is.null(time.pow)==TRUE & is.null(restricted)) time.pow<-rep(0,pg);
if (is.null(time.pow.test)==TRUE & is.null(restricted)) time.pow<-rep(0,pg);
silent <- c(silent,rep(0,ntimes-1));
if (!is.null(restricted)) time2 <- pmin(time2,restricted)
### print(restricted); print(time2); print(table(deltatau)); print(table(status));
### print(table(status)); print(times); print(summary(Gctimes)); print(summary(Gcx));
### important to start in 0 for linear model
est<-matrix(0,ntimes,ps+1)
if (model!="additive") est[,2] <- log(mean(time2))
betaS<-rep(0,ps);
if (model=="prop") betaS[1] <- log(mean(time2))
ordertime <- order(eventtime);
out<-.C("resmean",
as.double(times),as.integer(ntimes),as.double(time2),
as.integer(deltatau), as.integer(status),as.double(Gcx),
as.double(X),as.integer(n),as.integer(px),
as.integer(Nit), as.double(betaS), as.double(score),
as.double(hess), as.double(est), as.double(var),
as.integer(sim),as.integer(antsim),as.integer(rani),
as.double(test), as.double(testOBS), as.double(Ut),
as.double(simUt),as.integer(weighted),as.double(gamma),
as.double(var.gamma),as.integer(fixed),as.double(Z),
as.integer(pg),as.integer(trans),as.double(gamma2),
as.integer(cause),as.integer(line),as.integer(detail),
as.double(biid),as.double(gamiid),as.integer(resample.iid),
as.double(time.pow),as.integer(clusters),as.integer(antclust),
as.double(time.pow.test),as.integer(silent),
as.double(conv),as.double(tau),as.integer(estimator),
as.integer(cause),as.double(weights),as.double(Gctimes),
as.integer(ordertime-1),as.integer(conservative),as.integer(cens.code),
PACKAGE="timereg")
gamma<-matrix(out[[24]],pg,1);
var.gamma<-matrix(out[[25]],pg,pg);
gamma2<-matrix(out[[30]],ps,1);
rownames(gamma2)<-covnamesX;
conv <- list(convp=out[[41]],convd=out[[42]]);
if (fixed==0) gamma<-NULL;
if (resample.iid==1) {
biid<-matrix(out[[34]],ntimes,antclust*px);
if (fixed==1) gamiid<-matrix(out[[35]],antclust,pg) else gamiid<-NULL;
B.iid<-list();
for (i in (0:(antclust-1))*px) {
B.iid[[i/px+1]]<-matrix(biid[,i+(1:px)],ncol=px);
colnames(B.iid[[i/px+1]])<-covnamesX; }
if (fixed==1) colnames(gamiid)<-covnamesZ
} else B.iid<-gamiid<-NULL;
if (sim==1) {
simUt<-matrix(out[[22]],ntimes,50*ps); UIt<-list();
for (i in (0:49)*ps) UIt[[i/ps+1]]<-as.matrix(simUt[,i+(1:ps)]);
Ut<-matrix(out[[21]],ntimes,ps+1);
test<-matrix(out[[19]],antsim,3*ps); testOBS<-out[[20]];
supUtOBS<-apply(abs(as.matrix(Ut[,-1])),2,max);
p<-ps
for (i in 1:(3*p)) testval<-c(testval,pval(test[,i],testOBS[i]))
for (i in 1:p) unifCI<-as.vector(c(unifCI,percen(test[,i],0.95)));
pval.testBeq0<-as.vector(testval[1:p]);
pval.testBeqC<-as.vector(testval[(p+1):(2*p)]);
pval.testBeqC.is<-as.vector(testval[(2*p+1):(3*p)]);
obs.testBeq0<-as.vector(testOBS[1:p]);
obs.testBeqC<-as.vector(testOBS[(p+1):(2*p)]);
obs.testBeqC.is<-as.vector(testOBS[(2*p+1):(3*p)]);
sim.testBeq0<-as.matrix(test[,1:p]);
sim.testBeqC<-as.matrix(test[,(p+1):(2*p)]);
sim.testBeqC.is<-as.matrix(test[,(2*p+1):(3*p)]);
} else {test<-unifCI<-Ut<-UIt<-pval.testBeq0<-pval.testBeqC<-obs.testBeq0<-
obs.testBeqC<- sim.testBeq0<-sim.testBeqC<-
sim.testBeqC.is<- pval.testBeqC.is<-
obs.testBeqC.is<-NULL;
}
est<-matrix(out[[14]],ntimes,ps+1);
score<-matrix(out[[12]],ntimes,ps+1);
var<-matrix(out[[15]],ntimes,ps+1);
colnames(var)<-colnames(est)<-c("time",covnamesX);
if (sim>=1) {
colnames(Ut)<- c("time",covnamesX)
names(unifCI)<-names(pval.testBeq0)<- names(pval.testBeqC)<-
names(pval.testBeqC.is)<- names(obs.testBeq0)<- names(obs.testBeqC)<-
names(obs.testBeqC.is)<- colnames(sim.testBeq0)<- colnames(sim.testBeqC)<-
colnames(sim.testBeqC.is)<- covnamesX;
}
if (fixed==1) { rownames(gamma)<-c(covnamesZ);
colnames(var.gamma)<- rownames(var.gamma)<-c(covnamesZ); }
colnames(score)<-c("time",covnamesX);
if (is.na(sum(score))==TRUE) score<-NA else
if (sum(score[,-1])<0.00001) score<-sum(score[,-1]);
ud<-list(cum=est,var.cum=var,gamma=gamma,score=score,
gamma2=gamma2,var.gamma=var.gamma,robvar.gamma=var.gamma,
pval.testBeq0=pval.testBeq0,pval.testBeqC=pval.testBeqC,
obs.testBeq0=obs.testBeq0,
obs.testBeqC.is=obs.testBeqC.is,
obs.testBeqC=obs.testBeqC,pval.testBeqC.is=pval.testBeqC.is,
conf.band=unifCI,B.iid=B.iid,gamma.iid=gamiid,
test.procBeqC=Ut,sim.test.procBeqC=UIt,conv=conv,
cens.weights=Gcx,time=time2,delta.tau=deltatau,time2tau=time2tau)
ud$call<-match.call(); ud$model<-model; ud$n<-n;
ud$formula<-formula; class(ud)<-"resmean";
attr(ud, "Call") <- match.call()
attr(ud, "Formula") <- formula
attr(ud, "time.pow") <- time.pow
attr(ud, "cause") <- cause
attr(ud, "restricted") <- restricted
attr(ud, "times") <- times
attr(ud, "model") <- model
attr(ud, "estimator") <- estimator
return(ud);
} ## }}}
#' @export
print.resmean <- function (x,...) { ## {{{
object <- x; rm(x);
if (!inherits(object, 'resmean')) stop ("Must be an resmean object")
if (is.null(object$gamma)==TRUE) semi<-FALSE else semi<-TRUE
# We print information about object:
cat(paste("Residual mean model with",object$model,"\n"))
if (!is.null(attr(object,"restricted")))
cat(paste("Restricted at ",attr(object,"restricted","\n")))
cat("\n")
cat(" Nonparametric terms : ");
cat(colnames(object$cum)[-1]); cat(" \n");
if (semi) {
cat(" Parametric terms : ");
cat(rownames(object$gamma));
cat(" \n");
}
if (object$conv$convd>=1) {
cat("Warning problem with convergence for time points:\n")
cat(object$cum[object$conv$convp>0,1])
cat("\nReadjust analyses by removing points\n") }
cat(" \n");
} ## }}}
#' @export
coef.resmean <- function(object, digits=3,...) { ## {{{
coefBase(object,digits=digits)
} ## }}}
#' @export
summary.resmean <- function (object,digits = 3,ci=0, alpha=0.05,silent=0, ...) { ## {{{
if (!inherits(object, 'resmean')) stop ("Must be a resmean object")
if (is.null(object$gamma)==TRUE) semi<-FALSE else semi<-TRUE
if (silent==0) {
# We print information about object:
cat(paste("Residual mean model with",object$model,"\n"))
### cat("Residual mean model \n\n")
if (attr(object,"estimator")==1) cat("Residual mean life:")
if (attr(object,"estimator")==2) cat("Years Lost to cause given event:")
if (attr(object,"estimator")==3) cat("Years Lost to cause:")
if (!is.null(attr(object,"restricted")))
cat(paste("Restricted at ",attr(object,"restricted","\n")))
cat("\n"); cat("\n")
modelType<-object$model
#if (modelType=="additive" || modelType=="rcif")
if (sum(object$obs.testBeq0)==FALSE) cat("No test for non-parametric terms\n") else
timetest(object,digits=digits);
}
if (semi) { if (silent==0) cat("Parametric terms : \n");
out=coef.resmean(object);
out=signif(out,digits=digits)
if (silent==0) { print(out); cat(" \n"); }
} else { if (silent==0) cat("Non-Parametric terms : \n");
if (nrow(object$cum)==1) {
out=cbind(c(object$cum[,-1]),c(object$var.cum[,-1]^.5))
colnames(out) <- c("resmean","se")
out=signif(out,digits=digits)
} else {
se.cum <- object$var.cum[,-1]^.5
colnames(se.cum) <- paste("se",colnames(se.cum),paste="")
out=cbind(object$cum,se.cum)
out=signif(head(out),digits=digits)
}
if (silent==0) cat(" \n");
}
if (ci==1) { out <- round(cbind(out,out[,1]+qnorm(alpha/2)*out[,2],out[,1]+qnorm(1-alpha/2)*out[,2]),digits=digits);
nn <- ncol(out);
colnames(out)[(nn-1):nn] <- c("lower","upper")
}
if (silent==0) {
if (object$conv$convd>=1) {
cat("WARNING problem with convergence for time points:\n")
cat(object$cum[object$conv$convp>0,1])
cat("\nReadjust analyses by removing points\n\n")
}
cat(" Call: \n")
dput(attr(object, "Call"))
cat("\n")
}
out
} ## }}}
#' @export
vcov.resmean <- function(object, ...) {
rv <- object$robvar.gamma
if (!identical(rv, matrix(0, nrow = 1L, ncol = 1L))) rv # else return NULL
}
#' @export
plot.resmean <- function (x, pointwise.ci=1, hw.ci=0,
sim.ci=0, specific.comps=FALSE,level=0.05, start.time = 0,
stop.time = 0, add.to.plot=FALSE, mains=TRUE, xlab="Time",
ylab ="Coefficients",score=FALSE,...){ ## {{{
object <- x; rm(x);
if (!inherits(object,'resmean') ){ stop ("Must be output from res.mean function") }
if (score==FALSE) {
B<-object$cum;
V<-object$var.cum;
p<-dim(B)[[2]];
if (sum(specific.comps)==FALSE){
comp<-2:p
} else {
comp<-specific.comps+1
}
if (stop.time==0) {
stop.time<-max(B[,1]);
}
med<-B[,1]<=stop.time & B[,1]>=start.time
B<-B[med,];
V<-V[med,];
c.alpha<- qnorm(1-level/2)
for (v in comp) {
c.alpha<- qnorm(1-level/2)
est<-B[,v];
ul<-B[,v]+c.alpha*V[,v]^.5;
nl<-B[,v]-c.alpha*V[,v]^.5;
if (add.to.plot==FALSE) {
plot(B[,1],est,ylim=1.05*range(ul,nl),type="l",xlab=xlab,ylab=ylab)
if (mains==TRUE) title(main=colnames(B)[v]);
} else {
lines(B[,1],est,type="l");
}
if (pointwise.ci>=1) {
lines(B[,1],ul,lty=pointwise.ci,type="l");
lines(B[,1],nl,lty=pointwise.ci,type="l");
}
if (hw.ci>=1) {
if (level!=0.05){
cat("Hall-Wellner bands only 95 % \n");
}
tau<-length(B[,1])
nl<-B[,v]-1.27*V[tau,v]^.5*(1+V[,v]/V[tau,v])
ul<-B[,v]+1.27*V[tau,v]^.5*(1+V[,v]/V[tau,v])
lines(B[,1],ul,lty=hw.ci,type="l");
lines(B[,1],nl,lty=hw.ci,type="l");
}
if (sim.ci>=1) {
if (is.null(object$conf.band)==TRUE){
cat("Uniform simulation based bands only computed for n.sim> 0\n")
}
if (level!=0.05){
c.alpha<-percen(object$sim.testBeq0[,v-1],1-level)
} else {
c.alpha<-object$conf.band[v-1];
}
nl<-B[,v]-c.alpha*V[,v]^.5;
ul<-B[,v]+c.alpha*V[,v]^.5;
lines(B[,1],ul,lty=sim.ci,type="l");
lines(B[,1],nl,lty=sim.ci,type="l");
}
abline(h = 0)
}
} else {
# plot score proces
if (is.null(object$pval.testBeqC)==TRUE) {
cat("Simulations not done \n");
cat("To construct p-values and score processes under null n.sim>0 \n");
} else {
if (ylab=="Cumulative regression function"){
ylab<-"Test process";
}
dim1<-ncol(object$test.procBeqC)
if (sum(specific.comps)==FALSE){
comp<-2:dim1
} else {
comp<-specific.comps+1
}
for (i in comp){
ranyl<-range(object$test.procBeqC[,i]);
for (j in 1:50){
ranyl<-range(c(ranyl,(object$sim.test.procBeqC[[j]])[,i-1]));
}
mr<-max(abs(ranyl));
plot(object$test.procBeqC[,1],
object$test.procBeqC[,i],
ylim=c(-mr,mr),lwd=2,xlab=xlab,ylab=ylab,type="l")
if (mains==TRUE){
title(main=colnames(object$test.procBeqC)[i]);
}
for (j in 1:50){
lines(object$test.procBeqC[,1],
as.matrix(object$sim.test.procBeqC[[j]])[,i-1],col="grey",lwd=1,lty=1,type="l")
}
lines(object$test.procBeqC[,1],object$test.procBeqC[,i],lwd=2,type="l")
}
}
}
} ## }}}
scheike/timereg documentation built on Jan. 25, 2023, 3:43 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 vcov.resmean coef.resmean res.mean
Documented in coef.resmean res.mean
#' Residual mean life (restricted)
#'
#' Fits a semiparametric model for the residual life (estimator=1): \deqn{ E(
#' \min(Y,\tau) -t | Y>=t) = h_1( g(t,x,z) ) } or cause specific years lost of
#' Andersen (2012) (estimator=3) \deqn{ E( \tau- \min(Y_j,\tau) | Y>=0) =
#' \int_0^t (1-F_j(s)) ds = h_2( g(t,x,z) ) } where \eqn{Y_j = \sum_j Y
#' I(\epsilon=j) + \infty * I(\epsilon=0)} or (estimator=2) \deqn{ E( \tau-
#' \min(Y_j,\tau) | Y<\tau, \epsilon=j) = h_3( g(t,x,z) ) = h_2(g(t,x,z))
#' F_j(\tau,x,z) }{} where \eqn{F_j(s,x,z) = P(Y<\tau, \epsilon=j | x,z )} for a
#' known link-function \eqn{h()} and known prediction-function \eqn{g(t,x,z)}
#'
#' Uses the IPCW for the score equations based on \deqn{ w(t)
#' \Delta(\tau)/P(\Delta(\tau)=1| T,\epsilon,X,Z) ( Y(t) - h_1(t,X,Z)) } and
#' where \eqn{\Delta(\tau)}{} is the at-risk indicator given data and requires a
#' IPCW model.
#'
#' Since timereg version 1.8.4. the response must be specified with the
#' \code{\link{Event}} function instead of the \code{\link{Surv}} function and
#' the arguments.
#'
#' @param formula a formula object, with the response on the left of a '~'
#' operator, and the terms on the right. The response must be a survival object
#' as returned by the `Event' function. The status indicator is not important
#' here. Time-invariant regressors are specified by the wrapper const(), and
#' cluster variables (for computing robust variances) by the wrapper cluster().
#' @param data a data.frame with the variables.
#' @param cause For competing risk models specificies which cause we consider.
#' @param restricted gives a possible restriction times for means.
#' @param times specifies the times at which the estimator is considered.
#' Defaults to all the times where an event of interest occurs, with the first
#' 10 percent or max 20 jump points removed for numerical stability in
#' simulations.
#' @param Nit number of iterations for Newton-Raphson algorithm.
#' @param clusters specifies cluster structure, for backwards compability.
#' @param gamma starting value for constant effects.
#' @param n.sim number of simulations in resampling.
#' @param weighted Not implemented. To compute a variance weighted version of
#' the test-processes used for testing time-varying effects.
#' @param model "additive", "prop"ortional.
#' @param detail if 0 no details are printed during iterations, if 1 details
#' are given.
#' @param interval specifies that we only consider timepoints where the
#' Kaplan-Meier of the censoring distribution is larger than this value.
#' @param resample.iid to return the iid decomposition, that can be used to
#' construct confidence bands for predictions
#' @param cens.model specified which model to use for the ICPW, KM is
#' Kaplan-Meier alternatively it may be "cox" or "aalen" model for further
#' flexibility.
#' @param cens.formula specifies the regression terms used for the regression
#' model for chosen regression model. When cens.model is specified, the
#' default is to use the same design as specified for the competing risks
#' model. "KM","cox","aalen","weights". "weights" are user specified weights
#' given is cens.weight argument.
#' @param time.pow specifies that the power at which the time-arguments is
#' transformed, for each of the arguments of the const() terms, default is 1
#' for the additive model and 0 for the proportional model.
#' @param time.pow.test specifies that the power the time-arguments is
#' transformed for each of the arguments of the non-const() terms. This is
#' relevant for testing if a coefficient function is consistent with the
#' specified form A_l(t)=beta_l t^time.pow.test(l). Default is 1 for the
#' additive model and 0 for the proportional model.
#' @param silent if 0 information on convergence problems due to non-invertible
#' derviates of scores are printed.
#' @param conv gives convergence criterie in terms of sum of absolute change of
#' parameters of model
#' @param estimator specifies what that is estimated.
#' @param cens.weights censoring weights for estimating equations.
#' @param conservative for slightly conservative standard errors.
#' @param weights weights for estimating equations.
#' @return returns an object of type 'comprisk'. With the following arguments:
#' \item{cum}{cumulative timevarying regression coefficient estimates are
#' computed within the estimation interval.} \item{var.cum}{pointwise variances
#' estimates. } \item{gamma}{estimate of proportional odds parameters of
#' model.} \item{var.gamma}{variance for gamma. } \item{score}{sum of absolute
#' value of scores.} \item{gamma2}{estimate of constant effects based on the
#' non-parametric estimate. Used for testing of constant effects.}
#' \item{obs.testBeq0}{observed absolute value of supremum of cumulative
#' components scaled with the variance.} \item{pval.testBeq0}{p-value for
#' covariate effects based on supremum test.} \item{obs.testBeqC}{observed
#' absolute value of supremum of difference between observed cumulative process
#' and estimate under null of constant effect.} \item{pval.testBeqC}{p-value
#' based on resampling.} \item{obs.testBeqC.is}{observed integrated squared
#' differences between observed cumulative and estimate under null of constant
#' effect.} \item{pval.testBeqC.is}{p-value based on resampling.}
#' \item{conf.band}{resampling based constant to construct 95\% uniform
#' confidence bands.} \item{B.iid}{list of iid decomposition of non-parametric
#' effects.} \item{gamma.iid}{matrix of iid decomposition of parametric
#' effects.} \item{test.procBeqC}{observed test process for testing of
#' time-varying effects} \item{sim.test.procBeqC}{50 resample processes for for
#' testing of time-varying effects} \item{conv}{information on convergence for
#' time points used for estimation.}
#' @author Thomas Scheike
#' @references
#' Andersen (2013), Decomposition of number of years lost according
#' to causes of death, Statistics in Medicine, 5278-5285.
#'
#' Scheike, and Cortese (2015), Regression Modelling of Cause Specific Years Lost,
#'
#' Scheike, Cortese and Holmboe (2015), Regression Modelling of Restricted
#' Residual Mean with Delayed Entry,
#' @keywords survival
#' @examples
#'
#' data(bmt);
#' tau <- 100
#'
#' ### residual restricted mean life
#' out<-res.mean(Event(time,cause>=1)~factor(tcell)+factor(platelet),data=bmt,cause=1,
#' times=0,restricted=tau,n.sim=0,model="additive",estimator=1);
#' summary(out)
#'
#' out<-res.mean(Event(time,cause>=1)~factor(tcell)+factor(platelet),data=bmt,cause=1,
#' times=seq(0,90,5),restricted=tau,n.sim=0,model="additive",estimator=1);
#' par(mfrow=c(1,3))
#' plot(out)
#'
#' ### restricted years lost given death
#' out21<-res.mean(Event(time,cause)~factor(tcell)+factor(platelet),data=bmt,cause=1,
#' times=0,restricted=tau,n.sim=0,model="additive",estimator=2);
#' summary(out21)
#' out22<-res.mean(Event(time,cause)~factor(tcell)+factor(platelet),data=bmt,cause=2,
#' times=0,restricted=tau,n.sim=0,model="additive",estimator=2);
#' summary(out22)
#'
#'
#' ### total restricted years lost
#' out31<-res.mean(Event(time,cause)~factor(tcell)+factor(platelet),data=bmt,cause=1,
#' times=0,restricted=tau,n.sim=0,model="additive",estimator=3);
#' summary(out31)
#' out32<-res.mean(Event(time,cause)~factor(tcell)+factor(platelet),data=bmt,cause=2,
#' times=0,restricted=tau,n.sim=0,model="additive",estimator=3);
#' summary(out32)
#'
#'
#' ### delayed entry
#' nn <- nrow(bmt)
#' entrytime <- rbinom(nn,1,0.5)*(bmt$time*runif(nn))
#' bmt$entrytime <- entrytime
#'
#' bmtw <- prep.comp.risk(bmt,times=tau,time="time",entrytime="entrytime",cause="cause")
#'
#' out<-res.mean(Event(time,cause>=1)~factor(tcell)+factor(platelet),data=bmtw,cause=1,
#' times=0,restricted=tau,n.sim=0,model="additive",estimator=1,
#' cens.model="weights",weights=bmtw$cw,cens.weights=1/bmtw$weights);
#' summary(out)
#'
#' @export
res.mean<-function(formula,data=parent.frame(),cause=1,restricted=NULL,times=NULL,Nit=50,
clusters=NULL,gamma=0,n.sim=0,weighted=0,model="additive",detail=0,interval=0.01,resample.iid=1,
cens.model="KM",cens.formula=NULL,time.pow=NULL,time.pow.test=NULL,silent=1,conv=1e-6,estimator=1,cens.weights=NULL,
conservative=1,weights=NULL){
## {{{
# restricted residual mean life models, using IPCW
# two models, additive and proportional
# trans=1 E(min(Y,tau) - t | Y>= t) = ( x' b(b)+ z' gam (tau-t))
# trans=2 E(min(Y,tau) - t | Y>= t) = exp( x' b(b)+ z' gam (tau-t))
# unrestricted residual mean life models, using IPCW
# trans=1 E(Y - t | Y>= t) = ( x' b(b)+ z' gam )
# trans=2 E(Y - t | Y>= t) = exp( x' b(b)+ z' gam )
if (model=="additive") trans<-1;
if (model=="prop") trans<-2;
if (model=="additive") Nit=2; ### one for estimation and one for residuals
# estimator =1 E(min(Y,tau)-t | Y>=t) or E(Y - t | Y>=t)
# estimator =2 Year lost due to cause 1 up to time tau given event
# estimator =2 E( tau - min(T,tau) | T <= tau, epsilon=j)
# estimator =3 PKA Years lost due to cause 1 up to time tau
# estimator =3 E( tau - min(T_j,tau) | T>t ) = \int_t^tau (1-F_j(s)) ds / S(t)
# estimator =3 E( tau - min(T,tau) | T <= tau, epsilon=j, T>t) * F_j(tau)
cens.tau <- 1
line <- 0
m<-match.call(expand.dots = FALSE);
m$gamma<-m$times<-m$cause<-m$Nit<-m$weighted<-m$n.sim<-
### m$cens.tau <-
m$model<- m$detail<- m$cens.model<-m$time.pow<-m$silent<-
m$interval<- m$clusters<-m$resample.iid<-m$restricted <- m$weights <-
m$time.pow.test<-m$conv<-m$estimator <- m$cens.weights <-
m$conservative <- m$cens.formula <- NULL
special <- c("const","cluster")
if (missing(data)) {
Terms <- terms(formula, special)
} else {
Terms <- terms(formula, special, data = data)
}
m$formula <- Terms
m[[1]] <- as.name("model.frame")
m <- eval(m, parent.frame())
mt <- attr(m, "terms")
intercept <- attr(mt, "intercept")
Y <- model.extract(m, "response")
## {{{ Event stuff
cens.code <- attr(Y,"cens.code")
if (ncol(Y)==3) stop("Left-truncation, through weights \n");
time2 <- eventtime <- Y[,1]
status <- delta <- Y[,2]
event <- (status==cause)
entrytime <- rep(0,length(time2))
if (sum(event)==0) stop("No events of interest in data\n");
## }}}
if (n.sim==0) sim<-0 else sim<-1; antsim<-n.sim;
des<-read.design(m,Terms)
X<-des$X; Z<-des$Z; npar<-des$npar; px<-des$px; pz<-des$pz;
covnamesX<-des$covnamesX; covnamesZ<-des$covnamesZ;
if(is.null(clusters)){ clusters <- des$clusters}
if(is.null(clusters)){
clusters <- 0:(nrow(X) - 1)
antclust <- nrow(X)
} else {
clusters <- as.integer(factor(clusters))-1
antclust <- length(unique(clusters))
}
pxz <-px+pz;
### always survival case
status; ### cause==1 and cause==0 censoring
if (is.null(restricted)) tau <- max(time2)+0.1 else tau <- restricted
if (is.null(times)) {
times<-sort(unique(time2[status==cause]));
if (tau>0) times <- times[times<tau];
###times<-times[-c(1:5)];
} else times <- sort(times);
n<-nrow(X); ntimes<-length(times);
if (npar==TRUE) {Z<-matrix(0,n,1); pg<-1; fixed<-0;} else {fixed<-1;pg<-pz;}
### non-cens at restricted time
delta<-(status!=cens.code)
if (!is.null(restricted)) {
if (cens.tau==1) {
time2tau <- pmin(time2,restricted)
deltatau <- rep(1,n)
deltatau[status==cens.code & time2<restricted] <- 0
} else { time2tau <- time2; deltatau <- delta*1;}
} else { time2tau <- time2; deltatau <- delta*1;}
if (is.null(weights)==TRUE) weights <- rep(1,n);
### cens model ## {{{
if (cens.model=="KM") { ## {{{
ud.cens<-survfit(Surv(time2,status==cens.code)~+1);
Gfit<-cbind(ud.cens$time,ud.cens$surv)
Gfit<-rbind(c(0,1),Gfit);
Gcx<-Cpred(Gfit,time2tau)[,2];
Gtimes<-Cpred(Gfit,times)[,2];
Gctimes<-Cpred(Gfit,times)[,2];
Gctimes[Gctimes<=0] <- 1 ## }}}
} else if (cens.model=="cox") { ## {{{
if (npar==TRUE) XZ<-X[,-1] else XZ<-cbind(X,Z)[,-1];
ud.cens<-cox.aalen(Surv(time2,status==cens.code)~prop(XZ),n.sim=0,robust=0);
Gcx<-Cpred(ud.cens$cum,time2tau)[,2];
RR<-exp(XZ %*% ud.cens$gamma)
Gcx<-exp(-Gcx*RR)
Gfit<-rbind(c(0,1),cbind(time2,Gcx));
Gctimes<-Cpred(Gfit,times,strict=TRUE)[,2];
Gctimes[Gctimes<=0] <- 1 ## }}}
} else if (cens.model=="aalen") { ## {{{
if (npar==TRUE) XZ <- X[,-1] else XZ <- cbind(X,Z)[,-1];
ud.cens <- aalen(Surv(time2,status==cens.code)~XZ,n.sim=0,robust=0);
Gcx <- Cpred(ud.cens$cum,time2tau)[,-1];
Gcx1 <- Gcx
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(time2,Gcx));
Gctimes<-Cpred(Gfit,times)[,2];
Gctimes[Gctimes<=0] <- 1 ## }}}
} else if (cens.model=="weights") {
if (length(weights)!=n) stop(paste("Weights length=",length(weights),"do not have length \n",length(time2),"problems with missing values?\n"));
Gcx <- cens.weights
ord2 <- order(time2)
Gctimes <- Cpred(cbind(time2[ord2],Gcx[ord2]),times)
Gctimes[Gctimes<=0] <- 1
} else { stop('Unknown censoring model') }
## }}}
if (resample.iid == 1) {
biid <- double(ntimes* antclust * px);
gamiid<- double(antclust *pg);
} else {
gamiid <- biid <- NULL;
}
ps<-px;
hess<-matrix(0,ps,ps); var<-score<-matrix(0,ntimes,ps+1);
if (sum(gamma)==0) gamma<-rep(0,pg); gamma2<-rep(0,ps);
test<-matrix(0,antsim,3*ps); testOBS<-rep(0,3*ps); unifCI<-c();
testval<-c(); rani<--round(runif(1)*10000);
Ut<-matrix(0,ntimes,ps+1); simUt<-matrix(0,ntimes,50*ps);
var.gamma<-matrix(0,pg,pg);
pred.covs.sem<-0
if (is.null(time.pow)==TRUE & !is.null(restricted)) time.pow<-rep(1,pg);
if (is.null(time.pow.test)==TRUE & !is.null(restricted)) time.pow<-rep(1,pg);
### if (is.null(time.pow)==TRUE & model=="additive") time.pow<-rep(1,pg);
if (is.null(time.pow)==TRUE & is.null(restricted)) time.pow<-rep(0,pg);
if (is.null(time.pow.test)==TRUE & is.null(restricted)) time.pow<-rep(0,pg);
silent <- c(silent,rep(0,ntimes-1));
if (!is.null(restricted)) time2 <- pmin(time2,restricted)
### print(restricted); print(time2); print(table(deltatau)); print(table(status));
### print(table(status)); print(times); print(summary(Gctimes)); print(summary(Gcx));
### important to start in 0 for linear model
est<-matrix(0,ntimes,ps+1)
if (model!="additive") est[,2] <- log(mean(time2))
betaS<-rep(0,ps);
if (model=="prop") betaS[1] <- log(mean(time2))
ordertime <- order(eventtime);
out<-.C("resmean",
as.double(times),as.integer(ntimes),as.double(time2),
as.integer(deltatau), as.integer(status),as.double(Gcx),
as.double(X),as.integer(n),as.integer(px),
as.integer(Nit), as.double(betaS), as.double(score),
as.double(hess), as.double(est), as.double(var),
as.integer(sim),as.integer(antsim),as.integer(rani),
as.double(test), as.double(testOBS), as.double(Ut),
as.double(simUt),as.integer(weighted),as.double(gamma),
as.double(var.gamma),as.integer(fixed),as.double(Z),
as.integer(pg),as.integer(trans),as.double(gamma2),
as.integer(cause),as.integer(line),as.integer(detail),
as.double(biid),as.double(gamiid),as.integer(resample.iid),
as.double(time.pow),as.integer(clusters),as.integer(antclust),
as.double(time.pow.test),as.integer(silent),
as.double(conv),as.double(tau),as.integer(estimator),
as.integer(cause),as.double(weights),as.double(Gctimes),
as.integer(ordertime-1),as.integer(conservative),as.integer(cens.code),
PACKAGE="timereg")
gamma<-matrix(out[[24]],pg,1);
var.gamma<-matrix(out[[25]],pg,pg);
gamma2<-matrix(out[[30]],ps,1);
rownames(gamma2)<-covnamesX;
conv <- list(convp=out[[41]],convd=out[[42]]);
if (fixed==0) gamma<-NULL;
if (resample.iid==1) {
biid<-matrix(out[[34]],ntimes,antclust*px);
if (fixed==1) gamiid<-matrix(out[[35]],antclust,pg) else gamiid<-NULL;
B.iid<-list();
for (i in (0:(antclust-1))*px) {
B.iid[[i/px+1]]<-matrix(biid[,i+(1:px)],ncol=px);
colnames(B.iid[[i/px+1]])<-covnamesX; }
if (fixed==1) colnames(gamiid)<-covnamesZ
} else B.iid<-gamiid<-NULL;
if (sim==1) {
simUt<-matrix(out[[22]],ntimes,50*ps); UIt<-list();
for (i in (0:49)*ps) UIt[[i/ps+1]]<-as.matrix(simUt[,i+(1:ps)]);
Ut<-matrix(out[[21]],ntimes,ps+1);
test<-matrix(out[[19]],antsim,3*ps); testOBS<-out[[20]];
supUtOBS<-apply(abs(as.matrix(Ut[,-1])),2,max);
p<-ps
for (i in 1:(3*p)) testval<-c(testval,pval(test[,i],testOBS[i]))
for (i in 1:p) unifCI<-as.vector(c(unifCI,percen(test[,i],0.95)));
pval.testBeq0<-as.vector(testval[1:p]);
pval.testBeqC<-as.vector(testval[(p+1):(2*p)]);
pval.testBeqC.is<-as.vector(testval[(2*p+1):(3*p)]);
obs.testBeq0<-as.vector(testOBS[1:p]);
obs.testBeqC<-as.vector(testOBS[(p+1):(2*p)]);
obs.testBeqC.is<-as.vector(testOBS[(2*p+1):(3*p)]);
sim.testBeq0<-as.matrix(test[,1:p]);
sim.testBeqC<-as.matrix(test[,(p+1):(2*p)]);
sim.testBeqC.is<-as.matrix(test[,(2*p+1):(3*p)]);
} else {test<-unifCI<-Ut<-UIt<-pval.testBeq0<-pval.testBeqC<-obs.testBeq0<-
obs.testBeqC<- sim.testBeq0<-sim.testBeqC<-
sim.testBeqC.is<- pval.testBeqC.is<-
obs.testBeqC.is<-NULL;
}
est<-matrix(out[[14]],ntimes,ps+1);
score<-matrix(out[[12]],ntimes,ps+1);
var<-matrix(out[[15]],ntimes,ps+1);
colnames(var)<-colnames(est)<-c("time",covnamesX);
if (sim>=1) {
colnames(Ut)<- c("time",covnamesX)
names(unifCI)<-names(pval.testBeq0)<- names(pval.testBeqC)<-
names(pval.testBeqC.is)<- names(obs.testBeq0)<- names(obs.testBeqC)<-
names(obs.testBeqC.is)<- colnames(sim.testBeq0)<- colnames(sim.testBeqC)<-
colnames(sim.testBeqC.is)<- covnamesX;
}
if (fixed==1) { rownames(gamma)<-c(covnamesZ);
colnames(var.gamma)<- rownames(var.gamma)<-c(covnamesZ); }
colnames(score)<-c("time",covnamesX);
if (is.na(sum(score))==TRUE) score<-NA else
if (sum(score[,-1])<0.00001) score<-sum(score[,-1]);
ud<-list(cum=est,var.cum=var,gamma=gamma,score=score,
gamma2=gamma2,var.gamma=var.gamma,robvar.gamma=var.gamma,
pval.testBeq0=pval.testBeq0,pval.testBeqC=pval.testBeqC,
obs.testBeq0=obs.testBeq0,
obs.testBeqC.is=obs.testBeqC.is,
obs.testBeqC=obs.testBeqC,pval.testBeqC.is=pval.testBeqC.is,
conf.band=unifCI,B.iid=B.iid,gamma.iid=gamiid,
test.procBeqC=Ut,sim.test.procBeqC=UIt,conv=conv,
cens.weights=Gcx,time=time2,delta.tau=deltatau,time2tau=time2tau)
ud$call<-match.call(); ud$model<-model; ud$n<-n;
ud$formula<-formula; class(ud)<-"resmean";
attr(ud, "Call") <- match.call()
attr(ud, "Formula") <- formula
attr(ud, "time.pow") <- time.pow
attr(ud, "cause") <- cause
attr(ud, "restricted") <- restricted
attr(ud, "times") <- times
attr(ud, "model") <- model
attr(ud, "estimator") <- estimator
return(ud);
} ## }}}
#' @export
print.resmean <- function (x,...) { ## {{{
object <- x; rm(x);
if (!inherits(object, 'resmean')) stop ("Must be an resmean object")
if (is.null(object$gamma)==TRUE) semi<-FALSE else semi<-TRUE
# We print information about object:
cat(paste("Residual mean model with",object$model,"\n"))
if (!is.null(attr(object,"restricted")))
cat(paste("Restricted at ",attr(object,"restricted","\n")))
cat("\n")
cat(" Nonparametric terms : ");
cat(colnames(object$cum)[-1]); cat(" \n");
if (semi) {
cat(" Parametric terms : ");
cat(rownames(object$gamma));
cat(" \n");
}
if (object$conv$convd>=1) {
cat("Warning problem with convergence for time points:\n")
cat(object$cum[object$conv$convp>0,1])
cat("\nReadjust analyses by removing points\n") }
cat(" \n");
} ## }}}
#' @export
coef.resmean <- function(object, digits=3,...) { ## {{{
coefBase(object,digits=digits)
} ## }}}
#' @export
summary.resmean <- function (object,digits = 3,ci=0, alpha=0.05,silent=0, ...) { ## {{{
if (!inherits(object, 'resmean')) stop ("Must be a resmean object")
if (is.null(object$gamma)==TRUE) semi<-FALSE else semi<-TRUE
if (silent==0) {
# We print information about object:
cat(paste("Residual mean model with",object$model,"\n"))
### cat("Residual mean model \n\n")
if (attr(object,"estimator")==1) cat("Residual mean life:")
if (attr(object,"estimator")==2) cat("Years Lost to cause given event:")
if (attr(object,"estimator")==3) cat("Years Lost to cause:")
if (!is.null(attr(object,"restricted")))
cat(paste("Restricted at ",attr(object,"restricted","\n")))
cat("\n"); cat("\n")
modelType<-object$model
#if (modelType=="additive" || modelType=="rcif")
if (sum(object$obs.testBeq0)==FALSE) cat("No test for non-parametric terms\n") else
timetest(object,digits=digits);
}
if (semi) { if (silent==0) cat("Parametric terms : \n");
out=coef.resmean(object);
out=signif(out,digits=digits)
if (silent==0) { print(out); cat(" \n"); }
} else { if (silent==0) cat("Non-Parametric terms : \n");
if (nrow(object$cum)==1) {
out=cbind(c(object$cum[,-1]),c(object$var.cum[,-1]^.5))
colnames(out) <- c("resmean","se")
out=signif(out,digits=digits)
} else {
se.cum <- object$var.cum[,-1]^.5
colnames(se.cum) <- paste("se",colnames(se.cum),paste="")
out=cbind(object$cum,se.cum)
out=signif(head(out),digits=digits)
}
if (silent==0) cat(" \n");
}
if (ci==1) { out <- round(cbind(out,out[,1]+qnorm(alpha/2)*out[,2],out[,1]+qnorm(1-alpha/2)*out[,2]),digits=digits);
nn <- ncol(out);
colnames(out)[(nn-1):nn] <- c("lower","upper")
}
if (silent==0) {
if (object$conv$convd>=1) {
cat("WARNING problem with convergence for time points:\n")
cat(object$cum[object$conv$convp>0,1])
cat("\nReadjust analyses by removing points\n\n")
}
cat(" Call: \n")
dput(attr(object, "Call"))
cat("\n")
}
out
} ## }}}
#' @export
vcov.resmean <- function(object, ...) {
rv <- object$robvar.gamma
if (!identical(rv, matrix(0, nrow = 1L, ncol = 1L))) rv # else return NULL
}
#' @export
plot.resmean <- function (x, pointwise.ci=1, hw.ci=0,
sim.ci=0, specific.comps=FALSE,level=0.05, start.time = 0,
stop.time = 0, add.to.plot=FALSE, mains=TRUE, xlab="Time",
ylab ="Coefficients",score=FALSE,...){ ## {{{
object <- x; rm(x);
if (!inherits(object,'resmean') ){ stop ("Must be output from res.mean function") }
if (score==FALSE) {
B<-object$cum;
V<-object$var.cum;
p<-dim(B)[[2]];
if (sum(specific.comps)==FALSE){
comp<-2:p
} else {
comp<-specific.comps+1
}
if (stop.time==0) {
stop.time<-max(B[,1]);
}
med<-B[,1]<=stop.time & B[,1]>=start.time
B<-B[med,];
V<-V[med,];
c.alpha<- qnorm(1-level/2)
for (v in comp) {
c.alpha<- qnorm(1-level/2)
est<-B[,v];
ul<-B[,v]+c.alpha*V[,v]^.5;
nl<-B[,v]-c.alpha*V[,v]^.5;
if (add.to.plot==FALSE) {
plot(B[,1],est,ylim=1.05*range(ul,nl),type="l",xlab=xlab,ylab=ylab)
if (mains==TRUE) title(main=colnames(B)[v]);
} else {
lines(B[,1],est,type="l");
}
if (pointwise.ci>=1) {
lines(B[,1],ul,lty=pointwise.ci,type="l");
lines(B[,1],nl,lty=pointwise.ci,type="l");
}
if (hw.ci>=1) {
if (level!=0.05){
cat("Hall-Wellner bands only 95 % \n");
}
tau<-length(B[,1])
nl<-B[,v]-1.27*V[tau,v]^.5*(1+V[,v]/V[tau,v])
ul<-B[,v]+1.27*V[tau,v]^.5*(1+V[,v]/V[tau,v])
lines(B[,1],ul,lty=hw.ci,type="l");
lines(B[,1],nl,lty=hw.ci,type="l");
}
if (sim.ci>=1) {
if (is.null(object$conf.band)==TRUE){
cat("Uniform simulation based bands only computed for n.sim> 0\n")
}
if (level!=0.05){
c.alpha<-percen(object$sim.testBeq0[,v-1],1-level)
} else {
c.alpha<-object$conf.band[v-1];
}
nl<-B[,v]-c.alpha*V[,v]^.5;
ul<-B[,v]+c.alpha*V[,v]^.5;
lines(B[,1],ul,lty=sim.ci,type="l");
lines(B[,1],nl,lty=sim.ci,type="l");
}
abline(h = 0)
}
} else {
# plot score proces
if (is.null(object$pval.testBeqC)==TRUE) {
cat("Simulations not done \n");
cat("To construct p-values and score processes under null n.sim>0 \n");
} else {
if (ylab=="Cumulative regression function"){
ylab<-"Test process";
}
dim1<-ncol(object$test.procBeqC)
if (sum(specific.comps)==FALSE){
comp<-2:dim1
} else {
comp<-specific.comps+1
}
for (i in comp){
ranyl<-range(object$test.procBeqC[,i]);
for (j in 1:50){
ranyl<-range(c(ranyl,(object$sim.test.procBeqC[[j]])[,i-1]));
}
mr<-max(abs(ranyl));
plot(object$test.procBeqC[,1],
object$test.procBeqC[,i],
ylim=c(-mr,mr),lwd=2,xlab=xlab,ylab=ylab,type="l")
if (mains==TRUE){
title(main=colnames(object$test.procBeqC)[i]);
}
for (j in 1:50){
lines(object$test.procBeqC[,1],
as.matrix(object$sim.test.procBeqC[[j]])[,i-1],col="grey",lwd=1,lty=1,type="l")
}
lines(object$test.procBeqC[,1],object$test.procBeqC[,i],lwd=2,type="l")
}
}
}
} ## }}}
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.