Bayesian Population Analysis using WinBUGSのコードをStanに移植してみる（第7章その1）：Taglibro de H：SSブログ

	ブログをはじめるログイン

Bayesian Population ..｜Bayesian Population .. ブログトップ

Bayesian Population Analysis using WinBUGSのコードをStanに移植してみる（第7章その1）　[統計] [編集]

第6章の次は第7章 Cormack-Jolly-Seber (CJS) モデルによる捕獲再捕獲データのモデリング。これは、Stan Modeling Language User's Guide and Reference Manualの第11章11.3節に実装例があるので、それを参考にしました。

7.3節 Models with Constant Parameters

cjs_c_c_code <- "
functions {
  int first_capture(int[] y_i) {
    for (k in 1:size(y_i))
      if (y_i[k])
        return k;
    return 0;
  }

  int last_capture(int[] y_i) {
    for (k_rev in 0:(size(y_i) - 1)) {
      int k;
      k <- size(y_i) - k_rev;
      if (y_i[k])
        return k;
    }
    return 0;
  }

  matrix prob_uncaptured(int nind, int n_occasions,
                         matrix p, matrix phi) {
    matrix[nind, n_occasions] chi;

    for (i in 1:nind) {
      chi[i, n_occasions] <- 1.0;
      for (t in 1:(n_occasions - 1)) {
        int t_curr;
        int t_next;

        t_curr <- n_occasions - t;
        t_next <- t_curr + 1;
        chi[i, t_curr] <- (1 - phi[i, t_curr]) +
                          phi[i, t_curr] *
                          (1 - p[i, t_next - 1]) *
                          chi[i, t_next];
      }
    }
    return chi;
  }
}

data {
  int<lower=0> nind;
  int<lower=2> n_occasions;
  int<lower=0,upper=1> y[nind, n_occasions];
}

transformed data {
  int<lower=0,upper=n_occasions> first[nind];
  int<lower=0,upper=n_occasions> last[nind];

  for (i in 1:nind)
    first[i] <- first_capture(y[i]);
  for (i in 1:nind)
    last[i] <- last_capture(y[i]);
}

parameters {
  real<lower=0,upper=1> mean_phi;
  real<lower=0,upper=1> mean_p;
}

transformed parameters {
  matrix<lower=0,upper=1>[nind, n_occasions - 1] phi;
  matrix<lower=0,upper=1>[nind, n_occasions - 1] p;
  matrix<lower=0,upper=1>[nind, n_occasions] chi;

  // Constraints
  for (i in 1:nind) {
    for (t in 1:(first[i] - 1)) {
      phi[i, t] <- 0;
      p[i, t] <- 0;
    }
    for (t in first[i]:(n_occasions - 1)) {
      phi[i, t] <- mean_phi;
      p[i, t] <- mean_p;
    }
  }

  chi <- prob_uncaptured(nind, n_occasions, p, phi);
}

model {
  // Priors
  mean_phi ~ uniform(0, 1);    // Prior for mean survival
  mean_p ~ uniform(0, 1);      // Prior for mean recapture

  // Likelihood
  for (i in 1:nind) {
    if (first[i] > 0) {
      for (t in (first[i] + 1):last[i]) {
        1 ~ bernoulli(phi[i, t - 1]);
        y[i, t] ~ bernoulli(p[i, t - 1]);
      }
      1 ~ bernoulli(chi[i, last[i]]);
    }
  }
}
"

## Bundle data
stan_data <- list(y = CH,
                  nind = dim(CH)[1],
                  n_occasions = dim(CH)[2])

## Initial values
inits <- function(){list(mean_phi = runif(1, 0, 1),
                         mean_p = runif(1, 0, 1))}

## Parameters monitored
params <- c("mean_phi", "mean_p")

## MCMC settings
ni <- 2000
nt <- 1
nb <- 1000
nc <- 3

## Stan
model <- stan_model(model_code = cjs_c_c_code)
cjs_c_c <- sampling(model,
                data = stan_data, init = inits, pars = params,
                chains = nc, thin = nt, iter = ni, warmup = nb,
                seed = 2,
                open_progress = FALSE)

