simulations/plot_example_1.R

   1 source("RemembR/R/RemembeR.R")
   2 library(arrow)
   3 library(data.table)
   4 library(ggplot2)
   5
   6 df <- data.table(read_feather("example_1.feather"))
   7
   8 x.naive <- df[,.(N, m, Bxy, Bxy.est.naive, Bxy.ci.lower.naive, Bxy.ci.upper.naive)]
   9 x.naive <- x.naive[,':='(true.in.ci = as.integer((Bxy >= Bxy.ci.lower.naive) & (Bxy <= Bxy.ci.upper.naive)),
  10                          zero.in.ci = (0 >= Bxy.ci.lower.naive) & (0 <= Bxy.ci.upper.naive),
  11                          bias = Bxy - Bxy.est.naive,
  12                          sign.correct = as.integer(sign(Bxy) == sign(Bxy.est.naive)))]
  13
  14 x.naive.plot <- x.naive[,.(p.true.in.ci = mean(true.in.ci),
  15                            mean.bias = mean(bias),
  16                            p.sign.correct = mean(as.integer(sign.correct & (! zero.in.ci))),
  17                            mean.estimate=mean(Bxy.est.naive),
  18                            var.estimate=var(Bxy.est.naive),
  19                            Bxy=mean(Bxy),
  20                            variable='x',
  21                            method='Naive'
  22                            ),
  23                         by=c('N','m')]
  24
  25 x.amelia.full <- df[,.(N, m, Bxy, Bxy.est.true, Bxy.ci.lower.amelia.full, Bxy.ci.upper.amelia.full, Bxy.est.amelia.full)]
  26
  27 x.amelia.full <- x.amelia.full[,':='(true.in.ci = (Bxy.est.true >= Bxy.ci.lower.amelia.full) & (Bxy.est.true <= Bxy.ci.upper.amelia.full),
  28                                      zero.in.ci = (0 >= Bxy.ci.lower.amelia.full) & (0 <= Bxy.ci.upper.amelia.full),
  29                                      bias = Bxy.est.true - Bxy.est.amelia.full,
  30                                      sign.correct = sign(Bxy.est.true) == sign(Bxy.est.amelia.full))]
  31
  32 x.amelia.full.plot <- x.amelia.full[,.(p.true.in.ci = mean(as.integer(true.in.ci)),
  33                                        mean.bias = mean(bias),
  34                                        p.sign.correct = mean(as.integer(sign.correct & (! zero.in.ci))),
  35                                        variable='x',
  36                                        method='Multiple imputation'
  37                                        ),
  38                                     by=c('N','m')]
  39
  40
  41
  42 g.amelia.full <- df[,.(N, m, Bgy.est.true, Bgy.est.amelia.full, Bgy.ci.lower.amelia.full, Bgy.ci.upper.amelia.full)]
  43 g.amelia.full <- g.amelia.full[,':='(true.in.ci = (Bgy.est.true >= Bgy.ci.lower.amelia.full) & (Bgy.est.true <= Bgy.ci.upper.amelia.full),
  44                                      zero.in.ci = (0 >= Bgy.ci.lower.amelia.full) & (0 <= Bgy.ci.upper.amelia.full),
  45                                      bias =  Bgy.est.amelia.full - Bgy.est.true,
  46                                      sign.correct = sign(Bgy.est.true) == sign(Bgy.est.amelia.full))]
  47
  48 g.amelia.full.plot <- g.amelia.full[,.(p.true.in.ci = mean(as.integer(true.in.ci)),
  49                                        mean.bias = mean(bias),
  50                                        p.sign.correct = mean(as.integer(sign.correct & (! zero.in.ci))),
  51                                        variable='g',
  52                                        method='Multiple imputation'
  53                                        ),
  54                                     by=c('N','m')]
  55
  56 x.amelia.nok <- df[,.(N, m, Bxy.est.true, Bxy.est.amelia.nok, Bxy.ci.lower.amelia.nok, Bxy.ci.upper.amelia.nok)]
  57 x.amelia.nok <- x.amelia.nok[,':='(true.in.ci = (Bxy.est.true >= Bxy.ci.lower.amelia.nok) & (Bxy.est.true <= Bxy.ci.upper.amelia.nok),
  58                                      zero.in.ci = (0 >= Bxy.ci.lower.amelia.nok) & (0 <= Bxy.ci.upper.amelia.nok),
  59                                      bias =  Bxy.est.amelia.nok - Bxy.est.true,
  60                                      sign.correct = sign(Bxy.est.true) == sign(Bxy.est.amelia.nok))]
  61
  62  x.amelia.nok.plot <- x.amelia.nok[,.(p.true.in.ci = mean(as.integer(true.in.ci)),
  63                                        mean.bias = mean(bias),
  64                                      p.sign.correct = mean(as.integer(sign.correct & (! zero.in.ci))),
  65                                      variable='x',
  66                                      method='Multiple imputation (Classifier features unobserved)'
  67                                      ),
  68                                     by=c('N','m')]
  69
  70 g.amelia.nok <- df[,.(N, m, Bgy.est.true, Bgy.est.amelia.nok, Bgy.ci.lower.amelia.nok, Bgy.ci.upper.amelia.nok)]
  71 g.amelia.nok <- g.amelia.nok[,':='(true.in.ci = (Bgy.est.true >= Bgy.ci.lower.amelia.nok) & (Bgy.est.true <= Bgy.ci.upper.amelia.nok),
  72                                      zero.in.ci = (0 >= Bgy.ci.lower.amelia.nok) & (0 <= Bgy.ci.upper.amelia.nok),
  73                                      bias =  Bgy.est.amelia.nok - Bgy.est.true,
  74                                      sign.correct = sign(Bgy.est.true) == sign(Bgy.est.amelia.nok))]
  75
  76 g.amelia.nok.plot <- g.amelia.nok[,.(p.true.in.ci = mean(as.integer(true.in.ci)),
  77                                        mean.bias = mean(bias),
  78                                      p.sign.correct = mean(as.integer(sign.correct & (! zero.in.ci))),
  79                                      variable='g',
  80                                      method='Multiple imputation (Classifier features unobserved)'
  81                                      ),
  82                                     by=c('N','m')]
  83
  84
  85 x.mecor <- df[,.(N,m,Bxy.est.true, Bxy.est.mecor,Bxy.lower.mecor, Bxy.upper.mecor)]
  86
  87 x.mecor <- x.mecor[,':='(true.in.ci = (Bxy.est.true >= Bxy.lower.mecor) & (Bxy.est.true <= Bxy.upper.mecor),
  88                                      zero.in.ci = (0 >= Bxy.lower.mecor) & (0 <= Bxy.upper.mecor),
  89                                      bias =  Bxy.est.mecor - Bxy.est.true,
  90                                      sign.correct = sign(Bxy.est.true) == sign(Bxy.est.mecor))]
  91
  92 x.mecor.plot <- x.mecor[,.(p.true.in.ci = mean(as.integer(true.in.ci)),
  93                                        mean.bias = mean(bias),
  94                                      p.sign.correct = mean(as.integer(sign.correct & (! zero.in.ci))),
  95                                      variable='x',
  96                                      method='Regression Calibration'
  97                                      ),
  98                                     by=c('N','m')]
  99
 100 g.mecor <- df[,.(N,m,Bgy.est.true, Bgy.est.mecor,Bgy.lower.mecor, Bgy.upper.mecor)]
 101
 102 g.mecor <- g.mecor[,':='(true.in.ci = (Bgy.est.true >= Bgy.lower.mecor) & (Bgy.est.true <= Bgy.upper.mecor),
 103                                      zero.in.ci = (0 >= Bgy.lower.mecor) & (0 <= Bgy.upper.mecor),
 104                                      bias =  Bgy.est.mecor - Bgy.est.true,
 105                                      sign.correct = sign(Bgy.est.true) == sign(Bgy.est.mecor))]
 106
 107 g.mecor.plot <- g.mecor[,.(p.true.in.ci = mean(as.integer(true.in.ci)),
 108                                        mean.bias = mean(bias),
 109                                      p.sign.correct = mean(as.integer(sign.correct & (! zero.in.ci))),
 110                                      variable='g',
 111                                      method='Regression Calibration'
 112                                      ),
 113                                     by=c('N','m')]
 114
 115 ## x.mecor <- df[,.(N,m,Bgy.est.true, Bgy.est.mecor,Bgy.ci.lower.mecor, Bgy.ci.upper.mecor)]
 116
 117 ## x.mecor <- x.mecor[,':='(true.in.ci = (Bgy.est.true >= Bgy.ci.lower.mecor) & (Bgy.est.true <= Bgy.ci.upper.mecor),
 118 ##                                      zero.in.ci = (0 >= Bgy.ci.lower.mecor) & (0 <= Bgy.ci.upper.mecor),
 119 ##                                      bias =  Bgy.est.mecor - Bgy.est.true,
 120 ##                                      sign.correct = sign(Bgy.est.true) == sign(Bgy.est.mecor))]
 121
 122 ## x.mecor.plot <- x.mecor[,.(p.true.in.ci = mean(as.integer(true.in.ci)),
 123 ##                                        mean.bias = mean(bias),
 124 ##                                      p.sign.correct = mean(as.integer(sign.correct & (! zero.in.ci))),
 125 ##                                      variable='g',
 126 ##                                      method='Regression Calibration'
 127 ##                                      ),
 128 ##                                     by=c('N','m')]
 129
 130
 131 x.gmm <- df[,.(N,m,Bxy.est.true, Bxy.est.gmm,Bxy.ci.lower.gmm, Bxy.ci.upper.gmm)]
 132 x.gmm <- x.gmm[,':='(true.in.ci = (Bxy.est.true >= Bxy.ci.lower.gmm) & (Bxy.est.true <= Bxy.ci.upper.gmm),
 133                                      zero.in.ci = (0 >= Bxy.ci.lower.gmm) & (0 <= Bxy.ci.upper.gmm),
 134                                      bias =  Bxy.est.gmm - Bxy.est.true,
 135                                      sign.correct = sign(Bxy.est.true) == sign(Bxy.est.gmm))]
 136
 137 x.gmm.plot <- x.gmm[,.(p.true.in.ci = mean(as.integer(true.in.ci)),
 138                                        mean.bias = mean(bias),
 139                                      p.sign.correct = mean(as.integer(sign.correct & (! zero.in.ci))),
 140                                      variable='x',
 141                                      method='2SLS+gmm'
 142                                      ),
 143                                     by=c('N','m')]
 144
 145 g.gmm <- df[,.(N,m,Bgy.est.true, Bgy.est.gmm,Bgy.ci.lower.gmm, Bgy.ci.upper.gmm)]
 146 g.gmm <- g.gmm[,':='(true.in.ci = (Bgy.est.true >= Bgy.ci.lower.gmm) & (Bgy.est.true <= Bgy.ci.upper.gmm),
 147                                      zero.in.ci = (0 >= Bgy.ci.lower.gmm) & (0 <= Bgy.ci.upper.gmm),
 148                                      bias =  Bgy.est.gmm - Bgy.est.true,
 149                                      sign.correct = sign(Bgy.est.true) == sign(Bgy.est.gmm))]
 150
 151 g.gmm.plot <- g.gmm[,.(p.true.in.ci = mean(as.integer(true.in.ci)),
 152                                        mean.bias = mean(bias),
 153                                      p.sign.correct = mean(as.integer(sign.correct & (! zero.in.ci))),
 154                                      variable='g',
 155                                      method='2SLS+gmm'
 156                                      ),
 157                                     by=c('N','m')]
 158
 159 plot.df <- rbindlist(list(x.naive.plot,g.naive.plot,x.amelia.full.plot,g.amelia.full.plot,x.amelia.nok.plot,g.amelia.nok.plot, x.mecor.plot, g.mecor.plot, x.gmm.plot, g.gmm.plot))
 160
 161 remember(plot.df,'example.1.plot.df')
 162
 163
 164 ggplot(plot.df,aes(y=N,x=m,color=p.sign.correct)) + geom_point() + facet_grid(variable ~ method) + scale_color_viridis_c(option='D') + theme_minimal() + xlab("Number of gold standard labels") + ylab("Total sample size")
 165
 166 ggplot(plot.df,aes(y=N,x=m,color=abs(mean.bias))) + geom_point() + facet_grid(variable ~ method) + scale_color_viridis_c(option='D') + theme_minimal() + xlab("Number of gold standard labels") + ylab("Total sample size")
 167
 168 ests <- df[,.(Bxy.est.true = mean(Bxy.est.true),
 169               Bxy.est.naive = mean(Bxy.est.naive),
 170               Bxy.est.feasible = mean(Bxy.est.feasible),
 171               Bxy.est.amelia.full = mean(Bxy.est.amelia.full),
 172               Bxy.est.amelia.nok = mean(Bxy.est.amelia.nok),
 173               Bxy.est.mecor = mean(Bxy.est.mecor),
 174               Bxy.est.gmm = mean(Bxy.est.gmm)),
 175            by=c("N","m")]