我正在尝试使用已知参数模拟变异数据,以将其进一步用于测试回归函数。在此模拟中,我希望突变计数取决于变量:
mutations ~ intercept + beta_cancer + beta_gene + beta_int + offset(log(ntAtRisk)))
其中offset参数是理论上可以发生的最大计数。
用参数创建表
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
使用不同的函数制作,其中offset
变量要么作为普通变量包含,要么在第二种情况下指定为偏移量。但是,我正在做的测试没有通过任何一项。
我需要突变计数不大于ntAtRisk,但不幸的是事实并非如此。我在互联网上找不到如何将偏移量包括在模拟中的功能。我有什么选择?
ggplot(dat, aes(ntAtRisk, mutations+0.5)) +
geom_point() +
xlim(0, max(dat$ntAtRisk)) +
ylim(0, max(dat$ntAtRisk)) +
geom_abline(color = "red")
当您为带偏移量的泊松negbin拟合glm时,系数和截距的总和不能大于1,因为log(offset)是从log(response)中减去的,并且始终小于1,例如:
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作为平均值,然后指定一个恒定的theta值:
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。您可以编写代码来手动更正此问题,或者如果您的预测如此之高,我想您需要真正检查数据。
本文收集自互联网,转载请注明来源。
如有侵权,请联系 [email protected] 删除。
我来说两句