-### 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)
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")
#### 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]
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)
+}
projects / ml_measurement_error_public.git / blobdiff