R/plotGLM.R

Defines functions plotGLM

Documented in plotGLM

plotGLM <-
function(model = NULL, obs = NULL, pred = NULL, link = "logit",
         plot.values = TRUE, plot.digits = 3, xlab = "Logit (Y)",
         ylab = "Predicted probability", main = "Model plot", ...) {
  # version 1.9 (13 Apr 2016)

  model.provided <- ifelse(is.null(model), FALSE, TRUE)

  if (model.provided) {
    if(!("glm" %in% class(model) && model$family$family == "binomial" && model$family$link == "logit")) stop ("'model' must be an object of class 'glm' with 'binomial' family and 'logit' link.")
    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
    )
  
  if (any(pred < 0) | any(pred > 1)) warning("Some of your 'pred' values are outside the [0,1] interval; are you sure these represent probabilities? Unexpected or incorrect results may arise. Consider properly rescaling you 'pred' values, or obtaining real probabilities instead.")
  
  pred[pred == 0] <- 2e-16  # avoid log 0 below
  pred[pred == 1] <- 1 - 2e-16  # avoid division by 0 below
  
  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 <- round(Dsquared(model = model, adjust = FALSE), plot.digits)
    Rsq <- RsqGLM(model = model)
    CoxSnell <- round(Rsq$CoxSnell, plot.digits)
    McFadden <- round(Rsq$McFadden, plot.digits)
    Nagelkerke <- round(Rsq$Nagelkerke, plot.digits)
    Tjur <- round(Rsq$Tjur, plot.digits)

    #minX <- min(log(mod$fitted/(1 - mod$fitted)))
    
    #if (max(logit) > abs(min(logit)))  x.loc <- c(max(logit), 1)
    #else x.loc <- c(min(logit), 0)
    
    if (max(logit) > abs(min(logit)))  x.loc <- max(logit)
    else x.loc <- min(logit)

    text(x.loc, 0.95, substitute(paste(D^2 == a), list(a = Dsq)), adj = 0)
    text(x.loc, 0.8, substitute(paste(R[Cox-Snell]^2 == a), list(a = CoxSnell)), adj = 0)
    text(x.loc, 0.6, substitute(paste(R[McFadden]^2 == a), list(a = McFadden)), adj = 0)
    text(x.loc, 0.4, substitute(paste(R[Nagelkerke]^2 == a), list(a = Nagelkerke)), adj = 0)
    text(x.loc, 0.2, substitute(paste(R[Tjur]^2 == a), list(a = Tjur)), adj = 0)
    
    #text(x = x.loc[1], y = 0.6, adj = x.loc[2], labels = substitute(paste(D^2 == a), list(a = round(Dsq, plot.digits))))
    #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, plot.digits))))
    #}  # end if model provided
    
  }  # end if plot values
}

Try the modEvA package in your browser

Any scripts or data that you put into this service are public.

modEvA documentation built on May 29, 2017, 10:47 a.m.