Add q-values

[ataquette]

Hello!

I love the modelsummary table and would really love to use it but I need to add q-values (adjusted p-values). I am not sure how to do this and if the tidy.custom function would allow that. I have posted my question here

https://stackoverflow.com/questions/68648459/how-to-add-a-q-value-adjusted-p-value-to-a-modelsummary-table-after-pooling-th

I think it would really help to have an inbuilt function for that Thank you

[vincentarelbundock]

Hi,

I’m on vacation now without a computer, so I can’t give you working code, but it should be very easy to add that. In principle all you need is a function that returns a data.frame with two columns. For example:

tidy_custom.mipo = function(model, …) {
out = data.frame(
   term = names(coef(model))
   p.value = your_adjustment_function(model))
return(out)
}

Then all models should be adjusted automatically.

[vincentarelbundock]

@ataquette were you eventually able to get this working? I'm back so could take a look now. However, your SO question does not give me enough info, since the example is not fully reproducible (in particular, I don't have exposures objects).

[ataquette]

Hello!

Thanks for getting back to me and sorry to have bothered you on your holidays! Ok so no, not fully, because I work with lists...and I am slow!!

Here is the code - would be awesome if you get it to work!!

Thanks for much

# Setup

## packages
library(Hmisc)
library(pbapply)
library(mice)
library(data.table)
library(modelsummary)
library(nnet)
library(mlogit)

##  set global R options
pboptions(type = "timer", char = "=") # initialize progress bar

# Univariate Models 

##  set list of exposures to test
exposures <- c(
    Cs(
        age,
        ethnic_group,
        med_lastyr_binary_num,
        vit_lastyr_binary_num,
        maritalstatus,
        gender,
        bmi_group,
        physact_before,
        income_before,
        smoke_before,
        contra_current_groupall,
        cyclelength_before_group,
        period_lengh_before_group,
        heavyperiod_before,
        cycle_irreg_before,
        nbdeliveries,
        vaccine_type,
        vaccine_timing,
        covid_group,
        covid_tested,
        vaccinated,
        changesmenses_num,
        lifesatchanges,
        disease_before_endometriosis_num,
        disease_before_pcos_num,
        disease_before_underthyroid_num,
        disease_before_overthyroid_num,
        disease_before_uterpolyp_num,
        disease_before_uterfibroid_num,
        disease_before_intercystitis_num,
        disease_before_eatingdisord_num

    )
)

## model function
models <- function(x) {
    lapply(multfin, function(y)
        multinom(as.formula(
            paste0(
                "vaccine_changescycle ~ ",
                x
            )
        ), data = y, model = TRUE)
    )
}

## run XWAS models
models_univariate <- as.list(seq(1,length(exposures))) # create list to store model results
models_univariate <- pblapply(exposures, models)

# Pool 

pool_univariate <- as.list(seq(1,length(exposures))) # create list to store pool results

## run pool
for(j in seq_along(exposures)) {
    pool_univariate[[j]] <- pool(models_univariate[[j]])
}

# Q values 

library(data.table)
smry <-lapply(pool_univariate,
             summary,
              conf.int = TRUE,
            conf.level = 0.95)
all_models <- rbindlist(smry)

all_models$q.value <- p.adjust(all_models$p.value,method = "BH")

# Table

msuni<-modelsummary(pool_univariate,
                   group = model + term ~ y.level + y.level,
                    stars = c('*' = 0.05,'**' = 0.01,'***' = 0.001),
                    exponentiate = TRUE,estimate = "{estimate}{stars} [{conf.low}, {conf.high}]", 
                    statistic = NULL ,
                    gof_omit = ".*",title = 'Risk factors for vaccine related changes in menstrual cycles',
                    coef_omit = "Intercept", 
                    output = "kableExtra") 

[vincentarelbundock]

Still not reproducible:

models_univariate <- pblapply(exposures, models)
  |                                                  | 0 % ~calculating  Error in lapply(multfin, function(y) multinom(as.formula(paste0("vaccine_changescycle ~ ",  : 
  object 'multfin' not found

[ataquette]

Yes ok multfin is a list of 5 imputed datasets. Do you need the .Rdata file?

[vincentarelbundock]

Well, if you don't have a choice, the .Rdata would do, but really, the best approach would be to compose a truly minimal and self-contained example. You're also much more likely to get help on Stack Overflow if you supply a minimal reproducible example.

I found these two links very useful when learning how to create MWEs:

https://reprex.tidyverse.org/articles/reprex-dos-and-donts.html

https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example

The reprex package can be especially useful for this.

[ataquette]

Thanks so much - I will try that this evening - thanks!

[vincentarelbundock]

For a reason unrelated to this particular issue, I wrote a short case study showing how to adjust p values of a linear regression model. Posting here in case it’s useful.

[EDIT: Also copied a similar thing to StackOverflow in an answer]

The mechanism described in the previous section can also be used to add custom statistics to a table. For example, many researchers want to adjust the p values that they report to account for multiple comparisons. The p.adjust function in can calculate many types of corrected p values, adjusted following the methods of Bonferroni, Holm, and many others.

To adjust the p values of a linear regression model, we define a new tidy_custom.lm method which returns a data frame with one column called term and another column with the new statistic we wish to report. We can also add statistics to the bottom of the table by defining a glance_custom.lm method which returns a data frame with one row and one piece of information per column. For instance, if the analyst plans to conduct 10 tests with the drat coefficient, they could write:

library(modelsummary)
mod <- lm(mpg ~ hp + drat + vs + am, data = mtcars)

tidy_custom.lm <- function(x, ...) {
  out <- broom::tidy(x)
  out$bonferroni <- p.adjust(out$p.value, n = 10, method = "bonferroni")
  out$holm <- p.adjust(out$p.value, n = 10, method = "holm")
  return(out)
}

glance_custom.lm <- function(x, ...) {
  out <- data.frame(
    "Num.Comparisons" = "10",
    "Model" = class(x)[1]
  )
  return(out)
}

Then, we call modelsummary and use strings in the statistic argument to label the different p values. This code produces Table :

modelsummary(mod,
  statistic = c(
    "p = {p.value}",
    "p (Bonferroni) = {bonferroni}",
    "p (Holm) = {holm}"
  ),
  gof_omit = ".*",
  fmt = 5
)

	Model 1
(Intercept)	19.16936
	p = 0.00188
	p (Bonferroni) = 0.01875
	p (Holm) = 0.01688
hp	−0.04365
	p = 0.00045
	p (Bonferroni) = 0.00445
	p (Holm) = 0.00445
drat	1.25285
	p = 0.42027
	p (Bonferroni) = 1.00000
	p (Holm) = 1.00000
vs	2.32286
	p = 0.13511
	p (Bonferroni) = 1.00000
	p (Holm) = 0.94579
am	4.43462
	p = 0.00593
	p (Bonferroni) = 0.05926
	p (Holm) = 0.04741

Session info

``` r sessioninfo::session_info() #> ─ Session info ─────────────────────────────────────────────────────────────── #> setting value #> version R version 4.1.0 (2021-05-18) #> os macOS Big Sur 10.16 #> system x86_64, darwin17.0 #> ui X11 #> language (EN) #> collate en_CA.UTF-8 #> ctype en_CA.UTF-8 #> tz America/Toronto #> date 2021-08-09 #> #> ─ Packages ─────────────────────────────────────────────────────────────────── #> package * version date lib #> assertthat 0.2.1 2019-03-21 [1] #> backports 1.2.1 2020-12-09 [1] #> broom 0.7.9 2021-07-27 [1] #> checkmate 2.0.0 2020-02-06 [1] #> cli 3.0.1 2021-07-17 [1] #> colorspace 2.0-2 2021-06-24 [1] #> crayon 1.4.1 2021-02-08 [1] #> DBI 1.1.1 2021-01-15 [1] #> digest 0.6.27 2020-10-24 [1] #> dplyr 1.0.7 2021-06-18 [1] #> ellipsis 0.3.2 2021-04-29 [1] #> evaluate 0.14 2019-05-28 [1] #> fansi 0.5.0 2021-05-25 [1] #> fs 1.5.0 2020-07-31 [1] #> generics 0.1.0 2020-10-31 [1] #> glue 1.4.2 2020-08-27 [1] #> htmltools 0.5.1.1 2021-01-22 [1] #> httr 1.4.2 2020-07-20 [1] #> kableExtra 1.3.4 2021-02-20 [1] #> knitr 1.33 2021-04-24 [1] #> languageserver * 0.3.10 2021-04-20 [1] #> lifecycle 1.0.0 2021-02-15 [1] #> magrittr 2.0.1 2020-11-17 [1] #> modelsummary * 0.8.1.9000 2021-08-09 [1] #> munsell 0.5.0 2018-06-12 [1] #> pillar 1.6.2 2021-07-29 [1] #> pkgconfig 2.0.3 2019-09-22 [1] #> purrr 0.3.4 2020-04-17 [1] #> R.cache 0.15.0 2021-04-30 [1] #> R.methodsS3 1.8.1 2020-08-26 [1] #> R.oo 1.24.0 2020-08-26 [1] #> R.utils 2.10.1 2020-08-26 [1] #> R6 2.5.0 2020-10-28 [1] #> reprex 2.0.0 2021-04-02 [1] #> rlang 0.4.11 2021-04-30 [1] #> rmarkdown 2.10 2021-08-06 [1] #> rstudioapi 0.13 2020-11-12 [1] #> rvest 1.0.1 2021-07-26 [1] #> scales 1.1.1 2020-05-11 [1] #> sessioninfo 1.1.1 2018-11-05 [1] #> stringi 1.7.3 2021-07-16 [1] #> stringr 1.4.0 2019-02-10 [1] #> styler 1.4.1 2021-03-30 [1] #> svglite 2.0.0 2021-02-20 [1] #> systemfonts 1.0.2 2021-05-11 [1] #> tables 0.9.6 2020-09-22 [1] #> tibble 3.1.3 2021-07-23 [1] #> tidyr 1.1.3 2021-03-03 [1] #> tidyselect 1.1.1 2021-04-30 [1] #> utf8 1.2.2 2021-07-24 [1] #> vctrs 0.3.8 2021-04-29 [1] #> viridisLite 0.4.0 2021-04-13 [1] #> webshot 0.5.2 2019-11-22 [1] #> withr 2.4.2 2021-04-18 [1] #> xfun 0.25 2021-08-06 [1] #> xml2 1.3.2 2020-04-23 [1] #> yaml 2.2.1 2020-02-01 [1] #> source #> standard (@0.2.1) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> standard (@1.1.1) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> Github (vincentarelbundock/modelsummary@4467537) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> standard (@2.0.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> #> [1] /Library/Frameworks/R.framework/Versions/4.1/Resources/library ```

Session info

``` r sessioninfo::session_info() #> ─ Session info ─────────────────────────────────────────────────────────────── #> setting value #> version R version 4.1.0 (2021-05-18) #> os macOS Big Sur 10.16 #> system x86_64, darwin17.0 #> ui X11 #> language (EN) #> collate en_CA.UTF-8 #> ctype en_CA.UTF-8 #> tz America/Toronto #> date 2021-08-09 #> #> ─ Packages ─────────────────────────────────────────────────────────────────── #> package * version date lib #> assertthat 0.2.1 2019-03-21 [1] #> backports 1.2.1 2020-12-09 [1] #> broom 0.7.9 2021-07-27 [1] #> checkmate 2.0.0 2020-02-06 [1] #> cli 3.0.1 2021-07-17 [1] #> colorspace 2.0-2 2021-06-24 [1] #> crayon 1.4.1 2021-02-08 [1] #> DBI 1.1.1 2021-01-15 [1] #> digest 0.6.27 2020-10-24 [1] #> dplyr 1.0.7 2021-06-18 [1] #> ellipsis 0.3.2 2021-04-29 [1] #> evaluate 0.14 2019-05-28 [1] #> fansi 0.5.0 2021-05-25 [1] #> fs 1.5.0 2020-07-31 [1] #> generics 0.1.0 2020-10-31 [1] #> glue 1.4.2 2020-08-27 [1] #> htmltools 0.5.1.1 2021-01-22 [1] #> httr 1.4.2 2020-07-20 [1] #> kableExtra 1.3.4 2021-02-20 [1] #> knitr 1.33 2021-04-24 [1] #> languageserver * 0.3.10 2021-04-20 [1] #> lifecycle 1.0.0 2021-02-15 [1] #> magrittr 2.0.1 2020-11-17 [1] #> modelsummary * 0.8.1.9000 2021-08-09 [1] #> munsell 0.5.0 2018-06-12 [1] #> pillar 1.6.2 2021-07-29 [1] #> pkgconfig 2.0.3 2019-09-22 [1] #> purrr 0.3.4 2020-04-17 [1] #> R.cache 0.15.0 2021-04-30 [1] #> R.methodsS3 1.8.1 2020-08-26 [1] #> R.oo 1.24.0 2020-08-26 [1] #> R.utils 2.10.1 2020-08-26 [1] #> R6 2.5.0 2020-10-28 [1] #> reprex 2.0.0 2021-04-02 [1] #> rlang 0.4.11 2021-04-30 [1] #> rmarkdown 2.10 2021-08-06 [1] #> rstudioapi 0.13 2020-11-12 [1] #> rvest 1.0.1 2021-07-26 [1] #> scales 1.1.1 2020-05-11 [1] #> sessioninfo 1.1.1 2018-11-05 [1] #> stringi 1.7.3 2021-07-16 [1] #> stringr 1.4.0 2019-02-10 [1] #> styler 1.4.1 2021-03-30 [1] #> svglite 2.0.0 2021-02-20 [1] #> systemfonts 1.0.2 2021-05-11 [1] #> tables 0.9.6 2020-09-22 [1] #> tibble 3.1.3 2021-07-23 [1] #> tidyr 1.1.3 2021-03-03 [1] #> tidyselect 1.1.1 2021-04-30 [1] #> utf8 1.2.2 2021-07-24 [1] #> vctrs 0.3.8 2021-04-29 [1] #> viridisLite 0.4.0 2021-04-13 [1] #> webshot 0.5.2 2019-11-22 [1] #> withr 2.4.2 2021-04-18 [1] #> xfun 0.25 2021-08-06 [1] #> xml2 1.3.2 2020-04-23 [1] #> yaml 2.2.1 2020-02-01 [1] #> source #> standard (@0.2.1) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> standard (@1.1.1) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> Github (vincentarelbundock/modelsummary@4467537) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> standard (@2.0.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> CRAN (R 4.1.0) #> #> [1] /Library/Frameworks/R.framework/Versions/4.1/Resources/library ```

[ataquette]

OK so I am hoping this bit of code will be reproducible

## Create dataset ------------------------------------------------------------

set.seed(42)  ## for sake of reproducibility

multfin<-list() ## list of datasets
nbdatasets<-5  
n <- 240  ## size of each dataset

for (i in (1:nbdatasets)){ 

dat <- data.frame(
        age=sample(18:30, n, replace=TRUE),
        ethnic_group=rep(LETTERS[1:2], n/2),
        med_lastyr_binary_num=factor(rep(sample(c("Yes","No")), n/2)),
        vit_lastyr_binary_num=factor(rep(sample(c("Yes","No")), n/2)),
        maritalstatus=factor(rep(sample(c("Married","Other")), n/2)),
        gender=factor(rep(sample(c("Woman","Other"),prob=c(.95,.5), size=n, replace=TRUE))),
        bmi_group=factor(rep(sample(c("Normal Weigth","Overweight","Underweight")), n/3)),
        physact_before=factor(rep(sample(c("Yes","No")), n/2)),
        income_before=factor(rep(sample(c("Low","Middle","High")), n/3)),
        smoke_before=factor(rep(sample(c("Yes","No")), n/2)),
        contra_current_groupall=factor(rep(sample(c("Hormonal","Non-hormonal")), n/2)),
        cyclelength_before_group=factor(rep(sample(c("24-40","Other")), n/2)),
        period_lengh_before_group=factor(rep(sample(c("5/6","Other")), n/2)),
        heavyperiod_before=factor(rep(sample(c("Yes","No")), n/2)),
        cycle_irreg_before=factor(rep(sample(c("Yes","No")), n/2)),
        vaccine_type=factor(rep(sample(c("Pfizer","AstraZeneca"),prob=c(.55,.45), size=n, replace=TRUE))),
        vaccine_timing=factor(rep(sample(c("Before Jan 21","Feb 21","March 21", "April 21")), n/4)),
        covid_group=factor(rep(sample(c("No","Acute Covid","Long Covid")),prob=c(.75,.15,.10), size=n, replace=TRUE)),
        covid_tested=factor(rep(sample(c("No","Yes","False Negative"),prob=c(.75,.15,.10), size=n, replace=TRUE))),
        vaccinated=factor(rep(sample(c("One-shot","Two-shots")), n/2)),
        changesmenses_num=rnorm(n),
        lifesatchanges=rnorm(n),
        disease_before_endometriosis_num=factor(rep(sample(c("Yes","No")), n/2)),
        disease_before_pcos_num=factor(rep(sample(c("Yes","No")), n/2)),
        disease_before_underthyroid_num=factor(rep(sample(c("Yes","No")), n/2)),
        disease_before_overthyroid_num=factor(rep(sample(c("Yes","No")), n/2)),
        disease_before_uterpolyp_num=factor(rep(sample(c("Yes","No")), n/2)),
        disease_before_uterfibroid_num=factor(rep(sample(c("Yes","No")), n/2)),
        disease_before_intercystitis_num=factor(rep(sample(c("Yes","No")), n/2)),
        disease_before_eatingdisord_num=factor(rep(sample(c("Yes","No")), n/2)),
        id=1:n,
        vaccine_changescycle=factor(rep(sample(c("No","Yes, more disrupted","Yes, less disrupted","Other"),prob=c(.75,.15,.5,.5), size=n, replace=TRUE))))

# fill list
multfin[[i]]<-dat
}

## Univariate analysis ------------------------------------------------------------

# packages
library(Hmisc)
library(pbapply)
library(mice)
library(modelsummary)
library(nnet)
library(mlogit)

# set global R options
pboptions(type = "timer", char = "=") # initialize progress bar

# set list of exposures to test
exposures <- c(
    Cs(
        age,
        ethnic_group,
        med_lastyr_binary_num,
        vit_lastyr_binary_num,
        maritalstatus,
        gender,
        bmi_group,
        physact_before,
        income_before,
        smoke_before,
        contra_current_groupall,
        cyclelength_before_group,
        period_lengh_before_group,
        heavyperiod_before,
        cycle_irreg_before,
        vaccine_type,
        vaccine_timing,
        covid_group,
        covid_tested,
        vaccinated,
        changesmenses_num,
        lifesatchanges,
        disease_before_endometriosis_num,
        disease_before_pcos_num,
        disease_before_underthyroid_num,
        disease_before_overthyroid_num,
        disease_before_uterpolyp_num,
        disease_before_uterfibroid_num,
        disease_before_intercystitis_num,
        disease_before_eatingdisord_num

    )
)

## model function ------------------------------------------------------------
models <- function(x) {
    lapply(multfin, function(y)
        multinom(as.formula(
            paste0(
                "vaccine_changescycle ~ ",
                x
            )
        ), data = y, model = TRUE)
    )
}

## run  models
models_univariate <- as.list(seq(1,length(exposures))) # create list to store model results
models_univariate <- pblapply(exposures, models)

### Pool -------------------------------------------------------------------

pool_univariate <- as.list(seq(1,length(exposures))) # create list to store pool results

# run pool
for(j in seq_along(exposures)) {
    pool_univariate[[j]] <- pool(models_univariate[[j]])
}

However then when I try to extract a table for pool_univariate it does not work (but it used to before I set up the tidy.custom function)

library(modelsummary)
#> Warning: le package 'modelsummary' a été compilé avec la version R 4.0.5

tidy_custom.mipo <- function(x, ...) {
  out <- broom::tidy(x)
  out$bonferroni <- p.adjust(out$p.value, n = 30, method = "bonferroni")
  out$holm <- p.adjust(out$p.value, n = 30, method = "holm")
  return(out)
}

glance_custom.mipo<- function(x, ...) {
  out <- data.frame(
    "Num.Comparisons" = "30",
    "Model" = class(x)[1]
  )
  return(out)
}

## Model summary table output----------------------------------------------

msuni<-modelsummary(pool_univariate, 
                   group = model + term ~ y.level + y.level,
                    stars = c('*' = 0.05,'**' = 0.01,'***' = 0.001),
                    exponentiate = TRUE,estimate = "{estimate}{stars} [{conf.low}, {conf.high}]", 
                statistic = c(
                    "p = {p.value}",
                    "p (Holm) = {holm}"),
                    gof_omit = ".*",title = 'Risk factors for vaccine related changes in menstrual cycles',
                    coef_omit = "Intercept", 
                    output = "kableExtra")
#> Error in sanitize_models(models): objet 'pool_univariate' introuvable

[vincentarelbundock]

I see that you just copied my functions. But should we expect them to work as-is? Does the p.adjust function extract adjust p values from a mipo object? Probably not...

How would you normally adjust the p values from such objects?

[ataquette]

This is what I was using

## Get the q values 
library(data.table)
adjustment_pvalues<-function(model){ 

  smry <-lapply(pool_univariate,
       summary,
       conf.int = TRUE,
       conf.level = 0.95)
  all_models <- rbindlist(smry)
       # this is the FDR p values
  all_models$q.value <- p.adjust(all_models$p.value,method = "BH")
  return(all_models)
}  

This returns a data frame

[vincentarelbundock]

Right. So my guess is that code similar to this should go in your tidy_custom.mipo function.

[ataquette]

I did try that

tidy_custom.mipo = function(model) {
  out = data.frame(
    adjustment_pvalues(model),
   term = all_models$term,
   p.value = all_models$q.value)
return(out)
}

but I don't think that works....between the list of models and the dataframe with q. values....

[ataquette]

I'll try to get it to work - mustn't be far off

[vincentarelbundock]

Yeah, sorry. I don't know mipo enough to do this right now (would require me to dive in more than I can afford). The trick is that your tidy_custom.mipo function should take a SINGLE combined model as input, and return a SINGLE data frame with term and q.value column as output.

Since this is more about extract q values from mice objects, I'm going to close this issue, but feel free to keep the discussion going if you have specific questions or want to exchange ideas.