The modelsum function

Beth Atkinson, Ethan Heinzen, Pat Votruba, Jason Sinnwell, Shannon McDonnell and Greg Dougherty

Introduction

Very often we are asked to summarize model results from multiple fits into a nice table. The endpoint might be of different types (e.g., survival, case/control, continuous) and there may be several independent variables that we want to examine univariately or adjusted for certain variables such as age and sex. Locally at Mayo, the SAS macros %modelsum, %glmuniv, and %logisuni were written to create such summary tables. With the increasing interest in R, we have developed the function modelsum to create similar tables within the R environment.

In developing the modelsum function, the goal was to bring the best features of these macros into an R function. However, the task was not simply to duplicate all the functionality, but rather to make use of R’s strengths (modeling, method dispersion, flexibility in function definition and output format) and make a tool that fits the needs of R users. Additionally, the results needed to fit within the general reproducible research framework so the tables could be displayed within an R markdown report.

This report provides step-by-step directions for using the functions associated with modelsum. All functions presented here are available within the arsenal package. An assumption is made that users are somewhat familiar with R markdown documents. For those who are new to the topic, a good initial resource is available at rmarkdown.rstudio.com.

Simple Example

The first step when using the modelsum function is to load the arsenal package. All the examples in this report use a dataset called mockstudy made available by Paul Novotny which includes a variety of types of variables (character, numeric, factor, ordered factor, survival) to use as examples.

> require(arsenal)
> data(mockstudy) # load data
> dim(mockstudy)  # look at how many subjects and variables are in the dataset 
[1] 1499   14
> # help(mockstudy) # learn more about the dataset and variables
> str(mockstudy) # quick look at the data
'data.frame':   1499 obs. of  14 variables:
 $ case       : int  110754 99706 105271 105001 112263 86205 99508 90158 88989 90515 ...
 $ age        : int  67 74 50 71 69 56 50 57 51 63 ...
  ..- attr(*, "label")= chr "Age in Years"
 $ arm        : chr  "F: FOLFOX" "A: IFL" "A: IFL" "G: IROX" ...
  ..- attr(*, "label")= chr "Treatment Arm"
 $ sex        : Factor w/ 2 levels "Male","Female": 1 2 2 2 2 1 1 1 2 1 ...
 $ race       : chr  "Caucasian" "Caucasian" "Caucasian" "Caucasian" ...
  ..- attr(*, "label")= chr "Race"
 $ fu.time    : int  922 270 175 128 233 120 369 421 387 363 ...
 $ fu.stat    : int  2 2 2 2 2 2 2 2 2 2 ...
 $ ps         : int  0 1 1 1 0 0 0 0 1 1 ...
 $ hgb        : num  11.5 10.7 11.1 12.6 13 10.2 13.3 12.1 13.8 12.1 ...
 $ bmi        : num  25.1 19.5 NA 29.4 26.4 ...
  ..- attr(*, "label")= chr "Body Mass Index (kg/m^2)"
 $ alk.phos   : int  160 290 700 771 350 569 162 152 231 492 ...
 $ ast        : int  35 52 100 68 35 27 16 12 25 18 ...
 $ mdquality.s: int  NA 1 1 1 NA 1 1 1 1 1 ...
 $ age.ord    : Ord.factor w/ 8 levels "10-19"<"20-29"<..: 6 7 4 7 6 5 4 5 5 6 ...

To create a simple linear regression table (the default), use a formula statement to specify the variables that you want summarized. The example below predicts BMI with the variables sex and age.

> tab1 <- modelsum(bmi ~ sex + age, data=mockstudy)

If you want to take a quick look at the table, you can use summary on your modelsum object and the table will print out as text in your R console window. If you use summary without any options you will see a number of \(\&nbsp;\) statements which translates to “space” in HTML.

Pretty text version of table

If you want a nicer version in your console window then adding the text=TRUE option.

> summary(tab1, text=TRUE)


|             |estimate |std.error |p.value |adj.r.squared |Nmiss |
|:------------|:--------|:---------|:-------|:-------------|:-----|
|(Intercept)  |27.491   |0.181     |< 0.001 |0.004         |33    |
|sex Female   |-0.731   |0.290     |0.012   |              |      |
|(Intercept)  |26.424   |0.752     |< 0.001 |0.000         |33    |
|Age in Years |0.013    |0.012     |0.290   |              |      |

