Неправильные p-значения, скорректированные по Холму, в psych :: corr.test ()?

Question

При подготовке учебных материалов по скорректированным по Холму значениям p я решил внести коррективы как вручную, так и с помощью R 's psych :: corr.test (),Все 26 от руки и отрегулированные значения R p соответствуют, кроме двух.Я предполагаю, что это ошибка пользователя, но я не могу понять, что я делаю неправильно.

Код ниже будет:

1. Воспроизвести p-значенияпод вопросом (психологический пакет).
2. Ранжирование нескорректированных значений p (самое низкое = 1).
3. Вычислить p-значения Холма как p (unadj.) * Rank.
4. Сравнить скорректированные значения p.

library(psych)  
dat = state.x77  
R.out = corr.test(dat)  
R.out$p  
unadj.p = ifelse(lower.tri(R.out$p) == F, NA, R.out$p)  
p.ranks = 29 - rank(unadj.p, na.last = T)  
p.ranks = matrix(ifelse(p.ranks < 1, NA, p.ranks), 8, 8)  
myHolm = unadj.p * p.ranks  
myHolm = ifelse(myHolm > 1, 1, myHolm)  
myHolm = t(myHolm)  
round(myHolm, 4)  
round(R.out$p, 4)  
myHolm == R.out$p  

Результаты анализов приведены ниже.Скорректированные значения p в первой (моей) таблице соответствуют значениям во второй таблице - из corr.test () - за исключением двух в столбцах строки 1 7: 8.

Это мой первый пост здесь- извините за любые отклонения от нормы публикации.

>   round(myHolm,4)  
     [,1] [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]  
[1,]   NA    1 1.0000 1.0000 0.2328 1.0000 0.2398 0.8765  
[2,]   NA   NA 0.0286 0.2343 1.0000 0.0000 1.0000 0.1616  
[3,]   NA   NA     NA 0.0002 0.0000 0.0000 0.0000 1.0000  
[4,]   NA   NA     NA     NA 0.0000 0.0002 0.7918 1.0000  
[5,]   NA   NA     NA     NA     NA 0.0065 0.0011 1.0000  
[6,]   NA   NA     NA     NA     NA     NA 0.1583 0.2510  
[7,]   NA   NA     NA     NA     NA     NA     NA 1.0000  
[8,]   NA   NA     NA     NA     NA     NA     NA     NA  
>   round(R.out$p,4)  
           Population Income Illiteracy Life Exp Murder HS Grad  Frost   Area  
Population     0.0000 1.0000     1.0000   1.0000 0.2328  1.0000 0.2510 1.0000  
Income         0.1467 0.0000     0.0286   0.2343 1.0000  0.0000 1.0000 0.1616  
Illiteracy     0.4569 0.0015     0.0000   0.0002 0.0000  0.0000 0.0000 1.0000  
Life Exp       0.6387 0.0156     0.0000   0.0000 0.0000  0.0002 0.7918 1.0000  
Murder         0.0146 0.1080     0.0000   0.0000 0.0000  0.0065 0.0011 1.0000  
HS Grad        0.4962 0.0000     0.0000   0.0000 0.0003  0.0000 0.1583 0.2510  
Frost          0.0184 0.1141     0.0000   0.0660 0.0001  0.0088 0.0000 1.0000  
Area           0.8765 0.0095     0.5938   0.4581 0.1106  0.0179 0.6828 0.0000  
>   myHolm == R.out$p  
           Population Income Illiteracy Life Exp Murder HS Grad Frost  Area  
Population         NA   TRUE       TRUE     TRUE   TRUE    TRUE FALSE FALSE  
Income             NA     NA       TRUE     TRUE   TRUE    TRUE  TRUE  TRUE  
Illiteracy         NA     NA         NA     TRUE   TRUE    TRUE  TRUE  TRUE  
Life Exp           NA     NA         NA       NA   TRUE    TRUE  TRUE  TRUE  
Murder             NA     NA         NA       NA     NA    TRUE  TRUE  TRUE  
HS Grad            NA     NA         NA       NA     NA      NA  TRUE  TRUE  
Frost              NA     NA         NA       NA     NA      NA    NA  TRUE  
Area               NA     NA         NA       NA     NA      NA    NA    NA  

Answer 1

Добро пожаловать на SO, ваш первый пост выглядит великолепно, поэтому проблем нет.

С точки зрения ваших расчетов вы столкнулись с проблемой, потому что при умножении unadj.p * p.ranks порядок ваших рангов не сохраняется. Например, если вы проверяете значение столбца 7 строки 1, или [1,7], (у вас 0,2398), оно ниже значения столбца 8 строки 7 или [7,8] (0,2510). Этого не должно быть, поскольку ваша матрица p.ranks показывает, что они должны быть на 13 и 14 ранге соответственно (при этом [1,7] является более высоким из двух).

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

    library(psych)  
    dat = state.x77  
    R.out = corr.test(dat)  
    R.out$p  
    unadj.p = ifelse(lower.tri(R.out$p)==F,NA,R.out$p)  

    # convert into vector for ease of calculation
    p <- as.numeric(unadj.p)
    # remove missing values
    p <- p[!is.na(p)]

    # find the ranks of p
    pr <- rank(p)

    # put p in order
    po <- p[order(p)]

    # put ranks in order (1 is smallest)
    pro <- pr[order(pr, decreasing = T)]

    # now they are in order we can take the CUMULATIVE MAX to preserve order
    pcum <- cummax(po * pro)

    # now put back in our order and stick in our matrix
    myHolm <-unadj.p
    myHolm[!is.na(myHolm)] <- pcum[pr]


    myHolm = ifelse(myHolm>1,1,myHolm)  
    myHolm = t(myHolm)  

    myHolm == R.out$p  

#            Population Income Illiteracy Life Exp Murder HS Grad Frost Area
# Population         NA   TRUE       TRUE     TRUE   TRUE    TRUE  TRUE TRUE
# Income             NA     NA       TRUE     TRUE   TRUE    TRUE  TRUE TRUE
# Illiteracy         NA     NA         NA     TRUE   TRUE    TRUE  TRUE TRUE
# Life Exp           NA     NA         NA       NA   TRUE    TRUE  TRUE TRUE
# Murder             NA     NA         NA       NA     NA    TRUE  TRUE TRUE
# HS Grad            NA     NA         NA       NA     NA      NA  TRUE TRUE
# Frost              NA     NA         NA       NA     NA      NA    NA TRUE
# Area               NA     NA         NA       NA     NA      NA    NA   NA