Need help with discrete Poisson parameters

I have a model…

data{
    vector[54] power_Light;
    vector[54] logREMS_Light;
    vector[54] REMS_Light;
    vector[54] logNREMS_Light;
    vector[54] NREMS_Light;
    vector[54] power_Dark;
    vector[54] logREMS_Dark;
    vector[54] REMS_Dark;
    vector[54] logNREMS_Dark;
    vector[54] NREMS_Dark;
    int SWA;
    int dREMS;
    int nREMS;
    int dNREMS;
    int nNREMS;
    int CCL[54];
    int LPS[54];
    int Nobs;
    int n;
    vector[54] SWA_Light;
    vector[54] dREMS_Light;
    int nREMS_Light[54];
    vector[54] dNREMS_Light;
    int nNREMS_Light[54];
    vector[54] SWA_Dark;
    vector[54] dREMS_Dark;
    int nREMS_Dark[54];
    vector[54] dNREMS_Dark;
    int nNREMS_Dark[54];
    vector[54] Wake_Dark;
    int LPSxCCL[54];
    int Animal[54];
}
parameters{
    vector<lower=0>[18] a_nNREMS_Dark;
    real<lower=0> sigma_a_nNREMS_Dark;
    vector<lower=0>[18] a_dNREMS_Dark;
    real<lower=0> sigma_a_dNREMS_Dark;
    vector<lower=0>[18] a_nREMS_Dark;
    real<lower=0> sigma_a_nREMS_Dark;
    vector<lower=0>[18] a_dREMS_Dark;
    real<lower=0> sigma_a_dREMS_Dark;
    vector<lower=0>[18] a_SWA_Dark;
    real<lower=0> sigma_a_SWA_Dark;
    vector[4] a_SWA_Dark_LPSxCCL;
    vector[4] a_dREMS_Dark_LPSxCCL;
    vector[4] a_nREMS_Dark_LPSxCCL;
    vector[4] a_dNREMS_Dark_LPSxCCL;
    vector[4] a_nNREMS_Dark_LPSxCCL;
    corr_matrix[5] Rho_LPSxCCL_Dark;
    vector<lower=0>[5] Sigma_LPSxCCL_Dark;
    vector<lower=0>[18] a_nNREMS_Light;
    real<lower=0> sigma_a_nNREMS_Light;
    real b_nNREMS_DarkWake;
    real<lower=0> sigma_b_nNREMS_DarkWake;
    vector<lower=0>[18] a_dNREMS_Light;
    real<lower=0> sigma_a_dNREMS_Light;
    real b_dNREMS_DarkWake;
    real<lower=0> sigma_b_dNREMS_DarkWake;
    vector<lower=0>[18] a_nREMS_Light;
    real<lower=0> sigma_a_nREMS_Light;
    real b_nREMS_DarkWake;
    real<lower=0> sigma_b_nREMS_DarkWake;
    vector<lower=0>[18] a_dREMS_Light;
    real<lower=0> sigma_a_dREMS_Light;
    real b_dREMS_DarkWake;
    real<lower=0> sigma_b_dREMS_DarkWake;
    vector<lower=0>[18] a_SWA_Light;
    real<lower=0> sigma_a_SWA_Light;
    real b_SWA_DarkWake;
    real<lower=0> sigma_b_SWA_DarkWake;
    vector[4] a_SWA_Light_LPSxCCL;
    vector[4] a_dREMS_Light_LPSxCCL;
    vector[4] a_nREMS_Light_LPSxCCL;
    vector[4] a_dNREMS_Light_LPSxCCL;
    vector[4] a_nNREMS_Light_LPSxCCL;
    corr_matrix[5] Rho_LPSxCCL_Light;
    vector<lower=0>[5] Sigma_LPSxCCL_Light;
    corr_matrix[10] Rho;
    vector<lower=0>[10] Sigma;
}
model{
    vector[54] mu_nNREMS_Dark;
    vector[54] mu_dNREMS_Dark;
    vector[54] mu_nREMS_Dark;
    vector[54] mu_dREMS_Dark;
    vector[54] mu_SWA_Dark;
    vector[54] mu_nNREMS_Light;
    vector[54] mu_dNREMS_Light;
    vector[54] mu_nREMS_Light;
    vector[54] mu_dREMS_Light;
    vector[54] mu_SWA_Light;
    Sigma ~ exponential( 0.1 );
    Rho ~ lkj_corr( 2 );
    Sigma_LPSxCCL_Light ~ exponential( 0.1 );
    Rho_LPSxCCL_Light ~ lkj_corr( 2 );
    {
    vector[5] YY[4];
    vector[5] MU;
    MU = [ 0 , 0 , 0 , 0 , 0 ]';
    for ( j in 1:4 ) YY[j] = [ a_nNREMS_Light_LPSxCCL[j] , a_dNREMS_Light_LPSxCCL[j] , a_nREMS_Light_LPSxCCL[j] , a_dREMS_Light_LPSxCCL[j] , a_SWA_Light_LPSxCCL[j] ]';
    YY ~ multi_normal( MU , quad_form_diag(Rho_LPSxCCL_Light , Sigma_LPSxCCL_Light) );
    }
    sigma_b_SWA_DarkWake ~ exponential( 1 );
    b_SWA_DarkWake ~ normal( 0 , sigma_b_SWA_DarkWake );
    sigma_a_SWA_Light ~ exponential( 1 );
    a_SWA_Light ~ normal( 1.1 , sigma_a_SWA_Light );
    for ( i in 1:54 ) {
        mu_SWA_Light[i] = a_SWA_Light[Animal[i]] + a_SWA_Light_LPSxCCL[LPSxCCL[i]] + b_SWA_DarkWake * Wake_Dark[i];
    }
    sigma_b_dREMS_DarkWake ~ exponential( 0.1 );
    b_dREMS_DarkWake ~ normal( 0 , sigma_b_dREMS_DarkWake );
    sigma_a_dREMS_Light ~ exponential( 0.1 );
    a_dREMS_Light ~ normal( 90 , sigma_a_dREMS_Light );
    for ( i in 1:54 ) {
        mu_dREMS_Light[i] = a_dREMS_Light[Animal[i]] + a_dREMS_Light_LPSxCCL[LPSxCCL[i]] + b_dREMS_DarkWake * Wake_Dark[i];
    }
    sigma_b_nREMS_DarkWake ~ exponential( 0.2 );
    b_nREMS_DarkWake ~ normal( 0 , sigma_b_nREMS_DarkWake );
    sigma_a_nREMS_Light ~ exponential( 0.2 );
    a_nREMS_Light ~ normal( 48 , sigma_a_nREMS_Light );
    for ( i in 1:54 ) {
        mu_nREMS_Light[i] = a_nREMS_Light[Animal[i]] + a_nREMS_Light_LPSxCCL[LPSxCCL[i]] + b_nREMS_DarkWake * Wake_Dark[i];
    }
    sigma_b_dNREMS_DarkWake ~ exponential( 0.1 );
    b_dNREMS_DarkWake ~ normal( 0 , sigma_b_dNREMS_DarkWake );
    sigma_a_dNREMS_Light ~ exponential( 0.1 );
    a_dNREMS_Light ~ normal( 300 , sigma_a_dNREMS_Light );
    for ( i in 1:54 ) {
        mu_dNREMS_Light[i] = a_dNREMS_Light[Animal[i]] + a_dNREMS_Light_LPSxCCL[LPSxCCL[i]] + b_dNREMS_DarkWake * Wake_Dark[i];
    }
    sigma_b_nNREMS_DarkWake ~ exponential( 0.2 );
    b_nNREMS_DarkWake ~ normal( 0 , sigma_b_nNREMS_DarkWake );
    sigma_a_nNREMS_Light ~ exponential( 0.2 );
    a_nNREMS_Light ~ normal( 96 , sigma_a_nNREMS_Light );
    for ( i in 1:54 ) {
        mu_nNREMS_Light[i] = a_nNREMS_Light[Animal[i]] + a_nNREMS_Light_LPSxCCL[LPSxCCL[i]] + b_nNREMS_DarkWake * Wake_Dark[i];
    }
    Sigma_LPSxCCL_Dark ~ exponential( 0.1 );
    Rho_LPSxCCL_Dark ~ lkj_corr( 2 );
    {
    vector[5] YY[4];
    vector[5] MU;
    MU = [ 0 , 0 , 0 , 0 , 0 ]';
    for ( j in 1:4 ) YY[j] = [ a_nNREMS_Dark_LPSxCCL[j] , a_dNREMS_Dark_LPSxCCL[j] , a_nREMS_Dark_LPSxCCL[j] , a_dREMS_Dark_LPSxCCL[j] , a_SWA_Dark_LPSxCCL[j] ]';
    YY ~ multi_normal( MU , quad_form_diag(Rho_LPSxCCL_Dark , Sigma_LPSxCCL_Dark) );
    }
    sigma_a_SWA_Dark ~ exponential( 1 );
    a_SWA_Dark ~ normal( 0.9 , sigma_a_SWA_Dark );
    for ( i in 1:54 ) {
        mu_SWA_Dark[i] = a_SWA_Dark[Animal[i]] + a_SWA_Dark_LPSxCCL[LPSxCCL[i]];
    }
    sigma_a_dREMS_Dark ~ exponential( 0.1 );
    a_dREMS_Dark ~ normal( 90 , sigma_a_dREMS_Dark );
    for ( i in 1:54 ) {
        mu_dREMS_Dark[i] = a_dREMS_Dark[Animal[i]] + a_dREMS_Dark_LPSxCCL[LPSxCCL[i]];
    }
    sigma_a_nREMS_Dark ~ exponential( 0.2 );
    a_nREMS_Dark ~ normal( 48 , sigma_a_nREMS_Dark );
    for ( i in 1:54 ) {
        mu_nREMS_Dark[i] = a_nREMS_Dark[Animal[i]] + a_nREMS_Dark_LPSxCCL[LPSxCCL[i]];
    }
    sigma_a_dNREMS_Dark ~ exponential( 0.1 );
    a_dNREMS_Dark ~ normal( 300 , sigma_a_dNREMS_Dark );
    for ( i in 1:54 ) {
        mu_dNREMS_Dark[i] = a_dNREMS_Dark[Animal[i]] + a_dNREMS_Dark_LPSxCCL[LPSxCCL[i]];
    }
    sigma_a_nNREMS_Dark ~ exponential( 0.2 );
    a_nNREMS_Dark ~ normal( 96 , sigma_a_nNREMS_Dark );
    for ( i in 1:54 ) {
        mu_nNREMS_Dark[i] = a_nNREMS_Dark[Animal[i]] + a_nNREMS_Dark_LPSxCCL[LPSxCCL[i]];
    }
    {
    vector[10] YY[54];
    vector[10] MU[54];
    for ( j in 1:54 ) MU[j] = [ mu_nNREMS_Dark[j] , mu_dNREMS_Dark[j] , mu_nREMS_Dark[j] , mu_dREMS_Dark[j] , mu_SWA_Dark[j] , mu_nNREMS_Light[j] , mu_dNREMS_Light[j] , mu_nREMS_Light[j] , mu_dREMS_Light[j] , mu_SWA_Light[j] ]';
    for ( j in 1:54 ) YY[j] = [ nNREMS_Dark[j] , dNREMS_Dark[j] , nREMS_Dark[j] , dREMS_Dark[j] , SWA_Dark[j] , nNREMS_Light[j] , dNREMS_Light[j] , nREMS_Light[j] , dREMS_Light[j] , SWA_Light[j] ]';
    YY ~ multi_normal( MU , quad_form_diag(Rho , Sigma) );
    }
}

…where there are a few parameters that are currently a truncated normal but would like to make Poisson; these are a_nNREMS_Dark, a_nREMS_Dark, a_nNREMS_Light, and a_nREMS_Light. The problem is that I’m quite shaky with the whole marginalizing-out business. Can anybody help me with this?

Thanks!

2 Likes

The model as it stands now is much too complex for me to make any suggestions. I suggest you trim it down to the minimal working example so we can more clearly see what needs to be done. You also have a bunch of hard-coded stuff, like 54 everywhere. I suggest you pass that as data so your program is generalisable to any number of animals.

1 Like