1
$\begingroup$

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 the next 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?

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      
$\endgroup$

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.