I have a dataset for which I need to study both fixed and random effects. And I would love doing it with a LMM or a GLMM, but my data is not normal at all and I cannot transform by any means.

So... my only option left is Bayesian statistics (which I know they rock but I don't really understand)

My data looks like this

  d.Treatment d.Meters d.Plot d.percentageLoss
1      Control    -5.00     18         77.14286
2      Control     9.55     18         76.30058
3      Control     4.70     18         76.04790
4      Control    -5.00     27         81.06509
5      Control    20.00     27         80.58824
6      Control    -0.10     18         78.23529
7      Control    20.00     18         77.77778
8      Control   100.00     31         77.57576
9    ControlNM    20.00     37         84.07643
10   ControlNM     0.00     37         84.66258
11   ControlNM    50.00     38         85.23490
12   ControlNM    50.00     37         83.13953
13   ControlNM    20.00     40         83.13953
14   ControlNM    -5.00     39         84.21053
15   ControlNM     0.00     39         90.41916
16   ControlNM    20.00     39         84.93976
17   ControlNM     0.00     40         85.54217
18   ControlNM     5.00     38         82.69231
19   ControlNM    -5.00     38         85.29412
20   ControlNM     1.00     38         82.68156
21   ControlNM     1.00     40         81.33333
22   ControlNM     1.00     37         83.88889
23   ControlNM     5.00     37         82.60870
24   ControlNM    20.00     38         85.71429
25   ControlNM     1.00     39         80.66667
26   ControlNM    -5.00     37         83.33333
27          Mu     3.55     25         76.81159
28          Mu     5.90     20         75.14793
29          Mu    -0.10     25         77.71429
30          Mu     9.55     29         76.51007
31          Mu     4.70     25         82.60870
32          Mu     2.90     20         76.92308
33          Mu     3.55     24         77.30496
34          Mu     3.55     20         80.35714
35          Mu     0.55     24         68.15287
36          Mu     8.90     29         78.53107
37          Mu     8.90     25         77.05882
38          Mu     1.70     24         74.85714
39          Mu    -5.00     29         82.87293
40          Mu     1.70     29         77.10843...

So I want to know how much percentageLoss is affected by treatment + distance (fixed effects) with plot as a random effect.

For a GLMM it looks as follows

percentageLoss ~ Meters + Treatment + (1 | Plot)

I would run this on R and of course I'm not asking you to do it for me, but I have read two extensive tutorials so far and it feels like they're in chinese.

Could you please explain me how the bayesian MCMC works in this case, which priors I need to assume and why (how do I pick my priors if I have no information on similar data) and also what would the posterior mean.

I know it's a lot to ask for but there is not a single "guide for dummies" (statistiscally challenged students) on the topic... Any literature or similar R codes would already be of great help.

  • $\begingroup$ Check out the brms package in R. It lets you use the GLMM syntax to fit Bayesian mixed effects models. $\endgroup$ – Will Jun 30 '17 at 1:07

I'll try to brief and get some intuition across. In answer to your three questions:

  1. How does MCMC work? There are lots of packages which implement MCMC algorithms, which for the most part means you don't really need a deep understanding of their inner workings. The basic principle is: I try a bunch (thousands to hundreds of thousands) of different parameter values and pick the 'best' ones -- I pick a possible combination of parameters somehow (simplest case is randomly, as in Metropolis-Hastings). Then I compute the value of the posterior at those parameter values and compare it to the value of the posterior at the previous step. If the new value is higher than the previous, i.e. Ive found a better set of parameters, I accept the new values. If not, I reject it with a certain probability. So basically I iteratively find better and better parameters, but I allow some bad ones to sneak in too because I can't be completely sure that my data are any good. The rate at which I allow bad parameters to sneak in depends on how bad they are.

  2. Choice of priors. There are a few basic principles for choosing priors when you don't have substantive prior information. Andrew Gelman is probably the main guy for this, he has a great blog and a bunch of really good publications on what he calls 'weakly informative priors'. Basically something broad (large standard deviation) with a zero mean is generally OK. The idea is that, while I don't necessarily have any information about the precise distribution of my parameter a priori, I do know that it shouldn't be excessively large, so I can just use a prior that doesn't preclude any reasonable values but does preclude values which are unreasonably large, say in the thousands for a standardised predictor. So in short, broad zero-mean normal priors are generally good as a default. Bear in mind that 'broad' is relative to the scale of your predictors.

  3. What does the posterior mean? The posterior is just a distribution over plausible parameter values, just like you might have a normal distribution over the heights of females or whatever. Conceptually the posterior is slightly different in that it reflects a 'degree of belief' rather than the frequency at which we can expect to observe a given outcome -- the y axis is like a measure of credibility; given the data and my prior beliefs, how confident am I that the parameter lies in a given range -- this is what the posterior tells you.

EDIT: I realise I didn't really address the mixed-effects part of the question. Essentially in a Bayesian framework, you assume that the mixed effects are the sum of two components, $\theta_{ij} = \tilde{\theta_i} + \eta_{ij}$. The priors for the fixed effects $\tilde{\theta_i}$ are chosen according to the problem, either informative or weakly informative will do. The priors for the random effects $\eta_{ij}$ are the distribution you're assuming for the mixed effects, which will usually be zero-mean Gaussian, but can be anything. Because the variance of the random effect is unknown, we also need a prior on that -- a uniform prior might be OK, but Gelman advocates for half-cauchy priors, and that's what I'd use.

This tutorial is pretty practically-based and includes code etc.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.