함수 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)