Pretty Rmarkdown version of table

In order for the report to look nice within an R markdown (knitr) report, you just need to specify results="asis" when creating the r chunk. This changes the layout slightly (compresses it) and bolds the variable names.

> summary(tab1)
estimate std.error p.value adj.r.squared Nmiss
(Intercept) 27.491 0.181 < 0.001 0.004 33
sex Female -0.731 0.290 0.012
(Intercept) 26.424 0.752 < 0.001 0.000 33
Age in Years 0.013 0.012 0.290

Data frame version of table

If you want a data.frame version, simply use as.data.frame.

> as.data.frame(tab1)
  y.term                  y.label strata.term adjustment model        term
1    bmi Body Mass Index (kg/m^2)             unadjusted     1 (Intercept)
2    bmi Body Mass Index (kg/m^2)             unadjusted     1   sexFemale
3    bmi Body Mass Index (kg/m^2)             unadjusted     2 (Intercept)
4    bmi Body Mass Index (kg/m^2)             unadjusted     2         age
         label term.type    estimate  std.error       p.value adj.r.squared
1  (Intercept) Intercept 27.49147713 0.18134740  0.000000e+00  3.632258e-03
2   sex Female      Term -0.73105055 0.29032223  1.190605e-02  3.632258e-03
3  (Intercept) Intercept 26.42372272 0.75211474 1.279109e-196  8.354809e-05
4 Age in Years      Term  0.01304859 0.01231653  2.895753e-01  8.354809e-05
  Nmiss
1    33
2    33
3    33
4    33

Add an adjustor to the model

The argument adjust allows the user to indicate that all the variables should be adjusted for these terms. To adjust each model for age and sex (for instance), we use adjust = ~ age + sex:

> tab2 <- modelsum(alk.phos ~ arm + ps + hgb, adjust= ~age + sex, data=mockstudy)
> summary(tab2)
estimate std.error p.value adj.r.squared Nmiss
(Intercept) 175.548 20.587 < 0.001 -0.001 266
Treatment Arm F: FOLFOX -13.701 8.730 0.117
Treatment Arm G: IROX -2.245 9.860 0.820
Age in Years -0.017 0.319 0.956
sex Female 3.016 7.521 0.688
(Intercept) 148.391 19.585 < 0.001 0.045 266
ps 46.721 5.987 < 0.001
Age in Years -0.084 0.311 0.787
sex Female 1.169 7.343 0.874
(Intercept) 336.554 32.239 < 0.001 0.031 266
hgb -13.845 2.137 < 0.001
Age in Years 0.095 0.314 0.763
sex Female -5.980 7.516 0.426

Models for each endpoint type

To make sure the correct model is run you need to specify “family”. The options available right now are : gaussian, binomial, survival, and poisson. If there is enough interest, additional models can be added.

Gaussian

Fit and summarize linear regression model

Look at whether there is any evidence that AlkPhos values vary by study arm after adjusting for sex and age (assuming a linear age relationship).

> fit <- lm(alk.phos ~ arm + age + sex, data=mockstudy)
> summary(fit)

Call:
lm(formula = alk.phos ~ arm + age + sex, data = mockstudy)

Residuals:
    Min      1Q  Median      3Q     Max 
-168.80  -81.45  -47.17   37.39  853.56 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  175.54808   20.58665   8.527   <2e-16 ***
armF: FOLFOX -13.70062    8.72963  -1.569    0.117    
armG: IROX    -2.24498    9.86004  -0.228    0.820    
age           -0.01741    0.31878  -0.055    0.956    
sexFemale      3.01598    7.52097   0.401    0.688    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 128.5 on 1228 degrees of freedom
  (266 observations deleted due to missingness)
Multiple R-squared:  0.002552,  Adjusted R-squared:  -0.0006969 
F-statistic: 0.7855 on 4 and 1228 DF,  p-value: 0.5346
> plot(fit)

The results suggest that the endpoint may need to be transformed. Calculating the Box-Cox transformation suggests a log transformation.

> require(MASS)
> boxcox(fit)

> fit2 <- lm(log(alk.phos) ~ arm + age + sex, data=mockstudy)
> summary(fit2)

