The modelsum function

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

06 September, 2019

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.

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.

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.

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:

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).

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

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

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

Create a summary table using modelsum

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

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

Create a summary table using modelsum

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
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

Create a summary table using modelsum

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

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.

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\).

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

Create a summary table using modelsum

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
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)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.1188  -0.4041   0.3242   0.9727   4.3588  

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)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.1188  -0.4041   0.3242   0.9727   4.3588  

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

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)).

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.

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.

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.

If you want to add labels to other variables, there are a couple of options. First, you could add labels to the variables in your dataset.

estimate std.error p.value adj.r.squared Nmiss
(Intercept) 26.793 0.766 < 0.001 0.004 33
Age, yrs 0.012 0.012 0.348
sex Female -0.718 0.291 0.014

You can also use the built-in data.frame method for labels<-:

estimate std.error p.value adj.r.squared Nmiss
(Intercept) 26.793 0.766 < 0.001 0.004 33
Age, yrs 0.012 0.012 0.348
sex Female -0.718 0.291 0.014

Another option is to add labels after you have created the table

estimate std.error p.value adj.r.squared Nmiss
(Intercept) 26.793 0.766 < 0.001 0.004 33
Age, yrs 0.012 0.012 0.348
Female -0.718 0.291 0.014

Alternatively, you can check the variable labels and manipulate them with a function called labels, which works on the modelsum object.

estimate std.error p.value adj.r.squared Nmiss
(Intercept) 26.793 0.766 < 0.001 0.004 33
Baseline Age (yrs) 0.012 0.012 0.348
Female -0.718 0.291 0.014

3. Don’t show intercept values

estimate std.error p.value adj.r.squared Nmiss
mdquality.s -0.326 1.093 0.766 -0.001 252
sex Female -1.208 0.610 0.048 0.002 0

4. Don’t show results for adjustment variables

