sara_stats.Rmd

---
title: "Presentation Station Indignation"
author: "Sara Herbst"
date: "2025-09-02"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Cool Stats

```{r deps}
library(lme4)
library(boot) 
library(emmeans)
library(pbkrtest)
library(lmerTest)
library(ggplot2)
library(dplyr)
library(tidyr)
```

```{r read}

v1_df <- read.csv(file = "df_stats_v1_filtered_scaled.csv", sep = ";") |>
  mutate(
    total_count = count / scaled.count
  )
head(v1_df)
```

```{r analysis}
m <- lmer(log(count) ~ age * region + log(total_count) + (1|mouse), data=v1_df)
summary(m)
```

```{r emmeans}
# Estimated marginal means on log scale
emm <- emmeans(m, ~ age * region, type = "response")

# Back-transform to counts
summary(emm)
```

```{r plot}
age_levels <- levels(v1_df$age)
n_age <- length(age_levels)
dodge_width <- 0.4

emm_df <- emm |>
  as.data.frame() |>  
  group_by(region) |>
  mutate(
    x_num = as.numeric(region),
    age_idx = as.numeric(age),
    dodge_offset = (age_idx - (n_age + 1)/2) * dodge_width,
    x_pos = x_num + dodge_offset,
  )

# Needle half-width
needle_half_width <- 0.175

# Vertical expected counts plot
ggplot(emm_df) +
  # CI ribbons
  geom_linerange(aes(ymin = lower.CL, ymax = upper.CL,
                     x = x_pos, color = age),
                 size = 15, alpha = 0.3) +
  # Needles = point estimates
  geom_segment(aes(x = x_pos - needle_half_width,
                   xend = x_pos + needle_half_width,
                   y = response, yend = response,
                   color = age),
               size = 0.8) +
  scale_x_continuous(breaks = (1:length(levels(emm_df$region))) + 0.5,
                     labels = levels(emm_df$region)) +
  labs(title = "Cell Counts", y = "Expected count", x = NULL) +
  theme_minimal() +
  scale_color_manual(values = c("skyblue", "salmon")) +  
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, size = 16),
    axis.text.y = element_text(size = 16),
    axis.title = element_text(size = 16),
    plot.title = element_text(size = 20)
  )
```

```{r save1}
ggsave(
  "count.png",
  plot = last_plot(),
  width = 10,
  height = 8,
  dpi = 300
)
```

```{r fold}
fold_df <- contrast(
  emm,
  method = "revpairwise",
  by = "region",
  type = "response",
  adjust = "none"
) %>%
  summary(infer = TRUE, type = "response") %>%
  as.data.frame() %>%
  rename(
    fold_change = ratio,
    fold_lower  = lower.CL,
    fold_upper  = upper.CL
  ) %>%
  mutate(x = as.numeric(factor(region)))

# ----------------------------
# Plot fold-change
# ----------------------------
needle_half_width <- 0.26
ggplot(fold_df, aes(x = x)) +
  # Ribbon = 95% CI
  geom_linerange(aes(ymin = fold_lower, ymax = fold_upper),
                 size = 16, alpha = 0.3, color = "skyblue") +
  # Needle = point estimate
  geom_segment(aes(x = x - needle_half_width,
                   xend = x + needle_half_width,
                   y = fold_change, yend = fold_change),
               size = 0.8, color = "skyblue") +
  # Reference line at 1
  geom_hline(yintercept = 1, linetype = "dashed", color = "gray50") +
  # x-axis labels
  scale_x_continuous(breaks = fold_df$x, labels = fold_df$region) +
  labs(title = "Fold Change", y = "Fold-change (p7-13 / p1-7)", x = NULL) +
  theme_minimal() +  
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, size = 16),
    axis.text.y = element_text(size = 16),
    axis.title = element_text(size = 16),
    plot.title = element_text(size = 20)
  )
```

```{r save2}
ggsave(
  "folds.png",
  plot = last_plot(),
  width = 6,
  height = 8,
  dpi = 300
)
```