R/funs.R

In wx202/OptSurrogateSurv: PTE (proportion of treatment explained) estimate for optimally-transformed surrogate with censored primary outcome

#' Title
#'
#' @param zz
#' @param zi
#' @param bw
#'
#' @return
#' @export
#'
#' @examples
Kern.FUN <- function(zz,zi,bw)
{
  out = (VTM(zz,length(zi))- zi)/bw
  dnorm(out)/bw

}


#' Title
#'
#' @param vc
#' @param dm
#'
#' @return
#' @export
#'
#' @examples
VTM<-function(vc, dm){
  matrix(vc, ncol=length(vc), nrow=dm, byrow=T)
}


#' Title
#'
#' @param xob
#' @param deltaob
#' @param aob
#' @param n
#' @param t
#' @param t.0
#'
#' @return
#' @export
#'
#' @examples
WEIGHT<-function(xob,deltaob,aob,n,t,t.0){
  x=xob[aob==0]; delta=1-deltaob[aob==0]
  xsort=sort(x[delta==1])
  risk=VTM(x,length(xsort))>=xsort
  risk.n=apply(risk,1,sum)
  Lam=cumsum(1/risk.n)
  s0=exp(-Lam)
  sur0=data.frame("time"=xsort,"surv"=s0)
  x=xob[aob==1]; delta=1-deltaob[aob==1]
  xsort=sort(x[delta==1])
  risk=VTM(x,length(xsort))>=xsort
  risk.n=apply(risk,1,sum)
  Lam=cumsum(1/risk.n)
  s1=exp(-Lam)
  sur1=data.frame("time"=xsort,"surv"=s1)
  ind0=c(sapply(1:n, function(kk){which.min(abs(xob[kk]-sur0$time))}))
  G0=sur0$surv[ind0]
  ind1=c(sapply(1:n, function(kk){which.min(abs(xob[kk]-sur1$time))}))
  G1=sur1$surv[ind1]
  G=(1-aob)*G0+aob*G1
  G.t0=(1-aob)*G0[which.min(abs(t.0-xob))]+aob*G1[which.min(abs(t.0-xob))]
  G.t=(1-aob)*G0[which.min(abs(t-xob))]+aob*G1[which.min(abs(t-xob))]
  wei.t0=(xob<=t.0)*deltaob/G+(xob>t.0)/G.t0; #wei.t0[is.nan(wei.t0)]=max(wei.t0[!is.nan(wei.t0)])
  wei.t=(xob<=t)*deltaob/G+(xob>t)/G.t; #wei.t[is.nan(wei.t)]=max(wei.t[!is.nan(wei.t)])

  out=cbind(wei.t0,wei.t,G.t,G.t0,G1,G0,G)#,pte2
}


#' Title
#'
#' @param xob
#' @param deltaob
#' @param aob
#' @param n
#' @param v
#' @param t
#' @param t.0
#'
#' @return
#' @export
#'
#' @examples
WEIGHT.p<-function(xob,deltaob,aob,n,v,t,t.0){
  x=xob[aob==0];delta=1-deltaob[aob==0]
  xsort=sort(x[delta==1])
  risk=VTM(x,length(xsort))>=xsort
  risk.n=apply(risk*VTM(v[aob==0],length(xsort)),1,sum)
  N=(VTM(x,length(xsort))<=xsort)*VTM(delta,length(xsort))
  dN=rbind( N[1,],N[-1,]-N[-length(xsort),] )
  nu=apply(VTM(v[aob==0],length(xsort))*dN, 1, sum)
  Lam=cumsum(nu/risk.n)
  s0=exp(-Lam)
  sur0=data.frame("time"=xsort,"surv"=s0)
  x=xob[aob==1];delta=1-deltaob[aob==1]
  xsort=sort(x[delta==1])
  risk=VTM(x,length(xsort))>=xsort
  risk.n=apply(risk*VTM(v[aob==1],length(xsort)),1,sum)
  N=(VTM(x,length(xsort))<=xsort)*VTM(delta,length(xsort))
  dN=rbind( N[1,],N[-1,]-N[-length(xsort),] )
  nu=apply(dN*VTM(v[aob==1],length(xsort)), 1, sum)
  Lam=cumsum(nu/risk.n)
  s1=exp(-Lam)
  sur1=data.frame("time"=xsort,"surv"=s1)
  ind0=c(sapply(1:n, function(kk){which.min(abs(xob[kk]-sur0$time))}))
  G0=sur0$surv[ind0]
  ind1=c(sapply(1:n, function(kk){which.min(abs(xob[kk]-sur1$time))}))
  G1=sur1$surv[ind1]
  G=(1-aob)*G0+aob*G1
  G.t0=(1-aob)*G0[which.min(abs(t.0-xob))]+aob*G1[which.min(abs(t.0-xob))]
  G.t=(1-aob)*G0[which.min(abs(t-xob))]+aob*G1[which.min(abs(t-xob))]
  wei.t0=(xob<=t.0)*deltaob/G+(xob>t.0)/G.t0; #wei.t0[is.nan(wei.t0)]=max(wei.t0[!is.nan(wei.t0)])
  wei.t=(xob<=t)*deltaob/G+(xob>t)/G.t; #wei.t[is.nan(wei.t)]=max(wei.t[!is.nan(wei.t)])

  out=cbind(wei.t0,wei.t,G.t,G.t0,G1,G0,G)#,pte2
}


