So-net無料ブログ作成
検索選択

Bayesian Population Analysis using WinBUGSのコードをStanに移植してみる(第8章) [統計]

第7章がおわって、第8章 Estimation of Survival Using Mark-Recovery Data です。

8.2. The mark-recovery model as a state-space model

8.2.2. Analysis of a model with constant parameters

mr_ss_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;
  }

  real cell_prob(int n_occasions, row_vector s, row_vector r,
                   int first, int last) {
    vector[n_occasions] pr;

    pr[first] <- (1 - s[first]) * r[first];

    for (j in (first + 1):n_occasions - 1)
      pr[j] <- prod(segment(s, first, j - first)) * (1 - s[j]) * r[j];

    for (j in 1:(first - 1))
      pr[j] <- 0;

    for (t in 1:n_occasions - 1)
      pr[n_occasions] <- 1 - sum(head(pr, n_occasions - 1));

    return pr[last];
  }
}

data {
  int<lower=0> nind;
  int<lower=0> 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_s;
  real<lower=0,upper=1> mean_r;
}

transformed parameters {
  matrix<lower=0,upper=1>[nind, n_occasions - 1] s;
  matrix<lower=0,upper=1>[nind, n_occasions - 1] r;

  // Constraints
  for (i in 1:nind) {
    for (t in 1:(first[i] - 1)) {
      s[i, t] <- 0;
      r[i, t] <- 0;
    }
    for (t in first[i]:n_occasions - 1) {
      s[i, t] <- mean_s;
      r[i, t] <- mean_r;
    }
  }
}

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

  // Likelihood 
  for (i in 1:nind){
    real pr;

    if (first[i] > 0) {
      if (last[i] > first[i]) // Recovered
        pr <- cell_prob(n_occasions, s[i], r[i], first[i], last[i] - 1);
      else                    // Never Recovered
        pr <- cell_prob(n_occasions, s[i], r[i], first[i], n_occasions);
      1 ~ bernoulli(pr);
    }
  }
}
"

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

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

## Parameters monitored
params <- c("mean_s", "mean_r")

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

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

8.3. The mark-recovery model fitted with the multinomial likelihood

8.3.1. Constant parameters

mr_mnl_code <- "
data {
  int n_occasions;
  int marr[n_occasions, n_occasions + 1];
}

parameters {
  real<lower=0,upper=1> mean_s;
  real<lower=0,upper=1> mean_r;
}

transformed parameters {
  vector[n_occasions] s;
  vector[n_occasions] r;
  simplex[n_occasions + 1] pr[n_occasions];

  // Constraints
  for (t in 1:n_occasions) {
    s[t] <- mean_s;
    r[t] <- mean_r;
  }

  // Define the cell probabilities of the m-array
  // Main diagonal
  for (t in 1:n_occasions) {
    pr[t][t] <- (1 - s[t]) * r[t];

    // Above main diagonal
    for (j in (t + 1):n_occasions)
      pr[t][j] <- prod(segment(s, t, j - t)) * (1 - s[j]) * r[j];

    // Below main diagonal
    for (j in 1:(t - 1))
      pr[t][j] <- 0;
  } // t

  // Last column: probability of non-recovery
  for (t in 1:n_occasions)
    pr[t][n_occasions + 1] <- 1 - sum(head(pr[t], n_occasions));
}

model {
  // Priors
  mean_s ~ uniform(0, 1);           // Prior for mean survival
  mean_r ~ uniform(0, 1);           // Prior for mean recovery

  // Define the multinomial likelihood
  for (t in 1:n_occasions)
    marr[t] ~ multinomial(pr[t]);
}
"

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

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

## Parameters monitored
params <- c("mean_s", "mean_r")

## MCMC settings
ni <- 3000
nt <- 2
nb <- 1000
nc <- 4

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

8.3.2. Age-dependent parameters

mr_mnl_age_code <- "
data {
  int n_occasions;
  int marr_j[n_occasions, n_occasions + 1];
  int marr_a[n_occasions, n_occasions + 1];
}

parameters {
  real<lower=0,upper=1> mean_sj;
  real<lower=0,upper=1> mean_sa;
  real<lower=0,upper=1> mean_rj;
  real<lower=0,upper=1> mean_ra;
}

