Demographics

Row

Sex

Age

League

School

Sport

Sport Level

Row

RTL Summary

RTL Sex

RTL Age

RTL League

RTL School

RTL Sport

RTL Total

Sex

Age

League

RTL Group Summary

Row

RTP Summary

RTP Sex

RTP Age

RTP League

RTP School

RTP Sport

RTL Total

Sex

Age

League

Test One PCSS Summary Scores

Row

Total Symptom Score

Total Symptom Score Summary

Sex

Sex Summary

Age

Age Summary

Row

Headache-Migraine

Headache-Migraine Summary

Sex

Sex Summary

Age

Age Summary

Headache-Migraine Normalized

Headache-Migraine Summary Normalized

Sex Normalized

Sex Summary Normalized

Age Normalized

Age Summary Normalized

Row

Cognitive

Cognitive Summary

Sex

Sex Summary

Age

Age Summary

Cognitive Normalized

Cognitive Summary Normalized

Sex Normalized

Sex Summary Normalized

Age Normalized

Age Summary Normalized

Row

Anxiety-Mood

Anxiety-Mood Summary

Sex

Sex Summary

Age

Age Summary

Anxiety-Mood Normalized

Anxiety-Mood Summary Normalized

Sex Normalized

Sex Summary Normalized

Age Normalized

Age Summary Normalized

Row

Ocular-Motor

Ocular-Motor Summary

Sex

Sex Summary

Age

Age Summary

Ocular-Motor Normalized

Ocular-Motor Summary Normalized

Sex Normalized

Sex Summary Normalized

Age Normalized

Age Summary Normalized

Row

Vestibular

Vestibular Summary

Sex

Sex Summary

Age

Age Summary

Vestibular Normalized

Vestibular Summary Normalized

Sex Normalized

Sex Summary Normalized

Age Normalized

Age Summary Normalized

Row

Sleep

Sleep Summary

Sex

Sex Summary

Age

Age Summary

Sleep Normalized

Sleep Summary Normalized

Sex Normalized

Sex Summary Normalized

Age Normalized

Age Summary Normalized

Models

Row

Plot 1

Plot 2

Plot 3

Row

Test 1 Model


Call:
lm(formula = dys_btwn_onset_rtp_3 ~ total_symptom_score_post_injury_1, 
    data = simsimp)

Residuals:
     Min       1Q   Median       3Q      Max 
-11.2686  -5.0552  -0.9726   4.1376  17.1376 

Coefficients:
                                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)                       11.86243    0.28879  41.076   <2e-16 ***
total_symptom_score_post_injury_1  0.05508    0.01399   3.937    9e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.398 on 753 degrees of freedom
Multiple R-squared:  0.02017,   Adjusted R-squared:  0.01887 
F-statistic:  15.5 on 1 and 753 DF,  p-value: 8.998e-05

Sex Model


Call:
lm(formula = dys_btwn_onset_rtp_3 ~ gender, data = simsimp)

Residuals:
     Min       1Q   Median       3Q      Max 
-11.7553  -4.7553  -0.7553   4.2447  16.2447 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  12.7553     0.3847  33.152   <2e-16 ***
genderMale   -0.3515     0.4861  -0.723     0.47    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.461 on 753 degrees of freedom
Multiple R-squared:  0.000694,  Adjusted R-squared:  -0.0006331 
F-statistic: 0.5229 on 1 and 753 DF,  p-value: 0.4698

Age Model


Call:
lm(formula = dys_btwn_onset_rtp_3 ~ age, data = simsimp)

Residuals:
     Min       1Q   Median       3Q      Max 
-11.4951  -4.8490  -0.9225   4.1572  17.1510 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  11.0000     3.7192   2.958   0.0032 **
age14         1.9225     3.7583   0.512   0.6091   
age15         2.4951     3.7462   0.666   0.5056   
age16         0.8490     3.7482   0.226   0.8209   
age17         1.0307     3.7533   0.275   0.7837   
age18         0.8367     3.8314   0.218   0.8272   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.442 on 749 degrees of freedom
Multiple R-squared:  0.01189,   Adjusted R-squared:  0.005294 
F-statistic: 1.803 on 5 and 749 DF,  p-value: 0.11

Row

Additive Model: Sex and Test 1 Severity


Call:
lm(formula = dys_btwn_onset_rtp_3 ~ gender + total_symptom_score_post_injury_1, 
    data = simsimp)

Residuals:
    Min      1Q  Median      3Q     Max 
-11.223  -5.059  -0.950   4.091  17.091 

Coefficients:
                                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)                       11.90918    0.43950  27.097  < 2e-16 ***
genderMale                        -0.06879    0.48716  -0.141 0.887751    
total_symptom_score_post_injury_1  0.05478    0.01416   3.869 0.000119 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.402 on 752 degrees of freedom
Multiple R-squared:  0.0202,    Adjusted R-squared:  0.01759 
F-statistic: 7.752 on 2 and 752 DF,  p-value: 0.0004652

Additive Model: Age and Test 1 Severity


Call:
lm(formula = dys_btwn_onset_rtp_3 ~ age + total_symptom_score_post_injury_1, 
    data = simsimp)

Residuals:
     Min       1Q   Median       3Q      Max 
-11.3507  -4.6577  -0.8519   3.8259  16.2139 

Coefficients:
                                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)                       10.55657    3.68569   2.864   0.0043 ** 
age14                              1.68560    3.72315   0.453   0.6509    
age15                              2.22957    3.71130   0.601   0.5482    
age16                              0.75889    3.71269   0.204   0.8381    
age17                              0.80998    3.71813   0.218   0.8276    
age18                              0.12860    3.79934   0.034   0.9730    
total_symptom_score_post_injury_1  0.05543    0.01412   3.925 9.46e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.381 on 748 degrees of freedom
Multiple R-squared:  0.03183,   Adjusted R-squared:  0.02407 
F-statistic: 4.099 on 6 and 748 DF,  p-value: 0.0004671

Additive Model: Sex, Age and Test 1 Severity


Call:
lm(formula = dys_btwn_onset_rtp_3 ~ gender + age + total_symptom_score_post_injury_1, 
    data = simsimp)

Residuals:
     Min       1Q   Median       3Q      Max 
-11.3305  -4.6722  -0.8339   3.7926  16.1788 

Coefficients:
                                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)                       10.60862    3.72333   2.849 0.004503 ** 
genderMale                        -0.05006    0.49133  -0.102 0.918870    
age14                              1.66833    3.72947   0.447 0.654763    
age15                              2.21253    3.71752   0.595 0.551914    
age16                              0.73817    3.72071   0.198 0.842790    
age17                              0.79100    3.72524   0.212 0.831903    
age18                              0.12258    3.80231   0.032 0.974291    
total_symptom_score_post_injury_1  0.05518    0.01434   3.848 0.000129 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.385 on 747 degrees of freedom
Multiple R-squared:  0.03185,   Adjusted R-squared:  0.02277 
F-statistic:  3.51 on 7 and 747 DF,  p-value: 0.001027

