# This script was written for data analysis reported in the manuscript:
# Freshwater trematodes differ from marine trematodes in patterns connected with 
#   Division of Labor
# By Allison T Neal, Moira Stettner, Renytzabelle Ortega-Cotto, Daniel Dieringer,
#   and Lydia Reed
# Submitted to PeerJ in October 2023
# The script was written by Allison Neal and is provided to allow easy replication
#   and critical evaluation of the analysis performed.  The script author does not
#   claim to be an efficient or expert R coder, so please excuse any inefficiencies 
#   in the code.

# This script performs analysis to determine whether redial volume appears bimodal

### First, import data (see referenced script file for more information about data)
source('/Users/aneal1/OneDrive - Norwich University/2_Research/1_Trematode Research/Div of Labor/Analysis/Redial Sizes/DOLVTsizesImportCleanUp.R')

### Analyze Volume ###
# Approximate volume of redia assuming they are cylindrical: pi*r^2*h
data$vol<-pi*(data$width/2)^2*data$length
data$log.vol<-log(data$vol,10)
data$volR<-pi*(data$widthR/2)^2*data$lengthR
data$log.volR<-log(data$volR,10)


# Perform Shapiro-Wilk Test for normality on data from each snail
#   Write a function 
vol.test<-function(sample){
	# plot data
	layout(matrix(1:2,nrow=1,ncol=2))
  hist(data$log.volR[data$snail==sample],main="reliable")
  hist(data$log.vol[data$snail==sample],main="with unreliable")
	# test vs. normal distribution
	n.unr<-sum(!is.na(data$log.vol[data$snail==sample])) #sample size with unreliable measurements included
	sw.unr<-shapiro.test(data$log.vol[data$snail==sample]) # Shapiro Wilk test with unreliable measurements included
	n.r<-sum(!is.na(data$log.volR[data$snail==sample])) # sample size with only reliable measurements
	sw.r<-shapiro.test(data$log.volR[data$snail==sample]) # Shapiro Wilk test with only reliable measurements
	values<-as.numeric(c(sw.r[1],sw.r[2],n.r,sw.unr[1],sw.unr[2],n.unr)) # all numbers to be reported
	names(values)<-c("SWT W (rel)","SWT P (rel)","N (rel)","SWT W (unr)", "SWT P (unr)",'N (unr)') # labels for numbers to be reported
	return(values)
}


# This provides the results for all infections in one table
#   Create empty matrix to store results
vol.results<-matrix(data = NA, nrow = length(levels(data$snail)), ncol = 7)
#   Put snail numbers in first column of results matrix
vol.results[,1]<-levels(data$snail)
#   Returns and stores values from the statistical tests for each infection
for(i in c(1:4,6:length(levels(data$snail)))){ # excludes 5th infection, which is H109 (no data)
  vol.results[i,2:7] <- vol.test(levels(data$snail)[i])
}
# Prints results in console
vol.results

# Run function for each infection individually so that histograms can be examined
#   An examination of the histograms was used to classify each infection as bimodal/not
#   and also to identify where the break occurred for histograms that did appear
#   bimodal.
vol.test("H101")
vol.test("H102")
vol.test("H104")
vol.test("H107") # good bimodal example
vol.test("H109")
vol.test("H110")
vol.test("H127")
vol.test("H148")
vol.test("H156")
vol.test("H158")
vol.test("H70")
vol.test("H96")
vol.test("P107") # example that's not bimodal or normal
vol.test("P112")
vol.test("P113")
vol.test("P115")
vol.test("P116")
vol.test("P117") # not bimodal or normal
vol.test("P119")
vol.test("P120")
vol.test("P121")
vol.test("P122")
vol.test("P125")
vol.test("P126")
vol.test("P128")
vol.test("P130")
vol.test("P131")
vol.test("P135")
vol.test("P136")
vol.test("P137")
vol.test("P139")
vol.test("P140")
vol.test("P141")
vol.test("P150")
vol.test("P152") # bimodal, but with small peak
vol.test("P153") 
vol.test("P154")
vol.test("P155")
vol.test("P157")
vol.test("P160")
vol.test("P162")
vol.test("P163")
vol.test("P164")
vol.test("P165")
vol.test("P166")
vol.test("P170")
vol.test("P177")
vol.test("P178")
vol.test("Ph100")
vol.test("Ph107")
vol.test("Ph109")
vol.test("Ph112")
vol.test("Ph114")
vol.test("Ph115")
vol.test("Ph116")
vol.test("Ph126")
vol.test("Ph159")
vol.test("Ph95")
vol.test("Ph96")
vol.test("V141")
vol.test("V159")
vol.test("V168")
vol.test("V172")
vol.test("V182")
vol.test("V183")
vol.test("V184")


###################################
# Figure for publication
pdf(file = "/Users/aneal1/OneDrive - Norwich University/2_Research/1_Trematode Research/Div of Labor/DOL Drafts/Fig1Volumes",
    width=5,height=3.5)
layout(matrix(1:3,ncol=3),widths=1,heights=c(1,1,1))

# Bimodal examples
par(mar=c(4,3.5,2,0.1))
hist(data$log.volR[data$snail=="H107"], main = '',ylab='',xlab = '',cex.main=1.8)
title(main = "A", adj = 0, line = -0.5,cex.main = 1.8)
mtext("# rediae",side=2,at=4.5,line=2.2)
mtext("log(volume)",side = 1, at = 6.5, line = 2.2)
par(mar=c(4,2,2,0.1))
hist(data$log.volR[data$snail=="P152"],main="",ylab='',xlab="",cex.main=1.8)
mtext("log(volume)",side = 1, at = 7, line = 2.2)
title(main = "B", adj = 0.1, line = -0.5,cex.main = 1.8)
# not bimodal example
hist(data$log.volR[data$snail=="P117"],main="",ylab='',xlab="",cex.main=1.8)
title(main = "C", adj = 0, line = -0.5,cex.main = 1.8)
mtext("log(volume)",side = 1, at = 7, line = 2.2)

dev.off()