Widefield_10x_Ipsilateral_KO_Fish_Pdgfrb.qmd

---
title-block-banner: true
title: "Analysis of PDGFRβ mRNA (FISH) in the injured brain following cerebral ischemia"
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β
  - Flourescence in situ hybridization (FISH)
  - 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 PDGFRβ expression (by FISH) in the injured brain following cerebral ischemia. Previous results using PDGFRβ(Td-tomato from transgenic animals) suggest that KLF4 depletion reduces PDGFRβ expression. To confirm those results we evaluated PDGFRβ at mRNA levels. 

**Parent dataset:** PDGFRβ(FISH) stained ischemic hemispheres imaged at 10x (with stitching). We took regions of interest (ROIs) from three different areas in the injure cortex/striatum to quantify PDGFRβ expression. Samples are grouped at 0 (Sham), and 14 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_10x_ROIs_KO_Fish_Pdgfrb.zip`. Please note that we excluded mice `KL31` because of an inefficient staining (the section dried during the staining procedure).

**Working dataset**: The `Data_Raw/Widefield_10x_ROIs_KO_Fish_Pdgfrb/Widefield_10x_ROIs_KO_Fish_Pdgfrb_Pdgfrb_Objects.csv dataframes containing the raw output from CellProfiller [@stirling2021]. The CellProfiler pipeline is available at https://osf.io/yndua.

We perform scientific inference based on the number of cells, intensity and area. We also perform PPA and TDA to evaluate topological features of cells.

# 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", "reshape2", "spatstat", "scales", "tidybayes", "tidyverse", "viridis"))


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

# Analysis of the number of PDGFRβ cells

## Load and handle the datasets

We load the `Data_Raw/Widefield_10x_ROIs_KO_Fish_Pdgfrb/Widefield_10x_ROIs_KO_Fish_Pdgfrb_Image.csv` datasets containing raw cell counts from cellProfiler.

```{r}
#| label: tbl-Pdgfrb_Counts_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_df <- read.csv(file = "Data_Raw/Widefield_10x_ROIs_KO_Fish_Pdgfrb/Widefield_10x_ROIs_KO_Fish_Pdgfrb_Image.csv", header = TRUE)

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

We are interested in the `FileName_Pdgfrb_Raw`column containing the identification data for the images and the `Count_Pdgfrb_Objects` indicating the number of PDGFRβ+ cells. Next, we subset the data frame to select the columns of interest and give them meaningful names.

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

## We subset the relevant columns (cell number)
Pdgfrb_Counts <- subset(Pdgfrb_df, select = c("FileName_Pdgfrb_Raw", "Count_Pdgfrb_Objects"))

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

Pdgfrb_Counts <- subset(Pdgfrb_Counts, select = -c(FileName_Pdgfrb_Raw))

## We Rename the relevant columns 
colnames(Pdgfrb_Counts) <- c("Counts", "MouseID", "DPI", "Genotype", "Section")

## We set the factors
Pdgfrb_Counts$DPI <- factor(Pdgfrb_Counts$DPI, levels = c("0D", "14D"))
Pdgfrb_Counts$Genotype <- factor(Pdgfrb_Counts$Genotype, levels = c("WT", "KO"))
Pdgfrb_Counts$Section <- factor(Pdgfrb_Counts$Section, levels = c("S1", "S2", "S3"))

write.csv(Pdgfrb_Counts, "Data_Processed/Widefield_10x_ROIs_KO_Fish_Pdgfrb/Widefield_10x_ROIs_KO_Fish_Pdgfrb_Counts.csv", row.names = FALSE)

```

With the data handled, we perform exploratory data visualization.

## Exploratory data visualization

We plot the data to see its distribution and have a first hint for modeling strategies.

