Multivariate Analysis

mv.model = glm(hospitalizations ~ sex+cesd+age, data = HELPrct, family = "poisson")
summary(mv.model)
## 
## Call:
## glm(formula = hospitalizations ~ sex + cesd + age, family = "poisson", 
##     data = HELPrct)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.7643  -1.7065  -0.8024   0.1558  20.0901  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.697066   0.157546  -4.425 9.67e-06 ***
## sexmale     -0.072835   0.061790  -1.179    0.238    
## cesd         0.013541   0.002204   6.145 8.01e-10 ***
## age          0.038229   0.003238  11.808  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2261.9  on 452  degrees of freedom
## Residual deviance: 2077.8  on 449  degrees of freedom
## AIC: 3100.2
## 
## Number of Fisher Scoring iterations: 6
mv.model %>%
  tbl_regression(exponentiate = TRUE) %>%
  as_gt() %>%
  gt::tab_source_note(gt::md("*This is mv model*"))
Characteristic IRR1 95% CI1 p-value
sex
female — —
male 0.93 0.82, 1.05 0.2
CESD at baseline 1.01 1.01, 1.02 <0.001
age (years) 1.04 1.03, 1.05 <0.001
This is mv model
1 IRR = Incidence Rate Ratio, CI = Confidence Interval

Univariate Model after producing univariate model of our own, and it ’s IRR matches exp(beta)

uv.model = glm(hospitalizations ~ cesd, data = HELPrct, family = "poisson")
summary(uv.model)
## 
## Call:
## glm(formula = hospitalizations ~ cesd, family = "poisson", data = HELPrct)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.9009  -1.5347  -0.8500   0.1062  21.5052  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 0.601920   0.082267   7.317 2.54e-13 ***
## cesd        0.015181   0.002212   6.864 6.71e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2261.9  on 452  degrees of freedom
## Residual deviance: 2213.8  on 451  degrees of freedom
## AIC: 3232.2
## 
## Number of Fisher Scoring iterations: 6
uv.model %>%
  tbl_regression(exponentiate = TRUE) %>%
  as_gt() %>%
  gt::tab_source_note(gt::md("*This is uv model*"))
Characteristic IRR1 95% CI1 p-value
CESD at baseline 1.02 1.01, 1.02 <0.001
This is uv model
1 IRR = Incidence Rate Ratio, CI = Confidence Interval

Univariate model table for all variables in data put toghter ( it matches with above IRR)

HELPrct %>%
 dplyr:: select(hospitalizations,sex,cesd,age) %>%
  tbl_uvregression(
    method = glm,
    y = hospitalizations,
    method.args = list(family = poisson(link="log")),
    exponentiate = TRUE,
    pvalue_fun = ~ style_pvalue(.x, digits = 2)
  ) %>%
  add_global_p() %>% # add global p-value
  add_nevent() %>% # add number of events of the outcome
  add_q() %>% # adjusts global p-values for multiple testing
  bold_p() %>% # bold p-values under a given threshold (default 0.05)
  bold_p(t = 0.10, q = TRUE) %>% # now bold q-values under the threshold of 0.10
  bold_labels()
## add_q: Adjusting p-values with
## `stats::p.adjust(x$table_body$p.value, method = "fdr")`
Characteristic N Event N IRR1 95% CI1 p-value q-value2
sex 453 1386 0.004 0.004
female — —
male 0.84 0.75, 0.95
CESD at baseline 453 1386 1.02 1.01, 1.02 <0.001 <0.001
age (years) 453 1386 1.04 1.03, 1.05 <0.001 <0.001
1 IRR = Incidence Rate Ratio, CI = Confidence Interval
2 False discovery rate correction for multiple testing