turtles {bridgesampling}

Turtles Data from Janzen, Tucker, and Paukstis (2000)

Description

This data set contains information about 244 newborn turtles from 31 different clutches. For each turtle, the data set includes information about survival status (column y; 0 = died, 1 = survived), birth weight in grams (column x), and clutch (family) membership (column clutch; an integer between one and 31). The clutches have been ordered according to mean birth weight.

Usage

turtles

Format

A data.frame with 244 rows and 3 variables.

Source

Janzen, F. J., Tucker, J. K., & Paukstis, G. L. (2000). Experimental analysis of an early life-history stage: Selection on size of hatchling turtles. Ecology, 81(8), 2290-2304. doi: 10.2307/177115

Overstall, A. M., & Forster, J. J. (2010). Default Bayesian model determination methods for generalised linear mixed models. Computational Statistics & Data Analysis, 54, 3269-3288. doi: 10.1016/j.csda.2010.03.008

Sinharay, S., & Stern, H. S. (2005). An empirical comparison of methods for computing Bayes factors in generalized linear mixed models. Journal of Computational and Graphical Statistics, 14(2), 415-435. doi: 10.1198/106186005X47471

Examples

## Not run: 

################################################################################
# BAYESIAN GENERALIZED LINEAR MIXED MODEL (PROBIT REGRESSION)
################################################################################

library(bridgesampling)
library(rstan)

data("turtles")

#-------------------------------------------------------------------------------
# plot data
#-------------------------------------------------------------------------------

# reproduce Figure 1 from Sinharay & Stern (2005)
xticks <- pretty(turtles$clutch)
yticks <- pretty(turtles$x)

plot(1, type = "n", axes = FALSE, ylab = "", xlab = "", xlim = range(xticks),
     ylim =  range(yticks))
points(turtles$clutch, turtles$x, pch = ifelse(turtles$y, 21, 4), cex = 1.3,
       col = ifelse(turtles$y, "black", "darkred"), bg = "grey", lwd = 1.3)
axis(1, cex.axis = 1.4)
mtext("Clutch Identifier", side = 1, line = 2.9, cex = 1.8)
axis(2, las = 1, cex.axis = 1.4)
mtext("Birth Weight (Grams)", side = 2, line = 2.6, cex = 1.8)

#-------------------------------------------------------------------------------
# Analysis: Natural Selection Study (compute same BF as Sinharay & Stern, 2005)
#-------------------------------------------------------------------------------

### H0 (model without random intercepts) ###
H0_code <-
"data {
  int<lower = 1> N;
  int<lower = 0, upper = 1> y[N];
  real<lower = 0> x[N];
}
parameters {
  real alpha0_raw;
  real alpha1_raw;
}
transformed parameters {
  real alpha0 = sqrt(10.0) * alpha0_raw;
  real alpha1 = sqrt(10.0) * alpha1_raw;
}
model {
  // priors
  target += normal_lpdf(alpha0_raw | 0, 1);
  target += normal_lpdf(alpha1_raw | 0, 1);

  // likelihood
  for (i in 1:N) {
    target += bernoulli_lpmf(y[i] | Phi(alpha0 + alpha1 * x[i]));
  }
}"

### H1 (model with random intercepts) ###
H1_code <-
"data {
  int<lower = 1> N;
  int<lower = 0, upper = 1> y[N];
  real<lower = 0> x[N];
  int<lower = 1> C;
  int<lower = 1, upper = C> clutch[N];
}
parameters {
  real alpha0_raw;
  real alpha1_raw;
  vector[C] b_raw;
  real<lower = 0> sigma2;
}
transformed parameters {
  vector[C] b;
  real<lower = 0> sigma = sqrt(sigma2);
  real alpha0 = sqrt(10.0) * alpha0_raw;
  real alpha1 = sqrt(10.0) * alpha1_raw;
  b = sigma * b_raw;
}
model {
  // priors
  target += - 2 * log(1 + sigma2); // p(sigma2) = 1 / (1 + sigma2) ^ 2
  target += normal_lpdf(alpha0_raw | 0, 1);
  target += normal_lpdf(alpha1_raw | 0, 1);

  // random effects
  target += normal_lpdf(b_raw | 0, 1);

  // likelihood
  for (i in 1:N) {
    target += bernoulli_lpmf(y[i] | Phi(alpha0 + alpha1 * x[i] + b[clutch[i]]));
  }
}"

set.seed(1)
### fit models ###
stanfit_H0 <- stan(model_code = H0_code,
                   data = list(y = turtles$y, x = turtles$x, N = nrow(turtles)),
                   iter = 15500, warmup = 500, chains = 4, seed = 1)
stanfit_H1 <- stan(model_code = H1_code,
                   data = list(y = turtles$y, x = turtles$x, N = nrow(turtles),
                               C = max(turtles$clutch), clutch = turtles$clutch),
                   iter = 15500, warmup = 500, chains = 4, seed = 1)

set.seed(1)
### compute (log) marginal likelihoods ###
bridge_H0 <- bridge_sampler(stanfit_H0)
bridge_H1 <- bridge_sampler(stanfit_H1)

### compute approximate percentage errors ###
error_measures(bridge_H0)$percentage
error_measures(bridge_H1)$percentage

### summary ###
summary(bridge_H0)
summary(bridge_H1)

### compute Bayes factor ("true" value: BF01 = 1.273) ###
bf(bridge_H0, bridge_H1)


## End(Not run)