Chapter 12 Example 2.7

We want to apply the Region of Practical Equivalence approach to the runstitching differences. Specifically, we will apply this to the model fit.0.5. This assumes a prior of mean 0 and standard deviation .5 for the mean with a prior expectation of .75 and reasonably dispersed for the standard deviation.

Get the data:
> library(cfcdae)
> library(bcfcdae)
> data(RunStitch)
> differences <- RunStitch[,"Standard"] - RunStitch[,"Ergonomic"]
Fit the basic model from Example 2.6
> fit.0.5 <- bglmm(differences ~ 1, sigma0.mean=.75, sigma0.shape=3, 
+                  sigma.mean.int=.5, sigma.shape.int=20000, 
+                  quiet=TRUE,seed=10007)

Suppose that you want to advertise the ergonomic setup as being different from the standard. In that case, you might elect a very narrow region of equivalence, say (-.01,.01). Such a narrow ROPE could make sense for a customer who makes a lot of shirts. On the other hand, if you have to reequip your factory and retrain your employees, you might think that the two setups are equivalent unless these is a much larger change. Then you might have a ROPE of (-.2,.2).

To estimate the posterior probability of the ROPE, we need to extract the MCMC samples and see what fraction of them fall into the ROPE. We can extract all of the MCMC samples for the intercept via
> mcmc.samples <- fixef(fit.0.5,summary=FALSE)
There will be a column for each parameter, but here the parameter of interest is the intercept.
> mcmc.intercept <- mcmc.samples[,"(Intercept)"]
This is a vector of length 4,000, and we need to find the fractions of samples that fall into our two potential ROPEs.
> length(mcmc.intercept)
[1] 4000
> sum(mcmc.intercept >  -.01 & mcmc.intercept < .01)/length(mcmc.intercept)
[1] 0.02775
> sum(mcmc.intercept >  -.2 & mcmc.intercept < .2)/length(mcmc.intercept)
[1] 0.606

The probability in the narrow ROPE is estimated at .0225, so with this very stringent definition of equivalence, we would probably conclude against equivalence. On the other hand, the probability of the wider ROPE is about .6, so we would probably not decide in favor of either equivalence or non-equivalence.