# amelia says use normal distribution for binary variables.
-
- amelia.out.k <- amelia(df, m=200, p2s=0, idvars=c('y','ystar','w'))
- mod.amelia.k <- zelig(y.obs~x+z, model='ls', data=amelia.out.k$imputations, cite=FALSE)
- (coefse <- combine_coef_se(mod.amelia.k, messages=FALSE))
- est.x.mi <- coefse['x','Estimate']
- est.x.se <- coefse['x','Std.Error']
- result <- append(result,
- list(Bxy.est.amelia.full = est.x.mi,
+ amelia_result <- list(Bxy.est.amelia.full = NA,
+ Bxy.ci.upper.amelia.full = NA,
+ Bxy.ci.lower.amelia.full = NA,
+ Bzy.est.amelia.full = NA,
+ Bzy.ci.upper.amelia.full = NA,
+ Bzy.ci.lower.amelia.full = NA
+ )
+
+ tryCatch({
+ amelia.out.k <- amelia(df, m=200, p2s=0, idvars=c('y','ystar','w'))
+ mod.amelia.k <- zelig(y.obs~x+z, model='ls', data=amelia.out.k$imputations, cite=FALSE)
+ (coefse <- combine_coef_se(mod.amelia.k, messages=FALSE))
+ est.x.mi <- coefse['x','Estimate']
+ est.x.se <- coefse['x','Std.Error']
+
+ est.z.mi <- coefse['z','Estimate']
+ est.z.se <- coefse['z','Std.Error']
+ amelia_result <- list(Bxy.est.amelia.full = est.x.mi,
Bxy.ci.upper.amelia.full = est.x.mi + 1.96 * est.x.se,
- Bxy.ci.lower.amelia.full = est.x.mi - 1.96 * est.x.se
- ))
-
- est.z.mi <- coefse['z','Estimate']
- est.z.se <- coefse['z','Std.Error']
-
- result <- append(result,
- list(Bzy.est.amelia.full = est.z.mi,
+ Bxy.ci.lower.amelia.full = est.x.mi - 1.96 * est.x.se,
+ Bzy.est.amelia.full = est.z.mi,
Bzy.ci.upper.amelia.full = est.z.mi + 1.96 * est.z.se,
Bzy.ci.lower.amelia.full = est.z.mi - 1.96 * est.z.se
- ))
+ )
+ },
+ error = function(e){
+ result[['error']] <- e}
+ )
+ result <- append(result,amelia_result)
return(result)
Bxy.ci.lower.naive = naive.ci.Bxy[1],
Bzy.ci.upper.naive = naive.ci.Bzy[2],
Bzy.ci.lower.naive = naive.ci.Bzy[1]))
-
+ amelia_result <- list(
+ Bxy.est.amelia.full = NULL,
+ Bxy.ci.upper.amelia.full = NULL,
+ Bxy.ci.lower.amelia.full = NULL,
+ Bzy.est.amelia.full = NULL,
+ Bzy.ci.upper.amelia.full = NULL,
+ Bzy.ci.lower.amelia.full = NULL
+ )
+
+ tryCatch({
+ amelia.out.k <- amelia(df, m=200, p2s=0, idvars=c('x','w'))
+ mod.amelia.k <- zelig(y~x.obs+z, model='ls', data=amelia.out.k$imputations, cite=FALSE)
+ (coefse <- combine_coef_se(mod.amelia.k))
+
+ est.x.mi <- coefse['x.obs','Estimate']
+ est.x.se <- coefse['x.obs','Std.Error']
+ est.z.mi <- coefse['z','Estimate']
+ est.z.se <- coefse['z','Std.Error']
+
+ amelia_result <- list(Bxy.est.amelia.full = est.x.mi,
+ Bxy.ci.upper.amelia.full = est.x.mi + 1.96 * est.x.se,
+ Bxy.ci.lower.amelia.full = est.x.mi - 1.96 * est.x.se,
+ Bzy.est.amelia.full = est.z.mi,
+ Bzy.ci.upper.amelia.full = est.z.mi + 1.96 * est.z.se,
+ Bzy.ci.lower.amelia.full = est.z.mi - 1.96 * est.z.se
+ )
- amelia.out.k <- amelia(df, m=200, p2s=0, idvars=c('x','w'))
- mod.amelia.k <- zelig(y~x.obs+z, model='ls', data=amelia.out.k$imputations, cite=FALSE)
- (coefse <- combine_coef_se(mod.amelia.k, messages=FALSE))
+ },
- est.x.mi <- coefse['x.obs','Estimate']
- est.x.se <- coefse['x.obs','Std.Error']
- result <- append(result,
- list(Bxy.est.amelia.full = est.x.mi,
- Bxy.ci.upper.amelia.full = est.x.mi + 1.96 * est.x.se,
- Bxy.ci.lower.amelia.full = est.x.mi - 1.96 * est.x.se
- ))
+ error = function(e){
+ result[['error']] <- e}
+ )
- est.z.mi <- coefse['z','Estimate']
- est.z.se <- coefse['z','Std.Error']
- result <- append(result,
- list(Bzy.est.amelia.full = est.z.mi,
- Bzy.ci.upper.amelia.full = est.z.mi + 1.96 * est.z.se,
- Bzy.ci.lower.amelia.full = est.z.mi - 1.96 * est.z.se
- ))
+ result <- append(result, amelia_result)
+
+ mle_result <- list(Bxy.est.mle = NULL,
+ Bxy.ci.upper.mle = NULL,
+ Bxy.ci.lower.mle = NULL,
+ Bzy.est.mle = NULL,
+ Bzy.ci.upper.mle = NULL,
+ Bzy.ci.lower.mle = NULL)
+ tryCatch({
temp.df <- copy(df)
temp.df <- temp.df[,x:=x.obs]
mod.caroll.lik <- measerr_mle(temp.df, outcome_formula=outcome_formula, proxy_formula=proxy_formula, truth_formula=truth_formula)
-
- ## tryCatch({
- ## mod.calibrated.mle <- mecor(y ~ MeasError(w_pred, reference = x.obs) + z, df, B=400, method='efficient')
- ## (mod.calibrated.mle)
- ## (mecor.ci <- summary(mod.calibrated.mle)$c$ci['x.obs',])
- ## result <- append(result, list(
- ## Bxy.est.mecor = mecor.ci['Estimate'],
- ## Bxy.ci.upper.mecor = mecor.ci['UCI'],
- ## Bxy.ci.lower.mecor = mecor.ci['LCI'])
- ## )
-
-
-
- fischer.info <- NA
- ci.upper <- NA
- ci.lower <- NA
-
- tryCatch({fischer.info <- solve(mod.caroll.lik$hessian)
+ fischer.info <- solve(mod.caroll.lik$hessian)
+ coef <- mod.caroll.lik$par
ci.upper <- coef + sqrt(diag(fischer.info)) * 1.96
ci.lower <- coef - sqrt(diag(fischer.info)) * 1.96
+ mle_result <- list(Bxy.est.mle = coef['x'],
+ Bxy.ci.upper.mle = ci.upper['x'],
+ Bxy.ci.lower.mle = ci.lower['x'],
+ Bzy.est.mle = coef['z'],
+ Bzy.ci.upper.mle = ci.upper['z'],
+ Bzy.ci.lower.mle = ci.lower['z'])
},
error=function(e) {result[['error']] <- as.character(e)
})
- coef <- mod.caroll.lik$par
- result <- append(result,
- list(Bxy.est.mle = coef['x'],
- Bxy.ci.upper.mle = ci.upper['x'],
- Bxy.ci.lower.mle = ci.lower['x'],
- Bzy.est.mle = coef['z'],
- Bzy.ci.upper.mle = ci.upper['z'],
- Bzy.ci.lower.mle = ci.lower['z']))
-
- mod.zhang.lik <- zhang.mle.iv(df)
- coef <- coef(mod.zhang.lik)
- ci <- confint(mod.zhang.lik,method='quad')
- result <- append(result,
- list(Bxy.est.zhang = coef['Bxy'],
- Bxy.ci.upper.zhang = ci['Bxy','97.5 %'],
- Bxy.ci.lower.zhang = ci['Bxy','2.5 %'],
- Bzy.est.zhang = coef['Bzy'],
- Bzy.ci.upper.zhang = ci['Bzy','97.5 %'],
- Bzy.ci.lower.zhang = ci['Bzy','2.5 %']))
+ result <- append(result, mle_result)
+
+ mod.zhang.lik <- zhang.mle.iv(df)
+ coef <- coef(mod.zhang.lik)
+ ci <- confint(mod.zhang.lik,method='quad')
+ result <- append(result,
+ list(Bxy.est.zhang = coef['Bxy'],
+ Bxy.ci.upper.zhang = ci['Bxy','97.5 %'],
+ Bxy.ci.lower.zhang = ci['Bxy','2.5 %'],
+ Bzy.est.zhang = coef['Bzy'],
+ Bzy.ci.upper.zhang = ci['Bzy','97.5 %'],
+ Bzy.ci.lower.zhang = ci['Bzy','2.5 %']))
## What if we can't observe k -- most realistic scenario. We can't include all the ML features in a model.
## amelia.out.nok <- amelia(df, m=200, p2s=0, idvars=c("x","w_pred"), noms=noms)
projects / ml_measurement_error_public.git / blobdiff