Построение графика с размерами выборки и оценками мощности - PullRequest
Новая
       39

Построение графика с размерами выборки и оценками мощности

1 голос
/ 19 июня 2019

Я моделировал линейную модель 1000 раз, используя случайно сгенерированные значения роста и веса, и случайным образом назначал каждого участника лечению или не лечению (коэффициент 1 и 0). Допустим, модель была: 

lm(bmi~height + weight + treatment, data = df)

Я сейчас борюсь за следующее:

Теперь модели необходимо циклически проходить по размерам выборки от 300 до 500 с шагом 10 для каждой из 1000 репликаций и сохранять долю смоделированных экспериментов со значениями р менее 0,05 с целью оценки мощности, которая может обнаружить изменение 0,5 bmi между двумя группами лечения при уровне значимости 5%.

После выполнения вышеизложенного мне нужно создать фигуру, которая наилучшим образом отображает размеры выборки по оси x и предполагаемую мощность по оси y, а также отражает наименьший размер выборки для достижения оценки мощности 80%. отличным цветом.

Есть идеи, как и куда идти дальше?

Спасибо, Chris

1 Ответ

0 голосов
/ 19 июня 2019

Я бы сделал это примерно так: 

library(dplyr)
library(ggplot2)

# first, encapsulate the steps required to generate one sample of data
# at a given sample size, run the model, and extract the treatment p-value
do_simulate <- function(n) {
  # use assumed data generating process to simulate data and add error
  data <- tibble(height = rnorm(n, 69, 0.1), 
                 weight = rnorm(n, 197.8, 1.9), 
                 treatment = sample(c(0, 1), n, replace = TRUE),
                 error = rnorm(n, sd = 1.75),
                 bmi = 703 * weight / height^2 + 0.5 * treatment + error)

  # model the data
  mdl <- lm(bmi ~ height + weight + treatment, data = data)

  # extract p-value for treatment effect
  summary(mdl)[["coefficients"]]["treatment", "Pr(>|t|)"]
}

# second, wrap that single simulation in a replicate so that you can perform
# many simulations at a given sample size and estimate power as the proportion
# of simulations that achieve a significant p-value
simulate_power <- function(n, alpha = 0.05, r = 1000) {
  p_values <- replicate(r, do_simulate(n))
  power <- mean(p_values < alpha)
  return(c(n, power))
}

# third, estimate power at each of your desired 
# sample sizes and restructure that data for ggplot
mx <- vapply(seq(300, 500, 10), simulate_power, numeric(2))
plot_data <- tibble(n = mx[1, ], 
                    power = mx[2, ])

# fourth, make a note of the minimum sample size to achieve your desired power
plot_data %>% 
  filter(power > 0.80) %>% 
  top_n(-1, n) %>% 
  pull(n) -> min_n

# finally, construct the plot
ggplot(plot_data, aes(x = n, y = power)) + 
  geom_smooth(method = "loess", se = FALSE) + 
  geom_vline(xintercept = min_n)
...