Как вычислить значения P для увеличения размера выборки с помощью цикла?

Хотя я понимаю, что вы можете использовать циклы for, это хорошая возможность использовать sapply или lapply. Я продемонстрирую альтернативы, используя iris. Хотя я собираюсь использовать упрощенный тест «не равно 5» для iris$Sepal.Length для всех образцов, вам следует обновить alternative= и другие аргументы для ваших конкретных данных.

Вариант 1: если все, что вам когда-либо понадобится, это значение p, мы можем зафиксировать только это ... или мы можем захватить всю модель и выполнить поиск p-значений на втором этапе.

Вариант 2: мы можем использовать одну из функций *apply, которая хорошо читается (как только вы привыкнете к векторному коду R), или вы можете придерживаться цикла for. У первого варианта есть преимущества в читаемости, хотя вам может быть удобнее использовать цикл for, и в этом случае вам действительно следует предварительно выделить список/вектор. (Причина для предварительного определения длинного, но пустого списка/вектора: хотя вы можете легко объединить вектор out с out <- c(out, newstuff), повторение этого в долгосрочной перспективе будет очень неэффективно. Я настоятельно не рекомендую делать это в «крупном масштабе».)

Впереди некоторые примечания:

• Я использую set.seed(2) для каждого, чтобы результаты были идентичными. Вы не должны использовать его до тех пор, пока вам не потребуется строгая воспроизводимость. Обычно не требуется для производственных/академических отчетов.
• Я делаю seq_len вместо 2:length(...) из-за привычек: когда делаешь что-то программно, хорошо, когда он изящно терпит неудачу. Если по какой-то причине в будущем вы используете 1:length(nrow(x)), а x оказывается, что у него 0 строк, то 1:0 создает вектор длины 2, что нелогично (и почти наверняка нарушит последующий код). Вместо этого seq_len(0) создает вектор длины 2, что хорошо. Опять же, здесь менее критично, но это хорошая привычка. (Кстати: seq_along(0) по-прежнему выводит вектор длины 1, поэтому он тоже подвержен этой проблеме.)
• Я использовал seq_len(...)[-1], чтобы отбросить «1», потому что t-тест с одним данным невозможно выполнить. Можно было бы также сделать 1 + seq_len(nrow(x)-1).

1: for петля, только p-значение

set.seed(2)
out <- rep(NA, nrow(iris))
for (i in seq_len(nrow(iris))[-1]) {
  thisdat <- sample(iris$Sepal.Length, size = i)
  out[i] <- t.test(thisdat, mu = 5)$p.value
}
summary(out)
#      Min.   1st Qu.    Median      Mean   3rd Qu.      Max.      NA's 
# 0.0000000 0.0000000 0.0000000 0.0080013 0.0000001 0.4156151         1 

(Вы можете предположить, что out одинакова для всех последующих примеров, поэтому я не буду это показывать.)

2. *apply, только p-значение

set.seed(2)
out <- sapply(seq_len(nrow(iris))[-1], function(i) {
  thisdat <- sample(iris$Sepal.Length, size = i)
  t.test(thisdat, mu = 5)$p.value
})

sapply принимает вектор и обычно возвращает одно из:

• vector если все возвращаемые значения идеальной длины 1;
• matrix если все возвращаемые значения являются векторами одинаковой длины; или
• list в любое другое время.

Из-за этого некоторые программисты предпочитают lapply (который всегда возвращает list) или vapply (для которого вы должны объявить, какое возвращаемое значение вы ожидаете... и оно терпит неудачу, когда появляется что-то еще). Можно сделать:

set.seed(2)
out <- vapply(seq_len(nrow(iris))[-1], function(i) {
  thisdat <- sample(iris$Sepal.Length, size = i)
  t.test(thisdat, mu = 5)$p.value
}, numeric(1))

(Попробуйте изменить numeric(1) на numeric(2), и вы увидите ошибку values must be length 2, but FUN(X[[1]]) result is length 1.)

Для варианта lapply он очень похож на мой четвертый метод ниже.

Обратите внимание, что length(out) здесь будет nrow(iris)-1, потому что мы пропускаем его во входном векторе seq_len(nrow(iris))[-1]. Это означает, что технически будет разница в summary(out): не будет NA. В остальном все числа равны.

3. for петля, полная модель

Здесь нам нужно хранить гораздо больше, чем просто одно число, поэтому нам нужно хранить его в list.

