0

I have some categorical data showing plant family, fruit type, fruit colour, and seed dispersal, where the response variable (Dispersal) is either 0 for no or 1 for yes.

 test1.3
             FAMILY FRUIT_TYPE COLOUR_RF Dispersal
3   Erythroxylaceae          D         R         1
4         Lamiaceae          D         G         1
8        Clusiaceae          D         Y         1
12       Clusiaceae          D         Y         1
16        Myrtaceae          D         R         1
19        Rubiaceae          D         R         0
22    Anacardiaceae          D        Br         1
25  Melastomataceae          D         R         1
29         Moraceae          F         Y         1
32         Moraceae          F        Br         1
35         Fabaceae          C        Br         1
37        Lauraceae          D        PG         1
39      Sapindaceae          D        Br         1
41        Myrtaceae          D         R         1
43         Moraceae          D         G         1
45       Clusiaceae          D         Y         1
51         Moraceae          F         Y         1
52        Lauraceae          D         R         0
54        Rubiaceae          D         Y         0
57    Euphorbiaceae          D        PY         0
75  Dichapetalaceae          D         Y         1
83         Moraceae          F         R         1
86        Myrtaceae          D         Y         1
91       Solanaceae          D         Y         1
94      Myrsinaceae          D         Y         1
98      Connaraceae          D         O         1
101       Ochnaceae          C         R         1
104      Proteaceae          D        Br         0
107      Clusiaceae          D         R         1
114      Clusiaceae          D         P         1
116      Clusiaceae          D         P         1
119     Smilacaceae          D         R         1
124     Apocynaceae          D         Y         1
129     Apocynaceae          D        Br         1
138     Icacinaceae          D         Y         0
141        Moraceae          D         Y         1
144     Myrsinaceae          D         R         0
147  Pittosporaceae          D         O         1
150     Sapindaceae          C         O         1
154        Fabaceae          C         Y         1
158     Aphloiaceae          C         W         1
169    Celastraceae          D         O         1
176        Oleaceae          D         P         0
179       Rubiaceae          D         Y         0
182       Meliaceae          D         R         0
186     Apocynaceae          C        PY         1
188      Salicaceae          D         R         1
192     Burseraceae          D        Br         0
195         Araceae          D         Y         0
198       Rubiaceae          D         P         1
202        Rutaceae          D         O         1
206 Torrecilliaceae          D        PY         0
214       Arecaceae          D        PY         1
220        Rutaceae          D        PY         0
223       Rubiaceae          D        DR         0
230       Rubiaceae          D         B         0
237      Clusiaceae          D         Y         1
244     Myrsinaceae          D         R         1
247 Melastomataceae          D         R         0
250    Aquifoliacea          D         R         1
260       Marsaceae          D         W         1
263        Vitaceae          D        DR         0
266       Lamiaceae          D         B         0
278    Hypericaceae          D         O         1
281     Cannelaceae          D         Y         0
283       Rubiaceae          D         R         1
289      Sapotaceae          D        Br         1
293       Lauraceae          D         R         0
343     Sapindaceae          D         O         0
344       Lauraceae          D         P         0
362      Clusiaceae          D        Gr         1
366   Anacardiaceae          D        Br         1
370       Lauraceae          D         P         1
374       Lauraceae          D         R         0
399       Lauraceae          D         R         0
405       Lauraceae          D         R         0

I am performing binomial GLMs and am using the MuMIn package to 'dredge' a global model.

Global<-glm(Dispersal~FAMILY+COLOUR_RF+FRUIT_TYPE+COLOUR_RF*FRUIT_TYPE+COLOUR_RF*FAMILY, family=binomial(link="logit"), na.action="na.fail", data=test1.3)

When I try to plot the most significant model based on AICc weight (as I have done many times before with other iterations), I get the warning error:

library(MuMIn)
    x<-dredge(Global)

Fixed term is "(Intercept)"
Warning messages:
1: glm.fit: fitted probabilities numerically 0 or 1 occurred 
2: glm.fit: fitted probabilities numerically 0 or 1 occurred 
3: glm.fit: fitted probabilities numerically 0 or 1 occurred 
4: glm.fit: fitted probabilities numerically 0 or 1 occurred 
5: glm.fit: algorithm did not converge 
6: glm.fit: fitted probabilities numerically 0 or 1 occurred 
7: glm.fit: fitted probabilities numerically 0 or 1 occurred 
8: glm.fit: fitted probabilities numerically 0 or 1 occurred 
9: glm.fit: fitted probabilities numerically 0 or 1 occurred 

par(mfrow=c(1,1))
par(mar = c(5,5,14,5))
plot(subset(x, delta <2), labAsExpr = T, xlab=c("Predictor Variables for 'P. edwardsi' Dispersal"), ylab=c(expression(paste("Model Number by Cumulative Akaike Weight ", (omega) ))))

"Error in pal[, 1L + ((i - 1L)%%npal)] : incorrect number of dimensions"

Does anyone know why this is occurring? I have used the same code for many other iterations (i.e. changing the response variable) and have had no previous issues.

Ryan Rothman
  • 163
  • 1
  • 4
  • 12

