runehaubo / ordinal Goto Github PK

Hello,

Can you please clarify where the random effects covariance structure is specified in the source code? I would like to specify the variance-covariance structure of 6 random effects parameters as shown below. This is a modified version of the Toeplitz covariance structure that can be implemented in SAS PROC MIXED. I have spent several hours sorting through the source code but have been unsuccessful in modifying the code to permit this covariance structure. Any assistance would be greatly appreciated!

I can't seem to find in the documentation specifics on how missing values are handled in clmm(). My assumption was that the default is complete case, as it is for most functions, but my output indicates less values are being deleted than I would anticipate under CCA.

As such I am assuming clmm() is implementing FIML under the Missing at Random assumption as the default? This is a great if true given that this will produced unbiased estimates, I'm just having trouble finding documentation to support that this is indeed how the function is operating. I understand the documentation indicates that MLE is being used, but based on my understanding, MLE doesn't guarantee/imply FIML (and most functions where the outcome is categorical can't implement FIML, at least in R). I am hoping you can help clarify how these missing values are being handled within this function.

Thanks in advance!

Amanda

This happens whenver the random effects specification passed to clmm is longer than a certain limit - maybe 60 characters. The error generated is this one:

"Error in names(retrms) <- term.names :
'names' attribute [2] must be the same length as the vector [1]"

I suppose the numbers can change.

It seems to be due to this line:

term.names <- unlist(lapply(barlist, function(x) deparse(x)))

line 220 in the current version of clmm.R.
I have made a local copy of clmm.R in which this was solved by adding a width.cutoff argument:

term.names <- unlist(lapply(barlist, function(x) deparse(x,width.cutoff = 500)))

By the way, any plans to add a corStruct argument to clmm? Or any workaround you know of?

Thanks for the very useful package.

Hi,

I was wondering whether the monotonicity constraint could be enforced during MLE as a constraint. Would require switichg to an external library like nlopt(r) but might get rid of the nuicance of not being able to use a fit just because the constraint is violated.

I was comparing the results of the cumulative link model results from mgcv and ordinal package. For analogous models, every result (coefficients and AIC) seems to be almost the same, except for the fitted (or predicted) values. I found it when I was trying to calculate Somers Delta as a goodness of fit index for the models, using DescTools package, and I noted that comparing models with different loglikelihood in ordinal package, sometimes a model with low loglikelihood would have higher Somers Delta, which does not make so much sense to me. However, with mgcv package, the models with higher loglikelihood were more coherent having higher Somers Delta index. I will share here some code that I have been testing.

#load packages
library(mgcv)
library(ordinal)
library(DescTools)

#Simulate some ordered categorical data...
set.seed(3);n<-400
dat <- gamSim(1,n=n)
dat$f <- dat$f - mean(dat$f)
alpha <- c(-Inf,-1,0,5,Inf)
R <- length(alpha)-1
y <- dat$f
u <- runif(n)
u <- dat$f + log(u/(1-u))
for (i in 1:R) {
y[u > alpha[i]&u <= alpha[i+1]] <- i
}
dat$y <- y

#Fit with mgcv
b <- gam(y~x0+x1+x2+x3,family=ocat(R=R),data=dat)
summary(b)
AIC(b)
SomersDelta(x=b$fitted.values, y=dat$y, conf.level = 0.95) # 0.58

#Fit with ordinal
b2<-clm(as.factor(y)~x0+x1+x2+x3,data=dat,link="logit",threshold="flexible")
summary(b2)
AIC(b2)
SomersDelta(x=b2$fitted.values, y=dat$y, conf.level = 0.95) # -0.27

plot(x=b$fitted.values,y=dat$y)
plot(x=b2$fitted.values,y=dat$y)
plot(x=b$fitted.values,y=b2$fitted.values)

I would like to create a simple function that i) automates the process of evaluating whether particular independent variables violate the parallel regression assumption and then ii) subsequently returns a regression with the in which the independent variables that do violate the parallel regression assumption are evaluated accordingly. To that end, I've written the following code, where I use the nominal_test() function in the ordinal package to implement part i).

runOrdReg = function(dv,  vnames, dataset){
       require(ordinal)
       # run parallel regression
       parallelReg = clm(formula(paste0(dv, '~', paste(vnames, collapse = '+'))), 
                       data = dataset, link = 'logit')
       
       # implement nominal_test
       nominalTest = nominal_test(parallelReg)
      
      # select outcome and nominal variables based off of nominal test
        nomVars = rownames(nominalTest)[which(nominalTest$`Pr(>Chi)`<0.1)]      
        outcomeVars = setdiff(vnames, nomVars)

      # run non parallel regression
          nonParallelReg  <- clm(as.formula(paste0(dv, '~', paste(outcomeVars, collapse = '+'))),
                 nominal = as.formula(paste0( '~', paste(nomVars, collapse = '+'))) ,
                 data = dataset, link = 'logit')

       return(nonParallelReg)
}

However, I have found that for some reason, the nominal_test() function does not work as desired when using it inside another function. To illustrate, if you run the following code,

set.seed(2)
data = data.frame(y = factor(sample(1:4, 100, replace = TRUE)),
                  x1 = rnorm(100),
                  x2 = rnorm(100, 10, 5),
                  x3 = sample(-100:100, 100, replace = TRUE))

runOrdReg = function(dv,  vnames, dataset){
       require(ordinal)
       # run parallel regression
       parallelReg = clm(formula(paste0(dv, '~', paste(vnames, collapse = '+'))), 
                       data = dataset, link = 'logit')
}
model_func = runOrdReg(dv = 'y', vnames = paste0('x', 1:3), dataset = data)
nominal_test(model_func)

The nominal_test() function does not return any results. However, running the same model, but outside of the user-created function, nominal_test() DOES return results, as evidenced if you run the following:

model = clm(as.formula(paste0('y', '~', paste(paste0('x', 1:3), collapse = '+'))), data = data, link = 'logit')
nominal_test(model)

I'm guessing the issue lies somewhere in the source code for nominal_test(), would it be possible to look into this? Thanks very much!

I was trying to run a simple clm() until I got the above mentioned error.
I think this error is caused by the change in R 4.2 as explained in this post

I am using a multinomial DV measuring provision of social support. This variable was originally constructed as 0 = No support, 1 = One Support, 2 = Multiplex support. I factored "No Support", "One Support" and "Multiplex Support" before running the regressions with clmm(). However, I remember that in STATA one can establish a reference with baseoutcome(). How can I do this with the ordinal package?

