0

我有一些分类数据显示植物家族、果实类型、果实颜色和种子散布，其中响应变量（散布）为 0 表示否或 1 表示是。

 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

我正在执行二项式 GLM，并且正在使用 MuMIn 包来“疏通”全局模型。

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)

当我尝试根据 AICc 权重绘制最重要的模型时（正如我之前在其他迭代中所做的那样），我收到警告错误：

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"

有谁知道为什么会这样？我在许多其他迭代中使用了相同的代码（即更改响应变量）并且以前没有问题。

4

2 回答 2

1

之所以会出现此问题，是因为subset(x, delta < 2)只选择了一种情况（请参阅 参考资料x$delta）。显然，该函数plot.model.selection无法处理此问题。您可以通过增加子集范围或在绘图函数中插入一些代码来修复错误来解决此问题。这是我通过添加以下行来解决的问题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)
}

您还需要通过 将该Weights功能从隐藏环境中移除Weights = MuMIn:::Weights。然后这对我有用：

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) ))))
于 2016-08-31T16:16:28.613 回答
1

这是您要发送到的对象的 str() 调用顶部的输出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

我怀疑您只是进行了一个限制性调用，它没有提供绘图函数执行任何有用工作所需的尽可能多的信息。使用subset(x, delta <5)成功。

于 2016-08-31T16:19:15.610 回答