累积链接混合模型的概率预测 [英] Probability predictions with cumulative link mixed models

问题描述

我正在尝试使用ordinal包来拟合累积链接混合模型，但是对于获取预测概率我有些不了解.我使用ordinal包中的以下示例:

I am trying to fit cumulative link mixed models with the ordinal package but there is something I do not understand about obtaining the prediction probabilities. I use the following example from the ordinal package:

   library(ordinal)
data(soup)
## More manageable data set:
dat <- subset(soup, as.numeric(as.character(RESP)) <=  24)
dat$RESP <- dat$RESP[drop=TRUE]
m1 <- clmm2(SURENESS ~ PROD, random = RESP, data = dat, link="logistic",  Hess = TRUE,doFit=T)
summary(m1)
str(dat)

现在我正在尝试获取新数据集的概率预测

Now I am trying to get predictions of probabilities for a new dataset

newdata1=data.frame(PROD=factor(c("Ref", "Ref")), SURENESS=factor(c("6","6")))

使用

predict(m1, newdata=newdata1)

但是我遇到了以下错误

Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) : 
  contrasts can be applied only to factors with 2 or more levels

为什么会出现此错误? predict.clmm2()的语法中是否有错误?通常，predict.clmm2()会输出哪些概率? Pr(J<j)还是Pr(J=j)?有人可以指出我的信息(站点，书本)，这些材料专门用于将分类(序数)序数混合模型与R拟合.通过我在文献和网络中的搜索，大多数研究人员都将这些模型与SAS拟合.

Why am I getting this error? Is there something in the syntax of predict.clmm2() wrong? Generally which probabilities does does predict.clmm2() output? The Pr(J<j) or Pr(J=j)? Could someone point me to information (site, books) material regarding fitting categorical (ordinal) ordinal mixed models specifically with R. From my search in the literature and net, most researchers fit these kind of models with SAS.

推荐答案

您没有说出您要纠正的内容，但是当我使用它时，我没有收到任何错误提示:

You did not say what you corrected, but when I use this, I get no error:

newdata1=data.frame(PROD=factor(c("Test", "Test"), levels=levels(dat$PROD)), 
                    SURENESS=factor(c("1","1")) )
predict(m1, newdata=newdata1)

除非您将所有因子水平对齐，使它们与输入数据一致，否则带有newdata参数的predict.clmm2的输出将没有多大意义.

The output from predict.clmm2 with a newdata argument will not make much sense unless you get all the factor levels aligned so they are in the agreement with the input data:

> newdata1=data.frame(
                PROD=factor(c("Ref", "Test"), levels=levels(dat$PROD)), 
                SURENESS=factor(c("1","1")) )
> predict(m1, newdata=newdata1)
 [1] 1 1 1 1 1 1 1 1 1 1 1 1

不太有趣.该预测是针对只有一个级别的结果，处于该级别的概率为1. (虚假的预测.)但是重新创建原始有序结果的结构更有意义:

Not very interesting. The prediction is for an outcome with only one level to have a probability of 1 of being in that level. (A vacuous prediction.) But recreating the structure of the original ordered outcomes is more meaningful:

> newdata1=data.frame(
             PROD=factor(c("Ref", "Test"), levels=levels(dat$PROD)), 
             SURENESS=factor(c("1","1"), levels=levels(dat$SURENESS)) , )
> predict(m1, newdata=newdata1)
[1] 0.20336975 0.03875713

您可以通过汇总各个级别的所有预测来回答评论中的问题:

You can answer the question in the comments by assembling all the predictions for various levels:

> sapply(as.character(1:6), function(x){ newdata1=data.frame(PROD=factor(c("Ref", "Test"), levels=levels(dat$PROD)), SURENESS=factor(c(x,x), levels=levels(dat$SURENESS))  );predict(m1, newdata=newdata1)})
              1          2          3          4         5         6
[1,] 0.20336975 0.24282083 0.10997039 0.07010327 0.1553313 0.2184045
[2,] 0.03875713 0.07412618 0.05232823 0.04405965 0.1518367 0.6388921
> out <- .Last.value
> rowSums(out)
[1] 1 1

概率是Pr(J=j|X=x & Random=all).

这篇关于累积链接混合模型的概率预测的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！