Call:
lm(formula = log(alk.phos) ~ arm + age + sex, data = mockstudy)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.0098 -0.4470 -0.1065  0.4205  2.0620 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)   4.9692474  0.1025239  48.469   <2e-16 ***
armF: FOLFOX -0.0766798  0.0434746  -1.764    0.078 .  
armG: IROX   -0.0192828  0.0491041  -0.393    0.695    
age          -0.0004058  0.0015876  -0.256    0.798    
sexFemale     0.0179253  0.0374553   0.479    0.632    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6401 on 1228 degrees of freedom
  (266 observations deleted due to missingness)
Multiple R-squared:  0.003121,  Adjusted R-squared:  -0.0001258 
F-statistic: 0.9613 on 4 and 1228 DF,  p-value: 0.4278
> plot(fit2)

Finally, look to see whether there there is a non-linear relationship with age.

> require(splines)
Loading required package: splines
> fit3 <- lm(log(alk.phos) ~ arm + ns(age, df=2) + sex, data=mockstudy)
> 
> # test whether there is a difference between models 
> stats::anova(fit2,fit3)
Analysis of Variance Table

Model 1: log(alk.phos) ~ arm + age + sex
Model 2: log(alk.phos) ~ arm + ns(age, df = 2) + sex
  Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
1   1228 503.19                              
2   1227 502.07  1    1.1137 2.7218 0.09924 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> 
> # look at functional form of age
> termplot(fit3, term=2, se=T, rug=T)

In this instance it looks like there isn’t enough evidence to say that the relationship is non-linear.

Extract data using the broom package

The broom package makes it easy to extract information from the fit.

> tmp <- tidy(fit3) # coefficients, p-values
> class(tmp)
[1] "tbl_df"     "tbl"        "data.frame"
> tmp
# A tibble: 6 × 5
  term             estimate std.error statistic   p.value
  <chr>               <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)        4.76      0.141     33.8   1.93e-177
2 armF: FOLFOX      -0.0767    0.0434    -1.77  7.78e-  2
3 armG: IROX        -0.0195    0.0491    -0.396 6.92e-  1
4 ns(age, df = 2)1   0.330     0.260      1.27  2.04e-  1
5 ns(age, df = 2)2  -0.101     0.0935    -1.08  2.82e-  1
6 sexFemale          0.0183    0.0374     0.489 6.25e-  1
> 
> glance(fit3)
# A tibble: 1 × 12
  r.squared adj.r.squared sigma statistic p.value    df logLik   AIC   BIC
      <dbl>         <dbl> <dbl>     <dbl>   <dbl> <dbl>  <dbl> <dbl> <dbl>
1   0.00533       0.00127 0.640      1.31   0.255     5 -1196. 2405. 2441.
# ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

Create a summary table using modelsum

> ms.logy <- modelsum(log(alk.phos) ~ arm + ps + hgb, data=mockstudy, adjust= ~age + sex, 
+                     family=gaussian,  
+                     gaussian.stats=c("estimate","CI.lower.estimate","CI.upper.estimate","p.value"))
> summary(ms.logy)
estimate CI.lower.estimate CI.upper.estimate p.value
(Intercept) 4.969 4.768 5.170 < 0.001
Treatment Arm F: FOLFOX -0.077 -0.162 0.009 0.078
Treatment Arm G: IROX -0.019 -0.116 0.077 0.695
Age in Years -0.000 -0.004 0.003 0.798
sex Female 0.018 -0.056 0.091 0.632
(Intercept) 4.832 4.640 5.023 < 0.001
ps 0.226 0.167 0.284 < 0.001
Age in Years -0.001 -0.004 0.002 0.636
sex Female 0.009 -0.063 0.081 0.814
(Intercept) 5.765 5.450 6.080 < 0.001
hgb -0.069 -0.090 -0.048 < 0.001
Age in Years 0.000 -0.003 0.003 0.925
sex Female -0.027 -0.101 0.046 0.468

Binomial

Fit and summarize logistic regression model

> boxplot(age ~ mdquality.s, data=mockstudy, ylab=attr(mockstudy$age,'label'), xlab='mdquality.s')

> 
> fit <- glm(mdquality.s ~ age + sex, data=mockstudy, family=binomial)
> summary(fit)

Call:
glm(formula = mdquality.s ~ age + sex, family = binomial, data = mockstudy)

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  2.329442   0.514684   4.526 6.01e-06 ***
age         -0.002353   0.008256  -0.285    0.776    
sexFemale    0.039227   0.195330   0.201    0.841    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 807.68  on 1246  degrees of freedom
Residual deviance: 807.55  on 1244  degrees of freedom
  (252 observations deleted due to missingness)