```{r}
#| label: fig-Pdgfb_Counts_Expl
#| include: true
#| warning: false
#| message: false
#| results: false
#| fig-cap: Exploratory data visualization for the number of PDGFRβ+ cells
#| fig-width: 9
#| fig-height: 4

set.seed(8807)

# Boxplot
##################

Pdgfrb_Counts_box <- 
  ggplot(
    data  = Pdgfrb_Counts, 
    aes(x = DPI, 
        y = Counts,
        color = Genotype)
    ) +
  geom_boxplot() +

scale_y_continuous(name = "Number of PDGFRβ+ cells") +
scale_x_discrete(
  name   ="Days post-ischemia (DPI) ",
  breaks = c("0D", "14D")
  ) +
Plot_theme

#Plot the result
Pdgfrb_Counts_box
```

@fig-Pdgfb_Counts_Expl shows that cerebral ischemia comprises an increase in the number of PDGFRβ cells. The plots at 14D show a tiny difference between the means of the two genotypes.

## Statistical modeling for the number of PDGFRβ+ cells

Given that we took three different photos per animal, we'll fit a multilevel model to capture the both the population-level effects of Genotype and the group variability due to individual differences among the mice. We fit the following model:

 - **Pdgfrb_Counts_Mdl1:** Considers `Genotype` as a predictor and `MouseID` as random effects. The model takes the following notation:
 
 $$
 Counts_i = \beta_0 + \beta_1 \times Genotype_i + u_{MouseID[i]} + \epsilon_i
 $$
Where $β0$ is the intercept, $𝛽1$, represents the effect of Genotype on Counts. $𝑢MouseID$ is the group-level effect associated with the $𝑖$i-th mouse, and $𝜖$ denotes the residual error terms.
 .
By including (1 | MouseID), we accounts for the fact that observations from the same mouse are likely to be more similar to each other than to observations from different mice

### Fit the models

We employ `brms` to fit the model.

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

Pdgfrb_Counts_14D <- Pdgfrb_Counts[Pdgfrb_Counts$DPI == "14D",]

# Model 1: Genotype as a linear predictor
###########################################
Pdgfrb_Counts_Mdl1 <- bf(Counts ~ Genotype + (1 | MouseID))

get_prior(Pdgfrb_Counts_Mdl1, Pdgfrb_Counts_14D , family = negbinomial())

# Fit model 1
Pdgfrb_Counts_Fit1 <- 
  brm(
    family  = negbinomial(),
    data    = Pdgfrb_Counts_14D,
    formula = Pdgfrb_Counts_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_KO_Fish_Pdgfrb/Widefield_10x_ROIs_KO_Fish_Pdgfrb_Counts_Fit1.rds",
    file_refit = "never")
```

### Model diagnostics

To evaluate sample predictions, we perform the model diagnostics 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_Counts_Diagnostics
#| include: true
#| warning: false
#| message: false
#| results: false
#| fig-cap: Model diagnostics for the number of PDGFRβ+ cells
#| fig-height: 4
#| fig-width: 5

set.seed(8807)

Pdgfrb_Counts_Fit1_pp <- 
  brms::pp_check(Pdgfrb_Counts_Fit1, 
                 ndraws = 100) +
  labs(title = "Posterior predictive checks",
  subtitle = "Formula: Counts ~ Genotype + (1 | MouseID)") +
  Plot_theme  
  
Pdgfrb_Counts_Fit1_pp
```

We see that our model captures the trend of the observations.

## Model results

### Visualizaing the posterior

After validating the model, we plot the full posterior distributions for the population and group (Mouse) level effects.:

```{r}
#| label: fig-Pdgfrb_Counts_Posterior
#| include: true
#| warning: false
#| message: false
#| results: false
#| fig-cap: Posterior distribution for the number of of PDGFRβ+ cells
#| fig-width: 10
#| fig-height: 4

# For population-Level effects:

set.seed(8807)

Pdgfrb_Count_Draws <- Pdgfrb_Counts_14D %>%
  group_by(Genotype, MouseID) %>% 
  add_predicted_draws(Pdgfrb_Counts_Fit1, re_formula = NA) 

# Plotting
Pdgfrb_Count_fig <- ggplot(Pdgfrb_Count_Draws, aes(x = Genotype, y = .prediction)) +
  stat_halfeye() + # Visualize the posterior predictive distribution
  geom_jitter(data = Pdgfrb_Counts_14D, aes(x = Genotype, y = Counts),
              width = 0.15, alpha = 0.6, size = 2) + # Actual measurements
  labs(x = "Genotype") + 
  scale_y_continuous(name = "Number of PDGFRβ+ cells",
                     limits = c(0, 400),
                      breaks = seq(0, 400, 100) ) +
  Plot_theme +
  coord_flip()

ggsave(
  plot     = Pdgfrb_Count_fig, 
  filename = "Plots/Widefield_10x_ROIs_KO_Fish_Pdgfrb/Widefield_10x_ROIs_KO_Fish_Pdgfrb_Counts.png", 
  width    = 9, 
  height   = 9, 
  units    = "cm")

# For Group-Level effects

# Extracting group-level effects
posterior_samples <- Pdgfrb_Counts_Fit1 %>% 
  spread_draws(r_MouseID[MouseID,])  # Extract random effects for AnimalID

# Plotting the group-level effects
Pdgfrb_Counts_group <- ggplot(posterior_samples, 
                              aes(x = reorder(MouseID,r_MouseID), 
                                  y = r_MouseID)) +
  stat_halfeye(fill = "lightgray", color = "black") +
  labs(title = "Mouse-Level effects",
       x = "MouseID",
       y = "Deviation from mean") +
  coord_flip()+
    geom_hline(yintercept=0, size=1, linetype=2)+
  Plot_theme

ggsave(
  plot     = Pdgfrb_Counts_group, 
  filename = "Plots/Widefield_10x_ROIs_KO_Fish_Pdgfrb/Widefield_10x_ROIs_KO_Fish_Pdgfrb_Group.png", 
  width    = 9, 
  height   = 12, 
  units    = "cm")


Pdgfrb_Count_fig + Pdgfrb_Counts_group

```

@fig-Pdgfrb_Counts_Posterior shows the full posterior distribution for population and group (MouseID) level effects. We see that WT and KO mice does not differ meaningfully in the number of PDGFR-B cells. 


### Posterior summary

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

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

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

modelsummary(Pdgfrb_Counts_Fit1, 
             shape = term ~ model + statistic,
             centrality = "mean", 
             title = "Number of PDGFR-B+ nuclei in defined ROIs of the ipsilateral hemisphere",
             statistic = "conf.int",
             gof_omit = 'ELPD|ELDP s.e|LOOIC|LOOIC s.e|WAIC|RMSE',
             output = "Tables/html/Widefield_10x_ROIs_KO_Fish_Pdgfrb_Counts_Fit1_Table.html",
             )

Pdgfrb_Counts_Fit1_Table <- modelsummary(Pdgfrb_Counts_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_Counts_Fit1_Table, filename = "Tables/tex/Widefield_10x_ROIs_KO_Fish_Pdgfrb_Counts_Fit1_Table.tex")
```
The table displays the regression coefficients o the log scale.


# Analysis of PDGFR-B area and intensity

## Load and handle the dataset 

Now, we handle the `Widefield_10x_ROIs_KO_Fish_Pdgfrb_Pdgfrb_Objects.csv` to analyze the area and intensity of the PDGFR-B+ cells

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

Pdgfrb_AreaInt_df <- read.csv("Data_Raw/Widefield_10x_ROIs_KO_Fish_Pdgfrb/Widefield_10x_ROIs_KO_Fish_Pdgfrb_Pdgfrb_Objects.csv", header = TRUE)

gt::gt(Pdgfrb_AreaInt_df[1:10,])
```
We are interested in the `FileName_Pdgfrb_Raw`column containing the identification data for the images, the `AreaShape_Area` and `Intensity_IntegratedIntensity_Pdgfrb_Raw` indicating the area and intensity per cell, respectively. We also add `Location_Center_X` and `Location_Center_X` to get the cell's coordinates. Next, we subset the data frame to select the columns of interest and give them meaningful names.

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

