Selecting best binomial model. Case of flowers pollinated by bees

Question

I'm looking to build a binomial model in R for a fascinating experiment involving flower pollination at various temperatures. In this experiment, plants exhibit a dynamic range of flowers and floral stems. Some plants even boast different floral stems, each with a varying number of flowers.

As bees work their magic and successfully pollinate these flowers, they transform into fruits. My curiosity leads me to explore whether the probability of pollination, or the generation of fruits, is influenced by temperature. I'd love to delve deeper into this and hear your thoughts on the matter!

For this, I use the glmer function from the lme4 package.

Now, I'll describe how my database is structured.

library(Matrix)
library(lme4)


temperature <- rep(seq(25, length.out = 8), each = 5)

flowers <- c(7,8,14,9,9,14,10,5,7,8,13,10,10,9,7,14,10,3,11,8,6,8,6,7,7,4,12,9,6,13,10,8,12,12,11,11,11,9,8,8)

fruits <- c(1,0,5,0,0 ,3 ,2 ,0 ,1 ,0 ,0 ,3 ,0 ,5 ,3 ,1 ,2 ,0 ,4 ,0 ,5 ,0 ,6 ,7 ,7 ,0 ,5 ,8 ,4 ,1 ,2 ,4 ,5 ,6 ,9 ,7 ,3 ,5 ,0 ,8)

ID <-  c("plant 1",  "plant 2",  "plant 3",  "plant 3",  "plant 3",
         "plant 4",  "plant 4",  "plant 5",  "plant 5",  "plant 6",
         "plant 7",  "plant 8", "plant 9",  "plant 10", "plant 10", 
         "plant 11",  "plant 11",  "plant 11","plant 12","plant 12",
         "plant 13",  "plant 13", "plant 14", "plant 14","plant 14",
         "plant 15","plant 15",  "plant 16", "plant 16", "plant 16",
         "plant 17", "plant 17", "plant 18","plant 18", "plant 19",
         "plant 20", "plant 20", "plant 21", "plant 21","plant 22")

spike <- c(1,1,1,2,3,
           1,2,1,2,1,
           1,1,1,1,2,
           1,2,3,1,2,
           1,2,1,2,3,
           1,2,1,2,3,
           1,2,1,2,1,
           1,2,1,2,1)

df_flowers <- data.frame(temperature,flowers,fruits,ID, spike)

As you can see, the variable "spike" simply indicates, when applicable, the number corresponding to each floral stem. Since the floral stem depends on the same flower, I believe it should be nested within each plant.

As mentioned in some literature, one can use the cbind function to create a small matrix with proportions. This is done as follows:

response2<- cbind(df_flowers$flowers,df_flowers$flowers- df_flowers$fruits)

Now, let's create the models, starting with the simplest one, which does not consider the plant ID as a random effect. idk if thats correct, or if i need to begin with the complex one

model1<- glm(response2 ~ temperature, family = binomial, data = df_flowers)
summary(model1)

We can observe that the probability of pollination is indeed influenced by temperature. However, this analysis doesn't account for the previously described implications.

Hence, we'll create a second model, incorporating the plant ID as a random variable:

model2<-glmer(response2 ~ temperature + (1|ID), family =  binomial, data = df_flowers)
summary(model2)

I think we've progressed quite well up to this point, and I would appreciate your support at this stage. I'm not sure if I did the statistically correct thing by nesting 'spike.' It's a doubt I have about whether this is correct. Can you help/explain me?

model3<- glmer(response2 ~ temperature + (1|ID/spike), family =  binomial, data = df_flowers)
summary(model3)

Also, I have the issue that this is a simplified example, but in my real data, I have another variable that could explain pollination. It is the abundance of bees at each experimental site. I was thinking of including it as a covariate with the symbol +, but I really don't know how these approaches should be done.

Should I first try simple models and then make them more complex? Or does it depend on whether we find significant differences? It could also be started with the most complex and then remove variables that do not explain.

Thank you very much for your attention and your responses

Answer 1

There are quite a few separate topics to address here.

Response matrix

Per the glm documentation, you can provide the response for a binomial GLM "as a two-column matrix with the columns giving the numbers of successes and failures". Your response2 object is now programmed as the number of trials and the number of failures, so I'm quite sure this should be the following instead:

response2 <- cbind(df_flowers$fruits, df_flowers$flowers - df_flowers$fruits)

Random effects or not?

You seem to be unsure whether a standard GLM or a GLMM (mixed effects) is more appropriate. This is not usually a question you answer based on your observed data, but rather on a hierarchy in the probable data-generating mechanism.

Specifically, the standard GLM assumes that every observation - every flower of every plant - is independent. It is quite likely however that some plants are just 'naturally better' than others at making nice flowers, or attracting bees, or actually turning flowers into fruit, so assuming that each and every flower is described by the same parameter regardless of which plant they came from might be a strong assumption. If that turns out to be false the standard GLM will be misspecified, giving biased variance estimates and possibly biased parameter estimates as well.

As a rule of thumb it's usually better to have a random effect that may exist than not to have one that does exist. Testing whether you need a random effect is not particularly straightforward because a redundant one will be on the boundary of the parameter space (zero variance or $\sigma_R=0$), see for example this question for some more discussion. Based on your sampling scheme alone (plants forming natural 'clusters' of flowers) I'm quite convinced you should be including a random effect, and the observed fit ($\sigma_R\gg 0$) only confirms this.

Random effects structure

As just discussed, random effects model hierarchy in your data: in your case you have flowers that cluster within plants, and a random intercept per plant (1|ID) allows flowers within one plant to be 'more similar' than ones from different plants.

Extending this to an intercept for flower nested in plant (1|ID/spike) (actually plant nested in flower but your coding makes these the same) will try to fit a random intercept for every flower/plant combination, which in your dataset is for every observation. This random effect is fully confounded with the residual variance so suffers from technical issues in producing the fit, but more importantly it is not meaningful (I'm ignoring specific technical cases here where this might produce side effects like 'robustifying' variance estimates). If you had multiple observations within one flower this could have been another clustering level and a relevant extension, right now it makes very little sense.

Additional covariates

You mentioned you wanted to include the abundance of bees as a potential covariate. This makes a lot of sense, because the assumption that more bees = more pollination seems very reasonable. There's at least two ways to approach this problem depending on whether you're in a more exploratory or confirmatory setting:

• when exploring you're not sure which covariates might be most relevant to your outcome, and you want to test the impact of their in/exclusion. It sounds like this where you are. For one covariate a likelihood ratio test (fitting your models with maximum likelihood!) will do this for you. If you have more predictors a better way to perform selection is probably some kind of lasso/ridge regression, as you have to account for the interdependence between predictors which you won't capture by stepwise selection. Quite a lot has been written about this already, I suggest you search around.
• when confirming you would rather specify a model upfront and leave predictors in as nuisance parameters, regardless of their magnitude or statistical significance. For this you'd usually have some prior data to corroborate that bee abundance and pollination may indeed be related, and simply want to account for that in your confirmatory model.

Purely technically you can indeed include this numeric predictor as response2 ~ temperature + bee_abundance, where the parameter estimate will give you the log odds ratio for a one unit increase in whatever bee_abundance expresses.