Row

Interaction Model: Sex, Age and Test 1 Severity

Adjusted R-squared of 0.058


Call:
lm(formula = dys_btwn_onset_rtp_3 ~ gender * age * total_symptom_score_post_injury_1, 
    data = simsimp)

Residuals:
    Min      1Q  Median      3Q     Max 
-12.794  -4.461  -0.815   4.027  17.577 

Coefficients: (3 not defined because of singularities)
                                                    Estimate Std. Error t value
(Intercept)                                          1.66457    5.36191   0.310
genderMale                                           9.54937    4.00873   2.382
age14                                                8.66500    5.48506   1.580
age15                                                9.75859    5.44467   1.792
age16                                               10.38849    5.42846   1.914
age17                                               11.62258    5.45836   2.129
age18                                                2.55073    3.78956   0.673
total_symptom_score_post_injury_1                    0.11075    0.20952   0.529
genderMale:age14                                    -8.64881    4.26434  -2.028
genderMale:age15                                    -7.19646    4.17635  -1.723
genderMale:age16                                   -10.61545    4.15720  -2.554
genderMale:age17                                   -12.45559    4.20771  -2.960
genderMale:age18                                          NA         NA      NA
genderMale:total_symptom_score_post_injury_1        -0.13749    0.21292  -0.646
age14:total_symptom_score_post_injury_1              0.04036    0.21446   0.188
age15:total_symptom_score_post_injury_1             -0.03247    0.21211  -0.153
age16:total_symptom_score_post_injury_1             -0.10483    0.21342  -0.491
age17:total_symptom_score_post_injury_1             -0.04594    0.21384  -0.215
age18:total_symptom_score_post_injury_1                   NA         NA      NA
genderMale:age14:total_symptom_score_post_injury_1   0.16459    0.22591   0.729
genderMale:age15:total_symptom_score_post_injury_1   0.06747    0.21994   0.307
genderMale:age16:total_symptom_score_post_injury_1   0.21456    0.22180   0.967
genderMale:age17:total_symptom_score_post_injury_1   0.08820    0.22184   0.398
genderMale:age18:total_symptom_score_post_injury_1        NA         NA      NA
                                                   Pr(>|t|)   
(Intercept)                                         0.75631   
genderMale                                          0.01747 * 
age14                                               0.11460   
age15                                               0.07349 . 
age16                                               0.05605 . 
age17                                               0.03356 * 
age18                                               0.50110   
total_symptom_score_post_injury_1                   0.59725   
genderMale:age14                                    0.04290 * 
genderMale:age15                                    0.08528 . 
genderMale:age16                                    0.01087 * 
genderMale:age17                                    0.00317 **
genderMale:age18                                         NA   
genderMale:total_symptom_score_post_injury_1        0.51864   
age14:total_symptom_score_post_injury_1             0.85078   
age15:total_symptom_score_post_injury_1             0.87837   
age16:total_symptom_score_post_injury_1             0.62344   
age17:total_symptom_score_post_injury_1             0.82997   
age18:total_symptom_score_post_injury_1                  NA   
genderMale:age14:total_symptom_score_post_injury_1  0.46650   
genderMale:age15:total_symptom_score_post_injury_1  0.75911   
genderMale:age16:total_symptom_score_post_injury_1  0.33368   
genderMale:age17:total_symptom_score_post_injury_1  0.69105   
genderMale:age18:total_symptom_score_post_injury_1       NA   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.268 on 734 degrees of freedom
Multiple R-squared:  0.0833,    Adjusted R-squared:  0.05833 
F-statistic: 3.335 on 20 and 734 DF,  p-value: 1.342e-06

Interaction Model: Sex:Age plus Test 1 Severity

Adjusted R-squared of 0.047


Call:
lm(formula = dys_btwn_onset_rtp_3 ~ gender * age + total_symptom_score_post_injury_1, 
    data = simsimp)

Residuals:
     Min       1Q   Median       3Q      Max 
-11.7190  -4.5015  -0.8785   4.2035  17.1137 

Coefficients: (1 not defined because of singularities)
                                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)                         3.65748    4.31919   0.847  0.39738    
genderMale                          6.91721    2.32780   2.972  0.00306 ** 
age14                               8.30329    4.41232   1.882  0.06025 .  
age15                               8.22884    4.38704   1.876  0.06109 .  
age16                               7.83392    4.37830   1.789  0.07398 .  
age17                               9.81038    4.39417   2.233  0.02587 *  
age18                               1.42804    3.77885   0.378  0.70561    
total_symptom_score_post_injury_1   0.05316    0.01422   3.738  0.00020 ***
genderMale:age14                   -6.42936    2.57943  -2.493  0.01290 *  
genderMale:age15                   -5.45667    2.50853  -2.175  0.02993 *  
genderMale:age16                   -7.18383    2.50626  -2.866  0.00427 ** 
genderMale:age17                  -10.36710    2.54144  -4.079 5.01e-05 ***
genderMale:age18                         NA         NA      NA       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.304 on 743 degrees of freedom
Multiple R-squared:  0.06135,   Adjusted R-squared:  0.04746 
F-statistic: 4.415 on 11 and 743 DF,  p-value: 1.882e-06

Test 1 Severity ANOVA

Row

Test 1 Boxplot

Row

Two-Way ANOVA: Sex and Test 1 Severity

PWC

---
title: "HCAMP Paper Version 2"
output: 
  flexdashboard::flex_dashboard:
    orientation: rows
    social: menu
    source_code: embed
    vertical_layout: scroll
    theme: united
---

```{r setup, include=FALSE}
library(flexdashboard)
library(tidyverse)
library(here)
library(janitor)
library(rio)
library(colorblindr)
library(gghighlight)
library(forcats)
library(ggrepel)
library(knitr)
library(kableExtra)
library(reactable)
library(plotly)
library(glue)
library(fs)
library(rstatix)
library(ggpubr)
library(writexl) 
library(remotes)
library(profvis) 


# theme_fivethirtyeight <- function(base_size = 15, base_family = "") {
#   theme_grey(base_size = base_size, base_family = base_family) %+replace%
#     theme(
# 
#       # Base elements which are not used directly but inherited by others
#       line =              element_line(colour = '#DADADA', size = 0.75,
#                                        linetype = 1, lineend = "butt"),
#       rect =              element_rect(fill = "#F0F0F0", colour = "#F0F0F0",
#                                        size = 0.5, linetype = 1),
#       text =              element_text(family = base_family, face = "plain",
#                                        colour = "#656565", size = base_size,
#                                        hjust = 0.5, vjust = 0.5, angle = 0,
#                                        lineheight = 0.9),
# 
#       # Modified inheritance structure of text element
#       plot.title =        element_text(size = rel(1.5), family = '' ,
#                                        face = 'bold', hjust = -0.05,
#                                        vjust = 1.5, colour = '#3B3B3B'),
#       axis.title.x =      element_text(),
#       axis.title.y =      element_text(),
#       axis.text =         element_text(),
# 
#       # Modified inheritance structure of line element
#       axis.ticks =        element_line(),
#       panel.grid.major =  element_line(),
#       panel.grid.minor =  element_blank(),
# 
#       # Modified inheritance structure of rect element
#       plot.background =   element_rect(),
#       panel.background =  element_rect(),
#       legend.key =        element_rect(colour = '#DADADA'),
# 
#       # Modifiying legend.position
#       legend.position = 'none',
# 
#       complete = TRUE
#     )
# }
# 
# 
# theme_set(theme_fivethirtyeight())


theme_set(theme_minimal(15) +
            theme(legend.position = "bottom",
                  panel.grid.major.x = element_line(color = "gray60"),
                  panel.grid.minor.x = element_blank(),
                  panel.grid.major.y = element_blank())
          )
```