library(lme4)
fm1 <- lmer(Reaction ~ Days + (Days | Subject), data=sleepstudy)
ranef(fm1)
library(ordinal)
ranef(fm1)

The first call to ranef works, the second returns an error

Error in UseMethod("ranef") : 
  no applicable method for 'ranef' applied to an object of class "c('lmerMod', 'merMod')"

The reason for the error is that generic ranef function, rather than being imported to the ordinal package from nlme, is defined inside the package (at line 29 of clmm.ranef.R). Can you please delete these lines and instead add to NAMESPACE importFrom(nlme, ranef) (which is currently actually at lines 42-43 of the file, but it has been commented out).

Received an error when running the following command using the {modeldata} package's ames data.

> ordinal::clm(Overall_Cond ~ Lot_Frontage + Lot_Area, data = modeldata::ames)
Error in x$code == 0L || action == "silent" : 
  'length = 2' in coercion to 'logical(1)'

Looks like the error condition noted in the significant user changes to R v 4.3.0. Indeed when debugging I found that this model receives two warning codes. Below is a snippet of what x returns from the debugging frames.

...
# $code
# [1] 2 3
# 
# $cond.H
# [1] 23671919819
# 
# $messages
# [1] "Model is nearly unidentifiable: very large eigenvalue\n - Rescale variables?" 
# [2] "Model is nearly unidentifiable: large eigenvalue ratio\n - Rescale variables?"
...

I think this error could be fixed in the print.conv.check() function if on line 291 of convergence.R
if(x$code == 0L || action == "silent") return(invisible())

was changed to

if(x$code[[1]] == 0L || action == "silent") return(invisible())

It is my understanding that this code should only execute if the first element of the x$code vector is 0L which itself only occurs when there are no other convergence warning codes.

I would be happy to make this change and submit a pull request if this would be helpful.

My R environment:

> sessionInfo()
R version 4.3.1 (2023-06-16 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19045)

Matrix products: default


locale:
[1] LC_COLLATE=English_United States.utf8  LC_CTYPE=English_United States.utf8    LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C                           LC_TIME=English_United States.utf8    

time zone: America/Chicago
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
 [1] modeldata_1.2.0     utf8_1.2.3          R6_2.5.1            numDeriv_2016.8-1.1 tidyselect_1.2.0    Matrix_1.6-1       
 [7] lattice_0.21-8      magrittr_2.0.3      glue_1.6.2          tibble_3.2.1        ordinal_2022.11-16  ucminf_1.2.0       
[13] pkgconfig_2.0.3     dplyr_1.1.2         generics_0.1.3      lifecycle_1.0.3     cli_3.6.1           fansi_1.0.4        
[19] grid_4.3.1          vctrs_0.6.3         compiler_4.3.1      rstudioapi_0.15.0   tools_4.3.1         nlme_3.1-162       
[25] pillar_1.9.0        rlang_1.1.1         MASS_7.3-60 

Hi there,
Thank you for constructing this package! I really appreciate it.
I met a problem with its application.
When I was trying to fit a clmm model, I tried those two ways to write the formula:
Method 1:

m1 <- clmm(c_score ~ modalityID + (1 + modalityID|caseID) + (1 + modalityID|readerID) +  (1|readerID : caseID),
            data=data.temp, link="logit", control = ctrl)

This resulted in the warning: no. random effects (=1082) >= no. observations (=819).
Method 2:

m2 <- clmm(c_score ~ modalityID + (1|caseID) + (1 |readerID) +  (1|readerID : caseID) + (1|modalityID:caseID) + (1|modalityID:readerID),
            data=data.temp, link="logit", control = ctrl)

This one was like: no. random effects (=1062) >= no. observations (=819)
I was just wondering why there is a 20 difference of random effects in the model?
My dataset is not fully cross, it has Number of groups: readerID:caseID 803, modalityID:caseID 123, caseID 58, modalityID:readerID 43, readerID 35, which has a sum of 1062.
Could you help me with this? Why does method 1 have 20 more random effects than method 2?
Thank you very much. Let me know if there is any problem.

Sincerely,
Leona Wang

Hi, I am using the clm function to fit an model on ordinal rating data, and I keep run into "Hessian is not positive definite", or "Hessian is numerically singular" issue for 70% of my datasets. The model is quite simple, in that there is only one fixed effect, what will you recommend to fix this problem? By the way, is there a way to make it speed up the computation when using clmm? And it would be helpful to understand why hessian is computed when there is no random effect? Thank you very much!

Thanks for creating this package.

I'm attempting to use nominal_test to test the proportional odds assumption for some models with a single continuous predictor, which attempt to model the relationship between a mock exam percentage score and the final grade achieved by a student.

Sometimes a value for Pr(>Chi) is returned but other times none of the values for the second row (labelled MockScore in this case) are all missing. I haven't yet worked out specifically when this occurs.

Please see the reproducible examples below.

library(ordinal)


# this example produces expected output

data1 <- data.frame(
  MockScore = c(
    0.5, 0.52, 0.54, 0.56, 0.58, 0.6, 0.62, 0.64, 0.66, 0.68, 0.7,
    0.72, 0.74, 0.76, 0.78, 0.8, 0.82, 0.84, 0.86, 0.88, 0.9
  ),
  FinalGrade = factor(
    c("D", "D", "D", "C", "D", "C", "D", "C", "C", "B", "B", "C",
      "B", "B", "B", "A", "B", "A", "A", "A", "A"),
    levels = c("D", "C", "B", "A")
  )
)

mod1 <- clm(FinalGrade ~ MockScore, data = data1)
nominal_test(mod1)
#> Tests of nominal effects
#> 
#> formula: FinalGrade ~ MockScore
#>           Df  logLik    AIC     LRT Pr(>Chi)
#> <none>       -10.723 29.446                 
#> MockScore  2 -10.454 32.908 0.53863   0.7639


# the nominal_test for this example has a number of missing values

data2 <- data.frame(
  MockScore = c(0.404, 0.627, 0.543, 0.46, 0.557, 0.543, 0.376, 0.18, 0.446),
  FinalGrade = factor(
    c("C", "A*", "A", "A", "B", "A", "B", "C", "A*"),
    levels = c("C", "B", "A", "A*")
  )
)

mod2 <- clm(FinalGrade ~ MockScore, data = data2)
nominal_test(mod2)
#> Tests of nominal effects
#> 
#> formula: FinalGrade ~ MockScore
#>           Df  logLik    AIC LRT Pr(>Chi)
#> <none>       -9.4604 26.921             
#> MockScore

