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  2 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

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] 86
## [1] 90
## [1] 122
## [1] 118

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)
## DIAGNOSTIC(S) FROM PARSER:
## Info:
## Left-hand side of sampling statement (~) may contain a non-linear transform of a parameter or local variable.
## If it does, you need to include a target += statement with the log absolute determinant of the Jacobian of the transform.
## Left-hand-side of sampling statement:
##     YY[1, :, k] ~ multi_normal(...)
## Info:
## Left-hand side of sampling statement (~) may contain a non-linear transform of a parameter or local variable.
## If it does, you need to include a target += statement with the log absolute determinant of the Jacobian of the transform.
## Left-hand-side of sampling statement:
##     YY[i, :, k] ~ multi_normal(...)
## 
## 
## SAMPLING FOR MODEL 'iCARH' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.001959 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 19.59 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.009264 seconds (Warm-up)
## Chain 1:                0.003827 seconds (Sampling)
## Chain 1:                0.013091 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 > 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
}

plot of chunk unnamed-chunk-4 Treatments effects

if(!is.null(fit.sim)){
  gplot = iCARH.plotTreatmentEffect(fit.sim)
  gplot
}

plot of chunk unnamed-chunk-5

Pathway analysis

if(!is.null(fit.sim)){
  gplot = iCARH.plotPathwayPerturbation(fit.sim, path.names=names(data.sim$pathways))
  gplot
}

plot of chunk unnamed-chunk-6 Normality assumptions

if(!is.null(fit.sim)){
  par(mfrow=c(1,2))
  iCARH.checkNormality(fit.sim)
}

plot of chunk unnamed-chunk-7