```{r global, include=FALSE}
#all clean sims data
sims_concussion_data <- read_csv(here("data", "sims_concussion_data.csv"))

sims_concussion_data <- sims_concussion_data %>% 
  mutate(age = as.factor(age))

simsimp <- read_csv(here("data", "clean_impact_sims_data.csv"))

mod_data <- read_csv(here("data", "clean_impact_sims_data.csv"))

mod_datav2 <- read_csv(here("data", "clean_impact_sims_data.csv"))

str(simsimp)

simsimp <- simsimp %>% 
  mutate(dataset = as.factor(dataset),
         school_year = as.factor(school_year),
         school = as.factor(school),
         league = as.factor(league),
         gender = as.factor(gender),
         age = as.factor(age),
         sport = as.factor(sport),
         injury = as.factor(injury)) %>% 
  mutate_if(is.numeric, round, digits = 3)

mod_data <- mod_data %>% 
  mutate(dataset = as.factor(dataset),
         school_year = as.factor(school_year),
         school = as.factor(school),
         league = as.factor(league),
         gender = as.factor(gender),
         sport = as.factor(sport),
         injury = as.factor(injury)) %>% 
  mutate_if(is.numeric, round, digits = 3)

str(mod_data)

mod_datav2 <- mod_datav2 %>% 
  mutate(dataset = as.factor(dataset),
         school_year = as.factor(school_year),
         school = as.factor(school),
         league = as.factor(league),
         gender = as.factor(gender),
         sport = as.factor(sport),
         injury = as.factor(injury)) %>% 
  mutate_if(is.numeric, round, digits = 3)

str(mod_datav2)
```

```{r, include=FALSE}
# creating new data set for machine learning models

names(simsimp)
str(simsimp)
names(mod_data)

mod_data[mod_data$age <= 16, "age_group"] <- "13-16"
mod_data[mod_data$age >=17, "age_group"] <- "17-18"

names(mod_data)

mod_data %>% 
  count(age_group)

str(mod_data)

mod_data <- mod_data %>% 
  mutate(age_group = as.factor(age_group))

str(mod_data)

mod_data[mod_data$dys_btwn_onset_rtp_3 <= 6, "rtl_group"] <- "0-6"
mod_data[mod_data$dys_btwn_onset_rtp_3 >= 7 & mod_data$dys_btwn_onset_rtp_3 <= 12,
         "rtl_group"] <- "7-12"
mod_data[mod_data$dys_btwn_onset_rtp_3 >= 13 & mod_data$dys_btwn_onset_rtp_3 <= 18,
         "rtl_group"] <- "13-18"
mod_data[mod_data$dys_btwn_onset_rtp_3 >= 19 & mod_data$dys_btwn_onset_rtp_3 <= 24,
         "rtl_group"] <- "19-24"
mod_data[mod_data$dys_btwn_onset_rtp_3 >= 25 & mod_data$dys_btwn_onset_rtp_3 <= 30,
         "rtl_group"] <- "25-30"

str(mod_data)

mod_data <- mod_data %>% 
  mutate(rtl_group = as.factor(rtl_group))

mod_data %>% 
  count(rtl_group)

mod_data_2 <- mod_data

mod_data_2[mod_data_2$total_symptom_score_post_injury_1 <= 16,
           "test_1_pcss_group"] <- "0-16"
mod_data_2[mod_data_2$total_symptom_score_post_injury_1 >= 17 &
           mod_data_2$total_symptom_score_post_injury_1 <= 32,
           "test_1_pcss_group"] <- "17-32"
mod_data_2[mod_data_2$total_symptom_score_post_injury_1 >= 33 &
           mod_data_2$total_symptom_score_post_injury_1 <= 48,
           "test_1_pcss_group"] <- "33-48"
mod_data_2[mod_data_2$total_symptom_score_post_injury_1 >= 49 &
           mod_data_2$total_symptom_score_post_injury_1 <= 64,
           "test_1_pcss_group"] <- "49-64"
mod_data_2[mod_data_2$total_symptom_score_post_injury_1 >= 65 &
           mod_data_2$total_symptom_score_post_injury_1 <= 80,
           "test_1_pcss_group"] <- "65-80"
mod_data_2[mod_data_2$total_symptom_score_post_injury_1 >= 81,
           "test_1_pcss_group"] <- "81 or higher"

str(mod_data_2)

mod_data_2 <- mod_data_2 %>% 
  mutate(test_1_pcss_group = as.factor(test_1_pcss_group))

str(mod_data_2)

mod_data_2 %>% 
  count(test_1_pcss_group)


# per Tom's suggestion, divding RTL group into two levels - splitting around 12, which is the mean


mod_datav2[mod_datav2$age <= 16, "age_group"] <- "13-16"
mod_datav2[mod_datav2$age >=17, "age_group"] <- "17-18"

mod_datav2[mod_datav2$dys_btwn_onset_rtp_3 <= 12, "rtl_group"] <- "0-12"
mod_datav2[mod_datav2$dys_btwn_onset_rtp_3 >= 13, "rtl_group"] <- "13-30"

mod_datav2[mod_datav2$total_symptom_score_post_injury_1 <= 16,
           "test_1_pcss_group"] <- "0-16"
mod_datav2[mod_datav2$total_symptom_score_post_injury_1 >= 17 &
           mod_datav2$total_symptom_score_post_injury_1 <= 32,
           "test_1_pcss_group"] <- "17-32"
mod_datav2[mod_datav2$total_symptom_score_post_injury_1 >= 33 &
           mod_datav2$total_symptom_score_post_injury_1 <= 48,
           "test_1_pcss_group"] <- "33-48"
mod_datav2[mod_datav2$total_symptom_score_post_injury_1 >= 49 &
           mod_datav2$total_symptom_score_post_injury_1 <= 64,
           "test_1_pcss_group"] <- "49-64"
mod_datav2[mod_datav2$total_symptom_score_post_injury_1 >= 65 &
           mod_datav2$total_symptom_score_post_injury_1 <= 80,
           "test_1_pcss_group"] <- "65-80"
mod_datav2[mod_datav2$total_symptom_score_post_injury_1 >= 81,
           "test_1_pcss_group"] <- "81 or higher"


str(mod_datav2)

mod_datav2 <- mod_datav2 %>% 
  mutate(rtl_group = as.factor(rtl_group),
         test_1_pcss_group = as.factor(test_1_pcss_group),
         age_group = as.factor(age_group))

str(mod_datav2)

mod_datav2 %>% 
  count(rtl_group)
```