^{Created on 2019-02-26 by the reprex package (v0.2.1)}

devtools::session_info()
#> ─ Session info ──────────────────────────────────────────────────────────
#>  setting  value                       
#>  version  R version 3.5.2 (2018-12-20)
#>  os       macOS Sierra 10.12.6        
#>  system   x86_64, darwin15.6.0        
#>  ui       X11                         
#>  language (EN)                        
#>  collate  en_GB.UTF-8                 
#>  ctype    en_GB.UTF-8                 
#>  tz       Europe/London               
#>  date     2019-02-26                  
#> 
#> ─ Packages ──────────────────────────────────────────────────────────────
#>  package     * version   date       lib source        
#>  assertthat    0.2.0     2017-04-11 [1] CRAN (R 3.5.0)
#>  backports     1.1.3     2018-12-14 [1] CRAN (R 3.5.0)
#>  callr         3.1.1     2018-12-21 [1] CRAN (R 3.5.0)
#>  cli           1.0.1     2018-09-25 [1] CRAN (R 3.5.0)
#>  crayon        1.3.4     2017-09-16 [1] CRAN (R 3.5.0)
#>  desc          1.2.0     2018-05-01 [1] CRAN (R 3.5.0)
#>  devtools      2.0.1     2018-10-26 [1] CRAN (R 3.5.1)
#>  digest        0.6.18    2018-10-10 [1] CRAN (R 3.5.0)
#>  evaluate      0.13      2019-02-12 [1] CRAN (R 3.5.2)
#>  fs            1.2.6     2018-08-23 [1] CRAN (R 3.5.0)
#>  glue          1.3.0     2018-07-17 [1] CRAN (R 3.5.0)
#>  highr         0.7       2018-06-09 [1] CRAN (R 3.5.0)
#>  htmltools     0.3.6     2017-04-28 [1] CRAN (R 3.5.0)
#>  knitr         1.21      2018-12-10 [1] CRAN (R 3.5.1)
#>  lattice       0.20-38   2018-11-04 [1] CRAN (R 3.5.2)
#>  magrittr      1.5       2014-11-22 [1] CRAN (R 3.5.0)
#>  MASS          7.3-51.1  2018-11-01 [1] CRAN (R 3.5.2)
#>  Matrix        1.2-15    2018-11-01 [1] CRAN (R 3.5.2)
#>  memoise       1.1.0     2017-04-21 [1] CRAN (R 3.5.0)
#>  numDeriv      2016.8-1  2016-08-27 [1] CRAN (R 3.5.0)
#>  ordinal     * 2018.8-25 2018-08-25 [1] CRAN (R 3.5.0)
#>  pkgbuild      1.0.2     2018-10-16 [1] CRAN (R 3.5.0)
#>  pkgload       1.0.2     2018-10-29 [1] CRAN (R 3.5.0)
#>  prettyunits   1.0.2     2015-07-13 [1] CRAN (R 3.5.0)
#>  processx      3.2.1     2018-12-05 [1] CRAN (R 3.5.0)
#>  ps            1.3.0     2018-12-21 [1] CRAN (R 3.5.0)
#>  R6            2.4.0     2019-02-14 [1] CRAN (R 3.5.2)
#>  Rcpp          1.0.0     2018-11-07 [1] CRAN (R 3.5.0)
#>  remotes       2.0.2     2018-10-30 [1] CRAN (R 3.5.0)
#>  rlang         0.3.1     2019-01-08 [1] CRAN (R 3.5.2)
#>  rmarkdown     1.11      2018-12-08 [1] CRAN (R 3.5.0)
#>  rprojroot     1.3-2     2018-01-03 [1] CRAN (R 3.5.0)
#>  sessioninfo   1.1.1     2018-11-05 [1] CRAN (R 3.5.0)
#>  stringi       1.3.1     2019-02-13 [1] CRAN (R 3.5.2)
#>  stringr       1.4.0     2019-02-10 [1] CRAN (R 3.5.2)
#>  testthat      2.0.1     2018-10-13 [1] CRAN (R 3.5.0)
#>  ucminf        1.1-4     2016-08-18 [1] CRAN (R 3.5.0)
#>  usethis       1.4.0     2018-08-14 [1] CRAN (R 3.5.0)
#>  withr         2.1.2     2018-03-15 [1] CRAN (R 3.5.0)
#>  xfun          0.5       2019-02-20 [1] CRAN (R 3.5.2)
#>  yaml          2.2.0     2018-07-25 [1] CRAN (R 3.5.0)
#> 
#> [1] /Library/Frameworks/R.framework/Versions/3.5/Resources/library

Is there a predict method in development for the clmm model, can I use the clmm2 predict as a template for now?

Cheers

Derek

I received an inquiry -- rvlenth/emmeans#141 -- about support for clmm2 models in the emmeans package. These are not currently supported, and my impression is that clmm2 is an older generation of objects so that would not be a priority. But I could be wrong.

Recent changes in emmeans make it easy to add to or override support incorporated in the emmeans package. For example, if you were to provide recover_data.clmm and emm_basis.clmm methods in your package and register them via emmeans::.emm_register(), they would be used instead of the methods that currently exist in emmeans when the ordinal package is loaded. If you do this (and I encourage it), let me know and we can work out a plan for removing that support from emmeans.

And of course, if you were to do the same for clmm2 methods, those methods would become available.

Hey,

just a brief question, is the AIC computed from a clm fit comparable to the AIC obtained from a linear model fitted with lm? Would be nice to clarify that in the docs.

Best,

Kevin

Hi Rune,

Have you considered contributing to packages like lavaan or psych to enable ordered logistic functions in the SEM framework aka GSEM? This is possible in Stata and I guess MPlus, but not in R, AFAIK. Would be magnificent. Currently, polychoric matrices can instead be loaded into SEM.

venlig hilsen
Johan

Hi everybody,

I am using "clm" function (ordinal package) to explore the associations between a numeric ordinal dependent variable (5 levels, from 1 to 5) and several covariates. In particular, I am using the partial proportional odds model specifying (in "nominal =") the covariates that can vary (those variables that violated the parallel assumption) according to the DV levels.
I find some difficulties in interpreting the output and for this I kindly ask for support.
Variables whose coefficients of effect are constant predict the highest levels of DV. On the other hand, what do variables (those variable that can vary according to DV level) predict in the "nominal" output?
Judging by the values ​​it seems that they predict each lower level of DV. For example, between levels 1/2 they predict level 1 (?), between 2/3 they predict 2 (?), etc.
In VGAM package manual (p. 216) is stated that "By default, the cumulative probabilities used are P (Y≤1), P (Y≤2),..., P (Y≤J).". Does the same apply to "ordinal" package? I couldn't find information on this in the manual.

