simulations/example_2.R

   1 ### EXAMPLE 2: demonstrates how measurement error can lead to a type sign error in a covariate
   2 ### Even when you have a good predictor, if it's biased against a covariate you can get the wrong sign.
   3 ### Even when you include the proxy variable in the regression.
   4 ### But with some ground truth and multiple imputation, you can fix it.
   5
   6 library(argparser)
   7 library(mecor)
   8 library(ggplot2)
   9 library(data.table)
  10 library(filelock)
  11 library(arrow)
  12 library(Amelia)
  13 library(Zelig)
  14 library(predictionError)
  15 options(amelia.parallel="no",
  16         amelia.ncpus=1)
  17
  18 ## SETUP:
  19 ### we want to estimate x -> y; x is MAR
  20 ### we have x -> k; k -> w; x -> w is used to predict x via the model w.
  21 ### A realistic scenario is that we have an NLP model predicting something like "racial harassment" in social media comments
  22 ### The labels x are binary, but the model provides a continuous predictor
  23
  24 ### simulation:
  25 #### how much power do we get from the model in the first place? (sweeping N and m)
  26 ####
  27 source("simulation_base.R")
  28
  29 simulate_latent_cocause <- function(N, m, B0, Bxy, Bgy, Bkx, Bgx, seed){
  30     set.seed(seed)
  31
  32     ## the true value of x
  33
  34     g <- rbinom(N, 1, 0.5)
  35     k <- rnorm(N, 0, 1)
  36     xprime <- Bkx*k + Bgx * g + rnorm(N,0,1)
  37     xvec <- scale(xprime)
  38
  39     y <-  Bxy * xvec  + Bgy * g + rnorm(N, 0, 1) + B0
  40
  41     df <- data.table(x=xvec,k=k,y=y,g=g)
  42     names(df) <- c('x','k','y','g')
  43     if( m < N){
  44         df <- df[sample(nrow(df), m), x.obs := x]
  45     } else {
  46         df <- df[, x.obs := x]
  47     }
  48
  49     w.model <- lm(x ~ k,df)
  50     w <- predict(w.model,data.frame(k=k))
  51     w <- logistic(w + rnorm(N,0,sd(w)*0.1))
  52     ## y = B0 + B1x + e
  53
  54     df[,':='(w=w, w_pred = as.integer(w>0.5))]
  55     return(df)
  56 }
  57
  58
  59 parser <- arg_parser("Simulate data and fit corrected models")
  60 parser <- add_argument(parser, "--N", default=1000, help="number of observations of w")
  61 parser <- add_argument(parser, "--m", default=200, help="m the number of ground truth observations")
  62 parser <- add_argument(parser, "--seed", default=4321, help='seed for the rng')
  63 parser <- add_argument(parser, "--outfile", help='output file', default='example_2.feather')
  64 args <- parse_args(parser)
  65
  66 Ns <- c(1000, 10000, 1e6)
  67 ms <- c(100, 250, 500, 1000)
  68 B0 <- 0
  69 Bxy <- 0.2
  70 Bgy <- -0.2
  71 Bkx <- 2
  72 Bgx <- 3
  73
  74 outline <- run_simulation(simulate_latent_cocause(args$N, args$m, B0, Bxy, Bgy, Bkx, Bgx, args$seed)
  75                          ,list('N'=args$N,'m'=args$m,'B0'=B0,'Bxy'=Bxy,'Bgy'=Bgy, 'Bkx'=Bkx, 'Bgx'=Bgx, 'seed'=args$seed))
  76
  77 outfile_lock <- lock(paste0(args$outfile, '_lock'),exclusive=TRUE)
  78 if(file.exists(args$outfile)){
  79     logdata <- read_feather(args$outfile)
  80     logdata <- rbind(logdata,as.data.table(outline))
  81 } else {
  82     logdata <- as.data.table(outline)
  83 }
  84
  85 print(outline)
  86 write_feather(logdata, args$outfile)
  87 unlock(outfile_lock)