```{r, include=FALSE}
#helpful functions 

mean_2 <- function(x) {
  z <- na.omit(x)
  sum(z) / length(z)
}

my_mean <- function(x) {
  mean(x[x >= 0], na.rm = TRUE)
}

create_react_time <- function(df, var) {
    df %>% 
      summarize(Mean = mean({{var}}),
                SD = sd({{var}}),
                Min = min({{var}}),
                Max = max({{var}}),
                Total = length({{var}})) %>% 
      mutate_if(is.numeric, round, 2) %>% 
      reactable(columns = list(
        Mean = colDef(format = colFormat(separators = TRUE, suffix = " days")),
        SD = colDef(format = colFormat(separators = TRUE, suffix = " days")),
        Min = colDef(format = colFormat(separators = TRUE, suffix = " days")),
        Max = colDef(format = colFormat(separators = TRUE, suffix = " days")),
        Total = colDef(format = colFormat(separators = TRUE, suffix = " concussions"))
      ))
}

create_react_time2 <- function(df, x, var) {
    df %>% 
    group_by({{x}}) %>% 
      summarize(Mean = mean({{var}}),
                SD = sd({{var}}),
                Min = min({{var}}),
                Max = max({{var}}),
                Total = length({{var}})) %>% 
      mutate_if(is.numeric, round, 2) %>% 
      reactable(columns = list(
        Mean = colDef(format = colFormat(separators = TRUE, suffix = " days")),
        SD = colDef(format = colFormat(separators = TRUE, suffix = " days")),
        Min = colDef(format = colFormat(separators = TRUE, suffix = " days")),
        Max = colDef(format = colFormat(separators = TRUE, suffix = " days")),
        Total = colDef(format = colFormat(separators = TRUE, suffix = " concussions"))
      ))
}

create_react <- function(df, var) {
    df %>% 
      summarize(Mean = mean({{var}}),
                SD = sd({{var}}),
                Min = min({{var}}),
                Max = max({{var}}),
                Total = length({{var}})) %>% 
      mutate_if(is.numeric, round, 3) %>% 
      reactable(columns = list(
        Mean = colDef(format = colFormat(separators = TRUE)),
        SD = colDef(format = colFormat(separators = TRUE)),
        Min = colDef(format = colFormat(separators = TRUE)),
        Max = colDef(format = colFormat(separators = TRUE)),
        Total = colDef(format = colFormat(separators = TRUE, suffix = " concussions"))
      ))
}


create_react_age <- function(df, var) {
    df %>% 
    group_by(age) %>% 
      summarize(Mean = mean({{var}}),
                SD = sd({{var}}),
                Min = min({{var}}),
                Max = max({{var}}),
                Total = length({{var}})) %>% 
      mutate_if(is.numeric, round, 3) %>% 
      reactable(columns = list(
        Mean = colDef(format = colFormat(separators = TRUE)),
        SD = colDef(format = colFormat(separators = TRUE)),
        Min = colDef(format = colFormat(separators = TRUE)),
        Max = colDef(format = colFormat(separators = TRUE)),
        Total = colDef(format = colFormat(separators = TRUE, suffix = " concussions"))
      ))
}

create_react_gender <- function(df, var) {
    df %>% 
    group_by(gender) %>% 
      summarize(Mean = mean({{var}}),
                SD = sd({{var}}),
                Min = min({{var}}),
                Max = max({{var}}),
                Total = length({{var}})) %>% 
      mutate_if(is.numeric, round, 3) %>% 
      reactable(columns = list(
        Mean = colDef(format = colFormat(separators = TRUE)),
        SD = colDef(format = colFormat(separators = TRUE)),
        Min = colDef(format = colFormat(separators = TRUE)),
        Max = colDef(format = colFormat(separators = TRUE)),
        Total = colDef(format = colFormat(separators = TRUE, suffix = " concussions"))
      ))
}


my_mean(simsimp$dys_btwn_onset_test_4)
```

```{r, include=FALSE}
simsimp %>% 
  count(student_id)

length(unique(simsimp$student_id))
length(unique(simsimp$gender))



simsimp %>% 
  group_by(row, gender) %>% 
  count()
```

# Demographics

Sidebar {.sidebar}
------------

The **Sex** table displays the total number of injuries by sex used in the data set. The total number of injuries is 755 that can be utilized for analysis. Like the previous iteration of the paper, some individuals sustained multiple injuries that are tracked individually. This is a characteristic that one of the reviewers specified we describe more to better explain the sample. The tables displayed present data representing the total number of _injuries_, which include instances of repeat injuries. Data on the number of unique individuals is outlined here: 

  *  **Number of females:** 271
  * **Number of males:** 460
  
  * 260 females sustained one tracked injury
  * 447 males sustained one tracked injury 
  * 10 females sustained two tracked injuries
  * 12 males sustained two tracked injuries
  * 1 female sustained three tracked injuries
  * 1 male sustained three tracked injuries 

Row {.tabset}
-----------------------------------------------------------------------

### Sex 

```{r, include=TRUE}
simsimp %>% 
  group_by(gender) %>% 
  summarize(total = n()) %>% 
  arrange(desc(total)) %>% 
  reactable(
    columns = list(
      gender = colDef(name = "Sex",
                      align = "center"),
      total = colDef(name = "Total",
                     align = "center",
                     format = colFormat(suffix = " injuries"))),
    pagination = TRUE,
    striped = TRUE,
    outlined = TRUE,
    compact = TRUE,
    highlight = TRUE,
    bordered = TRUE
    )
```

```{r, include=FALSE}
sims_sex <- simsimp %>% 
  group_by(gender) %>% 
  summarize(total = n()) %>% 
  arrange(desc(total))

sims_sex_plot <- ggplot(sims_sex, aes(fct_reorder(gender, total), total)) +
  geom_col(fill = "blue",
           alpha = 0.7) +
  scale_y_continuous(limits = c(0, 600),
                     breaks = c(0, 200, 400, 600)) +
  coord_flip() +
  labs(x = "",
       y = "Total")
```

```{r, include=FALSE}
ggplotly(sims_sex_plot)
```

### Age