Thanks in advance for your help.

Best
Simone

Hello Rune,

I am using the clmm function which allows for more than one random effects (nested) variables to be included. However, then I also want to use the predict function to predict probabilities - is this functionality not available with clmm? Is there a way around?

Thanks,

Using R 4.0.2 and ordinal 2019.12-10, I get a warning when using a formula with long names (but not changing the underlying data or the formula structure).

See for example this reprex:

library(ordinal)
wine1 <- wine
# short version:
mm <- clmm(rating ~ temp + contact + (1|judge),
           data = wine1) # no warning
names(wine1)[2] <- "a_very_very_very_very_very_long_name"
# long version:
mm <- clmm(a_very_very_very_very_very_long_name ~ temp + contact + (1|judge),
           data = wine1) # warning
#> Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
#>   Consider formula(paste(x, collapse = " ")) instead.
# it's the length of the whole formula that matters, not the length
# of a single variable's name:
wine2 <- wine
names(wine2) <- lapply(names(wine), paste0, "_quite_long")
mm <- clmm(rating_quite_long ~ temp_quite_long + contact_quite_long + (1|judge_quite_long),
           data = wine2) # warning too
#> Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
#>   Consider formula(paste(x, collapse = " ")) instead.

^{Created on 2020-07-23 by the reprex package (v0.3.0)}

I believe this might be because of how the clmm.formulae() function is defined, and might very well need a fix similar to issue #13 and PR #30, but I have a limited understanding of the functions involved, so I could be completely wrong! 🙂

Dear Rune,

Thank you for developing and maintaining the R package “ordinal”.

I have the following dataset:

> df
   Fish_ID Tank Diet     Rank
1      148   A1   FM   Normal
2      149   A1   FM     Mild
3      160   A2   FM   Normal
4      161   A2   FM Moderate
5      162   A2   FM Moderate
6      178   A3   FM     Mild
7      180   A3   FM     Mild
8      196   A4   FM Moderate
9      197   A4   FM Moderate
10     198   A4   FM Moderate
11     145   B1   IM   Normal
12     146   B1   IM   Normal
13     147   B1   IM   Normal
14     199   B2   IM Moderate
15     200   B2   IM   Normal
16     201   B2   IM   Normal
17     205   B3   IM Moderate
18     206   B3   IM   Normal
19     207   B3   IM   Normal
20     214   B4   IM   Normal
21     215   B4   IM     Mild
22     216   B4   IM   Normal
> str(df)
'data.frame':	22 obs. of  4 variables:
 $ Fish_ID: chr  "148" "149" "160" "161" ...
 $ Tank   : Factor w/ 8 levels "A1","A2","A3",..: 1 1 2 2 2 3 3 4 4 4 ...
 $ Diet   : Factor w/ 2 levels "FM","IM": 1 1 1 1 1 1 1 1 1 1 ...
 $ Rank   : Ord.factor w/ 3 levels "Normal"<"Mild"<..: 1 2 1 3 3 2 2 3 3 3 ...

I ran the clmm and found a huge condition number of the Hessian. It was mentioned in the tutorial that high condition numbers of the Hessian indicate that the model is ill defined, that the model can be simplified, that possibly some parameters are not identifiable, and that optimization of the model can be difficult. How can I check if my dataset meets the model assumptions of clmm? What can I do in the next steps with a huge Hessian condition number?

> mod1 <- clmm(Rank ~ Diet + (1|Tank), data = df)
> summary(mod1)
Cumulative Link Mixed Model fitted with the Laplace approximation

formula: Rank ~ Diet + (1 | Tank)
data:    df

 link  threshold nobs logLik AIC   niter    max.grad cond.H 
 logit flexible  22   -19.44 46.88 266(288) 2.59e-08 5.9e+06

Random effects:
 Groups Name        Variance Std.Dev. 
 Tank   (Intercept) 5.83e-09 7.636e-05
Number of groups:  Tank 8 

Coefficients:
       Estimate Std. Error z value Pr(>|z|)  
DietIM  -2.1272     0.9199  -2.312   0.0208 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Threshold coefficients:
              Estimate Std. Error z value
Normal|Mild    -1.1039     0.6627  -1.666
Mild|Moderate  -0.1227     0.6026  -0.204

Addtionally, when I applied the more accurate approximation to the model, the adaptive Gauss-Hermite quadrature method, I got a warning message in the model summary and many summary statistics are missing. Why is that?

> mod2 <- clmm(Rank ~ Diet + (1|Tank), data = df, Hess = TRUE, nAGQ = 10)
> summary(mod2)
Cumulative Link Mixed Model fitted with the adaptive Gauss-Hermite quadrature approximation with 10 quadrature points 

formula: Rank ~ Diet + (1 | Tank)
data:    df

 link  threshold nobs logLik AIC   niter    max.grad cond.H
 logit flexible  22   -19.44 46.88 264(286) 2.84e-08 NaN   

Random effects:
 Groups Name        Variance  Std.Dev. 
 Tank   (Intercept) 5.431e-09 7.369e-05
Number of groups:  Tank 8 

Coefficients:
       Estimate Std. Error z value Pr(>|z|)
DietIM   -2.127         NA      NA       NA

Threshold coefficients:
              Estimate Std. Error z value
Normal|Mild    -1.1039         NA      NA
Mild|Moderate  -0.1227         NA      NA
Warning message:
In summary.clmm(mod2) :
  Variance-covariance matrix of the parameters is not defined

I'm trying to test the significance of a single (categorical) fixed effect, and I ran into this small problem when I speicified a null model with no fixed effects:

> mm1 <- clmm(SURENESS ~ PROD + (1|RESP) + (1|RESP:PROD), data = soup,
+             link = "probit", threshold = "equidistant")
> 
> mm0 <- clmm(SURENESS ~ (1|RESP) + (1|RESP:PROD), data = soup,
+             link = "probit", threshold = "equidistant")
Error in environment(fixedForm) <- environment(form) <- form.envir : 
  cannot set attribute on a symbol
> 
> mm0 <- clmm(SURENESS ~ 1 + (1|RESP) + (1|RESP:PROD), data = soup,
+             link = "probit", threshold = "equidistant")

Making the intercept explicit worked, but I thought I should let you know about this edge case.

Note that this also affects drop1:

