# R version: 4.2.0
library(tidyverse)
library(data.table)

################
### DATA
## load final data set for obstruent tokens
data_obs = read.table("Obs_V1_trackdata_sd_ws_zh_younger.txt")

###################
### FLMM
## data preparation for FLMM 4: MSG vs. ZSG
# exclude speakers from Vienna
data_msg_zsg <- data_obs %>%
  subset(variety != "Austrian_Standard")

# select, rename and convert relevant columns for FLMM
curve_info_msg_zsg = data_msg_zsg %>%
  select(c("bundle", "session", "word", "repetition", "scaledF0", "t.norm", "variety", "V1.quantity")) %>%
  rename(n_long = bundle, subject_long = session, word_long = word, combi_long = repetition,
         y_vec = scaledF0, t = t.norm,
         covariate.1 = variety, covariate.2 = V1.quantity) %>%
  mutate_at(c("n_long", "subject_long", "word_long"), as.factor)

curve_info_msg_zsg = curve_info_msg_zsg %>%
  mutate_at(c("n_long", "subject_long", "word_long"), as.integer)

# set long vowel tokens by MSG speakers as reference mean
curve_info_msg_zsg$covariate.1 <- if_else(curve_info_msg_zsg$covariate.1 == "German_Standard", 0, 1)
curve_info_msg_zsg$covariate.2 <- if_else(curve_info_msg_zsg$covariate.2 == "long", 0, 1)

# get mean and range of data points per f0 trajectory included in the analysis
n_msg_zsg = count(curve_info_msg_zsg, n_long)
mean(n_msg_zsg$n)
range(n_msg_zsg$n)

# convert data.frame to data.table
curve_info_msg_zsg <- as.data.table(curve_info_msg_zsg)
curve_info_msg_zsg = setDT(curve_info_msg_zsg)

## data preparation for FLMM 5: MSG vs. VSG
# exclude speakers from Zurich
data_msg_vsg <- data_obs %>%
  subset(variety != "Swiss_Standard")

# select, rename and convert relevant columns for FLMM
curve_info_msg_vsg = data_msg_vsg %>%
  select(c("bundle", "session", "word", "repetition", "scaledF0", "t.norm", "variety", "V1.quantity")) %>%
  rename(n_long = bundle, subject_long = session, word_long = word, combi_long = repetition,
         y_vec = scaledF0, t = t.norm,
         covariate.1 = variety, covariate.2 = V1.quantity) %>%
  mutate_at(c("n_long", "subject_long", "word_long"), as.factor)

curve_info_msg_vsg = curve_info_msg_vsg %>%
  mutate_at(c("n_long", "subject_long", "word_long"), as.integer)

# set long vowel tokens by MSG speakers as reference mean
curve_info_msg_vsg$covariate.1 <- if_else(curve_info_msg_vsg$covariate.1 == "German_Standard", 0, 1)
curve_info_msg_vsg$covariate.2 <- if_else(curve_info_msg_vsg$covariate.2 == "long", 0, 1)

# get mean and range of data points per f0 trajectory included in the analysis
n_msg_vsg = count(curve_info_msg_vsg, n_long)
mean(n_msg_vsg$n)
range(n_msg_vsg$n)

# convert data.frame to data.table
curve_info_msg_vsg <- as.data.table(curve_info_msg_vsg)
curve_info_msg_vsg = setDT(curve_info_msg_vsg)

## data preparation for FLMM 6: ZSG vs. VSG
# exclude speakers from Munich
data_zsg_vsg <- data_obs %>%
  subset(variety != "German_Standard")

# select, rename and convert relevant columns for FLMM
curve_info_zsg_vsg = data_zsg_vsg %>%
  select(c("bundle", "session", "word", "repetition", "scaledF0", "t.norm", "variety", "V1.quantity")) %>%
  rename(n_long = bundle, subject_long = session, word_long = word, combi_long = repetition,
         y_vec = scaledF0, t = t.norm,
         covariate.1 = variety, covariate.2 = V1.quantity) %>%
  mutate_at(c("n_long", "subject_long", "word_long"), as.factor)

curve_info_zsg_vsg = curve_info_zsg_vsg %>%
  mutate_at(c("n_long", "subject_long", "word_long"), as.integer)

# set long vowel tokens by ZSG speakers as reference mean
curve_info_zsg_vsg$covariate.1 <- if_else(curve_info_zsg_vsg$covariate.1 == "Swiss_Standard", 0, 1)
curve_info_zsg_vsg$covariate.2 <- if_else(curve_info_zsg_vsg$covariate.2 == "long", 0, 1)

