changes from klone

[ml_measurement_error_public.git] / simulations / measerr_methods.R

diff --git a/simulations/measerr_methods.R b/simulations/measerr_methods.R
index 63f8bc19df116b2f6d148c55cefbd89cf2dc7b85..8f8eb5db082ff140885aee2cdd467258d2fddfbd 100644 (file)
--- a/simulations/measerr_methods.R
+++ b/simulations/measerr_methods.R
@@ -15,18 +15,38 @@ library(bbmle)
 
 ### ideal formulas for example 2
 # test.fit.2 <- measerr_mle(df, y ~ x + z, gaussian(), w_pred ~ x + z + y + y:x, binomial(link='logit'), x ~ z)
 
 ### ideal formulas for example 2
 # test.fit.2 <- measerr_mle(df, y ~ x + z, gaussian(), w_pred ~ x + z + y + y:x, binomial(link='logit'), x ~ z)

-
+likelihood.logistic <- function(model.params, outcome, model.matrix){
+    ll <- vector(mode='numeric', length=length(outcome))
+    ll[outcome == 1] <- plogis(model.params %*% t(model.matrix[outcome==1,]), log=TRUE)
+    ll[outcome == 0] <- plogis(model.params %*% t(model.matrix[outcome==0,]), log=TRUE, lower.tail=FALSE)
+    return(ll)
+}

 
 ## outcome_formula <- y ~ x + z; proxy_formula <- w_pred ~ y + x + z + x:z + x:y + z:y 
 
 ## 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'),method='optim'){
+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'){
+    df.obs <- model.frame(outcome_formula, df)
+    proxy.model.matrix <- model.matrix(proxy_formula, df)
+    proxy.variable <- all.vars(proxy_formula)[1]
+
+    df.proxy.obs <- model.frame(proxy_formula,df)
+    proxy.obs <- with(df.proxy.obs, eval(parse(text=proxy.variable)))
+
+    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.obs)
+
+    df.unobs <- df[is.na(df[[response.var]])]
+    df.unobs.y1 <- copy(df.unobs)
+    df.unobs.y1[[response.var]] <- 1
+    df.unobs.y0 <- copy(df.unobs)
+    df.unobs.y0[[response.var]] <- 0
+
+    outcome.model.matrix.y1 <- model.matrix(outcome_formula, df.unobs.y1)
+    proxy.model.matrix.y1 <- model.matrix(proxy_formula, df.unobs.y1)
+    proxy.model.matrix.y0 <- model.matrix(proxy_formula, df.unobs.y0)
+    proxy.unobs <- with(df.unobs, eval(parse(text=proxy.variable)))

 
     nll <- function(params){
 
     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.obs)

 
         param.idx <- 1
         n.outcome.model.covars <- dim(outcome.model.matrix)[2]
 
         param.idx <- 1
         n.outcome.model.covars <- dim(outcome.model.matrix)[2]


@@ -39,12 +59,9 @@ measerr_mle_dv <- function(df, outcome_formula, outcome_family=binomial(link='lo
             ll.y.obs[y.obs==0] <- plogis(outcome.params %*% t(outcome.model.matrix[y.obs==0,]),log=TRUE,lower.tail=FALSE)
         }
 
             ll.y.obs[y.obs==0] <- plogis(outcome.params %*% t(outcome.model.matrix[y.obs==0,]),log=TRUE,lower.tail=FALSE)
         }
 

