Почему тесты гипотез на переделанных наборах данных слишком часто отклоняют нуль?

tl; dr: Начиная с набора данных, сгенерированного под нулевым значением, я повторно проанализировал случаи с заменой и провел проверку гипотезы для каждого повторно выбранного набора данных. Эти проверки гипотез отклоняют ноль более 5% времени.

В приведенном ниже очень простом моделировании я генерирую наборы данных с $X \sim N(0,1) \amalg Y \sim N(0,1)$ и подгоняю простую модель OLS для каждого. Затем для каждого набора данных я генерирую 1000 новых наборов данных путем повторной выборки строк исходного набора данных с заменой (алгоритм, специально описанный в классическом тексте Дэвисона и Хинкли как подходящий для линейной регрессии). Для каждого из них мне подходит одна и та же модель OLS. В конечном итоге, около 16% тестов гипотез в загрузочных выборках отклоняют ноль , тогда как мы должны получить 5% (как мы делаем в исходных наборах данных).

Я подозревал, что это как-то связано с повторными наблюдениями, вызывающими завышенные ассоциации, поэтому для сравнения я попробовал два других подхода в приведенном ниже коде (закомментировано). В методе 2 я исправляю $X$ , а затем заменяю $Y$ остатками с передискретизацией из модели OLS в исходном наборе данных. В методе 3 я рисую случайную подвыборку без замены. Обе эти альтернативы работают, то есть их тесты на гипотезы отклоняют ноль в 5% случаев.

Мой вопрос: прав ли я, что повторные наблюдения являются виновником? Если это так, учитывая, что это стандартный подход к начальной загрузке, где именно мы нарушаем стандартную теорию начальной загрузки?

Обновление № 1: Больше симуляций

$Y$

# note: simulation takes 5-10 min on my laptop; can reduce boot.reps
#  and n.sims.run if wanted
# set the number of cores: can change this to match your machine
library(doParallel)
registerDoParallel(cores=8)
boot.reps = 1000
n.sims.run = 1000

for ( j in 1:n.sims.run ) {

  # make initial dataset from which to bootstrap
  # generate under null
  d = data.frame( X1 = rnorm( n = 1000 ), Y1 = rnorm( n = 1000 ) )

  # fit OLS to original data
  mod.orig = lm( Y1 ~ X1, data = d )
  bhat = coef( mod.orig )[["X1"]]
  se = coef(summary(mod.orig))["X1",2]
  rej = coef(summary(mod.orig))["X1",4] < 0.05

  # run all bootstrap iterates
  parallel.time = system.time( {
    r = foreach( icount( boot.reps ), .combine=rbind ) %dopar% {

      # Algorithm 6.2: Resample entire cases - FAILS
      # residuals of this model are repeated, so not normal?
      ids = sample( 1:nrow(d), replace=TRUE )
      b = d[ ids, ]

      # # Method 2: Resample just the residuals themselves - WORKS
      # b = data.frame( X1 = d$X1, Y1 = sample(mod.orig$residuals, replace = TRUE) )

      # # Method 3: Subsampling without replacement - WORKS
      # ids = sample( 1:nrow(d), size = 500, replace=FALSE )
      # b = d[ ids, ]

      # save stats from bootstrap sample
      mod = lm( Y1 ~ X1, data = b ) 
      data.frame( bhat = coef( mod )[["X1"]],
                  se = coef(summary(mod))["X1",2],
                  rej = coef(summary(mod))["X1",4] < 0.05 )

    }
  } )[3]


  ###### Results for This Simulation Rep #####
  r = data.frame(r)
  names(r) = c( "bhat.bt", "se.bt", "rej.bt" )

  # return results of each bootstrap iterate
  new.rows = data.frame( bt.iterate = 1:boot.reps,
                         bhat.bt = r$bhat.bt,
                         se.bt = r$se.bt,
                         rej.bt = r$rej.bt )
  # along with results from original sample
  new.rows$bhat = bhat
  new.rows$se = se
  new.rows$rej = rej

  # add row to output file
  if ( j == 1 ) res = new.rows
  else res = rbind( res, new.rows )
  # res should have boot.reps rows per "j" in the for-loop

  # simulation rep counter
  d$sim.rep = j

}  # end loop over j simulation reps



##### Analyze results #####

# dataset with only one row per simulation
s = res[ res$bt.iterate == 1, ]

# prob of rejecting within each resample
# should be 0.05
mean(res$rej.bt); mean(s$rej)

Обновление № 2: ответ

$H_0$

При повторной выборке нулевого значения ожидаемое значение коэффициента регрессии равно нулю. При повторной выборке некоторого наблюдаемого набора данных ожидаемое значение является наблюдаемым коэффициентом для этих данных. Это не ошибка типа I, если P <= 0,05 при повторной выборке наблюдаемых данных. Фактически, это ошибка типа II, если у вас P> 0,05.

Вы можете получить некоторую интуицию, вычислив корреляцию между abs (b) и средним значением (P). Вот более простой код для копирования того, что вы сделали, плюс вычисление корреляции между ошибкой b и «типа I» в наборе симуляций

