Как подсчитать количество значимых p-значений в симуляции?
Допустим, у меня есть клиническое испытание, в котором 32 мыши инфицированы болезнью. Я лечу 16 из них своим лечением, а остальные 16 оставляю под контролем. Перед проведением эксперимента я запускаю моделирование. Я сделаю 32 наблюдения из равномерного (0,1) распределения. Я хотел бы подсчитать, сколько раз я достиг значения p <0,01 после завершения моделирования. Я чувствую, что почти у цели, но не знаю, как это собрать:
nSims <- 10000 #number of simulated experiments
p <-numeric(nSims) #set up empty container for all simulated p-values
sig<-0
for(i in 1:nSims){ #for each simulated experiment
#generating 32 observations total from uniform(0,1) distribution
control.year1 <- runif (16, min = 0, max = 1)
treat.year1 <- runif (16, min = 0, max = 1)
#Creating dichotomous variable for those get better/don't get better
control.respond <- ifelse(control.year1<=0.05,1,0)
treat.respond <- ifelse(treat.year1<=0.30,1,0)
#perform t-test
z <- t.test(control.respond,treat.respond)
p[i]<-z$p.value #get the p-value and store it
# want to count number of significant p-values - not sure how to do it
significance <- ifelse(z$p.value<= 0.01,sum(sig, 1),0)
}
Разве вы не можете решить это аналитически?
вас могут заинтересовать некоторые встроенные калькуляторы мощности в R (apropos("^power"))
@BenBolker, спасибо, это не последняя симуляция - я моделирую пятилетний эксперимент, который может прекратиться раньше, если лечение будет эффективным, я просто пытаюсь правильно понять цикл для первого года эксперимента.
@Marius Спасибо за внимание, вы правы, я собираюсь использовать односторонний точный тест рыболова.
@Mob Я мог бы, но для этой симуляции я имитирую эксперимент, который продлится 5 лет, но может прекратиться раньше, если лечение будет эффективным.
Ответы 1
Это решение не самое элегантное, но для обработки данных использует magrittr и dplyr. Сначала я создал матрицу для хранения ваших смоделированных данных:
library(magrittr)
library(dplyr)
n <- 100
control.years <- as.data.frame(matrix(runif (16*n, min=0, max=1),ncol=16))
treat.years <- as.data.frame(matrix(runif (16*n, min=0, max=1),ncol=16))
Затем я создал структуру данных для сбора p-значений для всех t-тестов:
for (i in 1:n) { p[i] <- t.test(control.years[i,],treat.years[i,])$p.value }
Вы можете отфильтровать p-значения в желаемом диапазоне:
> as.data.frame(p) %>% filter(p<0.05)
p
1 0.03173299
2 0.01652114
3 0.00471807
Или вы можете создать новую переменную, которая сообщает вам, является ли она значимой или нет:
> as.data.frame(p) %>% mutate(sig=ifelse(p<0.05,1,0))
p sig
1 0.65233254 0
2 0.50731231 0
3 0.11657045 0
...
29 0.03173299 1
Или вы можете узнать, сколько значимых p-значений:
> z <- as.data.frame(p) %>% mutate(sig=ifelse(p<0.05,1,0))
> table(z$sig)
0 1
97 3
спасибо, это именно то, о чем я думал, и спасибо за разные варианты.
Вы должны иметь возможность подсчитать количество значимых результатов в конце, вне цикла for, с помощью
sum(p < 0.01). Но вы проводите t-тест для двоичного результата? Я не уверен, имеет ли это смысл.