groceryheist/ml_measurement_error_public
1
0
Fork 0
Code

update simulation base from hyak

Browse Source
This commit is contained in:
Nathan TeBlunthuis 2023-02-28 16:28:35 -08:00
parent 69948cae1e
commit c1dbbfd0dd
1 changed files with 33 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

45
simulations/simulation_base.R
View File

@ -89,7 +89,7 @@ my.mle <- function(df){
return(mlefit) return(mlefit)
} }
run_simulation_depvar <- function(df, result, outcome_formula=y~x+z, proxy_formula=w_pred~y){ run_simulation_depvar <- function(df, result, outcome_formula=y~x+z, proxy_formula=w_pred~y, confint_method='quad'){
(accuracy <- df[,mean(w_pred==y)]) (accuracy <- df[,mean(w_pred==y)])
result <- append(result, list(accuracy=accuracy)) result <- append(result, list(accuracy=accuracy))
@ -150,11 +150,23 @@ run_simulation_depvar <- function(df, result, outcome_formula=y~x+z, proxy_formu
temp.df <- copy(df) temp.df <- copy(df)
temp.df[,y:=y.obs] temp.df[,y:=y.obs]
mod.caroll.lik <- measerr_mle_dv(temp.df, outcome_formula=outcome_formula, proxy_formula=proxy_formula)
fischer.info <- solve(mod.caroll.lik$hessian) if(confint_method=='quad'){
coef <- mod.caroll.lik$par mod.caroll.lik <- measerr_mle_dv(temp.df, outcome_formula=outcome_formula, proxy_formula=proxy_formula)
ci.upper <- coef + sqrt(diag(fischer.info)) * 1.96 fischer.info <- solve(mod.caroll.lik$hessian)
ci.lower <- coef - sqrt(diag(fischer.info)) * 1.96 coef <- mod.caroll.lik$par
ci.upper <- coef + sqrt(diag(fischer.info)) * 1.96
ci.lower <- coef - sqrt(diag(fischer.info)) * 1.96
}
else{ ## confint_method is 'profile'
mod.caroll.lik <- measerr_mle_dv(temp.df, outcome_formula=outcome_formula, proxy_formula=proxy_formula, method='bbmle')
coef <- coef(mod.caroll.lik)
ci <- confint(mod.caroll.lik, method='spline')
ci.lower <- ci[,'2.5 %']
ci.upper <- ci[,'97.5 %']
}
result <- append(result, result <- append(result,
list(Bxy.est.mle = coef['x'], list(Bxy.est.mle = coef['x'],
Bxy.ci.upper.mle = ci.upper['x'], Bxy.ci.upper.mle = ci.upper['x'],
@ -216,7 +228,7 @@ run_simulation_depvar <- function(df, result, outcome_formula=y~x+z, proxy_formu
## outcome_formula, proxy_formula, and truth_formula are passed to measerr_mle ## outcome_formula, proxy_formula, and truth_formula are passed to measerr_mle
run_simulation <- function(df, result, outcome_formula=y~x+z, proxy_formula=NULL, truth_formula=NULL){ run_simulation <- function(df, result, outcome_formula=y~x+z, proxy_formula=NULL, truth_formula=NULL, confint_method='quad'){
accuracy <- df[,mean(w_pred==x)] accuracy <- df[,mean(w_pred==x)]
accuracy.y0 <- df[y<=0,mean(w_pred==x)] accuracy.y0 <- df[y<=0,mean(w_pred==x)]
@ -328,11 +340,20 @@ run_simulation <- function(df, result, outcome_formula=y~x+z, proxy_formula=NUL
tryCatch({ tryCatch({
temp.df <- copy(df) temp.df <- copy(df)
temp.df <- temp.df[,x:=x.obs] 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) if(confint_method=='quad'){
fischer.info <- solve(mod.caroll.lik$hessian) mod.caroll.lik <- measerr_mle(temp.df, outcome_formula=outcome_formula, proxy_formula=proxy_formula, truth_formula=truth_formula, method='optim')
coef <- mod.caroll.lik$par fischer.info <- solve(mod.caroll.lik$hessian)
ci.upper <- coef + sqrt(diag(fischer.info)) * 1.96 coef <- mod.caroll.lik$par
ci.lower <- coef - sqrt(diag(fischer.info)) * 1.96 ci.upper <- coef + sqrt(diag(fischer.info)) * 1.96
ci.lower <- coef - sqrt(diag(fischer.info)) * 1.96
} else { # confint_method == 'bbmle'
mod.caroll.lik <- measerr_mle(temp.df, outcome_formula=outcome_formula, proxy_formula=proxy_formula, truth_formula=truth_formula, method='bbmle')
coef <- coef(mod.caroll.lik)
ci <- confint(mod.caroll.lik, method='spline')
ci.lower <- ci[,'2.5 %']
ci.upper <- ci[,'97.5 %']
}
mle_result <- list(Bxy.est.mle = coef['x'], mle_result <- list(Bxy.est.mle = coef['x'],
Bxy.ci.upper.mle = ci.upper['x'], Bxy.ci.upper.mle = ci.upper['x'],
Bxy.ci.lower.mle = ci.lower['x'], Bxy.ci.lower.mle = ci.lower['x'],