Вероятностные прогнозы с комбинированными смешанными моделями

Question

Я пытаюсь установить смешанные модели совокупной связи с пакетом 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() неправильно? Как правило, какие вероятности дает прогноз. Clmm2() output? Pr(J<j) или Pr(J=j)? Может ли кто-нибудь указать мне на информацию (сайт, книги) материал, касающийся установки категориальных (порядковых) ординарных смешанных моделей, в частности с R. Из моего поиска в литературе и в сети большинство исследователей подходят к этим типам моделей с SAS.

Answer 1

Вы не сказали, что вы исправить, но когда я использую это, я не получаю сообщение об ошибке:

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

Выход из predict.clmm2 с NewData аргумент не имеет особого смысла, если вы не получите все коэффициенты выравниваются так, что они находятся в согласии со своими входными данными:

> 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 быть на этом уровне. (A бессодержательным предсказание.) Но воссоздать структуру исходных упорядоченными исходами является более значимым:

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