Widefield_10x_ROIs_CD31-Pdgfrb_Coloc.qmd

---
title-block-banner: true
title: "Analysis of PDGFRβ attachment to vasculature (CD31+)"
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 vasculature
  - 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
    number-sections: true
    theme: spacelab
knitr:
  opts_chunk: 
    warning: false
    message: false
csl: science.csl
bibliography: references.bib
editor: 
  markdown: 
    wrap: 72
---

# Preview

This notebook reports the analysis of PDGFRβ localization in the
vasculature in defined ROIs (injured cortex and striatum, and
perilesional cortex).

**Parent dataset:** CD31 and PDGFRβ stained ROIs imaged at 10x using
widefield microscopy. These data set contain three major groups. First,
animals with cortico-striatal injuries grouped at 0 (sham), 3, 7, 14,
and 30 days post-ischemia (DPI). Second, animals with striatal injuries
at 14 and 30 DPI; and finally, a group of sham animals scarified at
during the time course to control for protein recombination after
tamoxifen injection. The raw images and pre-processing scripts (if
applicable) are available at the Zenodo repository
(10.5281/zenodo.10553084) under the name
`Widefield_10x_ROIs_CD31-Pdgfrb.zip`.

**Working dataset**: The
`Data_Raw/Raw_Widefield_10x_ROIs_CD31-Pdgfrb_Coloc.csv`data frame
containing the raw output from CellProfiller [@stirling2021]. The
pipeline used to perform PDGFRβ/CD31 colocalization is available at
https://osf.io/6ec89/.

We perform scientific inference based on the ratio of PDGFRβ cells
attached to the brain vasculature stained with CD31.

# 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", "emmeans", "ggplot","gtsummary",  "modelr", "modelsummary", "patchwork", "poorman", "tidybayes", "tidyverse", "viridis"))


library(bayesplot)
library(bayestestR)
library(brms)
library(dplyr)
library(easystats)
library(emmeans)
library(ggplot2)
library(gtsummary)
library(modelr)
library(modelsummary)
library(patchwork)
library(poorman)
library(plyr)
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
      ) 
```

# Load and handle the datasets

We load the `Data_Raw/Raw_Widefield_10x_ROIs_CD31-Pdgfrb_Coloc.csv` raw
data set to obtain the variables of interest.

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

Pdgfrb_CD31_Raw <- read.csv (file = "Data_Raw/Widefield_10x_ROIs_CD31-Pdgfrb/Raw_Widefield_10x_ROIs_CD31-Pdgfrb_Coloc.csv", header = TRUE)

# Eliminate unnecessary columns
Pdgfrb_CD31_Coloc <- subset(Pdgfrb_CD31_Raw, select = c(Count_Colocalization, Count_PDGFR_PrimaryObjects, FileName_PDGFR))

# Extract metadata information from image name
Pdgfrb_CD31_Coloc <- cbind(Pdgfrb_CD31_Coloc, do.call(rbind , strsplit(Pdgfrb_CD31_Coloc$FileName_PDGFR, "[_\\.]"))[,1:5])

# Eliminate File_Name column
Pdgfrb_CD31_Coloc <- subset(Pdgfrb_CD31_Coloc, select = -c(FileName_PDGFR))

# Change column names
colnames(Pdgfrb_CD31_Coloc) <- c("Pdgfrb_Perivascular", "Pdgfrb_Total", "AnimalID", "DPI", "Condition", "Lesion", "Region")
  
# Reordering the table
Pdgfrb_CD31_Coloc <- subset(Pdgfrb_CD31_Coloc, select = c("AnimalID", "DPI", "Condition", "Lesion", "Region", "Pdgfrb_Perivascular", "Pdgfrb_Total"))

# Create the variable for parenchymal cells
Pdgfrb_CD31_Coloc$Pdgfrb_Parenchymal <- Pdgfrb_CD31_Coloc$Pdgfrb_Total - Pdgfrb_CD31_Coloc$Pdgfrb_Perivascular

# Setting factors
Pdgfrb_CD31_Coloc$DPI <- factor(Pdgfrb_CD31_Coloc$DPI, levels = c("3D", "7D", "14D", "30D"))

Pdgfrb_CD31_Coloc$Region <- factor(Pdgfrb_CD31_Coloc$Region, levels = c("Peri", "Str", "Ctx"))

Pdgfrb_CD31_Coloc$Condition <- factor(Pdgfrb_CD31_Coloc$Condition, levels = c("Sham", "MCAO"))

Pdgfrb_CD31_Coloc$Lesion <- factor(Pdgfrb_CD31_Coloc$Lesion, levels = c("L0", "L1", "L2"))

# Create an additional DPI variable (numeric)

DPI_mapping <- c("0D" = "0", "3D" = "3", "7D" = "7", "14D" = "14", "30D" = "30")
Pdgfrb_CD31_Coloc$DPI2 <- as.numeric(DPI_mapping[as.character(Pdgfrb_CD31_Coloc$DPI)])

write.csv(Pdgfrb_CD31_Coloc, "Data_Processed/Widefield_10x_ROIs_CD31-Pdgfrb/Widefield_10x_ROIs_CD31-Pdgfrb_Coloc.csv", row.names = FALSE)
```

We save the `Wide10x_ROIs_CD31-Pdgfrb_Coloc.csv` containing the
following variables:

-   **AnimalID**: Unique animal identification
-   **DPI**: Days post-ischemia with factor levels (3D, 7D, 14D, 30D)
-   **Condition**: MCAO = Animals submitted to middle cerebral artery
    occlusion. Sham = Healthy animals. For Sham animals, DPI equalizes
    time of tamoxifen recombination to MCAO animals.
-   **Lesion**: Lesion type. L0 = No lesion. L1 = Cortico-striatal
    lesion. L2 = Striatal-only lesion.
-   **Region**: Brain region where imaging was performed. Ctx = Cortex
    (injured). Str = Striaum (injured). Peri = Perilesion (healthy).
-   **Pdgfr_Perivascular**: PDGFRβ+ cells attached (colocalizing) with
    CD31 (blood vessels)
-   **Pdgfr_Total**: Total number of PDGFRβ+ cells.
-   **Pdgfr_Parenchymal**: PDGFRβ+ cells not attached (colocalizing)
    with CD31 (blood vessels)
-   **DPI2**: Days-post ischemia as continuous (nemeric variable)

# Analysis of of PDGFRβ-CD31 colocalization in cortico-striatal lesions

First, our objective is to analyze the proportion of PDGFRβ/CD31
colocalization in animals with cortico-striatal lesions. To facilitate
the modeling and visualization of results using a binominal
distribution, we fit different models per brain region (Cortex, striatum
and perilesion).

## Analysis of PDGFRβ-CD31 colocalization in the cortex

We subset the dataset to exclude sham mice and animals with only
striatal injuries.

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

Pdgfrb_CD31_CtxMCAO <- filter(Pdgfrb_CD31_Coloc, Lesion == "L1", Region == "Ctx")
```

### Exploratory data visualization

We visualize the number of parenchymal PDGFRβ+ cells in the injured
cortex.

```{r}
#| label: fig-Pdgfrb_CtxMCAO_Exploratory
#| include: true
#| warning: false
#| message: false
#| results: false
#| fig-cap: Exploratory data visualization for PDGFRβ/CD31 colocalization
#| fig-width: 5
#| fig-height: 4

