rethinking R Package

When working with binary (0/1) predictors that contain missing values (NA), Stan cannot sample the discrete parameters directly. map2stan handles this by performing semi-automated marginalization. It constructs a mixture model where the likelihood is the sum of the log-likelihoods conditional on all possible combinations of the missing 0/1 values.

To use this, specify a distribution for the binary variable (e.g., bernoulli(phi)) and ensure you provide constraints for the hyperparameters (e.g., phi must be between 0 and 1).

If you need to recover the unobserved values, set do_discrete_imputation=TRUE. This will compute the posterior probability of each missing value being 1.

# Example: Binary predictor with missingness
N <- 100
N_miss <- 10
x <- rbinom( N , size=1 , prob=0.5 )
y <- rnorm( N , 2*x , 1 )
x[ sample(1:N,size=N_miss) ] <- NA

f6 <- alist(
    y ~ dnorm( mu , sigma ),
    mu <- a + b*x,
    x ~ bernoulli( phi ),
    a ~ dnorm( 0 , 100 ),
    b ~ dnorm( 0  , 10 ),
    phi ~ beta( 1 , 1 ),
    sigma ~ dcauchy(0,2)
)

# Use do_discrete_imputation=TRUE to get imputed probabilities
m6 <- map2stan( f6 , data=list(y=y,x=x) , constraints=list(phi="lower=0,upper=1") ,
      do_discrete_imputation=TRUE )

precis( m6 , depth=2 )

Gaussian processes can be implemented in ulam using the cov_GPL2 macro.

For a standard Gaussian process, you define a covariance matrix SIGMA using cov_GPL2( Dmat , etasq , rhosq , 0.01 ), where Dmat is a distance matrix.

For more complex (e.g., non-centered) Gaussian processes, use the <<- operator for direct assignment to avoid loops during Stan compilation. This allows you to perform linear algebra operations, such as multiplying a Cholesky factor by a vector of random effects, directly within the model definition.

# Simple Gaussian Process example
m_GP1 <- ulam(
    alist(
        y ~ poisson( mu ),
        log(mu) <- a + aj[society] + b*log_pop,
        a ~ normal(0,10),
        b ~ normal(0,1),
        vector[10]: aj ~ multi_normal( 0 , SIGMA ),
        matrix[10,10]: SIGMA <- cov_GPL2( Dmat , etasq , rhosq , 0.01 ),
        etasq ~ exponential(1),
        rhosq ~ exponential(1)
    ),
    data=dat )

The ulam function supports if-then logic and custom distribution assignments using the | (conditional) operator and the custom keyword. This is useful for coding mixture models like zero-inflated Poisson or models with discrete missing values.

When using custom, you are defining custom target updates in Stan. You can use standard Stan functions like log_mix, log1m, and poisson_lpmf within these blocks.

Example: A zero-inflated Poisson model where y|y==0 uses a mixture and y|y>0 uses a standard Poisson likelihood.

# zero-inflated poisson model
m_zip <- ulam(
    alist(
        y|y==0 ~ custom( log_mix( p , 0 , poisson_lpmf(0|lambda) ) ),
        y|y>0 ~ custom( log1m(p) + poisson_lpmf(y|lambda) ),
        logit(p) <- ap,
        log(lambda) <- al + bl*x,
        ap ~ dnorm(0,1),
        al ~ dnorm(0,10),
        bl ~ normal(0,1)
    ) ,
    data=list(y=y,x=x) )

The ulam() function allows for complex multilevel models, including varying intercepts and varying slopes. You can use explicit vector/matrix declarations to handle grouped effects.

Varying Intercepts Example:

m_glmm1 <- ulam(
    alist(
        admit ~ binomial(applications,p),
        logit(p) <- a[dept] + b*male,
        a[dept] ~ normal( abar , sigma ),
        abar ~ normal( 0 , 4 ),
        sigma ~ half_normal(0,1),
        b ~ normal(0,1)
    ), data=UCBadmit )

Varying Slopes with Vector Declaration: To declare a vector of length $N$ for each group, use the syntax vector[N]:name[group_index]. For example, vector[2]:v[dept] declares a vector of length 2 for each unique dept.

m_glmm3 <- ulam(
    alist(
        admit ~ binomial(applications,p),
        logit(p) <- v[dept,1] + v[dept,2]*male,
        vector[2]:v[dept] ~ multi_normal( c(abar,bbar) , Rho , sigma ),
        abar ~ normal( 0 , 4 ),
        bbar ~ normal(0,1),
        sigma ~ half_normal(0,1),
        Rho ~ lkjcorr(2)
    ), data=UCBadmit )