#' Title
#'
#' @param vv
#' @param t
#' @param t.0
#' @param data
#' @param data1
#' @param data2
#' @param s
#' @param nn
#'
#' @return
#' @export
#'
#' @examples
resam<- function(vv,t,t.0,data,data1,data2,s,nn,step){
  n1=nrow(data1); n2=nrow(data2)

  ################ pte2
  xob=data1[,1];yob=data1[,2];deltaob=data1[,3];s.ob=data1[,4];aob=data1[,5]; n=n1; v=vv[1:nrow(data1)]

  temp=WEIGHT.p(xob,deltaob,aob,n=n,v,t,t.0)
  wei.t0=temp[,1];wei.t=temp[,2]
  bw = 1.06*sd(s.ob)*n^(-1/5)/(n^0.11)
  kern = Kern.FUN(zz=s,zi=s.ob,bw)
  p0.t0.s.hat=apply(as.numeric(v)*(aob==0)*(xob>t.0)*kern*wei.t0,2,sum)/sum(as.numeric(v)*(aob==0)*wei.t0)
  p.t0.s.hat=apply(as.numeric(v)*(xob>t.0)*kern*wei.t0,2,mean)/mean(as.numeric(v)*wei.t0)
  m.s.hat=(apply(as.numeric(v)*(xob>t)*kern*wei.t,2,sum)/sum(as.numeric(v)*wei.t))/
    (apply(as.numeric(v)*(xob>t.0)*kern*wei.t0,2,sum)/sum(as.numeric(v)*wei.t0))
  kern2 = Kern.FUN(zz=s.ob,zi=s.ob,bw)
  m.sob.hat=(apply(as.numeric(v)*(xob>t)*kern2*wei.t,2,sum)/sum(as.numeric(v)*wei.t))/
    (apply(as.numeric(v)*(xob>t.0)*kern2*wei.t0,2,sum)/sum(as.numeric(v)*wei.t0))
  c.hat=sum(as.numeric(v)*(aob==0)*(xob>t.0)*(xob>t)*wei.t)/sum(as.numeric(v)*(aob==0)*wei.t)-
    sum(as.numeric(v)*(aob==0)*m.sob.hat*(xob>t.0)*wei.t0)/sum(as.numeric(v)*(aob==0)*wei.t0)
  integrand<-p0.t0.s.hat^2/p.t0.s.hat
  temp=(integrand[1] + integrand[nn+1] + 2*sum(integrand[seq(2,nn,by=2)]) + 4 *sum(integrand[seq(3,nn-1, by=2)]) )*step/3
  p0t=sum(as.numeric(v)*(aob==0)*(xob<=t.0)*wei.t0)/sum(as.numeric(v)*(aob==0)*wei.t0)
  pt=mean(as.numeric(v)*(xob<=t.0)*wei.t0)/mean(as.numeric(v)*wei.t0)
  g.s.hat=m.s.hat+p0.t0.s.hat/p.t0.s.hat *c.hat/(temp+p0t^2/pt)
  g1.hat=p0t/pt *c.hat/(temp+p0t^2/pt)
  g.s.hat.1=g.s.hat
  g1.hat.1=g1.hat

  ####
  xob=data2[,1];yob=data2[,2];deltaob=data2[,3];s.ob=data2[,4];aob=data2[,5];n=n2;v=vv[(nrow(data1)+1):nrow(data)]

  temp=WEIGHT.p(xob,deltaob,aob,n=n,v,t,t.0)
  wei.t0=temp[,1];wei.t=temp[,2]
  causal=sum(v*(xob>t)*aob*wei.t)/sum(v*aob*wei.t)-sum(v*(xob>t)*(1-aob)*wei.t)/sum(v*(1-aob)*wei.t)
  tempind=c(sapply(1:n, function(kk){which.min(abs(s.ob[kk]-s))}))
  causals=sum(((xob<=t.0)*g1.hat+(xob>t.0)*g.s.hat[tempind])*v*aob*wei.t0)/sum(v*aob*wei.t0)-
    sum(((xob<=t.0)*g1.hat+(xob>t.0)*g.s.hat[tempind])*(1-aob)*v*wei.t0)/sum(v*(1-aob)*wei.t0)
  pte2=causals/causal

  ################ pte1
  xob=data2[,1];yob=data2[,2];deltaob=data2[,3];s.ob=data2[,4];aob=data2[,5];n=n2; v=vv[(nrow(data1)+1):nrow(data)]

  temp=WEIGHT.p(xob,deltaob,aob,n=n,v,t,t.0)
  wei.t0=temp[,1];wei.t=temp[,2]
  bw = 1.06*sd(s.ob)*n^(-1/5)/(n^0.11)
  kern = Kern.FUN(zz=s,zi=s.ob,bw)
  p0.t0.s.hat=apply(as.numeric(v)*(aob==0)*(xob>t.0)*kern*wei.t0,2,sum)/sum(as.numeric(v)*(aob==0)*wei.t0)
  p.t0.s.hat=apply(as.numeric(v)*(xob>t.0)*kern*wei.t0,2,mean)/mean(as.numeric(v)*wei.t0)
  m.s.hat=(apply(as.numeric(v)*(xob>t)*kern*wei.t,2,sum)/sum(as.numeric(v)*wei.t))/
    (apply(as.numeric(v)*(xob>t.0)*kern*wei.t0,2,sum)/sum(as.numeric(v)*wei.t0))
  kern2 = Kern.FUN(zz=s.ob,zi=s.ob,bw)
  m.sob.hat=(apply(as.numeric(v)*(xob>t)*kern2*wei.t,2,sum)/sum(as.numeric(v)*wei.t))/
    (apply(as.numeric(v)*(xob>t.0)*kern2*wei.t0,2,sum)/sum(as.numeric(v)*wei.t0))
  c.hat=sum(as.numeric(v)*(aob==0)*(xob>t.0)*(xob>t)*wei.t)/sum(as.numeric(v)*(aob==0)*wei.t)-
    sum(as.numeric(v)*(aob==0)*m.sob.hat*(xob>t.0)*wei.t0)/sum(as.numeric(v)*(aob==0)*wei.t0)
  integrand<-p0.t0.s.hat^2/p.t0.s.hat
  temp=(integrand[1] + integrand[nn+1] + 2*sum(integrand[seq(2,nn,by=2)]) + 4 *sum(integrand[seq(3,nn-1, by=2)]) )*step/3
  p0t=sum(as.numeric(v)*(aob==0)*(xob<=t.0)*wei.t0)/sum(as.numeric(v)*(aob==0)*wei.t0)
  pt=mean(as.numeric(v)*(xob<=t.0)*wei.t0)/mean(as.numeric(v)*wei.t0)
  g.s.hat=m.s.hat+p0.t0.s.hat/p.t0.s.hat *c.hat/(temp+p0t^2/pt)
  g1.hat=p0t/pt *c.hat/(temp+p0t^2/pt)
  g.s.hat.2=g.s.hat
  g1.hat.2=g1.hat

  ####
  xob=data1[,1];yob=data1[,2];deltaob=data1[,3];s.ob=data1[,4];aob=data1[,5];n=n1;v=vv[1:nrow(data1)]

  temp=WEIGHT.p(xob,deltaob,aob,n=n,v,t,t.0)
  wei.t0=temp[,1];wei.t=temp[,2]
  causal=sum(v*(xob>t)*aob*wei.t)/sum(v*aob*wei.t)-sum(v*(xob>t)*(1-aob)*wei.t)/sum(v*(1-aob)*wei.t)
  tempind=c(sapply(1:n, function(kk){which.min(abs(s.ob[kk]-s))}))
  causals=sum(((xob<=t.0)*g1.hat+(xob>t.0)*g.s.hat[tempind])*v*aob*wei.t0)/sum(v*aob*wei.t0)-
    sum(((xob<=t.0)*g1.hat+(xob>t.0)*g.s.hat[tempind])*(1-aob)*v*wei.t0)/sum(v*(1-aob)*wei.t0)
  pte1=causals/causal

  pte=(pte1+pte2)/2
  g1.es=(g1.hat.1+g1.hat.2)/2
  gs.es=(g.s.hat.1+g.s.hat.2)/2

  out=c(pte2,pte1,pte,g1.es,gs.es)

}