By ECII


2013-07-05 14:43:01 8 Comments

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

with

predict(m1, newdata=newdata1)

but I am getting the following error

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

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.

1 comments

@42- 2013-07-05 16:03:57

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)

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

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

The probabilities are Pr(J=j|X=x & Random=all).

@ECII 2013-07-05 17:11:18

Thanks. I guess I missed the fact that its important to spacify matcing values and labels for the categorical regressors. Is this specific to predict.clmm2() ? Do you also happen to know what kind of probabilities are in the output of predict.clmm2()? Are they Pr(J<j) or Pr(J=j)?

@42- 2013-07-05 17:42:45

Not just the regressors, but also the outcomes.

@ECII 2013-07-06 07:37:40

Thanks a lot. Just to check, the fits are of the model log(odds)=a+bx right? I am asking because other programs tend to fit log(odds)=a-bx

@42- 2013-07-08 21:38:36

You should study the vignette: cran.r-project.org/web/packages/ordinal/vignettes/clm_intro.‌​pdf . It looks like the clm package uses follows the convention you attribute to "other programs". I think it lets the probabilities sum to unity.

Related Questions

Sponsored Content

1 Answered Questions

[SOLVED] Function predict() for object of segmented.lme()

1 Answered Questions

[SOLVED] Factor/level error in mixed model

  • 2016-12-17 22:47:59
  • elduderino260
  • 89 View
  • 0 Score
  • 1 Answer
  • Tags:   r glm mixed-models

1 Answered Questions

[SOLVED] R: Plotting predictions of MASS polr ordinal model

1 Answered Questions

[SOLVED] Adaboost in R: Predict for data that does not have dependent variable

  • 2015-10-26 12:16:22
  • Karthaveeryarjun Vinjamoori
  • 711 View
  • 2 Score
  • 1 Answer
  • Tags:   r boost predict adaboost

1 Answered Questions

Predicting probabilities in a logistic regression model

  • 2015-01-14 10:11:11
  • Tapiwa
  • 300 View
  • 0 Score
  • 1 Answer
  • Tags:   r

Sponsored Content