X-Git-Url: https://code.communitydata.science/ml_measurement_error_public.git/blobdiff_plain/47e9367ed5c61b721bdc17cddd76bced4f8ed621..979dc14b6861ae31f00d56392fd5b8cf69f17333:/simulations/measerr_methods.R diff --git a/simulations/measerr_methods.R b/simulations/measerr_methods.R index 6bf8c3f..00f1746 100644 --- a/simulations/measerr_methods.R +++ b/simulations/measerr_methods.R @@ -102,17 +102,211 @@ measerr_mle_dv <- function(df, outcome_formula, outcome_family=binomial(link='lo return(fit) } +## Experimental, and not necessary if errors are independent. +measerr_irr_mle <- function(df, outcome_formula, outcome_family=gaussian(), rater_formula, proxy_formula, proxy_family=binomial(link='logit'), truth_formula, truth_family=binomial(link='logit')){ + + ### in this scenario, the ground truth also has measurement error, but we have repeated measures for it. + + ## probability of y given observed data. + df.obs <- df[!is.na(x.obs.1)] + proxy.variable <- all.vars(proxy_formula)[1] + df.x.obs.1 <- copy(df.obs)[,x:=1] + df.x.obs.0 <- copy(df.obs)[,x:=0] + y.obs <- df.obs[,y] + + nll <- function(params){ + outcome.model.matrix.x.obs.0 <- model.matrix(outcome_formula, df.x.obs.0) + outcome.model.matrix.x.obs.1 <- model.matrix(outcome_formula, df.x.obs.1) + + param.idx <- 1 + n.outcome.model.covars <- dim(outcome.model.matrix.x.obs.0)[2] + outcome.params <- params[param.idx:n.outcome.model.covars] + param.idx <- param.idx + n.outcome.model.covars + + sigma.y <- params[param.idx] + param.idx <- param.idx + 1 + + ll.y.x.obs.0 <- dnorm(y.obs, outcome.params %*% t(outcome.model.matrix.x.obs.0),sd=sigma.y, log=TRUE) + ll.y.x.obs.1 <- dnorm(y.obs, outcome.params %*% t(outcome.model.matrix.x.obs.1),sd=sigma.y, log=TRUE) + + ## assume that the two coders are statistically independent conditional on x + ll.x.obs.0.x0 <- vector(mode='numeric', length=nrow(df.obs)) + ll.x.obs.1.x0 <- vector(mode='numeric', length=nrow(df.obs)) + ll.x.obs.0.x1 <- vector(mode='numeric', length=nrow(df.obs)) + ll.x.obs.1.x1 <- vector(mode='numeric', length=nrow(df.obs)) + + rater.model.matrix.x.obs.0 <- model.matrix(rater_formula, df.x.obs.0) + rater.model.matrix.x.obs.1 <- model.matrix(rater_formula, df.x.obs.1) + + n.rater.model.covars <- dim(rater.model.matrix.x.obs.0)[2] + rater.0.params <- params[param.idx:(n.rater.model.covars + param.idx - 1)] + param.idx <- param.idx + n.rater.model.covars + + rater.1.params <- params[param.idx:(n.rater.model.covars + param.idx - 1)] + param.idx <- param.idx + n.rater.model.covars + + # probability of rater 0 if x is 0 or 1 + ll.x.obs.0.x0[df.obs$x.obs.0==1] <- plogis(rater.0.params %*% t(rater.model.matrix.x.obs.0[df.obs$x.obs.0==1,]), log=TRUE) + ll.x.obs.0.x0[df.obs$x.obs.0==0] <- plogis(rater.0.params %*% t(rater.model.matrix.x.obs.0[df.obs$x.obs.0==0,]), log=TRUE, lower.tail=FALSE) + ll.x.obs.0.x1[df.obs$x.obs.0==1] <- plogis(rater.0.params %*% t(rater.model.matrix.x.obs.1[df.obs$x.obs.0==1,]), log=TRUE) + ll.x.obs.0.x1[df.obs$x.obs.0==0] <- plogis(rater.0.params %*% t(rater.model.matrix.x.obs.1[df.obs$x.obs.0==0,]), log=TRUE, lower.tail=FALSE) + + # probability of rater 1 if x is 0 or 1 + ll.x.obs.1.x0[df.obs$x.obs.1==1] <- plogis(rater.1.params %*% t(rater.model.matrix.x.obs.0[df.obs$x.obs.1==1,]), log=TRUE) + ll.x.obs.1.x0[df.obs$x.obs.1==0] <- plogis(rater.1.params %*% t(rater.model.matrix.x.obs.0[df.obs$x.obs.1==0,]), log=TRUE, lower.tail=FALSE) + ll.x.obs.1.x1[df.obs$x.obs.1==1] <- plogis(rater.1.params %*% t(rater.model.matrix.x.obs.1[df.obs$x.obs.1==1,]), log=TRUE) + ll.x.obs.1.x1[df.obs$x.obs.1==0] <- plogis(rater.1.params %*% t(rater.model.matrix.x.obs.1[df.obs$x.obs.1==0,]), log=TRUE, lower.tail=FALSE) + + proxy.model.matrix.x0 <- model.matrix(proxy_formula, df.x.obs.0) + proxy.model.matrix.x1 <- model.matrix(proxy_formula, df.x.obs.1) + + n.proxy.model.covars <- dim(proxy.model.matrix.x0)[2] + proxy.params <- params[param.idx:(n.proxy.model.covars+param.idx-1)] + param.idx <- param.idx + n.proxy.model.covars + + proxy.obs <- with(df.obs, eval(parse(text=proxy.variable))) + + if( (proxy_family$family=="binomial") & (proxy_family$link=='logit')){ + ll.w.obs.x0 <- vector(mode='numeric',length=dim(proxy.model.matrix.x0)[1]) + ll.w.obs.x1 <- vector(mode='numeric',length=dim(proxy.model.matrix.x1)[1]) + + # proxy_formula likelihood using logistic regression + ll.w.obs.x0[proxy.obs==1] <- plogis(proxy.params %*% t(proxy.model.matrix.x0[proxy.obs==1,]),log=TRUE) + ll.w.obs.x0[proxy.obs==0] <- plogis(proxy.params %*% t(proxy.model.matrix.x0[proxy.obs==0,]),log=TRUE, lower.tail=FALSE) + + ll.w.obs.x1[proxy.obs==1] <- plogis(proxy.params %*% t(proxy.model.matrix.x1[proxy.obs==1,]),log=TRUE) + ll.w.obs.x1[proxy.obs==0] <- plogis(proxy.params %*% t(proxy.model.matrix.x1[proxy.obs==0,]),log=TRUE, lower.tail=FALSE) + } + + ## assume that the probability of x is a logistic regression depending on z + truth.model.matrix.obs <- model.matrix(truth_formula, df.obs) + n.truth.params <- dim(truth.model.matrix.obs)[2] + truth.params <- params[param.idx:(n.truth.params + param.idx - 1)] + + ll.obs.x0 <- plogis(truth.params %*% t(truth.model.matrix.obs), log=TRUE) + ll.obs.x1 <- plogis(truth.params %*% t(truth.model.matrix.obs), log=TRUE, lower.tail=FALSE) + + ll.obs <- colLogSumExps(rbind(ll.y.x.obs.0 + ll.x.obs.0.x0 + ll.x.obs.1.x0 + ll.obs.x0 + ll.w.obs.x0, + ll.y.x.obs.1 + ll.x.obs.0.x1 + ll.x.obs.1.x1 + ll.obs.x1 + ll.w.obs.x1)) + + ### NOW FOR THE FUN PART. Likelihood of the unobserved data. + ### we have to integrate out x.obs.0, x.obs.1, and x. + + + ## THE OUTCOME + df.unobs <- df[is.na(x.obs)] + df.x.unobs.0 <- copy(df.unobs)[,x:=0] + df.x.unobs.1 <- copy(df.unobs)[,x:=1] + y.unobs <- df.unobs$y + + outcome.model.matrix.x.unobs.0 <- model.matrix(outcome_formula, df.x.unobs.0) + outcome.model.matrix.x.unobs.1 <- model.matrix(outcome_formula, df.x.unobs.1) + + ll.y.unobs.x0 <- dnorm(y.unobs, outcome.params %*% t(outcome.model.matrix.x.unobs.0), sd=sigma.y, log=TRUE) + ll.y.unobs.x1 <- dnorm(y.unobs, outcome.params %*% t(outcome.model.matrix.x.unobs.1), sd=sigma.y, log=TRUE) + + + ## THE UNLABELED DATA + + + ## assume that the two coders are statistically independent conditional on x + ll.x.unobs.0.x0 <- vector(mode='numeric', length=nrow(df.unobs)) + ll.x.unobs.1.x0 <- vector(mode='numeric', length=nrow(df.unobs)) + ll.x.unobs.0.x1 <- vector(mode='numeric', length=nrow(df.unobs)) + ll.x.unobs.1.x1 <- vector(mode='numeric', length=nrow(df.unobs)) + + df.x.unobs.0[,x.obs := 1] + df.x.unobs.1[,x.obs := 1] + + rater.model.matrix.x.unobs.0 <- model.matrix(rater_formula, df.x.unobs.0) + rater.model.matrix.x.unobs.1 <- model.matrix(rater_formula, df.x.unobs.1) + + + ## # probability of rater 0 if x is 0 or 1 + ## ll.x.unobs.0.x0 <- colLogSumExps(rbind(plogis(rater.0.params %*% t(rater.model.matrix.x.unobs.0), log=TRUE), + ## plogis(rater.0.params %*% t(rater.model.matrix.x.unobs.0), log=TRUE, lower.tail=TRUE))) + + ## ll.x.unobs.0.x1 <- colLogSumExps(rbind(plogis(rater.0.params %*% t(rater.model.matrix.x.unobs.1), log=TRUE), + ## plogis(rater.0.params %*% t(rater.model.matrix.x.unobs.1), log=TRUE, lower.tail=TRUE))) + + ## # probability of rater 1 if x is 0 or 1 + ## ll.x.unobs.1.x0 <- colLogSumExps(rbind(plogis(rater.1.params %*% t(rater.model.matrix.x.unobs.0), log=TRUE), + ## plogis(rater.1.params %*% t(rater.model.matrix.x.unobs.0), log=TRUE, lower.tail=TRUE))) + + ## ll.x.unobs.1.x1 <- colLogSumExps(rbind(plogis(rater.1.params %*% t(rater.model.matrix.x.unobs.1), log=TRUE), + ## plogis(rater.1.params %*% t(rater.model.matrix.x.unobs.1), log=TRUE, lower.tail=TRUE))) + + + proxy.unobs <- with(df.unobs, eval(parse(text=proxy.variable))) + proxy.model.matrix.x0.unobs <- model.matrix(proxy_formula, df.x.unobs.0) + proxy.model.matrix.x1.unobs <- model.matrix(proxy_formula, df.x.unobs.1) + + if( (proxy_family$family=="binomial") & (proxy_family$link=='logit')){ + ll.w.unobs.x0 <- vector(mode='numeric',length=dim(proxy.model.matrix.x0)[1]) + ll.w.unobs.x1 <- vector(mode='numeric',length=dim(proxy.model.matrix.x1)[1]) + + + # proxy_formula likelihood using logistic regression + ll.w.unobs.x0[proxy.unobs==1] <- plogis(proxy.params %*% t(proxy.model.matrix.x0.unobs[proxy.unobs==1,]),log=TRUE) + ll.w.unobs.x0[proxy.unobs==0] <- plogis(proxy.params %*% t(proxy.model.matrix.x0.unobs[proxy.unobs==0,]),log=TRUE, lower.tail=FALSE) + + ll.w.unobs.x1[proxy.unobs==1] <- plogis(proxy.params %*% t(proxy.model.matrix.x1.unobs[proxy.unobs==1,]),log=TRUE) + ll.w.unobs.x1[proxy.unobs==0] <- plogis(proxy.params %*% t(proxy.model.matrix.x1.unobs[proxy.unobs==0,]),log=TRUE, lower.tail=FALSE) + } + + truth.model.matrix.unobs <- model.matrix(truth_formula, df.unobs) + + ll.unobs.x0 <- plogis(truth.params %*% t(truth.model.matrix.unobs), log=TRUE) + ll.unobs.x1 <- plogis(truth.params %*% t(truth.model.matrix.unobs), log=TRUE, lower.tail=FALSE) + + ll.unobs <- colLogSumExps(rbind(ll.unobs.x0 + ll.w.unobs.x0 + ll.y.unobs.x0, + ll.unobs.x1 + ll.w.unobs.x1 + ll.y.unobs.x1)) + + return(-1 *( sum(ll.obs) + sum(ll.unobs))) + } + + outcome.params <- colnames(model.matrix(outcome_formula,df)) + lower <- rep(-Inf, length(outcome.params)) + + if(outcome_family$family=='gaussian'){ + params <- c(outcome.params, 'sigma_y') + lower <- c(lower, 0.00001) + } else { + params <- outcome.params + } + + rater.0.params <- colnames(model.matrix(rater_formula,df)) + params <- c(params, paste0('rater_0',rater.0.params)) + lower <- c(lower, rep(-Inf, length(rater.0.params))) + + rater.1.params <- colnames(model.matrix(rater_formula,df)) + params <- c(params, paste0('rater_1',rater.1.params)) + lower <- c(lower, rep(-Inf, length(rater.1.params))) + + proxy.params <- colnames(model.matrix(proxy_formula, df)) + params <- c(params, paste0('proxy_',proxy.params)) + lower <- c(lower, rep(-Inf, length(proxy.params))) + + truth.params <- colnames(model.matrix(truth_formula, df)) + params <- c(params, paste0('truth_', truth.params)) + lower <- c(lower, rep(-Inf, length(truth.params))) + start <- rep(0.1,length(params)) + names(start) <- params + + fit <- optim(start, fn = nll, lower=lower, method='L-BFGS-B', hessian=TRUE, control=list(maxit=1e6)) + return(fit) +} + + measerr_mle <- function(df, outcome_formula, outcome_family=gaussian(), proxy_formula, proxy_family=binomial(link='logit'), truth_formula, truth_family=binomial(link='logit')){ measrr_mle_nll <- function(params){ df.obs <- model.frame(outcome_formula, df) - proxy.variable <- all.vars(proxy_formula)[1] proxy.model.matrix <- model.matrix(proxy_formula, df) - response.var <- all.vars(outcome_formula)[1] y.obs <- with(df.obs,eval(parse(text=response.var))) - + outcome.model.matrix <- model.matrix(outcome_formula, df) param.idx <- 1 @@ -125,7 +319,7 @@ measerr_mle <- function(df, outcome_formula, outcome_family=gaussian(), proxy_fo sigma.y <- params[param.idx] param.idx <- param.idx + 1 - # outcome_formula likelihood using linear regression + # outcome_formula likelihood using linear regression ll.y.obs <- dnorm(y.obs, outcome.params %*% t(outcome.model.matrix),sd=sigma.y, log=TRUE) } @@ -138,7 +332,7 @@ measerr_mle <- function(df, outcome_formula, outcome_family=gaussian(), proxy_fo if( (proxy_family$family=="binomial") & (proxy_family$link=='logit')){ ll.w.obs <- vector(mode='numeric',length=dim(proxy.model.matrix)[1]) - # proxy_formula likelihood using logistic regression + # proxy_formula likelihood using logistic regression ll.w.obs[proxy.obs==1] <- plogis(proxy.params %*% t(proxy.model.matrix[proxy.obs==1,]),log=TRUE) ll.w.obs[proxy.obs==0] <- plogis(proxy.params %*% t(proxy.model.matrix[proxy.obs==0,]),log=TRUE, lower.tail=FALSE) } @@ -154,12 +348,12 @@ measerr_mle <- function(df, outcome_formula, outcome_family=gaussian(), proxy_fo if( (truth_family$family=="binomial") & (truth_family$link=='logit')){ ll.x.obs <- vector(mode='numeric',length=dim(truth.model.matrix)[1]) - # truth_formula likelihood using logistic regression + # truth_formula likelihood using logistic regression ll.x.obs[truth.obs==1] <- plogis(truth.params %*% t(truth.model.matrix[truth.obs==1,]),log=TRUE) ll.x.obs[truth.obs==0] <- plogis(truth.params %*% t(truth.model.matrix[truth.obs==0,]),log=TRUE, lower.tail=FALSE) } - # add the three likelihoods + # add the three likelihoods ll.obs <- sum(ll.y.obs + ll.w.obs + ll.x.obs) ## likelihood for the predicted data @@ -177,9 +371,9 @@ measerr_mle <- function(df, outcome_formula, outcome_family=gaussian(), proxy_fo outcome.model.matrix.x1 <- model.matrix(outcome_formula, df.unobs.x1) if(outcome_family$family=="gaussian"){ - # likelihood of outcome - ll.y.x0 <- dnorm(outcome.unobs, outcome.params %*% t(outcome.model.matrix.x0), sd=sigma.y, log=TRUE) - ll.y.x1 <- dnorm(outcome.unobs, outcome.params %*% t(outcome.model.matrix.x1), sd=sigma.y, log=TRUE) + # likelihood of outcome + ll.y.x0 <- dnorm(outcome.unobs, outcome.params %*% t(outcome.model.matrix.x0), sd=sigma.y, log=TRUE) + ll.y.x1 <- dnorm(outcome.unobs, outcome.params %*% t(outcome.model.matrix.x1), sd=sigma.y, log=TRUE) } if( (proxy_family$family=='binomial') & (proxy_family$link=='logit')){ @@ -190,7 +384,7 @@ measerr_mle <- function(df, outcome_formula, outcome_family=gaussian(), proxy_fo ll.w.x0 <- vector(mode='numeric', length=dim(df.unobs)[1]) ll.w.x1 <- vector(mode='numeric', length=dim(df.unobs)[1]) - # likelihood of proxy + # likelihood of proxy ll.w.x0[proxy.unobs==1] <- plogis(proxy.params %*% t(proxy.model.matrix.x0[proxy.unobs==1,]), log=TRUE) ll.w.x1[proxy.unobs==1] <- plogis(proxy.params %*% t(proxy.model.matrix.x1[proxy.unobs==1,]), log=TRUE) @@ -200,7 +394,7 @@ measerr_mle <- function(df, outcome_formula, outcome_family=gaussian(), proxy_fo if(truth_family$link=='logit'){ truth.model.matrix <- model.matrix(truth_formula, df.unobs.x0) - # likelihood of truth + # likelihood of truth ll.x.x1 <- plogis(truth.params %*% t(truth.model.matrix), log=TRUE) ll.x.x0 <- plogis(truth.params %*% t(truth.model.matrix), log=TRUE, lower.tail=FALSE) }