library(argparser) library(mecor) library(ggplot2) library(data.table) library(filelock) library(arrow) library(Amelia) library(Zelig) library(predictionError) options(amelia.parallel="no", amelia.ncpus=1) setDTthreads(40) source("irr_dv_simulation_base.R") ## one way to do it is by adding correlation to x.obs and y that isn't in w. ## in other words, the model is missing an important feature of x.obs that's related to y. simulate_data <- function(N, m, B0, Bxy, Bzy, seed, prediction_accuracy=0.73, coder_accuracy=0.8){ set.seed(seed) # make w and y dependent z <- rbinom(N, 1, 0.5) x <- rbinom(N, 1, 0.5) ystar <- Bzy * z + Bxy * x + B0 y <- rbinom(N,1,plogis(ystar)) # glm(y ~ x + z, family="binomial") df <- data.table(x=x,y=y,ystar=ystar,z=z) df <- df[sample(nrow(df), m), y.obs := y] coder.0.correct <- rbinom(m, 1, coder_accuracy) coder.1.correct <- rbinom(m, 1, coder_accuracy) df[!is.na(y.obs),y.obs.0 := as.numeric((.SD$y.obs & coder.0.correct) | (!.SD$y.obs & !coder.0.correct))] df[!is.na(y.obs),y.obs.1 := as.numeric((.SD$y.obs & coder.1.correct) | (!.SD$y.obs & !coder.1.correct))] odds.y1 <- qlogis(prediction_accuracy) odds.y0 <- qlogis(prediction_accuracy,lower.tail=F) df[y==0,w:=plogis(rlogis(.N,odds.y0))] df[y==1,w:=plogis(rlogis(.N,odds.y1))] df[,w_pred := as.integer(w > 0.5)] print(mean(df$y == df$y.obs.0,na.rm=T)) print(mean(df$y == df$y.obs.1,na.rm=T)) print(mean(df[x==0]$y == df[x==0]$w_pred)) print(mean(df[x==1]$y == df[x==1]$w_pred)) print(mean(df$w_pred == df$y)) return(df) } parser <- arg_parser("Simulate data and fit corrected models") parser <- add_argument(parser, "--N", default=5000, help="number of observations of w") parser <- add_argument(parser, "--m", default=500, help="m the number of ground truth observations") parser <- add_argument(parser, "--seed", default=16, help='seed for the rng') parser <- add_argument(parser, "--outfile", help='output file', default='example_2.feather') parser <- add_argument(parser, "--y_explained_variance", help='what proportion of the variance of y can be explained?', default=0.1) parser <- add_argument(parser, "--prediction_accuracy", help='how accurate is the predictive model?', default=0.73) ## parser <- add_argument(parser, "--x_bias_y1", help='how is the classifier biased when y = 1?', default=-0.75) ## parser <- add_argument(parser, "--x_bias_y0", help='how is the classifier biased when y = 0 ?', default=0.75) parser <- add_argument(parser, "--Bxy", help='coefficient of x on y', default=0.3) parser <- add_argument(parser, "--Bzy", help='coeffficient of z on y', default=-0.3) parser <- add_argument(parser, "--outcome_formula", help='formula for the outcome variable', default="y~x+z") parser <- add_argument(parser, "--proxy_formula", help='formula for the proxy variable', default="w_pred~y+y.obs.1+y.obs.0") parser <- add_argument(parser, "--coder_accuracy", help='How accurate are the coders?', default=0.8) args <- parse_args(parser) B0 <- 0 Bxy <- args$Bxy Bzy <- args$Bzy df <- simulate_data(args$N, args$m, B0, Bxy, Bzy, args$seed, args$prediction_accuracy, args$coder_accuracy) # result <- list('N'=args$N,'m'=args$m,'B0'=B0,'Bxy'=Bxy,'Bzy'=Bzy, 'seed'=args$seed, 'y_explained_variance'=args$y_explained_variance, 'prediction_accuracy'=args$prediction_accuracy, 'x_bias_y0'=args$x_bias_y0,'x_bias_y1'=args$x_bias_y1,'outcome_formula' = args$outcome_formula, 'proxy_formula' = args$proxy_formula) result <- list('N'=args$N,'m'=args$m,'B0'=B0,'Bxy'=Bxy,'Bzy'=Bzy, 'seed'=args$seed, 'y_explained_variance'=args$y_explained_variance, 'prediction_accuracy'=args$prediction_accuracy, 'outcome_formula' = args$outcome_formula, 'proxy_formula' = args$proxy_formula) outline <- run_simulation_depvar(df, result, outcome_formula = as.formula(args$outcome_formula), proxy_formula = as.formula(args$proxy_formula)) outfile_lock <- lock(paste0(args$outfile, '_lock'),exclusive=TRUE) if(file.exists(args$outfile)){ logdata <- read_feather(args$outfile) logdata <- rbind(logdata,as.data.table(outline),fill=TRUE) } else { logdata <- as.data.table(outline) } print(outline) write_feather(logdata, args$outfile) unlock(outfile_lock) warnings()