### 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
-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){
- 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]
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)]
-
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])
}
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
- 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])
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])
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)
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)
}
-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]
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]
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.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)]
- 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'){
- 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"){
-
# 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')){
- 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'){
- truth.model.matrix <- model.matrix(truth_formula, df.unobs.x0)
# 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)
}
}
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))
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)
## 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)
## likelihood of observed data
target <- -1 * sum(lls)
- print(target)
- print(params)
+# print(target)
+# print(params)
return(target)
}
}
projects / ml_measurement_error_public.git / blobdiff