## We subset the relevant columns (cell number)
Pdgfrb_AreaInt <- subset(Pdgfrb_AreaInt_df, select = c("FileName_Pdgfrb_Raw", "AreaShape_Area", "Intensity_IntegratedIntensity_Pdgfrb_Raw", "Location_Center_X", "Location_Center_Y"))

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

Pdgfrb_AreaInt <- subset(Pdgfrb_AreaInt, select = -c(FileName_Pdgfrb_Raw))

## We Rename the relevant columns 
colnames(Pdgfrb_AreaInt) <- c("Area", "Intensity","x-coord", "y-coord", "MouseID", "DPI", "Genotype", "Section")

## We set the factors
Pdgfrb_AreaInt$DPI <- factor(Pdgfrb_AreaInt$DPI, levels = c("0D", "14D"))
Pdgfrb_AreaInt$Genotype <- factor(Pdgfrb_AreaInt$Genotype, levels = c("WT", "KO"))
Pdgfrb_AreaInt$Section <- factor(Pdgfrb_AreaInt$Section, levels = c("S1", "S2", "S3"))

write.csv(Pdgfrb_AreaInt, "Data_Processed/Widefield_10x_ROIs_KO_Fish_Pdgfrb/Widefield_10x_ROIs_KO_Fish_Pdgfrb_AreaInt.csv", row.names = FALSE)

```

## Exploratory data visualization

Next, we visualize the raw data to guide the statistical modeling. 

```{r}
#| label: fig-Pdgfrb_AreaInt_Expl
#| include: true
#| warning: false
#| message: false
#| results: false
#| fig-cap: Exploratory data visualization for PDGFR-B area and intensity
#| fig-width: 9
#| fig-height: 4

set.seed(8807)

# Box plot for area
##################

Pdgfrb_Area_Box <- 
  ggplot(
    data  = Pdgfrb_AreaInt, 
    aes(x = DPI, 
        y = Area,
        color = Genotype)
    ) +
  geom_boxplot()+

Plot_theme

# Box plot for intensity
##################

Pdgfrb_Intensity_Box <- 
  ggplot(
    data  = Pdgfrb_AreaInt, 
    aes(x = DPI, 
        y = Intensity,
        color = Genotype)
    ) +
  geom_boxplot()+

Plot_theme

#Plot the result
Pdgfrb_Area_Box + Pdgfrb_Intensity_Box


Pdgfrb_Intensity_scr <- 
  ggplot(
    data  = Pdgfrb_AreaInt, 
    aes(x = Area, 
        y = Intensity,
        color = Genotype)
    ) +
  geom_point()

Pdgfrb_Intensity_scr 
```

@fig-Pdgfrb_AreaInt_Expl shows that both area and intensity of the cells change after MCAO. On the other hand, while the area seems similar for both genotypes at 14 DPI, there is an apparent change in the intensity of the cells that is worth modeling. 

## Statistical modeling for PDGFR-B area/intensity

To analyze the changes in PDGFRβ expression considering area and intensity, we utilized a mixed-effects model (family = student) with the interactions of Area and Genotype as predictors. We use MouseID as random term as previously specified. The model specification is as follows:

$$
Intensity_i = \beta_0 + \beta_1 \times Area_i + \beta_2 \times Genotype_i + \beta_3 \times (Area_i \times Genotype_i) + u_{MouseID[i]} + \epsilon_i
$$
Here, $\beta_0$ is the intercept (area and Genotype at reference levels); $\beta_1$ is the effect of area on intensity; $\beta_2$ is the effect of Genotype on Intensity; $\beta_3$ represents the interaction effect between Area and Genotype;$u_{MouseID[i]}$ accounts for the random effect associated with each mouse; $\epsilon_i$ is the error term.   

### Fit the model

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

# Model 1

Pdgfrb_AreaInt_14D <- Pdgfrb_AreaInt[Pdgfrb_AreaInt$DPI == "14D",]

Pdgfrb_AreaInt_Mdl1 <- bf(Intensity ~ Area * Genotype + (1 | MouseID))

get_prior(Pdgfrb_AreaInt_Mdl1, Pdgfrb_AreaInt_14D, family = student)

# Fit model 1
Pdgfrb_AreaInt_Fit1 <- 
  brm(
    family  = student,
    data    = Pdgfrb_AreaInt_14D,
    formula = Pdgfrb_AreaInt_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_KO_Fish_Pdgfrb/Widefield_10x_ROIs_KO_Fish_Pdgfrb_AreaInt_Fit1.rds",
    file_refit = "never")
```

