моделировать отрицательное биномиальное распределение с помощью переменной смещения

Question

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

mutations ~ intercept + beta_cancer + beta_gene + beta_int + offset(log(ntAtRisk)))

, где параметр смещения - максимальное число подсчетов, которое теоретически может произойти.

Создание таблицы с параметрами

ncancers <- 20
ngenes <- 20

beta <- CJ(cancer = as.factor(0:ncancers), gene =  as.factor(0:ngenes))
beta[, beta_cancer := rnorm(n = (ncancers+1), sd = 1)[cancer]]
beta[, beta_gene := rnorm(n = (ngenes+1), sd = 1)[gene]]
beta[, beta_int := rnorm(n = (ngenes+1)*(ncancers+1), sd = 1.5)]
beta[, ntAtRisk := abs(round(rnorm(n = (ngenes+1)*(ncancers+1), mean = 5000, sd  = 2000), digits = 0))[gene]]
beta[, intercept := rnorm(n = (ngenes+1)*(ncancers+1), mean = 2, sd = 1)[gene]]

beta[cancer == "0", c("beta_cancer", "beta_int") := 0] # reference cancer type
beta[gene == "0", c("beta_gene", "beta_int") := 0] # reference gene

Имитация счетчиков мутаций

beta[, mu := exp(intercept + beta_cancer + beta_gene + beta_int + log(ntAtRisk))]
setkey(beta, cancer, gene)

dat <- beta
setkey(dat, cancer, gene)
dat[, mutations := rnbinom(n = nrow(dat), mu = mu, size = 1.5)]
dat[, mutations2 := MASS::rnegbin(n = nrow(dat), 
                                  mu = exp(intercept + beta_cancer + beta_gene + 
                                           beta_int + offset(log(ntAtRisk))), 
                                  theta = 1.5)]

mutations и mutations2 выполняются с использованием различных функций, где * Переменная 1020 * либо включена как обычная переменная, либо, во втором случае, указана как смещение. Тем не менее, тест, который я делаю, не проходит ни одного из них.

Мне нужно, чтобы число мутаций не превышало ntAtRisk, но, к сожалению, это не так. В inte rnet я не смог найти, как я могу включить смещение в симуляцию. Какие у меня варианты?

ggplot(dat, aes(ntAtRisk, mutations+0.5)) +
  geom_point() +
  xlim(0, max(dat$ntAtRisk)) + 
  ylim(0, max(dat$ntAtRisk)) + 
  geom_abline(color = "red") 

Answer 1

Когда вы подбираете glm для Пуассона, Негбина со смещением, сумма ваших коэффициентов и перехватов не может быть больше 1, потому что log (смещение) вычитается из log (ответа), и это всегда <1, для пример: </p>

n=seq(100,1000,by=100)
mu = n/5
y = rnbinom(n = 10,size =1.5,mu=mu)
glm.nb(y~1+offset(log(n)))

Call:  glm.nb(formula = y ~ 1 + offset(log(n)), init.theta = 1.217692649, 
    link = log)

Coefficients:
(Intercept)  
     -1.424 

Это очень сложная симуляция для настройки из-за ограничений, в вашем случае, я предлагаю установить перехват как нечто очень низкое, так как, скорее всего, мутации (если я правильно понял), в любом случае не так часто:

set.seed(222)
beta <- CJ(cancer = as.factor(0:ncancers), gene =  as.factor(0:ngenes))
beta[, beta_cancer := rnorm(n = (ncancers+1))[cancer]]
beta[, beta_gene := rnorm(n = (ngenes+1))[gene]]
beta[, beta_int := rnorm(n = (ngenes+1)*(ncancers+1))]
beta[, ntAtRisk := abs(round(rnorm(n = (ngenes+1)*(ncancers+1), mean = 5000, sd  = 2000), digits = 0))[gene]]
beta[, intercept := runif(n = (ngenes+1)*(ncancers+1),min=-5,max=-3)[gene]]
beta[cancer == "0", c("beta_cancer", "beta_int") := 0] # reference cancer type
beta[gene == "0", c("beta_gene", "beta_int") := 0] # reference gene

На этом этапе вы будете учитывать смещение, добавляя лог-термин, нет необходимости добавлять смещение позже:

beta[, mu := exp(intercept + beta_cancer + beta_gene + beta_int + log(ntAtRisk))]
setkey(beta, cancer, gene)

Теперь мы моделируем данные, предоставляя среднее значение как mu, и вы задаете постоянное значение тета:

dat <- beta
setkey(dat, cancer, gene)
dat[, mutations := rnbinom(n = nrow(dat), mu = mu, size = 1.5)]

ggplot(dat, aes(ntAtRisk, mutations+0.5)) +
  geom_point() +
  xlim(0, max(dat$ntAtRisk)) + 
  ylim(0, max(dat$ntAtRisk)) + 
  geom_abline(color = "red") 

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