Nothing
make.jagsboralnullmodel <- function (family, num.X = 0, X.ind = NULL, num.traits = 0, which.traits = NULL,
row.eff = "none", row.ids = NULL, ranef.ids = NULL,
offset = NULL, trial.size = 1, n, p, model.name = NULL,
prior.control = list(type = c("normal","normal","normal","uniform"), hypparams = c(10, 10, 10, 30),
ssvs.index = -1, ssvs.g = 1e-6, ssvs.traitsindex = -1)) {
check_which_traits(num.traits = num.traits, which.traits = which.traits, num.X = num.X, makejagsboralfile_messages = TRUE)
if(is.null(which.traits)) {
which.traits <- vector("list",num.X+1)
for(k in 1:(num.X+1))
which.traits[[k]] <- 0
}
complete_trial_size <- check_trial_size(family = family, trial.size = trial.size, create.complete.trial.size = TRUE, y = matrix(NA,nrow=1,ncol=p))
complete_family <- check_family(family = family, y = matrix(1,nrow=1,ncol=p), traits = NULL) ## Done just to produce the complete_family vector
if(row.eff != "none" && is.null(row.ids)) {
row.ids <- matrix(1:n, ncol = 1)
message("row.ids assumed to be matrix with one column and elements 1,2,...n i.e., a row-specific intercept.")
}
if(!is.null(row.ids)) {
row.ids <- as.matrix(row.ids)
if(nrow(row.ids) != n)
stop("Number of rows in the matrix row.ids should be equal to n.")
if(is.null(colnames(row.ids)))
colnames(row.ids) <- paste0("ID", 1:ncol(row.ids))
}
ranef.ids <- check_ranef_ids(ranef.ids = ranef.ids, y = matrix(0,n,p))
if(!is.null(offset)) {
if(!is.matrix(offset))
stop("offset could be a matrix with the same dimensions as y.")
}
prior.control <- fillin_prior_control(x = prior.control)
check_prior_control(prior.control = prior.control)
if(length(prior.control$ssvs.index) == 1 & num.X > 0)
prior.control$ssvs.index <- rep(prior.control$ssvs.index, num.X)
if(num.traits > 0) {
if(!is.list(prior.control$ssvs.traitsindex)) {
prior.control$ssvs.traitsindex <- vector("list",num.X+1)
for(k in 1:(num.X+1))
prior.control$ssvs.traitsindex[[k]] <- rep(-1,length(which.traits[[k]]))
}
if(is.list(prior.control$ssvs.traitsindex))
{
check_ssvstraits(prior.control$ssvs.traitsindex, which.traits)
}
}
X.ind <- check_X_ind(X.ind = X.ind, p = p, num.X = num.X, prior.control = prior.control)
index.ord.cols <- which(complete_family == "ordinal")
index.tweed.cols <- which(complete_family == "tweedie")
##-------------------
## Checks done; starting writing JAGS script!
##-------------------
model_script <- paste0("## JAGS model written for boral version ", packageDescription("boral")$Version, " on ", as.character(Sys.time()), " ##\n\n model {")
model_script <- c(model_script, "\t ## Data Level ## \n\t for(i in 1:n) {")
write.resp.script <- setup_respfamilies(p = p, complete.family = complete_family, num.lv = 0, row.eff = row.eff, row.ids = row.ids, ranef.ids = ranef.ids,
offset = offset, num.X = num.X, complete.trial.size = complete_trial_size, index.tweed.cols = index.tweed.cols, index.ord.cols = index.ord.cols)
model_script <- c(model_script, write.resp.script)
model_script <- c(model_script, paste0("\t\t } \n\n\t ## Process level and priors ##"))
rm(write.resp.script)
## Build prior strings for all priors distributions
prior.strings <- construct_prior_strings(x = prior.control)
## Code for response-specific intercept. Note this is set up different to how X variables are set up to save some coding space!
## No traits or traits included but not regressed against intercept
if(num.traits == 0 || (num.traits > 0 & which.traits[[1]][1] == 0)) {
## Not ordinal columns, then as per usual
if(length(index.ord.cols) == 0)
model_script <- c(model_script, paste0("\t for(j in 1:p) { lv.coefs[j,1] ~ ", prior.strings$p1, " } ## Separate response intercepts"))
## If 1 ordinal column, then intercept for this column equal 0
if(length(index.ord.cols) == 1) {
model_script <- c(model_script, paste0("\t lv.coefs[",index.ord.cols, ",1] <- 0 ## Single ordinal response intercept"))
for(j in (1:p)[-index.ord.cols])
model_script <- c(model_script, paste0("\t lv.coefs[", j, ",1] ~ ", prior.strings$p1, "All other response intercepts"))
}
## More than 1 ordinal column, then set up random intercept for this response
if(length(index.ord.cols) > 1) {
if(length(index.ord.cols) == p)
model_script <- c(model_script, paste0("\t for(j in 1:p) { lv.coefs[j,1] ~ dnorm(0,pow(ordinal.sigma,-2)) } ## Random intercept for all ordinal response"))
else
{
for(j in index.ord.cols)
model_script <- c(model_script, paste0("\t lv.coefs[",j, ",1] ~ dnorm(0,pow(ordinal.sigma,-2)) ## Random intercept for all ordinal response"))
for(j in (1:p)[-index.ord.cols])
model_script <- c(model_script, paste0("\t lv.coefs[", j, ",1] ~ ", prior.strings$p1, "All other response intercepts"))
}
model_script <- c(model_script, paste0("\t ordinal.sigma ~ ", prior.strings$p4))
}
if((num.traits > 0 & which.traits[[1]][1] == 0)) {
model_script <- c(model_script, paste0("\t traits.int[1] <- 0; for(l in 1:num.traits) { traits.coefs[1,l] <- 0 } \n\t trait.sigma[1] <- 0 ## Traits not used for intercept"))
}
}
## Traits included and intercept regressed against them
if(num.traits > 0 & all(which.traits[[1]] > 0)) {
## If there are 0 or > 1 ordinal columns, then regress all intercepts against traits
if(length(index.ord.cols) != 1) {
model_script <- c(model_script, paste0("\t for(j in 1:p) { lv.coefs[j,1] ~ dnorm(traits.int[1] + inprod(traits[j,],traits.coefs[1,1:num.traits]),pow(trait.sigma[1],-2)) } ## Response intercepts regressed against traits"))
}
## If there is 1 ordinal column, do not regress this intercept against trait
if(length(index.ord.cols) == 1) {
model_script <- c(model_script, paste0("\t lv.coefs[",index.ord.cols, ",1] <- 0 ## Ordinal response intercept"))
for(j in (1:p)[-index.ord.cols])
model_script <- c(model_script, paste0("\t lv.coefs[", j, ",1] ~ dnorm(traits.int[1] + inprod(traits[",j,",],traits.coefs[1,1:num.traits]),pow(trait.sigma[1],-2)) ## All other intercepts"))
}
for(l in which.traits[[1]]) {
if(prior.control$ssvs.traitsindex[[1]][which.traits[[1]] == l] == -1)
model_script <- c(model_script, paste0("\t traits.coefs[",1,",",l,"] ~ ", prior.strings$p1, " ## Traits used for intercept"))
if(prior.control$ssvs.traitsindex[[1]][which.traits[[1]] == l] == 0)
{
ssvs.prior.string <- paste0("dnorm(0,pow(", prior.control$hypparams[1], "*((1-ssvs.traitscoefs1",l,")*", prior.control$ssvs.g, " + ssvs.traitscoefs1", l, "),-1)); ssvs.traitscoefs1", l, " ~ dbern(0.5)")
model_script <- c(model_script, paste0("\t traits.coefs[",1,",",l,"] ~ ", ssvs.prior.string, " ## Traits used for intercept"))
}
}
if(length((1:num.traits)[-which.traits[[1]]]) > 0) {
for(l in (1:num.traits)[-which.traits[[1]]])
{
model_script <- c(model_script, paste0("\t traits.coefs[",1,",",l,"] <- 0 ## Traits not used for intercept"))
}
}
model_script <- c(model_script, paste0("\t traits.int[1] ~ ", prior.strings$p1))
model_script <- c(model_script, paste0("\t trait.sigma[1] ~ ", prior.strings$p4))
}
if(any(complete_family == "tweedie"))
model_script <- c(model_script, paste0("\t powerparam ~ dunif(1,2) ## Tweedie power parameter"))
if(any(complete_family == "ordinal")) {
model_script <- c(model_script, paste0("\t for(k in 1:(num.ord.levels-1)) { cutoffs0[k] ~ ", prior.strings$p1, " }"))
model_script <- c(model_script, paste0("\t cutoffs[1:(num.ord.levels-1)] <- sort(cutoffs0) ## Ordinal cutoffs"))
}
## Priors on row effects
if(row.eff == "fixed") {
for(k in 1:ncol(row.ids)) {
model_script <- c(model_script, paste0("\n\t row.coefs.ID", k, "[1] <- 0"))
model_script <- c(model_script, paste0("\n\t for(i in 2:n.ID[", k, "]) { row.coefs.ID", k, "[i] ~ ", prior.strings$p1, " } "))
}
}
if(row.eff == "random") {
for(k in 1:ncol(row.ids)) {
model_script <- c(model_script, paste0("\n\t for(i in 1:n.ID[", k, "]) { row.coefs.ID", k, "[i] ~ dnorm(0, pow(row.sigma.ID", k, ",-2)) } "))
model_script <- c(model_script, paste0("\t row.sigma.ID", k, " ~ ",prior.strings$p4))
}
#model_script <- c(model_script, paste0("\t row.ranef.mean ~ ", prior.strings$p1))
}
## Priors on response-specific random intercepts
if(!is.null(ranef.ids)) {
for(k0 in 1:ncol(ranef.ids)) {
model_script <- c(model_script, paste0("\n\t for(j in 1:p) { for(i in 1:n.ranefID[", k0, "]) { ranef.coefs.ID", k0, "[j,i] ~ dnorm(0, pow(ranef.sigma.ID", k0, "[j],-2)) } }"))
model_script <- c(model_script, paste0("\t for(j in 1:p) { ranef.sigma.ID", k0, "[j] ~ ",prior.strings$p4, " }"))
}
}
## Prior for X-coefficients, controlled by prior.control$hypparams[3]
if(num.X > 0) {
model_script <- c(model_script, paste0("\n"))
## Traits not included in model
if(num.traits == 0) {
for(i in 1:length(prior.control$ssvs.index)) {
if(prior.control$ssvs.index[i] == -1) {
if(!is.null(X.ind))
model_script <- c(model_script, paste0("\t for(j in 1:p) { X.coefs[j,", i, "] ~ ", prior.strings$p3, "I(-X.ind[j,",i,"],X.ind[j,",i,"]) } "))
if(is.null(X.ind))
model_script <- c(model_script, paste0("\t for(j in 1:p) { X.coefs[j,", i, "] ~ ", prior.strings$p3, " } "))
}
if(prior.control$ssvs.index[i] == 0) {
ssvs.prior.string <- paste0("dnorm(0,pow(", prior.control$hypparams[3], "*((1-ssvs.indX", i, "[j])*", prior.control$ssvs.g, " + ssvs.indX", i, "[j]),-1)); ssvs.indX", i, "[j] ~ dbern(0.5)")
model_script <- c(model_script, paste0("\t for(j in 1:p) { X.coefs[j,", i, "] ~ ", ssvs.prior.string, " }"))
}
if(prior.control$ssvs.index[i] > 0) {
ssvs.prior.string <- paste0("dnorm(0,pow(", prior.control$hypparams[3], "*((1-ssvs.gp", prior.control$ssvs.index[i], ")*", prior.control$ssvs.g, " + ssvs.gp", prior.control$ssvs.index[i], "),-1))")
model_script <- c(model_script, paste0("\t for(j in 1:p) { X.coefs[j,", i, "] ~ ", ssvs.prior.string, " } "))
}
}
}
if(num.traits > 0) { for(i in 1:num.X) {
## Traits included but X coefs not regressed against them
if(all(which.traits[[i+1]] == 0)) {
model_script <- c(model_script, paste0("\t for(j in 1:p) { X.coefs[j,", i, "] ~ ", prior.strings$p3, " } ## Coefficient not regressed against any traits"))
model_script <- c(model_script, paste0("\t traits.int[",i+1,"] <- 0; trait.sigma[",i+1,"] <- 0; for(l in 1:num.traits) { traits.coefs[",i+1,",l] <- 0 } \n"))
}
## Traits included and X coefs regressed against some of them
if(all(which.traits[[i+1]] > 0)) {
model_script <- c(model_script, paste0("\t for(j in 1:p) { X.coefs[j,", i, "] ~ dnorm(traits.int[",i+1,"] + inprod(traits[j,],traits.coefs[",i+1,",1:num.traits]),pow(trait.sigma[",i+1,"],-2)) } "))
for(l in which.traits[[i+1]]) {
if(prior.control$ssvs.traitsindex[[i+1]][which.traits[[i+1]] == l] == -1)
model_script <- c(model_script, paste0("\t traits.coefs[",i+1,",",l,"] ~ ", prior.strings$p3, " ## Traits used for this X.coefs"))
if(prior.control$ssvs.traitsindex[[i+1]][which.traits[[i+1]] == l] == 0)
{
ssvs.prior.string <- paste0("dnorm(0,pow(", prior.control$hypparams[3], "*((1-ssvs.traitscoefs",i+1,l,")*", prior.control$ssvs.g, " + ssvs.traitscoefs", i+1, l, "),-1)); ssvs.traitscoefs", i+1, l, " ~ dbern(0.5)")
model_script <- c(model_script, paste0("\t traits.coefs[",i+1,",",l,"] ~ ", ssvs.prior.string, " ## Traits used for this X.coefs"))
}
}
if(length((1:num.traits)[-which.traits[[i+1]]]) > 0) {
for(l in (1:num.traits)[-which.traits[[i+1]]]) {
model_script <- c(model_script, paste0("\t traits.coefs[",i+1,",",l,"] <- 0 ## traits not used for this X.coefs"))
}
}
model_script <- c(model_script, paste0("\t traits.int[",i+1,"] ~ ", prior.strings$p3))
model_script <- c(model_script, paste0("\t trait.sigma[",i+1,"] ~ ", prior.strings$p4, "\n"))
}
} }
model_script <- c(model_script, paste0(""))
if(any(prior.control$ssvs.index > 0)) {
for(i in unique(prior.control$ssvs.index[prior.control$ssvs.index > 0]))
model_script <- c(model_script, paste0("\t ssvs.gp", i, " ~ dbern(0.5)"))
}
}
# if(num.X > 0 & any(family == "multinom")) {
# model_script <- c(model_script, paste0("\t for(j in 1:",length(index_multinom_cols),") { for(i in 1:num.X) { X.multinom.params[j,i,1] <- 0 } } "))
# model_script <- c(model_script, paste0("\t for(j in 1:",length(index_multinom_cols),") { for(i in 1:num.X) { for(k in 2:num.multinom.levels) { X.multinom.params[j,i,k] ~ dnorm(0,",1/prior.control$hypparams[1],") } } } "))
# }
## Prior on dispersion parameters, controlled by prior.control$hypparams[4]
if(!all(complete_family %in% c("poisson", "ztpoisson", "binomial", "ordinal", "multinom", "exponential"))) {
model_script <- c(model_script, paste0("\t for(j in 1:p) { lv.coefs[j,num.lv+2] ~ ", prior.strings$p4, " } ## Dispersion parameters"))
}
model_script <- c(model_script, "\n\t }")
if(!is.null(model.name)) { write(model_script, file = model.name) }
if(is.null(model.name)) { write(model_script, file = "jagsboralmodel.txt") }
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.