AIC: 813.55

Number of Fisher Scoring iterations: 4
> 
> # create Odd's ratio w/ confidence intervals
> tmp <- data.frame(summary(fit)$coef)
> tmp
                Estimate  Std..Error    z.value     Pr...z..
(Intercept)  2.329441734 0.514683688  4.5259677 6.011977e-06
age         -0.002353404 0.008255814 -0.2850602 7.755980e-01
sexFemale    0.039227292 0.195330166  0.2008256 8.408350e-01
> 
> tmp$OR <- round(exp(tmp[,1]),2)
> tmp$lower.CI <- round(exp(tmp[,1] - 1.96* tmp[,2]),2)
> tmp$upper.CI <- round(exp(tmp[,1] + 1.96* tmp[,2]),2)
> names(tmp)[4] <- 'P-value'
> 
> kable(tmp[,c('OR','lower.CI','upper.CI','P-value')])
OR lower.CI upper.CI P-value
(Intercept) 10.27 3.75 28.17 0.000006
age 1.00 0.98 1.01 0.775598
sexFemale 1.04 0.71 1.53 0.840835 
> 
> # Assess the predictive ability of the model
> 
> # code using the pROC package
> require(pROC)
> pred <- predict(fit, type='response')
> tmp <- pROC::roc(mockstudy$mdquality.s[!is.na(mockstudy$mdquality.s)]~ pred, plot=TRUE, percent=TRUE)
Setting levels: control = 0, case = 1
Setting direction: controls < cases

> tmp$auc
Area under the curve: 50.69%

Extract data using broom package

The broom package makes it easy to extract information from the fit.

> tidy(fit, exp=T, conf.int=T) # coefficients, p-values, conf.intervals
# A tibble: 3 × 7
  term        estimate std.error statistic    p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>      <dbl>    <dbl>     <dbl>
1 (Intercept)   10.3     0.515       4.53  0.00000601    3.83      28.9 
2 age            0.998   0.00826    -0.285 0.776         0.981      1.01
3 sexFemale      1.04    0.195       0.201 0.841         0.712      1.53
> 
> glance(fit) # model summary statistics
# A tibble: 1 × 8
  null.deviance df.null logLik   AIC   BIC deviance df.residual  nobs
          <dbl>   <int>  <dbl> <dbl> <dbl>    <dbl>       <int> <int>
1          808.    1246  -404.  814.  829.     808.        1244  1247

Create a summary table using modelsum

> summary(modelsum(mdquality.s ~ age + bmi, data=mockstudy, adjust=~sex, family=binomial))
OR CI.lower.OR CI.upper.OR p.value concordance Nmiss
(Intercept) 10.272 3.831 28.876 < 0.001 0.507 252
Age in Years 0.998 0.981 1.014 0.776
sex Female 1.040 0.712 1.534 0.841
(Intercept) 4.814 1.709 13.221 0.003 0.550 273
Body Mass Index (kg/m^2) 1.023 0.987 1.063 0.220
sex Female 1.053 0.717 1.561 0.794 
> 
> fitall <- modelsum(mdquality.s ~ age, data=mockstudy, family=binomial,
+                    binomial.stats=c("Nmiss2","OR","p.value"))
> summary(fitall)
OR p.value Nmiss2
(Intercept) 10.493 < 0.001 252
Age in Years 0.998 0.766

Survival

Fit and summarize a Cox regression model

> require(survival)
Loading required package: survival
> 
> # multivariable model with all 3 terms
> fit  <- coxph(Surv(fu.time, fu.stat) ~ age + sex + arm, data=mockstudy)
> summary(fit)
Call:
coxph(formula = Surv(fu.time, fu.stat) ~ age + sex + arm, data = mockstudy)

  n= 1499, number of events= 1356 

                  coef exp(coef)  se(coef)      z Pr(>|z|)    
age           0.004600  1.004611  0.002501  1.839   0.0659 .  
sexFemale     0.039893  1.040699  0.056039  0.712   0.4765    
armF: FOLFOX -0.454650  0.634670  0.064878 -7.008 2.42e-12 ***
armG: IROX   -0.140785  0.868676  0.072760 -1.935   0.0530 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

             exp(coef) exp(-coef) lower .95 upper .95
