R/runPCModel_final.R
In cammiller/imagingPC: Analysis of Imaging Mass Spectrometry Data using a Process Convolution Approach
## code to create data for model and run model
#' Run the model
#'
#' The \code{runPCModel} function runs the model created by the
#' \code{\link{writePCModel}} function, using the information generated by
#' the \code{\link{createPCData}} function.
#'
#' @param modelObj An object of class \code{PCModelList}. This object is the
#' result of using the \code{\link{writePCModel}} function.
#' @param PCDataObj An object of class \code{PCDataList}. This object is the
#' result of using the \code{\link{createPCData}} function.
#' @param nIter An integer specifying he total number of iterations for the
#' MCMC, pre-thinning. This includes the burn-in.
#' @param nBurnin An integer specifying the number of iterations
#' (pre-thinning) to be discarded.
#' @param nThin An integer specifying the thinning interval.
#' @param slideVar An optional argument specifying a slide variable. This
#' option is only useful for data with zeros and when those zeros are assumed
#' to be censored observations. In this case, the censoring threshold can
#' vary by each run of the mass spectrometer. When this argument is used,
#' the censoring threshold will be allowed to vary by slide.
#' @param keepModel A TRUE/FALSE argument specifying whether or not to keep
#' the model (found in the modelObj).
#' @param keepModelObj A TRUE/FALSE argument specifying whether the
#' compiled NIMBLE model object should be kept after running the model. This
#' object can be large, but it is required for restarting the MCMC.
#' @param monitorCoefOnly A TRUE/FALSE argument stating whether to monitor the
#' the model coefficients only (TRUE) or all stochastic nodes (FALSE).
#' Regardless of the value of this argument, if any categorical covariates are
#' included in the model, the differences between all subgroups of the
#' categorical covariate(s) will be monitored in the MCMC with assigned nodes.
#' @param returnSample A TRUE/FALSE argument specifying whether to return the
#' entire MCMC sample. If this argument is FALSE, only the summary measures
#' will be returned.
#' @param returnAllSummary A TRUE/FALSE argument specifying whether to return
#' summary measures for all monitored nodes.
#' @param justReturnData A TRUE/FALSE argument. If justReturnData=TRUE, then
#' no MCMC will be run. Instead, the data and constants lists required to run
#' the model will be returned. This is useful for modifying the model to
#' personal specifications.
#' @param convergenceTol A numeric value indicating the highest acceptable
#' value for the Brooks-Gelman-Rubin (BGR) statistic to consider model
#' nodes converged.
#' @param maxRuns A numeric value indicating the maximum number of times a
#' model will be restarted to try to achieve convergence. After running the
#' MCMC for nIter iterations, the model is checked for convergence. A model
#' is considered converged if the model coefficients converge according to
#' the BGR statistic. For the marginalized two-part model, only the
#' coefficients for the marginalized part are used. If the model does not
#' converge, then the MCMC will be restarted and run for another
#' iterIncrement iterations (post-thinning). Convergence will then be
#' rechecked. The maxRuns argument specifies how many times the MCMC can
#' be restarted after the initial run of nIter iterations.
#' @param iterWindow A numeric value specifying the post-thinning number of
#' iterations (per chain) that will be used to make inference, if the model
#' does not converge with the first nIter iterations.
#' @param sliceSamplers A TRUE/FALSE argument specifying whether or not to
#' use slice samplers for all stochastic nodes. If
#' \code{sliceSamplers=TRUE} then NIMBLE's onlySlice option in the
#' \code{\link[nimble]{configureMCMC}} is set to TRUE. If
#' \code{sliceSamplers=FALSE} then NIMBLE's default MCMC settings are used.
#' Using slice samplers is useful for reducing autocorrelation in the MCMC,
#' though it significantly increases computation time.
#'
#' @return A list of class \code{PCResults} containing the results of running
#' the MCMC. The list will contain additional objects and information
#' specified by the arguments.
#'
#' \describe{
#' \item{\code{results}}{A matrix containing summary statistics (mean
#' and the 2.5%, 25%, 50%, 75%, and 97.5% percentiles) for each model
#' coefficient. For ease of reading, the names of the model coefficients
#' (beta1, beta2,...) are replaced with the names of their corresponding
#' covariates.}
#' \item{\code{model}}{The model specified from the
#' \code{\link{writePCModel}} function. The model is output if
#' \code{keepModel=TRUE}.}
#' \item{\code{modelObj}}{The compiled NIMBLE model object after running
#' until covergence or maxRuns is reached. This object is output only if
#' \code{keepModelObj=TRUE}. This object can be large (>100MB) for even
#' small imaging mass spectrometry datasets, so some thought should be
#' put into deciding whether or not to keep this object.}
#' \item{\code{fullSummary}}{A matrix of summary information for all
#' nodes in the model, This will be returned if
#' \code{returnAllSummary=TRUE}.}
#' \item{\code{sample}}{An \code{mcmc.list} object. See the \code{coda}
#' package for more details on such objects. This is a list of two
#' matrices, one for each chain. Each matrix has rows equal to the
#' number of thinned iterations and columns equal to the number of nodes
#' monitored. The sample is returned if \code{returnSample=TRUE}.}
#' \item{\code{dataList}}{The data list used to create a compiled
#' NIMBLE model. For all data this list will include the matrix
#' generated from the smoothing kernel function. For data with zeros
#' (censored and true), the model is specified using the zeros trick. In
#' this case \code{dataList} also has a vector of zeros.}
#' \item{\code{constantsList}}{The list of constants used to create
#' a compiled NIMBLE model. This list is typically composed of
#' covariates and other vectors needed to navigate the covariate
#' information. However, for data with zeros (censored and true),
#' the model is specified using the zeros trick, and so the outcome
#' data is also included in this list.}
#' \item{\code{convergenceTol}}{The limit on the BGR statistic used to
#' determine convergence. The default is 1.10, so if the BGR values for
#' all model coefficients are <=1.10, then the model is considered
#' converged.}
#' }
#'
#' @examples
#' data("TAMdata")
#' # The dataset is trimmed only for the speed of the example
#' TAMdata <- TAMdata[TAMdata$subject < 3, ]
#' TAMdata <- rScale(TAMdata, subjectVar = 'subject', sampleVar = 'ROI',
#' xCoord = 'x', yCoord = 'y')
#' rangs <- estRange(TAMdata, outcome = 'X1282.auc', spatialVar = 'TAM',
#' semivEst = 'modulus', logTransform = TRUE)
#' structs <- chooseStructures(rangs)
#' PCdat <- createPCData(structs, trimData = FALSE,
#' covariates = c("secondary", "TAM", "secTAM"),
#' covariateTypes = c("binary", "binary", "binary"),
#' covariateLevels = c("sample", "raster", "raster"))
#' PCmod <- writePCModel(PCdat, multiSampsPerSubj = TRUE, typeOfZero = "censored")
#' PCresults <- runPCModel(modelObj = PCmod, PCDataObj = PCdat, slideVar='slide',
#' monitorCoefOnly = FALSE,
#' nBurnin = 15000, nIter = 40000, nThin = 25)
#'
#' @references de Valpine, P., D. Turek, C.J. Paciorek, C. Anderson-Bergman,
#' D. Temple Lang, and R. Bodik. 2017. Programming with models: writing
#' statistical algorithms for general model structures with NIMBLE.
#' \emph{Journal of Computational and Graphical Statistics}. 26: 403-413.
runPCModel <- function(modelObj, PCDataObj, nIter = 50000, nBurnin = 10000, nThin = 5, slideVar = NULL, keepModel = TRUE, keepModelObj = FALSE,
monitorCoefOnly = TRUE, returnSample = TRUE, returnAllSummary = FALSE, justReturnData = FALSE, convergenceTol = 1.1, maxRuns = 5, iterIncrement = 500,
iterWindow = 1000, sliceSamplers=FALSE) {
if (class(modelObj) != "PCModelList")
stop("modelObj must of class PCModelList This object must be created using the writePCModel function.")
if (class(PCDataObj) != "PCDataList")
stop("PCDataObj must of class PCDataList. This object must be created using the createPCData function.")
if (is.numeric(nIter) == FALSE)
stop("nIter must be numeric.")
if ((nIter%%1 == 0) == FALSE)
stop("nIter must be an integer.")
if (nIter <= 0)
stop("nIter must be greater than 0.")
if (is.numeric(nBurnin) == FALSE)
stop("nBurnin must be numeric.")
if ((nBurnin%%1 == 0) == FALSE)
stop("nBurnin must be an integer.")
if (nBurnin < 0)
stop("nBurnin must be greater than or equal to 0.")
if (is.numeric(nThin) == FALSE)
stop("nThin must be numeric.")
if ((nThin%%1 == 0) == FALSE)
stop("nThin must be an integer.")
if (nThin <= 0)
stop("nThin must be greater than 0.")
if (is.numeric(convergenceTol) == FALSE)
stop("convergenceTol must be numeric.")
if (convergenceTol <= 1)
stop("convergenceTol should be greater than 1.")
if (is.numeric(maxRuns) == FALSE)
stop("maxRuns must be numeric.")
if ((maxRuns%%1 == 0) == FALSE)
stop("maxRuns must be an integer.")
if (maxRuns < 0)
stop("maxRuns cannot be negative.")
if (is.numeric(iterIncrement) == FALSE)
stop("iterIncrement must be numeric.")
if ((iterIncrement%%1 == 0) == FALSE)
stop("iterIncrement must be an integer.")
if (iterIncrement <= 0)
stop("iterIncrement must be greater than 0.")
if (is.numeric(iterWindow) == FALSE)
stop("iterWindow must be numeric.")
if ((iterWindow%%1 == 0) == FALSE)
stop("iterWindow must be an integer.")
if (iterWindow <= 0)
stop("iterWindow must be greater than 0.")
if (iterWindow > (((nIter - nBurnin)/nThin) + iterIncrement))
stop("iterWindow must be smaller than ((nIter-nBurnin)/nThin)+iterIncrement")
covs <- PCDataObj[["covs"]]
subjectVar <- PCDataObj[["subjectVar"]]
spatialVar <- PCDataObj[["spatialVar"]]
outcome <- PCDataObj[["outcome"]]
recStructures<-PCDataObj[['recStructures']]
model <- modelObj[["model"]]
covariates <- modelObj[["covariates"]]
covariateLevels<-modelObj[['covriateLevels']]
covariateTypes <- modelObj[["covariateTypes"]]
covariatesForBinary<-modelObj[['covariatesForBinary']]
rastpervar<-PCDataObj[['rastersPerVar']]
if(is.null(spatialVar)==FALSE){
#PCDataObj[['data']] <- PCDataObj[['data']][order(PCDataObj[['data']][[spatialVar]]), ]
spatvarlevels <- sort(unique(PCDataObj[['data']][[spatialVar]]))
nvarlevels <- length(spatvarlevels)
}
covstokeep <- rep(NA, length(covs))
for (i in 1:length(covs)) {
covstokeep[i] <- ifelse(covs[[i]]$covariate %in% covariates, 1, 0)
}
covstokeep <- which(covstokeep == 1)
covsb <- covs[covstokeep]
ordcovnames <- rep(NA, length(covs))
for (i in 1:length(covs)) {
ordcovnames[i] <- covs[[i]]$covariate
}
covariates <- factor(covariates, levels = ordcovnames)
covariates <- sort(covariates)
covariates <- as.character(covariates)
covariatesb <- covariates[covstokeep]
covariateTypesb <- covariateTypes[covstokeep]
covariateLevelsb <- covariateLevels[covstokeep]
nvars <- length(covsb)
covarnamesb <- list()
for (i in 1:nvars) {
if (covariateTypesb[i] %in% c("binary", "continuous")) {
covarnamesb[[i]] <- paste0(covariates[i])
}
if (covariateTypesb[i] == "categorical") {
covarnamesb[[i]] <- paste0(covariates[i], "_", colnames(covsb[[i]]$mapping))
}
}
var.namesb <- paste(unlist(covarnamesb))
if (nrow(PCDataObj[["data"]][PCDataObj[["data"]][[outcome]] == 0, ]) > 0 & modelObj[["typeOfZero"]] == "true") {
if (is.null(covariatesForBinary) == TRUE) {
covsa <- covsb
covariatesForBinary <- covariatesb
}
if (is.null(covariatesForBinary) == FALSE) {
covstokeep <- rep(NA, length(covs))
for (i in 1:length(covs)) {
covstokeep[i] <- ifelse(covs[[i]]$covariate %in% covariatesForBinary, 1, 0)
}
covstokeep <- which(covstokeep == 1)
covsa <- covs[covstokeep]
}
# need to get correct order of covariatesa
covariatesForBinary <- factor(covariatesForBinary, levels = ordcovnames)
covariatesForBinary <- sort(covariatesForBinary)
covariatesForBinary <- as.character(covariatesForBinary)
covariatesa <- covariatesForBinary[covstokeep]
covariateTypesa <- covariateTypes[covstokeep]
covariateLevelsa <- covariateLevels[covstokeep]
nvarsa <- length(covsa)
covarnamesa <- list()
for (i in 1:nvarsa) {
if (covariateTypesa[i] %in% c("binary", "continuous")) {
covarnamesa[[i]] <- paste0(covariatesForBinary[i])
}
if (covariateTypesa[i] == "categorical") {
covarnamesa[[i]] <- paste0(covariatesForBinary[i], "_", colnames(covsa[[i]]$mapping))
}
}
var.namesa <- paste(unlist(covarnamesa))
var.namesall <- unique(c(var.namesa, var.namesb))
}
if(is.null(spatialVar)==FALSE & any(recStructures==0)){
for(i in 1:length(unique(rastpervar$cov.level))){
assign(paste0('rastpervar',i), rastpervar[rastpervar$cov.level==spatvarlevels[i],])
}
}
############################################## Create data and constant lists #####
# create constants common to all models
const <- list()
for (i in 1:length(covsb)) {
if (covsb[[i]]$type == "categorical") {
for (j in 1:ncol(covsb[[i]]$mapping)) {
const[[paste0(covsb[[i]]$covariate, "_", colnames(covsb[[i]]$mapping)[j])]] <- covsb[[i]]$mapping[, j]
}
}
if (covsb[[i]]$type != "categorical") {
if (covsb[[i]]$level != "raster") {
const[[paste0(covsb[[i]]$covariate)]] <- covsb[[i]]$info[, 2]
}
if (covsb[[i]]$level == "raster") {
const[[paste0(covsb[[i]]$covariate)]] <- covsb[[i]]$info
}
}
}
if (nrow(PCDataObj[["data"]][PCDataObj[["data"]][[outcome]] == 0, ]) > 0 & modelObj[["typeOfZero"]] == "true") {
for (i in 1:length(covsa)) {
if (covsa[[i]]$type == "categorical") {
for (j in 1:ncol(covsa[[i]]$mapping)) {
const[[paste0(covsa[[i]]$covariate, "_", colnames(covsa[[i]]$mapping)[j])]] <- covsa[[i]]$mapping[, j]
}
}
if (covsa[[i]]$type != "categorical") {
if (covsa[[i]]$level != "raster") {
const[[paste0(covsa[[i]]$covariate)]] <- covsa[[i]]$info[, 2]
}
if (covsa[[i]]$level == "raster") {
const[[paste0(covsa[[i]]$covariate)]] <- covsa[[i]]$info
}
}
}
}
if (is.null(subjectVar) == FALSE) {
const[["Subj"]] <- PCDataObj[["nSubjs"]]
const[["Sample"]] <- PCDataObj[["cNSampsPerSubj"]]
const[["Raster"]] <- PCDataObj[["cNRastPerSamp"]]
}
if (is.null(subjectVar) == TRUE) {
const[["Sample"]] <- PCDataObj[["nSamps"]]
const[["Raster"]] <- PCDataObj[["cNRastPerSamp"]]
}
if (is.null(spatialVar) == FALSE) {
if(sum(PCDataObj[['recStructures']])==length(PCDataObj[['recStructures']])){
for (i in 1:PCDataObj$nVarLevels) {
const[[paste0("nsup", i)]] <- PCDataObj[["nSupportSites"]][, i]
}
}
if(any(PCDataObj[['recStructures']]==0)){
for (i in 1:PCDataObj$nVarLevels){
if(PCDataObj[['recStructures']][i]==1){
const[[paste0("nsup", i)]] <- PCDataObj[["nSupportSites"]][, i]
}
if(PCDataObj[['recStructures']][i]==0){
#create vector of start and stop
startstop<-rep(NA,(nrow(eval(parse(text = paste0('rastpervar',i))))*2))
for(j in 1:nrow(eval(parse(text = paste0('rastpervar',i))))){
startstop[(j-1)*2+1]<-eval(parse(text = paste0('rastpervar',i)))$start[j]
startstop[(j-1)*2+2]<-eval(parse(text = paste0('rastpervar',i)))$stop[j]
}
const[[paste0('samprast',i)]]<-startstop
}
}
}
}
if (is.null(spatialVar) == TRUE & all(PCDataObj[['recStructures']]==1)) {
const[["nsup"]] <- PCDataObj[["nSupportSites"]]
}
datl <- list()
# lognormal likelihood
if (nrow(PCDataObj[["data"]][PCDataObj[["data"]][[outcome]] == 0, ]) == 0) {
datl[["y"]] <- PCDataObj[["data"]][[outcome]]
}
# censored and true models
if (nrow(PCDataObj[["data"]][PCDataObj[["data"]][[outcome]] == 0, ]) > 0) {
const[["C"]] <- 10000
PCDataObj[["data"]]$zeros <- 0
datl[["zeros"]] <- PCDataObj[["data"]]$zeros
datl[["Kmat"]] <- PCDataObj[["KMat"]]
if (modelObj[["typeOfZero"]] == "censored")
{
PCDataObj[["data"]]$ind <- NA
for (i in 1:nrow(PCDataObj[["data"]])) {
PCDataObj[["data"]]$ind[i] <- ifelse(PCDataObj[["data"]][[outcome]][i] > 0, 1, 0)
}
const[["obs"]] <- PCDataObj[["data"]]$ind
if (is.null(slideVar) == TRUE) {
tstar <- min(PCDataObj[["data"]][PCDataObj[["data"]][[outcome]] > 0, outcome])
for (i in 1:nrow(PCDataObj[["data"]])) {
if (PCDataObj[["data"]][[outcome]][i] == 0) {
PCDataObj[["data"]][[outcome]][i] <- tstar
}
}
}
if (is.null(slideVar) == FALSE) {
tstarmat <- matrix(NA, nrow = length(unique(PCDataObj[["data"]][[slideVar]])), ncol = 2)
tstarmat <- as.data.frame(tstarmat)
colnames(tstarmat) <- c("slide", "tstar")
tstarmat$slide <- unique(PCDataObj[["data"]][[slideVar]])
for (i in 1:nrow(tstarmat)) {
tstarmat$tstar[i] <- min(PCDataObj[["data"]][[outcome]][PCDataObj[["data"]][[outcome]] > 0 & PCDataObj[["data"]][[slideVar]] == tstarmat$slide[i]])
}
for (i in 1:nrow(PCDataObj[["data"]])) {
if (PCDataObj[["data"]][[outcome]][i] == 0) {
PCDataObj[["data"]][[outcome]][i] <- tstarmat$tstar[tstarmat$slide == PCDataObj[["data"]][[slideVar]][i]]
}
}
}
const[["y"]] <- PCDataObj[["data"]][[outcome]]
# llcal.censored <- nimbleFunction(run = function(y = double(0), obs = double(0), sigma = double(0), mu = double(0)) {
# returnType(double(0))
# if (obs != 0) {
# return(log(dnorm((log(y) - mu), sd = sigma)))
# }
# if (obs == 0) {
# return(log(pnorm((log(y) - mu), sd = sigma)))
# }
# })
assign('llcal.censored',
eval(parse(text = c('nimbleFunction(run = function(y = double(0), obs = double(0), sigma = double(0), mu = double(0)) {',
'returnType(double(0))',
'if (obs != 0) {',
'return(log(dnorm((log(y) - mu), sd = sigma)))',
'}',
'if (obs == 0) {',
'return(log(pnorm((log(y) - mu), sd = sigma)))',
'}',
'})'))), envir=.GlobalEnv)
} # censored model
if (modelObj[["typeOfZero"]] == "true")
{
const[["y"]] <- PCDataObj[["data"]][[outcome]]
llcal.true <- nimbleFunction(run = function(y = double(0), binprob = double(0), sigma = double(0), mu = double(0)) {
returnType(double(0))
if (y != 0) {
return((log(binprob) - log(y) - 0.5 * log(2 * 3.14159265359) - log(sigma) - (1/(2 * pow(sigma, 2))) * pow((log(y) - mu), 2)))
}
if (y == 0) {
return(log(1 - binprob))
}
})
} # true zeros model
}
################################# Create inits list #####
inits <- list()
# inits common to all models
inits[["beta0"]] <- mean(log(PCDataObj[["data"]][[outcome]][PCDataObj[["data"]][[outcome]] > 0]))
betanames <- list()
for (i in 1:nvars) {
if (covariateTypesb[i] %in% c("binary", "continuous")) {
betanames[[i]] <- paste0("beta", i)
}
if (covariateTypesb[i] == "categorical") {
betanames[[i]] <- paste0("beta", i, "_", seq(1:ncol(covsb[[i]]$mapping)))
}
}
betanames <- paste(unlist(betanames))
for (i in 1:length(betanames)) {
inits[[paste0(betanames[i])]] <- 0
}
if (modelObj[["coefPrior"]] == "sdunif") {
inits[["sdb0"]] <- 0.1
sdbnames <- list()
for (i in 1:nvars) {
if (covariateTypesb[i] %in% c("binary", "continuous")) {
sdbnames[[i]] <- paste0("sdb", i)
}
if (covariateTypesb[i] == "categorical") {
sdbnames[[i]] <- paste0("sdb", i, "_", seq(1:ncol(covsb[[i]]$mapping)))
}
}
sdbnames <- paste(unlist(sdbnames))
for (i in 1:length(sdbnames)) {
inits[[paste0(sdbnames[i])]] <- 0.1
}
}
if (modelObj[["errorVarianceLevel"]] == "overall") {
inits[["sigma"]] <- 0.1
}
if (modelObj[["errorVarianceLevel"]] == "subject") {
inits[["sigma"]] <- rep(0.1, PCDataObj[["nSubjs"]])
}
if (modelObj[["errorVarianceLevel"]] == "sample") {
inits[["sigma"]] <- rep(0.1, PCDataObj[["nSamps"]])
}
if (modelObj[["multiSampsPerSubj"]] == TRUE & nrow(PCDataObj[["data"]][PCDataObj[["data"]][[outcome]] == 0, ]) > 0 & modelObj[["typeOfZero"]] ==
"censored") {
inits[["bi"]] <- rep(0, PCDataObj[["nSubjs"]])
inits[["sdbi"]] <- 0.1
}
if (modelObj[["multiSampsPerSubj"]] == TRUE & nrow(PCDataObj[["data"]][PCDataObj[["data"]][[outcome]] == 0, ]) == 0) {
inits[["bi"]] <- rep(0, PCDataObj[["nSubjs"]])
inits[["sdbi"]] <- 0.1
}
if (modelObj[["multiSampsPerSubj"]] == TRUE & nrow(PCDataObj[["data"]][PCDataObj[["data"]][[outcome]] == 0, ]) > 0 & modelObj[["typeOfZero"]] ==
"true") {
inits[["bi1"]] <- rep(0, PCDataObj[["nSubjs"]])
inits[["bi2"]] <- rep(0, PCDataObj[["nSubjs"]])
inits[["rho"]] <- 0
}
if (nrow(PCDataObj[["data"]][PCDataObj[["data"]][[outcome]] == 0, ]) > 0 & modelObj[["typeOfZero"]] == "true") {
alphanames <- list()
for (i in 1:nvarsa) {
if (covariateTypesa[i] %in% c("binary", "continuous")) {
alphanames[[i]] <- paste0("alpha", i)
}
if (covariateTypesa[i] == "categorical") {
alphanames[[i]] <- paste0("alpha", i, "_", seq(1:ncol(covsa[[i]]$mapping)))
}
}
alphanames <- paste(unlist(alphanames))
for (i in 1:length(alphanames)) {
inits[[paste0(alphanames[i])]] <- 0
}
if (modelObj[["coefPrior"]] == "sdunif") {
inits[["sda0"]] <- 0.1
sdanames <- list()
for (i in 1:nvars) {
if (covariateTypesa[i] %in% c("binary", "continuous")) {
sdanames[[i]] <- paste0("sda", i)
}
if (covariateTypesa[i] == "categorical") {
sdanames[[i]] <- paste0("sda", i, "_", seq(1:ncol(covsa[[i]]$mapping)))
}
}
sdanames <- paste(unlist(sdanames))
for (i in 1:length(sdanames)) {
inits[[paste0(sdanames[i])]] <- 0.1
}
}
inits[["v"]] <- rep(0, PCDataObj[["nSamps"]])
inits[["sdv"]] <- 0.1
}
if (is.null(spatialVar) == TRUE & all(PCDataObj[['recStructures']]==1)) {
xc <- matrix(nrow = PCDataObj[["nSamps"]], ncol = max(PCDataObj[["nSupportSites"]]))
for (i in 1:PCDataObj[["nSamps"]]) {
xc[i, 1:PCDataObj[["nSupportSites"]][i]] <- 0
}
inits[["xc"]] <- xc
if (modelObj[["latentVarianceLevel"]] == "overall") {
inits[["sdxc"]] <- 0.1
}
if (modelObj[["latentVarianceLevel"]] == "subject") {
inits[["sdxc"]] <- rep(0.1, PCDataObj[["nSubjs"]])
}
if (modelObj[["latentVarianceLevel"]] == "sample") {
inits[["sdxc"]] <- rep(0.1, PCDataObj[["nSamps"]])
}
}
if (is.null(spatialVar) == FALSE) {
for (j in 1:PCDataObj[["nVarLevels"]]) {
if(PCDataObj[['recStructures']][j]==1){
assign(paste0("xc", j), matrix(NA, nrow = length(PCDataObj[["nSupportSites"]][[paste0("nsups.", j)]]), ncol = max(PCDataObj[["nSupportSites"]][[paste0("nsups.",
j)]])))
eval(parse(text = c(paste0("for(i in 1:length(PCDataObj[[\"gT0SupportSites\"]][[", j, "]])){"), paste0("xc", j, "[PCDataObj[[\"gT0SupportSites\"]][[",
j, "]][i],1:PCDataObj[[\"nSupportSites\"]]$nsups.", j, "[PCDataObj[[\"gT0SupportSites\"]][[1]][i]]]<-0"), "}")))
eval(parse(text = paste0("inits[[\"xc", j, "\"]]<-xc", j)))
if (modelObj[["latentVarianceLevel"]] == "overall") {
assign(paste0("sdxc", j), 0.1)
eval(parse(text = paste0("inits[[\"sdxc", j, "\"]]<-sdxc", j)))
}
if (modelObj[["latentVarianceLevel"]] == "subject") {
assign(paste0("sdxc", j), rep(0.1, PCDataObj[["nSubjs"]]))
eval(parse(text = paste0("inits[[\"sdxc", j, "\"]]<-sdxc", j)))
}
if (modelObj[["latentVarianceLevel"]] == "sample") {
assign(paste0("sdxc", j), rep(0.1, PCDataObj[["nSamps"]]))
eval(parse(text = paste0("inits[[\"sdxc", j, "\"]]<-sdxc", j)))
}
}
if(PCDataObj[['recStructures']][j]==0){
eval(parse(text = paste0('inits[["rk',j,'"]]<-rep(NA,max(rastpervar$stop))')))
eval(parse(text = c(paste0('for(i in 1:nrow(rastpervar',j,')){'),paste0('inits[["rk',j,'"]][rastpervar',j,'$start[i]:rastpervar',j,'$stop[i]]<-0'),'}')))
eval(parse(text = paste0('inits[["sdrk',j,'"]]<-0.1')))
}
}
}
######################## run PC model #####
if (justReturnData == TRUE) {
outlist <- list()
if (keepModel == TRUE) {
outlist[["model"]] <- model
}
outlist[["dataList"]] <- datl
outlist[["constantsList"]] <- const
outlist[['initsList']]<-inits
outlist[["convergenceTol"]] <- convergenceTol
return(outlist)
}
if (justReturnData == FALSE) {
if ("categorical" %in% covariateTypesb) {
addbetanodes <- list()
nodetemp <- list()
for (i in 1:nvars) {
if (covariateTypesb[i] == "categorical") {
# to provide names, need to get mapping matrix, column of betas and column of corresponding name
betanames2 <- paste0("beta", i, "_", seq(1:ncol(covsb[[i]]$mapping)))
covarnames2 <- paste0(covariates[i], "_", colnames(covsb[[i]]$mapping))
namemap <- data.frame(betanames2, covarnames2)
nodetemp[[i]] <- paste0("beta", i, "_", seq(1:ncol(covsb[[i]]$mapping)))
nodemat <- expand.grid(paste(nodetemp[[i]]), paste(nodetemp[[i]]))
nodemat <- nodemat[nodemat[, 1] != nodemat[, 2], ]
nodemat[, 3] <- 0
for (j in 1:nrow(nodemat)) {
if (as.numeric(strsplit(as.character(nodemat[j, 1]), "_")[[1]][2]) > as.numeric(strsplit(as.character(nodemat[j, 2]), "_")[[1]][2])) {
nodemat[j, 3] <- 1
}
}
nodemat <- nodemat[nodemat[, 3] == 1, ]
nodemat <- nodemat[order(nodemat[, 1]), ]
nodemat <- nodemat[, c(1, 2)]
nodemat <- as.data.frame(nodemat)
colnames(nodemat) <- c("beta1", "beta2")
nodemat$cov1 <- NA
nodemat$cov2 <- NA
for (j in 1:nrow(nodemat)) {
nodemat$cov1[j] <- as.character(namemap$covarnames2[namemap$betanames2 == nodemat$beta1[j]])
nodemat$cov2[j] <- as.character(namemap$covarnames2[namemap$betanames2 == nodemat$beta2[j]])
}
# need to subtract columns
addbetanodes[[i]] <- rep(NA, nrow(nodemat))
for (j in 1:nrow(nodemat)) {
addbetanodes[[i]][j] <- paste0(nodemat$cov1[j], "_minus_", nodemat$cov2[j])
}
}
}
addbetanodes <- paste(unlist(addbetanodes))
}
# vector of betas
betas <- betanames
if (nrow(PCDataObj[["data"]][PCDataObj[["data"]][[outcome]] == 0, ]) > 0 & modelObj[["typeOfZero"]] == "true") {
betas<-c(betanames,alphanames)
}
nimmod <- nimbleModel(code = model, constants = const, data = datl, inits = inits)
cnimmod <- compileNimble(nimmod)
if (monitorCoefOnly == TRUE) {
confmod <- configureMCMC(nimmod, monitors = betas, print = FALSE, onlySlice = sliceSamplers)
}
if (monitorCoefOnly == FALSE) {
confmod <- configureMCMC(nimmod, monitors = cnimmod$getNodeNames(stochOnly = TRUE, includeData = FALSE), print = FALSE, onlySlice = sliceSamplers)
}
if (any(covariateTypesb == "categorical")) {
confmod$addMonitors(addbetanodes)
}
bmod <- buildMCMC(confmod)
modMCMC <- compileNimble(bmod, project = nimmod)
samp <- runMCMC(modMCMC, niter = nIter, nburnin = nBurnin, nchains = 2, samplesAsCodaMCMC = TRUE, thin = nThin)
samp[[1]] <- samp[[1]][, is.na(samp[[1]][1, ]) == FALSE]
samp[[2]] <- samp[[2]][, is.na(samp[[2]][1, ]) == FALSE]
samp <- as.mcmc.list(samp)
bgrs <- gelman.diag(samp, autoburnin = FALSE, multivariate = FALSE)
bgrtemp <- bgrs$psrf[betas, 1]
if (max(bgrtemp) <= convergenceTol | maxRuns==0) {
if(max(bgrtemp) <= convergenceTol){
cat('The model coefficients converged.\n')
} else {
cat('The model coefficients did not converge.\n')
}
}
if (max(bgrtemp) > convergenceTol) {
cat("With a burnin of", nBurnin, "and a sample of", nrow(samp[[1]]), "(post-thinning), the model coefficients did not converge. Will keep running to try to reach convergence.\n")
iter <- 1
flag <- 1
while (flag) {
sampnew <- runMCMC(modMCMC, niter = iterIncrement * nThin, nburnin = 0, nchains = 2, samplesAsCodaMCMC = TRUE, thin = nThin)
sampnew[[1]] <- sampnew[[1]][, is.na(sampnew[[1]][1, ]) == FALSE]
sampnew[[2]] <- sampnew[[2]][, is.na(sampnew[[2]][1, ]) == FALSE]
sampnew <- as.mcmc.list(sampnew)
samp[[1]] <- rbind(samp[[1]], sampnew[[1]])
samp[[2]] <- rbind(samp[[2]], sampnew[[2]])
samp[[1]] <- samp[[1]][(nrow(samp[[1]]) - (iterWindow - 1)):nrow(samp[[1]]), ]
samp[[2]] <- samp[[2]][(nrow(samp[[2]]) - (iterWindow - 1)):nrow(samp[[2]]), ]
samp[[1]] <- as.mcmc(samp[[1]])
samp[[2]] <- as.mcmc(samp[[2]])
samp <- as.mcmc.list(samp)
bgrs <- gelman.diag(samp, autoburnin = FALSE, multivariate = FALSE)
bgrtemp <- bgrs$psrf[betas, 1]
if (max(bgrtemp) <= convergenceTol | iter >= maxRuns) {
flag <- 0
}
if (max(bgrtemp) > convergenceTol) {
cat("Run", iter, "complete without convergence.")
}
if (iter == maxRuns) {
cat("Max number of runs reached.")
}
iter <- iter + 1
}
}
sampres <- summary(samp)
sampstat <- sampres$statistics
sampquant <- sampres$quantiles
ressum <- cbind(sampstat[, 1], sampquant)
colnames(ressum) <- c("Mean", "2.5%", "25%", "50%", "75%", "97.5%")
# row.names(ressum) which(row.names(ressum)==betas[1])
# reassign row names
for (i in 1:length(betanames)) {
row.names(ressum)[which(row.names(ressum) == betas[i])] <- var.namesb[i]
}
if (nrow(PCDataObj[["data"]][PCDataObj[["data"]][[outcome]] == 0, ]) > 0 & modelObj[["typeOfZero"]] == "true") {
for (i in 1:length(alphanames)) {
row.names(ressum)[which(row.names(ressum) == alphanames[i])] <- paste0(var.namesa[i], "_binary")
}
}
# create small results matrix
if (nrow(PCDataObj[["data"]][PCDataObj[["data"]][[outcome]] == 0, ]) == 0) {
ressum.slim <- ressum[which(row.names(ressum) %in% var.namesb), ]
if("categorical" %in% covariateTypesb){
ressum.slim <- ressum[which(row.names(ressum) %in% c(var.namesb, addbetanodes)), ]
}
}
if (nrow(PCDataObj[["data"]][PCDataObj[["data"]][[outcome]] == 0, ]) > 0 & modelObj[["typeOfZero"]] == "true") {
ressum.slim <- ressum[which(row.names(ressum) %in% c(var.namesb, paste0(var.namesa, "_binary"))), ]
if("categorical" %in% covariateTypesb){
ressum.slim <- ressum[which(row.names(ressum) %in% c(var.namesb, paste0(var.namesa, "_binary"), addbetanodes)), ]
}
}
# create mapping for coefficient names
coefnamemap<-matrix(nrow=length(var.namesb), ncol=2)
coefnamemap<-as.data.frame(coefnamemap)
colnames(coefnamemap)<-c('coefficient','variableName')
coefnamemap[,1]<-betanames
coefnamemap[,2]<-var.namesb
if("categorical" %in% covariateTypesb){
coefnametemp<-matrix(nrow=length(addbetanodes), ncol=2)
coefnametemp<-as.data.frame(coefnametemp)
colnames(coefnametemp)<-c('coefficient','variableName')
coefnametemp[,1]<-addbetanodes
coefnametemp[,2]<-addbetanodes
coefnamemap<-rbind(coefnamemap,coefnametemp)
}
if (nrow(PCDataObj[["data"]][PCDataObj[["data"]][[outcome]] == 0, ]) > 0 & modelObj[["typeOfZero"]] == "true"){
coefalphatemp<-matrix(nrow=length(var.namesa), ncol=2)
coefalphatemp<-as.data.frame(coefalphatemp)
colnames(coefalphatemp)<-c('coefficient','variableName')
coefalphatemp[,1]<-alphanames
coefalphatemp[,2]<-var.namesa
coefnamemap<-rbind(coefnamemap,coefalphatemp)
}
rm(llcal.censored, envir = .GlobalEnv)
outlist <- list()
outlist[["results"]] <- ressum.slim
if (keepModel == TRUE) {
outlist[["model"]] <- model
}
if (keepModelObj == TRUE) {
outlist[["modelObj"]] <- modMCMC
}
if (returnAllSummary == TRUE) {
outlist[["fullSummary"]] <- ressum
}
if (returnSample == TRUE) {
outlist[["sample"]] <- samp
}
outlist[['coefNameMap']]<-coefnamemap
outlist[["dataList"]] <- datl
outlist[["constantsList"]] <- const
outlist[['initsList']]<-inits
outlist[["convergenceTol"]] <- convergenceTol
class(outlist) <- 'PCResults'
return(outlist)
}
}
cammiller/imagingPC documentation built on June 28, 2019, 12:04 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
## code to create data for model and run model
#' Run the model
#'
#' The \code{runPCModel} function runs the model created by the
#' \code{\link{writePCModel}} function, using the information generated by
#' the \code{\link{createPCData}} function.
#'
#' @param modelObj An object of class \code{PCModelList}. This object is the
#' result of using the \code{\link{writePCModel}} function.
#' @param PCDataObj An object of class \code{PCDataList}. This object is the
#' result of using the \code{\link{createPCData}} function.
#' @param nIter An integer specifying he total number of iterations for the
#' MCMC, pre-thinning. This includes the burn-in.
#' @param nBurnin An integer specifying the number of iterations
#' (pre-thinning) to be discarded.
#' @param nThin An integer specifying the thinning interval.
#' @param slideVar An optional argument specifying a slide variable. This
#' option is only useful for data with zeros and when those zeros are assumed
#' to be censored observations. In this case, the censoring threshold can
#' vary by each run of the mass spectrometer. When this argument is used,
#' the censoring threshold will be allowed to vary by slide.
#' @param keepModel A TRUE/FALSE argument specifying whether or not to keep
#' the model (found in the modelObj).
#' @param keepModelObj A TRUE/FALSE argument specifying whether the
#' compiled NIMBLE model object should be kept after running the model. This
#' object can be large, but it is required for restarting the MCMC.
#' @param monitorCoefOnly A TRUE/FALSE argument stating whether to monitor the
#' the model coefficients only (TRUE) or all stochastic nodes (FALSE).
#' Regardless of the value of this argument, if any categorical covariates are
#' included in the model, the differences between all subgroups of the
#' categorical covariate(s) will be monitored in the MCMC with assigned nodes.
#' @param returnSample A TRUE/FALSE argument specifying whether to return the
#' entire MCMC sample. If this argument is FALSE, only the summary measures
#' will be returned.
#' @param returnAllSummary A TRUE/FALSE argument specifying whether to return
#' summary measures for all monitored nodes.
#' @param justReturnData A TRUE/FALSE argument. If justReturnData=TRUE, then
#' no MCMC will be run. Instead, the data and constants lists required to run
#' the model will be returned. This is useful for modifying the model to
#' personal specifications.
#' @param convergenceTol A numeric value indicating the highest acceptable
#' value for the Brooks-Gelman-Rubin (BGR) statistic to consider model
#' nodes converged.
#' @param maxRuns A numeric value indicating the maximum number of times a
#' model will be restarted to try to achieve convergence. After running the
#' MCMC for nIter iterations, the model is checked for convergence. A model
#' is considered converged if the model coefficients converge according to
#' the BGR statistic. For the marginalized two-part model, only the
#' coefficients for the marginalized part are used. If the model does not
#' converge, then the MCMC will be restarted and run for another
#' iterIncrement iterations (post-thinning). Convergence will then be
#' rechecked. The maxRuns argument specifies how many times the MCMC can
#' be restarted after the initial run of nIter iterations.
#' @param iterWindow A numeric value specifying the post-thinning number of
#' iterations (per chain) that will be used to make inference, if the model
#' does not converge with the first nIter iterations.
#' @param sliceSamplers A TRUE/FALSE argument specifying whether or not to
#' use slice samplers for all stochastic nodes. If
#' \code{sliceSamplers=TRUE} then NIMBLE's onlySlice option in the
#' \code{\link[nimble]{configureMCMC}} is set to TRUE. If
#' \code{sliceSamplers=FALSE} then NIMBLE's default MCMC settings are used.
#' Using slice samplers is useful for reducing autocorrelation in the MCMC,
#' though it significantly increases computation time.
#'
#' @return A list of class \code{PCResults} containing the results of running
#' the MCMC. The list will contain additional objects and information
#' specified by the arguments.
#'
#' \describe{
#' \item{\code{results}}{A matrix containing summary statistics (mean
#' and the 2.5%, 25%, 50%, 75%, and 97.5% percentiles) for each model
#' coefficient. For ease of reading, the names of the model coefficients
#' (beta1, beta2,...) are replaced with the names of their corresponding
#' covariates.}
#' \item{\code{model}}{The model specified from the
#' \code{\link{writePCModel}} function. The model is output if
#' \code{keepModel=TRUE}.}
#' \item{\code{modelObj}}{The compiled NIMBLE model object after running
#' until covergence or maxRuns is reached. This object is output only if
#' \code{keepModelObj=TRUE}. This object can be large (>100MB) for even
#' small imaging mass spectrometry datasets, so some thought should be
#' put into deciding whether or not to keep this object.}
#' \item{\code{fullSummary}}{A matrix of summary information for all
#' nodes in the model, This will be returned if
#' \code{returnAllSummary=TRUE}.}
#' \item{\code{sample}}{An \code{mcmc.list} object. See the \code{coda}
#' package for more details on such objects. This is a list of two
#' matrices, one for each chain. Each matrix has rows equal to the
#' number of thinned iterations and columns equal to the number of nodes
#' monitored. The sample is returned if \code{returnSample=TRUE}.}
#' \item{\code{dataList}}{The data list used to create a compiled
#' NIMBLE model. For all data this list will include the matrix
#' generated from the smoothing kernel function. For data with zeros
#' (censored and true), the model is specified using the zeros trick. In
#' this case \code{dataList} also has a vector of zeros.}
#' \item{\code{constantsList}}{The list of constants used to create
#' a compiled NIMBLE model. This list is typically composed of
#' covariates and other vectors needed to navigate the covariate
#' information. However, for data with zeros (censored and true),
#' the model is specified using the zeros trick, and so the outcome
#' data is also included in this list.}
#' \item{\code{convergenceTol}}{The limit on the BGR statistic used to
#' determine convergence. The default is 1.10, so if the BGR values for
#' all model coefficients are <=1.10, then the model is considered
#' converged.}
#' }
#'
#' @examples
#' data("TAMdata")
#' # The dataset is trimmed only for the speed of the example
#' TAMdata <- TAMdata[TAMdata$subject < 3, ]
#' TAMdata <- rScale(TAMdata, subjectVar = 'subject', sampleVar = 'ROI',
#' xCoord = 'x', yCoord = 'y')
#' rangs <- estRange(TAMdata, outcome = 'X1282.auc', spatialVar = 'TAM',
#' semivEst = 'modulus', logTransform = TRUE)
#' structs <- chooseStructures(rangs)
#' PCdat <- createPCData(structs, trimData = FALSE,
#' covariates = c("secondary", "TAM", "secTAM"),
#' covariateTypes = c("binary", "binary", "binary"),
#' covariateLevels = c("sample", "raster", "raster"))
#' PCmod <- writePCModel(PCdat, multiSampsPerSubj = TRUE, typeOfZero = "censored")
#' PCresults <- runPCModel(modelObj = PCmod, PCDataObj = PCdat, slideVar='slide',
#' monitorCoefOnly = FALSE,
#' nBurnin = 15000, nIter = 40000, nThin = 25)
#'
#' @references de Valpine, P., D. Turek, C.J. Paciorek, C. Anderson-Bergman,
#' D. Temple Lang, and R. Bodik. 2017. Programming with models: writing
#' statistical algorithms for general model structures with NIMBLE.
#' \emph{Journal of Computational and Graphical Statistics}. 26: 403-413.
runPCModel <- function(modelObj, PCDataObj, nIter = 50000, nBurnin = 10000, nThin = 5, slideVar = NULL, keepModel = TRUE, keepModelObj = FALSE,
monitorCoefOnly = TRUE, returnSample = TRUE, returnAllSummary = FALSE, justReturnData = FALSE, convergenceTol = 1.1, maxRuns = 5, iterIncrement = 500,
iterWindow = 1000, sliceSamplers=FALSE) {
if (class(modelObj) != "PCModelList")
stop("modelObj must of class PCModelList This object must be created using the writePCModel function.")
if (class(PCDataObj) != "PCDataList")
stop("PCDataObj must of class PCDataList. This object must be created using the createPCData function.")
if (is.numeric(nIter) == FALSE)
stop("nIter must be numeric.")
if ((nIter%%1 == 0) == FALSE)
stop("nIter must be an integer.")
if (nIter <= 0)
stop("nIter must be greater than 0.")
if (is.numeric(nBurnin) == FALSE)
stop("nBurnin must be numeric.")
if ((nBurnin%%1 == 0) == FALSE)
stop("nBurnin must be an integer.")
if (nBurnin < 0)
stop("nBurnin must be greater than or equal to 0.")
if (is.numeric(nThin) == FALSE)
stop("nThin must be numeric.")
if ((nThin%%1 == 0) == FALSE)
stop("nThin must be an integer.")
if (nThin <= 0)
stop("nThin must be greater than 0.")
if (is.numeric(convergenceTol) == FALSE)
stop("convergenceTol must be numeric.")
if (convergenceTol <= 1)
stop("convergenceTol should be greater than 1.")
if (is.numeric(maxRuns) == FALSE)
stop("maxRuns must be numeric.")
if ((maxRuns%%1 == 0) == FALSE)
stop("maxRuns must be an integer.")
if (maxRuns < 0)
stop("maxRuns cannot be negative.")
if (is.numeric(iterIncrement) == FALSE)
stop("iterIncrement must be numeric.")
if ((iterIncrement%%1 == 0) == FALSE)
stop("iterIncrement must be an integer.")
if (iterIncrement <= 0)
stop("iterIncrement must be greater than 0.")
if (is.numeric(iterWindow) == FALSE)
stop("iterWindow must be numeric.")
if ((iterWindow%%1 == 0) == FALSE)
stop("iterWindow must be an integer.")
if (iterWindow <= 0)
stop("iterWindow must be greater than 0.")
if (iterWindow > (((nIter - nBurnin)/nThin) + iterIncrement))
stop("iterWindow must be smaller than ((nIter-nBurnin)/nThin)+iterIncrement")
covs <- PCDataObj[["covs"]]
subjectVar <- PCDataObj[["subjectVar"]]
spatialVar <- PCDataObj[["spatialVar"]]
outcome <- PCDataObj[["outcome"]]
recStructures<-PCDataObj[['recStructures']]
model <- modelObj[["model"]]
covariates <- modelObj[["covariates"]]
covariateLevels<-modelObj[['covriateLevels']]
covariateTypes <- modelObj[["covariateTypes"]]
covariatesForBinary<-modelObj[['covariatesForBinary']]
rastpervar<-PCDataObj[['rastersPerVar']]
if(is.null(spatialVar)==FALSE){
#PCDataObj[['data']] <- PCDataObj[['data']][order(PCDataObj[['data']][[spatialVar]]), ]
spatvarlevels <- sort(unique(PCDataObj[['data']][[spatialVar]]))
nvarlevels <- length(spatvarlevels)
}
covstokeep <- rep(NA, length(covs))
for (i in 1:length(covs)) {
covstokeep[i] <- ifelse(covs[[i]]$covariate %in% covariates, 1, 0)
}
covstokeep <- which(covstokeep == 1)
covsb <- covs[covstokeep]
ordcovnames <- rep(NA, length(covs))
for (i in 1:length(covs)) {
ordcovnames[i] <- covs[[i]]$covariate
}
covariates <- factor(covariates, levels = ordcovnames)
covariates <- sort(covariates)
covariates <- as.character(covariates)
covariatesb <- covariates[covstokeep]
covariateTypesb <- covariateTypes[covstokeep]
covariateLevelsb <- covariateLevels[covstokeep]
nvars <- length(covsb)
covarnamesb <- list()
for (i in 1:nvars) {
if (covariateTypesb[i] %in% c("binary", "continuous")) {
covarnamesb[[i]] <- paste0(covariates[i])
}
if (covariateTypesb[i] == "categorical") {
covarnamesb[[i]] <- paste0(covariates[i], "_", colnames(covsb[[i]]$mapping))
}
}
var.namesb <- paste(unlist(covarnamesb))
if (nrow(PCDataObj[["data"]][PCDataObj[["data"]][[outcome]] == 0, ]) > 0 & modelObj[["typeOfZero"]] == "true") {
if (is.null(covariatesForBinary) == TRUE) {
covsa <- covsb
covariatesForBinary <- covariatesb
}
if (is.null(covariatesForBinary) == FALSE) {
covstokeep <- rep(NA, length(covs))
for (i in 1:length(covs)) {
covstokeep[i] <- ifelse(covs[[i]]$covariate %in% covariatesForBinary, 1, 0)
}
covstokeep <- which(covstokeep == 1)
covsa <- covs[covstokeep]
}
# need to get correct order of covariatesa
covariatesForBinary <- factor(covariatesForBinary, levels = ordcovnames)
covariatesForBinary <- sort(covariatesForBinary)
covariatesForBinary <- as.character(covariatesForBinary)
covariatesa <- covariatesForBinary[covstokeep]
covariateTypesa <- covariateTypes[covstokeep]
covariateLevelsa <- covariateLevels[covstokeep]
nvarsa <- length(covsa)
covarnamesa <- list()
for (i in 1:nvarsa) {
if (covariateTypesa[i] %in% c("binary", "continuous")) {
covarnamesa[[i]] <- paste0(covariatesForBinary[i])
}
if (covariateTypesa[i] == "categorical") {
covarnamesa[[i]] <- paste0(covariatesForBinary[i], "_", colnames(covsa[[i]]$mapping))
}
}
var.namesa <- paste(unlist(covarnamesa))
var.namesall <- unique(c(var.namesa, var.namesb))
}
if(is.null(spatialVar)==FALSE & any(recStructures==0)){
for(i in 1:length(unique(rastpervar$cov.level))){
assign(paste0('rastpervar',i), rastpervar[rastpervar$cov.level==spatvarlevels[i],])
}
}
############################################## Create data and constant lists #####
# create constants common to all models
const <- list()
for (i in 1:length(covsb)) {
if (covsb[[i]]$type == "categorical") {
for (j in 1:ncol(covsb[[i]]$mapping)) {
const[[paste0(covsb[[i]]$covariate, "_", colnames(covsb[[i]]$mapping)[j])]] <- covsb[[i]]$mapping[, j]
}
}
if (covsb[[i]]$type != "categorical") {
if (covsb[[i]]$level != "raster") {
const[[paste0(covsb[[i]]$covariate)]] <- covsb[[i]]$info[, 2]
}
if (covsb[[i]]$level == "raster") {
const[[paste0(covsb[[i]]$covariate)]] <- covsb[[i]]$info
}
}
}
if (nrow(PCDataObj[["data"]][PCDataObj[["data"]][[outcome]] == 0, ]) > 0 & modelObj[["typeOfZero"]] == "true") {
for (i in 1:length(covsa)) {
if (covsa[[i]]$type == "categorical") {
for (j in 1:ncol(covsa[[i]]$mapping)) {
const[[paste0(covsa[[i]]$covariate, "_", colnames(covsa[[i]]$mapping)[j])]] <- covsa[[i]]$mapping[, j]
}
}
if (covsa[[i]]$type != "categorical") {
if (covsa[[i]]$level != "raster") {
const[[paste0(covsa[[i]]$covariate)]] <- covsa[[i]]$info[, 2]
}
if (covsa[[i]]$level == "raster") {
const[[paste0(covsa[[i]]$covariate)]] <- covsa[[i]]$info
}
}
}
}
if (is.null(subjectVar) == FALSE) {
const[["Subj"]] <- PCDataObj[["nSubjs"]]
const[["Sample"]] <- PCDataObj[["cNSampsPerSubj"]]
const[["Raster"]] <- PCDataObj[["cNRastPerSamp"]]
}
if (is.null(subjectVar) == TRUE) {
const[["Sample"]] <- PCDataObj[["nSamps"]]
const[["Raster"]] <- PCDataObj[["cNRastPerSamp"]]
}
if (is.null(spatialVar) == FALSE) {
if(sum(PCDataObj[['recStructures']])==length(PCDataObj[['recStructures']])){
for (i in 1:PCDataObj$nVarLevels) {
const[[paste0("nsup", i)]] <- PCDataObj[["nSupportSites"]][, i]
}
}
if(any(PCDataObj[['recStructures']]==0)){
for (i in 1:PCDataObj$nVarLevels){
if(PCDataObj[['recStructures']][i]==1){
const[[paste0("nsup", i)]] <- PCDataObj[["nSupportSites"]][, i]
}
if(PCDataObj[['recStructures']][i]==0){
#create vector of start and stop
startstop<-rep(NA,(nrow(eval(parse(text = paste0('rastpervar',i))))*2))
for(j in 1:nrow(eval(parse(text = paste0('rastpervar',i))))){
startstop[(j-1)*2+1]<-eval(parse(text = paste0('rastpervar',i)))$start[j]
startstop[(j-1)*2+2]<-eval(parse(text = paste0('rastpervar',i)))$stop[j]
}
const[[paste0('samprast',i)]]<-startstop
}
}
}
}
if (is.null(spatialVar) == TRUE & all(PCDataObj[['recStructures']]==1)) {
const[["nsup"]] <- PCDataObj[["nSupportSites"]]
}
datl <- list()
# lognormal likelihood
if (nrow(PCDataObj[["data"]][PCDataObj[["data"]][[outcome]] == 0, ]) == 0) {
datl[["y"]] <- PCDataObj[["data"]][[outcome]]
}
# censored and true models
if (nrow(PCDataObj[["data"]][PCDataObj[["data"]][[outcome]] == 0, ]) > 0) {
const[["C"]] <- 10000
PCDataObj[["data"]]$zeros <- 0
datl[["zeros"]] <- PCDataObj[["data"]]$zeros
datl[["Kmat"]] <- PCDataObj[["KMat"]]
if (modelObj[["typeOfZero"]] == "censored")
{
PCDataObj[["data"]]$ind <- NA
for (i in 1:nrow(PCDataObj[["data"]])) {
PCDataObj[["data"]]$ind[i] <- ifelse(PCDataObj[["data"]][[outcome]][i] > 0, 1, 0)
}
const[["obs"]] <- PCDataObj[["data"]]$ind
if (is.null(slideVar) == TRUE) {
tstar <- min(PCDataObj[["data"]][PCDataObj[["data"]][[outcome]] > 0, outcome])
for (i in 1:nrow(PCDataObj[["data"]])) {
if (PCDataObj[["data"]][[outcome]][i] == 0) {
PCDataObj[["data"]][[outcome]][i] <- tstar
}
}
}
if (is.null(slideVar) == FALSE) {
tstarmat <- matrix(NA, nrow = length(unique(PCDataObj[["data"]][[slideVar]])), ncol = 2)
tstarmat <- as.data.frame(tstarmat)
colnames(tstarmat) <- c("slide", "tstar")
tstarmat$slide <- unique(PCDataObj[["data"]][[slideVar]])
for (i in 1:nrow(tstarmat)) {
tstarmat$tstar[i] <- min(PCDataObj[["data"]][[outcome]][PCDataObj[["data"]][[outcome]] > 0 & PCDataObj[["data"]][[slideVar]] == tstarmat$slide[i]])
}
for (i in 1:nrow(PCDataObj[["data"]])) {
if (PCDataObj[["data"]][[outcome]][i] == 0) {
PCDataObj[["data"]][[outcome]][i] <- tstarmat$tstar[tstarmat$slide == PCDataObj[["data"]][[slideVar]][i]]
}
}
}
const[["y"]] <- PCDataObj[["data"]][[outcome]]
# llcal.censored <- nimbleFunction(run = function(y = double(0), obs = double(0), sigma = double(0), mu = double(0)) {
# returnType(double(0))
# if (obs != 0) {
# return(log(dnorm((log(y) - mu), sd = sigma)))
# }
# if (obs == 0) {
# return(log(pnorm((log(y) - mu), sd = sigma)))
# }
# })
assign('llcal.censored',
eval(parse(text = c('nimbleFunction(run = function(y = double(0), obs = double(0), sigma = double(0), mu = double(0)) {',
'returnType(double(0))',
'if (obs != 0) {',
'return(log(dnorm((log(y) - mu), sd = sigma)))',
'}',
'if (obs == 0) {',
'return(log(pnorm((log(y) - mu), sd = sigma)))',
'}',
'})'))), envir=.GlobalEnv)
} # censored model
if (modelObj[["typeOfZero"]] == "true")
{
const[["y"]] <- PCDataObj[["data"]][[outcome]]
llcal.true <- nimbleFunction(run = function(y = double(0), binprob = double(0), sigma = double(0), mu = double(0)) {
returnType(double(0))
if (y != 0) {
return((log(binprob) - log(y) - 0.5 * log(2 * 3.14159265359) - log(sigma) - (1/(2 * pow(sigma, 2))) * pow((log(y) - mu), 2)))
}
if (y == 0) {
return(log(1 - binprob))
}
})
} # true zeros model
}
################################# Create inits list #####
inits <- list()
# inits common to all models
inits[["beta0"]] <- mean(log(PCDataObj[["data"]][[outcome]][PCDataObj[["data"]][[outcome]] > 0]))
betanames <- list()
for (i in 1:nvars) {
if (covariateTypesb[i] %in% c("binary", "continuous")) {
betanames[[i]] <- paste0("beta", i)
}
if (covariateTypesb[i] == "categorical") {
betanames[[i]] <- paste0("beta", i, "_", seq(1:ncol(covsb[[i]]$mapping)))
}
}
betanames <- paste(unlist(betanames))
for (i in 1:length(betanames)) {
inits[[paste0(betanames[i])]] <- 0
}
if (modelObj[["coefPrior"]] == "sdunif") {
inits[["sdb0"]] <- 0.1
sdbnames <- list()
for (i in 1:nvars) {
if (covariateTypesb[i] %in% c("binary", "continuous")) {
sdbnames[[i]] <- paste0("sdb", i)
}
if (covariateTypesb[i] == "categorical") {
sdbnames[[i]] <- paste0("sdb", i, "_", seq(1:ncol(covsb[[i]]$mapping)))
}
}
sdbnames <- paste(unlist(sdbnames))
for (i in 1:length(sdbnames)) {
inits[[paste0(sdbnames[i])]] <- 0.1
}
}
if (modelObj[["errorVarianceLevel"]] == "overall") {
inits[["sigma"]] <- 0.1
}
if (modelObj[["errorVarianceLevel"]] == "subject") {
inits[["sigma"]] <- rep(0.1, PCDataObj[["nSubjs"]])
}
if (modelObj[["errorVarianceLevel"]] == "sample") {
inits[["sigma"]] <- rep(0.1, PCDataObj[["nSamps"]])
}
if (modelObj[["multiSampsPerSubj"]] == TRUE & nrow(PCDataObj[["data"]][PCDataObj[["data"]][[outcome]] == 0, ]) > 0 & modelObj[["typeOfZero"]] ==
"censored") {
inits[["bi"]] <- rep(0, PCDataObj[["nSubjs"]])
inits[["sdbi"]] <- 0.1
}
if (modelObj[["multiSampsPerSubj"]] == TRUE & nrow(PCDataObj[["data"]][PCDataObj[["data"]][[outcome]] == 0, ]) == 0) {
inits[["bi"]] <- rep(0, PCDataObj[["nSubjs"]])
inits[["sdbi"]] <- 0.1
}
if (modelObj[["multiSampsPerSubj"]] == TRUE & nrow(PCDataObj[["data"]][PCDataObj[["data"]][[outcome]] == 0, ]) > 0 & modelObj[["typeOfZero"]] ==
"true") {
inits[["bi1"]] <- rep(0, PCDataObj[["nSubjs"]])
inits[["bi2"]] <- rep(0, PCDataObj[["nSubjs"]])
inits[["rho"]] <- 0
}
if (nrow(PCDataObj[["data"]][PCDataObj[["data"]][[outcome]] == 0, ]) > 0 & modelObj[["typeOfZero"]] == "true") {
alphanames <- list()
for (i in 1:nvarsa) {
if (covariateTypesa[i] %in% c("binary", "continuous")) {
alphanames[[i]] <- paste0("alpha", i)
}
if (covariateTypesa[i] == "categorical") {
alphanames[[i]] <- paste0("alpha", i, "_", seq(1:ncol(covsa[[i]]$mapping)))
}
}
alphanames <- paste(unlist(alphanames))
for (i in 1:length(alphanames)) {
inits[[paste0(alphanames[i])]] <- 0
}
if (modelObj[["coefPrior"]] == "sdunif") {
inits[["sda0"]] <- 0.1
sdanames <- list()
for (i in 1:nvars) {
if (covariateTypesa[i] %in% c("binary", "continuous")) {
sdanames[[i]] <- paste0("sda", i)
}
if (covariateTypesa[i] == "categorical") {
sdanames[[i]] <- paste0("sda", i, "_", seq(1:ncol(covsa[[i]]$mapping)))
}
}
sdanames <- paste(unlist(sdanames))
for (i in 1:length(sdanames)) {
inits[[paste0(sdanames[i])]] <- 0.1
}
}
inits[["v"]] <- rep(0, PCDataObj[["nSamps"]])
inits[["sdv"]] <- 0.1
}
if (is.null(spatialVar) == TRUE & all(PCDataObj[['recStructures']]==1)) {
xc <- matrix(nrow = PCDataObj[["nSamps"]], ncol = max(PCDataObj[["nSupportSites"]]))
for (i in 1:PCDataObj[["nSamps"]]) {
xc[i, 1:PCDataObj[["nSupportSites"]][i]] <- 0
}
inits[["xc"]] <- xc
if (modelObj[["latentVarianceLevel"]] == "overall") {
inits[["sdxc"]] <- 0.1
}
if (modelObj[["latentVarianceLevel"]] == "subject") {
inits[["sdxc"]] <- rep(0.1, PCDataObj[["nSubjs"]])
}
if (modelObj[["latentVarianceLevel"]] == "sample") {
inits[["sdxc"]] <- rep(0.1, PCDataObj[["nSamps"]])
}
}
if (is.null(spatialVar) == FALSE) {
for (j in 1:PCDataObj[["nVarLevels"]]) {
if(PCDataObj[['recStructures']][j]==1){
assign(paste0("xc", j), matrix(NA, nrow = length(PCDataObj[["nSupportSites"]][[paste0("nsups.", j)]]), ncol = max(PCDataObj[["nSupportSites"]][[paste0("nsups.",
j)]])))
eval(parse(text = c(paste0("for(i in 1:length(PCDataObj[[\"gT0SupportSites\"]][[", j, "]])){"), paste0("xc", j, "[PCDataObj[[\"gT0SupportSites\"]][[",
j, "]][i],1:PCDataObj[[\"nSupportSites\"]]$nsups.", j, "[PCDataObj[[\"gT0SupportSites\"]][[1]][i]]]<-0"), "}")))
eval(parse(text = paste0("inits[[\"xc", j, "\"]]<-xc", j)))
if (modelObj[["latentVarianceLevel"]] == "overall") {
assign(paste0("sdxc", j), 0.1)
eval(parse(text = paste0("inits[[\"sdxc", j, "\"]]<-sdxc", j)))
}
if (modelObj[["latentVarianceLevel"]] == "subject") {
assign(paste0("sdxc", j), rep(0.1, PCDataObj[["nSubjs"]]))
eval(parse(text = paste0("inits[[\"sdxc", j, "\"]]<-sdxc", j)))
}
if (modelObj[["latentVarianceLevel"]] == "sample") {
assign(paste0("sdxc", j), rep(0.1, PCDataObj[["nSamps"]]))
eval(parse(text = paste0("inits[[\"sdxc", j, "\"]]<-sdxc", j)))
}
}
if(PCDataObj[['recStructures']][j]==0){
eval(parse(text = paste0('inits[["rk',j,'"]]<-rep(NA,max(rastpervar$stop))')))
eval(parse(text = c(paste0('for(i in 1:nrow(rastpervar',j,')){'),paste0('inits[["rk',j,'"]][rastpervar',j,'$start[i]:rastpervar',j,'$stop[i]]<-0'),'}')))
eval(parse(text = paste0('inits[["sdrk',j,'"]]<-0.1')))
}
}
}
######################## run PC model #####
if (justReturnData == TRUE) {
outlist <- list()
if (keepModel == TRUE) {
outlist[["model"]] <- model
}
outlist[["dataList"]] <- datl
outlist[["constantsList"]] <- const
outlist[['initsList']]<-inits
outlist[["convergenceTol"]] <- convergenceTol
return(outlist)
}
if (justReturnData == FALSE) {
if ("categorical" %in% covariateTypesb) {
addbetanodes <- list()
nodetemp <- list()
for (i in 1:nvars) {
if (covariateTypesb[i] == "categorical") {
# to provide names, need to get mapping matrix, column of betas and column of corresponding name
betanames2 <- paste0("beta", i, "_", seq(1:ncol(covsb[[i]]$mapping)))
covarnames2 <- paste0(covariates[i], "_", colnames(covsb[[i]]$mapping))
namemap <- data.frame(betanames2, covarnames2)
nodetemp[[i]] <- paste0("beta", i, "_", seq(1:ncol(covsb[[i]]$mapping)))
nodemat <- expand.grid(paste(nodetemp[[i]]), paste(nodetemp[[i]]))
nodemat <- nodemat[nodemat[, 1] != nodemat[, 2], ]
nodemat[, 3] <- 0
for (j in 1:nrow(nodemat)) {
if (as.numeric(strsplit(as.character(nodemat[j, 1]), "_")[[1]][2]) > as.numeric(strsplit(as.character(nodemat[j, 2]), "_")[[1]][2])) {
nodemat[j, 3] <- 1
}
}
nodemat <- nodemat[nodemat[, 3] == 1, ]
nodemat <- nodemat[order(nodemat[, 1]), ]
nodemat <- nodemat[, c(1, 2)]
nodemat <- as.data.frame(nodemat)
colnames(nodemat) <- c("beta1", "beta2")
nodemat$cov1 <- NA
nodemat$cov2 <- NA
for (j in 1:nrow(nodemat)) {
nodemat$cov1[j] <- as.character(namemap$covarnames2[namemap$betanames2 == nodemat$beta1[j]])
nodemat$cov2[j] <- as.character(namemap$covarnames2[namemap$betanames2 == nodemat$beta2[j]])
}
# need to subtract columns
addbetanodes[[i]] <- rep(NA, nrow(nodemat))
for (j in 1:nrow(nodemat)) {
addbetanodes[[i]][j] <- paste0(nodemat$cov1[j], "_minus_", nodemat$cov2[j])
}
}
}
addbetanodes <- paste(unlist(addbetanodes))
}
# vector of betas
betas <- betanames
if (nrow(PCDataObj[["data"]][PCDataObj[["data"]][[outcome]] == 0, ]) > 0 & modelObj[["typeOfZero"]] == "true") {
betas<-c(betanames,alphanames)
}
nimmod <- nimbleModel(code = model, constants = const, data = datl, inits = inits)
cnimmod <- compileNimble(nimmod)
if (monitorCoefOnly == TRUE) {
confmod <- configureMCMC(nimmod, monitors = betas, print = FALSE, onlySlice = sliceSamplers)
}
if (monitorCoefOnly == FALSE) {
confmod <- configureMCMC(nimmod, monitors = cnimmod$getNodeNames(stochOnly = TRUE, includeData = FALSE), print = FALSE, onlySlice = sliceSamplers)
}
if (any(covariateTypesb == "categorical")) {
confmod$addMonitors(addbetanodes)
}
bmod <- buildMCMC(confmod)
modMCMC <- compileNimble(bmod, project = nimmod)
samp <- runMCMC(modMCMC, niter = nIter, nburnin = nBurnin, nchains = 2, samplesAsCodaMCMC = TRUE, thin = nThin)
samp[[1]] <- samp[[1]][, is.na(samp[[1]][1, ]) == FALSE]
samp[[2]] <- samp[[2]][, is.na(samp[[2]][1, ]) == FALSE]
samp <- as.mcmc.list(samp)
bgrs <- gelman.diag(samp, autoburnin = FALSE, multivariate = FALSE)
bgrtemp <- bgrs$psrf[betas, 1]
if (max(bgrtemp) <= convergenceTol | maxRuns==0) {
if(max(bgrtemp) <= convergenceTol){
cat('The model coefficients converged.\n')
} else {
cat('The model coefficients did not converge.\n')
}
}
if (max(bgrtemp) > convergenceTol) {
cat("With a burnin of", nBurnin, "and a sample of", nrow(samp[[1]]), "(post-thinning), the model coefficients did not converge. Will keep running to try to reach convergence.\n")
iter <- 1
flag <- 1
while (flag) {
sampnew <- runMCMC(modMCMC, niter = iterIncrement * nThin, nburnin = 0, nchains = 2, samplesAsCodaMCMC = TRUE, thin = nThin)
sampnew[[1]] <- sampnew[[1]][, is.na(sampnew[[1]][1, ]) == FALSE]
sampnew[[2]] <- sampnew[[2]][, is.na(sampnew[[2]][1, ]) == FALSE]
sampnew <- as.mcmc.list(sampnew)
samp[[1]] <- rbind(samp[[1]], sampnew[[1]])
samp[[2]] <- rbind(samp[[2]], sampnew[[2]])
samp[[1]] <- samp[[1]][(nrow(samp[[1]]) - (iterWindow - 1)):nrow(samp[[1]]), ]
samp[[2]] <- samp[[2]][(nrow(samp[[2]]) - (iterWindow - 1)):nrow(samp[[2]]), ]
samp[[1]] <- as.mcmc(samp[[1]])
samp[[2]] <- as.mcmc(samp[[2]])
samp <- as.mcmc.list(samp)
bgrs <- gelman.diag(samp, autoburnin = FALSE, multivariate = FALSE)
bgrtemp <- bgrs$psrf[betas, 1]
if (max(bgrtemp) <= convergenceTol | iter >= maxRuns) {
flag <- 0
}
if (max(bgrtemp) > convergenceTol) {
cat("Run", iter, "complete without convergence.")
}
if (iter == maxRuns) {
cat("Max number of runs reached.")
}
iter <- iter + 1
}
}
sampres <- summary(samp)
sampstat <- sampres$statistics
sampquant <- sampres$quantiles
ressum <- cbind(sampstat[, 1], sampquant)
colnames(ressum) <- c("Mean", "2.5%", "25%", "50%", "75%", "97.5%")
# row.names(ressum) which(row.names(ressum)==betas[1])
# reassign row names
for (i in 1:length(betanames)) {
row.names(ressum)[which(row.names(ressum) == betas[i])] <- var.namesb[i]
}
if (nrow(PCDataObj[["data"]][PCDataObj[["data"]][[outcome]] == 0, ]) > 0 & modelObj[["typeOfZero"]] == "true") {
for (i in 1:length(alphanames)) {
row.names(ressum)[which(row.names(ressum) == alphanames[i])] <- paste0(var.namesa[i], "_binary")
}
}
# create small results matrix
if (nrow(PCDataObj[["data"]][PCDataObj[["data"]][[outcome]] == 0, ]) == 0) {
ressum.slim <- ressum[which(row.names(ressum) %in% var.namesb), ]
if("categorical" %in% covariateTypesb){
ressum.slim <- ressum[which(row.names(ressum) %in% c(var.namesb, addbetanodes)), ]
}
}
if (nrow(PCDataObj[["data"]][PCDataObj[["data"]][[outcome]] == 0, ]) > 0 & modelObj[["typeOfZero"]] == "true") {
ressum.slim <- ressum[which(row.names(ressum) %in% c(var.namesb, paste0(var.namesa, "_binary"))), ]
if("categorical" %in% covariateTypesb){
ressum.slim <- ressum[which(row.names(ressum) %in% c(var.namesb, paste0(var.namesa, "_binary"), addbetanodes)), ]
}
}
# create mapping for coefficient names
coefnamemap<-matrix(nrow=length(var.namesb), ncol=2)
coefnamemap<-as.data.frame(coefnamemap)
colnames(coefnamemap)<-c('coefficient','variableName')
coefnamemap[,1]<-betanames
coefnamemap[,2]<-var.namesb
if("categorical" %in% covariateTypesb){
coefnametemp<-matrix(nrow=length(addbetanodes), ncol=2)
coefnametemp<-as.data.frame(coefnametemp)
colnames(coefnametemp)<-c('coefficient','variableName')
coefnametemp[,1]<-addbetanodes
coefnametemp[,2]<-addbetanodes
coefnamemap<-rbind(coefnamemap,coefnametemp)
}
if (nrow(PCDataObj[["data"]][PCDataObj[["data"]][[outcome]] == 0, ]) > 0 & modelObj[["typeOfZero"]] == "true"){
coefalphatemp<-matrix(nrow=length(var.namesa), ncol=2)
coefalphatemp<-as.data.frame(coefalphatemp)
colnames(coefalphatemp)<-c('coefficient','variableName')
coefalphatemp[,1]<-alphanames
coefalphatemp[,2]<-var.namesa
coefnamemap<-rbind(coefnamemap,coefalphatemp)
}
rm(llcal.censored, envir = .GlobalEnv)
outlist <- list()
outlist[["results"]] <- ressum.slim
if (keepModel == TRUE) {
outlist[["model"]] <- model
}
if (keepModelObj == TRUE) {
outlist[["modelObj"]] <- modMCMC
}
if (returnAllSummary == TRUE) {
outlist[["fullSummary"]] <- ressum
}
if (returnSample == TRUE) {
outlist[["sample"]] <- samp
}
outlist[['coefNameMap']]<-coefnamemap
outlist[["dataList"]] <- datl
outlist[["constantsList"]] <- const
outlist[['initsList']]<-inits
outlist[["convergenceTol"]] <- convergenceTol
class(outlist) <- 'PCResults'
return(outlist)
}
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.