Этот вопрос является ответом на ответ @Greg Snow на вопрос, который я задал относительно анализа мощности с помощью логистической регрессии и SAS Proc GLMPOWER.

Если я планирую эксперимент и проанализирую результаты в факторной логистической регрессии, как я могу использовать симуляцию (и здесь ) для анализа мощности?

Вот простой пример, где есть две переменные, первая принимает три возможных значения {0,03, 0,06, 0,09}, а вторая - фиктивный индикатор {0,1}. Для каждого мы оцениваем уровень ответов для каждой комбинации (количество респондентов / количество людей, продаваемых на рынок). Кроме того, мы хотим иметь в 3 раза больше первой комбинации факторов, чем другие (которые можно считать равными), потому что эта первая комбинация является нашей проверенной верной версией. Это установка, подобная приведенной в курсе SAS, упомянутом в связанном вопросе.

Модель, которая будет использоваться для анализа результатов, будет логистической регрессией с основными эффектами и взаимодействием (ответ 0 или 1).

mod <- glm(response ~ Var1 + Var2 + I(Var1*Var2))

Как я могу имитировать набор данных для использования с этой моделью для анализа мощности?

Когда я запускаю это через SAS Proc GLMPOWER(использование STDDEV =0.05486016 которого соответствует sqrt(p(1-p))где p - средневзвешенное значение показанных ответов):

data exemplar;
  input Var1 $ Var2 $ response weight;
  datalines;
    3 0 0.0025  3
    3 1 0.00395 1
    6 0 0.003   1
    6 1 0.0042  1
    9 0 0.0035  1
    9 1 0.002   1;
run;

proc glmpower data=exemplar;
  weight weight;
  class Var1 Var2;
  model response = Var1 | Var2;
  power
    power=0.8
    ntotal=.
    stddev=0.05486016;
run;

Примечание: GLMPOWERбудут использоваться только переменные класса (именные), поэтому 3, 6, 9, приведенные выше, обрабатываются как символы и могут иметь значения low, mid и high или любые другие три строки. Когда будет проведен реальный анализ, для расчета любой кривизны будет использоваться Var1 с числовым значением (и мы будем включать полиномиальный термин Var1 * Var1).

Выход из SAS

Таким образом, мы видим, что нам нужно 762,112 как размер нашей выборки (основной эффект Var2 труднее всего оценить) с мощностью, равной 0,80, и альфа, равной 0,05. Мы бы распределили их так, чтобы базовая комбинация была в 3 раза больше (т. Е. 0,375 * 762112), а остальные просто в равной степени попадали в другие 5 комбинаций.

Отборочные:

• $N$$ES$$\alpha$

  • $N$$\alpha=.05$
  • $N$

• В дополнение к превосходному посту @ GregSnow, здесь можно найти еще одно действительно замечательное руководство по анализу мощности на основе моделирования на основе CV: расчет статистической мощности . Подводя итог основным идеям:

  1. выяснить эффект, который вы хотите иметь возможность обнаружить
  2. генерировать N данных из этого возможного мира
  3. выполнить анализ, который вы намереваетесь провести над этими ложными данными
  4. Сохраните, являются ли результаты «значительными» в соответствии с выбранной вами альфа
  5. $B$$N$
  6. $N$

• $p$$p$$B$$p$$B$

• В R основным способом генерирования двоичных данных с заданной вероятностью «успеха» является ? Rbinom

  • Например, чтобы получить количество успехов из 10 испытаний Бернулли с вероятностью p, код будет rbinom(n=10, size=1, prob=p)(вы, вероятно, захотите присвоить результат переменной для хранения)
  • Вы также можете генерировать такие данные менее элегантно, используя ? runif , например,ifelse(runif(1)<=p, 1, 0)
  • если вы считаете, что результаты опосредованы скрытой гауссовой переменной, вы можете сгенерировать скрытую переменную как функцию ваших ковариат с помощью команды ? rnorm , а затем преобразовать их в вероятности pnorm()и использовать их в своем rbinom()коде.

• $var1^2$$var1*var2$$var1^2*var2$

• Хотя мой ответ здесь написан в контексте другого вопроса, мой ответ здесь: Разница между логит-моделями и пробит-моделями содержит много базовой информации об этих типах моделей.

• Так же , как существуют различные виды коэффициентов ошибок типа I , когда есть несколько гипотез (например, за отличие частоты ошибок , коэффициент ошибок familywise , и за семью частоты ошибок ), поэтому существуют различные виды питания * (например, для а один предварительно заданный эффект , для любого эффекта и для всех эффектов ). Вы также можете искать возможности для обнаружения определенной комбинации эффектов или возможности одновременной проверки модели в целом. Из вашего описания вашего кода SAS я думаю, что он ищет последний. Однако, исходя из вашего описания вашей ситуации, я предполагаю, что вы хотите обнаружить эффекты взаимодействия как минимум.

  • * ссылка: Maxwell, SE (2004). Постоянство слабых исследований в психологических исследованиях: причины, последствия и средства правовой защиты. Психологические методы , 9 , 2 , с. 147-163.
  • Ваши эффекты довольно малы (не путать с низкой частотой откликов), поэтому нам будет трудно добиться хорошей мощи.
  • Обратите внимание на то, что, хотя все они звучат довольно схожими, они во многом не совпадают (например, очень возможно получить значимую модель без существенных эффектов - обсуждается здесь: как регрессия может быть существенной, но все предикторы не значительный? или значительный эффект, но там, где модель несущественна - обсуждается здесь: Значение коэффициентов в линейной регрессии: значимый t-критерий против незначимой F-статистики ), что будет проиллюстрировано ниже.

• Чтобы по-другому взглянуть на проблемы, связанные с властью, см. Мой ответ здесь: Как сообщить общую точность оценки корреляций в контексте обоснования размера выборки.