## Model diagnostics

Next, we perform model diagnostics:

```{r}
#| label: fig-Pdgfrb_Diagnostics
#| include: true
#| warning: false
#| message: false
#| results: false
#| fig-cap: Model diagnostics for PDGFR-B expression in individual cells
#| fig-height: 5
#| fig-width: 10

set.seed(8807)

Pdgfrb_AreaInt_Fit1_pp <- 
  brms::pp_check(Pdgfrb_AreaInt_Fit1, 
                 ndraws = 100) +
  labs(title = "Posterior predictive checks",
  subtitle = "Formula: Intensity ~ Area * Genotype + (1 | MouseID)") +
  scale_x_continuous(limits=c(0, 20)) +
  Plot_theme  

Pdgfrb_AreaInt_Fit1_pp 
```

We see that this model does not align with the observed data. Therefore, we fit an additional model with scaled predictors.

### Fit the new model

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

# Model 2

Pdgfrb_AreaInt_14D_Scaled <- data.frame(scale(Pdgfrb_AreaInt_14D[,c("Intensity","Area")], scale=T, center=T), MouseID=Pdgfrb_AreaInt_14D$MouseID, DPI=Pdgfrb_AreaInt_14D$DPI, Genotype = Pdgfrb_AreaInt_14D$Genotype, Section = Pdgfrb_AreaInt_14D$Section)


Pdgfrb_AreaInt_Mdl2 <- bf(Intensity ~ Area * Genotype + (1 | MouseID))

get_prior(Pdgfrb_AreaInt_Mdl1, Pdgfrb_AreaInt_14D_Scaled, family = student)

Pdgfrb_AreaInt_Prior <- 
  c(prior(normal(0,5), class = b),
    prior(normal(0,5), class = Intercept, lb= 0))

# Fit model 2
Pdgfrb_AreaInt_Fit2 <- 
  brm(
    family  = student,
    data    = Pdgfrb_AreaInt_14D_Scaled,
    formula = Pdgfrb_AreaInt_Mdl2,
    prior   = Pdgfrb_AreaInt_Prior,
    chains  = 4,
    cores   = 4,
    warmup  = 2500, 
    iter    = 5000, 
    seed    = 8807,
    control = list(adapt_delta = 0.99, max_treedepth = 15),
    file    = "Models/Widefield_10x_ROIs_KO_Fish_Pdgfrb/Widefield_10x_ROIs_KO_Fish_Pdgfrb_AreaInt_Fit2.rds",
    file_refit = "never")
```
```{r}
#| label: fig-Pdgfrb_Diagnostics2
#| include: true
#| warning: false
#| message: false
#| results: false
#| fig-cap: Model diagnostics for PDGFR-B expression in individual cells
#| fig-height: 5
#| fig-width: 10

set.seed(8807)

Pdgfrb_AreaInt_Fit1_pp <- 
  brms::pp_check(Pdgfrb_AreaInt_Fit2, 
                 ndraws = 100) +
  labs(title = "Posterior predictive checks",
  subtitle = "Formula: Intensity ~ Area * Genotype + (1 | MouseID)") +
  scale_x_continuous(limits=c(0, 20)) +
  Plot_theme  

Pdgfrb_AreaInt_Fit1_pp 
```

We see that the model captures the overall trend but deviates slightly. We verify further the model:

```{r}
#| label: fig-Pdgfrb_Diagnostics3
#| include: true
#| warning: false
#| message: false
#| results: false
#| fig-cap: Model diagnostics for PDGFR-B expression in individual cells
#| fig-height: 5
#| fig-width: 10

