Hot answers tagged

8

It is actually possible to have missing values in the predictors, although you need to put a prior distribution on any missing values (or all of the values, if this is easier). So if you have some missing values in peace[i] then you need an additional line in the code as follows: peace[i] ~ dnorm(0, 10^-6) # Or whatever would be an appropriate prior ...


8

One draw from binomial distribution generally is enough. But it depends of the data you have. If you have the data of how many heads in the individual coin flips have been seen in total, then binomial distribution is enough, no need for detailed model with N bernoulli flips. However, if you have data on results of individual coin flips and you need to ...


8

It might be late... but, Please note 2 things: Adding data points is not advised as it would change degrees of freedom. Mean estimates of fixed effect could be well estimated, but all inference should be avoided with such models. It is hard to "let the data speaks" if you change it. Of course it only works with integer-valued weights (you cannot duplicate ...


8

I don't know for sure what the trick is, but this is my guess. Using JAGS syntax to specify $\xi \sim \mathcal D(\alpha)$, you would normally do something like this: xi ~ dirichlet(alpha[]) JAGS would then not allow you to assign a prior to $\alpha = (\alpha_1, \ldots, \alpha_J)$. Instead, let $\xi^\star_j \sim \mbox{Gamma}(\alpha_j, 1)$. Then it can be ...


7

Both models will give the exact same results. Why? The Likelihood principle. RJags is an R package that uses the software JAGS to conduct Bayesian inference, and any fully Bayesian procedure, one where inference proceeds from the posterior distribution, will satisfy the Likelihood principle. Essentially, the Likelihood principle states that if two ...


7

I was asked to re-post this answer here from my comment at http://doingbayesiandataanalysis.blogspot.com/2012/01/complete-example-of-right-censoring-in.html The specifics of this answer relate to the model in that comment, but the concepts apply to the topic here. The core of the JAGS model for censored data is this: isCensored[i] ~ dinterval( y[i] , ...


7

Following the suggestion from user777, it looks like the answer to my first question is "use Stan." After rewriting the model in Stan, here are the trajectories (4 chains x 5000 iterations after burn-in): And the autocorrelation plots: Much better! For completeness, here's the Stan code: data { // Data: Exogenously given ...


7

When using Markov chain Monte Carlo (MCMC) algorithms in Bayesian analysis, often the goal is to sample from the posterior distribution. We resort to MCMC when other independent sampling techniques are not possible (like rejection sampling). The problem however with MCMC is that the resulting samples are correlated. This is because each subsequent sample is ...


6

Other answered is correct about BUGS but his answer does not apply to JAGS (at least, not to rjags, R2jags might be different). I haven't used JAGS directly, but the writer of rjags is the creator of JAGS so I would guess they use the same convention. In rjags, the jags.model object keeps track of the number of iterations that the chain(s) have been run. ...


6

This an interesting problem termed 'species-sampling', that has received a lot of attention over the years, and encompasses many other estimation problems (such as mark-recapture). Suffice it to say, JAGS will not help you in this case--JAGS cannot handle Markov chains with a variable dimension across iterations. One must recourse to an MCMC scheme designed ...


5

(most likely) in one of these lines: denom <- 1 + sum(exp(p)) D1 <- (exp(p[1]))/denom D2 <- (exp(p[2]))/denom D3 <- (exp(p[3]))/denom When the code is running those, part of the contents of p isn't something you want to pass to an exp function. It says "invalid vector argument". Can exp take a vector argument? If it can, maybe that's a message ...


5

The logarithm of the binomial coefficient can be implemented in WinBUGS/JAGS by using the function logfact which is the logarithm of the factorial: $\ln(x!)$. Alternatively, we could use the function of the logarithm of the gamma function $\ln(\Gamma(x))$ loggam. The binomial coefficient is defined as: $$ \binom{n}{k} = \frac{n!}{k!(n-k)!} $$ Using the ...


5

I had exactly the same problem! After reading the R2jags::jags source I came to the same conclusion as @Hao Ye, I will expand on more detail: First, I noticed that this happens only for some models: if you remove the dgamma from your example and use dunif instead, i.e. if you modify your example this way: ... inits <- function(){list(beta0=rnorm(1), ...


5

Just to clarify your model Let $y_{ij}$ be reaction time for participant $i$ on trial $j$. $$y_{ij} \sim N(\mu_i, \sigma^2_i)$$ And then you model $\mu_i$ and $\sigma^2_i$ as coming from some other distribution with hyperparameters. You ask if I would want to compare two subjects what distributions should I compare? So for example, if you wanted to ...


5

In BUGS/JAGS, the order in which statements are written does not matter. In Stan statements execute in the order in which they are written (see Stan 2.2.0 Reference Manual, pg. 405). Thus your last statement is in the wrong place: lambda is sampled from a gamma distribution, but that happens after the previous statements, so it's sampled 'in the air'. ...


5

Your model estimate would be a useful prior. I have applied the following approach in LeBauer et al 2013, and have adapted code from priors_demo.Rmd below. To parameterize this prior using simulation, consider your model $$ \textrm{logLC}_{50} = b_0 X+b_1$$ Assume $b_0 \sim N(0.94, 0.03)$ and $b_1 \sim N(1.33, 0.1)$; $\textrm{Lkow}$ is known (a ...


5

The problem is not that sum(p[]) < 1 but that you are modelling an incompletely observed multinomial distribution. As you have found, one of the limitations of JAGS is that multivariate nodes cannot be incompletely observed. There is one possibility to rescue your model that depends on a particular choice of prior distribution for N[i]. If this has a ...


5

Generally, in Bayesian model you do predictions on new data the same way as you do with non-Bayesian models. As your example is complicated I will provide a simplified one to make things easier to illustrate. Say you want to estimate linear regression model $$ y_i = \beta_0 + \beta_1 x_i + \varepsilon_i $$ and based on the model you want to predict $y_\...


5

You want a distribution for each quarter (given a state), each state (given a region), and each region. That means you'll need at least some state parameters indexed by s (in your model b0, b1, b2), region parameters indexed by r (which I'll call c0, c1, c2), and global parameters (which I'll call d0, d1, d2). Your model should look something like this:...


4

I don't know if this is a solution for you, but since the lme4 glmer function can provide random intercept posterior median estimates and their conditional variance - and under the assumption of normality (for random effect), posterior median = posterior mode - wouldn't it be valid to do a parametric bootstrap repeatedly drawing from the estimated posterior ...


4

C <- 10000 #Constant 1/0 trick # Likelihood: for ( i in 1:ny ) { #Likelihood of the count model component LikCountModel[i] <- pow(mu[i],y[i])/y_fact[i]*exp(-mu[i]) #Count model component eta[i] <- bet0 + inprod( beta[] , B[i,] ) mu[i] <- exp(eta[i]) #ZI Component zeta[i] <- gamm0 + inprod( gamma[] , G[i,] ) w[i] <- exp(zeta[i])/(1+exp(...


4

First, it's worth pointing out thatglm does not perform bayesian regression. The 'weights' parameter is basically a short hand for "proportion of observations," which can be replaced with up-sampling your dataset appropriately. For example: x=1:10 y=jitter(10*x) w=sample(x,10) augmented.x=NULL augmented.y=NULL for(i in 1:length(x)){ augmented.x=c(...


4

Reading through the comments on the othe answers, I believe the correct answer to the question that was intended to be asked is "they don't", in general. As has been mentioned, they construct a DAG and look at the Markov blanket and then (roughly) do the following. If the Markov blanket around a node correspond to a full conditional that is in a lookup ...


4

I figured out the answer, with help from Martyn Plummer. My code uses the inverse link for the gamma model (and no inverse of the predictors). Also, this code requires the 'glm' module for JAGS. model{ # For the ones trick C <- 10000 # for every observation for(i in 1:N){ # define the logistic regression model, where w is the probability ...


4

In the answers there (if I understood correctly) I learned that within-subject variance does not effect inferences made about group means and it is ok to simply take the averages of averages to calculate group mean, then calculate within-group variance and use that to perform significance tests. Let me develop this idea here. The model for the ...


4

Here is one example of implementing a basic macro substitution system for JAGS scripts. Explanation of the system Define a function that takes as arguments any optional elements of the script. For any aspects of the script that vary across argument values, record a macro token. This should be some unique text. Starting and ending with some symbols may ...


4

Regardless of whether the given covariance matrix correctly models the covariance matrix of an AR(1) process (or an MA(1) process) or not, the sum of all the entries in a covariance matrix is the variance of the sum of the $n$ random variables. Since this variance must be nonnegative, we get that in order for your matrix to be a valid covariance matrix, it ...


4

JAGS model notation is almost exactly the same as would you describe this model mathematically: $$ \alpha \sim \mathrm{Normal}(0, .001) \\ \beta \sim \mathrm{Normal}(0, .001) \\ \sigma_y \sim \mathrm{Uniform}(0, 100) \\ \tau_y = 1 / \sigma_y^2 \\ \tau_x \sim \mathrm{Uniform}(.03, .05) \\ x_{0i} \sim \mathrm{Normal}(0, .04) \\ x_i \sim \mathrm{Normal}(x_{0i}...


4

The use of vague or informative prior depends on the amount of knowledge that you have for the parameters that you want to assign the prior. I consider the following cases: No experts information and big/small data set, in that case, a vague prior would do the job Experts information and big/small data set, in that case, a more informative prior based on ...


3

I suggest that instead of treating the $m_i$ as missing data, you integrate them out, since this marginalization leaves the form of the data distribution unchanged (in the sense that it's still a multivariate normal distribution). Marginally of the $m_i$, $\mathrm{E}(y_i)=\mu$ and $\mathrm{Cov}\left(y_{i},y_{j}\right)=\begin{cases}2\sigma_{m}^{2}+\sigma_{y}^...


Only top voted, non community-wiki answers of a minimum length are eligible