> drop1(mm1, ~PROD)
Error in environment(fixedForm) <- environment(form) <- form.envir : 
  cannot set attribute on a symbol

Dear @runehaubo,

Thank you for the very useful package ordinal!

I’m using the GitHub version 2020.8-22 (not CRAN version 2019.12-10 because of issue #34) in R 4.2.1 on Windows.

My question concerns the following warning that I get when running clmm() models: length(x) = 3 > 1' in coercion to 'logical(1). This might have to do with this issue.

The syntax of my model looks like this:

model = clmm(X  ~  Y1 + Y2 + Y3 + Y4 + Y1:Y2 + Y1:Y4 + (1+Y1 | subject), data = dataset, link = "logit")

I get the same warning for models with only one interaction term with Y1.

The warning only seems to appear for models with random slopes. Also, the warning disappears when I remove the 1+ part in (1+Y1 | subject). When fitting a model with two by-subject slopes, the warning appears again (also when excluding the 1+ part).

Could you tell me whether I’m doing something wrong, or whether this is a bug?

Thank you very much!

Hi,

Thank you for developing and maintaining the ordinal package.

I've run clmm for my dataset but how can I check the model assumption? Nominal_test works for clm object but not for clmm object. Any suggestions on how to do it in R?

Hi,
thank you very much for your package. Thanks to its ease-of-use I convinced colleagues working on ordinal variables to go from classical LMs to CLMs!
This is not really to raise an issue but rather asking for a new functionality. Would you consider implementing a model.matrix() method for clmm objects, as it already exists for clm objects?
Thanks!
Maxime

The dgumbel and pgumbel functions are only symmetric around 0. With a non-zero location parameter, the distribution should reflect around the location parameter.

> dgumbel(10, 10, 1, max=TRUE)
[1] 0.3678794
> dgumbel(10, 10, 1, max=FALSE)
[1] 0

Also these should match:

> pgumbel(10, 10, 1, lower.tail=TRUE, max=FALSE)
[1] 1
> pgumbel(10, 10, 1, lower.tail=FALSE, max=TRUE)
[1] 0.6321206

qgumbel, rgumbel give correct results; however, you may want to check ggumbel as well.

Hi @runehaubo,

I believe this issue is similar to the #34 which relates to a warning regarding formula length. However, this is occurring for the ordinal::clm() function instead. Another user has already raised this in #34 but the issue has since been closed.

Due to this warning, it's causing me issues with rendering tables using another package.

Can this function be fixed as well?

Thanks.

See reproducible example below:

# load packages
library(ordinal)
library(dplyr)

# data set
df <- mtcars %>% 
  mutate_all(.funs = ~ factor(.)) %>% 
  mutate(very_very_very_long_variable1 = cyl,
         very_very_very_long_variable2 = vs,
         very_very_very_long_variable3 = gear)

# this produces a warning
clm(formula = very_very_very_long_variable1 ~ very_very_very_long_variable2 + very_very_very_long_variable3,
    data = df)
#> Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
#>   Consider formula(paste(x, collapse = " ")) instead.

# this does not result in a warning
clm(formula = cyl ~ vs + gear,
    data = df)
#> formula: cyl ~ vs + gear
#> data:    df
#> 
#>  link  threshold nobs logLik AIC   niter max.grad cond.H 
#>  logit flexible  32   -15.00 39.99 5(0)  4.85e-07 2.8e+01
#> 
#> Coefficients:
#>    vs1  gear4  gear5 
#> -4.611 -3.180 -3.476 
#> 
#> Threshold coefficients:
#>    4|6    6|8 
#> -5.852 -2.891

^{Created on 2021-08-17 by the reprex package (v2.0.1)}

Hi @runehaubo, I have used your ordinal package intensively lately and created the following functions meanwhile, so that one clm() model can accept multiple different threshold types across the intercept and variables in the nominal effects with possibly user-specified constraints.

With these alternations, clm() can function like VGAM::vglm() where any coefficient can be constrain. I prefer the former as clm() also allows scale effects. My custom functions enable a much more parsimonious model specification and avoid numerical problems (i.e., extremely large coefficient) often associated with relaxing the parallel line assumption. All these functions work very well in my model fitting. As my functions inevitably alter many important ordinal package routines, I had to update several interval functions, such as predict.clm() which contains a small mistake in interval probability mass estimate.

I also expanded the functionality of post-fitting functions such as those for scale- and nominal-effect test. Now they correctly code interaction terms and show likelihood ratio test results when either adding (new) or removing (existing) terms from an expanded scope that may or may not include the variables already in the model. This appears useful and necessary, as one variable may have a coefficient estimate at 0 and thus is not retained in a fitted model; the reason could be that its nominal effects are -b, 0, and b across three thresholds, making the location-effect-only coefficient nonsignificant but some nominal-effect coefficients significant. Therefore, testing terms not in the fitted model is beneficial. I also made these test functions indicate warning messages of numerical and convergence problems. The nominal test function also allows the user to specify the threshold type for the added nominal-effect terms.

I can contribute by sharing my source codes if needed.

Related to emmeans support in Issue 22; if clmm2 is older and/or outdated, will support for scale or nominal arguments be supported in clmm?

After I ran the following model, the software-R crashed.
I am trying to run the model that includes one random effect, I got this issue. Variables names are perfectly coded.

m4<-clmm(Satis ~ rePS + reTL+ reMXL + (1|MJA), data=data1, link = "logit", threshold = "equidistant")

Thanks in advance

The model.matrix() method does not work correctly for clmm objects:

> wine.clmm <- clmm(rating ~ temp * contact + (1|judge), data = wine)
> model.matrix(wine.clmm)
Error in eval(predvars, data, env) : object 'rating' not found

This is important to me (i.e., the emmeans package) because if one happens to have a rank deficiency, we need the model matrix to extract the basis for non-estimable functions of the regression coefficients.

Closely related... In such rank-deficient cases, coef(<clm model>) returns NAs for dropped columns, whereas these are not included with coef(<clmm model>). That's inconsistent and requires special handling on my part.

Hi,

I am fitting a random effects model of the form clmm(OrdinalVar ~ 1 + (1|Group)). The model fits correctly, but every time I try to set the nAGQ argument to any number, I get the error in the subject.

I also tried using clmm2(OrdinalVar ~ 1, random = Group) and got an error: Error in eval(substitute(expr), data, enclos = parent.frame()) : NA/NaN/Inf in foreign function call (arg 6).

This is the call stack for clmm:

3: stop(gettextf("non-finite likelihood at starting value (%g)", 
       init.val), call. = FALSE)
2: clmm.fit.ssr(rho, control = control$optCtrl, method = control$method, 
       Hess)
1: clmm(AGradeC ~ 1 + (1 | ALangC), data = Datd, nAGQ = 7)

This is the call stack from clmm2:

10: .C("NRalgv3", as.integer(ctrl$trace), as.integer(ctrl$maxIter), 
        as.double(ctrl$gradTol), as.integer(ctrl$maxLineIter), as.integer(grFac), 
        as.double(stDev), as.double(o1), as.double(o2), as.double(eta1Fix), 
        as.double(eta2Fix), as.double(sigma), as.integer(linkInt), 
        as.double(weights), u = as.double(uStart), pr = as.double(pr), 
        funValue = double(1), gradValues = as.double(uStart), hessValues = as.double(rep(1, 
            length(uStart))), length(pr), length(uStart), maxGrad = double(1), 
        conv = 0L, as.double(lambda), Niter = as.integer(Niter))
9: eval(substitute(expr), data, enclos = parent.frame())
8: eval(substitute(expr), data, enclos = parent.frame())
7: with.default(rho, .C("NRalgv3", as.integer(ctrl$trace), as.integer(ctrl$maxIter), 
       as.double(ctrl$gradTol), as.integer(ctrl$maxLineIter), as.integer(grFac), 
       as.double(stDev), as.double(o1), as.double(o2), as.double(eta1Fix), 
       as.double(eta2Fix), as.double(sigma), as.integer(linkInt), 
       as.double(weights), u = as.double(uStart), pr = as.double(pr), 
       funValue = double(1), gradValues = as.double(uStart), hessValues = as.double(rep(1, 
           length(uStart))), length(pr), length(uStart), maxGrad = double(1), 
       conv = 0L, as.double(lambda), Niter = as.integer(Niter))[c("u", 
       "funValue", "gradValues", "hessValues", "maxGrad", "conv", 
       "Niter")])
6: with(rho, .C("NRalgv3", as.integer(ctrl$trace), as.integer(ctrl$maxIter), 
       as.double(ctrl$gradTol), as.integer(ctrl$maxLineIter), as.integer(grFac), 
       as.double(stDev), as.double(o1), as.double(o2), as.double(eta1Fix), 
       as.double(eta2Fix), as.double(sigma), as.integer(linkInt), 
       as.double(weights), u = as.double(uStart), pr = as.double(pr), 
       funValue = double(1), gradValues = as.double(uStart), hessValues = as.double(rep(1, 
           length(uStart))), length(pr), length(uStart), maxGrad = double(1), 
       conv = 0L, as.double(lambda), Niter = as.integer(Niter))[c("u", 
       "funValue", "gradValues", "hessValues", "maxGrad", "conv", 
       "Niter")])
5: rho$updateU(rho)
4: ObjFun(rhoM, x)
3: fn(.x, ...)
2: ucminf(rhoM$start, function(x) ObjFun(rhoM, x), control = rhoM$optCtrl)
1: clmm2(AGradeC ~ 1, random = ALangC, data = Datd, link = "logistic", 
       Hess = T, nAGQ = 7)

The variable ALangC is an unordered factor with 20 levels, but I get an error even if I use a character vector instead.

Thank you for any help,

k.

Inside a function, anova.clmm() doesn’t work properly. I think it looks for the objects in the wrong environment (the global one), so it uses the wrong objects (if similarly named global objects exist) or exists with an error (if they don’t). Here’s a simple reprex:

First I simulate five datasets:

library(ordinal)

# Simulate a simple data sett
simulate = function(n = 100) {
  y = ordered(sample(1:5, n, replace = TRUE))
  x = runif(n)
  gr = sample(1:10, n, replace =TRUE)
  data.frame(y, x, gr)
}

# Simulate 5 datasets
set.seed(123)
l_data = replicate(5, simulate(), simplify = FALSE)

Analysing one of the datasets works fine:

# Analyse the first dataset manually
# and extract the p value of the x coefficient
d = l_data[[1]]
mod0 = clmm(y ~ 1 + (1|gr), data = d)
mod1 = update(mod0, . ~ x + (1|gr))
p = anova(mod0, mod1)$`Pr(>Chisq)`[2]
p # OK
#> [1] 0.2188232

But I want to easily analyse all of them, so I write a function:

get_p_val = function(d) {
  mod0 = clmm(y ~ 1 + (1|gr), data = d)
  mod1 = update(mod0, . ~ x + (1|gr))
  p = anova(mod0, mod1)$`Pr(>Chisq)`[2]
  p
}

Calling the function seemingly works fine for the first dataset. But in fact, it uses the model objects from the global environment instead of the ones fitted inside the function. So we get the same P-value for all datasets:

get_p_val(l_data[[1]]) # *Looks* OK (but isn’t)
#> [1] 0.2188232
get_p_val(l_data[[2]]) # Wrong P-value (same as above)
#> [1] 0.2188232
get_p_val(l_data[[3]]) # Still the wrong P-value + strange warning
#> Warning in update.uC(rho): Non finite negative log-likelihood
#>   at iteration 26
#> [1] 0.2188232

Note that for l_data[[3]], there is an additional warning (which I guess is from fitting one of the models for the third data set).

I would typically use the get_p_val() function like this:

lapply(l_data, get_p_val) # Identical P-values
#> Warning in update.uC(rho): Non finite negative log-likelihood
#>   at iteration 26
#> [[1]]
#> [1] 0.2188232
#> 
#> [[2]]
#> [1] 0.2188232
#> 
#> [[3]]
#> [1] 0.2188232
#> 
#> [[4]]
#> [1] 0.2188232
#> 
#> [[5]]
#> [1] 0.2188232

It is easy to see that the function uses the global variables, since we get an error message if we remove one of them:

rm(mod1)
get_p_val(l_data[[3]])
#> Warning in update.uC(rho): Non finite negative log-likelihood
#>   at iteration 26
#> Error in (function (object, ..., type = c("I", "II", "III", "1", "2", : object 'mod1' not found

Also note that the anova() function for clm models does not have this problem. You can easily test this by just replacing clmm with clm in the code above.

I am interested in understanding the relationship between SurveyGroup(levels=1,2) and food availability (ordinal variable: None < Inconstant < Constant). I've run the following model with confounding variables are predictors:

mod<- clm(food_availability~ SurveyGroup + Age + education + gender + PropLifeFeeding + 
observation + centrality_mean + temp_seasonal, 
data = subset(econmix2, econmix2$season=="spring"), Hess = TRUE)

But then when I use nominal_test() and scale_test(), SurveyGroup does not meet the assumptions. So, I specify it in the nominal argument, but then I don't have a coefficient estimated and it is instead in the thresholds section of the summary output.

mod<- clm(food_availability~ gender + PropLifeFeeding + 
observation + centrality_mean + temp_seasonal, 
nominal = ~ SurveyGroup + Age + education,
data = subset(econmix2, econmix2$season=="spring"), Hess = TRUE)

Coefficients:
                Estimate Std. Error z value Pr(>|z|)    
genderMale      0.015329   0.149843   0.102 0.918520    
PropLifeFeeding 0.188023   0.311829   0.603 0.546530    
observation     0.036880   0.010640   3.466 0.000528 ***
centrality_mean 0.209999   0.052579   3.994  6.5e-05 ***
temp_seasonal   0.000742   0.007174   0.103 0.917617    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Threshold coefficients:
                                  Estimate Std. Error z value
None|Inconstant.(Intercept)      -3.321715   0.680552  -4.881
Inconstant|Constant.(Intercept)   0.627151   0.482889   1.299
None|Inconstant.SurveyGroup2      1.074639   0.340853   3.153
Inconstant|Constant.SurveyGroup2  0.096961   0.169354   0.573

How do I assess the relationships between SurveyGroup and food availability then?

Hi @runehaubo, I used your ordinal package intensively lately. I noticed that ordinal:::prob.predict.clm() that generate predicted probability masses may have incorrect components, as its confidence interval is built from only a logistic distribution around a linear predictor while clm() allows a few other distributions. Given the link function, this should be an easy fix.

I've generated a very simple multilevel ordinal model using HLM software and using the clmm function from the ordinal package in R. The threshold values and coefficients are comparable, however the signs for the coefficients (and only the coefficients; not the thresholds) are opposite (e.g., positive in HLM, negative in clmm). This then impacts the odds ratios for the predictors. HLM reports odds greater than one and clmm reports odds less than one. I anticipate there is a very easy reason for what's happening, and it's just escaping me. Can someone please advise on what I've overlooked?

The following example produces warnings due to variable name incompleteness:

library(ordinal)

options(warnPartialMatchDollar=TRUE) # warns when partial matching occurs

mm1 <- clmm(SURENESS ~ PROD + (1|RESP) + (1|RESP:PROD), data = soup,
                 link = "probit", threshold = "equidistant")

warnings()

mm2 <- clmm(SURENESS ~ PROD + (1|RESP), data = soup,
            link = "probit", threshold = "equidistant")

warnings()

Session info:

R version 4.1.2 (2021-11-01)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Big Sur 11.7.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] ordinal_2022.11-16 colorout_1.2-2    

