Hi

I’m trying to fit a bayesian model to do an analysis similar to Fusaroli et al. (2012): *Coming to Terms: Quantifying the Benefits of Linguistic Coordination.*

I’ve got data from a similar experiment where dyads complete a forced-choice task individually and together at various stimulus intensity levels, and with various stimuli (different vowels in my case).

In the paper above, they fitted psychometric models for individuals and dyads separately, and compared the estimated slopes to quantify the benefit they got from working together as a dyad compared to working individually.

Finally, they find that collective benefit is predicted by various features of the conversation between the dyad partners.

I’d like to combine all of this in stan in order to carry the uncertainty all the way through.

So - I’ve got two identical hierachical (dyad and stimulus) psychometric models, one fitting on the individual data, another on the dyad data.

Then in the transformed parameters block, I have basically, for each dyad:

collective_benefit = inv_logit(dyad slope) / inv_logit(individual slope);

Finally, I’d like to put a linear model on collective_benefit.

collective_benefit ~ normal(mu, sigma);

mu = alpha + beta * data;

[…]

It always gives me Jacobian warnings, but I can get as far as modelling it like a linear regression with just an intercept.

But as soon as I try to add (z-scaled) predictors, it spits out hundreds of errors and contains no samples.

Rejecting initial value: Error evaluating the log probability at the initial value.

Reading through the mailing list led me to the Change of variables vs transformations section of the reference manual, and as far as I can tell I’m doing a change of variables on collective benefit.

But I’m not sure I’m reading it right, and even then I have no idea how to get to a “log absolute determinant of the Jacobian of the transform” (in over my head here).

I also wondered if there was simply too much overlap in the collective benefit for the different dyads to support putting a model on it, but based on the posterior of the intercept-only model (see dropbox at the bottom), I don’t think that’s a problem.

Help, please?

Model below:

data{

// number of data points, groups, and vowels

int<lower=1> N;

int<lower=1> N_group;

int<lower=1> N_vowel;// binomial data for the two psychometric functions

int k_ind[N];

int k_gro[N];

int n_ind[N];

int n_gro[N];// input data

int group[N];

int intensity[N];

int vowel[N];// collective benefit model

real local_confidence[N_group];

}parameters{

// specify var-covar matrices

vector[2] gv_group[N_group];

vector[2] gv_vowel[N_vowel];

vector[2] iv_group[N_group];

vector[2] iv_vowel[N_vowel];corr_matrix[2] gRho_group;

corr_matrix[2] gRho_vowel;

corr_matrix[2] iRho_group;

corr_matrix[2] iRho_vowel;// scale parameters for the random effects

vector<lower=0>[2] gsigma_group;

vector<lower=0>[2] gsigma_vowel;

vector<lower=0>[2] isigma_group;

vector<lower=0>[2] isigma_vowel;// main effects

real ga;

real gb;

real ia;

real ib;// =========================================================

// parameters for the collective benefit model

real<lower=0> collective_sigma;

real col_a;

real col_b_confidence;

}transformed parameters{

real<lower=0> collective_benefit[N_group];

for ( k in 1:N_group ) {

collective_benefit[k] =

inv_logit(gv_group[k,2] * gsigma_group[2]) // group slope

/ inv_logit(iv_group[k,2] * isigma_group[2]); // individual slope

}

}model{

// psychometric model for individuals and groups

// logit(theta) = A + B * stimulus intensity

vector[N] gtheta;

vector[N] gA;

vector[N] gB;vector[N] itheta;

vector[N] iA;

vector[N] iB;real collective_mu[N];

// priors for the “main effects” in the psychometric models

ia ~ normal(0, 1);

ib ~ normal(0, 1);

ga ~ normal(0, 1);

gb ~ normal(0, 1);// correlations of the random effects

gv_group ~ multi_normal(rep_vector(0,2), gRho_group);

gv_vowel ~ multi_normal(rep_vector(0,2), gRho_vowel);

iv_group ~ multi_normal(rep_vector(0,2), iRho_group);

iv_vowel ~ multi_normal(rep_vector(0,2), iRho_vowel);gRho_group ~ lkj_corr( 4 );

gRho_vowel ~ lkj_corr( 4 );

iRho_group ~ lkj_corr( 4 );

iRho_vowel ~ lkj_corr( 4 );// scale of the random effects

gsigma_group ~ cauchy( 0 , 2 );

gsigma_vowel ~ cauchy( 0 , 2 );

isigma_group ~ cauchy( 0 , 2 );

isigma_vowel ~ cauchy( 0 , 2 );// psychometric function for both individuals and groups

for ( i in 1:N ) {

gA[i] = ga + gv_vowel[vowel[i],1] * gsigma_vowel[1] + gv_group[group[i],1] * gsigma_group[1];

gB[i] = gb + gv_vowel[vowel[i],2] * gsigma_vowel[2] + gv_group[group[i],2] * gsigma_group[2];

gtheta[i] = gA[i] + gB[i] * intensity[i];`iA[i] = ia + iv_vowel[vowel[i],1] * isigma_vowel[1] + iv_group[group[i],1] * isigma_group[1]; iB[i] = gb + iv_vowel[vowel[i],2] * isigma_vowel[2] + iv_group[group[i],2] * isigma_group[2]; itheta[i] = iA[i] + iB[i] * intensity[i];`

}

k_gro ~ binomial_logit(n_gro, gtheta);

k_ind ~ binomial_logit(n_ind, itheta);// here comes the collective benefit stuff

for (k in 1:N_group) {

collective_mu[k] = col_a + col_b_confidence * local_confidence[k];

}col_a ~ normal(1,1);

col_b_confidence ~ normal(0,1);

collective_sigma ~ cauchy(0,2);collective_benefit ~ normal(collective_mu, collective_sigma);

// Specifying only an intercept works, but …

// collective_benefit ~ normal(col_a, collective_sigma);

}

data on dropbox

R version 3.3.2 (2016-10-31)

Platform: x86_64-pc-linux-gnu (64-bit)

StanHeaders_2.15.0-1

rstan_2.15.1

clang version 4.0.0-1ubuntu1

Thanks in advance

- Malte