Using ulam to run a simple social relations model

Modeling
1

Trying to run my first Bayesian model with a toy example

The model is a simple modification from the social relations model in the statistical rethinking book.

The original model in the book

f_dyad <- alist(
    GAB ~ poisson( lambdaAB ),
    GBA ~ poisson( lambdaBA ),
    log(lambdaAB) <- a + T[D,1] ,
    log(lambdaBA) <- a + T[D,2] ,
    a ~ normal(0,1),

    ## dyad effects - non-centered
    transpars> matrix[N_dyads,2]:T <-
            compose_noncentered( rep_vector(sigma_T,2) , L_Rho_T , Z ),
    matrix[2,N_dyads]:Z ~ normal( 0 , 1 ),
    cholesky_factor_corr[2]:L_Rho_T ~ lkj_corr_cholesky( 2 ),
    sigma_T ~ exponential(1),

    ## compute correlation matrix for dyads
    gq> matrix[2,2]:Rho_T <<- Chol_to_Corr( L_Rho_T )
)

mGD <- ulam( f_dyad , data=sim_data , chains=4 , cores=4 , iter=2000 )

I want to try a linear version of it

f_dyad <- alist(
  GAB ~ normal(mAB,e),
  GBA ~ normal(mBA,e),
  mAB <- T[D,1] ,
  mBA <- T[D,2] ,
  e ~ exponential(1),
  
  ## dyad effects - non-centered
  transpars> matrix[N_dyads,2]:T <-
    compose_noncentered( rep_vector(sigma_T,2) , L_Rho_T , Z ),
  matrix[2,N_dyads]:Z ~ normal( 0 , 1 ),
  cholesky_factor_corr[2]:L_Rho_T ~ lkj_corr_cholesky( 2 ),
  sigma_T ~ exponential(1),
  
  ## compute correlation matrix for dyads
  gq> matrix[2,2]:Rho_T <<- Chol_to_Corr( L_Rho_T )
)

mGD <- ulam( f_dyad , data=sim_data , chains=4 , cores=4 , iter=2000 )

But the model won’t converge and I got bad effective sample size and rhat. What I am doing wrong here?

2

Sorry, @BayesLearner, but it looks like ulam is a package from Richard McElreath (author of Statistical Rethinking), not a part of the Stan project. Here’s their GitHub page:

So I’d suggest raising an issue on their GitHub or figuring out if they have a mailing list. Good luck!

3

The Poisson version is kind of a fudge in some ways because the dyadic effects model the extra-Poisson variation. That is, the dyadic effects are analogous to the observation-level random effects that are sometimes used to model overdispersion in count data, where in this case the random effects are modeled as correlated within dyads.

When doing a Social Relations Model for Gaussian outcomes, by contrast, you’ll want to use the residuals in the same way that a Seemingly Unrelated Regression model would. For instance, your dyadic effects should be distributed with the same residual variance (e) that you’re using for the main likelihood function. It’s the addition of the extraneous residuals and variance that is seemingly competing within the likelihood function that probably results in the poor mixing.

Outside of Stan, there are a couple of software options for fitting Gaussian Social Relations Models using MCMC estimation, including:

Stat-JR as described by Koster et. al: https://journals.sagepub.com/doi/pdf/10.1177/1525822X19889011

amen is an R package by Peter Hoff: [1506.08237] Dyadic data analysis with amen