Introduce jstable

Jinseob Kim

2024-10-20

Regression Tables from ‘GLM’, ‘GEE’, ‘GLMM’, ‘Cox’ and ‘survey’ Results for Publication.

Install

remotes::install_github("jinseob2kim/jstable")
library(jstable)

GLM Table

## Gaussian
glm_gaussian <- glm(mpg ~ cyl + disp, data = mtcars)
glmshow.display(glm_gaussian, decimal = 2)
#> $first.line
#> [1] "Linear regression predicting mpg\n"
#> 
#> $table
#>      crude coeff.(95%CI)   crude P value adj. coeff.(95%CI)    adj. P value
#> cyl  "-2.88 (-3.51,-2.24)" "< 0.001"     "-1.59 (-2.98,-0.19)" "0.034"     
#> disp "-0.04 (-0.05,-0.03)" "< 0.001"     "-0.02 (-0.04,0)"     "0.054"     
#> 
#> $last.lines
#> [1] "No. of observations = 32\nR-squared = 0.7596\nAIC value = 167.1456\n\n"
#> 
#> attr(,"class")
#> [1] "display" "list"

## Binomial
glm_binomial <- glm(vs ~ cyl + disp, data = mtcars, family = binomial)
glmshow.display(glm_binomial, decimal = 2)
#> $first.line
#> [1] "Logistic regression predicting vs\n"
#> 
#> $table
#>      crude OR.(95%CI)   crude P value adj. OR.(95%CI)    adj. P value
#> cyl  "0.2 (0.08,0.56)"  "0.002"       "0.15 (0.02,1.02)" "0.053"     
#> disp "0.98 (0.97,0.99)" "0.002"       "1 (0.98,1.03)"    "0.715"     
#> 
#> $last.lines
#> [1] "No. of observations = 32\nAIC value = 23.8304\n\n"
#> 
#> attr(,"class")
#> [1] "display" "list"

GEE Table: from geeglm object from geepack package

library(geepack) ## for dietox data
data(dietox)
dietox$Cu <- as.factor(dietox$Cu)
dietox$ddn <- as.numeric(rnorm(nrow(dietox)) > 0)
gee01 <- geeglm(Weight ~ Time + Cu, id = Pig, data = dietox, family = gaussian, corstr = "ex")
geeglm.display(gee01)
#> $caption
#> [1] "GEE(gaussian) predicting Weight by Time, Cu - Group Pig"
#> 
#> $table
#>                crude coeff(95%CI)   crude P value adj. coeff(95%CI)  
#> Time           "6.94 (6.79,7.1)"    "< 0.001"     "6.94 (6.79,7.1)"  
#> Cu: ref.=Cu000 NA                   NA            NA                 
#>       035      "-0.59 (-3.73,2.54)" "0.711"       "-0.84 (-3.9,2.23)"
#>       175      "1.9 (-1.87,5.66)"   "0.324"       "1.77 (-1.9,5.45)" 
#>                adj. P value
#> Time           "< 0.001"   
#> Cu: ref.=Cu000 NA          
#>       035      "0.593"     
#>       175      "0.345"     
#> 
#> $metric
#>                                  crude coeff(95%CI) crude P value
#>                                  NA                 NA           
#> Estimated correlation parameters "0.775"            NA           
#> No. of clusters                  "72"               NA           
#> No. of observations              "861"              NA           
#>                                  adj. coeff(95%CI) adj. P value
#>                                  NA                NA          
#> Estimated correlation parameters NA                NA          
#> No. of clusters                  NA                NA          
#> No. of observations              NA                NA

gee02 <- geeglm(ddn ~ Time + Cu, id = Pig, data = dietox, family = binomial, corstr = "ex")
geeglm.display(gee02)
#> $caption
#> [1] "GEE(binomial) predicting ddn by Time, Cu - Group Pig"
#> 
#> $table
#>                crude OR(95%CI)    crude P value adj. OR(95%CI)     adj. P value
#> Time           "0.94 (0.91,0.98)" "0.005"       "0.94 (0.91,0.98)" "0.005"     
#> Cu: ref.=Cu000 NA                 NA            NA                 NA          
#>       035      "1.04 (0.71,1.53)" "0.842"       "1.04 (0.71,1.54)" "0.833"     
#>       175      "1.03 (0.71,1.5)"  "0.887"       "1.03 (0.7,1.5)"   "0.882"     
#> 
#> $metric
#>                                  crude OR(95%CI) crude P value adj. OR(95%CI)
#>                                  NA              NA            NA            
#> Estimated correlation parameters "0.018"         NA            NA            
#> No. of clusters                  "72"            NA            NA            
#> No. of observations              "861"           NA            NA            
#>                                  adj. P value
#>                                  NA          
#> Estimated correlation parameters NA          
#> No. of clusters                  NA          
#> No. of observations              NA

