119 lines
4.6 KiB
R
119 lines
4.6 KiB
R
options(amelia.parallel="no",
|
|
amelia.ncpus=1)
|
|
library(matrixStats) # for numerically stable logsumexps
|
|
source("../misclassificationmodels/R/likelihoods.R")
|
|
source("../misclassificationmodels/R/glm_fixit.R")
|
|
|
|
## ... are the parts of the formular.
|
|
run_simulation <- function(..., family=gaussian(),
|
|
proxy_family = binomial(link='logit'),
|
|
truth_family = binomial(link='logit'),
|
|
data,
|
|
data2,
|
|
result
|
|
){
|
|
|
|
accuracy <- df[,mean(w_pred==x)]
|
|
accuracy.x.y0 <- df[y<=0,mean(w_pred==x)]
|
|
accuracy.x.y1 <- df[y>=0,mean(w_pred==x)]
|
|
accuracy.z.y0 <- df[y<=0,mean(a_pred==z)]
|
|
accuracy.z.y1 <- df[y>=0,mean(a_pred==z)]
|
|
cor.y.xi <- cor(df$x - df$w_pred, df$y)
|
|
cor.y.zi <- cor(df$z - df$a_pred, df$y)
|
|
|
|
fnr <- df[w_pred==0,mean(w_pred!=x)]
|
|
fnr.y0 <- df[(w_pred==0) & (y<=0),mean(w_pred!=x)]
|
|
fnr.y1 <- df[(w_pred==0) & (y>=0),mean(w_pred!=x)]
|
|
|
|
fpr <- df[w_pred==1,mean(w_pred!=x)]
|
|
fpr.y0 <- df[(w_pred==1) & (y<=0),mean(w_pred!=x)]
|
|
fpr.y1 <- df[(w_pred==1) & (y>=0),mean(w_pred!=x)]
|
|
cor.resid.w_pred <- cor(resid(lm(y~x+z,df)),df$w_pred)
|
|
|
|
result <- append(result, list(accuracy=accuracy,
|
|
accuracy.x.y0=accuracy.x.y0,
|
|
accuracy.x.y1=accuracy.x.y1,
|
|
accuracy.z.y0=accuracy.z.y0,
|
|
accuracy.z.y1=accuracy.z.y1,
|
|
cor.y.xi=cor.y.xi,
|
|
fnr=fnr,
|
|
fnr.y0=fnr.y0,
|
|
fnr.y1=fnr.y1,
|
|
fpr=fpr,
|
|
fpr.y0=fpr.y0,
|
|
fpr.y1=fpr.y1,
|
|
cor.resid.w_pred=cor.resid.w_pred
|
|
))
|
|
|
|
result <- append(result, list(cor.xz=cor(df$x,df$z)))
|
|
(model.true <- lm(y ~ x + z, data=df))
|
|
true.ci.Bxy <- confint(model.true)['x',]
|
|
true.ci.Bzy <- confint(model.true)['z',]
|
|
|
|
result <- append(result, list(Bxy.est.true=coef(model.true)['x'],
|
|
Bzy.est.true=coef(model.true)['z'],
|
|
Bxy.ci.upper.true = true.ci.Bxy[2],
|
|
Bxy.ci.lower.true = true.ci.Bxy[1],
|
|
Bzy.ci.upper.true = true.ci.Bzy[2],
|
|
Bzy.ci.lower.true = true.ci.Bzy[1]))
|
|
|
|
## mle_result <- list(
|
|
## Bxy.est.naive = NULL,
|
|
## Bxy.ci.upper.naive = NULL,
|
|
## Bxy.ci.lower.naive = NULL,
|
|
## Bzy.est.naive = NULL,
|
|
## Bzy.ci.upper.naive = NULL,
|
|
## Bzy.ci.lower.naive = NULL,
|
|
## Bxy.est.feasible = NULL,
|
|
## Bxy.ci.upper.feasible = NULL,
|
|
## Bxy.ci.lower.feasible = NULL,
|
|
## Bzy.est.feasible = NULL,
|
|
## Bzy.ci.upper.feasible = NULL,
|
|
## Bzy.ci.lower.feasible = NULL,
|
|
## Bxy.est.mle = NULL,
|
|
## Bxy.ci.upper.mle = NULL,
|
|
## Bxy.ci.lower.mle = NULL,
|
|
## Bzy.est.mle = NULL,
|
|
## Bzy.ci.upper.mle = NULL,
|
|
## Bzy.ci.lower.mle = NULL
|
|
## )
|
|
## tryCatch({
|
|
mod <- glm_fixit(..., family=family, data=data, data2=data2)
|
|
|
|
coef_corrected <- coef(mod, which_model = "corrected")
|
|
ci_corrected <- confint(mod, which_model = "corrected")
|
|
|
|
coef_feasible <- coef(mod, which_model = "feasible")
|
|
ci_feasible <- confint(mod, which_model = "feasible")
|
|
|
|
coef_naive <- coef(mod, which_model = "naive")
|
|
ci_naive <- confint(mod, which_model = "naive")
|
|
|
|
mle_result <- list(
|
|
Bxy.est.mle = coef_corrected["x"],
|
|
Bxy.ci.upper.mle = ci_corrected["x", 2],
|
|
Bxy.ci.lower.mle = ci_corrected["x", 1],
|
|
Bzy.est.mle = coef_corrected["z"],
|
|
Bzy.ci.upper.mle = ci_corrected["z", 2],
|
|
Bzy.ci.lower.mle = ci_corrected["z", 1],
|
|
Bxy.est.feasible = coef_feasible["x"],
|
|
Bxy.ci.upper.feasible = ci_feasible["x", 2],
|
|
Bxy.ci.lower.feasible = ci_feasible["x",1],
|
|
Bzy.est.feasible = coef_feasible["z"],
|
|
Bzy.ci.upper.feasible = ci_feasible["z",2],
|
|
Bzy.ci.lower.feasible = ci_feasible["z",1],
|
|
Bxy.est.naive = coef_naive["w"],
|
|
Bxy.ci.upper.naive = ci_naive["w",2],
|
|
Bxy.ci.lower.naive = ci_naive["w",1],
|
|
Bzy.est.naive = coef_naive["a"],
|
|
Bzy.ci.upper.naive = ci_naive["a",2],
|
|
Bzy.ci.lower.naive = ci_naive["a",1]
|
|
)
|
|
## print(mle_result)
|
|
## },
|
|
## error=function(e) {result[['error']] <- as.character(e)
|
|
## })
|
|
result <- append(result, mle_result)
|
|
return(result)
|
|
}
|