simulations/pl_methods.R

   1 library(stats4)
   2 library(bbmle)
   3 library(matrixStats)
   4
   5 zhang.mle.dv <- function(df){
   6     df.obs <- df[!is.na(y.obs)]
   7     df.unobs <- df[is.na(y.obs)]
   8
   9     fp <- df.obs[(w_pred==1) & (y.obs != w_pred),.N]
  10     tn <- df.obs[(w_pred == 0) & (y.obs == w_pred),.N]
  11     fpr <- fp / (fp+tn)
  12
  13     fn <- df.obs[(w_pred==0) & (y.obs != w_pred), .N]
  14     tp <- df.obs[(w_pred==1) & (y.obs == w_pred),.N]
  15     fnr <- fn / (fn+tp)
  16
  17     nll <- function(B0=0, Bxy=0, Bzy=0){
  18
  19
  20         ## observed case
  21         ll.y.obs <- vector(mode='numeric', length=nrow(df.obs))
  22         ll.y.obs[df.obs$y.obs==1] <- with(df.obs[y.obs==1], plogis(B0 + Bxy * x + Bzy * z,log=T))
  23         ll.y.obs[df.obs$y.obs==0] <- with(df.obs[y.obs==0], plogis(B0 + Bxy * x + Bzy * z,log=T,lower.tail=FALSE))
  24
  25         ll <- sum(ll.y.obs)
  26
  27         pi.y.1 <- with(df.unobs,plogis(B0 + Bxy * x + Bzy*z, log=T))
  28         #pi.y.0 <- with(df.unobs,plogis(B0 + Bxy * x + Bzy*z, log=T,lower.tail=FALSE))
  29
  30         lls <- with(df.unobs, colLogSumExps(rbind(w_pred * colLogSumExps(rbind(log(fpr), log(1 - fnr - fpr)+pi.y.1)),
  31         (1-w_pred) * (log(1-fpr) - exp(log(1-fnr-fpr)+pi.y.1)))))
  32
  33         ll <- ll + sum(lls)
  34 #        print(paste0(B0,Bxy,Bzy))
  35 #        print(ll)
  36         return(-ll)
  37     }
  38     mlefit <- mle2(minuslogl = nll, control=list(maxit=1e6),method='L-BFGS-B',lower=c(B0=-Inf, Bxy=-Inf, Bzy=-Inf),
  39                    upper=c(B0=Inf, Bxy=Inf, Bzy=Inf))
  40     return(mlefit)
  41 }
  42
  43
  44 ## model from Zhang's arxiv paper, with predictions for y
  45 ## Zhang got this model from Hausman 1998
  46 zhang.mle.iv <- function(df){
  47     df.obs <- df[!is.na(x.obs)]
  48     df.unobs <- df[is.na(x.obs)]
  49
  50     tn <- df.obs[(w_pred == 0) & (x.obs == w_pred),.N]
  51     fn <- df.obs[(w_pred==0) & (x.obs==1), .N]
  52     npv <- tn / (tn + fn)
  53
  54
  55     tp <- df.obs[(w_pred==1) & (x.obs == w_pred),.N]
  56     fp <- df.obs[(w_pred==1) & (x.obs == 0),.N]
  57     ppv <- tp / (tp + fp)
  58
  59     nll <- function(B0=0, Bxy=0, Bzy=0, sigma_y=9){
  60
  61     ## fpr = 1 - TNR
  62     ### Problem: accounting for uncertainty in ppv / npv
  63
  64     ## fnr = 1 - TPR
  65     ll.y.obs <- with(df.obs, dnorm(y, B0 + Bxy * x + Bzy * z, sd=sigma_y,log=T))
  66
  67     ll <- sum(ll.y.obs)
  68     # unobserved case; integrate out x
  69     ll.x.1 <- with(df.unobs, dnorm(y, B0 + Bxy + Bzy * z, sd = sigma_y, log=T))
  70     ll.x.0 <- with(df.unobs, dnorm(y, B0 + Bzy * z, sd = sigma_y,log=T))
  71
  72     ## case x == 1
  73     lls.x.1 <- colLogSumExps(rbind(log(ppv) + ll.x.1, log(1-ppv) + ll.x.0))
  74
  75     ## case x == 0
  76     lls.x.0 <- colLogSumExps(rbind(log(1-npv) + ll.x.1, log(npv) + ll.x.0))
  77
  78     lls <- colLogSumExps(rbind(df.unobs$w_pred * lls.x.1, (1-df.unobs$w_pred) * lls.x.0))
  79
  80     ll <- ll + sum(lls)
  81
  82     }
  83     mlefit <- mle2(minuslogl = nll, control=list(maxit=1e6), lower=list(sigma_y=0.00001, B0=-Inf, Bxy=-Inf, Bzy=-Inf),
  84                    upper=list(sigma_y=Inf, B0=Inf, Bxy=Inf, Bzy=Inf),method='L-BFGS-B')
  85     return(mlefit)
  86 }