loaded via a namespace (and not attached):
[1] MASS_7.3-55         compiler_4.1.2      Matrix_1.5-1       
[4] tools_4.1.2         nlme_3.1-155        ucminf_1.1-4.1     
[7] grid_4.1.2          numDeriv_2016.8-1.1 lattice_0.20-45  

Hello,

I ran a clmm with two main factors and one random factor which is "participants".
One of the main factors has 6 levels that can be either categorical or numerical.
Since it is my first time running this type of model, I do not understand why I'm getting this warning.

Do you know how could I resolve this problem?

Thank you in advance for your help and support.

I think we could probably do with testing whether coef outputs fall within confint outputs to test for multiple comparison, aka whether beta_1 (the coefficient of x_1) is significantly different from beta_2 (the coefficient of x_2). But such a mechanism probably already exists somewhere that supports clm objects, no?
I checked multcomp package and some pairwise.* functions, but none can directly deal with the clm output. Maybe there needs to be some S3 method implemented for clm before it can be supported by those other packages?

Hi,
I am using the package ordinal in R to run ordinal logistic regression mixed models on ordinal data with 4 categories. I have to use clmm2 because one of my explanatory variables does not satisfy the proportional odds assumption and the nominal effects are not available in clmm. However, I do not know how to extract the confidence interval for the odds ratio of the fixed effects with clmm2.

