simulations/measerr_methods.R

   1 library(formula.tools)
   2 library(matrixStats)
   3 library(optimx)
   4 library(bbmle)
   5 ## df: dataframe to model
   6 ## outcome_formula: formula for y | x, z
   7 ## outcome_family: family for y | x, z
   8 ## proxy_formula: formula for w | x, z, y
   9 ## proxy_family: family for w | x, z, y
  10 ## truth_formula: formula for x | z
  11 ## truth_family: family for x | z
  12
  13 ### ideal formulas for example 1
  14 # test.fit.1 <- measerr_mle(df, y ~ x + z, gaussian(), w_pred ~ x, binomial(link='logit'), x ~ z)
  15
  16 ### ideal formulas for example 2
  17 # test.fit.2 <- measerr_mle(df, y ~ x + z, gaussian(), w_pred ~ x + z + y + y:x, binomial(link='logit'), x ~ z)
  18
  19
  20 ## outcome_formula <- y ~ x + z; proxy_formula <- w_pred ~ y + x + z + x:z + x:y + z:y
  21 measerr_mle_dv <- function(df, outcome_formula, outcome_family=binomial(link='logit'), proxy_formula, proxy_family=binomial(link='logit'),method='optim'){
  22
  23     nll <- function(params){
  24         df.obs <- model.frame(outcome_formula, df)
  25         proxy.variable <- all.vars(proxy_formula)[1]
  26         proxy.model.matrix <- model.matrix(proxy_formula, df)
  27         response.var <- all.vars(outcome_formula)[1]
  28         y.obs <- with(df.obs,eval(parse(text=response.var)))
  29         outcome.model.matrix <- model.matrix(outcome_formula, df.obs)
  30
  31         param.idx <- 1
  32         n.outcome.model.covars <- dim(outcome.model.matrix)[2]
  33         outcome.params <- params[param.idx:n.outcome.model.covars]
  34         param.idx <- param.idx + n.outcome.model.covars
  35
  36         if((outcome_family$family == "binomial") & (outcome_family$link == 'logit')){
  37             ll.y.obs <- vector(mode='numeric', length=length(y.obs))
  38             ll.y.obs[y.obs==1] <- plogis(outcome.params %*% t(outcome.model.matrix[y.obs==1,]),log=TRUE)
  39             ll.y.obs[y.obs==0] <- plogis(outcome.params %*% t(outcome.model.matrix[y.obs==0,]),log=TRUE,lower.tail=FALSE)
  40         }
  41
  42         df.obs <- model.frame(proxy_formula,df)
  43         n.proxy.model.covars <- dim(proxy.model.matrix)[2]
  44         proxy.params <- params[param.idx:(n.proxy.model.covars+param.idx-1)]
  45
  46         param.idx <- param.idx + n.proxy.model.covars
  47         proxy.obs <- with(df.obs, eval(parse(text=proxy.variable)))
  48
  49         if( (proxy_family$family=="binomial") & (proxy_family$link=='logit')){
  50             ll.w.obs <- vector(mode='numeric',length=dim(proxy.model.matrix)[1])
  51             ll.w.obs[proxy.obs==1] <- plogis(proxy.params %*% t(proxy.model.matrix[proxy.obs==1,]),log=TRUE)
  52             ll.w.obs[proxy.obs==0] <- plogis(proxy.params %*% t(proxy.model.matrix[proxy.obs==0,]),log=TRUE, lower.tail=FALSE)
  53         }
  54
  55         ll.obs <- sum(ll.y.obs + ll.w.obs)
  56
  57         df.unobs <- df[is.na(df[[response.var]])]
  58         df.unobs.y1 <- copy(df.unobs)
  59         df.unobs.y1[[response.var]] <- 1
  60         df.unobs.y0 <- copy(df.unobs)
  61         df.unobs.y0[[response.var]] <- 0
  62
  63         ## integrate out y
  64         outcome.model.matrix.y1 <- model.matrix(outcome_formula, df.unobs.y1)
  65
  66         if((outcome_family$family == "binomial") & (outcome_family$link == 'logit')){
  67             ll.y.unobs.1 <- vector(mode='numeric', length=dim(outcome.model.matrix.y1)[1])
  68             ll.y.unobs.0 <- vector(mode='numeric', length=dim(outcome.model.matrix.y1)[1])
  69             ll.y.unobs.1 <- plogis(outcome.params %*% t(outcome.model.matrix.y1),log=TRUE)
  70             ll.y.unobs.0 <- plogis(outcome.params %*% t(outcome.model.matrix.y1),log=TRUE,lower.tail=FALSE)
  71         }
  72
  73         proxy.model.matrix.y1 <- model.matrix(proxy_formula, df.unobs.y1)
  74         proxy.model.matrix.y0 <- model.matrix(proxy_formula, df.unobs.y0)
  75         proxy.unobs <- with(df.unobs, eval(parse(text=proxy.variable)))
  76
  77         if( (proxy_family$family=="binomial") & (proxy_family$link=='logit')){
  78             ll.w.unobs.1 <- vector(mode='numeric',length=dim(proxy.model.matrix.y1)[1])
  79             ll.w.unobs.0 <- vector(mode='numeric',length=dim(proxy.model.matrix.y0)[1])
  80             ll.w.unobs.1[proxy.unobs==1] <- plogis(proxy.params %*% t(proxy.model.matrix.y1[proxy.unobs==1,]),log=TRUE)
  81             ll.w.unobs.1[proxy.unobs==0] <- plogis(proxy.params %*% t(proxy.model.matrix.y1[proxy.unobs==0,]),log=TRUE, lower.tail=FALSE)
  82
  83             ll.w.unobs.0[proxy.unobs==1] <- plogis(proxy.params %*% t(proxy.model.matrix.y0[proxy.unobs==1,]),log=TRUE)
  84             ll.w.unobs.0[proxy.unobs==0] <- plogis(proxy.params %*% t(proxy.model.matrix.y0[proxy.unobs==0,]),log=TRUE, lower.tail=FALSE)
  85         }
  86
  87         ll.unobs.1 <- ll.y.unobs.1 + ll.w.unobs.1
  88         ll.unobs.0 <- ll.y.unobs.0 + ll.w.unobs.0
  89         ll.unobs <- sum(colLogSumExps(rbind(ll.unobs.1,ll.unobs.0)))
  90         ll <- ll.unobs + ll.obs
  91         return(-ll)
  92     }
  93
  94     params <- colnames(model.matrix(outcome_formula,df))
  95     lower <- rep(-Inf, length(params))
  96     proxy.params <- colnames(model.matrix(proxy_formula, df))
  97     params <- c(params, paste0('proxy_',proxy.params))
  98     lower <- c(lower, rep(-Inf, length(proxy.params)))
  99     start <- rep(0.1,length(params))
 100     names(start) <- params
 101
 102     if(method=='optim'){
 103         fit <- optim(start, fn = nll, lower=lower, method='L-BFGS-B', hessian=TRUE, control=list(maxit=1e6))
 104     } else {
 105         quoted.names <- gsub("[\\(\\)]",'',names(start))
 106         print(quoted.names)
 107         text <- paste("function(", paste0(quoted.names,'=',start,collapse=','),"){params<-c(",paste0(quoted.names,collapse=','),");return(nll(params))}")
 108
 109         measerr_mle_nll <- eval(parse(text=text))
 110         names(start) <- quoted.names
 111         names(lower) <- quoted.names
 112         fit <- mle2(minuslogl=measerr_mle_nll, start=start, lower=lower, parnames=params,control=list(maxit=1e6),method='L-BFGS-B')
 113     }
 114     return(fit)
 115 }
 116
 117
 118 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'){
 119
 120     df.obs <- model.frame(outcome_formula, df)
 121     response.var <- all.vars(outcome_formula)[1]
 122     proxy.variable <- all.vars(proxy_formula)[1]
 123     truth.variable <- all.vars(truth_formula)[1]
 124     outcome.model.matrix <- model.matrix(outcome_formula, df)
 125     proxy.model.matrix <- model.matrix(proxy_formula, df)
 126     y.obs <- with(df.obs,eval(parse(text=response.var)))
 127
 128     measerr_mle_nll <- function(params){
 129         param.idx <- 1
 130         n.outcome.model.covars <- dim(outcome.model.matrix)[2]
 131         outcome.params <- params[param.idx:n.outcome.model.covars]
 132         param.idx <- param.idx + n.outcome.model.covars
 133
 134         ## likelihood for the fully observed data
 135         if(outcome_family$family == "gaussian"){
 136             sigma.y <- params[param.idx]
 137             param.idx <- param.idx + 1
 138                                         #  outcome_formula likelihood using linear regression
 139             ll.y.obs <- dnorm(y.obs, outcome.params %*% t(outcome.model.matrix),sd=sigma.y, log=TRUE)
 140         }
 141
 142         df.obs <- model.frame(proxy_formula,df)
 143         n.proxy.model.covars <- dim(proxy.model.matrix)[2]
 144         proxy.params <- params[param.idx:(n.proxy.model.covars+param.idx-1)]
 145         param.idx <- param.idx + n.proxy.model.covars
 146         proxy.obs <- with(df.obs, eval(parse(text=proxy.variable)))
 147
 148         if( (proxy_family$family=="binomial") & (proxy_family$link=='logit')){
 149             ll.w.obs <- vector(mode='numeric',length=dim(proxy.model.matrix)[1])
 150
 151                                         # proxy_formula likelihood using logistic regression
 152             ll.w.obs[proxy.obs==1] <- plogis(proxy.params %*% t(proxy.model.matrix[proxy.obs==1,]),log=TRUE)
 153             ll.w.obs[proxy.obs==0] <- plogis(proxy.params %*% t(proxy.model.matrix[proxy.obs==0,]),log=TRUE, lower.tail=FALSE)
 154         }
 155
 156         df.obs <- model.frame(truth_formula, df)
 157
 158         truth.obs <- with(df.obs, eval(parse(text=truth.variable)))
 159         truth.model.matrix <- model.matrix(truth_formula,df)
 160         n.truth.model.covars <- dim(truth.model.matrix)[2]
 161
 162         truth.params <- params[param.idx:(n.truth.model.covars + param.idx - 1)]
 163
 164         if( (truth_family$family=="binomial") & (truth_family$link=='logit')){
 165             ll.x.obs <- vector(mode='numeric',length=dim(truth.model.matrix)[1])
 166
 167                                         # truth_formula likelihood using logistic regression
 168             ll.x.obs[truth.obs==1] <- plogis(truth.params %*% t(truth.model.matrix[truth.obs==1,]),log=TRUE)
 169             ll.x.obs[truth.obs==0] <- plogis(truth.params %*% t(truth.model.matrix[truth.obs==0,]),log=TRUE, lower.tail=FALSE)
 170         }
 171
 172                                         # add the three likelihoods
 173         ll.obs <- sum(ll.y.obs + ll.w.obs + ll.x.obs)
 174
 175         ## likelihood for the predicted data
 176         ## integrate out the "truth" variable.
 177
 178         if(truth_family$family=='binomial'){
 179             df.unobs <- df[is.na(eval(parse(text=truth.variable)))]
 180             df.unobs.x1 <- copy(df.unobs)
 181             df.unobs.x1[,'x'] <- 1
 182             df.unobs.x0 <- copy(df.unobs)
 183             df.unobs.x0[,'x'] <- 0
 184             outcome.unobs <- with(df.unobs, eval(parse(text=response.var)))
 185
 186             outcome.model.matrix.x0 <- model.matrix(outcome_formula, df.unobs.x0)
 187             outcome.model.matrix.x1 <- model.matrix(outcome_formula, df.unobs.x1)
 188             if(outcome_family$family=="gaussian"){
 189
 190                                         # likelihood of outcome
 191                 ll.y.x0 <- dnorm(outcome.unobs, outcome.params %*% t(outcome.model.matrix.x0), sd=sigma.y, log=TRUE)
 192                 ll.y.x1 <- dnorm(outcome.unobs, outcome.params %*% t(outcome.model.matrix.x1), sd=sigma.y, log=TRUE)
 193             }
 194
 195             if( (proxy_family$family=='binomial') & (proxy_family$link=='logit')){
 196
 197                 proxy.model.matrix.x0 <- model.matrix(proxy_formula, df.unobs.x0)
 198                 proxy.model.matrix.x1 <- model.matrix(proxy_formula, df.unobs.x1)
 199                 proxy.unobs <- df.unobs[[proxy.variable]]
 200                 ll.w.x0 <- vector(mode='numeric', length=dim(df.unobs)[1])
 201                 ll.w.x1 <- vector(mode='numeric', length=dim(df.unobs)[1])
 202
 203                                         # likelihood of proxy
 204                 ll.w.x0[proxy.unobs==1] <- plogis(proxy.params %*% t(proxy.model.matrix.x0[proxy.unobs==1,]), log=TRUE)
 205                 ll.w.x1[proxy.unobs==1] <- plogis(proxy.params %*% t(proxy.model.matrix.x1[proxy.unobs==1,]), log=TRUE)
 206
 207                 ll.w.x0[proxy.unobs==0] <- plogis(proxy.params %*% t(proxy.model.matrix.x0[proxy.unobs==0,]), log=TRUE,lower.tail=FALSE)
 208                 ll.w.x1[proxy.unobs==0] <- plogis(proxy.params %*% t(proxy.model.matrix.x1[proxy.unobs==0,]), log=TRUE,lower.tail=FALSE)
 209             }
 210
 211             if(truth_family$link=='logit'){
 212                 truth.model.matrix <- model.matrix(truth_formula, df.unobs.x0)
 213                                         # likelihood of truth
 214                 ll.x.x1 <- plogis(truth.params %*% t(truth.model.matrix), log=TRUE)
 215                 ll.x.x0 <- plogis(truth.params %*% t(truth.model.matrix), log=TRUE, lower.tail=FALSE)
 216             }
 217         }
 218
 219         ll.x0 <- ll.y.x0 + ll.w.x0 + ll.x.x0
 220         ll.x1 <- ll.y.x1 + ll.w.x1 + ll.x.x1
 221         ll.unobs <- sum(colLogSumExps(rbind(ll.x0, ll.x1)))
 222         return(-(ll.unobs + ll.obs))
 223     }
 224
 225     outcome.params <- colnames(model.matrix(outcome_formula,df))
 226     lower <- rep(-Inf, length(outcome.params))
 227
 228     if(outcome_family$family=='gaussian'){
 229         params <- c(outcome.params, 'sigma_y')
 230         lower <- c(lower, 0.00001)
 231     } else {
 232         params <- outcome.params
 233     }
 234
 235     proxy.params <- colnames(model.matrix(proxy_formula, df))
 236     params <- c(params, paste0('proxy_',proxy.params))
 237     lower <- c(lower, rep(-Inf, length(proxy.params)))
 238
 239     truth.params <- colnames(model.matrix(truth_formula, df))
 240     params <- c(params, paste0('truth_', truth.params))
 241     lower <- c(lower, rep(-Inf, length(truth.params)))
 242     start <- rep(0.1,length(params))
 243     names(start) <- params
 244
 245     if(method=='optim'){
 246         fit <- optim(start, fn = measerr_mle_nll, lower=lower, method='L-BFGS-B', hessian=TRUE, control=list(maxit=1e6))
 247     } else { # method='mle2'
 248
 249         quoted.names <- gsub("[\\(\\)]",'',names(start))
 250
 251         text <- paste("function(", paste0(quoted.names,'=',start,collapse=','),"){params<-c(",paste0(quoted.names,collapse=','),");return(measerr_mle_nll(params))}")
 252
 253         measerr_mle_nll_mle <- eval(parse(text=text))
 254         names(start) <- quoted.names
 255         names(lower) <- quoted.names
 256         fit <- mle2(minuslogl=measerr_mle_nll_mle, start=start, lower=lower, parnames=params,control=list(maxit=1e6),method='L-BFGS-B')
 257     }
 258
 259     return(fit)
 260 }
 261
 262 ## Experimental, but probably works.
 263 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'){
 264
 265     ### in this scenario, the ground truth also has measurement error, but we have repeated measures for it.
 266     # this time we never get to observe the true X
 267     outcome.model.matrix <- model.matrix(outcome_formula, df)
 268     proxy.model.matrix <- model.matrix(proxy_formula, df)
 269     response.var <- all.vars(outcome_formula)[1]
 270     proxy.var <- all.vars(proxy_formula)[1]
 271     param.var <- all.vars(truth_formula)[1]
 272     truth.var<- all.vars(truth_formula)[1]
 273     y <- with(df,eval(parse(text=response.var)))
 274
 275     nll <- function(params){
 276         param.idx <- 1
 277         n.outcome.model.covars <- dim(outcome.model.matrix)[2]
 278         outcome.params <- params[param.idx:n.outcome.model.covars]
 279         param.idx <- param.idx + n.outcome.model.covars
 280
 281
 282         if(outcome_family$family == "gaussian"){
 283             sigma.y <- params[param.idx]
 284             param.idx <- param.idx + 1
 285         }
 286
 287
 288         n.proxy.model.covars <- dim(proxy.model.matrix)[2]
 289         proxy.params <- params[param.idx:(n.proxy.model.covars+param.idx-1)]
 290         param.idx <- param.idx + n.proxy.model.covars
 291
 292         df.temp <- copy(df)
 293
 294         if((truth_family$family == "binomial")
 295            & (truth_family$link=='logit')){
 296             integrate.grid <- expand.grid(replicate(1 + length(coder_formulas), c(0,1), simplify=FALSE))
 297             ll.parts <- matrix(nrow=nrow(df),ncol=nrow(integrate.grid))
 298             for(i in 1:nrow(integrate.grid)){
 299                 # setup the dataframe for this row
 300                 row <- integrate.grid[i,]
 301
 302                 df.temp[[param.var]] <- row[[1]]
 303                 ci <- 2
 304                 for(coder_formula in coder_formulas){
 305                     coder.var <- all.vars(coder_formula)[1]
 306                     df.temp[[coder.var]] <- row[[ci]]
 307                     ci <- ci + 1
 308                 }
 309
 310                 outcome.model.matrix <- model.matrix(outcome_formula, df.temp)
 311                 if(outcome_family$family == "gaussian"){
 312                     ll.y <- dnorm(y, outcome.params %*% t(outcome.model.matrix), sd=sigma.y, log=TRUE)
 313                 }
 314
 315                 if(proxy_family$family=="binomial" & (proxy_family$link=='logit')){
 316                     proxy.model.matrix <- model.matrix(proxy_formula, df.temp)
 317                     ll.w <- vector(mode='numeric', length=dim(proxy.model.matrix)[1])
 318                     proxyvar <- with(df.temp,eval(parse(text=proxy.var)))
 319                     ll.w[proxyvar==1] <- plogis(proxy.params %*% t(proxy.model.matrix[proxyvar==1,]),log=TRUE)
 320                     ll.w[proxyvar==0] <- plogis(proxy.params %*% t(proxy.model.matrix[proxyvar==0,]),log=TRUE,lower.tail=FALSE)
 321                 }
 322
 323                 ## probability of the coded variables
 324                 coder.lls <- matrix(nrow=nrow(df.temp),ncol=length(coder_formulas))
 325                 ci <- 1
 326                 for(coder_formula in coder_formulas){
 327                     coder.model.matrix <- model.matrix(coder_formula, df.temp)
 328                     n.coder.model.covars <- dim(coder.model.matrix)[2]
 329                     coder.params <- params[param.idx:(n.coder.model.covars + param.idx - 1)]
 330                     param.idx <- param.idx + n.coder.model.covars
 331                     coder.var <- all.vars(coder_formula)[1]
 332                     x.obs <- with(df.temp, eval(parse(text=coder.var)))
 333                     true.codervar <- df[[all.vars(coder_formula)[1]]]
 334
 335                     ll.coder <- vector(mode='numeric', length=dim(coder.model.matrix)[1])
 336                     ll.coder[x.obs==1] <- plogis(coder.params %*% t(coder.model.matrix[x.obs==1,]),log=TRUE)
 337                     ll.coder[x.obs==0] <- plogis(coder.params %*% t(coder.model.matrix[x.obs==0,]),log=TRUE,lower.tail=FALSE)
 338
 339                     # don't count when we know the observed value, unless we're accounting for observed value
 340                     ll.coder[(!is.na(true.codervar)) & (true.codervar != x.obs)] <- NA
 341                     coder.lls[,ci] <- ll.coder
 342                     ci <- ci + 1
 343                 }
 344
 345                 truth.model.matrix <- model.matrix(truth_formula, df.temp)
 346                 n.truth.model.covars <- dim(truth.model.matrix)[2]
 347                 truth.params <- params[param.idx:(n.truth.model.covars + param.idx - 1)]
 348
 349                 for(coder_formula in coder_formulas){
 350                     coder.model.matrix <- model.matrix(coder_formula, df.temp)
 351                     n.coder.model.covars <- dim(coder.model.matrix)[2]
 352                     param.idx <- param.idx - n.coder.model.covars
 353                 }
 354
 355                 x <- with(df.temp, eval(parse(text=truth.var)))
 356                 ll.truth <- vector(mode='numeric', length=dim(truth.model.matrix)[1])
 357                 ll.truth[x==1] <- plogis(truth.params %*% t(truth.model.matrix[x==1,]), log=TRUE)
 358                 ll.truth[x==0] <- plogis(truth.params %*% t(truth.model.matrix[x==0,]), log=TRUE, lower.tail=FALSE)
 359
 360                 true.truthvar <- df[[all.vars(truth_formula)[1]]]
 361
 362                 if(!is.null(true.truthvar)){
 363                                         # ll.truth[(!is.na(true.truthvar)) & (true.truthvar != truthvar)] <- -Inf
 364                     # ll.truth[(!is.na(true.truthvar)) & (true.truthvar == truthvar)] <- 0
 365                 }
 366                 ll.parts[,i] <- ll.y + ll.w + apply(coder.lls,1,sum) + ll.truth
 367
 368             }
 369
 370             lls <- rowLogSumExps(ll.parts,na.rm=TRUE)
 371
 372         ## likelihood of observed data
 373             target <- -1 * sum(lls)
 374             return(target)
 375         }
 376     }
 377
 378     outcome.params <- colnames(model.matrix(outcome_formula,df))
 379     lower <- rep(-Inf, length(outcome.params))
 380
 381     if(outcome_family$family=='gaussian'){
 382         params <- c(outcome.params, 'sigma_y')
 383         lower <- c(lower, 0.00001)
 384     } else {
 385         params <- outcome.params
 386     }
 387
 388     proxy.params <- colnames(model.matrix(proxy_formula, df))
 389     params <- c(params, paste0('proxy_',proxy.params))
 390     positive.params <- paste0('proxy_',truth.var)
 391     lower <- c(lower, rep(-Inf, length(proxy.params)))
 392     names(lower) <- params
 393     lower[positive.params] <- 0.01
 394     ci <- 0
 395
 396     for(coder_formula in coder_formulas){
 397         coder.params <- colnames(model.matrix(coder_formula,df))
 398         params <- c(params, paste0('coder_',ci,coder.params))
 399         positive.params <- paste0('coder_', ci, truth.var)
 400         ci <- ci + 1
 401         lower <- c(lower, rep(-Inf, length(coder.params)))
 402         names(lower) <- params
 403         lower[positive.params] <- 0.01
 404     }
 405
 406     truth.params <- colnames(model.matrix(truth_formula, df))
 407     params <- c(params, paste0('truth_', truth.params))
 408     lower <- c(lower, rep(-Inf, length(truth.params)))
 409     start <- rep(0.1,length(params))
 410     names(start) <- params
 411     names(lower) <- params
 412
 413     if(method=='optim'){
 414         print(start)
 415         fit <- optim(start, fn = nll, lower=lower, method='L-BFGS-B', hessian=TRUE, control=list(maxit=1e6))
 416     } else {
 417
 418         quoted.names <- gsub("[\\(\\)]",'',names(start))
 419         print(quoted.names)
 420         text <- paste("function(", paste0(quoted.names,'=',start,collapse=','),"){params<-c(",paste0(quoted.names,collapse=','),");return(nll(params))}")
 421
 422         measerr_mle_nll <- eval(parse(text=text))
 423         names(start) <- quoted.names
 424         names(lower) <- quoted.names
 425         fit <- mle2(minuslogl=measerr_mle_nll, start=start, lower=lower, method='L-BFGS-B',control=list(maxit=1e6))
 426     }
 427
 428     return(fit)
 429 }
 430
 431 ## Experimental, and does not work.
 432 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'){
 433     integrate.grid <- expand.grid(replicate(1 + length(coder_formulas), c(0,1), simplify=FALSE))
 434     print(integrate.grid)
 435
 436
 437     outcome.model.matrix <- model.matrix(outcome_formula, df)
 438     n.outcome.model.covars <- dim(outcome.model.matrix)[2]
 439
 440
 441     ### in this scenario, the ground truth also has measurement error, but we have repeated measures for it.
 442     # this time we never get to observe the true X
 443     nll <- function(params){
 444         param.idx <- 1
 445         outcome.params <- params[param.idx:n.outcome.model.covars]
 446         param.idx <- param.idx + n.outcome.model.covars
 447         proxy.model.matrix <- model.matrix(proxy_formula, df)
 448         n.proxy.model.covars <- dim(proxy.model.matrix)[2]
 449         response.var <- all.vars(outcome_formula)[1]
 450
 451         if(outcome_family$family == "gaussian"){
 452             sigma.y <- params[param.idx]
 453             param.idx <- param.idx + 1
 454         }
 455
 456         proxy.params <- params[param.idx:(n.proxy.model.covars+param.idx-1)]
 457         param.idx <- param.idx + n.proxy.model.covars
 458
 459         df.temp <- copy(df)
 460
 461         if((outcome_family$family == "binomial")
 462            & (outcome_family$link=='logit')){
 463             ll.parts <- matrix(nrow=nrow(df),ncol=nrow(integrate.grid))
 464             for(i in 1:nrow(integrate.grid)){
 465                 # setup the dataframe for this row
 466                 row <- integrate.grid[i,]
 467
 468                 df.temp[[response.var]] <- row[[1]]
 469                 ci <- 2
 470                 for(coder_formula in coder_formulas){
 471                     codervar <- all.vars(coder_formula)[1]
 472                     df.temp[[codervar]] <- row[[ci]]
 473                     ci <- ci + 1
 474                 }
 475
 476                 outcome.model.matrix <- model.matrix(outcome_formula, df.temp)
 477                 if(outcome_family$family == "gaussian"){
 478                     ll.y <- dnorm(df.temp[[response.var]], outcome.params %*% t(outcome.model.matrix), sd=sigma.y, log=T)
 479                 }
 480
 481                 if(outcome_family$family == "binomial" & (outcome_family$link=='logit')){
 482                     ll.y <- vector(mode='numeric',length=nrow(df.temp))
 483                     ll.y[df.temp[[response.var]]==1] <- plogis( outcome.params %*% t(outcome.model.matrix), log=TRUE)
 484                     ll.y[df.temp[[response.var]]==0] <- plogis( outcome.params %*% t(outcome.model.matrix), log=TRUE,lower.tail=FALSE)
 485                 }
 486
 487                 if(proxy_family$family=="binomial" & (proxy_family$link=='logit')){
 488                     proxy.model.matrix <- model.matrix(proxy_formula, df.temp)
 489                     ll.w <- vector(mode='numeric', length=dim(proxy.model.matrix)[1])
 490                     proxyvar <- with(df.temp,eval(parse(text=all.vars(proxy_formula)[1])))
 491                     ll.w[proxyvar==1] <- plogis(proxy.params %*% t(proxy.model.matrix[proxyvar==1,]),log=TRUE)
 492                     ll.w[proxyvar==0] <- plogis(proxy.params %*% t(proxy.model.matrix[proxyvar==0,]),log=TRUE,lower.tail=FALSE)
 493                 }
 494
 495                 ## probability of the coded variables
 496                 coder.lls <- matrix(nrow=nrow(df.temp),ncol=length(coder_formulas))
 497                 ci <- 1
 498                 for(coder_formula in coder_formulas){
 499                     coder.model.matrix <- model.matrix(coder_formula, df.temp)
 500                     n.coder.model.covars <- dim(coder.model.matrix)[2]
 501                     coder.params <- params[param.idx:(n.coder.model.covars + param.idx - 1)]
 502                     param.idx <- param.idx + n.coder.model.covars
 503                     codervar <- with(df.temp, eval(parse(text=all.vars(coder_formula)[1])))
 504                     true.codervar <- df[[all.vars(coder_formula)[1]]]
 505
 506                     ll.coder <- vector(mode='numeric', length=dim(coder.model.matrix)[1])
 507                     ll.coder[codervar==1] <- plogis(coder.params %*% t(coder.model.matrix[codervar==1,]),log=TRUE)
 508                     ll.coder[codervar==0] <- plogis(coder.params %*% t(coder.model.matrix[codervar==0,]),log=TRUE,lower.tail=FALSE)
 509
 510                     # don't count when we know the observed value, unless we're accounting for observed value
 511                     ll.coder[(!is.na(true.codervar)) & (true.codervar != codervar)] <- NA
 512                     coder.lls[,ci] <- ll.coder
 513                     ci <- ci + 1
 514                 }
 515
 516                 for(coder_formula in coder_formulas){
 517                     coder.model.matrix <- model.matrix(coder_formula, df.temp)
 518                     n.coder.model.covars <- dim(coder.model.matrix)[2]
 519                     param.idx <- param.idx - n.coder.model.covars
 520                 }
 521
 522                 ll.parts[,i] <- ll.y + ll.w + apply(coder.lls,1,function(x) sum(x))
 523
 524             }
 525
 526             lls <- rowLogSumExps(ll.parts,na.rm=TRUE)
 527
 528             ## likelihood of observed data
 529             target <- -1 * sum(lls)
 530             print(target)
 531             print(params)
 532             return(target)
 533         }
 534     }
 535
 536     outcome.params <- colnames(model.matrix(outcome_formula,df))
 537     response.var <- all.vars(outcome_formula)[1]
 538     lower <- rep(-Inf, length(outcome.params))
 539
 540     if(outcome_family$family=='gaussian'){
 541         params <- c(outcome.params, 'sigma_y')
 542         lower <- c(lower, 0.00001)
 543     } else {
 544         params <- outcome.params
 545     }
 546
 547     ## constrain the model of the coder and proxy vars
 548     ## this is to ensure identifiability
 549     ## it is a safe assumption because the coders aren't hostile (wrong more often than right)
 550     ## so we can assume that y ~Bw, B is positive
 551     proxy.params <- colnames(model.matrix(proxy_formula, df))
 552     positive.params <- paste0('proxy_',response.var)
 553     params <- c(params, paste0('proxy_',proxy.params))
 554     lower <- c(lower, rep(-Inf, length(proxy.params)))
 555     names(lower) <- params
 556     lower[positive.params] <- 0.001
 557
 558     ci <- 0
 559     for(coder_formula in coder_formulas){
 560         coder.params <- colnames(model.matrix(coder_formula,df))
 561         latent.coder.params <- coder.params %in% response.var
 562         params <- c(params, paste0('coder_',ci,coder.params))
 563         positive.params <- paste0('coder_',ci,response.var)
 564         ci <- ci + 1
 565         lower <- c(lower, rep(-Inf, length(coder.params)))
 566         names(lower) <-params
 567         lower[positive.params] <- 0.001
 568     }
 569
 570     ## init by using the "loco model"
 571     temp.df <- copy(df)
 572     temp.df <- temp.df[y.obs.1 == y.obs.0, y:=y.obs.1]
 573     loco.model <- glm(outcome_formula, temp.df, family=outcome_family)
 574
 575     start <- rep(1,length(params))
 576     names(start) <- params
 577     start[names(coef(loco.model))] <- coef(loco.model)
 578     names(lower) <- params
 579     if(method=='optim'){
 580         print(lower)
 581         fit <- optim(start, fn = nll, lower=lower, method='L-BFGS-B', hessian=TRUE,control=list(maxit=1e6))
 582     } else {
 583
 584         quoted.names <- gsub("[\\(\\)]",'',names(start))
 585         print(quoted.names)
 586         text <- paste("function(", paste0(quoted.names,'=',start,collapse=','),"){params<-c(",paste0(quoted.names,collapse=','),");return(nll(params))}")
 587
 588         measerr_mle_nll <- eval(parse(text=text))
 589         names(start) <- quoted.names
 590         names(lower) <- quoted.names
 591         fit <- mle2(minuslogl=measerr_mle_nll, start=start, lower=lower, parnames=params,control=list(maxit=1e6),method='L-BFGS-B')
 592     }
 593
 594     return(fit)
 595 }
 596