[R-sig-ME] mixture DIC model

[Ben Bolker]

We'd like to help you, but your question is very hard to understand.

On Fri, Apr 14, 2017 at 3:40 AM, Mohammadreza Kordbagheri via
R-sig-mixed-models <r-sig-mixed-models at r-project.org> wrote:
> hello
> my name is alireza kordbagheri of shahid beheshti university tehran.
>
> i have two question of you:
>
> can you help me for my question:
> one question: i read your article tittle:
>
> i dont know, how draw heatmap ANOVA GAMMA.

   What does this mean?

A heatmap is an image plot drawn with colours drawn from a "heat
spectrum" (red to white)

ANOVA is analysis of variance (abused in R's anova() function as a
general-purpose tool for comparing nested models)

Gamma is a distribution (and a function)

>
> two question:i writing code model with JAGS softwar.
> i used of function GB2.

  What does that mean?
>
> I have problem, what DIC estimation?

  What is your question?

> This article is open access at
>
> https://projecteuclid.org/euclid.ba/1340370933.
>
> library(rjags)
>
> cat(file = "model.txt",
> "
> data {
>   C <- 100
>   for (i in 1:N) {
>     zeros[i] <- 0
>   }
> }
> model {
>   for (i in 1:N) {
>     zeros[i] ~ dpois(loglik[i] + C)
>     loglik[i] <- -log(pi[1] * l1[i] + pi[2] * l2[i])
>     l1[i] <- abs(a[1]) * pow(y[i], a[1] * p[1] - 1) / (pow(b1[i], a[1] * p[1]) * exp(loggam(p[1]) + loggam(q[1]) - loggam(p[1] + q[1])) * pow(1 + pow(y[i] / b1[i], a[1]), p[1] + q[1]))
>     l2[i] <- abs(a[2]) * pow(y[i], a[2] * p[2] - 1) / (pow(b2[i], a[2] * p[2]) * exp(loggam(p[2]) + loggam(q[2]) - loggam(p[2] + q[2])) * pow(1 + pow(y[i] / b2[i], a[2]), p[2] + q[2]))
>     b1[i]<- exp(mu[1] + alpha[1] * f[i] + beta1[1] * col1[i] + beta2[1] * col2[i] + beta3[1] * col3[i] + loggam(p[1]) + loggam(q[1]) - loggam(p[1] + 1 / a[1]) - loggam(q[1] - 1 / a[1]))
>     b2[i]<- exp(mu[2] + alpha[2] * f[i] + beta1[2] * col1[i] + beta2[2] * col2[i] + beta3[2] * col3[i] + loggam(p[2]) + loggam(q[2]) - loggam(p[2] + 1 / a[2]) - loggam(q[2] - 1 / a[2]))
>   }
>   pi[1] ~ dunif(0, 1)
>   pi[2] <- 1 - pi[1]
>   for(k in 1:2) {
>     mu[k] ~ dnorm(0, .01)
>     alpha[k] ~ dnorm(0, .01)
>     beta1[k] ~ dnorm(0, .01)
>     beta2[k] ~ dnorm(0, .01)
>     beta3[k] ~ dnorm(0, .01)
>     a[k] ~ dnorm(0, .01)
>     p[k] ~ dgamma(.01, .01)I(0, 90)
>     q[k] ~ dgamma(.01, .001)I(1, 3.5)
>   }
> }
> ")
>
> data <- matrix(c(
> 3323,  5, 1, 0, 0,
> 8332,  8, 0, 1, 1,
> 9572,  5, 1, 0, 0,
> 10172, 10, 1, 0, 1,
> 7631,  2, 1, 0, 0,
> 3855,  9, 0, 1, 1,
> 3252,  6, 1, 0, 1,
> 4433,  8, 0, 1, 1,
> 2188,  7, 1, 0, 0,
>   333,  4, 1, 0, 0,
>   199,  3, 0, 1, 1,
>   692,  9, 0, 1, 0,
>   311, 12, 1, 0, 0,
> 0.01,  5, 1, 0, 1,
>   405,  6, 0, 1, 0,
>   293,  6, 0, 1, 0,
>   76,  7, 1, 0, 1,
>   14,  9, 1, 0, 0,
> 3785, 21, 0, 1, 1,
> 10342, 11, 0, 1, 1), ncol = 5, byrow = TRUE)
>
> bugs.data <- list(N = nrow(data),
>                   y = data[, 1],
>                   f = data[, 2],
>                   col1 = data[, 3],
>                   col2 = data[, 4],
>                   col3 = data[, 5])
> set.seed(123)
> inits <- lapply(1:3, function(i) {
>     list(mu = rnorm(2, 0, 1),
>         alpha = rnorm(2, 0, 1),
>         beta1 = rnorm(2, 0, 1),
>         beta2 = rnorm(2, 0, 1),
>         beta3 = rnorm(2, 0, 1),
>         pi = c(runif(1, 0, 1), NA),
>         a = runif(2, -4, -2),
>         p = runif(2, 1, 2),
>         q = runif(2, 1, 2),
>         .RNG.name = "base::Mersenne-Twister",
>         .RNG.seed = i)
>     })
> pars <- c("mu", "alpha", "beta1", "beta2", "beta3", "pi",
>           "a", "p", "q")
> model <- jags.model("model.txt",
>                     data = bugs.data,
>                     inits = inits,
>                     n.chains = 3,
>                     n.adapt = 1000)
> update(model, 3000)
> samp <- coda.samples(model,
>                     variable.names = pars,
>                     n.iter = 4000, thin = 4)
> gelman.diag(samp, multivariate = FALSE)
> traceplot(samp[, 15]) # pi[1]
> traceplot(samp[, 16]) # pi[2]
>
> can you help me.
>
> Regards,
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models