# get mean and range of data points per f0 trajectory included in the analysis
n_zsg_vsg = count(curve_info_zsg_vsg, n_long)
mean(n_zsg_vsg$n)
range(n_zsg_vsg$n)

# convert data.frame to data.table
curve_info_zsg_vsg <- as.data.table(curve_info_zsg_vsg)
curve_info_zsg_vsg = setDT(curve_info_zsg_vsg)

## unload add-on packages (necessary for all FLMM functions to work properly)
invisible(lapply(paste0("package:", names(sessionInfo()$otherPkgs)),
                 detach,
                 character.only = TRUE, unload = TRUE))

## model obstruent data
library(sparseFLMM) # v.0.4.2

# specify the covariates' form and the symmetric interaction matrix (identical for all three models below)
covariate = TRUE
num_covariates = 2
covariate_form = rep("by",2)
interaction = TRUE
which_interaction = matrix(c(FALSE, TRUE, TRUE, FALSE), byrow = TRUE, nrow = 2, ncol = 2)

## FLMM 4: MSG vs. ZSG
# ATTENTION: model computation takes approx. 20 minutes!
results_msg_zsg <- sparseFLMM(curve_info = curve_info_msg_zsg,
                      covariate = covariate,
                      num_covariates = num_covariates, covariate_form = covariate_form,
                      interaction = interaction, which_interaction = which_interaction,
                      bf_covs = c(5,5,5), m_covs = list(c(2,3),c(2,3),c(2,3)),
                      use_famm = TRUE)

## FLMM 5: MSG vs. VSG
# model computation takes approx. 15 minutes
results_msg_vsg <- sparseFLMM(curve_info = curve_info_msg_vsg,
                              covariate = covariate,
                              num_covariates = num_covariates, covariate_form = covariate_form,
                              interaction = interaction, which_interaction = which_interaction,
                              bf_covs = c(5,5,5), m_covs = list(c(2,3),c(2,3),c(2,3)),
                              use_famm = TRUE)

## FLMM 6: ZSG vs. VSG
# ATTENTION: model computation takes approx. 20 to 30 minutes!
results_zsg_vsg <- sparseFLMM(curve_info = curve_info_zsg_vsg,
                              covariate = covariate,
                              num_covariates = num_covariates, covariate_form = covariate_form,
                              interaction = interaction, which_interaction = which_interaction,
                              bf_covs = c(5,5,5), m_covs = list(c(2,3),c(2,3),c(2,3)),
                              use_famm = TRUE)

## plot obstruent results
# set plot device to two rows, two columns
par(mfrow = c(2,2))

# define scaling of the y-axis for summed effects plot (identical for Figures 7 to 9)
miny <- -1
maxy <- 1

# save each plot with jpeg() as JPEG, height = 960, width = 960

## Figure 7: MSG vs. ZSG
# define reference mean and plot grid
intercept_msg_zsg <- results_msg_zsg$fpc_famm_hat_tri_constr$intercept
y_mean_msg_zsg <- results_msg_zsg$fpc_famm_hat_tri_constr$famm_cb_mean$value
my_grid_msg_zsg <- results_msg_zsg$my_grid

## plot effect of covariate.1 (Variety)
y_cov1_msg_zsg <- results_msg_zsg$fpc_famm_hat_tri_constr$famm_cb_covariate.1$value
se_cov1_msg_zsg <- results_msg_zsg$fpc_famm_hat_tri_constr$famm_cb_covariate.1$se

plot(my_grid_msg_zsg, y_cov1_msg_zsg, t = "l", lwd = 3, xlab = "Normalized time",
     ylab = "", ylim = c(-1,1), cex.axis = 1.4,
     main = expression(paste("a. Variety ", f[1](t))), cex.main = 2.0, cex.lab = 1.6)
title(ylab = expression(paste(Delta, " Scaled f0")), line = 2.5, cex.lab = 1.6)

# add confidence bands
lines(x = my_grid_msg_zsg, y = (y_cov1_msg_zsg + 2*se_cov1_msg_zsg), lty = 2, lwd = 3)
lines(x = my_grid_msg_zsg, y = (y_cov1_msg_zsg - 2*se_cov1_msg_zsg), lty = 2, lwd = 3)

# add zero reference line
abline(h = 0, lty = 3)

