X-Git-Url: https://code.communitydata.science/ml_measurement_error_public.git/blobdiff_plain/47e9367ed5c61b721bdc17cddd76bced4f8ed621..d9d3e47a44ddead1cdf7a649bc0e9849c2219498:/simulations/plot_example.R?ds=sidebyside diff --git a/simulations/plot_example.R b/simulations/plot_example.R index 7a853b7..4da045b 100644 --- a/simulations/plot_example.R +++ b/simulations/plot_example.R @@ -5,52 +5,58 @@ library(ggplot2) library(filelock) library(argparser) +source("summarize_estimator.R") + + parser <- arg_parser("Simulate data and fit corrected models.") -parser <- add_argument(parser, "--infile", default="", help="name of the file to read.") +parser <- add_argument(parser, "--infile", default="robustness_2.feather", help="name of the file to read.") +parser <- add_argument(parser, "--remember-file", default="remembr.RDS", help="name of the remember file.") parser <- add_argument(parser, "--name", default="", help="The name to safe the data to in the remember file.") args <- parse_args(parser) -summarize.estimator <- function(df, suffix='naive', coefname='x'){ - - part <- df[,c('N', - 'm', - 'Bxy', - paste0('B',coefname,'y.est.',suffix), - paste0('B',coefname,'y.ci.lower.',suffix), - paste0('B',coefname,'y.ci.upper.',suffix), - 'y_explained_variance', - 'Bzx', - 'Bzy', - 'accuracy_imbalance_difference' - ), - with=FALSE] + + +## summarize.estimator <- function(df, suffix='naive', coefname='x'){ + +## part <- df[,c('N', +## 'm', +## 'Bxy', +## paste0('B',coefname,'y.est.',suffix), +## paste0('B',coefname,'y.ci.lower.',suffix), +## paste0('B',coefname,'y.ci.upper.',suffix), +## 'y_explained_variance', +## 'Bzx', +## 'Bzy', +## 'accuracy_imbalance_difference' +## ), +## with=FALSE] - true.in.ci <- as.integer((part$Bxy >= part[[paste0('B',coefname,'y.ci.lower.',suffix)]]) & (part$Bxy <= part[[paste0('B',coefname,'y.ci.upper.',suffix)]])) - zero.in.ci <- as.integer(0 >= part[[paste0('B',coefname,'y.ci.lower.',suffix)]]) & (0 <= part[[paste0('B',coefname,'y.ci.upper.',suffix)]]) - bias <- part$Bxy - part[[paste0('B',coefname,'y.est.',suffix)]] - sign.correct <- as.integer(sign(part$Bxy) == sign(part[[paste0('B',coefname,'y.est.',suffix)]])) - - part <- part[,':='(true.in.ci = true.in.ci, - zero.in.ci = zero.in.ci, - bias=bias, - sign.correct =sign.correct)] - - part.plot <- part[, .(p.true.in.ci = mean(true.in.ci), - mean.bias = mean(bias), - mean.est = mean(.SD[[paste0('B',coefname,'y.est.',suffix)]]), - var.est = var(.SD[[paste0('B',coefname,'y.est.',suffix)]]), - est.upper.95 = quantile(.SD[[paste0('B',coefname,'y.est.',suffix)]],0.95,na.rm=T), - est.lower.95 = quantile(.SD[[paste0('B',coefname,'y.est.',suffix)]],0.05,na.rm=T), - N.sims = .N, - p.sign.correct = mean(as.integer(sign.correct & (! zero.in.ci))), - variable=coefname, - method=suffix - ), - by=c("N","m",'y_explained_variance','Bzx', 'Bzy', 'accuracy_imbalance_difference') - ] +## true.in.ci <- as.integer((part$Bxy >= part[[paste0('B',coefname,'y.ci.lower.',suffix)]]) & (part$Bxy <= part[[paste0('B',coefname,'y.ci.upper.',suffix)]])) +## zero.in.ci <- as.integer(0 >= part[[paste0('B',coefname,'y.ci.lower.',suffix)]]) & (0 <= part[[paste0('B',coefname,'y.ci.upper.',suffix)]]) +## bias <- part$Bxy - part[[paste0('B',coefname,'y.est.',suffix)]] +## sign.correct <- as.integer(sign(part$Bxy) == sign(part[[paste0('B',coefname,'y.est.',suffix)]])) + +## part <- part[,':='(true.in.ci = true.in.ci, +## zero.in.ci = zero.in.ci, +## bias=bias, +## sign.correct =sign.correct)] + +## part.plot <- part[, .(p.true.in.ci = mean(true.in.ci), +## mean.bias = mean(bias), +## mean.est = mean(.SD[[paste0('B',coefname,'y.est.',suffix)]]), +## var.est = var(.SD[[paste0('B',coefname,'y.est.',suffix)]]), +## est.upper.95 = quantile(.SD[[paste0('B',coefname,'y.est.',suffix)]],0.95,na.rm=T), +## est.lower.95 = quantile(.SD[[paste0('B',coefname,'y.est.',suffix)]],0.05,na.rm=T), +## N.sims = . +## p.sign.correct = mean(as.integer(sign.correct & (! zero.in.ci))), +## variable=coefname, +## method=suffix +## ), +## by=c("N","m",'y_explained_variance','Bzx', 'Bzy', 'accuracy_imbalance_difference') +## ] - return(part.plot) -} +## return(part.plot) +## } build_plot_dataset <- function(df){ @@ -70,13 +76,13 @@ build_plot_dataset <- function(df){ z.amelia.full <- summarize.estimator(df, 'amelia.full', 'z') - x.mecor <- summarize.estimator(df, 'mecor', 'x') + ## x.mecor <- summarize.estimator(df, 'mecor', 'x') - z.mecor <- summarize.estimator(df, 'mecor', 'z') + ## z.mecor <- summarize.estimator(df, 'mecor', 'z') - x.mecor <- summarize.estimator(df, 'mecor', 'x') + ## x.mecor <- summarize.estimator(df, 'mecor', 'x') - z.mecor <- summarize.estimator(df, 'mecor', 'z') + ## z.mecor <- summarize.estimator(df, 'mecor', 'z') x.mle <- summarize.estimator(df, 'mle', 'x') @@ -91,31 +97,48 @@ build_plot_dataset <- function(df){ z.gmm <- summarize.estimator(df, 'gmm', 'z') accuracy <- df[,mean(accuracy)] - plot.df <- rbindlist(list(x.true,z.true,x.naive,z.naive,x.amelia.full,z.amelia.full,x.mecor, z.mecor, x.gmm, z.gmm, x.feasible, z.feasible,z.mle, x.mle, x.zhang, z.zhang, x.gmm, z.gmm),use.names=T) + plot.df <- rbindlist(list(x.true,z.true,x.naive,z.naive,x.amelia.full,z.amelia.full,x.gmm, z.gmm, x.feasible, z.feasible,z.mle, x.mle, x.zhang, z.zhang, x.gmm, z.gmm),use.names=T) plot.df[,accuracy := accuracy] plot.df <- plot.df[,":="(sd.est=sqrt(var.est)/N.sims)] return(plot.df) } -plot.df <- read_feather(args$infile) -print(unique(plot.df$N)) +sims.df <- read_feather(args$infile) +unique(sims.df[,.N,by=.(N,m)]) +print(unique(sims.df$N)) # df <- df[apply(df,1,function(x) !any(is.na(x)))] -if(!('Bzx' %in% names(plot.df))) - plot.df[,Bzx:=NA] +if(!('Bzx' %in% names(sims.df))) + sims.df[,Bzx:=NA] -if(!('accuracy_imbalance_difference' %in% names(plot.df))) - plot.df[,accuracy_imbalance_difference:=NA] +if(!('accuracy_imbalance_difference' %in% names(sims.df))) + sims.df[,accuracy_imbalance_difference:=NA] -unique(plot.df[,'accuracy_imbalance_difference']) +unique(sims.df[,'accuracy_imbalance_difference']) +change.remember.file(args$remember_file, clear=TRUE) #plot.df <- build_plot_dataset(df[accuracy_imbalance_difference==0.1][N==700]) -plot.df <- build_plot_dataset(plot.df) +plot.df <- build_plot_dataset(sims.df) remember(plot.df,args$name) +set.remember.prefix(gsub("plot.df.","",args$name)) + +remember(median(sims.df$cor.xz),'med.cor.xz') +remember(median(sims.df$accuracy),'med.accuracy') +remember(median(sims.df$accuracy.y0),'med.accuracy.y0') +remember(median(sims.df$accuracy.y1),'med.accuracy.y1') +remember(median(sims.df$fpr),'med.fpr') +remember(median(sims.df$fpr.y0),'med.fpr.y0') +remember(median(sims.df$fpr.y1),'med.fpr.y1') +remember(median(sims.df$fnr),'med.fnr') +remember(median(sims.df$fnr.y0),'med.fnr.y0') +remember(median(sims.df$fnr.y1),'med.fnr.y1') + +remember(median(sims.df$cor.resid.w_pred),'cor.resid.w_pred') + #ggplot(df,aes(x=Bxy.est.mle)) + geom_histogram() + facet_grid(accuracy_imbalance_difference ~ Bzy) ## ## ## df[gmm.ER_pval<0.05]