R/03a.compLTP_binary.R

In btrm: Bayesian Treed Regression Model for Personalized Prediction and Precision Diagnostics

#postb=function(b){
#  z.j.mat=cbind(1,z.j)
#  bhat.mat=matrix(c(b[1],b[2]),nrow=2,ncol=1)
#  zb.j=c(z.j.mat%*%bhat.mat)
#  loglik=sum(log(exp(zb.j)/(1+exp(zb.j)))[which(y.j==1)])+sum(log(1/(1+exp(zb.j)))[which(y.j==0)])

#  a=matrix(c(0,0),nrow=2,ncol=1)
#  A=matrix(c(100,0,0,100),2,2)
#  logprior=-log(2*pi)-(1/2)*log(det(A))-(1/2)*t(bhat.mat-a) %*% ginv(A) %*% (bhat.mat-a)

#  logpost=loglik+logprior
#  return(logpost)
#}
#optim(c(1,1),postb,control=list(fnscale=-1),hessian=FALSE)

#LTP: likelihood, tree probabiity, posterior
compLTP_binary=function(ET,base,power){

  ###################################################
  ###1. log likelihood
  ###################################################
  loglik=0
  bhat=list()
  yhat=NA #yhat
  for(j in ET$terminal){
    wh.j=which(ET$node.hat==j) #idx for subgroup j
    y.j=ET$y[wh.j] #y for subgroup j
    z.j=ET$z[wh.j,ET$marker[j]] #selected marker for subgroup j
    depth.j=floor(log2(j))+1

    #bayesian logistic wt prior: biv indep normal
    ssglm.j=arm::bayesglm(y.j~z.j,family=binomial,
                          prior.mean = 0, prior.scale = 100, prior.df = Inf, #normal prior for beta1
                          prior.mean.for.intercept = 0, prior.scale.for.intercept = 100, prior.df.for.intercept = Inf, #normal prior for beta0
                          scaled=FALSE)


    #compute predictive posterior (Laplace approximation)
    #(a) loglik
    z.j.mat=cbind(1,z.j)
    bhat.mat=matrix(ssglm.j$coefficients,nrow=2,ncol=1)
    zb.j=c(z.j.mat%*%bhat.mat)
    loglik.j=sum(log(exp(zb.j)/(1+exp(zb.j)))[which(y.j==1)])+sum(log(1/(1+exp(zb.j)))[which(y.j==0)])

    #(b) logprior
    a=matrix(c(0,0),nrow=2,ncol=1) #prior
    A=matrix(c(100,0,0,100),2,2)
    A.inv=ginv(A)
    logprior.j=-log(2*pi)-(1/2)*log(det(A))-(1/2)*t(bhat.mat-a)%*%ginv(A)%*%(bhat.mat-a)

    #(c) marginal loglik
    p.j=exp(zb.j)/(1+exp(zb.j))
    p.jj=p.j*(1-p.j)
    n.j=length(p.jj)
    H=A.inv
    for(i in 1:n.j)
      H=H+p.jj[i]*z.j.mat[i,]%*%t(z.j.mat[i,])
    loglik=loglik+loglik.j+logprior.j+log(2*pi*(det(H)^(1/2))) #Laplace approximation

    bhat[[j]]=ssglm.j$coefficients
    yhat[wh.j]=ssglm.j$fitted.values
  }
  ET$bhat=bhat
  ET$loglik=loglik
  ET$yhat=yhat

  ###################################################
  ##2 log tree prob
  ###################################################
  logTreeProb=0
  n.terminal=ET$terminal
  for(m in ET$terminal){ #each terminal node
    depth=floor(log2(m))
    logTreeProb = logTreeProb+log(1-prior.tree(base,power,depth))+log(ET$dir.marker[ET$marker[m]])
                              #1-splitting prob                  # marker sel prob
  }

  if(ET$numNodes>1){ #each internal node
    for(m in ET$internal){
      depth=floor(log2(m))
      logTreeProb=logTreeProb+log(prior.tree(base,power,depth))+log(ET$dir.predictor[ET$splitVariable[m]])+log(ET$eta[m])
                              #splitting prob                   # feature sel prob                         #cutoff sel prob
      #eta[m] is the number of unique values for mth predction selected for the mth internal node
    }
  }
  ET$logTreeProb=logTreeProb

  ###################################################
  ###3. log posterior
  ###################################################
  ET$logPosterior=loglik+logTreeProb #proportion of posterior since noming constant is ignored.

  return(ET)
}