age             1.0046     0.9954    0.9997    1.0095
sexFemale       1.0407     0.9609    0.9324    1.1615
armF: FOLFOX    0.6347     1.5756    0.5589    0.7207
armG: IROX      0.8687     1.1512    0.7532    1.0018

Concordance= 0.563  (se = 0.009 )
Likelihood ratio test= 56.21  on 4 df,   p=2e-11
Wald test            = 56.26  on 4 df,   p=2e-11
Score (logrank) test = 56.96  on 4 df,   p=1e-11
> 
> # check proportional hazards assumption
> fit.z <- cox.zph(fit)
> fit.z
       chisq df    p
age     1.41  1 0.24
sex     1.08  1 0.30
arm     1.80  2 0.41
GLOBAL  4.68  4 0.32
> plot(fit.z[1], resid=FALSE) # makes for a cleaner picture in this case
> abline(h=coef(fit)[1], col='red')

> 
> # check functional form for age using pspline (penalized spline)
> # results are returned for the linear and non-linear components
> fit2 <- coxph(Surv(fu.time, fu.stat) ~ pspline(age) + sex + arm, data=mockstudy)
> fit2
Call:
coxph(formula = Surv(fu.time, fu.stat) ~ pspline(age) + sex + 
    arm, data = mockstudy)

                         coef se(coef)      se2    Chisq   DF       p
pspline(age), linear  0.00443  0.00237  0.00237  3.48989 1.00  0.0617
pspline(age), nonlin                            13.11270 3.08  0.0047
sexFemale             0.03993  0.05610  0.05607  0.50663 1.00  0.4766
armF: FOLFOX         -0.46240  0.06494  0.06493 50.69608 1.00 1.1e-12
armG: IROX           -0.15243  0.07301  0.07299  4.35876 1.00  0.0368

Iterations: 6 outer, 16 Newton-Raphson
     Theta= 0.954 
Degrees of freedom for terms= 4.1 1.0 2.0 
Likelihood ratio test=70.1  on 7.08 df, p=2e-12
n= 1499, number of events= 1356 
> 
> # plot smoothed age to visualize why significant
> termplot(fit2, se=T, terms=1)
> abline(h=0)

> 
> # The c-statistic comes out in the summary of the fit
> summary(fit2)$concordance
          C       se(C) 
0.568432549 0.008487495 
> 
> # It can also be calculated using the survConcordance function
> survConcordance(Surv(fu.time, fu.stat) ~ predict(fit2), data=mockstudy)
Warning in survConcordance(Surv(fu.time, fu.stat) ~ predict(fit2), data = mockstudy): 'survConcordance' is deprecated.
Use 'concordance' instead.
See help("Deprecated")
Warning in survConcordance.fit(Y, x, strat, wt): 'survConcordance.fit' is deprecated.
Use 'concordancefit' instead.
See help("Deprecated")
$concordance
concordant 
 0.5684325 

$stats
concordant discordant  tied.risk  tied.time   std(c-d) 
 620221.00  470282.00    5021.00     766.00   19235.49 

$n
[1] 1499

$std.err
   std(c-d) 
0.008779125 

$call
survConcordance(formula = Surv(fu.time, fu.stat) ~ predict(fit2), 
    data = mockstudy)

attr(,"class")
[1] "survConcordance"

Extract data using broom package

The broom package makes it easy to extract information from the fit.

> tidy(fit) # coefficients, p-values
# A tibble: 4 × 5
  term         estimate std.error statistic  p.value
  <chr>           <dbl>     <dbl>     <dbl>    <dbl>
1 age           0.00460   0.00250     1.84  6.59e- 2
2 sexFemale     0.0399    0.0560      0.712 4.77e- 1
3 armF: FOLFOX -0.455     0.0649     -7.01  2.42e-12
4 armG: IROX   -0.141     0.0728     -1.93  5.30e- 2
> 
> glance(fit) # model summary statistics
# A tibble: 1 × 18
      n nevent statistic.log p.value.log statistic.sc p.value.sc statistic.wald
  <int>  <dbl>         <dbl>       <dbl>        <dbl>      <dbl>          <dbl>
