simulations/irr_simulation_base.R

   1 library(matrixStats) # for numerically stable logsumexps
   2
   3 options(amelia.parallel="no",
   4         amelia.ncpus=1)
   5 library(Amelia)
   6 source("measerr_methods.R")
   7 source("pl_methods.R")
   8
   9 run_simulation <- function(df, result, outcome_formula = y ~ x + z, proxy_formula = w_pred ~ x, coder_formulas=list(x.obs.1~x, x.obs.0~x), truth_formula = x ~ z){
  10
  11     accuracy <- df[,mean(w_pred==x)]
  12     result <- append(result, list(accuracy=accuracy))
  13
  14     (model.true <- lm(y ~ x + z, data=df))
  15     true.ci.Bxy <- confint(model.true)['x',]
  16     true.ci.Bzy <- confint(model.true)['z',]
  17
  18     result <- append(result, list(Bxy.est.true=coef(model.true)['x'],
  19                                   Bzy.est.true=coef(model.true)['z'],
  20                                   Bxy.ci.upper.true = true.ci.Bxy[2],
  21                                   Bxy.ci.lower.true = true.ci.Bxy[1],
  22                                   Bzy.ci.upper.true = true.ci.Bzy[2],
  23                                   Bzy.ci.lower.true = true.ci.Bzy[1]))
  24
  25
  26
  27
  28
  29     loa0.feasible <- lm(y ~ x.obs.0 + z, data = df[!(is.na(x.obs.1))])
  30
  31     loa0.ci.Bxy <- confint(loa0.feasible)['x.obs.0',]
  32     loa0.ci.Bzy <- confint(loa0.feasible)['z',]
  33
  34     result <- append(result, list(Bxy.est.loa0.feasible=coef(loa0.feasible)['x.obs.0'],
  35                                   Bzy.est.loa0.feasible=coef(loa0.feasible)['z'],
  36                                   Bxy.ci.upper.loa0.feasible = loa0.ci.Bxy[2],
  37                                   Bxy.ci.lower.loa0.feasible = loa0.ci.Bxy[1],
  38                                   Bzy.ci.upper.loa0.feasible = loa0.ci.Bzy[2],
  39                                   Bzy.ci.lower.loa0.feasible = loa0.ci.Bzy[1]))
  40     print("fitting loa0 model")
  41
  42     df.loa0.mle <- copy(df)
  43     df.loa0.mle[,x:=x.obs.0]
  44     loa0.mle <- measerr_mle(df.loa0.mle, outcome_formula=outcome_formula, proxy_formula=proxy_formula, truth_formula=truth_formula)
  45     fisher.info <- solve(loa0.mle$hessian)
  46     coef <- loa0.mle$par
  47     ci.upper <- coef + sqrt(diag(fisher.info)) * 1.96
  48     ci.lower <- coef - sqrt(diag(fisher.info)) * 1.96
  49
  50     result <- append(result, list(Bxy.est.loa0.mle=coef['x'],
  51                                   Bzy.est.loa0.mle=coef['z'],
  52                                   Bxy.ci.upper.loa0.mle = ci.upper['x'],
  53                                   Bxy.ci.lower.loa0.mle = ci.lower['x'],
  54                                   Bzy.ci.upper.loa0.mle = ci.upper['z'],
  55                                   Bzy.ci.lower.loa0.mle = ci.upper['z']))
  56
  57
  58
  59     loco.feasible <- lm(y ~ x.obs.1 + z, data = df[(!is.na(x.obs.1)) & (x.obs.1 == x.obs.0)])
  60
  61
  62     loco.feasible.ci.Bxy <- confint(loco.feasible)['x.obs.1',]
  63     loco.feasible.ci.Bzy <- confint(loco.feasible)['z',]
  64
  65     result <- append(result, list(Bxy.est.loco.feasible=coef(loco.feasible)['x.obs.1'],
  66                                   Bzy.est.loco.feasible=coef(loco.feasible)['z'],
  67                                   Bxy.ci.upper.loco.feasible = loco.feasible.ci.Bxy[2],
  68                                   Bxy.ci.lower.loco.feasible = loco.feasible.ci.Bxy[1],
  69                                   Bzy.ci.upper.loco.feasible = loco.feasible.ci.Bzy[2],
  70                                   Bzy.ci.lower.loco.feasible = loco.feasible.ci.Bzy[1]))
  71
  72
  73     (model.naive <- lm(y~w_pred+z, data=df))
  74
  75     naive.ci.Bxy <- confint(model.naive)['w_pred',]
  76     naive.ci.Bzy <- confint(model.naive)['z',]
  77
  78     result <- append(result, list(Bxy.est.naive=coef(model.naive)['w_pred'],
  79                                   Bzy.est.naive=coef(model.naive)['z'],
  80                                   Bxy.ci.upper.naive = naive.ci.Bxy[2],
  81                                   Bxy.ci.lower.naive = naive.ci.Bxy[1],
  82                                   Bzy.ci.upper.naive = naive.ci.Bzy[2],
  83                                   Bzy.ci.lower.naive = naive.ci.Bzy[1]))
  84
  85     print("fitting loco model")
  86
  87     df.loco.mle <- copy(df)
  88     df.loco.mle[,x.obs:=NA]
  89     df.loco.mle[(x.obs.0)==(x.obs.1),x.obs:=x.obs.0]
  90     df.loco.mle[,x.true:=x]
  91     df.loco.mle[,x:=x.obs]
  92     print(df.loco.mle[!is.na(x.obs.1),mean(x.true==x,na.rm=TRUE)])
  93     loco.accuracy <- df.loco.mle[(!is.na(x.obs.1)) & (x.obs.1 == x.obs.0),mean(x.obs.1 == x.true)]
  94     loco.mle <- measerr_mle(df.loco.mle, outcome_formula=outcome_formula, proxy_formula=proxy_formula, truth_formula=truth_formula)
  95     fisher.info <- solve(loco.mle$hessian)
  96     coef <- loco.mle$par
  97     ci.upper <- coef + sqrt(diag(fisher.info)) * 1.96
  98     ci.lower <- coef - sqrt(diag(fisher.info)) * 1.96
  99
 100     result <- append(result, list(loco.accuracy=loco.accuracy,
 101                                   Bxy.est.loco.mle=coef['x'],
 102                                   Bzy.est.loco.mle=coef['z'],
 103                                   Bxy.ci.upper.loco.mle = ci.upper['x'],
 104                                   Bxy.ci.lower.loco.mle = ci.lower['x'],
 105                                   Bzy.ci.upper.loco.mle = ci.upper['z'],
 106                                   Bzy.ci.lower.loco.mle = ci.lower['z']))
 107
 108     df.double.proxy.mle <- copy(df)
 109     df.double.proxy.mle[,x.obs:=NA]
 110     print("fitting double proxy model")
 111
 112     double.proxy.mle <- measerr_irr_mle(df.double.proxy.mle, outcome_formula=outcome_formula, proxy_formula=proxy_formula, coder_formulas=coder_formulas[1], truth_formula=truth_formula)
 113     fisher.info <- solve(double.proxy.mle$hessian)
 114     coef <- double.proxy.mle$par
 115     ci.upper <- coef + sqrt(diag(fisher.info)) * 1.96
 116     ci.lower <- coef - sqrt(diag(fisher.info)) * 1.96
 117
 118     result <- append(result, list(
 119                                   Bxy.est.double.proxy=coef['x'],
 120                                   Bzy.est.double.proxy=coef['z'],
 121                                   Bxy.ci.upper.double.proxy = ci.upper['x'],
 122                                   Bxy.ci.lower.double.proxy = ci.lower['x'],
 123                                   Bzy.ci.upper.double.proxy = ci.upper['z'],
 124                                   Bzy.ci.lower.double.proxy = ci.lower['z']))
 125
 126     df.triple.proxy.mle <- copy(df)
 127     df.triple.proxy.mle[,x.obs:=NA]
 128
 129     print("fitting triple proxy model")
 130     triple.proxy.mle <- measerr_irr_mle(df.triple.proxy.mle, outcome_formula=outcome_formula, proxy_formula=proxy_formula, coder_formulas=coder_formulas, truth_formula=truth_formula)
 131     fisher.info <- solve(triple.proxy.mle$hessian)
 132     coef <- triple.proxy.mle$par
 133     ci.upper <- coef + sqrt(diag(fisher.info)) * 1.96
 134     ci.lower <- coef - sqrt(diag(fisher.info)) * 1.96
 135
 136     result <- append(result, list(
 137                                   Bxy.est.triple.proxy=coef['x'],
 138                                   Bzy.est.triple.proxy=coef['z'],
 139                                   Bxy.ci.upper.triple.proxy = ci.upper['x'],
 140                                   Bxy.ci.lower.triple.proxy = ci.lower['x'],
 141                                   Bzy.ci.upper.triple.proxy = ci.upper['z'],
 142                                   Bzy.ci.lower.triple.proxy = ci.lower['z']))
 143     tryCatch({
 144     amelia.out.k <- amelia(df.loco.mle, m=200, p2s=0, idvars=c('x.true','w','x.obs.1','x.obs.0','x'))
 145     mod.amelia.k <- zelig(y~x.obs+z, model='ls', data=amelia.out.k$imputations, cite=FALSE)
 146     (coefse <- combine_coef_se(mod.amelia.k, messages=FALSE))
 147
 148     est.x.mi <- coefse['x.obs','Estimate']
 149     est.x.se <- coefse['x.obs','Std.Error']
 150     result <- append(result,
 151                      list(Bxy.est.amelia.full = est.x.mi,
 152                           Bxy.ci.upper.amelia.full = est.x.mi + 1.96 * est.x.se,
 153                           Bxy.ci.lower.amelia.full = est.x.mi - 1.96 * est.x.se
 154                           ))
 155
 156     est.z.mi <- coefse['z','Estimate']
 157     est.z.se <- coefse['z','Std.Error']
 158
 159     result <- append(result,
 160                      list(Bzy.est.amelia.full = est.z.mi,
 161                           Bzy.ci.upper.amelia.full = est.z.mi + 1.96 * est.z.se,
 162                           Bzy.ci.lower.amelia.full = est.z.mi - 1.96 * est.z.se
 163                           ))
 164
 165     },
 166     error = function(e){
 167         message("An error occurred:\n",e)
 168         result$error <-paste0(result$error,'\n', e)
 169     }
 170     )
 171
 172     tryCatch({
 173
 174         mod.zhang.lik <- zhang.mle.iv(df.loco.mle)
 175         coef <- coef(mod.zhang.lik)
 176         ci <- confint(mod.zhang.lik,method='quad')
 177         result <- append(result,
 178                          list(Bxy.est.zhang = coef['Bxy'],
 179                               Bxy.ci.upper.zhang = ci['Bxy','97.5 %'],
 180                               Bxy.ci.lower.zhang = ci['Bxy','2.5 %'],
 181                               Bzy.est.zhang = coef['Bzy'],
 182                               Bzy.ci.upper.zhang = ci['Bzy','97.5 %'],
 183                               Bzy.ci.lower.zhang = ci['Bzy','2.5 %']))
 184     },
 185
 186     error = function(e){
 187         message("An error occurred:\n",e)
 188         result$error <- paste0(result$error,'\n', e)
 189     })
 190
 191     df <- df.loco.mle
 192     N <- nrow(df)
 193     m <- nrow(df[!is.na(x.obs)])
 194     p <- v <- train <- rep(0,N)
 195     M <- m
 196     p[(M+1):(N)] <- 1
 197     v[1:(M)] <- 1
 198     df <- df[order(x.obs)]
 199     y <- df[,y]
 200     x <- df[,x.obs]
 201     z <- df[,z]
 202     w <- df[,w_pred]
 203     # gmm gets pretty close
 204     (gmm.res <- predicted_covariates(y, x, z, w, v, train, p, max_iter=100, verbose=TRUE))
 205
 206     result <- append(result,
 207                      list(Bxy.est.gmm = gmm.res$beta[1,1],
 208                           Bxy.ci.upper.gmm = gmm.res$confint[1,2],
 209                           Bxy.ci.lower.gmm = gmm.res$confint[1,1],
 210                           gmm.ER_pval = gmm.res$ER_pval
 211                           ))
 212
 213     result <- append(result,
 214                      list(Bzy.est.gmm = gmm.res$beta[2,1],
 215                           Bzy.ci.upper.gmm = gmm.res$confint[2,2],
 216                           Bzy.ci.lower.gmm = gmm.res$confint[2,1]))
 217
 218
 219
 220     return(result)
 221
 222 }