commit all: what a mess!

.gitignore

@@ -52,3 +52,4 @@ po/*~
 *.zip
 *.pdf
 *.csv
+/*.csv

00 master.R

@@ -10,9 +10,9 @@
 # 
 # Focus has been on universal functions, and the amount of flags has accumulated...
 # 
-# Only retransformation is performed on independent variables (these are not marked in the exported tables). 
+# No transformation is performed. Instead a sensitivity analysis is performed using splines.
-# Coefs and CIs are transformed according to this: 
+# 
-# https://stats.stackexchange.com/questions/93089/reporting-regression-statistics-after-logarithmic-transformation to allow for interpretation.
+# Robust CIs are calculated for all analyses.
 
 ## =============================================================================
 ## Table 1
@@ -24,6 +24,7 @@ source("table_1.R")
 
 tbl_1
 
+
 ## =============================================================================
 ## Primary regression analysis
 ## =============================================================================
@@ -40,15 +41,7 @@ biv_mul <- TRUE # Sets flag for both bivariate and multivariate or only multivar
 
 strat_var <- NULL # Variable to stratify by. Only one variable(!)
 
-trans_vars <- TRUE # Transform variables? T/F
+robust_ci <- TRUE
-
-sqrt_vars<-"pase_0" # Variables to sqrt-transfom
-log1p_vars<-"nihss_0" # Variables to log1p-transform, not outcome
-log_vars <- NULL # Variables to log-transform
-
-log1p_vars_all<-c(log1p_vars,outs) # All variables to log1p-transform, incl outcome
-
-trans_back <- TRUE # Back transform variables? T/F
 
 print_tbl <- TRUE # Print tables? T/F
 
@@ -85,15 +78,7 @@ biv_mul <- TRUE # Sets flag for both bivariate and multivariate or only multivar
 
 strat_var <- NULL # Variable to stratify by. Only one variable(!)
 
-trans_vars <- TRUE # Transform variables? T/F
+robust_ci  <- TRUE
-
-sqrt_vars<-"pase_0" # Variables to sqrt-transfom
-log1p_vars<-"nihss_0" # Variables to log1p-transform, not outcome
-log_vars <- NULL # Variables to log-transform
-
-log1p_vars_all<-c(log1p_vars,outs) # All variables to log1p-transform, incl outcome
-
-trans_back <- TRUE # Back transform variables? T/F
 
 print_tbl <- TRUE # Print tables? T/F
 
@@ -120,6 +105,21 @@ source("biv_mul_man.R")
 
 source("data_format.R")
 
+vars<-c("active_treat", # vars() is defined in data_format.R, here "
+        "pase_0", 
+        "female",
+        "age",
+        "cohab",
+        "ever_smoker",
+        "diabetes", 
+        "hypertension", 
+        "afli", 
+        "ami", 
+        "tci",
+        "pad",
+        "nihss_0",
+        "any_reperf")
+
 outs<-c("mdi_1_enr","mdi_6_newobs_enr")
 
 lbl_x <- "_inter" # Extra label in file name
@@ -128,17 +128,9 @@ inter_reg <- c("active_treat","pase_0") # Interaction variables to include (only
 
 biv_mul <- TRUE # Sets flag for both bivariate and multivariate or only multivariate analysis 
 
-trans_vars <- TRUE # Transform variables? T/F
-
 strat_var <- NULL # Variable to stratify by. Only one variable(!)
 
-sqrt_vars<-"pase_0" # Variables to sqrt-transfom
+robust_ci <- TRUE
-log1p_vars<-"nihss_0" # Variables to log1p-transform, not outcome
-log_vars <- NULL # Variables to log-transform
-
-log1p_vars_all<-c(log1p_vars,outs) # All variables to log1p-transform, incl outcome
-
-trans_back <- TRUE # Back transform variables? T/F
 
 print_tbl <- TRUE # Print tables? T/F
 
@@ -173,19 +165,11 @@ lbl_x <- "_nihss-pase" # Extra label in file name
 
 inter_reg <- NULL # Interaction variables to include (only multivariate)
 
-biv_mul <- FALSE # Sets flag for both bivariate and multivariate or only multivariate analysis 
+biv_mul <- TRUE # Sets flag for both bivariate and multivariate or only multivariate analysis 
-
-trans_vars <- TRUE # Transform variables? T/F
 
 strat_var <- NULL # Variable to stratify by. Only one variable(!)
 
-sqrt_vars<-"pase_0" # Variables to sqrt-transfom
+robust_ci <- TRUE
-log1p_vars<-NULL # Variables to log1p-transform, not outcome
-log_vars <- NULL # Variables to log-transform
-
-log1p_vars_all<-c(outs) # All variables to log1p-transform, incl outcome
-
-trans_back <- TRUE # Back transform variables? T/F
 
 print_tbl <- TRUE # Print tables? T/F
 

Archive/biv_mul_man.R

@@ -0,0 +1,118 @@
+## 
+## NOTES
+## 
+## Leaving gtsummary as this is slow and a little inflexible.
+## 
+
+## =============================================================================
+## Requirements
+## =============================================================================
+
+source("function_trans_cols.R")
+
+if (trans_vars==TRUE){
+  # If trans_vars flag is TRUE, transform specified variables
+  dta<-trans_cols(dta_backup,sqrts=sqrt_vars,log1ps = log1p_vars_all)
+} else {dta<-dta_backup}
+
+library(dplyr)
+library(gtsummary)
+
+source("function_back_trans.R")
+
+## =============================================================================
+## Loop
+## =============================================================================
+
+dec <- 3
+bm_list<-list()
+
+if (biv_mul){
+  do_biv=TRUE
+  do_mul=TRUE
+} else {
+  do_biv=FALSE
+  do_mul=TRUE
+}
+
+for (i in 1:length(outs)){
+
+  if (!is.null(strat_var)){
+    ls<-split(dta,dta[strat_var])
+    strat_list<-list()
+    
+    for (j in 1:length(ls)){
+      X<-ls[[j]][,vars] # Choosing the dataset for the given strata and stated variables
+      y <- ls[[j]][,outs[i]] # Choosing the outcome data for the given strata and the current outcome variable
+      
+      # Flagging factors and continous variables
+      # 
+      sel_f<-X|>
+        sapply(is.factor)
+      f.names <- colnames(X)[sel_f]
+      
+      source("function_reg_table.R")
+      rt<-reg_table(X,y,m.biv=do_biv,m.mul=do_mul,trans.var = trans_vars,sqrt.vars=sqrt_vars,inter.add=inter_reg)
+      
+      rt<-rt|>mutate(across(matches('co|lo|hi|pv'),as.numeric))
+      
+      rt[,grepl('co|lo|hi|pv',names(rt))] <- sapply( round( rt[,grepl('co|lo|hi|pv',names(rt))], 3 ),formatC, format='f', digits=3 )
+      
+      # Saving the data.frame in indexed list
+      strat_list[[j]] <- rt
+      
+      # Using just the first four characters for shorter colnames in "do.call" later.
+      names(strat_list)[j] <- paste0(substr(strat_var,1,4), 
+                                     ".",
+                                     levels(dta[,strat_var])[j])
+      
+    }
+    
+    bm_list[[i]] <- do.call(cbind, strat_list) # cbinds all elements in list, independent of number of strata
+    names(bm_list)[i] <- outs[i]
+    
+  } else {
+    X<-dta[,vars]
+    y <- dta[,outs[i]]
+    
+    # Flagging factors and continous variables
+    # 
+    sel_f<-X|>
+      sapply(is.factor)
+    f.names <- colnames(X)[sel_f]
+    
+    source("function_reg_table.R")
+    rt<-reg_table(X,y,m.biv=do_biv,m.mul=do_mul,trans.var = trans_vars,sqrt.vars=sqrt_vars,inter.add=inter_reg)
+    
+    rt<-rt|>mutate(across(matches('co|lo|hi|pv'),as.numeric))
+    
+    rt[,grepl('co|lo|hi|pv',names(rt))] <- sapply( round( rt[,grepl('co|lo|hi|pv',names(rt))], 3 ),formatC, format='f', digits=3 )
+    
+    bm_list[[i]] <- rt
+    names(bm_list)[i] <- outs[i]
+    
+  }
+  
+  }
+## =============================================================================
+## Print
+## =============================================================================
+
+# File name depending onstratification, transformation and back transformation
+if (print_tbl==TRUE){
+  fnm<-paste0("bm_16_tbl",lbl_x)
+  if (!is.null(strat_var)){fnm<-paste0(fnm,"_strat")}
+  if (trans_vars==TRUE){fnm<-paste0(fnm,"_trans")}
+  # if (trans_back==TRUE){fnm<-paste0(fnm,"_back")}
+  
+  export<-data.frame()
+  for (i in 1:length(bm_list)){
+    export<-rbind(export,paste(names(bm_list)[i]),bm_list[[i]])
+  }
+  
+  write.csv(export,paste0(fnm,".csv"))
+  
+}
+
+
+

Data skim.Rmd

@@ -0,0 +1,22 @@
+---
+title: "Data overview"
+author: "AGDamsbo"
+date: "2022-10-03"
+output:
+  html_document:
+    df_print: paged
+---
+
+```{r setup, include=FALSE}
+knitr::opts_chunk$set(echo = TRUE)
+```
+
+## Data skim
+
+
+```{r message=FALSE}
+source("data_format.R")
+
+dta |> skimr::skim(c(vars, "mdi_1_enr","mdi_6_newobs_enr", "active_treat"))
+```
+

NOK.Rmd

@@ -0,0 +1,108 @@
+---
+title: "NOK"
+author: "AGDamsbo"
+date: "`r Sys.Date()`"
+output: html_document
+---
+
+```{r setup, include=FALSE}
+knitr::opts_chunk$set(echo = TRUE)
+```
+
+## PASE by next of kin
+
+Hermed analyse og oversigt over TALOS patienter med PASE der ikke er udfyldt af dem selv.
+
+## Oversigt efter PASE gruppe
+
+Variablen "own_write" angiver om det af PASE-skemaet fremgår hvorvidt pt selv har udfyldt skemaet. Er TRUE/FALSE og kun TRUE vises.
+Variablen "own_answer" er min bedste vurdering udfra den registrerede kommentar, om hvorvidt det er pts eget svar (alternativt vurdering ved pårørende). Er TRUE/FALSE og kun TRUE vises.
+Variablen "afasi" angiver grad af afasi ved den akutte NIHSS.
+
+```{r echo=FALSE, message=FALSE, warning=FALSE}
+source("api_dataset.R")
+source("data_format.R")
+
+df <- dta|> left_join(df_own,by=c("rnumb"="record_id"))
+
+tbl_vars<-unique(c("own_write","own_answer","afasi","active_treat",vars,"inc_time","mdi_1_enr","mdi_6_newobs_enr","mdi_1","mdi_6_newobs"))
+```
+
+```{r}
+df |> 
+  tbl_summary(by = "pase_0_q",
+              missing = "ifany",
+              include = all_of(tbl_vars),
+              missing_text="(Missing)"#,
+              # label = lab_sel(labels_all,tbl1_vars)
+  )|>
+  add_n()|>
+  add_overall() |> 
+    gtsummary::italicize_levels() |> 
+  as_gt() |>
+  # modify with gt functions
+  gt::tab_header("Baseline Characteristics") |>
+  gt::tab_options(
+    table.font.size = "small",
+    data_row.padding = gt::px(1))
+```
+
+## Oversigt efter "own_answer" til sammenligning
+
+```{r}
+df |> #mutate(own_answer = own_answer+1) |> 
+  tbl_summary(by = "own_answer",
+              missing = "ifany",
+              include = all_of(tbl_vars),
+              missing_text="(Missing)"#,
+              # label = lab_sel(labels_all,tbl1_vars)
+  )|>
+  add_n()|>
+  add_overall() |> 
+  add_p() |> 
+  gtsummary::bold_p() |> 
+  gtsummary::italicize_levels() |> 
+  as_gt() |>
+  # modify with gt functions
+  gt::tab_header("PASE by next of kin") |>
+  gt::tab_options(
+    table.font.size = "small",
+    data_row.padding = gt::px(1))
+```
+
+## Sensititvitetsanalyse
+
+Begrænset til kun at inkludere ptt med eget svar
+
+```{r}
+source("data_format.R")
+
+dta <- inner_join(dta,df_own |> filter(own_answer) |> select(record_id),by=c("rnumb"="record_id"))
+
+outs<-c("mdi_1_enr","mdi_6_newobs_enr")
+
+lbl_x <- NULL # Extra label in file name
+
+inter_reg <- NULL # Interaction variables to include (only multivariate)
+
+biv_mul <- TRUE # Sets flag for both bivariate and multivariate or only multivariate analysis 
+
+print_tbl <- FALSE # Print tables? T/F
+
+robust_ci = TRUE
+
+# source("biv_mul.R")
+
+source("biv_mul_man.R")
+
+export<-data.frame()
+  for (i in 1:length(bm_list)){
+    export<-rbind(export,paste(names(bm_list)[i]),bm_list[[i]])
+  }
+
+export |> gt::gt()|>
+  gt::tab_options(
+    table.font.size = "small",
+    data_row.padding = gt::px(1))
+```
+

api_dataset.R

@@ -0,0 +1,38 @@
+# keyring::key_set("talos_api")
+
+
+df <-
+  REDCapR::redcap_read_oneshot(
+    "https://redcap.au.dk/api/",
+    keyring::key_get("talos_api"),forms = c("basis","reg","nihss_0","pase_0")
+  )$data |> 
+  dplyr::select(c(
+  "record_id",
+  "talos_nihss13_0",
+  "talos_nihss16_0",
+  "talos_pase12_0",
+  "talos_pase12x_0"
+)) 
+
+# For inspection of a few missing.
+# ds <- df[df$talos_pase12_0=="Not available",c("record_id","talos_pase10_0","talos_pase12_0")]
+
+
+library(tidyverse)
+
+# Grov kategorisering
+df |>
+  mutate(own = grepl("(interview)|(hjælp)|(adspurgt)", talos_pase12x_0)) |> 
+  arrange(desc(own)) |> readODS::write_ods("pase.ods")
+
+# Manual sorting
+
+df_own <- readODS::read_ods("pase_edit.ods") 
+
+df_own <- df_own |>
+  mutate(
+    own_write = talos_pase12_0=="1. Ja",
+    own_answer = ifelse(own_write, TRUE, own),
+    afasi = factor(talos_nihss13_0,labels=c("0:ingen","1:mild","2:svær","3:stum",NA)) 
+  ) |>
+  select(c("record_id", "own_write", "own_answer", "afasi"))

biv_mul_man.R

@@ -8,26 +8,21 @@
 ## Requirements
 ## =============================================================================
 
-source("function_trans_cols.R")
+source("function_reg_table_robust.R")
 
-if (trans_vars==TRUE){
+# dta<-dta_backup
-  # If trans_vars flag is TRUE, transform specified variables
-  dta<-trans_cols(dta_backup,sqrts=sqrt_vars,log1ps = log1p_vars_all)
-} else {dta<-dta_backup}
 
 library(dplyr)
 library(gtsummary)
+library(sandwich)
+library(lmtest)
 
-source("function_back_trans.R")
-
-vars_all<-c("active_treat",vars)
 
 ## =============================================================================
 ## Loop
 ## =============================================================================
 
 dec <- 3
-bm_list<-list()
 
 if (biv_mul){
   do_biv=TRUE
@@ -37,60 +32,28 @@ if (biv_mul){
   do_mul=TRUE
 }
 
-for (i in 1:length(outs)){
 
-  if (!is.null(strat_var)){
+bm_list <- lapply(seq_along(outs), function(i){
-    ls<-split(dta,dta[strat_var])
+  X<-dta[,vars]
-    strat_list<-list()
-    
-    for (j in 1:length(ls)){
-      X<-ls[[i]][,vars]
-      y <- ls[[i]][,outs[i]]
-      
-      # Flagging factors and continous variables
-      # 
-      sel_f<-X|>
-        sapply(is.factor)
-      f.names <- colnames(X)[sel_f]
-      
-      source("function_reg_table.R")
-      rt<-reg_table(X,y,m.biv=do_biv,m.mul=do_mul,trans.var = trans_vars,sqrt.vars=sqrt_vars,inter.add=inter_reg)
-      
-      rt<-rt|>mutate(across(matches('co|lo|hi|pv'),as.numeric))
-      
-      rt[,grepl('co|lo|hi|pv',names(rt))] <- sapply( round( rt[,grepl('co|lo|hi|pv',names(rt))], 3 ),formatC, format='f', digits=3 )
-      
-      strat_list[[j]] <- rt
-      names(strat_list)[j] <- levels(dta[,strat_var])[j]
-      
-    }
-    
-    bm_list[[i]] <- cbind(names(strat_list)[1],strat_list[[1]],names(strat_list)[2],strat_list[[2]])
-    names(bm_list)[i] <- outs[i]
-    
-  } else {
-    X<-dta[,vars_all]
   y <- dta[,outs[i]]
   
   # Flagging factors and continous variables
   # 
-    sel_f<-X|>
+  # sel_f <- X |>
-      sapply(is.factor)
+  #   sapply(is.factor)
-    f.names <- colnames(X)[sel_f]
+  # f.names <- colnames(X)[sel_f]
   
-    source("function_reg_table.R")
+  rt<-reg_table(X,y,m.biv=do_biv,m.mul=do_mul,inter.add=inter_reg,robust=robust_ci)
-    rt<-reg_table(X,y,m.biv=do_biv,m.mul=do_mul,trans.var = trans_vars,sqrt.vars=sqrt_vars,inter.add=inter_reg)
   
   rt<-rt|>mutate(across(matches('co|lo|hi|pv'),as.numeric))
   
   rt[,grepl('co|lo|hi|pv',names(rt))] <- sapply( round( rt[,grepl('co|lo|hi|pv',names(rt))], 3 ),formatC, format='f', digits=3 )
   
-    bm_list[[i]] <- rt
+  rt
-    names(bm_list)[i] <- outs[i]
+})
 
-  }
+names(bm_list) <- outs
 
-  }
 ## =============================================================================
 ## Print
 ## =============================================================================
@@ -98,9 +61,7 @@ for (i in 1:length(outs)){
 # File name depending onstratification, transformation and back transformation
 if (print_tbl==TRUE){
   fnm<-paste0("bm_16_tbl",lbl_x)
-  if (!is.null(strat_var)){fnm<-paste0(fnm,"_strat")}
+  if (robust_ci==TRUE){fnm<-paste0(fnm,"_robust")}
-  if (trans_vars==TRUE){fnm<-paste0(fnm,"_trans")}
-  # if (trans_back==TRUE){fnm<-paste0(fnm,"_back")}
   
   export<-data.frame()
   for (i in 1:length(bm_list)){

bm_16_tbl_inter_trans.csv

@@ -1,33 +0,0 @@
-"","name","pred","co","lo","hi","pv","co.p","lo.p","hi.p"
-"1","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr"
-"2","active_treat","active_treat.yes","0.234","-0.232","0.699","0.325","0.263","-0.207","1.013"
-"3","pase_0","pase_0","-0.000","-0.002","0.000","0.266","-0.000","-0.002","0.000"
-"4","female","female.yes","0.180","0.019","0.341","0.028","0.197","0.020","0.406"
-"5","age","age","-0.006","-0.013","0.000","0.059","-0.006","-0.013","0.000"
-"6","cohab","cohab.yes","0.039","-0.122","0.199","0.637","0.039","-0.115","0.220"
-"7","ever_smoker","ever_smoker.yes","-0.131","-0.288","0.027","0.103","-0.123","-0.250","0.027"
-"8","diabetes","diabetes.yes","0.218","-0.011","0.448","0.062","0.244","-0.011","0.565"
-"9","hypertension","hypertension.yes","0.062","-0.091","0.214","0.428","0.064","-0.087","0.239"
-"10","afli","afli.yes","-0.152","-0.353","0.050","0.140","-0.141","-0.298","0.051"
-"11","ami","ami.yes","0.175","-0.105","0.456","0.220","0.192","-0.100","0.577"
-"12","tci","tci.yes","0.075","-0.366","0.515","0.739","0.078","-0.306","0.674"
-"13","pad","pad.yes","0.150","-0.237","0.537","0.447","0.162","-0.211","0.712"
-"14","nihss_0","nihss_0","0.116","0.009","0.223","0.034","0.116","0.009","0.223"
-"15","any_reperf","any_reperf.yes","-0.011","-0.172","0.151","0.896","-0.011","-0.158","0.163"
-"16","active_treat*pase_0+","active_treat*pase_0+","-0.011","-0.049","0.027","0.570","-0.011","-0.049","0.027"
-"17","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr"
-"18","active_treat","active_treat.yes","-0.056","-0.624","0.511","0.845","-0.055","-0.464","0.667"
-"19","pase_0","pase_0","-0.001","-0.004","0.000","0.060","-0.001","-0.004","0.000"
-"20","female","female.yes","0.264","0.076","0.452","0.006","0.302","0.079","0.571"
-"21","age","age","-0.007","-0.015","0.000","0.063","-0.007","-0.015","0.000"
-"22","cohab","cohab.yes","-0.008","-0.198","0.182","0.933","-0.008","-0.179","0.199"
-"23","ever_smoker","ever_smoker.yes","-0.001","-0.183","0.182","0.996","-0.001","-0.168","0.200"
-"24","diabetes","diabetes.yes","-0.017","-0.299","0.265","0.906","-0.017","-0.259","0.304"
-"25","hypertension","hypertension.yes","0.100","-0.079","0.278","0.273","0.105","-0.076","0.321"
-"26","afli","afli.yes","-0.063","-0.301","0.175","0.604","-0.061","-0.260","0.191"
-"27","ami","ami.yes","0.303","-0.032","0.638","0.076","0.354","-0.032","0.894"
-"28","tci","tci.yes","-0.696","-1.211","-0.180","0.008","-0.501","-0.702","-0.164"
-"29","pad","pad.yes","0.304","-0.155","0.763","0.194","0.355","-0.144","1.144"
-"30","nihss_0","nihss_0","0.133","0.006","0.261","0.040","0.133","0.006","0.261"
-"31","any_reperf","any_reperf.yes","-0.035","-0.223","0.153","0.718","-0.034","-0.200","0.166"
-"32","active_treat*pase_0+","active_treat*pase_0+","-0.001","-0.047","0.045","0.965","-0.001","-0.047","0.045"

bm_16_tbl_nihss-pase_trans.csv

@@ -1,14 +0,0 @@
-"","name","pred","co","lo","hi","pv","co.p","lo.p","hi.p"
-"1","nihss_0","nihss_0","nihss_0","nihss_0","nihss_0","nihss_0","nihss_0","nihss_0","nihss_0"
-"2","active_treat","active_treat.yes","0.002","-0.119","0.123","0.975","0.002","-0.112","0.130"
-"3","pase_0","pase_0","-0.000","-0.001","0.000","0.067","-0.000","-0.001","0.000"
-"4","female","female.yes","0.086","-0.050","0.222","0.213","0.090","-0.049","0.249"
-"5","age","age","0.004","-0.001","0.010","0.127","0.004","-0.001","0.010"
-"6","cohab","cohab.yes","0.085","-0.049","0.220","0.215","0.089","-0.048","0.246"
-"7","ever_smoker","ever_smoker.yes","0.008","-0.124","0.140","0.908","0.008","-0.117","0.150"
-"8","diabetes","diabetes.yes","-0.000","-0.192","0.192","0.999","-0.000","-0.174","0.211"
-"9","hypertension","hypertension.yes","-0.149","-0.275","-0.022","0.021","-0.138","-0.240","-0.022"
-"10","afli","afli.yes","0.229","0.063","0.395","0.007","0.258","0.065","0.485"
-"11","ami","ami.yes","0.034","-0.195","0.263","0.771","0.035","-0.177","0.301"
-"12","tci","tci.yes","-0.170","-0.535","0.195","0.361","-0.156","-0.414","0.216"
-"13","pad","pad.yes","-0.273","-0.608","0.062","0.110","-0.239","-0.455","0.063"

bm_16_tbl_strat_trans.csv

@@ -1,29 +0,0 @@
-"","names(strat_list)[1]","name","pred","biv_co","biv_lo","biv_hi","biv_pv","biv_co.p","biv_lo.p","biv_hi.p","mul_co","mul_lo","mul_hi","mul_pv","mul_co.p","mul_lo.p","mul_hi.p","names(strat_list)[2]","name","pred","biv_co","biv_lo","biv_hi","biv_pv","biv_co.p","biv_lo.p","biv_hi.p","mul_co","mul_lo","mul_hi","mul_pv","mul_co.p","mul_lo.p","mul_hi.p"
-"1","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr"
-"2","no","pase_0","pase_0","-0.000","-0.002","0.000","0.226","-0.000","-0.002","0.000","-0.000","-0.002","0.000","0.364","-0.000","-0.002","0.000","yes","pase_0","pase_0","-0.000","-0.002","0.000","0.226","-0.000","-0.002","0.000","-0.000","-0.002","0.000","0.364","-0.000","-0.002","0.000"
-"3","no","female","female.yes","-0.032","-0.245","0.180","0.765","-0.032","-0.217","0.198","0.025","-0.204","0.255","0.830","0.025","-0.185","0.290","yes","female","female.yes","-0.032","-0.245","0.180","0.765","-0.032","-0.217","0.198","0.025","-0.204","0.255","0.830","0.025","-0.185","0.290"
-"4","no","age","age","0.001","-0.007","0.009","0.847","0.001","-0.007","0.009","-0.001","-0.010","0.008","0.768","-0.001","-0.010","0.008","yes","age","age","0.001","-0.007","0.009","0.847","0.001","-0.007","0.009","-0.001","-0.010","0.008","0.768","-0.001","-0.010","0.008"
-"5","no","cohab","cohab.yes","-0.028","-0.235","0.180","0.792","-0.027","-0.210","0.197","0.039","-0.188","0.266","0.736","0.040","-0.171","0.304","yes","cohab","cohab.yes","-0.028","-0.235","0.180","0.792","-0.027","-0.210","0.197","0.039","-0.188","0.266","0.736","0.040","-0.171","0.304"
-"6","no","ever_smoker","ever_smoker.yes","-0.115","-0.333","0.103","0.301","-0.109","-0.283","0.109","-0.114","-0.341","0.113","0.323","-0.108","-0.289","0.119","yes","ever_smoker","ever_smoker.yes","-0.115","-0.333","0.103","0.301","-0.109","-0.283","0.109","-0.114","-0.341","0.113","0.323","-0.108","-0.289","0.119"
-"7","no","diabetes","diabetes.yes","0.125","-0.162","0.412","0.392","0.133","-0.149","0.509","0.176","-0.130","0.482","0.259","0.192","-0.122","0.619","yes","diabetes","diabetes.yes","0.125","-0.162","0.412","0.392","0.133","-0.149","0.509","0.176","-0.130","0.482","0.259","0.192","-0.122","0.619"
-"8","no","hypertension","hypertension.yes","0.064","-0.134","0.262","0.525","0.066","-0.126","0.300","-0.054","-0.272","0.164","0.625","-0.053","-0.238","0.178","yes","hypertension","hypertension.yes","0.064","-0.134","0.262","0.525","0.066","-0.126","0.300","-0.054","-0.272","0.164","0.625","-0.053","-0.238","0.178"
-"9","no","afli","afli.yes","-0.061","-0.337","0.216","0.665","-0.059","-0.286","0.241","-0.126","-0.418","0.167","0.397","-0.118","-0.342","0.181","yes","afli","afli.yes","-0.061","-0.337","0.216","0.665","-0.059","-0.286","0.241","-0.126","-0.418","0.167","0.397","-0.118","-0.342","0.181"
-"10","no","ami","ami.yes","0.230","-0.116","0.575","0.192","0.258","-0.110","0.778","0.348","-0.044","0.741","0.082","0.416","-0.043","1.097","yes","ami","ami.yes","0.230","-0.116","0.575","0.192","0.258","-0.110","0.778","0.348","-0.044","0.741","0.082","0.416","-0.043","1.097"
-"11","no","tci","tci.yes","-0.371","-1.065","0.324","0.295","-0.310","-0.655","0.383","-0.356","-1.044","0.332","0.309","-0.300","-0.648","0.393","yes","tci","tci.yes","-0.371","-1.065","0.324","0.295","-0.310","-0.655","0.383","-0.356","-1.044","0.332","0.309","-0.300","-0.648","0.393"
-"12","no","pad","pad.yes","0.321","-0.111","0.752","0.145","0.378","-0.105","1.122","0.287","-0.195","0.770","0.241","0.333","-0.177","1.159","yes","pad","pad.yes","0.321","-0.111","0.752","0.145","0.378","-0.105","1.122","0.287","-0.195","0.770","0.241","0.333","-0.177","1.159"
-"13","no","nihss_0","nihss_0","0.153","0.010","0.297","0.036","0.153","0.010","0.297","0.238","0.083","0.393","0.003","0.238","0.083","0.393","yes","nihss_0","nihss_0","0.153","0.010","0.297","0.036","0.153","0.010","0.297","0.238","0.083","0.393","0.003","0.238","0.083","0.393"
-"14","no","any_reperf","any_reperf.yes","-0.087","-0.289","0.114","0.396","-0.083","-0.251","0.121","-0.194","-0.417","0.029","0.089","-0.176","-0.341","0.030","yes","any_reperf","any_reperf.yes","-0.087","-0.289","0.114","0.396","-0.083","-0.251","0.121","-0.194","-0.417","0.029","0.089","-0.176","-0.341","0.030"
-"15","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr"
-"16","no","pase_0","pase_0","-0.002","-0.006","-0.000","0.007","-0.002","-0.006","-0.000","-0.002","-0.006","0.000","0.059","-0.002","-0.006","0.000","yes","pase_0","pase_0","-0.002","-0.006","-0.000","0.007","-0.002","-0.006","-0.000","-0.002","-0.006","0.000","0.059","-0.002","-0.006","0.000"
-"17","no","female","female.yes","0.305","0.061","0.549","0.015","0.356","0.063","0.731","0.274","0.002","0.546","0.048","0.315","0.002","0.726","yes","female","female.yes","0.305","0.061","0.549","0.015","0.356","0.063","0.731","0.274","0.002","0.546","0.048","0.315","0.002","0.726"
-"18","no","age","age","0.000","-0.009","0.009","0.977","0.000","-0.009","0.009","-0.008","-0.019","0.003","0.134","-0.008","-0.019","0.003","yes","age","age","0.000","-0.009","0.009","0.977","0.000","-0.009","0.009","-0.008","-0.019","0.003","0.134","-0.008","-0.019","0.003"
-"19","no","cohab","cohab.yes","-0.067","-0.325","0.190","0.607","-0.065","-0.277","0.209","-0.102","-0.381","0.178","0.475","-0.097","-0.317","0.195","yes","cohab","cohab.yes","-0.067","-0.325","0.190","0.607","-0.065","-0.277","0.209","-0.102","-0.381","0.178","0.475","-0.097","-0.317","0.195"
-"20","no","ever_smoker","ever_smoker.yes","0.103","-0.145","0.351","0.413","0.109","-0.135","0.420","0.098","-0.163","0.359","0.460","0.103","-0.150","0.432","yes","ever_smoker","ever_smoker.yes","0.103","-0.145","0.351","0.413","0.109","-0.135","0.420","0.098","-0.163","0.359","0.460","0.103","-0.150","0.432"
-"21","no","diabetes","diabetes.yes","0.080","-0.339","0.498","0.709","0.083","-0.287","0.645","0.079","-0.343","0.500","0.713","0.082","-0.290","0.649","yes","diabetes","diabetes.yes","0.080","-0.339","0.498","0.709","0.083","-0.287","0.645","0.079","-0.343","0.500","0.713","0.082","-0.290","0.649"
-"22","no","hypertension","hypertension.yes","0.285","0.047","0.523","0.019","0.330","0.049","0.686","0.158","-0.103","0.419","0.235","0.171","-0.098","0.520","yes","hypertension","hypertension.yes","0.285","0.047","0.523","0.019","0.330","0.049","0.686","0.158","-0.103","0.419","0.235","0.171","-0.098","0.520"
-"23","no","afli","afli.yes","-0.079","-0.404","0.245","0.631","-0.076","-0.332","0.278","-0.076","-0.411","0.259","0.655","-0.073","-0.337","0.296","yes","afli","afli.yes","-0.079","-0.404","0.245","0.631","-0.076","-0.332","0.278","-0.076","-0.411","0.259","0.655","-0.073","-0.337","0.296"
-"24","no","ami","ami.yes","0.203","-0.254","0.661","0.382","0.225","-0.225","0.936","0.178","-0.326","0.682","0.488","0.195","-0.278","0.977","yes","ami","ami.yes","0.203","-0.254","0.661","0.382","0.225","-0.225","0.936","0.178","-0.326","0.682","0.488","0.195","-0.278","0.977"
-"25","no","tci","tci.yes","-0.718","-1.431","-0.005","0.049","-0.512","-0.761","-0.005","-0.806","-1.521","-0.092","0.027","-0.553","-0.781","-0.088","yes","tci","tci.yes","-0.718","-1.431","-0.005","0.049","-0.512","-0.761","-0.005","-0.806","-1.521","-0.092","0.027","-0.553","-0.781","-0.088"
-"26","no","pad","pad.yes","0.339","-0.439","1.118","0.391","0.404","-0.355","2.058","0.293","-0.575","1.161","0.507","0.340","-0.437","2.194","yes","pad","pad.yes","0.339","-0.439","1.118","0.391","0.404","-0.355","2.058","0.293","-0.575","1.161","0.507","0.340","-0.437","2.194"
-"27","no","nihss_0","nihss_0","0.028","-0.136","0.193","0.734","0.028","-0.136","0.193","-0.056","-0.239","0.126","0.544","-0.056","-0.239","0.126","yes","nihss_0","nihss_0","0.028","-0.136","0.193","0.734","0.028","-0.136","0.193","-0.056","-0.239","0.126","0.544","-0.056","-0.239","0.126"
-"28","no","any_reperf","any_reperf.yes","0.037","-0.211","0.285","0.768","0.038","-0.190","0.330","0.183","-0.097","0.462","0.199","0.201","-0.092","0.588","yes","any_reperf","any_reperf.yes","0.037","-0.211","0.285","0.768","0.038","-0.190","0.330","0.183","-0.097","0.462","0.199","0.201","-0.092","0.588"

bm_16_tbl_trans.csv

@@ -1,31 +1,33 @@
-"","name","pred","biv_co","biv_lo","biv_hi","biv_pv","biv_co.p","biv_lo.p","biv_hi.p","mul_co","mul_lo","mul_hi","mul_pv","mul_co.p","mul_lo.p","mul_hi.p"
+"","name","co.x","biv_lo","biv_hi","biv_pv","co.y","mul_lo","mul_hi","mul_pv"
-"1","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr"
+"1","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr","mdi_1_enr"
-"2","active_treat","active_treat.yes","0.053","-0.090","0.196","0.464","0.055","-0.086","0.217","0.105","-0.039","0.250","0.151","0.111","-0.038","0.284"
+"2","pase_0_qq_2","-0.589","-2.679","1.500","0.580","-0.036","-2.199","2.127","0.974"
-"3","pase_0","pase_0","-0.001","-0.002","-0.000","0.006","-0.001","-0.002","-0.000","-0.000","-0.002","0.000","0.059","-0.000","-0.002","0.000"
+"3","pase_0_qq_3","-2.279","-4.151","-0.406","0.017","-2.361","-4.189","-0.532","0.011"
-"4","female","female.yes","0.145","-0.005","0.295","0.058","0.156","-0.005","0.343","0.183","0.023","0.344","0.025","0.201","0.023","0.410"
+"4","pase_0_qq_4","-3.019","-4.683","-1.355","0.000","-2.382","-4.285","-0.480","0.014"
-"5","age","age","-0.001","-0.007","0.004","0.688","-0.001","-0.007","0.004","-0.006","-0.012","0.000","0.060","-0.006","-0.012","0.000"
+"5","femaleyes","1.392","-0.017","2.801","0.053","1.680","0.211","3.150","0.025"
-"6","cohab","cohab.yes","0.019","-0.132","0.169","0.807","0.019","-0.123","0.184","0.040","-0.120","0.201","0.620","0.041","-0.113","0.222"
+"6","age","0.003","-0.050","0.056","0.903","-0.059","-0.119","0.001","0.055"
-"7","ever_smoker","ever_smoker.yes","-0.088","-0.240","0.064","0.256","-0.084","-0.214","0.066","-0.134","-0.291","0.023","0.095","-0.125","-0.252","0.023"
+"7","cohabyes","-0.001","-1.393","1.390","0.998","0.415","-1.044","1.874","0.576"
-"8","diabetes","diabetes.yes","0.159","-0.065","0.383","0.163","0.172","-0.063","0.466","0.217","-0.013","0.446","0.064","0.242","-0.012","0.563"
+"8","ever_smokeryes","-0.904","-2.286","0.479","0.200","-0.798","-2.129","0.534","0.240"
-"9","hypertension","hypertension.yes","0.138","-0.004","0.281","0.058","0.148","-0.004","0.324","0.061","-0.091","0.214","0.430","0.063","-0.087","0.238"
+"9","diabetesyes","0.949","-1.098","2.996","0.363","0.959","-1.041","2.960","0.347"
-"10","afli","afli.yes","-0.146","-0.342","0.050","0.143","-0.136","-0.290","0.051","-0.149","-0.351","0.052","0.146","-0.139","-0.296","0.053"
+"10","hypertensionyes","1.395","0.089","2.702","0.036","0.636","-0.677","1.948","0.342"
-"11","ami","ami.yes","0.147","-0.112","0.405","0.266","0.158","-0.106","0.499","0.177","-0.103","0.457","0.216","0.194","-0.098","0.580"
+"11","afliyes","-1.047","-2.746","0.652","0.227","-1.185","-2.807","0.438","0.152"
-"12","tci","tci.yes","0.138","-0.316","0.591","0.550","0.148","-0.271","0.807","0.078","-0.362","0.518","0.729","0.081","-0.304","0.678"
+"12","amiyes","2.527","-0.491","5.545","0.101","2.150","-0.941","5.240","0.172"
-"13","pad","pad.yes","0.139","-0.219","0.497","0.447","0.149","-0.197","0.643","0.150","-0.237","0.537","0.446","0.162","-0.211","0.711"
+"13","tciyes","2.151","-4.206","8.508","0.507","2.137","-4.254","8.528","0.512"
-"14","nihss_0","nihss_0","0.106","0.007","0.204","0.036","0.106","0.007","0.204","0.117","0.010","0.224","0.032","0.117","0.010","0.224"
+"14","padyes","3.275","-1.509","8.058","0.179","4.002","-1.489","9.493","0.153"
-"15","any_reperf","any_reperf.yes","0.021","-0.127","0.168","0.785","0.021","-0.119","0.183","-0.013","-0.174","0.149","0.879","-0.012","-0.160","0.160"
+"15","nihss_0","0.176","0.032","0.320","0.017","0.182","0.032","0.332","0.018"
-"16","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr"
+"16","any_reperfyes","0.151","-1.191","1.493","0.825","-0.236","-1.643","1.171","0.742"
-"17","active_treat","active_treat.yes","-0.084","-0.249","0.081","0.316","-0.081","-0.220","0.084","-0.068","-0.237","0.100","0.426","-0.066","-0.211","0.105"
+"17","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr","mdi_6_newobs_enr"
-"18","pase_0","pase_0","-0.001","-0.003","-0.000","0.002","-0.001","-0.003","-0.000","-0.001","-0.003","-0.000","0.017","-0.001","-0.003","-0.000"
+"18","pase_0_qq_2","-1.939","-4.202","0.325","0.093","-1.469","-3.817","0.880","0.220"
-"19","female","female.yes","0.246","0.074","0.418","0.005","0.279","0.077","0.519","0.264","0.077","0.451","0.006","0.302","0.080","0.570"
+"19","pase_0_qq_3","-3.334","-5.389","-1.279","0.002","-3.346","-5.484","-1.207","0.002"
-"20","age","age","0.001","-0.005","0.008","0.658","0.001","-0.005","0.008","-0.007","-0.015","0.000","0.062","-0.007","-0.015","0.000"
+"20","pase_0_qq_4","-3.218","-5.354","-1.083","0.003","-2.767","-5.190","-0.344","0.025"
-"21","cohab","cohab.yes","-0.094","-0.269","0.080","0.290","-0.090","-0.236","0.084","-0.008","-0.197","0.181","0.935","-0.008","-0.179","0.199"
+"21","femaleyes","2.350","0.768","3.932","0.004","2.266","0.547","3.986","0.010"
-"22","ever_smoker","ever_smoker.yes","-0.021","-0.197","0.155","0.817","-0.020","-0.178","0.168","-0.001","-0.183","0.182","0.994","-0.001","-0.167","0.199"
+"22","age","0.009","-0.054","0.072","0.779","-0.054","-0.122","0.013","0.115"
-"23","diabetes","diabetes.yes","0.012","-0.260","0.284","0.932","0.012","-0.229","0.329","-0.017","-0.299","0.265","0.905","-0.017","-0.259","0.303"
+"23","cohabyes","-1.162","-2.738","0.414","0.148","-0.409","-2.032","1.213","0.620"
-"24","hypertension","hypertension.yes","0.232","0.068","0.396","0.006","0.261","0.070","0.487","0.100","-0.078","0.278","0.272","0.105","-0.075","0.321"
+"24","ever_smokeryes","0.011","-1.502","1.525","0.988","0.241","-1.259","1.741","0.753"
-"25","afli","afli.yes","-0.025","-0.255","0.206","0.835","-0.024","-0.225","0.229","-0.063","-0.300","0.175","0.605","-0.061","-0.259","0.191"
+"25","diabetesyes","-0.343","-2.379","1.694","0.741","-0.544","-2.598","1.511","0.603"
-"26","ami","ami.yes","0.200","-0.103","0.502","0.195","0.221","-0.098","0.652","0.303","-0.031","0.638","0.076","0.354","-0.031","0.893"
+"26","hypertensionyes","1.357","-0.047","2.761","0.058","0.592","-0.914","2.099","0.440"
-"27","tci","tci.yes","-0.660","-1.180","-0.140","0.013","-0.483","-0.693","-0.130","-0.695","-1.210","-0.180","0.008","-0.501","-0.702","-0.165"
+"27","afliyes","-0.092","-2.074","1.889","0.927","-0.682","-2.616","1.252","0.488"
-"28","pad","pad.yes","0.416","0.001","0.831","0.050","0.516","0.001","1.296","0.304","-0.155","0.762","0.194","0.355","-0.143","1.143"
+"28","amiyes","1.483","-1.174","4.141","0.273","2.282","-0.716","5.279","0.135"
-"29","nihss_0","nihss_0","0.142","0.024","0.260","0.018","0.142","0.024","0.260","0.134","0.007","0.260","0.039","0.134","0.007","0.260"
+"29","tciyes","-4.472","-6.447","-2.498","0.000","-4.458","-6.453","-2.462","0.000"
-"30","any_reperf","any_reperf.yes","-0.028","-0.198","0.142","0.749","-0.027","-0.179","0.153","-0.035","-0.222","0.153","0.715","-0.034","-0.199","0.165"
+"30","padyes","2.129","-1.510","5.769","0.251","1.461","-2.690","5.611","0.490"
+"31","nihss_0","0.225","0.045","0.404","0.014","0.208","0.009","0.407","0.040"
+"32","any_reperfyes","-0.036","-1.480","1.409","0.961","-0.369","-1.839","1.101","0.622"

data_format.R

@@ -10,7 +10,7 @@ library(dplyr)
 library(tidyr)
 
 
-dta_all<-read.csv("/Volumes/Data/depression/dep_dataset.csv")
+dta_all<-read.csv("/Volumes/Data 1/depression/dep_dataset.csv")
 
 
 # Defining patients to include for analysis
@@ -62,9 +62,9 @@ dta$pase_0_bin<-factor(dta$pase_0_bin,levels=c("lower","higher"))
 
 dta$nihss_0_isna<-is.na(dta$nihss_0)
 
+dta$pase_0_q<-factor(dta$pase_0_q)
 
-
+vars<-c("pase_0_q",
-vars<-c("pase_0",
         "female",
         "age",
         "cohab",
@@ -76,11 +76,13 @@ vars<-c("pase_0",
         "tci",
         "pad",
         "nihss_0",
-        "any_reperf")
+        "any_reperf",
+        "active_treat")
 
 # tbl1_vars<-c("thrombolysis", "thrombechtomy","inc_time")
 
 labels_all<-list(active_treat~"Active trial treatment",
+                 pase_0_q~"PASE score quartile",
                  pase_0~"PASE score",
                  age~"Age",
                  female~"Female sex",

experiments_splines.R

@@ -0,0 +1,89 @@
+# Experiments
+
+library("sandwich")
+library("lmtest")
+
+vars<-c("pase_0", # New variables for analysis
+        "female",
+        "age",
+        "cohab",
+        "ever_smoker",
+        "diabetes", 
+        "hypertension", 
+        "afli", 
+        "ami", 
+        "tci",
+        "pad")
+
+m<-lm(dta$mdi_1_enr~., data = dta[,vars])
+summary(m)
+coeftest(m, vcov = sandwich)
+coefci(m, vcov = sandwich)
+
+source("function_reg_table_robust.R")
+
+reg_table(X = dta[,vars], y = dta[,"mdi_1_enr"]) |> gt()
+
+## Så virker funktionen. Indbygges i script som tidligere.
+
+
+vars<-c("pase_0", # New variables for analysis
+        "female",
+        "age",
+        "cohab",
+        "ever_smoker",
+        "diabetes", 
+        "hypertension", 
+        "afli", 
+        "ami", 
+        "tci",
+        "pad")
+quants <- quantile(dta$pase_0,na.rm=TRUE)
+
+m<-lm(dta$mdi_1_enr~splines::ns(pase_0, knots = quants[2:4])+., data = dta[,vars])
+
+coeftest(m, vcov = sandwich) |> tbl_regression() |> bold_p()
+
+
+p1 <- ggplot(dta) + 
+  aes(pase_0,mdi_1_enr) +
+  geom_point(na.rm = TRUE) +
+  geom_smooth(method='lm', formula= y~x, na.rm = TRUE)+
+  labs(title="Linear fit", x = "PASE score", y = "MDI score 1 month")
+
+p2 <- ggplot() +
+  aes(m$model$pase_0, predict(m)) +
+  geom_point(color="blue") +
+  geom_smooth(method='lm',formula = y~splines::ns(x, knots = quants[2:4])) +
+  geom_vline(xintercept=quants) +
+  ylim(0,max(dta$mdi_1_enr)) +
+  labs(title="Qubic splines with predicted data", x = "PASE score", y = "MDI score 1 month")
+
+library(patchwork)
+p1+p2
+
+# p2 <- ggplot() +
+#   aes(m$model$pase_0, predict(m)) +
+#   geom_point(color="maroon",shape=25) +
+#   geom_smooth(method='lm',formula = y~splines::ns(x, knots = quants[2:4])) +
+#   geom_vline(xintercept=quants) +
+#   # ylim(0,max(dta$mdi_1_enr)) +
+#   labs(title="Qubic splines with predicted data", x = "PASE score", y = "MDI score 1 month")+
+#   theme_bw(14)
+
+library(npreg)
+
+dta_c <- dta[,c("mdi_1_enr","pase_0")] |> na.omit()
+
+x <- dta_c$pase_0
+y <- dta_c$mdi_1_enr
+
+m1<-lm(dta$mdi_1_enr~dta$pase_0)
+
+
+plot(ss(x, y, knots=quants), xlab = "PASE score", ylab = "MDI score 1 month", ylim=c(min(y),max(y)))
+# rug(x)
+points(x,y)
+abline(coef(m1), lty = 2)
+legend("topright", legend = c("ss", "lm"), lty = 1:2, bty = "n")
+

function_reg_table.R

@@ -42,6 +42,7 @@ for (j in colnames(X)){
   }
   
   df_b<-df_b|>data.frame()|>mutate(across(matches('co|lo|hi'),as.numeric))
+  colnames(df_b)[3:ncol(df_b)]<-paste0("biv_",colnames(df_b)[3:ncol(df_b)])
 }
   
   if (m.mul){
@@ -70,15 +71,11 @@ for (j in colnames(X)){
       }
       
       df_m<-df_m|>data.frame()|>mutate(across(matches('co|lo|hi'),as.numeric))
+      colnames(df_m)[3:ncol(df_m)]<-paste0("mul_",colnames(df_m)[3:ncol(df_m)])
       
   }
 
-  
-  
   if (all(m.biv,m.mul)){
-    colnames(df_b)[3:ncol(df_b)]<-paste0("biv_",colnames(df_b)[3:ncol(df_b)])
-    colnames(df_m)[3:ncol(df_m)]<-paste0("mul_",colnames(df_m)[3:ncol(df_m)])
-    
     return(merge(df_b,df_m,by=c("name","pred"),sort=FALSE))
   } else if (m.biv) {
     return(df_b)

function_reg_table_robust.R

@@ -0,0 +1,74 @@
+reg_table <- function(X,y,m.biv=TRUE,m.mul=TRUE,inter.add=NULL,robust=TRUE){
+  # method One of biv, mul, biv_mul
+  
+  # X = dta[,vars]
+  # y = dta[,"mdi_1_enr"]
+  
+  cols<-c("name", "co", "lo", "hi", "pv")
+  
+  if(!is.null(inter.add)){
+    form_add<-paste0(paste0(inter.add,collapse = "*"),"+")
+    m.biv <- FALSE # If interaction term is added, only multivariate is performed.
+  } else {
+    form_add=NULL
+  }
+  
+if (m.biv){
+  df_b <- data.frame(matrix(NA, ncol = length(cols)))
+  names(df_b)<-cols
+  
+  
+for (j in colnames(X)){
+  m<-lm(formula(paste0("y~",j)),X)
+  if(robust){
+    m <- coeftest(m, vcov = sandwich)
+    pv <- m[-1,4] # Avoids dependency
+  } else {
+      pv <- summary(m)$coefficients[-1,4] # Avoids dependency
+      }
+  
+  ci <- confint(m)
+  lo <- ci[-1, 1]
+  hi <- ci[-1, 2]
+  co <- coef(m)[-1]
+
+  df_b <- rbind(df_b, cbind(name=names(co), co, lo, hi, pv))
+  
+}
+  
+  df_b <- df_b[-1,]
+  
+  df_b<-df_b|>data.frame()|>mutate(across(matches('co|lo|hi'),as.numeric))
+  colnames(df_b)[3:ncol(df_b)]<-paste0("biv_",colnames(df_b)[3:ncol(df_b)])
+}
+  
+  if (m.mul){
+    m<-lm(formula(paste0("y~",form_add,".")),X)
+    if(robust){
+      m <- coeftest(m, vcov = sandwich)
+      pv <- m[-1,4] # Avoids dependency
+    } else {
+      pv <- summary(m)$coefficients[-1,4] # Avoids dependency
+      }
+      
+      ci <- confint(m)
+      lo <- ci[-1, 1]
+      hi <- ci[-1, 2]
+      co <- coef(m)[-1]
+
+      df_m <- cbind(name=c(names(co)), co, lo, hi, pv)
+
+      df_m<-df_m|>data.frame()|>mutate(across(matches('co|lo|hi'),as.numeric))
+      colnames(df_m)[3:ncol(df_m)]<-paste0("mul_",colnames(df_m)[3:ncol(df_m)])
+      
+  }
+
+  if (all(m.biv,m.mul)){
+    merge(df_b,df_m,by=c("name"),sort=FALSE)
+  } else if (m.biv) {
+    return(df_b)
+  } else {
+      df_m
+    }
+  
+}

master.Rmd

@@ -0,0 +1,362 @@
+---
+title: "Master file for analysis"
+author: "Andreas Gammelgaard Damsbo"
+date: "`r format(Sys.time(), '%d.%b.%Y')`"
+toc: true
+output:
+  html_document:
+    df_print: paged
+---
+
+# Noter
+
+Ny analyse efter snak med Jan 14.10.22
+
+Plan:
+
+- Analyser med rebust variance estimering (sandwich eller clubSandwich) for at komme omkring transformation af outcome.
+- Opdel primær exposure (PASE0) i fire gupper
+- Suppler med spline-k analyse med fire knuder
+  - Behøver visuel tolkning. Jan bruger natural splines. splines::ns()
+- Supplerende analyse med 
+
+
+# Overview
+```{r}
+source("data_format.R")
+
+dta |> skimr::skim(c(vars, "mdi_1_enr","mdi_6_newobs_enr", "active_treat", "pase_0"))
+```
+
+
+# Table 1
+
+```{r}
+source("data_format.R")
+
+source("table_1.R")
+
+tbl_1
+
+write.csv(tbl_1[["_data"]][c(3,4,6:10)],"table_1_quarts.csv")
+```
+
+
+
+# Primære analyse
+
+```{r}
+## =============================================================================
+## Primary regression analysis
+## =============================================================================
+
+source("data_format.R")
+
+outs<-c("mdi_1_enr","mdi_6_newobs_enr")
+
+lbl_x <- NULL # Extra label in file name
+
+inter_reg <- NULL # Interaction variables to include (only multivariate)
+
+biv_mul <- TRUE # Sets flag for both bivariate and multivariate or only multivariate analysis 
+
+robust_ci <- TRUE
+
+print_tbl <- TRUE # Print tables? T/F
+
+source("biv_mul_man.R")
+
+export
+```
+
+
+
+# Spline analyse
+
+## One month
+
+```{r message=FALSE, warning=FALSE}
+library("sandwich")
+library("lmtest")
+
+source("data_format.R")
+
+vars<-c("pase_0", # New variables for analysis
+        "female",
+        "age",
+        "cohab",
+        "ever_smoker",
+        "diabetes", 
+        "hypertension", 
+        "afli", 
+        "ami", 
+        "tci",
+        "pad",
+        "active_treat")
+quants <- quantile(dta$pase_0,na.rm=TRUE)
+
+m<-lm(dta$mdi_1_enr~splines::ns(pase_0, knots = quants[2:4])+., data = dta[,vars])
+
+coeftest(m, vcov = sandwich) |> tbl_regression() |> bold_p()
+```
+
+
+## Six months
+
+```{r}
+m<-lm(dta$mdi_6_newobs_enr~splines::ns(pase_0, knots = quants[2:4])+., data = dta[,vars])
+
+coeftest(m, vcov = sandwich) |> tbl_regression() |> bold_p()
+```
+
+## Ploting
+
+```{r}
+library(ggplot2)
+# dta[,vars]
+spline_plot <- function(ds=dta, 
+                        x.n="pase_0", 
+                        y.n, 
+                        t.lab = "Model fit", 
+                        x.lab = "PASE score", 
+                        y.lab="MDI score", 
+l.size = .7,
+l.color = "black",
+l.type = c("dotdash", "solid"),
+l.lab = c("Linear", "Cubic \nsplines"),
+p.color = "grey50",
+p.shape=20,
+p.size=.5,
+s.knots=quants,
+th.size=12, 
+m.var=NULL){
+  p <- ggplot() + 
+  aes(ds[,x.n],ds[,y.n]) +
+  geom_point(na.rm = TRUE,color=p.color,shape=p.shape, size=p.size) +
+  geom_smooth(method='lm', formula= y~x, na.rm = TRUE, se = FALSE, aes(linetype = "A"),size=l.size, color=l.color) +
+  geom_vline(xintercept=s.knots[2:(length(quants)-1)], linetype = "dashed") +
+  labs(title=t.lab, x = x.lab, y = y.lab)+
+  theme_bw(th.size)
+  
+  if (!is.null(m.var)){
+    m <- lm(ds[,y.n]~splines::ns(ds[,x.n], knots = s.knots[2:(length(quants)-1)])+., data = ds[,m.var])
+    
+    df <- cbind(ds[x.n],stats::predict(m, newdata=ds, interval="confidence", method="lm"))
+    
+    p + geom_smooth(data=df, aes(x=df[,x.n], y=fit, ymin=lwr, ymax=upr, linetype = "B"), size = l.size, colour = l.color, se = TRUE, stat = "smooth", na.rm=TRUE)+
+  scale_linetype_manual(name="Fit lines", values=l.type, labels = l.lab)
+  } else {
+    p+
+  geom_smooth(method='lm',formula = y~splines::ns(x, knots = s.knots[2:(length(quants)-1)]), na.rm = TRUE, aes(linetype = "B"),size=l.size, color=l.color)+
+  scale_linetype_manual(name="Fit lines", values=l.type, labels = l.lab)
+  }
+}
+```
+
+```{r}
+ds=dta
+x.n="pase_0"
+y.n="mdi_1_enr"
+t.lab = "Model fit"
+x.lab = "PASE score"
+y.lab="MDI score"
+l.size = .7
+l.color = "black"
+l.type = c("dotdash", "solid")
+l.lab = c("Linear", "Cubic \nsplines")
+p.color = "grey50"
+p.shape=20
+p.size=.5
+s.knots=quants
+th.size=12
+m.var=c(vars,"mdi_1_enr")
+```
+
+
+```{r out.width="100%"}
+vars_all <- c(vars,"mdi_1_enr")
+
+p1 <- dta |> select(c("mdi_1_enr",vars)) |> spline_plot(y.n="mdi_1_enr",t.lab="One month")
+p2 <- spline_plot(y.n="mdi_6_newobs_enr",t.lab="Six months")
+
+dta |> select(c("mdi_1_enr",vars)) |> spline_plot(y.n="mdi_1_enr",t.lab="One month",m.var=c(vars,"mdi_1_enr"))
+
+library(patchwork)
+p1/p2 + 
+  plot_layout(guides = "collect") +
+  plot_annotation(tag_levels = 'A',
+                  caption="Vertical dashed lines mark PASE score quartiles.")
+
+ggsave("spline_plot.png", plot = last_plot(), device = NULL, path = NULL,
+       scale = 1, width = 150, height = 200, dpi = 450, limitsize = TRUE,
+       units = "mm")
+```
+
+## Multivariate
+
+```{r}
+df <- dta[,c(vars,"mdi_6_newobs_enr")]
+model <- lm(mdi_6_newobs_enr~splines::ns(pase_0, knots = quants[2:4])+., data = df)
+# model <- coeftest(model, vcov = sandwich)
+
+df <- cbind(df,stats::predict(model, newdata=dta, interval="confidence"))
+
+
+g <- ggplot(df, aes(x=pase_0))
+g <- g + geom_point(aes(y = fit), size = 2, colour = "blue")
+g <- g + geom_smooth(data=df, aes(y=fit, ymin=lwr, ymax=upr), size = 1.5, 
+            colour = "red", se = TRUE, stat = "smooth")
+```
+
+
+# Kun indsamlede MDI
+
+```{r}
+## =============================================================================
+## Sensitivity regression analysis
+## =============================================================================
+
+source("data_format.R")
+
+outs<-c("mdi_1","mdi_6")
+
+lbl_x <- "_sens" # Extra label in file name
+
+inter_reg <- NULL # Interaction variables to include (only multivariate)
+
+biv_mul <- TRUE # Sets flag for both bivariate and multivariate or only multivariate analysis 
+
+robust_ci <- TRUE
+
+print_tbl <- TRUE # Print tables? T/F
+
+source("biv_mul_man.R")
+
+export
+
+# Notes
+# - Eliminate anything on stratification to slim the functions
+# - Move whole setup over from "00 master.R", easy block execution in .Rmd
+```
+
+# Interaktion
+
+```{r}
+## =============================================================================
+## Interaction regression analysis
+## =============================================================================
+
+source("data_format.R")
+
+vars<-c("active_treat", # vars() is defined in data_format.R, here "
+        "pase_0", 
+        "female",
+        "age",
+        "cohab",
+        "ever_smoker",
+        "diabetes", 
+        "hypertension", 
+        "afli", 
+        "ami", 
+        "tci",
+        "pad",
+        "nihss_0",
+        "any_reperf")
+
+outs<-c("mdi_1_enr","mdi_6_newobs_enr")
+
+lbl_x <- "_inter" # Extra label in file name
+
+inter_reg <- c("active_treat","pase_0") # Interaction variables to include (only multivariate)
+
+biv_mul <- TRUE # Sets flag for both bivariate and multivariate or only multivariate analysis 
+
+print_tbl <- TRUE # Print tables? T/F
+
+
+# source("regression_interaction.R")
+
+source("biv_mul_man.R")
+
+export
+```
+
+
+# NIHSS-MDI
+
+```{r}
+## =============================================================================
+## NIHSS~PASE0 regression analysis
+## =============================================================================
+
+source("data_format.R")
+
+vars<-c("pase_0_q", # New variables for analysis
+        "female",
+        "age",
+        "cohab",
+        "ever_smoker",
+        "diabetes", 
+        "hypertension", 
+        "afli", 
+        "ami", 
+        "tci",
+        "pad")
+
+outs<-c("nihss_0")
+
+lbl_x <- "_nihss-pase" # Extra label in file name
+
+inter_reg <- NULL # Interaction variables to include (only multivariate)
+
+biv_mul <- TRUE # Sets flag for both bivariate and multivariate or only multivariate analysis 
+
+print_tbl <- TRUE # Print tables? T/F
+
+# source("regression_nihss-pase.R")
+
+source("biv_mul_man.R")
+
+export
+```
+
+```{r message=FALSE, warning=FALSE}
+vars<-c("pase_0", # New variables for analysis
+        "female",
+        "age",
+        "cohab",
+        "ever_smoker",
+        "diabetes", 
+        "hypertension", 
+        "afli", 
+        "ami", 
+        "tci",
+        "pad")
+quants <- quantile(dta$pase_0,na.rm=TRUE)
+
+m<-lm(dta$nihss_0~splines::ns(pase_0, knots = quants[2:4])+., data = dta[,vars])
+
+coeftest(m, vcov = sandwich) |> tbl_regression() |> bold_p()
+```
+
+```{r}
+library(ggplot2)
+p1 <- ggplot(dta) + 
+  aes(pase_0,nihss_0) +
+  geom_point(na.rm = TRUE,color="#008080",shape=23) +
+  geom_smooth(method='lm', formula= y~x, na.rm = TRUE, se = FALSE)+
+  labs(title="Linear fit", x = "PASE score", y = "Acute NIHSS")+
+  theme_bw(14)
+
+p2 <- ggplot() +
+  aes(m$model$pase_0, predict(m)) +
+  geom_point(color="maroon",shape=25) +
+  geom_smooth(method='lm',formula = y~splines::ns(x, knots = quants[2:4])) +
+  geom_vline(xintercept=quants) +
+  # ylim(0,max(dta$mdi_1_enr)) +
+  labs(title="Qubic splines with predicted data", x = "PASE score", y = "Acute NIHSS")+
+  theme_bw(14)
+
+library(patchwork)
+p1+p2
+```

table_1.R

@@ -1,15 +1,18 @@
-tbl1_vars<-c("active_treat",vars,"inc_time")
+tbl1_vars<-c("active_treat",vars,"inc_time","mdi_1_enr","mdi_6_newobs_enr","mdi_1","mdi_6_newobs")
 
 tbl_1<-dta|>
-  tbl_summary(missing = "ifany",
+  tbl_summary(by = "pase_0_q",
+              missing = "ifany",
               include = all_of(tbl1_vars),
-              missing_text="(Missing)",
+              missing_text="(Missing)"#,
-              label = lab_sel(labels_all,tbl1_vars)
+              # label = lab_sel(labels_all,tbl1_vars)
   )|>
   add_n()|>
+  add_overall() |> 
   as_gt() |>
   # modify with gt functions
   gt::tab_header("Baseline Characteristics") |>
   gt::tab_options(
     table.font.size = "small",
     data_row.padding = gt::px(1))
+