1  1499   1356          56.2    1.81e-11         57.0   1.26e-11           56.3
# ℹ 11 more variables: p.value.wald <dbl>, statistic.robust <dbl>,
#   p.value.robust <dbl>, r.squared <dbl>, r.squared.max <dbl>,
#   concordance <dbl>, std.error.concordance <dbl>, logLik <dbl>, AIC <dbl>,
#   BIC <dbl>, nobs <dbl>

Create a summary table using modelsum

> ##Note: You must use quotes when specifying family="survival" 
> ##      family=survival will not work
> summary(modelsum(Surv(fu.time, fu.stat) ~ arm, 
+                  adjust=~age + sex, data=mockstudy, family="survival"))
HR CI.lower.HR CI.upper.HR p.value concordance
Treatment Arm F: FOLFOX 0.635 0.559 0.721 < 0.001 0.563
Treatment Arm G: IROX 0.869 0.753 1.002 0.053
Age in Years 1.005 1.000 1.010 0.066
sex Female 1.041 0.932 1.162 0.477 
> 
> ##Note: the pspline term is not working yet
> #summary(modelsum(Surv(fu.time, fu.stat) ~ arm, 
> #                adjust=~pspline(age) + sex, data=mockstudy, family='survival'))

Poisson

Poisson regression is useful when predicting an outcome variable representing counts. It can also be useful when looking at survival data. Cox models and Poisson models are very closely related and survival data can be summarized using Poisson regression. If you have overdispersion (see if the residual deviance is much larger than degrees of freedom), you may want to use quasipoisson() instead of poisson(). Some of these diagnostics need to be done outside of modelsum.

Example 1: fit and summarize a Poisson regression model

For the first example, use the solder dataset available in the rpart package. The endpoint skips has a definite Poisson look.

> require(rpart) ##just to get access to solder dataset
> data(solder)
> hist(solder$skips)

> 
> fit <- glm(skips ~ Opening + Solder + Mask , data=solder, family=poisson)
> stats::anova(fit, test='Chi')
Analysis of Deviance Table

Model: poisson, link: log

Response: skips

Terms added sequentially (first to last)

        Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                      899     8788.2              
Opening  2   2920.5       897     5867.7 < 2.2e-16 ***
Solder   1   1168.4       896     4699.3 < 2.2e-16 ***
Mask     4   2015.7       892     2683.7 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> summary(fit)

Call:
glm(formula = skips ~ Opening + Solder + Mask, family = poisson, 
    data = solder)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.12220    0.07742  -14.50  < 2e-16 ***
OpeningM     0.57161    0.05707   10.02  < 2e-16 ***
OpeningS     1.81475    0.05044   35.98  < 2e-16 ***
SolderThin   0.84682    0.03327   25.45  < 2e-16 ***
MaskA3       0.51315    0.07098    7.23 4.83e-13 ***
MaskA6       1.81103    0.06609   27.40  < 2e-16 ***
MaskB3       1.20225    0.06697   17.95  < 2e-16 ***
MaskB6       1.86648    0.06310   29.58  < 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: 8788.2  on 899  degrees of freedom
Residual deviance: 2683.7  on 892  degrees of freedom
AIC: 4802.2

Number of Fisher Scoring iterations: 5

Overdispersion is when the Residual deviance is larger than the degrees of freedom. This can be tested, approximately using the following code. The goal is to have a p-value that is \(>0.05\).

> 1-pchisq(fit$deviance, fit$df.residual)
[1] 0

One possible solution is to use the quasipoisson family instead of the poisson family. This adjusts for the overdispersion.

> fit2 <- glm(skips ~ Opening + Solder + Mask, data=solder, family=quasipoisson)
> summary(fit2)

Call:
glm(formula = skips ~ Opening + Solder + Mask, family = quasipoisson, 
    data = solder)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.12220    0.13483  -8.323 3.19e-16 ***
OpeningM     0.57161    0.09939   5.751 1.22e-08 ***
OpeningS     1.81475    0.08784  20.660  < 2e-16 ***
SolderThin   0.84682    0.05794  14.615  < 2e-16 ***
MaskA3       0.51315    0.12361   4.151 3.62e-05 ***
MaskA6       1.81103    0.11510  15.735  < 2e-16 ***
MaskB3       1.20225    0.11663  10.308  < 2e-16 ***
MaskB6       1.86648    0.10989  16.984  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for quasipoisson family taken to be 3.033198)

    Null deviance: 8788.2  on 899  degrees of freedom
Residual deviance: 2683.7  on 892  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 5

