update simulation and mle code

[ml_measurement_error_public.git] / simulations / simulation_base.R

diff --git a/simulations/simulation_base.R b/simulations/simulation_base.R
index 0f03276c432f257ede87dcc8c0d78fca6834557b..ee46ec6e6d303462ff71c9b62c132e82752a76fb 100644 (file)
--- a/simulations/simulation_base.R
+++ b/simulations/simulation_base.R
@@ -41,21 +41,26 @@ my.pseudo.mle <- function(df){
 ## Zhang got this model from Hausman 1998
 ### I think this is actually eqivalent to the pseudo.mle method
 zhang.mle.iv <- function(df){
 ## Zhang got this model from Hausman 1998
 ### I think this is actually eqivalent to the pseudo.mle method
 zhang.mle.iv <- function(df){

-    nll <- function(B0=0, Bxy=0, Bzy=0, sigma_y=0.1, ppv=0.9, npv=0.9){

     df.obs <- df[!is.na(x.obs)]
     df.unobs <- df[is.na(x.obs)]
 
     df.obs <- df[!is.na(x.obs)]
     df.unobs <- df[is.na(x.obs)]
 

+    tn <- df.obs[(w_pred == 0) & (x.obs == w_pred),.N]
+    pn <- df.obs[(w_pred==0), .N]
+    npv <- tn / pn
+
+    tp <- df.obs[(w_pred==1) & (x.obs == w_pred),.N]
+    pp <- df.obs[(w_pred==1),.N]
+    ppv <- tp / pp
+
+    nll <- function(B0=0, Bxy=0, Bzy=0, sigma_y=0.1){
+

     ## fpr = 1 - TNR
     ### Problem: accounting for uncertainty in ppv / npv
     ## fpr = 1 - TNR
     ### Problem: accounting for uncertainty in ppv / npv

-    
-    ll.w1x1.obs <- with(df.obs[(w_pred==1)], dbinom(x.obs,size=1,prob=ppv,log=T))
-    ll.w0x0.obs <- with(df.obs[(w_pred==0)], dbinom(1-x.obs,size=1,prob=npv,log=T))

 
     ## fnr = 1 - TPR
     ll.y.obs <- with(df.obs, dnorm(y, B0 + Bxy * x + Bzy * z, sd=sigma_y,log=T))
     ll <- sum(ll.y.obs)
 
     ## fnr = 1 - TPR
     ll.y.obs <- with(df.obs, dnorm(y, B0 + Bxy * x + Bzy * z, sd=sigma_y,log=T))
     ll <- sum(ll.y.obs)

-    ll <- ll + sum(ll.w1x1.obs) + sum(ll.w0x0.obs)
-
+    

     # unobserved case; integrate out x
     ll.x.1 <- with(df.unobs, dnorm(y, B0 + Bxy + Bzy * z, sd = sigma_y, log=T))
     ll.x.0 <- with(df.unobs, dnorm(y, B0 + Bzy * z, sd = sigma_y,log=T))
     # unobserved case; integrate out x
     ll.x.1 <- with(df.unobs, dnorm(y, B0 + Bxy + Bzy * z, sd = sigma_y, log=T))
     ll.x.0 <- with(df.unobs, dnorm(y, B0 + Bzy * z, sd = sigma_y,log=T))


@@ -66,55 +71,90 @@ zhang.mle.iv <- function(df){
     ## case x == 0
     lls.x.0 <- colLogSumExps(rbind(log(1-npv) + ll.x.1, log(npv) + ll.x.0))
 
     ## case x == 0
     lls.x.0 <- colLogSumExps(rbind(log(1-npv) + ll.x.1, log(npv) + ll.x.0))
 

-    lls <- colLogSumExps(rbind(lls.x.1, lls.x.0))
+    lls <- colLogSumExps(rbind(df.unobs$w_pred * lls.x.1, (1-df.unobs$w_pred) * lls.x.0))

     ll <- ll + sum(lls)
     return(-ll)
     }    
     ll <- ll + sum(lls)
     return(-ll)
     }    

-    mlefit <- mle2(minuslogl = nll, control=list(maxit=1e6), lower=list(sigma_y=0.0001, B0=-Inf, Bxy=-Inf, Bzy=-Inf,ppv=0.00001, npv=0.00001),
-                   upper=list(sigma_y=Inf, B0=Inf, Bxy=Inf, Bzy=Inf, ppv=0.99999,npv=0.99999),method='L-BFGS-B')
+    mlefit <- mle2(minuslogl = nll, control=list(maxit=1e6), lower=list(sigma_y=0.0001, B0=-Inf, Bxy=-Inf, Bzy=-Inf),
+                   upper=list(sigma_y=Inf, B0=Inf, Bxy=Inf, Bzy=Inf),method='L-BFGS-B')

     return(mlefit)
 }
 
     return(mlefit)
 }
 

-## this is equivalent to the pseudo-liklihood model from Carolla
-zhang.mle.dv <- function(df){
+## this is equivalent to the pseudo-liklihood model from Caroll
+## zhang.mle.dv <- function(df){

 
 

-    nll <- function(B0=0, Bxy=0, Bzy=0, ppv=0.9, npv=0.9){
-    df.obs <- df[!is.na(y.obs)]
+##     nll <- function(B0=0, Bxy=0, Bzy=0, ppv=0.9, npv=0.9){
+##     df.obs <- df[!is.na(y.obs)]

 
 

-    ## fpr = 1 - TNR
-    ll.w0y0 <- with(df.obs[y.obs==0],dbinom(1-w_pred,1,npv,log=TRUE))
-    ll.w1y1 <- with(df.obs[y.obs==1],dbinom(w_pred,1,ppv,log=TRUE))
-
-    # observed case
-    ll.y.obs <- vector(mode='numeric', length=nrow(df.obs))
-    ll.y.obs[df.obs$y.obs==1] <- with(df.obs[y.obs==1], plogis(B0 + Bxy * x + Bzy * z,log=T))
-    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))
-
-    ll <- sum(ll.y.obs) + sum(ll.w0y0) + sum(ll.w1y1)
-
-    # unobserved case; integrate out y
-    ## case y = 1
-    ll.y.1 <- vector(mode='numeric', length=nrow(df))
-    pi.y.1 <- with(df,plogis(B0 + Bxy * x + Bzy*z, log=T))
-    ## P(w=1| y=1)P(y=1) + P(w=0|y=1)P(y=1) = P(w=1,y=1) + P(w=0,y=1)
-    lls.y.1 <- colLogSumExps(rbind(log(ppv) + pi.y.1, log(1-ppv) + pi.y.1))
+##     ## fpr = 1 - TNR
+##     ll.w0y0 <- with(df.obs[y.obs==0],dbinom(1-w_pred,1,npv,log=TRUE))
+##     ll.w1y1 <- with(df.obs[y.obs==1],dbinom(w_pred,1,ppv,log=TRUE))
+
+##     # observed case
+##     ll.y.obs <- vector(mode='numeric', length=nrow(df.obs))
+##     ll.y.obs[df.obs$y.obs==1] <- with(df.obs[y.obs==1], plogis(B0 + Bxy * x + Bzy * z,log=T))
+##     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))
+
+##     ll <- sum(ll.y.obs) + sum(ll.w0y0) + sum(ll.w1y1)
+
+##     # unobserved case; integrate out y
+##     ## case y = 1
+##     ll.y.1 <- vector(mode='numeric', length=nrow(df))
+##     pi.y.1 <- with(df,plogis(B0 + Bxy * x + Bzy*z, log=T))
+##     ## P(w=1| y=1)P(y=1) + P(w=0|y=1)P(y=1) = P(w=1,y=1) + P(w=0,y=1)
+##     lls.y.1 <- colLogSumExps(rbind(log(ppv) + pi.y.1, log(1-ppv) + pi.y.1))

     
     

-    ## case y = 0
-    ll.y.0 <- vector(mode='numeric', length=nrow(df))
-    pi.y.0 <- with(df,plogis(B0 + Bxy * x + Bzy*z, log=T,lower.tail=FALSE))
+##     ## case y = 0
+##     ll.y.0 <- vector(mode='numeric', length=nrow(df))
+##     pi.y.0 <- with(df,plogis(B0 + Bxy * x + Bzy*z, log=T,lower.tail=FALSE))
+
+##     ## P(w=1 | y=0)P(y=0) + P(w=0|y=0)P(y=0) = P(w=1,y=0) + P(w=0,y=0)
+##     lls.y.0 <- colLogSumExps(rbind(log(npv) + pi.y.0, log(1-npv) + pi.y.0))
+
+##     lls <- colLogSumExps(rbind(lls.y.1, lls.y.0))
+##     ll <- ll + sum(lls)
+##     return(-ll)
+##     }    
+##     mlefit <- mle2(minuslogl = nll, control=list(maxit=1e6),method='L-BFGS-B',lower=list(B0=-Inf, Bxy=-Inf, Bzy=-Inf, ppv=0.001,npv=0.001),
+##                    upper=list(B0=Inf, Bxy=Inf, Bzy=Inf,ppv=0.999,npv=0.999))
+##     return(mlefit)
+## }

 
 

-    ## P(w=1 | y=0)P(y=0) + P(w=0|y=0)P(y=0) = P(w=1,y=0) + P(w=0,y=0)
-    lls.y.0 <- colLogSumExps(rbind(log(npv) + pi.y.0, log(1-npv) + pi.y.0))
+zhang.mle.dv <- function(df){
+    df.obs <- df[!is.na(y.obs)]
+    df.unobs <- df[is.na(y.obs)]

 
 

-    lls <- colLogSumExps(rbind(lls.y.1, lls.y.0))
-    ll <- ll + sum(lls)
-    return(-ll)
+    fp <- df.obs[(w_pred==1) & (y.obs != w_pred),.N]
+    p <- df.obs[(w_pred==1),.N]
+    fpr <- fp / p
+    fn <- df.obs[(w_pred==0) & (y.obs != w_pred), .N]
+    n <- df.obs[(w_pred==0),.N]
+    fnr <- fn / n
+
+    nll <- function(B0=0, Bxy=0, Bzy=0){
+
+
+        ## observed case
+        ll.y.obs <- vector(mode='numeric', length=nrow(df.obs))
+        ll.y.obs[df.obs$y.obs==1] <- with(df.obs[y.obs==1], plogis(B0 + Bxy * x + Bzy * z,log=T))
+        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))
+
+        ll <- sum(ll.y.obs)
+
+        pi.y.1 <- with(df,plogis(B0 + Bxy * x + Bzy*z, log=T))
+        pi.y.0 <- with(df,plogis(B0 + Bxy * x + Bzy*z, log=T,lower.tail=FALSE))
+
+        lls <- with(df.unobs, colLogSumExps(rbind(w_pred * colLogSumExps(rbind(log(fpr), log(1 - fnr - fpr)+pi.y.1)),
+        (1-w_pred) * colLogSumExps(rbind(log(1-fpr), log(1 - fnr - fpr)+pi.y.0)))))
+    
+        ll <- ll + sum(lls)
+        return(-ll)

     }    
     }    