```{r, include=TRUE}
simsimp %>% 
  group_by(age) %>% 
  summarize(total = n()) %>% 
  reactable(
    columns = list(
      age = colDef(name = "Age",
                      align = "center"),
      total = colDef(name = "Total",
                     align = "center",
                     format = colFormat(suffix = " injuries"))),
    pagination = TRUE,
    striped = TRUE,
    outlined = TRUE,
    compact = TRUE,
    highlight = TRUE,
    bordered = TRUE
    )
```

```{r, include=FALSE}
sims_age <- simsimp %>% 
  mutate(age = as.factor(age)) %>% 
  group_by(age) %>% 
  summarize(total = n()) %>% 
  arrange(desc(total))



sims_age_plot <- ggplot(sims_age, aes(fct_reorder(age, total), total)) +
  geom_col(fill = "blue",
           alpha = 0.7) +
  coord_flip() +
  labs(x = "Age",
       y = "Total")
```

```{r, include=FALSE}
ggplotly(sims_age_plot)
```

### League

```{r, include=TRUE}
simsimp %>% 
  group_by(league) %>% 
  summarize(total = n()) %>% 
  arrange(desc(total)) %>% 
  reactable(
    columns = list(
      league = colDef(name = "League",
                      align = "center"),
      total = colDef(name = "Total",
                     align = "center",
                     format = colFormat(suffix = " injuries"))),
    pagination = TRUE,
    striped = TRUE,
    outlined = TRUE,
    compact = TRUE,
    highlight = TRUE,
    bordered = TRUE
    )
```

### School

```{r}
simsimp %>% 
  group_by(school) %>% 
  summarize(total = n()) %>% 
  arrange(desc(total)) %>% 
  reactable(
    columns = list(
      school = colDef(name = "School",
                      align = "center"),
      total = colDef(name = "Total",
                     align = "center",
                     format = colFormat(suffix = " injuries"))),
    pagination = TRUE,
    striped = TRUE,
    outlined = TRUE,
    compact = TRUE,
    highlight = TRUE,
    bordered = TRUE,
    searchable = TRUE
    )
```

### Sport

```{r}
simsimp %>% 
  group_by(sport) %>% 
  summarize(total = n()) %>% 
  arrange(desc(total)) %>% 
  reactable(
    columns = list(
      sport = colDef(name = "Sport",
                      align = "center"),
      total = colDef(name = "Total",
                     align = "center",
                     format = colFormat(suffix = " injuries"))),
    pagination = TRUE,
    striped = TRUE,
    outlined = TRUE,
    compact = TRUE,
    highlight = TRUE,
    bordered = TRUE,
    searchable = TRUE
    )
```

### Sport Level

```{r, include=FALSE}
simsimp %>% 
  group_by(level) %>% 
  summarize(total = n()) %>% 
  arrange(desc(total)) %>% 
  reactable(
    columns = list(
      level = colDef(name = "Level",
                      align = "center"),
      total = colDef(name = "Total",
                     align = "center",
                     format = colFormat(suffix = " injuries"))),
    pagination = TRUE,
    striped = TRUE,
    outlined = TRUE,
    compact = TRUE,
    highlight = TRUE,
    bordered = TRUE
    )
```

Row {.tabset}
-----------------------------------------------------------------------

### RTL Summary 

```{r, include=TRUE}
create_react_time(simsimp, dys_btwn_onset_rtp_3)
```

### RTL Sex

```{r, include=TRUE}
create_react_time2(simsimp, gender, dys_btwn_onset_rtp_3)
```

### RTL Age

```{r, include=TRUE}
create_react_time2(simsimp, age, dys_btwn_onset_rtp_3)
```

### RTL League

```{r, include=TRUE}
create_react_time2(simsimp, league, dys_btwn_onset_rtp_3)
```

### RTL School

```{r, include=TRUE}
create_react_time2(simsimp, school, dys_btwn_onset_rtp_3)
```

### RTL Sport

```{r, include=TRUE}
create_react_time2(simsimp, sport, dys_btwn_onset_rtp_3)
```

```{r, include=FALSE}
rtl_smry_plot <- ggplot(simsimp, aes(dys_btwn_onset_rtp_3)) +
  geom_histogram(fill = "#56B4E9",
                color = "white", 
                alpha = 0.9,
                bins = 10) +
  labs(x = "Days to Complete RTL",
       y = "Number of Injuries")

rtp_smry_plot <- ggplot(simsimp, aes(dys_btwn_onset_rtp_7)) +
  geom_histogram(fill = "#56B4E9",
                color = "white", 
                alpha = 0.9,
                bins = 10) +
  labs(x = "Days to Complete RTP",
       y = "Number of Injuries")


rtl_smry_plot2 <- function(df, x, y) {
  p <- ggplot(df, aes({{x}})) +
  geom_histogram(fill = "#56B4E9",
                color = "white", 
                alpha = 0.9,
                bins = 10)
  p + facet_wrap(vars({{y}})) +
  labs(x = "Days to Complete RTL",
       y = "Number of Injuries")
}

rtp_smry_plot2 <- function(df, x, y) {
  p <- ggplot(df, aes({{x}})) +
  geom_histogram(fill = "#56B4E9",
                color = "white", 
                alpha = 0.9,
                bins = 10)
  p + facet_wrap(vars({{y}})) +
  labs(x = "Days to Complete RTP",
       y = "Number of Injuries")
}

rtl_smry_plot2(simsimp, dys_btwn_onset_rtp_3, gender)
```

### RTL Total 

```{r, include=TRUE}
ggplotly(rtl_smry_plot)
```

### Sex

```{r, include=TRUE}
ggplotly(rtl_smry_plot2(simsimp, dys_btwn_onset_rtp_3, gender))
```

### Age

```{r, include=TRUE}
ggplotly(rtl_smry_plot2(simsimp, dys_btwn_onset_rtp_3, age))
```

### League

```{r, include=TRUE}
ggplotly(rtl_smry_plot2(simsimp, dys_btwn_onset_rtp_3, league))
```

### RTL Group Summary

```{r, include=TRUE}
mod_data %>% 
  group_by(gender, age_group, rtl_group) %>% 
  summarize(total = n()) %>% 
  reactable(
    columns = list(
      gender = colDef(name = "Sex",
                      align = "center"),
      age_group = colDef(name = "Age Group",
                         align = "center"),
      rtl_group = colDef(name = "RTL Group",
                         align = "center",
                         format = colFormat(suffix = " days")),
      total = colDef(name = "Total",
                     align = "center",
                     format = colFormat(suffix = " injuries"))),
    pagination = TRUE,
    striped = TRUE,
    outlined = TRUE,
    compact = TRUE,
    highlight = TRUE,
    bordered = TRUE    
  )
```

Row {.tabset}
-----------------------------------------------------------------------

### RTP Summary 

