From c1dbbfd0dd88defca0ce00425910757e436284ad Mon Sep 17 00:00:00 2001 From: Nathan TeBlunthuis Date: Tue, 28 Feb 2023 16:28:35 -0800 Subject: [PATCH] update simulation base from hyak --- simulations/simulation_base.R | 45 +++++++++++++++++++++++++---------- 1 file changed, 33 insertions(+), 12 deletions(-) diff --git a/simulations/simulation_base.R b/simulations/simulation_base.R index 73544e9..bafd7d3 100644 --- a/simulations/simulation_base.R +++ b/simulations/simulation_base.R @@ -89,7 +89,7 @@ my.mle <- function(df){ 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)]) 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[,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) - coef <- mod.caroll.lik$par - ci.upper <- coef + sqrt(diag(fischer.info)) * 1.96 - ci.lower <- coef - sqrt(diag(fischer.info)) * 1.96 + + if(confint_method=='quad'){ + mod.caroll.lik <- measerr_mle_dv(temp.df, outcome_formula=outcome_formula, proxy_formula=proxy_formula) + 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 + } + 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, list(Bxy.est.mle = coef['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 -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.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({ 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) - 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 + if(confint_method=='quad'){ + mod.caroll.lik <- measerr_mle(temp.df, outcome_formula=outcome_formula, proxy_formula=proxy_formula, truth_formula=truth_formula, method='optim') + 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 + } 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'], Bxy.ci.upper.mle = ci.upper['x'], Bxy.ci.lower.mle = ci.lower['x'], -- 2.39.2