-    mlefit <- mle2(minuslogl = nll, control=list(maxit=1e6),method='L-BFGS-B',lower=list(B0=-Inf, Bxy=-Inf, Bzy=-Inf, ppv=0.001,npv=0.001),
-                   upper=list(B0=Inf, Bxy=Inf, Bzy=Inf,ppv=0.999,npv=0.999))
+    mlefit <- mle2(minuslogl = nll, control=list(maxit=1e6),method='L-BFGS-B',lower=c(B0=-Inf, Bxy=-Inf, Bzy=-Inf),
+                   upper=c(B0=Inf, Bxy=Inf, Bzy=Inf))

     return(mlefit)
 }
     return(mlefit)
 }

-
+ 

 ## This uses the likelihood approach from Carroll page 353.
 ## assumes that we have a good measurement error model
 my.mle <- function(df){
 ## This uses the likelihood approach from Carroll page 353.
 ## assumes that we have a good measurement error model
 my.mle <- function(df){


@@ -211,7 +251,7 @@ run_simulation_depvar <- function(df, result, outcome_formula=y~x+z, proxy_formu
     naivecont.ci.Bxy <- confint(model.naive.cont)['x',]
     naivecont.ci.Bzy <- confint(model.naive.cont)['z',]
 
     naivecont.ci.Bxy <- confint(model.naive.cont)['x',]
     naivecont.ci.Bzy <- confint(model.naive.cont)['z',]
 

-    ## my implementatoin of liklihood based correction
+    ## my implementation of liklihood based correction

 
     temp.df <- copy(df)
     temp.df[,y:=y.obs]
 
     temp.df <- copy(df)
     temp.df[,y:=y.obs]


@@ -241,7 +281,8 @@ run_simulation_depvar <- function(df, result, outcome_formula=y~x+z, proxy_formu
                           Bzy.est.zhang = coef['Bzy'],
                           Bzy.ci.upper.zhang = ci['Bzy','97.5 %'],
                           Bzy.ci.lower.zhang = ci['Bzy','2.5 %']))
                           Bzy.est.zhang = coef['Bzy'],
                           Bzy.ci.upper.zhang = ci['Bzy','97.5 %'],
                           Bzy.ci.lower.zhang = ci['Bzy','2.5 %']))