## plot effect of covariate.2 (Syllable structure)
y_cov2_msg_zsg <- results_msg_zsg$fpc_famm_hat_tri_constr$famm_cb_covariate.2$value
se_cov2_msg_zsg <- results_msg_zsg$fpc_famm_hat_tri_constr$famm_cb_covariate.2$se

plot(my_grid_msg_zsg, y_cov2_msg_zsg, t = "l",lwd = 3, xlab = "Normalized time",
     ylab = "", ylim = c(-1,1), cex.axis = 1.4,
     main = expression(paste("b. Syllable structure ", f[2](t))), cex.main = 2.0, cex.lab = 1.6)
title(ylab = expression(paste(Delta, " Scaled f0")), line = 2.5, cex.lab = 1.6)

# add confidence bands
lines(x = my_grid_msg_zsg, y = (y_cov2_msg_zsg + 2*se_cov2_msg_zsg), lty = 2, lwd = 3)
lines(x = my_grid_msg_zsg, y = (y_cov2_msg_zsg - 2*se_cov2_msg_zsg), lty = 2, lwd = 3)

# add zero reference line
abline(h = 0, lty = 3)

## plot interaction effect of both covariates (Variety * Syllable structure)
y_inter_1_2_msg_zsg <- results_msg_zsg$fpc_famm_hat_tri_constr$famm_cb_inter_1_2$value
se_inter_1_2_msg_zsg <- results_msg_zsg$fpc_famm_hat_tri_constr$famm_cb_inter_1_2$se

plot(my_grid_msg_zsg, y_inter_1_2_msg_zsg, t = "l", lwd = 3, xlab = "Normalized time",
     ylab = "", ylim = c(-1,1), cex.axis = 1.4,
     main = expression(paste("c. Variety * Syllable structure ", f[3](t))),
     cex.main = 2.0, cex.lab = 1.6)
title(ylab = expression(paste(Delta, " Scaled f0")), line = 2.5, cex.lab = 1.6)

# add confidence bands
lines(x = my_grid_msg_zsg, y = (y_inter_1_2_msg_zsg + 2*se_inter_1_2_msg_zsg),
      lty = 2,lwd = 3)
lines(x = my_grid_msg_zsg, y = (y_inter_1_2_msg_zsg - 2*se_inter_1_2_msg_zsg),
      lty = 2,lwd = 3)

# add zero reference line
abline(h = 0, lty = 3)

## plot summed effects curves: effects of Variety and Syllable structure on scaled f0
# define summed effects
none_msg_zsg <- intercept_msg_zsg + y_mean_msg_zsg
Variety_msg_zsg <- intercept_msg_zsg + y_mean_msg_zsg + y_cov1_msg_zsg
V1.quantity_msg_zsg <- intercept_msg_zsg + y_mean_msg_zsg + y_cov2_msg_zsg
Variety_V1.quantity_msg_zsg <- intercept_msg_zsg + y_mean_msg_zsg + y_cov1_msg_zsg + y_cov2_msg_zsg + y_inter_1_2_msg_zsg

plot(my_grid_msg_zsg, none_msg_zsg, col = "black", lty = "solid",lwd = 4,t = "l",
     ylim = c(miny, maxy), ylab = "", xlab = "", xaxt="n", yaxt = "n")
par(new = TRUE)
plot(my_grid_msg_zsg, Variety_msg_zsg, col = "black", lty = "dashed", lwd = 4,t = "l",
     ylim = c(miny, maxy), ylab = "", xlab = "", xaxt = "n", yaxt = "n")
par(new = TRUE)
plot(my_grid_msg_zsg, V1.quantity_msg_zsg, col = "azure4", lty = "solid",lwd = 3, t = "l",
     ylim = c(miny, maxy), ylab = "", xlab = "", xaxt = "n", yaxt = "n")
par(new = TRUE)
plot(my_grid_msg_zsg, Variety_V1.quantity_msg_zsg, col = "azure4", lty = "dashed",lwd = 3,
     ylim = c(miny, maxy),
     ylab = "",
     xlab = "Normalized time",
     main = expression(paste("d. Vowels preceding obstruents")), t = "l", cex.axis = 1.4,
     cex.main = 2.0, cex.lab = 1.6)
title(ylab = "Summed effects: scaled f0", line = 2.5, cex.lab = 1.6)

# add zero reference line
abline(h = 0, lty = 3)

# add legend text
leg.text <- c(expression(paste("Munich Standard German, open/long (baseline)")),
              expression(paste("Zurich Standard German, open/long")),
              expression(paste("Munich Standard German, closed/short")),
              expression(paste("Zurich Standard German, closed/short")))
