Расхождение в результатах epiDisplay ordinal.or.display значение p и рассчитанное p для polr из пакета MASS?

Я пытаюсь рассчитать отношения шансов и связанные с ними p-значения для набора биомаркеров, используя порядковую логистическую регрессию (в частности, polr из пакета MASS). Я использовал функцию 'ordinal.or.display ()' из пакета epiDisplay для просмотра результатов регрессии, но заметил несоответствие между отображаемым p-значением и тем, что я вычисляю вручную ... это примерно вдвое больше большой, когда я рассчитываю его по нормальному распределению. Мне не хватает чего-то особенного для порядковой логистической регрессии, или это проблема с функцией epiDisplay?

Я попытался посмотреть документацию для пакета epiDisplay (https://cran.r-project.org/web/packages/epiDisplay/epiDisplay.pdf), но не нашел ничего, что объясняло бы, как там рассчитывается p-значение. Любая помощь или дополнительные знания приветствуются!

#the model using polr from MASS
#generating artifical outcome var from mt cars

mtcars <- mtcars
mtcars$outcome <- round(runif(nrow(mtcars))*5)
myMod <-polr(ordered(outcome) ~ factor(am)+ factor(carb)+ wt,
             data = mtcars,
             Hess = TRUE)

summary <- summary(myMod)
(ctable <- coef(summary(myMod)))
p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2

## combined table: p = 0.624
(ctable <- cbind(ctable, "p value" = p))


ordinal.or.display(myMod) #p value = 0.312

Я ожидаю, что значения p будут одинаковыми - возможно, epiDisplay не умножается на два?


r ordinal logistic-regression statistics mass
person svanalsten    schedule 21.06.2019    source источник

Ответы (1)


arrow_upward
0
arrow_downward

Пожалуйста, устанавливайте начальное число случайных чисел всякий раз, когда вы имитируете данные. В противном случае мы не сможем воспроизвести ваши результаты.

set.seed(123)

Вы использовали нормальное распределение для получения p-значения, когда вам следует использовать t-распределение. Итак, ваш код для вычисления p-значения должен быть:

df <- summary(myMod)$df.residual
p <- pt(abs(ctable[, "t value"]), df=df, lower.tail=FALSE) * 2

Таким же образом функция ordinal.or.display из пакета epiDisplay выполняет вычисления. Вы можете проверить это, посмотрев код примерно в строке 9.

library(epiDisplay)
ordinal.or.display

Теперь проверьте:

data(mtcars)
mtcars$outcome <- round(runif(nrow(mtcars))*5)

myMod <- polr(ordered(outcome) ~ factor(am)+ factor(carb)+ wt,
              data = mtcars,
              Hess = TRUE)

ctable <- coef(summary(myMod))

df <- summary(myMod)$df.residual
p <- pt(abs(ctable[, "t value"]), df=df, lower.tail=FALSE) * 2

options(digits=2)

(ctable <- cbind(ctable, "p value:" = p))

       Value Std. Error t value p value:
am1   -2.335       1.19  -1.962    0.064
carb2  0.045       0.92   0.049    0.961
carb3  1.840       1.40   1.319    0.202
carb4 -0.618       1.16  -0.534    0.599
carb6  1.721       1.85   0.933    0.362
carb8  2.418       2.10   1.152    0.263
wt    -0.875       0.68  -1.295    0.210

ordinal.or.display(myMod)
#Waiting for profiling to be done...
      Ordinal OR lower95ci upper95ci P value
am1   0.097      0.009     0.971     0.064  
carb2 1.046      0.166     6.477     0.961  
carb3 6.298      0.428     112.827   0.202  
carb4 0.539      0.051     5.065     0.599  
carb6 5.59       0.121     273.826   0.362  
carb8 11.218     0.162     864.225   0.263  
wt    0.417      0.106     1.573     0.210

Для меня это то же самое.

person Edward    schedule 27.06.2020