transformed parameters {
  vector[n_occasions] sj;
  vector[n_occasions] sa;
  vector[n_occasions] rj;
  vector[n_occasions] ra;
  simplex[n_occasions + 1] pr_a[n_occasions];
  simplex[n_occasions + 1] pr_j[n_occasions];

  // Constraints
  for (t in 1:n_occasions) {
    sj[t] <- mean_sj;
    sa[t] <- mean_sa;
    rj[t] <- mean_rj;
    ra[t] <- mean_ra;
  }

  // Define the cell probabilities of the juvenile m-array
  // Main diagonal
  for (t in 1:n_occasions) {
    pr_j[t, t] <- (1 - sj[t]) * rj[t];

    // Further above main diagonal
    for (j in (t + 2):n_occasions)
      pr_j[t, j] <- sj[t] * prod(segment(sa, t + 1, j - t - 1)) *
                    (1 - sa[j]) * ra[j];

    // Below main diagonal
    for (j in 1:(t - 1))
      pr_j[t, j] <- 0;
  } //t

  for (t in 1:(n_occasions - 1))
    // One above main diagonal
    pr_j[t, t + 1] <- sj[t] * (1 - sa[t + 1]) * ra[t + 1];

  // Last column: probability of non-recovery
  for (t in 1:n_occasions)
    pr_j[t, n_occasions + 1] <- 1 - sum(head(pr_j[t], n_occasions));

  // Define the cell probabilities of the adult m-array
  // Main diagonal
  for (t in 1:n_occasions) {
    pr_a[t,t] <- (1 - sa[t]) * ra[t];
    // Above main diagonal
    for (j in (t + 1):n_occasions)
      pr_a[t, j] <- prod(segment(sa, t, j - t)) * (1 - sa[j]) * ra[j];

    // Below main diagonal
    for (j in 1:(t - 1))
      pr_a[t, j] <- 0;
  } //t

  // Last column: probability of non-recovery
  for (t in 1:n_occasions)
   pr_a[t, n_occasions + 1] <- 1 - sum(head(pr_a[t], n_occasions));
}

model {
  // Priors
  mean_sj ~ uniform(0, 1);         // Prior for mean juv. survival
  mean_sa ~ uniform(0, 1);         // Prior for mean ad. survival
  mean_rj ~ uniform(0, 1);         // Prior for mean juv. recovery
  mean_ra ~ uniform(0, 1);         // Prior for mean ad. recovery

  // Define the multinomial likelihood
  for (t in 1:n_occasions) {
    increment_log_prob(multinomial_log(marr_j[t], pr_j[t]));
    increment_log_prob(multinomial_log(marr_a[t], pr_a[t]));
  }
}
"

## Bundle data
stan_data <- list(marr_j = marr.j,
                  marr_a = marr.a,
                  n_occasions = dim(marr.j)[2]-1)

## Initial values
inits <- function(){list(mean_sj = runif(1, 0, 1),
                         mean_sa = runif(1, 0, 1),
                         mean_rj = runif(1, 0, 1),
                         mean_ra = runif(1, 0, 1))}  

## Parameters monitored
params <- c("mean_sj", "mean_rj", "mean_sa", "mean_ra")

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

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

8.4. Real data example: age-dependent survival in Swiss red kites

mr_mnl_age3_code <- "
data {
  int n_age;
  int marr_j[n_age + 1];
  int marr_a[n_age + 1];
}

parameters {
  real<lower=0,upper=1> sjuv;
  real<lower=0,upper=1> ssub;
  real<lower=0,upper=1> sad;
  real<lower=0,upper=1> rjuv;
  real<lower=0,upper=1> rad;
}

transformed parameters {
  simplex[n_age + 1] pr_a;
  simplex[n_age + 1] pr_j;

  // Define the cell probabilities of the juvenile m-array
  // First element
  pr_j[1] <- (1 - sjuv) * rjuv;

  // Second element
  pr_j[2] <- sjuv * (1 - ssub) * rad;

  // Third and further elements
  for (t in 3:n_age)
    pr_j[t] <- sjuv * ssub * pow(sad, (t - 3)) * (1 - sad) * rad;

  // Probability of non-recovery
  pr_j[n_age + 1] <- 1 - sum(head(pr_j, n_age));

  // Define the cell probabilities of the adult m-array
  // All elements
  for (t in 1:n_age)
    pr_a[t] <- pow(sad, (t - 1)) * (1 - sad) * rad;

  // Probability of non-recovery
  pr_a[n_age + 1] <- 1 - sum(head(pr_a, n_age));
}

model {
  // Priors
  sjuv ~ beta(4.2, 2.8); // Informative prior for juv. survival: Analysis A
//sjuv ~ uniform(0, 1);  // Non-informative for juv. survival prior: Analysis B
  ssub ~ uniform(0, 1);  // Prior for subad. survival
  sad ~ uniform(0, 1);   // Prior for ad. survival
  rjuv ~ uniform(0, 1);  // Prior for juv. recovery
  rad ~ uniform(0, 1);   // Prior for ad. recovery

  // Define the multinomial likelihood
  increment_log_prob(multinomial_log(marr_j, pr_j));
  increment_log_prob(multinomial_log(marr_a, pr_a));
}
"

## Bundle data
stan_data <- list(marr_j = marray.juv,
                  marr_a = marray.ad,
                  n_age = length(marray.juv) - 1)

## Initial values
inits <- function(){list(sjuv = runif(1, 0, 1),
                         ssub = runif(1, 0, 1),
                         sad = runif(1, 0, 1),
                         rjuv = runif(1, 0, 1),
                         rad = runif(1, 0, 1))}  

## Parameters monitored
params <- c("sjuv", "ssub", "sad", "rjuv", "rad")

## MCMC settings
ni <- 4000
nt <- 2
nb <- 2000
nc <- 4

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


タグ:STAn
nice!(1)  コメント(0)  トラックバック(1) 
共通テーマ:日記・雑感

nice! 1

コメント 0

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

Facebook コメント

トラックバック 1