legend("bottomleft", legend = leg.text, cex = 1.3, bty = "n",
       col = c("black","black","azure4","azure4"),
       lty = c("solid","dashed","solid","dashed"), 
       lwd = c(2, 2, 1, 1))

## Figure 8: MSG vs. VSG
# define reference mean and plot grid
intercept_msg_vsg <- results_msg_vsg$fpc_famm_hat_tri_constr$intercept
y_mean_msg_vsg <- results_msg_vsg$fpc_famm_hat_tri_constr$famm_cb_mean$value
my_grid_msg_vsg <- results_msg_vsg$my_grid

## plot effect of covariate.1 (Variety)
y_cov1_msg_vsg <- results_msg_vsg$fpc_famm_hat_tri_constr$famm_cb_covariate.1$value
se_cov1_msg_vsg <- results_msg_vsg$fpc_famm_hat_tri_constr$famm_cb_covariate.1$se

plot(my_grid_msg_vsg, y_cov1_msg_vsg, t = "l", lwd = 3, xlab = "Normalized time",
     ylab = "", ylim = c(-1,1), cex.axis = 1.4,
     main = expression(paste("a. Variety ", f[1](t))), cex.main = 2.0, cex.lab = 1.6)
title(ylab = expression(paste(Delta, " Scaled f0")), line = 2.5, cex.lab = 1.6)

# add confidence bands
lines(x = my_grid_msg_vsg, y = (y_cov1_msg_vsg + 2*se_cov1_msg_vsg), lty = 2, lwd = 3)
lines(x = my_grid_msg_vsg, y = (y_cov1_msg_vsg - 2*se_cov1_msg_vsg), lty = 2, lwd = 3)

# add zero reference line
abline(h = 0, lty = 3)

## plot effect of covariate.2 (Syllable structure)
y_cov2_msg_vsg <- results_msg_vsg$fpc_famm_hat_tri_constr$famm_cb_covariate.2$value
se_cov2_msg_vsg <- results_msg_vsg$fpc_famm_hat_tri_constr$famm_cb_covariate.2$se

plot(my_grid_msg_vsg, y_cov2_msg_vsg, t = "l",lwd = 3, xlab = "Normalized time",
     ylab = "", ylim = c(-1,1.4), cex.axis = 1.4,
     main = expression(paste("b. Syllable structure ", f[2](t))), cex.main = 2.0, cex.lab = 1.6)
title(ylab = expression(paste(Delta, " Scaled f0")), line = 2.5, cex.lab = 1.6)

# add confidence bands
lines(x = my_grid_msg_vsg, y = (y_cov2_msg_vsg + 2*se_cov2_msg_vsg), lty = 2, lwd = 3)
lines(x = my_grid_msg_vsg, y = (y_cov2_msg_vsg - 2*se_cov2_msg_vsg), lty = 2, lwd = 3)

# add zero reference line
abline(h = 0, lty = 3)

## plot interaction effect of both covariates (Variety * Syllable structure)
y_inter_1_2_msg_vsg <- results_msg_vsg$fpc_famm_hat_tri_constr$famm_cb_inter_1_2$value
se_inter_1_2_msg_vsg <- results_msg_vsg$fpc_famm_hat_tri_constr$famm_cb_inter_1_2$se

plot(my_grid_msg_vsg, y_inter_1_2_msg_vsg, t = "l", lwd = 3, xlab = "Normalized time",
     ylab = "", ylim = c(-1,1), cex.axis = 1.4,
     main = expression(paste("c. Variety * Syllable structure ", f[3](t))),
     cex.main = 2.0, cex.lab = 1.6)
title(ylab = expression(paste(Delta, " Scaled f0")), line = 2.5, cex.lab = 1.6)

# add confidence bands
lines(x = my_grid_msg_vsg, y = (y_inter_1_2_msg_vsg + 2*se_inter_1_2_msg_vsg),
      lty = 2,lwd = 3)
lines(x = my_grid_msg_vsg, y = (y_inter_1_2_msg_vsg - 2*se_inter_1_2_msg_vsg),
      lty = 2,lwd = 3)

# add zero reference line
abline(h = 0, lty = 3)