-        df.obs <- model.frame(proxy_formula,df)

         n.proxy.model.covars <- dim(proxy.model.matrix)[2]
         proxy.params <- params[param.idx:(n.proxy.model.covars+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
         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 <- vector(mode='numeric',length=dim(proxy.model.matrix)[1])
 
         if( (proxy_family$family=="binomial") & (proxy_family$link=='logit')){
             ll.w.obs <- vector(mode='numeric',length=dim(proxy.model.matrix)[1])


@@ -53,15 +70,8 @@ measerr_mle_dv <- function(df, outcome_formula, outcome_family=binomial(link='lo
         }
 
         ll.obs <- sum(ll.y.obs + ll.w.obs)
         }
 
         ll.obs <- sum(ll.y.obs + ll.w.obs)

-
-        df.unobs <- df[is.na(df[[response.var]])]
-        df.unobs.y1 <- copy(df.unobs)
-        df.unobs.y1[[response.var]] <- 1
-        df.unobs.y0 <- copy(df.unobs)
-        df.unobs.y0[[response.var]] <- 0

         
         ## integrate out y
         
         ## integrate out y

-        outcome.model.matrix.y1 <- model.matrix(outcome_formula, df.unobs.y1)

 
         if((outcome_family$family == "binomial") & (outcome_family$link == 'logit')){
             ll.y.unobs.1 <- vector(mode='numeric', length=dim(outcome.model.matrix.y1)[1])
 
         if((outcome_family$family == "binomial") & (outcome_family$link == 'logit')){
             ll.y.unobs.1 <- vector(mode='numeric', length=dim(outcome.model.matrix.y1)[1])


@@ -70,10 +80,6 @@ measerr_mle_dv <- function(df, outcome_formula, outcome_family=binomial(link='lo
             ll.y.unobs.0 <- plogis(outcome.params %*% t(outcome.model.matrix.y1),log=TRUE,lower.tail=FALSE)
         }
 
             ll.y.unobs.0 <- plogis(outcome.params %*% t(outcome.model.matrix.y1),log=TRUE,lower.tail=FALSE)
         }
 

-        proxy.model.matrix.y1 <- model.matrix(proxy_formula, df.unobs.y1)
-        proxy.model.matrix.y0 <- model.matrix(proxy_formula, df.unobs.y0)
-        proxy.unobs <- with(df.unobs, eval(parse(text=proxy.variable)))
-

         if( (proxy_family$family=="binomial") & (proxy_family$link=='logit')){
             ll.w.unobs.1 <- vector(mode='numeric',length=dim(proxy.model.matrix.y1)[1])
             ll.w.unobs.0 <- vector(mode='numeric',length=dim(proxy.model.matrix.y0)[1])
         if( (proxy_family$family=="binomial") & (proxy_family$link=='logit')){
             ll.w.unobs.1 <- vector(mode='numeric',length=dim(proxy.model.matrix.y1)[1])
             ll.w.unobs.0 <- vector(mode='numeric',length=dim(proxy.model.matrix.y0)[1])


@@ -100,7 +106,7 @@ measerr_mle_dv <- function(df, outcome_formula, outcome_family=binomial(link='lo
     names(start) <- params
     
     if(method=='optim'){
     names(start) <- params
     
     if(method=='optim'){

-        fit <- optim(start, fn = nll, lower=lower, method='L-BFGS-B', hessian=TRUE, control=list(maxit=1e6))
+        fit <- optim(start, fn = nll, lower=lower, method=optim_method, hessian=TRUE, control=list(maxit=maxit))

     } else {
         quoted.names <- gsub("[\\(\\)]",'',names(start))
         print(quoted.names)
     } else {
         quoted.names <- gsub("[\\(\\)]",'',names(start))
         print(quoted.names)


@@ -109,13 +115,13 @@ measerr_mle_dv <- function(df, outcome_formula, outcome_family=binomial(link='lo
         measerr_mle_nll <- eval(parse(text=text))
         names(start) <- quoted.names
         names(lower) <- quoted.names
         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')
+        fit <- mle2(minuslogl=measerr_mle_nll, start=start, lower=lower, parnames=params,control=list(maxit=maxit),method=optim_method)

     }
     return(fit)
 }
 
 
     }
     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'),method='optim'){
+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'){

 
     df.obs <- model.frame(outcome_formula, df)
     response.var <- all.vars(outcome_formula)[1]
 
     df.obs <- model.frame(outcome_formula, df)
     response.var <- all.vars(outcome_formula)[1]


@@ -125,6 +131,31 @@ measerr_mle <- function(df, outcome_formula, outcome_family=gaussian(), proxy_fo
     proxy.model.matrix <- model.matrix(proxy_formula, df)
     y.obs <- with(df.obs,eval(parse(text=response.var)))
 
     proxy.model.matrix <- model.matrix(proxy_formula, df)
     y.obs <- with(df.obs,eval(parse(text=response.var)))
 

+    df.proxy.obs <- model.frame(proxy_formula,df)
+    proxy.obs <- with(df.proxy.obs, eval(parse(text=proxy.variable)))
+    n.proxy.model.covars <- dim(proxy.model.matrix)[2]
+
+    df.truth.obs <- model.frame(truth_formula, df)
+    truth.obs <- with(df.truth.obs, eval(parse(text=truth.variable)))
+    truth.model.matrix <- model.matrix(truth_formula,df.truth.obs)
+    n.truth.model.covars <- dim(truth.model.matrix)[2]
+
+    df.unobs <- df[is.na(eval(parse(text=truth.variable)))]
+    df.unobs.x1 <- copy(df.unobs)
+    df.unobs.x1[,truth.variable] <- 1
+    df.unobs.x0 <- copy(df.unobs)
+    df.unobs.x0[,truth.variable] <- 0
+    outcome.unobs <- with(df.unobs, eval(parse(text=response.var)))
+            
+    outcome.model.matrix.x0 <- model.matrix(outcome_formula, df.unobs.x0)
+    outcome.model.matrix.x1 <- model.matrix(outcome_formula, df.unobs.x1)
+        
+    proxy.model.matrix.x0 <- model.matrix(proxy_formula, df.unobs.x0)
+    proxy.model.matrix.x1 <- model.matrix(proxy_formula, df.unobs.x1)
+    proxy.unobs <- df.unobs[[proxy.variable]]
+
+    truth.model.matrix.unobs <- model.matrix(truth_formula, df.unobs.x0)
+

     measerr_mle_nll <- function(params){
         param.idx <- 1
         n.outcome.model.covars <- dim(outcome.model.matrix)[2]
     measerr_mle_nll <- function(params){
         param.idx <- 1
         n.outcome.model.covars <- dim(outcome.model.matrix)[2]


@@ -137,82 +168,48 @@ measerr_mle <- function(df, outcome_formula, outcome_family=gaussian(), proxy_fo
             param.idx <- param.idx + 1
                                         #  outcome_formula likelihood using linear regression
             ll.y.obs <- dnorm(y.obs, outcome.params %*% t(outcome.model.matrix),sd=sigma.y, log=TRUE)
             param.idx <- param.idx + 1
                                         #  outcome_formula likelihood using linear regression
             ll.y.obs <- dnorm(y.obs, outcome.params %*% t(outcome.model.matrix),sd=sigma.y, log=TRUE)

-        }
-        
-        df.obs <- model.frame(proxy_formula,df)
-        n.proxy.model.covars <- dim(proxy.model.matrix)[2]
+        } else if( (outcome_family$family == "binomial") & (outcome_family$link == "logit") )
+            ll.y.obs <- likelihood.logistic(outcome.params, y.obs, outcome.model.matrix)
+
+    

         proxy.params <- params[param.idx:(n.proxy.model.covars+param.idx-1)]
         param.idx <- param.idx + n.proxy.model.covars
         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 <- vector(mode='numeric',length=dim(proxy.model.matrix)[1])

 
 

-                                        # 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)
-        }
+        if( (proxy_family$family=="binomial") & (proxy_family$link=='logit'))
+            ll.w.obs <- likelihood.logistic(proxy.params, proxy.obs, proxy.model.matrix)

 
 

-        df.obs <- model.frame(truth_formula, df)
-
-        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]
-        

         truth.params <- params[param.idx:(n.truth.model.covars + param.idx - 1)]
 
         truth.params <- params[param.idx:(n.truth.model.covars + param.idx - 1)]
 

-        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
-            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)
-        }
+        if( (truth_family$family=="binomial") & (truth_family$link=='logit'))
+            ll.x.obs <- likelihood.logistic(truth.params, truth.obs, truth.model.matrix)

         
         

-                                        # 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
         ## integrate out the "truth" variable. 
         
         if(truth_family$family=='binomial'){
         ll.obs <- sum(ll.y.obs + ll.w.obs + ll.x.obs)
 
         ## likelihood for the predicted data
         ## integrate out the "truth" variable. 
         
         if(truth_family$family=='binomial'){

-            df.unobs <- df[is.na(eval(parse(text=truth.variable)))]
-            df.unobs.x1 <- copy(df.unobs)
-            df.unobs.x1[,'x'] <- 1
-            df.unobs.x0 <- copy(df.unobs)
-            df.unobs.x0[,'x'] <- 0
-            outcome.unobs <- with(df.unobs, eval(parse(text=response.var)))
-            
-            outcome.model.matrix.x0 <- model.matrix(outcome_formula, df.unobs.x0)
-            outcome.model.matrix.x1 <- model.matrix(outcome_formula, df.unobs.x1)

             if(outcome_family$family=="gaussian"){
             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)

+            } else if( (outcome_family$family == "binomial") & (outcome_family$link == "logit") ){
+                ll.y.x1 <- likelihood.logistic(outcome.params, outcome.unobs, outcome.model.matrix.x1)
+                ll.y.x0 <- likelihood.logistic(outcome.params, outcome.unobs, outcome.model.matrix.x0)

             }
 
             if( (proxy_family$family=='binomial') & (proxy_family$link=='logit')){
 
             }
 
             if( (proxy_family$family=='binomial') & (proxy_family$link=='logit')){
 

-                proxy.model.matrix.x0 <- model.matrix(proxy_formula, df.unobs.x0)
-                proxy.model.matrix.x1 <- model.matrix(proxy_formula, df.unobs.x1)
-                proxy.unobs <- df.unobs[[proxy.variable]]
-                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
-                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)
+                ll.w.x0 <- likelihood.logistic(proxy.params, proxy.unobs, proxy.model.matrix.x0)
+                ll.w.x1 <- likelihood.logistic(proxy.params, proxy.unobs, proxy.model.matrix.x1)

 
 

-                ll.w.x0[proxy.unobs==0] <- plogis(proxy.params %*% t(proxy.model.matrix.x0[proxy.unobs==0,]), log=TRUE,lower.tail=FALSE)
-                ll.w.x1[proxy.unobs==0] <- plogis(proxy.params %*% t(proxy.model.matrix.x1[proxy.unobs==0,]), log=TRUE,lower.tail=FALSE)

             }
 
             if(truth_family$link=='logit'){
             }
 
             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)
+                ll.x.x1 <- plogis(truth.params %*% t(truth.model.matrix.unobs), log=TRUE)
+                ll.x.x0 <- plogis(truth.params %*% t(truth.model.matrix.unobs), log=TRUE, lower.tail=FALSE)

             }
         }
 
             }
         }
 


@@ -243,7 +240,7 @@ measerr_mle <- function(df, outcome_formula, outcome_family=gaussian(), proxy_fo
     names(start) <- params
 
     if(method=='optim'){
     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))
+        fit <- optim(start, fn = measerr_mle_nll, lower=lower, method=optim_method, hessian=TRUE, control=list(maxit=maxit))

     } else { # method='mle2'
                 
         quoted.names <- gsub("[\\(\\)]",'',names(start))
     } else { # method='mle2'
                 
         quoted.names <- gsub("[\\(\\)]",'',names(start))


@@ -253,7 +250,7 @@ measerr_mle <- function(df, outcome_formula, outcome_family=gaussian(), proxy_fo
         measerr_mle_nll_mle <- eval(parse(text=text))
         names(start) <- quoted.names
         names(lower) <- quoted.names
         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')
+        fit <- mle2(minuslogl=measerr_mle_nll_mle, start=start, lower=lower, parnames=params,control=list(maxit=maxit),method=optim_method)

     }
 
     return(fit)
     }
 
     return(fit)


@@ -431,7 +428,7 @@ measerr_irr_mle <- function(df, outcome_formula, outcome_family=gaussian(), code
 ## 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))
 ## 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)
+#    print(integrate.grid)

 
 
     outcome.model.matrix <- model.matrix(outcome_formula, df)
 
 
     outcome.model.matrix <- model.matrix(outcome_formula, df)


@@ -527,8 +524,8 @@ measerr_irr_mle_dv <- function(df, outcome_formula, outcome_family=binomial(link
 
             ## likelihood of observed data 
             target <- -1 * sum(lls)
 
             ## likelihood of observed data 
             target <- -1 * sum(lls)

-            print(target)
-            print(params)
+#            print(target)
+#            print(params)

             return(target)
         }
     }
             return(target)
         }
     }