---
title: "Quantitative Mapping of Human Hair Greying and Reversal in Relation to Life Stress"
output: 
  flexdashboard::flex_dashboard:
    orientation: columns
    vertical_layout: fill
    social: menu
    source_code: embed
  html_document:
    df_print: paged
runtime: shiny
---

```{r setup, include=FALSE}
library(flexdashboard)
library(mvtnorm)
library(tidyverse)
library(zoo)
library(plotly)
library(DT)
# Creates ? Help Popups
helpPopup <- function(title) {
  tagList(
    tags$a(
      href='#', title=title, 
      `data-toggle`='popover', `data-animation`=TRUE,
      `data-placement`='right', `data-trigger`='hover',
      
      tags$i(class='glyphicon glyphicon-question-sign')
    )
  )
} # END helpPopup

#function for generating the simulation data
sim_data_generate = function(beta1 = 1,beta2 = 1,sigma1 = 10,sigma2 = 11,
                             sigma3 = 12,sigmae = 1,rho12 = 0.5,
                             rho13 = 0.55,rho23 = 0.6,step = 0.1,
                             window_width = 5,stress_start = c(20,20),
                             stress_length = c(0.1,10), stress_strength = c(100,10),
                             thres = 1900,n=100,
                             seed_num = 2019){
  set.seed(seed_num)
  #Correlation matrix
  R <- matrix(numeric(3*3), nrow = 3)   # will be the correlation matrix
  diag(R) <- 1                        # set diagonal to 1
  R[upper.tri(R)] <- c(rho12, rho13, rho23)       # fill in to upper right
  R[lower.tri(R)] <- c(rho12, rho13, rho23)       # fill in to lower left
  #Covairance matrix
  G = diag(c(sigma1,sigma2,sigma3)) %*% R %*% diag(c(sigma1,sigma2,sigma3))
  #Stress data for one hair
  stress_data = tibble(age = seq(from = step, to = 100, by = step))
  for(i in 1:length(stress_strength)){#loop over different stress period
    stress_data_v = tibble(stress_raw = c(rep(0,stress_start[i]/step), 
                                        rep(stress_strength[i],stress_length[i]/step),#stress period
                                        rep(0,(100 - stress_start[i] - stress_length[i])/step)
                         )) %>% 
      #Cumulative sum of a window using window function rollapply
      mutate(stress_window = rollapply(stress_raw, width = window_width, 
                                       by = 1, FUN = sum, align = "right", fill = 0)) %>% 
      select(!!paste("stress_window",i,sep = "") := stress_window)
    stress_data = bind_cols(stress_data, stress_data_v)
  }

  stress_data = stress_data %>% 
    mutate(stress_window = rowSums(.[2:ncol(stress_data)])) %>% select(age, stress_window)
    
                           
 
  
  #Generate data for 100,000 * 20 = 2,000,000
  time_point = length(seq(from = step, to = 100, by = step)) # number of time points
  stress_data_all = do.call("rbind", replicate(n, stress_data, simplify = FALSE)) %>% 
    bind_cols(id = rep(1:n, each = time_point))#stress data for all hair
  sim_data = as_tibble(data.frame(id = 1:n, rmvnorm(n, rep(0,3), G)))
  names(sim_data) = c("id","b0", "b1", "b2")
  sim_data = left_join(sim_data, stress_data_all)
  dat = tibble(id = rep(1:n, each = time_point), erorr = rnorm(n * time_point,0,sigmae), 
                   age = rep(seq(from = step, to = 100, by = step), n),
                   beta1 = rep(beta1,n * time_point), beta2 = rep(beta2,n * time_point)) %>% 
    left_join(sim_data) %>% 
    mutate(white_no_stress = b0 + (abs(b1) + beta1) * age + erorr, 
           white = b0 + (abs(b1) + beta1) * age + (abs(b2) + beta2) * stress_window + erorr)
  # thres = dat %>% filter(age==70) %>% select(white_no_stress) %>% 
  #   pull() %>% quantile(probs = quantile_percent)
  dat_binay_white = dat %>% 
    mutate(white_binary_no_stress = white_no_stress > thres, white_binary_stress = white > thres)
  dat_summary = dat_binay_white %>% group_by(age) %>% 
    summarise(percentage_white = sum(white_binary_stress)/n)
  return(list(full_dat = dat,percentage_dat = dat_summary,
              thres = thres,stress_dat = stress_data))
}


#Preserve the data after the action button 'runsim'



dat_sim = eventReactive(input$runsim,{
  #get the stress input
  stress_start = as.numeric(str_split(input$stress_start,",")[[1]])
  stress_length = as.numeric(str_split(input$stress_length,",")[[1]])
  stress_strength = as.numeric(str_split(input$stress_strength,",")[[1]])
  #simulate the data
  sim_data_generate(beta1 = input$beta1,beta2 = input$beta2,
                    sigma1 = input$sigma1,sigma2 = input$sigma2,
                    sigma3 = input$sigma3, sigmae = input$sigmae,
                    rho12 = input$rho12,rho13 = input$rho13,
                    rho23 = input$rho23,step = 0.1, 
                    window_width = input$window_width,
                    stress_start = stress_start,
                    stress_strength = stress_strength,
                    stress_length = stress_length,
                    thres = input$thres,# quantile_percent = input$quantile_percent/100
                    n = input$n, seed_num = input$seed_num)
}, ignoreNULL=FALSE)



```
Simulation Plots
=======================================================================



