This is the code to make the graph of predicted values at a single time point. You’ll have to modify it, because you have time as an additional predictor, but the logic should be the same.

Get data

library(brms)
## Loading required package: Rcpp
## Loading 'brms' package (version 2.13.5). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').
## 
## Attaching package: 'brms'
## The following object is masked from 'package:stats':
## 
##     ar
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.2     ✓ purrr   0.3.4
## ✓ tibble  3.0.4     ✓ dplyr   1.0.2
## ✓ tidyr   1.1.2     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(tidybayes)
## 
## Attaching package: 'tidybayes'
## The following objects are masked from 'package:brms':
## 
##     dstudent_t, pstudent_t, qstudent_t, rstudent_t
DOD_42mo <- read.csv("Longitudinal/DOD_42mo.csv")
DOD_42mo <-DOD_42mo  %>% 
 mutate(
    Prime = ifelse(Prime.Type==.5, yes="DOD", no="POD"), 
    Overlap = ifelse(Verb.match==.5, yes="Match", no="No Match")
    )

Plot condition means and SDs using raw data. One of the error bars is wrong here but I’m not sure why. But I think you get the logic.

(p1 <- DOD_42mo %>% 
  group_by(Prime, Overlap) %>%
  summarise(mean = mean(RESPONSE.CODE, na.rm=TRUE), 
            sd = sd(RESPONSE.CODE, na.rm=TRUE)) %>%
  
ggplot(aes(y=mean, x=Prime, fill=Overlap)) +
  geom_bar(stat="identity",position="dodge") +
 geom_errorbar(aes(ymin=mean, ymax=mean+sd), width=.2,
                 position=position_dodge(.9)) +
  ylim(0, .6)
)
## `summarise()` regrouping output by 'Prime' (override with `.groups` argument)

logit_prior = prior(normal(0,2), class=sd)

Now, fit the model.

m2 <- brm(RESPONSE.CODE ~ 1 + Prime.Type*Verb.match +   (1 +  Prime.Type*Verb.match  || Part..ID) + (1 +  Prime.Type*Verb.match  || Target.verb), logit_prior,family="bernoulli", cores=4, data=DOD_42mo)
## Warning: Rows containing NAs were excluded from the model.
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## Start sampling

Save the posterior samples from this model as an object.

post <- posterior_samples(m2)

If we wanted to look at the posterior distribution (modal point + highest posterior density interval) for specific parameters we could do so as follows. I think this code should be pretty straightforward.

post %>%
select(b_Intercept, b_Prime.Type, b_Verb.match, `b_Prime.Type:Verb.match`) %>% 
gather(parameter) %>% # make long
group_by(parameter) %>% # group by parameter
mean_hdi(value) # highest posterior desnity for values within each parameter. 

If you want to get condition means, you just combined those parameters to calculated the predicted values for each condition as well. This is how you would get means and hdis for each condition, convert them back to probabilities and plot them.

# Note that in transmute() I've used inv_logit_scaled(). This function transforms from logits to probabilities. To plot this stuff on the logit scale, just get rid of that function. 
p2 <- post %>%
select(b_Intercept, b_Prime.Type, b_Verb.match, `b_Prime.Type:Verb.match`) %>%
transmute(
  DOD_Match = inv_logit_scaled(b_Intercept +  .5*b_Prime.Type + .5*b_Verb.match + .25*`b_Prime.Type:Verb.match`), # use the regression equation for each condition to predict its value. 
  DOD_NoMatch = inv_logit_scaled(b_Intercept + .5*b_Prime.Type - .5*b_Verb.match - .25*`b_Prime.Type:Verb.match`), 
  POD_Match = inv_logit_scaled(b_Intercept -  .5*b_Prime.Type + .5*b_Verb.match - .25*`b_Prime.Type:Verb.match`), 
  POD_NoMatch = inv_logit_scaled(b_Intercept -  .5*b_Prime.Type - .5*b_Verb.match + .25*`b_Prime.Type:Verb.match`)
) %>%
gather(parameter) %>%
group_by(parameter) %>%
mean_hdi(value)  %>%
  mutate(
    Prime = ifelse(parameter=="DOD_Match" | parameter=="DOD_NoMatch", "DOD", "POD"),
    Overlap= ifelse(parameter=="DOD_Match" | parameter=="POD_Match", "Match", "UnMatch")
  ) %>%
  ggplot(aes(y=value, x=Prime, fill=Overlap)) +
  geom_bar(stat="identity",position="dodge") +
 geom_errorbar(aes(ymin=.lower, ymax=.upper), width=.2,
                 position=position_dodge(.9)) +
  ylim(0, .65)

library(patchwork)
(p1 | p2) + ylim(0, .6)
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.