R/isoph.ti.R

Defines functions isoph.ti

isoph.ti=function(TIME, STATUS, Z, X, shape, K, maxdec, maxiter, eps){
  #the comp speed can be further improved by separating shape, inc or dec
  
  #sorted by z
  n=length(STATUS)
  data=data.frame(TIME,STATUS,Z)
  if(!is.null(X)){
    data=data.frame(TIME,STATUS,Z,X)
    nparm=1; if(!is.vector(X)) nparm=ncol(X)
    if(nparm>1){; NR.ft=NR2.ft
    }else{;       NR.ft=NR1.ft
    }
  }
  data=data[order(data$Z),]
  
  if(shape=='increasing'){
    RPA_ti=RPA_ti_inc;
    BTF=BTF_inc;
  }else{
    RPA_ti=RPA_ti_dec;
    BTF=BTF_dec;
  }
  
  #remove subjects whose cov is less than z^*_(1)
  n.total=n
  if(shape=='increasing'){
    if(sum(data$STATUS[data$Z==data$Z[1]])==0){
      idx=which(data$STATUS==1)[1]
      idx2=idx:n;      n=length(idx2)
    }
  }else{
    if(sum(data$STATUS[data$Z==data$Z[n]])==0){
      idx=tail(which(data$STATUS==1),1)
      idx2=1:idx;      n=length(idx2)
    }
  }
  if(n.total>n)  data=data[idx2,]
  
  #re-define variable
  t=data$TIME
  status=data$STATUS
  z=data$Z
  if(!is.null(X)){
    if(nparm==1){ x=as.numeric(data[,4])
    }else{;       x=as.matrix(data[,-c(1,2,3)])
    }
  }
  
  z.obs=unique(z[status==1])
  m=length(z.obs)
  
  t.obs=sort(unique(t))
  nt=length(t.obs)
  
  #anchor
  k=sum(z.obs<K)
  if(k==0) k=1 #for the case when min(z.obs) < K
  zk=z.obs[k]
  
  #counting process
  Y=dN=matrix(0,n,nt)       #row is the subj, col is the time corresponding z_(i);
  for(i in 1:n){
    rank.t=which(t[i]==t.obs)
    Y[i,][1:rank.t]=1
    if(status[i]==1) dN[i,][rank.t]=1
  }
  
  #initial
  zk=z.obs[k]
  Z.BAR=Z-zk
  
  if(is.null(X)){ #no trt group
    options(warn=2)
    cox.warning<-try(cox.fit<-coxph(Surv(TIME,STATUS)~Z.BAR), silent=TRUE)
    options(warn=0) 

    if(class(cox.warning)=="try-error"){
      beta.hat=0.01
    }else{
      beta.hat=cox.fit$coefficients
      if(shape=='increasing'){;       psi= abs(beta.hat)*(z.obs-zk)
      }else if(shape=='decreasing'){; psi=-abs(beta.hat)*(z.obs-zk);  }
    }
  }else{
    options(warn=2)
    cox.warning<-try(cox.fit<-coxph(Surv(TIME,STATUS)~Z.BAR+X), silent=TRUE)
    options(warn=0)
    
    if(class(cox.warning)=="try-error"){
      beta.hat=0.01; beta=rep(0,nparm);
    }else{
      beta.hat=cox.fit$coefficients
      beta=beta.hat[-1]
      if(shape=='increasing'){;       psi= abs(beta.hat[1])*(z.obs-zk)
      }else if(shape=='decreasing'){; psi=-abs(beta.hat[1])*(z.obs-zk);  }
    }
  }
  
  #picm & beta for newton raphson algo
  if(is.null(X)){ #no trt group
    #RPA (with interval), Y, dN
    rpa.Y=RPA_ti(n, nt, m, z, z.obs, Y, dN)
    Y2=rpa.Y$Y2
    dN2=rpa.Y$dN2
    
    #picm
    dNsum=colSums(dN2)
    Delta=rowSums(dN2)
    
    dist=0; exp.beta=NA
    picm=picm.ft(psi,m,z.obs,zk,k, dN2,Y2,dNsum,Delta, eps,maxiter, shape)
    psi.new=picm$psi.new
    
    psi.full=BTF(m, n, z,z.obs, psi.new)
    if(picm$conv==0) stop
  }else{ #trt group
    iter=0;  dist=1; beta.new=rep(NA,nparm)
    while(dist>=eps){
      iter=iter+1
      if(iter>maxiter) break
      
      #RPA (with interval), Y, dN
      Yest=matrix(NA,n,nt)
      if(nparm>1){ for(j in 1:nt) Yest[,j]=Y[,j]*exp(x%*%beta)
      }else{       for(j in 1:nt) Yest[,j]=Y[,j]*exp(x*beta)
      }
      
      rpa.Y=RPA_ti(n, nt, m, z, z.obs, Yest, dN)
      Y2=rpa.Y$Y2
      dN2=rpa.Y$dN2
      
      #picm
      dNsum=colSums(dN2)
      Delta=rowSums(dN2)
      
      #estimate psi
      picm=picm.ft(psi,m,z.obs,zk,k, dN2,Y2,dNsum,Delta, eps,maxiter, shape)
      psi.new=picm$psi.new
      
      psi.full=BTF(m, n, z,z.obs, psi.new)
      if(picm$conv==0) stop
      
      #estimate beta (Y1&w1 or Y2&w2 should be the same);
      beta.new=NR.ft(x,beta,psi.full,n,nt,Y,dN, maxiter,eps)
      
      #update;
      dist=sqrt(sum((psi.new-psi)^2))+sqrt(sum((beta.new-beta)^2))
      #this can be reduced to: (1-picm$conv)+sqrt(sum((beta.new-beta)^2))
      
      psi=psi.new
      beta=beta.new
    }
    exp.beta=round(exp(beta.new), maxdec)
    exp.beta=formatC( exp.beta, format='f', digits=maxdec)
  }
  
  ###
  #picm result
  ###
  conv="converged"
  if(dist>=eps)    conv="not converged"
  if(picm$conv==0) conv="not converged"
  
  #back to full rank (later)
  psi.obs=round(psi.new, maxdec)
  
  #level sets
  psi.uniq=unique(psi.obs)
  n.lv=length(psi.uniq)
  
  lv.sets=c()
  zmin=formatC( round(min(z),maxdec), format='f', digits=maxdec)
  zmax=formatC( round(max(z),maxdec), format='f', digits=maxdec)
  
  if(n.lv==1){ #only one level sets
    lv.sets[1]=paste('[',zmin,',',zmax,']', sep='')
  }else{
    lv=list()
    if(shape=='increasing'){
      for(i in 1:n.lv)
        lv[[i]]=formatC( round(z.obs[which(psi.obs==psi.uniq[i])],maxdec)[1], format='f', digits=maxdec)
      for(i in 1:n.lv){
        if(i==n.lv){; lv.sets[i]=paste('[',lv[i][1],', ',zmax,']', sep='')
        }else{;       lv.sets[i]=paste('[',lv[i],', ',lv[i+1],')', sep='')
        }
      }
    }else{
      for(i in 1:n.lv)
        lv[i]=formatC( tail(round(z.obs[which(psi.obs==psi.uniq[i])],maxdec),1), format='f', digits=maxdec)
      for(i in 1:n.lv){
        if(i==1){;   lv.sets[i]=paste('[',zmin,', ',lv[i][1],']', sep='')
        }else{;      lv.sets[i]=paste('(',lv[i-1],', ',lv[i],']', sep='')
        };};};};
  
  psi.hat=formatC( unique(psi.obs), format='f', digits=maxdec)
  HR.hat=formatC( unique(exp(psi.obs)), format='f', digits=maxdec)
  #recover if status[1]=0 associated with z(1)
  if(n.total>n){
    n=n.total
    if(shape=='increasing'){
      psi.hat=c(-Inf,psi.hat)
      HR.hat=c(0,HR.hat)
      zmin0=formatC( round(min(Z),maxdec), format='f', digits=maxdec)
      lv.sets=c(paste('[',zmin0,', ',zmin,')', sep=''), lv.sets)
    }else{
      psi.hat=c(psi.hat,-Inf)
      HR.hat=c(HR.hat,0)
      zmax0=formatC( round(max(Z),maxdec), format='f', digits=maxdec)
      lv.sets=c(lv.sets, paste('(',zmax,', ',zmax0,']', sep=''))
    }
  }
  
  est=data.frame(psi.hat=psi.hat, HR.hat=HR.hat, lv.set=lv.sets)
  names(est)=c("psi.hat","exp(psi.hat)","level set of psi.hat")
  
  z.range=range(Z)
  
  #temprary to compute psi.full
  psi.full2=rep(NA,n.total)
  z.full=sort(Z)
  for(i in 1:n){
    psi.full2[which(z.full%in%z[i])]=psi.full[i]
  }
  psi.full2[is.na(psi.full2)]=-Inf
  psi.full=psi.full2
  
  return(list(est=est, exp.beta=exp.beta, conv=conv,
              psi=psi.obs, z=z.obs, z.range=z.range, K=K, shape=shape, n=n, nevent=sum(STATUS), njump=m,
              psi.full=psi.full, z.full=z.full))
}

Try the isoph package in your browser

Any scripts or data that you put into this service are public.

isoph documentation built on May 1, 2019, 10:15 p.m.