Mixed model Table: lmerMod or glmerMod object from lme4 package

library(lme4)
l1 <- lmer(Weight ~ Time + Cu + (1 | Pig), data = dietox)
lmer.display(l1, ci.ranef = T)
#> $table
#>                      crude coeff(95%CI) crude P value adj. coeff(95%CI)
#> Time                   6.94 (6.88,7.01)     0.0000000  6.94 (6.88,7.01)
#> Cu: ref.=Cu000                     <NA>            NA              <NA>
#>       035            -0.58 (-4.67,3.51)     0.7811327 -0.84 (-4.47,2.8)
#>       175              1.9 (-2.23,6.04)     0.3670740  1.77 (-1.9,5.45)
#> Random effects                     <NA>            NA              <NA>
#> Pig                 40.34 (28.08,54.95)            NA              <NA>
#> Residual             11.37 (10.3,12.55)            NA              <NA>
#> Metrics                            <NA>            NA              <NA>
#> No. of groups (Pig)                  72            NA              <NA>
#> No. of observations                 861            NA              <NA>
#> Log-likelihood                  -2400.8            NA              <NA>
#> AIC value                        4801.6            NA              <NA>
#>                     adj. P value
#> Time                   0.0000000
#> Cu: ref.=Cu000                NA
#>       035              0.6527264
#>       175              0.3442309
#> Random effects                NA
#> Pig                           NA
#> Residual                      NA
#> Metrics                       NA
#> No. of groups (Pig)           NA
#> No. of observations           NA
#> Log-likelihood                NA
#> AIC value                     NA
#> 
#> $caption
#> [1] "Linear mixed model fit by REML : Weight ~ Time + Cu + (1 | Pig)"

l2 <- glmer(ddn ~ Time + Cu + (1 | Pig), data = dietox, family = "binomial")
lmer.display(l2)
#> $table
#>                      crude OR(95%CI) crude P value   adj. OR(95%CI)
#> Time                0.94 (0.91,0.98)   0.004070315 0.94 (0.91,0.98)
#> Cu: ref.=Cu000                  <NA>            NA             <NA>
#>       035           1.04 (0.73,1.49)   0.828431852  1.04 (0.72,1.5)
#>       175           1.03 (0.71,1.48)   0.880735113 1.03 (0.71,1.49)
#> Random effects                  <NA>            NA             <NA>
#> Pig                             0.07            NA             <NA>
#> Metrics                         <NA>            NA             <NA>
#> No. of groups (Pig)               72            NA             <NA>
#> No. of observations              861            NA             <NA>
#> Log-likelihood               -590.41            NA             <NA>
#> AIC value                    1190.83            NA             <NA>
#>                     adj. P value
#> Time                 0.004059061
#> Cu: ref.=Cu000                NA
#>       035            0.818591808
#>       175            0.875496753
#> Random effects                NA
#> Pig                           NA
#> Metrics                       NA
#> No. of groups (Pig)           NA
#> No. of observations           NA
#> Log-likelihood                NA
#> AIC value                     NA
#> 
#> $caption
#> [1] "Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) : ddn ~ Time + Cu + (1 | Pig)"

Cox model with frailty or cluster options