7.4.1節 Fixed Time Effects

cjs_temp_fixeff_code <- "
functions {
  int first_capture(int[] y_i) {
    for (k in 1:size(y_i))
      if (y_i[k])
        return k;
    return 0;
  }

  int last_capture(int[] y_i) {
    for (k_rev in 0:(size(y_i) - 1)) {
      int k;
      k <- size(y_i) - k_rev;
      if (y_i[k])
        return k;
    }
    return 0;
  }

  matrix prob_uncaptured(int nind, int n_occasions,
                         matrix p, matrix phi) {
    matrix[nind, n_occasions] chi;

    for (i in 1:nind) {
      chi[i, n_occasions] <- 1.0;
      for (t in 1:(n_occasions - 1)) {
        int t_curr;
        int t_next;

        t_curr <- n_occasions - t;
        t_next <- t_curr + 1;
        chi[i, t_curr] <- (1 - phi[i, t_curr]) +
                          phi[i, t_curr] *
                          (1 - p[i, t_next - 1]) *
                          chi[i, t_next];
      }
    }
    return chi;
  }
}

data {
  int<lower=0> nind;
  int<lower=2> n_occasions;
  int<lower=0,upper=1> y[nind, n_occasions];
}

transformed data {
  int<lower=0,upper=n_occasions> first[nind];
  int<lower=0,upper=n_occasions> last[nind];

  for (i in 1:nind)
    first[i] <- first_capture(y[i]);
  for (i in 1:nind)
    last[i] <- last_capture(y[i]);
}

parameters {
  real<lower=0,upper=1> alpha[n_occasions - 1];
  real<lower=0,upper=1> beta[n_occasions - 1];
}

transformed parameters {
  matrix<lower=0,upper=1>[nind, n_occasions - 1] phi;
  matrix<lower=0,upper=1>[nind, n_occasions - 1] p;
  matrix<lower=0,upper=1>[nind, n_occasions] chi;

  // Constraints
  for (i in 1:nind) {
    for (t in 1:(first[i] - 1)) {
      phi[i, t] <- 0;
      p[i, t] <- 0;
    }
    for (t in first[i]:(n_occasions - 1)) {
      phi[i, t] <- alpha[t];
      p[i, t] <- beta[t];
    }
  }

  chi <- prob_uncaptured(nind, n_occasions, p, phi);
}

model {
  // Priors
  alpha ~ uniform(0, 1);    // Prior for mean survival
  beta ~ uniform(0, 1);     // Prior for mean recapture

  // Likelihood
  for (i in 1:nind) {
    if (first[i] > 0) {
      for (t in (first[i] + 1):last[i]) {
        1 ~ bernoulli(phi[i, t - 1]);
        y[i, t] ~ bernoulli(p[i, t - 1]);
      }
      1 ~ bernoulli(chi[i, last[i]]);
    }
  }
}
"

## Bundle data
stan_data <- list(y = CH,
                  nind = dim(CH)[1],
                  n_occasions = dim(CH)[2])

## Initial values
inits <- function(){list(mean_phi = runif(1, 0, 1),
                         mean_p = runif(1, 0, 1))}

## Parameters monitored
params <- c("alpha", "beta")

## MCMC settings
ni <- 4000
nt <- 3
nb <- 1000
nc <- 3

## Stan
model <- stan_model(model_code = cjs_temp_fixeff_code)
cjs_temp_fixeff <- sampling(model,
                            data = stan_data, init = inits, pars = params,
                            chains = nc, thin = nt, iter = ni, warmup = nb,
                            seed = 2,
                            open_progress = FALSE)

7.4.2節 Random Time Effects

