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?

Thanks in advance,

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):
  [1] DBI_1.2.3           gld_2.6.6           readxl_1.4.3
  [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
 [52] withr_3.0.1         askpass_1.2.0       evaluate_0.24.0
 [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