Widefield_5x_Ipsilateral_Pdgfrb_LowHigh_Str.qmd

---
title-block-banner: true
title: "Analysis of PDGFR-β-low and PDGFR-β high+ cells in the ipsilateral hemisphere (striatal-only lesions)"
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-β
  - Brain injury
  - 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

In this notebook, we analyze the proportion and distribution of PDGFR-β_Low and PDGFR-β_High+ in brains with only striatal lesions. We performed a separate image pre-processing, cell detection and data analysis for these animals given that the experiments were carried out at a different moment. We performed automatic cell detection and classification using QuPath [@bankhead2017]. We also create point patterns using the coordinates of detecting cells to perform Point Pattern analysis.

**Parent dataset:** PDGFR-β and GFAP-stained ischemic hemispheres imaged at 5x (with stitching). Samples are grouped at 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`. Individual cells were detected and classified into PDGFR-β^low^ (Pdgfrb_NonReact) and PDGFR-β^high^ (Pdgfrb_React) using QuPath [@bankhead2017].The complete QuPath project, including classifiers and output data as .tsv files is available at https://osf.io/ty4z5.

**Working dataset**: The `Widefield_5x_Ipsilateral_Gfap-Pdgfrb_Inten_Str.csv`data frame containing the number of PDGFR-β^low^ (Pdgfrb_NonReact) and PDGFR-β^high^ (Pdgfrb_React) cells in the ischemic hemisphere. Here, we analyze the proportion and distribution of these populations.

# 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("devtools")
#library(devtools)

#install.packages(c("bayesplot", "bayestestR", "brms","broom.mixed", "dplyr", "easystats", "distributional", "ggplot","gtsummary", "modelbased", "modelr", "modelsummary", "patchwork", "poorman", "tidybayes", "tidyverse"))

library(bayesplot)
library(bayestestR)
library(brms)
library(broom.mixed)
library(dplyr)
library(easystats)
library(distributional)
library(ggplot2)
library(gtsummary)
library(modelbased)
library(modelr)
library(modelsummary)
library(patchwork)
library(poorman)
library(tidybayes)
library(tidyverse)
```

# 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

```{r}
#| label: Pdgfrb_Striatum_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_Striatum <- read.csv(file = "Data_Processed/Widefield_5x_Ipsilateral_Gfap-Pdgfrb/Widefield_5x_Ipsilateral_Gfap-Pdgfrb_Inten_Str.csv", header = TRUE)

gt::gt(Pdgfrb_Striatum[1:10,])
```

For this analysis, we focus on `DPI` (Days post-ischemia), `Pdgfrb_NonReact`, `Pdgfrb_React` variables to analyze these cells proportions in the ischemic brain. Next, we visualize the raw data to guide the statistical modeling. We plot the response variables as a scatter plot fitting a linear model, given that this data set comprises just two time points.

```{r}
#| label: fig-Pdgfrb_LowHigh_Str_Exploratory
#| include: true
#| warning: false
#| message: false
#| results: false
#| fig-cap: Exploratory data visualization for PDGFR-β^low^ and PDGFR-β^high^ (striatal lesions)
#| fig-width: 9
#| fig-height: 4


set.seed(8807)

# PDGFR-β^low^
###################

Striatum_Pdgfr_Low_Sctr <- 
  ggplot(
    data  = Pdgfrb_Striatum, 
    aes(x = DPI, 
        y = Pdgfrb_NonReact)) +
geom_smooth(
  method = "lm", 
  se     = TRUE,
  color  = "black") +
geom_jitter(
  width = 0.5, 
  shape = 1, 
  size  = 1.5, 
  color = "black") +
  scale_y_continuous(name= expression("Number of PDGFR-β"^low)) +
  scale_x_continuous(name="Days post-ischemia (DPI) ",
                     breaks=c(14,30)) +
  Plot_theme
  

# PDGFR-β^low^
#######################

Striatum_Pdgfr_High_Sctr <- 
  ggplot(
    data  = Pdgfrb_Striatum, 
    aes(x = DPI, 
        y = Pdgfrb_React)) +
geom_smooth(
  method = "lm", 
  se     = TRUE,
  color  = "black") +
geom_jitter(
  width = 0.5, 
  shape = 1, 
  size  = 1.5, 
  color = "black") +
  scale_y_continuous(name= expression("Number of PDGFR-β"^high)) +
  scale_x_continuous(name="Days post-ischemia (DPI) ",
                     breaks=c(14,30)) +
  Plot_theme
  
Striatum_Pdgfr_Low_Sctr | Striatum_Pdgfr_High_Sctr
```

@Pdgfrb_Striatum_LowHigh_Exploratory show that PDGFR-β^low^ (non-reactive cells) show a wide dispersion give the reduced number of observations at 14 DPI. Meanwhile, PDGFR-β^high^ display an increasing trend from 14 to 30 DPI. As done for brain with cortico-striatal lesions, we fit the model using a binomial family distribution. Please refer to the `Widefield_5x_Ipsilateral_Pdgfrb_LowHigh.qmd` notebook for additional details of this distribution model. 

# Statistical modeling 

## Fit the model

We fit a linear model using a binomial distribution `DPI` as a linear predictor for the probability of PDGFR-𝛽^high^:

$$
\log\left(\frac{p_{i}}{1 - p_{i}}\right) = \alpha + \beta_{1} \times DPI_{i}
$$ 
This model uses weakly informative priors based on the information gathered from the model for corticostriatal lesions. The prior take the following notation:  
 
