update real data examples code and rerun project.

[ml_measurement_error_public.git] / simulations / simulation_base.R

diff --git a/simulations/simulation_base.R b/simulations/simulation_base.R
index 82b17a737ae05c9a98109da6164a9c75a10aebc3..e715edfaf61b9b1f423001808187ad66828a43b0 100644 (file)
--- a/simulations/simulation_base.R
+++ b/simulations/simulation_base.R
@@ -151,10 +151,10 @@ run_simulation_depvar <- function(df, result, outcome_formula=y~x+z, proxy_formu
     temp.df <- copy(df)
     temp.df[,y:=y.obs]
     mod.caroll.lik <- measerr_mle_dv(temp.df, outcome_formula=outcome_formula, proxy_formula=proxy_formula)
     temp.df <- copy(df)
     temp.df[,y:=y.obs]
     mod.caroll.lik <- measerr_mle_dv(temp.df, outcome_formula=outcome_formula, proxy_formula=proxy_formula)

-    fisher.info <- solve(mod.caroll.lik$hessian)
+    fischer.info <- solve(mod.caroll.lik$hessian)

     coef <- mod.caroll.lik$par
     coef <- mod.caroll.lik$par

-    ci.upper <- coef + sqrt(diag(fisher.info)) * 1.96
-    ci.lower <- coef - sqrt(diag(fisher.info)) * 1.96
+    ci.upper <- coef + sqrt(diag(fischer.info)) * 1.96
+    ci.lower <- coef - sqrt(diag(fischer.info)) * 1.96

     result <- append(result,
                      list(Bxy.est.mle = coef['x'],
                           Bxy.ci.upper.mle = ci.upper['x'],
     result <- append(result,
                      list(Bxy.est.mle = coef['x'],
                           Bxy.ci.upper.mle = ci.upper['x'],


@@ -180,26 +180,35 @@ run_simulation_depvar <- function(df, result, outcome_formula=y~x+z, proxy_formu
     
 
     # amelia says use normal distribution for binary variables.
     
 
     # amelia says use normal distribution for binary variables.

-
-    amelia.out.k <- amelia(df, m=200, p2s=0, idvars=c('y','ystar','w'))
-    mod.amelia.k <- zelig(y.obs~x+z, model='ls', data=amelia.out.k$imputations, cite=FALSE)
-    (coefse <- combine_coef_se(mod.amelia.k, messages=FALSE))
-    est.x.mi <- coefse['x','Estimate']
-    est.x.se <- coefse['x','Std.Error']
-    result <- append(result,
-                     list(Bxy.est.amelia.full = est.x.mi,
+    amelia_result <- list(Bxy.est.amelia.full = NA,
+                          Bxy.ci.upper.amelia.full = NA,
+                          Bxy.ci.lower.amelia.full = NA,
+                          Bzy.est.amelia.full = NA,
+                          Bzy.ci.upper.amelia.full = NA,
+                          Bzy.ci.lower.amelia.full = NA
+                          )
+
+    tryCatch({
+        amelia.out.k <- amelia(df, m=200, p2s=0, idvars=c('y','ystar','w'))
+        mod.amelia.k <- zelig(y.obs~x+z, model='ls', data=amelia.out.k$imputations, cite=FALSE)
+        (coefse <- combine_coef_se(mod.amelia.k, messages=FALSE))
+        est.x.mi <- coefse['x','Estimate']
+        est.x.se <- coefse['x','Std.Error']
+
+        est.z.mi <- coefse['z','Estimate']
+        est.z.se <- coefse['z','Std.Error']
+        amelia_result <- list(Bxy.est.amelia.full = est.x.mi,

                           Bxy.ci.upper.amelia.full = est.x.mi + 1.96 * est.x.se,
                           Bxy.ci.upper.amelia.full = est.x.mi + 1.96 * est.x.se,

-                          Bxy.ci.lower.amelia.full = est.x.mi - 1.96 * est.x.se
-                          ))
-
-    est.z.mi <- coefse['z','Estimate']
-    est.z.se <- coefse['z','Std.Error']
-
-    result <- append(result,
-                     list(Bzy.est.amelia.full = est.z.mi,
+                          Bxy.ci.lower.amelia.full = est.x.mi - 1.96 * est.x.se,
+                          Bzy.est.amelia.full = est.z.mi,

                           Bzy.ci.upper.amelia.full = est.z.mi + 1.96 * est.z.se,
                           Bzy.ci.lower.amelia.full = est.z.mi - 1.96 * est.z.se
                           Bzy.ci.upper.amelia.full = est.z.mi + 1.96 * est.z.se,
                           Bzy.ci.lower.amelia.full = est.z.mi - 1.96 * est.z.se

-                          ))
+                          )
+    },
+    error = function(e){
+    result[['error']] <- e}
+    )
+    result <- append(result,amelia_result)

 
     return(result)
 
 
     return(result)
 


@@ -299,11 +308,32 @@ run_simulation <-  function(df, result, outcome_formula=y~x+z, proxy_formula=NUL
         temp.df <- copy(df)
         temp.df <- temp.df[,x:=x.obs]
         mod.caroll.lik <- measerr_mle(temp.df, outcome_formula=outcome_formula, proxy_formula=proxy_formula, truth_formula=truth_formula)
         temp.df <- copy(df)
         temp.df <- temp.df[,x:=x.obs]
         mod.caroll.lik <- measerr_mle(temp.df, outcome_formula=outcome_formula, proxy_formula=proxy_formula, truth_formula=truth_formula)

-        fisher.info <- solve(mod.caroll.lik$hessian)
-        coef <- mod.caroll.lik$par
-        ci.upper <- coef + sqrt(diag(fisher.info)) * 1.96
-        ci.lower <- coef - sqrt(diag(fisher.info)) * 1.96
-        
+
+    ## tryCatch({
+    ## mod.calibrated.mle <- mecor(y ~ MeasError(w_pred, reference = x.obs) + z, df, B=400, method='efficient')
+    ## (mod.calibrated.mle)
+    ## (mecor.ci <- summary(mod.calibrated.mle)$c$ci['x.obs',])
+    ## result <- append(result, list(
+    ##                              Bxy.est.mecor = mecor.ci['Estimate'],
+    ##                              Bxy.ci.upper.mecor = mecor.ci['UCI'],
+    ##                              Bxy.ci.lower.mecor = mecor.ci['LCI'])
+    ##                  )
+
+
+
+    fischer.info <- NA
+    ci.upper <- NA
+    ci.lower <- NA
+
+    tryCatch({fischer.info <- solve(mod.caroll.lik$hessian)
+        ci.upper <- coef + sqrt(diag(fischer.info)) * 1.96
+        ci.lower <- coef - sqrt(diag(fischer.info)) * 1.96
+    },
+
+    error=function(e) {result[['error']] <- as.character(e)
+    })
+
+    coef <- mod.caroll.lik$par

         
         result <- append(result,
                          list(Bxy.est.mle = coef['x'],
         
         result <- append(result,
                          list(Bxy.est.mle = coef['x'],