Bug Report: Missing Confidence Intervals and P-Values in tbl_regression for Survey Design with Many Independent Variables

[zhaohongxin0]

Hello,

I’d like to report a bug in the tbl_regression function when creating regression tables for survey designs with many independent variables. The issue is that the confidence intervals and p-values are missing from the table.

my R data attached:
data_frame.zip

Below is a simplified version of my analysis code:

library(gtsummary)
library(survey)

data_frame <- data_to_save

# Constructing a complex survey design object
svy_design <- survey::svydesign(strata = ~SDMVSTRA, id = ~SDMVPSU, weights = ~WTINT2YR, nest = TRUE, data = data_frame)

# Fitting a survey-weighted logistic regression model
result <- survey::svyglm(RIAGENDR ~ DMDHRGND + DMDCITZN + SIALANG + DMDHRAGZ + RIDRETH3 + DMDHRMAZ + INDFMPIR, family = binomial(), design = svy_design)

tbl_regression(result, exponentiate = TRUE)

Here’s an example of the output, where the confidence intervals and p-values are missing:

Interestingly, when we examine the model summary, the standard errors are present, but the p-values are missing:

# Print model summary with statistical measures
summary(result)

According to the svyglm function documentation, specifying the degrees of freedom brings back the p-values:

# Print model summary with degrees of freedom adjustment (df = degf)
summary(result, df.resid = degf(result$survey.design))

# Print model summary with degrees of freedom adjustment (df = Inf)
summary(result, df.resid = Inf)

In summary, when the number of independent variables is large, the tbl_regression output lacks confidence intervals and p-values. I suspect that the missing confidence intervals are due to tbl_regression extracting profile likelihood CI, and perhaps switching to Wald CI might resolve this. The missing p-values seem related to the degrees of freedom setting, which works fine when using df.resid = Inf or df.resid = degf(result$survey.design) in the summary() function. The question is, how can this setting be passed to tbl_regression?

Could you please look into this issue? It's becoming a significant limitation as most of my analyses involve many independent variables, making the tbl_regression results unusable.

I've attached my R dataset for your reference.

Thank you!

[larmarange]

Hi. I will look at it on Tuesday.

In the meantime, could you explore broom::tidy() results? Can you pass df.resid?

Alternatively, could you check broom.helpers::tidy_parameters?

Gtsummary::tbl_regression() uses broom.helpers::tidy_plus_plus() who uses broom::tidy()

Any additional parameters passed to tbl_regression() will be passed through to broom::tidy().

So the question lies with broom::tidy(). Can we pass the number of degrees of freedom?

[larmarange]

After exploration, it appears that the number of degrees of freedom of your model residuals is 0, resulting in a problem in computing p-values and confidence intervals by default. I guess that this is due to an insufficient number of observations in regard to your survey design and the number of coefficients to estimate in your model. However, my stats skills are insufficient to answer that specific point and to know if such model is still valid.

Now, as seen below, df.resid cannot be passed to the tidy() method. Eventually you can pass ddf to fix confidence intervals, but with no effect on p-values. This issue has been mentioned here tidymodels/broom#1096 although not fixed.

A quick fix, manually change mod$df.residual.

library(gtsummary)
library(survey)
#> Le chargement a nécessité le package : grid
#> Le chargement a nécessité le package : Matrix
#> Le chargement a nécessité le package : survival
#> 
#> Attachement du package : 'survey'
#> L'objet suivant est masqué depuis 'package:graphics':
#> 
#>     dotchart

load("./data_frame.RData")

data_frame <- data_to_save

# Constructing a complex survey design object
svy_design <- survey::svydesign(
  strata = ~SDMVSTRA, 
  id = ~SDMVPSU, 
  weights = ~WTINT2YR, 
  nest = TRUE, 
  data = data_frame
)

# Fitting a survey-weighted logistic regression model
mod <- survey::svyglm(
  RIAGENDR ~ DMDHRGND + DMDCITZN + SIALANG + 
    DMDHRAGZ + RIDRETH3 + DMDHRMAZ + INDFMPIR, 
  family = binomial(), 
  design = svy_design
)
#> Warning in eval(family$initialize): nombre de succès non entier dans un glm
#> binomial !

