Widefield_20x_ROIs_Pdgfrb_Haralick.qmd

---
title-block-banner: true
title: "Analysis of Haralick features from PDGFRβ in defined ROIs"
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: 
  - Haralick features
  - Image analysis
  - 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
---

# Preview

This notebook reports the analysis of the area and Haralick features for PDGFRβ in defined ROIs of the ipsilateral hemisphere.

**Parent dataset:** PDGFRβ (td-tomato) Z-plane ROIs imaged at 20x. Images were mas intensity projected to obtain a single plane image. The images are grouped at 0, 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_20x_ROIs_Pdgfrb.zip`.

**Working dataset**: This notebook works with two different datasets. First, the `Data_Raw/Widefield_20x_ROIs_Pdgfrb_Cells/CellProfiler/Image.csv`containing the quantification of the area covered by the detected cells after processing using CellProfiler [@stirling2021]. The pipeline is available at https://osf.io/qmyk4.

On the other hand, we load the `Data_Raw/Widefield_20x_ROIs_Pdgfrb_Cells/Raw_Widefield_20x_ROIs_Pdgfrb_Haralick.csv` containing the computed Haralick features using the `Widefield_20x_ROIs_Pdgfrb_Haralick_BatchScript` notebook.

# 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","gtsummary", "modelbased", "modelr","modelsummary", "patchwork", "poorman","plyr", "spatstat", "tidybayes", "tidyverse"))

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

# Analysis PDGFRβ+ covered area

## Load and handle the datasets

We load the `Data_Raw/Widefield_20x_ROIs_Pdgfrb/CellProfiler/Image.csv` dataset containing the total objects area exported by CellProfiler.

```{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_Area <- read.csv(file = "Data_Raw/Widefield_20x_ROIs_Pdgfrb/CellProfiler/Image.csv", header = TRUE)

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

From the table, we are interested in the `AreaOccupied_AreaOccupied_Cells`and the `FileName_Raw`columns containing the measured area and metadata information from the image names. Next, we subset the dataset to select the columns of interest and give them meaningful names. Please note that to facilitate the reading of results, we downscale the pixel area count by thousand.

When performing pre-processing, we noticed that the CellProfiler pipeline yielded 0 cell count for three images because an unknown error. We exclude this images with biased counts.

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

## We subset the relevant columns (cell number)
Pdgfrb_Area <- subset(Pdgfrb_Area, select = c("FileName_Raw", "AreaOccupied_AreaOccupied_Cells"))

## And extract metadata from the image name
Pdgfrb_Area  <- cbind(Pdgfrb_Area, do.call(rbind , strsplit(Pdgfrb_Area$FileName_Raw, "[_\\.]"))[,1:5])

Pdgfrb_Area <- subset(Pdgfrb_Area, select = -c(FileName_Raw))

## We Rename the relevant columns 
colnames(Pdgfrb_Area) <- c("Area", "MouseID", "DPI", "Condition", "Lesion", "Region")

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

Pdgfrb_Area$Region <- factor(Pdgfrb_Area$Region, levels = c("Peri", "Str", "Ctx"), ordered = TRUE)

Pdgfrb_Area$Condition <- factor(Pdgfrb_Area$Condition, levels = c("SHAM", "MCAO"))

Pdgfrb_Area$Lesion <- factor(Pdgfrb_Area$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_Area$DPI_cont <- as.numeric(DPI_mapping[as.character(Pdgfrb_Area$DPI)])

## To facilitate reading
Pdgfrb_Area$Area <- Pdgfrb_Area$Area/1000

# Exclued error-derived 0 counts
Pdgfrb_Area <- filter(Pdgfrb_Area, Area != "0")

write.csv(Pdgfrb_Area, "Data_Processed/Widefield_20x_ROIs_Pdgfrb/Widefield_20x_ROIs_Pdgfrb_Area.csv", row.names = FALSE)

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

Next, we perform exploratory data visualization.

### Exploratory data visualization

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

```{r}
#| label: fig-Pdgfrb_Exploratory
#| include: true
#| warning: false
#| message: false
#| results: false
#| fig-cap: Exploratory data visualization for PDGFRβ area
#| fig-width: 12
#| fig-height: 5

set.seed(8807)

Pdgfrb_Area_EDV <- 
  ggplot(
    data  = Pdgfrb_Area, 
    aes(x = DPI_cont, 
        y = Area)) +
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 +
  facet_wrap(~ Region)

Pdgfrb_Area_EDV
```

@fig-Pdgfrb_Exploratory shows the distribution of area per DPI and regions. As expected, the cortical regions present the most prominent increase in area. The differences in the fit show little variation for linear and non linear models.Given that our objective is to capture the main differences in the trends between regions we will develop a linear model.

# Statistical modeling for PDGFRβ area

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 model with a student_t distribution:

-   **Pdgfrb_Area_Mdl1:** We use `DPI` and Region as a linear predictors for PDGFRβ Area. The model takes the following notation:

$$
Area = \beta_{0} + \beta_{1} * DPI + \beta_{2} * Region + \beta_{3} * DPI * Region + \epsilon_{i}
$$

This model uses the following weakly informative priors:

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

## Fit the models for PDGFRβ area

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

# Model 1: DPI and Region as a linear predictors
###########################################
Pdgfrb_Area_Mdl1 <- bf(Area ~ DPI_cont * Region)

get_prior(Pdgfrb_Area_Mdl1, Pdgfrb_Area)

Pdgfrb_Area_Prior1 <- 
  c(prior(normal(0,100), class = b),
    prior(student_t(3, 150,50), class = Intercept, lb= 0),
    prior(student_t(3, 0, 100), class = sigma, lb=0))

# Fit model 1
Pdgfrb_Area_Fit1 <- 
  brm(
    family  = student,  
    data    = Pdgfrb_Area,
    formula = Pdgfrb_Area_Mdl1,
    prior   = Pdgfrb_Area_Prior1,
    chains  = 4,
    cores   = 4,
    warmup  = 2500, 
    iter    = 5000, 
    seed    = 8807,
    control = list(adapt_delta = 0.99, max_treedepth = 15),
    file    = "Models/Widefield_20x_ROIs_Pdgfrb_Haralick/Widefield_20x_ROIs_Pdgfrb_Area_Fit1.rds",
    file_refit = "never") 
                     
# Add loo for model comparison
Pdgfrb_Area_Fit1 <- 
  add_criterion(Pdgfrb_Area_Fit1, c("loo", "waic", "bayes_R2"))
```

## Model diagnostics

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

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

set.seed(8807)

color_scheme_set("darkgray")

# Model 4
Pdgfrb_Area_Mdl1_pp <- 
  brms::pp_check(Pdgfrb_Area_Fit1, 
                 ndraws = 100) +
  labs(title = "Posterior predictive checks (model 4)",
  subtitle = "Formula: Area ~ DPI * Region") +
  scale_x_continuous(limits=c(-200, 1000)) +
  Plot_theme 

Pdgfrb_Area_Mdl1_pp
```

@fig-Pdgfrb_Area_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_Area_Fit1)
```

# Model results

### Visualizing the results

We visualize the results using the `conditional_effects` function form `brms`.

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

set.seed(8807)


# We create the graph for convex hull
Pdgfrb_Area_DPI <- 
  conditional_effects(Pdgfrb_Area_Fit1)

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

Pdgfrb_Area_fig <- Pdgfrb_Area_DPI  + 
  scale_y_continuous(name = expression ("PDGFRβ area (pixels)")) +
  scale_x_continuous(name="DPI"                   ,
                     breaks = c(0, 3, 7, 14, 30),
                     labels = c("0", "3", "7", "14", "30")) +
  
  geom_point(data=Pdgfrb_Area, 
               aes(y = Area, 
                   x = DPI_cont, 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.4, 1), 
        legend.direction = "horizontal")

ggsave(
  plot     = Pdgfrb_Area_fig, 
  filename = "Plots/Widefield_20x_ROIs_Pdgfrb_Haralick/Widefield_20x_ROIs_Pdgfrb_Area.png", 
  width    = 10, 
  height   = 9, 
  units    = "cm")

Pdgfrb_Area_fig
```

@fig-Pdgfrb_Area_CondEff shows the fit for PDGFRβ area. We can visualize an increasing trend, specially for cortex.

## Posterior summary

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

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

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

modelsummary(Pdgfrb_Area_Fit1, 
             shape = term ~ model + statistic,
             centrali2ty = "mean", 
             title = "PDGFRβ are in defined ROIs following MCAO",
             statistic = "conf.int",
             gof_omit = 'ELPD|ELDP s.e|LOOIC|LOOIC s.e|WAIC|RMSE',
             output = "Tables/html/Widefield_20x_ROIs_Pdgfrb_Area_Fit1_Table.html",
             )

Pdgfrb_Area_Fit1_Table <- modelsummary(Pdgfrb_Area_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_Area_Fit1_Table, 
            filename = "Tables/tex/Widefield_20x_ROIs_Pdgfrb_Area_Fit1_Table.tex")
```

With this summary we can visualize estimates and the uncertainty for each parameter. Specifically, we can observe a difference of 14.7 (CI95% 12.6 - 16.6) for the slope of cortex, denoting a marked difference compared to perilesion and striatum.

# Analysis of Haralick features for PDGFRβ (Texture analysis)

## Load the dataset

We load and handle the `Raw_Widefield_20x_ROIs_Pdgfrb_Haralick.csv` dataset containing the Haralick features calculated using the `mahotas` library for Python [@mahotas]. Here, We subset the data and make it usable.

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


# Load the data
Pdgfrb_Haralick <- read.csv(file = 'Data_Raw/Widefield_20x_ROIs_Pdgfrb/Raw_Widefield_20x_ROIs_Pdgfrb_Haralick.csv', header = TRUE)

# Subset the data frame (Eliminate filename column)
Pdgfrb_Haralick <- subset(Pdgfrb_Haralick, select = -c(Filename))

Pdgfrb_Haralick $DPI <- factor(Pdgfrb_Haralick$Day, levels = c("0D", "3D", "7D", "14D", "30D"))

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

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

Here, we define the Haralick features for reference:

-   **Angular Second Moment (ASM) / Energy:** Measures the uniformity of an image. A higher value indicates that the image has more uniform textures or constant regions.

-   **Contrast:** Represents the difference between the highest and the lowest intensity value in the co-occurrence matrix. It measures the amount of local variations present in an image.

-   **Correlation:** Measures the joint probability occurrence of the specified pixel pairs. It provides information about the linear dependency of gray levels in the neighboring pixels.

-   **Sum of Squares / Variance:** It provides a measure of the squared differences from the mean intensity value.

-   **Inverse Difference Moment (IDM) / Homogeneity:** Represents the local homogeneity of an image. The values are high when the local textures are consistent or homogeneous.

-   **Sum Average:** Represents the average intensity value of the co-occurrence matrix.

-   **Sum Variance:** Measures the variance of the sum of the intensity values from the average value in the co-occurrence matrix.

-   **Sum Entropy:** Represents the randomness or complexity in the sum of the intensity values of the co-occurrence matrix.

-   **Entropy:** Provides a measure of the randomness or complexity in the image. Higher values indicate more complex textures.

-   **Difference Variance:** Represents the variance in the differences between the intensity values of pairs of pixels.

-   **Difference Entropy:** Measures the randomness or complexity in the differences between the intensity values of pairs of pixels.

-   **Informational Correlation 1 (Info Corr 1):** Represents the correlation between the occurrence of the specified pixel pairs and their average intensity values.

-   **Informational Correlation 2 (Info Corr 2):** Provides another measure of the correlation between the occurrence of the specified pixel pairs and their average intensity values. It's typically more sensitive to changes than Info Corr 1.

## Principal component analysis (PCA)

Next, we perform principal component analysis (PCA) as a dimensionality reduction method to for multiple variables. We will use the `factoMineR` [@lê2008; @FactoMineR] and `factoextra` [@factoextra] packages to have access to the PCA and visualization functions.

PCA uses an orthogonal transformation to convert a set of possibly correlated observations into a set of values of linearly uncorrelated variables called principal components. The number of principal components is generally less than the number of original variables. Here, the first principal component accounts for the higher variability in the data. This allow the transformation of a large set of variables into a smaller inputs containing all the information.

### Define the PCA variables

The dimensionality reduction is applicable exclusively to numerical. Here, we select the relevant columns. We specify `center=T` and `scale=T` arguments into the `scale` function for scaling of all variables.

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

# Subset variables from the data set
Pdgfrb_Haralick_Variables <- data.frame(scale(Pdgfrb_Haralick[, c(5:17)], center=T, scale=T))
```

### Obtain a correlation matrix

In this section, we generate a correlation matrix using to `ggcorrplot` package [@ggcorrplot] to reveal the proximity between different variables. The higher the value, the higher the correlation between two variables. Conversely, the closer is to -1, the most negatively correlated they are.

```{r}
#| label: fig-Pdgfrb_Haralick_Matrix
#| include: true
#| warning: false
#| message: false
#| #| fig-width: 6
#| fig-height: 6

# Create the correlation matrix
Pdgfrb_Haralick_corr <- cor(Pdgfrb_Haralick_Variables)

# Plot and save the correlation matrix
Pdgfrb_Haralick_corr_fig <- ggcorrplot(Pdgfrb_Haralick_corr)

ggsave(
  plot     = Pdgfrb_Haralick_corr_fig, 
  filename = "Plots/Widefield_20x_ROIs_Pdgfrb_Haralick/Widefield_20x_ROIs_Pdgfrb_Haralick_corr_fig.png", 
  width    = 12, 
  height   = 10, 
  units    = "cm")

Pdgfrb_Haralick_corr_fig 
```

### Perfrom the PCA

Now, we execute the principal component analysis using the `princomp` function. Then, we print graphical and numerical summaries.We use the `fviz_eig` function from the `factoextra`package to generate the Scree plot.

```{r}
#| label: fig-Pdgfrb_Haralick_PCA
#| include: true
#| warning: false
#| message: false
#| fig-cap: Scree plot for Haralick features
#| fig-width: 4
#| fig-height: 4

# Perform the PCA
Pdgfrb_Haralick_PCA <- PCA(Pdgfrb_Haralick_Variables)

# Plot the Cree plot
fviz_eig(Pdgfrb_Haralick_PCA, addlabels = TRUE)
```

The Scree plot shows that the first component alone comprises 67.5% of the explained variance.

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

# Print the summary
summary(Pdgfrb_Haralick_PCA)

# Print the summary
Pdgfrb_Haralick_PCA$loadings[, 1:2]
```

The results of `summary(Pdgfrb_Haralick_PCA)` show 10 principal components (Comp.1 - Comp.10) that represent the number of variables in our data frame. Each of this components explains the percentage of total variance in the data.

### Biplot of the components attributes

Now, we generate a biplot to witness the similarities between the variables and its impact on each component. All the variables grouped are positively correlated. Also, this plot shows the weight of each variable as distance to the point of origin. Finally, negative correlated variables are displayed in the opposite sides of the Cartesian plane. The color depicts the representation of each variable, with green are higher represented variables.

```{r}
#| label: fig_Pdgfrb_Haralick_Biplot
#| include: true
#| warning: false
#| message: false
#| #| fig-cap: Biplot for Haralick features
#| fig-width: 4
#| fig-height: 4

Pdgfrb_Haralick_Biplot_fig <- fviz_pca_biplot(Pdgfrb_Haralick_PCA,
             geom.ind = "point", # show points only (but not "text")
             col.ind = Pdgfrb_Haralick$Region, # color by groups
             palette = c("blue", "darkgreen", "red"),
             addEllipses = TRUE, # Concentration ellipses
             legend.title = "Region"
             ) +
  labs(tile = "") +
  Plot_theme +
  
  scale_color_manual(name = "Region",   
                     values=c("#0048BA", "red", "darkgreen"),
                     labels=c("Perilesion", "Striatum", "Cortex")) +
                     
  scale_fill_manual(name = "Region", 
                    values=c("#0048BA", "red", "darkgreen"),
                    labels=c("Perilesion", "Striatum", "Cortex")) +
  
  scale_shape(name = "Region", 
                    labels=c("Perilesion", "Striatum", "Cortex")) +
  
  theme(legend.position = "top", legend.direction = "horizontal")
  
  

ggsave(
  plot     = Pdgfrb_Haralick_Biplot_fig, 
  filename = "Plots/Widefield_20x_ROIs_Pdgfrb_Haralick/Widefield_20x_ROIs_Pdgfrb_Haralick_Biplot.png", 
  width    = 12, 
  height   = 12, 
  units    = "cm")

Pdgfrb_Haralick_Biplot_fig
```

Altogether, the PCA analysis suggest that all the variables have a similar loading for the components. This means, we cannot detect more relevant features. Given that we dealing with images that change in intensity, distribution, and density in the course of injury, we judge that the variations in the texture, uniformity, and the complexity of the patterns are features of interest. Based on this, we decided to rely in the following features:

**Contrast:** Given the changes in intensity and distribution, the contrast reveal the amount of local variations in the image. This features will capture the differences as the image changes intensity during injury conditions.

**Entropy:** Entropy provides a measure of the randomness or complexity in the image. Considering that injury results in more complex cellular distributions or more erratic patterns, we expect a higher entropy. Conversely, if cells become more uniform or the patterns become simpler post-injury, entropy will be lower.

**Inverse Difference Moment (IDM) / Homogeneity:** Given our interest in the distribution and density of cells, homogeneity reveals cells cluster or dispersion during the injury course.

### Sypder plot

We create a spyder plot to summarize the select features. For this, we need to prepare the data in tiny format and summarize the variables.

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

# Subset the variables of interest
Pdgfrb_Spyder_Variables <- subset(Pdgfrb_Haralick_Variables, select = c(contrast, inverse_difference_moment, entropy))

# Select columns containing metadata
metadata <- Pdgfrb_Haralick[,c(1:3)]

# Merge the metadata and the variables of interest
Pdgfrb_Haralick_Spyder <- cbind(metadata, Pdgfrb_Spyder_Variables)

# Group by class (cell type), summarize and reshape the date set
Pdgfrb_Haralick_Spyder_Sum <- Pdgfrb_Haralick_Spyder %>%
   group_by(Region) %>%
   summarise(across(4:6, mean)) %>%
  reshape_longer(c("contrast", "inverse_difference_moment",  "entropy"))

# Set variables as factors and rename them
Pdgfrb_Haralick_Spyder_Sum$name <- 
  factor(Pdgfrb_Haralick_Spyder_Sum$name, 
  levels = c("contrast", "inverse_difference_moment",  "entropy"), 
  labels = c("Contrast", "IDM",  "Entropy"))

Pdgfrb_Haralick_Spyder_Sum$Region <- 
  factor(Pdgfrb_Haralick_Spyder_Sum$Region, 
   levels = c("Peri", "Str",  "Ctx"), 
   labels = c("Perilesion", "Striatum", "Cortex"))

# Create the graph
Pdgfrb_Haralick_Spyder_fig <- Pdgfrb_Haralick_Spyder_Sum %>%
  ggplot(aes(
    x = name,
    y = value,
    color =Region,
    group = Region,
    fill = Region
  )) +
  geom_polygon(linewidth = 2, alpha = 0.3) +
  scale_x_discrete(name = "") +
  scale_y_continuous(name = "Value (scaled)") +
  coord_radar() +
  theme_radar() +
  theme(
    legend.position="none",
    legend.direction="horizontal",
    legend.title = element_text(colour="black", face="bold", size=10),
    legend.text = element_text(colour="black", size=8),
    legend.margin=margin(t = 0, unit='cm'))+
  
  guides(fill=guide_legend(nrow=2,byrow=TRUE)) +
  guides(color=guide_legend(nrow=2,byrow=TRUE)) 
   
 
# Save the graph 
ggsave(
  plot     = Pdgfrb_Haralick_Spyder_fig, 
  filename = "Plots/Widefield_20x_ROIs_Pdgfrb_Haralick/Widefield_20x_ROIs_Pdgfrb_Haralick_Spyder.png", 
  width    = 9, 
  height   = 15, 
  units    = "cm")

Pdgfrb_Haralick_Spyder_fig
```

## Statistical modeling

In this case, we perform statistical modeling using a logistic regression in `brms`. The response variable "Region" is modeled as a function of Entropy, Contrast and IDM. In this case, we set the family to `categorical` for multinomial logistic regression model with a logit link function.

$$
Region_{i} = logit ^{-1} (\beta_{0} + \beta_{1} * Entropy_{i} + \beta_{2} * Contrast_{i} + \beta_{3} * IDM_{i} + \epsilon_{i})
$$ This model uses the following weakly-informative priors:

$$
\begin{align}
\beta_{1} \sim Normal(0, 1) \\
\alpha \sim normal(0, 1) \\
\end{align}
$$ \



## Fit the model

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

set.seed(8807)

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

Pdgfrb_Haralick_Mdl1 <- bf(Region ~ entropy + contrast + inverse_difference_moment)

get_prior(Pdgfrb_Haralick_Mdl1, data = Pdgfrb_Haralick_Spyder)

Pdgfrb_Haralick_prior  <- 
  c(prior(normal(0, 1), class = "Intercept"),
    prior(normal(0, 1), class = b))

# Fit model 1
Pdgfrb_Haralick_Fit1 <- 
  brm(
    data    = Pdgfrb_Haralick_Spyder,
    family  = "categorical",
    formula = Pdgfrb_Haralick_Mdl1,
    prior   = Pdgfrb_Haralick_prior,
    chains  = 4,
    cores   = 4,
    warmup  = 2500, 
    iter    = 5000, 
    seed    = 8807,
    control = list(adapt_delta = 0.99, max_treedepth = 15),
    file    = "Models/Widefield_20x_ROIs_Pdgfrb_Haralick/Widefield_20x_ROIs_Pdgfrb_Haralick_Fit1.rds",
    file_refit = "never") 
```

## Model diagnostics

We check the model fitting with `pp_checks`

```{r}
#| label: fig-Pdgfrb_Haralick_Diagnostics
#| include: true
#| warning: false
#| message: false
#| fig-cap: Posterior predictive checks for PDGFRβ+ morphology
#| fig-height: 4
#| fig-width: 5

set.seed(8807)

color_scheme_set("darkgray")

# Model 1
Pdgfrb_Haralick_Mdl1_pp <- 
  brms::pp_check(Pdgfrb_Haralick_Fit1, 
                 ndraws = 100) +
  labs(title = "Posterior predictive checks (model 1)",
  subtitle = "Formula: Region ~ entropy + contrast + inverse_difference_moment") +
  #scale_y_continuous(limits=c(0, 0.15)) +
  Plot_theme  

Pdgfrb_Haralick_Mdl1_pp
```

@fig-Pdgfrb_Haralick_Diagnostics displays a good fit between predictions and data. We can explore further the model using `shinystan`.

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

#launch_shinystan(Pdgfrb_Haralick_Fit1)
```

## Model results

### Visualize conditional effects

We plot the model using `conditional_effects` from the `brm`s package.

```{r}
#| label: fig-Pdgfrb_Haralick_CE
#| include: true
#| warning: false
#| message: false
#| fig-cap: Posterior predictive checks for PDGFRβ+ Haralick features
#| fig-height: 5
#| fig-width: 12

# We create the graph for convex hull
Pdgfrb_Haralick_Class <- 
  conditional_effects(Pdgfrb_Haralick_Fit1, categorical = TRUE)

Pdgfrb_Haralick_Convex <- plot(Pdgfrb_Haralick_Class, 
       plot = FALSE)[[1]]

Pdgfrb_Haralick_Convex_fig <- Pdgfrb_Haralick_Convex  + 
  scale_y_continuous(name = expression ("P(Region)")) +
  scale_x_continuous(name="Entropy (scaled)") +
  Plot_theme +
  theme(legend.position = "top", legend.direction = "horizontal")

ggsave(
  plot     = Pdgfrb_Haralick_Convex_fig, 
  filename = "Plots/Widefield_20x_ROIs_Pdgfrb_Haralick/Widefield_20x_ROIs_Pdgfrb_Haralick_Fit1_Entropy.png", 
  width    = 12, 
  height   = 10, 
  units    = "cm")


# We create the graph for Branch length
Pdgfrb_Haralick_Branch <- plot(Pdgfrb_Haralick_Class, 
       plot = FALSE)[[2]]

Pdgfrb_Haralick_Branch_fig <- Pdgfrb_Haralick_Branch  + 
  scale_y_continuous(name = expression ("P(Region)")) +
  scale_x_continuous(name="Contrast (scaled)") +
  Plot_theme +
  theme(legend.position = "top", legend.direction = "horizontal")

ggsave(
  plot     = Pdgfrb_Haralick_Branch_fig, 
  filename = "Plots/Widefield_20x_ROIs_Pdgfrb_Haralick/Widefield_20x_ROIs_Pdgfrb_Haralick_Fit1_Contrast.png", 
  width    = 12, 
  height   = 10, 
  units    = "cm")


# We create the graph for Intensity
Pdgfrb_Haralick_Inten <- plot(Pdgfrb_Haralick_Class, 
       plot = FALSE)[[3]]

Pdgfrb_Haralick_Inten_fig <- Pdgfrb_Haralick_Inten  + 
  scale_y_continuous(name = expression ("Pr(Region)")) +
  scale_x_continuous(name="IDM (scaled)") +
  Plot_theme +
  theme(legend.position = "top", legend.direction = "horizontal")

ggsave(
  plot     = Pdgfrb_Haralick_Inten_fig, 
  filename = "Plots/Widefield_20x_ROIs_Pdgfrb_Haralick/Widefield_20x_ROIs_Pdgfrb_Haralick_Fit1_IDM.png", 
  width    = 12, 
  height   = 10, 
  units    = "cm")

Pdgfrb_Haralick_Convex_fig | Pdgfrb_Haralick_Branch_fig | Pdgfrb_Haralick_Inten_fig
```

### Posterior summary

Next, we plot the posterior summary for both models using the `describe_posterior` function:

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

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

modelsummary(Pdgfrb_Haralick_Fit1, 
             shape = term ~ model + statistic,
             centrali2ty = "mean", 
             title = "Probability of Region for specified Haralick features following MCAO",
             statistic = "conf.int",
             gof_omit = 'ELPD|ELDP s.e|LOOIC|LOOIC s.e|WAIC|RMSE',
             output = "Tables/html/Widefield_20x_ROIs_Pdgfrb_Haralick_Fit1_Table.html",
             )

Pdgfrb_Haralick_Fit1_Table <- modelsummary(Pdgfrb_Haralick_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_Haralick_Fit1_Table, 
            filename = "Tables/tex/Widefield_20x_ROIs_Pdgfrb_Haralick_Fit1_Table.tex")

```

# References

::: {#refs}
:::

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