X-Git-Url: https://code.communitydata.science/ml_measurement_error_public.git/blobdiff_plain/46e2d1fe4876a9ed906b723f9e5f74fcc949e339..d0c5766bdf867a81a2477d2cac1d40812110af90:/simulations/simulation_base.R?ds=inline diff --git a/simulations/simulation_base.R b/simulations/simulation_base.R index 0f03276..27f0276 100644 --- a/simulations/simulation_base.R +++ b/simulations/simulation_base.R @@ -41,21 +41,26 @@ my.pseudo.mle <- function(df){ ## Zhang got this model from Hausman 1998 ### I think this is actually eqivalent to the pseudo.mle method zhang.mle.iv <- function(df){ - nll <- function(B0=0, Bxy=0, Bzy=0, sigma_y=0.1, ppv=0.9, npv=0.9){ df.obs <- df[!is.na(x.obs)] df.unobs <- df[is.na(x.obs)] + tn <- df.obs[(w_pred == 0) & (x.obs == w_pred),.N] + pn <- df.obs[(w_pred==0), .N] + npv <- tn / pn + + tp <- df.obs[(w_pred==1) & (x.obs == w_pred),.N] + pp <- df.obs[(w_pred==1),.N] + ppv <- tp / pp + + nll <- function(B0=0, Bxy=0, Bzy=0, sigma_y=0.1){ + ## fpr = 1 - TNR ### Problem: accounting for uncertainty in ppv / npv - - ll.w1x1.obs <- with(df.obs[(w_pred==1)], dbinom(x.obs,size=1,prob=ppv,log=T)) - ll.w0x0.obs <- with(df.obs[(w_pred==0)], dbinom(1-x.obs,size=1,prob=npv,log=T)) ## fnr = 1 - TPR ll.y.obs <- with(df.obs, dnorm(y, B0 + Bxy * x + Bzy * z, sd=sigma_y,log=T)) ll <- sum(ll.y.obs) - ll <- ll + sum(ll.w1x1.obs) + sum(ll.w0x0.obs) - + # unobserved case; integrate out x ll.x.1 <- with(df.unobs, dnorm(y, B0 + Bxy + Bzy * z, sd = sigma_y, log=T)) ll.x.0 <- with(df.unobs, dnorm(y, B0 + Bzy * z, sd = sigma_y,log=T)) @@ -66,55 +71,90 @@ zhang.mle.iv <- function(df){ ## case x == 0 lls.x.0 <- colLogSumExps(rbind(log(1-npv) + ll.x.1, log(npv) + ll.x.0)) - lls <- colLogSumExps(rbind(lls.x.1, lls.x.0)) + lls <- colLogSumExps(rbind(df.unobs$w_pred * lls.x.1, (1-df.unobs$w_pred) * lls.x.0)) ll <- ll + sum(lls) return(-ll) } - mlefit <- mle2(minuslogl = nll, control=list(maxit=1e6), lower=list(sigma_y=0.0001, B0=-Inf, Bxy=-Inf, Bzy=-Inf,ppv=0.00001, npv=0.00001), - upper=list(sigma_y=Inf, B0=Inf, Bxy=Inf, Bzy=Inf, ppv=0.99999,npv=0.99999),method='L-BFGS-B') + mlefit <- mle2(minuslogl = nll, control=list(maxit=1e6), lower=list(sigma_y=0.0001, B0=-Inf, Bxy=-Inf, Bzy=-Inf), + upper=list(sigma_y=Inf, B0=Inf, Bxy=Inf, Bzy=Inf),method='L-BFGS-B') return(mlefit) } -## this is equivalent to the pseudo-liklihood model from Carolla -zhang.mle.dv <- function(df){ +## this is equivalent to the pseudo-liklihood model from Caroll +## zhang.mle.dv <- function(df){ - nll <- function(B0=0, Bxy=0, Bzy=0, ppv=0.9, npv=0.9){ - df.obs <- df[!is.na(y.obs)] +## nll <- function(B0=0, Bxy=0, Bzy=0, ppv=0.9, npv=0.9){ +## df.obs <- df[!is.na(y.obs)] - ## fpr = 1 - TNR - ll.w0y0 <- with(df.obs[y.obs==0],dbinom(1-w_pred,1,npv,log=TRUE)) - ll.w1y1 <- with(df.obs[y.obs==1],dbinom(w_pred,1,ppv,log=TRUE)) - - # observed case - ll.y.obs <- vector(mode='numeric', length=nrow(df.obs)) - ll.y.obs[df.obs$y.obs==1] <- with(df.obs[y.obs==1], plogis(B0 + Bxy * x + Bzy * z,log=T)) - ll.y.obs[df.obs$y.obs==0] <- with(df.obs[y.obs==0], plogis(B0 + Bxy * x + Bzy * z,log=T,lower.tail=FALSE)) - - ll <- sum(ll.y.obs) + sum(ll.w0y0) + sum(ll.w1y1) - - # unobserved case; integrate out y - ## case y = 1 - ll.y.1 <- vector(mode='numeric', length=nrow(df)) - pi.y.1 <- with(df,plogis(B0 + Bxy * x + Bzy*z, log=T)) - ## P(w=1| y=1)P(y=1) + P(w=0|y=1)P(y=1) = P(w=1,y=1) + P(w=0,y=1) - lls.y.1 <- colLogSumExps(rbind(log(ppv) + pi.y.1, log(1-ppv) + pi.y.1)) +## ## fpr = 1 - TNR +## ll.w0y0 <- with(df.obs[y.obs==0],dbinom(1-w_pred,1,npv,log=TRUE)) +## ll.w1y1 <- with(df.obs[y.obs==1],dbinom(w_pred,1,ppv,log=TRUE)) + +## # observed case +## ll.y.obs <- vector(mode='numeric', length=nrow(df.obs)) +## ll.y.obs[df.obs$y.obs==1] <- with(df.obs[y.obs==1], plogis(B0 + Bxy * x + Bzy * z,log=T)) +## ll.y.obs[df.obs$y.obs==0] <- with(df.obs[y.obs==0], plogis(B0 + Bxy * x + Bzy * z,log=T,lower.tail=FALSE)) + +## ll <- sum(ll.y.obs) + sum(ll.w0y0) + sum(ll.w1y1) + +## # unobserved case; integrate out y +## ## case y = 1 +## ll.y.1 <- vector(mode='numeric', length=nrow(df)) +## pi.y.1 <- with(df,plogis(B0 + Bxy * x + Bzy*z, log=T)) +## ## P(w=1| y=1)P(y=1) + P(w=0|y=1)P(y=1) = P(w=1,y=1) + P(w=0,y=1) +## lls.y.1 <- colLogSumExps(rbind(log(ppv) + pi.y.1, log(1-ppv) + pi.y.1)) - ## case y = 0 - ll.y.0 <- vector(mode='numeric', length=nrow(df)) - pi.y.0 <- with(df,plogis(B0 + Bxy * x + Bzy*z, log=T,lower.tail=FALSE)) +## ## case y = 0 +## ll.y.0 <- vector(mode='numeric', length=nrow(df)) +## pi.y.0 <- with(df,plogis(B0 + Bxy * x + Bzy*z, log=T,lower.tail=FALSE)) + +## ## P(w=1 | y=0)P(y=0) + P(w=0|y=0)P(y=0) = P(w=1,y=0) + P(w=0,y=0) +## lls.y.0 <- colLogSumExps(rbind(log(npv) + pi.y.0, log(1-npv) + pi.y.0)) + +## lls <- colLogSumExps(rbind(lls.y.1, lls.y.0)) +## ll <- ll + sum(lls) +## return(-ll) +## } +## mlefit <- mle2(minuslogl = nll, control=list(maxit=1e6),method='L-BFGS-B',lower=list(B0=-Inf, Bxy=-Inf, Bzy=-Inf, ppv=0.001,npv=0.001), +## upper=list(B0=Inf, Bxy=Inf, Bzy=Inf,ppv=0.999,npv=0.999)) +## return(mlefit) +## } - ## P(w=1 | y=0)P(y=0) + P(w=0|y=0)P(y=0) = P(w=1,y=0) + P(w=0,y=0) - lls.y.0 <- colLogSumExps(rbind(log(npv) + pi.y.0, log(1-npv) + pi.y.0)) +zhang.mle.dv <- function(df){ + df.obs <- df[!is.na(y.obs)] + df.unobs <- df[is.na(y.obs)] - lls <- colLogSumExps(rbind(lls.y.1, lls.y.0)) - ll <- ll + sum(lls) - return(-ll) + fp <- df.obs[(w_pred==1) & (y.obs != w_pred),.N] + p <- df.obs[(w_pred==1),.N] + fpr <- fp / p + fn <- df.obs[(w_pred==0) & (y.obs != w_pred), .N] + n <- df.obs[(w_pred==0),.N] + fnr <- fn / n + + nll <- function(B0=0, Bxy=0, Bzy=0){ + + + ## observed case + ll.y.obs <- vector(mode='numeric', length=nrow(df.obs)) + ll.y.obs[df.obs$y.obs==1] <- with(df.obs[y.obs==1], plogis(B0 + Bxy * x + Bzy * z,log=T)) + ll.y.obs[df.obs$y.obs==0] <- with(df.obs[y.obs==0], plogis(B0 + Bxy * x + Bzy * z,log=T,lower.tail=FALSE)) + + ll <- sum(ll.y.obs) + + pi.y.1 <- with(df,plogis(B0 + Bxy * x + Bzy*z, log=T)) + pi.y.0 <- with(df,plogis(B0 + Bxy * x + Bzy*z, log=T,lower.tail=FALSE)) + + lls <- with(df.unobs, colLogSumExps(rbind(w_pred * colLogSumExps(rbind(log(fpr), log(1 - fnr - fpr)+pi.y.1)), + (1-w_pred) * colLogSumExps(rbind(log(1-fpr), log(1 - fnr - fpr)+pi.y.0))))) + + ll <- ll + sum(lls) + return(-ll) } - mlefit <- mle2(minuslogl = nll, control=list(maxit=1e6),method='L-BFGS-B',lower=list(B0=-Inf, Bxy=-Inf, Bzy=-Inf, ppv=0.001,npv=0.001), - upper=list(B0=Inf, Bxy=Inf, Bzy=Inf,ppv=0.999,npv=0.999)) + mlefit <- mle2(minuslogl = nll, control=list(maxit=1e6),method='L-BFGS-B',lower=c(B0=-Inf, Bxy=-Inf, Bzy=-Inf), + upper=c(B0=Inf, Bxy=Inf, Bzy=Inf)) return(mlefit) } - + ## This uses the likelihood approach from Carroll page 353. ## assumes that we have a good measurement error model my.mle <- function(df){ @@ -170,11 +210,19 @@ run_simulation_depvar <- function(df, result, outcome_formula=y~x+z, proxy_formu accuracy <- df[,mean(w_pred==y)] result <- append(result, list(accuracy=accuracy)) + error.cor.x <- cor(df$x, df$w - df$x) + result <- append(result, list(error.cor.x = error.cor.x)) + model.null <- glm(y~1, data=df,family=binomial(link='logit')) (model.true <- glm(y ~ x + z, data=df,family=binomial(link='logit'))) + (lik.ratio <- exp(logLik(model.true) - logLik(model.null))) + true.ci.Bxy <- confint(model.true)['x',] true.ci.Bzy <- confint(model.true)['z',] + + result <- append(result, list(lik.ratio=lik.ratio)) + result <- append(result, list(Bxy.est.true=coef(model.true)['x'], Bzy.est.true=coef(model.true)['z'], Bxy.ci.upper.true = true.ci.Bxy[2], @@ -211,7 +259,7 @@ run_simulation_depvar <- function(df, result, outcome_formula=y~x+z, proxy_formu naivecont.ci.Bxy <- confint(model.naive.cont)['x',] naivecont.ci.Bzy <- confint(model.naive.cont)['z',] - ## my implementatoin of liklihood based correction + ## my implementation of liklihood based correction temp.df <- copy(df) temp.df[,y:=y.obs] @@ -241,7 +289,8 @@ run_simulation_depvar <- function(df, result, outcome_formula=y~x+z, proxy_formu Bzy.est.zhang = coef['Bzy'], Bzy.ci.upper.zhang = ci['Bzy','97.5 %'], Bzy.ci.lower.zhang = ci['Bzy','2.5 %'])) - + + # amelia says use normal distribution for binary variables. tryCatch({ @@ -278,11 +327,36 @@ 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=w_pred~x, truth_formula=x~z){ +run_simulation <- function(df, result, outcome_formula=y~x+z, proxy_formula=NULL, truth_formula=NULL){ accuracy <- df[,mean(w_pred==x)] - result <- append(result, list(accuracy=accuracy)) - + accuracy.y0 <- df[y<=0,mean(w_pred==x)] + accuracy.y1 <- df[y>=0,mean(w_pred==x)] + cor.y.xi <- cor(df$x - df$w_pred, df$y) + + fnr <- df[w_pred==0,mean(w_pred!=x)] + fnr.y0 <- df[(w_pred==0) & (y<=0),mean(w_pred!=x)] + fnr.y1 <- df[(w_pred==0) & (y>=0),mean(w_pred!=x)] + + fpr <- df[w_pred==1,mean(w_pred!=x)] + fpr.y0 <- df[(w_pred==1) & (y<=0),mean(w_pred!=x)] + fpr.y1 <- df[(w_pred==1) & (y>=0),mean(w_pred!=x)] + cor.resid.w_pred <- cor(resid(lm(y~x+z,df)),df$w_pred) + + result <- append(result, list(accuracy=accuracy, + accuracy.y0=accuracy.y0, + accuracy.y1=accuracy.y1, + cor.y.xi=cor.y.xi, + fnr=fnr, + fnr.y0=fnr.y0, + fnr.y1=fnr.y1, + fpr=fpr, + fpr.y0=fpr.y0, + fpr.y1=fpr.y1, + cor.resid.w_pred=cor.resid.w_pred + )) + + result <- append(result, list(cor.xz=cor(df$x,df$z))) (model.true <- lm(y ~ x + z, data=df)) true.ci.Bxy <- confint(model.true)['x',] true.ci.Bzy <- confint(model.true)['z',] @@ -320,7 +394,7 @@ run_simulation <- function(df, result, outcome_formula=y~x+z, proxy_formula=w_p tryCatch({ - amelia.out.k <- amelia(df, m=200, p2s=0, idvars=c('x','w_pred')) + 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))