Plotting effects in Dirichlet regression
0
0
Entering edit mode
3 days ago
Juan • 0

Hi!

I performed different regression analyses, including Dirichlet regression, on compositional data. I want to create a plot similar to the one obtained when you use the package effects.

I have tried this piece of code:

sputum_cell_prop_col<-c("MS_DIFF_MAKROS", "MS_DIFF_NEUT", "MS_DIFF_EOS", "MS_DIFF_LYM",
"MS_DIFF_FLIMMEREPITHEL", "MS_DIFF_MONOS")

meta2$Y <- DR_data(mutate_all(meta2[sputum_cell_prop_col], function(x) as.numeric(as.character(x)))) res2<-DirichReg(Y ~ Eotaxin.3 + G.CSF + IFN.g + IL.10+ IL.13 + IL.17 + IL.1alpha + IL.1F7 + IL.24 + IL.33 + IL.4 + IL.5+ IL.8 + Periostin + SCGB1A.1 + TNF.alpha + D_SEX+D_AGE + D_SmokingStatus+year+OCS_YN + ICS_YN, meta2, control=list(iterlim = 50000)) new = data.frame( Eotaxin.3 = mean(meta2$Eotaxin.3),
G.CSF = seq(min(meta2$G.CSF) - 1,max(meta2$G.CSF) + 1,length.out=280),
IFN.g = mean(meta2$IFN.g), IL.10 = mean(meta2$IL.10),
IL.13 = mean(meta2$IL.13), IL.17 = mean(meta2$IL.17),
IL.1alpha = mean(meta2$IL.1alpha), IL.1F7 = mean(meta2$IL.1F7),
IL.24 = mean(meta2$IL.24), IL.33 = mean(meta2$IL.33),
IL.4 = mean(meta2$IL.4), IL.5 = mean(meta2$IL.5),
IL.8 = mean(meta2$IL.8), Periostin = mean(meta2$Periostin),
SCGB1A.1 = mean(meta2$SCGB1A.1), TNF.alpha = mean(meta2$TNF.alpha),
year = 2018,
D_SEX = meta2$D_SEX, D_AGE = mean(meta2$D_AGE),
D_SmokingStatus = meta2$D_SmokingStatus, OCS_YN = meta2$OCS_YN,
ICS_YN = meta2$ICS_YN ) new$year <- as.factor(new$year) new$D_SEX <- as.factor(new$D_SEX) new$D_SmokingStatus <- as.factor(new$D_SmokingStatus) new$OCS_YN <- as.factor(new$OCS_YN) new$ICS_YN <- as.factor(new\$ICS_YN)

a <- predict(res2,newdata = new)


This piece of code has the purpose of calculating the effects exclusively for the variable G.CSF. However, in this case, I should use the whole set of factor levels instead of fixing only one as the constant value.

Questions:

Is there a way to create a plot for specific covariate effects similar to the plots provided by the function Effects when I use other regression models like linear or multimodal? If so, is the previous piece of code suitable for this task?

Kind regards,

Juan

Session info

R version 4.4.1 (2024-06-14 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 22631)

Matrix products: default

locale:
[1] LC_COLLATE=German_Germany.utf8  LC_CTYPE=German_Germany.utf8
[3] LC_MONETARY=German_Germany.utf8 LC_NUMERIC=C
[5] LC_TIME=German_Germany.utf8

time zone: Europe/Berlin
tzcode source: internal

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods
[8] base

other attached packages:
[1] betareg_3.2-0         ggsignif_0.6.4        sciplot_1.2-0
[4] gridExtra_2.3         DescTools_0.99.55     MASS_7.3-61
[7] AER_1.2-12            survival_3.7-0        sandwich_3.1-0
[10] lmtest_0.9-40         zoo_1.8-12            car_3.1-2
[13] effects_4.2-2         carData_3.0-5         ordinal_2023.12-4.1
[16] Rtsne_0.17            M3C_1.26.0            ComplexHeatmap_2.20.0
[19] missRanger_2.6.0      limma_3.60.4          openxlsx_4.2.6.1
[22] RColorBrewer_1.1-3    ggplot2_3.5.1         dplyr_1.1.4
[25] DirichletReg_0.7-1    Formula_1.2-5

loaded via a namespace (and not attached):
[4] rlang_1.1.4         magrittr_2.0.3      clue_0.3-65
[7] GetoptLong_1.0.5    e1071_1.7-14        matrixStats_1.3.0
[10] compiler_4.4.1      flexmix_2.3-19      png_0.1-8
[13] vctrs_0.6.5         pkgconfig_2.0.3     shape_1.4.6.1
[16] crayon_1.5.3        fastmap_1.2.0       labeling_0.4.3
[19] utf8_1.2.4          rmarkdown_2.28      nloptr_2.1.1
[22] miscTools_0.6-28    modeltools_0.2-23   xfun_0.47
[25] jsonlite_1.8.8      parallel_4.4.1      cluster_2.1.6
[28] R6_2.5.1            stringi_1.8.4       reticulate_1.38.0
[31] boot_1.3-30         estimability_1.5.1  cellranger_1.1.0
[34] numDeriv_2016.8-1.1 Rcpp_1.0.13         iterators_1.0.14
[37] knitr_1.48          snow_0.4-4          IRanges_2.38.1
[40] Matrix_1.7-0        splines_4.4.1       nnet_7.3-19
[43] tidyselect_1.2.1    rstudioapi_0.16.0   abind_1.4-5
[46] yaml_2.3.10         doParallel_1.0.17   maxLik_1.5-2.1
[49] codetools_0.2-20    lattice_0.22-6      tibble_3.2.1
[55] proxy_0.4-27        survey_4.4-2        zip_2.3.1
[58] circlize_0.4.16     pillar_1.9.0        foreach_1.5.2
[61] stats4_4.4.1        insight_0.20.3      generics_0.1.3
[64] S4Vectors_0.42.1    rootSolve_1.8.2.4   munsell_0.5.1
[67] scales_1.3.0        minqa_1.2.8         class_7.3-22
[70] glue_1.7.0          lmom_3.0            tools_4.4.1
[73] data.table_1.15.4   lme4_1.1-35.5       RSpectra_0.16-2
[76] Exact_3.3           mvtnorm_1.2-6       mitools_2.4
[79] matrixcalc_1.0-6    umap_0.2.10.0       colorspace_2.1-1
[82] nlme_3.1-166        cli_3.6.3           expm_1.0-0
[85] fansi_1.0.6         corpcor_1.6.10      doSNOW_1.0.20
[88] gtable_0.3.5        digest_0.6.37       BiocGenerics_0.50.0
[91] ucminf_1.2.2        farver_2.1.2        rjson_0.2.22
[94] htmlwidgets_1.6.4   htmltools_0.5.8.1   lifecycle_1.0.4
[97] httr_1.4.7          GlobalOptions_0.1.2 statmod_1.5.0
[100] openssl_2.2.1

R compositional-data dirichlet-regression • 164 views