X-Git-Url: https://code.communitydata.science/ml_measurement_error_public.git/blobdiff_plain/979dc14b6861ae31f00d56392fd5b8cf69f17333..b52b4f7daaba8a877b041ddb24c8f36b466ddc5b:/simulations/measerr_methods.R

diff --git a/simulations/measerr_methods.R b/simulations/measerr_methods.R
index 00f1746..087c608 100644
--- a/simulations/measerr_methods.R
+++ b/simulations/measerr_methods.R
@@ -1,6 +1,6 @@
 library(formula.tools)
 library(matrixStats)
-
+library(bbmle)
 ## df: dataframe to model
 ## outcome_formula: formula for y | x, z
 ## outcome_family: family for y | x, z
@@ -17,7 +17,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,12 +98,23 @@ 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)
 }
 
 ## 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')){
+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'),method='optim'){
 
     ### in this scenario, the ground truth also has measurement error, but we have repeated measures for it. 
 
@@ -293,14 +304,28 @@ measerr_irr_mle <- function(df, outcome_formula, outcome_family=gaussian(), rate
     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')){
+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'){
 
-    measrr_mle_nll <- function(params){
+    measerr_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)
@@ -425,8 +450,21 @@ 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
-    
-    fit <- optim(start, fn = measrr_mle_nll, lower=lower, method='L-BFGS-B', hessian=TRUE, control=list(maxit=1e6))
+
+    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)
 }
+