cjs_temp_raneff_code <- "
functions {
  int first_capture(int[] y_i) {
    for (k in 1:size(y_i))
      if (y_i[k])
        return k;
    return 0;
  }

  int last_capture(int[] y_i) {
    for (k_rev in 0:(size(y_i) - 1)) {
      int k;
      k <- size(y_i) - k_rev;
      if (y_i[k])
        return k;
    }
    return 0;
  }

  matrix prob_uncaptured(int nind, int n_occasions,
                         matrix p, matrix phi) {
    matrix[nind, n_occasions] chi;

    for (i in 1:nind) {
      chi[i, n_occasions] <- 1.0;
      for (t in 1:(n_occasions - 1)) {
        int t_curr;
        int t_next;

        t_curr <- n_occasions - t;
        t_next <- t_curr + 1;
        chi[i, t_curr] <- (1 - phi[i, t_curr]) +
                          phi[i, t_curr] *
                          (1 - p[i, t_next - 1]) *
                          chi[i, t_next];
      }
    }
    return chi;
  }
}

data {
  int<lower=0> nind;
  int<lower=2> n_occasions;
  int<lower=0,upper=1> y[nind, n_occasions];
}

transformed data {
  int<lower=0,upper=n_occasions> first[nind];
  int<lower=0,upper=n_occasions> last[nind];

  for (i in 1:nind)
    first[i] <- first_capture(y[i]);
  for (i in 1:nind)
    last[i] <- last_capture(y[i]);
}

parameters {
  real<lower=0,upper=1> mean_phi;
  real<lower=0,upper=1> mean_p;
  real<lower=0> sigma;
  real epsilon[n_occasions - 1];
}

transformed parameters {
  matrix<lower=0,upper=1>[nind, n_occasions - 1] phi;
  matrix<lower=0,upper=1>[nind, n_occasions - 1] p;
  matrix<lower=0,upper=1>[nind, n_occasions] chi;
  real mu;

  mu <- logit(mean_phi);
  // Constraints
  for (i in 1:nind) {
    for (t in 1:(first[i] - 1)) {
      phi[i, t] <- 0;
      p[i, t] <- 0;
    }
    for (t in first[i]:(n_occasions - 1)) {
      phi[i, t] <- inv_logit(mu + epsilon[t]);
      p[i, t] <- mean_p;
    }
  }

  chi <- prob_uncaptured(nind, n_occasions, p, phi);
}

model {
  // Priors
  mean_phi ~ uniform(0, 1);    // Prior for mean survival
  mean_p ~ uniform(0, 1);      // Prior for mean recapture
  epsilon ~ normal(0, sigma);
  sigma ~ uniform(0, 10);

  for (i in 1:nind) {
    if (first[i] > 0) {
      for (t in (first[i] + 1):last[i]) {
        1 ~ bernoulli(phi[i, t - 1]);
        y[i, t] ~ bernoulli(p[i, t - 1]);
      }
      1 ~ bernoulli(chi[i, last[i]]);
    }
  }
}

generated quantities {
  real<lower=0> sigma2;

  sigma2 <- sigma^2;
}
"

## Bundle data
stan_data <- list(y = CH,
                  nind = dim(CH)[1],
                  n_occasions = dim(CH)[2])

## Initial values
inits <- function(){list(mean_phi = runif(1, 0, 1),
                         sigma = runif(1, 0, 10),
                         mean_p = runif(1, 0, 1))}  

## Parameters monitored
params <- c("mean_phi", "mean_p", "sigma2")

## MCMC settings
ni <- 4000
nt <- 3
nb <- 1000
nc <- 3

model <- stan_model(model_code = cjs_temp_raneff_code)
cjs_ran <- sampling(model,
                data = stan_data, init = inits, pars = params,
                chains = nc, thin = nt, iter = ni, warmup = nb,
                seed = 2,
                open_progress = FALSE)

7.4.3節 Temporal Covariates

cjs_cov_raneff_code <- "
functions {
  int first_capture(int[] y_i) {
    for (k in 1:size(y_i))
      if (y_i[k])
        return k;
    return 0;
  }

  int last_capture(int[] y_i) {
    for (k_rev in 0:(size(y_i) - 1)) {
      int k;
      k <- size(y_i) - k_rev;
      if (y_i[k])
        return k;
    }
    return 0;
  }

  matrix prob_uncaptured(int nind, int n_occasions,
                         matrix p, matrix phi) {
    matrix[nind, n_occasions] chi;

    for (i in 1:nind) {
      chi[i, n_occasions] <- 1.0;
      for (t in 1:(n_occasions - 1)) {
        int t_curr;
        int t_next;

        t_curr <- n_occasions - t;
        t_next <- t_curr + 1;
        chi[i, t_curr] <- (1 - phi[i, t_curr]) +
                          phi[i, t_curr] *
                          (1 - p[i, t_next - 1]) *
                          chi[i, t_next];
      }
    }
    return chi;
  }
}

