X-Git-Url: https://code.communitydata.science/ml_measurement_error_public.git/blobdiff_plain/2cd447c327744263d5f94b20e1146cdf31b2ec2c..47e9367ed5c61b721bdc17cddd76bced4f8ed621:/simulations/measerr_methods.R

diff --git a/simulations/measerr_methods.R b/simulations/measerr_methods.R
index ab87d71..6bf8c3f 100644
--- a/simulations/measerr_methods.R
+++ b/simulations/measerr_methods.R
@@ -57,7 +57,7 @@ measerr_mle_dv <- function(df, outcome_formula, outcome_family=binomial(link='lo
         df.unobs.y1 <- copy(df.unobs)
         df.unobs.y1[[response.var]] <- 1
         df.unobs.y0 <- copy(df.unobs)
-        df.unobs.y0[[response.var]] <- 1
+        df.unobs.y0[[response.var]] <- 0
         
         ## integrate out y
         outcome.model.matrix.y1 <- model.matrix(outcome_formula, df.unobs.y1)
@@ -124,6 +124,8 @@ 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
             ll.y.obs <- dnorm(y.obs, outcome.params %*% t(outcome.model.matrix),sd=sigma.y, log=TRUE)
         }
         
@@ -135,6 +137,8 @@ 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
             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)
         }
@@ -149,10 +153,13 @@ 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
             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
         ll.obs <- sum(ll.y.obs + ll.w.obs + ll.x.obs)
 
         ## likelihood for the predicted data
@@ -169,6 +176,8 @@ measerr_mle <- function(df, outcome_formula, outcome_family=gaussian(), proxy_fo
             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"){
+
+                # 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)
             }
@@ -181,6 +190,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
                 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)
 
@@ -190,8 +200,9 @@ 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)
-                ll.x.x0 <- plogis(truth.params %*% t(truth.model.matrix), log=TRUE)
-                ll.x.x1 <- plogis(truth.params %*% t(truth.model.matrix), log=TRUE, lower.tail=FALSE)
+                # 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)
             }
         }