X-Git-Url: https://code.communitydata.science/ml_measurement_error_public.git/blobdiff_plain/47e9367ed5c61b721bdc17cddd76bced4f8ed621..d0c5766bdf867a81a2477d2cac1d40812110af90:/simulations/simulation_base.R

diff --git a/simulations/simulation_base.R b/simulations/simulation_base.R
index ee46ec6..27f0276 100644
--- a/simulations/simulation_base.R
+++ b/simulations/simulation_base.R
@@ -210,11 +210,19 @@ run_simulation_depvar <- function(df, result, outcome_formula=y~x+z, proxy_formu
 
     accuracy <- df[,mean(w_pred==y)]
     result <- append(result, list(accuracy=accuracy))
+    error.cor.x <- cor(df$x, df$w - df$x)
+    result <- append(result, list(error.cor.x = error.cor.x))
 
+    model.null <- glm(y~1, data=df,family=binomial(link='logit'))
     (model.true <- glm(y ~ x + z, data=df,family=binomial(link='logit')))
+    (lik.ratio <- exp(logLik(model.true) - logLik(model.null)))
+
     true.ci.Bxy <- confint(model.true)['x',]
     true.ci.Bzy <- confint(model.true)['z',]
 
+
+    result <- append(result, list(lik.ratio=lik.ratio))
+
     result <- append(result, list(Bxy.est.true=coef(model.true)['x'],
                                   Bzy.est.true=coef(model.true)['z'],
                                   Bxy.ci.upper.true = true.ci.Bxy[2],
@@ -322,8 +330,33 @@ run_simulation_depvar <- function(df, result, outcome_formula=y~x+z, proxy_formu
 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.y0 <- df[y<=0,mean(w_pred==x)]
+    accuracy.y1 <- df[y>=0,mean(w_pred==x)]
+    cor.y.xi <- cor(df$x - df$w_pred, df$y)
+
+    fnr <- df[w_pred==0,mean(w_pred!=x)]
+    fnr.y0 <- df[(w_pred==0) & (y<=0),mean(w_pred!=x)]
+    fnr.y1 <- df[(w_pred==0) & (y>=0),mean(w_pred!=x)]
+
+    fpr <- df[w_pred==1,mean(w_pred!=x)]
+    fpr.y0 <- df[(w_pred==1) & (y<=0),mean(w_pred!=x)]
+    fpr.y1 <- df[(w_pred==1) & (y>=0),mean(w_pred!=x)]
+    cor.resid.w_pred <- cor(resid(lm(y~x+z,df)),df$w_pred)
+
+    result <- append(result, list(accuracy=accuracy,
+                                  accuracy.y0=accuracy.y0,
+                                  accuracy.y1=accuracy.y1,
+                                  cor.y.xi=cor.y.xi,
+                                  fnr=fnr,
+                                  fnr.y0=fnr.y0,
+                                  fnr.y1=fnr.y1,
+                                  fpr=fpr,
+                                  fpr.y0=fpr.y0,
+                                  fpr.y1=fpr.y1,
+                                  cor.resid.w_pred=cor.resid.w_pred
+                                  ))
+
+    result <- append(result, list(cor.xz=cor(df$x,df$z)))
     (model.true <- lm(y ~ x + z, data=df))
     true.ci.Bxy <- confint(model.true)['x',]
     true.ci.Bzy <- confint(model.true)['z',]