| title | author | date | output | editor_options | ||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
gaps |
Julie Jung |
Jan 26th, 2022 |
|
|
Presented here is the data analysis from our gaps experiment, conducted 2013 to 2017 and published in Animal Cognition in 2022. If any of this code and/or analysis proves helpful to you, please cite our paper:
Julie Jung, Ming Guo, Mark E. Crovella, J. Gregory McDaniel, and Karen M. Warkentin. Frog embryos use multiple levels of temporal pattern in risk assessment for vibration-cued escape hatching. Animal Cognition (2022). https://doi.org/10.1007/s10071-022-01634-4
library(here)
library(openxlsx)
library(janitor)
library(ggplot2)
library(ggbeeswarm)
library(cowplot)
library(MASS)
library(multcomp)
library(stargazer)
library(knitr)
library(tidyr)
library(dplyr)
library(sciplot)
library(aod)
library(lme4)
library(dunn.test)
library(DescTools) #dunnett's test
library(R.matlab)
library(car)
#FFTs
These frequences are in Hz and from what Ming called "stim 16" and what I call the "base pattern" stimulus.
FFT.df <- read.xlsx(xlsxFile="raw_data.xlsx", sheet="FFTs", stringsAsFactors = TRUE)
widerFFTs <- FFT.df %>%
pivot_longer(-Frequencies, names_to = "Treatment", values_to = "Amps")
str(widerFFTs)
max(widerFFTs$Amps)
StimColors <- c("black", "red")
Stim_fig <- ggplot(data=widerFFTs, aes(x=Frequencies, y=Amps, color=Treatment)) +
geom_line(size=1) +
scale_y_continuous(limits=c(50, 85),
breaks=c(84.58859-40, 84.58859-30, 84.58859-20, 84.58859-10, 84.58859),
labels=c("-40", " ", "-20", " ", "0"))+
scale_x_continuous(limits=c(0, 200),
breaks=c(0, 50, 100, 150, 200),
labels=c("0", " ", "100", " ", "200"))+
labs(x = "Frequency (Hz)",
y = "Relative amplitude (dB)") +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank()) +
theme_cowplot(font_size = 12, line_size = 1, font_family="Helvetica") +
theme(axis.text.y=element_text(size=12, colour= "black")) +
theme(axis.text.x=element_text(size=12, colour= "black")) +
theme(axis.title.x=element_text(size=12, colour = "black")) +
theme(axis.title.y=element_text(size=12, colour = "black")) +
theme(legend.justification = c(1, 1), legend.position = c(0.9, 0.7), legend.text=element_text(size=12)) +
scale_colour_manual(name=Treatment,
values=StimColors,
aesthetics = "colour",
guide="legend") +
theme(legend.position = "NONE") +
annotate("text", x = 160, y = 62, size=4, color="black", label = "Tines")+
annotate("text", x = 130, y = 55, size=4, color="red", label = "Tray")
Stim_fig
ggsave("~/Desktop/StimFig_v2.pdf",
plot = Stim_fig, # or give ggplot object name as in myPlot,
width = 18, height = 10,
units = "cm", # other options c("in", "cm", "mm"),
dpi = 300,
useDingbats=FALSE)
Note: The stimulus names in this file are slightly confusing and I apologize for that. It's because Ming Guo had different names for the stimuli initially (stim16 for base pattern, stim21 for primer only, etc.), and at first I renamed the stimuli letters (A-D for Experiment I) before deciding in the paper to have more descriptive names... I probably could have changed all the file names here but for fear of mucking something up accidentally I shall leave that up to you good folks.
library(tuneR)
stimA<- readWave('stimulus_21-primer_only.wav')
stimB<-readWave('stimulus_17-little_bursts.wav')
stimC<- readWave('stimulus_16-base_pattern.wav')
stimD<- readWave('stimulus_22-30s_gap.wav')
str(stimA) # the wav file has one channel (@left) containing 132300 sample points each
str(stimB) #6725250
str(stimC) #6615000
str(stimD) #7375725
# extract signal
sA<- stimA@left #we’ll select and work only with one of the channels from now onwards
sB <- stimB@left
sC <- stimC@left
sD <- stimD@left
# determine duration
#can use str() to calculate duration = sample points / sample rate
length(sA)/[email protected] #6 s
length(sB)/[email protected] #305 s
length(sC)/[email protected] #300 s
length(sD)/[email protected] #334.5 s
#If our wav file has a 16-bit depth (train_audio@bit), this means that the sound pressure values are mapped to integer values that can range from -2^¹⁵ to (2^¹⁵)-1.
#convert our sound array to floating point values ranging from -1 to 1 as follows:
sA <- sA / 2^(stimA@bit -1)
sB <- sB / 2^(stimB@bit -1)
sC <- sC / 2^(stimC@bit -1)
sD <- sD / 2^(stimD@bit -1)
#Create an array containing the time points first:
sA.timeArray <- (0:(length(sA)-1)) / [email protected]
sB.timeArray <- (0:(length(sB)-1)) / [email protected]
sC.timeArray <- (0:(length(sC)-1)) / [email protected]
sD.timeArray <- (0:(length(sD)-1)) / [email protected]
sA.df<- data.frame(sA.timeArray, sA)
sB.df<- data.frame(sB.timeArray, sB)
sC.df<- data.frame(sC.timeArray, sC)
sD.df<- data.frame(sD.timeArray, sD)
A time representation of the sound can be obtained by plotting the pressure values against the time axis.
# demean to remove DC offset
sA = sA - mean(sA)
sB = sB - mean(sB)
sC = sC - mean(sC)
sD = sD - mean(sD)
# plot waveform
sAfig<- plot(sA, type = 'l', xlab = 'Samples', ylab = 'Amplitude', xlim=c(0,length(sA)))
sBfig<- plot(sB, type = 'l', xlab = 'Samples', ylab = 'Amplitude', xlim=c(0,length(sB)))
sCfig<- plot(sC, type = 'l', xlab = 'Samples', ylab = 'Amplitude', xlim=c(0,length(sC)))
sDfig<- plot(sD, type = 'l', xlab = 'Samples', ylab = 'Amplitude', xlim=c(0,length(sD)))
sA.sub<-sA.df[seq(1, nrow(sA.df), 500), ] #take every 500 rows
sB.sub<-sB.df[seq(1, nrow(sB.df), 500), ] #take every 500 rows
sC.sub<-sC.df[seq(1, nrow(sC.df), 500), ] #take every 500 rows
sD.sub<-sD.df[seq(1, nrow(sD.df), 500), ] #take every 500 rows
sAfig <- ggplot(data=sA.sub, aes(x=sA.timeArray, y=sA)) +
geom_line(size=0.25) +
scale_x_continuous(limits=c(0, 360),
breaks=c(0, 60, 120, 180, 240, 300, 360),
labels=c("0", "60", "120", "180", "240", "300", "360"))+
scale_y_continuous(labels=NULL) +
labs(x = "Time (s)",
y = "Amplitude") +
theme_cowplot(font_size = 16, line_size = 0.75, font_family="Helvetica") +
theme(axis.text.y=element_text(size=16, colour= "black")) +
theme(axis.text.x=element_text(size=16, colour= "black")) +
theme(axis.title.x=element_text(size=16, colour = "black")) +
theme(axis.title.y=element_text(size=16, colour = "black"))+
theme(legend.position = "NONE")
sAfig
sBfig<- ggplot(data=sB.sub, aes(x=sB.timeArray, y=sB)) +
geom_line(size=0.25) +
scale_x_continuous(limits=c(0, 360),
breaks=c(0, 60, 120, 180, 240, 300, 360),
labels=c("0", "60", "120", "180", "240", "300", "360"))+
scale_y_continuous(labels=NULL) +
labs(x = "Time (s)",
y = "Amplitude") +
theme_cowplot(font_size = 16, line_size = 0.75, font_family="Helvetica") +
theme(axis.text.y=element_text(size=16, colour= "black")) +
theme(axis.text.x=element_text(size=16, colour= "black")) +
theme(axis.title.x=element_text(size=16, colour = "black")) +
theme(axis.title.y=element_text(size=16, colour = "black"))+
theme(legend.position = "NONE")
sBfig
sCfig <- ggplot(data=sC.sub, aes(x=sC.timeArray, y=sC)) +
geom_line(size=0.25) +
scale_x_continuous(limits=c(0, 360),
breaks=c(0, 60, 120, 180, 240, 300, 360),
labels=c("0", "60", "120", "180", "240", "300", "360"))+
scale_y_continuous(labels=NULL) +
labs(x = "Time (s)",
y = "Amplitude") +
theme_cowplot(font_size = 16, line_size = 0.75, font_family="Helvetica") +
theme(axis.text.y=element_text(size=16, colour= "black")) +
theme(axis.text.x=element_text(size=16, colour= "black")) +
theme(axis.title.x=element_text(size=16, colour = "black")) +
theme(axis.title.y=element_text(size=16, colour = "black"))+
theme(legend.position = "NONE")
sCfig
sDfig <- ggplot(data=sD.sub, aes(x=sD.timeArray, y=sD)) +
geom_line(size=0.25) +
scale_x_continuous(limits=c(0, 360),
breaks=c(0, 60, 120, 180, 240, 300, 360),
labels=c("0", "60", "120", "180", "240", "300", "360"))+
scale_y_continuous(labels=NULL) +
labs(x = "Time (s)",
y = "Amplitude") +
theme_cowplot(font_size = 16, line_size = 0.75, font_family="Helvetica") +
theme(axis.text.y=element_text(size=16, colour= "black")) +
theme(axis.text.x=element_text(size=16, colour= "black")) +
theme(axis.title.x=element_text(size=16, colour = "black")) +
theme(axis.title.y=element_text(size=16, colour = "black"))+
theme(legend.position = "NONE")
sDfig
combine plots
Exp1_stimuli <- plot_grid(sAfig, sBfig, sCfig, sDfig, labels = "AUTO", axis = "br", ncol = 1, nrow=4, align = 'v')
ggsave("~/Desktop/Exp1_stimuli_v2.pdf",
plot = Exp1_stimuli, # or give ggplot object name as in myPlot,
width = 18, height = 20,
units = "cm", # other options c("in", "cm", "mm"),
dpi = 300,
useDingbats=FALSE)
#stimulus fig (exp2)
stimC<- readWave('stimulus_16-base_pattern.wav')
stimE<- readWave('stimulus_27-3bursts.wav')
stimF<-readWave('stimulus_28-10bursts.wav')
str(stimC) # the wav file has one channel (@left) containing 6615000 sample points each
str(stimE) #6945750
str(stimF) #6824475
# extract signal
sC<- stimC@left #we’ll select and work only with one of the channels from now onwards
sE <- stimE@left
sF <- stimF@left
# determine duration
#can use str() to calculate duration = sample points / sample rate
length(sC)/[email protected] #300 s
length(sE)/[email protected] #315 s
length(sF)/[email protected] #309.5 s
#If our wav file has a 16-bit depth (train_audio@bit), this means that the sound pressure values are mapped to integer values that can range from -2^¹⁵ to (2^¹⁵)-1.
#convert our sound array to floating point values ranging from -1 to 1 as follows:
sC <- sC / 2^(stimC@bit -1)
sE <- sE / 2^(stimE@bit -1)
sF <- sF / 2^(stimF@bit -1)
#Create an array containing the time points first:
sC.timeArray <- (0:(length(sC)-1)) / [email protected]
sE.timeArray <- (0:(length(sE)-1)) / [email protected]
sF.timeArray <- (0:(length(sF)-1)) / [email protected]
sC.df<- data.frame(sC.timeArray, sC)
sE.df<- data.frame(sE.timeArray, sE)
sF.df<- data.frame(sF.timeArray, sF)
A time representation of the sound can be obtained by plotting the pressure values against the time axis.
# demean to remove DC offset
sC = sC - mean(sC)
sE = sE - mean(sE)
sF = sF - mean(sF)
sC.sub<-sC.df[seq(1, nrow(sC.df), 500), ] #take every 500 rows
sE.sub<-sE.df[seq(1, nrow(sE.df), 500), ] #take every 500 rows
sF.sub<-sF.df[seq(1, nrow(sF.df), 500), ] #take every 500 rows
sCfig <- ggplot(data=sC.sub, aes(x=sC.timeArray, y=sC)) +
geom_line(size=0.25) +
scale_x_continuous(limits=c(0, 360),
breaks=c(0, 60, 120, 180, 240, 300, 360),
labels=c("0", "60", "120", "180", "240", "300", "360"))+
scale_y_continuous(labels=NULL) +
labs(x = "Time (s)",
y = "Amplitude") +
theme_cowplot(font_size = 16, line_size = 0.75, font_family="Helvetica") +
theme(axis.text.y=element_text(size=16, colour= "black")) +
theme(axis.text.x=element_text(size=16, colour= "black")) +
theme(axis.title.x=element_text(size=16, colour = "black")) +
theme(axis.title.y=element_text(size=16, colour = "black"))+
theme(legend.position = "NONE")
sCfig
sEfig <- ggplot(data=sE.sub, aes(x=sE.timeArray, y=sE)) +
geom_line(size=0.25) +
scale_x_continuous(limits=c(0, 360),
breaks=c(0, 60, 120, 180, 240, 300, 360),
labels=c("0", "60", "120", "180", "240", "300", "360"))+
scale_y_continuous(labels=NULL) +
labs(x = "Time (s)",
y = "Amplitude") +
theme_cowplot(font_size = 16, line_size = 0.75, font_family="Helvetica") +
theme(axis.text.y=element_text(size=16, colour= "black")) +
theme(axis.text.x=element_text(size=16, colour= "black")) +
theme(axis.title.x=element_text(size=16, colour = "black")) +
theme(axis.title.y=element_text(size=16, colour = "black"))+
theme(legend.position = "NONE")
sEfig
sFfig<- ggplot(data=sF.sub, aes(x=sF.timeArray, y=sF)) +
geom_line(size=0.25) +
scale_x_continuous(limits=c(0, 360),
breaks=c(0, 60, 120, 180, 240, 300, 360),
labels=c("0", "60", "120", "180", "240", "300", "360"))+
scale_y_continuous(labels=NULL) +
labs(x = "Time (s)",
y = "Amplitude") +
theme_cowplot(font_size = 16, line_size = 0.75, font_family="Helvetica") +
theme(axis.text.y=element_text(size=16, colour= "black")) +
theme(axis.text.x=element_text(size=16, colour= "black")) +
theme(axis.title.x=element_text(size=16, colour = "black")) +
theme(axis.title.y=element_text(size=16, colour = "black"))+
theme(legend.position = "NONE")
sFfig
combine plots
Exp2_stimuli <- plot_grid(sCfig, sEfig, sFfig, labels = "AUTO", axis = "br", ncol = 1, nrow=3, align = 'v')
ggsave("~/Desktop/Exp2_stimuli_v2.pdf",
plot = Exp2_stimuli, # or give ggplot object name as in myPlot,
width = 18, height = 20,
units = "cm", # other options c("in", "cm", "mm"),
dpi = 300,
useDingbats=FALSE)
#stimulus fig (exp3)
#stimC<- readWave('stimulus_16-base_pattern.wav')
stimD_after3p <- readWave('stimulus_30.wav')
stimG<- readWave('stimulus_31-after_17_pulses.wav')
stimH<-readWave('stimulus_32-after_43_pulses.wav')
str(stimD_after3p) #7210350
str(stimG) #7210350
str(stimH) #7210350
# extract signal
sD3 <- stimD_after3p@left
sG <- stimG@left
sH <- stimH@left
# determine duration
#can use str() to calculate duration = sample points / sample rate
length(sD3)/[email protected] #327
length(sG)/[email protected] #327 s
length(sH)/[email protected] #327 s
#If our wav file has a 16-bit depth (train_audio@bit), this means that the sound pressure values are mapped to integer values that can range from -2^¹⁵ to (2^¹⁵)-1.
#convert our sound array to floating point values ranging from -1 to 1 as follows:
sG <- sG / 2^(stimG@bit -1)
sH <- sH / 2^(stimH@bit -1)
#Create an array containing the time points first:
sG.timeArray <- (0:(length(sG)-1)) / [email protected]
sH.timeArray <- (0:(length(sH)-1)) / [email protected]
sG.df<- data.frame(sG.timeArray, sG)
sH.df<- data.frame(sH.timeArray, sH)
A time representation of the sound can be obtained by plotting the pressure values against the time axis.
# demean to remove DC offset
sG = sG - mean(sG)
sH = sH - mean(sH)
#sC.sub<-sC.df[seq(1, nrow(sC.df), 500), ] #take every 500 rows
sG.sub<-sG.df[seq(1, nrow(sG.df), 500), ] #take every 500 rows
sH.sub<-sH.df[seq(1, nrow(sH.df), 500), ] #take every 500 rows
sGfig <- ggplot(data=sG.sub, aes(x=sG.timeArray, y=sG)) +
geom_line(size=0.25) +
scale_x_continuous(limits=c(0, 360),
breaks=c(0, 60, 120, 180, 240, 300, 360),
labels=c("0", "60", "120", "180", "240", "300", "360"))+
scale_y_continuous(labels=NULL) +
labs(x = "Time (s)",
y = "Amplitude") +
theme_cowplot(font_size = 16, line_size = 0.75, font_family="Helvetica") +
theme(axis.text.y=element_text(size=16, colour= "black")) +
theme(axis.text.x=element_text(size=16, colour= "black")) +
theme(axis.title.x=element_text(size=16, colour = "black")) +
theme(axis.title.y=element_text(size=16, colour = "black"))+
theme(legend.position = "NONE")
sGfig
sHfig<- ggplot(data=sH.sub, aes(x=sH.timeArray, y=sH)) +
geom_line(size=0.25) +
scale_x_continuous(limits=c(0, 360),
breaks=c(0, 60, 120, 180, 240, 300, 360),
labels=c("0", "60", "120", "180", "240", "300", "360"))+
scale_y_continuous(labels=NULL) +
labs(x = "Time (s)",
y = "Amplitude") +
theme_cowplot(font_size = 16, line_size = 0.75, font_family="Helvetica") +
theme(axis.text.y=element_text(size=16, colour= "black")) +
theme(axis.text.x=element_text(size=16, colour= "black")) +
theme(axis.title.x=element_text(size=16, colour = "black")) +
theme(axis.title.y=element_text(size=16, colour = "black"))+
theme(legend.position = "NONE")
sHfig
combine plots
Exp3_stimuli <- plot_grid(sCfig, sDfig, sGfig, sHfig, labels = "AUTO", axis = "br", ncol = 1, nrow=4, align = 'v')
ggsave("~/Desktop/Exp3_stimuli_v1.pdf",
plot = Exp3_stimuli, # or give ggplot object name as in myPlot,
width = 18, height = 20,
units = "cm", # other options c("in", "cm", "mm"),
dpi = 300,
useDingbats=FALSE)
#stimulus fig (exp4)
#stimC<- readWave('stimulus_16-base_pattern.wav')
#stimD<- readWave('---')
stimI<- readWave('gapped45sec.wav')
stimJ<-readWave('gapped1min.wav')
str(stimI) #7574175
str(stimJ) #7904925
# extract signal
sI <- stimI@left
sJ <- stimJ@left
# determine duration
#can use str() to calculate duration = sample points / sample rate
length(sI)/[email protected] #343.5 s
length(sJ)/[email protected] #358.5 s
#If our wav file has a 16-bit depth (train_audio@bit), this means that the sound pressure values are mapped to integer values that can range from -2^¹⁵ to (2^¹⁵)-1.
#convert our sound array to floating point values ranging from -1 to 1 as follows:
sI <- sI / 2^(stimI@bit -1)
sJ <- sJ / 2^(stimJ@bit -1)
#Create an array containing the time points first:
sI.timeArray <- (0:(length(sI)-1)) / [email protected]
sJ.timeArray <- (0:(length(sJ)-1)) / [email protected]
sI.df<- data.frame(sI.timeArray, sI)
sJ.df<- data.frame(sJ.timeArray, sJ)
sI.sub<-sI.df[seq(1, nrow(sI.df), 600), ]
sJ.sub<-sJ.df[seq(1, nrow(sJ.df), 700), ]
A time representation of the sound can be obtained by plotting the pressure values against the time axis.
# demean to remove DC offset
sI = sI - mean(sI)
sJ = sJ - mean(sJ)
sIfig <- ggplot(data=sI.sub, aes(x=sI.timeArray, y=sI)) +
geom_line(size=0.25) +
scale_x_continuous(limits=c(0, 360),
breaks=c(0, 60, 120, 180, 240, 300, 360),
labels=c("0", "60", "120", "180", "240", "300", "360"))+
scale_y_continuous(labels=NULL) +
labs(x = "Time (s)",
y = "Amplitude") +
theme_cowplot(font_size = 16, line_size = 0.75, font_family="Helvetica") +
theme(axis.text.y=element_text(size=16, colour= "black")) +
theme(axis.text.x=element_text(size=16, colour= "black")) +
theme(axis.title.x=element_text(size=16, colour = "black")) +
theme(axis.title.y=element_text(size=16, colour = "black"))+
theme(legend.position = "NONE")
sIfig
sJfig<- ggplot(data=sJ.sub, aes(x=sJ.timeArray, y=sJ)) +
geom_line(size=0.25) +
scale_x_continuous(limits=c(0, 360),
breaks=c(0, 60, 120, 180, 240, 300, 360),
labels=c("0", "60", "120", "180", "240", "300", "360"))+
scale_y_continuous(labels=NULL) +
labs(x = "Time (s)",
y = "Amplitude") +
theme_cowplot(font_size = 16, line_size = 0.75, font_family="Helvetica") +
theme(axis.text.y=element_text(size=16, colour= "black")) +
theme(axis.text.x=element_text(size=16, colour= "black")) +
theme(axis.title.x=element_text(size=16, colour = "black")) +
theme(axis.title.y=element_text(size=16, colour = "black"))+
theme(legend.position = "NONE")
sJfig
combine plots
Exp4_stimuli <- plot_grid(sCfig, sDfig, sIfig, sJfig, labels = "AUTO", axis = "br", ncol = 1, nrow=4, align = 'v')
ggsave("~/Desktop/Exp4_stimuli_v3.pdf",
plot = Exp4_stimuli, # or give ggplot object name as in myPlot,
width = 18, height = 20,
units = "cm", # other options c("in", "cm", "mm"),
dpi = 300,
useDingbats=FALSE)
#snake sample fig
snake_sample<- readWave('ImantodesSample-2003-C631b[2303].wav')
str(snake_sample) #330728
# extract signal
snake <- snake_sample@left
# determine duration
#can use str() to calculate duration = sample points / sample rate
length(snake)/[email protected] #15 s
#If our wav file has a 16-bit depth (train_audio@bit), this means that the sound pressure values are mapped to integer values that can range from -2^¹⁵ to (2^¹⁵)-1.
#convert our sound array to floating point values ranging from -1 to 1 as follows:
snake <- snake / 2^(snake_sample@bit -1)
#Create an array containing the time points first:
snake.timeArray <- (0:(length(snake)-1)) / [email protected]
snake.df<- data.frame(snake.timeArray, snake)
snake.sub<-snake.df[seq(1, nrow(snake.df), 500), ] #take every 600 rows
A time representation of the sound can be obtained by plotting the pressure values against the time axis.
# demean to remove DC offset
snake = snake - mean(snake)
snakefig <- ggplot(data=snake.sub, aes(x=snake.timeArray, y=snake)) +
geom_line(size=0.25) +
scale_x_continuous(limits=c(0, 15),
breaks=c(0, 5, 10, 15),
labels=c("0", "5", "10", "15"))+
scale_y_continuous(labels=NULL) +
labs(x = "Time (s)",
y = "Amplitude") +
theme_cowplot(font_size = 16, line_size = 0.75, font_family="Helvetica") +
theme(axis.text.y=element_text(size=16, colour= "black")) +
theme(axis.text.x=element_text(size=16, colour= "black")) +
theme(axis.title.x=element_text(size=16, colour = "black")) +
theme(axis.title.y=element_text(size=16, colour = "black"))+
theme(legend.position = "NONE")
snakefig
combine plots
ggsave("~/Desktop/snakefig2.pdf",
plot = snakefig, # or give ggplot object name as in myPlot,
width = 18, height = 13,
units = "cm", # other options c("in", "cm", "mm"),
dpi = 300,
useDingbats=FALSE)
#Hatching Data:
ALL exps I-IV
data <- read.xlsx(xlsxFile="raw_data.xlsx", sheet="data") %>%
clean_names() %>%
remove_empty() %>%
mutate(age = age_d + set_up_time_h/24 + set_up_time_m/(60*24), na.rm=T) %>% #seconds not included here
mutate(mean_stage = (stage_1+stage_2+stage_3)/3) %>% #only relevant for experiment 4
mutate(proportion_hatched_trays = number_hatched/n_test_eggs) %>% #only relevant for experiment 4
mutate(proportion_hatched_tines = hatched_during_playback/(hatched_during_playback+competent_and_no_hatch_during_playback))%>% #only relevant for experiments 1-3
mutate(diff_last_and_first = last_hatched_s - first_hatched_s) %>%
mutate(diff_75_percent_and_25_percent = x75_percent_hatched_s - x25_percent_hatched_s) %>%
mutate(diff_90_percent_and_10_percent = x90_percent_hatched_s - x10_percent_hatched_s) %>%
fill(stimulus, age, experiment, clutch) #fills stimulus name, embryo age, experiment, etc.
data$experiment<-as.factor(data$experiment)
data$stimulus<-as.factor(data$stimulus)
data$clutch<-as.factor(data$clutch)
data$set<-as.factor(data$set)
data$hatching_response<-as.factor(data$hatching_response)
Age here we've calculated as age in days of set up time
quantile(data$age, na.rm=T)
mean(data$age, na.rm=T) #reported in Methods
se(data$age) #reported in Methods
range(data$age, na.rm=T) #reported in Methods
prefixonly stimulus in experiment 1 vs. exp 2
prefix_only <- subset(data, data$stimulus=="prefixonly", na.rm=T)
glmer0 <- glmer(hatching_response ~ (1|clutch), family=binomial(logit), data=primer_only)
glmer1 <- glmer(hatching_response ~ as.factor(experiment) + (1|clutch), family=binomial(logit), data=primer_only)
Anova(glmer1) #no significant diffs between hatching response to primer only control in E1 vs. E2.
anova(glmer0, glmer1) #no significant diffs between hatching response to primer only control in E1 vs. E2.
gapsonly stimulus in experiment 1 vs. exp 2
gaps_only <- subset(data, data$stimulus=="gapsonly", na.rm=T)
glmer0 <- glmer(hatching_response ~ (1|clutch), family=binomial(logit), data=gaps_only)
glmer1 <- glmer(hatching_response ~ as.factor(experiment) + (1|clutch), family=binomial(logit), data=gaps_only)
Anova(glmer1) #no significant diffs between hatching response to gaps only control in E1 vs. E2.
anova(glmer0, glmer1) #no significant diffs between hatching response to gaps only control in E1 vs. E2.
base pattern stimulus in experiments 1 - 3
base_pattern <- subset(data, data$stimulus=="basepattern", na.rm=T)
glmer0 <- glmer(hatching_response ~ (1|clutch), family=binomial(logit), data=base_pattern)
glmer1 <- glmer(hatching_response ~ as.factor(experiment) + (1|clutch), family=binomial(logit), data=base_pattern)
Anova(glmer1) #no significant diffs between hatching response to base pattern control in E1 vs. E2. vs. E3
anova(glmer0, glmer1) #no significant diffs between hatching response to base pattern control in E1 vs. E2. vs. E3
prefix plus gap stimulus in experiment 1 vs exp 3
prefix_plus_gap <- subset(data, data$stimulus=="prefixplusgap", na.rm=T)
glmer0 <- glmer(hatching_response ~ (1|clutch), family=binomial(logit), data=prefix_plus_gap)
glmer1 <- glmer(hatching_response ~ as.factor(experiment) + (1|clutch), family=binomial(logit), data=prefix_plus_gap)
Anova(glmer1) #no significant diffs between hatching response to prefix plus gap control in E1 vs. E3
anova(glmer0, glmer1) #no significant diffs between hatching response to prefix plus gap control in E1 vs. E3
prefixonly stimulus in experiment 1 vs. exp 2
prefix_only <- subset(data, data$stimulus=="prefixonly", na.rm=T)
glm0 <- glm(latency_to_hatch_s ~ 1, family=Gamma(inverse), data=prefix_only)
glm1 <- glm(latency_to_hatch_s ~ as.factor(experiment), family=Gamma(inverse), data=prefix_only) #took out clutch as a random effect here before it was overfit -- got a singularity warning...
Anova(glm1) #no significant diffs between hatching timing to primer only control in E1 vs. E2.
anova(glm0, glm1) #no significant diffs between hatching timing to primer only control in E1 vs. E2.
gapsonly stimulus in experiment 1 vs. exp 2
gaps_only <- subset(data, data$stimulus=="gapsonly", na.rm=T)
glm0 <- glm(latency_to_hatch_s ~ 1, family=Gamma(inverse), data=gaps_only)
glm1 <- glm(latency_to_hatch_s ~ as.factor(experiment), family=Gamma(inverse), data=gaps_only)
Anova(glm1) # yes significant diffs between hatching timing to gaps only control in E1 vs. E2.
anova(glm0, glm1) # this doesn't work though
base pattern stimulus in experiments 1 - 3
base_pattern <- subset(data, data$stimulus=="basepattern", na.rm=T)
glm0 <- glm(latency_to_hatch_s ~ 1, family=Gamma(inverse), data=base_pattern)
glm1 <- glm(latency_to_hatch_s ~ as.factor(experiment), family=Gamma(inverse), data=base_pattern)
Anova(glm1) #no significant diffs between hatching response to base pattern control in E1 vs. E2. vs. E3
anova(glm0, glm1) # this doesn't work either
prefix plus gap stimulus in experiment 1 vs exp 3
prefix_plus_gap <- subset(data, data$stimulus=="prefixplusgap", na.rm=T)
glm0 <- glm(latency_to_hatch_s ~ 1, family=Gamma(inverse), data=prefix_plus_gap)
glm1 <- glm(latency_to_hatch_s ~ as.factor(experiment), family=Gamma(inverse), data=prefix_plus_gap)
Anova(glm1) #no significant diffs between hatching response to prefix plus gap control in E1 vs. E3
anova(glm0, glm1) #
#Exp 1
####fig
exp1_combined <- subset(data, data$experiment==1)
summary_e1<-
exp1_combined %>%
group_by(stimulus) %>%
summarize(count = n(), #number of individuals per stimulus REPORTED
mean = mean(proportion_hatched_tines, na.rm=T),
SE = se(proportion_hatched_tines, na.rm=T), #NA bc only 1 per clutch
SD = sd(proportion_hatched_tines, na.rm=T)
)
kable(summary_e1,title="Mean & SE",digits=4)
N_clutches_e1<-
exp1_combined %>%
group_by(stimulus, clutch) %>%
summarize(count = n(), #number of individuals per clutch (each row here is a clutch)
mean = mean(proportion_hatched_tines, na.rm=T),
SE = se(proportion_hatched_tines, na.rm=T) #NA bc only 1 per clutch
)
kable(N_clutches_e1,title="Mean & SE",digits=4)
N_e1<-
N_clutches_e1 %>%
group_by(stimulus) %>%
summarize(count = n(), #number of clutches (REPORTED IN FIG 4)
meanProportionHatched = mean(mean, na.rm=T),
SEProportionHatched = se(mean, na.rm=T)
)
kable(N_e1,title="Mean & SE",digits=4)
Since the data we are visualizing here is pretty small, why not forgo the boxplot?
# Function to produce summary statistics (mean and +/- se)
data_summary <- function(x) {
m <- mean(x, na.rm=T)
ymin <- m-se(x, na.rm=T)
ymax <- m+se(x, na.rm=T)
return(c(y=m,ymin=ymin,ymax=ymax))
}
exp1_combined$stimulus <- factor(exp1_combined$stimulus , levels=c("basepattern", "prefixplusgap", "prefixonly", "gapsonly"))
fig_e1 <- exp1_combined %>% ggplot() +
aes(x=as.factor(stimulus), y=proportion_hatched_tines) +
geom_dotplot(dotsize=0.5, binaxis='y', stackdir='center', position=position_dodge(), fill='white', stroke=2, binwidth=1/30) +
stat_summary(fun.data=data_summary, color='red')+
scale_y_continuous(limits=c(0, 1),
breaks=c(0, 0.25, 0.5, 0.75, 1),
labels=c("0", " ", "0.5", " ", "1")
)+
scale_x_discrete(labels=c("base pattern", "prefix + gap", "prefix only", "gaps only")) +
labs(x = "Stimulus",
y = "Proportion hatched") +
theme(panel.background = element_blank())+
#theme_cowplot(font_size = 15, line_size = 0.75, font_family="Helvetica") +
theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
#theme(axis.text.y=element_text(colour= "black")) +
#theme(axis.text.x=element_text(colour= "black")) +
#theme(axis.title.x=element_text(colour = "black")) +
#theme(axis.title.y=element_text(colour = "black")) +
theme(legend.position="none")+
theme(axis.ticks.length=unit(-0.25, "cm"), axis.text.x = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm")), axis.text.y = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm")) )
fig_e1
ggsave("~/Desktop/fig_e1.pdf",
plot = fig_e1, # or give ggplot object name as in myPlot,
width = 16, height = 16,
units = "cm", # other options c("in", "cm", "mm"),
dpi = 400,
useDingbats=FALSE)
compare "prefixonly" vs. "gapsonly" vs. "basepattern" vs. "prefixplusgap"
shapiro.test(exp1_combined$proportion_hatched_tines) #P<0.05, so not normal
kruskal.test(proportion_hatched_tines ~ stimulus, data=exp1_combined) #reported in results exp 1
test<-pairwise.wilcox.test(exp1_combined$proportion_hatched_tines,
exp1_combined$stimulus,
p.adjust.method = "none")
library(multcompView)
test$p.value
library(reshape)
(a <- melt(test$p.value))
a.cc <- na.omit(a)
a.pvals <- a.cc[, 3]
names(a.pvals) <- paste(a.cc[, 1], a.cc[, 2], sep="-")
a.pvals
multcompLetters(a.pvals)
##midpoint hatching timing
nogaps_BPvsPPG <- subset(data, data$experiment==1 & (data$stimulus=="basepattern" | data$stimulus=="prefixplusgap") , na.rm=T)
kruskal.test(x50_percent_hatched_s ~ stimulus, data=nogaps_BPvsPPG)
pairwise.wilcox.test(nogaps_BPvsPPG$x50_percent_hatched_s,
nogaps_BPvsPPG$gaps_removed,
p.adjust.method = "none")
##midrange hatching timing
pairwise.wilcox.test(nogaps_BPvsPPG$diff_75_percent_and_25_percent,
nogaps_BPvsPPG$gaps_removed,
p.adjust.method = "none")
#Exp 2
####fig
exp2_combined <- subset(data, data$experiment==2)
N_summary_e2<-
exp2_combined %>%
group_by(stimulus) %>%
summarize(count = n(), #number of individuals per stimulus
mean = mean(proportion_hatched_tines, na.rm=T),
SE = se(proportion_hatched_tines, na.rm=T)
)
kable(N_summary_e2,title="Mean & SE",digits=4)
N_clutches_e2<-
exp2_combined %>%
group_by(stimulus, clutch) %>%
summarize(count = n(), #number of individuals per clutch (each row here is a clutch)
mean = mean(proportion_hatched_tines, na.rm=T),
SE = se(proportion_hatched_tines, na.rm=T) #NA bc only 1 per clutch
)
kable(N_clutches_e2,title="Mean & SE",digits=4)
N_e2<-
N_clutches_e2 %>%
group_by(stimulus) %>%
summarize(count = n(), #number of clutches
meanProportionHatched = mean(mean, na.rm=T),
SEProportionHatched = se(mean, na.rm=T)
)
kable(N_e2,title="Mean & SE",digits=4)
with dotplot
exp2_combined$stimulus <- factor(exp2_combined$stimulus , levels=c("basepattern", "tenpulsegroups", "threepulsegroups"))
fig_e2 <- exp2_combined %>% ggplot() +
aes(x=as.factor(stimulus), y=proportion_hatched_tines) +
geom_dotplot(dotsize=0.5, binaxis='y', stackdir='center', position=position_dodge(), fill='white', stroke=2, binwidth=1/30) +
stat_summary(fun.data=data_summary, color='red')+
scale_y_continuous(limits=c(0, 1),
breaks=c(0, 0.25, 0.5, 0.75, 1),
labels=c("0", " ", " ", " ", "1")
)+
scale_x_discrete(labels=c("base pattern", "10-pulse groups", "3-pulse groups")) +
labs(x = NULL,
y = "Proportion of clutch hatched") +
theme(panel.background = element_blank())+
#theme_cowplot(font_size = 15, line_size = 0.75, font_family="Helvetica") +
theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
#theme(axis.text.y=element_text(colour= "black")) +
#theme(axis.text.x=element_text(colour= "black")) +
#theme(axis.title.x=element_text(colour = "black")) +
#theme(axis.title.y=element_text(colour = "black")) +
theme(legend.position="none")+
theme(axis.ticks.length=unit(-0.25, "cm"), axis.text.x = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm")), axis.text.y = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm")) )
fig_e2
ggsave("~/Desktop/fig_e2.pdf",
plot = fig_e2, # or give ggplot object name as in myPlot,
width = 16, height = 16,
units = "cm", # other options c("in", "cm", "mm"),
dpi = 400,
useDingbats=FALSE)
shapiro.test(exp2_combined$proportion_hatched_tines) #P<0.05, so not normal!
kruskal.test(proportion_hatched_tines ~ stimulus, data=exp2_combined)
test<-pairwise.wilcox.test(exp2_combined$proportion_hatched_tines,
exp2_combined$stimulus,
p.adjust.method = "none")
library(multcompView)
test$p.value
library(reshape)
(a <- melt(test$p.value))
a.cc <- na.omit(a)
a.pvals <- a.cc[, 3]
names(a.pvals) <- paste(a.cc[, 1], a.cc[, 2], sep="-")
a.pvals
multcompLetters(a.pvals)
If you consider a per-pulse basis, i.e. if you calculate proportion hatched by the number of pulses in each stimulus (PH/150 or PH/30, etc) is the hatching response actually relatively higher for the stimulus with fewer pulses? Probably yes? even though N is low.
exp2_combined <- exp2_combined %>%
mutate(perpulse = case_when(stimulus=="basepattern" ~ proportion_hatched_tines/150,
stimulus=="tenpulsegroups" ~ proportion_hatched_tines/70,
stimulus=="threepulsegroups" ~ proportion_hatched_tines/30)) %>%
mutate(normalized_perpulse = (perpulse-min(perpulse, na.rm=T))/(max(perpulse, na.rm=T)-min(perpulse, na.rm=T)))
shapiro.test(exp2_combined$perpulse) #P<0.05, so not normal
kruskal.test(perpulse ~ stimulus, data=exp2_combined)
test<-pairwise.wilcox.test(exp2_combined$perpulse,
exp2_combined$stimulus,
p.adjust.method = "none")
#library(multcompView)
test$p.value
#library(reshape)
(a <- melt(test$p.value))
a.cc <- na.omit(a)
a.pvals <- a.cc[, 3]
names(a.pvals) <- paste(a.cc[, 1], a.cc[, 2], sep="-")
a.pvals
multcompLetters(a.pvals)
Revise figure
exp2_combined$stimulus <- factor(exp2_combined$stimulus , levels=c("basepattern", "tenpulsegroups", "threepulsegroups"))
range(exp2_combined$normalized_perpulse, na.rm=T)
fig_e2_perpulse <- exp2_combined %>% ggplot() +
aes(x=as.factor(stimulus), y=normalized_perpulse) +
geom_dotplot(dotsize=0.5, binaxis='y', stackdir='center', position=position_dodge(), fill='white', stroke=2, binwidth=1/30) +
stat_summary(fun.data=data_summary, color='red')+
geom_dotplot(dotsize=0.5, binaxis='y', stackdir='center', position=position_dodge(), fill='white', stroke=2, binwidth=1/30) +
scale_y_continuous(limits=c(0, 1),
breaks=c(0, 0.25, 0.5, 0.75, 1),
labels=c("0", " ", " ", " ", "1")
)+
scale_x_discrete(labels=c("base pattern", "10-pulse groups", "3-pulse groups")) +
labs(x = NULL,
y = "Hatching response (normalized, per pulse)") +
theme(panel.background = element_blank())+
theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
theme(legend.position="none")+
theme(axis.ticks.length=unit(-0.25, "cm"), axis.text.x = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm")), axis.text.y = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm")) )
fig_e2_perpulse
ggsave("~/Desktop/fig_e2_perpulse.pdf",
plot = fig_e2_perpulse, # or give ggplot object name as in myPlot,
width = 16, height = 16,
units = "cm", # other options c("in", "cm", "mm"),
dpi = 400,
useDingbats=FALSE)
##midpoint hatching timing
nogaps_3vs10vsBP <- subset(data, data$experiment==2 & (data$stimulus=="basepattern" | data$stimulus=="threepulsegroups"| data$stimulus=="tenpulsegroups") , na.rm=T)
kruskal.test(x50_percent_hatched_s ~ stimulus, data=nogaps_3vs10vsBP)
pairwise.wilcox.test(nogaps_3vs10vsBP$x50_percent_hatched_s,
nogaps_3vs10vsBP$gaps_removed,
p.adjust.method = "none")
##midrange hatching timing
pairwise.wilcox.test(nogaps_3vs10vsBP$diff_75_percent_and_25_percent,
nogaps_3vs10vsBP$gaps_removed,
p.adjust.method = "none")
exp3_combined <- subset(data, data$experiment==3)
summary_e3<-
exp3_combined %>%
group_by(stimulus) %>%
summarize(count = n(), #number of individuals per stimulus
mean = mean(proportion_hatched_tines, na.rm=T),
SE = se(proportion_hatched_tines, na.rm=T), #NA bc only 1 per clutch
SD = sd(proportion_hatched_tines, na.rm=T)
)
kable(summary_e3,title="Mean & SE",digits=4)
N_clutches_e3<-
exp3_combined %>%
group_by(stimulus, clutch) %>%
summarize(count = n(), #number of individuals per clutch (each row here is a clutch)
mean = mean(proportion_hatched_tines, na.rm=T),
SE = se(proportion_hatched_tines, na.rm=T) #NA bc only 1 per clutch
)
kable(N_clutches_e3,title="Mean & SE",digits=4)
N_e3<-
N_clutches_e3 %>%
group_by(stimulus) %>%
summarize(count = n(), #number of clutches
meanProportionHatched = mean(mean, na.rm=T),
SEProportionHatched = se(mean, na.rm=T)
)
kable(N_e3,title="Mean & SE",digits=4)
with dotplot
exp3_combined$stimulus <- factor(exp3_combined$stimulus , levels=c("basepattern", "prefixplusgap", "gapafterseventeen", "gapafterfortythree"))
fig_e3 <- exp3_combined %>% ggplot() +
aes(x=as.factor(stimulus), y=proportion_hatched_tines) +
geom_dotplot(dotsize=0.5, binaxis='y', stackdir='center', position=position_dodge(), fill='white', stroke=2, binwidth=1/30) +
stat_summary(fun.data=data_summary, color='red')+
scale_y_continuous(limits=c(0, 1),
breaks=c(0, 0.25, 0.5, 0.75, 1),
labels=c("0", " ", " ", " ", "1")
)+
scale_x_discrete(labels=c("base pattern", "assessment (gap at 3 pulses)", "peak hatching (gap at 17 pulses)", "post-peak hatching (gap at 43 pulses)")) +
labs(x = NULL,
y = "Proportion of clutch hatched") +
theme(panel.background = element_blank())+
#theme_cowplot(font_size = 15, line_size = 0.75, font_family="Helvetica") +
theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
#theme(axis.text.y=element_text(colour= "black")) +
#theme(axis.text.x=element_text(colour= "black")) +
#theme(axis.title.x=element_text(colour = "black")) +
#theme(axis.title.y=element_text(colour = "black")) +
theme(legend.position="none")+
theme(axis.ticks.length=unit(-0.25, "cm"), axis.text.x = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm")), axis.text.y = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm")) )
fig_e3
ggsave("~/Desktop/fig_e3.pdf",
plot = fig_e3, # or give ggplot object name as in myPlot,
width = 16, height = 16,
units = "cm", # other options c("in", "cm", "mm"),
dpi = 400,
useDingbats=FALSE)
shapiro.test(exp3_combined$proportion_hatched_tines) #P.0.05, so normal! BUT I'LL stay consistent with others
kruskal.test(proportion_hatched_tines ~ stimulus, data=exp3_combined)
test<-pairwise.wilcox.test(exp3_combined$proportion_hatched_tines,
exp3_combined$stimulus,
p.adjust.method = "none")
library(multcompView)
test$p.value
library(reshape)
(a <- melt(test$p.value))
a.cc <- na.omit(a)
a.pvals <- a.cc[, 3]
names(a.pvals) <- paste(a.cc[, 1], a.cc[, 2], sep="-")
a.pvals
multcompLetters(a.pvals)
##midpoint hatching timing
nogaps_e3allstims <- subset(data, data$experiment==3, na.rm=T)
kruskal.test(x50_percent_hatched_s ~ stimulus, data=nogaps_e3allstims)
pairwise.wilcox.test(nogaps_e3allstims$x50_percent_hatched_s,
nogaps_e3allstims$gaps_removed,
p.adjust.method = "none")
##midrange hatching timing
pairwise.wilcox.test(nogaps_e3allstims$diff_75_percent_and_25_percent,
nogaps_e3allstims$gaps_removed,
p.adjust.method = "none")
####fig
exp4_combined <- subset(data, data$experiment==4)
To get number of clutches
N_summary_e4<-
exp4_combined %>%
group_by(stimulus) %>%
summarize(count = n(), #number of individuals per stimulus
mean = mean(proportion_hatched_trays, na.rm=T),
SE = se(proportion_hatched_trays, na.rm=T)
)
kable(N_summary_e4,title="Mean & SE",digits=4)
N_clutches_e4<-
exp4_combined %>%
group_by(stimulus, clutch) %>%
summarize(count = n(), #number of individuals per clutch (each row here is a clutch)
mean = mean(proportion_hatched_trays, na.rm=T),
SE = se(proportion_hatched_trays, na.rm=T) #sometimes multiple per clutch!!
)
kable(N_clutches_e4,title="Mean & SE",digits=4)
# number of trays per "no gap" = N clutches + 2
# number of trays per "45 s gap" = N clutches + 1
N_e4<-
N_clutches_e4 %>%
group_by(stimulus) %>%
summarize(count = n(), #number of clutches (REPORTED)
meanProportionHatched = mean(mean, na.rm=T),
SEProportionHatched = se(mean, na.rm=T)
)
kable(N_e4,title="Mean & SE",digits=4)
# "no gap" = N = 35
# "45 s gap" N = 17
exp4_combined$stimulus <- factor(exp4_combined$stimulus, levels=c("nogap", "thirtysgap", "fortyfivesgap", "sixtysgap"))
fig_e4 <- exp4_combined %>% ggplot() +
aes(x=as.factor(stimulus), y=proportion_hatched_trays) +
geom_dotplot(dotsize=0.5, binaxis='y', stackdir='center', position=position_dodge(), fill='white', stroke=2, binwidth=1/30) +
stat_summary(fun.data=data_summary, color='red')+
scale_y_continuous(limits=c(0, 1),
breaks=c(0, 0.25, 0.5, 0.75, 1),
labels=c("0", " ", " ", " ", "1")
)+
scale_x_discrete(labels=c("base pattern", "30 s gap", "45 s gap", "60 s gap")) +
labs(x = NULL,
y = "Proportion of tray hatched") +
theme(panel.background = element_blank())+
theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
theme(legend.position="none")+
theme(axis.ticks.length=unit(-0.25, "cm"), axis.text.x = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm")), axis.text.y = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm")) )
fig_e4
ggsave("~/Desktop/fig_e4.pdf",
plot = fig_e4, # or give ggplot object name as in myPlot,
width = 16, height = 16,
units = "cm", # other options c("in", "cm", "mm"),
dpi = 400,
useDingbats=FALSE)
shapiro.test(exp4_combined$proportion_hatched_trays) #P<0.05, so not normal
kruskal.test(proportion_hatched_trays ~ stimulus, data=exp4_combined)
test<-pairwise.wilcox.test(exp4_combined$proportion_hatched_trays,
exp4_combined$stimulus,
p.adjust.method = "none")
#library(multcompView)
test$p.value
#library(reshape)
(a <- melt(test$p.value))
a.cc <- na.omit(a)
a.pvals <- a.cc[, 3]
names(a.pvals) <- paste(a.cc[, 1], a.cc[, 2], sep="-")
a.pvals
multcompLetters(a.pvals)
N_summary_e4_ltoh<- exp4_combined %>%
filter(!is.na(latency_to_hatch_s)) %>% #to filter out the NAs
group_by(stimulus) %>%
summarize(count = n(), #number of individuals per stimulus
mean = mean(latency_to_hatch_s, na.rm=T),
SE = se(latency_to_hatch_s, na.rm=T) #sometimes multiple per clutch!! and sometimes no latency available bc none hatched
)
kable(N_summary_e4_ltoh,title="Mean & SE",digits=4)
N_clutches_e4_ltoh<-
exp4_combined %>%
filter(!is.na(latency_to_hatch_s)) %>% #to filter out the NAs
group_by(stimulus, clutch) %>%
summarize(count = n(), #number of individuals per clutch (each row here is a clutch)--- i should count these up to make sure it lines up with what's reported in the figures right now.
mean = mean(latency_to_hatch_s, na.rm=T),
SE = se(latency_to_hatch_s, na.rm=T) #sometimes multiple per clutch!! and sometimes no latency available bc none hatched
)
kable(N_clutches_e4_ltoh,title="Mean & SE",digits=4)
#45s gap : c249 has 2 latencies
#no gap: c93 has 2 latencies
N_e4_ltoh<-
N_clutches_e4_ltoh %>%
group_by(stimulus) %>%
summarize(count = n(), #number of CLUTCHES
meanLtoH = mean(mean, na.rm=T),
SELtoH = se(mean, na.rm=T)
)
kable(N_e4_ltoh,title="Mean & SE",digits=4)
#45s gap: N trays = 16 + 1 = 17
#no gap: N trays = 29 + 1 = 30
#right now in figure 8b:
# big font : N of trays (subscript = N embryos)
# 8c) N of trays tested for latencies (subscript = N embryos for which we got latencies - difference is number of trays where no embryos hatched)
# 8d) N of trays tested for latencies (subscript = N embryos for which we got latencies - difference is number of trays where no embryos hatched)
range(exp4_combined$latency_to_hatch_s, na.rm = T)
fig_e4_ltoh <- exp4_combined %>% ggplot() +
aes(x=as.factor(stimulus), y=latency_to_hatch_s) +
#geom_dotplot(dotsize=0.5, binaxis='y', stackdir='center', position=position_dodge(), fill='white', stroke=2, binwidth=1/30) +
geom_dotplot(dotsize=0.5, binaxis='y', stackdir='center', position=position_dodge(), fill='white', stroke=2, binwidth=3.5) + #use binwidth=3.8 or 4 to make a little bigger
stat_summary(fun.data=data_summary, color='red')+
scale_y_continuous(limits=c(0, 140),
breaks=c(0, 30, 60, 90, 120),
labels=c("0", " ", "60", " ", "120")
)+
scale_x_discrete(labels=c("base pattern", "30 s gap", "45 s gap", "60 s gap")) +
labs(x = NULL,
y = "Latency to hatch (s)") +
theme(panel.background = element_blank())+
theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
theme(legend.position="none")+
theme(axis.ticks.length=unit(-0.25, "cm"), axis.text.x = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm")), axis.text.y = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm")) )
fig_e4_ltoh
ggsave("~/Desktop/e4_ltoh_plots.pdf",
plot = fig_e4_ltoh, # or give ggplot object name as in myPlot,
width = 16, height = 16,
units = "cm", # other options c("in", "cm", "mm"),
dpi = 400,
useDingbats=FALSE)
shapiro.test(exp4_combined$latency_to_hatch_s) #P<0.05, so not normal
kruskal.test(latency_to_hatch_s ~ stimulus, data=exp4_combined)
test<-pairwise.wilcox.test(exp4_combined$latency_to_hatch_s, exp4_combined$stimulus,
p.adjust.method = "none")
#library(multcompView)
test$p.value
#library(reshape)
(a <- melt(test$p.value))
a.cc <- na.omit(a)
a.pvals <- a.cc[, 3]
names(a.pvals) <- paste(a.cc[, 1], a.cc[, 2], sep="-")
a.pvals
multcompLetters(a.pvals)
###subtract gaps
Using conditional formatting in dplyr
#if tines_trays == stim C --> no change to lto_h
#if tines_trays == stim D and lto_h > 0.5 (30 sec), then subtract 0.5 from lto_h
#if tines_trays == stim I and lto_h > 0.75 (45 sec), then subtract 0.75 from lto_h
#if tines_trays == stim J and lto_h > 1 (min), then subtract 1 from lto_h
exp4_combined <- exp4_combined %>%
mutate(trunc_ltoh = case_when(stimulus=="nogap" ~ latency_to_hatch_s,
stimulus=="thirtysgap" & latency_to_hatch_s>30 ~ latency_to_hatch_s-30,
stimulus=="fortyfivesgap" & latency_to_hatch_s>45 ~ latency_to_hatch_s-45,
stimulus=="sixtysgap" & latency_to_hatch_s>60 ~ latency_to_hatch_s-60,
TRUE ~ as.numeric(as.character(latency_to_hatch_s))
))
Values that are never matched by a logical test get a default replacement value: NA. To keep these values from getting NAs, include a final catch-all test and replacement in the last line.
how many hatchlings do we have latency data for? (which which we got latency)
N_summary_e4_truncltoh<- exp4_combined %>%
filter(!is.na(trunc_ltoh)) %>% #to filter out the NAs
group_by(stimulus) %>%
summarize(count = n(), #number of individuals per stimulus
mean = mean(trunc_ltoh, na.rm=T),
SE = se(trunc_ltoh, na.rm=T) #sometimes multiple per clutch!! and sometimes no latency available bc none hatched
)
kable(N_summary_e4_truncltoh,title="Mean & SE",digits=4)
N_clutches_e4_ltoh<-
exp4_combined %>%
filter(!is.na(trunc_ltoh)) %>% #to filter out the NAs
group_by(stimulus, clutch) %>%
summarize(count = n(), #number of individuals per clutch (each row here is a clutch)--- i should count these up to make sure it lines up with what's reported in the figures right now.
mean = mean(trunc_ltoh, na.rm=T),
SE = se(trunc_ltoh, na.rm=T) #sometimes multiple per clutch!! and sometimes no latency available bc none hatched
)
kable(N_clutches_e4_ltoh,title="Mean & SE",digits=4)
N_e4_ltoh<-
N_clutches_e4_ltoh %>%
group_by(stimulus) %>%
summarize(count = n(), #number of CLUTCHES
meanLtoH = mean(mean, na.rm=T),
SELtoH = se(mean, na.rm=T)
)
kable(N_e4_ltoh,title="Mean & SE",digits=4)
latency_clutches_e4<-
exp4_combined %>%
group_by(stimulus, clutch) %>%
summarize(num = n(), #sample sizes for individuals per clutch
meanlat = mean(trunc_ltoh, na.rm=TRUE)
)
kable(latency_clutches_e4, digits=3)
N_latency_clutches_e4 <- latency_clutches_e4 %>%
group_by(stimulus) %>%
summarise(clutches= n(), #number of clutches
meanlatency = mean(meanlat, na.rm=TRUE),
medianlatency = median(meanlat, na.rm=TRUE),
SEnumlatency = sd(meanlat, na.rm=TRUE)/sqrt(n())
)
kable(N_latency_clutches_e4, digits=3)
range(exp4_combined$trunc_ltoh, na.rm=T)
fig_e4_ltoh_trunc <- exp4_combined %>% ggplot() +
aes(x=as.factor(stimulus), y=trunc_ltoh) +
#geom_dotplot(dotsize=0.5, binaxis='y', stackdir='center', position=position_dodge(), fill='white', stroke=2, binwidth=1/30) +
geom_dotplot(dotsize=0.5, binaxis='y', stackdir='center', position=position_dodge(), fill='white', stroke=2, binwidth=3.5) + #use binwidth=3.8 or 4 to make a little bigger
stat_summary(fun.data=data_summary, color='red')+
scale_y_continuous(limits=c(0, 140),
breaks=c(0, 30, 60, 90, 120),
labels=c("0", " ", "30", " ", "60") #converting from seconds to # of vibration pulses bc 1 pulse every 2 seconds
)+
scale_x_discrete(labels=c("base pattern", "30 s gap", "45 s gap", "60 s gap")) +
labs(x = NULL,
y = "Latency (# of vibration pulses)") +
theme(panel.background = element_blank())+
theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
theme(legend.position="none")+
theme(axis.ticks.length=unit(-0.25, "cm"), axis.text.x = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm")), axis.text.y = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm")) )
fig_e4_ltoh_trunc
ggsave("~/Desktop/e4_ltoh_boxplots_trunc.pdf",
plot = fig_e4_ltoh_trunc, # or give ggplot object name as in myPlot,
width = 16, height = 16,
units = "cm", # other options c("in", "cm", "mm"),
dpi = 400,
useDingbats=FALSE)
shapiro.test(exp4_combined$trunc_ltoh) #P<0.05, so not normal
kruskal.test(trunc_ltoh ~ stimulus, data=exp4_combined)
test<-pairwise.wilcox.test(exp4_combined$trunc_ltoh,
exp4_combined$stimulus,
p.adjust.method = "none")
#library(multcompView)
test$p.value
#library(reshape)
(a <- melt(test$p.value))
a.cc <- na.omit(a)
a.pvals <- a.cc[, 3]
names(a.pvals) <- paste(a.cc[, 1], a.cc[, 2], sep="-")
a.pvals
multcompLetters(a.pvals)
what fraction of total hatched embryos hatched during or after these last 3 "extra" pulses in the "prefix + gap" stimulus (stimulus D or stimulus 22 for Ming)?
#hatching times for stimulus 22 (stim D, primer + gaps), experiment 1
hatching_times_D <- readMat("hatching_times_22.mat")
#stimulus length = 6 + 28.5 + 300 = 334.5
#334.5 s
##334.5 s-6s = 328.5
# matrix of 12 lists --> 12 trials of this stimulus (12 clutches)
# each list = hatching times for any embryos that hatched
df<-data.frame(unlist(hatching_times_D), nrow=length(hatching_times_D), byrow=TRUE)
#any embryos that hatched and the TIMES at which they hatched (in seconds)
#how many hatched? i.e.how many elements in this matrix?
length(df$unlist.hatching_times_D.) #340 individuals hatched
#how much hatching in the LAST 3 pulses (6 seconds)?
subsetdf<- subset(df, unlist.hatching_times_D.>=328.5)
length(subsetdf$unlist.hatching_times_D.) #2 individuals hatched in the last 6 seconds
2/340 #2 hatched in the last 6 seconds out of 340
how many hatching did the gaps only stimulus (only 10 pulses) elicit, compared to hatching from the first 10 pulses of other stimuli?
could provide evidence that this is probs not simply a reflection of the low total number of pulses in this stimulus - i.e. if the more rapid accumulation of 10 pulses in other stimuli elicited more hatching...
SO How much hatching occurred in response to the first 10 pulses of all the other stimuli?
gaps only stimulus from exp 1
#hatching times for stimulus 17_1 (stim B, gapsonly), experiment 1
hatching_times_17_1 <- readMat("hatching_times_17_1.mat")
# matrix of 10 lists --> 10 trials of this stimulus (10 clutches)
# each list = hatching times for any embryos that hatched
df<-data.frame(unlist(hatching_times_17_1), nrow=length(hatching_times_17_1), byrow=TRUE)
#any embryos that hatched and the TIMES at which they hatched (in seconds)
#how many hatched? i.e.how many elements in this matrix?
length(df$unlist.hatching_times_17_1.) #13 individuals hatched
#how much hatching in the first 10 pulses?
subsetdf<- subset(df, unlist.hatching_times_17_1.<=(2*10))
length(subsetdf$unlist.hatching_times_17_1.) #2 individuals hatched in the first 10 pulses (20 seconds)
2/13 #only 2 out of 13 hatched in the first 10 pulses (20 seconds)