```{r, include=TRUE}
create_react_time(simsimp, dys_btwn_onset_rtp_7)
```

### RTP Sex

```{r, include=TRUE}
create_react_time2(simsimp, gender, dys_btwn_onset_rtp_7)
```

### RTP Age

```{r, include=TRUE}
create_react_time2(simsimp, age, dys_btwn_onset_rtp_7)
```

### RTP League

```{r, include=TRUE}
create_react_time2(simsimp, league, dys_btwn_onset_rtp_7)
```

### RTP School

```{r, include=TRUE}
create_react_time2(simsimp, school, dys_btwn_onset_rtp_7)
```

### RTP Sport

```{r, include=TRUE}
create_react_time2(simsimp, sport, dys_btwn_onset_rtp_7)
```

### RTL Total 

```{r, include=TRUE}
ggplotly(rtp_smry_plot)
```

### Sex

```{r, include=TRUE}
ggplotly(rtp_smry_plot2(simsimp, dys_btwn_onset_rtp_7, gender))
```

### Age

```{r, include=TRUE}
ggplotly(rtp_smry_plot2(simsimp, dys_btwn_onset_rtp_7, age))
```

### League

```{r, include=TRUE}
ggplotly(rtp_smry_plot2(simsimp, dys_btwn_onset_rtp_7, league))
```


# Test One PCSS Summary Scores

Row {.tabset}
-----------------------------------------------------------------------

```{r, include=FALSE}
score_hist <- function(df, x) {
  ggplot(df, aes({{x}})) +
  geom_histogram(fill = "#56B4E9",
                color = "white", 
                alpha = 0.9,
                bins = 25) +
    labs(x = "Symptom Severity",
         y = "Number of Injuries")
}

gender_hist <- function(df, x) {
  ggplot(df, aes({{x}})) +
  geom_histogram(fill = "#56B4E9",
                color = "white", 
                alpha = 0.9,
                bins = 25) +
    facet_wrap(~gender) +
    labs(x = "Symptom Severity",
         y = "Number of Injuries")
}

age_hist <- function(df, x) {
  ggplot(df, aes({{x}})) +
  geom_histogram(fill = "#56B4E9",
                color = "white", 
                alpha = 0.9,
                bins = 25) +
    facet_wrap(~age) +
    labs(x = "Symptom Severity",
         y = "Number of Injuries")
}

names(simsimp)
```

### Total Symptom Score 

```{r, include=TRUE}
ggplotly(score_hist(simsimp, total_symptom_score_post_injury_1))
```

### Total Symptom Score Summary

```{r, include=TRUE}
create_react(simsimp, total_symptom_score_post_injury_1)
```

### Sex

```{r, include=TRUE}
ggplotly(gender_hist(simsimp, total_symptom_score_post_injury_1))
```

### Sex Summary

```{r, include=TRUE}
create_react_gender(simsimp, total_symptom_score_post_injury_1)
```

### Age

```{r, include=TRUE}
ggplotly(age_hist(simsimp, total_symptom_score_post_injury_1))
```

### Age Summary

```{r, include=TRUE}
create_react_age(simsimp, total_symptom_score_post_injury_1)
```


Row {.tabset}
-----------------------------------------------------------------------

### Headache-Migraine 

```{r, include=TRUE}
ggplotly(score_hist(simsimp, headache_migraine_cluster_score_post_injury_1))
```

### Headache-Migraine Summary

```{r, include=TRUE}
create_react(simsimp, headache_migraine_cluster_score_post_injury_1)
```

### Sex

```{r, include=TRUE}
ggplotly(gender_hist(simsimp, headache_migraine_cluster_score_post_injury_1))
```

### Sex Summary

```{r, include=TRUE}
create_react_gender(simsimp, headache_migraine_cluster_score_post_injury_1)
```

### Age

```{r, include=TRUE}
ggplotly(age_hist(simsimp, headache_migraine_cluster_score_post_injury_1))
```

### Age Summary

```{r, include=TRUE}
create_react_age(simsimp, headache_migraine_cluster_score_post_injury_1)
```

### Headache-Migraine Normalized

```{r, include=TRUE}
ggplotly(score_hist(simsimp, headache_migraine_test_1))
```

### Headache-Migraine Summary Normalized

```{r, include=TRUE}
create_react(simsimp, headache_migraine_test_1)
```

### Sex Normalized

```{r, include=TRUE}
ggplotly(gender_hist(simsimp, headache_migraine_test_1))
```

### Sex Summary Normalized

```{r, include=TRUE}
create_react_gender(simsimp, headache_migraine_test_1)
```

### Age Normalized 

```{r, include=TRUE}
ggplotly(age_hist(simsimp, headache_migraine_test_1))
```

### Age Summary Normalized

```{r, include=TRUE}
create_react_age(simsimp, headache_migraine_test_1)
```

Row {.tabset}
-----------------------------------------------------------------------

### Cognitive

```{r, include=TRUE}
ggplotly(score_hist(simsimp, cognitive_cluster_score_post_injury_1))
```

### Cognitive Summary

```{r, include=TRUE}
create_react(simsimp, cognitive_cluster_score_post_injury_1)
```

### Sex

```{r, include=TRUE}
ggplotly(gender_hist(simsimp, cognitive_cluster_score_post_injury_1))
```

### Sex Summary

```{r, include=TRUE}
create_react_gender(simsimp, cognitive_cluster_score_post_injury_1)
```

### Age

```{r, include=TRUE}
ggplotly(age_hist(simsimp, cognitive_cluster_score_post_injury_1))
```

### Age Summary

```{r, include=TRUE}
create_react_age(simsimp, cognitive_cluster_score_post_injury_1)
```

### Cognitive Normalized

```{r, include=TRUE}
ggplotly(score_hist(simsimp, cognitive_test_1))
```

### Cognitive Summary Normalized

```{r, include=TRUE}
create_react(simsimp, cognitive_test_1)
```

### Sex Normalized

```{r, include=TRUE}
ggplotly(gender_hist(simsimp, cognitive_test_1))
```

### Sex Summary Normalized

```{r, include=TRUE}
create_react_gender(simsimp, cognitive_test_1)
```

### Age Normalized 

```{r, include=TRUE}
ggplotly(age_hist(simsimp, cognitive_test_1))
```

### Age Summary Normalized

```{r, include=TRUE}
create_react_age(simsimp, cognitive_test_1)
```


Row {.tabset}
-----------------------------------------------------------------------

### Anxiety-Mood

```{r, include=TRUE}
ggplotly(score_hist(simsimp, anxiety_mood_cluster_score_post_injury_1))
```

### Anxiety-Mood Summary

```{r, include=TRUE}
create_react(simsimp, anxiety_mood_cluster_score_post_injury_1)
```

### Sex

```{r, include=TRUE}
ggplotly(gender_hist(simsimp, anxiety_mood_cluster_score_post_injury_1))
```

