The data is the observation from your werewolf detector that the mailman is a werewolf. That is information you went out and collected, presumably to test the hypothesis that the mailman is a werewolf.
The likelihood is the probability of the data given the model. In this case, it’s the probability that the data point (the werewolf detector observing that the mailman is a werewolf) would be observed, given that the model is true (the model being the hypothesis that the mailman is a werewolf). As we know the detector correctly identifies werewolves 99% of the time, it’s 0.99.
The prior is the information from your field guide that 1 in every 1000 people is a werewolf. This is knowledge that you have before you set out to test whether anyone is a werewolf.
The posterior is the final result: if we observe a positive werewolf test on an individual, given our prior knowledge that 1/1000 people are werewolves and the fact that our test is 99% accurate for both positive and negative results, the posterior probability that the mailman is a werewolf is 9%.
Incidentally, this seems like a very low probability but it can be demonstrated using Bayes’ Theorem. See the worked example on the Wikipedia page on Bayes’ Theorem. It uses drug testing as an example, which is not as fun as werewolf testing, but the math is the same.
Trace plot A resembles a “hairy caterpillar,” representing good model convergence. All Markov chains have arrived at the correct solution and are moving back and forth around it. The amount they go back and forth is proportional to the uncertainty we have about the value of the parameter.
In trace plot B, some of the Markov chains show trends in one direction or another represents poor model convergence. The chains are still moving around trying to find the correct solution. They are not mixing well, meaning they are not moving back and forth with the same stationary distribution.
In the first case, we can tell the researcher they can use the posterior distribution to make some inferences. They can do things like calculate the median and quantiles of the distribution, and make predictions.
In the second case, we would recommend the researcher to either run the MCMC for a greater number of iterations, or change the prior distributions to constrain the chains to explore parameter space more narrowly so they are more likely to converge on the correct solution.
Here is how you would proceed if you wanted to add two-way interaction terms to the model. This adds the interaction terms to the last model we fit that has random slopes for all main effects.
fit_withinteractions <- brm(yield ~ 1 + N + P + K + N:P + N:K + P:K + (1 + N + P + K | site) + (1 | block:site),
data = stvincent,
prior = c(prior(normal(0, 10), class = b)),
chains = 4,
iter = 2000,
warmup = 1000,
seed = 5,
file = 'fits/fit_2wayinteractions')Then, to examine the model diagnostics, look at the trace plots and the posterior predictive check plot (output not shown here). We should see “hairy caterpillar” traceplots and a posterior predictive check plot where the fitted values taken from the posterior have a similar distribution to the data.
plot(fit_withinteractions)
pp_check(fit_withinteractions)Next, look at the model summary. This will show you that the R-hat statistics are close to 1 to further confirm the model converged. That means we can make inferences based on the parameter estimates.
summary(fit_withinteractions)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: yield ~ 1 + N + P + K + N:P + N:K + P:K + (1 + N + P + K | site) + (1 | block:site)
## Data: stvincent (Number of observations: 324)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Multilevel Hyperparameters:
## ~block:site (Number of levels: 36)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 1.45 0.33 0.86 2.17 1.00 1555 2597
##
## ~site (Number of levels: 9)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 4.45 1.37 2.48 7.67 1.00 2394 2489
## sd(N1N) 0.86 0.70 0.03 2.61 1.00 2891 2330
## sd(N2N) 3.42 1.04 1.89 5.93 1.00 2612 2967
## sd(N3N) 3.54 1.33 1.45 6.66 1.00 2678 2117
## sd(P1P) 1.14 0.82 0.06 3.05 1.00 2192 2095
## sd(P2P) 1.69 0.91 0.18 3.74 1.00 1774 1347
## sd(P3P) 1.02 0.86 0.04 3.12 1.00 2354 1825
## sd(K1K) 1.45 0.94 0.09 3.54 1.00 1746 2132
## sd(K2K) 1.46 0.84 0.15 3.35 1.00 2165 1825
## sd(K3K) 2.06 1.21 0.15 4.82 1.00 1947 1759
## cor(Intercept,N1N) -0.00 0.30 -0.58 0.55 1.00 7187 2782
## cor(Intercept,N2N) 0.01 0.25 -0.46 0.49 1.00 3787 3215
## cor(N1N,N2N) 0.07 0.31 -0.54 0.63 1.00 1691 2330
## cor(Intercept,N3N) -0.02 0.25 -0.51 0.47 1.00 5388 3472
## cor(N1N,N3N) 0.09 0.31 -0.52 0.65 1.00 2319 2830
## cor(N2N,N3N) 0.39 0.24 -0.12 0.80 1.00 4073 3625
## cor(Intercept,P1P) -0.03 0.29 -0.58 0.53 1.00 7128 2962
## cor(N1N,P1P) -0.03 0.31 -0.61 0.56 1.00 4124 3131
## cor(N2N,P1P) 0.20 0.30 -0.43 0.72 1.00 4254 2946
## cor(N3N,P1P) 0.11 0.30 -0.48 0.65 1.00 5123 3311
## cor(Intercept,P2P) 0.17 0.27 -0.38 0.67 1.00 6294 3279
## cor(N1N,P2P) 0.05 0.31 -0.54 0.62 1.00 4254 3333
## cor(N2N,P2P) -0.07 0.27 -0.58 0.47 1.00 4757 3532
## cor(N3N,P2P) -0.01 0.27 -0.52 0.49 1.00 4061 3338
## cor(P1P,P2P) 0.06 0.30 -0.52 0.62 1.00 3366 3177
## cor(Intercept,P3P) 0.02 0.29 -0.53 0.55 1.00 8832 3217
## cor(N1N,P3P) -0.02 0.31 -0.60 0.58 1.00 5392 3037
## cor(N2N,P3P) 0.01 0.30 -0.59 0.58 1.00 6713 2976
## cor(N3N,P3P) 0.01 0.30 -0.57 0.56 1.00 5593 3342
## cor(P1P,P3P) 0.04 0.30 -0.53 0.59 1.00 3909 3325
## cor(P2P,P3P) 0.00 0.30 -0.58 0.58 1.00 3541 3398
## cor(Intercept,K1K) 0.13 0.28 -0.43 0.65 1.00 6422 3103
## cor(N1N,K1K) -0.07 0.31 -0.63 0.53 1.00 3937 3094
## cor(N2N,K1K) 0.15 0.29 -0.44 0.66 1.00 5357 3261
## cor(N3N,K1K) 0.09 0.29 -0.47 0.62 1.00 4345 3209
## cor(P1P,K1K) -0.01 0.30 -0.58 0.57 1.00 3522 3587
## cor(P2P,K1K) 0.17 0.30 -0.46 0.69 1.00 2952 3091
## cor(P3P,K1K) -0.04 0.30 -0.60 0.54 1.00 3256 3526
## cor(Intercept,K2K) -0.07 0.27 -0.58 0.48 1.00 6778 3106
## cor(N1N,K2K) -0.00 0.30 -0.59 0.56 1.00 4158 3220
## cor(N2N,K2K) -0.23 0.28 -0.71 0.35 1.00 5393 3217
## cor(N3N,K2K) -0.16 0.28 -0.65 0.39 1.00 4523 3424
## cor(P1P,K2K) -0.08 0.29 -0.62 0.50 1.00 3415 3259
## cor(P2P,K2K) -0.09 0.29 -0.61 0.49 1.00 3927 3490
## cor(P3P,K2K) 0.03 0.30 -0.55 0.59 1.00 2770 3204
## cor(K1K,K2K) -0.05 0.29 -0.59 0.52 1.00 3205 3525
## cor(Intercept,K3K) -0.20 0.28 -0.71 0.38 1.00 6213 2895
## cor(N1N,K3K) -0.04 0.30 -0.59 0.55 1.00 4396 3548
## cor(N2N,K3K) 0.16 0.28 -0.39 0.67 1.00 4790 3328
## cor(N3N,K3K) 0.07 0.28 -0.48 0.59 1.00 5066 3165
## cor(P1P,K3K) -0.02 0.30 -0.58 0.54 1.00 3567 3729
## cor(P2P,K3K) -0.12 0.29 -0.64 0.48 1.00 3350 3362
## cor(P3P,K3K) 0.01 0.30 -0.58 0.59 1.00 3120 3420
## cor(K1K,K3K) 0.01 0.30 -0.57 0.56 1.00 2997 3289
## cor(K2K,K3K) -0.04 0.29 -0.57 0.53 1.00 2917 3375
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 13.11 1.61 9.93 16.35 1.00 1776 2368
## N1N 0.67 8.15 -15.46 16.55 1.00 7928 2823
## N2N 4.40 1.45 1.50 7.22 1.00 3200 2975
## N3N 0.68 8.48 -16.27 16.96 1.00 8068 2893
## P1P 1.23 8.21 -14.68 16.67 1.00 7034 2985
## P2P -0.16 1.08 -2.41 1.91 1.00 4289 3059
## P3P 0.38 8.75 -16.74 17.63 1.00 7914 3105
## K1K 1.02 8.22 -15.05 17.06 1.00 6040 2486
## K2K -0.89 1.06 -2.93 1.20 1.00 4051 3186
## K3K 0.56 8.86 -17.09 17.91 1.00 8814 2895
## N1N:P1P 0.28 8.31 -16.03 16.64 1.00 6838 3004
## N2N:P1P 0.09 10.12 -19.64 19.82 1.00 8100 2409
## N3N:P1P 0.74 8.57 -15.95 17.50 1.00 6920 2888
## N1N:P2P -0.02 9.82 -19.24 19.43 1.00 12890 2874
## N2N:P2P -0.19 1.02 -2.18 1.84 1.00 7300 3043
## N3N:P2P 0.02 9.90 -19.24 19.68 1.00 10381 2906
## N1N:P3P 0.24 9.15 -17.62 17.81 1.00 8732 2769
## N2N:P3P 0.19 9.86 -18.82 19.46 1.00 10948 3157
## N3N:P3P 0.28 10.05 -19.30 19.99 1.00 10874 2579
## N1N:K1K -0.12 8.32 -16.44 16.15 1.00 7457 3031
## N2N:K1K 0.11 10.06 -19.63 19.55 1.00 12203 2812
## N3N:K1K 0.74 8.93 -16.53 18.16 1.00 9565 3130
## N1N:K2K -0.07 10.06 -19.52 19.24 1.00 9188 2461
## N2N:K2K -0.90 1.03 -2.90 1.08 1.00 7541 3006
## N3N:K2K -0.17 10.13 -19.91 19.52 1.00 12657 2745
## N1N:K3K 0.70 8.83 -16.94 18.00 1.00 7751 2861
## N2N:K3K 0.05 10.03 -19.50 19.22 1.00 11066 3168
## N3N:K3K -0.06 9.82 -18.66 19.41 1.00 10669 2321
## P1P:K1K 0.45 8.04 -14.69 16.22 1.00 7378 3301
## P2P:K1K 0.08 9.62 -18.90 18.88 1.00 9148 2672
## P3P:K1K 0.14 8.81 -17.01 17.47 1.00 8336 3133
## P1P:K2K -0.06 10.11 -19.57 19.85 1.00 11943 2690
## P2P:K2K 1.96 1.03 -0.08 3.97 1.00 5780 3261
## P3P:K2K -0.00 10.26 -20.42 19.62 1.00 12314 2810
## P1P:K3K 0.58 9.05 -16.60 18.28 1.00 8575 2194
## P2P:K3K -0.10 9.73 -18.95 18.47 1.00 8617 2748
## P3P:K3K 0.01 10.02 -19.05 19.59 1.00 9871 2875
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 3.04 0.14 2.78 3.33 1.00 4073 3139
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).We can see here that the point estimates of the interaction slopes are mostly close to zero with wide credible intervals, indicating weak to nonexistent interactions.
Finally, we can use LOO cross-validation to compare the model to others. In this case I will just compare the model with both interaction terms to the random-slope model with only the three main effects and no interaction terms that we fit during the main part of the lesson.
fit_withinteractions <- add_criterion(fit_withinteractions, 'loo')
loo_compare(fit_withinteractions, fit_randomNPKslopes)## elpd_diff se_diff
## fit_randomNPKslopes 0.0 0.0
## fit_withinteractions -0.9 2.3This shows us that the fit that includes the interaction terms does not perform any better in the LOO cross-validation than the fit that ignores interactions. This is further evidence that interactions between fertilization treatments are not very strong in this particular dataset.