[R-sig-ME] predict() with cumulative link mixed models fitted with clmm2

[Rune Haubo]

On 30 November 2016 at 15:25, Miina Jari (Luke) <jari.miina at luke.fi> wrote:
> Dear list,
>
> I would like to predict() in the ordinal package using a model fitted with clmm2. The aim is to calculate the classification table of the ordered logistic regression model.
>
> Using the wine data and clm the classes can be predicted and the classification efficiency can be calculated.
>> m1 <- clm(rating ~ temp + contact, data = wine)
>> m1
> formula: rating ~ temp + contact
> data:    wine
>
> link  threshold nobs logLik AIC    niter max.grad cond.H
>  logit flexible  72   -86.49 184.98 6(0)  4.01e-12 2.7e+01
>
> Coefficients:
>   tempwarm contactyes
>      2.503      1.528
>
> Threshold coefficients:
>    1|2    2|3    3|4    4|5
> -1.344  1.251  3.467  5.006
>> predict(m1, newdata=wine, type="class")
> $fit
> [1] 2 2 3 3 3 3 4 4 2 2 3 3 3 3 4 4 2 2 3 3 3 3 4 4 2 2 3 3 3 3 4 4 2 2 3
> [36] 3 3 3 4 4 2 2 3 3 3 3 4 4 2 2 3 3 3 3 4 4 2 2 3 3 3 3 4 4 2 2 3 3 3 3
> [71] 4 4
> Levels: 1 2 3 4 5
>
> But when a model is fitted using clmm2 so that there are a nominal predictor (temp) and random effect (judge), predict() doesn't give classes but probabilities. The probabilities are predicted for the classes of rating indicated in the data, aren't they?

Yes, that is correct. The predict method for clmm2 dispatches to
predict.clm2, which does not have a type argument which explains why
it has no effect in the example below. However, if you look at
?predict.clm2 there is a section in the examples explaining how to
compute the class predictions (referencing predict.polr).  You should
be able to achieve the same with clmm2.

>> m4 <- clmm2(rating ~ contact, random = judge, nominal = ~ temp, data = wine, Hess = TRUE)
>> m4
> Cumulative Link Mixed Model fitted with the Laplace approximation
>
> Call:
> clmm2(location = rating ~ contact, nominal = ~temp, random = judge,
>     data = wine, Hess = TRUE)
>
> Random effects:
>            Var Std.Dev
> judge 1.197886 1.09448
>
> Location coefficients:
> contactyes
>   1.784829
>
> No Scale coefficients
>
> Threshold coefficients:
>                    1|2       2|3       3|4       4|5
> (Intercept)  -1.530285  1.355902  4.450757  21.96606
> tempwarm    -26.859452 -2.683005 -3.352416 -19.09974
>
> log-likelihood: -80.39439
> AIC: 180.7888
>> head(predict(m4, newdata=wine, type="class"))
> [1] 0.61714084 0.19337206 0.54060344 0.06501408 0.19620705 0.19620705
>
> Note that probabilities are predicted for an average judge by including the data used to fit the model in the newdata argument of predict.
>
> My questions:
> (1) How do I predict() the probabilities for all the classes (1-5) of rating to determine the classification efficiency?

Please see above.

> (2) It seems that the sign of the predictor will change when it appears as a nominal predictor (see coefficients for tempwarm in m1 and m4). However, according to the equations (1) and (3) in the tutorial on fitting CLMs with ordinal (https://cran.r-project.org/web/packages/ordinal/vignettes/clm_tutorial.pdf), the sign should not change. To me, the sign of the coefficients of nominal predictors and threshold parameters should be the same, or have I missed something?

That is my mistake - sorry, and thanks for pointing it out. The
nominal effects are implemented as effects on the thresholds (theta_j)
 - not the regression parameters (beta) and therefore the sign is
different from what is reflected in eq 3 for the nominal effects. The
right-hand side of eq 3 should have read theta_1j + theta_2j(contact)
- beta(temp) to be true to the implementation.

Hope this helps
Rune

>
> Regards,
> Jari Miina
>
>
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models