-                          
+
+    

 
     # amelia says use normal distribution for binary variables.
     tryCatch({
 
     # amelia says use normal distribution for binary variables.
     tryCatch({


@@ -278,7 +319,7 @@ run_simulation_depvar <- function(df, result, outcome_formula=y~x+z, proxy_formu
 
 
 ## outcome_formula, proxy_formula, and truth_formula are passed to measerr_mle 
 
 
 ## outcome_formula, proxy_formula, and truth_formula are passed to measerr_mle 

-run_simulation <-  function(df, result, outcome_formula=y~x+z, proxy_formula=w_pred~x, truth_formula=x~z){
+run_simulation <-  function(df, result, outcome_formula=y~x+z, proxy_formula=NULL, truth_formula=NULL){

 
     accuracy <- df[,mean(w_pred==x)]
     result <- append(result, list(accuracy=accuracy))
 
     accuracy <- df[,mean(w_pred==x)]
     result <- append(result, list(accuracy=accuracy))


@@ -320,7 +361,7 @@ run_simulation <-  function(df, result, outcome_formula=y~x+z, proxy_formula=w_p
                                   
 
     tryCatch({
                                   
 
     tryCatch({

-    amelia.out.k <- amelia(df, m=200, p2s=0, idvars=c('x','w_pred'))
+    amelia.out.k <- amelia(df, m=200, p2s=0, idvars=c('x','w'))

     mod.amelia.k <- zelig(y~x.obs+z, model='ls', data=amelia.out.k$imputations, cite=FALSE)
     (coefse <- combine_coef_se(mod.amelia.k, messages=FALSE))
 
     mod.amelia.k <- zelig(y~x.obs+z, model='ls', data=amelia.out.k$imputations, cite=FALSE)
     (coefse <- combine_coef_se(mod.amelia.k, messages=FALSE))