set.seed(8807)

Pdgfrb_Parenchymal_10x <- 
  ggplot(
    data  = Pdgfrb_CD31_CtxMCAO, 
    aes(x = DPI2, 
        y = Pdgfrb_Parenchymal)) +
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= expression("Number of parenchymal PDGFRβ cells")) +
  scale_x_continuous(name="Days post-ischemia (DPI) ",
                     breaks=c(0, 3, 7,14,30)) +
  Plot_theme 

Pdgfrb_Parenchymal_10x
```

@fig-Pdgfrb_CtxMCAO_Exploratory shows that the number of parenchymal
cells strongly increases in the cortex after 3 Days post ischemia and
seem the remain constant over time. However, to gain a more
comprehensive picture, we will model the number of parenchymal PDGFRβ
cells conditional on the total number of cells.

### Statistical modeling

We will fit a linear and a non-linear Bayesian model using a binomial
distribution to predict the proportion of parenchymal PDGFRβ cells by
DPI.

-   **Pdgfr_CtxMCAO_Mdl1** DPI as a linear predictor of PDGFRβ, with
    the following notation:

$$ 
P(Parenchymal | Total) \sim Binomial(n = Total, p) \\
\text{logit}(P(Y -1)) = \beta*0 +* \beta{DPI} \times DPI 
$$

Where $Y$ represents the occurrence of parenchymal cells, $P(Y=1)$ is
the probability of observing parenchymal cells, $\beta_0$ is the
intercept, and $\beta_{DPI}$ is the coefficient for the effect of DPI on
the log-odds of observing parenchymal cells.

Next, we incorporate a smooth term for DPI, allowing a non-linear
relationship between DPI and the log-odds of observing parenchymal cells
within the total number fo cells. The use of a smoothing function with k
= 4 represents a flexible, spline-based curve to model this
relationship:

-   **Pdgfr_CtxMCAO_Mdl2** Splines model on DPI, with the following
    notation:

$$ 
Parenchymal | Total \sim Binomial(p) \\
\text{logit}(p) = s(DPI2, k = 4) 
$$

Both models use default `brms` flat priors.

#### Fit the models

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

# Model 1: DPI as a single predictor
Pdgfrb_CtxMCAO_Mdl1 <- bf(Pdgfrb_Parenchymal | trials(Pdgfrb_Total) ~ DPI2)

get_prior(Pdgfrb_CtxMCAO_Mdl1, Pdgfrb_CD31_CtxMCAO, family = binomial())

# Fit model 1
Pdgfrb_CtxMCAO_Fit1 <- 
  brm(
    data    = Pdgfrb_CD31_CtxMCAO,
    family  = binomial(), 
    formula = Pdgfrb_CtxMCAO_Mdl1,
    chains  = 4,
    cores   = 4,
    warmup  = 2500, 
    iter    = 5000, 
    seed    = 8807,
    control = list(adapt_delta = 0.99, max_treedepth = 15),
    file    = "Models/Widefield_10x_ROIs_CD31-Pdgfrb_Coloc/Widefield_10x_ROIs_CD31-Pdgfrb_CtxMCAO_Fit1.rds",
    file_refit = "never") 
                     
# Add loo for model comparison
Pdgfrb_CtxMCAO_Fit1 <- 
  add_criterion(Pdgfrb_CtxMCAO_Fit1, c("loo", "waic", "bayes_R2"))


# Model 2: DPI with splines

Pdgfrb_CtxMCAO_Mdl2 <- bf(Pdgfrb_Parenchymal | trials(Pdgfrb_Total) ~ s(DPI2, k = 4))

get_prior(Pdgfrb_CtxMCAO_Mdl2, Pdgfrb_CD31_CtxMCAO, family = binomial())

# Fit model 2
Pdgfrb_CtxMCAO_Fit2 <- 
  brm(
    data    = Pdgfrb_CD31_CtxMCAO,
    family  = binomial(), 
    formula = Pdgfrb_CtxMCAO_Mdl2,
    knots   = list(DPI = c(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_10x_ROIs_CD31-Pdgfrb_Coloc/Widefield_10x_ROIs_CD31-Pdgfrb_CtxMCAO_Fit2.rds",
    file_refit = "never") 
                     
# Add loo for model comparison
Pdgfrb_CtxMCAO_Fit2 <- 
  add_criterion(Pdgfrb_CtxMCAO_Fit2, c("loo", "waic", "bayes_R2"))
```

#### Model comparison

We compare the fitted models using WAIC.

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

Pdgfrb_CtxMCAO_Comp <- 
  compare_performance(
    Pdgfrb_CtxMCAO_Fit1, 
    Pdgfrb_CtxMCAO_Fit2 
    )

Pdgfrb_CtxMCAO_Comp
```

Let's see it in graphical terms:

```{r}
#| label: fig-Pdgfrb_CtxMCAO_Compare
#| include: true
#| warning: false
#| message: false
#| results: false
#| #| fig-cap: Model coparison for PDGFRβ/CD31 colocalization
#| fig-width: 5
#| fig-height: 4

Pdgfrb_CtxMCAO_W <- 
loo_compare(
  Pdgfrb_CtxMCAO_Fit1, 
  Pdgfrb_CtxMCAO_Fit2, 
  criterion = "waic")

# Generate WAIC graph
Pdgfrb_CtxMCAO_WAIC <- 
  Pdgfrb_CtxMCAO_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_CtxMCAO_Fit1", 
             "Pdgfrb_CtxMCAO_Fit2"), 
             
    labels=c("Mdl1", 
             "Mdl2")) +
    
  coord_flip() +
  labs(x = "", 
       y = "WAIC (score)",
       title = "") +
  Plot_theme

Pdgfrb_CtxMCAO_WAIC
```

The graph shows the model with splines is far less penalized that the
linear model. This offer us good support to continue inference using the
model with splines.

#### Model diagnostics

We check the model fitting using posterior predictive checks

```{r}
#| label: Pdgfrb_CtxMCAO_Diagnistics
#| include: true
#| warning: false
#| message: false
#| results: false
#| fig-cap: Model diagnostics for PDGFRβ/CD31 colocalization (Cortex)
#| fig-width: 5
#| fig-height: 4

set.seed(8807)

Pdgfrb_CtxMCAO_Fit2_pp <- 
  brms::pp_check(Pdgfrb_CtxMCAO_Fit2, 
                 ndraws = 100) +
  labs(title = "Posterior predictive checks",
  subtitle = "Formula: Formula: PDGFR_Parenchymal | PDGFR_Total ~ s(DPI, k = 4)") +
  Plot_theme  
  
Pdgfrb_CtxMCAO_Fit2_pp
```

The predictions follow the same pattern that the observed data. However,
we can appreciate a moderate deviation in the density of the peak about
200. Nonetheless, considering the R2 = 0.83 from this model, we believe
it has a good prediction accuracy.

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


#launch_shinystan(Pdgfrb_CtxMCAO_Fit2)
```

### Model results

#### Visualization of conditional effects

We use the `conditiona_effects` function to see the posterior
distribution:

