Regression Tables from ‘GLM’, ‘GEE’, ‘GLMM’, ‘Cox’ and ‘survey’ Results for Publication.
## 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"geeglm object from
geepack packagelibrary(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 "1 (0.96,1.03)" "0.896" "1 (0.96,1.03)" "0.896"
#> Cu: ref.=Cu000 NA NA NA NA
#> 035 "1.03 (0.77,1.38)" "0.832" "1.03 (0.77,1.38)" "0.832"
#> 175 "1.16 (0.86,1.56)" "0.327" "1.16 (0.86,1.56)" "0.327"
#>
#> $metric
#> crude OR(95%CI) crude P value adj. OR(95%CI)
#> NA NA NA
#> Estimated correlation parameters "-0.016" 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 NAlmerMod or glmerMod
object from lme4 packagelibrary(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 1 (0.96,1.04) 0.9016730 1 (0.96,1.04)
#> Cu: ref.=Cu000 <NA> NA <NA>
#> 035 1.03 (0.74,1.43) 0.8494035 1.03 (0.74,1.43)
#> 175 1.16 (0.83,1.62) 0.3778824 1.16 (0.83,1.62)
#> Random effects <NA> NA <NA>
#> Pig 0 NA <NA>
#> Metrics <NA> NA <NA>
#> No. of groups (Pig) 72 NA <NA>
#> No. of observations 861 NA <NA>
#> Log-likelihood -596.35 NA <NA>
#> AIC value 1202.71 NA <NA>
#> adj. P value
#> Time 0.9009070
#> Cu: ref.=Cu000 NA
#> 035 0.8489406
#> 175 0.3777180
#> 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)"frailty or cluster
optionslibrary(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 "1.61 (1.25,2.08)" "< 0.001" "1.56 (1.22,2)" "< 0.001"
#> age "1.02 (1.01,1.03)" "0.007" "1.01 (1,1.02)" "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.64 (1.31,2.05)" "< 0.001" "1.58 (1.26,1.99)" "< 0.001"
#> age "1.02 (1,1.04)" "0.041" "1.01 (0.99,1.03)" "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"coxme object from
coxme packagelibrary(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"svyglm object from
survey packagelibrary(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"svycoxph object from
survey packagedata(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')"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)