Widefield_5x_Ipsilateral_Pdgfrb_IntDen.qmd

---
title-block-banner: true
title: "Analysis of PDGFRβ reactivity (Integrated density) in the ipsilateral hemisphere"
subtitle: "Data analysis notebook"
date: today
date-format: full
author: 
  - name: "Daniel Manrique-Castano"
    orcid: 0000-0002-1912-1764
    degrees:
      - PhD
    affiliation: 
      - name: Univerisity Laval 
        department: Psychiatry and Neuroscience
        group: Laboratory of neurovascular interactions 
note: "GitHub: https://daniel-manrique.github.io/"

keywords: 
  - PDGFRβ
  - Scar formation
  - Brain injury
  - MCAO
  - Bayesian modeling 
   
license: "CC BY"

format:
   pdf: 
    toc: true
    number-sections: true
    colorlinks: true
   html:
    code-fold: true
    embed-resources: true
    toc: true
    toc-depth: 2
    toc-location: left
    theme: spacelab

knitr:
  opts_chunk: 
    warning: false
    message: false

csl: science.csl
bibliography: references.bib
---

# Preview

This notebook reports the analysis for PDGFRβ expression using the FIJI-derived "Integrated density measurement".

**Parent dataset:** PDGFRβ and GFAP-stained ischemic hemispheres imaged at 5x (with stitching). Samples are grouped at 0 (Sham), 3, 7, 14, and 30 days post-ischemia (DPI). The raw images and pre-processing scripts (if applicable) are available at the Zenodo repository (10.5281/zenodo.10553084) under the name `Widefield_5x_Ipsilateral_Gfap-Pdgfrb.zip`.

