Skip to content

Commit

Permalink
Add Probit model
Browse files Browse the repository at this point in the history
  • Loading branch information
xsswang committed Apr 19, 2022
1 parent fc9d04f commit 17929e4
Show file tree
Hide file tree
Showing 19 changed files with 1,603 additions and 72 deletions.
2 changes: 2 additions & 0 deletions NAMESPACE
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@ S3method(coef,summary.remiod)
S3method(print,summary.remiod)
S3method(summary,remiod)
export(extract_MIdata)
export(get_subset)
export(list.models)
export(mcmcplot)
export(miAnalyze)
Expand Down Expand Up @@ -31,5 +32,6 @@ importFrom(future,sequential)
importFrom(mcmcse,mcse)
importFrom(rjags,coda.samples)
importFrom(rjags,jags.model)
importFrom(utils,getAnywhere)
importFrom(utils,getFromNamespace)
importFrom(utils,globalVariables)
134 changes: 134 additions & 0 deletions R/JAGSmodel_opm.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,134 @@

paste_phi = function (nr, env = parent.frame()) {
vn <- env$.info$varname
ind <- env$.index
gk <- env$.isgk
paste0("phi( c_",vn,"[", nr, "] - eta_",vn,"[",ind, "] ) ")
}

write_opm = function (info, index, isgk = FALSE, indent = 4L) {
.info <- info
.index <- index
.isgk <- isgk
probs <- cvapply(2L:(info$ncat - 1L), function(k, .info = info, .index = index, .isgk = isgk) {
paste0(tab(indent), paste_p(k), " <- ", (if (isTRUE(.info$rev)) {
paste0(paste_phi(k - 1L), " - ", paste_phi(k))
}
else {
paste0(paste_phi(k), " - ", paste_phi(k - 1L))
}))
})

paste0(tab(indent), paste_p(1L), " <- ", if (isTRUE(info$rev)) {
1 - paste_phi(1L)
} else {
paste_phi(1L)
}, "\n",
paste(probs, collapse = "\n"), "\n",

tab(indent), paste_p(info$ncat), " <- ", if (isTRUE(info$rev)) {
paste_phi(info$ncat)
} else {
paste0("1 - ", paste_phi(info$ncat - 1L))
})
}

write_priors_opm = function (info) {

# sign <- if (isTRUE(info$rev)) {
# " + "
# } else {
# " - "
# }
gammas <- cvapply(1L:(info$ncat - 1L), function(k) {
if (k == 1L) {
paste0(tab(), "gamma_", info$varname, "[",
k, "] ~ dnorm(mu_delta_ordinal, tau_delta_ordinal)")
}
else {
paste0(tab(), "gamma_", info$varname, "[",
k, "] ~ dgamma (0.001, 0.001)")
}
})

deltas <- cvapply(1L:(info$ncat - 1L), function(k) {
paste0(tab(), "c_", info$varname, "[",
k, "] <- sum( gamma_", info$varname, "[1:",k,"])")
})

paste0(c(gammas, "", deltas), collapse = "\n")
}




jagsmodel_opm = function (info) {
if (info$ncat < 3L) {
errormsg("A cumulative logit mixed model is supposed to be fitted for the\n
variable %s but %s only has %s categories.",
dQuote(info$varname), dQuote(info$varname), info$ncat)
}
if (!is.null(info$hc_list)) {
errormsg("I found a random effects structure. Did you mean to use %s\n
instead of %s?",
dQuote("clmm"), dQuote("clm"))
}
indent <- 4L + 4L + nchar(info$varname) + 7L
index <- info$index[gsub("M_", "", info$resp_mat)]
linpred <- if (length(info$lp[[info$resp_mat]]) > 0L) {
paste_linpred(info$parname, info$parelmts[[info$resp_mat]],
matnam = info$resp_mat, index = index, cols = info$lp[[info$resp_mat]],
scale_pars = info$scale_pars[[info$resp_mat]])
}
else {
"0"
}
linpred_nonprop <- if (!is.null(attr(info$parelmts[[info$resp_mat]],
"nonprop"))) {
rhs <- cvapply(attr(info$parelmts[[info$resp_mat]], "nonprop"),
function(par_elmts) {
add_linebreaks(paste_linpred(info$parname, par_elmts,
matnam = info$resp_mat, index = index, cols = attr(info$lp,
"nonprop")[[info$resp_mat]], scale_pars = info$scale_pars[[info$resp_mat]]),
indent = indent)
})
paste0("\n\n", paste0(tab(4L), "eta_", info$varname,
"_", seq_along(rhs), "[", index, "] <- ",
rhs, collapse = "\n"))
}
dummies <- if (!is.null(info$dummy_cols)) {
paste0("\n", paste0(paste_dummies(resp_mat = info$resp_mat,
resp_col = info$resp_col, dummy_cols = info$dummy_cols,
index = index, refs = info$refs), collapse = "\n"),
"\n")
}
paste0("\r", tab(), add_dashes(paste0("# Cumulative logit model for ",
info$varname)), "\n", tab(), "for (", index,
" in 1:", info$N[gsub("M_", "", info$resp_mat)],
") {", "\n", tab(4L), info$resp_mat, "[",
index, ", ", info$resp_col, "] ~ dcat(p_",
info$varname, "[", index, ", 1:", info$ncat,
"])", "\n", tab(4L), "eta_", info$varname,
"[", index, "] <- ", add_linebreaks(linpred,
indent = indent), linpred_nonprop, "\n\n",
write_opm(info, index), "\n\n",
#write_logits(info,index, nonprop = !is.null(linpred_nonprop)), "\n",
dummies,
info$trafos, tab(), "}", "\n\n",
tab(), "# Priors for the model for ", info$varname,
"\n", if (!is.null(info$lp[[info$resp_mat]])) {
paste0(tab(), "for (k in ", min(unlist(c(info$parelmts,
lapply(info$parelmts, attr, "nonprop")))),
":", max(unlist(c(info$parelmts, lapply(info$parelmts,
attr, "nonprop")))), ") {", "\n",
get_priordistr(info$shrinkage, type = "ordinal",
parname = info$parname), tab(), "}",
"\n\n")
}, write_priors_opm(info))
}






17 changes: 12 additions & 5 deletions R/get_MI_RB.R
Original file line number Diff line number Diff line change
Expand Up @@ -61,18 +61,25 @@ get_MI_RB <- function(object, treatment, method=c("MAR","J2R","CR","delta"), del

# get a summary of the relevant characteristics of the imputed variables
varinfo <- lapply(object$info_list[vars], function(x) {
#if (x$modeltype == 'opm') x$modeltype == 'clm'
data.frame(varname = x$varname,
modeltype = x$modeltype,
family = ifelse(!is.null(x$family), x$family, NA),
stringsAsFactors = FALSE)
})

if (varinfo[[1]]$modeltype == 'clm'){
mcUpdateFun = switch(method,
'MAR' = prep_MCMC,
'J2R' = clm_MI_J2R,
'CR' = clm_MI_CR,
'delta' = clm_MI_delta)
mcUpdateFun = switch(method,
'MAR' = prep_MCMC,
'J2R' = clm_MI_J2R,
'CR' = clm_MI_CR,
'delta' = clm_MI_delta)
} else if (varinfo[[1]]$modeltype == 'opm') {
mcUpdateFun = switch(method,
'MAR' = prep_MCMC,
'J2R' = opm_MI_J2R,
'CR' = opm_MI_CR,
'delta' = opm_MI_delta)
} else {
mcUpdateFun = switch(method,
'MAR' = prep_MCMC,
Expand Down
Loading

0 comments on commit 17929e4

Please sign in to comment.