Socks example with die score

What is the effect of the die on success in throwing the sock? Modify the code to include the die score as well as distance and hand. You should centre the die scores as we did distance by subtracting the mean.

Model answer

We need to add the die score and its coefficient and the prior for the coefficient to the JAGS model:

# Save this in the file "sox_DHd.jags"
model{
   # Likelihood
   for(i in 1:N) {
     logit(p[i]) <- b0 + bDist * distance[i] + 
         bHand * hand[i] + bDie * die[i]
     y[i] ~ dbern(p[i])
   }
   # Priors
   p0 ~ dbeta(1, 1)
   b0 <- logit(p0)
   bDist ~ dunif(-5, 5)
   bHand ~ dunif(-5, 5)
   bDie ~ dunif(-5, 5)
 }

We follow on from the code in the main page. We centre the data and prepare the list for JAGS:

# Centre the die scores
dieC <- sox$die - mean(sox$die)

jagsData <- list(N=N, y = sox$result, distance = distanceC,
    hand = hand, die = dieC)
str(jagsData)
# List of 5
# $ N       : int 840
# $ y       : int [1:840] 1 1 0 0 0 0 0 0 1 1 ...
# $ distance: num [1:840] -1.5 -1.5 -0.5 -0.5 0.5 0.5 1.5 1.5 ...
# $ hand    : num [1:840] 1 0 1 0 1 0 1 0 1 0 ...
# $ die     : num [1:840] -3.414 0.586 -3.414 0.586 1.586 ...

wanted <- c("p0", "bDist", "bHand", "bDie")

Now we are ready to run the model and look at the results:

# Takes 20 secs
( soxDHd <- jags(jagsData, NULL, wanted, "sox_DHd.jags",
  n.chains=3, n.iter=10000, parallel=TRUE, DIC=FALSE) )
#         mean    sd   2.5%    50%  97.5%  Rhat n.eff
# p0     0.385 0.026  0.334  0.385  0.438     1 18844
# bDist -0.940 0.077 -1.091 -0.939 -0.793     1 19927
# bHand  0.403 0.157  0.096  0.402  0.713     1 30000
# bDie   0.025 0.028 -0.028  0.025  0.080     1 30000

diagPlot(soxDHd)

mco <- mcmcOutput(soxDHd)
plot(mco, "bDie", compVal=0)

The results for p0, bDist and bHand are almost exactly the same as before, so we’ll only show the posterior for bDie:

The first thing to note is that the posterior mean of bDie is very small: getting a die score which is higher by 1 unit is equivalent to stepping forward by just 3cm. Secondly, we are about 80% certain that it’s positive.

So what do you think? Does the die score effect the probability of getting the sock in the box? Remember that the die was rolled after the throw, and I can’t think of a biological mechanism for the die to have an effect.

We include the die roll in the activity to show the effect of a covariate which actually has no effect on the response, but which can show an effect when you run the analysis. Sometimes when we do the activity, the die effect appears to be definitely positive or negative, but that doesn’t mean that we should conclude that we have made a major scientific discovery!

The main lesson is, “Only include covariates in your model if they make biological sense!”

One thought on “Socks example with die score”

Leave a Reply

Your email address will not be published. Required fields are marked *

The maximum upload file size: 1 MB. You can upload: image, document, text, archive.