Bayes@Lund2017 examples

bayes_Lund_2017.R

# script for examples in Bayes@Lund2017 presentation

# joint model example
notes_max <- 10
rate_max <- 5
pm <- matrix( NA , nrow=rate_max+1 , ncol=notes_max+1 )
for ( i in 1:(rate_max+1) )
    for ( j in 1:(notes_max+1) )
        pm[i,j] <- dpois( j-1 , lambda=i ) * (1/(rate_max+1))

# basic case
N_days <- 7
alpha <- 20
beta <- 10
cat <- sample( 0:1 , size=N_days , prob=c(0.5,0.5) , replace=TRUE )
notes <- rpois( N_days , lambda=(1-cat)*alpha + cat*beta )

library(rethinking) # needs Experimental branch (as of 19 April 2017)

dat_list <- list(notes=notes,cat=cat)

m1 <- map2stan(
    alist(
        notes ~ poisson(lambda),
        lambda <- (1-cat)*alpha + cat*beta,
        alpha ~ exponential(0.1),
        beta ~ exponential(0.1)
    ),
    constraints=list(alpha="lower=0",beta="lower=0"),
    data=dat_list )

# GLMM

dat_list_glmm <- dat_list
dat_list_glmm$id <- c(1,1,2,3,3,4,5)

mx <- map2stan(
    alist(
        notes ~ poisson(lambda),
        lambda <- (1-cat)*alpha[id] + cat*beta[id],
        alpha[id] ~ exponential("1.0/alpha_bar"),
        beta[id] ~ exponential("1.0/beta_bar"),
        alpha_bar ~ exponential(0.1),
        beta_bar ~ exponential(0.1)
    ),
    constraints=list(alpha="lower=0",beta="lower=0",alpha_bar="lower=0",beta_bar="lower=0"),
    data=dat_list_glmm , start=list(alpha=rep(1,5),beta=rep(1,5)) , sample=TRUE )

# missing cat data
cat_obs <- cat
cat_obs[c(1,4)] <- NA
dat_list2 <- list( notes=notes , cat=cat_obs )

m2 <- map2stan(
    alist(
        notes ~ poisson(lambda),
        lambda <- (1-cat)*alpha + cat*beta,
        cat ~ bernoulli(kappa),
        kappa ~ beta(4,4),
        alpha ~ exponential(0.1),
        beta ~ exponential(0.1)
    ),
    constraints=list(alpha="lower=0",beta="lower=0",kappa="lower=0,upper=1"),
    data=dat_list2)

# stan version
mcatmiss_code <- "
data{
    int N;
    int notes[N];
    int cat[N];
}
parameters{
    real<lower=0,upper=1> kappa;
    real<lower=0> beta;
    real<lower=0> alpha;
}
model{
    beta ~ exponential( 0.1 );
    alpha ~ exponential( 0.1 );
    kappa ~ beta( 4 , 4 );
    for ( i in 1:N ) {
        if ( cat[i]==-1 ) { // cat missing
            target += log_mix( kappa ,
                    poisson_lpmf( notes[i] | beta ),
                    poisson_lpmf( notes[i] | alpha )
                );
        } else { // cat not missing
            cat[i] ~ bernoulli(kappa);
            notes[i] ~ poisson( (1-cat[i])*alpha + cat[i]*beta );
        } 
    }//i
}//model
generated quantities{
    vector[N] cat_impute;
    for ( i in 1:N ) {
        real logPxy;
        real logPy;
        if ( cat[i]==-1 ) {
            // need P(x|y)
            // P(x|y) = P(x,y)/P(y)
            // P(x,y) = P(x)P(y|x)
            // P(y) = P(x==1)P(y|x==1) + P(x==0)P(y|x==0)
            logPxy = log(kappa) + poisson_lpmf(notes[i]|beta);
            logPy = log_mix( kappa ,
                    poisson_lpmf( notes[i] | beta ),
                    poisson_lpmf( notes[i] | alpha ) );
            cat_impute[i] = exp( logPxy - logPy );
        } else {
            cat_impute[i] = cat[i];
        }
    }//i
}//gq
"

cat_obs <- cat
cat_obs[c(1,4)] <- -1
dat_list4 <- list( notes=notes , cat=cat_obs , N=N_days )

m_miss <- stan( model_code=mcatmiss_code , data=dat_list4 , chains=1 )

# measurement error (occupancy model)
# half of time, present cat not observed by data collector, but seen by bird
cat_obs <- ifelse( cat==0 , 0 , ifelse(runif(N_days)<0.5,0,1) )

dat_list3 <- list(
    notes=notes,
    cat_true=ifelse(cat_obs==1,1,NA),
    cat_obs=cat_obs
)

# doesn't work in map2stan, but would look like this
m3 <- map2stan(
    alist(
        notes ~ poisson(lambda),
        lambda <- (1-cat_true)*alpha + cat_true*beta,
        cat_obs ~ bernoulli( cat_true*delta ),
        cat_true ~ bernoulli(kappa),
        kappa ~ beta(4,4),
        delta ~ beta(4,4),
        alpha ~ exponential(0.1),
        beta ~ exponential(0.1)
    ),
    constraints=list(alpha="lower=0",beta="lower=0"),
    data=dat_list3 )

# working Stan version

m3stan_code <- "
// static occupancy model
data {
    int N;
    int notes[N];
    int cat[N];
}
parameters {
    real<lower=0> beta;
    real<lower=0> alpha;
    real<lower=0,upper=1> kappa;    // prob cat present
    real<lower=0,upper=1> delta;    // prob of detecting cat
}
model {
    beta ~ exponential( 0.1 );
    alpha ~ exponential( 0.1 );
    kappa ~ beta(4,4);
    delta ~ beta(4,4);

    for ( i in 1:N ) {
        if ( cat[i]==1 )
            // cat present and detected
            target += log(kappa) + log(delta) + poisson_lpmf( notes[i] | beta );
        if ( cat[i]==0 ) {
            // cat not observed, but cannot be sure not there
            // marginalize over unknown cat state:
            // (1) cat present and not detected
            // (2) cat absent
            target += log_sum_exp(
                    log(kappa) + log1m(delta) + poisson_lpmf( notes[i] | beta ),
                    log1m(kappa) + poisson_lpmf( notes[i] | alpha ) );
        }//cat==0
    }//i
}
"

dat_list3 <- dat_list
dat_list3$N <- N_days

m3stan <- stan( model_code=m3stan_code , data=dat_list3 , chains=1 )