add mle methods from carroll

[ml_measurement_error_public.git] / simulations / 01_two_covariates.R

diff --git a/simulations/01_two_covariates.R b/simulations/01_two_covariates.R
index 419403d20fd415bc55772c6460961f5a6ef91be4..7b8e12ee90f53f855782561209f6b83505967fc7 100644 (file)
--- a/simulations/01_two_covariates.R
+++ b/simulations/01_two_covariates.R
@@ -50,8 +50,8 @@ simulate_data <- function(N, m, B0=0, Bxy=0.2, Bgy=-0.2, Bgx=0.2, y_explained_va
     }
 
     df <- df[,w_pred:=x]
-
     df <- df[sample(1:N,(1-prediction_accuracy)*N),w_pred:=(w_pred-1)**2]
+    w <- predict(glm(x ~ w_pred,data=df,family=binomial(link='logit')),type='response')
     df <- df[,':='(w=w, w_pred = w_pred)]
     return(df)
 }
@@ -61,15 +61,20 @@ parser <- add_argument(parser, "--N", default=500, help="number of observations
 parser <- add_argument(parser, "--m", default=100, help="m the number of ground truth observations")
 parser <- add_argument(parser, "--seed", default=4321, help='seed for the rng')
 parser <- add_argument(parser, "--outfile", help='output file', default='example_1.feather')
+parser <- add_argument(parser, "--y_explained_variance", help='what proportion of the variance of y can be explained?', default=0.005)
+parser <- add_argument(parser, "--gx_explained_variance", help='what proportion of the variance of x can be explained by g?', default=0.15)
+parser <- add_argument(parser, "--prediction_accuracy", help='how accurate is the predictive model?', default=0.73)
+
 args <- parse_args(parser)
 
 B0 <- 0
 Bxy <- 0.2
 Bgy <- -0.2
-Bgx <- 0.5
+Bgx <- 0.4
+
+df <- simulate_data(args$N, args$m, B0, Bxy, Bgy, Bgx, seed=args$seed, y_explained_variance = args$y_explained_variance, gx_explained_variance = args$gx_explained_variance, prediction_accuracy=args$prediction_accuracy)
 
-df <- simulate_data(args$N, args$m, B0, Bxy, Bgy, Bgx, seed=args$seed, y_explained_variance = 0.025, gx_explained_variance = 0.15)
-result <- list('N'=args$N,'m'=args$m,'B0'=B0,'Bxy'=Bxy,'Bgy'=Bgy, 'Bgx'=Bgx, 'seed'=args$seed)
+result <- list('N'=args$N,'m'=args$m,'B0'=B0,'Bxy'=Bxy,'Bgy'=Bgy, 'Bgx'=Bgx, 'seed'=args$seed, 'y_explained_variance' = args$y_explained_variance, 'gx_explained_variance' = args$gx_explained_variance, "prediction_accuracy"=args$prediction_accuracy)
 outline <- run_simulation(df, result)
 
 outfile_lock <- lock(paste0(args$outfile, '_lock'),exclusive=TRUE)