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