OR CI.lower.OR CI.upper.OR p.value concordance Nmiss
(Intercept) 10.272 3.831 28.876 < 0.001 0.507 252
Age, yrs 0.998 0.981 1.014 0.776
(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

5. Summarize multiple variables without typing them out

Often one wants to summarize a number of variables. Instead of typing by hand each individual variable, an alternative approach is to create a formula using the paste command with the collapse="+" option.

[1] “ps+hgb+bmi+alk.phos+ast”

mdquality.s ~ ps + hgb + bmi + alk.phos + ast

OR CI.lower.OR CI.upper.OR p.value concordance Nmiss
(Intercept) 14.628 10.755 20.399 < 0.001 0.620 460
ps 0.461 0.332 0.639 < 0.001
(Intercept) 1.236 0.272 5.560 0.783 0.573 460
hgb 1.176 1.040 1.334 0.011
(Intercept) 4.963 1.818 13.292 0.002 0.549 273
Body Mass Index (kg/m^2) 1.023 0.987 1.062 0.225
(Intercept) 10.622 7.687 14.794 < 0.001 0.552 460
alk.phos 0.999 0.998 1.000 0.159
(Intercept) 10.936 7.912 15.232 < 0.001 0.545 460
ast 0.995 0.988 1.001 0.099

These steps can also be done using the formulize function.

mdquality.s ~ ps + hgb + bmi

OR CI.lower.OR CI.upper.OR p.value concordance Nmiss
(Intercept) 14.628 10.755 20.399 < 0.001 0.620 460
ps 0.461 0.332 0.639 < 0.001
(Intercept) 1.236 0.272 5.560 0.783 0.573 460
hgb 1.176 1.040 1.334 0.011
(Intercept) 4.963 1.818 13.292 0.002 0.549 273
Body Mass Index (kg/m^2) 1.023 0.987 1.062 0.225

6. Subset the dataset used in the analysis

Here are two ways to get the same result (limit the analysis to subjects age>50 and in the F: FOLFOX treatment group).

estimate std.error p.value adj.r.squared Nmiss
(Intercept) 122.577 46.924 0.009 -0.001 108
age 0.619 0.719 0.390
(Intercept) 164.814 7.673 < 0.001 -0.002 108
sex Female -5.497 12.118 0.650
(Intercept) 238.658 33.705 < 0.001 0.010 119
bmi -2.776 1.207 0.022
estimate std.error p.value adj.r.squared Nmiss
(Intercept) 4.872 0.039 < 0.001 -0.002 108
sex Female -0.005 0.062 0.931
(Intercept) 4.770 0.040 < 0.001 0.027 108
ps 0.183 0.050 < 0.001
(Intercept) 5.207 0.172 < 0.001 0.007 119
Body Mass Index (kg/m^2) -0.012 0.006 0.044
estimate std.error p.value adj.r.squared Nmiss
(Intercept) 178.812 14.550 < 0.001 0.007 77
ps 20.834 13.440 0.122
sex Female -17.542 16.656 0.293
(Intercept) 373.008 104.272 < 0.001 0.009 77
Body Mass Index (kg/m^2) -8.239 4.727 0.083
sex Female -24.058 16.855 0.155
estimate std.error p.value adj.r.squared Nmiss
(Intercept) 169.112 57.013 0.006 0.294 0
ps 254.901 68.100 < 0.001
sex Female 49.566 67.643 0.470
(Intercept) 453.070 200.651 0.033 -0.049 1
Body Mass Index (kg/m^2) -5.993 7.408 0.426
sex Female -22.308 79.776 0.782

7. Create combinations of variables on the fly

estimate std.error p.value adj.r.squared Nmiss
(Intercept) 59.714 1.314 < 0.001 0.003 252
interaction(mdquality.s, sex) 1.Male 0.730 1.385 0.598
interaction(mdquality.s, sex) 0.Female 0.988 2.134 0.643
interaction(mdquality.s, sex) 1.Female -1.021 1.425 0.474

8. Transform variables on the fly

Certain transformations need to be surrounded by I() so that R knows to treat it as a variable transformation and not some special model feature. If the transformation includes any of the symbols / - + ^ * then surround the new variable by I().

OR CI.lower.OR CI.upper.OR p.value concordance Nmiss
(Intercept) 0.656 0.382 1.124 0.126 0.514 0
Age, yrs 1.045 0.957 1.142 0.326
(Intercept) 0.633 0.108 3.698 0.611 0.508 33
Body Mass Index (kg/m^2) 1.092 0.638 1.867 0.748
(Intercept) 0.722 0.503 1.029 0.074 0.502 252
mdquality.s 1.045 0.719 1.527 0.819

9. Change the ordering of the variables or delete a variable

estimate std.error p.value adj.r.squared Nmiss
(Intercept) 26.424 0.752 < 0.001 0.000 33
Age, yrs 0.013 0.012 0.290
(Intercept) 27.491 0.181 < 0.001 0.004 33
sex Female -0.731 0.290 0.012
(Intercept) 27.944 0.253 < 0.001 0.011 294
alk.phos -0.005 0.001 < 0.001
estimate std.error p.value adj.r.squared Nmiss
(Intercept) 26.424 0.752 < 0.001 0.000 33
Age, yrs 0.013 0.012 0.290
(Intercept) 27.491 0.181 < 0.001 0.004 33
sex Female -0.731 0.290 0.012
estimate std.error p.value adj.r.squared Nmiss
(Intercept) 26.424 0.752 < 0.001 0.000 33
Age, yrs 0.013 0.012 0.290
(Intercept) 27.491 0.181 < 0.001 0.004 33
sex Female -0.731 0.290 0.012

10. Merge two modelsum objects together

It is possible to merge two modelsum objects so that they print out together, however you need to pay attention to the columns that are being displayed. It is sometimes easier to combine two models of the same family (such as two sets of linear models). Overlapping y-variables will have their x-variables concatenated, and (if all=TRUE) non-overlapping y-variables will have their tables printed separately.

[1] “modelsum” “arsenal_table”

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, yrs 0.013 0.012 0.290
OR CI.lower.OR CI.upper.OR p.value concordance Nmiss
(Intercept) 1.236 0.272 5.560 0.783 0.573 460
hgb 1.176 1.040 1.334 0.011
(Intercept) 10.622 7.687 14.794 < 0.001 0.552 460
alk.phos 0.999 0.998 1.000 0.159

11. Add a title to the table

When creating a pdf the tables are automatically numbered and the title appears below the table. In Word and HTML, the titles appear un-numbered and above the table.

Demographics
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, yrs 0.013 0.012 0.290

12. Modify how missing values are treated

Depending on the report you are writing you have the following options:

estimate N
(Intercept) 27.331 1205
ast -0.005
(Intercept) 26.424 1466
Age, yrs 0.013
estimate Nmiss2
(Intercept) 27.331 294
ast -0.005
(Intercept) 26.424 33
Age, yrs 0.013
estimate Nmiss
(Intercept) 27.331 294
ast -0.005
(Intercept) 26.424 33
Age, yrs 0.013
estimate
(Intercept) 27.331
ast -0.005
(Intercept) 26.424
Age, yrs 0.013

13. Modify the number of digits used

Within modelsum.control function there are 3 options for controlling the number of significant digits shown.

estimate std.error p.value adj.r.squared Nmiss
(Intercept) 27.4915 0.1813 < 0.001 0.0036 33
sex Female -0.7311 0.2903 0.012
(Intercept) 26.4237 0.7521 < 0.001 0.0001 33
Age, yrs 0.0130 0.0123 0.290
(Intercept) 26.4937 0.2447 < 0.001 0.0079 33
fu.time 0.0011 0.0003 < 0.001

14. Use case-weights in the models

Occasionally it is of interest to fit models using case weights. The modelsum function allows you to pass on the weights to the models and it will do the appropriate fit.

No Case Weights used
OR CI.lower.OR CI.upper.OR p.value concordance Nmiss
(Intercept) 0.590 0.473 0.735 < 0.001 0.550 210
ast 1.003 0.998 1.008 0.258
(Intercept) 0.578 0.306 1.093 0.091 0.500 29
Body Mass Index (kg/m^2) 1.003 0.980 1.026 0.808
(Intercept) 1.006 0.386 2.631 0.990 0.514 210
hgb 0.965 0.894 1.043 0.372
Case Weights used
OR CI.lower.OR CI.upper.OR p.value concordance Nmiss
(Intercept) 0.956 0.837 1.091 0.504 0.550 210
ast 1.003 1.000 1.006 0.068
(Intercept) 0.957 0.658 1.393 0.820 0.500 29
Body Mass Index (kg/m^2) 1.002 0.988 1.016 0.780
(Intercept) 1.829 1.031 3.248 0.039 0.514 210
hgb 0.956 0.913 1.001 0.058

15. Use modelsum within an Sweave document

For those users who wish to create tables within an Sweave document, the following code seems to work.

\documentclass{article}

\usepackage{longtable}
\usepackage{pdfpages}

\begin{document}

\section{Read in Data}
<<echo=TRUE>>=
require(arsenal)
require(knitr)
require(rmarkdown)
data(mockstudy)

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

\section{Convert Summary.modelsum to LaTeX}
<<echo=TRUE, results='hide', message=FALSE>>=
capture.output(summary(tab1), file="Test.md")

## Convert R Markdown Table to LaTeX
render("Test.md", pdf_document(keep_tex=TRUE))
@ 

\includepdf{Test.pdf}

\end{document}

16. Export modelsum results to a .CSV file

When looking at multiple variables it is sometimes useful to export the results to a csv file. The as.data.frame function creates a data frame object that can be exported or further manipulated within R.

17. Write modelsum object to a separate Word or HTML file

18. Use modelsum in R Shiny

The easiest way to output a modelsum() object in an R Shiny app is to use the tableOutput() UI in combination with the renderTable() server function and as.data.frame(summary(modelsum())):

This can be especially powerful if you feed the selections from a selectInput(multiple = TRUE) into formulize() to make the table dynamic!

23. Use modelsum in bookdown

Since the backbone of modelsum() is knitr::kable(), tables still render well in bookdown. However, print.summary.modelsum() doesn’t use the caption= argument of kable(), so some tables may not have a properly numbered caption. To fix this, use the method described on the bookdown site to give the table a tag/ID.

24. Model multiple endpoints

You can now use list() on the left-hand side of modelsum() to give multiple endpoints. Note that only one “family” can be specified this way (use merge() instead if you want multiple families).

estimate std.error p.value adj.r.squared Nmiss
(Intercept) 58.053 1.614 < 0.001 -0.001 33
Body Mass Index (kg/m^2) 0.059 0.055 0.289
Treatment Arm F: FOLFOX 0.593 0.718 0.408
Treatment Arm G: IROX 0.171 0.819 0.834
(Intercept) 60.108 0.597 < 0.001 0.001 0
sex Female -1.232 0.611 0.044
Treatment Arm F: FOLFOX 0.693 0.709 0.329
Treatment Arm G: IROX 0.148 0.812 0.855
estimate std.error p.value adj.r.squared Nmiss
(Intercept) 11.565 0.267 < 0.001 0.005 294
Body Mass Index (kg/m^2) 0.028 0.009 0.003
Treatment Arm F: FOLFOX 0.046 0.118 0.699
Treatment Arm G: IROX 0.065 0.133 0.624
(Intercept) 12.505 0.096 < 0.001 0.032 266
sex Female -0.642 0.099 < 0.001
Treatment Arm F: FOLFOX 0.131 0.115 0.256
Treatment Arm G: IROX 0.131 0.130 0.313

To avoid confusion about which table is which endpoint, you can set term.name=TRUE in summary(). This takes the labels for each endpoint and puts them in the top-left of the table.

Age, yrs estimate std.error p.value adj.r.squared Nmiss
(Intercept) 58.053 1.614 < 0.001 -0.001 33
Body Mass Index (kg/m^2) 0.059 0.055 0.289
Treatment Arm F: FOLFOX 0.593 0.718 0.408
Treatment Arm G: IROX 0.171 0.819 0.834
(Intercept) 60.108 0.597 < 0.001 0.001 0
sex Female -1.232 0.611 0.044
Treatment Arm F: FOLFOX 0.693 0.709 0.329
Treatment Arm G: IROX 0.148 0.812 0.855
hgb estimate std.error p.value adj.r.squared Nmiss
(Intercept) 11.565 0.267 < 0.001 0.005 294
Body Mass Index (kg/m^2) 0.028 0.009 0.003
Treatment Arm F: FOLFOX 0.046 0.118 0.699
Treatment Arm G: IROX 0.065 0.133 0.624
(Intercept) 12.505 0.096 < 0.001 0.032 266
sex Female -0.642 0.099 < 0.001
Treatment Arm F: FOLFOX 0.131 0.115 0.256
Treatment Arm G: IROX 0.131 0.130 0.313

25. Model data by a non-test group (strata)

You can also specify a grouping variable that doesn’t get tested (but instead separates results): a strata variable.

Treatment Arm estimate std.error p.value adj.r.squared Nmiss
A: IFL (Intercept) 59.147 2.783 < 0.001 -0.002 9
Body Mass Index (kg/m^2) 0.019 0.100 0.851
(Intercept) 59.903 0.683 < 0.001 -0.002 0
sex Female -0.651 1.151 0.572
F: FOLFOX (Intercept) 57.194 2.407 < 0.001 0.001 20
Body Mass Index (kg/m^2) 0.112 0.087 0.197
(Intercept) 60.691 0.574 < 0.001 0.000 0
sex Female -0.962 0.901 0.286
G: IROX (Intercept) 59.188 2.873 < 0.001 -0.003 4
Body Mass Index (kg/m^2) 0.023 0.104 0.822
(Intercept) 60.702 0.759 < 0.001 0.007 0
sex Female -2.346 1.200 0.051
Treatment Arm estimate std.error p.value adj.r.squared Nmiss
A: IFL (Intercept) 11.247 0.459 < 0.001 0.013 77
Body Mass Index (kg/m^2) 0.039 0.017 0.018
(Intercept) 12.527 0.109 < 0.001 0.037 69
sex Female -0.703 0.182 < 0.001
F: FOLFOX (Intercept) 11.661 0.414 < 0.001 0.004 157
Body Mass Index (kg/m^2) 0.026 0.015 0.085
(Intercept) 12.661 0.095 < 0.001 0.037 141
sex Female -0.707 0.151 < 0.001
G: IROX (Intercept) 11.874 0.457 < 0.001 0.001 60
Body Mass Index (kg/m^2) 0.019 0.017 0.264
(Intercept) 12.565 0.121 < 0.001 0.016 56
sex Female -0.470 0.188 0.013

26. Add multiple sets of adjustors to the model

By putting multiple formulas into a list, you can use multiple sets of adjustors. Use ~ 1 for an “unadjusted” model. By using the adjustment.names=TRUE argument and giving names to your adjustor sets in the list, you can name the various analyses.

> multi.adjust <- modelsum(list(age, bmi) ~ fu.time + ast, adjust = list(Unadjusted = ~ 1, "Adjusted for Arm" = ~ arm), data = mockstudy)
> summary(multi.adjust, adjustment.names = TRUE)


|adjustment       |                            |estimate |std.error |p.value |adj.r.squared |Nmiss |
|:----------------|:---------------------------|:--------|:---------|:-------|:-------------|:-----|
|Unadjusted       |(Intercept)                 |60.766   |0.512     |< 0.001 |0.002         |0     |
|                 |**fu.time**                 |-0.001   |0.001     |0.061   |              |      |
|Adjusted for Arm |(Intercept)                 |60.420   |0.663     |< 0.001 |0.002         |0     |
|                 |**fu.time**                 |-0.001   |0.001     |0.039   |              |      |
|                 |**Treatment Arm F: FOLFOX** |0.868    |0.717     |0.227   |              |      |
|                 |**Treatment Arm G: IROX**   |0.163    |0.812     |0.841   |              |      |
|Unadjusted       |(Intercept)                 |61.343   |0.547     |< 0.001 |0.004         |266   |
|                 |**ast**                     |-0.030   |0.012     |0.014   |              |      |
|Adjusted for Arm |(Intercept)                 |61.236   |0.757     |< 0.001 |0.005         |266   |
|                 |**ast**                     |-0.030   |0.012     |0.015   |              |      |
|                 |**Treatment Arm F: FOLFOX** |0.653    |0.779     |0.403   |              |      |
|                 |**Treatment Arm G: IROX**   |-0.728   |0.880     |0.408   |              |      |


|adjustment       |                            |estimate |std.error |p.value |adj.r.squared |Nmiss |
|:----------------|:---------------------------|:--------|:---------|:-------|:-------------|:-----|
|Unadjusted       |(Intercept)                 |26.494   |0.245     |< 0.001 |0.008         |33    |
|                 |**fu.time**                 |0.001    |0.000     |< 0.001 |              |      |
|Adjusted for Arm |(Intercept)                 |26.658   |0.317     |< 0.001 |0.007         |33    |
|                 |**fu.time**                 |0.001    |0.000     |< 0.001 |              |      |
|                 |**Treatment Arm F: FOLFOX** |-0.280   |0.341     |0.413   |              |      |
|                 |**Treatment Arm G: IROX**   |-0.237   |0.385     |0.538   |              |      |
|Unadjusted       |(Intercept)                 |27.331   |0.259     |< 0.001 |-0.000        |294   |
|                 |**ast**                     |-0.005   |0.006     |0.433   |              |      |
|Adjusted for Arm |(Intercept)                 |27.291   |0.356     |< 0.001 |-0.001        |294   |
|                 |**ast**                     |-0.004   |0.006     |0.440   |              |      |
|                 |**Treatment Arm F: FOLFOX** |0.181    |0.368     |0.623   |              |      |
|                 |**Treatment Arm G: IROX**   |-0.161   |0.414     |0.698   |              |      |
> summary(multi.adjust, adjustment.names = TRUE, show.intercept = FALSE, show.adjust = FALSE)


|adjustment       |            |estimate |std.error |p.value |adj.r.squared |Nmiss |
|:----------------|:-----------|:--------|:---------|:-------|:-------------|:-----|
|Unadjusted       |**fu.time** |-0.001   |0.001     |0.061   |0.002         |0     |
|Adjusted for Arm |**fu.time** |-0.001   |0.001     |0.039   |0.002         |0     |
|Unadjusted       |**ast**     |-0.030   |0.012     |0.014   |0.004         |266   |
|Adjusted for Arm |**ast**     |-0.030   |0.012     |0.015   |0.005         |266   |


|adjustment       |            |estimate |std.error |p.value |adj.r.squared |Nmiss |
|:----------------|:-----------|:--------|:---------|:-------|:-------------|:-----|
|Unadjusted       |**fu.time** |0.001    |0.000     |< 0.001 |0.008         |33    |
|Adjusted for Arm |**fu.time** |0.001    |0.000     |< 0.001 |0.007         |33    |
|Unadjusted       |**ast**     |-0.005   |0.006     |0.433   |-0.000        |294   |
|Adjusted for Arm |**ast**     |-0.004   |0.006     |0.440   |-0.001        |294   |

Available Function Options

Summary statistics

The available summary statistics, by varible type, are:

The full description of these parameters that can be shown for models include:

modelsum.control settings

A quick way to see what arguments are possible to utilize in a function is to use the args() command. Settings involving the number of digits can be set in modelsum.control or in summary.modelsum.

summary.modelsum settings

The summary.modelsum function has options that modify how the table appears (such as adding a title or modifying labels).