set.seed(2)
out <- vector("list", nrow(iris))
for (i in seq_len(nrow(iris))[-1]) {
  thisdat <- sample(iris$Sepal.Length, size = i)
  out[[i]] <- t.test(thisdat, mu = 5)
}
str(out[1:3])
# List of 3
#  $ : NULL
#  $ :List of 9
#   ..$ statistic  : Named num 1.31
#   .. ..- attr(*, "names")= chr "t"
#   ..$ parameter  : Named num 1
#   .. ..- attr(*, "names")= chr "df"
#   ..$ p.value    : num 0.416
#   ..$ conf.int   : num [1:2] -2.41 14.11
#   .. ..- attr(*, "conf.level")= num 0.95
#   ..$ estimate   : Named num 5.85
#   .. ..- attr(*, "names")= chr "mean of x"
#   ..$ null.value : Named num 5
#   .. ..- attr(*, "names")= chr "mean"
#   ..$ alternative: chr "two.sided"
#   ..$ method     : chr "One Sample t-test"
#   ..$ data.name  : chr "thisdat"
#   ..- attr(*, "class")= chr "htest"
#  $ :List of 9
#   ..$ statistic  : Named num 1.76
#   .. ..- attr(*, "names")= chr "t"
#   ..$ parameter  : Named num 2
#   .. ..- attr(*, "names")= chr "df"
#   ..$ p.value    : num 0.22
#   ..$ conf.int   : num [1:2] 3.61 8.33
#   .. ..- attr(*, "conf.level")= num 0.95
#   ..$ estimate   : Named num 5.97
#   .. ..- attr(*, "names")= chr "mean of x"
#   ..$ null.value : Named num 5
#   .. ..- attr(*, "names")= chr "mean"
#   ..$ alternative: chr "two.sided"
#   ..$ method     : chr "One Sample t-test"
#   ..$ data.name  : chr "thisdat"
#   ..- attr(*, "class")= chr "htest"

Список довольно длинный, но вы можете видеть, что (1) первый элемент пуст, что неудивительно, поскольку мы пропускаем i из 1; и (2) каждый из элементов после этого содержит все, что вы ожидаете от модели.

Хорошо, давайте пройдемся по этому. Сначала мы выделяем полный список, а затем запускаем цикл for, как и раньше. Единственная разница в цикле заключается в том, что мы сохраняем всю модель (нужно out[[i]] вместо out[i]), а не только $p.value. Теперь, чтобы получить значение p, мы можем использовать цикл for или sapply, я продемонстрирую последний:

head(sapply(out[-1], `[[`, "p.value"))
# [1] 0.41561507 0.22019340 0.05766889 0.08544124 0.03243253 0.09059092

# more verbose, same thing though, showing the "anonymous-function" definition
head(sapply(out[-1], function(m) m$p.value))

Я использовал out[-1], потому что мы знаем, что первый пуст. Мы могли бы легко сделать out <- out[-1] сразу после цикла for выше.

Из моделей можно получить любое другое свойство, например коэффициент модели, используя определение «анонимной функции», которое я продемонстрировал выше.

4. *sapply, полная модель

Это может вас не сильно удивить.

set.seed(2)
out <- lapply(seq_len(nrow(iris))[-1], function(i) {
  thisdat <- sample(iris$Sepal.Length, size = i)
  out[[i]] <- t.test(thisdat, mu = 5)
})

Если вы посмотрите на них, первый элемент не пуст (аналогично приведенному выше примеру sapply), потому что мы даже не запускали и не выделяли его предварительно.

Затем вы можете делать все, что хотите, с отдельными элементами списка:

out[[1]]$p.value
# [1] 0.4156151

str(out[[17]])
# List of 9
#  $ statistic  : Named num 3.98
#   ..- attr(*, "names")= chr "t"
#  $ parameter  : Named num 17
#   ..- attr(*, "names")= chr "df"
#  $ p.value    : num 0.000974
#  $ conf.int   : num [1:2] 5.48 6.57
#   ..- attr(*, "conf.level")= num 0.95
#  $ estimate   : Named num 6.03
#   ..- attr(*, "names")= chr "mean of x"
#  $ null.value : Named num 5
#   ..- attr(*, "names")= chr "mean"
#  $ alternative: chr "two.sided"
#  $ method     : chr "One Sample t-test"
#  $ data.name  : chr "thisdat"
#  - attr(*, "class")= chr "htest"

out[[19]]$statistic
#        t 
# 3.420489 

Если вы хотите получить всю тестовую статистику, аналогичную получению p-значений, вы можете просто сделать:

head(sapply(out, `[[`, "statistic"))
#        t        t        t        t        t        t 
# 1.307692 1.761625 3.000000 2.273030 2.935307 2.014477 

У @ r2evans есть отличный ответ. Я просто сосредоточусь на вашем коде и попытаюсь вывести его на график.

Улучшения включают в себя:

1. Синтаксис for таков: for (i in seq_along(n)) перебирать каждый i. В вашем случае вы действительно хотите сделать for (i in 2:200), потому что i==1 не сможет вычислить p.value.
2. Образец ваших данных должен быть назначен переменной. Как есть, ничего не происходит. Кроме того, вы можете напрямую поместить оператор sample в вызов t.test().
3. Вы хотите сохранить результат каждого цикла в pvalue. Если бы это работало как есть, pvalue закончилось бы последним значением вашего цикла.

Мне нравится серия apply, потому что вам не нужно заранее ничего явно выделять.

set.seed(1)
n <- 50
results <- sapply(seq(2, n)
                  , function(n) {
                   t.test(sample(iris$Sepal.Length, n, replace = T), mu = 5.5, alternative = 'greater')$p.value
                  })

plot(y = results, x = seq(2, n))

Теоретически все, что вам нужно сделать, это заменить iris$Sepal.Length на data$column1 и использовать любой n, который вы предпочитаете.