broom::tidy(mod, conf.int = TRUE)
#> # A tibble: 16 × 7
#>    term                  estimate std.error statistic p.value conf.low conf.high
#>    <chr>                    <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
#>  1 (Intercept)             0.462     0.0980     4.71      NaN      NaN       NaN
#>  2 DMDHRGNDMale           -0.979     0.0643   -15.2       NaN      NaN       NaN
#>  3 DMDCITZNDon't Know      0.127     0.720      0.177     NaN      NaN       NaN
#>  4 DMDCITZNNot a citize…   0.0695    0.150      0.464     NaN      NaN       NaN
#>  5 DMDCITZNRefused         1.15      0.785      1.46      NaN      NaN       NaN
#>  6 SIALANGSpanish         -0.0331    0.107     -0.308     NaN      NaN       NaN
#>  7 DMDHRAGZ40-59 years     0.0523    0.0806     0.649     NaN      NaN       NaN
#>  8 DMDHRAGZ60+ years       0.271     0.0763     3.55      NaN      NaN       NaN
#>  9 RIDRETH3Non-Hispanic…   0.159     0.0717     2.22      NaN      NaN       NaN
#> 10 RIDRETH3Non-Hispanic…   0.0324    0.0724     0.447     NaN      NaN       NaN
#> 11 RIDRETH3Non-Hispanic…  -0.0242    0.0676    -0.359     NaN      NaN       NaN
#> 12 RIDRETH3Other Hispan…  -0.116     0.0909    -1.28      NaN      NaN       NaN
#> 13 RIDRETH3Other Race -…  -0.130     0.147     -0.883     NaN      NaN       NaN
#> 14 DMDHRMAZNever Married  -0.0941    0.0761    -1.24      NaN      NaN       NaN
#> 15 DMDHRMAZWidowed/Divo…   0.0337    0.0775     0.436     NaN      NaN       NaN
#> 16 INDFMPIR                0.0124    0.0220     0.563     NaN      NaN       NaN

broom::tidy(mod, conf.int = TRUE, df.resid = Inf)
#> Error in `broom_confint_terms()`:
#> ! Arguments in `...` must be used.
#> ✖ Problematic argument:
#> • df.resid = Inf
#> ℹ Did you misspell an argument name?

broom.helpers::tidy_parameters(mod, df.resid = Inf)
#> Warning: Model has zero degrees of freedom!
#> Warning: Model has zero degrees of freedom!
#> Warning in stats::pt(-abs(statistic$Statistic), df = dof): Production de NaN
#> Warning: Model has zero degrees of freedom!
#>                                           term    estimate  std.error
#> 1                                  (Intercept)  0.46217286 0.09804752
#> 2                                 DMDHRGNDMale -0.97886910 0.06432897
#> 3                           DMDCITZNDon't Know  0.12734062 0.72046327
#> 4              DMDCITZNNot a citizen of the US  0.06951323 0.14984462
#> 5                              DMDCITZNRefused  1.14990120 0.78521016
#> 6                               SIALANGSpanish -0.03307634 0.10734405
#> 7                          DMDHRAGZ40-59 years  0.05234068 0.08064448
#> 8                            DMDHRAGZ60+ years  0.27056091 0.07629910
#> 9                   RIDRETH3Non-Hispanic Asian  0.15921068 0.07168616
#> 10                  RIDRETH3Non-Hispanic Black  0.03236700 0.07238431
#> 11                  RIDRETH3Non-Hispanic White -0.02424546 0.06761533
#> 12                      RIDRETH3Other Hispanic -0.11625534 0.09091682
#> 13 RIDRETH3Other Race - Including Multi-Racial -0.12976618 0.14700891
#> 14                       DMDHRMAZNever Married -0.09408349 0.07612769
#> 15          DMDHRMAZWidowed/Divorced/Separated  0.03373998 0.07746102
#> 16                                    INDFMPIR  0.01236432 0.02195053
#>    conf.level conf.low conf.high   statistic df.error p.value
#> 1        0.95      NaN       NaN   4.7137640        0     NaN
#> 2        0.95      NaN       NaN -15.2166131        0     NaN
#> 3        0.95      NaN       NaN   0.1767482        0     NaN
#> 4        0.95      NaN       NaN   0.4639021        0     NaN
#> 5        0.95      NaN       NaN   1.4644502        0     NaN
#> 6        0.95      NaN       NaN  -0.3081339        0     NaN
#> 7        0.95      NaN       NaN   0.6490299        0     NaN
#> 8        0.95      NaN       NaN   3.5460563        0     NaN
#> 9        0.95      NaN       NaN   2.2209402        0     NaN
#> 10       0.95      NaN       NaN   0.4471550        0     NaN
#> 11       0.95      NaN       NaN  -0.3585794        0     NaN
#> 12       0.95      NaN       NaN  -1.2787000        0     NaN
#> 13       0.95      NaN       NaN  -0.8827096        0     NaN
#> 14       0.95      NaN       NaN  -1.2358642        0     NaN
#> 15       0.95      NaN       NaN   0.4355737        0     NaN
#> 16       0.95      NaN       NaN   0.5632809        0     NaN