### Sex Summary

```{r, include=TRUE}
create_react_gender(simsimp, anxiety_mood_cluster_score_post_injury_1)
```

### Age

```{r, include=TRUE}
ggplotly(age_hist(simsimp, anxiety_mood_cluster_score_post_injury_1))
```

### Age Summary

```{r, include=TRUE}
create_react_age(simsimp, anxiety_mood_cluster_score_post_injury_1)
```

### Anxiety-Mood Normalized

```{r, include=TRUE}
ggplotly(score_hist(simsimp, anxiety_mood_test_1))
```

### Anxiety-Mood Summary Normalized

```{r, include=TRUE}
create_react(simsimp, anxiety_mood_test_1)
```

### Sex Normalized

```{r, include=TRUE}
ggplotly(gender_hist(simsimp, anxiety_mood_test_1))
```

### Sex Summary Normalized

```{r, include=TRUE}
create_react_gender(simsimp, anxiety_mood_test_1)
```

### Age Normalized 

```{r, include=TRUE}
ggplotly(age_hist(simsimp, anxiety_mood_test_1))
```

### Age Summary Normalized

```{r, include=TRUE}
create_react_age(simsimp, anxiety_mood_test_1)
```

Row {.tabset}
-----------------------------------------------------------------------

### Ocular-Motor 

```{r, include=TRUE}
ggplotly(score_hist(simsimp, ocular_motor_cluster_score_post_injury_1))
```

### Ocular-Motor Summary

```{r, include=TRUE}
create_react(simsimp, ocular_motor_cluster_score_post_injury_1)
```

### Sex

```{r, include=TRUE}
ggplotly(gender_hist(simsimp, ocular_motor_cluster_score_post_injury_1))
```

### Sex Summary

```{r, include=TRUE}
create_react_gender(simsimp, ocular_motor_cluster_score_post_injury_1)
```

### Age

```{r, include=TRUE}
ggplotly(age_hist(simsimp, ocular_motor_cluster_score_post_injury_1))
```

### Age Summary

```{r, include=TRUE}
create_react_age(simsimp, ocular_motor_cluster_score_post_injury_1)
```

### Ocular-Motor Normalized

```{r, include=TRUE}
ggplotly(score_hist(simsimp, ocular_motor_test_1))
```

### Ocular-Motor Summary Normalized

```{r, include=TRUE}
create_react(simsimp, ocular_motor_test_1)
```

### Sex Normalized

```{r, include=TRUE}
ggplotly(gender_hist(simsimp, ocular_motor_test_1))
```

### Sex Summary Normalized

```{r, include=TRUE}
create_react_gender(simsimp, ocular_motor_test_1)
```

### Age Normalized 

```{r, include=TRUE}
ggplotly(age_hist(simsimp, ocular_motor_test_1))
```

### Age Summary Normalized

```{r, include=TRUE}
create_react_age(simsimp, ocular_motor_test_1)
```


Row {.tabset}
-----------------------------------------------------------------------

### Vestibular 

```{r, include=TRUE}
ggplotly(score_hist(simsimp, vestibular_cluster_score_post_injury_1))
```

### Vestibular Summary

```{r, include=TRUE}
create_react(simsimp, vestibular_cluster_score_post_injury_1)
```

### Sex

```{r, include=TRUE}
ggplotly(gender_hist(simsimp, vestibular_cluster_score_post_injury_1))
```

### Sex Summary

```{r, include=TRUE}
create_react_gender(simsimp, vestibular_cluster_score_post_injury_1)
```

### Age

```{r, include=TRUE}
ggplotly(age_hist(simsimp, vestibular_cluster_score_post_injury_1))
```

### Age Summary

```{r, include=TRUE}
create_react_age(simsimp, vestibular_cluster_score_post_injury_1)
```

### Vestibular Normalized

```{r, include=TRUE}
ggplotly(score_hist(simsimp, vestibular_test_1))
```

### Vestibular Summary Normalized

```{r, include=TRUE}
create_react(simsimp, vestibular_test_1)
```

### Sex Normalized

```{r, include=TRUE}
ggplotly(gender_hist(simsimp, vestibular_test_1))
```

### Sex Summary Normalized

```{r, include=TRUE}
create_react_gender(simsimp, vestibular_test_1)
```

### Age Normalized 

```{r, include=TRUE}
ggplotly(age_hist(simsimp, vestibular_test_1))
```

### Age Summary Normalized

```{r, include=TRUE}
create_react_age(simsimp, vestibular_test_1)
```


Row {.tabset}
-----------------------------------------------------------------------

### Sleep

```{r, include=TRUE}
ggplotly(score_hist(simsimp, sleep_cluster_score_post_injury_1))
```

### Sleep Summary

```{r, include=TRUE}
create_react(simsimp, sleep_cluster_score_post_injury_1)
```

### Sex

```{r, include=TRUE}
ggplotly(gender_hist(simsimp, sleep_cluster_score_post_injury_1))
```

### Sex Summary

```{r, include=TRUE}
create_react_gender(simsimp, sleep_cluster_score_post_injury_1)
```

### Age

```{r, include=TRUE}
ggplotly(age_hist(simsimp, sleep_cluster_score_post_injury_1))
```

### Age Summary

```{r, include=TRUE}
create_react_age(simsimp, sleep_cluster_score_post_injury_1)
```

### Sleep Normalized

```{r, include=TRUE}
ggplotly(score_hist(simsimp, sleep_test_1))
```

### Sleep Summary Normalized

```{r, include=TRUE}
create_react(simsimp, sleep_test_1)
```

### Sex Normalized

```{r, include=TRUE}
ggplotly(gender_hist(simsimp, sleep_test_1))
```

### Sex Summary Normalized

```{r, include=TRUE}
create_react_gender(simsimp, sleep_test_1)
```

### Age Normalized 

```{r, include=TRUE}
ggplotly(age_hist(simsimp, sleep_test_1))
```

### Age Summary Normalized

```{r, include=TRUE}
create_react_age(simsimp, sleep_test_1)
```


# Models

Sidebar {.sidebar}
------------

The interaction models are the strongest relative models, approaching adjusted *R-*squared values of 0.05; however, the amount of variance accounted for within these models is still very small. Considering the plots at the top of the page, it is evident there is not a linear relationship between the hypothesized predictor variables (age, sex, test 1 PCSS score) and RTL duration time. 

Row {.tabset}
-----------------------------------------------------------------------

### Plot 1

```{r, include=FALSE}
p1 <- ggplot(simsimp, aes(dys_btwn_onset_rtp_3, total_symptom_score_post_injury_1)) +
  geom_point(color = "gray70") +
  geom_smooth() +
  geom_smooth(method = "lm",
              color = "magenta") +
  labs(x = "Days to Complete RTL",
       y = "Test 1 Total Symptom Severity Score")
```

```{r, include=TRUE}
ggplotly(p1)
```