Простая последующая мощность для логистической регрессии в R:

Допустим, ваши положительные ответы отражают реальную ситуацию в мире и что вы разослали 10 000 писем. Какова сила, чтобы обнаружить эти эффекты? (Обратите внимание, что я известен тем, что написал «комично неэффективный» код, ниже предполагается, что за ним проще следить, а не оптимизировать для эффективности; на самом деле он довольно медленный.)

set.seed(1)

repetitions = 1000
N = 10000
n = N/8
var1  = c(   .03,    .03,    .03,    .03,    .06,    .06,    .09,   .09)
var2  = c(     0,      0,      0,      1,      0,      1,      0,     1)
rates = c(0.0025, 0.0025, 0.0025, 0.00395, 0.003, 0.0042, 0.0035, 0.002)

var1    = rep(var1, times=n)
var2    = rep(var2, times=n)
var12   = var1**2
var1x2  = var1 *var2
var12x2 = var12*var2

significant = matrix(nrow=repetitions, ncol=7)

startT = proc.time()[3]
for(i in 1:repetitions){
  responses          = rbinom(n=N, size=1, prob=rates)
  model              = glm(responses~var1+var2+var12+var1x2+var12x2, 
                           family=binomial(link="logit"))
  significant[i,1:5] = (summary(model)$coefficients[2:6,4]<.05)
  significant[i,6]   = sum(significant[i,1:5])
  modelDev           = model$null.deviance-model$deviance
  significant[i,7]   = (1-pchisq(modelDev, 5))<.05
}
endT = proc.time()[3]
endT-startT

sum(significant[,1])/repetitions      # pre-specified effect power for var1
[1] 0.042
sum(significant[,2])/repetitions      # pre-specified effect power for var2
[1] 0.017
sum(significant[,3])/repetitions      # pre-specified effect power for var12
[1] 0.035
sum(significant[,4])/repetitions      # pre-specified effect power for var1X2
[1] 0.019
sum(significant[,5])/repetitions      # pre-specified effect power for var12X2
[1] 0.022
sum(significant[,7])/repetitions      # power for likelihood ratio test of model
[1] 0.168
sum(significant[,6]==5)/repetitions   # all effects power
[1] 0.001
sum(significant[,6]>0)/repetitions    # any effect power
[1] 0.065
sum(significant[,4]&significant[,5])/repetitions   # power for interaction terms
[1] 0.017

Таким образом, мы видим, что 10000 писем на самом деле не достигают 80% мощности (любого рода) для обнаружения этих ответов. (Я не совсем уверен в том, что делает код SAS, чтобы объяснить резкое несоответствие между этими подходами, но этот код концептуально прост - если медленен - ​​и я потратил некоторое время на его проверку, и я думаю, что они результаты разумны.)

Основанная на моделировании априорная мощность для логистической регрессии:

$N$$N$$N$$N$

$N$$N$

sum(significant[,1])/repetitions      # pre-specified effect power for var1
[1] 0.115
sum(significant[,2])/repetitions      # pre-specified effect power for var2
[1] 0.091
sum(significant[,3])/repetitions      # pre-specified effect power for var12
[1] 0.059
sum(significant[,4])/repetitions      # pre-specified effect power for var1X2
[1] 0.606
sum(significant[,5])/repetitions      # pre-specified effect power for var12X2
[1] 0.913
sum(significant[,7])/repetitions      # power for likelihood ratio test of model
[1] 1
sum(significant[,6]==5)/repetitions   # all effects power
[1] 0.005
sum(significant[,6]>0)/repetitions    # any effect power
[1] 0.96
sum(significant[,4]&significant[,5])/repetitions   # power for interaction terms
[1] 0.606

$var1^2$significant

$N$

Ответ @ Gung отлично подходит для понимания. Вот подход, который я бы использовал:

mydat <- data.frame( v1 = rep( c(3,6,9), each=2 ),
    v2 = rep( 0:1, 3 ), 
    resp=c(0.0025, 0.00395, 0.003, 0.0042, 0.0035, 0.002) )

fit0 <- glm( resp ~ poly(v1, 2, raw=TRUE)*v2, data=mydat,
    weight=rep(100000,6), family=binomial)
b0 <- coef(fit0)


simfunc <- function( beta=b0, n=10000 ) {
    w <- sample(1:6, n, replace=TRUE, prob=c(3, rep(1,5)))
    mydat2 <- mydat[w, 1:2]
    eta <- with(mydat2,  cbind( 1, v1, 
                v1^2, v2,
                v1*v2,
                v1^2*v2 ) %*% beta )
    p <- exp(eta)/(1+exp(eta))
    mydat2$resp <- rbinom(n, 1, p)

    fit1 <- glm( resp ~ poly(v1, 2)*v2, data=mydat2,
        family=binomial)
    fit2 <- update(fit1, .~ poly(v1,2) )
    anova(fit1,fit2, test='Chisq')[2,5]
}

out <- replicate(100, simfunc(b0, 10000))
mean( out <= 0.05 )
hist(out)
abline(v=0.05, col='lightgrey')

Эта функция проверяет общий эффект v2, модели могут быть изменены, чтобы посмотреть на другие типы тестов. Мне нравится писать это как функцию, так что, когда я хочу протестировать что-то другое, я могу просто изменить аргументы функции. Вы также можете добавить индикатор выполнения или использовать параллельный пакет для ускорения процесса.

Здесь я только что сделал 100 повторений, я обычно начинаю с этого уровня, чтобы найти приблизительный размер выборки, затем увеличиваю количество итераций, когда нахожусь в правильном парке мячей (нет необходимости тратить время на 10000 итераций, когда у вас есть 20% мощности).