Extract data using broom package

The broom package makes it easy to extract information from the fit.

> tidy(fit) # coefficients, p-values
# A tibble: 8 × 5
  term        estimate std.error statistic   p.value
  <chr>          <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)   -1.12     0.0774    -14.5  1.29e- 47
2 OpeningM       0.572    0.0571     10.0  1.29e- 23
3 OpeningS       1.81     0.0504     36.0  1.66e-283
4 SolderThin     0.847    0.0333     25.5  6.47e-143
5 MaskA3         0.513    0.0710      7.23 4.83e- 13
6 MaskA6         1.81     0.0661     27.4  2.45e-165
7 MaskB3         1.20     0.0670     18.0  4.55e- 72
8 MaskB6         1.87     0.0631     29.6  2.71e-192
> 
> glance(fit) # model summary statistics
# A tibble: 1 × 8
  null.deviance df.null logLik   AIC   BIC deviance df.residual  nobs
          <dbl>   <int>  <dbl> <dbl> <dbl>    <dbl>       <int> <int>
1         8788.     899 -2393. 4802. 4841.    2684.         892   900

Create a summary table using modelsum

> summary(modelsum(skips~Opening + Solder + Mask, data=solder, family="quasipoisson"))
RR CI.lower.RR CI.upper.RR p.value
(Intercept) 1.533 1.179 1.952 < 0.001
Opening M 2.328 1.733 3.167 < 0.001
Opening S 7.491 5.780 9.888 < 0.001
(Intercept) 2.904 2.423 3.446 < 0.001
Solder Thin 2.808 2.295 3.458 < 0.001
(Intercept) 1.611 1.135 2.204 0.005
Mask A3 1.469 0.995 2.214 0.059
Mask A6 8.331 5.839 12.222 < 0.001
Mask B3 3.328 2.309 4.920 < 0.001
Mask B6 6.466 4.598 9.378 < 0.001 
> summary(modelsum(skips~Opening + Solder + Mask, data=solder, family="poisson"))
RR CI.lower.RR CI.upper.RR p.value
(Intercept) 1.533 1.397 1.678 < 0.001
Opening M 2.328 2.089 2.599 < 0.001
Opening S 7.491 6.805 8.267 < 0.001
(Intercept) 2.904 2.750 3.065 < 0.001
Solder Thin 2.808 2.637 2.992 < 0.001
(Intercept) 1.611 1.433 1.804 < 0.001
Mask A3 1.469 1.280 1.690 < 0.001
Mask A6 8.331 7.341 9.487 < 0.001
Mask B3 3.328 2.923 3.800 < 0.001
Mask B6 6.466 5.724 7.331 < 0.001

Example 2: fit and summarize a Poisson regression model

This second example uses the survival endpoint available in the mockstudy dataset. There is a close relationship between survival and Poisson models, and often it is easier to fit the model using Poisson regression, especially if you want to present absolute risk.

> # add .01 to the follow-up time (.01*1 day) in order to keep everyone in the analysis
> fit <- glm(fu.stat ~ offset(log(fu.time+.01)) + age + sex + arm, data=mockstudy, family=poisson)
> summary(fit)

Call:
glm(formula = fu.stat ~ offset(log(fu.time + 0.01)) + age + sex + 
    arm, family = poisson, data = mockstudy)

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -5.875627   0.108984 -53.913  < 2e-16 ***
age           0.003724   0.001705   2.184   0.0290 *  
sexFemale     0.027321   0.038575   0.708   0.4788    
armF: FOLFOX -0.335141   0.044600  -7.514 5.72e-14 ***
armG: IROX   -0.107776   0.050643  -2.128   0.0333 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 2113.5  on 1498  degrees of freedom
Residual deviance: 2048.0  on 1494  degrees of freedom
AIC: 5888.2

Number of Fisher Scoring iterations: 5
> 1-pchisq(fit$deviance, fit$df.residual)
[1] 0
> 
> coef(coxph(Surv(fu.time,fu.stat) ~ age + sex + arm, data=mockstudy))
         age    sexFemale armF: FOLFOX   armG: IROX 
 0.004600011  0.039892735 -0.454650445 -0.140784996 
> coef(fit)[-1]
         age    sexFemale armF: FOLFOX   armG: IROX 
 0.003723763  0.027320917 -0.335141090 -0.107775577 
