update plotting code and makefile

[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 c52a3dc3c45e7d1769455adf17225d1b16a2aba3..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[,w_pred:=x]

-

     df <- df[sample(1:N,(1-prediction_accuracy)*N),w_pred:=(w_pred-1)**2]
     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)
 }
     df <- df[,':='(w=w, w_pred = w_pred)]
     return(df)
 }


@@ -60,16 +60,21 @@ parser <- arg_parser("Simulate data and fit corrected models")
 parser <- add_argument(parser, "--N", default=500, help="number of observations of w")
 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, "--N", default=500, help="number of observations of w")
 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_2_B.feather')
+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
 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)
 outline <- run_simulation(df, result)
 
 outfile_lock <- lock(paste0(args$outfile, '_lock'),exclusive=TRUE)