Chapter 12 声学分析项目案例
#install.packages("plotfunctions")
library(tidyverse)
# 统计建模
library(mgcv)
library(lme4)
library(car)
library(lsmeans)
library(rstatix)
require(nlme)
## 作图
library(itsadug)
# plot the comparison
library(cowplot)
# remove scientific notation
options(scipen=999)
# 计算置信区间
library(rcompanion)案例12.1 越南语南北方言声调曲拱对比
12.0.3 数据处理
## 数据导入及整理
contour_data_all = read.csv("data/ch12/nv.sv.data.2024-08-24.csv")
# modelling NV difference
nv.data = contour_data_all%>%
ungroup()%>%
filter(tone %in% c("NV415", "NV214") )%>%
#filter(subject != 214)%>%
rename(Time = Timepoint)%>%
select(-language)%>%
mutate(id = as.factor(id),
tone = as.factor(tone),
subject = as.factor(subject))
#sort data per individual trajectory (for autocorrelation)
nv.data = data.frame(nv.data)
nv.data$id = as.factor(nv.data$id)
nv.data = start_event(nv.data, event = c("subject", "id"))12.0.3.1 数据建模
### 越南语北部方言声调建模
# preliminary model without autocorrelation correction
model.nv = bam(normF0 ~ tone + s(Time, by=tone)+
s(Time, subject, by = tone, bs = "fs", m = 1),
data= nv.data)
summary(model.nv)
nvacf = acf_resid(model.nv)
#获取参数
nvacf[2]
## with autocorrelation correction
model.nv.b = bam(normF0 ~ tone + s(Time, by=tone)+
s(Time, subject, by = tone, bs = "fs", m = 1),
rho = nvacf[2],
AR.start = nv.data$start.event,
discrete = T,
nthreads = 2,
data= nv.data)
summary(model.nv.b)
# plot
plot_smooth(model.nv.b, ylab = "Normalised f0",
xlab = "Normalised time",
view="Time", cond=list(tone=c("NV415")),
ylim=c(-0.5,1), rm.ranef=T,
main="Northen Vietnamese",
rug=FALSE, col="blue",
hide.label = "fitted values")
plot_smooth(model.nv.b, view="Time",
cond=list(tone=c("NV214")),
v0 = 6, rug=FALSE,rm.ranef=T,
col="green",add=T)# 验证越南语北部方言是否存在这个声调对立
# （一）模型比较法
# 模型1区分两个对立
model.nv = bam(normF0 ~ tone + s(Time, by=tone)+
s(Time, subject, by = tone, bs = "fs", m = 1),
data= nv.data)
# 模型2不区分两个对立
model.nv0 = bam(normF0 ~ tone + s(Time)+
s(Time, subject, by = tone, bs = "fs", m = 1),
data= nv.data)
#比较两个模型
compareML(model.nv, model.nv0)
#此方法缺点是计算量比较大。但是对于不是很复杂的模型没有影响。# （二）差值比较法
# 通过修改模型的设定，将表示两个原始平滑函数之间差异的函数包含在模型中。接下来，发现这个差异平滑函数是显著的，这就表明有显著差异。为了拟合这个新模型，我们首先需要创建一个新的二元变量，该变量在一个水平上等于0，而在另一个水平上等于1。我们现在创建一个名为IS214的变量，该变量在单词‘NV214’上为1，而在单词‘NV415’上为0。
nv.data$IS214 = (nv.data$tone == "NV214")*1
model.nv.binary = bam(normF0 ~ s(Time) + s(Time, by=IS214)+
s(Time, subject, by = tone, bs = "fs", m = 1),
rho = nvacf[2],
AR.start = nv.data$start.event,
discrete = T,
nthreads = 2,
data= nv.data)
summary(model.nv.binary)# （三）
nv.data$toneO = as.ordered(nv.data$tone)
contrasts(nv.data$toneO) = "contr.treatment"
model.nv.ord = bam(normF0 ~ toneO+ s(Time) + s(Time, by=toneO)+
s(Time, subject, by = tone, bs = "fs", m = 1),
rho = nvacf[2],
AR.start = nv.data$start.event,
discrete = T,
nthreads = 2,
data= nv.data)
summary(model.nv.ord )### 越南语南部方言声调建模
sv.data = contour_data_all%>%
ungroup()%>%
filter(tone %in% c("SV415", "SV214") )%>%
rename(Time = Timepoint)%>%
select(-language)%>%
mutate(id = as.factor(id),
tone = as.factor(tone),
subject = as.factor(subject))
#sort data per individual trajectory (for autocorrelation)
sv.data = start_event(sv.data, event = c("subject", "id"))
# preliminary model without autocorrelation correction
model.sv = bam(normF0 ~ tone + s(Time, by=tone)+
s(Time, subject, by = tone, bs = "fs", m = 1),
data= sv.data)
summary(model.sv)
svacf = acf_resid(model.sv)
svacf[2]
## with autocorrelation correction
model.sv.b = bam(normF0 ~ tone + s(Time, by=tone)+
s(Time, subject, by = tone, bs = "fs", m = 1),
rho = svacf[2],
AR.start = sv.data$start.event,
discrete = T,
nthreads = 2,
data= sv.data)
summary(model.sv.b)
# plot
plot_smooth(model.sv.b, ylab = "Normalised f0",
xlab = "Normalised time",
view="Time", cond=list(tone=c("SV415")),
ylim=c(-0.5,1), rm.ranef=T,
main="Southen Vietnamese",
rug=FALSE, col="blue",
hide.label = "fitted values")
plot_smooth(model.sv.b, view="Time",
cond=list(tone=c("SV214")),
v0 = 6, rug=FALSE,rm.ranef=T,
col="green",add=T)
model.sv = bam(normF0 ~ tone + s(Time, by=tone)+
s(Time, subject, by = tone, bs = "fs", m = 1),
data= sv.data)
model.sv0 = bam(normF0 ~ tone + s(Time)+
s(Time, subject, by = tone, bs = "fs", m = 1),
data= sv.data)
compareML(model.sv.b, model.sv.b0)
# solution 2
sv.data$IS214 = (sv.data$tone == "SV214")*1
model.sv.binary = bam(normF0 ~ s(Time) + s(Time, by=IS214)+
s(Time, subject, by = tone, bs = "fs", m = 1),
rho = svacf[2],
AR.start = sv.data$start.event,
discrete = T,
nthreads = 2,
data= sv.data)
summary(model.sv.binary)
# soluation 3
sv.data$toneO = as.ordered(sv.data$tone)
contrasts(sv.data$toneO) = "contr.treatment"
model.sv.ord = bam(normF0 ~ toneO+ s(Time) + s(Time, by=toneO)+
s(Time, subject, by = tone, bs = "fs", m = 1),
rho = svacf[2],
AR.start = sv.data$start.event,
discrete = T,
nthreads = 2,
data= sv.data)
summary(model.sv.ord )# png(paste(dir.fig, paste("gam.mdl",Sys.Date(),"png",sep = "."),sep = ""), units="in", width=14, height=10, res=600)
par(mfrow=c(1,2))
# Northen Vietnamese
plot_smooth(model.nv.b, ylab = "Normalised f0",
xlab = "Normalised time",
view="Time", cond=list(tone=c("NV415")),
ylim=c(-0.5,0.6), rm.ranef=T,
main="Northern Vietnamese",
rug=FALSE, col="black",
hide.label = "fitted values")
plot_smooth(model.nv.b, view="Time",
cond=list(tone=c("NV214")),
v0 = 6, rug=FALSE,rm.ranef=T,
col="gray",add=T)
# Southern vietnamese
plot_smooth(model.sv.b, ylab = "Normalised f0",
xlab = "Normalised time",
view="Time", cond=list(tone=c("SV415")),
ylim=c(-0.5,0.6), rm.ranef=T,
main="Southern Vietnamese",
rug=FALSE, col="black",
hide.label = "fitted values")
plot_smooth(model.sv.b, view="Time",
cond=list(tone=c("SV214")),
v0 = 6, rug=FALSE,rm.ranef=T,
col="gray",add=T)
# dev.off()
# write.csv(nv.data,
# file = paste("data/ch12/",
# paste("nv.sv.data",
# Sys.Date(),
# "csv",
# sep = "."),
# sep = "") ,row.names = FALSE)案例12.2 汉语母语者模仿泰语声调研究
详细背景请看原文献
研究方法
数据处理
##### import eprime data
ep_md <- read.delim("data/ch12/md_ep.txt")
#data cleaning
ep_md_clean <- ep_md%>%
filter(Procedure.Block. != "prac")%>%
mutate(isi = dplyr::recode(ExperimentName,
imit_test = "short",
imit_test2000ms = "long"),
tone.Trial. = as.factor(tone.Trial.))%>%
select(isi, subject = Subject, session = Session,
block = Procedure.Block., informant = speaker.Trial.,
syllable = syllable.Trial., token = token.Trial.,
tone = tone.Trial., order = Block)%>%
mutate(cgload = dplyr::recode(block, block1 = "sssv", block2 = "sssv",
block3 = "ssdv", block4 = "ssdv",
block5 = "dssv", block6 = "dssv",
block7 = "dsdv", block8 = "dsdv"))%>%
mutate(speaker = substr(cgload, 1, 2),
vowel_vari = substr(cgload, 3, 4))%>%
unite(col = "bundle", c(subject, session), sep ="_", remove = FALSE)%>%
group_by(bundle)%>%
mutate(id = 1:n(),
language = "md")
##################
##Imitation data##
##################
# ProsodyPro outputs acoustic measures in the **means** files.
# We only used duration in these files that measured the duration of the syllable.
# For other F0 related measures, we calculated from the pitch contour (10%-90%).
#duration
#set working directory
duration_md_file = list.files("data/ch12/md_mean",full.names=T)
# duration data
#mandarin data
df_md_dur = data.frame()
bin = data.frame()
for (i in 1:length(duration_md_file)){
bin = read.delim(duration_md_file[i])
bin = mutate(bin, bundle = str_extract(duration_md_file[i], "[0-9]{3}_[1-5]") )
#bin = mutate(bin, id = row_number())
df_md_dur = rbind(df_md_dur, bin)
}
dur_md = df_md_dur%>%
group_by(bundle)%>%
mutate(id = 1: n())%>% # add id for each recording
left_join(., ep_md_clean, by = c("id","bundle"))%>%
filter(., rowLabel != "na")
#missing data
missing_md = df_md_dur%>%
group_by(bundle)%>%
mutate(id = 1: n())%>% # add id for each recording
left_join(., ep_md_clean, by = c("id","bundle"))%>%
filter(., rowLabel == "na")%>%
mutate(count = 1)%>%
group_by(subject, syllable, tone, speaker, vowel_vari)%>%
summarise(n = sum(count))
# write.csv(missing_md,
# file = paste(dir.processed,
# paste("missing_md",
# Sys.Date(),
# "csv",
# sep = "."),
# sep = "") ,row.names = FALSE)
#
#loading 10 point f0 data
filename_md = list.files("data/ch12/md_norm_f0",full.names=T)
#mandarin
df_md_f0 = data.frame()
bin = data.frame()
for (i in 1:length(filename_md)){
bin = read.delim(filename_md[i])
bin = mutate(bin, bundle = str_extract(filename_md[i], "[0-9]{3}_[1-5]") )
#bin = mutate(bin, id = row_number())
df_md_f0 = rbind(df_md_f0, bin)
}
#attaching e-prime data
#mandarin data
f0_md = df_md_f0%>%
group_by(bundle)%>%
mutate(id = rep(c(1:80), each = 10))%>% #segement 80 tokens per session per speaker
mutate(temp = as.factor(cut(NormalizedTime, 80)))%>%
group_by(bundle, temp)%>%
mutate(timepoint = 1:10)%>%
filter(., rowLabel != "na")%>%
left_join(., ep_md_clean, by = c("id","bundle"))
# write.csv(f0_md,
# file = paste(dir.processed,
# paste("f0.md.10pnt",
# Sys.Date(),
# "csv",
# sep = "."),
# sep = "") ,row.names = FALSE)
#
# Labonov normalisation
# 8 points selection
md_imit_8points = f0_md %>%
group_by(subject)%>%
mutate(F0norm = (F0 - mean(F0))/sqrt(sum(F0^2)/length(F0)),
language = "md")%>%
filter(timepoint>1 & timepoint < 10)%>%
group_by(bundle, id)%>%
mutate(timepoint = 1:n())%>% # recode timepoint
select("language","subject","F0","isi","id","tone","speaker",
"vowel_vari","syllable", "token",
"timepoint","F0norm","bundle","cgload","informant")%>%
ungroup()
md_imit_discrete = md_imit_8points %>%
group_by(bundle, id, language, isi,
speaker, tone, vowel_vari,syllable, token, cgload,informant,subject)%>%
summarise(excursion = max(F0norm)-min(F0norm),
meanf0 = mean(F0norm),
maxloc = which.max(F0norm)/8)%>%
ungroup()
# write.csv(md_imit_8points,
# file = paste(dir.processed,
# paste("md_imit_8points",
# Sys.Date(),
# "csv",
# sep = "."),
# sep = "") ,row.names = FALSE)
# md_all = left_join(md, imit_excursion, by = c("bundle", "id"))
#####################
##stimuli data#######
#####################
#stimulus syllabe duration
stim.duration <- read.delim("data/ch12/stim_duration.txt", na.strings="--undefined--", stringsAsFactors=FALSE)
colnames(stim.duration)= c("rowLabel", "stim_duration")
stim.duration = stim.duration%>%
mutate(tone = str_extract(rowLabel, "[0-9]{2,3}"),
informant = str_extract(rowLabel, "F[1-4]"),
token = str_extract(rowLabel, "[1-6]{1}$"),
syllable = str_extract(rowLabel, "m[a,i]"),
#normalize duration
dur_stim_norm = (stim_duration - mean(stim_duration))/sqrt(sum(stim_duration^2)/length(stim_duration)))
#write.csv(table.data, file = "table1.csv")
#Thai 10 point
raw.10points = read.delim("data/ch12/stim_normf0.txt")
#get excursion
stim.8point = raw.10points%>%
mutate(., tone = str_extract(Filename, "[0-9]{2,3}"),
informant = str_extract(Filename, "F[1-4]"),
token = str_extract(Filename, "[1-6]{1}$"),
syllable = str_extract(Filename, "m[a-z]"))%>%
gather(timepoint, f0, 2:11)%>%
group_by(informant )%>%
mutate(F0norm.stim = (f0- mean(f0))/sqrt(sum(f0^2)/length(f0)))%>%
filter(timepoint%in%c("X2","X3","X4","X5","X6","X7","X8","X9"))%>%
group_by(informant, tone, syllable,token)%>%
mutate(timepoint = 1:n())%>%
select("f0","tone","informant","syllable","token","timepoint","F0norm.stim")
write.csv(stim.8point,
file = paste(dir.processed,
paste("stim.8point",
Sys.Date(),
"csv",
sep = "."),
sep = "") ,row.names = FALSE)
stim.discret = stim.8point%>%
group_by(informant, tone, token, syllable)%>%
summarise(mean.stim=mean(F0norm.stim),
excursion.stim = max(F0norm.stim)-min(F0norm.stim),
maxloc.stim = which.max(F0norm.stim)/8)
stimuli.data.all = stim.discret%>%
left_join(., stim.duration, by = c("informant", "tone", "token", "syllable"))%>%
select(informant,
tone,
token, syllable,
stim_duration,
mean.stim,
excursion.stim,
maxloc.stim,
dur_stim_norm)%>%
ungroup(informant)
####################
# difference scores#
####################
# duration difference score
ds.md.dur = dur_md %>%
mutate(., informant = as.character(informant),
syllable = str_extract(syllable, "m[a-z]{1}"),
syllable = as.character(syllable),
token = as.character(token),
tone = as.character(tone))%>%
group_by(subject)%>%
# normalize duration
mutate(dur_norm = (duration- mean( duration))/sqrt(sum( duration^2)/length( duration)))%>%
ungroup()%>%
#join stimulus data
left_join(., stimuli.data.all,
by = c("informant","syllable","token","tone"))%>%
mutate(df_dur= dur_norm - dur_stim_norm)
# write.csv(ds_md_dur, file = "data/paper3/processed/ds_md_dur.csv", row.names = FALSE)
write.csv(ds_md_dur,
file = paste(dir.processed,
paste("ds.md.dur",
Sys.Date(),
"csv",
sep = "."),
sep = "") ,row.names = FALSE)
# discrete difference score
ds.md.discrete = md_imit_discrete %>%
mutate(., informant = as.character(informant),
syllable = str_extract(syllable, "m[a-z]{1}"),
syllable = as.character(syllable),
token = as.character(token),
tone = as.character(tone))%>%
#join stimulus data
left_join(., stimuli.data.all,
by = c("informant","syllable","token","tone"))%>%
mutate(., df_meanf0 = meanf0 - mean.stim,
df_excursion = excursion - excursion.stim,
df_maxloc = maxloc - maxloc.stim)
# # for plot difference score contour
# ds_md_8point = md_imit_8points %>%
# mutate(., informant = as.character(informant),
# syllable = str_extract(syllable, "m[a-z]{1}"),
# syllable = as.character(syllable),
# token = as.character(token),
# tone = as.character(tone))%>%
# #join stimulus data
# left_join(., stim.8point,
# by = c("informant","syllable","token","tone","timepoint"))%>%
# mutate(., df_meanf0 = F0norm - F0norm.stim)
#
#
# save the data
write.csv(ds_md_discrete,
file = paste(dir.processed,
paste("ds.md.discrete",
Sys.Date(),
"csv",
sep = "."),
sep = "") ,row.names = FALSE)研究结果
数据建模
#import data
# this is the f0-related measures
ds.md.discrete <- read.csv(paste(dir.processed,
paste("ds.md.discrete",
"2020-08-29",
"csv",
sep = "."),
sep = ""))
# this is for duration data
ds.md.dur <- read.csv(paste(dir.processed,
paste("ds.md.dur",
"2020-08-29",
"csv",
sep = "."),
sep = ""))
#### Experiemnt 1 ####
# convert tone into a categorical variable
ds.md.dur$tone = as.factor(ds.md.dur$tone)
ds.md.dur$subject = as.factor(ds.md.dur$subject)
ds.md.discrete$tone = as.factor(ds.md.discrete$tone)
ds.md.discrete$subject = as.factor(ds.md.discrete$subject)
#### md duration modelling #####
imit.md.dur = lmer(df_dur ~ isi*speaker*vowel_vari* tone +(1|subject),
data = ds.md.dur)
summary(imit.md.dur)
imit.md.dur.smry = data.frame(coef(summary(imit.md.dur)))
imit.md.dur.smry$variable = rownames(imit.md.dur.smry)
write.csv(imit.md.dur.smry,
file = paste(dir.tbl,
paste("imit.md.dur.smry",
Sys.Date(),
"csv",
sep = "."),
sep = "") ,row.names = FALSE)
appendix0a.new = Anova(imit.md.dur, test = "F")
emm_options(pbkrtest.limit = 5100) # set this larger than the data points
imit.md.dur.post1 = lsmeans(imit.md.dur, pairwise ~ tone)
imit.md.dur.post2 = lsmeans(imit.md.dur, pairwise ~ tone:speaker)
# imit.md.dur.post3 = lsmeans(imit.md.dur, pairwise ~ tone:isi)
# calculating effect size
# imit.md.dur.post1.eff = data.frame (eff_size(imit.md.dur.post1,
# sigma = sigma(imit.md.dur), edf = Inf))
# imit.md.dur.post2.eff = data.frame (eff_size(imit.md.dur.post2,
# sigma = sigma(imit.md.dur), edf = Inf))
imit.md.dur.post1.fnl = data.frame(imit.md.dur.post1$contrasts)%>%
bind_cols(imit.md.dur.post1.eff)%>%
mutate(type = "duration",
contrast = contrast...1,
contrast2 = contrast...7,
d = round(effect.size, 3),
estimate = round(estimate, 3),
SE = round(SE...3, 3),
df = round(df...4, 0),
t = round(t.ratio, 2),
p = round(p.value,3))%>%
select(type, contrast, contrast2, d, estimate, SE, df, t, p)
imit.md.dur.post2.fnl = data.frame(imit.md.dur.post2$contrasts)%>%
bind_cols(imit.md.dur.post2.eff)%>%
mutate(type = "duration",
contrast = contrast...1,
contrast2 = contrast...7,
d = round(effect.size, 3),
estimate = round(estimate, 3),
SE = round(SE...3, 3),
df = round(df...4, 0),
t = round(t.ratio, 2),
p = round(p.value,3))%>%
select(type, contrast, contrast2, d, estimate, SE, df,t, p)
write.csv(appendix0a,
file = "data/paper3/tables/appendix0a.csv", row.names = FALSE)
write.csv(imit.md.dur.post1$contrasts,
file = "data/paper3/tables/ap_dur.a.csv", row.names = FALSE)
write.csv(imit.md.dur.post2$contrasts,
file = "data/paper3/tables/ap_dur.b.csv", row.names = FALSE)
########## md f0 mean modelling ################
# ds_md_discrete$df_meanf0.abs = abs(ds_md_discrete$df_meanf0)
imit.md.f0mean = lmer(df_meanf0 ~ isi*speaker*vowel_vari* tone +(1|subject),
data = ds.md.discrete)
# imit.md.f0mean = lmer(df_meanf0 ~ isi*speaker*vowel_vari* tone +
# (1 + speaker|subject)+(1 + vowel_vari|subject),
# data = ds.md.discrete)
summary(imit.md.f0mean)
imit.md.f0mean.smry = data.frame(coef(summary(imit.md.f0mean)))
imit.md.f0mean.smry$variable = rownames(imit.md.f0mean.smry)
write.csv(imit.md.f0mean.smry,
file = paste(dir.tbl,
paste("imit.md.f0mean.smry",
Sys.Date(),
"csv",
sep = "."),
sep = "") ,row.names = FALSE)
#Caculating P for fixed factors
appendix1a = Anova(imit.md.f0mean, test = "F")
emm_options(pbkrtest.limit = 5100) # set this larger than your data points
imit.md.f0mean.post1 = lsmeans(imit.md.f0mean, pairwise ~ tone)
imit.md.f0mean.post2 = lsmeans(imit.md.f0mean, pairwise ~ tone:speaker)
imit.md.f0mean.post3 = lsmeans(imit.md.f0mean, pairwise ~ tone:isi)
ds.md.241 = ds.md.discrete%>%
filter(tone == 241)%>%
select(speaker, df_meanf0)
ds.md.241 = ds.all.discrete%>%
#filter(tone == 241)%>%
groupwiseMean(df_meanf0 ~ language + tone + isi + speaker,
data = .,
conf = 0.95,
digits = 3)
ds.all.discrete%>%
groupwiseMean(df_meanf0 ~ language + speaker,
data = .,
conf = 0.95,
digits = 3)
## effect size
# imit.md.f0mean.post1.eff = data.frame (eff_size(imit.md.f0mean.post1,
# sigma = sigma(imit.md.f0mean), edf = Inf))
# imit.md.f0mean.post2.eff = data.frame (eff_size(imit.md.f0mean.post2,
# sigma = sigma(imit.md.f0mean), edf = Inf))
# imit.md.f0mean.post3.eff = data.frame (eff_size(imit.md.f0mean.post3,
# sigma = sigma(imit.md.f0mean), edf = Inf))
#
# imit.md.f0mean.post1.fnl = data.frame(imit.md.f0mean.post1$contrasts)%>%
# bind_cols(imit.md.f0mean.post1.eff)%>%
# mutate(type = "f0mean",
# contrast = contrast...1,
# contrast2 = contrast...7,
# d = round(effect.size, 3),
# estimate = round(estimate, 3),
# SE = round(SE...3, 3),
# df = round(df...4, 0),
# t = round(t.ratio, 2),
# p = round(p.value,3))%>%
# select(type, contrast, contrast2, d, estimate, SE, df,t, p)
#
# imit.md.f0mean.post2.fnl = data.frame(imit.md.f0mean.post2$contrasts)%>%
# bind_cols(imit.md.f0mean.post2.eff)%>%
# mutate(type = "f0mean",
# contrast = contrast...1,
# contrast2 = contrast...7,
# d = round(effect.size, 3),
# estimate = round(estimate, 3),
# SE = round(SE...3, 3),
# df = round(df...4, 0),
# t = round(t.ratio, 2),
# p = round(p.value,3))%>%
# select(type, contrast, contrast2, d, estimate, SE, df, t, p)
#
# imit.md.f0mean.post3.fnl = data.frame(imit.md.f0mean.post3$contrasts)%>%
# bind_cols(imit.md.f0mean.post3.eff)%>%
# mutate(type = "f0mean",
# contrast = contrast...1,
# contrast2 = contrast...7,
# d = round(effect.size, 3),
# estimate = round(estimate, 3),
# SE = round(SE...3, 3),
# df = round(df...4, 0),
# t = round(t.ratio, 2),
# p = round(p.value,3))%>%
# select(type, contrast, contrast2, d, estimate, SE, df,t, p)
#
# write.csv(imit.md.f0mean.post1$contrasts, file = "data/paper3/tables/ap_meanf0.a.csv", row.names = FALSE)
# write.csv(imit.md.f0mean.post2$contrasts, file = "data/paper3/tables/ap_meanf0.b.csv", row.names = FALSE)
# write.csv(imit.md.f0mean.post3$contrasts, file = "data/paper3/tables/ap_meanf0.c.csv", row.names = FALSE)
# absolute value of deviation scores
# ds.md.discrete$df_meanf0.abs = abs(ds.md.discrete$df_meanf0)
#
# imit.md.f0mean.abs = lmer(df_meanf0.abs ~ isi*speaker*vowel_vari* tone +(1|subject),
# data = ds.md.discrete)
#
# #Caculating P for fixed factors
#
# appendix1a.abs = Anova(imit.md.f0mean.abs, test = "F")
##### md f0 excursion modeling #####
imit.md.excursion = lmer(df_excursion ~ isi*speaker*vowel_vari* tone +(1|subject),
data = ds.md.discrete)
summary(imit.md.excursion)
imit.md.excursion.smry = data.frame(coef(summary(imit.md.excursion)))
imit.md.excursion.smry$variable = rownames(imit.md.excursion.smry)
write.csv(imit.md.excursion.smry,
file = paste(dir.tbl,
paste("imit.md.excursion.smry",
Sys.Date(),
"csv",
sep = "."),
sep = "") ,row.names = FALSE)
appendix1b = Anova(imit.md.excursion, test = "F")
imit.md.excursion.post1 = lsmeans(imit.md.excursion, pairwise ~ tone)
imit.md.excursion.post2 = lsmeans(imit.md.excursion, pairwise ~ tone:speaker)
imit.md.excursion.post3 = lsmeans(imit.md.excursion, pairwise ~ tone:isi)
## effect size
imit.md.excursion.post1.eff = data.frame (eff_size(imit.md.excursion.post1,
sigma = sigma(imit.md.excursion), edf = Inf))
imit.md.excursion.post2.eff = data.frame (eff_size(imit.md.excursion.post2,
sigma = sigma(imit.md.excursion), edf = Inf))
imit.md.excursion.post3.eff = data.frame (eff_size(imit.md.excursion.post3,
sigma = sigma(imit.md.excursion), edf = Inf))
imit.md.excursion.post1.fnl = data.frame(imit.md.excursion.post1$contrasts)%>%
bind_cols(imit.md.excursion.post1.eff)%>%
mutate(type = "excursion",
contrast = contrast...1,
contrast2 = contrast...7,
d = round(effect.size, 3),
estimate = round(estimate, 3),
SE = round(SE...3, 3),
df = round(df...4, 0),
t = round(t.ratio, 2),
p = round(p.value,3))%>%
select(type, contrast, contrast2, d, estimate, SE, df, t,p)
imit.md.excursion.post2.fnl = data.frame(imit.md.excursion.post2$contrasts)%>%
bind_cols(imit.md.excursion.post2.eff)%>%
mutate(type = "excursion",
contrast = contrast...1,
contrast2 = contrast...7,
d = round(effect.size, 3),
estimate = round(estimate, 3),
SE = round(SE...3, 3),
df = round(df...4, 0),
t = round(t.ratio, 2),
p = round(p.value,3))%>%
select(type, contrast, contrast2, d, estimate, SE, df, t, p)
imit.md.excursion.post3.fnl = data.frame(imit.md.excursion.post3$contrasts)%>%
bind_cols(imit.md.excursion.post3.eff)%>%
mutate(type = "excursion",
contrast = contrast...1,
contrast2 = contrast...7,
d = round(effect.size, 3),
estimate = round(estimate, 3),
SE = round(SE...3, 3),
df = round(df...4, 0),
t = round(t.ratio, 2),
p = round(p.value,3))%>%
select(type, contrast, contrast2, d, estimate, SE, df, t, p)
write.csv(imit.md.excursion.post1$contrasts,
file = "data/paper3/tables/ap_excursion.a.csv", row.names = FALSE)
write.csv(imit.md.excursion.post2$contrasts,
file = "data/paper3/tables/ap_excursion.b.csv", row.names = FALSE)
write.csv(imit.md.excursion.post3$contrasts,
file = "data/paper3/tables/ap_excursion.c.csv", row.names = FALSE)
ds.md.discrete$df_excursion.abs = abs(ds.md.discrete$df_excursion)
imit.md.excursion.abs = lmer(df_excursion.abs ~ isi*speaker*vowel_vari* tone +(1|subject),
data = ds.md.discrete)
#Caculating P for fixed factors
#appendix1b.abs = Anova(imit.md.excursion.abs, test = "F")
##### md f0 maxloc modeling ####
imit.md.maxloc = lmer(df_maxloc ~ isi*speaker*vowel_vari* tone +(1|subject),
data = ds.md.discrete)
summary(imit.md.maxloc)
imit.md.maxloc.smry = data.frame(coef(summary(imit.md.maxloc)))
imit.md.maxloc.smry$variable = rownames(imit.md.maxloc.smry)
write.csv(imit.md.maxloc.smry,
file = paste(dir.tbl,
paste("imit.md.maxloc.smry",
Sys.Date(),
"csv",
sep = "."),
sep = "") ,row.names = FALSE)
appendix1c = Anova(imit.md.maxloc, test = "F")
imit.md.maxloc.post1 = lsmeans(imit.md.maxloc, pairwise ~ tone)
imit.md.maxloc.post2 = lsmeans(imit.md.maxloc, pairwise ~ tone:speaker)
imit.md.maxloc.post3 = lsmeans(imit.md.maxloc, pairwise ~ tone:isi)
lsmeans(imit.md.maxloc, pairwise ~ speaker)
## effect size
# imit.md.maxloc.post1.eff = data.frame (eff_size(imit.md.maxloc.post1,
# sigma = sigma(imit.md.maxloc), edf = Inf))
# imit.md.maxloc.post2.eff = data.frame (eff_size(imit.md.maxloc.post2,
# sigma = sigma(imit.md.maxloc), edf = Inf))
# imit.md.maxloc.post3.eff = data.frame (eff_size(imit.md.maxloc.post3,
# sigma = sigma(imit.md.maxloc), edf = Inf))
imit.md.maxloc.post1.fnl = data.frame(imit.md.maxloc.post1$contrasts)%>%
bind_cols(imit.md.maxloc.post1.eff)%>%
mutate(type = "maxloc",
contrast = contrast...1,
contrast2 = contrast...7,
d = round(effect.size, 3),
estimate = round(estimate, 3),
SE = round(SE...3, 3),
df = round(df...4, 0),
t = round(t.ratio, 2),
p = round(p.value,3))%>%
select(type, contrast, contrast2, d, estimate, SE, df,t, p)
imit.md.maxloc.post2.fnl = data.frame(imit.md.maxloc.post2$contrasts)%>%
bind_cols(imit.md.maxloc.post2.eff)%>%
mutate(type = "maxloc",
contrast = contrast...1,
contrast2 = contrast...7,
d = round(effect.size, 3),
estimate = round(estimate, 3),
SE = round(SE...3, 3),
df = round(df...4, 0),
t = round(t.ratio, 2),
p = round(p.value,3))%>%
select(type, contrast, contrast2, d, estimate, SE, df,t, p)
imit.md.maxloc.post3.fnl = data.frame(imit.md.maxloc.post3$contrasts)%>%
bind_cols(imit.md.maxloc.post3.eff)%>%
mutate(type = "maxloc",
contrast = contrast...1,
contrast2 = contrast...7,
d = round(effect.size, 3),
estimate = round(estimate, 3),
SE = round(SE...3, 3),
df = round(df...4, 0),
t = round(t.ratio, 2),
p = round(p.value,3))%>%
select(type, contrast, contrast2, d, estimate, SE, df,t, p)
write.csv(imit.md.maxloc.post1$contrasts, file = "data/paper3/tables/ap_maxloc.a.csv", row.names = FALSE)
write.csv(imit.md.maxloc.post2$contrasts, file = "data/paper3/tables/ap_maxloc.b.csv", row.names = FALSE)
write.csv(imit.md.maxloc.post3$contrasts, file = "data/paper3/tables/ap_maxloc.c.csv", row.names = FALSE)
ds.md.discrete$df_maxloc.abs = abs(ds.md.discrete$df_maxloc)
imit.md.maxloc.abs = lmer(df_maxloc.abs ~ isi*speaker*vowel_vari* tone +(1|subject),
data = ds.md.discrete)
#Caculating P for fixed factors
# appendix1c.abs = Anova(imit.md.maxloc.abs, test = "F")
# tone main effect
#imit.md.dur.post1 = data.frame(imit.md.dur.post1$contrasts)
#imit.md.f0mean.post1 = data.frame(imit.md.f0mean.post1$contrasts)
#imit.md.excursion.post1 = data.frame(imit.md.excursion.post1$contrasts)
#imit.md.maxloc.post1 = data.frame(imit.md.maxloc.post1$contrasts)
md_tone_main = rbind(imit.md.dur.post1.fnl,
imit.md.f0mean.post1.fnl,
imit.md.excursion.post1.fnl,
imit.md.maxloc.post1.fnl)
write.csv(md_tone_main,
file = paste(dir.tbl,
paste("md_tone_main",
Sys.Date(),
"csv",
sep = "."),
sep = "") ,row.names = FALSE)
# tone by talker interaction
md.tone.talker = rbind(imit.md.dur.post2.fnl,
imit.md.f0mean.post2.fnl,
imit.md.excursion.post2.fnl,
imit.md.maxloc.post2.fnl)%>%
mutate(t1 = str_extract(contrast, "[0-9]{2,3}"),
t2 = str_extract(contrast, "- [0-9]{2,3}"),
t2 = str_extract(t2, "[0-9]{2,3}"))%>%
filter(t1 == t2)%>%
select(-t1, -t2)
write.csv(md.tone.talker,
file = paste(dir.tbl,
paste("md.tone.talker",
Sys.Date(),
"csv",
sep = "."),
sep = "") ,row.names = FALSE)
# tone by memory load interaction
md.tone.isi = rbind( imit.md.f0mean.post3.fnl,
imit.md.excursion.post3.fnl,
imit.md.maxloc.post3.fnl)%>%
mutate(t1 = str_extract(contrast, "[0-9]{2,3}"),
t2 = str_extract(contrast, "- [0-9]{2,3}"),
t2 = str_extract(t2, "[0-9]{2,3}"))%>%
filter(t1 == t2)%>%
select(-t1, -t2)
write.csv(md.tone.isi,
file = paste(dir.tbl,
paste("md.tone.isi",
Sys.Date(),
"csv",
sep = "."),
sep = "") ,row.names = FALSE)
##### batch analysis #######
appendix1a.old = Anova(imit.md.f0mean, test = "F")
appendix1b.old = Anova(imit.md.excursion, test = "F")
appendix1c.old = Anova(imit.md.maxloc, test = "F")
appendix1a.new = Anova(imit.md.f0mean, test = "F")
appendix1b.new = Anova(imit.md.excursion, test = "F")
appendix1c.new = Anova(imit.md.maxloc, test = "F")
appendix0a.new = Anova(imit.md.dur, test = "F")
md_model.p3 = cbind(appendix1a,appendix1b,appendix1c)
# write.csv(md_model, file = "data/paper3/tables/md_model.csv")12.1 课程项目孵化器
声学分析可以使用IPA开放数据 UCLA 语音学实验室档案 http://archive.phonetics.ucla.edu/ UCLA 语音学实验室档案。收录了世界各地的语言录音样本。这些数据集对于语音研究、音系学研究，都非常有帮助，涉及到了200个语言，并且配有相应的国际音标转写，以及田野调查的一些笔记。 Sanders, N. C., & Chin, S. B. (2009). Phonological Distance Measures. Journal of Quantitative Linguistics, 16(1), 96–114. Vakulenko, M. (2021). Calculation of Phonetic Distances between Speech Sounds. Journal of Quantitative Linguistics, 28(3), 223–236. Vakulenko, M. O. (2023). Unified Parametrization of Phonetic Features and Numerical Calculation of Phonetic Distances between Speech Sounds. Journal of Quantitative Linguistics, 30(1), 67–85. Bauer, L., Warren, P., Bardsley, D., Kennedy, M., & Major, G. (2007). New Zealand English. Journal of the International Phonetic Association, 37(01), 97. https://doi.org/10.1017/S0025100306002830 Chen, Y., & Gussenhoven, C. (2015). Shanghai Chinese. Journal of the International Phonetic Association, 45(3), 321–337. https://doi.org/10.1017/S0025100315000043 Cox, F., & Palethorpe, S. (2007). Australian English. Journal of the International Phonetic Association, 37(3), 341–350. https://doi.org/10.1017/S0025100307003192 Deterding, D. (1997). The Formants of Monophthong Vowels in Standard Southern British English Pronunciation. Journal of the International Phonetic Association, 27(1–2), 47–55. https://doi.org/10.1017/S0025100300005417 Deterding, D. (2006). The North Wind versus a Wolf: Short texts for the description and measurement of English pronunciation. Journal of the International Phonetic Association, 36(2), 187–196. https://doi.org/10.1017/S0025100306002544 Ferragne, E., & Pellegrino, F. (2010). Formant frequencies of vowels in 13 accents of the British Isles. Journal of the International Phonetic Association, 40(1), 1–34. https://doi.org/10.1017/S0025100309990247 Harrington, J., Palethorpe, S., & Watson, C. (2000). Monophthongal vowel changes in Received Pronunciation: An acoustic analysis of the Queen’s Christmas broadcasts. Journal of the International Phonetic Association, 30(1–2), 63–78. https://doi.org/10.1017/S0025100300006666 Hawkins, S., & Midgley, J. (2005). Formant frequencies of RP monophthongs in four age groups of speakers. Journal of the International Phonetic Association, 35(2), 183–199. https://doi.org/10.1017/S0025100305002124 Hillenbrand, J. M. (2003). American English: Southern Michigan. Journal of the International Phonetic Association, 33(1), 121–126. https://doi.org/10.1017/S0025100303001221 Lee, W.-S., & Zee, E. (2003). Standard Chinese (Beijing). Journal of the International Phonetic Association, 33(1), 109–112. https://doi.org/10.1017/S0025100303001208 Li, Q., Chen, Y., & Xiong, Z. (2019). Tianjin Mandarin. Journal of the International Phonetic Association, 49(1), 109–128. https://doi.org/10.1017/S0025100317000287 Roach, P. (2004). British English: Received Pronunciation. Journal of the International Phonetic Association, 34(2), 239–245. https://doi.org/10.1017/S0025100304001768 Shi, M., & Chen, Y. (undefined/ed). Lili Wu Chinese. Journal of the International Phonetic Association, 1–23. https://doi.org/10.1017/S0025100320000092 Villarreal, D., Clark, L., Hay, J., & Watson, K. (2020). From categories to gradience: Auto-coding sociophonetic variation with random forests. Laboratory Phonology: Journal of the Association for Laboratory Phonology, 11(1), Article 1. https://doi.org/10.5334/labphon.216 Watson, K. (2007). Liverpool English. Journal of the International Phonetic Association, 37(3), 351–360. https://doi.org/10.1017/S0025100307003180 Zee, E. (1991). Chinese (Hong Kong Cantonese). Journal of the International Phonetic Association, 21(1), 46–48.