函數lipsitz.test {generalhoslem}不適用於對象clm {ordinal}

Question

我現在正在嘗試使用lipsitz.test {generalhoslem}來測試ordianl模型的擬合度。根據document，函數可以處理polr和clm。但是，當我嘗試在lipsitz.test函數中使用clm時，會發生錯誤。這裏是一個例子

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. 

有沒有解決方案來解決這個問題？非常感謝。

Answer 1

我不確定這是不是主題，應該在CrossValidated上。這部分是測試代碼的問題，部分是測試本身的統計數據。

有兩個問題。我在使用clm時發現了代碼中的一個錯誤，並且會將修復推向CRAN（下面的更正代碼）。

然而，示例數據似乎確實存在一個更基本的問題。基本上，Lipsitz測試需要使用分組的虛擬變量來擬合新模型。當用這個例子擬合新模型時，模型失敗，因此一些係數不計算。如果使用polr，則新模型會收到警告，說明它是等級不足的;如果使用clm，則新模型會得到一條消息，即由於奇點而導致兩個係數未擬合。我認爲這個示例數據集不適合這種分析。

已更正的代碼在下方，我使用了一個更大的示例數據集，在該數據集上運行測試。

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)