library(survival) library(rootSolve) library(mvtnorm) library(np) library(MASS) ################################################################################################################### ### CNSFull() ### Description : Function to implement the conditional score method with one time-dependent immune biomarker X(u) on full cohort data. ### Reference : Tsiatis and Davidian, 2001 ### Note : 1. The sample size of full cohort data is N. ### 2. The maximum number of immune response measurements per subject is mk. ### 3. Details for the argument "model" which specifies the form of the Cox regression ### If model=="margin", exp{X(u)*beta} ### If model=="adj", exp{X(u)*beta+z^Teta} ### If model=="inter", exp{X(u)*beta+z^Teta+X(u)*z[,1]*gamma}, where z[,1] is a binary vector ### Arguments : ### t : N*mk matrix of time points when immune response values are measured. ### w : N*mk matrix of values of immune response at time points given by matrix t. ### z : If model=="margin", by default = NULL. If model=="adj", N*p matrix of time-independent variables. If method=="inter", N*p matrix of time-independent ### and the first column must be a vector of binary 0/1 indicators. ### time : N-dimensional vector of observed follow up time. ### event : N-dimensional vector of observed event status, 0=non-event, 1=event. ### model : character string specifies the form of Cox regression. See Note 3 above. ### para0 : vector of initial parameter values for (beta,eta^T,gamma)^T. Default is calculated from coxph() on observed immune responses. ### ftfun : function specifies the immune response trajectory, X(u)=alpha^Tftfun(u). Default is linear in time, i.e. ftfun=function(u) c(1,u), ### and X(u)=alpha1 + alpha2*u ### Value : the function returns a list including ### coef : estimated Cox model coefficients. ### covmtx : estimated covariance matrix for coef. ### merr : estimated variance of measurement error and its standard error. ### convergence : convergence parameters. Output from function multiroot() in Package rootSolve. CNSFull <- function(t,w,z=NULL,time,event,model,para0=NULL,ftfun=function(u) c(1,u)){ ## set up some measurements N <- nrow(t) # sample size mk <- ncol(t) # maximum number of measurement time points per subject q <- length(ftfun(1)) # number of parameters in the random effects model J <- apply(w,1,function(wi){sum(!is.na(wi))}) # number of measurements per subjects if(all(J<=q)) stop("Inadequate number of measurements to fit this random effects model.") if((model=="adj" | model=="inter") & is.null(z)) stop("Matrix z is missing.") ## reform z to be a matrix if(!is.null(z) & class(z)!="matrix") z <- as.matrix(z) p <- ifelse(model=="margin",0,ncol(z)) # number of time-independent covariates in the Cox regression model ## the first column of z is to be included in the interaction term if(model=="inter") trt <- z[,1] np <- ifelse(model=="margin",1,ifelse(model=="adj",1+p,2+p)) # number of parameters in the Cox model ## obtain initial parameter values from coxph() if not specified by the user if(is.null(para0)){ lt1 <- c(t(t)) lw <- c(t(w)) lid <- rep(1:N,each=mk) ltime <- rep(time,each=mk) lt2 <- c(lt1[-1],tail(ltime,1)) ldata <- data.frame(id=lid,time1=lt1,time2=lt2,w=lw,event=0,ltime) ldata <- ldata[complete.cases(ldata),] for(i in unique(ldata$id)){ idx <- which(ldata$id==i) ldata$time2[tail(idx,1)] <- time[i] } ldata <- ldata[ldata$time1 != ldata$time2,] for(i in which(event==1)){ idx <- which(ldata$id==i) if(length(idx)>0) ldata$event[tail(idx,1)] <- 1 } if(model=="margin"){ para0 <- coef(coxph(Surv(time1,time2,event)~w,data=ldata)) } if(model=="adj"){ tempz <- data.frame(id=1:N,z=z) ldata.z <- merge(ldata,tempz,by="id") names(ldata.z)[-(1:ncol(ldata))] <- paste("z",1:p,sep="") para0 <- coef(coxph(as.formula(paste("Surv(time1,time2,event)~w+", paste("z",1:p,collapse="+",sep=""))),data=ldata.z)) } if(model=="inter"){ tempz <- data.frame(id=1:N,z=z) ldata.z <- merge(ldata,tempz,by="id") names(ldata.z)[-(1:ncol(ldata))] <- paste("z",1:p,sep="") ldata.z$inter.w.z1 <- ldata.z$w*ldata.z$z1 para0 <- coef(coxph(as.formula(paste("Surv(time1,time2,event)~w+", paste("z",1:p,collapse="+",sep=""),"+inter.w.z1")),data=ldata.z)) } } ## estimate the variance of measurement error sigmae^2 resi <- rep(0,N) resi[which(J>=q)] <- sapply(which(J>=q),function(i){ ki <- !is.na(w[i,])&!is.na(t[i,]) di <- t(sapply(t[i,ki],ftfun)) xls <- drop(solve(crossprod(di))%*%crossprod(di,w[i,ki])) crossprod(w[i,ki]-di%*%xls) }) sigmae <- sqrt(sum(resi*(J>=q))/sum((J-q)*(J>=q))) ## some functions # J(u), the number of avlaialbe measurement time points up to time u Jiu <- function(i,u) sum(t[i,]<=u & !is.na(t[i,]) & !is.na(w[i,])) Ju <- function(u) sapply(1:N,Jiu,u) # Y(u), at risk process Yiu <- function(i,u) 1*(Jiu(i,u) >= q & time[i]>=u) Yu <- function(u) 1*(Ju(u)>=q & time>=u) # dN(u), event increment process diu <- function(i,u) 1*(Jiu(i,u) >= q & event[i]==1 & time[i]==u) du <- function(u) 1*(Ju(u)>=q & event==1 & time==u) ## functions to get least squares estimates \hat{\alpha} and SigmaR, which are later used to calculate the complete and ## sufficient statistics Q(u,theta,sigmae^2) in function Qiu() # for subject i at time u SigmaRxlsiu <- function(i,u){ kiu <- which(t[i,]<=u & !is.na(t[i,]) & !is.na(w[i,])) #(i,u) must satisfy length(Jiu)>=q ti <- t[i,]; wi <- w[i,]; f <- ftfun(u) diu <- t(sapply(ti[kiu],ftfun)) diuinv <- solve(crossprod(diu)) ss <- drop(sigmae^2*(t(f)%*%diuinv%*%f)) xx <- drop(diuinv %*% crossprod(diu,wi[kiu])) # least square estimates for random effects alpha return(list(SigmaR=ss,xls=xx)) } # for all N subjects at time u SigmaRxlsu <- function(u){ xx <- matrix(0,N,q) ss <-rep(0,N) for(i in which(Yu(u)==1)){ temp <- SigmaRxlsiu(i,u) ss[i] <- temp$SigmaR xx[i,] <- temp$xls } return(list(SigmaR=ss,xls=xx)) } # for all N subjects at all observed event time points orderu <- unique(sort(time[event==1])) xxss <- lapply(orderu,SigmaRxlsu) xlsall <- SigmaRall <- NULL for(i in 1:length(orderu)){ xlsall[[i]] <- xxss[[i]][[2]] SigmaRall[[i]] <- xxss[[i]][[1]] } ## functions for Q(u,theta,sigmae^2), estimating equations, and variance ## the functions are in similar structure for each model if(model=="margin"){ # Q_i(u,theta,sigmae^2) Qiu <- function(i,u,beta){ A <- SigmaRxlsiu(i,u) drop(crossprod(ftfun(u),A$xls) + A$SigmaR*beta*diu(i,u)) } # E_F^{1}(u,theta,sigmae^2)/E_F^{0}(u,theta,sigmae^2) Eu <- function(u,beta){ YY <- Yu(u) # at risk indicator RR <- YY if(sum(RR)==0) return(rep(NA,np)) ii <- which(orderu==u) ss <- SigmaRall[[ii]] xx <- xlsall[[ii]] dN <- du(u) QQ <- RR*drop(xx%*%ftfun(u) + ss*beta*dN) e0 <- RR*exp(QQ*beta-0.5*beta^2*ss) e1 <- RR*QQ*e0 sum(e1,na.rm=T) / sum(e0,na.rm=T) } # estimating equations eefun <- function(beta){ ps <- sapply(which(event==1&J>=q),function(i){ Qiu(i,time[i],beta)-Eu(time[i],beta) }) sum(ps,na.rm=T) } # variance estimates Varfun <- function(beta,sigmae){ Eu_VAR <- function(u,beta){ YY <- Yu(u) if(sum(YY)==0) return(list(E0=NA,E1=NA,dE0dpara=NA,dE1dpara=NA, dE0dsigmae2=NA,dE1dsigmae2=NA)) ii <- which(orderu==u) ss <- SigmaRall[[ii]] xx <- xlsall[[ii]] dN <- du(u) xhat <- drop(xx%*%ftfun(u)) # \hat X(u) QQ <- YY*(xhat + ss*beta*dN) # Q_i(u,theta,sigmae^2) e0 <- YY*exp(QQ*beta-0.5*beta^2*ss) # E_i^{0}(u,theta,sigmae^2) e1 <- YY*e0*QQ # E_i^{1}(u,theta,sigmae^2) E0 <- sum(e0,na.rm=T) # E_F^{0}(u,theta,sigmae^2) E1 <- sum(e1,na.rm=T) # E_F^{1}(u,theta,sigmae^2) # some derivatives de0dpara <- e0*(xhat+(2*dN-1)*ss*beta) de1dpara <- QQ*de0dpara + e0*dN*ss de0dsigmae2 <- e0*(dN-0.5)*beta^2*ss/sigmae^2 de1dsigmae2 <- QQ*de0dsigmae2 + e0*ss*beta*dN/sigmae^2 dE0dpara <- sum(de0dpara,na.rm=T) dE1dpara <- sum(de1dpara,na.rm=T) dE0dsigmae2 <- sum(de0dsigmae2,na.rm=T) dE1dsigmae2 <- sum(de1dsigmae2,na.rm=T) return(list(E0=E0,E1=E1,dE0dpara=dE0dpara,dE1dpara=dE1dpara, dE0dsigmae2=dE0dsigmae2,dE1dsigmae2=dE1dsigmae2)) } temps <- lapply(orderu,Eu_VAR,beta) # calculate the "meat" part of the sandwich variance estimate lambda0 <- rep(NA,length(orderu)) ee1 <- matrix(0,N,np) for(j in 1:length(orderu)){ u <- orderu[j] YY <- Yu(u) if(sum(YY)==0) next E0 <- temps[[j]]$E0 E1 <- temps[[j]]$E1 ss <- SigmaRall[[j]] xx <- xlsall[[j]] dN <- du(u) QQ <- YY*drop(xx%*%ftfun(u) + ss*beta*dN) e0 <- YY*exp(QQ*beta-0.5*beta^2*ss) lambda0[j] <- sum(dN)/E0 # baseline hazard tempee1 <- (QQ - E1/E0)*(dN - YY*e0*lambda0[j]) tempee1[is.na(tempee1)] <- 0 ee1 <- ee1 + tempee1 } ee2 <- (resi-(J-q)*sigmae^2)*(J>=q) ee <- cbind(ee1,ee2) # calculate the "bread" part of the sandwich variance estimate dee1 <- matrix(0,np,np+1) for(i in which(event==1&J>=q)){ j <- which(orderu==time[i]) YY <- Yu(time[i]) if(sum(YY)==0) next E0 <- temps[[j]]$E0 E1 <- temps[[j]]$E1 dE0dpara <- temps[[j]]$dE0dpara dE1dpara <- temps[[j]]$dE1dpara dE0dsigmae2 <- temps[[j]]$dE0dsigmae2 dE1dsigmae2 <- temps[[j]]$dE1dsigmae2 si <- SigmaRall[[j]][i] dee1dpara <- -(E0*dE1dpara - tcrossprod(E1,dE0dpara))/E0^2 + si dee1dsigmae2 <- -(E0*dE1dsigmae2 - E1*dE0dsigmae2)/E0^2 + si*beta/sigmae^2 tempdee1 <- cbind(dee1dpara,dee1dsigmae2) tempdee1[is.na(tempdee1)] <- 0 dee1 <- dee1 + tempdee1 } dee2dsigmae2 <- -sum((J-q)*(J>=q),na.rm=T) dee2 <- c(rep(0,np),dee2dsigmae2) dee <- rbind(dee1,dee2) # get the sandwich variance estimate bread <- dee/N meat <- crossprod(ee)/N bbinverse <- try(solve(bread),T) if(class(bbinverse)=="try-error") varest <- matrix(NA,np+1,np+1) else varest <- bbinverse%*%meat%*%t(bbinverse)/N colnames(varest) <- rownames(varest) <- c("X(u)","sigmae2") # get the baseline hazard estimate lambda0 <- cbind(time=orderu,lambda0=lambda0) return(list(varest=varest,lambda0=lambda0)) } } if(model=="adj"){ Qiu <- function(i,u,beta){ A <- SigmaRxlsiu(i,u) drop(crossprod(ftfun(u),A$xls) + A$SigmaR*beta*diu(i,u)) } Eu <- function(u,beta,eta){ YY <- Yu(u) RR <- YY if(sum(RR)==0) return(rep(NA,np)) ii <- which(orderu==u) ss <- SigmaRall[[ii]] xx <- xlsall[[ii]] dN <- du(u) ez <- drop(exp(z%*%eta)) QQ <- RR*drop(xx%*%ftfun(u) + ss*beta*dN) e0 <- RR*exp(QQ*beta-0.5*beta^2*ss)*ez e1 <- RR*e0*cbind(QQ,z) colSums(e1,na.rm=T) / sum(e0,na.rm=T) } eefun <- function(para){ beta <- para[1] eta <- para[-1] ps <- sapply(which(event==1&J>=q),function(i){ c(Qiu(i,time[i],beta),z[i,])-Eu(time[i],beta,eta) }) rowSums(ps,na.rm=T) } Varfun <- function(para,sigmae){ beta <- para[1] eta <- para[-1] Eu_VAR <- function(u,beta,eta){ YY <- Yu(u) if(sum(YY)==0) return(list(E0=NA,E1=rep(NA,np),dE0dpara=rep(NA,np), dE1dpara=rep(NA,np^2), dE0dsigmae2=NA,dE1dsigmae2=rep(NA,np))) ii <- which(orderu==u) ss <- SigmaRall[[ii]] xx <- xlsall[[ii]] dN <- du(u) xhat <- drop(xx%*%ftfun(u)) QQ <- YY*(xhat + ss*beta*dN) Qz <- YY*cbind(QQ,z) ez <- drop(exp(z%*%eta)) e0 <- YY*exp(QQ*beta-0.5*beta^2*ss)*ez e1 <- YY*e0*Qz E0 <- sum(e0,na.rm=T) E1 <- colSums(e1,na.rm=T) de0dpara <- e0*cbind(xhat+(2*dN-1)*ss*beta,z) de1dpara <- t(sapply(1:N,function(i){ a <- tcrossprod(Qz[i,],de0dpara[i,]) a[1,1] <- a[1,1] + e0[i]*dN[i]*ss[i] a })) de0dsigmae2 <- e0*(dN-0.5)*beta^2*ss/sigmae^2 de1dsigmae2 <- cbind(Qz[,1]*de0dsigmae2 + e0*ss*beta*dN/sigmae^2, Qz[,-1]*de0dsigmae2) dE0dpara <- colSums(de0dpara,na.rm=T) dE1dpara <- colSums(de1dpara,na.rm=T) dE0dsigmae2 <- sum(de0dsigmae2,na.rm=T) dE1dsigmae2 <- colSums(de1dsigmae2,na.rm=T) return(list(E0=E0,E1=E1,dE0dpara=dE0dpara,dE1dpara=dE1dpara, dE0dsigmae2=dE0dsigmae2,dE1dsigmae2=dE1dsigmae2)) } temps <- lapply(orderu,Eu_VAR,beta,eta) lambda0 <- rep(NA,length(orderu)) ee1 <- matrix(0,N,np) for(j in 1:length(orderu)){ u <- orderu[j] YY <- Yu(u) if(sum(YY)==0) next E0 <- temps[[j]]$E0 E1 <- temps[[j]]$E1 ss <- SigmaRall[[j]] xx <- xlsall[[j]] dN <- du(u) QQ <- YY*drop(xx%*%ftfun(u) + ss*beta*dN) Qz <- YY*cbind(QQ,z) ez <- drop(exp(z%*%eta)) e0 <- YY*exp(QQ*beta-0.5*beta^2*ss)*ez lambda0[j] <- sum(dN)/E0 tempee1 <- sweep(Qz,MARGIN=2,E1/E0,"-")*(dN - YY*e0*lambda0[j]) tempee1[is.na(tempee1)] <- 0 ee1 <- ee1 + tempee1 } ee2 <- (resi-(J-q)*sigmae^2)*(J>=q) ee <- cbind(ee1,ee2) dee1 <- matrix(0,np,np+1) for(i in which(event==1)){ j <- which(orderu==time[i]) YY <- Yu(time[i]) if(sum(YY)==0) next E0 <- temps[[j]]$E0 E1 <- temps[[j]]$E1 dE0dpara <- temps[[j]]$dE0dpara dE1dpara <- matrix(temps[[j]]$dE1dpara,np,np) dE0dsigmae2 <- temps[[j]]$dE0dsigmae2 dE1dsigmae2 <- temps[[j]]$dE1dsigmae2 si <- SigmaRall[[j]][i] dee1dpara <- -(E0*dE1dpara - tcrossprod(E1,dE0dpara))/E0^2 dee1dpara[1,1] <- dee1dpara[1,1] + si dee1dsigmae2 <- -(E0*dE1dsigmae2 - E1*dE0dsigmae2)/E0^2 dee1dsigmae2[1] <- dee1dsigmae2[1] + si*beta/sigmae^2 tempdee1 <- cbind(dee1dpara,dee1dsigmae2) tempdee1[is.na(tempdee1)] <- 0 dee1 <- dee1 + tempdee1 } dee2dsigmae2 <- -sum((J-q)*(J>=q),na.rm=T) dee2 <- c(rep(0,np),dee2dsigmae2) dee <- rbind(dee1,dee2) bread <- dee/N meat <- crossprod(ee)/N bbinverse <- try(solve(bread),T) if(class(bbinverse)=="try-error") varest <- matrix(NA,np+1,np+1) else varest <- bbinverse%*%meat%*%t(bbinverse)/N colnames(varest) <- rownames(varest) <- c("X(u)",paste("z",1:p,sep=""),"sigmae2") lambda0 <- cbind(time=orderu,lambda0=lambda0) return(list(varest=varest,lambda0=lambda0)) } } if(model=="inter"){ Qiu <- function(i,u,beta,gamma){ A <- SigmaRxlsiu(i,u) drop(crossprod(ftfun(u),A$xls) + A$SigmaR*(beta+gamma*trt[i])*diu(i,u)) } Eu <- function(u,beta,eta,gamma){ YY <- Yu(u) RR <- YY if(sum(RR)==0) return(rep(NA,np)) ii <- which(orderu==u) ss <- SigmaRall[[ii]] xx <- xlsall[[ii]] bg <- (beta+gamma*trt) dN <- du(u) ez <- drop(exp(z%*%eta)) QQ <- RR*drop(xx%*%ftfun(u) + ss*bg*dN) e0 <- RR*exp(QQ*bg-0.5*bg^2*ss)*ez e1 <- RR*e0*cbind(QQ,z,QQ*trt) colSums(e1,na.rm=T) / sum(e0,na.rm=T) } eefun <- function(para){ beta <- para[1] eta <- para[2:(np-1)] gamma <- para[np] ps <- sapply(which(event==1&J>=q),function(i){ qi <- Qiu(i,time[i],beta,gamma) c(qi,z[i,],qi*trt[i])-Eu(time[i],beta,eta,gamma) }) rowSums(ps,na.rm=T) } Varfun <- function(para,sigmae){ beta <- para[1] eta <- para[2:(np-1)] gamma <- para[np] bg <- (beta+gamma*trt) Eu_VAR <- function(u,beta,eta,gamma){ YY <- Yu(u) if(sum(YY)==0) return(list(E0=NA,E1=rep(NA,np),dE0dpara=rep(NA,np), dE1dpara=rep(NA,np^2), dE0dsigmae2=NA,dE1dsigmae2=rep(NA,np))) ii <- which(orderu==u) ss <- SigmaRall[[ii]] xx <- xlsall[[ii]] bg <- (beta+gamma*trt) dN <- du(u) xhat <- drop(xx%*%ftfun(u)) QQ <- YY*(xhat + ss*bg*dN) Qz <- YY*cbind(QQ,z,QQ*trt) ez <- drop(exp(z%*%eta)) e0 <- YY*exp(QQ*bg-0.5*bg^2*ss)*ez e1 <- YY*e0*Qz E0 <- sum(e0,na.rm=T) E1 <- colSums(e1,na.rm=T) de0dpara <- e0*cbind(xhat+(2*dN-1)*ss*bg, z, xhat*trt+(2*dN-1)*ss*bg*trt) de1dpara <- t(sapply(1:N,function(i){ tcrossprod(Qz[i,],de0dpara[i,]) })) de1dpara <- de1dpara + e0*ss*dN*cbind(1,matrix(0,N,p), trt,matrix(0,N,p*np), trt,matrix(0,N,p), trt^2) de0dsigmae2 <- e0*(dN-0.5)*bg^2*ss/sigmae^2 de1dsigmae2 <- Qz*de0dsigmae2 + e0*ss*bg*dN/sigmae^2*cbind(1,matrix(0,N,p),trt) dE0dpara <- colSums(de0dpara,na.rm=T) dE1dpara <- colSums(de1dpara,na.rm=T) dE0dsigmae2 <- sum(de0dsigmae2,na.rm=T) dE1dsigmae2 <- colSums(de1dsigmae2,na.rm=T) return(list(E0=E0,E1=E1,dE0dpara=dE0dpara,dE1dpara=dE1dpara, dE0dsigmae2=dE0dsigmae2,dE1dsigmae2=dE1dsigmae2)) } temps <- lapply(orderu,Eu_VAR,beta,eta,gamma) # calculate estimating equation values for each subject lambda0 <- rep(NA,length(orderu)) ee1 <- matrix(0,N,np) for(j in 1:length(orderu)){ u <- orderu[j] YY <- Yu(u) if(sum(YY)==0) next E0 <- temps[[j]]$E0 E1 <- temps[[j]]$E1 ss <- SigmaRall[[j]] xx <- xlsall[[j]] dN <- du(u) QQ <- YY*drop(xx%*%ftfun(u) + ss*bg*dN) Qz <- YY*cbind(QQ,z,QQ*trt) ez <- drop(exp(z%*%eta)) e0 <- YY*exp(QQ*bg-0.5*bg^2*ss)*ez lambda0[j] <- sum(dN)/E0 tempee1 <- sweep(Qz,MARGIN=2,E1/E0,"-")*(dN - YY*e0*lambda0[j]) tempee1[is.na(tempee1)] <- 0 ee1 <- ee1 + tempee1 } ee2 <- (resi-(J-q)*sigmae^2)*(J>=q) ee <- cbind(ee1,ee2) dee1 <- matrix(0,np,np+1) for(i in which(event==1)){ j <- which(orderu==time[i]) YY <- Yu(time[i]) if(sum(YY)==0) next E0 <- temps[[j]]$E0 E1 <- temps[[j]]$E1 dE0dpara <- temps[[j]]$dE0dpara dE1dpara <- matrix(temps[[j]]$dE1dpara,np,np) dE0dsigmae2 <- temps[[j]]$dE0dsigmae2 dE1dsigmae2 <- temps[[j]]$dE1dsigmae2 si <- SigmaRall[[j]][i] dee1dpara <- -(E0*dE1dpara - tcrossprod(E1,dE0dpara))/E0^2 dqzidpara <- c(si,rep(0,p),trt[i]*si) dee1dpara[1,] <- dee1dpara[1,] + dqzidpara dee1dpara[p+2,] <- dee1dpara[p+2,] + trt[i]*dqzidpara dee1dsigmae2 <- -(E0*dE1dsigmae2 - E1*dE0dsigmae2)/E0^2 + c(1,rep(0,p),trt[i])*si*(bg[i])/sigmae^2 tempdee1 <- cbind(dee1dpara,dee1dsigmae2) tempdee1[is.na(tempdee1)] <- 0 dee1 <- dee1 + tempdee1 } dee2dsigmae2 <- -sum((J-q)*(J>=q),na.rm=T) dee2 <- c(rep(0,2+p),dee2dsigmae2) dee <- rbind(dee1,dee2) bread <- dee/N meat <- crossprod(ee)/N bbinverse <- try(solve(bread),T) if(class(bbinverse)=="try-error") varest <- matrix(NA,np+1,np+1) else varest <- bbinverse%*%meat%*%t(bbinverse)/N colnames(varest) <- rownames(varest) <- c("X(u)",paste("z",1:p,sep=""),"X(u)*z1","sigmae2") lambda0 <- cbind(time=orderu,lambda0=lambda0) return(list(varest=varest,lambda0=lambda0)) } } ## numeric solution from estimating equations for Cox regression parameters (apply all default settings in multiroot() function) paraest <- multiroot(eefun,start=para0) ## sandwich variance estimate for estimated Cox regression parameters vest <- Varfun(paraest$root,sigmae) ## summary function output covmtx <- vest$varest[1:np,1:np] coef <- c(paraest$root) names(coef) <- colnames(covmtx) #bhazard <- vest$lambda0 error <- c(sigmae^2,sqrt(vest$varest[np+1,np+1])) names(error) <- c("Est.","Std.Err.") convergence <- list(f.root=paraest$f.root,iter=paraest$iter,estim.precis=paraest$estim.precis) names(convergence$f.root) <- names(coef) return(list(coef=coef,covmtx=covmtx,merr=error, convergence=convergence)) } ########################################################################################## ### CNSBernoulliIPW() ### Description : Function to implement IPW conditional score estimator on two-phase data from Bernoulli sampling. ### Note : The sample size of full cohort data is N and the total number of subjects sampled ### into the phaseII sample is n. The number of rows of the matrices ### and the length of vectors used in following arguments must all be N. Values for subjects not sampled for ### immune response measurements are NA. The maximum number of immune response measurements per subject is mk. ### Arguments : ### t : N*mk matrix of time points to make immune response measurements. ### w : N*mk matrix of immune response values. ### z : If model=="margin", by default = NULL. If model=="adj", N*p matrix of time-independent variables. If method=="inter", ### N*p matrix of time-independent and the first column must be a vector of binary 0/1 indicators. ### time : N-dimensional vector of follow up time. ### event : N-dimensional vector of event status, 0=non-event, 1=event. ### model : which form of Cox regression is fitted. If "margin", exp(X(u)*beta). If "adj", exp(X(u)*beta+z^Teta). If "inter", ### exp(X(u)*beta+z^Teta+X(u)*z[,1]*gamma), and the first column of z must be 0/1 treatment indicators. ### Pi : N-dimensional vector of pre-specified sampling probabilities, if ==NULL, estimated from X (see below). ### phaseII : N-dimensional vector of phaseII sample indicators, 0=not a phaseII sample, 1= a phaseII sample. ### X : matrix of variables to estimate the sampling probabilities if Pi is Null. If simple random sampling, then ### leave X=NULL or X=matrix(1,N,1), and sampling probabilites are estimated as n/N. ### para0 : vector of initial parameters for (beta,eta^T,gamma)^T. Default is calculated from ### coxph() on observed immune responses. ### ftfun : function specifies the immune response trajectory, X(u) = alpha^Tftfun(u). ### Default is linear in time, i.e. ftfun = function(u) c(1,u). ### Value : the function returns a list including ### coef : estimated Cox model coefficients. ### covmtx : estimated covariance matrix for coef. ### merr : estimated variance of measurement error and its standard error. ### convergence : convergence parameters. Output from multiroot(). CNSBernoulliIPW <- function(t,w,z=NULL,time,event,model,Pi=NULL,phaseII,X=NULL, para0=NULL,ftfun=function(u) c(1,u)){ N <- nrow(t) n <- sum(phaseII) mk <- ncol(t) ## number of parameters in the random effects model q <- length(ftfun(1)) ## number of measurements per subjects J <- apply(w,1,function(wi){sum(!is.na(wi))}) if(all(J<=q)) stop("Inadequate number of measurements to fit this random effects model.") if((model=="adj" | model=="inter") & is.null(z)) stop("Matrix z is missing.") p <- ifelse(model=="margin",0,ncol(z)) ## number of parameters in the Cox model np <- ifelse(model=="margin",1,ifelse(model=="adj",1+p,2+p)) ## use only phaseII subjects t <- t[phaseII==1,] w <- w[phaseII==1,] event <- event[phaseII==1] time <- time[phaseII==1] J <- J[phaseII==1] if(!is.null(z) & class(z)!="matrix") z <- as.matrix(z) z <- as.matrix(z[phaseII==1,]) if(model=="inter") trt <- z[,1] ## estimate the sampling probs rho <- NULL if(is.null(Pi)){ ## random sampling if(is.null(X)|all(X==1)){ rho <- log(mean(phaseII)/(1-mean(phaseII))) rhox <- matrix(1,N) rhox2 <- matrix(1,n) Pi <- rep(mean(phaseII),N) }else{ Pimod <- glm(phaseII~X,"binomial") rho <- as.numeric(Pimod$coef) rhox <- model.matrix(Pimod) rhox2 <- rhox[phaseII==1,] Pi <- as.numeric(predict(Pimod,type="response")) } } wt2 <- 1/Pi[phaseII==1] ## obtain initial parameters from coxph() if(is.null(para0)){ lt1 <- c(t(t)) lw <- c(t(w)) lwt <- rep(wt2,each=mk) lid <- rep(1:n,each=mk) ltime <- rep(time,each=mk) lt2 <- c(lt1[-1],tail(ltime,1)) ldata <- data.frame(id=lid,wt=lwt,time1=lt1,time2=lt2,w=lw,event=0,ltime) ldata <- ldata[complete.cases(ldata),] for(i in unique(ldata$id)){ idx <- which(ldata$id==i) ldata$time2[tail(idx,1)] <- time[i] } ldata <- ldata[ldata$time1 != ldata$time2,] for(i in which(event==1)){ idx <- which(ldata$id==i) if(length(idx)>0) ldata$event[tail(idx,1)] <- 1 } if(model=="margin"){ para0 <- coef(coxph(Surv(time1,time2,event)~w,weights=wt,data=ldata)) } if(model=="adj"){ tempz <- data.frame(id=1:n,z=z) ldata.z <- merge(ldata,tempz,by="id") names(ldata.z)[-(1:ncol(ldata))] <- paste("z",1:p,sep="") para0 <- coef(coxph(as.formula(paste("Surv(time1,time2,event)~w+", paste("z",1:p,collapse="+",sep=""))), weights=wt,data=ldata.z)) } if(model=="inter"){ tempz <- data.frame(id=1:n,z=z) ldata.z <- merge(ldata,tempz,by="id") names(ldata.z)[-(1:ncol(ldata))] <- paste("z",1:p,sep="") ldata.z$inter.w.z1 <- ldata.z$w*ldata.z$z1 para0 <- coef(coxph(as.formula(paste("Surv(time1,time2,event)~w+", paste("z",1:p,collapse="+",sep=""), "+inter.w.z1")), weights=wt,data=ldata.z)) } } ## estimate the variance of measurement error sigmae^2 resi <- rep(0,n) resi[which(J>=q)] <- sapply(which(J>=q),function(i){ ki <- !is.na(w[i,])&!is.na(t[i,]) design <- t(sapply(t[i,ki],ftfun)) xls <- drop(solve(crossprod(design))%*%crossprod(design,w[i,ki])) crossprod(w[i,ki]-design%*%xls) }) sigmae <- sqrt(sum(wt2*resi*(J>=q))/sum(wt2*(J-q)*(J>=q))) ## J(u) Jiu <- function(i,u) length(which(t[i,]<=u & !is.na(t[i,]) & !is.na(w[i,]))) Ju <- function(u) sapply(1:n,Jiu,u) ## Y(u) Yiu <- function(i,u) 1*(Jiu(i,u) >= q & time[i]>=u) Yu <- function(u) 1*(Ju(u)>=q & time>=u) ## dN(u) diu <- function(i,u) 1*(Jiu(i,u) >= q & event[i]==1 & time[i]==u) du <- function(u) 1*(Ju(u)>=q & event==1 & time==u) ## least squares estimates and SigmaR SigmaRxlsiu <- function(i,u){ kiu <- which(t[i,]<=u & !is.na(t[i,]) & !is.na(w[i,])) #(i,u) must satisfy length(kiu)>=q ti <- t[i,]; wi <- w[i,]; f <- ftfun(u) diu <- t(sapply(ti[kiu],ftfun)) diuinv <- solve(crossprod(diu)) ss <- drop(sigmae^2*(t(f)%*%diuinv%*%f)) xx <- drop(diuinv %*% crossprod(diu,wi[kiu])) return(list(SigmaR=ss,xls=xx)) } SigmaRxlsu <- function(u){ xx <- matrix(0,n,q) ss <-rep(0,n) for(i in which(Yu(u)==1)){ temp <- SigmaRxlsiu(i,u) ss[i] <- temp$SigmaR xx[i,] <- temp$xls } return(list(SigmaR=ss,xls=xx)) } ## for each subject, get xls and SigmaR at each observed event time orderu <- unique(sort(time[event==1])) xxss <- lapply(orderu,SigmaRxlsu) xlsall <- SigmaRall <- NULL for(i in 1:length(orderu)){ xlsall[[i]] <- xxss[[i]][[2]] SigmaRall[[i]] <- xxss[[i]][[1]] } if(model=="margin"){ ## Q_i(u,theta,sigmae^2) Qiu <- function(i,u,beta){ A <- SigmaRxlsiu(i,u) drop(crossprod(ftfun(u),A$xls) + A$SigmaR*beta*diu(i,u)) } ## E_F^(1)(u,theta,sigmae^2)/E_F^(0)(u,theta,sigmae^2) Eu <- function(u,beta){ YY <- Yu(u) if(sum(YY)==0) return(rep(NA,np)) ii <- which(orderu==u) ss <- SigmaRall[[ii]] xx <- xlsall[[ii]] dN <- du(u) QQ <- YY*drop(xx%*%ftfun(u) + ss*beta*dN) e0 <- YY*exp(QQ*beta-0.5*beta^2*ss) e1 <- YY*e0*QQ sum(e1*wt2,na.rm=T) / sum(e0*wt2,na.rm=T) } ## estimating equations eefun <- function(beta){ ps <- sapply(which(event==1&J>=q),function(i){ wt2[i]*(c(Qiu(i,time[i],beta))-Eu(time[i],beta)) }) sum(ps,na.rm=T) } ## variance Varfun <- function(beta,sigmae,rho){ Eu_VAR <- function(u,beta){ YY <- Yu(u) if(sum(YY)==0) return(list(E0=NA,E1=NA,dE0dpara=NA,dE1dpara=NA, dE0dsigmae2=NA,dE1dsigmae2=NA)) ii <- which(orderu==u) ss <- SigmaRall[[ii]] xx <- xlsall[[ii]] dN <- du(u) xhat <- drop(xx%*%ftfun(u)) QQ <- YY*(xhat + ss*beta*dN) E0 <- YY*exp(QQ*beta-0.5*beta^2*ss) E1 <- YY*E0*QQ dE0dpara <- E0*(xhat+(2*dN-1)*beta*ss) dE1dpara <- QQ*dE0dpara + dN*E0*ss dE0dsigmae2 <- E0*(dN-0.5)*beta^2*ss/sigmae^2 dE1dsigmae2 <- QQ*dE0dsigmae2 + beta*ss*dN*E0/sigmae^2 return(list(E0=sum(wt2*E0,na.rm=T),E1=sum(wt2*E1,na.rm=T), dE0dpara=sum(wt2*dE0dpara,na.rm=T),dE1dpara=sum(wt2*dE1dpara,na.rm=T), dE0dsigmae2=sum(wt2*dE0dsigmae2,na.rm=T),dE1dsigmae2=sum(wt2*dE1dsigmae2,na.rm=T))) } temps <- lapply(orderu,Eu_VAR,beta) # calculate the "meat" part of the sandwich variance estimate lambda0 <- rep(NA,length(orderu)) Mi1 <- matrix(0,n,np) for(j in 1:length(orderu)){ u <- orderu[j] YY <- Yu(u) if(sum(YY)==0) next E0 <- temps[[j]]$E0 E1 <- temps[[j]]$E1 ss <- SigmaRall[[j]] xx <- xlsall[[j]] dN <- du(u) QQ <- YY*drop(xx%*%ftfun(u) + ss*beta*dN) e0 <- YY*exp(QQ*beta-0.5*beta^2*ss) lambda0[j] <- sum(wt2*dN)/E0 tempMi1 <- (QQ - E1/E0)*(dN - YY*e0*lambda0[j]) tempMi1[is.na(tempMi1)] <- 0 Mi1 <- Mi1 + tempMi1 } Mi <- cbind(Mi1,(resi-(J-q)*sigmae^2)*(J>=q)) if(!is.null(rho)) Si <- (phaseII - Pi)*rhox # calculate the "bread" part of the sandwich variance estimate dU1dparadsigmae2 <- matrix(0,np,np+1) for(i in which(event==1)){ if(sum(Yu(time[i]))==0) next j <- which(orderu==time[i]) E0 <- temps[[j]]$E0 E1 <- temps[[j]]$E1 dE0dpara <- temps[[j]]$dE0dpara dE1dpara <- temps[[j]]$dE1dpara dE0dsigmae2 <- temps[[j]]$dE0dsigmae2 dE1dsigmae2 <- temps[[j]]$dE1dsigmae2 si <- SigmaRall[[j]][i] dU1dpara <- -(E0*dE1dpara - tcrossprod(E1,dE0dpara))/E0^2 + si dU1dsigmae2 <- -(E0*dE1dsigmae2 - E1*dE0dsigmae2)/E0^2 + si*beta/sigmae^2 tempdU1 <- cbind(dU1dpara,dU1dsigmae2) tempdU1[is.na(tempdU1)] <- 0 dU1dparadsigmae2 <- dU1dparadsigmae2 + wt2[i]*tempdU1 } dU2dsigmae2 <- -sum(wt2*(J-q)*(J>=q),na.rm=T) dUdparadsigmae2 <- rbind(dU1dparadsigmae2,c(rep(0,np),dU2dsigmae2)) if(is.null(rho)){ bread <- dUdparadsigmae2/N meat <- crossprod(wt2*Mi)/N }else{ bread <- dUdparadsigmae2/N dinvpidrho <- as.matrix(-rhox2/drop(exp(rhox2%*%rho))) EdUdrho <- crossprod(Mi,dinvpidrho)/N EdSdrho <- matrix(apply(-Pi*(1-Pi)*t(apply(rhox,1,tcrossprod)),2,mean,na.rm=T),length(rho),length(rho)) Ri <- matrix(0,N,np+1) Ri[phaseII==1,] <- wt2*Mi bread <- rbind(cbind(bread,EdUdrho),cbind(matrix(0,length(rho),np+1),EdSdrho)) meat <- crossprod(cbind(Ri,Si))/N } bbinverse <- try(solve(bread),T) if(class(bbinverse)=="try-error") varest <- matrix(NA,np+1,np+1) else varest <- (bbinverse%*%meat%*%t(bbinverse)/N)[1:(np+1),1:(np+1)] colnames(varest) <- rownames(varest) <- c("X(u)","sigmae2") lambda0 <- cbind(time=orderu,lambda0=lambda0) return(list(varest=varest,lambda0=lambda0)) } } if(model=="adj"){ Qiu <- function(i,u,beta){ A <- SigmaRxlsiu(i,u) drop(crossprod(ftfun(u),A$xls) + A$SigmaR*beta*diu(i,u)) } Eu <- function(u,beta,eta){ YY <- Yu(u) if(sum(YY)==0) return(rep(NA,np)) ii <- which(orderu==u) ss <- SigmaRall[[ii]] xx <- xlsall[[ii]] dN <- du(u) ez <- drop(exp(z%*%eta)) QQ <- YY*drop(xx%*%ftfun(u) + ss*beta*dN) e0 <- YY*exp(QQ*beta-0.5*beta^2*ss)*ez e1 <- e0*cbind(QQ,z) colSums(e1*wt2,na.rm=T) / sum(e0*wt2,na.rm=T) } eefun <- function(para){ beta <- para[1] eta <- para[-1] ps <- sapply(which(event==1&J>=q),function(i){ wt2[i]*(c(Qiu(i,time[i],beta),z[i,])-Eu(time[i],beta,eta)) }) rowSums(ps,na.rm=T) } Varfun <- function(para,sigmae,rho){ beta <- para[1] eta <- para[-1] Eu_VAR <- function(u,beta,eta){ YY <- Yu(u) if(sum(YY)==0) return(list(E0=NA,E1=rep(NA,np),dE0dpara=rep(NA,np), dE1dpara=rep(NA,np^2), dE0dsigmae2=NA,dE1dsigmae2=rep(NA,np))) ii <- which(orderu==u) ss <- SigmaRall[[ii]] xx <- xlsall[[ii]] dN <- du(u) xhat <- drop(xx%*%ftfun(u)) QQ <- YY*(xhat + ss*beta*dN) Qz <- YY*cbind(QQ,z) ez <- drop(exp(z%*%eta)) E0 <- YY*exp(QQ*beta-0.5*beta^2*ss)*ez E1 <- YY*E0*Qz dE0dpara <- E0*cbind(xhat+(2*dN-1)*ss*beta,z) dE1dpara <- t(sapply(1:n,function(i){ a <- tcrossprod(c(Qz[i,]),dE0dpara[i,]) a[1,1] <- a[1,1] + dN[i]*E0[i]*ss[i] a })) dE0dsigmae2 <- E0*(dN-0.5)*beta^2*ss/sigmae^2 dE1dsigmae2 <- Qz*dE0dsigmae2 + cbind(ss*beta*dN*E0/sigmae^2,matrix(0,n,p)) return(list(E0=sum(wt2*E0,na.rm=T),E1=colSums(wt2*E1,na.rm=T), dE0dpara=colSums(wt2*dE0dpara,na.rm=T),dE1dpara=colSums(wt2*dE1dpara,na.rm=T), dE0dsigmae2=sum(wt2*dE0dsigmae2,na.rm=T),dE1dsigmae2=colSums(wt2*dE1dsigmae2,na.rm=T))) } temps <- lapply(orderu,Eu_VAR,beta,eta) lambda0 <- rep(NA,length(orderu)) Mi1 <- matrix(0,n,np) for(j in 1:length(orderu)){ u <- orderu[j] YY <- Yu(u) if(sum(YY)==0) next E0 <- temps[[j]]$E0 E1 <- temps[[j]]$E1 ss <- SigmaRall[[j]] xx <- xlsall[[j]] dN <- du(u) QQ <- YY*drop(xx%*%ftfun(u) + ss*beta*dN) Qz <- YY*cbind(QQ,z) ez <- drop(exp(z%*%eta)) e0 <- YY*exp(QQ*beta-0.5*beta^2*ss)*ez lambda0[j] <- sum(wt2*dN)/E0 tempMi1 <- sweep(Qz,MARGIN=2,E1/E0,"-")*(dN - YY*e0*lambda0[j]) tempMi1[is.na(tempMi1)] <- 0 Mi1 <- Mi1 + tempMi1 } Mi <- cbind(Mi1,(resi-(J-q)*sigmae^2)*(J>=q)) if(!is.null(rho)) Si <- (phaseII - Pi)*rhox # calculate the derivative of estimating equations w.r.t. (beta,eta,sigmae2,rho) dU1dparadsigmae2 <- matrix(0,np,np+1) for(i in which(event==1)){ if(sum(Yu(time[i]))==0) next j <- which(orderu==time[i]) E0 <- temps[[j]]$E0 E1 <- temps[[j]]$E1 dE0dpara <- temps[[j]]$dE0dpara dE1dpara <- matrix(temps[[j]]$dE1dpara,np,np) dE0dsigmae2 <- temps[[j]]$dE0dsigmae2 dE1dsigmae2 <- temps[[j]]$dE1dsigmae2 si <- SigmaRall[[j]][i] dU1dpara <- -(E0*dE1dpara - tcrossprod(E1,dE0dpara))/E0^2 dU1dpara[1,1] <- dU1dpara[1,1] + si dU1dsigmae2 <- -(E0*dE1dsigmae2 - E1*dE0dsigmae2)/E0^2 dU1dsigmae2[1] <- dU1dsigmae2[1] + si*beta/sigmae^2 tempdU1 <- cbind(dU1dpara,dU1dsigmae2) tempdU1[is.na(tempdU1)] <- 0 dU1dparadsigmae2 <- dU1dparadsigmae2 + wt2[i]*tempdU1 } dU2dsigmae2 <- -sum(wt2*(J-q)*(J>=q),na.rm=T) dUdparadsigmae2 <- rbind(dU1dparadsigmae2,c(rep(0,np),dU2dsigmae2)) if(is.null(rho)){ bread <- dUdparadsigmae2/N meat <- crossprod(wt2*Mi)/N }else{ bread <- dUdparadsigmae2/N dinvpidrho <- as.matrix(-rhox2/drop(exp(rhox2%*%rho))) EdUdrho <- crossprod(Mi,dinvpidrho)/N EdSdrho <- matrix(apply(-Pi*(1-Pi)*t(apply(rhox,1,tcrossprod)),2,mean,na.rm=T),length(rho),length(rho)) Ri <- matrix(0,N,np+1) Ri[phaseII==1,] <- wt2*Mi bread <- rbind(cbind(bread,EdUdrho),cbind(matrix(0,length(rho),np+1),EdSdrho)) meat <- crossprod(cbind(Ri,Si))/N } bbinverse <- try(solve(bread),T) if(class(bbinverse)=="try-error") varest <- matrix(NA,np+1,np+1) else varest <- (bbinverse%*%meat%*%t(bbinverse)/N)[1:(np+1),1:(np+1)] colnames(varest) <- rownames(varest) <- c("X(u)",paste("z",1:p,sep=""),"sigmae2") lambda0 <- cbind(time=orderu,lambda0=lambda0) return(list(varest=varest,lambda0=lambda0)) } } if(model=="inter"){ Qiu <- function(i,u,beta,gamma){ A <- SigmaRxlsiu(i,u) drop(crossprod(ftfun(u),A$xls) + A$SigmaR*(beta+gamma*trt[i])*diu(i,u)) } Eu <- function(u,beta,eta,gamma){ YY <- Yu(u) if(sum(YY)==0) return(rep(NA,np)) ii <- which(orderu==u) ss <- SigmaRall[[ii]] xx <- xlsall[[ii]] bg <- (beta+gamma*trt) dN <- du(u) ez <- drop(exp(z%*%eta)) QQ <- YY*drop(xx%*%ftfun(u) + ss*bg*dN) e0 <- YY*exp(QQ*bg-0.5*bg^2*ss)*ez e1 <- YY*e0*cbind(QQ,z,QQ*trt) colSums(e1*wt2,na.rm=T) / sum(e0*wt2,na.rm=T) } eefun <- function(para){ beta <- para[1] eta <- para[2:(np-1)] gamma <- para[np] ps <- sapply(which(event==1&J>=q),function(i){ qi <- Qiu(i,time[i],beta,gamma) wt2[i]*(c(qi,z[i,],qi*trt[i])-Eu(time[i],beta,eta,gamma)) }) rowSums(ps,na.rm=T) } Varfun <- function(para,sigmae,rho){ beta <- para[1] eta <- para[2:(p+1)] gamma <- para[p+2] bg <- (beta+gamma*trt) Eu_VAR <- function(u,beta,eta,gamma){ YY <- Yu(u) if(sum(YY)==0) return(list(E0=NA,E1=rep(NA,np),dE0dpara=rep(NA,np), dE1dpara=rep(NA,np^2), dE0dsigmae2=NA,dE1dsigmae2=rep(NA,np))) ii <- which(orderu==u) ss <- SigmaRall[[ii]] xx <- xlsall[[ii]] dN <- du(u) bg <- (beta+gamma*trt) xhat <- drop(xx%*%ftfun(u)) QQ <- YY*(xhat + ss*bg*dN) Qz <- YY*cbind(QQ,z,QQ*trt) ez <- drop(exp(z%*%eta)) E0 <- YY*exp(QQ*bg-0.5*bg^2*ss)*ez E1 <- YY*E0*Qz dE0dpara <- E0*cbind(xhat+(2*dN-1)*ss*bg, z, xhat*trt+(2*dN-1)*ss*bg*trt) dE1dpara <- t(sapply(1:n,function(i){ tcrossprod(Qz[i,],dE0dpara[i,]) })) dE1dpara<- dE1dpara + E0*ss*dN*cbind(1,matrix(0,n,p),trt,matrix(0,n,np*p),trt,matrix(0,n,p),trt^2) dE0dsigmae2 <- E0*(dN-0.5)*bg^2*ss/sigmae^2 dE1dsigmae2 <- Qz*dE0dsigmae2 + E0*ss*bg*dN/sigmae^2*cbind(1,matrix(0,n,p),trt) return(list(E0=sum(wt2*E0,na.rm=T),E1=colSums(wt2*E1,na.rm=T), dE0dpara=colSums(wt2*dE0dpara,na.rm=T),dE1dpara=colSums(wt2*dE1dpara,na.rm=T), dE0dsigmae2=sum(wt2*dE0dsigmae2,na.rm=T),dE1dsigmae2=colSums(wt2*dE1dsigmae2,na.rm=T))) } temps <- lapply(orderu,Eu_VAR,beta,eta,gamma) lambda0 <- rep(NA,length(orderu)) Mi1 <- matrix(0,n,np) for(j in 1:length(orderu)){ u <- orderu[j] YY <- Yu(u) if(sum(YY)==0) next E0 <- temps[[j]]$E0 E1 <- temps[[j]]$E1 ss <- SigmaRall[[j]] xx <- xlsall[[j]] dN <- du(u) QQ <- YY*drop(xx%*%ftfun(u) + ss*bg*dN) Qz <- YY*cbind(QQ,z,QQ*trt) ez <- drop(exp(z%*%eta)) e0 <- YY*exp(QQ*bg-0.5*bg^2*ss)*ez lambda0[j] <- sum(dN)/E0 tempMi1 <- sweep(Qz,MARGIN=2,E1/E0,"-")*(dN - YY*e0*lambda0[j]) tempMi1[is.na(tempMi1)] <- 0 Mi1 <- Mi1 + tempMi1 } Mi <- cbind(Mi1,(resi-(J-q)*sigmae^2)*(J>=q)) if(!is.null(rho)) Si <- (phaseII - Pi)*rhox dU1dparadsigmae2 <- matrix(0,np,np+1) for(i in which(event==1)){ if(sum(Yu(time[i]))==0) next j <- which(orderu==time[i]) E0 <- temps[[j]]$E0 E1 <- temps[[j]]$E1 dE0dpara <- temps[[j]]$dE0dpara dE1dpara <- matrix(temps[[j]]$dE1dpara,np,np) dE0dsigmae2 <- temps[[j]]$dE0dsigmae2 dE1dsigmae2 <- temps[[j]]$dE1dsigmae2 si <- SigmaRall[[j]][i] dU1dpara <- -(E0*dE1dpara - tcrossprod(E1,dE0dpara))/E0^2 dqzidpara <- c(si,rep(0,p),trt[i]*si) dU1dpara[1,] <- dU1dpara[1,] + dqzidpara dU1dpara[p+2,] <- dU1dpara[p+2,] + trt[i]*dqzidpara dU1dsigmae2 <- -(E0*dE1dsigmae2 - E1*dE0dsigmae2)/E0^2 + c(1,rep(0,p),trt[i])*si*(bg[i])/sigmae^2 tempdU1 <- cbind(dU1dpara,dU1dsigmae2) tempdU1[is.na(tempdU1)] <- 0 dU1dparadsigmae2 <- dU1dparadsigmae2 + wt2[i]*tempdU1 } dU2dsigmae2 <- -sum(wt2*(J-q)*(J>=q),na.rm=T) dUdparadsigmae2 <- rbind(dU1dparadsigmae2,c(rep(0,np),dU2dsigmae2)) if(is.null(rho)){ bread <- dUdparadsigmae2/N meat <- crossprod(wt2*Mi)/N }else{ bread <- dUdparadsigmae2/N dinvpidrho <- as.matrix(-rhox2/drop(exp(rhox2%*%rho))) EdUdrho <- crossprod(Mi,dinvpidrho)/N EdSdrho <- matrix(apply(-Pi*(1-Pi)*t(apply(rhox,1,tcrossprod)),2,mean,na.rm=T),length(rho),length(rho)) Ri <- matrix(0,N,np+1) Ri[phaseII==1,] <- wt2*Mi bread <- rbind(cbind(bread,EdUdrho),cbind(matrix(0,length(rho),np+1),EdSdrho)) meat <- crossprod(cbind(Ri,Si))/N } bbinverse <- try(solve(bread),T) if(class(bbinverse)=="try-error") varest <- matrix(NA,np+1,np+1) else varest <- (bbinverse%*%meat%*%t(bbinverse)/N)[1:(np+1),1:(np+1)] colnames(varest) <- rownames(varest) <- c("X(u)",paste("z",1:p,sep=""),"X(u)*z1","sigmae2") lambda0 <- cbind(time=orderu,lambda0=lambda0) return(list(varest=varest,lambda0=lambda0)) } } ## estimate the parameters paraest <- multiroot(eefun,start=para0) ## estimate the variance matrix vest <- Varfun(paraest$root,sigmae,rho) ## output covmtx <- vest$varest[1:np,1:np] coef <- c(paraest$root) names(coef) <- colnames(covmtx) #bhazard <- vest$lambda0 error <- c(sigmae^2,sqrt(vest$varest[np+1,np+1])) names(error) <- c("Est.","Std.Err.") convergence <- list(f.root=paraest$f.root,iter=paraest$iter,estim.precis=paraest$estim.precis) names(convergence$f.root) <- names(coef) return(list(coef=coef,covmtx=covmtx,merr=error, convergence=convergence)) } ########################################################################################## ### CNSBernoulliAIPWnp() ### Description : Function to implement AIPW conditional score estimator on two-phase data from Bernoulli sampling. ### Note : The sample size of full cohort data is N and the total number of subjects sampled ### into the phaseII sample is n. The number of rows of the matrices ### and the length of vectors used in following arguments must all be N. Values for subjects not sampled for ### immune response measurements are NA. The maximum number of immune response measurements per subject is mk. ### Arguments : ### t : N*mk matrix of time points to make immune response measurements. ### w : N*mk matrix of immune response values. ### z : If model=="margin", by default = NULL. If model=="adj", N*p matrix of time-independent variables. If method=="inter", ### N*p matrix of time-independent and the first column must be a vector of binary 0/1 indicators. ### time : N-dimensional vector of follow up time. ### event : N-dimensional vector of event status, 0=non-event, 1=event. ### model : which form of Cox regression is fitted. If "margin", exp(X(u)*beta). If "adj", exp(X(u)*beta+z^Teta). If "inter", ### exp(X(u)*beta+z^Teta+X(u)*z[,1]*gamma), and the first column of z must be 0/1 treatment indicators. ### Pi : N-dimensional vector of pre-specified sampling probabilities, if ==NULL, estimated from X (see below). ### phaseII : N-dimensional vector of phaseII sample indicators, 0=not a phaseII sample, 1= a phaseII sample. ### X : matrix of variables to estimate the sampling probabilities if Pi is Null. If simple random sampling, then ### leave X=NULL or X=matrix(1,N,1), and sampling probabilites are estimated as n/N. ### W : matrix of variables used in the augmentation term to estimate the conditional expectation E[.|W] ### H : diagonal bandwidth matrix ### para0 : vector of initial parameters for (beta,eta^T,gamma)^T. Default is calculated from ### coxph() on observed immune responses. ### ftfun : function specifies the immune response trajectory, X(u) = alpha^Tftfun(u). ### Default is linear in time, i.e. ftfun = function(u) c(1,u). ### maxiter : maximal number of iterations allowed, see multiroot(). ### Value : the function returns a list including ### coef : estimated Cox model coefficients. ### covmtx : estimated covariance matrix for coef. ### merr : estimated variance of measurement error and its standard error. ### convergence : convergence parameters. Output from multiroot(). CNSBernoulliAIPWnp <- function(t,w,z=NULL,time,event,model=c("margin","adj","inter"),phaseII,Pi=NULL,W, X=NULL,H=NULL,para0=NULL,ftfun=NULL,maxiter=50){ if(is.null(ftfun)) ftfun <- function(uu) return(c(1,uu)) q <- length(ftfun(1)) N <- length(phaseII) if(is.null(W)) stop("Missing W!") W <- as.matrix(W) #if(!is.null(z)) W <- cbind(W,z) else W <- as.matrix(W) sd2 <- function(x) sqrt(var(x,na.rm=T)*(length(x)-1)/length(x)) nW <- dim(W)[2] if(is.null(H)){ h <- sapply(1:nW,function(i) 0.75*sd2(W[,i])*N^(-1/3)) H <- diag(h,nW) } t[phaseII==0,] <- NA w[phaseII==0,] <- NA J <- apply(w,1,function(wi){sum(!is.na(wi))}) if(all(J<=q)==1) stop("Need to have at least q measurements to fit the model!") p <- ifelse(model=="margin",0,ncol(z)) if(model=="inter") trt <- z[,1] np <- ifelse(model=="margin",1,ifelse(model=="adj",1+p,2+p)) # unknown sampling prob if(is.null(Pi)){ if(is.null(X)) stop("Missing X and Pi!") Pimod <- glm(phaseII~X,"binomial") rho <- as.numeric(Pimod$coef) rhox <- model.matrix(Pimod) Pi <- as.numeric(predict(Pimod,type="response")) }else{ rho <- NULL } if(is.null(para0)){ mk <- ncol(t) lt <- c(t(t)) lw <- c(t(w)) lwt <- rep(1/Pi,each=mk) lid <- rep(1:N,each=mk) lphase2 <- rep(phaseII,each=mk) ltime <- rep(time,each=mk) ldata <- data.frame(lid,lphase2,lwt,lt,lw,ltime,levent=0) ldata <- ldata[complete.cases(ldata),] ldata <- ldata[ldata$lphase2==1,] lt2 <- c(ldata[-1,"lt"],ldata[dim(ldata)[1],"ltime"]) ldata$lt2 <- lt2 for(i in 1:N){ idx <- which(ldata[,"lid"]==i) if(length(idx)>0) ldata[tail(idx,1),"lt2"] <- time[i] } ldata <- ldata[ldata$lt != ldata$lt2,] eid <- which(event==1) for(i in eid){ idx <- which(ldata[,"lid"]==i) if(length(idx)>0) ldata[tail(idx,1),"levent"] <- 1 } if(model=="margin"){ para0 <- coef(coxph(Surv(lt,lt2,levent)~lw,data=ldata,weights=lwt)) } if(model=="adj"){ wz <- data.frame(lid=1:N,z=z) ldata.z <- merge(ldata,wz,by="lid") names(ldata.z)[-(1:8)] <- paste("z.",1:p,sep="") para0 <- coef(coxph(as.formula(paste("Surv(lt,lt2,levent)~lw+",paste("z.",1:p,collapse="+",sep=""))), data=ldata.z,weights=lwt)) } if(model=="inter"){ wz <- data.frame(lid=1:N,z=z) ldata.z <- merge(ldata,wz,by="lid") names(ldata.z)[-(1:8)] <- paste("z.",1:p,sep="") ldata.z$inter <- ldata.z$lw*ldata.z$z.1 para0 <- coef(coxph(as.formula(paste("Surv(lt,lt2,levent)~lw+",paste("z.",1:p,collapse="+",sep=""),"+inter")), data=ldata.z,weights=lwt)) } } ## estimate sigmae, phaseII smaple resi <- rep(0,N) resi[which(J>=q&phaseII==1)] <- sapply(which(J>=q&phaseII==1),function(i){ ji <- !is.na(w[i,])&!is.na(t[i,]) design <- t(sapply(t[i,ji],ftfun)) xls <- drop(solve(crossprod(design))%*%crossprod(design,w[i,ji])) drop(crossprod(w[i,ji]-design%*%xls)) }) kq <- J>=q ai <- kq*resi bi <- kq*(J-q) Eai <- npreg(tydat=ai[phaseII==1],txdat=W[phaseII==1,],exdat=W,bws=diag(H))$mean Ebi <- npreg(tydat=bi[phaseII==1],txdat=W[phaseII==1,],exdat=W,bws=diag(H))$mean sigmae <- sqrt(sum(phaseII/Pi*ai + (1-phaseII/Pi)*Eai) / sum(phaseII/Pi*bi + (1-phaseII/Pi)*Ebi)) sei <- phaseII/Pi*(ai-bi*sigmae^2) + (1-phaseII/Pi)*(Eai-Ebi*sigmae^2) ## number of measurements up to u Jiu <- function(i,u) length(which(t[i,]<=u & !is.na(t[i,]) & !is.na(w[i,]))) Ju <- function(u) sapply(1:N,Jiu,u) ## at risk at u Yiu <- function(i,u) 1*(Jiu(i,u) >= q & time[i]>=u) Yu <- function(u) 1*(Ju(u)>=q & time>=u) ## event process at u diu <- function(i,u) 1*(Jiu(i,u) >= q & event[i]==1 & time[i]==u) du <- function(u) 1*(Ju(u)>=q & event==1 & time==u) ## preparation: xls and SigmaR SigmaRxlsiu <- function(i,u){ Jiu <- which(t[i,]<=u & !is.na(t[i,]) & !is.na(w[i,])) #(i,u) must satisfy length(Jiu)>=q ti <- t[i,]; wi <- w[i,]; f <- ftfun(u) diu <- t(sapply(ti[Jiu],ftfun)) diuinv <- solve(crossprod(diu)) ss <- drop(sigmae^2*(t(f)%*%diuinv%*%f)) xx <- drop(diuinv %*% crossprod(diu,wi[Jiu])) return(list(SigmaR=ss,xls=xx)) } SigmaRxlsu <- function(u){ xx <- matrix(0,N,q) ss <- rep(0,N) for(i in which(Ju(u)>=q&phaseII==1)){ temp <- SigmaRxlsiu(i,u) ss[i] <- temp$SigmaR xx[i,] <- temp$xls } return(list(SigmaR=ss,xls=xx)) } ## orderuc orderu <- unique(sort(time[which(event==1)])) xxss <- lapply(orderu,SigmaRxlsu) xlsall <- SigmaRall <- NULL for(i in 1:length(orderu)){ xlsall[[i]] <- xxss[[i]][[2]] SigmaRall[[i]] <- xxss[[i]][[1]] } ker <- sapply(which(phaseII==1),function(i){ dmvnorm(W, mean=W[i,],sigma=H) }) if(model=="margin"){ Qiu <- function(i,u,beta){ A <- SigmaRxlsiu(i,u) drop(crossprod(ftfun(u),A$xls) + A$SigmaR*beta*diu(i,u)) } Ebar_aipw <- function(u,beta){ YY <- Yu(u) KK <- Ju(u) #RR <- 1*(KJ>=q) RR <- 1*(YY==1&phaseII==1) if(sum(RR)==0) return(list(NA)) ii <- which(orderu==u) ss <- SigmaRall[[ii]] xx <- xlsall[[ii]] dN <- du(u) QQ <- RR*drop(xx%*%ftfun(u) + ss*beta*dN) e0q <- RR*exp(QQ*beta-0.5*beta^2*ss) e1q <- RR*QQ*e0q use <- 1*(e0q!=Inf) ker0 <- ker[,RR[phaseII==1]==1] yobs <- cbind(QQ,e0q,e1q)[RR==1,] jj <- !is.na(yobs[,2]) & yobs[,2]!=Inf if(all(!jj)) return(list(NA)) #jj <- rep(T,sum(phaseII)) ## conditional expectation given complete data num <- ker0[,jj]%*%yobs[jj,] den <- rowSums(ker0[,jj],na.rm=T) newy <- num/den newy[is.na(newy)] <- 0 Eq <- newy[,1] Ee0q <- newy[,2] Ee1q <- newy[,3] e0 <- (phaseII/Pi*e0q + (1-phaseII/Pi)*Ee0q)*(time>=u)*use e1 <- (phaseII/Pi*e1q + (1-phaseII/Pi)*Ee1q)*(time>=u)*use qz <- (phaseII/Pi*QQ + (1-phaseII/Pi)*Eq)*(time>=u)*use return(list(ebar=sum(e1,na.rm=T)/sum(e0,na.rm=T),qz=qz)) } estfun <- function(beta){ ps <- sapply(which(event==1),function(i){ temp <- Ebar_aipw(time[i],beta) if(all(is.na(temp[[1]]))) return(0) (temp$qz[i] - temp$ebar) }) return(sum(ps,na.rm=T)) } } if(model=="adj"){ Qiu <- function(i,u,beta){ A <- SigmaRxlsiu(i,u) drop(crossprod(ftfun(u),A$xls) + A$SigmaR*beta*diu(i,u)) } Ebar_aipw <- function(u,beta,eta){ YY <- Yu(u) KK <- Ju(u) RR <- 1*(YY==1 & phaseII==1) if(sum(RR)==0) return(list(NA)) ii <- which(orderu==u) ss <- SigmaRall[[ii]] xx <- xlsall[[ii]] dN <- du(u) QQ <- RR*drop(xx%*%ftfun(u) + ss*beta*dN) e0q <- RR*exp(QQ*beta-0.5*beta^2*ss) e1q <- RR*QQ*e0q use <- 1*(e0q!=Inf) ker0 <- ker[,RR[phaseII==1]==1] yobs <- cbind(QQ,e0q,e1q)[RR==1,] jj <- !is.na(yobs[,2]) & yobs[,2]!=Inf if(all(!jj)) return(list(NA)) #jj <- rep(T,sum(phaseII)) ## conditional expectation given complete data num <- ker0[,jj]%*%yobs[jj,] den <- rowSums(ker0[,jj],na.rm=T) newy <- num/den newy[is.na(newy)] <- 0 Eq <- newy[,1] Ee0q <- newy[,2] Ee1q <- newy[,3] ez <- drop(exp(z%*%eta)) e0q2 <- (phaseII/Pi*e0q + (1-phaseII/Pi)*Ee0q)*(time>=u)*use e1q2 <- (phaseII/Pi*e1q + (1-phaseII/Pi)*Ee1q)*(time>=u)*use q2 <- (phaseII/Pi*QQ + (1-phaseII/Pi)*Eq)*(time>=u)*use e0 <- e0q2*ez e1 <- cbind(e1q2,z*e0q2)*ez qz <- cbind(q2,z) return(list(ebar=colSums(e1,na.rm=T)/sum(e0,na.rm=T),qz=qz)) } estfun <- function(para){ beta <- para[1] eta <- para[-1] ps <- sapply(which(event==1),function(i){ temp <- Ebar_aipw(time[i],beta,eta) if(all(is.na(temp[[1]]))) return(rep(0,np)) (temp$qz[i,] - temp$ebar) }) rowSums(ps,na.rm=T) } } if(model=="inter"){ Qiu <- function(i,u,beta,gamma){ A <- SigmaRxlsiu(i,u) drop(crossprod(ftfun(u),A$xls) + A$SigmaR*(beta+gamma*trt[i])*diu(i,u)) } Ebar_aipw <- function(u,beta,eta,gamma){ YY <- Yu(u) KK <- Ju(u) RR <- 1*(YY==1 & phaseII==1) if(sum(RR)==0) return(list(NA)) ii <- which(orderu==u) ss <- SigmaRall[[ii]] xx <- xlsall[[ii]] dN <- du(u) bg <- (beta+gamma*trt) QQ <- RR*drop(xx%*%ftfun(u) + ss*bg*dN) e0q <- RR*exp(QQ*bg-0.5*bg^2*ss) e1q <- RR*QQ*e0q use <- 1*(e0q!=Inf) ker0 <- ker[,RR[phaseII==1]==1] yobs <- cbind(QQ,e0q,e1q)[RR==1,] jj <- !is.na(yobs[,2]) & yobs[,2]!=Inf if(all(!jj)) return(list(NA)) #jj <- rep(T,sum(phaseII)) ## conditional expectation given complete data num <- ker0[,jj]%*%yobs[jj,] den <- rowSums(ker0[,jj],na.rm=T) newy <- num/den newy[is.na(newy)] <- 0 Eq <- newy[,1] Ee0q <- newy[,2] Ee1q <- newy[,3] ez <- drop(exp(z%*%eta)) e0q2 <- (phaseII/Pi*e0q + (1-phaseII/Pi)*Ee0q)*(time>=u)*use e1q2 <- (phaseII/Pi*e1q + (1-phaseII/Pi)*Ee1q)*(time>=u)*use q2 <- (phaseII/Pi*QQ + (1-phaseII/Pi)*Eq)*(time>=u)*use e0 <- e0q2*ez e1 <- cbind(e1q2,z*e0q2,e1q2*trt)*ez qz <- cbind(q2,z,q2*trt) return(list(ebar=colSums(e1,na.rm=T)/sum(e0,na.rm=T),qz=qz)) } estfun <- function(para){ beta <- para[1] eta <- para[2:(p+1)] gamma <- para[np] ps <- sapply(which(event==1),function(i){ temp <- Ebar_aipw(time[i],beta,eta,gamma) if(all(is.na(temp[[1]]))) return(rep(0,np)) (temp$qz[i,] - temp$ebar) }) rowSums(ps,na.rm=T) } } Varfun <- function(para,sigmae,rho){ if(model=="margin"){ beta <- para bg <- beta } if(model=="adj"){ beta <- para[1] eta <- para[-1] bg <- beta } if(model=="inter"){ beta <- para[1] eta <- para[2:(p+1)] gamma <- para[np] bg <- beta+gamma*trt } Ebar_VAR <- function(u){ YY <- Yu(u) KK <- Ju(u) RR <- 1*(YY==1 & phaseII==1) if(sum(RR)==0) return(list(NA)) ii <- which(orderu==u) ss <- SigmaRall[[ii]] xx <- xlsall[[ii]] dN <- du(u) xhat <- drop(xx%*%ftfun(u)) QQ <- RR*(xhat + ss*bg*dN) e0q <- RR*exp(QQ*bg-0.5*bg^2*ss) e1q <- RR*QQ*e0q dqdbeta <- ss*dN dqdsigmae2 <- ss*bg*dN/sigmae^2 de0qdbeta <- e0q*(xhat+(2*dN-1)*ss*bg) de0qdsigmae2 <- e0q*ss*bg^2*(dN-0.5)/sigmae^2 de1qdbeta <- QQ*de0qdbeta + dqdbeta*e0q de1qdsigmae2 <- QQ*de0qdsigmae2 + e0q*dqdsigmae2 yobs <- cbind(QQ,e0q,e1q,RR,de0qdbeta,de1qdbeta, de0qdsigmae2,de1qdsigmae2, dqdbeta,dqdsigmae2)[RR==1,] ker0 <- ker[,RR[phaseII==1]==1] jj <- !is.na(yobs[,2]) & yobs[,2]!=Inf if(all(!jj)) return(list(NA)) #jj <- rep(T,sum(phaseII)) ## conditional expectation given complete data num <- ker0[,jj]%*%yobs[jj,] den <- rowSums(ker0[,jj],na.rm=T) newy <- num/den newy[is.na(newy)] <- 0 Eq <- newy[,1] Ee0q <- newy[,2] Ee1q <- newy[,3] Ek <- newy[,4] Ede0qdbeta <- newy[,5] Ede1qdbeta <- newy[,6] Ede0qdsigmae2 <- newy[,7] Ede1qdsigmae2 <- newy[,8] Edqdbeta <- newy[,9] Edqdsigmae2 <- newy[,10] use <- 1*(e0q!=Inf) q2 <- (phaseII/Pi*QQ + (1-phaseII/Pi)*Eq)*(time>=u)*use #rr <- (phaseII/Pi*RR + (1-phaseII/Pi)*Ek)*(time>=u)*use e0q2 <- (phaseII/Pi*e0q + (1-phaseII/Pi)*Ee0q)*(time>=u)*use e1q2 <- (phaseII/Pi*e1q + (1-phaseII/Pi)*Ee1q)*(time>=u)*use de0qdbeta2 <- (phaseII/Pi*de0qdbeta + (1-phaseII/Pi)*Ede0qdbeta)*(time>=u)*use de1qdbeta2 <- (phaseII/Pi*de1qdbeta + (1-phaseII/Pi)*Ede1qdbeta)*(time>=u)*use de0qdsigmae22 <- (phaseII/Pi*de0qdsigmae2 + (1-phaseII/Pi)*Ede0qdsigmae2)*(time>=u)*use de1qdsigmae22 <- (phaseII/Pi*de1qdsigmae2 + (1-phaseII/Pi)*Ede1qdsigmae2)*(time>=u)*use dqdbeta2 <- (phaseII/Pi*dqdbeta + (1-phaseII/Pi)*Edqdbeta)*(time>=u)*use dqdsigmae22 <- (phaseII/Pi*dqdsigmae2 + (1-phaseII/Pi)*Edqdsigmae2)*(time>=u)*use if(model=="margin"){ qz <- cbind(q2) e0 <- e0q2 e1 <- cbind(e1q2) de0dpara <- sum(de0qdbeta2,na.rm=T) de1dpara <- sum(de1qdbeta2,na.rm=T) de0dsigmae2 <- sum(de0qdsigmae22,na.rm=T) de1dsigmae2 <- sum(de1qdsigmae22,na.rm=T) } if(model=="adj"){ ez <- drop(exp(z%*%eta)) qz <- cbind(q2,z) e0 <- e0q2*ez e1 <- cbind(e1q2,z*e0q2)*ez if(p==1) z2 <- z*z else z2 <- t(apply(z,1,tcrossprod)) de0dpara <- colSums(ez*cbind(de0qdbeta2,e0q2*z),na.rm=T) de1dpara11 <- sum(ez*de1qdbeta2,na.rm=T) # 1*1 de1dpara12 <- matrix(colSums(ez*e1q2*z,na.rm=T),1,p,byrow=T) # 1*p de1dpara21 <- matrix(colSums(ez*de0qdbeta2*z,na.rm=T),p,1) # p*1 de1dpara22 <- matrix(colSums(ez*e0q2*z2,na.rm=T),p,p) #p*p de1dpara <- cbind(rbind(de1dpara11,de1dpara21),rbind(de1dpara12,de1dpara22)) de0dsigmae2 <- sum(ez*de0qdsigmae22,na.rm=T) de1dsigmae2 <- colSums(ez*cbind(de1qdsigmae22,z*de0qdsigmae22),na.rm=T) } if(model=="inter"){ ez <- drop(exp(z%*%eta)) qz <- cbind(q2,z,q2*trt) e0 <- e0q2*ez e1 <- cbind(e1q2,z*e0q2,e1q2*trt)*ez if(p==1) z2 <- z*z else z2 <- t(apply(z,1,tcrossprod)) de0dpara <- colSums(ez*cbind(de0qdbeta2,e0q2*z,de0qdbeta2*trt),na.rm=T) de1dpara11 <- sum(ez*de1qdbeta2,na.rm=T) # 1*1 de1dpara12 <- matrix(colSums(ez*e1q2*z,na.rm=T),1,p,byrow=T) # 1*p de1dpara13 <- sum(ez*de1qdbeta2*trt,na.rm=T) # 1*1 de1dpara21 <- matrix(colSums(ez*de0qdbeta2*z,na.rm=T),p,1) # p*1 de1dpara22 <- matrix(colSums(ez*e0q2*z2,na.rm=T),p,p) #p*p de1dpara23 <- matrix(colSums(ez*de0qdbeta2*z*trt,na.rm=T),p,1) #p*1 de1dpara31 <- sum(ez*de1qdbeta2*trt,na.rm=T) #1*1 de1dpara32 <- matrix(colSums(ez*e1q2*trt*z,na.rm=T),1,p,byrow=T) #1*p de1dpara33 <- sum(ez*de1qdbeta2*trt*trt,na.rm=T) #1*1 de1dpara <- cbind(rbind(de1dpara11,de1dpara21,de1dpara31),rbind(de1dpara12,de1dpara22,de1dpara32), rbind(de1dpara13,de1dpara23,de1dpara33)) de0dsigmae2 <- sum(ez*de0qdsigmae22,na.rm=T) de1dsigmae2 <- colSums(ez*cbind(de1qdsigmae22,z*de0qdsigmae22,trt*de1qdsigmae22),na.rm=T) } return(list(use=use,ebar=colMeans(e1,na.rm=T)/mean(e0,na.rm=T),e1=e1,e0=e0,qz=qz, de0dpara=de0dpara,de1dpara=de1dpara, de0dsigmae2=de0dsigmae2,de1dsigmae2=de1dsigmae2, dqdbeta=dqdbeta2,dqdsigmae2=dqdsigmae22)) } temps <- lapply(orderu,Ebar_VAR) # baseline hazard lambda0 <- rep(NA,length(orderu)) # calculate M1,M2 Mi1 <- matrix(0,N,np) #totalrr <- rep(0,N) for(j in 1:length(orderu)){ u <- orderu[j] YY <- Yu(u) YYstar <- 1*(time>=u) dN <- du(u) dNstar <- 1*(time==u&event==1) KK <- Ju(u) RR <- 1*(YY==1 & phaseII==1) if(sum(RR)==0) next ss <- SigmaRall[[j]] xx <- xlsall[[j]] tj <- temps[[j]] use <- tj$use==1 lambda0u <- mean(dN)/mean(tj$e0[use],na.rm=T) lambda0[j] <- lambda0u #print(sum((QQ-tj$ebar)*(dN-YY*lambda0u*exp(QQ*bg-0.5*bg^2*ss))>0)) dN2 <- phaseII/Pi*dN + (1-phaseII/Pi)*(dNstar) if(any(is.na(tj$ebar[[1]]))) next Mi1[use,] <- Mi1[use,] + (dNstar*tj$qz - YYstar*lambda0u*tj$e1 - matrix(tj$ebar,N,np,byrow=T)*(dN2 - YYstar*lambda0u*tj$e0))[use,] #totalrr <- totalrr + use if(any(is.na(Mi1))) stop } Mi <- cbind(Mi1,sei) # calculate the dM1dparadsigmae2 dM1dparadsigmae2 <- matrix(0,np,np+1) for(i in which(event==1)){ if(sum(Yu(time[i]))==0) next j <- which(orderu==time[i]) e0 <- sum(temps[[j]]$e0,na.rm=T) e1 <- colSums(temps[[j]]$e1,na.rm=T) de0dpara <- temps[[j]]$de0dpara de1dpara <- temps[[j]]$de1dpara de0dsigmae2 <- temps[[j]]$de0dsigmae2 de1dsigmae2 <- temps[[j]]$de1dsigmae2 dqdbeta <- temps[[j]]$dqdbeta dqdsigmae2 <- temps[[j]]$dqdsigmae2 si <- SigmaRall[[j]][i] debardpara <- (e0*de1dpara - tcrossprod(e1,de0dpara))/e0^2 debardsigmae2 <- (e0*de1dsigmae2 - e1*de0dsigmae2)/e0^2 #if(is.na(-debardpara)) print(i) if(model=="margin"){ dM1dpara <- -debardpara + dqdbeta[i] dM1dsigmae2 <- -debardsigmae2 + dqdsigmae2[i] vvname <- c("X(u)","sigmae2") } if(model=="adj"){ dM1dpara <- -debardpara dM1dpara[1,1] <- dM1dpara[1,1] + dqdbeta[i] dM1dsigmae2 <- -debardsigmae2 dM1dsigmae2[1] <- dM1dsigmae2[1] + dqdsigmae2[i] vvname <- c("X(u)",paste("z",1:p,sep=""),"sigmae2") } if(model=="inter"){ dM1dpara <- -debardpara dqidpara <- c(dqdbeta[i],rep(0,p),dqdbeta[i]*trt[i]) dM1dpara[1,] <- dM1dpara[1,] + dqidpara dM1dpara[p+2,] <- dM1dpara[p+2,] + trt[i]*dqidpara dM1dsigmae2 <- -debardsigmae2 dM1dsigmae2 <- dM1dsigmae2 + c(dqdsigmae2[i],rep(0,p),trt[i]*dqdsigmae2[i]) vvname <- c("X(u)",paste("z",1:p,sep=""),"X(u)*z1","sigmae2") } tempdM1 <- cbind(dM1dpara,dM1dsigmae2) tempdM1[is.na(tempdM1)] <- 0 dM1dparadsigmae2 <- dM1dparadsigmae2 + tempdM1 } dM2dsigmae2 <- -sum(phaseII/Pi*bi+(1-phaseII/Pi)*Ebi,na.rm=T) dMdparadsigmae2 <- rbind(dM1dparadsigmae2,c(rep(0,np),dM2dsigmae2)) bread <- dMdparadsigmae2/N meat <- crossprod(Mi)/N bbinverse <- try(solve(bread),T) #diag(sqrt(bbinverse%*%meat%*%t(bbinverse)/N)) if(class(bbinverse)=="try-error") varest <- matrix(NA,np+1,np+1) else varest=(bbinverse%*%meat%*%t(bbinverse)/N)[1:(np+1),1:(np+1)] colnames(varest) <- rownames(varest) <- vvname lambda0 <- cbind(time=orderu,lambda0=lambda0) return(list(varest=varest,lambda0=lambda0)) } ## estimate the parameters paraest <- multiroot(estfun,start=para0,maxiter=maxiter) ## estimate the variance matrix vest <- Varfun(paraest$root,sigmae,rho) ## output covmtx <- vest$varest[1:np,1:np] coef <- c(paraest$root) names(coef) <- colnames(covmtx) #bhazard <- vest$lambda0 error <- c(sigmae^2,sqrt(vest$varest[np+1,np+1])) names(error) <- c("Est.","Std.Err.") convergence <- list(f.root=paraest$f.root,iter=paraest$iter,estim.precis=paraest$estim.precis) names(convergence$f.root) <- names(coef) return(list(coef=coef,covmtx=covmtx,merr=error, convergence=convergence)) }