data {
  int<lower=0> nind;
  int<lower=2> n_occasions;
  int<lower=0,upper=1> y[nind, n_occasions];
  real x[n_occasions - 1];
}

transformed data {
  int<lower=0,upper=n_occasions> first[nind];
  int<lower=0,upper=n_occasions> last[nind];
  vector<lower=0,upper=nind>[n_occasions] n_captured;

  for (i in 1:nind)
    first[i] <- first_capture(y[i]);
  for (i in 1:nind)
    last[i] <- last_capture(y[i]);
  n_captured <- rep_vector(0, n_occasions);
  for (t in 1:n_occasions)
    for (i in 1:nind)
      if (y[i, t])
        n_captured[t] <- n_captured[t] + 1;
}

parameters {
  real beta;
  real<lower=0,upper=1> mean_phi;
  real<lower=0,upper=1> mean_p;
  real epsilon[n_occasions - 1];
  real<lower=0> sigma;
}

transformed parameters {
  matrix<lower=0,upper=1>[nind, n_occasions - 1] phi;
  matrix<lower=0,upper=1>[nind, n_occasions - 1] p;
  matrix<lower=0,upper=1>[nind, n_occasions] chi;
  real mu;

  mu <- logit(mean_phi);
  // Constraints
  for (i in 1:nind) {
    for (t in 1:(first[i] - 1)) {
      phi[i, t] <- 0;
      p[i, t] <- 0;
    }
    for (t in first[i]:(n_occasions - 1)) {
      phi[i, t] <- inv_logit(mu + beta * x[t] + epsilon[t]);
      p[i, t] <- mean_p;
    }
  }

  chi <- prob_uncaptured(nind, n_occasions, p, phi);
}

model {
  // Priors
  mean_phi ~ uniform(0, 1);     // Prior for mean survival
  mean_p ~ uniform(0, 1);       // Prior for mean recapture
  beta ~ normal(0, 100);        // Prior for slope parameter
  sigma ~ uniform(0, 10);
  epsilon ~ normal(0, sigma);

  for (i in 1:nind) {
    if (first[i] > 0) {
      for (t in (first[i] + 1):last[i]) {
        1 ~ bernoulli(phi[i, t - 1]);
        y[i, t] ~ bernoulli(p[i, t - 1]);
      }
      1 ~ bernoulli(chi[i, last[i]]);
    }
  }
}

generated quantities {
  real<lower=0> sigma2;
  real<lower=0,upper=1> phi_est[n_occasions - 1];

  sigma2 <- sigma^2;
  for (t in 1:(n_occasions - 1))
    phi_est[t] <- inv_logit(mu + beta * x[t] + epsilon[t]); // Yearly survival
}
"

## Bundle data
stan_data <- list(y = CH,
                  nind = dim(CH)[1],
                  n_occasions = dim(CH)[2],
                  x = winter)

## Initial values
inits <- function(){list(mu = rnorm(1),
                         sigma = runif(1, 0, 5),
                         beta = runif(1, -5, 5),
                         mean_p = runif(1, 0, 1))}  

## Parameters monitored
params <- c("mean_phi", "mean_p", "phi_est", "sigma2", "beta")

# MCMC settings
ni <- 4000
nt <- 3
nb <- 1000
nc <- 3

## Stan
model <- stan_model(model_code = cjs_cov_raneff_code)
cjs_cov <- sampling(model,
                    data = stan_data, init = inits, pars = params,
                    chains = nc, thin = nt, iter = ni, warmup = nb,
                    seed = 2,
                    open_progress = FALSE)

最後のモデルは、かたむき(beta)がゆるく、ランダム効果が大きく推定される傾向にあるかもしれません。

タグ：STAn