Расчет мощности путем симуляции для glmer в R - PullRequest
Новая
       65

Расчет мощности путем симуляции для glmer в R

0 голосов
/ 08 марта 2019

Я пытаюсь выполнить расчет мощности путем моделирования для следующей ситуации.

  1. Доступны следующие воспроизводимые данные:

structure(list(id = c("151", "151", "151", "151", "151", "158", "158", "158", "158", "166", "166", "166", "166", "166", "173", "173", "173", "173", "173", "173", "173", "176", "176", "176", "176", "176", "176", "176", "176", "176", "176", "201", "201", "201", "201", "201", "213", "213", "213", "213", "213", "213", "213", "219", "219", "219", "219", "219", "219", "219", "219", "219", "220", "220", "220", "220", "220", "220", "220", "221", "221", "221", "221", "221", "221", "221", "221", "227", "227", "227", "227", "227", "227", "227", "227", "227", "227", "228", "228", "228", "228", "228", "228", "231", "231", "231", "231", "231", "231", "234", "234", "234", "234", "234", "234", "234", "234", "234", "246", "246", "246", "246", "246", "246", "246", "246", "246", "246", "247", "247", "247", "247", "247", "247", "247", "247", "247", "247", "261", "261", "261", "261", "261", "261", "266", "266", "266", "266", "266", "266", "266", "266", "266", "266", "273", "273", "273", "273", "273", "273", "273", "273", "273", "276", "276", "276", "276", "276", "276", "276", "276", "276", "276", "287", "287", "287", "287", "287", "287", "287", "287", "304", "304", "304", "304", "304", "304", "304", "304", "310", "310", "310", "310", "310", "310", "310", "312", "312", "312", "312", "312", "312", "312", "312", "312", "312", "318", "318", "318", "318", "318", "318", "318", "318", "327", "327", "327", "327", "327", "327", "327", "327", "332", "332", "332", "332", "332", "332", "332", "332", "332"), arm = structure(c(1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("C", "I"), class = "factor"), timepoint = structure(c(1L, 1L, 1L, 3L, 3L, 1L, 3L, 2L, 2L, 1L, 1L, 1L, 3L, 3L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 2L, 2L, 1L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 3L, 3L, 3L, 2L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 2L, 2L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 3L, 3L, 2L, 2L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 2L, 2L, 1L, 1L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 3L, 3L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 2L, 2L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 2L, 2L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 2L, 2L, 1L, 1L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 2L, 2L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 2L, 2L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 2L, 2L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 2L, 2L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 3L, 3L, 2L, 2L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 2L), .Label = c("b", "i", "f"), class = "factor"), social = c(1, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2 ), deno = c(2L, 1L, 1L, 2L, 3L, 3L, 1L, 2L, 1L, 3L, 1L, 1L, 6L, 3L, 1L, 1L, 3L, 9L, 3L, 3L, 3L, 2L, 5L, 3L, 2L, 2L, 6L, 3L, 4L, 1L, 3L, 1L, 11L, 1L, 5L, 2L, 2L, 7L, 1L, 1L, 3L, 4L, 6L, 4L, 1L, 6L, 3L, 19L, 3L, 3L, 10L, 12L, 1L, 4L, 1L, 1L, 2L, 11L, 1L, 5L, 1L, 1L, 4L, 9L, 7L, 1L, 6L, 5L, 2L, 9L, 3L, 2L, 6L, 2L, 4L, 2L, 14L, 1L, 2L, 1L, 7L, 1L, 1L, 1L, 1L, 2L, 3L, 1L, 3L, 13L, 5L, 14L, 17L, 8L, 8L, 4L, 6L, 7L, 10L, 2L, 8L, 6L, 4L, 4L, 2L, 1L, 7L, 11L, 9L, 1L, 2L, 1L, 4L, 2L, 1L, 5L, 1L, 6L, 3L, 4L, 1L, 9L, 1L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 12L, 1L, 2L, 3L, 1L, 4L, 4L, 1L, 7L, 2L, 1L, 9L, 2L, 1L, 3L, 8L, 3L, 3L, 4L, 3L, 3L, 5L, 4L, 5L, 10L, 2L, 2L, 3L, 2L, 3L, 4L, 2L, 3L, 2L, 6L, 3L, 1L, 7L, 7L, 1L, 4L, 5L, 2L, 4L, 2L, 3L, 6L, 4L, 2L, 4L, 1L, 2L, 4L, 2L, 4L, 8L, 3L, 3L, 2L, 2L, 1L, 4L, 2L, 2L, 1L, 1L, 3L, 1L, 2L, 1L, 8L, 2L, 9L, 4L, 2L, 4L, 2L, 1L, 3L, 1L, 2L), nume = c(0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 2L, 0L, 7L, 2L, 1L, 6L, 8L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 2L, 3L, 1L, 0L, 3L, 2L, 0L, 3L, 1L, 1L, 4L, 1L, 3L, 0L, 7L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 3L, 2L, 2L, 1L, 0L, 2L, 1L, 5L, 1L, 2L, 1L, 0L, 1L, 1L, 0L, 1L, 6L, 3L, 0L, 0L, 1L, 0L, 1L, 0L, 2L, 0L, 3L, 0L, 0L, 0L, 2L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 2L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 1L, 3L, 1L, 1L, 1L, 0L, 2L, 3L, 0L, 0L, 3L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 2L, 1L, 1L, 0L, 4L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 2L, 3L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 2L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L)), class = "data.frame", row.names = c(NA, -211L))

Существует несколько наблюдений за "цифрой" и "дено" в разные "моменты времени" для каждого из лиц под столбцом "id".«социальный» - ковариата, не особенно интересная для этого вопроса.

Я смоделировал отношение число / deno со смешанной (по id) регрессии Пуассона, как

model <- glmer(nume ~ arm*timepoint + (1|id) , data=data, subset=timepoint%in%c("b","f"), offset=log(deno), family='poisson')

Существует смещение "deno«для учета воздействия каждого наблюдения« числа », а подмножество данных берется только в два момента времени« b »(базовый уровень) и« f »(продолжение), потому что вы можете заметить, что не всеу людей есть «i» (промежуточный) момент времени.

Это кажется довольно разумным.

Теперь я хотел бы рассчитать мощность этого объекта, чтобы определить отношение шансов для взаимодействия, скажем, равное 0.5, 1, 1.5.Я хочу знать, в порядке ли размер моей выборки 27 особей, и что было бы (если так) лучше.

Я настроил следующий смоделированный набор данных, учитывая то же числоиндивидуумы реального случая и одинаковые пропорции в двух руках, одинаковые пропорции в gender и social.Я предполагаю 4 измерений в каждый момент времени для каждого человека, которых у меня нет на самом деле.Фоновая структура набора данных является фиксированной и не входит в симуляцию (то есть, например, какой человек является мужчиной / женщиной и т. Д.) 

n.ind       <- 27
prop.arms   <- 0.444
prop.gender <- 0.852
prop.social <- 0.556
howmany.b   <- 4
howmany.i   <- 4
howmany.f   <- 4

id     <- 1:n.ind
arm    <- sample(rep(c('C','I'), c(round(n.ind*prop.arms,0), n.ind-  
          round(n.ind*prop.arms,0) ) ), n.docs, replace=F)
gender <- sample(rep(c('M','F'), c(round(n.ind*prop.gender,0), n.ind-
          round(n.ind*prop.gender,0) ) ), n.ind, replace=F)
social <- sample(rep(c(1,2), c(round(n.ind*prop.social,0), n.ind-
          round(n.ind*prop.social,0) ) ), n.ind, replace=F)
stage.i <- sapply(arm, function(x) ifelse(x%in%"C", 0, 4) )

Сначала я настроил фрейм данных с функциями

ind.data <- data.frame(id=id, arm=arm, gender=gender, social=social, 
            stage.b=rep(howmany.b, n.docs), stage.i=stage.i, 
            stage.f=rep(howmany.f, n.docs))
ind.data$howmany <- 
            apply(doc.data[,c('stage.b','stage.i','stage.f')],1,sum)

ind    <- rep(1: n.ind, ind.data $howmany)
arm    <- rep(doc.data$arm, ind.data $howmany)
gender <- rep(doc.data$gender, ind.data $howmany )
social <- rep(doc.data$social, ind.data $howmany)
when   <- rep(rep(c('b','i','f'),n.ind), 
                as.vector(t(as.matrix(
                doc.data[,c('stage.b','stage.i','stage.f')]))) )

Затем я установил «фиксированную» часть того, что будет моделируемым набором данных (год, пол, социальный статус и, когда измеряется, для каждого человека): 

data.sim <- data.frame(ind=ind, arm=arm, gender=gender, 
            social=social, when=when)

## denominators are random in reality
data.sim$deno  <- sample( c(1:14,17,19), nrow(data.sim), 
                  prob=prop.table(table(data$deno)), replace=T)

Здесь запускается функция симуляции: 

sim1 <- function(bArm, bwhen, bint, b0, Vind, Verror ){

#random effect on individuals   
D.re <- rnorm(1:n.ind, 0, sqrt(Vind))
# residuals
eps  <- rnorm(nrow(data.sim), 0, sqrt(Verror))

# simulation from model  
data.sim$nume <- rpois(nrow(data.sim), 
             data.sim$deno*exp(D.re[ind] + bArm*(arm=='I') + 
             bwhen*(when=='f') + bint*(arm=='I')*(when=='f') + eps))

# fit the model as done with real data
fit1 <- glmer(nume ~ arm*when + (1|ind), data=data.sim, 
        subset=when%in%c("b","f"), offset=log(deno), family='poisson')

# check p-val of interaction
summary(fit1)$coef[4,4]
}

Теперь запустите симуляции: 

n.simul = 50
library(tcltk)
pb <- tkProgressBar(max=n.simul)
out1 <- replicate(n.simul, {setTkProgressBar(pb, 
        getTkProgressBar(pb)+1); sim1(bArm=0, bwhen=0, bint=0, b0=0, 
        Vdoc=1, Verror=1)})

mean(out1<0.05)

Этот код взят из адаптации к моему случаю this , который, как мне показалось, был довольно близок к моему случаю.

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

Мой вопрос прост: я не знаю, правильно ли я это делаю, и я действительно не знаю, какие параметры следует ввести в симуляцию.

Я понимаю, что, если меня интересует соотношение шансов взаимодействия, мой bint в симуляции должен быть 0.5, 1, 1.5, как указано выше, но как насчет других?Могу ли я также объяснить влияние на отдельные ковариаты?

Я знаю, что это довольно длинный вопрос, но я был бы очень признателен, если бы вы могли помочь мне с этим немного.

Лучше всего,EM

...