> 
> # results from the Poisson model can then be described as risk ratios (similar to the hazard ratio)
> exp(coef(fit)[-1])
         age    sexFemale armF: FOLFOX   armG: IROX 
   1.0037307    1.0276976    0.7152372    0.8978291 
> 
> # As before, we can model the dispersion which alters the standard error
> fit2 <- glm(fu.stat ~ offset(log(fu.time+.01)) + age + sex + arm, 
+             data=mockstudy, family=quasipoisson)
> summary(fit2)

Call:
glm(formula = fu.stat ~ offset(log(fu.time + 0.01)) + age + sex + 
    arm, family = quasipoisson, data = mockstudy)

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -5.875627   0.566666 -10.369   <2e-16 ***
age           0.003724   0.008867   0.420    0.675    
sexFemale     0.027321   0.200572   0.136    0.892    
armF: FOLFOX -0.335141   0.231899  -1.445    0.149    
armG: IROX   -0.107776   0.263318  -0.409    0.682    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for quasipoisson family taken to be 27.03493)

    Null deviance: 2113.5  on 1498  degrees of freedom
Residual deviance: 2048.0  on 1494  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 5

Extract data using broom package

The broom package makes it easy to extract information from the fit.

> tidy(fit) ##coefficients, p-values
# A tibble: 5 × 5
  term         estimate std.error statistic  p.value
  <chr>           <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)  -5.88      0.109     -53.9   0       
2 age           0.00372   0.00171     2.18  2.90e- 2
3 sexFemale     0.0273    0.0386      0.708 4.79e- 1
4 armF: FOLFOX -0.335     0.0446     -7.51  5.72e-14
5 armG: IROX   -0.108     0.0506     -2.13  3.33e- 2
> 
> glance(fit) ##model summary statistics
# A tibble: 1 × 8
  null.deviance df.null logLik   AIC   BIC deviance df.residual  nobs
          <dbl>   <int>  <dbl> <dbl> <dbl>    <dbl>       <int> <int>
1         2114.    1498 -2939. 5888. 5915.    2048.        1494  1499

Create a summary table using modelsum

Remember that the result from modelsum is different from the fit above. The modelsum summary shows the results for age + offset(log(fu.time+.01)) then sex + offset(log(fu.time+.01)) instead of age + sex + arm + offset(log(fu.time+.01)).

> summary(modelsum(fu.stat ~ age, adjust=~offset(log(fu.time+.01))+ sex + arm, 
+                  data=mockstudy, family=poisson))
RR CI.lower.RR CI.upper.RR p.value
(Intercept) 0.003 0.002 0.003 < 0.001
Age in Years 1.004 1.000 1.007 0.029
sex Female 1.028 0.953 1.108 0.479
Treatment Arm F: FOLFOX 0.715 0.656 0.781 < 0.001
Treatment Arm G: IROX 0.898 0.813 0.991 0.033

Additional Examples

Here are multiple examples showing how to use some of the different options.

1. Change summary statistics globally

There are standard settings for each type of model regarding what information is summarized in the table. This behavior can be modified using the modelsum.control function. In fact, you can save your standard settings and use that for future tables.

> mycontrols  <- modelsum.control(gaussian.stats=c("estimate","std.error","adj.r.squared","Nmiss"),
+                                 show.adjust=FALSE, show.intercept=FALSE)                            
> tab2 <- modelsum(bmi ~ age, adjust=~sex, data=mockstudy, control=mycontrols)
> summary(tab2)
estimate std.error adj.r.squared Nmiss
Age in Years 0.012 0.012 0.004 33

You can also change these settings directly in the modelsum call.

> tab3 <- modelsum(bmi ~  age, adjust=~sex, data=mockstudy,
+                  gaussian.stats=c("estimate","std.error","adj.r.squared","Nmiss"), 
+                  show.intercept=FALSE, show.adjust=FALSE)
> summary(tab3)
estimate std.error adj.r.squared Nmiss
Age in Years 0.012 0.012 0.004 33

2. Add labels to independent variables

In the above example, age is shown with a label (Age in Years), but sex is listed “as is”. This is because the data was created in SAS and in the SAS dataset, age had a label but sex did not. The label is stored as an attribute within R.

> ## Look at one variable's label
> attr(mockstudy$age,'label')
[1] "Age in Years"
> 
> ## See all the variables with a label