Column {.sidebar}
-----------------------------------------------------------------------

```{r}

actionButton("runsim", "Run simulation",icon("paper-plane"), 
    style="color: #fff; background-color: #337ab7; border-color: #2e6da4")

sliderInput("beta1", 
            tagList("\\( \\beta_{1} \\): Fixed rate of increase in aging factor (slope)"), 
            min = 0, max = 20, 
            value = 16)

sliderInput("beta2", 
            tagList("\\( \\beta_{2} \\): Fixed rate of hair sensitivity to stress (slope)"), 
            min = 0, max = 20, 
            value = 0)

sliderInput("sigma1", 
            tagList("\\( \\sigma_{0} \\): Standard deviation of initial aging factor value at age=0 "), 
            min = 0, max = 50, 
            value = 10)

sliderInput("sigma2", 
            tagList("\\( \\sigma_{1} \\): Standard deviation of aging factor rate across the lifespan "), 
            min = 0, max = 50, 
            value = 13)

sliderInput("sigma3", 
            tagList("\\( \\sigma_{2} \\): Standard deviation of hair sensitivity to stress rate "), 
            min = 0, max = 50, 
            value = 0)

sliderInput("sigmae", 
            tagList("\\( \\sigma_{e} \\): Standard deviation of measurement error"), 
            min = 0, max = 50, 
            value = 10)

sliderInput("rho12", 
            tagList("\\( \\rho_{01} \\): Correlation between initial aging factor value at age 0 (\\( b_i0 \\)) and aging factor rate across the lifespan (\\( b_i1 \\))"), 
            min = 0, max = 1, 
            value = 0)

sliderInput("rho13", 
            tagList("\\( \\rho_{02} \\): Correlation between initial aging factor value at age 0 (\\( b_i0 \\)) and hair sensitivity to stress rate (\\( b_i2 \\))"), 
            min = 0, max = 1, 
            value = 0)

sliderInput("rho23", 
            tagList("\\( \\rho_{12} \\): Correlation between aging factor rate across the lifespan (\\( b_i1 \\)) and hair sensitivity to stress rate (\\( b_i2 \\)) "), 
            min = 0, max = 1, 
            value = 0)

sliderInput("window_width", 
            tagList("Window width: Pre-set number of years accumulating stress before age of evaluation"),
            min = 0, max = 10, 
            value = 0)

textInput("stress_start", 
          tagList("Stress start time: The starting time for stress. For multiple stressors, please separate the different start time with ‘,’"), 
          "0")

textInput("stress_strength", 
          tagList("Stress strenght: please separate the different stress strength with ‘,’"), 
          "0")

textInput("stress_length", 
          tagList("Stress duration: please separate the different stress duration with ‘,’"), 
          "0")
# sliderInput("stress_length1", "Stress length1", min = 0, max = 20, 
#             value = 0.1, step = 0.1)
# sliderInput("stress_strength1", "Stress strength1", min = 0, max = 500, 
#             value = 100)
# 
# sliderInput("stress_length2", "Stress length2", min = 0, max = 10, 
#             value = 10, step = 0.1)
# sliderInput("stress_strength2", "Stress strength2", min = 0, max = 500, 
#             value = 50)




numericInput('thres', tagList("Greying threshold"), 
             value = 1920)
# sliderInput("quantile_percent", 
#             tagList("Threshold", helpPopup("The whiteloading threshold for hair turning white")), 
#             min = 1, max = 99, 
#             value = 75)
sliderInput("n", 
            tagList("N: Number of hairs in simulation"), 
            min = 100, max = 1000, 
            value = 1000)
numericInput("seed_num", 
             tagList("Simulation seed number"),
             value = 2019)
sliderInput("per_hair_show",
            tagList("Percent of hairs to show in the white loading plot"),
            min = 1, max = 100, value=10)



```