For clmm, I understand that confidence intervals of the odds ratio are obtained by transforming the limits of confidence intervals of the parameter estimates, provided by the confint method as shown below using the wine dataset:

require(ordinal)
fm1 <- clmm(rating ~ contact + temp + (1|judge), data = wine)
round(exp(confint(fm1)), 2)

However, for clmm2, the confint method only provides confidence intervals for the standard deviation of the random effects, and not for the fixed effects (here contact and temp, with the proportional odds assumption relaxed for contact in our example).

fm2 <- clmm2(rating ~ temp, nominal = ~ contact, random = judge, data = wine, Hess = TRUE)
pr2 <- profile(fm2)
confint(pr2)

Do you know any other method to obtain the confidence intervals for the fixed effects with clmm2? I have not been able to find the answer to this question anywhere.

Dear author,

I have another several questions regarding to clmm/clm functions. Now I am applying mainly clmm to my datasets, fitting a mixed model. and I keep encountering error messages like "Hessian is not positive definite" or "Hessian is numerically singular", and errors related to eigen values such as "Error in eigen(H, symmetric = TRUE, only.values = TRUE) : infinite or missing values in 'x'".

My model is
model <- clmm(ordered(VALUE) ~ fixed + (1|random), data = data, threshold="flexible")

One of my datasets run into the last error as mentioned above, where fixed has 8 levels, random has 25 levels, response has 200 observations in total. Histogram of response VALUE is as follows:

However, some of the errors could be solved by changing the link function from logit to probit. So I am wondering what is the preference between those two link functions when fitting a mixed model using clmm? What will you suggest when model is having hessian matrix not positive definite errors? And also why clmm does not treat effects with less or equal to 2 levels as random effect? Thank you very much for your help!

Hi,

Could you please tell me how to set two random variables in clmm2?
Like:
clmm2(Y ~ X1 + X2, random = c(X3,X4) , data = data_new, Hess=TRUE, link='logit')

Thanks!

Hi,

I want to add a matrix of 1000*1000 covariance matrix as a random effect? Could you please provide any suggestions on how could I do it? Is it okay to convert to the list and add? otherwise, how could I add whole matrix as a random effect?

Thanks

Hi,

I ran the nominal test for multiple clm objects using the map() function from the purrr package but it returned null values. The nominal test can be run for the individual dataset repeatedly, returning relevant statistics. I open an issue under the github repository of purrr package. The maintainer of the purrr pakcage suggested that this could be a bug in the ordinal pakcage which assumes the data in the global environment.

