#Updated 7.19.23

# STATISTICS
#Variable: AR Grains Above Back Per Cell
#Linear mixed model with sample nested within subjectID and treated as a random factor
#Fixed factors are Motor Nucleus (lumbar, brachial) and Species (Sp, Rp, Xl) and their interaction
library (lme4)
library (lmerTest)
library (performance)
require (mosaic)

model1 = lmer (GrAbBackPerCell ~ MotorNuc*Species + (1|SubjectID/Sample), data =ARData)
summary (model1)
anova (model1)

#Post hoc pairwise comparisons using Tukey method
library (emmeans)
emmeans (model1, specs=pairwise ~ MotorNuc:Species)



#Load data on AR Ratio of lumbar:brachial within individuals
#ANOVA to test for species differences in AR Ratio L:B

favstats (ARRatio ~ Species, data=ARRatioData)

onewayanova <- aov(ARRatio ~ Species, data = ARRatioData)
summary (onewayanova)



#Variable: MN Count
#Linear mixed model with sample nested within subjectID and treated as a random factor
#Fixed factors are Motor Nucleus (lumbar, brachial) and Species (Sp, Rp, Xl) and their interaction
library (lme4)
library (lmerTest)
library (performance)
require (mosaic)

model2 = lmer (MNCount ~ MotorNuc*Species + (1|SubjectID/Sample), data =ARData)
summary (model2)
anova (model2)


#Post hoc pairwise comparisons using Tukey method
library (emmeans)
emmeans (model2, specs=pairwise ~ MotorNuc:Species)



#Variable: IN Count
#Linear mixed model with sample nested within subjectID and treated as a random factor
#Fixed factors are Motor Nucleus (lumbar, brachial) and Species (Sp, Rp, Xl) and their interaction
library (lme4)
library (lmerTest)
library (performance)
require (mosaic)

model3 = lmer (INCount ~ MotorNuc*Species + (1|SubjectID/Sample), data =ARData)
summary (model3)
anova (model3)

#Post hoc pairwise comparisons using Tukey method
library (emmeans)
emmeans (model3, specs=pairwise ~ MotorNuc:Species)



#FIGURES
#Run this code first: Load packages
library (ggplot2)
library (dplyr)
library (wesanderson)


# Grains Above Back per Cell, Grouped By Species, Boxplots with jittered points
ggplot (ARData, aes (x=factor (Species, level=c('Sp', 'Rp', 'XL')), y=GrAbBackPerCell, fill=MotorNuc)) +
  geom_boxplot (alpha=0.5) +
  geom_point(position=position_jitterdodge(jitter.width=0.2), alpha=0.5) +
  scale_fill_grey () +
  theme_classic ()

#Same graph as above with individuals denoted by different shapes

ggplot (ARData, aes (x=factor (Species, level=c('Sp', 'Rp', 'XL')), y=GrAbBackPerCell), fill=MotorNuc) +
  geom_boxplot (aes (fill=MotorNuc), position=position_dodge(width=0.85),
  alpha=0.5, outlier.shape=NA) +
  geom_point (aes (shape=SubjectID, group=MotorNuc),
  position=position_jitterdodge(dodge.width=0.85, jitter.width=0.2),
  alpha=0.5) +
  scale_shape_manual (values = c(Sp24=15,Sp29=17,Sp30=19,Rp9=15,Rp11=17,Rp10=19,XL10=15,XL9=17,XL12=19)) +
  theme_classic () +
  scale_fill_grey () +
  guides (shape=FALSE)
  
  



#Plot of Interneuron, Motorneuron Counts
#Run this code first: Order the levels of factor Species the way that I want them in dataframe ARData
ARData$Species <- factor (ARData$Species, levels=c("Sp", "Rp", "XL"), ordered=TRUE)

#Plot Interneuron Counts
#Grouped by Motor Nucleus
ggplot (ARData, aes (x=MotorNuc, y=INCount, fill=Species)) +
  geom_boxplot (alpha=0.5) +
  geom_point (position=position_jitterdodge(jitter.width=0.15), alpha=0.5) +
  scale_fill_manual (values=wes_palette (n=3, name="FantasticFox1")) +
  theme_classic ()

#Same graph as above with individuals denoted by different shapes

ggplot (ARData, aes (x=MotorNuc, y=INCount, fill=Species)) +
  geom_boxplot (aes (fill=Species), position=position_dodge(width=0.85),
  alpha=0.5, outlier.shape=NA) +
  geom_point (aes (shape=SubjectID, group=Species),
  position=position_jitterdodge(dodge.width=0.85, jitter.width=0.2),
  alpha=0.5) +
  scale_shape_manual (values = c(Sp24=15,Sp29=17,Sp30=19,Rp9=15,Rp11=17,Rp10=19,XL10=15,XL9=17,XL12=19)) +
  scale_fill_manual (values=wes_palette (n=3, name="FantasticFox1")) +
  theme_classic () +
  guides (shape=FALSE)


#Plot Motorneuron Counts
#Grouped by Motor Nucleus
ggplot (ARData, aes (x=MotorNuc, y=MNCount, fill=Species)) +
  geom_boxplot (alpha=0.5) +
  ylim (0, 30) +
  geom_point (position=position_jitterdodge(jitter.width=0.15), alpha=0.5) +
  scale_fill_manual (values=wes_palette (n=3, name="FantasticFox1")) +
  theme_classic ()
  
#Same graph as above with individuals denoted by different shapes

ggplot (ARData, aes (x=MotorNuc, y=MNCount, fill=Species)) +
  geom_boxplot (aes (fill=Species), position=position_dodge(width=0.85),
  alpha=0.5, outlier.shape=NA) +
  ylim (0,30) +
  geom_point (aes (shape=SubjectID, group=Species),
  position=position_jitterdodge(dodge.width=0.85, jitter.width=0.2),
  alpha=0.5) +
  scale_shape_manual (values = c(Sp24=15,Sp29=17,Sp30=19,Rp9=15,Rp11=17,Rp10=19,XL10=15,XL9=17,XL12=19)) +
  scale_fill_manual (values=wes_palette (n=3, name="FantasticFox1")) +
  theme_classic () +
  guides (shape=FALSE)