CARlasso is the main interface to work with the model, for details see the reference.
This is the case when data is Gaussian or can be transformed to Gaussian, for example, (log) biomass of trees. We will provide an example with simulated data. The use of the
Carlasso is similar to
lm, we could supply a formula and a dataframe with both responses and predictors.
First, we simulate data under a 5-node AR1 model:
set.seed(42) <- simu_AR1(n=100,k=5, rho=0.7) dt head(dt) #> y1 y2 y3 y4 y5 x1 x2 #> 1 1.2705558 0.94046103 -0.9356549 -0.5202384 -0.1928792 1.3709584 1.2009654 #> 2 0.6027324 0.06671259 0.9096216 1.3740787 0.9584424 -0.5646982 1.0447511 #> 3 -0.1139447 -1.61460997 -0.1140354 -0.8013874 0.4673093 0.3631284 -1.0032086 #> 4 0.5912523 2.13482891 2.2626483 1.8791249 2.3825324 0.6328626 1.8484819 #> 5 -2.1837420 -3.68267385 -2.3069358 -3.6243897 -2.2644276 0.4042683 -0.6667734 #> 6 -2.2367855 -1.92498967 -3.2242525 -1.5495283 -1.4191343 -0.1061245 0.1055138 #> x3 x4 x5 #> 1 -2.0009292 -0.004620768 1.33491259 #> 2 0.3337772 0.760242168 -0.86927176 #> 3 1.1713251 0.038990913 0.05548695 #> 4 2.0595392 0.735072142 0.04906691 #> 5 -1.3768616 -0.146472627 -0.57835573 #> 6 -1.1508556 -0.057887335 -0.99873866
To use the Normal version, we should set
link="identity" which is the default. In this case, we are setting
adaptive=TRUE to use the adaptive version of CAR-LASSO (for more details, see the paper):
<- CARlasso(y1+y2+y3+y4+y5~x1+x2+x3+x4+x5, data = dt, adaptive = TRUE) car_res plot(car_res,tol = 0.05)
The color of the edge represents the type of correlation (negative=blue, positive=red) and the width of the edge corresponds to the effect size. Response nodes are represented by circles and predictor nodes are represented by triangles.
We can have a more formal horseshoe inference on the structure of the network which will update the
# with horseshoe inference <- horseshoe(car_res) car_res plot(car_res)
This is common in the case of microbe-related studies and some ecological applications with relative abundances. For instance, microbe relative abundance come from sequencing and in this case, the sum of “abundance” is determined by the sequence depth rather than the real total abundance. The data is usually described as “compositional”. In
CARlasso, this type of data are modeled as Logit-Normal-multinomial. In this case, we need to have a “reference level” taxa and all others are “relative” to this taxa.
First, we take a look at the data which is still a dataframe with all predictors and responses
1:5,1:7] mgp154[#> id Alistipes Bacteroides Barnesiella Blautia Butyrivibrio Caloramator #> 1 EM245 42 2101 0 334 61 156 #> 2 EM109 2910 6080 40 597 178 175 #> 3 EM248 2850 5392 0 327 46 87 #> 4 EM274 1559 7210 43 1206 50 428 #> 5 EM330 1916 5067 0 1827 64 5847
To run the composition model, we need to set
<- CARlasso(Alistipes+Bacteroides+ gut_res +Parabacteroides+all_others~ Eubacterium+Age+Gender+Stratum, BMIdata = mgp154,link = "logit", adaptive = TRUE, n_iter = 2000, n_burn_in = 1000, thin_by = 2)
Note that in this case
all_others (an existing column in our data), i.e. the last one in the left hand side of the formula will be the reference level.
We can update the network inference by a horseshoe method to determine when edges will be considered non-existent. More details can be found in the paper.
# horseshoe will take a while, as it needs to sample the latent normal too <- horseshoe(gut_res) gut_res plot(gut_res)
This is common in a lot of ecological applications. For instance, number of seedlings within a site. The responses are counts rather than continuous. In
CARlasso, it is modeled as Poisson with log-Normal rate.
We will use the same compositional data as before to illustrate the counts model. However, it is important to note that relative abundances should not be considered as counts. To distinguish between compositional and count data, one can ask the question: is the sum decided by us? If yes, we want to use compositional models.
To run the count model, we need to set
<- CARlasso(Alistipes+Bacteroides+ gut_res +Parabacteroides+all_others~ Eubacterium+Age+Gender+Stratum, BMIdata = mgp154,link = "log", adaptive = TRUE, r_beta = 0.1, # default sometimes cause singularity in Poisson model due to exponential transformation, slightly change can fix it. n_iter = 2000, n_burn_in = 1000, thin_by = 2) # horseshoe will take a while, as it's currently implemented in R rather than C++ <- horseshoe(gut_res) gut_res plot(gut_res)