PS: I also tried for loop or lapply() function, both returned null values as I got from the map().

Regards,
Yanxian

Hi Rune,

I've got a somewhat complicated dataset with repeated measures within 8 collections. I think, for the most part, I am beginning to understand how to form my model -- but I'm very confused as to how I show clmm2 that I have repeated measures. I'm thinking the correlation of repeated measures within each collection may differ from each other a bit, so I need to try an unstructured variance-covariance structure; I was thinking compound symmetry. Can clmm2 handle this? And if so, how?? I.e., what does the syntax look like?

Any help is much appreciated -- thank you so much for maintaining this package!

Kindly,
Charity

Hi! I really need help with this. I run the following model in R:

 clmm_br<-clmm(Grado_amenaza~Life_Form  + size_max_cm + leaf_length_mean  + petals_length_mean + 
 silicua_length_mean + bloom_length + categ_color+ (1|Genero) , data=brasic1)

I didn't get any warning or error but when I run summary(clmm_br) I can't get the p-values:

  summary(clmm_br)

 Cumulative Link Mixed Model fitted with the Laplace approximation

 formula: Grado_amenaza ~ Life_Form + size_max_cm + leaf_length_mean +  
    petals_length_mean + silicua_length_mean + bloom_length +  
    categ_color + (1 | Genero)
 data:    brasic1

 link  threshold nobs logLik AIC    niter      max.grad cond.H
 logit flexible  76   -64.18 160.36 1807(1458) 1.50e-03 NaN   

 Random effects:

  Groups Name        Variance       Std.Dev. 
  Genero (Intercept) 0.000000008505 0.00009222

 Number of groups:  Genero 39 

 Coefficients:
                      Estimate Std. Error z value Pr(>|z|)
Life_Form[T.G]        2.233338         NA      NA       NA
Life_Form[T.Hem]      0.577112         NA      NA       NA
Life_Form[T.Hyd]    -22.632916         NA      NA       NA
Life_Form[T.Th]      -1.227512         NA      NA       NA
size_max_cm           0.006442         NA      NA       NA
leaf_length_mean      0.008491         NA      NA       NA
petals_length_mean    0.091623         NA      NA       NA
silicua_length_mean  -0.036001         NA      NA       NA
bloom_length         -0.844697         NA      NA       NA
categ_color[T.2]     -2.420793         NA      NA       NA
categ_color[T.3]      1.268585         NA      NA       NA
categ_color[T.4]      1.049953         NA      NA       NA

Threshold coefficients:
    Estimate Std. Error z value
1|3   -1.171         NA      NA
3|4    1.266         NA      NA
4|5    1.800         NA      NA
(4 observations deleted due to missingness)

I tried with no random effects and excluding the rows with NAs but it's the same

The structure of my data here:

str(brasic1)
tibble[,13] [80 x 13] (S3: tbl_df/tbl/data.frame)
 $ ID                 : num [1:80] 135 137 142 145 287 295 585 593 646 656 ...
 $ Genero             : chr [1:80] "Alyssum" "Alyssum" "Alyssum" "Alyssum" ...
 $ Cons.stat          : chr [1:80] "LC" "VU" "VU" "LC" ...
 $ Amenazada          : num [1:80] 0 1 1 0 1 0 0 1 0 0 ...
 $ Grado_amenaza      : Factor w/ 5 levels "1","3","4","5",..: 1 2 2 1 4 1 1 2 1 1 ...
 $ Life_Form          : chr [1:80] "Th" "Hem" "Hem" "Th" ...
 $ size_max_cm        : num [1:80] 12 6 7 15 20 27 60 62 50 60 ...
 $ leaf_length_mean   : num [1:80] 7.5 7 11 14.5 31.5 45 90 65 65 39 ...
 $ petals_length_mean : num [1:80] 2.2 3.5 5.5 2.55 6 8 10.5 9.5 9.5 2.9 ...
 $ silicua_length_mean: num [1:80] 3.5 4 3.5 4 25 47.5 37.5 41.5 17.5 2.9 ...
 $ X2n                : num [1:80] 32 NA 16 16 NA NA 20 20 18 14 ...
 $ bloom_length       : num [1:80] 2 1 2 2 2 2 2 2 11 2 ...
 $ categ_color        : chr [1:80] "1" "4" "4" "4" ...

Thank you so much for your help:)

As per https://orcid.org/0000-0002-4494-3399 this line:

https://github.com/runehaubo/ordinal/blob/f38646c44d149b91aa043ce24550800c2c92002b/inst/CITATION#L8C42-L8C60

should be:

comment = c(ORCID = "0000-0002-4494-3399")),

Hi,

I'm trying to fit a model with 6 random slopes using the clmm-function, but I keep getting the same error message.

Removal one random slope solves this error, but is there a way to include all random slopes?

Thanks very much for developing and maintaining this package. I am finding it very useful for modelling enteric infections and duration of illness.

Is there a way in clm to specify comparisons for the response levels in a nominal model? Specifically, I would like to compare all levels of the response to a reference level, as 0|1, 0|2, 0|3 etc.

Reading through the vignettes, I didn't find any examples of this.

Thanks!!

When there is a nominal predictor, confint() produces an error. It seems related to the fact that in such case, coef() returns a term with NA.

mod <- ordinal::clm(rating ~ temp * contact, data = ordinal::wine)
coef(mod)
#>                 1|2                 2|3                 3|4                 4|5 
#>          -1.4112620           1.1435537           3.3770825           4.9419823 
#>            tempwarm          contactyes tempwarm:contactyes 
#>           2.3211843           1.3474604           0.3595489
confint(mod)
#>                           2.5 %   97.5 %
#> tempwarm             0.99435182 3.761793
#> contactyes           0.08378091 2.694828
#> tempwarm:contactyes -1.45985126 2.180286

mod <- ordinal::clm(rating ~ temp * contact, data = ordinal::wine, nominal = ~contact)
coef(mod)
#>     1|2.(Intercept)     2|3.(Intercept)     3|4.(Intercept)     4|5.(Intercept) 
#>          -1.3784288           1.1187987           3.3521730           4.4521888 
#>      1|2.contactyes      2|3.contactyes      3|4.contactyes      4|5.contactyes 
#>          -1.5171328          -1.3164389          -1.2997110          -0.6150591 
#>            tempwarm          contactyes tempwarm:contactyes 
#>           2.2619528                  NA           0.5388898
confint(mod)
#> Error: Cannot profile model where standard errors are NA

^{Created on 2020-11-02 by the reprex package (v0.3.0)}