library(survival)
fit1 <- coxph(Surv(time, status) ~ ph.ecog + age, cluster = inst, lung, model = T) ## model = T: to extract original data
fit2 <- coxph(Surv(time, status) ~ ph.ecog + age + frailty(inst), lung, model = T)
cox2.display(fit1)
#> $table
#>         crude HR(95%CI)       crude P value adj. HR(95%CI)     adj. P value
#> ph.ecog "38.34 (38.32,38.36)" "< 0.001"     "33.03 (33,33.06)" "< 0.001"   
#> age     "14.67 (14.46,14.88)" "0.007"       "5.59 (4.73,6.61)" "0.085"     
#> 
#> $ranef
#>         [,1] [,2] [,3] [,4]
#> cluster   NA   NA   NA   NA
#> inst      NA   NA   NA   NA
#> 
#> $metric
#>                         [,1] [,2] [,3] [,4]
#> <NA>                      NA   NA   NA   NA
#> No. of observations  226.000   NA   NA   NA
#> No. of events        163.000   NA   NA   NA
#> AIC                 1463.797   NA   NA   NA
#> 
#> $caption
#> [1] "Marginal Cox model on time ('time') to event ('status') - Group inst"
cox2.display(fit2)
#> $table
#>         crude HR(95%CI) crude P value adj. HR(95%CI)     adj. P value
#> ph.ecog "1 (0.38,2.63)" "< 0.001"     "1 (0.41,2.47)"    "< 0.001"   
#> age     "1.04 (1,1.08)" "0.041"       "1.25 (1.22,1.28)" "0.225"     
#> 
#> $ranef
#>         [,1] [,2] [,3] [,4]
#> frailty   NA   NA   NA   NA
#> inst      NA   NA   NA   NA
#> 
#> $metric
#>                         [,1] [,2] [,3] [,4]
#> <NA>                      NA   NA   NA   NA
#> No. of observations  226.000   NA   NA   NA
#> No. of events        163.000   NA   NA   NA
#> AIC                 1463.223   NA   NA   NA
#> 
#> $caption
#> [1] "Frailty Cox model on time ('time') to event ('status') - Group inst"

Cox mixed effect model Table: coxme object from coxme package

library(coxme)
fit <- coxme(Surv(time, status) ~ ph.ecog + age + (1 | inst), lung)
coxme.display(fit)
#> $table
#>         crude HR(95%CI)    crude P value adj. HR(95%CI)     adj. P value
#> ph.ecog "1.66 (1.32,2.09)" "< 0.001"     "1.61 (1.27,2.03)" "< 0.001"   
#> age     "1.02 (1,1.04)"    "0.043"       "1.01 (0.99,1.03)" "0.227"     
#> 
#> $ranef
#>                 [,1] [,2] [,3] [,4]
#> Random effect     NA   NA   NA   NA
#> inst(Intercept) 0.02   NA   NA   NA
#> 
#> $metric
#>                     [,1] [,2] [,3] [,4]
#> <NA>                  NA   NA   NA   NA
#> No. of groups(inst)   18   NA   NA   NA
#> No. of observations  226   NA   NA   NA
#> No. of events        163   NA   NA   NA
#> 
#> $caption
#> [1] "Mixed effects Cox model on time ('time') to event ('status') - Group inst"

GLM for survey data : svyglm object from survey package

library(survey)
data(api)
apistrat$tt <- c(rep(1, 20), rep(0, nrow(apistrat) - 20))
apistrat$tt2 <- factor(c(rep(0, 40), rep(1, nrow(apistrat) - 40)))

dstrat <- svydesign(id = ~1, strata = ~stype, weights = ~pw, data = apistrat, fpc = ~fpc)
ds <- svyglm(api00 ~ ell + meals + tt2, design = dstrat)
ds2 <- svyglm(tt ~ ell + meals + tt2, design = dstrat, family = quasibinomial())
svyregress.display(ds)
#> $first.line
#> [1] "Linear regression predicting api00- weighted data\n"
#> 
#> $table
#>             crude coeff.(95%CI)    crude P value adj. coeff.(95%CI)   
#> ell         "-3.73 (-4.36,-3.1)"   "< 0.001"     "-0.51 (-1.27,0.26)" 
#> meals       "-3.38 (-3.71,-3.05)"  "< 0.001"     "-3.11 (-3.65,-2.57)"
#> tt2: 1 vs 0 "10.98 (-34.44,56.39)" "0.634"       "6.24 (-17.83,30.32)"
#>             adj. P value
#> ell         "0.195"     
#> meals       "< 0.001"   
#> tt2: 1 vs 0 "0.61"      
#> 
#> $last.lines
#> [1] "No. of observations = 200\nAIC value = 2308.0628\n\n"
#> 
#> attr(,"class")
#> [1] "display" "list"
svyregress.display(ds2)
#> $first.line
#> [1] "Logistic regression predicting tt- weighted data\n"
#> 
#> $table
#>             crude OR.(95%CI)   crude P value adj. OR.(95%CI)    adj. P value
#> ell         "1.02 (1,1.05)"    "0.047"       "1.11 (1.02,1.21)" "0.02"      
#> meals       "1.01 (0.99,1.03)" "0.255"       "0.97 (0.93,1.01)" "0.151"     
#> tt2: 1 vs 0 "0 (0,0)"          "< 0.001"     "0 (0,0)"          "< 0.001"   
#> 
#> $last.lines
#> [1] "No. of observations = 200\n\n"
#> 
#> attr(,"class")
#> [1] "display" "list"

