1

我正在尝试将累积链接混合模型与ordinal包相匹配，但对于获得预测概率，我有一些不明白的地方。我使用ordinal包中的以下示例：

   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)

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

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 拟合。

4

1 回答 1

3

你没有说你更正了什么，但是当我使用它时，我没有得到任何错误：

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

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

> 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。（一个空洞的预测。）但是重新创建原始有序结果的结构更有意义：

> 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

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

> 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)。

于 2013-07-05T16:03:57.263 回答