【问题标题】:Probability predictions with model averaged Cumulative Link Mixed Models fitted with clmm in ordinal package在序数包中安装 clmm 的模型平均累积链接混合模型的概率预测
【发布时间】:2017-06-28 06:08:57
【问题描述】:

我发现predict 函数目前未在使用ordinal R 包中的clmm 函数拟合的累积链接混合模型中实现。虽然predict 是在同一个包中为clmm2 实现的,但我选择应用clmm,因为后者允许多个随机效果。此外,我还安装了几个clmm 模型,并使用MuMIn 包中的model.avg 函数执行模型平均。理想情况下,我想使用平均模型来预测概率。然而,虽然MuMIn 支持clmm 模型,但predict 也不适用于普通模型。

有没有办法破解predict 函数,以便该函数不仅可以预测来自clmm 模型的概率,还可以使用来自clmm 的模型平均系数进行预测(即“平均”类的对象) ?例如:

require(ordinal)
require(MuMIn)

mm1 <- clmm(SURENESS ~ PROD + (1|RESP) + (1|RESP:PROD), data = soup,
        link = "probit", threshold = "equidistant")

## test random effect:
mm2 <- clmm(SURENESS ~ PROD + (1|RESP) + (1|RESP:PROD), data = soup,
        link = "logistic", threshold = "equidistant")

#create a model selection object
mm.sel<-model.sel(mm1,mm2)

##perform a model average
mm.avg<-model.avg(mm.sel)


#create new data and predict
new.data<-soup

##predict with indivindual model
predict(mm1, new.data)

我收到以下错误消息: 在 UseMethod("predict") 中: 没有适用于 predict 的方法应用于“clmm”类的对象

 ##predict with model average
 predict(mm.avg, new.data)

返回另一个错误: predict.averaging(mm.avg, new.data) 中的错误: predict 模型“mm1”和“mm2”导致错误

【问题讨论】:

  • 为什么这个问题不直接问软件包作者?这似乎极有可能“过于宽泛”,因为它需要理论和实施的努力才能以有原则的方式进行。

标签: r predict mixed-models ordinal


【解决方案1】:

我找到了一个潜在的解决方案(粘贴在下面),但无法处理我的数据。

这里的解决方案:https://gist.github.com/mainambui/c803aaf857e54a5c9089ea05f91473bc

我认为问题在于我使用的系数数量,但经验不足,无法弄清楚。希望这可以帮助某人。

这是我正在使用的模型和新数据,尽管它实际上是模型平均版本。但相同的预测因子。

ma10 <- clmm(Location3 ~ Sex * Grass3 + Sex * Forb3 + (1|Tag_ID), data = 
IP_all_dunes)
ma_1 <- model.avg(ma10, ma8, ma5)##top 3 models
new_ma<- data.frame(Sex = c("m","f","m","f","m","f","m","f"),
                Grass3 = c("1","1","1","1","0","0","0","0"),
                Forb3 = c("0","0","1","1","0","0","1","1"))


# Arguments:
#  - model = a clmm model
#  - modelAvg = a clmm model average (object of class averaging)
#  - newdata = a dataframe of new data to apply the model to
# Returns a dataframe of predicted probabilities for each row and response level
fake.predict.clmm <- function(modelAvg, newdata) {
  # Actual prediction function
  pred <- function(eta, theta, cat = 1:(length(theta) + 1), inv.link = plogis) {
    Theta <- c(-1000, theta, 1000)
    sapply(cat, function(j) inv.link(Theta[j + 1] - eta) - inv.link(Theta[j] - 
eta))
  }

  # Multiply each row by the coefficients
  #coefs <- c(model$beta, unlist(model$ST))##turn off if a model average is used
  beta <- modelAvg$coefficients[2,3:12]
  coefs <- c(beta, unlist(modelAvg$ST))

  xbetas <- sweep(newdata, MARGIN=2, coefs, `*`)

  # Make predictions
  Theta<-modelAvg$coefficients[2,1:2]
  #pred.mat <- data.frame(pred(eta=rowSums(xbetas), theta=model$Theta))
  pred.mat <- data.frame(pred(eta=rowSums(xbetas), theta=Theta))
  #colnames(pred.mat) <- levels(model$model[,1])
  a<-attr(modelAvg, "modelList")
  colnames(pred.mat) <- levels(a[[1]]$model[,1])

  pred.mat
}

【讨论】:

    【解决方案2】:

    我也一直在使用 clmm,是的,我确认 predict.clmm 没有(还没有?)实现。我还没有检查 fake.predict.clmm 的源代码。它可能会起作用。如果没有,您将无法手动操作或使用 predict.clmm2。

    【讨论】:

    • 正如目前所写,您的答案尚不清楚。请edit 添加其他详细信息,以帮助其他人了解这如何解决所提出的问题。你可以找到更多关于如何写好答案的信息in the help center
    猜你喜欢
    • 1970-01-01
    • 1970-01-01
    • 2013-06-11
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2021-05-31
    • 1970-01-01
    相关资源
    最近更新 更多