broom.helpers::tidy_parameters(mod, ddf = Inf)
#> Warning: Model has zero degrees of freedom!
#> Warning: Model has zero degrees of freedom!
#> Warning in stats::pt(-abs(statistic$Statistic), df = dof): Production de NaN
#> Warning: Model has zero degrees of freedom!
#>                                           term    estimate  std.error
#> 1                                  (Intercept)  0.46217286 0.09804752
#> 2                                 DMDHRGNDMale -0.97886910 0.06432897
#> 3                           DMDCITZNDon't Know  0.12734062 0.72046327
#> 4              DMDCITZNNot a citizen of the US  0.06951323 0.14984462
#> 5                              DMDCITZNRefused  1.14990120 0.78521016
#> 6                               SIALANGSpanish -0.03307634 0.10734405
#> 7                          DMDHRAGZ40-59 years  0.05234068 0.08064448
#> 8                            DMDHRAGZ60+ years  0.27056091 0.07629910
#> 9                   RIDRETH3Non-Hispanic Asian  0.15921068 0.07168616
#> 10                  RIDRETH3Non-Hispanic Black  0.03236700 0.07238431
#> 11                  RIDRETH3Non-Hispanic White -0.02424546 0.06761533
#> 12                      RIDRETH3Other Hispanic -0.11625534 0.09091682
#> 13 RIDRETH3Other Race - Including Multi-Racial -0.12976618 0.14700891
#> 14                       DMDHRMAZNever Married -0.09408349 0.07612769
#> 15          DMDHRMAZWidowed/Divorced/Separated  0.03373998 0.07746102
#> 16                                    INDFMPIR  0.01236432 0.02195053
#>    conf.level conf.low conf.high   statistic df.error p.value
#> 1        0.95      NaN       NaN   4.7137640        0     NaN
#> 2        0.95      NaN       NaN -15.2166131        0     NaN
#> 3        0.95      NaN       NaN   0.1767482        0     NaN
#> 4        0.95      NaN       NaN   0.4639021        0     NaN
#> 5        0.95      NaN       NaN   1.4644502        0     NaN
#> 6        0.95      NaN       NaN  -0.3081339        0     NaN
#> 7        0.95      NaN       NaN   0.6490299        0     NaN
#> 8        0.95      NaN       NaN   3.5460563        0     NaN
#> 9        0.95      NaN       NaN   2.2209402        0     NaN
#> 10       0.95      NaN       NaN   0.4471550        0     NaN
#> 11       0.95      NaN       NaN  -0.3585794        0     NaN
#> 12       0.95      NaN       NaN  -1.2787000        0     NaN
#> 13       0.95      NaN       NaN  -0.8827096        0     NaN
#> 14       0.95      NaN       NaN  -1.2358642        0     NaN
#> 15       0.95      NaN       NaN   0.4355737        0     NaN
#> 16       0.95      NaN       NaN   0.5632809        0     NaN