# We obtain the posterior
Pdgfrb_AreaInt_Posterior <- posterior_predict (Pdgfrb_AreaInt_Fit2, ndraws = 1000)

# We plot the mean per Genotype
Pdgfrb_mean_Prod <- ppc_stat_grouped(y=Pdgfrb_AreaInt_14D_Scaled$Intensity,yrep=Pdgfrb_AreaInt_Posterior, group = Pdgfrb_AreaInt_14D_Scaled$Genotype, stat = "mean", binwidth = 0.01)+
xlab("Mean intensity")+
#scale_x_continuous(limits=c(10, 40)) +
Plot_theme 

Pdgfrb_mean_Prod
```
We see this models captures better the mean for WT than for KO.

## Model results

After model diagnostics, we plot the full posterior distributions:

### Visualizing the results

```{r}
#| label: fig-Pdgfrb_AreaInt_Posterior
#| include: true
#| warning: false
#| message: false
#| results: false
#| fig-cap: Posterior distribution for PDGFR-B expression
#| fig-width: 5
#| fig-height: 4

set.seed(8807)

Pdgfrb_AreaInt_DPI <- 
  conditional_effects(Pdgfrb_AreaInt_Fit2)

Pdgfrb_AreaInt_DPI <- plot(Pdgfrb_AreaInt_DPI, 
       plot = FALSE, points = TRUE, point_args = list(size =0.1))[[3]]

Pdgfrb_AreaInt_fig <- Pdgfrb_AreaInt_DPI + 
  scale_y_continuous(name = expression ("PDGFRβ intensity"),
                     limits = c(0, 10)) +
  scale_x_continuous(name="Nuclei area") +
  
  scale_color_manual(name = "",
                       values = c("#4DB6D0", "#D9717D"),
                       labels = c("Wildtype", expression ("PDGFRβ"^KLF4_KO))) +
 scale_fill_manual(name = "",
                     values = c("#4DB6D0", "#D9717D"),
   labels = c("Wildtype", expression ("PDGFRβ"^KLF4_KO))) +
  
  Plot_theme +
  theme(legend.position = "top",
        legend.direction = "horizontal")

  ggsave(
  plot     = Pdgfrb_AreaInt_fig, 
  filename = "Plots/Widefield_10x_ROIs_KO_Fish_Pdgfrb/Widefield_10x_ROIs_Fish_Pdgfrb_IntArea.png",   width    = 9, 
  height   = 9, 
  units    = "cm")

Pdgfrb_AreaInt_fig 
```

@fig-Pdgfrb_AreaInt_Posterior shows the posterior distribution of PDGFR-B expression by cell, accouting for the nuclei size.


### Posterior summary

Next, we plot the posterior summary:
```{r}
#| label: Pdgfrb_AreaInt_DescribePosterior
#| include: true
#| warning: false
#| message: false
#| results: false
#| cache: true

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

modelsummary(Pdgfrb_AreaInt_Fit2, 
             shape = term ~ model + statistic,
             centrality = "mean", 
             title = "PDGFR-B+ LABELING INTENSITY in defined ROIs of the ipsilateral hemisphere",
             statistic = "conf.int",
             gof_omit = 'ELPD|ELDP s.e|LOOIC|LOOIC s.e|WAIC|RMSE',
             output = "Tables/html/Widefield_10x_ROIs_KO_Fish_Pdgfrb_AreaInt_Fit1_Table.html",
             )

Pdgfrb_AreaInt_Fit2_Table <- modelsummary(Pdgfrb_AreaInt_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_AreaInt_Fit2_Table, filename = "Tables/tex/Widefield_10x_ROIs_KO_Fish_Pdgfrb_AreaInt_Fit1_Table.tex")
```
The table displays the regression coefficients without a baseline (intercept value).It also shows he coefficients for individual animals included as group-level effects.

# References

::: {#refs}
:::

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