Cox model for survey data :svycoxph object from survey package

data(pbc, package = "survival")
pbc$sex <- factor(pbc$sex)
pbc$stage <- factor(pbc$stage)
pbc$randomized <- with(pbc, !is.na(trt) & trt > 0)
biasmodel <- glm(randomized ~ age * edema, data = pbc, family = binomial)
pbc$randprob <- fitted(biasmodel)

if (is.null(pbc$albumin)) pbc$albumin <- pbc$alb ## pre2.9.0

dpbc <- svydesign(id = ~1, prob = ~randprob, strata = ~edema, data = subset(pbc, randomized))

model <- svycoxph(Surv(time, status > 0) ~ sex + protime + albumin + stage, design = dpbc)
svycox.display(model)
#> Stratified Independent Sampling design (with replacement)
#> svydesign(id = ~1, prob = ~randprob, strata = ~edema, data = subset(pbc, 
#>     randomized))
#> Stratified Independent Sampling design (with replacement)
#> svydesign(id = ~1, prob = ~randprob, strata = ~edema, data = subset(pbc, 
#>     randomized))
#> Stratified Independent Sampling design (with replacement)
#> svydesign(id = ~1, prob = ~randprob, strata = ~edema, data = subset(pbc, 
#>     randomized))
#> Stratified Independent Sampling design (with replacement)
#> svydesign(id = ~1, prob = ~randprob, strata = ~edema, data = subset(pbc, 
#>     randomized))
#> Stratified Independent Sampling design (with replacement)
#> svydesign(id = ~1, prob = ~randprob, strata = ~edema, data = subset(pbc, 
#>     randomized))
#> $table
#>               crude HR(95%CI)      crude P value adj. HR(95%CI)       
#> sex: f vs m   "0.62 (0.4,0.97)"    "0.038"       "0.55 (0.33,0.9)"    
#> protime       "1.37 (1.09,1.72)"   "0.006"       "1.52 (1.2,1.91)"    
#> albumin       "0.2 (0.14,0.29)"    "< 0.001"     "0.31 (0.2,0.47)"    
#> stage: ref.=1 NA                   NA            NA                   
#>    2          "5.67 (0.77,41.78)"  "0.089"       "10.94 (1.01,118.54)"
#>    3          "9.78 (1.37,69.94)"  "0.023"       "17.03 (1.69,171.6)" 
#>    4          "22.89 (3.2,163.47)" "0.002"       "22.56 (2.25,226.41)"
#>               adj. P value
#> sex: f vs m   "0.017"     
#> protime       "< 0.001"   
#> albumin       "< 0.001"   
#> stage: ref.=1 NA          
#>    2          "0.049"     
#>    3          "0.016"     
#>    4          "0.008"     
#> 
#> $metric
#>                        [,1] [,2] [,3] [,4]
#> <NA>                     NA   NA   NA   NA
#> No. of observations  312.00   NA   NA   NA
#> No. of events        144.00   NA   NA   NA
#> AIC                 1483.12   NA   NA   NA
#> 
#> $caption
#> [1] "Survey cox model on time ('time') to event ('status > 0')"

Sub-group analysis for Cox/svycox model

library(dplyr)
lung %>%
  mutate(
    status = as.integer(status == 1),
    sex = factor(sex),
    kk = factor(as.integer(pat.karno >= 70)),
    kk1 = factor(as.integer(pat.karno >= 60))
  ) -> lung

# TableSubgroupMultiCox(Surv(time, status) ~ sex, var_subgroups = c("kk", "kk1"), data=lung, line = T)

## Survey data
library(survey)
data.design <- svydesign(id = ~1, data = lung, weights = ~1)
# TableSubgroupMultiCox(Surv(time, status) ~ sex, var_subgroups = c("kk", "kk1"), data = data.design, line = F)