X-Git-Url: https://code.communitydata.science/ml_measurement_error_public.git/blobdiff_plain/588bdd7ed74cf8fe8fd0f15df58a6a40c26ebae5..2cd447c327744263d5f94b20e1146cdf31b2ec2c:/simulations/01_two_covariates.R?ds=sidebyside diff --git a/simulations/01_two_covariates.R b/simulations/01_two_covariates.R index 7b8e12e..73e8939 100644 --- a/simulations/01_two_covariates.R +++ b/simulations/01_two_covariates.R @@ -1,8 +1,10 @@ -### EXAMPLE 2_b: demonstrates how measurement error can lead to a type sign error in a covariate -### This is the same as example 2, only instead of x->k we have k->x. -### Even when you have a good predictor, if it's biased against a covariate you can get the wrong sign. -### Even when you include the proxy variable in the regression. -### But with some ground truth and multiple imputation, you can fix it. +### EXAMPLE 2_b: demonstrates how measurement error can lead to a type +### sign error in a covariate This is the same as example 2, only +### instead of x->k we have k->x. Even when you have a good +### predictor, if it's biased against a covariate you can get the +### wrong sign. Even when you include the proxy variable in the +### regression. But with some ground truth and multiple imputation, +### you can fix it. library(argparser) library(mecor) @@ -12,9 +14,9 @@ library(filelock) library(arrow) library(Amelia) library(Zelig) + library(predictionError) -options(amelia.parallel="no", - amelia.ncpus=1) +options(amelia.parallel="no", amelia.ncpus=1) source("simulation_base.R") @@ -28,20 +30,18 @@ source("simulation_base.R") #### how much power do we get from the model in the first place? (sweeping N and m) #### -simulate_data <- function(N, m, B0=0, Bxy=0.2, Bgy=-0.2, Bgx=0.2, y_explained_variance=0.025, gx_explained_variance=0.15, prediction_accuracy=0.73, seed=1){ +simulate_data <- function(N, m, B0=0, Bxy=0.2, Bzy=-0.2, Bzx=0.2, y_explained_variance=0.025, prediction_accuracy=0.73, seed=1){ set.seed(seed) - g <- rbinom(N, 1, 0.5) - - x.var.epsilon <- var(Bgx *g) * ((1-gx_explained_variance)/gx_explained_variance) - x.epsilon <- rnorm(N,sd=sqrt(x.var.epsilon)) - xprime <- Bgx * g + x.epsilon - x <- as.integer(logistic(scale(xprime)) > 0.5) + z <- rbinom(N, 1, 0.5) + # x.var.epsilon <- var(Bzx *z) * ((1-zx_explained_variance)/zx_explained_variance) + xprime <- Bzx * z #+ x.var.epsilon + x <- rbinom(N,1,plogis(xprime)) - y.var.epsilon <- (var(Bgy * g) + var(Bxy *x) + 2*cov(Bxy*x,Bgy*g)) * ((1-y_explained_variance)/y_explained_variance) + y.var.epsilon <- (var(Bzy * z) + var(Bxy *x) + 2*cov(Bxy*x,Bzy*z)) * ((1-y_explained_variance)/y_explained_variance) y.epsilon <- rnorm(N, sd = sqrt(y.var.epsilon)) - y <- Bgy * g + Bxy * x + y.epsilon + y <- Bzy * z + Bxy * x + y.epsilon - df <- data.table(x=x,xprime=xprime,y=y,g=g) + df <- data.table(x=x,y=y,z=z) if(m < N){ df <- df[sample(nrow(df), m), x.obs := x] @@ -49,42 +49,53 @@ simulate_data <- function(N, m, B0=0, Bxy=0.2, Bgy=-0.2, Bgx=0.2, y_explained_va df <- df[, x.obs := x] } - 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)] + ## how can you make a model with a specific accuracy? + w0 =(1-x)**2 + (-1)**(1-x) * prediction_accuracy + + ## how can you make a model with a specific accuracy, with a continuous latent variable. + # now it makes the same amount of mistake to each point, probably + # add mean0 noise to the odds. + + w.noisey.odds = rlogis(N,qlogis(w0)) + df[,w := plogis(w.noisey.odds)] + df[,w_pred:=as.integer(w > 0.5)] + (mean(df$x==df$w_pred)) return(df) } 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=1000, help="number of observations of w") +parser <- add_argument(parser, "--m", default=200, help="m the number of ground truth observations") +parser <- add_argument(parser, "--seed", default=57, 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, "--y_explained_variance", help='what proportion of the variance of y can be explained?', default=0.05) +# parser <- add_argument(parser, "--zx_explained_variance", help='what proportion of the variance of x can be explained by z?', default=0.3) parser <- add_argument(parser, "--prediction_accuracy", help='how accurate is the predictive model?', default=0.73) - +parser <- add_argument(parser, "--Bzx", help='coefficient of z on x?', default=1) args <- parse_args(parser) B0 <- 0 -Bxy <- 0.2 -Bgy <- -0.2 -Bgx <- 0.4 +Bxy <- 0.3 +Bzy <- -0.3 +Bzx <- args$Bzx -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) +if (args$m < args$N){ -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) + df <- simulate_data(args$N, args$m, B0, Bxy, Bzy, Bzx, seed=args$seed + 500, y_explained_variance = args$y_explained_variance, prediction_accuracy=args$prediction_accuracy) -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)) -} else { - logdata <- as.data.table(outline) -} + result <- list('N'=args$N,'m'=args$m,'B0'=B0,'Bxy'=Bxy,'Bzy'=Bzy, 'Bzx'=Bzx, 'seed'=args$seed, 'y_explained_variance' = args$y_explained_variance, 'zx_explained_variance' = args$zx_explained_variance, "prediction_accuracy"=args$prediction_accuracy, "error"="") -print(outline) -write_feather(logdata, args$outfile) -unlock(outfile_lock) + outline <- run_simulation(df, result) + + 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) +}