2 Answers2

1

The issue is caused because subset(x, delta < 2) only selects one case (see x$delta). Apparently, the function plot.model.selection cannot handle this. You can fix this by increasing the subset range, or by inserting some code in the plot function to fix the error. Here's my fix by adding the line pal <- as.matrix(pal):

plot0 = function (x, ylab = NULL, xlab = NULL, labels = attr(x, "terms"), 
               labAsExpr = FALSE, col = c("SlateGray", "SlateGray2"), col2 = "white", 
               border = par("col"), par.lab = NULL, par.vlab = NULL, axes = TRUE, 
               ann = TRUE, ...) 
{
  if (is.null(xlab)) 
    xlab <- NA
  if (is.null(ylab)) 
    ylab <- expression("Cumulative Akaike weight" ~ ~(omega))
  op <- par(..., no.readonly = TRUE)
  on.exit(par(op))
  cumweight <- cumsum(weight <- Weights(x))
  stdweight <- weight/max(weight)
  n <- nrow(x)
  m <- length(attr(x, "terms"))
  plot.new()
  plot.window(xlim = c(0, m), ylim = c(1, 0), xaxs = "i", 
              yaxs = "i")
  pal <- if (is.na(col2)) 
    rbind(col)
  else vapply(col, function(x) grDevices::rgb((grDevices::colorRamp(c(col2, 
                                                                      x)))(stdweight), maxColorValue = 255), character(n))
  pal <- as.matrix(pal) # the only new line
  npal <- ncol(pal)

  for (i in 1L:m) rect(i - 1, c(0, cumweight), i, 
                       c(cumweight, 1), 
                       col = ifelse(is.na(x[, i]), NA, pal[, 1L + ((i - 1L)%%npal)]), 
                       border = border)
  if (ann) {
    labCommonArg <- list(col = par("col.axis"), font = par("font.axis"), 
                         cex = par("cex.axis"))
    if (labAsExpr) {
      labels <- gsub(":", "%*%", labels, perl = TRUE)
      labels <- gsub("\\B_?(\\d+)(?![\\w\\._])", "[\\1]", 
                     labels, perl = TRUE)
      labels <- parse(text = labels)
    }
    arg <- c(list(side = 3L, padj = 0.5, line = 1, las = 2), 
             labCommonArg)
    for (i in names(par.lab)) arg[i] <- par.lab[i]
    if (is.expression(labels)) {
      if (length(labels) != m) 
        stop("length of 'labels' is not equal to number of terms")
      for (i in 1L:m) do.call("mtext", c(list(text = as.expression(labels[[i]]), 
                                              at = i - 0.5), arg))
    }
    else if (!is.null(labels) && !is.na(labels)) {
      if (length(labels) != m) 
        stop("length of 'labels' is not equal to number of terms")
      do.call("mtext", c(list(text = labels, at = 1L:m - 
                                0.5), arg))
    }
    arg <- c(list(side = 4, las = 2, line = 1, adj = 1), 
             labCommonArg)
    for (i in names(par.vlab)) arg[i] <- par.vlab[i]
    ss <- weight > -(1.2 * strheight("I", cex = arg$cex))
    arg[["at"]] <- (c(0, cumweight[-n]) + cumweight)[ss]/2
    arg[["text"]] <- rownames(x)[ss]
    arg$line <- arg$line + max(strwidth(arg[["text"]], cex = arg$cex, 
                                        units = "in"))/par("mai")[4L] * par("mar")[4L]
    do.call(mtext, arg)
    title(ylab = ylab, xlab = xlab)
  }
  if (axes) {
    axis(2L, col = border, col.ticks = border)
    box(col = border)
  }
  invisible(x)
}

You will also need to take the Weights function out of the hidden environment via Weights = MuMIn:::Weights. Then this works for me:

plot0(subset(x, delta < 2), labAsExpr = T, 
      xlab=c("Predictor Variables for 'P. edwardsi' Dispersal"), 
      ylab=c(expression(paste("Model Number by Cumulative Akaike Weight ", (omega) ))))
Vandenman
  • 3,046
  • 20
  • 33
1

This is the output of the top of the str() call on the object you are sending to plot():

> xsub <- subset(x, delta <2)
> str(xsub)
Classes ‘model.selection’ and 'data.frame': 1 obs. of  11 variables:
 $ (Intercept)         : num 17.6
 $ COLOUR_RF           : Factor w/ 1 level "+": NA
 $ FAMILY              : Factor w/ 1 level "+": NA
 $ FRUIT_TYPE          : Factor w/ 1 level "+": 1
 $ COLOUR_RF:FAMILY    : Factor w/ 1 level "+": NA
 $ COLOUR_RF:FRUIT_TYPE: Factor w/ 1 level "+": NA
 $ df                  : int 3
 $ logLik              : num -44.3
 $ AICc                : num 94.8
 $ delta               : num 0
 $ weight              : num 1

I suspect you have simply made a restrictive call that doesn't provide as much information as needed for the plotting function to do any useful work. Using subset(x, delta <5) succeeds.

IRTFM
  • 258,963
  • 21
  • 364
  • 487