Как рассчитать статистическую мощность на величину эффекта регрессии Кокса в R?

Question

Мне очень жаль, что это очень общий вопрос, но мои навыки работы с Google полностью меня подводят.

Я бы хотел изобразить размер выборки по мощности для каждого размера эффекта (или не менее 0,1) в соответствии с тестом Кокса PH.Коллега (теперь слева) произвел анализ мощности линейной модели в R как заполнитель для регрессии Кокса.Есть пакеты (pwr), но ни один из них не касается регрессии Кокса.

Мне нужно сделать это в R, чтобы я мог наложить результаты анализа Кокса поверх этого анализа мощности, чтобы показать, насколько точны мои результаты на самом деле.

Это скриншот того, что мне представили, и, по сути, того, что я хотел бы воссоздать, не выдумывая.

Answer 1

Мощность определяется как условная вероятность Pr(reject H0 | H1 is true). Можно показать , что мощность простого одностороннего теста t может быть аппроксимирована следующей функцией

B <- function(beta, n, sd = 6) {
    sem <- sd / sqrt(n)
    1 - pnorm(1.64 - beta / sem)
}

Здесь beta - эторазмер эффекта, и в этом случае соответствует среднему значению для популяции (или на языке линейных моделей бета соответствует перехвату).n - это размер выборки, а sd - «настраиваемый» параметр, описывающий стандартное отклонение (таким образом, неопределенность, которую вы имеете в измерениях).Я выбрал это значение таким образом, чтобы мы (более или менее) воспроизводили результаты из вашей фигуры.

Теперь мы можем рассчитать мощность B для различных значений размера выборки n и величины эффекта beta,Мы выбираем те же значения, что и на графике, который вы показываете.

n <- c(100, 500, 2500, 10000, 50000, 250000)
beta <- c(0.1, 0.2, 0.5, 1.0)

library(tidyverse)
outer(beta, n, B) %>%
    data.frame(row.names = beta) %>%
    setNames(n) %>%
    rownames_to_column("beta") %>%
    gather(n, power, -beta) %>%
    mutate(
        n = as.numeric(n),
        beta = factor(beta, unique(beta))) %>%
    ggplot(aes(n, power, colour = beta, group = beta)) +
    geom_line() +
    geom_point() +
    scale_y_continuous(labels = scales::percent) +
    scale_x_log10() +
    labs(x = "sample size (N)", y = "power (%)")

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