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