Multiple Linear Regression Models with Bonferroni corrected p-values

Question

Please see my sample data below:


dad <- data.frame(type = c("new", "new", "old", "new", "old", "new", "old", "old", "new", "new"),
                   outcome = c(68,76,57,67,89,98,99,120,99,67),
                  loan1=c(98000 ,56668,67000,87999,87945, 89768, 65738, 87547, 78937, 20000),
                  loan2 =c(56768,45679, 23453, 87673, 56749, 45783, 67836, 54673, 45379, 78483),
                  house=c("semi", "detached", "town", "semi", "town", "town", "semi", "detached", "semi", "town"))

I am trying to run separate multivariate linear regression models for each house and would like to correct for multiple testing by adjusting the p-values with Bonferroni correction. In my actual dataset I am running 7 backward stepwise selection models (initially starting with the same independent variables) to come to the final multivariate regression models. Every model has the same outcome as well. I am wondering if there is a way to get the adjusted p-value with the below code at all? Any help would be appreciated, thank you!

#Example 
full <- lm(outcome ~loan1+loan2,data = subset(dad, house=="semi"))
summary(full)
tab_model(full)

Try `p.adjust` p.adjust(p, method = p.adjust.methods, n = length(p)) p.adjust.methods # c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY", # "fdr", "none") — GuedesBF, Mar 07 '23 at 23:06
Thanks for your response! I have tried this with little luck. Do you happen to know what the 'p' stands for? I assume n refers to the number of models I am running? — T K, Mar 07 '23 at 23:18
The 'p' stands for the 'p.value' element from the model. Se can get It with `model$p.value` or `model[["p.value]]` — GuedesBF, Mar 07 '23 at 23:42

score 1 · Answer 1 · answered Mar 08 '23 at 00:08

The tab_model() function has an argument for adjusted p-values.

library(sjPlot)
dad <- data.frame(type = c("new", "new", "old", "new", "old", "new", "old", "old", "new", "new"),
                  outcome = c(68,76,57,67,89,98,99,120,99,67),
                  loan1=c(98000 ,56668,67000,87999,87945, 89768, 65738, 87547, 78937, 20000),
                  loan2 =c(56768,45679, 23453, 87673, 56749, 45783, 67836, 54673, 45379, 78483),
                  house=c("semi", "detached", "town", "semi", "town", "town", "semi", "detached", "semi", "town"))

full <- lapply(unique(dad$house), \(h){
  args <- list(
    formula = outcome ~ loan1 + loan2, 
    data= subset(dad, house == h))
  do.call(lm, args)
})
tab_model(full[[1]], p.adjust="bonferroni")
tab_model(full[[2]], p.adjust="bonferroni")
tab_model(full[[3]], p.adjust="bonferroni")
})

GuedesBF · Answer 2 · 2023-03-08T01:20:24.503

0

Try p.adjust(). Call it with the argument method = 'bonferroni'

library(purrr)

Unnadjusted p.values:

full |>
    summary() |>
    pluck(coefficients) |>
    (\(x) x[-1, 4])()

    loan1     loan2 
0.0814972 0.1389506

adjusted p.values:

full |>
    summary() |>
    pluck(coefficients) |>
    (\(x) x[-1, 4])() |>
    p.adjust(method = 'bonferroni')

    loan1     loan2 
0.1629944 0.2779012

edited Mar 08 '23 at 01:20

answered Mar 07 '23 at 23:11

GuedesBF

8,409
5
19
37

Thanks for this, quick question. What does the following code do? "(\(x) x[-1, 4])() |>" – T K Mar 08 '23 at 16:02

Multiple Linear Regression Models with Bonferroni corrected p-values

2 Answers2