可以与BKMR包联用,估计环境混合物的影响及环境化学物暴露间的交互作用,BKMR包可探讨化学物间交互作用,但无法提供交互作用的PIP值,NLinteracton包可弥补这一缺陷。
安装
library(devtools)
install_github(repo = "jantonelli111/NLinteraction")
library(NLinteraction)
构建模拟数据集
n = 100
p = 10
pc = 1
X = matrix(rnorm(n*p), n, p)
C = matrix(rnorm(n*pc), nrow=n)
TrueH = function(X) {
return(1.5*(X[,2]*X[,3]) - 1.6*(X[,4]^2 * X[,5]))
}
Y = 5 + C + TrueH(X) + rnorm(n)
数据集中有10种暴露,一个协变量及结局变量,X与Y为非线性。
MCMC方程
接下来构建模型,首先确定模型中自由度,这里选了自由度为2和3(模型中的ns参数),nlter为MCMC迭代次数,其他参数也可以改变,建议用标准参数除非有强烈的先验信息及完全知道数据的先验分布 。
NLmod2 = NLint(Y=Y, X=X, C=C, nIter=10000, nBurn=5000, thin=5, nChains=2, ns=2)
NLmod3 = NLint(Y=Y, X=X, C=C, nIter=10000, nBurn=5000, thin=5, nChains=2, ns=3)
如果C或者X是数据库的话,需要转换为矩阵:
NLmod2 = NLint(Y=Y, X=as.matrix(X), C=as.matrix(C),
nIter=10000, nBurn=5000, thin=5, nChains=2, ns=2)
比较模型的WAIC
print(c(NLmod2$waic,NLmod3$waic))
自由度为2时WAIC最小,推荐比较自由度为1-5是WAIC的值,这里只展示了2和3。现在选择NLmod2为最终选择的model.
NLmod = NLmod2
Posterior inclusion probabilities
计算PIP值
NLmod$MainPIP
可视化
barplot(NLmod$MainPIP)
二阶交互作用的PIP值
NLmod$InteractionPIP
可视化
plotInt(NLmod = NLmod)
可计算任意化学物间交互作用
InteractionProb(NLmod=NLmod, Xsub=c(4,5,9))
X4,X5与X9的交互作用的PIP值为0。该方程可用于计算任意交互作用。
可视化暴露响应曲线
以下为X4与X5的交互作用的示例,我们将在X5的不同水平上绘制X4与结局的关联。
par(mfrow=c(1,3), pty='s')
plotSurface1d(NLmod = NLmod, X=X, C=C, j1=4, j2=5,
gridLength=30, quantile_j2=0.2, quantile_rest = 0.5,
xlab="X4", ylab="Posterior predictive distribution",
main="20th quantile of X5")
plotSurface1d(NLmod = NLmod, X=X, C=C, j1=4, j2=5,
gridLength=30, quantile_j2=0.5, quantile_rest = 0.5,
xlab="X4", ylab="Posterior predictive distribution",
main="50th quantile of X5")
plotSurface1d(NLmod = NLmod, X=X, C=C, j1=4, j2=5,
gridLength=30, quantile_j2=0.8, quantile_rest = 0.5,
xlab="X4", ylab="Posterior predictive distribution",
main="80th quantile of X5")
我们看到,曲线在X5的不同水平时看起来截然不同,因为它们之间具有很强的交互作用。
References
Joseph Antonelli, Maitreyi Mazumdar, David Bellinger, David C. Christiani, Robert Wright, Brent A. Coull. Estimating the health effects of environmental mixtures using Bayesian semiparametric regression and sparsity inducing priors. 2019. The Annals of Applied Statistics. To appear.