diff --git a/DESCRIPTION b/DESCRIPTION index 3dbe7a4..3ccbe3f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -12,4 +12,5 @@ LazyData: true RoxygenNote: 7.3.2 Imports: vcrpart, - rjags + rjags, + dclone diff --git a/NAMESPACE b/NAMESPACE index 914eeea..a765abf 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -4,5 +4,6 @@ export(bpcm) export(pcm) export(res_ij) export(residif) +import(dclone) import(rjags) import(vcrpart) diff --git a/R/bpcm.R b/R/bpcm.R index 8acca22..43b4435 100644 --- a/R/bpcm.R +++ b/R/bpcm.R @@ -14,6 +14,7 @@ #' @param diag.plots boolean indicating whether the JAGS diagnosis plots should be displayed #' @return A data.frame containing various model outputs #' @import rjags +#' @import dclone #' @export bpcm <- function(df=NULL,grp=NULL,is.dif=NULL,is.unif=NULL,priors=NULL,param=list(),verbose=T,diag.plots=F) { @@ -121,7 +122,7 @@ bpcm <- function(df=NULL,grp=NULL,is.dif=NULL,is.unif=NULL,priors=NULL,param=lis if (verbose) { cat('Initializing Chains...') } - dclone::parJagsModel(cl,name="mod.jags",data = data,inits = NULL,file = "/home/corentin/Documents/These/Recherche/SPT/src/bpcm",n.chains=n.chains) + dclone::parJagsModel(cl,name="mod.jags",data = data,inits = NULL,file = system.file("jags","bpcm",package="SPT"),n.chains=n.chains) if (verbose) { cat('DONE \n') cat('Applying burn-in iterations...') @@ -231,7 +232,7 @@ bpcm <- function(df=NULL,grp=NULL,is.dif=NULL,is.unif=NULL,priors=NULL,param=lis if (verbose) { cat('Initializing Chains...') } - dclone::parJagsModel(cl,name="mod.jags",data = data,inits = NULL,file = "/home/corentin/Documents/These/Recherche/SPT/src/bpcm_beta",n.chains=n.chains) + dclone::parJagsModel(cl,name="mod.jags",data = data,inits = NULL,file = system.file("jags","bpcm_beta",package="SPT"),n.chains=n.chains) if (verbose) { cat('DONE \n') cat('Applying burn-in iterations...') @@ -356,7 +357,7 @@ bpcm <- function(df=NULL,grp=NULL,is.dif=NULL,is.unif=NULL,priors=NULL,param=lis if (verbose) { cat('Initializing Chains...') } - dclone::parJagsModel(cl,name="mod.jags",data = data,inits = NULL,file = "/home/corentin/Documents/These/Recherche/SPT/src/bpcm_dif",n.chains=n.chains) + dclone::parJagsModel(cl,name="mod.jags",data = data,inits = NULL,file = system.file("jags","bpcm_dif",package="SPT"),n.chains=n.chains) if (verbose) { cat('DONE \n') cat('Applying burn-in iterations...') diff --git a/inst/jags/bpcm b/inst/jags/bpcm new file mode 100644 index 0000000..8706406 --- /dev/null +++ b/inst/jags/bpcm @@ -0,0 +1,23 @@ +model { + for (i in 1:n) { + for (j in 1:p) { + Y[i,j] ~ dcat(prob[i,j,1:K[j]]) + } + theta[i] ~ dnorm(0,1) + } + for (i in 1:n) { + for (j in 1:p) { + for (k in 1:K[j] ) { + eta[i,j,k] <- theta[i] - delta[j,k] + psum[i,j,k] <- sum(eta[i,j,1:k]) + exp.psum[i,j,k] <- exp(psum[i,j,k]) + prob[i,j,k] <- exp.psum[i,j,k] / sum(exp.psum[i,j,1:K[j]]) + } } } + for (j in 1:p) { + delta[j,1] <- 0.0 + for (k in 2:K[j]) { + delta[j,k] ~ dnorm(m.delta, pr.delta) + } + } + pr.delta <- pow(s.delta, -2) +} diff --git a/inst/jags/bpcm_beta b/inst/jags/bpcm_beta new file mode 100644 index 0000000..ffc441e --- /dev/null +++ b/inst/jags/bpcm_beta @@ -0,0 +1,25 @@ +model { + for (i in 1:n) { + for (j in 1:p) { + Y[i,j] ~ dcat(prob[i,j,1:K[j]]) + } + theta[i] ~ dnorm(0,1) + } + for (i in 1:n) { + for (j in 1:p) { + for (k in 1:K[j] ) { + eta[i,j,k] <- theta[i] + (beta * Z[i,1]) - delta[j,k] + psum[i,j,k] <- sum(eta[i,j,1:k]) + exp.psum[i,j,k] <- exp(psum[i,j,k]) + prob[i,j,k] <- exp.psum[i,j,k] / sum(exp.psum[i,j,1:K[j]]) + } } } + for (j in 1:p) { + delta[j,1] <- 0.0 + for (k in 2:K[j]) { + delta[j,k] ~ dnorm(m.delta, pr.delta) + } + } + beta ~ dnorm(m.beta,pr.beta) + pr.delta <- pow(s.delta, -2) + pr.beta <- pow(s.beta,-2) +} diff --git a/inst/jags/bpcm_dif b/inst/jags/bpcm_dif new file mode 100644 index 0000000..14a13c2 --- /dev/null +++ b/inst/jags/bpcm_dif @@ -0,0 +1,46 @@ +model { + for (i in 1:n) { + for (j in 1:p) { + Y[i,j] ~ dcat(prob[i,j,1:K[j]]) + } + theta[i] ~ dnorm(0, 1) + } + for (i in 1:n) { + for (j in 1:p) { + for (k in 1:K[j] ) { + eta[i,j,k] <- theta[i] + (beta * Z[i,1]) - ( delta[j,k] + Z[i,1]*ifelse(k==1,0,gamma[j,k]) ) + psum[i,j,k] <- sum(eta[i,j,1:k]) + exp.psum[i,j,k] <- exp(psum[i,j,k]) + prob[i,j,k] <- exp.psum[i,j,k] / sum(exp.psum[i,j,1:K[j]]) + } } } + for (j in 1:pnodif) { + delta[j,1] <- 0.0 + gamma[j,1] <- 0.0 + for (k in 2:K[j]) { + delta[j,k] ~ dnorm(m.delta, pr.delta) + gamma[j,k] <- 0.0 + } + } + for (j in pnodif1:pnounif) { + delta[j,1] <- 0.0 + gamma[j,1] <- 0.0 + for (k in 2:K[j]) { + delta[j,k] ~ dnorm(m.delta, pr.delta) + gamma[j,k] ~ dnorm(m.gamma,s.gamma) + } + } + for (j in pnounif1:p) { + delta[j,1] <- 0.0 + gamma[j,1] <- 0.0 + gamma[j,2] ~ dnorm(m.gamma,pr.gamma) + delta[j,2] ~ dnorm(m.delta, pr.delta) + for (k in 3:K[j]) { + delta[j,k] ~ dnorm(m.delta, pr.delta) + gamma[j,k] <- gamma[j,2] + } + } + beta ~ dnorm(m.beta,pr.beta) + pr.delta <- pow(s.delta, -2) + pr.beta <- pow(s.beta,-2) + pr.gamma <- pow(s.gamma,-2) +}