Plot a generalized linear model

The plotGLM function, now included in the modEvA package (Barbosa et al. 2014), plots the observed (presence/absence) data and the predicted (probability) values of a generalized linear model against the y regression equation (logit) values. For now, it works only for logistic regression (binomial response, logit link). Optionally (if plot.values = TRUE) it can display in the plot the proportion of deviance explained (Weisberg 1980, Guisan & Zimmermann 2000), but this requires that the Dsquared function is loaded, and that a model object is provided (rather than just the observed and predicted values) if you want the adjusted deviance as well.

plotGLM <-
function(obs = NULL, pred = NULL, model = NULL, link = "logit", 
                    plot.values = FALSE, xlab = "Logit (y)", 
                    ylab = "Predicted probability", main = "Model plot", ...) {
  # version 1.7 (28 Oct 2014)
  # obs: presence/absence or other binary (1-0) observed data
  # pred: values predicted by a GLM of the binary observed data
  # model: instead of (and overriding) obs & pred, you can provide a model object of class "glm"
  # link: link function of the GLM; only 'logit' is implemented
  # plot.values: logical, whether to report the values of deviance and adjusted deviance in the plot (requires the 'model' argument and the 'Dsquared' function)
  # ...: arguments to pass to the 'plot' function
 
  model.provided <- ifelse(is.null(model), FALSE, TRUE)
 
  if (model.provided) {
    if (!is.null(pred)) message("Argument 'pred' ignored in favour of 'model'.")
    if (!is.null(obs)) message("Argument 'obs' ignored in favour of 'model'.")
    obs <- model$y
    pred <- model$fitted.values
  }  else { # if model not provided
    if (is.null(obs) | is.null(pred)) stop("You must provide either 'obs' and 'pred', or a 'model' object of class 'glm'")
  }  # end if model
 
  stopifnot(
    length(obs) == length(pred),
    obs %in% c(0, 1),
    pred >= 0,
    pred <= 1)
 
  logit <- log(pred / (1 - pred))
  if (!model.provided) model <- glm(obs ~ logit, family = "binomial")
 
  plot(pred ~ logit, ylim = c(0, 1), type = "n", xlab = xlab, ylab = ylab, 
       main = main, ...)
  points(obs ~ logit, pch = 1, col = "darkgrey")
  abline(v = 0, lty = 5, col = "grey")  # y of P = 0.5
  #abline(h = pred[which.min(abs(logit))], col = "lightgrey", lty = 2)
  points(pred ~ logit, pch = 20, cex = 0.6)
 
  if (plot.values) {
    Dsq <- Dsquared(model = model, adjust = FALSE)
    if (max(logit) > abs(min(logit)))  x.loc <- c(max(logit), 1)
    else x.loc <- c(min(logit), 0)
    text(x = x.loc[1], y = 0.6, adj = x.loc[2], labels = substitute(paste(D^2 == a), list(a = round(Dsq, 3))))
    if (model.provided) {  # adjDsq needs n parameters in original model, not just our model created from obs~logit
      adjDsq <- Dsquared(model = model, adjust = TRUE)
      text(x = x.loc[1], y = 0.4, adj = x.loc[2], labels = substitute(paste('D'['adj']^2) == a, list(a = round(adjDsq, 3))))  # http://stackoverflow.com/questions/10156417/subscripts-in-plots-in-r, plus a mail on substitute(paste) in Silwod-R
 
    }
  }  # end if plot.val
}

Created by Pretty R at inside-R.org

[presented with Pretty R]

For example, if you have a table called mydata with your species’ presence/absence in column 3 and your predicted data in column 20, paste in R the plotGLM function and then this command:

plotGLM(obs = mydata[ , 3], pred = mydata[ , 20])

Instead of obs and pred, you can provide a model object:

# mymodel <- glm(myspecies ~ var1 + var2 ……)

plotGLM(model = mymodel, plot.values = TRUE)

Rplot03

Grey circles are the presences (0) and absences (1); black dots are the corresponding predicted values. The vertical line at y = 0 corresponds to a predicted value of 0.5 (inflection point). You can add a legend and also alternative or additional titles and labels to the plot, e.g with arguments legend, ylab, xlab and main (see ?plot).

References

Barbosa A.M., Brown J.A. & Real R. (2014) modEvA – an R package for model evaluation and analysis. R package, version 0.1.

Guisan, A. & Zimmermann, N.E. (2000). Predictive habitat distribution models in ecologyEcological Modelling 135: 147-186

Weisberg, S. (1980). Applied Linear Regression. Wiley, New York

Advertisements

Comment

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s