```{r}
#| label: fig-Pdgfrb_CtxMCAO_CE
#| include: true
#| warning: false
#| message: false
#| results: false
#| fig-cap: Posterior distribution for PDGFRβ/CD31 colocalization
#| fig-height: 4
#| fig-width: 5

set.seed(8807)

# We create the graph for convex hull
Pdgfrb_CtxMCAO_DPI <- 
  conditional_effects(Pdgfrb_CtxMCAO_Fit2)

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

Pdgfrb_CtxMCAO_fig <- Pdgfrb_CtxMCAO_DPI  + 
  scale_y_continuous(name = expression ("(P) parenchymal PDGFRβ")) +
  scale_x_continuous(name="DPI"                   ,
                     breaks = c(3, 10, 20, 30),
                     labels = c("3", "10", "20", "30")) +

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

ggsave(
  plot     = Pdgfrb_CtxMCAO_fig, 
  filename = "Plots/Widefield_10x_ROIs_CD31-Pdgfrb_Coloc/Widefield_10x_ROIs_CD31-Pdgfrb_CtxMCAO_Fit2.png", 
  width    = 9, 
  height   = 9, 
  units    = "cm")

Pdgfrb_CtxMCAO_fig
```

@fig-Pdgfrb_CtxMCAO_CE show an increasing probability for parenchymal
PDGFRβ+ cells with a peak during the second week post injury, and a
higher uncertainty thereafter.

#### Posterior summary

Next, We plot the posterior summary using the `describe_posterior`
function:

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

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

modelsummary(Pdgfrb_CtxMCAO_Fit2, 
             shape = term ~ model + statistic,
             centrali2ty = "mean", 
             title = "PDGFRβ+ parenchymal cells following MCAO",
             statistic = "conf.int",
             gof_omit = 'ELPD|ELDP s.e|LOOIC|LOOIC s.e|WAIC|RMSE',
             output = "Tables/html/Widefield_10x_ROIs_CD31-Pdgfrb_Ctx_Fit2_Table.html",
             )

Pdgfrb_CtxMCAO_Fit2_Table <- modelsummary(Pdgfrb_CtxMCAO_Fit2, 
             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_CtxMCAO_Fit2_Table, 
            filename = "Tables/tex/Widefield_10x_ROIs_CD31-Pdgfrb_Ctx_Fit2_Table.tex")
```

Please note that we are unable to calculate derivatives for binomial
models.

## Analysis of PDGFRβ-CD31 colocalization in the striatum

We perform the same analysis for the injured striatum. As done before,
we filter the data set to select relevant ROIs.

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

Pdgfrb_CD31_StrMCAO <- filter(Pdgfrb_CD31_Coloc, Lesion == "L1", Region == "Str")
```

### Exploratory data visualization

We visualize the number of parenchymal PDGFRβ+ cells in the injured
striatum.

```{r}
#| label: fig-Pdgfrb_StrMCAO_Exploratory
#| include: true
#| warning: false
#| message: false
#| results: false
#| fig-cap: Exploratory data visualization for PDGFRβ/CD31 colocalization
#| fig-width: 5
#| fig-height: 4

set.seed(8807)

Pdgfrb_Parenchymal_10x <- 
  ggplot(
    data  = Pdgfrb_CD31_StrMCAO, 
    aes(x = DPI2, 
        y = Pdgfrb_Parenchymal)) +
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= expression("Number of parenchymal PDGFRβ cells")) +
  scale_x_continuous(name="Days post-ischemia (DPI) ",
                     breaks=c(0, 3, 7,14,30)) +
  Plot_theme 

Pdgfrb_Parenchymal_10x
```

In this case, we see that linear and non-linear models are close.
However, we fit a splines model for consistency.

### Statistical modeling

We reproduce the same same approach employed with the cortex. This time
we focus exclusively in the model with splines.

#### Fit the model

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

set.seed(8807)

# Model 2: DPI with splines
Pdgfrb_StrMCAO_Mdl2 <- bf(Pdgfrb_Parenchymal | trials(Pdgfrb_Total) ~ s(DPI2, k = 4))

get_prior(Pdgfrb_StrMCAO_Mdl2, Pdgfrb_CD31_StrMCAO, family = binomial())

