New scripts for revised regression with gtsummary. Nice, but slow!
This commit is contained in:
AG Damsbo
parent
2d7f0fc186
commit
152668f778
6
.gitignore
vendored
6
.gitignore
vendored
@ -45,3 +45,9 @@ docs/
|
||||
# translation temp files
|
||||
po/*~
|
||||
|
||||
# Customs
|
||||
*.RTF
|
||||
*.png
|
||||
*.html
|
||||
*.zip
|
||||
*.pdf
|
||||
|
||||
58
biv_mul.R
Normal file
58
biv_mul.R
Normal file
@ -0,0 +1,58 @@
|
||||
## =============================================================================
|
||||
## Requirements
|
||||
## =============================================================================
|
||||
|
||||
library(gtsummary)
|
||||
|
||||
## =============================================================================
|
||||
## Loop
|
||||
## =============================================================================
|
||||
|
||||
bm_list<-list()
|
||||
|
||||
for (i in 1:length(outs)){
|
||||
## Bivariate
|
||||
biv<-dta|>
|
||||
select(all_of(c("active_treat",vars,outs[i])))|>
|
||||
tbl_uvregression(data=_,
|
||||
y=outs[i],
|
||||
method=lm,
|
||||
label = lab_sel(labels_all,vars)
|
||||
) |>
|
||||
add_global_p()|>
|
||||
bold_p() |>
|
||||
bold_labels() |>
|
||||
italicize_levels()
|
||||
|
||||
## Multivariate
|
||||
mul<-dta |>
|
||||
select(all_of(c("active_treat",vars,outs[i])))|>
|
||||
lm(formula(paste(c(outs[i],"."),collapse="~")),
|
||||
data = _) |>
|
||||
tbl_regression(label = lab_sel(labels_all,vars)
|
||||
)|>
|
||||
add_n() |>
|
||||
add_global_p() |>
|
||||
bold_p() |>
|
||||
bold_labels() |>
|
||||
italicize_levels()
|
||||
|
||||
## Merge
|
||||
biv_mul<-tbl_merge(
|
||||
tbls = list(biv, mul),
|
||||
tab_spanner = c("**Bivariate linear regression**",
|
||||
"**Multivariate linear regression**")
|
||||
)
|
||||
|
||||
bm_list[[i]]<-biv_mul
|
||||
}
|
||||
|
||||
## =============================================================================
|
||||
## Big merge
|
||||
## =============================================================================
|
||||
|
||||
bm_16_tbl<-tbl_merge(
|
||||
tbls = bm_list,
|
||||
tab_spanner = c("**One month follow up**",
|
||||
"**Six months follow up**")
|
||||
)
|
||||
70
biv_mul_strat.R
Normal file
70
biv_mul_strat.R
Normal file
@ -0,0 +1,70 @@
|
||||
## =============================================================================
|
||||
## Requirements
|
||||
## =============================================================================
|
||||
|
||||
library(gtsummary)
|
||||
|
||||
## =============================================================================
|
||||
## Loop
|
||||
## =============================================================================
|
||||
|
||||
bm_list_strat<-list()
|
||||
|
||||
for (i in 1:length(outs)){
|
||||
## Bivariate
|
||||
biv<-dta |>
|
||||
select(all_of(c("active_treat",vars,outs[i]))) |>
|
||||
tbl_strata(
|
||||
strata = active_treat,
|
||||
.tbl_fun =
|
||||
~ .x %>%
|
||||
tbl_uvregression(data=.,
|
||||
y=outs[i],
|
||||
method=lm,
|
||||
label = lab_sel(labels_all,vars)
|
||||
)|>
|
||||
add_global_p()|>
|
||||
bold_p() |>
|
||||
bold_labels() |>
|
||||
italicize_levels(),
|
||||
.header = "**{strata}**, N = {n}"
|
||||
)
|
||||
|
||||
## Multivariate
|
||||
mul<-dta |>
|
||||
select(all_of(c("active_treat",vars,outs[i]))) |>
|
||||
tbl_strata(
|
||||
strata = active_treat,
|
||||
.tbl_fun =
|
||||
~ .x %>%
|
||||
lm(formula(paste(c(outs[i],"."),collapse="~")),
|
||||
data = .) |>
|
||||
tbl_regression(label = lab_sel(labels_all,vars)
|
||||
)|>
|
||||
add_n() |>
|
||||
add_global_p() |>
|
||||
bold_p() |>
|
||||
bold_labels() |>
|
||||
italicize_levels(),
|
||||
.header = "**{strata}**, N = {n}"
|
||||
)
|
||||
|
||||
## Merge
|
||||
biv_mul_strat<-tbl_merge(
|
||||
tbls = list(biv, mul),
|
||||
tab_spanner = c("**Bivariate linear regression**",
|
||||
"**Multivariate linear regression**")
|
||||
)
|
||||
|
||||
bm_list_strat[[i]]<-biv_mul_strat
|
||||
}
|
||||
|
||||
## =============================================================================
|
||||
## Big merge
|
||||
## =============================================================================
|
||||
|
||||
bm_16_tbl_strat<-tbl_merge(
|
||||
tbls = bm_list_strat,
|
||||
tab_spanner = c("**One month follow up**",
|
||||
"**Six months follow up**")
|
||||
)
|
||||
248
regression_rev.Rmd
Normal file
248
regression_rev.Rmd
Normal file
@ -0,0 +1,248 @@
|
||||
---
|
||||
title: "revised statistics"
|
||||
author: "AGDamsbo"
|
||||
date: "8/1/2022"
|
||||
output: html_document
|
||||
toc: true
|
||||
---
|
||||
|
||||
```{r setup, include=FALSE}
|
||||
knitr::opts_chunk$set(echo = TRUE)
|
||||
```
|
||||
|
||||
```{r}
|
||||
library(REDCapR)
|
||||
library(gtsummary)
|
||||
theme_gtsummary_compact()
|
||||
library(REDCapR)
|
||||
library(gt)
|
||||
|
||||
library(lubridate)
|
||||
library(dplyr)
|
||||
library(tidyr)
|
||||
|
||||
```
|
||||
|
||||
|
||||
|
||||
```{r}
|
||||
# token_talos<-read.csv("/Users/au301842/talos_redcap_token.csv",colClasses = "character")|>
|
||||
# names()|>
|
||||
# (\(x){ ## Shorthand for "anonymous lambda function"
|
||||
# substr(x,2,33)})()|>
|
||||
# suppressWarnings()
|
||||
#
|
||||
# library(REDCapR)
|
||||
# dta <- redcap_read_oneshot(
|
||||
# redcap_uri = "https://redcap.au.dk/api/",
|
||||
# token = token_talos
|
||||
# )$data|>
|
||||
# select(-c("cpr"))
|
||||
```
|
||||
|
||||
## 03aug22: Currently awaiting resolution to error message: "ERROR: You do not have permissions to use the API"
|
||||
## 03aug22: now working. No idea how...
|
||||
## Data stored in redcap is raw and in no way enriched. This process has to be put in a script.
|
||||
## Using locally stored data for now.
|
||||
|
||||
# Import
|
||||
|
||||
|
||||
```{r}
|
||||
dta_all<-read.csv("/Volumes/Data/depression/dep_dataset.csv")
|
||||
```
|
||||
|
||||
# Defining patients to include for analysis
|
||||
Only including cases with complete pase_0 and MDI at 1 & 6 months
|
||||
```{r}
|
||||
dta<-dta_all[!is.na(dta_all$pase_0),]
|
||||
# &!is.na(dta$mdi_1)&!is.na(dta$mdi_6)
|
||||
```
|
||||
|
||||
## Formatting
|
||||
```{r echo=FALSE}
|
||||
dta$diabetes<-factor(dta$diabetes)
|
||||
dta$pad<-factor(dta$pad)
|
||||
|
||||
dta$cohab<-ifelse(dta$civil=="partner","yes","no")|>
|
||||
factor()
|
||||
|
||||
dta$hypertension<-factor(dta$hypertension)
|
||||
|
||||
dta$afli[dta$afli=="unknown"]<-NA
|
||||
dta$afli<-factor(dta$afli)
|
||||
|
||||
dta$ever_smoker<-ifelse(dta$smoke_ever=="ever","yes","no")|>
|
||||
factor()
|
||||
|
||||
dta$ami<-factor(dta$ami)
|
||||
dta$tci<-factor(dta$tci)
|
||||
dta$thrombolysis<-factor(dta$thrombolysis)
|
||||
dta$thrombechtomy<-factor(dta$thrombechtomy)
|
||||
|
||||
dta$any_reperf<-ifelse(dta$rep_any=="rep","yes","no")|>
|
||||
factor()
|
||||
|
||||
dta$pad<-factor(dta$pad)
|
||||
dta$nihss_0<-as.numeric(dta$nihss_0)
|
||||
dta$age<-as.numeric(dta$age)
|
||||
|
||||
dta$active_treat<-ifelse(dta$rtreat=="Active","yes","no")|>
|
||||
factor()
|
||||
# dta$rtreat<-factor(dta$rtreat)
|
||||
|
||||
dta$female<-ifelse(dta$sex=="female","yes","no")|>
|
||||
factor()
|
||||
|
||||
dta$pase_0<-as.numeric(dta$pase_0)
|
||||
dta$pase_6<-as.numeric(dta$pase_6)
|
||||
dta$bmi<-as.numeric(dta$bmi)
|
||||
dta$mdi_6<-as.numeric(dta$mdi_6)
|
||||
dta$pase_0_bin<-factor(dta$pase_0_bin,levels=c("lower","higher"))
|
||||
|
||||
dta$nihss_0_isna<-is.na(dta$nihss_0)
|
||||
```
|
||||
|
||||
```{r}
|
||||
vars<-c("pase_0",
|
||||
"female",
|
||||
"age",
|
||||
"cohab",
|
||||
"ever_smoker",
|
||||
"diabetes",
|
||||
"hypertension",
|
||||
"afli",
|
||||
"ami",
|
||||
"tci",
|
||||
"pad",
|
||||
"nihss_0",
|
||||
"any_reperf")
|
||||
|
||||
# tbl1_vars<-c("thrombolysis", "thrombechtomy","inc_time")
|
||||
|
||||
labels_all<-list(active_treat~"Active trial treatment",
|
||||
pase_0~"PASE score",
|
||||
age~"Age",
|
||||
female~"Female sex",
|
||||
ever_smoker~"History of smoking",
|
||||
cohab~"Cohabitation",
|
||||
diabetes~"Known diabetes",
|
||||
hypertension~"Known hypertension",
|
||||
afli~"Known Atrialfibrillation",
|
||||
ami~"Previos myocardial infarction",
|
||||
tci~"Previos TIA",
|
||||
pad~"Known peripheral arteriotic disease",
|
||||
nihss_0~"NIHSS score",
|
||||
thrombolysis~"Thrombolytic therapy",
|
||||
thrombechtomy~"Endovascular treatment",
|
||||
any_reperf~"Any reperfusion treatment",
|
||||
inc_time~"Study inclusion time")
|
||||
|
||||
lab_sel<-function(label_list,variables_vector){
|
||||
## Helper function to select labels for gtsummary function from list of all labels based on selected variables.
|
||||
## Long names in try to ease reading.
|
||||
include_index<-c()
|
||||
for (i in 1:length(label_list)) {
|
||||
include_index[i]<-as.character(label_list[[i]])[2] %in% variables_vector
|
||||
}
|
||||
return(label_list[include_index])
|
||||
}
|
||||
```
|
||||
|
||||
# Table 1
|
||||
```{r}
|
||||
tbl1_vars<-c("active_treat",vars,"inc_time")
|
||||
|
||||
dta|>
|
||||
tbl_summary(missing = "ifany",
|
||||
include = all_of(tbl1_vars),
|
||||
missing_text="(Missing)",
|
||||
label = lab_sel(labels_all,tbl1_vars)
|
||||
)|>
|
||||
add_n()|>
|
||||
as_gt() |>
|
||||
# modify with gt functions
|
||||
gt::tab_header("Baseline Characteristics") |>
|
||||
gt::tab_options(
|
||||
table.font.size = "small",
|
||||
data_row.padding = gt::px(1))
|
||||
```
|
||||
|
||||
|
||||
# Regression - all
|
||||
|
||||
```{r}
|
||||
outs<-c("mdi_1_enr","mdi_6_newobs_enr")
|
||||
```
|
||||
|
||||
## Non-stratified
|
||||
|
||||
```{r}
|
||||
source("biv_mul.R")
|
||||
|
||||
bm_16_tbl_rtf <- file("bm_16_tbl.RTF", "w")
|
||||
writeLines(bm_16_tbl%>%as_gt()%>%as_rtf(), bm_16_tbl_rtf)
|
||||
close(bm_16_tbl_rtf)
|
||||
|
||||
bm_16_tbl %>% # build gtsummary table
|
||||
as_gt() %>% # convert to gt table
|
||||
gt::gtsave( # save table as image
|
||||
filename = "bm_16_tbl.png"
|
||||
)
|
||||
```
|
||||
|
||||
|
||||
## Stratified by treatment
|
||||
|
||||
```{r}
|
||||
source("biv_mul_strat.R")
|
||||
|
||||
bm_16_tbl_strat_rtf <- file("bm_16_tbl_strat.RTF", "w")
|
||||
writeLines(bm_16_tbl_strat%>%as_gt()%>%as_rtf(), bm_16_tbl_strat_rtf)
|
||||
close(bm_16_tbl_strat_rtf)
|
||||
|
||||
bm_16_tbl_strat %>% # build gtsummary table
|
||||
as_gt() %>% # convert to gt table
|
||||
gt::gtsave( # save table as image
|
||||
filename = "bm_16_tbl_strat.png"
|
||||
)
|
||||
```
|
||||
|
||||
# Regression - sensitivity
|
||||
|
||||
```{r}
|
||||
outs<-c("mdi_1","mdi_6")
|
||||
```
|
||||
|
||||
## Non-stratified
|
||||
|
||||
```{r}
|
||||
source("biv_mul.R")
|
||||
|
||||
bm_16_sens_tbl_rtf <- file("bm_16_sens_tbl.RTF", "w")
|
||||
writeLines(bm_16_tbl%>%as_gt()%>%as_rtf(), bm_16_sens_tbl_rtf)
|
||||
close(bm_16_sens_tbl_rtf)
|
||||
|
||||
bm_16_tbl %>% # build gtsummary table
|
||||
as_gt() %>% # convert to gt table
|
||||
gt::gtsave( # save table as image
|
||||
filename = "bm_16_tbl_sens.png"
|
||||
)
|
||||
```
|
||||
|
||||
|
||||
## Stratified by treatment
|
||||
|
||||
```{r}
|
||||
source("biv_mul_strat.R")
|
||||
|
||||
bm_16_tbl_strat_sens_rtf <- file("bm_16_tbl_strat_sens.RTF", "w")
|
||||
writeLines(bm_16_tbl_strat%>%as_gt()%>%as_rtf(), bm_16_tbl_strat_sens_rtf)
|
||||
close(bm_16_tbl_strat_sens_rtf)
|
||||
|
||||
bm_16_tbl_strat %>% # build gtsummary table
|
||||
as_gt() %>% # convert to gt table
|
||||
gt::gtsave( # save table as image
|
||||
filename = "bm_16_tbl_strat_sens.png"
|
||||
)
|
||||
```
|
||||
32
table_one.R
Normal file
32
table_one.R
Normal file
@ -0,0 +1,32 @@
|
||||
## =============================================================================
|
||||
## table_one.R
|
||||
##
|
||||
## Script by agdamsbo based on gtsummary.
|
||||
##
|
||||
## =============================================================================
|
||||
|
||||
library(readxl)
|
||||
library(dplyr)
|
||||
|
||||
## =============================================================================
|
||||
## Formatting
|
||||
## =============================================================================
|
||||
|
||||
# dta$Atrialfibrillation<-factor(ifelse(dta$AtrieflimrenDAP==1,1,0))
|
||||
# dta$Diabetes<-factor(ifelse(dta$DiabetesDAP==1,1,0))
|
||||
# dta$Hypertension<-factor(ifelse(dta$HypertensionDAP==1,1,0))
|
||||
# dta$Prev.TIA<-factor(ifelse(dta$TidlTci==1,1,0))
|
||||
# dta$Prev.AIS<-factor(ifelse(dta$TidlApo2==1,1,0))
|
||||
# dta$Smoker<-factor(ifelse(dta$Rygning==1,1,0))
|
||||
|
||||
|
||||
## =============================================================================
|
||||
## Table
|
||||
## =============================================================================
|
||||
|
||||
tbl1_tbl<-dta|>
|
||||
tbl_summary(missing = "ifany",
|
||||
include = all_of(tbl1_vars),
|
||||
missing_text="(Missing)"
|
||||
)|>
|
||||
add_n()
|
||||
Loading…
Reference in New Issue
Block a user