6  Exercise: Arabidopsis

To practice we will use an adapted data set on Arabidopsis from Banta, Stevens and Pigliucci (2010 Oikos 119(2), 359-369).

If you have gone through the exercise on linear mixed models, you will know this data set. This growth chamber experiment is about the impact of removing the apical meristem and nutrient level on rosette growth and fruit formation.

You can download the data set here. It is an R object. Store it on your computer and load it via this (after adapting path):

load("data/arabidopsis.obj")

In this exercise we will now model de number of fruits instead of on the rosette area.

The task is:

If you haven’t done this before:

• Find out from the data how the experiment is set up.
• Comment on the set-up and the contribution of the factors

For the number of fruits:

• Find out what family would be suitable
• Fit an appropriate mixed effect model (or several)
• Provide a meaningful summary of the analysis
• Make conclusions about the main goal(s) of the experiment

Some of the code is hidden behind the buttons below. Don’t reveal the solutions unless you really stuck

6.1 Discover the data

library(tidyverse)

str(arabidopsis)

arabidopsis <- arabidopsis %>% mutate(rack = as.factor(rack))

ggplot(arabidopsis,aes(x=total.fruits)) + 
  geom_histogram(bins=50) + 
  facet_grid(clipping~nutrient)

ggplot(arabidopsis,aes(x=total.fruits)) + 
  geom_histogram(bins=50) + 
  facet_grid(rack~region)

6.2 Model the data

library(brms)
library(rstan)
library(tidybayes)
rstan_options(auto_write=TRUE)

mhg.brm <- brm(data=arabidopsis, 
              bf(total.fruits ~ nutrient*clipping + (1|region/population)  + (1|rack),
                 hu ~ nutrient*clipping),
              family = hurdle_gamma(link="log"),
              chains = 4,cores = 4, iter=2000, warmup = 500,
              control = list(max_treedepth = 10,adapt_delta = .97))

summary(mhg.brm)

mhln.brm <- brm(data=arabidopsis, 
              bf(total.fruits ~ nutrient*clipping + (1|region/population)  + (1|rack),
                 hu ~ nutrient*clipping),
              family = hurdle_lognormal(),
              chains = 4,cores = 4, iter=2000, warmup = 500,
              control = list(max_treedepth = 10,adapt_delta = .97))

summary(mhln.brm)

mhg2.brm <- brm(data=arabidopsis, 
              bf(total.fruits ~ nutrient*clipping + (1|region/population)  + (1|rack),
                 hu ~ nutrient*clipping + (1|region)  + (1|rack)),
              family = hurdle_gamma(),
              chains = 4,cores = 4, iter=2000, warmup = 500,
              control = list(max_treedepth = 10,adapt_delta = .97))

summary(mhg2.brm)

loo(mhg.brm,mhln.brm,mhg2.brm)

6.3 Output on fixed effects

predframe <- arabidopsis %>% select(clipping,nutrient) %>% distinct()

hg2.pred <- bind_cols(predframe,predict(mhg2.brm,newdata=predframe,re_formula = NA,scale="response"))

# fitted values for full hurdle model
hg2.fitted <- bind_cols(predframe,fitted(mhg2.brm,newdata=predframe,re_formula = NA,scale="response"))
hg2.fitted
# fitted values for both components (the hurdle part and the gamma part)
hg2.hu <- bind_cols(predframe,fitted(mhg2.brm,newdata=predframe,re_formula = NA,dpar="hu",scale="response"))
hg2.gm <- bind_cols(predframe,fitted(mhg2.brm,newdata=predframe,re_formula = NA,dpar="mu",scale="response"))
hg2.hu
hg2.gm