0

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
)
Alison
  • 11
  • 3
  • I'm trying to understand what you're trying to do here. Do you, for example, have multiple measurements for each tree and each tree is a group? Or just species as a group and each species has multiple observed tree measurements? – ssp3nc3r Aug 02 '20 at 14:39
  • I have many measurements for each tree, and each tree is associated with a species and a site (which represents a location for a group of trees). There are multiple trees in each species category (two species total) and multiple trees of each species at each site (six sites total). – Alison Aug 03 '20 at 15:08

0 Answers0