Function lipsitz.test {generalhoslem} is not working for object clm{ordinal}

Question

I am now tring to test the goodness of fit of an ordianl model using lipsitz.test {generalhoslem}. According to the document, the function can deal with both polr and clm. However, when I try to use clm in the lipsitz.testfunction, an error occurs. Here is an example

library("ordinal")
library(generalhoslem)
data("wine")
fm1 <- clm(rating ~ temp * contact, data = wine)
lipsitz.test(fm1)

Error in names(LRstat) <- "LR statistic" : 
'names' attribute [1] must be the same length as the vector [0]
In addition: Warning message:
In lipsitz.test(fm1) :
n/5c < 6. Running this test when n/5c < 6 is not recommended.

Is there any solution to fix this? Thanks a lot.

Answer 1

I'm not sure if this is off-topic and should be on CrossValidated. It's partly a problem with the coding of the test and partly about the statistics of the test itself.

There are two problems. I've just spotted a bug in the code when using clm and will push a fix to CRAN (corrected code below).

There does however appear to be a more fundamental problem with the example data. Basically, the Lipsitz test requires fitting a new model with dummy variables of the groupings. When fitting the new model with this example, the model fails and so some of the coefficients are not calculated. If using polr, the new model gets the warning that it is rank-deficient; if using clm, the new model gets a message that two coefficients are not fitted due to singularities. I think this example data set is just unsuitable for this kind of analysis.

The corrected code is below and I have used a larger example dataset on which the test runs.

lipsitz.test <- function (model, g = NULL)  {
  oldmodel <- model
  if (class(oldmodel) == "polr") {
    yhat <- as.data.frame(fitted(oldmodel))
  } else if (class(oldmodel) == "clm") {
    predprob <- oldmodel$model[, 2:ncol(oldmodel$model)]
    yhat <- predict(oldmodel, newdata = predprob, type = "prob")$fit
  } else warning("Model is not of class polr or clm. Test may fail.")
  formula <- formula(oldmodel$terms)
  DNAME <- paste("formula: ", deparse(formula))
  METHOD <- "Lipsitz goodness of fit test for ordinal response models"
  obs <- oldmodel$model[1]
  if (is.null(g)) {
    g <- round(nrow(obs)/(5 * ncol(yhat)))
    if (g < 6) 
      warning("n/5c < 6. Running this test when n/5c < 6 is not recommended.")
  }
  qq <- unique(quantile(1 - yhat[, 1], probs = seq(0, 1, 1/g)))
  cutyhats <- cut(1 - yhat[, 1], breaks = qq, include.lowest = TRUE)
  dfobs <- data.frame(obs, cutyhats)
  dfobsmelt <- melt(dfobs, id.vars = 2)
  observed <- cast(dfobsmelt, cutyhats ~ value, length)
  if (g != nrow(observed)) {
    warning(paste("Not possible to compute", g, "rows. There might be too few observations."))
  }
  oldmodel$model <- cbind(oldmodel$model, cutyhats = dfobs$cutyhats)
  oldmodel$model$grp <- as.factor(vapply(oldmodel$model$cutyhats, 
                                         function(x) which(observed[, 1] == x), 1))
  newmodel <- update(oldmodel, . ~ . + grp, data = oldmodel$model)
  if (class(oldmodel) == "polr") {
    LRstat <- oldmodel$deviance - newmodel$deviance
  } else if (class(oldmodel) == "clm") {
    LRstat <- abs(-2 * (newmodel$logLik - oldmodel$logLik))
  }
  PARAMETER <- g - 1
  PVAL <- 1 - pchisq(LRstat, PARAMETER)
  names(LRstat) <- "LR statistic"
  names(PARAMETER) <- "df"
  structure(list(statistic = LRstat, parameter = PARAMETER, 
                 p.value = PVAL, method = METHOD, data.name = DNAME, newmoddata = oldmodel$model, 
                 predictedprobs = yhat), class = "htest")
}

library(foreign)
dt <- read.dta("http://www.ats.ucla.edu/stat/data/hsbdemo.dta")
fm3 <- clm(ses ~ female + read + write, data = dt)
lipsitz.test(fm3)
fm4 <- polr(ses ~ female + read + write, data = dt)
lipsitz.test(fm4)