# Fit model 2
Pdgfrb_StrMCAO_Fit2 <- 
  brm(
    data    = Pdgfrb_CD31_StrMCAO,
    family  = binomial(), 
    formula = Pdgfrb_StrMCAO_Mdl2,
    knots   = list(DPI = c(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_10x_ROIs_CD31-Pdgfrb_Coloc/Widefield_10x_ROIs_CD31-Pdgfrb_StrMCAO_Fit1.rds",
    file_refit = "never") 
                     
# Add loo for model comparison
Pdgfrb_StrMCAO_Fit2 <- 
  add_criterion(Pdgfrb_StrMCAO_Fit2, c("loo", "waic", "bayes_R2"))
```

#### Model diagnostics

We check the model fitting using posterior predictive checks

```{r}
#| label: fig-Pdgfrb_StrMCAO_Diagnistics
#| include: true
#| warning: false
#| message: false
#| results: false
#| fig-cap: Model diagnostics for PDGFRβ/CD31 colocalization (Striatum)
#| fig-height: 4
#| fig-width: 5

set.seed(8807)

Pdgfrb_StrMCAO_Fit2_pp <- 
  brms::pp_check(Pdgfrb_StrMCAO_Fit2, 
                 ndraws = 100) +
  labs(title = "Posterior predictive checks",
  subtitle = "Formula: Formula: PDGFR_Parenchymal | PDGFR_Total ~ s(DPI, k = 4)") +
  Plot_theme  
  
Pdgfrb_StrMCAO_Fit2_pp
```

We observe no significant deviations from the data. We can explore
further the model using `shinystan`.

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

#launch_shinystan(Pdgfrb_StrMCAO_Fit2)
```

### Model results

#### Visualization of conditional effects

```{r}
#| label: fig-Pdgfrb_StrMCAO_CE
#| include: true
#| warning: false
#| message: false
#| results: false
#| fig-cap: Posterior distribution for PDGFRβ/CD31 colocalization
#| fig-height: 4
#| fig-width: 5

set.seed(8807)

# We create the graph for convex hull
Pdgfrb_StrMCAO_DPI <- 
  conditional_effects(Pdgfrb_StrMCAO_Fit2)

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

Pdgfrb_StrMCAO_fig <- Pdgfrb_StrMCAO_DPI  + 
  scale_y_continuous(name = expression ("(p) parenchymal PDGFRβ")) +
    scale_x_continuous(name="DPI"                   ,
                     breaks = c(3, 10, 20, 30),
                     labels = c("3", "10", "20", "30")) +

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

ggsave(
  plot     = Pdgfrb_StrMCAO_fig, 
  filename = "Plots/Widefield_10x_ROIs_CD31-Pdgfrb_Coloc/Widefield_10x_ROIs_CD31-Pdgfrb_StrMCAO_Fit2.png", 
  width    = 9, 
  height   = 9, 
  units    = "cm")

Pdgfrb_StrMCAO_fig
```

@fig-Pdgfrb_StrMCAO_CE show an increasing probability for parenchymal
PDGFRβ+ cells with a peak during the second weeks post injury in 0.25.

#### Posterior summary

Next, We plot the posterior summary using the `describe_posterior`
function:

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

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

modelsummary(Pdgfrb_StrMCAO_Fit2, 
             shape = term ~ model + statistic,
             centrali2ty = "mean", 
             title = "PDGFRβ+ parenchymal cells following MCAO",
             statistic = "conf.int",
             gof_omit = 'ELPD|ELDP s.e|LOOIC|LOOIC s.e|WAIC|RMSE',
             output = "Tables/html/Widefield_10x_ROIs_CD31-Pdgfrb_Str_Fit2_Table.html",
             )

Pdgfrb_StrMCAO_Fit2_Table <- modelsummary(Pdgfrb_StrMCAO_Fit2, 
             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_StrMCAO_Fit2_Table, 
            filename = "Tables/tex/Widefield_10x_ROIs_CD31-Pdgfrb_Str_Fit2_Table.tex")
```

## Analysis of PDGFRβ-CD31 colocalization in the perilesion

We analyze the perilesion using the same approach.We begin by filtering
the perilesion

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

Pdgfrb_CD31_PeriMCAO <- filter(Pdgfrb_CD31_Coloc, Lesion == "L1", Region == "Peri")
```

### Exploratory data visualization

We visualize the number of parenchymal PDGFRβ+ cells in the perilesion.

```{r}
#| label: fig-Pdgfrb_PeriMCAO_Exploratory
#| include: true
#| warning: false
#| message: false
#| results: false
#| fig-cap: Exploratory data visualization for PDGFRβ/CD31 colocalization
#| fig-width: 5
#| fig-height: 4

set.seed(8807)

Pdgfrb_Parenchymal_10x <- 
  ggplot(
    data  = Pdgfrb_CD31_PeriMCAO, 
    aes(x = DPI2, 
        y = Pdgfrb_Parenchymal)) +
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= expression("Number of parenchymal PDGFRβ cells")) +
  scale_x_continuous(name="Days post-ischemia (DPI) ",
                     breaks=c(0, 3, 7,14,30)) +
  Plot_theme 

Pdgfrb_Parenchymal_10x
```

@fig-Pdgfrb_PeriMCAO_Exploratory reveals that the number of parenchymal
PDGFRβ cells tends to remain constant.

### Statistical modeling

#### Fit the model

We fit a similar statistical model with splines per DPI.

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

set.seed(8807)

# Model 1: DPI with splines

Pdgfrb_PeriMCAO_Mdl2 <- bf(Pdgfrb_Parenchymal | trials(Pdgfrb_Total) ~ s(DPI2, k = 4))

get_prior(Pdgfrb_PeriMCAO_Mdl2, Pdgfrb_CD31_PeriMCAO, family = binomial())

# Fit model 2
Pdgfrb_PeriMCAO_Fit2 <- 
  brm(
    data    = Pdgfrb_CD31_PeriMCAO,
    family  = binomial(), 
    formula = Pdgfrb_PeriMCAO_Mdl2,
    knots   = list(DPI = c(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_10x_ROIs_CD31-Pdgfrb_Coloc/Widefield_10x_ROIs_CD31-Pdgfrb_PeriMCAO_Fit2.rds",
    file_refit = "never") 
                     
# Add loo for model comparison
Pdgfrb_PeriMCAO_Fit2 <- 
  add_criterion(Pdgfrb_PeriMCAO_Fit2, c("loo", "waic", "bayes_R2"))
```

#### Model diagnostics

We plot posterior predictive checks

```{r}
#| label: fig-Pdgfrb_PeriMCAO_Diagnistics
#| include: true
#| warning: false
#| message: false
#| fig-cap: Model diagnostics for PDGFRβ/CD31 colocalization (Perilesion)
#| fig-width: 5
#| fig-height: 4

set.seed(8807)

Pdgfrb_PeriMCAO_Fit2_pp <- 
  brms::pp_check(Pdgfrb_PeriMCAO_Fit2, 
                 ndraws = 100) +
  labs(title = "Posterior predictive checks",
  subtitle = "Formula: Formula: PDGFR_Parenchymal | PDGFR_Total ~ s(DPI, k = 4)") +
  Plot_theme  
  
Pdgfrb_PeriMCAO_Fit2_pp
```

We observe similar trend and no significant deviations from the
observations.

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

#launch_shinystan(Pdgfrb_PeriMCAO_Fit2)
```

### Model results

#### Visualization of conditional effects

```{r}
#| label: fig-Pdgfrb_PeriMCAO_CE
#| include: true
#| warning: false
#| message: false
#| fig-cap: Model diagnostics for PDGFRβ/CD31 colocalization (Perilesion)
#| fig-width: 5
#| fig-height: 4

set.seed(8807)

# We create the graph for convex hull
Pdgfrb_PeriMCAO_DPI <- 
  conditional_effects(Pdgfrb_PeriMCAO_Fit2)

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

Pdgfrb_PeriMCAO_fig <- Pdgfrb_PeriMCAO_DPI  + 
  scale_y_continuous(name = expression ("(p) parenchymal PDGFRβ")) +
    scale_x_continuous(name="DPI"                   ,
                     breaks = c(3, 10, 20, 30),
                     labels = c("3", "10", "20", "30")) +
  Plot_theme +
  theme(legend.position = "top", legend.direction = "horizontal")

ggsave(
  plot     = Pdgfrb_PeriMCAO_fig, 
  filename = "Plots/Widefield_10x_ROIs_CD31-Pdgfrb_Coloc/Widefield_10x_ROIs_CD31-Pdgfrb_PeriMCAO_Fit2.png", 
  width    = 9, 
  height   = 9, 
  units    = "cm")

Pdgfrb_PeriMCAO_fig
```

@fig-Pdgfrb_PeriMCAO_CE We observe that the probability for perenchymal
PDGFRβ+ cells fluctuates little from 0.9 to 0.3 in perilesional
regions.

#### Posterior summary

Next, We plot the posterior summary using the `describe_posterior`
function:

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

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

modelsummary(Pdgfrb_PeriMCAO_Fit2, 
             shape = term ~ model + statistic,
             centrali2ty = "mean", 
             title = "PDGFRβ+ parenchymal cells following MCAO",
             statistic = "conf.int",
             gof_omit = 'ELPD|ELDP s.e|LOOIC|LOOIC s.e|WAIC|RMSE',
             output = "Tables/html/Widefield_10x_ROIs_CD31-Pdgfrb_Peri_Fit2_Table.html",
             )

Pdgfrb_PeriMCAO_Fit2_Table <- modelsummary(Pdgfrb_PeriMCAO_Fit2, 
             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_PeriMCAO_Fit2_Table, 
            filename = "Tables/tex/Widefield_10x_ROIs_CD31-Pdgfrb_Peri_Fit2_Table.tex")
```

## Analysis of PDGFRβ cells in all regions

Here, we model together not the probability, but the number of PDGFRβ
cells in all regions to complement our inferences. As done, previously,
we filter the dataset to select rows with L1 animals.

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

Pdgfrb_CD31_L1MCAO <- filter(Pdgfrb_CD31_Coloc, Lesion == "L1")

Pdgfrb_CD31_L1MCAO$Region <- factor(Pdgfrb_CD31_L1MCAO$Region, levels = c("Peri", "Str", "Ctx"))
```

### Fit the statistical model

To model the number of PDGFRβ across regions and DPIs, we use a
hurdle_lognormal distribution. This model is suitable for data
exhibiting zeros (indicating no occurrence of the event) and positive
outcomes. Although the number of PDGFRβ (as object counts) is commonly
modeled with Poisson distribution, we want to account for overdispersion
with the hurdle_lognormal, considering the variance across regions. The
model takes the following formulation:

Hurdle Component for Zero vs. Positive Outcomes: 

$$
\text{logit}(Pr(Y > 0)) = \alpha_0 + \alpha_1 \cdot DPI + \alpha_2 \cdot Region + \alpha_3 \cdot (DPI \times Region)
$$ 
and the continuous Component for Positive Outcomes: 
$$
\log(Y) \sim Normal(\mu, \sigma) \quad \text{for} \quad Y > 0
$$ 
where $\mu$ (the mean on the log scale) is modeled as:

$$
\mu = \beta_0 + \beta_1 \cdot DPI + \beta_2 \cdot Region + \beta_3 \cdot (DPI \times Region)
$$ The model uses the default `brms` flat priors.

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

set.seed(8807)

# Model 1

Pdgfrb_L1MCAO_Mdl1 <- bf(Pdgfrb_Parenchymal ~ DPI * Region)
                  
get_prior(Pdgfrb_L1MCAO_Mdl1, Pdgfrb_CD31_L1MCAO, family = hurdle_lognormal())

# Fit model 1
Pdgfrb_L1MCAO_Fit1 <- 
  brm(
    data    = Pdgfrb_CD31_L1MCAO,
    family  = hurdle_lognormal(), 
    formula = Pdgfrb_L1MCAO_Mdl1,
    chains  = 4,
    cores   = 4,
    warmup  = 2500, 
    iter    = 5000, 
    seed    = 8807,
    control = list(adapt_delta = 0.99, max_treedepth = 15),
    file    = "Models/Widefield_10x_ROIs_CD31-Pdgfrb_Coloc/Widefield_10x_ROIs_CD31-Pdgfrb_L1MCAO_Fit1.rds",
    file_refit = "never") 
                     
# Add loo for model comparison
Pdgfrb_L1MCAO_Fit1 <- 
  add_criterion(Pdgfrb_L1MCAO_Fit1, c("loo", "waic", "bayes_R2"))

```

### Model diagnostics

We plot posterior predictive checks\>

```{r}
#| label: fig-Pdgfrb_L1MCAO_Diagnistics
#| include: true
#| warning: false
#| message: false
#| results: false
#| fig-cap: Model diagnostics for PDGFRβ/CD31 colocalization (Striatum)
#| fig-height: 4
#| fig-width: 5

set.seed(8807)

Pdgfrb_L1MCAO_Fit1_pp <- 
  brms::pp_check(Pdgfrb_L1MCAO_Fit1, 
                 ndraws = 100) +
  labs(title = "Posterior predictive checks",
  subtitle = "Formula: Formula: Pdgfrb_Parenchymal ~ DPI * Region") +
  scale_x_continuous(limits = c(0, 500)) +
  Plot_theme  
  
Pdgfrb_L1MCAO_Fit1_pp
```

@fig-Pdgfrb_L1MCAO_Diagnistics shows no the model predictions align with
the observed data. We can explore further the model using `shinystan`.

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

#launch_shinystan(Pdgfrb_L1MCAO_Fit1)
```

### Model results

#### Visualization of conditional effects

```{r}
#| label: fig-Pdgfrb_L1MCAO_CE 
#| include: true
#| warning: false
#| message: false
#| fig-cap: Posterior for the number of PDGFRβ+ cells in different brain regions
#| fig-height: 4
#| fig-width: 5

set.seed(8807)

# We create the graph for convex hull
Pdgfrb_L1MCAO_DPI <- 
  conditional_effects(Pdgfrb_L1MCAO_Fit1)

Pdgfrb_L1MCAO_DPI <- plot(Pdgfrb_L1MCAO_DPI, 
       plot = FALSE)[[3]]

Pdgfrb_L1MCAO_fig <- Pdgfrb_L1MCAO_DPI  + 
  scale_y_continuous(name = expression ("Parenchymal PDGFRβ")) +
  scale_x_discrete(name="DPI") +
  
    geom_point(data=Pdgfrb_CD31_L1MCAO, 
               aes(y = Pdgfrb_Parenchymal, 
                   x = DPI, colour=Region),
             inherit.aes=FALSE, 
             alpha=0.5,
             size = 1,
             position=position_jitter(h=0, w=0.07)) +
  
  scale_color_manual(
    values = c("#0048BA", "red", "darkgreen"),
    labels = c("Perilesion", "Striatum", "Cortex"),
    name="Region"
    ) +
  scale_fill_manual(
    values = c("#0048BA", "red", "darkgreen"),
    labels = c("Perilesion", "Striatum", "Cortex"),
    name="Region"
    ) +
  
  Plot_theme +
  theme(legend.position = c(0.2, 0.8), 
        legend.direction = "vertical")

ggsave(
  plot     = Pdgfrb_L1MCAO_fig, 
  filename = "Plots/Widefield_10x_ROIs_CD31-Pdgfrb_Coloc/Widefield_10x_ROIs_CD31-Pdgfrb_Coloc_L1MCAO_Fit1.png", 
  width    = 10, 
  height   = 9, 
  units    = "cm")

Pdgfrb_L1MCAO_fig
```

@fig-Pdgfrb_L1MCAO_CE shows that parenchymal cells are present mainly in
the cortex.

#### Calculate the contrast

Additionally, we can calculate the contrast between regions per time
point, focusing on Cortex-Striatum.

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

# We generate a data frame with the contrast
Pdgfrb_Region_Contrast <- Pdgfrb_L1MCAO_Fit1%>%
emmeans(~ Region + DPI, var = "Pdgfrb_Parenchymal", epred = TRUE) %>%
contrast(method = "revpairwise") %>%
gather_emmeans_draws() %>% sample_n(100)

# We select contrast of interest 
Pdgfrb_Region_Contrast_Sub <- Pdgfrb_Region_Contrast[
  (Pdgfrb_Region_Contrast$contrast=="Ctx 3D - Str 3D"|
   Pdgfrb_Region_Contrast$contrast=="Ctx 7D - Str 7D"|
     Pdgfrb_Region_Contrast$contrast=="Ctx 14D - Str 14D"|
   Pdgfrb_Region_Contrast$contrast=="Ctx 30D - Str 30D"),]

Pdgfrb_Region_Contrast_Sub$contrast <- 
  factor(Pdgfrb_Region_Contrast_Sub$contrast, 
         levels = c("Ctx 3D - Str 3D", "Ctx 7D - Str 7D", "Ctx 14D - Str 14D", "Ctx 30D - Str 30D")) 
```

Next, we plot the contrast using ggplot

```{r}
#| label: Pdgfrb_Gfap_PlotRegions2
#| include: true
#| warning: false
#| message: false
#| fig-cap: Contrast for the number of PDGFRβ+ cells in different brain regions
#| fig-height: 4
#| fig-width: 5

Pdgfrb_Region_Contrast_Fig <- 
  Pdgfrb_Region_Contrast_Sub %>%
  ggplot(
    aes(x    = .value, 
        y    = contrast)) +
  
  stat_slab() +
  
  stat_pointinterval(
    point_interval = mode_hdi, 
    position = position_dodge(width = .95, preserve = "single")) +

  scale_y_discrete(
    name= "Cortex-Striatum",
    labels = c("3D","7D",  "14D", "30D")) +
  
  scale_x_continuous(
   name="PDGFRβ (contrast)",
   limits=c(-100, 600), 
   breaks=seq(-100,600, 200)) +

  Plot_theme 
  
ggsave(
  plot     = Pdgfrb_Region_Contrast_Fig, 
  filename = "Plots/Widefield_10x_ROIs_CD31-Pdgfrb_Coloc/Widefiled_10x_ROIs_CD31-Pdgfrb_RegionContrast.png", 
  width    = 9, 
  height   = 9, 
  units    = "cm")

Pdgfrb_Region_Contrast_Fig

```

#### Posterior summary

We plot the posterior summary using the `describe_posterior` function:

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

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

modelsummary(Pdgfrb_L1MCAO_Fit1, 
             shape = term ~ model + statistic,
             centrali2ty = "mean", 
             title = "PDGFRβ+ parenchymal cells following MCAO",
             statistic = "conf.int",
             gof_omit = 'ELPD|ELDP s.e|LOOIC|LOOIC s.e|WAIC|RMSE',
             output = "Tables/html/Widefield_10x_ROIs_CD31-Pdgfrb_Regions_Fit1_Table.html",
             )

Pdgfrb_L1MCAO_Fit1_Table <- modelsummary(Pdgfrb_L1MCAO_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_L1MCAO_Fit1_Table, 
            filename = "Tables/tex/Widefield_10x_ROIs_CD31-Pdgfrb_Regions_Fit1_Table.tex")
```

# Analysis of of PDGFRβ/CD31 colocalization in striatal lesions

In this section, we analyze the PDGFRβ/CD31are colocalization in brain
with striatal lesions to evaluate the impact of cortical injures in the
appearance of PDGFRβ+ parenchymal cells. To achieve this goal, we
subset the "L1" and "L2" lesions in our data set at 14 and 30 DPI.

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

Pdgfrb_CD31_Striatum <- subset(Pdgfrb_CD31_Coloc, Lesion != "L0" & Region == "Str" & DPI %in% c("14D", "30D"))

Pdgfrb_CD31_Striatum$Lesion <- factor(Pdgfrb_CD31_Striatum$Lesion, levels = c("L2", "L1"))

Pdgfrb_CD31_Striatum$DPI <- factor(Pdgfrb_CD31_Striatum$DPI, levels = c("14D", "30D"))
```

## Exploratory data visualization

We perform exploratory data visualization to appreciate data trends

```{r}
#| label: fig-Pdgfrb_Striatum_Exploratory
#| include: true
#| warning: false
#| message: false
#| fig-cap: Exploratory data visualization for the number of PDGFRβ+ cells in striatal lesions
#| fig-height: 4
#| fig-width: 5


set.seed(8807)

### For perenchymal cells
Pdgfrb_Coloc_Str_fig <- 
  ggplot(
    data  = Pdgfrb_CD31_Striatum, 
    aes(x = DPI, 
        y = Pdgfrb_Parenchymal,
        color = Lesion)) +
  geom_boxplot()+
  scale_y_continuous(name= expression("Parenchymal PDGFRβ+")) +
  scale_x_discrete(name="Days post-ischemia (DPI) ",
                     breaks=c("14D","30D")) +
  Plot_theme 

Pdgfrb_Coloc_Str_fig
```

In particular, we observe that animals with cortico-striatal injuries
present extreme values, reflecting a major variability probably related
to infarct size.

## Statistical modeling

Given that we count on two data points, we will fit linear models.
First, we fit a model for the total number of PDGFRβ cells using a
student-t distribution, having the interaction of DPI and Lesion as
predictors variables. This model takes the following notation:

$$
Cells = \beta_{0} + \beta_{1} * DPI + \beta_{2} * Lesion + \beta_{3} * DPI * Lesion
$$ This model uses default `brms` flat priors.

Secondly, we fit a binominal (logit) model to predict the proportion of
parenchymal PDGFRβ cells conditional on the total number of PDGFRβ
cells. Similarly, this model have the interaction of DPI and Lesion as a
predictor variables. The model takes the following notation:

$$
Y | n ~ Binominal(n, \mu) \\ 
logit(\mu) = \beta_{0} + \beta_{1} * DPI + \beta_{2} * Lesion + \beta_{3} * DPI * Lesion
$$

Where:

Y represents the response variable "Pdgfr_Parenchymal".

n represent the total number of trials or observations corresponding to
the response variable "Pdgfr_Total".

$\mu$ represent the probability of success in each trial..

$\beta_{0}$ represents the intercept term.

$\beta_{1}$ represents the coefficient for the predictor variable "DPI".

$\beta_{2}$ represents the coefficient for the predictor variable
"Lesion".

$\beta_{3}$ represents the coefficient for the interaction term between
"DPI" and "Lesion".

This model uses beta(2,5) priors for the intercept and the class "b"
variables based on the information obtained from the striatum fit for
cortio-striatal injuries in this same notebook.

### Fit the models

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

# Model 1
Pdgfrb_Striatum_Mdl1 <- bf(Pdgfrb_Total ~ DPI * Lesion,
                           sigma ~ DPI * Lesion)

get_prior(Pdgfrb_Striatum_Mdl1, Pdgfrb_CD31_Striatum, family = student())

# Fit model 1
Pdgfrb_Striatum_Fit1 <- 
  brm(
    data    = Pdgfrb_CD31_Striatum,
    family  = student(), 
    formula = Pdgfrb_Striatum_Mdl1,
    chains  = 4,
    cores   = 4,
    warmup  = 2500, 
    iter    = 5000, 
    seed    = 8807,
    control = list(adapt_delta = 0.99, max_treedepth = 15),
    file    = "Models/Widefield_10x_ROIs_CD31-Pdgfrb_Coloc/Widefield_10x_ROIs_CD31-Pdgfrb_Striatum_Fit1.rds",
    file_refit = "never") 
                     
# Add loo for model comparison
Pdgfrb_Striatum_Fit1 <- 
  add_criterion(Pdgfrb_Striatum_Fit1, c("loo", "waic", "bayes_R2"))


# Model 2
Pdgfrb_Striatum_Mdl2 <- bf(Pdgfrb_Parenchymal | trials(Pdgfrb_Total) ~ DPI * Lesion)

get_prior(Pdgfrb_Striatum_Mdl2, Pdgfrb_CD31_Striatum, family = binomial())

Pdgfrb_Striatum_Prior1 <- 
  c(prior(beta(2, 5), class = b))

# Fit model 2
Pdgfrb_Striatum_Fit2 <- 
  brm(
    data    = Pdgfrb_CD31_Striatum,
    family  = binomial(), 
    formula = Pdgfrb_Striatum_Mdl2,
    #prior   = Pdgfrb_Striatum_Prior1,
    chains  = 4,
    cores   = 4,
    warmup  = 2500, 
    iter    = 5000, 
    seed    = 8807,
    control = list(adapt_delta = 0.99, max_treedepth = 15),
    file    = "Models/Widefield_10x_ROIs_CD31-Pdgfrb_Coloc/Widefield_10x_ROIs_CD31-Pdgfrb_Striatum_Fit2.rds",
    file_refit = "never") 
                     
# Add loo for model comparison
Pdgfrb_Striatum_Fit2 <- 
  add_criterion(Pdgfrb_Striatum_Fit2, c("loo", "waic", "bayes_R2"))
```

### Model diagnostics

We plot posterior predictive checks

```{r}
#| label: fig-Pdgfrb_Striatum 
#| include: true
#| warning: false
#| message: false
#| fig-cap: Model diagnostics for the number of PDGFRβ+ cells in striatal lesions
#| fig-height: 5
#| fig-width: 12

set.seed(8807)

# For model 1
Pdgfrb_Striatum_Fit1_pp <- 
  brms::pp_check(Pdgfrb_Striatum_Fit1, 
                 ndraws = 100) +
  labs(title = "Posterior predictive checks",
  subtitle = "Formula: Pdgfrb_Parenchymal ~ DPI * Region, sigma ~ Region") +
  scale_x_continuous(limits = c(0, 1000)) +
  Plot_theme  
  

# For model 2
Pdgfrb_Striatum_Fit2_pp <- 
  brms::pp_check(Pdgfrb_Striatum_Fit2, 
                 ndraws = 100) +
  labs(title = "Posterior predictive checks",
  subtitle = "Formula: Pdgfrb_Parenchymal | trials(Pdgfrb_Total) ~ DPI2 * Lesion") +
  scale_x_continuous(limits = c(0, 200)) +
  Plot_theme  
  
Pdgfrb_Striatum_Fit1_pp | Pdgfrb_Striatum_Fit2_pp
```

The predictions follow the overall trend and do not show mayor
deviations from the data. However, please note that the predictions of
the second model have a moderate deviation from the offset. We can
explore further the model using `shinystan`.

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

#launch_shinystan(Pdgfr_Striatum_Fit1)
#launch_shinystan(Pdgfr_Striatum_Fit2)
```

## Model results

### Visualization of conditional effects

```{r}
#| label: fig-Pdgfrb_Striatum_CE 
#| include: true
#| warning: false
#| message: false
#| fig-cap: Posterior for the number of PDGFRβ+ cells in striatal lesions
#| fig-height: 5
#| fig-width: 12

set.seed(8807)

# We create the graph for the total number of cells
Pdgfrb_Striatum_DPI <- 
  conditional_effects(Pdgfrb_Striatum_Fit1)

Pdgfrb_Striatum_DPI <- plot(Pdgfrb_Striatum_DPI, 
       plot = FALSE)[[3]]

Pdgfrb_Striatum_fig1 <- Pdgfrb_Striatum_DPI  + 
  scale_y_continuous(name = expression ("PDGFRβ cells")) +
  scale_x_discrete(name="DPI") +
  scale_color_manual(
    values = c("#0048BA", "red"),
    labels = c("Striatal", "Cortico-Striatal"),
    name="Lesion"
    ) +
  scale_fill_manual(
    values = c("#0048BA", "red"),
    labels = c("Striatal", "Cortico-Striatal"),
    name="Lesion"
    ) +
  
  Plot_theme +
  theme(legend.position = "top", legend.direction = "horizontal")

ggsave(
  plot     = Pdgfrb_Striatum_fig1, 
  filename = "Plots/Widefield_10x_ROIs_CD31-Pdgfrb_Coloc/Widefield_10x_ROIs_CD31-Pdgfrb_Striatum_Fit1.png", 
  width    = 10, 
  height   = 9, 
  units    = "cm")


# We create the graph for sigma
Pdgfrb_Striatum_Sigma <- 
  conditional_effects(Pdgfrb_Striatum_Fit1, dpar = "sigma")

Pdgfrb_Striatum_Sigma <- plot(Pdgfrb_Striatum_Sigma, 
       plot = FALSE)[[3]]

Pdgfrb_StriatumSigma_fig1 <- Pdgfrb_Striatum_Sigma  + 
  scale_y_continuous(name = expression ("PDGFRβ cells (Sigma)")) +
  scale_x_discrete(name="DPI") +
  scale_color_manual(
    values = c("#0048BA", "red"),
    labels = c("Striatal", "Cortico-Striatal"),
    name="Lesion"
    ) +
  scale_fill_manual(
    values = c("#0048BA", "red"),
    labels = c("Striatal", "Cortico-Striatal"),
    name="Lesion"
    ) +
  
  Plot_theme +
  theme(legend.position = "top", legend.direction = "horizontal")

ggsave(
  plot     = Pdgfrb_StriatumSigma_fig1, 
  filename = "Plots/Widefield_10x_ROIs_CD31-Pdgfrb_Coloc/Widefield_10x_ROIs_CD31-Pdgfrb_Striatum_Fit1_Sigma.png", 
  width    = 10, 
  height   = 9, 
  units    = "cm")


# We create the graph for the probability of parenchymal
Pdgfrb_Striatum_DPI <- 
  conditional_effects(Pdgfrb_Striatum_Fit2)

Pdgfrb_Striatum_DPI <- plot(Pdgfrb_Striatum_DPI, 
       plot = FALSE)[[3]]

Pdgfrb_Striatum_fig2 <- Pdgfrb_Striatum_DPI  + 
  scale_y_continuous(name = expression ("(p) parenchymal PDGFRβ")) +
  scale_x_discrete(name="DPI",
                     breaks = c("14D", "30D"),
                     labels = c("14D", "30D")) +
  scale_color_manual(
    values = c("#0048BA", "red"),
    labels = c("Striatal", "Cortico-Striatal"),
    name="Lesion"
    ) +
  scale_fill_manual(
    values = c("#0048BA", "red"),
    labels = c("Striatal", "Cortico-Striatal"),
    name="Lesion"
    ) +
  
  Plot_theme +
  theme(legend.position = "top", legend.direction = "horizontal")

ggsave(
  plot     = Pdgfrb_Striatum_fig2, 
  filename = "Plots/Widefield_10x_ROIs_CD31-Pdgfrb_Coloc/Widefield_10x_ROIs_CD31-Pdgfrb_Striatum_Fit2.png", 
  width    = 10, 
  height   = 9, 
  units    = "cm")

Pdgfrb_Striatum_fig1 | Pdgfrb_StriatumSigma_fig1 | Pdgfrb_Striatum_fig2
```

@fig-Pdgfrb_Striatum_CE shows a substantial overlap in the total number
of PDGFRβ in both conditions. The decrease in the uncertainty at 30 DPI
for cortico-striatal lesions is notable Interestingly, cortico-striatal
lesion display a lower probability for parenchymal PDGFRβ cells at 30
DPI. This contrast sharply with the increasing trend showed by striatal
lesion. We speculate this reduction in the proportion of parenchymal
cells may be guided by increased cell mortality in cortico-striatal
lesions.

Next, we visualize the posterior summaries

### Posterior summary

Next, We plot the posterior summary using the `describe_posterior`
function:

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

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

modelsummary(Pdgfrb_Striatum_Fit1, 
             shape = term ~ model + statistic,
             centrali2ty = "mean", 
             title = "PDGFRβ+ cells following MCAO in striatal-only regions",
             statistic = "conf.int",
             gof_omit = 'ELPD|ELDP s.e|LOOIC|LOOIC s.e|WAIC|RMSE',
             output = "Tables/html/Widefield_10x_ROIs_CD31-Pdgfrb_Striatum_Fit1_Table.html",
             )

Pdgfrb_Striatum_Fit1_Table <- modelsummary(Pdgfrb_Striatum_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_Striatum_Fit1_Table, 
            filename = "Tables/tex/Widefield_10x_ROIs_CD31-Pdgfrb_Striatum_Fit1_Table.tex")


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

modelsummary(Pdgfrb_Striatum_Fit2, 
             shape = term ~ model + statistic,
             centrali2ty = "mean", 
             title = "Probability of parenchymal PDGFRβ+ cells following MCAO",
             statistic = "conf.int",
             gof_omit = 'ELPD|ELDP s.e|LOOIC|LOOIC s.e|WAIC|RMSE',
             output = "Tables/html/Widefield_10x_ROIs_CD31-Pdgfrb_Striatum_Fit2_Table.html",
             )

Pdgfrb_Striatum_Fit2_Table <- modelsummary(Pdgfrb_Striatum_Fit2, 
             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_Striatum_Fit2_Table, 
            filename = "Tables/tex/Widefield_10x_ROIs_CD31-Pdgfrb_Striatum_Fit2_Table.tex")
```

# Analisys of Sham/MCAO animals (effect of tamoxifen recombination)

To estimate the effects of tamoxifen recombination in the number of
parenchymal PDGFR-B cells, we analyze the cortical regions from MCAO and
Sham animals in a single Bayesian model.

## Subsetting the data

We filter the data set to select rows with cortical ROIs for MCAO and
Sham animals.

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

Pdgfrb_CD31_ShamMCAO <- filter(Pdgfrb_CD31_Coloc, Lesion != "L2", Region == "Peri")

Pdgfrb_CD31_ShamMCAO$Condition <- factor(Pdgfrb_CD31_ShamMCAO$Condition, levels = c("Sham", "MCAO"))
```

## Statistical modeling

As we want to evaluate the effects of recombination in each time point,
we fit a linear model using a Poisson distribution. Here, we fit the
model for the total number of parenchymal PDGFRβ cells conditioning on
the interaction between DPI and condition (sham or MCAO).The model takes
the following notation:

$$
Pdgfrb\_Parenchymal \sim Poisson(\lambda) \\
\log(\lambda) = \beta_0 + \beta_1 \cdot DPI + \beta_2 \cdot Condition + \beta_3 \cdot (DPI \times Condition)
$$

Here, $\log$ is the natural logarithm, serving as the link function that
relates the linear predictors to the mean of the Poisson distribution.
This log link ensures that the predictions for $\lambda$ are always
positive, as required for count data.

The model takes the following weakly informative priors:

$$
\beta_{0} \sim Normal(5,2) \\
\beta_{2} \sim Normal(0,2.5) \\
$$ 

### Fit the model

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

set.seed(8807)

Pdgfrb_ShamMCAO_Mdl1 <- bf(Pdgfrb_Parenchymal ~ DPI * Condition)
                          
get_prior(Pdgfrb_ShamMCAO_Mdl1, Pdgfrb_CD31_ShamMCAO, family = student)

# Model prior
Pdgfrb_ShamMCAO_Prior1 <- 
  c(prior(normal(5, 2), class = "Intercept", lb = 0),
    prior(normal(0, 2), class = b))
  

# Fit model 1
Pdgfrb_ShamMCAO_Fit1 <- 
  brm(
    data    = Pdgfrb_CD31_ShamMCAO,
    family  = poisson, 
    formula = Pdgfrb_ShamMCAO_Mdl1,
    prior   = Pdgfrb_ShamMCAO_Prior1,
    chains  = 4,
    cores   = 4,
    warmup  = 2500, 
    iter    = 5000, 
    seed    = 8807,
    control = list(adapt_delta = 0.99, max_treedepth = 15),
    file    = "Models/Widefield_10x_ROIs_CD31-Pdgfrb_Coloc/Widefield_10x_ROIs_CD31-Pdgfrb_ShamMCAO_Fit1",
    file_refit = "never") 
                     
# Add loo for model comparison
Pdgfrb_ShamMCAO_Fit1 <- 
  add_criterion(Pdgfrb_ShamMCAO_Fit1, c("loo", "waic", "bayes_R2"))
```

### Model diagnostics

We plot posterior predictive checks

```{r}
#| label: Pdgfrb_ShamMCAO_Diagnistics
#| include: true
#| warning: false
#| message: false
#| fig-cap: Model diagnostics for PDGFRβ in sham and MCAO animals
#| fig-width: 5
#| fig-height: 4

set.seed(8807)

Pdgfrb_ShamMCAO_Fit1_pp <- 
  brms::pp_check(Pdgfrb_ShamMCAO_Fit1, 
                 ndraws = 100) +
  labs(title = "Posterior predictive checks",
  subtitle = "Formula: Formula: Pdgfrb_Parenchymal ~ DPI * Condition,
                           sigma ~ DPI * Condition") +
  Plot_theme  
  
Pdgfrb_ShamMCAO_Fit1_pp
```

We observe that the model predictions are in agreement with the data
points. We can explore further the model using `shinystan`.

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

#launch_shinystan(Pdgfrb_ShamMCAO_Fit1)
```

## Model results

### Visualizing the results

```{r}
#| label: fig-Pdgfrb_ShamMCAO_CE
#| include: true
#| warning: false
#| message: false
#| results: false
#| fig-cap: Posterior distribution for PDGFRβ in sham and MCAO animals
#| fig-width: 5
#| fig-height: 4

set.seed(8807)

# We create the graph for convex hull
Pdgfrb_ShamMCAO_DPI <- 
  conditional_effects(Pdgfrb_ShamMCAO_Fit1)

Pdgfrb_ShamMCAO_DPI <- plot(Pdgfrb_ShamMCAO_DPI, 
       plot = FALSE)[[3]]

Pdgfrb_ShamMCAO_fig <- Pdgfrb_ShamMCAO_DPI  + 
  scale_y_continuous(name = expression ("Parenchymal PDGFRβ")) +
  scale_x_discrete(name="DPI",
                     breaks = c("3D", "7D", "14D", "30D")) +
                   
  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_ShamMCAO_fig, 
  filename = "Plots/Widefield_10x_ROIs_CD31-Pdgfrb_Coloc/Widefield_10x_ROIs_CD31-Pdgfrb_ShamMCAO_fig.png", 
  width    = 9, 
  height   = 9, 
  units    = "cm")

Pdgfrb_ShamMCAO_fig
```

@fig-Pdgfrb_ShamMCAO_CE reveals that time has no no effect on
recombination, and therefore, the appearance of parenchymal PDGFRβ
cells.

### Posterior summary

Next, We plot the posterior summary using the `describe_posterior`
function:

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

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

modelsummary(Pdgfrb_ShamMCAO_Fit1, 
             shape = term ~ model + statistic,
             centrali2ty = "mean", 
             title = "Parenchymal PDGFRβ+ cells in MCAO and sham animals",
             statistic = "conf.int",
             gof_omit = 'ELPD|ELDP s.e|LOOIC|LOOIC s.e|WAIC|RMSE',
             output = "Tables/html/Widefield_10x_ROIs_CD31-Pdgfrb_ShamMCAO_Fit1_Table.html",
             )

Pdgfrb_ShamMCAO_Fit1_Table <- modelsummary(Pdgfrb_ShamMCAO_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_ShamMCAO_Fit1_Table, 
            filename = "Tables/tex/Widefield_10x_ROIs_CD31-Pdgfrb_ShamMCAO_Fit1_Table.tex")
```

# References

::: {#refs}
:::

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