Nothing
##' Internal function.
##' Calculates Inverse Probability of Censoring
##' Weights (IPCW) and adds them to a data.frame
##'
##' @title Inverse Probability of Censoring Weights
##' @param formula Formula specifying the censoring model
##' @param data data frame
##' @param cluster clustering variable
##' @param same.cens For clustered data, should same censoring be assumed (bivariate probability calculated as mininum of the marginal probabilities)
##' @param obs.only Return data with uncensored observations only
##' @param weight.name Name of weight variable in the new data.frame
##' @param trunc.prob If TRUE truncation probabilities are also calculated and stored in 'weight.name2' (based on Clayton-Oakes gamma frailty model)
##' @param weight.name2 Name of truncation probabilities
##' @param cens.model Censoring model (default Aalens additive model)
##' @param pairs For paired data (e.g. twins) only the complete pairs are returned (With pairs=TRUE)
##' @param theta.formula Model for the dependence parameter in the Clayton-Oakes model (truncation only)
##' @param ... Additional arguments to censoring model
##' @author Klaus K. Holst
##' @examples
##' \dontrun{
##' data(prt)
##' prtw <- ipw(Surv(time,status==0)~country, data=prt[sample(nrow(prt),5000),],
##' cluster="id",weight.name="w")
##' plot(0,type="n",xlim=range(prtw$time),ylim=c(0,1),xlab="Age",ylab="Probability")
##' count <- 0
##' for (l in unique(prtw$country)) {
##' count <- count+1
##' prtw <- prtw[order(prtw$time),]
##' with(subset(prtw,country==l),
##' lines(time,w,col=count,lwd=2))
##' }
##' legend("topright",legend=unique(prtw$country),col=1:4,pch=-1,lty=1)
##' }
ipw <- function(formula,data,cluster,
same.cens=FALSE,obs.only=TRUE,
weight.name="w",
trunc.prob=FALSE,weight.name2="wt",
cens.model="aalen", pairs=FALSE,
theta.formula=~1, ...) {
##iid=TRUE,
##cens.args <- c(list(formula,n.sim=0,robust=0,data=eval(data)),list(...))
if (tolower(cens.model)%in%c("weibull","phreg.par","phreg.weibull")) {
ud.cens <- phreg.par(formula,data,...)
pr <- predict(ud.cens)
noncens <- which(!ud.cens$status)
} else {
m <- match.call(expand.dots = FALSE)
m <- m[match(c("", "formula", "data", "subset", "na.action"),
names(m), nomatch = 0)]
special <- c("strata", "factor", "NN", "cluster", "dummy")
if (missing(data))
Terms <- terms(formula)
else Terms <- terms(formula, data = data)
m$formula <- Terms
m[[1]] <- as.name("model.frame")
M <- eval(m, parent.frame())
censtime <- model.extract(M, "response")
if (ncol(censtime)==3) {
status <- censtime[,3]
otimes <- censtime[,2]
ltimes <- censtime[,1]
} else {
status <- censtime[,2]
otimes <- censtime[,1]
}
noncens <- !status
if (is.null(attr(terms(formula,"prop"),"specials")$prop)) {
ud.cens <- aalen(formula,n.sim=0,robust=0,data=data,...)
XZ <- model.matrix(formula,data)
Gcx <- ud.cens$cum[prodlim::sindex(ud.cens$cum[,1],otimes),-1]
##Gcx<-Cpred(ud.cens$cum,otimes)[,-1];
Gcx<-exp(-apply(Gcx*XZ,1,sum))
} else {
ud.cens <- cox.aalen(formula,n.sim=0,robust=0,data=data,...)
XZ <- model.matrix(formula,data)
Gcx<-Cpred(ud.cens$cum,otimes)[,-1];
Gcx<-exp(-apply(Gcx*XZ,1,sum))
}
##ud.cens <- do.call(cens.model,cens.args)
Gcx[Gcx>1]<-1; Gcx[Gcx<0]<-0
pr <- Gcx
}
if (trunc.prob & ncol(censtime)==3) { ## truncation
data$truncsurv <- Surv(ltimes,otimes,noncens)
trunc.formula <- update(formula,truncsurv~.)
ud.trunc <- aalen(trunc.formula,data=data,robust=0,n.sim=0,residuals=0,silent=1,max.clust=NULL,
clusters=data[,cluster], ...)
dependX0 <- model.matrix(theta.formula,data)
twostage.fit <-two.stage(ud.trunc,
data=data,robust=0,detail=0,
theta.des=dependX0)#,Nit=20,step=1.0,notaylor=1)
X <- model.matrix(trunc.formula,data)
Xnam <- colnames(X)
ww <- fast.reshape(cbind(X,".num"=seq(nrow(X)),".lefttime"=ltimes),varying=c(".num",".lefttime"),id=data[,cluster])
dependX <- model.matrix(theta.formula,ww)
Prob <- predict.two.stage(twostage.fit,X=ww[,Xnam],
times=ww[,".lefttime1"],times2=ww[,".lefttime2"],
theta.des=dependX)
P0 <- numeric(nrow(X))
P0[ww[,".num1"]] <- Prob$St1t2
P0[ww[,".num2"]] <- Prob$St1t2
data[,weight.name2] <- P0
}
data[,weight.name] <- pr
if (same.cens & !missing(cluster)) {
message("Minimum weights...")
myord <- order(data[,cluster])
data <- data[myord,,drop=FALSE]
id <- table(data[,cluster])
if (pairs) {
gem <- data[,cluster]%in%(names(id)[id==2])
id <- id[id==2]
data <- data[gem,]
}
d0 <- subset(data,select=c(cluster,weight.name))
noncens <- with(data,!eval(terms(formula)[[2]][[3]]))
d0[,"observed."] <- noncens
timevar <- paste("_",cluster,weight.name,sep="")
d0[,timevar] <- unlist(lapply(id,seq))
Wide <- reshape(d0,direction="wide",timevar=timevar,idvar=cluster)
W <- apply(Wide[,paste(weight.name,1:2,sep=".")],1,
function(x) min(x,na.rm=TRUE))
Wmarg <- d0[,weight.name]
data[,weight.name] <- 1/Wmarg
Wmin <- rep(W,id)
obs1only <- rep(with(Wide, observed..1 & (is.na(observed..2) | !observed..2)),id)
obs2only <- rep(with(Wide, observed..2 & (is.na(observed..1) | !observed..1)),id)
obsOne <- which(na.omit(obs1only|obs2only))
obsBoth <- rep(with(Wide, !is.na(observed..1) & !is.na(observed..2) & observed..2 & observed..1),id)
data[obsBoth,weight.name] <-
ifelse(noncens[obsBoth],1/Wmin[obsBoth],0)
data[obsOne,weight.name] <-
ifelse(noncens[obsOne],1/Wmarg[obsOne],0)
}
if (obs.only)
data <- data[noncens,,drop=FALSE]
return(data)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.