There is a lot of stuff online about eliciting p-values for lme4 linear mixed-effects regression model coefficients (just google it). I opted for MCMCglmm because I like MCMC , and because lme4 has deprecated its own MCMC implementation.
[Here]'s the MCMCglmm quick start guide.
The main idea is simulating the regression model and seeing how much of the posterior density lies in regions beyond a certain value (typically 0). This is the Bayesian version of what p-values are supposed to be doing anyway.
The expression we use to form a model is a bit different from lme4.
model <- lmer(score ~ nstates_amp_and_freq_n + (nstates_amp_and_freq_n | id) + condition:phase, data=all, REML=F)
mc <- MCMCglmm(score ~ nstates_amp_and_freq_n + condition:phase, random=~us(nstates_amp_and_freq_n):id, data=all, nitt = 50000, thin=100, pr=TRUE)
where ```pr = TRUE``` indicates we want to keep samples of random effects as well. The samples can be accessed using ```mc$Sol``` and ```summary(mc)``` produces a summary that includes p-values.