R/GeneraterKcode.R

In femoerman/PBPGM: Parallellised Bayesian Population Growth Models

#' GeneraterKcode function
#'
#' Generates the stan code to be used for the models, basd on the prior information provided. This allows to only compile one model for all data, saving a lot of computing time.
#' @param dd Dataframe with columns for population size data (colname=popsize), time data (colname=time) and unique identifiers for each population (colname=ident)
#' @param K.prior Prior value for the mean carrying capacity (K). Must be on log scale and numeric
#' @param Ksd.prior Prior value for the standard deviation on carrying capacity (K). Must be on log scale and numeric
#' @param r0.prior Prior value for the mean intrinsic rate of growth (r0). Must be on log scale and numeric
#' @param r0sd.prior Prior value for the standard deviation on intrinsic rate of growth (r0). Must be on log scale and numeric
#' @param N0.prior Prior value for the mean starting population size (N0). Must be on log scale and numeric
#' @param N0sd.prior Prior value for the standard deviation on starting population size (N0). Must be on log scale and numeric
#' @param sdev.prior Prior value for the standard deviation for the model fitting. Must be numeric. Defaults to 1.
#' @keywords Beverton Holt Stan code
#' @export
#' @examples
#' GeneratrKHcode()
GeneraterKcode <- function(K.prior, Ksd.prior, r0.prior, r0sd.prior, N0.prior, N0sd.prior, sdev.prior){
  part1 <- '
  // function that caclulates r-alpha logistic population growth
  functions{
    real[] odemodel(real t, real[] N, real[] p, real[] x_r, int[] x_i){
    // p[1]=r0, p[2]=K
    real dNdt[1];
      dNdt[1] = p[1] * (1 - (N[1]/p[2])) * N[1];
    return dNdt;
    }
  }

  data{
    int n;
    real log_N0;
    real log_N[n];
    real t0;
    real t[n];
  }

  transformed data {
    // not used here
    real x_r[0];
    int x_i[0];
  }

  parameters{
  real log_r;
  real log_K;
  real log_N0sim;
  real<lower=0> sdev;
  }

transformed parameters{
  // all of this was in the model section previously
  // I moved it here to be able to get waic because Nsim needs to be accessible in "generated quantities {}"

  real p[2];
  real Nsim[n,1]; // simulated values, matrix. dim1 = time, dim2 = dim_ODE = 1
  real N0sim_dummy[1]; // just a dummy because the ODE solver requires real[] instead of real

  // parameters for integrator
  p[1] = exp(log_r);
  p[2] = exp(log_K);
  N0sim_dummy[1] = exp(log_N0sim);  // see above

  // integrate ODE
  Nsim = integrate_ode_rk45(odemodel,N0sim_dummy,t0,t,p,x_r,x_i);
}

  model{
  // priors
  // note: it can be VERY helpful to estimate parameters on logscale,
  // especially if they have different orders of magnitude. here it works on regular scale, though.'
  part2 <- paste("log_r ~ normal(", r0.prior, ", ", r0sd.prior,");", sep="")
  part3 <- paste("log_K ~ normal(", K.prior, ", ", Ksd.prior,");", sep="")
  part4 <- paste("log_N0sim ~ normal(", N0.prior, ", ", N0sd.prior,");", sep="")
  part5 <- paste("sdev ~ cauchy(0, ", sdev.prior,");", sep="")
  part6 <-
  '
  // likelihood, normal (maybe lognormal helpful)
  log_N0 ~ normal(log_N0sim,sdev);
  for (i in 1:n){
  log_N[i] ~ normal(log(Nsim[i,1]),sdev);
  }
  }

  // cacluate log lik to get waic
  // from loo R package description
  generated quantities {
  real log_lik[n];
  for (nn in 1:n){
  log_lik[nn] = normal_lpdf(log_N[nn] | log(Nsim[nn,1]), sdev);
  }
  }
  '

  rK.code <- paste(part1, part2, part3, part4, part5, part6, sep="\n")
  return(rK.code)
}