Treatments effectslibrary(iCARH)
## Loading required package: rstan
## Loading required package: StanHeaders
## Loading required package: ggplot2
## rstan (Version 2.19.3, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## Loading required package: MASS
## Loading required package: glue
library(abind)
Tp=4L # timepoints
N=10L # number of samples
J=14L # number of metabolites
K=2L # number of bacteria species
P=8L # number of pathways
set.seed(12473)
For real data Build pathway matrices using iCARH.getPathwaysMat. Elements in KEGG id list may contain multiple KEGG ids per metabolite. If KEGG id unknown : “Unk[number]”.
keggid = list("Unk1", "C03299","Unk2","Unk3",
c("C08363", "C00712"), # allowing multiple ids per metabolite
)
pathways = iCARH.GetPathwaysMat(keggid, "rno")
To simulate data use iCARH.simulate
# Manually choose pathways
path.names = c("path:map00564","path:map00590","path:map00061","path:map00591",
"path:map00592","path:map00600","path:map01040","path:map00563")
data.sim = iCARH.simulate(Tp, N, J, P, K, path.names=path.names, Zgroupeff=c(0,4),
beta.val=c(1,-1,0.5, -0.5))
## Warning in iCARH.simulate(Tp, N, J, P, K, path.names = path.names, Zgroupeff = c(0, : Number of pathways reduced to 1 due to random selection of metabolites
## in the intially specified pathways.
XX = data.sim$XX
Y = data.sim$Y
Z = data.sim$Z
pathways = data.sim$pathways
XX[2,2,2] = NA #missing value example
Check inaccuracies between covariance and design matrices
pathways.bin = lapply(pathways, function(x) { y=1/(x+1); diag(y)=0; y})
adjmat = rowSums(abind::abind(pathways.bin, along = 3), dims=2)
cor.thresh = 0.7
# check number of metabolites in same pathway but not correlated
for(i in 1:Tp) print(sum(abs(cor(XX[i,,])[which(adjmat>0)])<cor.thresh))
## [1] 108
## [1] NA
## [1] 146
## [1] 124
Run iCARH model.
rstan::rstan_options(auto_write = TRUE)
options(mc.cores = 2)
# For testing, does not converge
fit.sim = iCARH.model(XX, Y, Z, groups=rep(c(0,1), each=5), pathways = pathways, control = list(max_treedepth=8),
iter = 4, chains = 1)
##
## SAMPLING FOR MODEL 'iCARH' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.001256 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 12.56 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: WARNING: No variance estimation is
## Chain 1: performed for num_warmup < 20
## Chain 1:
## Chain 1: Iteration: 1 / 4 [ 25%] (Warmup)
## Chain 1: Iteration: 2 / 4 [ 50%] (Warmup)
## Chain 1: Iteration: 3 / 4 [ 75%] (Sampling)
## Chain 1: Iteration: 4 / 4 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.005718 seconds (Warm-up)
## Chain 1: 0.003633 seconds (Sampling)
## Chain 1: 0.009351 seconds (Total)
## Chain 1:
## Warning: There were 2 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: Examine the pairs() plot to diagnose sampling problems
## Warning: The largest R-hat is NA, indicating chains have not mixed.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#r-hat
## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#tail-ess
# Not run
# fit.sim = iCARH.model(XX, Y, Z, pathways, control = list(adapt_delta = 0.99, max_treedepth=10),
# iter = 2000, chains = 2, pars=c("YY","Xmis","Ymis"), include=F)
Check convergence
if(!is.null(fit.sim)){
rhats = iCARH.checkRhats(fit.sim)
}
## Warning in iCARH.checkRhats(fit.sim): Some Rhats are not < 1.1 ! The model has
## not converged! Results may not be valid!
Processing results. Bacteria effects.
if(!is.null(fit.sim)){
gplot = iCARH.plotBeta(fit.sim)
gplot
}
Treatments effects
if(!is.null(fit.sim)){
gplot = iCARH.plotTreatmentEffect(fit.sim)
gplot
}
Pathway analysis
if(!is.null(fit.sim)){
gplot = iCARH.plotPathwayPerturbation(fit.sim, path.names=names(data.sim$pathways))
gplot
}
Normality assumptions
if(!is.null(fit.sim)){
par(mfrow=c(1,2))
iCARH.checkNormality(fit.sim)
}
## , , 1
##
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] -3.055853 -3.3793145 -0.5286002 -1.39049058 -0.2794554 0.09269145
## [2,] 2.361803 -3.8433470 3.0848752 2.19363869 -2.0087848 1.37631605
## [3,] 2.278252 1.7680472 -2.3757978 -0.10974800 2.7299226 -5.89012830
## [4,] -3.680416 0.4470671 -3.3951388 -0.08894923 4.1024955 -3.45549081
## [,7] [,8] [,9] [,10]
## [1,] -2.21933667 1.8124772 1.1893381 -5.6098692
## [2,] 1.70995590 -0.7735045 -3.7636450 -7.5312222
## [3,] 4.71779073 0.2556656 0.7580608 -4.6209620
## [4,] -0.01425826 3.2195493 1.5398419 0.3794799
##
## , , 2
##
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] -0.3284334 -1.9540838 0.8247847 -0.02307358 1.814267 -0.72370586
## [2,] 3.1572937 -1.6099841 2.5858064 -0.08797553 2.133968 -0.13824912
## [3,] -1.2004942 -0.6951897 0.6836295 1.64590589 1.991602 -2.55640859
## [4,] -2.8735361 2.4936193 -0.3133960 1.34831580 1.627402 -0.01595692
## [,7] [,8] [,9] [,10]
## [1,] 0.5834214 -0.7467125 1.235927 -2.769734
## [2,] 1.3636518 -0.2330343 -2.093739 -3.444191
## [3,] 3.1397456 1.3426825 -3.628814 -4.167833
## [4,] 1.2399843 -0.1533040 -1.410770 -5.507195
##
## , , 3
##
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] -2.99179870 0.64263857 -0.4959946 1.1008134 -2.5639959 -3.4356580
## [2,] -0.05841406 -0.94975705 -0.4718328 -1.7518027 1.5125396 -0.8472650
## [3,] -1.23102514 0.03043126 1.4076742 0.1712632 1.5140814 -1.6556234
## [4,] -0.64810474 -0.98184616 -1.6586561 -0.3184771 0.6027405 0.3235535
## [,7] [,8] [,9] [,10]
## [1,] 1.7536733 1.3998068 -2.6058169 -1.6157539
## [2,] -0.3670534 0.1279264 -0.7609899 -3.1654945
## [3,] 2.3453030 -0.6442833 -2.8676734 -0.5861365
## [4,] 1.8554132 1.8938348 1.8330476 0.3404682
##
## , , 4
##
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] -2.639715 -0.272807 -0.7516999 2.2486232 -3.8385865 -2.657051 0.6526438
## [2,] 2.868181 -3.358887 -0.8492988 0.9049273 0.5681844 1.922930 -3.3795687
## [3,] 2.507224 0.525538 1.5661607 -1.6293610 0.3297031 -3.684169 2.9127362
## [4,] 1.495347 -3.142945 -4.3462281 -2.8936067 -0.4824623 1.870963 -2.0038077
## [,8] [,9] [,10]
## [1,] 0.7830353 -0.850093 -1.7966030
## [2,] -2.4492945 1.907672 -5.2937683
## [3,] -0.9249798 -2.157835 0.4809685
## [4,] 2.2763935 1.822896 0.7619912
##
## , , 5
##
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] -4.8535371 3.485492 -1.92237352 -0.653504 0.6946749 -7.112814 0.5202526
## [2,] -0.6759193 1.132748 -2.76937397 -2.557627 3.4236638 -6.476971 -2.3195725
## [3,] -1.7773845 1.849471 -0.05848731 -2.530423 4.6300235 -6.572182 0.5543444
## [4,] -2.3735927 1.279223 -4.11537429 -2.112815 2.2722127 -3.478474 -0.6405798
## [,8] [,9] [,10]
## [1,] 0.8431982 -4.005770 -0.5035193
## [2,] -1.1783780 -4.189976 -3.4143809
## [3,] 0.4387649 -7.828129 -1.4432192
## [4,] 2.3705797 -5.069171 -0.9991446
##
## , , 6
##
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 2.8131901 -1.7016565 -0.8069688 -2.3809484 0.1348349 -0.14956354
## [2,] 1.2039484 -0.3153951 -0.5434478 -1.2008039 0.6651418 -0.78570678
## [3,] -0.9281675 -1.3019125 -1.1034234 0.7221181 1.1649595 -0.00336566
## [4,] -0.1861503 -0.4087768 1.5102522 2.0038806 0.4711048 0.88340780
## [,7] [,8] [,9] [,10]
## [1,] 1.3477207 1.013277 2.4380218 0.5129650
## [2,] 0.5392288 2.894000 1.5753261 2.0453682
## [3,] 0.5778530 1.233550 0.4436044 -0.2319781
## [4,] 1.1756630 -1.355127 2.4854631 -1.5476348
##
## , , 7
##
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] -0.2886281 -2.1736983 2.0019729 -1.6660257 0.2657682 -1.1912396 1.675534
## [2,] 2.2237208 -1.0804494 2.9629182 1.2863938 0.6277873 -0.3845123 0.475804
## [3,] -0.5222361 -1.3013527 0.4073982 0.8154761 2.0971199 -2.4726658 2.905872
## [4,] 0.4513626 0.9492984 -0.1109712 0.9847377 1.2744248 0.1854789 1.053684
## [,8] [,9] [,10]
## [1,] 1.44909730 1.5038480 -1.0614343
## [2,] -0.18677738 -0.7497199 0.5469377
## [3,] 1.29259014 -0.5054480 -1.0147067
## [4,] 0.06309161 2.3659821 1.9335704
##
## , , 8
##
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] -0.6991714 1.8233268 0.1107428 1.953416 -1.826508 -4.18121370
## [2,] 3.7911242 0.2218168 -3.3722381 -2.745609 5.912001 -2.26647607
## [3,] -2.1835732 -7.1703241 5.2185615 1.745051 4.719060 -0.01418309
## [4,] 3.7206698 -5.9284613 5.6741455 3.814411 -2.594932 6.50764830
## [,7] [,8] [,9] [,10]
## [1,] 4.5679480 1.701571 2.582798 -0.2829431
## [2,] -1.1504802 4.113455 6.411411 -0.8810533
## [3,] -4.4214606 4.978817 -2.637754 -2.2513353
## [4,] -0.8745123 -1.634676 -0.304199 -7.0798690
##
## , , 9
##
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] -0.5132092 1.701684 -1.711147 -0.5537785 -1.023229 0.2670966 6.2822362
## [2,] 1.1150431 1.415893 -3.097090 -4.0976480 3.186414 1.2497587 -0.4353590
## [3,] -5.1363191 -5.912558 4.203246 3.8613854 -2.101172 1.8809731 -0.7745889
## [4,] 2.4410089 1.342428 2.483303 1.8800734 -2.836810 4.4435987 3.6296483
## [,8] [,9] [,10]
## [1,] 1.366799 3.714978 2.5299142
## [2,] 5.717048 1.968557 4.3530160
## [3,] 4.284836 -1.556802 0.2549909
## [4,] -6.018854 3.095187 -3.5715283
##
## , , 10
##
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 2.888437 -2.0700012 -0.83170224 -2.8843816 2.144631 3.7891380 0.8799502
## [2,] 2.411762 -1.0020884 -0.02639399 -0.8186279 -1.684528 2.1618507 4.4822369
## [3,] 1.248299 -0.9532143 -2.03803852 0.1578336 -2.542552 1.7070866 2.5132356
## [4,] -1.838578 1.7971372 0.45497134 2.5540840 0.673773 -0.2679568 3.0758991
## [,8] [,9] [,10]
## [1,] -0.69660642 3.4979302 -0.4955946
## [2,] 2.64647024 0.4767633 3.9824721
## [3,] 2.46764872 2.3735845 3.0693980
## [4,] -0.06643391 1.4485614 1.6721279
##
## , , 11
##
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 0.8591543 -4.0113612 -0.308759 -1.9910027 1.0984663 -0.57266697 -0.7603881
## [2,] 2.9637231 -2.5020090 1.258663 1.8416461 0.5492347 -0.08309946 -1.4673534
## [3,] 2.3405730 -0.3372740 -3.159389 -0.3718894 2.2736436 -3.03106118 1.5748163
## [4,] 0.8946844 0.2781008 -3.728561 -0.2527563 1.2674482 -1.32187130 -0.5770917
## [,8] [,9] [,10]
## [1,] 1.6899552 0.90609948 -2.608252073
## [2,] -1.7794114 -3.83683355 -0.005965877
## [3,] 0.9622921 -2.09286679 -1.964947871
## [4,] 0.5123571 0.03824114 1.294292112
##
## , , 12
##
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 1.6465032 -2.4618811 -0.9409611 -1.948637 0.6140882 5.7351772 2.684787
## [2,] -1.0092234 0.6047679 -2.5363562 -2.474505 -2.3243741 4.0815955 7.637091
## [3,] -0.7811522 -1.4780116 -3.3653845 -1.044052 -3.2434457 4.1680197 3.011250
## [4,] -2.8702234 1.9582153 1.5057414 1.717147 -0.5126841 0.6566439 3.745728
## [,8] [,9] [,10]
## [1,] 2.053220 4.328277 1.426316
## [2,] 5.048069 1.228534 5.860397
## [3,] 4.420899 5.692705 4.115184
## [4,] 2.218755 3.816032 1.345206
##
## , , 13
##
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 0.64790696 1.687257 -0.76737227 -1.2504416 1.4970839 1.7233941 1.608752
## [2,] -0.37912792 2.839368 0.07955434 0.1177429 -1.6761170 0.8425606 4.652170
## [3,] 0.02143126 2.572597 -1.36322560 -0.2490572 -0.6095207 1.3725989 4.136028
## [4,] -2.47097763 4.236956 -1.22220999 0.7351078 0.6898370 -0.5316373 3.682320
## [,8] [,9] [,10]
## [1,] -0.47263324 3.801549 -1.1415585
## [2,] 0.14798414 1.678268 -0.1236237
## [3,] -0.78995076 3.402221 0.1291980
## [4,] 0.03969458 1.255692 -0.4745612
##
## , , 14
##
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] -2.09501622 -0.01114191 0.7218449 0.3034668 -1.1761003 -4.685661
## [2,] -0.28932768 -1.98861065 -0.2170103 -0.6930691 2.0964527 -1.222784
## [3,] -0.23781697 -0.81190883 2.3882106 -0.1581033 0.5454562 -2.447034
## [4,] -0.02981654 -2.68713759 -1.3982798 -1.0323053 -1.4014882 -1.019845
## [,7] [,8] [,9] [,10]
## [1,] -0.3419983 -1.219465 -2.291252 -1.0927290
## [2,] -1.9052025 -1.792254 -1.651508 -0.3946558
## [3,] -0.0754639 -3.073194 -3.229984 -0.6625487
## [4,] -1.3401919 -1.844466 -1.042655 0.5533183
WAIC
if(!is.null(fit.sim)){
waic = iCARH.waic(fit.sim)
waic
}
## [1] 4834.561
Posterior predictive checks MAD : mean absolute deviation between covariance matrices
if(!is.null(fit.sim)){
mad = iCARH.mad(fit.sim)
quantile(mad)
}
## 0% 25% 50% 75% 100%
## 0.1758039 1.8310247 4.0298413 7.8455358 57.3236375
Get missing data
if(!is.null(fit.sim)){
imp = iCARH.getDataImputation(fit.sim)
}