boot.reps = 1000
n.sims.run = 10
n <- 1000
b <- matrix(NA, nrow=boot.reps, ncol=n.sims.run)
p <- matrix(NA, nrow=boot.reps, ncol=n.sims.run)
for(sim_j in 1:n.sims.run){
  x <- rnorm(n)
  y <- rnorm(n)
  inc <- 1:n
  for(boot_i in 1:boot.reps){
    fit <- lm(y[inc] ~ x[inc])
    b[boot_i, sim_j] <- abs(coefficients(summary(fit))['x[inc]', 'Estimate'])
    p[boot_i, sim_j] <- coefficients(summary(fit))['x[inc]', 'Pr(>|t|)']
    inc <- sample(1:n, replace=TRUE)
  }
}
# note this is not really a type I error but whatever
type1 <- apply(p, 2, function(x) sum(x <= 0.05))/boot.reps
# correlation between b and "type I"
cor(b[1, ], type1)

обновить ответ по grand_chat не является причиной, частота P <= 0,05 составляет> 0,05. Ответ очень прост, и то, что я сказал выше - ожидаемое значение среднего значения для каждого повторного образца является исходным, наблюдаемым средним значением. Это вся основа начальной загрузки, которая была разработана для генерации стандартных ошибок / пределов достоверности для наблюдаемого среднего значения, а не для проверки гипотез. Поскольку ожидание не равно нулю, конечно, «ошибка типа I» будет больше, чем альфа. И именно поэтому будет существовать корреляция между величиной коэффициента (как далеко от нуля) и величиной отклонения «ошибки I типа» от альфа.

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

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

Основным следствием является то, что вывод, который вы делаете на основе нормальной теории, возвращает неверные значения при применении к образцу начальной загрузки. В частности, нормальная теория дает антиконсервативные правила принятия решений, потому что ваша выборка при начальной загрузке будет производить статистики, знаменатели которых меньше, чем можно было бы ожидать в обычной теории, из-за наличия дубликатов. В результате проверка гипотезы нормальной теории в конечном итоге отвергает нулевую гипотезу больше, чем ожидалось.$p$$t$

Я полностью согласен с ответом @ JWalker.

Есть еще один аспект этой проблемы. Это в процессе повторной выборки. Вы ожидаете, что коэффициент регрессии будет сосредоточен вокруг нуля, потому что вы предполагаете Xи Yнезависимы. Тем не менее, в вашей выборке вы делаете

ids = sample( 1:nrow(d), replace=TRUE )
  b = d[ ids, ]

который создает корреляцию, потому что вы делаете выборку Xи Yвместе. Например, скажем, что первая строка набора данных d: « (x1, y1)В наборе данных с передискретизацией P(Y = y1|X = x1) = 1, если Xи если Yони независимы, то P(Y|X = x1)следует нормальному распределению.

Так что еще один способ исправить это использовать

b = data.frame( X1 = rnorm( n = 1000 ), Y1 = rnorm( n = 1000 ) )

тот же код, который вы используете для генерации d, чтобы сделать X и Y независимыми друг от друга.

Эта же причина объясняет, почему он работает с остаточной передискретизацией (потому что Xне зависит от новой Y).

Самая большая проблема здесь заключается в том, что результаты модели являются ложными и, следовательно, крайне нестабильными, потому что модель просто соответствует шуму. В очень буквальном смысле. Y1 не является зависимой переменной из-за того, как были сгенерированы данные выборки.

Отредактируйте в ответ на комментарии. Позвольте мне еще раз попытаться объяснить мое мышление.

С OLS общая цель состоит в том, чтобы обнаружить и количественно определить основные взаимосвязи в данных. С реальными данными мы обычно не знаем их точно.

Но это искусственная проверка ситуации. Мы знаем механизм генерирования данных EXACT, он прямо в коде, размещенном OP.

X1 = rnorm (n = 1000), Y1 = rnorm (n = 1000)

Если мы выражаем это в знакомой форме регрессии OLS, то есть

Y1 = перехват + Beta1 * X1 + ошибка,
которая становится
Y1 = среднее (X1) + 0 (X1) + ошибка

Так что, на мой взгляд, это модель, выраженная в линейной ФОРМЕ, но на самом деле это НЕ линейная зависимость / модель, потому что нет наклона. Beta1 = 0,000000.

Когда мы создадим 1000 случайных точек данных, диаграмма рассеяния будет выглядеть как классический круговой спрей с дробовиком. Может быть некоторая корреляция между X1 и Y1 в конкретной выборке из 1000 случайных точек, которая была сгенерирована, но если это так, то это случайная случайность. Если OLS действительно находит корреляцию, то есть отвергает нулевую гипотезу об отсутствии наклона, поскольку мы точно знаем, что между этими двумя переменными действительно нет никакой связи, тогда OLS буквально нашла шаблон в Компоненте ошибок. Я охарактеризовал это как «подходящие по шуму» и «ложные».

Кроме того, одно из стандартных допущений / требований OLS заключается в том, что (+/-) «модель линейной регрессии является« линейной по параметрам ». Учитывая данные, я полагаю, что мы не выполняем это предположение. Следовательно, базовые статистические данные о значимости являются неточными. Я считаю, что нарушение предположения о линейности является прямой причиной неинтуитивных результатов начальной загрузки.

Когда я впервые прочитал эту проблему, она не утонула в том, что OP намеревался провести тестирование при нулевом [гипотезе].

Но случились бы такие же неинтуитивные результаты, если бы набор данных был сгенерирован как

X1 = rnorm (n = 1000), Y1 = X1 * .4 + rnorm (n = 1000)?