Column {data-width=650}
-----------------------------------------------------------------------

### Hair greying trajectories

```{r 1}
renderPlotly({
  
  ggplotly(
    dat_sim()$percentage_dat %>% 
      ggplot(aes(x = age, y = percentage_white)) + geom_col(color = "light blue") +
      theme_bw()  + scale_y_continuous(labels = scales::percent, limits = c(0,1)) + 
      labs(y = "Proportion of grey hairs", x = "Age (years)")
  )

})
```

Column {data-width=350}
-----------------------------------------------------------------------

### Accumulation of aging factor in individual hairs

```{r 2}
renderPlot({
  
  
  dat_sim()$full_dat %>% 
    #sample certain amount of hair to display
    filter(id %in% sample(1:input$n, round(input$n*(input$per_hair_show/100),0)) ) %>% 
    ggplot(aes(x = age, y = white, group = id)) + 
    geom_line(color = "light blue", alpha = 0.5)  +
    geom_hline(yintercept=dat_sim()[[3]], 
               color = "red") + theme_bw() + labs(y = "Aging Factor (A.U.)", x = "Age (years)")


})
```

### Window of stress

```{r 3}
renderPlot({
  
  dat_sim()$stress_dat %>% 
    ggplot(aes(x = age, y = stress_window)) + geom_line(color = "light blue") +
    theme_bw() + labs(y = "Stress Window (A.U.)", x = "Age (years)")

})


```

Download data
=======================================================================



Column {data-width=200}
-----------------------------------------------------------------------
### Percentage of white hair 
```{r}
downloadButton("downloadPer", "Download")
downloadHandler(
    filename = function() {
      paste("percentage_data.csv")
    },
    content = function(file) {
      write.csv(dat_sim()$percentage_dat, file, row.names = F)
    }
  )
renderDataTable({
  datatable(dat_sim()$percentage_dat,rownames = F)
})

```


Column {data-width=200}
-----------------------------------------------------------------------
### Stress data
```{r}
downloadButton("downloadStress", "Download")
downloadHandler(
    filename = function() {
      paste("stress_data.csv")
    },
    content = function(file) {
      write.csv(dat_sim()$stress_dat, file, row.names = F)
    }
  )
renderDataTable({
  datatable(dat_sim()$stress_dat,rownames = F)
})
```

Column {data-width=400}
----------------------------------------------------------------------
### Full simulation data
```{r}
downloadButton("downloadFull", "Download")
downloadHandler(
    filename = function() {
      paste("full_data.csv")
    },
    content = function(file) {
      write.csv(dat_sim()$full_dat, file, row.names = F)
    }
  )
renderDataTable({
  datatable(
    (dat_sim()$full_dat %>% 
    select(-white_no_stress) %>% 
    rename(stress = stress_window)),
    rownames = F
  ) 
})

```




