X-Git-Url: https://code.communitydata.science/ml_measurement_error_public.git/blobdiff_plain/47e9367ed5c61b721bdc17cddd76bced4f8ed621..c42b94110b18264fdd66ada100ee05232b7b81bb:/simulations/measerr_methods.R?ds=sidebyside diff --git a/simulations/measerr_methods.R b/simulations/measerr_methods.R index 6bf8c3f..63f8bc1 100644 --- a/simulations/measerr_methods.R +++ b/simulations/measerr_methods.R @@ -1,6 +1,7 @@ library(formula.tools) library(matrixStats) - +library(optimx) +library(bbmle) ## df: dataframe to model ## outcome_formula: formula for y | x, z ## outcome_family: family for y | x, z @@ -17,7 +18,7 @@ library(matrixStats) ## outcome_formula <- y ~ x + z; proxy_formula <- w_pred ~ y + x + z + x:z + x:y + z:y -measerr_mle_dv <- function(df, outcome_formula, outcome_family=binomial(link='logit'), proxy_formula, proxy_family=binomial(link='logit')){ +measerr_mle_dv <- function(df, outcome_formula, outcome_family=binomial(link='logit'), proxy_formula, proxy_family=binomial(link='logit'),method='optim'){ nll <- function(params){ df.obs <- model.frame(outcome_formula, df) @@ -98,23 +99,33 @@ measerr_mle_dv <- function(df, outcome_formula, outcome_family=binomial(link='lo 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)) + if(method=='optim'){ + fit <- optim(start, fn = nll, lower=lower, method='L-BFGS-B', hessian=TRUE, control=list(maxit=1e6)) + } else { + quoted.names <- gsub("[\\(\\)]",'',names(start)) + print(quoted.names) + text <- paste("function(", paste0(quoted.names,'=',start,collapse=','),"){params<-c(",paste0(quoted.names,collapse=','),");return(nll(params))}") + + measerr_mle_nll <- eval(parse(text=text)) + names(start) <- quoted.names + names(lower) <- quoted.names + fit <- mle2(minuslogl=measerr_mle_nll, start=start, lower=lower, parnames=params,control=list(maxit=1e6),method='L-BFGS-B') + } 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) +measerr_mle <- function(df, outcome_formula, outcome_family=gaussian(), proxy_formula, proxy_family=binomial(link='logit'), truth_formula, truth_family=binomial(link='logit'),method='optim'){ - 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) + df.obs <- model.frame(outcome_formula, df) + response.var <- all.vars(outcome_formula)[1] + proxy.variable <- all.vars(proxy_formula)[1] + truth.variable <- all.vars(truth_formula)[1] + outcome.model.matrix <- model.matrix(outcome_formula, df) + proxy.model.matrix <- model.matrix(proxy_formula, df) + y.obs <- with(df.obs,eval(parse(text=response.var))) + measerr_mle_nll <- function(params){ param.idx <- 1 n.outcome.model.covars <- dim(outcome.model.matrix)[2] outcome.params <- params[param.idx:n.outcome.model.covars] @@ -124,8 +135,7 @@ measerr_mle <- function(df, outcome_formula, outcome_family=gaussian(), proxy_fo if(outcome_family$family == "gaussian"){ 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,13 +148,13 @@ 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) } df.obs <- model.frame(truth_formula, df) - truth.variable <- all.vars(truth_formula)[1] + truth.obs <- with(df.obs, eval(parse(text=truth.variable))) truth.model.matrix <- model.matrix(truth_formula,df) n.truth.model.covars <- dim(truth.model.matrix)[2] @@ -154,12 +164,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 +187,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 +200,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 +210,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) } @@ -231,8 +241,356 @@ measerr_mle <- function(df, outcome_formula, outcome_family=gaussian(), proxy_fo lower <- c(lower, rep(-Inf, length(truth.params))) start <- rep(0.1,length(params)) names(start) <- params + + if(method=='optim'){ + fit <- optim(start, fn = measerr_mle_nll, lower=lower, method='L-BFGS-B', hessian=TRUE, control=list(maxit=1e6)) + } else { # method='mle2' + + quoted.names <- gsub("[\\(\\)]",'',names(start)) + + text <- paste("function(", paste0(quoted.names,'=',start,collapse=','),"){params<-c(",paste0(quoted.names,collapse=','),");return(measerr_mle_nll(params))}") + + measerr_mle_nll_mle <- eval(parse(text=text)) + names(start) <- quoted.names + names(lower) <- quoted.names + fit <- mle2(minuslogl=measerr_mle_nll_mle, start=start, lower=lower, parnames=params,control=list(maxit=1e6),method='L-BFGS-B') + } + + return(fit) +} + +## Experimental, but probably works. +measerr_irr_mle <- function(df, outcome_formula, outcome_family=gaussian(), coder_formulas=list(x.obs.0~x, x.obs.1~x), proxy_formula, proxy_family=binomial(link='logit'), truth_formula, truth_family=binomial(link='logit'),method='optim'){ + + ### in this scenario, the ground truth also has measurement error, but we have repeated measures for it. + # this time we never get to observe the true X + outcome.model.matrix <- model.matrix(outcome_formula, df) + proxy.model.matrix <- model.matrix(proxy_formula, df) + response.var <- all.vars(outcome_formula)[1] + proxy.var <- all.vars(proxy_formula)[1] + param.var <- all.vars(truth_formula)[1] + truth.var<- all.vars(truth_formula)[1] + y <- with(df,eval(parse(text=response.var))) + + nll <- function(params){ + param.idx <- 1 + n.outcome.model.covars <- dim(outcome.model.matrix)[2] + outcome.params <- params[param.idx:n.outcome.model.covars] + param.idx <- param.idx + n.outcome.model.covars + + + if(outcome_family$family == "gaussian"){ + sigma.y <- params[param.idx] + param.idx <- param.idx + 1 + } + + + n.proxy.model.covars <- dim(proxy.model.matrix)[2] + proxy.params <- params[param.idx:(n.proxy.model.covars+param.idx-1)] + param.idx <- param.idx + n.proxy.model.covars + + df.temp <- copy(df) + + if((truth_family$family == "binomial") + & (truth_family$link=='logit')){ + integrate.grid <- expand.grid(replicate(1 + length(coder_formulas), c(0,1), simplify=FALSE)) + ll.parts <- matrix(nrow=nrow(df),ncol=nrow(integrate.grid)) + for(i in 1:nrow(integrate.grid)){ + # setup the dataframe for this row + row <- integrate.grid[i,] + + df.temp[[param.var]] <- row[[1]] + ci <- 2 + for(coder_formula in coder_formulas){ + coder.var <- all.vars(coder_formula)[1] + df.temp[[coder.var]] <- row[[ci]] + ci <- ci + 1 + } + + outcome.model.matrix <- model.matrix(outcome_formula, df.temp) + if(outcome_family$family == "gaussian"){ + ll.y <- dnorm(y, outcome.params %*% t(outcome.model.matrix), sd=sigma.y, log=TRUE) + } + + if(proxy_family$family=="binomial" & (proxy_family$link=='logit')){ + proxy.model.matrix <- model.matrix(proxy_formula, df.temp) + ll.w <- vector(mode='numeric', length=dim(proxy.model.matrix)[1]) + proxyvar <- with(df.temp,eval(parse(text=proxy.var))) + ll.w[proxyvar==1] <- plogis(proxy.params %*% t(proxy.model.matrix[proxyvar==1,]),log=TRUE) + ll.w[proxyvar==0] <- plogis(proxy.params %*% t(proxy.model.matrix[proxyvar==0,]),log=TRUE,lower.tail=FALSE) + } + + ## probability of the coded variables + coder.lls <- matrix(nrow=nrow(df.temp),ncol=length(coder_formulas)) + ci <- 1 + for(coder_formula in coder_formulas){ + coder.model.matrix <- model.matrix(coder_formula, df.temp) + n.coder.model.covars <- dim(coder.model.matrix)[2] + coder.params <- params[param.idx:(n.coder.model.covars + param.idx - 1)] + param.idx <- param.idx + n.coder.model.covars + coder.var <- all.vars(coder_formula)[1] + x.obs <- with(df.temp, eval(parse(text=coder.var))) + true.codervar <- df[[all.vars(coder_formula)[1]]] + + ll.coder <- vector(mode='numeric', length=dim(coder.model.matrix)[1]) + ll.coder[x.obs==1] <- plogis(coder.params %*% t(coder.model.matrix[x.obs==1,]),log=TRUE) + ll.coder[x.obs==0] <- plogis(coder.params %*% t(coder.model.matrix[x.obs==0,]),log=TRUE,lower.tail=FALSE) + + # don't count when we know the observed value, unless we're accounting for observed value + ll.coder[(!is.na(true.codervar)) & (true.codervar != x.obs)] <- NA + coder.lls[,ci] <- ll.coder + ci <- ci + 1 + } + + truth.model.matrix <- model.matrix(truth_formula, df.temp) + n.truth.model.covars <- dim(truth.model.matrix)[2] + truth.params <- params[param.idx:(n.truth.model.covars + param.idx - 1)] + + for(coder_formula in coder_formulas){ + coder.model.matrix <- model.matrix(coder_formula, df.temp) + n.coder.model.covars <- dim(coder.model.matrix)[2] + param.idx <- param.idx - n.coder.model.covars + } + + x <- with(df.temp, eval(parse(text=truth.var))) + ll.truth <- vector(mode='numeric', length=dim(truth.model.matrix)[1]) + ll.truth[x==1] <- plogis(truth.params %*% t(truth.model.matrix[x==1,]), log=TRUE) + ll.truth[x==0] <- plogis(truth.params %*% t(truth.model.matrix[x==0,]), log=TRUE, lower.tail=FALSE) + + true.truthvar <- df[[all.vars(truth_formula)[1]]] + + if(!is.null(true.truthvar)){ + # ll.truth[(!is.na(true.truthvar)) & (true.truthvar != truthvar)] <- -Inf + # ll.truth[(!is.na(true.truthvar)) & (true.truthvar == truthvar)] <- 0 + } + ll.parts[,i] <- ll.y + ll.w + apply(coder.lls,1,sum) + ll.truth + + } + + lls <- rowLogSumExps(ll.parts,na.rm=TRUE) + + ## likelihood of observed data + target <- -1 * sum(lls) + return(target) + } + } + + 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 + } + + proxy.params <- colnames(model.matrix(proxy_formula, df)) + params <- c(params, paste0('proxy_',proxy.params)) + positive.params <- paste0('proxy_',truth.var) + lower <- c(lower, rep(-Inf, length(proxy.params))) + names(lower) <- params + lower[positive.params] <- 0.01 + ci <- 0 + + for(coder_formula in coder_formulas){ + coder.params <- colnames(model.matrix(coder_formula,df)) + params <- c(params, paste0('coder_',ci,coder.params)) + positive.params <- paste0('coder_', ci, truth.var) + ci <- ci + 1 + lower <- c(lower, rep(-Inf, length(coder.params))) + names(lower) <- params + lower[positive.params] <- 0.01 + } + + 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 + names(lower) <- params - fit <- optim(start, fn = measrr_mle_nll, lower=lower, method='L-BFGS-B', hessian=TRUE, control=list(maxit=1e6)) + if(method=='optim'){ + print(start) + fit <- optim(start, fn = nll, lower=lower, method='L-BFGS-B', hessian=TRUE, control=list(maxit=1e6)) + } else { + + quoted.names <- gsub("[\\(\\)]",'',names(start)) + print(quoted.names) + text <- paste("function(", paste0(quoted.names,'=',start,collapse=','),"){params<-c(",paste0(quoted.names,collapse=','),");return(nll(params))}") + + measerr_mle_nll <- eval(parse(text=text)) + names(start) <- quoted.names + names(lower) <- quoted.names + fit <- mle2(minuslogl=measerr_mle_nll, start=start, lower=lower, method='L-BFGS-B',control=list(maxit=1e6)) + } return(fit) } + +## Experimental, and does not work. +measerr_irr_mle_dv <- function(df, outcome_formula, outcome_family=binomial(link='logit'), coder_formulas=list(y.obs.0~y+w_pred+y.obs.1,y.obs.1~y+w_pred+y.obs.0), proxy_formula=w_pred~y, proxy_family=binomial(link='logit'),method='optim'){ + integrate.grid <- expand.grid(replicate(1 + length(coder_formulas), c(0,1), simplify=FALSE)) + print(integrate.grid) + + + outcome.model.matrix <- model.matrix(outcome_formula, df) + n.outcome.model.covars <- dim(outcome.model.matrix)[2] + + + ### in this scenario, the ground truth also has measurement error, but we have repeated measures for it. + # this time we never get to observe the true X + nll <- function(params){ + param.idx <- 1 + outcome.params <- params[param.idx:n.outcome.model.covars] + param.idx <- param.idx + n.outcome.model.covars + proxy.model.matrix <- model.matrix(proxy_formula, df) + n.proxy.model.covars <- dim(proxy.model.matrix)[2] + response.var <- all.vars(outcome_formula)[1] + + if(outcome_family$family == "gaussian"){ + sigma.y <- params[param.idx] + param.idx <- param.idx + 1 + } + + proxy.params <- params[param.idx:(n.proxy.model.covars+param.idx-1)] + param.idx <- param.idx + n.proxy.model.covars + + df.temp <- copy(df) + + if((outcome_family$family == "binomial") + & (outcome_family$link=='logit')){ + ll.parts <- matrix(nrow=nrow(df),ncol=nrow(integrate.grid)) + for(i in 1:nrow(integrate.grid)){ + # setup the dataframe for this row + row <- integrate.grid[i,] + + df.temp[[response.var]] <- row[[1]] + ci <- 2 + for(coder_formula in coder_formulas){ + codervar <- all.vars(coder_formula)[1] + df.temp[[codervar]] <- row[[ci]] + ci <- ci + 1 + } + + outcome.model.matrix <- model.matrix(outcome_formula, df.temp) + if(outcome_family$family == "gaussian"){ + ll.y <- dnorm(df.temp[[response.var]], outcome.params %*% t(outcome.model.matrix), sd=sigma.y, log=T) + } + + if(outcome_family$family == "binomial" & (outcome_family$link=='logit')){ + ll.y <- vector(mode='numeric',length=nrow(df.temp)) + ll.y[df.temp[[response.var]]==1] <- plogis( outcome.params %*% t(outcome.model.matrix), log=TRUE) + ll.y[df.temp[[response.var]]==0] <- plogis( outcome.params %*% t(outcome.model.matrix), log=TRUE,lower.tail=FALSE) + } + + if(proxy_family$family=="binomial" & (proxy_family$link=='logit')){ + proxy.model.matrix <- model.matrix(proxy_formula, df.temp) + ll.w <- vector(mode='numeric', length=dim(proxy.model.matrix)[1]) + proxyvar <- with(df.temp,eval(parse(text=all.vars(proxy_formula)[1]))) + ll.w[proxyvar==1] <- plogis(proxy.params %*% t(proxy.model.matrix[proxyvar==1,]),log=TRUE) + ll.w[proxyvar==0] <- plogis(proxy.params %*% t(proxy.model.matrix[proxyvar==0,]),log=TRUE,lower.tail=FALSE) + } + + ## probability of the coded variables + coder.lls <- matrix(nrow=nrow(df.temp),ncol=length(coder_formulas)) + ci <- 1 + for(coder_formula in coder_formulas){ + coder.model.matrix <- model.matrix(coder_formula, df.temp) + n.coder.model.covars <- dim(coder.model.matrix)[2] + coder.params <- params[param.idx:(n.coder.model.covars + param.idx - 1)] + param.idx <- param.idx + n.coder.model.covars + codervar <- with(df.temp, eval(parse(text=all.vars(coder_formula)[1]))) + true.codervar <- df[[all.vars(coder_formula)[1]]] + + ll.coder <- vector(mode='numeric', length=dim(coder.model.matrix)[1]) + ll.coder[codervar==1] <- plogis(coder.params %*% t(coder.model.matrix[codervar==1,]),log=TRUE) + ll.coder[codervar==0] <- plogis(coder.params %*% t(coder.model.matrix[codervar==0,]),log=TRUE,lower.tail=FALSE) + + # don't count when we know the observed value, unless we're accounting for observed value + ll.coder[(!is.na(true.codervar)) & (true.codervar != codervar)] <- NA + coder.lls[,ci] <- ll.coder + ci <- ci + 1 + } + + for(coder_formula in coder_formulas){ + coder.model.matrix <- model.matrix(coder_formula, df.temp) + n.coder.model.covars <- dim(coder.model.matrix)[2] + param.idx <- param.idx - n.coder.model.covars + } + + ll.parts[,i] <- ll.y + ll.w + apply(coder.lls,1,function(x) sum(x)) + + } + + lls <- rowLogSumExps(ll.parts,na.rm=TRUE) + + ## likelihood of observed data + target <- -1 * sum(lls) + print(target) + print(params) + return(target) + } + } + + outcome.params <- colnames(model.matrix(outcome_formula,df)) + response.var <- all.vars(outcome_formula)[1] + 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 + } + + ## constrain the model of the coder and proxy vars + ## this is to ensure identifiability + ## it is a safe assumption because the coders aren't hostile (wrong more often than right) + ## so we can assume that y ~Bw, B is positive + proxy.params <- colnames(model.matrix(proxy_formula, df)) + positive.params <- paste0('proxy_',response.var) + params <- c(params, paste0('proxy_',proxy.params)) + lower <- c(lower, rep(-Inf, length(proxy.params))) + names(lower) <- params + lower[positive.params] <- 0.001 + + ci <- 0 + for(coder_formula in coder_formulas){ + coder.params <- colnames(model.matrix(coder_formula,df)) + latent.coder.params <- coder.params %in% response.var + params <- c(params, paste0('coder_',ci,coder.params)) + positive.params <- paste0('coder_',ci,response.var) + ci <- ci + 1 + lower <- c(lower, rep(-Inf, length(coder.params))) + names(lower) <-params + lower[positive.params] <- 0.001 + } + + ## init by using the "loco model" + temp.df <- copy(df) + temp.df <- temp.df[y.obs.1 == y.obs.0, y:=y.obs.1] + loco.model <- glm(outcome_formula, temp.df, family=outcome_family) + + start <- rep(1,length(params)) + names(start) <- params + start[names(coef(loco.model))] <- coef(loco.model) + names(lower) <- params + if(method=='optim'){ + print(lower) + fit <- optim(start, fn = nll, lower=lower, method='L-BFGS-B', hessian=TRUE,control=list(maxit=1e6)) + } else { + + quoted.names <- gsub("[\\(\\)]",'',names(start)) + print(quoted.names) + text <- paste("function(", paste0(quoted.names,'=',start,collapse=','),"){params<-c(",paste0(quoted.names,collapse=','),");return(nll(params))}") + + measerr_mle_nll <- eval(parse(text=text)) + names(start) <- quoted.names + names(lower) <- quoted.names + fit <- mle2(minuslogl=measerr_mle_nll, start=start, lower=lower, parnames=params,control=list(maxit=1e6),method='L-BFGS-B') + } + + return(fit) +} +