## plot summed effects curves: effects of Variety and Syllable structure on scaled f0
# define summed effects
none_msg_vsg <- intercept_msg_vsg + y_mean_msg_vsg
Variety_msg_vsg <- intercept_msg_vsg + y_mean_msg_vsg + y_cov1_msg_vsg
V1.quantity_msg_vsg <- intercept_msg_vsg + y_mean_msg_vsg + y_cov2_msg_vsg
Variety_V1.quantity_msg_vsg <- intercept_msg_vsg + y_mean_msg_vsg + y_cov1_msg_vsg + y_cov2_msg_vsg + y_inter_1_2_msg_vsg

plot(my_grid_msg_vsg, none_msg_vsg, col = "black", lty = "solid",lwd = 4,t = "l",
     ylim = c(miny, maxy), ylab = "", xlab = "", xaxt="n", yaxt = "n")
par(new = TRUE)
plot(my_grid_msg_vsg, Variety_msg_vsg, col = "black", lty = "dashed", lwd = 4,t = "l",
     ylim = c(miny, maxy), ylab = "", xlab = "", xaxt = "n", yaxt = "n")
par(new = TRUE)
plot(my_grid_msg_vsg, V1.quantity_msg_vsg, col = "azure4", lty = "solid",lwd = 3, t = "l",
     ylim = c(miny, maxy), ylab = "", xlab = "", xaxt = "n", yaxt = "n")
par(new = TRUE)
plot(my_grid_msg_vsg, Variety_V1.quantity_msg_vsg, col = "azure4", lty = "dashed",lwd = 3,
     ylim = c(miny, maxy),
     ylab = "",
     xlab = "Normalized time",
     main = expression(paste("d. Vowels preceding obstruents")), t = "l", cex.axis = 1.4,
     cex.main = 2.0, cex.lab = 1.6)
title(ylab = "Summed effects: scaled f0", line = 2.5, cex.lab = 1.6)

# add zero reference line
abline(h = 0, lty = 3)

# add legend text
leg.text <- c(expression(paste("Munich Standard German, open/long (baseline)")),
              expression(paste("Vienna Standard German, open/long")),
              expression(paste("Munich Standard German, closed/short")),
              expression(paste("Vienna Standard German, closed/short")))
legend("bottomleft", legend = leg.text, cex = 1.3, bty = "n",
       col = c("black","black","azure4","azure4"),
       lty = c("solid","dotdash","solid","dotdash"), 
       lwd = c(2, 2, 1, 1))

## Figure 9: ZSG vs. VSG
# define reference mean and plot grid
intercept_zsg_vsg <- results_zsg_vsg$fpc_famm_hat_tri_constr$intercept
y_mean_zsg_vsg <- results_zsg_vsg$fpc_famm_hat_tri_constr$famm_cb_mean$value
my_grid_zsg_vsg <- results_zsg_vsg$my_grid

## plot effect of covariate.1 (Variety)
y_cov1_zsg_vsg <- results_zsg_vsg$fpc_famm_hat_tri_constr$famm_cb_covariate.1$value
se_cov1_zsg_vsg <- results_zsg_vsg$fpc_famm_hat_tri_constr$famm_cb_covariate.1$se

plot(my_grid_zsg_vsg, y_cov1_zsg_vsg, t = "l", lwd = 3, xlab = "Normalized time",
     ylab = "", ylim = c(-1,1), cex.axis = 1.4,
     main = expression(paste("a. Variety ", f[1](t))), cex.main = 2.0, cex.lab = 1.6)
title(ylab = expression(paste(Delta, " Scaled f0")), line = 2.5, cex.lab = 1.6)

# add confidence bands
lines(x = my_grid_zsg_vsg, y = (y_cov1_zsg_vsg + 2*se_cov1_zsg_vsg), lty = 2, lwd = 3)
lines(x = my_grid_zsg_vsg, y = (y_cov1_zsg_vsg - 2*se_cov1_zsg_vsg), lty = 2, lwd = 3)

# add zero reference line
abline(h = 0, lty = 3)

## plot effect of covariate.2 (Syllable structure)
y_cov2_zsg_vsg <- results_zsg_vsg$fpc_famm_hat_tri_constr$famm_cb_covariate.2$value
se_cov2_zsg_vsg <- results_zsg_vsg$fpc_famm_hat_tri_constr$famm_cb_covariate.2$se

plot(my_grid_zsg_vsg, y_cov2_zsg_vsg, t = "l",lwd = 3, xlab = "Normalized time",
     ylab = "", ylim = c(-1,1.6), cex.axis = 1.4,
     main = expression(paste("b. Syllable structure ", f[2](t))), cex.main = 2.0, cex.lab = 1.6)