## , , 1
## 
##            [,1]      [,2]      [,3]       [,4]       [,5]      [,6]      [,7]
## [1,] -10.213576  6.445569 -1.321549  7.8361808 -10.524550 -5.138802 10.053276
## [2,]   2.387296 -3.699570 -1.247321 -0.2126838   4.184093 15.116497 -6.919978
## [3,]   1.426918  4.176082 10.230606  5.0129742  -1.016624  3.161045 16.384401
## [4,]   7.630790 -1.800949 -9.513854 -8.5968441  -4.562519 16.110259  7.103706
##           [,8]      [,9]     [,10]
## [1,] 13.375901  1.524589 10.161447
## [2,]  1.455638 14.052005 -3.265036
## [3,]  3.491758 -5.201962  7.262404
## [4,] 11.728180 18.939164 17.182640
## 
## , , 2
## 
##            [,1]       [,2]      [,3]       [,4]       [,5]       [,6]
## [1,]  -5.088844  11.296419 -5.753211   7.306866  -8.057794 -12.937826
## [2,]   5.316184   1.607084 -7.177147 -13.198891  13.983388   9.955749
## [3,] -14.693839 -12.952884 20.013231  18.689553  -6.759622   5.167952
## [4,]  19.035867   8.210732 -6.406639  -6.896277 -13.193537  19.125974
##            [,7]      [,8]       [,9]      [,10]
## [1,]  13.917539  5.132316   1.704708  9.8751136
## [2,] -14.770455  5.533125  15.619601  3.0008821
## [3,]   5.045686  6.111875 -19.209248 -0.1314982
## [4,]  11.099443 -8.553219  16.421326 -2.3303916
## 
## , , 3
## 
##            [,1]      [,2]      [,3]       [,4]      [,5]      [,6]      [,7]
## [1,] -9.9504584  5.536070  2.235182  8.2883930 -6.992285 -1.310369  4.521084
## [2,] -0.7212141 -2.868168  1.646839  3.8826050  4.581448  8.959376 -2.013617
## [3,]  5.0468391  9.982679  7.218313  0.5175146  3.275421  3.369841 13.229531
## [4,]  2.1886473 -2.906611 -7.998087 -6.3579664  1.106099  9.311516  3.403617
##            [,8]        [,9]     [,10]
## [1,]  9.6654524 -1.94796676  7.460431
## [2,] -1.5987622  7.47553158 -2.974674
## [3,] -0.7280409  0.04128297  8.264015
## [4,] 15.0237564 10.79428097 17.086543
## 
## , , 4
## 
##            [,1]       [,2]      [,3]       [,4]      [,5]      [,6]       [,7]
## [1,] -5.6735893  2.9297849  1.817473  4.1115929 -6.349314 0.7618296  0.3147593
## [2,]  0.3899674 -1.9534183  2.490417  4.1720778 -2.736602 6.1923921 -2.4976122
## [3,]  4.5171129  6.5835503  4.037349 -0.9771058 -2.265911 0.7375205  7.2116565
## [4,]  0.1315630 -0.1985016 -5.501402 -4.4908088 -2.152796 5.0778110 -0.1693815
##             [,8]       [,9]     [,10]
## [1,]  4.94891090 -0.6700991  4.919355
## [2,] -1.55232897  3.3241391 -1.882332
## [3,]  0.03592352 -0.3325591  5.797383
## [4,]  9.10582993  4.3809158 10.193718
## 
## , , 5
## 
##           [,1]      [,2]      [,3]       [,4]       [,5]       [,6]      [,7]
## [1,] -3.984139  1.973527 -0.835078  3.1263158 -4.3867157 -1.6508347  2.256336
## [2,]  2.151434 -3.365746 -1.816938  0.3058599  0.4070566  2.1516427 -1.438986
## [3,]  3.616026  1.726672  2.356753 -1.3642883  0.5639876  0.1665699  5.190923
## [4,]  1.946208 -3.608431 -5.521780 -3.1060434 -1.5431763  3.6159837  0.783973
##            [,8]       [,9]     [,10]
## [1,]  4.5513520  0.7152833  4.074731
## [2,] -0.1702616  3.7426897 -1.132238
## [3,]  1.4680056 -0.7773321  4.121814
## [4,]  7.6920099  3.8710559  6.808840
## 
## , , 6
## 
##            [,1]      [,2]      [,3]       [,4]       [,5]       [,6]       [,7]
## [1,] -2.7383169  1.349171 -1.715640 -0.1733047 -5.4515381 0.00459552  2.2114041
## [2,]  0.9218435 -2.747078 -1.389818 -0.6874606  2.0572222 4.81556382 -2.3090983
## [3,]  2.0550494  4.407013  1.581830 -1.2781606 -0.1060202 1.42493353  6.7654836
## [4,]  1.4458444 -3.156923 -5.439342 -4.6184580 -1.3965580 5.27352367  0.4298045
##            [,8]       [,9]     [,10]
## [1,]  4.3998298  1.7325981  2.526150
## [2,] -0.5220565  5.6440897 -2.238925
## [3,] -0.2004671 -0.1082785  2.285320
## [4,]  5.3376058  7.5388413  5.494403
## 
## , , 7
## 
##            [,1]       [,2]       [,3]       [,4]       [,5]       [,6]
## [1,] -0.9304713  0.7619735  2.0733312 -0.5172805 -2.4249085 -0.5914261
## [2,]  1.7818442 -0.4231280  2.7001516  0.2930573 -2.3805807  0.6036866
## [3,]  2.2142434  2.5867173  1.5698787 -1.8561728 -1.1519809 -1.9932881
## [4,]  0.5245623  1.0687898 -0.7861162 -1.8544903 -0.7860856 -1.3389662
##            [,7]     [,8]      [,9]     [,10]
## [1,]  0.2044529 5.024195 1.1285320 0.1547737
## [2,] -0.3243634 3.411138 0.2144481 0.3852090
## [3,]  2.6710771 3.242820 1.3903077 0.9014257
## [4,] -0.1955173 4.874748 3.4419574 4.1483907
## 
## , , 8
## 
##            [,1]      [,2]       [,3]      [,4]       [,5]       [,6]      [,7]
## [1,] -4.1898564  2.152373  0.4204488  3.197670 -3.2849761 -0.6690823 4.7768449
## [2,]  1.2924799 -2.552609 -2.1316888  0.553516  4.1197509  3.7893594 0.9626007
## [3,]  0.3673145  0.456389  3.9086604  0.441065  1.8028190  1.6667007 5.5568017
## [4,]  3.3202784 -3.100891 -2.9153943 -1.892768 -0.6317664  5.9742256 3.0484901
##          [,8]       [,9]      [,10]
## [1,] 4.925343  2.3308563 4.11507373
## [2,] 2.177393  5.6761545 0.09809086
## [3,] 1.947129 -0.5151416 4.15243217
## [4,] 5.296060  5.3230370 5.70659568
## 
## , , 9
## 
##           [,1]      [,2]       [,3]      [,4]       [,5]      [,6]      [,7]
## [1,] -3.684052 -1.231501  0.2465122  1.151183 -2.3853800 1.7189077 0.3431556
## [2,] -1.681838 -3.351753  0.8912591  2.953810 -1.5102622 3.6270555 0.2309453
## [3,]  3.907627  3.158380  0.1239338 -1.523693 -0.6546571 0.8758389 5.7957651
## [4,] -1.364651 -2.594029 -3.3797259 -3.426324  0.9075057 1.1567036 0.5096630
##            [,8]      [,9]     [,10]
## [1,]  2.7400698 0.2869390  1.017397
## [2,] -0.8497008 0.2328078 -2.168019
## [3,] -0.6964438 3.2930774  1.135869
## [4,]  4.9903251 4.3043737  6.523448
## 
## , , 10
## 
##            [,1]      [,2]       [,3]       [,4]       [,5]      [,6]
## [1,] -12.927259  5.610130  0.1690126 11.1387697 -10.325396 -1.974889
## [2,]  -0.938883 -5.558102 -3.8215340  4.8683248   7.853790 11.787889
## [3,]   4.526204  6.649827  9.0366584  0.8475988   5.474417  6.756651
## [4,]   6.364211 -8.385255 -9.2160464 -9.1362982  -1.405948 16.339015
##             [,7]      [,8]      [,9]     [,10]
## [1,] 12.12481699 13.819728  1.400411  7.769596
## [2,]  0.09014148  3.043734 14.612502 -1.772361
## [3,] 15.78398463  2.988015  2.604647  7.772622
## [4,]  7.64536299 16.475761 14.402570 18.138213
## 
## , , 11
## 
##            [,1]       [,2]        [,3]      [,4]        [,5]      [,6]
## [1,] -0.2222178 -0.5631642 -0.41554197 0.3512547  0.39245499 1.0609386
## [2,] -0.1054625 -0.1384337 -0.05688655 0.4973455  0.57020966 1.2432559
## [3,]  0.2202642  0.5459460 -0.39119208 0.2880131 -0.08071313 0.2296872
## [4,] -0.5474087  0.2446263 -0.46352358 0.3429322 -0.38170819 0.1078989
##            [,7]      [,8]       [,9]     [,10]
## [1,] 0.51702245 1.5023415  1.8506450 0.3026201
## [2,] 0.08394884 0.5499894 -0.4490961 1.7217244
## [3,] 0.66197486 1.3940372  0.4154258 0.8776565
## [4,] 0.52452508 0.4610118  0.7739062 0.9082293
## 
## , , 12
## 
##            [,1]       [,2]        [,3]       [,4]       [,5]       [,6]
## [1,] -1.4767296 -1.9133019  0.84907171  1.2575403 -1.3145768  0.8741727
## [2,] -0.4942174 -2.0911269 -0.20228412  1.2994867 -0.7977411  1.4743672
## [3,]  3.6796563  1.9287542 -0.08235122 -0.6950127 -1.0366625 -0.1160956
## [4,] -0.7781045 -0.2432105 -1.95937220 -1.3454678  0.8171749 -1.0356078
##            [,7]       [,8]       [,9]      [,10]
## [1,] -1.0331139  0.8899643 -1.1700695 -0.5002178
## [2,]  0.7054288 -2.3175828 -1.7026772 -1.9161104
## [3,]  2.4941362 -2.1573259  0.4764347  0.1036484
## [4,] -0.5914276  2.7837869  0.8318956  2.4889030
## 
## , , 13
## 
##           [,1]      [,2]       [,3]       [,4]       [,5]      [,6]       [,7]
## [1,] -8.241621  3.922680  2.1529961  7.9321801 -6.4601546 0.0112584  6.4937715
## [2,]  2.234654 -4.782778  0.1160969  2.5358499  5.7681792 9.3824649 -1.4259504
## [3,]  4.835346  8.440078  6.7862946 -0.8044641  1.3768450 2.0525969 13.3341220
## [4,]  2.217657 -6.329729 -9.2774488 -6.5376260 -0.8984314 8.7524607  0.8877577
##            [,8]      [,9]     [,10]
## [1,]  9.0672784 1.4711123  6.189653
## [2,] -3.0813288 8.8085614 -3.771866
## [3,]  0.2648297 0.0411132  7.491355
## [4,] 11.9946514 9.3088694 12.385778
## 
## , , 14
## 
##            [,1]      [,2]      [,3]       [,4]       [,5]      [,6]       [,7]
## [1,] -0.6514879 -3.940836  3.921054 -0.9067297  0.8825901  6.738020 -6.8910077
## [2,] -2.8835196 -6.323699  8.268292  9.2942983 -7.0380196  5.957326  2.3629362
## [3,]  9.4333325  8.268532 -3.709367 -2.7268288 -2.1554031  1.402208  8.1322588
## [4,] -1.6403038  2.920595 -6.660719 -4.8565572  6.6357988 -6.430271  0.5170961
##           [,8]       [,9]      [,10]
## [1,]  1.138751  0.3137539 -3.0128319
## [2,] -1.978896 -6.0642900 -5.3957343
## [3,] -6.620977  6.8554359 -0.9983167
## [4,]  7.352914  4.1839656 10.7302697

WAIC

if(!is.null(fit.sim)){
  waic = iCARH.waic(fit.sim)
  waic
}
## [1] Inf

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.06433024   1.88518671   4.88946693  16.95681449 618.64202615

Get missing data

if(!is.null(fit.sim)){
  imp = iCARH.getDataImputation(fit.sim)
}