random slopes by two groups (species and site)?

I have a dataset of 20 bird species counted at 8 sites, and I'd like to estimate the species-specific change in migration timing at each site. Roughly half of those species are recorded at >= 5 sites; 7-15 species are recorded depending on the site (gray tiles indicate sites where a species is not counted or there is insufficient data):

If I was just focused on a single site, I think it would be sensible to fit a random slopes model (partially pooling by species) to obtain the species-specific change in migration timing (days per year) at that particular site. In this model, 'doy' is the day-of-year, and year is 1985, 1986, etc.:

m1 <- lmer(doy ~ year + (1 + year|sppID), data = df)


Now I'm wondering whether I could 'scale-up' this approach to also incorporate variation in siteID? If I run a linear model for each species x site combination, I not only see that different species are changing their migration timing in different directions, but that a single species can exhibit advancements/delays/no change in timing depending on the site.

My question is whether there's a coherent way to fit all the data with a single hierarchical model to obtain species-specific changes in timing that can also vary by site, and what the correct model specification would look like? I'm assuming too that some species would not be counted across enough sites to include in a model like this? It would be ok with me to drop those rare species in the dataset.

Thanks,

Jay

• Can you explain what the variables doy and year are, and also what "Years Counted" in the figure is ? – Robert Long Feb 5 at 11:37

Why not? :

m1 <- lmer(doy ~ year + siteID + year:siteID + (1 + year|sppID), data = df[sppID=="spp_1"])

If the species have very dissimilar responses and you are not interested in comparing sites, an interesting option would be to perform one model for each species and add sites as random term:

m2 <- lmer(doy ~ year + (1 + year|siteID), data = df[sppID=="spp_1"])

If you want to directly compare phenology trends between species, this model might also be usefull (more thought should be given to the random part specifications):

m3 <- lmer(doy ~ year + sppID + year:sspID + (1 + year:sspID|siteID), data = df)

If you are more interested in differences between sites than differences between species, I believe siteID and sspID could be switch in the m3 model:

m4 <- lmer(doy ~ year + siteID + year:siteID + (1 + year:siteID|sspID), data = df)