ASKurz / Statistical_Rethinking_with_brms_ggplot2_and_the_tidyverse_2_ed

The bookdown version lives here:
https://bookdown.org/content/4857/
Creative Commons Zero v1.0 Universal
125 stars 37 forks source link

Model 16.2 #17

Open ASKurz opened 3 years ago

ASKurz commented 3 years ago

In Section 16.3.2, McElreath fit a mixture of unobserved strategies. On page 534, we see his statistical model was:

Screen Shot 2020-11-26 at 9 49 16 AM

He fit the model using rstan::stan() and already has his Stan model code saved in the rethinking package. Here's how to fit his model:

# load
library(rethinking)

# data
data(Boxes)
data(Boxes_model)

# prep data 
dat_list <- list(
  N = nrow(Boxes),
  y = Boxes$y,
  majority_first = Boxes$majority_first )

# run the sampler
m16.2 <- stan( model_code=Boxes_model , data=dat_list , chains=3 , cores=3 )

If you execute cat(Boxes_model), you can get McElreath's Stan code:

data{
    int N;
    int y[N];
    int majority_first[N];
}
parameters{
    simplex[5] p;
}
model{
    vector[5] phi;

    // prior
    p ~ dirichlet( rep_vector(4,5) );

    // probability of data
    for ( i in 1:N ) {
        if ( y[i]==2 ) phi[1]=1; else phi[1]=0; // majority
        if ( y[i]==3 ) phi[2]=1; else phi[2]=0; // minority
        if ( y[i]==1 ) phi[3]=1; else phi[3]=0; // maverick
        phi[4]=1.0/3.0;                         // random
        if ( majority_first[i]==1 )             // follow first
            if ( y[i]==2 ) phi[5]=1; else phi[5]=0;
        else
            if ( y[i]==3 ) phi[5]=1; else phi[5]=0;

        // compute log( p_s * Pr(y_i|s )
        for ( j in 1:5 ) phi[j] = log(p[j]) + log(phi[j]);
        // compute average log-probability of y_i
        target += log_sum_exp( phi );
    }
}

To my eye, it looks like one might be able to fit a model like this with brms if you make a custom likelihood (see here). However, that would be well beyond my current skill set.