#########
# Description
#########
# This R code creates average f0 contours based on the data output by the Praat script MeasureIntensityDurationF0minF0maxF0contourpoints.praat.
# The code below assumes that your R dataframe with the measurements is called "Data".
# Your dependent variables will have different names. For example, I usually code condition, item, participant etc. in the filename 
# and then split up the filename Excel before I import the data into R because that's easiest most of the time.
# The code below refers only to one independent variable, which is called COND. You can search-and-replace that with the name of your variable,
# or adjust the code to include additional independent variables in the plot.
# I'm assuming that you have a variable called "participant" and a within-participant design, although that should be easy enough to change.
# 
# The column coding labels of the units for which the Praat script performed measurements is called "Vowel".
# If you chose to record labels of larger units that contain these intervals, they are identified in columns called "Word" and "Syllable".
# Below, I'm showing how to plot average contours for vowels by words and conditions.
#
#
# Anja Arnhold
# Version July 02, 2025
######### 

##
# Call necessary packages
##
library(ggplot2) # for plotting
library(reshape2)# needed for average contours (converting dataframe to long form)
library(Rmisc)# needed for average contours (summarySE)

##
# Convert measurements to semitones
##
# (if you want to plot in Hz, comment these lines out and adjust variable names below)
Data$F0AtPoint1 <- 12 * log2(Data$F0AtPoint1/100)
Data$F0AtPoint2 <- 12 * log2(Data$F0AtPoint2/100)
Data$F0AtPoint3 <- 12 * log2(Data$F0AtPoint3/100)
Data$F0AtPoint4 <- 12 * log2(Data$F0AtPoint4/100)
Data$F0AtPoint5 <- 12 * log2(Data$F0AtPoint5/100)
Data$F0AtPoint6 <- 12 * log2(Data$F0AtPoint6/100)
Data$F0AtPoint7 <- 12 * log2(Data$F0AtPoint7/100)
Data$F0AtPoint8 <- 12 * log2(Data$F0AtPoint8/100)
Data$F0AtPoint9 <- 12 * log2(Data$F0AtPoint9/100)
Data$F0AtPoint10 <- 12 * log2(Data$F0AtPoint10/100)

##
# Select relevant subset of the whole dataframe
##
colnames(Data) # see all names and select
# (Make sure to include all independent variables here, and any ID variables,
#  Vowel and the f0 columns have to be included, but it's ok to include more columns than necessary at this stage)
forf0 <- Data[,c("participant","COND","Word","Vowel","F0AtPoint1","F0AtPoint2",
                       "F0AtPoint3","F0AtPoint4","F0AtPoint5","F0AtPoint6", "F0AtPoint7", "F0AtPoint8", "F0AtPoint9", "F0AtPoint10")]

##
# Convert dataframe to long form
##
f0points <- melt(forf0,id=c("participant","COND","Word","Vowel"))
summary(f0points)
rm(forf0)

##
# Get averages and standard deviations
##
# (I use summarySEwithin because I usually have within-subject designs)
# (It is important to only include no more variables than you want to show in your plots here, 
# otherwise the plots will look weird.)
f0 <- summarySEwithin(f0points, measurevar="value",
                      withinvars=c("COND", "Vowel","Word","variable"),
                      idvar=c("participant"),
                      na.rm = T)
					  
##	
#Plot
##
f0plot <- ggplot(f0, aes(x=variable, y=value, colour=COND, group=COND, shape=COND)) +
  geom_errorbar(aes(ymin=value-se, ymax=value+se), width=.2) +
  geom_line(size=2) +
  geom_point(size=5.5) + facet_grid(Word ~ Vowel) + theme_bw() + ylab("F0 (in st re 100Hz)")