Details
=======================================================================
To simulate the greying process for a fixed person, we built a linear mixed model for the $i$th ($i=1,\ldots,\ n$) hair with two fixed effects ($\beta_1$ aging factor rate and $\beta_2$ hair sensitivity to stress rate) and three random effects ( $b_{i0}$ for initial aging factor value at age 0, $b_{i1}$ for aging factor rate across the lifespan and $b_{i2}$ for hair sensitivity to stress rate). We are taking the absolute value of these random effects because the slopes should be positive for age and accumulating stress. 
$$AgingFactorValue_{\left\{i,\ age\right\}}=\left|b_{i0}\right|+\left(\left|b_{i1}\right|+\beta_1\right)age+\left(\left|b_{i2}\right|+\beta_2\right)AccumulatingStress_{\left\{age\right\}}+e_i$$

Where AccumulatingStress is defined as:
$$
AccumulatingStress_{age}=\sum_{a=age-WindowWidth}^{age}{stress_a}
$$

The three random effects follow a multivariate normal:
$$\left(b_{i0},b_{i1},b_{i2}\right)\sim N\left(0,G\right)$$
with mean zero and a covariance structure:
$$
G=\left[\begin{matrix}\sigma_0^2&\rho_{01}\sigma_0\sigma_1&\rho_{02}\sigma_0\sigma_2\\\rho_{01}\sigma_1\sigma_0&\sigma_1^2&\rho_{12}\sigma_1\sigma_2\\\rho_{02}\sigma_2\sigma_0&\rho_{12}\sigma_2\sigma_1&\sigma_2^2\\\end{matrix}\right]
$$
All the correlation $\rho_{01},\rho_{02},\rho_{12}$ in the simulation are set to be positive.  
When the AgingFactorValue for ith hair reach a predefined threshold, the ith hair will turn grey. 
  
  
**Description for the parameters that could be changed (i is the fix person we are simulating)**  
1. $\beta_1$:  Fixed rate of increase in aging factor (slope)  
2. $\beta_2$:  Fixed rate of hair sensitivity to stress (slope)  
3. $\sigma_0$:  Standard deviation of initial aging factor value at age 0 for $b_{i0}$  
4. $\sigma_1$:  Standard deviation of aging factor rate across the lifespan for $b_{i1}$  
5. $\sigma_2$:  Standard deviation of hair sensitivity to stress rate for $b_{i2}$  
6. $\sigma_e$:  Standard deviation of measurement error for $e_i$  
7. $\rho_{01}$: Correlation between initial aging factor value at age 0 ($b_{i0}$) and aging factor rate across the lifespan ($b_{i1}$)  
8. $\rho_{02}$: Correlation between initial aging factor value at age 0 ($b_{i0}$) and hair sensitivity to stress rate ($b_{i2}$)  
9. $\rho_{12}$: Correlation between aging factor rate across the lifespan ($b_{i1}$) and hair sensitivity to stress rate ($b_{i2}$)  
10. Window width:  Pre-set number of years accumulating stress before age of evaluation.  
11. Stress start time: The starting time for stress. For multiple stressors, please separate the different start time with ','    
12. Stress strength: please separate the different stress strength with ','  
13. Stress duration: please separate the different stress duration with ','  
14. Threshold: The whiteloading threshold for hair turning white  
15. N: Number of hairs in simulation  
16. Seed number: Simulation seed number


About
=======================================================================
This app was created by *Junting Ren and Todd Ogden*.  
For more information please refer to the manuscript: Quantitative Mapping of Human Hair Greying and Reversal in Relation to Life Stress  
Ayelet Rosenberg, Shannon Rausser, Junting Ren, Eugene Mosharov, R Todd Ogden, Clay Lacefield, Gabriel Sturm, Ralf Paus, Martin Picard