update simulation and mle code

[ml_measurement_error_public.git] / simulations / summarize_estimator.R

diff --git a/simulations/summarize_estimator.R b/simulations/summarize_estimator.R
new file mode 100644 (file)
index 0000000..8199c06
--- /dev/null
+++ b/simulations/summarize_estimator.R
@@ -0,0 +1,42 @@
+
+summarize.estimator <- function(df, suffix='naive', coefname='x'){
+
+    part <- df[,c('N',
+                  'm',
+                  'Bxy',
+                  paste0('B',coefname,'y.est.',suffix),
+                  paste0('B',coefname,'y.ci.lower.',suffix),
+                  paste0('B',coefname,'y.ci.upper.',suffix),
+                  'y_explained_variance',
+                  'Bzx',
+                  'Bzy',
+                  'accuracy_imbalance_difference'
+                  ),
+               with=FALSE]
+    
+    true.in.ci <- as.integer((part$Bxy >= part[[paste0('B',coefname,'y.ci.lower.',suffix)]]) & (part$Bxy <= part[[paste0('B',coefname,'y.ci.upper.',suffix)]]))
+    zero.in.ci <- as.integer(0 >= part[[paste0('B',coefname,'y.ci.lower.',suffix)]]) & (0 <= part[[paste0('B',coefname,'y.ci.upper.',suffix)]])
+    bias <- part$Bxy - part[[paste0('B',coefname,'y.est.',suffix)]]
+    sign.correct <- as.integer(sign(part$Bxy) == sign(part[[paste0('B',coefname,'y.est.',suffix)]]))
+
+    part <- part[,':='(true.in.ci = true.in.ci,
+                       zero.in.ci = zero.in.ci,
+                       bias=bias,
+                       sign.correct =sign.correct)]
+
+    part.plot <- part[, .(p.true.in.ci = mean(true.in.ci),
+                          mean.bias = mean(bias),
+                          mean.est = mean(.SD[[paste0('B',coefname,'y.est.',suffix)]]),
+                          var.est = var(.SD[[paste0('B',coefname,'y.est.',suffix)]]),
+                          est.upper.95 = quantile(.SD[[paste0('B',coefname,'y.est.',suffix)]],0.95,na.rm=T),
+                          est.lower.95 = quantile(.SD[[paste0('B',coefname,'y.est.',suffix)]],0.05,na.rm=T),
+                          N.sims = .N,
+                          p.sign.correct = mean(as.integer(sign.correct & (! zero.in.ci))),
+                          variable=coefname,
+                          method=suffix
+                          ),
+                      by=c("N","m",'y_explained_variance','Bzx', 'Bzy', 'accuracy_imbalance_difference')
+                      ]
+    
+    return(part.plot)
+}