8
$\begingroup$

I am thinking about building a model (glm) where the response variable (y) is the cover (in per cent) of a plant species in a defined area, dependant of environmental variables. However, I don't think it is a binomial link, as the output is not a probability between 0 and 1. What would you suggest to use as a link function?

$\endgroup$
2
  • 3
    $\begingroup$ Sometimes theory will suggest a link. See stats.stackexchange.com/a/64039/919 for an example. $\endgroup$
    – whuber
    Feb 5, 2021 at 17:20
  • 2
    $\begingroup$ @whuber: Bill, we are like two peas in the same pod - we naturally gravitate towards the same questions! 😁 $\endgroup$ Feb 5, 2021 at 18:07

1 Answer 1

6
$\begingroup$

Nice question!

You can start your modelling efforts by converting your cover variable so that it is expressed as a proportion rather than a percentage (i.e., simply divide the values of your original cover variable by 100).

Once you get the converted variable, look to see how many 0 and 1 values you have, if any, in addition to values falling in the (0,1) interval.

If you don't have any 0 or 1 values, you can use beta regression modelling. A beta regression model is not a glm model but it has similarities to a glm model. The beta regression model has a logit link function; this means that you are modelling the logit-transformed mean cover as a function of various predictor variables. (Recall that cover is now expressed as a proportion.)

If you have a substantial number of 0 values but no 1 values, you will need to use a zero-inflated beta regression model for your modelling.

If you have a substantial number of 1 values but no 0 values, you will need to use a one-inflated beta regression model.

If you have a substantial number of 0 values and a substantial number of 1 values, you will need to use a zero-and-one-inflated beta regression model.

Do you use R? If yes, the gamlss package is your best bet for implementing these kinds of models - see, for instance, page 107 of this document:

https://www.gamlss.com/wp-content/uploads/2013/01/book-2010-Athens1.pdf.

The inflated versions of the beta regression models have multiple links because they simultaneously model multiple parameters. Some of these links can be logit links, some can be log links, etc. In gamlss, the family distribution for these models can be one of:

  1. BE for beta regression;

  2. BEINF0 for zero-inflated beta regression;

  3. BEINF1 for one-inflated beta regression;

  4. BEINF for zero-and-one inflated beta regression.

Addendum

Both betareg() and gamlss() afford the ability to fit a beta regression model to your data and then use the model for prediction purposes. The R code below provides an example for how you can do this. Note that there are some differences in output between the model summaries produced by the two functions from the same dataset.

#===========================================================
# Load Gasoline data into R 
#===========================================================

library(betareg)

data("GasolineYield", package = "betareg")

#===========================================================
# Fit beta regression model to Gasoline data using betareg()
#===========================================================

model_betareg <- betareg(yield ~ temp, 
                         link = "logit", 
                         link.phi = "log", 
                         data = GasolineYield)

 summary(model_betareg)

#===========================================================
# Fit beta regression model to Gasoline data using gamlss()
#===========================================================

library(gamlss) 

model_gamlss <- gamlss(yield ~ temp, 
                       family = BE(mu.link = "logit", 
                                   sigma.link = "log"), 
                       data = GasolineYield)

summary(model_gamlss)

#===========================================================
# Predict from model_betareg for a new temp value
#===========================================================

predict_betareg <- predict(model_betareg, 
                           newdata = data.frame(temp = 349), 
                           type="response")
predict_betareg

#===========================================================
# Predict from model_gamlss for a new temp value
#===========================================================

predict_gamlss <- predict(model_gamlss, 
                          newdata = data.frame(temp = 349), 
                          type="response")
predict_gamlss

The yield value predicted by model_betareg is 0.2046289 and the yield value produced by model_gamlss is 0.204634. These are very close, as expected.

Note that the prediction is performed on the scale of the response variable yield, so that the predicted values are expressed as proportions in the interval (0,1).

The newdata must be a dataframe which lists the values of all of the predictor variables included in the right-hand side of the beta regression model. The type of these values must match the type of the predictors variables as seen by R via the str() command:

str(Gasoline)

Here's the output for this str() command:

str(Gasoline)
Classes ‘nfnGroupedData’, ‘nfGroupedData’, ‘groupedData’ and 'data.frame':      
32 obs. of  6 variables:
$ yield   : num  6.9 14.4 7.4 8.5 8 2.8 5 12.2 10 15.2 ...
$ endpoint: num  235 307 212 365 218 235 285 205 267 300 ...
$ Sample  : Ord.factor w/ 10 levels "1"<"2"<"7"<"9"<..: 8 9 5 1 3 4 7 10 6 8 
...
$ API     : num  38.4 40.3 40 31.8 40.8 41.3 38.1 50.8 32.2 38.4 ...
$ vapor   : num  6.1 4.8 6.1 0.2 3.5 1.8 1.2 8.6 5.2 6.1 ...
$ ASTM    : num  220 231 217 316 210 267 274 190 236 220 ...

If temp were listed as a factor in your dataset (e.g., "High", "Medium", "Low"), you would specify it as such in your newdata:

predict_gamlss <- predict(model_gamlss, 
                          newdata = data.frame(temp = c("High","Low")), 
                          type="response")
predict_gamlss
$\endgroup$
13
  • 2
    $\begingroup$ Thank you very much! I am not an expert in statistical analysis and maybe that's the reason I have not come across beta regression in vegetation ecological studies. I am afraid to be faced with incomprehension when I come up with beta regression models. In a blog (plantecology.syr.edu/fridley/bio793/glm.html), I read about arsine square root transformation for the dependent variable and using linear models. They cited Crawley. I tried this out, however, I also get negative values in some cases, what is nonsense, because it is not possible to have a negative cover. $\endgroup$
    – parallax
    Feb 5, 2021 at 7:52
  • 1
    $\begingroup$ Maybe you didn’t know where to look? 😜 Here’s one reference on Estimating plant abundance using inflated beta distributions: Applied learnings from a lichen–caribou ecosystem by Keim et al. that will come in handy: onlinelibrary.wiley.com/doi/pdf/10.1002/ece3.2625. $\endgroup$ Feb 5, 2021 at 12:47
  • 1
    $\begingroup$ This is an interesting article and some quite convincing points, thank you again! I'll give beta regression a try;-) $\endgroup$
    – parallax
    Feb 5, 2021 at 13:49
  • 1
    $\begingroup$ I've tried beta regression using gamlss, however, I also need to predict to a new data set but this seems not yet implemented for gamlss. With the betareg package ist is possible. You write that the inverse logit transformation would be suitable for my purpose. Do you mean to transform (y/100) and then do the modelling, or am I wrong? I think it can't be done via a link function. $\endgroup$
    – parallax
    Feb 11, 2021 at 10:59
  • 2
    $\begingroup$ @parallax: You are mistaken that gamlss can't accommodate model predictions. Please see the Addendum to my original answer for R code that proves my point. $\endgroup$ Feb 11, 2021 at 19:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.