library(lme4)
library(rjags)
library(R2jags)

## **************************************************
## Repeatability of a sexual signal trait
## **************************************************

## --------------------------------------------------
## Read and examine the data

## 1.
ff <- read.csv('https://www.sfu.ca/~lmgonigl/materials-qm/data/flycatcher.csv',
               stringsAsFactors=TRUE,
               strip.white=TRUE)

## --------------------------------------------------
## fit linear model using lme4

out.lme4 <- lmer(patch ~ 1 + (1|bird), data=ff)
summ.lme4 <- summary(out.lme4)
summ.lme4

## --------------------------------------------------
## fit linear model using JAGS

## 1. construct data
##
## data
my.data <- list(patch=ff$patch,
                bird=ff$bird,
                n=nrow(ff),
                nbird=length(unique(ff$bird)))

## JAGS model
my.model <- function() {

  ## PRIORS
  ##
  ## intercept
  b.intercept ~ dnorm(0,0.001)

  ## residuals
  sigma.resid ~ dunif(0,20)
  tau.resid <- 1/(sigma.resid*sigma.resid)
  
  for(i in 1:n) {
    mu[i] <- b.intercept
    patch[i] ~ dnorm(mu[i], tau.resid)
  }
}

## inits
my.inits <- function() {
  list()
}

## params to track
my.params <- c('b.intercept',
               'sigma.resid')

res <- jags(data=my.data,
            inits=my.inits,
            parameters.to.save=my.params,
            model.file=my.model,
            n.chains=3,
            n.iter=11000,
            n.burnin=1000,
            n.thin=100,
            working.directory=NULL)

summ.jags <- res$BUGSoutput$summary[,c('mean','sd','2.5%','97.5%')]

summ.jags
summ.lme4