Оценка силы теста нормальности (в R)

9

Я хочу оценить точность тестов нормальности для разных размеров выборки в R (я понимаю, что тесты нормальности могут вводить в заблуждение ). Например, чтобы посмотреть на тест Шапиро-Уилка, я провожу следующую симуляцию (а также нанесение на график результатов) и ожидаю, что с увеличением размера выборки вероятность отклонения нуля уменьшается:

n <- 1000
pvalue_mat <- matrix(NA, ncol = 1, nrow = n)

for(i in 10:n){
    x1 <- rnorm(i, mean = 0, sd = 1)
    pvalue_mat[i,] <- shapiro.test(x1)$p.value
}   

plot(pvalue_mat)

Я бы подумал, что по мере роста размера выборки частота отбраковки должна быть ниже, однако она кажется довольно равномерной. Я думаю, что я неправильно понимаю это - любые мысли приветствуются.

user94759
источник
2
Возможно, вы захотите взглянуть на: stats.stackexchange.com/questions/2492/…
nico

Ответы:

7

Вы моделируете по нулевой гипотезе (нормальное распределение), поэтому уровень отклонения будет стремиться к уровню значимости, как и ожидалось. Чтобы оценить мощность, вам нужно смоделировать при любом ненормальном распределении. Есть бесконечные возможности / сценарии (например, гамма-распределения с увеличением асимметрии, t-распределение с уменьшением df и т. Д.) На выбор, в зависимости от объема вашего исследования.

Майкл М
источник
Спасибо за ответ. Когда я симулирую ненормальные распределения, я наблюдаю выпуклую картину относительно начала координат - то есть, когда размер выборки становится больше для любого ненормального распределения, вероятность отклонения нулевого значения нормальности увеличивается. Тем не менее, я не понимаю, почему это не наоборот при рисовании из нормального распределения: почему вероятность отклонения нулевого значения не уменьшается при увеличении размера выборки? Спасибо
user94759
3
Потому что вероятность совершения такой ошибки типа 1 по определению равна уровню значимости, который является постоянным. Или, иначе говоря, значения p равномерно распределены под нулем. Кстати, рекомендуется выполнять много симуляций для каждого параметра, включая выбор n, а не только одного, как в вашем коде.
Майкл М
7

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


αα

Второй критерий требует, чтобы мы указали, каким образом и насколько ноль не соответствует действительности. В случае с учебниками это легко, потому что альтернативы ограничены по объему и четко определены. С тестами распределения, такими как Shapiro-Wilk, альтернатива гораздо более расплывчата: они «ненормальные». Таким образом, при выборе среди тестов распределения аналитик, вероятно, должен будет провести свое собственное одноразовое исследование мощности, чтобы оценить, насколько хорошо тесты работают в отношении более конкретных альтернативных гипотез, которые имеют отношение к рассматриваемой проблеме.

ν1ν

αR

  • rdistимя функции для получения случайной выборки из некоторого распределения

  • nРазмер образцов на заказ rdist

  • n.iterколичество таких образцов для получения

  • ...rdistν

Остальные параметры управляют отображением результатов; они включены в основном для удобства создания цифр в этом ответе.

sim <- function(rdist, n, n.iter, prefix="",
                breaks=seq(0, 1, length.out=20), alpha=0.05,
                plot=TRUE, ...) {

  # The simulated P-values.
  # NB: The optional arguments "..." are passed to `rdist` to specify
  #     its parameters (if any).
  x <- apply(matrix(rdist(n*n.iter, ...), ncol=n.iter), 2, 
             function(y) shapiro.test(y)$p.value)

  # The histogram of P-values, if requested.
  if (plot) {
    power <- mean(x <= alpha)
    round.n <- 1+ceiling(log(1 + n.iter * power * (1-power), base=10) / 2)
    hist(x[x <= max(breaks)], xlab=paste("P value (n=", n, ")", sep=""), 
         breaks=breaks, 
         main=paste(prefix, "(power=", format(power, digits=round.n), ")", sep=""))
    # Specially color the "significant" part of the histogram
    hist(x[x <= alpha], breaks=breaks, col="#e0404080", add=TRUE)
  }

  # Return the array of P-values for any further processing.
  return(x)
}

5100.n520.

n.iter <- 10^5                 # Number of samples to generate
n.spec <- c(5, 10, 20)         # Sample sizes to study
par(mfrow=c(1,length(n.spec))) # Organize subsequent plots into a tableau
system.time(
  invisible(sapply(n.spec, function(n) sim(rnorm, n, n.iter, prefix="DF = Inf ")))
)

После указания параметров этот код также занимает всего одну строку. Это дает следующий вывод:

Гистограммы для нуля

01α=0.05,.04810.0499

10.2

νν=100ν=11001000), что не занимает много времени. Теперь в коде требуется двойной цикл (а в более сложных ситуациях нам часто нужны тройные или четырехкратные циклы, чтобы учесть все аспекты, которые нам нужно варьировать): один для изучения того, как мощность изменяется в зависимости от размера выборки, и другой для изучения того, как она меняется в зависимости от степени свободы. Тем не менее, еще раз, все делается в одной строке кода (третья и последняя):

df.spec <- c(64, 16, 4, 2, 1)
par(mfrow=c(length(n.spec), length(df.spec)))
for (n in n.spec) 
  for (df in df.spec)
    tmp <- sim(rt, n, n.iter, prefix=paste("DF =", df, ""), df=df)

Гистограммы для альтернатив

Небольшое изучение этой таблицы дает хорошее представление о силе. Я хотел бы обратить внимание на его наиболее важные и полезные аспекты:

  • ν=64ν=1

  • n=5n=20

  • 1

  • 201110086.57=13%205%95%

  • αα=0.10

    α0.05020%α=0.01α=0.005breakssim


Забавно, что так много можно почерпнуть из того, что фактически составляет три строки кода: одну для имитации выборок iid из указанного дистрибутива, одну, чтобы применить это к массиву нулевых распределений, и третью, чтобы применить его к массив альтернативных распределений. Это три шага, которые входят в любой анализ мощности: остальное - это просто суммирование и интерпретация результатов.

Whuber
источник
1

(Больше, чем комментарий, возможно, не полный ответ)

[I] ожидал бы, что с увеличением размера выборки вероятность отклонения нуля уменьшается

Оставляя в стороне соображения о предвзятых тестах (которые не являются чем-то необычным с точки зрения пригодности, поэтому стоит упомянуть), есть три ситуации, связанные с уровнем отклонения, которые можно рассмотреть:

1) коэффициент отклонения при моделировании от нуля (как вы, кажется, делаете в своем вопросе)

α

2) коэффициент отклонения при моделировании из некоторой альтернативы

Здесь показатель отклонения должен увеличиваться с увеличением n .

3) коэффициент отклонения для некоторого сбора реальных данных

Практически, ноль никогда не бывает истинным, и реальные данные будут иметь некоторую смесь величин ненормальности (как измерено статистикой теста). Если степень ненормальности не связана с размером выборки, коэффициент отклонения должен увеличиваться с увеличением n .

Так что на самом деле ни в одной из этих ситуаций мы не должны видеть снижение количества браков с размером выборки.

Glen_b - Восстановить Монику
источник