title(ylab = expression(paste(Delta, " Scaled f0")), line = 2.5, cex.lab = 1.6)

# add confidence bands
lines(x = my_grid_zsg_vsg, y = (y_cov2_zsg_vsg + 2*se_cov2_zsg_vsg), lty = 2, lwd = 3)
lines(x = my_grid_zsg_vsg, y = (y_cov2_zsg_vsg - 2*se_cov2_zsg_vsg), lty = 2, lwd = 3)

# add zero reference line
abline(h = 0, lty = 3)

## plot interaction effect of both covariates (Variety * Syllable structure)
y_inter_1_2_zsg_vsg <- results_zsg_vsg$fpc_famm_hat_tri_constr$famm_cb_inter_1_2$value
se_inter_1_2_zsg_vsg <- results_zsg_vsg$fpc_famm_hat_tri_constr$famm_cb_inter_1_2$se

plot(my_grid_zsg_vsg, y_inter_1_2_zsg_vsg, t = "l", lwd = 3, xlab = "Normalized time",
     ylab = "", ylim = c(-1,1), cex.axis = 1.4,
     main = expression(paste("c. Variety * Syllable structure ", f[3](t))),
     cex.main = 2.0, cex.lab = 1.6)
title(ylab = expression(paste(Delta, " Scaled f0")), line = 2.5, cex.lab = 1.6)

# add confidence bands
lines(x = my_grid_zsg_vsg, y = (y_inter_1_2_zsg_vsg + 2*se_inter_1_2_zsg_vsg),
      lty = 2,lwd = 3)
lines(x = my_grid_zsg_vsg, y = (y_inter_1_2_zsg_vsg - 2*se_inter_1_2_zsg_vsg),
      lty = 2,lwd = 3)

# add zero reference line
abline(h = 0, lty = 3)

## plot summed effects curves: effects of Variety and Syllable structure on scaled f0
# define summed effects
none_zsg_vsg <- intercept_zsg_vsg + y_mean_zsg_vsg
Variety_zsg_vsg <- intercept_zsg_vsg + y_mean_zsg_vsg + y_cov1_zsg_vsg
V1.quantity_zsg_vsg <- intercept_zsg_vsg + y_mean_zsg_vsg + y_cov2_zsg_vsg
Variety_V1.quantity_zsg_vsg <- intercept_zsg_vsg + y_mean_zsg_vsg + y_cov1_zsg_vsg + y_cov2_zsg_vsg + y_inter_1_2_zsg_vsg

plot(my_grid_zsg_vsg, none_zsg_vsg, col = "black", lty = "solid",lwd = 4,t = "l",
     ylim = c(miny, maxy), ylab = "", xlab = "", xaxt="n", yaxt = "n")
par(new = TRUE)
plot(my_grid_zsg_vsg, Variety_zsg_vsg, col = "black", lty = "dashed", lwd = 4,t = "l",
     ylim = c(miny, maxy), ylab = "", xlab = "", xaxt = "n", yaxt = "n")
par(new = TRUE)
plot(my_grid_zsg_vsg, V1.quantity_zsg_vsg, col = "azure4", lty = "solid",lwd = 3, t = "l",
     ylim = c(miny, maxy), ylab = "", xlab = "", xaxt = "n", yaxt = "n")
par(new = TRUE)
plot(my_grid_zsg_vsg, Variety_V1.quantity_zsg_vsg, col = "azure4", lty = "dashed",lwd = 3,
     ylim = c(miny, maxy),
     ylab = "",
     xlab = "Normalized time",
     main = expression(paste("d. Vowels preceding obstruents")), t = "l", cex.axis = 1.4,
     cex.main = 2.0, cex.lab = 1.6)
title(ylab = "Summed effects: scaled f0", line = 2.5, cex.lab = 1.6)

# add zero reference line
abline(h = 0, lty = 3)

# add legend text
leg.text <- c(expression(paste("Zurich Standard German, open/long (baseline)")),
              expression(paste("Vienna Standard German, open/long")),
              expression(paste("Zurich Standard German, closed/short")),
              expression(paste("Vienna Standard German, closed/short")))
legend("bottomleft", legend = leg.text, cex = 1.3, bty = "n",
       col = c("black","black","azure4","azure4"),
       lty = c("solid","dotdash","solid","dotdash"), 
       lwd = c(2, 2, 1, 1))

# set plot device back to default: one row, one column
par(mfrow = c(1,1))
dev.off()