Getting the most of mix models with random slopes

I use mix models as a way to find general patterns integrating different levels of information (i.e. the random effects). Sometimes you only want to focus on the general effects, but others the variation among levels is also of interest. If this is the case, using a random slope model is pretty cool, but making sense of lmer output is not trivial. I provide here code to get the random slopes and CI’s of such models using the iris dataset in R (mainly, because I am sure I will need to check this entry in the future myself)

#This shows how to get the random slopes and CI's for each level in a hierarchical model

Is there a general relationship between petal and sepal width? and how it differs by species?

plot(iris$Sepal.Width ~ iris$Petal.Width, col = iris$Species, las =1
#Our model with random slope and intercept
m2 <- lmer(data = iris, Sepal.Width ~ Petal.Width + (1 + Petal.Width|Species))
#extract fixed effects
#extract random effects
b=ranef(m2, condVar=TRUE)
# Extract the variances of the random effects
qq <- attr(b[[1]], "postVar")
e=e[2,2,] #here we want to access the Petal.Weigth, which is stored in column 2 in b[[1]], that's why I use the [,2,2]
#calculate CI's
#Plot betas and its errors
dotchart(mean_$Petal.Width, labels = rownames(mean_), cex = 0.5,         xlim = c(0.4,1.4), xlab = "betas")
#add CI's...
for (i in 1:nrow(mean_)){
     lines(x = c(liminf[i,1], limsup[i,1]), y = c(i,i)) 
#make final plot
plot(iris$Sepal.Width ~ iris$Petal.Width, col = iris$Species, las = 1)
#and plot each random slope
abline(a = b[[1]][1,1]+a[1], b= mean_$Petal.Width[1], col = "black")
abline(a = b[[1]][2,1]+a[1], b= mean_$Petal.Width[2], col = "red")
abline(a = b[[1]][3,1]+a[1], b= mean_$Petal.Width[3], col = "green")
#and general response
abline(a, lty = 2)

The code is fast and dirty, I know I can be more consistent, and the plot can be nicer (clip the slopes, add SE’s), but will work to get the idea.

& thanks to Lucas Garibaldi for sharing much of this knowledge with me.

And the gist:

Can functional diversity indices predict ecosystem function?

That is exactly the question Vesna and I had after reading a few papers using FD indices as proxies of function. Other than a few papers on plants (mostly experimental and with biomass production as the function studied) we couldn’t find the answer. Then is when we started a side project to gather the data necessary to answer our question. The problem is that there are many FD indexes (see previous post here) and many functions to test. So after realising that this was not a side project anymore we gathered data on as many functions as we can find (a.k.a. calling past collaborators and assaulting ecologists for data in the darker corridors of the conferences we attended) and developed a serious analysis.

The answer? Well, I wouldn’t use FD metrics to predict a single function (i.e. pollination to a single plant species or predation of a single aphid species) because their predictive power are usually below 0.5. However, it worth nothing that non-experimental single functions are usually better described by a few traits than by classic taxonomic diversity indices. It seems it depends on the function studied, but functional identity and to some degree weighted functional diversity metrics can improve our mechanistic understanding of the BEF effects, and this is pretty cool.

Plus, this is the first paper with two lead authors and lots of GitHub commits I did, and it worked really well!

Enjoy the paper here.