0

lipsitz.test{generalhoslem}を使用して順序モデルの適合度をテストしようとしています。ドキュメントによると、関数はpolrとclmの両方を処理できます。しかし、関数clm内で使用しようとするlipsitz.testと、エラーが発生します。ここに例があります

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.

これを修正する解決策はありますか？どうもありがとう。

4

1 に答える 1

1

これがトピックから外れているかどうかはわかりませんが、CrossValidated にする必要があります。部分的にはテストのコーディングの問題であり、部分的にはテスト自体の統計に関する問題です。

2 つの問題があります。使用時にコードのバグを見つけたのでclm、修正を CRAN にプッシュします (以下の修正済みコード)。

ただし、サンプル データには、より根本的な問題があるようです。基本的に、リプシッツ検定では、グループ化のダミー変数を使用して新しいモデルを適合させる必要があります。この例で新しいモデルを適合させると、モデルが失敗するため、一部の係数が計算されません。を使用するpolrと、新しいモデルはランク不足であるという警告を受け取ります。を使用するclmと、新しいモデルは、特異点のために 2 つの係数が適合しないというメッセージを受け取ります。このサンプル データ セットは、この種の分析にはまったく適していないと思います。

修正されたコードを以下に示します。テストを実行する大きなサンプル データセットを使用しました。

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)
于 2016-07-31T20:15:00.833 に答える