### Plot 2

```{r, include=FALSE}
p2 <- ggplot(simsimp, aes(dys_btwn_onset_rtp_3, total_symptom_score_post_injury_1)) +
  geom_point(color = "gray70") +
  geom_smooth(aes(color = gender),
              method = "lm") +
  labs(x = "Days to Complete RTL",
       y = "Test 1 Total Symptom Severity Score")
```

```{r, include=TRUE}
ggplotly(p2)
```


### Plot 3

```{r, include=FALSE}
p3 <- ggplot(simsimp, aes(dys_btwn_onset_rtp_3, total_symptom_score_post_injury_1)) +
  geom_point(color = "gray70") +
  geom_smooth(aes(color = age),
              method = "lm") +
  labs(x = "Days to Complete RTL",
       y = "Test 1 Total Symptom Severity Score")
```

```{r, include=TRUE}
ggplotly(p3)
```


Row {.tabset}
-----------------------------------------------------------------------

### Test 1 Model

```{r, include=FALSE}
# modeling age:total symptom score model

names(simsimp)

sex_test_1_mod <- lm(dys_btwn_onset_rtp_3 ~ gender*total_symptom_score_post_injury_1, 
                     data = simsimp)


summary(sex_test_1_mod)
confint(sex_test_1_mod)

sex_age_mod <- lm(dys_btwn_onset_rtp_3 ~ gender*age, data = simsimp)

summary(sex_age_mod)

sex_age_test_1_mod <- lm(dys_btwn_onset_rtp_3 ~ 
                            gender*age*total_symptom_score_post_injury_1, 
                          data = simsimp)

summary(sex_age_test_1_mod)


# LM examples 

test_1_sev_mod <- lm(dys_btwn_onset_rtp_3 ~ total_symptom_score_post_injury_1,
                     data = simsimp)

summary(test_1_sev_mod)
fitted(test_1_sev_mod)


```

```{r, include=TRUE}
summary(test_1_sev_mod)
```

### Sex Model

```{r, include=FALSE}
sex_mod <- lm(dys_btwn_onset_rtp_3 ~ gender, data = simsimp)
```

```{r, include=TRUE}
summary(sex_mod)
```

### Age Model

```{r, include=FALSE}
age_mod <- lm(dys_btwn_onset_rtp_3 ~ age, data = simsimp)
```

```{r, include=TRUE}
summary(age_mod)
```


Row {.tabset}
-----------------------------------------------------------------------

### Additive Model: Sex and Test 1 Severity 

```{r, include=FALSE}
sex_test_add_mod <- lm(dys_btwn_onset_rtp_3 ~ gender + total_symptom_score_post_injury_1, 
                     data = simsimp)
```

```{r, include=TRUE}
summary(sex_test_add_mod)
```

### Additive Model: Age and Test 1 Severity 

```{r, include=FALSE}
age_test_add_mod <- lm(dys_btwn_onset_rtp_3 ~ age + total_symptom_score_post_injury_1, 
                     data = simsimp)
```

```{r, include=TRUE}
summary(age_test_add_mod)
```


### Additive Model: Sex, Age and Test 1 Severity 

```{r, include=FALSE}
sex_age_test_add_mod <- lm(dys_btwn_onset_rtp_3 ~ gender + age + 
                             total_symptom_score_post_injury_1, data = simsimp)
```

```{r, include=TRUE}
summary(sex_age_test_add_mod)
```


Row {.tabset}
-----------------------------------------------------------------------

### Interaction Model: Sex, Age and Test 1 Severity 

Adjusted *R-*squared of 0.058

```{r, include=FALSE}
sex_age_test_int_mod <- lm(dys_btwn_onset_rtp_3 ~ gender*age*total_symptom_score_post_injury_1, data = simsimp)
```

```{r, include=TRUE}
summary(sex_age_test_int_mod)
```

### Interaction Model: Sex:Age plus Test 1 Severity 

Adjusted *R-*squared of 0.047

```{r, include=FALSE}
sex_age_int_plus_test <- lm(dys_btwn_onset_rtp_3 ~ gender*age + total_symptom_score_post_injury_1, data = simsimp)
```

```{r, include=TRUE}
summary(sex_age_int_plus_test)
```

```{r, include=FALSE}
#plot(sex_age_int_plus_test)
```

# Test 1 Severity ANOVA

Sidebar {.sidebar}
------------

The ANOVA generated results similar to the findings from the first iteration of the paper. Females report higher symptom severity than males across most clusters. The headache-migraine, sleep, and cognitive clusters are rated with the highest symptom severity. 

```{r, include=FALSE}
simsimp_long <- simsimp %>% 
  pivot_longer(
    cols = c(60:65),
    names_to = "symptom_cluster",
    values_to = "score_test_1",
    names_pattern = "(.*)_test_1"
  )
```

```{r, include=FALSE}
cluster_bxp <- function(df, x, y) {
  ggplot(df, aes({{x}}, {{y}}, fill = gender)) +
  geom_boxplot() + 
  scale_fill_OkabeIto() +
  coord_flip() +
  labs(x = "Symptom Cluster",
       y = "Scaled Severity Score") +
  theme(legend.position = "bottom")
}

cluster_bxp(simsimp_long, symptom_cluster, score_test_1)
```

Row {.tabset}
-----------------------------------------------------------------------

### Test 1 Boxplot

```{r, include=TRUE, fig.width=10}
cluster_bxp(simsimp_long, symptom_cluster, score_test_1)
```

Row {.tabset}
-----------------------------------------------------------------------

### Two-Way ANOVA: Sex and Test 1 Severity

```{r, include=FALSE}
anova_react <- function(df) {
  df %>% 
    reactable(
      defaultColDef = colDef(align = "center",
                             format = colFormat(digits = 2, separators = TRUE)),
      pagination = TRUE,
      striped = TRUE,
      outlined = TRUE,
      compact = TRUE,
      highlight = TRUE,
      bordered = TRUE
    )
}

pwc_react <- function(df) {
  df %>% 
    reactable(
      defaultColDef = colDef(align = "center",
                             format = colFormat(digits = 3, separators = TRUE)),
      pagination = TRUE,
      striped = TRUE,
      outlined = TRUE,
      compact = TRUE,
      highlight = TRUE,
      bordered = TRUE,
      searchable = TRUE
    )  
}
```

```{r, include=FALSE}
aov_res <- aov(score_test_1 ~ symptom_cluster * gender,
               data = simsimp_long)

summary(aov_res)

aov_res

Anova(aov_res, type = "III")

TukeyHSD(aov_res)

aov_res2 <- Anova(aov_res, type = "III")

tukey_res <- TukeyHSD(aov_res)

```

```{r, include=TRUE}
anova_react(aov_res2)
```

### PWC

```{r, include=TRUE}
tidy(tukey_res) %>% 
  pwc_react()
```