**Working dataset**: The `Data_Processed/Widefield_5x_Ipsilateral_Gfap-Pdgfrb/Widefield_5x_Ipsilateral_Gfap-Pdgfrb_Inten.csv`data frame containing the integrated density of PDGFRβ measured in the ipsilateral hemisphere (`Gfap_IntDen`) using a FIJI [script](https://osf.io/mqukx) [@schindelin2012]. We perform scientific inference on this variable as a basis to understand PDGFRβ reactivity in our model of cerebral ischemia. The current knowledge suggest a close to linear increase in GFAP expression from the injury onset. We anticipate that PDGFRβ expression will increase in ischemic regions after the injury onset.


# Install and load required packages

Install and load all required packages. Please uncomment (delete #) the line code if installation is required. Load the installed libraries each time you start a new R session.

```{r}
#| label: Install_Packages
#| include: true
#| warning: false
#| message: false

#install.packages("devtools")
#library(devtools)

#install.packages(c("bayesplot", "bayestestR", "brms","dplyr", "easystats", "ggplot","gtsumary", "modelbased", "modelr","modelsummary", "patchwork", "poorman", "tidybayes", "tidyverse", "viridis"))


library(bayesplot)
library(bayestestR)
library(brms)
library(dplyr)
library(easystats)
library(ggplot2)
library(gtsummary)
library(modelbased)
library(modelr)
library(modelsummary)
library(patchwork)
library(poorman)
library(tidybayes)
library(tidyverse)
library(viridis)
```

# Visual themes

We create a visual theme to use in our plots (ggplots).

```{r}
#| label: Plot_Theme
#| include: true
#| warning: false
#| message: false
  
Plot_theme <- theme_classic() +
  theme(
      plot.title = element_text(size=18, hjust = 0.5, face="bold"),
      plot.subtitle = element_text(size = 10, color = "black"),
      plot.caption = element_text(size = 12, color = "black"),
      axis.line = element_line(colour = "black", size = 1.5, linetype = "solid"),
      axis.ticks.length=unit(7,"pt"),
     
      axis.title.x = element_text(colour = "black", size = 16),
      axis.text.x = element_text(colour = "black", size = 16, angle = 0, hjust = 0.5),
      axis.ticks.x = element_line(colour = "black", size = 1),
      
      axis.title.y = element_text(colour = "black", size = 16),
      axis.text.y = element_text(colour = "black", size = 16),
      axis.ticks.y = element_line(colour = "black", size = 1),
      
      legend.position="right",
      legend.direction="vertical",
      legend.title = element_text(colour="black", face="bold", size=12),
      legend.text = element_text(colour="black", size=10),
      
      plot.margin = margin(t = 10,  # Top margin
                             r = 2,  # Right margin
                             b = 10,  # Bottom margin
                             l = 10) # Left margin
      ) 
```

# Exploratory data visualization

We load the `ProData/Pdgfrb_Gfap_Ipsi_5x_Inten.csv` dataset to very its content.

```{r}
#| label: tbl-Pdgfrb_Table
#| include: true
#| warning: false
#| message: false
#| tbl-cap: "Data set"

# We load the dataset in case is not present in the R environment
Pdgfrb_Summary <- read.csv(file = "Data_Processed/Widefield_5x_Ipsilateral_Gfap-Pdgfrb/Widefield_5x_Ipsilateral_Gfap-Pdgfrb_Inten.csv", header = TRUE)

gt::gt(Pdgfrb_Summary[1:10,])

```

From this table, we focus on `DPI` (Days post-ischemia) and `Pdgfrb_IntDen` variables.`Pdgfrb_IntDen` is the product of the stained area and the mean intensity as calculated in FIJI using the script available at https://osf.io/mqukx. Additionally, we'll use the `Tissue_Area` variable to fit a model relating PDGFRβ reactivity and hemispheric shrinkage.

Next, we visualize the raw data to guide the statistical modeling. We plot the response variable as a density and a scatter plot (per DPI). In the scatter plot, we fit lines for a lineal (black), 2-degree (red), and 3-degree (green) polynomial models.

```{r}
#| label: fig-Pdgfrb_IntDen_Expl
#| include: true
#| warning: false
#| message: false
#| results: false
#| fig-cap: Exploratory data visualization for PDGFRβ expression
#| fig-width: 9
#| fig-height: 4


set.seed(8807)


# Density plot
################

Pdgfrb_IntDen_Dens <- 
  ggplot(
    data  = Pdgfrb_Summary, 
    aes(x = Pdgfrb_IntDen)
    ) +
  geom_density(size = 1.5) +
  geom_rug(size = 1) +
  scale_x_continuous(name ="PDGFRβ (IntDen)") +
  scale_y_continuous(name = "Density") +
  Plot_theme


# Scatter plot
##################

Pdgfrb_IntDen_Sctr <- 
  ggplot(
    data  = Pdgfrb_Summary, 
    aes(x = DPI, 
        y = Pdgfrb_IntDen)
    ) +
geom_smooth(
  method = "lm", 
  se     = TRUE,
  color  = "black") +
geom_smooth(
  method  = "lm", 
  se      = TRUE, 
  formula = y ~ poly(x, 2), 
  color   = "darkred") +
geom_smooth(
  method  = "lm", 
  se      = TRUE, 
  formula = y ~ poly(x, 3), 
  color   = "darkgreen") +
geom_jitter(
  width = 0.5, 
  shape = 1, 
  size  = 1.5, 
  color = "black") +
scale_y_continuous(name = "PDGFRβ (integrated density)") +
scale_x_continuous(
  name   ="Days post-ischemia (DPI) ",
  breaks =c(0, 3,7,14,30)
  ) +
Plot_theme

#Plot the result
Pdgfrb_IntDen_Dens |  Pdgfrb_IntDen_Sctr 
```
@fig-Pdgfrb_IntDen_Expl shows a range of \~ 0 - 1500. Without accounting for days-post ischemia (DPI), the response variable exhibits a single peak far from normal distributions. On the other hand, the scatter shows that there are no substantial differences between the fitted lines. Still, non-linear models capture between the deviation at 15 DPI.

# Statistical modeling for PDGFRβ reactivity

We use the `brms` package [@brms; @burkner2017] to run Bayesian modeling. We employ weakly informative priors in all the cases to optimize the Markov chain Monte Carlo (MCMC) and the exploration of the parameter space. We fit the following models:

-   **Pdgfrb_IntDen_Mdl1:** We use `DPI` as a linear predictor for PDGFRβ integrated density. The model takes the following notation:

$$
IntDen_{i} = \alpha + \beta_{1} DPI_{i} + \epsilon_{i}
$$

This model uses the following user-defined priors:

$$
\begin{align}
\beta_{1} \sim Normal(0, 1) \\
\alpha \sim Student-t(3, 500, 300) \\
\sigma \sim Student-t(3, 0, 100), \sigma > 0 \\
\end{align}
$$

-   **Pdgfrb_IntDen_Mdl2:** We use smooth terms (splines) to model the influence of `DPI`:

$$
IntDen_{i} = f_{s(DPI)}(DPI_{i})
$$

This model uses the same priors as model 1.

-   **Pdgfrbb_IntDen_Mdl3:** We fit a splines model with heteroscedasticity (sigma as additional response variable) keeping `DPI` the unique predictor:

$$
IntDen_{i} = f_{s(DPI)}(DPI_{i}) + \epsilon_{i} \\
sigma_{i} = f_{s(DPI)}(DPI_{i}) + \epsilon_{i}
$$

This model uses the following used-defined priors:
$$
\begin{align}
b \sim \mathcal{N}(100, 50) \\
\text{Intercept} \sim t_{3}(500, 300), \text{lb}=0 \\
\sigma \sim t_{3}(0, 100), \text{class} = b, \text{lb}=0 \\
\end{align}
$$

-   **Pdgfrb_IntDen_Mdl4** Finally, we investigate the correlation between brain shrinkage and PDGFRβ integrated density. We do so by fitting a spline model with interaction using the `t2` function, which denote a tensor product smooth of the predictor variables `DPI` and `Tissue_Area` to model the response variable (`Pdgfr_IntDen`).

$$ 
IntDen = \beta_0 + \sum_{i=1}^{k_1} \sum_{j=1}^{k_2} \beta_{ij} B_{1i}(DPI) B_{2j}(Tissue\_Area) + \epsilon
$$ 
Where: - $ \beta_0$ is the intercept. - $B_{1i}(DPI)$ and $B\_{2j}(Tissue_Area)$ are basis functions for the `DPI` and `Tissue_Area` predictors. - $ k_1$ and $ k_2 $ are the number of basis functions for `DPI` and `Tissue_Area`, respectively. - $\beta_{ij}$ are the coefficients for the tensor product of the $i^{th}$ basis function of `DPI` and the $j^{th}$ basis function of `Tissue_Area`. Finally, $\epsilon$ is the residual error.

This model uses the brms-default priors.

## Fit the models for PDGFRβ integrated density

```{r}
#| label: PdgfrbIntDen_Modeling_Ipsi5x
#| include: true
#| warning: false
#| message: false
#| results: false
#| cache: true

# Model 1: DPI as a linear predictor
###########################################
Pdgfrb_IntDen_Mdl1 <- bf(Pdgfrb_IntDen ~ DPI)

get_prior(Pdgfrb_IntDen_Mdl1, Pdgfrb_Summary)

Pdgfrb_IntDen_Prior1 <- 
  c(prior(normal(0,1), class = b),
    prior(normal(100,50), class = Intercept, lb= 0),
    prior(student_t(1, 0, 50), class = sigma, lb=0))

# Fit model 1
Pdgfrb_IntDen_Fit1 <- 
  brm(
    data    = Pdgfrb_Summary,
    formula = Pdgfrb_IntDen_Mdl1,
    prior   = Pdgfrb_IntDen_Prior1,
    chains  = 4,
    cores   = 4,
    warmup  = 2500, 
    iter    = 5000, 
    seed    = 8807,
    control = list(adapt_delta = 0.99, max_treedepth = 15),
    file    = "Models/Widefield_5x_Ipsilateral_Gfap-Pdgfrb_IntDen/Widefield_5x_Ipsilateral_Pdgfrb_IntDen_Fit1.rds",
    file_refit = "never") 
                     
# Add loo for model comparison
Pdgfrb_IntDen_Fit1 <- 
  add_criterion(Pdgfrb_IntDen_Fit1, c("loo", "waic", "bayes_R2"))


# Model 2: Model with splines / DPI
########################################################

Pdgfrb_IntDen_Mdl2 <- bf(Pdgfrb_IntDen ~ s(DPI, k = 5))

get_prior(Pdgfrb_IntDen_Mdl2, Pdgfrb_Summary)

# Fit model 2
Pdgfrb_IntDen_Fit2 <- 
  brm(
    data    = Pdgfrb_Summary,
    formula = Pdgfrb_IntDen_Mdl2,
    prior   = Pdgfrb_IntDen_Prior1,
    knots   = list(DPI = c(0, 3, 7, 14, 30)),
    chains  = 4,
    cores   = 4,
    warmup  = 2500, 
    iter    = 5000, 
    seed    = 8807,
    control = list(adapt_delta = 0.99, max_treedepth = 15),
    file    = "Models/Widefield_5x_Ipsilateral_Gfap-Pdgfrb_IntDen/Widefield_5x_Ipsilateral_Pdgfrb_IntDen_Fit2.rds",
    file_refit = "never") 
                     
# Add loo for model comparison
Pdgfrb_IntDen_Fit2 <- 
  add_criterion(Pdgfrb_IntDen_Fit2, c("loo", "waic", "bayes_R2"))


# Model 3: Model with splines / DPI and sigma
########################################################

Pdgfrb_IntDen_Mdl3 <- 
  bf(Pdgfrb_IntDen ~ s(DPI, k=5),
     sigma ~  s(DPI, k=5)
     )

get_prior(Pdgfrb_IntDen_Mdl3, Pdgfrb_Summary)

Pdgfrb_IntDen_Prior2 <- 
  c(prior(normal(0, 1), class = b),
    prior(student_t(1, 100, 50), class = Intercept, lb= 0),
    prior(student_t(1, 0, 50), class = b, dpar=sigma, lb=0))

# Fit model 3
Pdgfrb_IntDen_Fit3 <- 
  brm(
    data    = Pdgfrb_Summary,
    formula = Pdgfrb_IntDen_Mdl3,
    prior   = Pdgfrb_IntDen_Prior2,
    knots   = list(DPI = c(0, 3, 7, 14, 30)),
    chains  = 4,
    cores   = 4,
    warmup  = 5000, 
    iter    = 10000, 
    seed    = 8807,
    control = list(adapt_delta = 0.999, max_treedepth = 15),
    file    = "Models/Widefield_5x_Ipsilateral_Gfap-Pdgfrb_IntDen/Widefield_5x_Ipsilateral_Pdgfrb_IntDen_Fit3.rds",
    file_refit = "never") 
                     
# Add loo for model comparison
Pdgfrb_IntDen_Fit3 <- 
  add_criterion(Pdgfrb_IntDen_Fit3, c("loo", "waic", "bayes_R2"))


# Model 4: Spline interaction between DPI and hemispheric area
########################################################

Pdgfrb_IntDen_Mdl4 <- 
  bf(Pdgfrb_IntDen ~ t2(DPI,Tissue_Area))

get_prior(Pdgfrb_IntDen_Mdl4, Pdgfrb_Summary)

# Fit model 4
Pdgfrb_IntDen_Fit4 <- 
  brm(
    data    = Pdgfrb_Summary,
    formula = Pdgfrb_IntDen_Mdl4,
    chains  = 4,
    cores   = 4,
    warmup  = 2500, 
    iter    = 5000, 
    seed    = 8807,
    control = list(adapt_delta = 0.99, max_treedepth = 15),
    file    = "Models/Widefield_5x_Ipsilateral_Gfap-Pdgfrb_IntDen/Widefield_5x_Ipsilateral_Pdgfrb_IntDen_Fit4.rds",
    file_refit = "never") 
                     
# Add loo for model comparison
Pdgfrb_IntDen_Fit4 <- 
  add_criterion(Pdgfrb_IntDen_Fit4, c("loo", "waic", "bayes_R2"))
```

## Model comparison

We perform model comparison using the WAIC criteria for estimating the out-of-sample prediction error. By considering both the log-likelihood of the observed data and the effective number of parameters, it provides a balance between model fit and complexity. Unlike some other criteria, WAIC inherently accounts for the posterior distribution of the parameters rather than relying on point estimates, making it particularly suited for Bayesian analyses. In this case, we exclude model 4 (Mdl4) for the comparison because it fulfills another purpose.

Given a data set and a Bayesian model, the WAIC is computed as:

$$
\text{WAIC} = -2 \times \left( \text{LLPD} - p_{\text{WAIC}} \right)
$$

Where: $\text{LLPD}$ is the log pointwise predictive density, calculated as the average log likelihood for each observed data point across the posterior samples. $\text{WAIC}$ is the effective number of parameters, computed as the difference between the average of the log likelihoods and the log likelihood of the averages across posterior samples.

We compare the models using the `loo` function and plot it as a graph:

```{r}
#| label: fig-PdgfrbIntDen_Compare
#| include: true
#| warning: false
#| message: false
#| results: false
#| fig-cap: Model camparison by WAIC
#| fig-height: 4
#| fig-width: 5

Pdgfrb_IntDen_W <- 
loo_compare(
  Pdgfrb_IntDen_Fit1, 
  Pdgfrb_IntDen_Fit2, 
  Pdgfrb_IntDen_Fit3, 
  criterion = "waic")

# Generate WAIC graph
Pdgfrb_IntDen_WAIC <- 
  Pdgfrb_IntDen_W[, 7:8] %>% 
  data.frame() %>% 
  rownames_to_column(var = "model_name") %>% 
  
ggplot(
  aes(x    = model_name, 
      y    = waic, 
      ymin = waic - se_waic, 
      ymax = waic + se_waic)
  ) +
  geom_pointrange(shape = 21) +
  scale_x_discrete(
    breaks=c("Pdgfrb_IntDen_Fit1", 
             "Pdgfrb_IntDen_Fit2", 
             "Pdgfrb_IntDen_Fit3"), 
            
    labels=c("Mdl1", 
             "Mdl2", 
             "Mdl3") 
             
    ) +
  coord_flip() +
  labs(x = "", 
       y = "WAIC (score)",
       title = "") +
  Plot_theme

Pdgfrb_IntDen_WAIC
```

We see that the model three is the less penalized. Therefore, we continue our statistical inference with this model. 

## Model diagnostics

To evaluate sample predictions, we perform the diagnostics for model 3 (Mdl3) using the `pp_check` (posterior predictive checks) function from `brms`. In the graph, 𝘺 shows the data and y\~ the simulated data.

```{r}
#| label: fig-PdgfrbIntDen_Diagnostics
#| include: true
#| warning: false
#| message: false
#| results: false
#| fig-cap: Model diagnostics for PDGFRβ expression
#| fig-height: 4
#| fig-width: 5

set.seed(8807)

# Model 4
Pdgfrb_IntDen_Mdl3_pp <- 
  brms::pp_check(Pdgfrb_IntDen_Fit3, 
                 ndraws = 100) +
  labs(title = "Posterior predictive checks (model 3)",
  subtitle = "Formula: Pdgfrb_IntDen ~ s(DPI, K=5), sigma ~ s(DPI, K=5)") +
  #scale_y_continuous(limits=c(0, 0.15)) +
  Plot_theme 

Pdgfrb_IntDen_Mdl3_pp
```

@fig-PdgfrbIntDen_Diagnostics shows that the model predictions do not deviate substantially from the data. We can explore further the model using `shinystan`:

```{r}
#| label: PdgfrbIntDen_Shiny
#| include: true
#| warning: false
#| message: false
#| results: false

#launch_shinystan(Pdgfrb_IntDen_Fit3)
```

# Model results

### Visualizing the results

First, we plot the full posterior distributions:

```{r}
#| label: fig-PdgfrbIntDen_fit
#| include: true
#| warning: false
#| message: false
#| results: false
#| fig-cap: Posterior for PDGFRβ IntDen
#| fig-width: 5
#| fig-height: 4

set.seed(8807)

# We define the predictor grid
nd <- Pdgfrb_Summary %>% 
  distinct(DPI)

# compute and save the fitted draws
Pdgfrb_Post <- fitted(
  Pdgfrb_IntDen_Fit3, 
  newdata = nd,
  # this line returns the full set of HMC draws
  summary = FALSE) %>% 
  data.frame() %>% 
  set_names(nd %>% pull(DPI)) 
  
Pdgfrb_Post %>% 
  # convert to the long format
  pivot_longer(everything()) %>% 
  mutate(name = fct_relevel(name, "0", "3", "7", "14", "30")) %>% 
  
  ggplot(aes(x = value, y = factor(name))) +
  # adjust the point and interval summary settings as desired
  stat_halfeye(aes(fill = after_stat(level)), .width = c(0.50,0.80,0.90)) +
  scale_fill_brewer() +
  scale_y_discrete(name = "DPI") +
                   
  scale_x_continuous(name = "PDGFRβ (IntDen)",
                     limits = c(0, 1200)) +
  Plot_theme
```
We see a constant increase in PDGFRβ reactivity up yo 30 DPI. Subsequently, we visualize the results using the `conditional_effects` function from `brms` and add graph aesthetics using ggplot.

The `conditional_effects` function display the effects of predictors for brms-fitted models. It computes conditional effects of specified predictors while marginalizing over the remaining predictors. Thereby, it allows to visualization of single variable effects to make the data more interpretable. Mathematically, let $y = f(x_1, x_2, \ldots, x_n)$ denote the modeled relationship between the variable $ y $ and the predictor variables $ x_1, x_2, \ldots, x_n $. The `conditional_effects` function calculates the conditional expectations $E[y | x_i]$ for a given $ x_i $, marginalizing over other predictors. Formally, it can be represented in this manner:

$$
 E[y | x_i] = \int f(x_1, x_2, \ldots, x_n) p(x_1, x_2, \ldots, x_{i-1}, x_{i+1}, \ldots, x_n) dx_1 dx_2 \ldots dx_{i-1} dx_{i+1} \ldots dx_n
$$
We use this function to graph the effect of DPI in PDGFRβ IntDen and the corresponding sigma:

```{r}
#| label: fig-PdgfrbIntDen_CondEff
#| include: true
#| warning: false
#| message: false
#| results: false
#| fig-cap: Posterior for PDGFRβ expression
#| fig-width: 9
#| fig-height: 4

set.seed(8807)


# Conditional effects for Pdgfrb_IntDen
#####################################

# We convert the estimates to a data frame
Pdgfrb_IntDen_DPI <- 
  conditional_effects(Pdgfrb_IntDen_Fit3)

Pdgfrb_IntDen_DPIdf <- 
  as.data.frame(Pdgfrb_IntDen_DPI$DPI)

Pdgfrb_IntDen_pred <- 
  Pdgfrb_Summary %>%
  data_grid(DPI, n = 10) %>%
  add_predicted_draws(
    Pdgfrb_IntDen_Fit3, 
    ndraws = 50)

Pdgfrb_IntDen_DPI_Fig <-  
  ggplot(
    data  = Pdgfrb_IntDen_DPIdf,
    aes(x = DPI,
        y = estimate__)
    )+
  geom_jitter(
    data  = Pdgfrb_Summary, 
    aes(x = DPI, y=Pdgfrb_IntDen),
    width = 0.5, 
    color = "black", 
    shape = 18, 
    size  = 2.5, 
    alpha = 2/3) +
  geom_ribbon(
    aes(ymin  = lower__, 
        ymax  = upper__),
        alpha = 0.2)+
  geom_line(
    size     = 1,
    position = position_dodge(0.05), 
    color    = "blue") +
  scale_y_continuous(
    name = "PDGFRβ (IntDen)",
    limits = c(0, 1500),
    breaks = seq(0, 1500, 500)
    ) +
  scale_x_continuous(name="DPI",
                     breaks = c(0, 3, 7, 14, 30)) +
  Plot_theme

ggsave(
  plot     = Pdgfrb_IntDen_DPI_Fig, 
  filename = "Plots/Widefield_5x_Ipsilateral_Gfap-Pdgfrb_IntDen/Widefield_5x_Ipsilateral_Pdgfrb_IntDen_Fit3.png", 
  width    = 9, 
  height   = 9, 
  units    = "cm")



# Conditional smooths for sigma
###############################

Pdgfrb_IntDen_Smooths <- 
  conditional_smooths(Pdgfrb_IntDen_Fit3)

Pdgfrb_IntDen_Sigma_df <- 
  as.data.frame(Pdgfrb_IntDen_Smooths$`sigma: s(DPI,k=5)`)

Pdgfrb_IntDen_Sigma_fig <-  
  ggplot(Pdgfrb_IntDen_Sigma_df,
         aes(x = DPI,y = estimate__))+
  geom_ribbon(
    aes(ymin  =lower__, 
        ymax = upper__), 
    alpha = 0.2)+
  geom_line(size  = 1,
            color ="blue",
            position = position_dodge(0.05)) +
  scale_y_continuous(name= expression("PDGFRβ (log)"[sigma])) +
  scale_x_continuous(name="DPI") +
  Plot_theme
          
ggsave(plot     = Pdgfrb_IntDen_Sigma_fig, 
       filename = "Plots/Widefield_5x_Ipsilateral_Gfap-Pdgfrb_IntDen/Widefield_5x_Ipsilateral_Pdgfrb_IntDen_Sigma_Fit3.png", 
       width    = 9, 
       height   = 9, 
       units    = "cm")

Pdgfrb_IntDen_DPI_Fig | Pdgfrb_IntDen_Sigma_fig
```

@fig-PdgfrbIntDen_CondEff shows the fit for PDGFRβ integrated density and sigma. We can visualize an increasing trend during the first two weeks after injury. In contrast to GFAP reactivity, our model predicts a *plateau* phase for PDGFRβ cells. This implies that PDGFRβ integrated density remains constant despite brain shrinkage. We ponder that these dynamics are featured by an exchange between area and intensity, given that PDGFRβ appear to remain (and likely loss) in regions prone to shrinkage.

## Posterior summary

We plot the posterior summary using the `describe_posterior` function from `bayestestR` package [@bayestestR; @makowski2019].

```{r}
#| label: PdgfrbIntDen_DescribePosterior
#| include: true
#| warning: false
#| message: false
#| results: false
#| cache: true

describe_posterior(
  Pdgfrb_IntDen_Fit3,
  effects = "all",
  test = c("p_direction", "rope"),
  component = "all",
  centrality = "median")

modelsummary(Pdgfrb_IntDen_Fit3, 
             shape = term ~ model + statistic,
             centrali2ty = "mean", 
             title = "PDGFRβ integrated density following MCAO",
             statistic = "conf.int",
             gof_omit = 'ELPD|ELDP s.e|LOOIC|LOOIC s.e|WAIC|RMSE',
             output = "Tables/html/Widefield_5x_Ipsilateral_Pdgfrb_IntDen_Fit3_Table.html",
             )

Pdgfrb_IntDen_Fit3_Table <- modelsummary(Pdgfrb_IntDen_Fit3, 
             shape = term ~ model + statistic,
             centrality = "mean", 
             statistic = "conf.int",
             gof_omit = 'ELPD|ELDP s.e|LOOIC|LOOIC s.e|WAIC|RMSE',
             output = "gt")
gt::gtsave (Pdgfrb_IntDen_Fit3_Table, filename = "Tables/tex/Widefield_5x_Ipsilateral_Pdgfrb_IntDen_Fit3_Table.tex")
```

## Estimate derivative

We estimate the derivative to separate positive from negative trends in the fitted curve and obtain more clear linear fits. In statistical modeling, derivatives provide a measure of how a predictor variable influences the response variable, offering a more nuanced interpretation beyond the sign and magnitude of parameter estimates.

In mathematical terms, the derivative of a smooth function $f(x)$ is defined as:

$$ 
f'(x) = \lim\_{{\Delta x \to 0}} \frac{f(x + \Delta x) - f(x)}{\Delta x} 
$$ {#eq-derivative}

For a Bayesian model where $y = f(x)$, the derivative $\frac{dy}{dx}$ or $f'(x)$ represents the rate of change in the response variable $y$ for a unit change in the predictor variable $x$. If the model is non-linear, this derivative is a function of $x$ itself, allowing us to visualize how this rate of change varies across different values of $x$.

First, we estimate point-to point derivatives:

```{r}
#| label: PdgfrbIntDen_EstimateRelation
#| include: true
#| warning: false
#| message: false
#| results: false
#| cache: true

Pdgfrb_IntDen_est <- 
  estimate_relation(
    Pdgfrb_IntDen_Fit3, 
    at = "DPI", 
    ci = c(0.95))

Pdgfrb_IntDen_est
```
The output depict the distinct "fragments" from the curve with predicted point estimates and their uncertainty. This table is more informative that the raw posterior summary from splines. However, we summarize further using the `estimate_smooth`function.

```{r}
#| label: PdgfrIntDen_EstimateSmooth
#| include: true
#| warning: false
#| message: false
#| results: false
#| cache: true
 
Pdgfrb_IntDen_Derivative <- estimate_smooth(Pdgfrb_IntDen_est, x = "DPI")
Pdgfrb_IntDen_Derivative

gt::gt(round(Pdgfrb_IntDen_Derivative, digits = 2)) %>%
gt::gtsave(filename = "Tables/html/Widefield_5x_Ipsilateral_Pdgfrb_IntDen_Fit3_Deriv.html")

gt::gt(round(Pdgfrb_IntDen_Derivative, digits = 2)) %>%
gt::gtsave(filename = "Tables/tex/Widefield_5x_Ipsilateral_Pdgfrb_IntDen_Fit3_Deriv.tex")
```
The output denote a sustained increase (slope = 36) for PDGFRβ integrated density up to \~ 23 DPI (Change of 841, 95%CI = 831 - 1222). This is followed by a slightly decrease in slope (-4.4) at 30 DPI.

## Visualize derivatives

Finally, we plot the derivatives:

```{r}
#| label: fig-Pdgfrb_EstimateSlopes
#| include: true
#| warning: false
#| message: false
#| results: false
#| fig-cap: Derivatives for PDGFRβ expression
#| fig-height: 4
#| fig-width: 5

Pdgfrb_IntDen_deriv <- 
  estimate_slopes(
    Pdgfrb_IntDen_Fit3,
    trend = "DPI",
     at = "DPI",
    length = 30
    )

# Visualize
Pdgfrb_IntDen_deriv_Fig <- 
  plot(Pdgfrb_IntDen_deriv) + 
  labs(x = "DPI", 
       y = "Effect of DPI",
       title = NULL
       ) +
  Plot_theme +
  theme(legend.position= "none")

ggsave(
  plot     = Pdgfrb_IntDen_deriv_Fig, 
  filename = "Plots/Widefield_5x_Ipsilateral_Gfap-Pdgfrb_IntDen/Widefield_5x_Ipsilateral_Pdgfrb_IntDen_Deriv.png", 
  width    = 9, 
  height   = 9, 
  units    = "cm")

Pdgfrb_IntDen_deriv_Fig 
```

@fig-Pdgfrb_EstimateSlopes shows that in the chronic phase (after two weeks post-ischemia), there is negligible change in PDGFRβ integrated density.

# PDGFRβ relation to hemispheric area

Previously, we fitted a model (Pdgfrb_IntDen \~ t2(DPI,Tissue_Area)) to investigate the relation between PDGFRβ integrated density and hemispheric area,as we did with GFAP.

## Model diagnostics

We run model diagnostics with `pp_check` to verify the predictions. In the graph, 𝘺 shows the data and y\~ the simulated data.

```{r}
#| label: fig-GfapIntDen_Hem_Diagnostics
#| include: true
#| warning: false
#| message: false
#| results: false
#| fig-cap: Model diagnostics for PDGFRβ expression / hemispheric area
#| fig-height: 4
#| fig-width: 5

set.seed(8807)

# Model 3
Pdgfrb_IntDen_Mdl4_pp <- 
  brms::pp_check(Pdgfrb_IntDen_Fit4, 
                 ndraws = 100) +
  labs(title = "Posterior predictive checks (model 4)",
  subtitle = "Formula: Pdgfrb_IntDen ~ t2(DPI,Tissue_Area)") +
  #scale_y_continuous(limits=c(0, 0.15)) +
  Plot_theme  
  
Pdgfrb_IntDen_Mdl4_pp
```
We can appreciate that the fitted model does not imply major deviations from the observed data. 


## Visualizaing the results

Now, we plot the results using the `conditional_smooths` function from `brms`.
```{r}
#| label: fig-PdgfrbIntDen_Hemis
#| include: true
#| warning: false
#| message: false
#| results: false
#| fig-cap: PDGFRβ relation to hemispheric area
#| fig-width: 5
#| fig-height: 4

set.seed(8807)

# We convert the estimates to a data frames
Pdgfrb_IntDen_Area <- conditional_smooths(Pdgfrb_IntDen_Fit4)

Pdgfrb_IntDen_Area_Fig  <- 
  plot(Pdgfrb_IntDen_Area, 
       plot = FALSE)[[1]] +
  
  scale_y_continuous(name = expression ("Hemispheric area (mm)"^2)) +
  scale_x_continuous(name="DPI") +
    guides(colour = guide_colourbar(direction = "horizontal")) +
  guides(colour = guide_colourbar(barheight = unit(2, "cm"))) +
  scale_color_gradientn(
    name = "IntDen", 
    colors = viridis(6)) +
  Plot_theme +
  theme(
    legend.position = c(0.17, 0.8))
  
ggsave(
  plot     = Pdgfrb_IntDen_Area_Fig , 
  filename = "Plots/Widefield_5x_Ipsilateral_Gfap-Pdgfrb_IntDen/Widefield_5x_Ipsilateral_Pdgfrb_IntDen-Area.png", 
  width    = 9, 
  height   = 9, 
  units    = "cm")

Pdgfrb_IntDen_Area_Fig
```

@fig-PdgfrbIntDen_Hemis estimates that PDGFRβ reactivity increases with DPI and brain shrinkage. This might imply that brain shrinkage is associated with PDGFRβ expression, an aspect that can be further explored with causal inference in our stroke model.

## Posterior summary

We plot the posterior summary as previously done:

```{r}
#| label: Pdgfrb_IntDenArea_DescribePosterior
#| include: true
#| warning: false
#| message: false
#| results: false
#| cache: true

describe_posterior(
  Pdgfrb_IntDen_Fit4,
  effects = "all",
  test = c("p_direction", "rope"),
  component = "all",
  centrality = "median")

modelsummary(Pdgfrb_IntDen_Fit4, 
             shape = term ~ model + statistic,
             centrali2ty = "mean", 
             title = "PDGFRβ integrated density following MCAO",
             statistic = "conf.int",
             gof_omit = 'ELPD|ELDP s.e|LOOIC|LOOIC s.e|WAIC|RMSE',
             output = "Tables/html/Widefield_5x_Ipsilateral_Pdgfrb_IntDen_Fit4_Table.html",
             )

Pdgfrb_IntDen_Fit4_Table <- modelsummary(Pdgfrb_IntDen_Fit4, 
             shape = term ~ model + statistic,
             centrality = "mean", 
             statistic = "conf.int",
             gof_omit = 'ELPD|ELDP s.e|LOOIC|LOOIC s.e|WAIC|RMSE',
             output = "gt")
gt::gtsave (Pdgfrb_IntDen_Fit4_Table, filename = "Tables/tex/Widefield_5x_Ipsilateral_Pdgfrb_IntDen_Fit4_Table.tex")
```

# References

::: {#refs}
:::


```{r}
sessionInfo()
```