We will start by loading the codalm R package, in addition to the ggtern package, which we will use to access the data.
library(codalm) library(ggtern) #> Loading required package: ggplot2 #> Registered S3 methods overwritten by 'ggtern': #> method from #> grid.draw.ggplot ggplot2 #> plot.ggplot ggplot2 #> print.ggplot ggplot2 #> -- #> Remember to cite, run citation(package = 'ggtern') for further info. #> -- #> #> Attaching package: 'ggtern' #> The following objects are masked from 'package:ggplot2': #> #> aes, annotate, ggplot, ggplot_build, ggplot_gtable, ggplotGrob, #> ggsave, layer_data, theme_bw, theme_classic, theme_dark, #> theme_gray, theme_light, theme_linedraw, theme_minimal, theme_void
We will now load in the data from the ggtern package. We will be analyzing how two different methods (image analysis or microscopic inspection) estimate the composition of 30 white blood cells. The format that we need is for both compositions to be in matrices, with one row per observation. We will also normalize the rows of these matrices to ensure that they sum to 1, although the codalm function would also take care of this for us.
data("WhiteCells", package = 'ggtern') image <- subset(WhiteCells, Experiment == "ImageAnalysis") image_mat <- as.matrix(image[,c("G", "L", "M")]) microscopic <- subset(WhiteCells, Experiment == "MicroscopicInspection") microscopic_mat <- as.matrix(microscopic[,c("G", "L", "M")]) image_mat <- image_mat / rowSums(image_mat) microscopic_mat <- microscopic_mat / rowSums(microscopic_mat)
To estimate the coefficient matrix B, we can use the codalm function.
B_est <- codalm(y = microscopic_mat, x = image_mat) B_est #> [,1] [,2] [,3] #> [1,] 9.742582e-01 0.02292755 2.814241e-03 #> [2,] 3.267605e-09 1.00000000 3.620972e-10 #> [3,] 9.484137e-18 0.04198815 9.580118e-01
To see the interpretation of this matrix, please see Fiksel et al. (2020). If all the rows of B_est are exactly the same, it is recommended to set accelerate = FALSE as a sensitivity check.
We can also use the bootstrap to estimate 95% confidence intervals. We will only use 50 bootstrap iterations as an example (nboot = 50), but is recommended to do more.
B_ci <- codalm_ci(y = microscopic_mat, x = image_mat, nboot = 50, conf = .95) B_ci$ci_L #> [,1] [,2] [,3] #> [1,] 9.576089e-01 0.01399186 1.156448e-03 #> [2,] 2.512029e-15 0.97220006 9.764192e-17 #> [3,] 6.197244e-21 0.01165546 9.145955e-01 B_ci$ci_U #> [,1] [,2] [,3] #> [1,] 9.837854e-01 0.03981767 4.560504e-03 #> [2,] 2.779953e-02 1.00000000 4.198096e-07 #> [3,] 3.143104e-06 0.08540455 9.868145e-01
These matrices given the lower and upper bounds for the confidence interval for each element of the coefficient matrix.
You can also take advantage of parallelization, if you have multiple cores available.
ncores <- 2 Sys.setenv(R_FUTURE_SUPPORTSMULTICORE_UNSTABLE = "quiet") B_ci_parallel <- codalm_ci(y = microscopic_mat, x = image_mat, nboot = 50, conf = .95, parallel = TRUE, ncpus = ncores, strategy = 'multisession') identical(B_ci$ci_L, B_ci_parallel$ci_L) #> [1] TRUE identical(B_ci$ci_U, B_ci_parallel$ci_U) #> [1] TRUE
Finally, we will do a permutation test for linear independence. Again, we will only do 50 permutations as an example, but in practice this number should be higher. For demonstration purposes, we will generate the compositional outcome independently of the compositional predictor
set.seed(123) x <- gtools::rdirichlet(100, rep(1, 3)) y <- gtools::rdirichlet(100, rep(1, 3)) indep_test_pval <- codalm_indep_test(y = y, x = x, nperms = 100, init.seed = 123) indep_test_pval #> [1] 0.28
This function can also be parallelized. Unlike the bootstrapping, there is no need to differentiate between whether the user is using a Windows or Unix system.
indep_test_pval_parallel <- codalm_indep_test(y = y, x = x, nperms = 100, init.seed = 123, parallel = TRUE, ncpus = ncores, strategy = 'multisession') indep_test_pval_parallel #> [1] 0.28