$$
\begin{align}
b \sim \mathcal{N}(1, 0.2),  \text{lb}=0 \\
\text{Intercept} \sim student_t(3, 0, 2), \text{lb}=0 \\
\end{align}
$$

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

Pdgfrb_Prop_Str_Mdl1 <- bf(Pdgfrb_React | trials(Pdgfrb_Total) ~ DPI)

get_prior(Pdgfrb_Prop_Str_Mdl1, Pdgfrb_Striatum, family = binomial())

Pdgfrb_Prop_Str_Prior1 <- 
  c(prior(normal(1,0.2), class = b, lb = 0),
    prior(student_t(3, 0, 2), class = Intercept, lb = 0))

# Fit model 1
Pdgfrb_Prop_Str_Fit1 <- 
  brm(
    data    = Pdgfrb_Striatum,
    family  = binomial(), 
    formula = Pdgfrb_Prop_Str_Mdl1,
    prior   = Pdgfrb_Prop_Str_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_Pdgfrb_LowHigh/Widefield_5x_Ipsilateral_Str_LowHigh_Fit1.rds",
    file_refit = "never") 
                     
# Add loo for model comparison
Pdgfrb_Prop_Strl_Fit1 <- 
  add_criterion(Pdgfrb_Prop_Str_Fit1, c("loo", "waic", "bayes_R2"))
```
## Model diagnostics

We perform model checks using the `pp_check` function:

```{r}
#| label: fig-Pdgfrb_LowHigh_Str_Diagnistics 
#| include: true
#| warning: false
#| message: false
#| results: false
#| cache: true
#| fig-cap: Model diagnostics for the probability of PDGFR-𝛽^high^ (striatal lesions)
#| fig-height: 4
#| fig-width: 5

set.seed(8807)

Pdgfrb_Prop_Str_Mdl1_pp <- 
  brms::pp_check(Pdgfrb_Prop_Str_Fit1, 
                 ndraws = 100) +
    labs(title = expression("Posterior predictive checks (PDGFR-β)"),
    subtitle = "Formula: PDGFR_React | PDGFR_Total ~ DPI") +
  Plot_theme  
 
Pdgfrb_Prop_Str_Mdl1_pp
```
Similar to our cortico-striatal injuries, our predictions exhibit a moderate deviations from the data. We further diagnose calculating Leave-One-Out Cross-Validation (Loo) using the `loo` function from `brms`.

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

loo(Pdgfrb_Prop_Str_Fit1)
```

The results show that most observations are problematic. In particular, 48% of observations belonging to the very bad category, suggest that the importance sampling weights for these observations are highly variable. This model does not have a good fit and must include additional variables. 

# Model results

## Visualization of conditional effects

```{r}
#| label: fig-Pdgfr_LowHigh_Str_CondEff
#| include: true
#| warning: false
#| message: false
#| results: false
#| fig-cap: Posterior for PDGFR-𝛽^low^ and PDGFR-𝛽^high^ (striatal lesions)
#| fig-width: 4
#| fig-height: 3

set.seed(8807)

# We create the graph for convex hull
Pdgfrb_Prop_Str_DPI <- 
  conditional_effects(Pdgfrb_Prop_Str_Fit1)

Pdgfrb_Prop_Str_DPI <- plot(Pdgfrb_Prop_Str_DPI, 
       plot = FALSE)[[1]]

Pdgfrb_Prop_Str_fig <- Pdgfrb_Prop_Str_DPI  + 
  scale_y_continuous(name = expression ("(p) Pdgfrb-β"^high)) +
  scale_x_continuous(name="DPI") +
  scale_color_manual(
    values = c("#0048BA", "red"),
    name="Condition"
    ) +
  scale_fill_manual(
    values = c("#0048BA", "red"),
    name="Condition"
    ) +

  Plot_theme +
  theme(legend.position = "top", legend.direction = "horizontal")

ggsave(
  plot     = Pdgfrb_Prop_Str_fig, 
  filename = "Plots/Widefield_5x_Ipsilateral_Pdgfrb_LowHigh/Widfield_5x_Ipsilateral_Pdgfrb_LowHigh_Str.png", 
  width    = 9, 
  height   = 9, 
  units    = "cm")

Pdgfrb_Prop_Str_fig
```

@fig-Pdgfr_LowHigh_Str_CondEff show an increasing probability for reactive PDGFR-β with a low uncertainty. This shows a clear contrast with cortico-striatal lesions which exhibit twice the probability at 14 and 30 DPI.

## Posterior summary

Here, we visualize the results in numerical terms using the `describe_posterior` function:

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

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

modelsummary(Pdgfrb_Prop_Str_Fit1, 
             shape = term ~ model + statistic,
             centrali2ty = "mean", 
             title = "PDGFR-β low and high-intensity populations in striatal-only lesions following MCAO",
             statistic = "conf.int",
             gof_omit = 'ELPD|ELDP s.e|LOOIC|LOOIC s.e|WAIC|RMSE',
             output = "Tables/html/Widefield_5x_Ipsilateral_Pdgfrb_LowHigh_Str_Table.html",
             )

Pdgfrb_Prop_Str_Fit1_Table <- modelsummary(Pdgfrb_Prop_Str_Fit1, 
             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_Prop_Str_Fit1_Table, filename = "Tables/tex/Widefield_5x_Ipsilateral_Pdgfrb_LowHigh_Str_Table.tex")
```


::: {#refs}
:::

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