What is the general process of choosing, confirming, and supporting the distribution used in a generalized linear model?

Question

I am trying to better understand the process of choosing and validating a distribution for a generalized linear model (glm). I understand that for the most part, you can narrow it down to a few distributions based on your knowledge of how the data was collected or what processes produced the data you collected. In a less concrete manner (in my current opinion), you can narrow it down based on certain characteristics of the data (e.g. if your response variable is yes/no or success/fail your probably going to use a binomial distribution). What is less clear to me is how we can support or confirm the use of a chosen distribution, and how to decide between more than one candidate distribution when you are not sure.

For instance, consider this data on plant diversity in response to a fully crossed treatments of fertilizer and light in grassland systems:

>dput(plants)
structure(list(Fertilizer = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("F-", "F+"
), class = "factor"), Light = structure(c(1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("L-", 
"L+"), class = "factor"), FL = structure(c(1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L), .Label = c("F-L-", 
"F-L+", "F+L-", "F+L+"), class = "factor"), LF = structure(c(1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L), .Label = c("L-F-", 
"L-F+", "L+F-", "L+F+"), class = "factor"), Diversity = c(6L, 
7L, 10L, 9L, 5L, 9L, 7L, 6L, 10L, 9L, 9L, 11L, 9L, 9L, 7L, 7L, 
4L, 4L, 6L, 4L, 5L, 5L, 4L, 4L, 7L, 8L, 9L, 9L, 10L, 10L, 10L, 
7L)), class = "data.frame", row.names = c(NA, -32L))

In these grassland systems, there is typically a loss of species diversity in response to fertilization, which may be due to light competition. My goal is to estimate whether the loss of species can be prevented by restoring light to these areas. There were a total of 32 plots, 8 received fertilizer and light(F+L+), 8 received fertilizer and no light (F+L-), 8 received no fertilizer and light (F-L+), and 8 received no fertilizer and no light(F-L-). The number of species was counted in each plot.

So with this being count data, I would think to use the Poisson distribution: PlantMod <- glm(diversity~FL, data = plants, family=Poisson(link="log)) note that is the same thing as: glm(diversity~Fertilizer+Light+Fertilizer*Light), because of the way the data is set up.

So now how I know the Poisson distribution was the best (or a good) choice, or if it would be better to switch to another distribution for count data? Do people typically just try all the models they can think of and see which one fits best? What if they were better off defining a new distribution?

Answer 1

I would try fitting some models with different families and looking at the outputs (z's, likelihoods etc). You could plot residuals vs fitted and observed vs fitted, and assess visually which models work best.

In addition, make sure that your predictions always make scientific sense i.e. if your model can return predictions of negative diversity, that wouldn't make scientific sense, so don't use that model.