shiverer.Rmd

```{r include=FALSE}
library(tidyverse)
library(readxl)
library(lme4)
library(cowplot)
library(philentropy)
library(dabestr)
library(lqmm)
```

# Goals
We have a dataset contain axon diameters of neurons from the optic nerve of control and mutant zebrafish. We'd like to know if the mean axon diameter, or the distribution of axon diameters, differs between groups. We want to implement tests that take account of within and between animal variance.

Here, we focus on the shiverer data set.

# Load and format data
```{r message = FALSE}
df <- read_excel("Shiverer Axon Diameter.xlsx", range = "A2:J293") %>%
    pivot_longer(cols = 1:10) %>%
    drop_na() %>%
    mutate(group = as_factor(ifelse(name %in% c("149-4555", "149-4567", "149-4591", "149-4593"), "control", "mutant")))
```

# Plot the data
Focus here on plot of individual mice, colour coded by group.
```{r}
(plot_by_id <- ggplot(data = df, aes(name, value)) +
     geom_violin(aes(colour = group), draw_quantiles = c(0.25, 0.5, 0.75)) +
     theme_cowplot(font_size = 14) +
     theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
    labs(x = 'Animal', y = 'Diameter (\u00B5 m)', colour = "Group"))

ggsave('Plots/violins_shiv.jpeg', plot_by_id)
```

```{r}
(plot_by_id <- ggplot(data = df, aes(name, log(value))) +
     geom_violin(aes(colour = group), draw_quantiles = c(0.25, 0.5, 0.75)) +
     theme_cowplot(font_size = 14) +
     theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
    labs(x = 'Animal', y = 'Diameter (\u00B5 m)', colour = "Group"))

ggsave('Plots/violins_log_shiv.jpeg', plot_by_id)
```


# Tests for differences in means
```{r}
mm_t <- lmer(value ~ group + (1 | name), data = df)
mm_t_null <- lmer(value ~ (1 | name), data = df)
summary(mm_t)
anova(mm_t, mm_t_null)
```

# Same test but with log transformed data
```{r}
mm_t <- lmer(log(value) ~ group + (1 | name), data = df)
mm_t_null <- lmer(log(value) ~ (1 | name), data = df)
summary(mm_t)
anova(mm_t, mm_t_null)
```



# Generate histograms for all animals
```{r make histograms}
### Make histograms for each animal
### Use log transformed data
names = unique(df$name)
### If submean = 1 then will substract means before making histograms 
histfun <- function(name, df, submean = 0) {
    sub <- ifelse(submean == 1, mean(df$value[df$name==name]), 0)
    hist(log(df$value[df$name==name]-sub), seq(-2,2,0.1), plot = FALSE)
}
hists <- sapply(names, histfun, df, submean=0)

### Convert results to tibble for use with tidyverse functions
rns <- rownames(hists)
hists_tib <- as_tibble(hists) %>%
    rownames_to_column(var = "rowname") %>%
    pivot_longer(-rowname, names_to = "column", values_to = "value") %>%
    pivot_wider(names_from = rowname, values_from = value)
colnames(hists_tib) <- c("animal", rns)

    
ggconvfun <- function(mids, density) {
    gghist <- cbind(mids, density)
    ggplot(gghist, aes(mids, density)) +
        geom_col() +
        theme_cowplot()
}

gghists <- map2(hists_tib$mids, hists_tib$density, ggconvfun)

plot_grid(plotlist = gghists)
```


# Evaluate distributions
Use package lqmm (https://www.jstatsoft.org/article/view/v057i13) to evaluate potential differences in the distributions.
Use Nelder-Mead optimization (derivative-free method) as gives fits with much narrower bounds on the intercept.
```{r}
fit.lqmm <- lqmm(fixed = value ~ group, random = ~1, group = name, tau = c(0.25,0.5, 0.75), nK = 7, type = "normal", data = df, control = c(method = "df", LP_max_iter = 5000))
summary(fit.lqmm, R = 500, seed = 2)
```

```{r}
fit.log.lqmm <- lqmm(fixed = log(value) ~ group, random = ~1, group = name, tau = c(0.25,0.5, 0.75), nK = 7, type = "normal", data = df, control = c(method = "df", LP_max_iter = 5000))
summary(fit.log.lqmm, R = 500, seed = 2)
```