# Trick, manually changing the df number for residuals

mod$df.residual <- Inf

broom::tidy(mod, conf.int = TRUE)
#> # A tibble: 16 × 7
#>    term                 estimate std.error statistic  p.value conf.low conf.high
#>    <chr>                   <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
#>  1 (Intercept)            0.462     0.0980     4.71  2.43e- 6   0.270     0.654 
#>  2 DMDHRGNDMale          -0.979     0.0643   -15.2   2.74e-52  -1.10     -0.853 
#>  3 DMDCITZNDon't Know     0.127     0.720      0.177 8.60e- 1  -1.28      1.54  
#>  4 DMDCITZNNot a citiz…   0.0695    0.150      0.464 6.43e- 1  -0.224     0.363 
#>  5 DMDCITZNRefused        1.15      0.785      1.46  1.43e- 1  -0.389     2.69  
#>  6 SIALANGSpanish        -0.0331    0.107     -0.308 7.58e- 1  -0.243     0.177 
#>  7 DMDHRAGZ40-59 years    0.0523    0.0806     0.649 5.16e- 1  -0.106     0.210 
#>  8 DMDHRAGZ60+ years      0.271     0.0763     3.55  3.91e- 4   0.121     0.420 
#>  9 RIDRETH3Non-Hispani…   0.159     0.0717     2.22  2.64e- 2   0.0187    0.300 
#> 10 RIDRETH3Non-Hispani…   0.0324    0.0724     0.447 6.55e- 1  -0.110     0.174 
#> 11 RIDRETH3Non-Hispani…  -0.0242    0.0676    -0.359 7.20e- 1  -0.157     0.108 
#> 12 RIDRETH3Other Hispa…  -0.116     0.0909    -1.28  2.01e- 1  -0.294     0.0619
#> 13 RIDRETH3Other Race …  -0.130     0.147     -0.883 3.77e- 1  -0.418     0.158 
#> 14 DMDHRMAZNever Marri…  -0.0941    0.0761    -1.24  2.17e- 1  -0.243     0.0551
#> 15 DMDHRMAZWidowed/Div…   0.0337    0.0775     0.436 6.63e- 1  -0.118     0.186 
#> 16 INDFMPIR               0.0124    0.0220     0.563 5.73e- 1  -0.0307    0.0554

mod |> 
  tbl_regression(
    exponentiate = TRUE
  ) |> 
  as_kable()

Characteristic	OR	95% CI	p-value
DMDHRGND
Female	—	—
Male	0.38	0.33, 0.43	<0.001
DMDCITZN
Citizen by birth or naturalization	—	—
Don’t Know	1.14	0.28, 4.66	0.9
Not a citizen of the US	1.07	0.80, 1.44	0.6
Refused	3.16	0.68, 14.7	0.14
SIALANG
English	—	—
Spanish	0.97	0.78, 1.19	0.8
DMDHRAGZ
20-39 years	—	—
40-59 years	1.05	0.90, 1.23	0.5
60+ years	1.31	1.13, 1.52	<0.001
RIDRETH3
Mexican American	—	—
Non-Hispanic Asian	1.17	1.02, 1.35	0.026
Non-Hispanic Black	1.03	0.90, 1.19	0.7
Non-Hispanic White	0.98	0.85, 1.11	0.7
Other Hispanic	0.89	0.74, 1.06	0.2
Other Race - Including Multi-Racial	0.88	0.66, 1.17	0.4
DMDHRMAZ
Married/Living with partner	—	—
Never Married	0.91	0.78, 1.06	0.2
Widowed/Divorced/Separated	1.03	0.89, 1.20	0.7
INDFMPIR	1.01	0.97, 1.06	0.6

^{Created on 2024-08-20 with reprex v2.1.1}

[ddsjoberg]

Thanks both! Hopefully broom will allow the update!