in R rethinking package, what is syntax for nested indices?

137 Views Asked by At

I'm running the following model using Richard McElreath's "rethinking" package in R. Each species (represented by species_dummy) has several trees associated with it (represented by the index [tree]). I would like to use nested indices or somehow use the species dummy variable to allow two different values of r_sigma_species, one for each species. I can't figure out how to make the indexing work. I don't believe any example like this is covered in the Statistical Rethinking book by McElreath.

I've tried the following:

  • r_tree[tree_index] ~ dnorm(0, r_sigma_species[species_dummy]) ... but this does not work because rethinking will not accept indices to the right of the ~

  • defining r_sigma_oak and r_sigma_juniper separately (i.e. the two species) and then using the dummy variable: r_sigma_tree[tree_index] ~ r_sigma_juniper * (1-species_dummy) + r_sigma_oak * species_dummy .... but rethinking won't interpret this because it doesn't automatically filter the data for species_dummy based on [tree_index], so there is a mismatch in the number of lines of data.

  • trying to combine the indices as in r_sigma_tree[tree_index][species_dummy] ~ r_sigma_species[species_dummy] ... rethinking doesn't seem to like combined indices.

Since I already have bunch of versions of the model using the "rethinking" package, I do not want to switch to rstan or any other syntax at this point.

model <- map2stan(
  alist(
    avg_v_cent ~ dnorm(mu,sigma),
    mu <- a * exp(r * day_cent),
    a <- a_base + a_species*species_dummy + a_tree[tree_index],
    r <- r_base + r_species*species_dummy + r_tree[tree_index][species_dummy],
    
    a_base ~ dnorm(0,1),
    a_species ~ dnorm(0,1),
    a_tree[tree_index] ~ dnorm(0, a_sigma_tree),
    
    r_base ~ dnorm(0,1),
    r_species ~ dnorm(0,1),
    r_tree[tree_index][species_dummy] ~ dnorm(0, r_sigma_species[species_dummy]),
    
    sigma ~ dcauchy(0,1), 
    a_sigma_tree ~ dcauchy(0,1), 
    r_sigma_species[species_dummy] ~ dcauchy(0,1)
  ),
  data=as.data.frame(dataset_2017_2.trim.test), 
  iter=2000,
  warmup = 1000, 
  control = list(max_treedepth = 13, adapt_delta = 0.9),
  chains = 1
)
0

There are 0 best solutions below