Miyazawa’s Pukiwiki
R Note/音響解析データの統計解析
はすでに存在します。
開始行:
*音響解析（音声コーパス）データの統計解析 [#aa6f3f70]
#contents
**準備 [#x5be6c17]
***R と R Studio のインストール [#eeb2823d]
-&pgid(,R Note);
***サンプルデータ [#ce3453b5]
#ref(CSJ_MFCC26_short.zip);
→ &color(red){''これは話者数の少ないデータなので、下記の...
//-[[宮澤幸希，白勢彩子，馬塚れい子，菊池英明，”知識獲得...
-CSJの音声データファイルの冒頭120秒をMFCC解析した結果。各...
filename CSJのファイル名。12ファイルある。
register レジスタ名。"APS"は学会講演録音。"Re-Read"は学...
speakerID 話者名。APSとRe-Readについて同じ話者がいるため...
gender 話者の性別。
seg 音素名。SclSは<cl>、SsvSは<sv>。音素ラベルのないフ...
segID 音素ID。各ファイルについて音素ごとに数値をつけて...
frameID フレームID。各ファイルについて0.01秒ごとに数値...
beg_time フレームの開始時間。
MFCC1～MFCC12 MFCC12次元の解析結果。
C0 パワーの解析結果。
dMFCC1～dMFCC12 デルタMFCC12次元の解析結果。
dC0 デルタパワーの解析結果。
--ここでの「フレーム」は[[音声解析の分析窓:http://speechr...
--解析の詳細は &pgid(,MFCC解析のツール);/[[HTK(Hidden Mar...
***パッケージインストール [#k8de8a47]
-Rを起動して以下のコマンドを実行します。初回に一回だけ実...
>install.packages("plyr") #データテーブルの集計用パッ...
>install.packages("ggplot2") #視覚化のためのパッケージ
>install.packages("reshape") #データテーブルの整形用パッ...
>install.packages("lme4") #一般化線形モデルのためのパ...
>install.packages("lmerTest") #一般化線形モデルの検定のた...
>install.packages("pwr") #検定力解析（パワーアナリシ...
**データテーブルの読み込み [#a560c213]
-ダウンロードした &ref(CSJ_MFCC26_short.zip); を解凍し、...
-Rを起動して、以下を実行((R Studioであれば、右上の「Envir...
>data <- read.table("E:/Users/m-kouki/Desktop/CSJ_MFCC26....
--パスの区切り文字は「/」または「\\」
--今回のデータの1行目はヘッダ（各列のデータの詳細を説明し...
--参考：[[R-Source 40. ファイルからデータを読み込む:http:...
-読み込まれたデータは変数 data に格納されます。以下のよう...
>data #データテーブルの中身を見る
>head(data) #データテーブルの冒頭数行を見る
>names(data) #変数名の一覧を得る
-
[1] "filename" "register" "speakerID" "gender" "se...
[8] "beg_time" "MFCC1" "MFCC2" "MFCC3" "MF...
[15] "MFCC7" "MFCC8" "MFCC9" "MFCC10" "MF...
[22] "dMFCC1" "dMFCC2" "dMFCC3" "dMFCC4" "dM...
[29] "dMFCC8" "dMFCC9" "dMFCC10" "dMFCC11" "dM...
-データの一部分を指定することができます。
>data[ , 1:4] #1列目から4列目までの全ての行
>data[ 1:4 ,] #1行目から4行目までの全ての列
**データの前処理 [#l688e65c]
***要約統計量を得る [#u5c296b9]
-以下を実行します。
>summary(data)
--[[データの尺度:http://speechresearch.fiw-web.net/114.ht...
---'''連続型変数 Area, Slope, Soil.pH, Worm.density は6つ...
--カテゴリカル変数のつもりで作成したデータが、summaryを表...
speakerID
Min. :1.000
1st Qu.:2.000
Median :4.000
Mean :3.596
3rd Qu.:5.000
Max. :6.000
--実際には、speakerIDの数値はID番号なので、平均値や中央値...
>data$speakerID <- factor(data$speakerID)
>data$segID <- factor(data$segID)
>data$frameID <- factor(data$frameID)
--もう一度summaryを確認します。
>summary(data)
-データの総数は以下で求められます。
>nrow(data)
***特定のデータを抜き出す [#lacf4d9c]
-単に特定の列だけ抜き出したい場合、例えば以下のように書け...
>data[ c("speakerID","MFCC1") ]
-条件に該当する行だけ抜き出したい場合、例えば以下のような...
>data[ C0 > 20 & gender == "F" , 1:9] # C0の値が20以上の...
-subset関数を使うやり方が便利。[[Rデータフレームから一部...
>subset(data, (C0 > 20) & (gender == "F"), select = 1:9 )
--ヘッダ名を指定して取り出すこともできる。
>subset(data, (C0 > 20) & (gender == "F"), select = c(spe...
-特定のデータのみを抜き出して新しいデータテーブルを作るこ...
>dataF <- subset(data, (gender == "F")) # 女性のみを取り...
>summary(dataF)
-欠損値を除外したい場合は na.omit 関数を使う。
>data <- na.omit(data)
--どれか一つの列に欠損値を含んでいた場合、その行ごと削除...
-複数の条件で抜き出したデータテーブルを結合する場合は、rb...
>dataV <- rbind( subset(data, seg=="i"), subset(data, se...
>summary(dataV)
--特定のデータを抜き出した後のデータテーブルをsummaryする...
seg
a :3519
o :3451
u :3304
i :1719
e :1514
aH : 0
(Other): 0
--再度factor関数にかけてやると、カウント数0のデータは消え...
>dataV$seg <- factor(dataV$seg)
>summary(dataV)
***特定のデータの計算結果を追記する [#e3780abf]
-MFCC1～MFCC12の平均を求めて、データテーブルに"MFCC.Mean"...
>data$MFCC.Mean <- (data$MFCC1 + data$MFCC2 + data$MFCC3 ...
--apply関数を使うこともできます。
>data$MFCC.Mean <- apply(data[, 9:20], 1, mean)
---各行について集計したい場合は第二引数に1、各列について...
--値のどこかにNaNを含む可能性があって、その値を除外して平...
>data$MFCC.Mean <- apply(data[, 9:20], 1, mean, na.rm=TRUE)
-dMFCC1～dMFCC12の二乗和((この値は音響的特徴の変動量の大...
>data$dMFCC.SquareSum <- data$dMFCC1^2 + data$dMFCC2^2 + ...
>summary(data)
-対数変換したデータを追記する
>data$dMFCC.SquareSum.log = log(data$dMFCC.SquareSum)
-C0を平均値が0、分散が1になるように標準化したC0.scale列を...
>data$C0.scale = scale(data$C0)
>summary(data)
--列名は C0.scale.V1 になっている。→ &color(red){''理由は...
-【発展】条件を指定してデータを追記するには、transform 関...
>data = transform(data, seg2 = ifelse( ( (seg == "i") | (...
--さらに条件を複雑にして、「短母音5種類には"sV"、長母音5...
>data = transform(data, seg3 = ifelse( ( (seg == "i") | (...
-【発展】特定のデータを含む行に隣接する行のデータを得る
--例えば、「短母音の行の直前と直後の行を探して、該当する...
--データテーブルの基本的な考え方は、「''1行の中に解析に必...
---まずは、seg列のみを取り出して新しいデータテーブルdata....
>data.seg = data[ c("seg") ]
---空欄を意味する1行1列のデータテーブルを作ります。
>seg.nothing = data.frame(seg=c("Nothing"))
---data.segを前に1列ずらしたデータテーブル data.seg.prev ...
>data.seg.prev = rbind(seg.nothing, data.seg)
>data.seg.next = rbind(data.seg, seg.nothing)
---続いて、data.seg.prev の最初の列、 data.seg.next の最...
>data.seg.prev = data.seg.prev[1:(nrow(data.seg.prev)-1) ...
>data.seg.next = data.seg.next[2:nrow(data.seg.next) , ]
---元のデータテーブル data に data.seg.prev と data.seg.n...
> data$seg.prev = data.seg.prev
> data$seg.next = data.seg.next
> head(data)
---これで、必要な情報は全て1行の中に集約されました。
--transform関数を使って該当条件のラベルを付加します。
>data = transform(data, seg.neighbor = ifelse( ( (seg.pre...
--このサンプルデータは各行が「フレーム（分析窓）」単位で...
***特定のデータの集計値を得る [#p0a9dcd9]
-ddply関数が大変便利です((浅井君、教えてくれてありがとう)...
-[[plyrパッケージのインストール:http://speechresearch.fiw...
>library(plyr)
-様々な実例は [[「plyrパッケージで君も前処理スタ☆」改め「...
-各話者・各レジスタごとに、C0の平均値と標準偏差を求めて、...
>dataSumm <- ddply(data, c('speakerID', 'register'), summ...
>dataSumm
--指定した条件ごとの集計値が出力されます。
speakerID register C0mean C0sd
1 1 APS 15.69648 2.949087
2 1 Re-Read 13.24777 2.565976
3 2 APS 15.75885 2.861640
4 2 Re-Read 13.95841 2.957040
5 3 APS 14.53344 2.584729
6 3 Re-Read 14.81330 2.511646
7 4 APS 13.91550 2.801826
8 4 Re-Read 13.18621 3.143895
9 5 APS 15.02620 2.385534
10 5 Re-Read 14.91980 3.005093
11 6 APS 15.10699 2.306960
12 6 Re-Read 14.23184 2.976182
-各話者・各レジスタごとに、C0の平均値と標準偏差を求めて、...
>data <- ddply(data, c('speakerID', 'register'), transfor...
>data
--元のデータに列が追加されています。
-各話者について、C0の値の平均値が0、標準偏差が1になるよう...
>data <- ddply(data, c('speakerID'), transform, C0.z = (C...
>head(data)
-【発展】特定条件の集計値を全行に追加する
--例えば「レジスタAPSに限定して、話者ごとにC0の平均値を求...
>#レジスタAPS以外のC0値をNAにした新しい列「C0_APS」を作り...
>data = transform(data, C0_APS = ifelse( register == "AP...
>#各話者のAPSのC0の平均値を求めて、すべての行に追記する
>data = ddply(data, c('speakerID'), transform, C0_APSmean...
-【発展】フレーム単位のデータテーブルから音素単位のデータ...
--サンプルデータ data の各行は、窓長0.01秒の時間フレーム...
--segIDについて集計値を求めます。segIDは異なる話者で重複...
> data2 <- ddply(data, c('speakerID', 'segID'), summarize...
> nrow(data2)
--
[1] 7153
---他のfactor変数を条件に加えても構いません。データの総行...
> data2 <- ddply(data, c('filename', 'register', 'speaker...
> nrow(data2)
---
[1] 7153
--よく行われる処理は、各音素について時間フレームの中心（...
> data2 <- ddply(data, c('speakerID', 'segID'), summarize...
---length(MFCC1) は各条件で抽出されたMFCC1のベクトルサイ...
---round(length(MFCC1)/2) は各音素のフレーム長の半分を四...
---MFCC1[round(length(MFCC1)/2)] は各音素の時間的中心点の...
***外れ値を除外する [#ma6808aa]
-ddply関数を使って、「話者＆レジスタごとにC0の平均値と標...
-まず、平均値＋2.0SD を上回る値は1、そうでない値は0をつけ...
>data <- ddply(data, c('speakerID', 'register'), transfor...
-続いて、平均値－2.0SD を下回る値は1、そうでない値は0をつ...
>data <- ddply(data, c('speakerID', 'register'), transfor...
-C0out1またはC0out2が1になっている行ははずれ値として削除...
>data2 = subset(data, C0out1=="0")
>data2 = subset(data2, C0out2=="0")
>nrow(data)
>nrow(data2)
***ランダムサンプリング [#x9fff2b3]
-ddply関数を使えば、指定した各条件について、該当するデー...
-話者＆レジスタ＆各音素カテゴリーごとに、C0とdC0をそれぞ...
>sampling = 10
>dataRand = ddply(data, c('speakerID', 'register', 'seg')...
>dataRand
--サンプリングしたC0とdC0の行番号は一致していない（それぞ...
--"replace = T"は、データが重複してサンプリングされること...
**データの視覚化 [#qcf6ad0a]
-データの特性を理解するために、視覚化はとても重要です。
***準備 [#add37ba8]
-ggplot2関数が大変多機能で、便利です((浅井君、教えてくれ...
-[[ggplot2パッケージのインストール:http://speechresearch....
>library(ggplot2)
>library(plyr)
>library(reshape)
-サンプルデータ &ref(CSJ_MFCC26_short.zip); をダウンロー...
-以後の説明では、短母音のみのデータテーブルを使います。
>data <- read.table("E:/Users/m-kouki/Desktop/CSJ_MFCC26....
>data$speakerID = factor(data$speakerID)
>data$segID = factor(data$segID)
>data$frameID = factor(data$frameID)
>dataV = subset(data, (seg=="i")|(seg=="e")|(seg=="o")|(s...
>summary(dataV)
***散布図 [#madd8730]
-単にplot関数に１つのデータ列を与えると、散布図になります。
>plot(dataV$C0)
#ref(Rplot_1.png);
--縦軸がC0の値、横軸は左からデータの行数です。
--ここで、分析に使う各パラメータの散布図を確認しておくと...
-plot関数に２つのデータ列を与えてみます。
>plot(dataV$MFCC2, dataV$MFCC4)
#ref(Rplot_2.png);
--ここでは、縦軸がMFCC4の値、横軸がMFCC2の値です。
--女性話者のデータに限定してみます。
>plot(subset(dataV, gender=="F")$MFCC2, subset(dataV, gen...
-ggplotを使うと、条件ごとの色分けも簡単にできます。以下は...
>ggplot(dataV) +
> geom_point(aes(x=MFCC2, y=MFCC4, color=seg)) +
> labs(title="All")
#ref(Rplot_3.png);
--女性話者のデータに限定してみます。後でグラフに機能を追...
>p = ggplot(subset(dataV, gender=="F")) +
> geom_point(aes(x=MFCC2, y=MFCC4, color=seg)) +
> labs(title="Female")
>p
--グラフの軸のサイズを指定することもできます。上記の変数 ...
>p = p + xlim(-30, 30) + ylim(-30, 30)
>p
#ref(Rplot_4.png);
-ddplyと組み合わせると色々と便利です。
--母音ごと・レジスタごとに平均値を求めて、レジスタで色分...
>dataSumm = ddply(dataV, c('seg', 'register'), summarize,...
>ggplot(dataSumm) +
> geom_text(aes(x=MFCC2mean, y=MFCC4mean, label=seg, col...
---散布図の各点を文字にしたい場合は、geom_textを使います...
#ref(Rplot_5.png);
--話者ごと・母音ごと・レジスタごとにMFCC2, MFCC4, MFCC6, ...
>dataSumm = ddply(dataV, c('speakerID', 'seg', 'register'...
>pairs(~dataSumm$MFCC2mean + dataSumm$MFCC4mean + dataSum...
#ref(Rplot_6.png);
---これはggplotを使っていません。((ggplotで同様のことはpl...
***ヒストグラム [#qcff188b]
-分析の最初のチェックとして、ヒストグラムの形状は非常に重...
--分析に用いる各要因の各水準ごとにヒストグラムをプロット...
--単峰性でも、ヒストグラムの形状がゆがんでいる場合、注意...
-全てのデータのMFCC1のヒストグラムを、通常の方法でプロット
>hist(dataV$MFCC1)
-ggplotを使用した方法
>ggplot( dataV, aes(x = MFCC1) ) + geom_histogram()
#ref(Rplot_7.png);
--全データのMFCC1は２峰性の分布になっています。各要因・各...
-ggplotとsubplotを組み合わせると、各要因・各水準のヒスト...
--話者別にヒストグラムをプロット
>ggplot(dataV) +
> geom_histogram( aes(x = MFCC1) ) +
> facet_wrap(~speakerID)
--レジスタ・話者・母音別にヒストグラムをプロット
>ggplot(dataV) +
> geom_histogram( aes(x = MFCC1) ) +
> facet_wrap(register~speakerID~seg)
--レジスタで色分けして、話者を縦グリッド・母音を横グリッ...
>ggplot(dataV) +
> geom_histogram( aes( x = MFCC1, fill=factor(register) ...
> facet_grid(speakerID~seg)
#ref(Rplot_8.png);
---/u/では全ての話者で二峰性の分布になっていることが分か...
-多次元ヒストグラム
--[[R generate 2D histogram from raw data:http://stackove...
--[[二変量データのヒストグラムと密度推定:http://www.okada...
***箱ヒゲ図（ボックスプロット）(([[統計学：Rを用いた入門...
-t検定や分散分析の前に、データの分布や有意傾向を視覚的に...
-全てのデータのMFCC1の箱ヒゲ図を、通常の方法でプロット
>boxplot(subset(dataV, select = c(MFCC1)))
--複数の列の値を同時にプロットすることもできます。
>boxplot(subset(dataV, select = c(MFCC1, MFCC2, MFCC3) ) )
#ref(Rplot_9.png);
--箱ヒゲ図の見方は以下の通りです。
(([[R+ 箱ひげ図の読み方:http://wb-nahce.info/rplus/statse...
---中央の横線：各データの全体の中央値（第2クォンタイル((...
---箱の上端の線：全体の中央値より大きいデータの中央値（第...
---箱の下端の線：全体の中央値より小さいデータの中央値（第...
---箱の中に、データ全体の約50%が含まれます。上下のヒゲは...
---ヒゲの外側にある値は、はずれ値として扱われます。
-全てのデータのMFCC1の箱ヒゲ図を、ggplotでは以下のように...
>ggplot(dataV) +
> geom_boxplot(aes(y = MFCC1, x = "MFCC1"))
-ggplotは、boxplotと特性が異なります。boxplotは異なる列の...
--レジスタで色分けして、母音ごとにMFCC1の値をプロット
>ggplot(dataV) +
> geom_boxplot(aes(y = MFCC1, x = seg, fill = register))
---ノッチを加えてプロット(('''この図にてへこんだ部分がノ...
>ggplot(dataV) +
> geom_boxplot(aes(y = MFCC1, x = seg, fill = register),...
#ref(Rplot_10.png);
--箱ヒゲ図は分布のばらつきを見るのに適していますが、''分...
>ggplot(dataV) +
> geom_violin(aes(y = MFCC1, x = seg, fill = register))
#ref(Rplot_15.png);
--レジスタで色分けして、話者ごとにMFCC1の値をプロットし、...
>ggplot(dataV) +
> geom_boxplot(aes(y = MFCC1, x = speakerID, fill = regi...
> facet_wrap( ~ seg)
#ref(Rplot_11.png);
-ggplotで、異なる列のデータを並べてプロットしたい場合は、...
>dataV.melt <- melt(subset(dataV, select = c(register, sp...
>head(dataV.melt)
--赤字のメッセージが出ますが問題ありません。dataV.meltの...
--レジスタで色分けして、MFCC1,MFCC2,MFCC3の値をプロットし...
>ggplot(dataV.melt) +
> geom_boxplot(aes(y = value, x = variable, fill = regis...
> facet_grid(speakerID ~ seg)
#ref(Rplot_12.png);
---この例では細かすぎて見づらくなっていますが、違いが出る...
***棒グラフ [#i5d11244]
-barplot関数を使った方法：[[バイオスタティスティクス 棒グ...
-ggplotを使う場合、まずはddplyを使ってデータを集計し、平...
--今回は、「各話者・各母音・各レジスター」別にMFCC2とMFCC...
--「各話者・各母音・各レジスター」別にMFCC2とMFCC4の平均...
>dataVSumm1 = ddply(dataV, c('speakerID', 'seg', 'registe...
>dataVSumm1
--reshapeパッケージのmelt関数を使って、データを並び替えま...
>dataVSumm1M = melt(dataVSumm1)
>dataVSumm1M
---%%エラーバーを表示しないなら、この時点でプロットできま...
>qplot(register, weight=value, geom="bar", data=dataVSumm...
---エラーバーを表示しないなら、その値の「各母音・各レジス...
>dataVSumm2 = ddply(dataVSumm1, c('seg', 'register'), sum...
>dataVSumm2
>dataVSumm2M = melt(dataVSumm2)
>dataVSumm2M
>
>qplot(register, weight=value, geom="bar", data=dataVSumm...
#ref(Rplot_13.png);
--エラーバーを表示する場合、その値の「各母音・各レジスタ...
>dataVSumm3 = ddply(dataVSumm1M, c('seg', 'register', 'va...
>dataVSumm3
---length(value)は、集計するデータサイズ。ここでは話者数...
--ggplot関数を使って棒グラフとエラーバーを表示します（参...
>ggplot(dataVSumm3, aes(x=variable, y=mean, fill=register...
> geom_bar(stat="identity", position="dodge") +
> geom_errorbar(aes(ymax=mean+se, ymin=mean-se), positio...
> facet_grid( ~ seg)
---"aes"の設定はgeom_barでもgeom_errorbarでも共通して使う...
#ref(Rplot_14.png);
---fillとfacetの設定を変えてプロットの形式を変更します。
--
>ggplot(dataVSumm3, aes(x=variable, y=mean, fill=seg)) +
> geom_bar(stat="identity", position="dodge") +
> geom_errorbar(aes(ymax=mean+se, ymin=mean-se), positio...
> facet_grid( ~ register)
#ref(Rplot_14a.png);
**統計解析 [#n2090f8b]
***準備 [#n9f75954]
-[[lme4パッケージとlmerTestパッケージのインストール:http:...
>library(lme4)
>library(lmerTest)
>library(plyr)
>library(reshape)
>library(pwr)
-サンプルデータ &ref(CSJ_MFCC26_short.zip); をダウンロー...
>data <- read.table("E:/Users/m-kouki/Desktop/CSJ_MFCC26....
>data$speakerID = factor(data$speakerID)
>data$segID = factor(data$segID)
>data$frameID = factor(data$frameID)
***t検定とパワーアナリシス [#l8aa50bb]
-対応のないt検定
--全てのデータのMFCC5の平均値とMFCC6の平均値の間に、有意...
--まずは箱ヒゲ図を書いてみます。
>boxplot(subset(data, select = c(MFCC5, MFCC6) ) )
#ref(Rplot_16.png);
---ほぼボックスが重複しているように見えますが、ノッチが極...
--F検定で等分散性を見てみます（Rではt.test関数で自動的に...
>var.test(data$MFCC5, data$MFCC6)
--
F test to compare two variances
data: data$MFCC5 and data$MFCC6
F = 1.017, num df = 40137, denom df = 40137, p-value = 0...
(略)
---p値は0.09なので、ぎりぎりですがMFCC5とMFCC6は「等分散...
--t検定を実施します。ここでは等分散であることを明記します。
>t.test(data$MFCC5, data$MFCC6, var.equal=TRUE)
--
Two Sample t-test
data: data$MFCC5 and data$MFCC6
t = 5.3208, df = 80274, p-value = 1.036e-07
(略)
---p値は0.001より小さいので、MFCC5とMFCC6の平均値は有意に...
--var.equalオプションなしで実行した場合、等分散であればt...
>t.test(data$MFCC4, data$MFCC9)
--
Welch Two Sample t-test
data: data$MFCC4 and data$MFCC9
t = 1.3549, df = 77047.3, p-value = 0.1754
(略)
---p値が0.05より大きいので、MFCC4とMFCC9の平均値は有意に...
-【注】t検定で比較する2群の標本数が大きく違う場合は、ウェ...
--上のデータの検定力を求めてみます。参考：[[効果量・検定...
--今回のt検定で使用した２つのデータサイズ確認
>length(data$MFCC4)
--
[1] 40138
--
>length(data$MFCC9)
--
[1] 40138
--検定力を求める
>pwr.t2n.test(n1=40138, n2=40138, sig.level=0.05, d=0.8)
---n1=40138 : データ群1のデータサイズ（サンプル数）です。
---n2=40138 : データ群2のデータサイズ（サンプル数）です。...
---d=0.8 : Cohenの効果量d には 0.8 を指定しています。これ...
---sig.level=0.05 : 有意水準0.05で統計を行う場合を想定し...
--結果は以下の通りです。
t test power calculation
n1 = 40138
n2 = 40138
d = 0.8
sig.level = 0.05
power = 1
alternative = two.sided
---power は 1 になりました。すなわち100回中統計を実施した...
---power の値は 0.8 以上であれば（100回中80回同じ結果にな...
--有意水準と検定力の関係について：[[例数設計の基礎 第８回...
-対応のあるt検定
--上でMFCC4とMFCC9の平均値に有意差がありませんでしたが、...
--そこで、話者ごとにMFCC4とMFCC9の平均値を求めて、その値...
> dataSumm1 <- ddply(data, c('speakerID'), summarize, MFC...
> dataSumm1
--
speakerID MFCC4 MFCC9
1 1 -0.36694394 -0.67979122
2 2 -0.02150908 0.28991697
3 3 0.23361956 -1.49124097
4 4 0.12591752 0.05997551
5 5 -0.97832807 0.45119489
6 6 -0.09088042 -0.02595399
--このデータは、各行が同じ話者から得られたデータなので、...
> t.test(dataSumm1$MFCC4, dataSumm1$MFCC9, paired = TRUE)
--
Paired t-test
data: dataSumm1$MFCC4 and dataSumm1$MFCC9
t = 0.1191, df = 5, p-value = 0.9098
---p値が0.05より大きいので、MFCC4とMFCC9の平均値は有意に...
--ただし、今回のデータは、''集計前よりもあまりにも数が少...
---ある数のデータセットから得られた統計の結果が信頼のでき...
--もう少し対応のあるデータを増やします。レジスタ・音素・...
> dataSumm2 <- ddply(data, c('speakerID', 'register', 'se...
> dataSumm2
--これでデータ数は454個になりました。対応のあるt検定を実...
> t.test(dataSumm2$MFCC4, dataSumm2$MFCC9, paired = TRUE)
--
Paired t-test
data: dataSumm2$MFCC4 and dataSumm2$MFCC9
t = 2.1648, df = 453, p-value = 0.03093
---p値は0.001より小さくなりました。したがって、MFCC4とMFC...
--このように、特に自然発話データの場合、どのような集計を...
***単回帰分析 [#tc65193f]
-お題：''各音素の持続時間（フレーム数）''と、''各音素のdM...
--音素が長いほど、発音に使える時間に余裕があるので、各音...
--分析に使用するデータテーブル dataSumm.SS を作ります。各...
>data$dMFCC.SquareSum <- data$dMFCC1^2 + data$dMFCC2^2 + ...
>dataSumm.SS <- ddply(data, c('speakerID', 'segID'), summ...
>head(dataSumm.SS)
--
speakerID segID seg.length dMFCC.SquareSum
1 1 4627 4 54.75966
2 1 4628 2 55.05001
3 1 4629 3 45.30628
4 1 4630 6 46.35851
5 1 4631 6 41.52738
6 1 4632 13 39.62246
-ここで、''目的変数（従属変数）''は dMFCC.SquareSum で、'...
--単回帰分析では、目的・説明の両変数は正規分布する連続変...
>par(mfrow=c(1,2))
>hist(dataSumm.SS$seg.length)
>hist(dataSumm.SS$dMFCC.SquareSum)
#ref(Rplot_17.png);
--どちらも、L字型の分布になっています。持続時間や二乗和の...
>dataSumm.SS$seg.length.log = log(dataSumm.SS$seg.length)
>dataSumm.SS$dMFCC.SquareSum.log = log(dataSumm.SS$dMFCC....
>par(mfrow=c(1,2))
>hist(dataSumm.SS$seg.length.log)
>hist(dataSumm.SS$dMFCC.SquareSum.log)
#ref(Rplot_18.png);
---これで正規分布に近づきました。
--続いて、散布図を描いてみます。
>par(mfrow=c(1,1))
>plot( dataSumm.SS$seg.length.log, dataSumm.SS$dMFCC.Squa...
#ref(Rplot_19.png);
--分布が縞々になってしまいました。これは、seg.length が整...
--相関分析にかけたところ、有意なある程度の負の相関（p-val...
>cor.test( dataSumm.SS$seg.length.log, dataSumm.SS$dMFCC....
--lm関数を使って回帰分析にかけます。
>ss.lm <- lm(dMFCC.SquareSum.log ~ seg.length.log, data =...
>summary(ss.lm)
---lm関数の引数で回帰式を指定します。
y(目的変数) = a(切片) + b(傾き)x(説明変数) → y~x
y = bx → y~-1+x
---結果は以下の通りです。
--
Call:
lm(formula = dMFCC.SquareSum.log ~ seg.length.log, data ...
Residuals:
Min 1Q Median 3Q Max
-3.00788 -0.29871 0.07117 0.36034 1.36548
---Residuals（残差）は、データと回帰モデルが予測したの値...
--
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.882726 0.014611 265.74 <2e-16 ***
seg.length.log -0.251529 0.008559 -29.39 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1...
Residual standard error: 0.512 on 6693 degrees of freedom
---(Intercept) は切片a、seg.length.log は説明変数の傾き（...
---目的変数に対して、説明変数 seg.length.log は有意な負の...
--
Multiple R-squared: 0.1143, Adjusted R-squared: 0.1142
F-statistic: 863.7 on 1 and 6693 DF, p-value: < 2.2e-16
--さきほどの図に回帰直線を追記します。
>abline(ss.lm)
--回帰診断図を描きます。
>par(mfrow=c(2,2))
>plot(ss.lm)
***一元配置分散分析＋多重比較法 [#i12035e7]
&color(red){''この節の手順は誤りがある可能性があります。...
-お題：''dMFCC1～dMFCC12の二乗和は、母音・ソノラント子音...
-３群以上の間の有意差の比較は、[[t検定:http://speechresea...
--dMFCC1～dMFCC12の二乗和を求め、[[ここ:http://speechrese...
>data$dMFCC.SquareSum <- data$dMFCC1^2 + data$dMFCC2^2 + ...
>data$dMFCC.SquareSum.log <- log(data$dMFCC.SquareSum)
>data = transform(data, seg2 = ifelse( ( (seg == "i") | (...
>data$seg2 <- factor(data$seg2)
>summary(data)
--箱ヒゲ図を表示してみます。
>data.melt.aov1 <- melt(subset(data, select = c(seg2, dMF...
>library(ggplot2)
>ggplot(data.melt.aov1) + geom_boxplot(aes(y = value, x =...
#ref(Rplot_20.png);
-対応のない一元配置分散分析
--関数aovを使って分散分析を実行します。
>ss.aov1 <- aov(dMFCC.SquareSum.log ~ seg2, data = data)
>ss.aov1
--
Call:
aov(formula = dMFCC.SquareSum.log ~ seg2, data = data)
Terms:
seg2 Residuals
Sum of Squares 1032.00 33207.52
Deg. of Freedom 2 40135
Residual standard error: 0.9096128
Estimated effects may be unbalanced
--分散分析表を見るには、summaryもしくはanova関数を使いま...
>summary(ss.aov1)
--
Df Sum Sq Mean Sq F value Pr(>F)
seg2 2 1032 516.0 623.6 <2e-16 ***
Residuals 40135 33208 0.8
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1...
---意味は以下の通りです。((分散分析表の枠は[[一元配置の分...
-------------------------------------------------------
要因 自由度 平方和 平均平方和 F-値 P-値
-------------------------------------------------------
F 2 1032 516.0 623.6 <2e-16 ***
誤差 40135 33208 0.8
-------------------------------------------------------
---P値が0.05より小さいため、説明変数（seg2）の各水準（vow...
--この結果だけでは、どの水準とどの水準に差があったのかは...
>TukeyHSD(ss.aov1)
--
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = dMFCC.SquareSum.log ~ seg2, data = da...
$seg2
diff lwr upr p...
sonorant-non-sonorant -0.1880719 -0.2295339 -0.1466100 ...
vowel-non-sonorant -0.3335195 -0.3556718 -0.3113673 ...
vowel-sonorant -0.1454476 -0.1870469 -0.1038483 ...
---diff : 平均値の差、lwr : 信頼区間の下限、upr : 信頼区...
---p値がゼロになってしまった。→ &color(red){''理由は？''};
--[[R言語で統計解析入門： くり返しのある一元配置分散分析...
>pairwise.t.test( data$dMFCC.SquareSum.log, data$seg2, p....
--
non-sonorant sonorant
sonorant <2e-16 -
vowel <2e-16 <2e-16
---各組み合わせについて有意差がある。
-対応のある一元配置分散分析
--レジスタ・話者別、また音素の種類ラベル（説明変数）別に...
> library(plyr)
> dataSumm.aov1p <- ddply(data, c('speakerID', 'register'...
> dataSumm.aov1p
--話者について対応をとり分散分析にかける場合、以下のよう...
> ss.aov1p.speaker <- aov(dMFCC.SquareSum.log ~ seg2 + sp...
> summary(ss.aov1p.speaker)
--
Df Sum Sq Mean Sq F value Pr(>F)
seg2 2 0.6355 0.3177 24.398 7.34e-07 ***
speakerID 5 0.1842 0.0368 2.829 0.0344 *
Residuals 28 0.3646 0.0130
(略)
---speakerIDの有意な影響が認められます。
---多重比較を行います。
--
> TukeyHSD(ss.aov1p.speaker)
--
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = dMFCC.SquareSum.log ~ seg2 + speakerI...
$seg2
diff lwr upr ...
sonorant-non-sonorant -0.1950974 -0.3103730 -0.07982184 ...
vowel-non-sonorant -0.3231252 -0.4384008 -0.20784963 ...
vowel-sonorant -0.1280278 -0.2433034 -0.01275223 ...
(略)
---各組み合わせについて有意差がある。
--レジスタについて対応をとり分散分析にかける場合、以下の...
> ss.aov1p.register <- aov(dMFCC.SquareSum.log ~ seg2 + s...
> summary(ss.aov1p.register)
--
Df Sum Sq Mean Sq F value Pr(>F)
seg2 2 0.6355 0.3177 19.276 3.21e-06 ***
register 1 0.0214 0.0214 1.299 0.263
Residuals 32 0.5275 0.0165
---registerの有意な影響は認められません。
---多重比較を行います。
--
>TukeyHSD(ss.aov1p.register)
--
(略)
$seg2
diff lwr up...
sonorant-non-sonorant -0.1950974 -0.3238975 -0.066297298...
vowel-non-sonorant -0.3231252 -0.4519253 -0.194325088...
vowel-sonorant -0.1280278 -0.2568279 0.000772313...
(略)
---registerの有意な影響がないにもかかわらず、レジスタにつ...
--speakerIDとregisterを結合して一つのラベル（speakerIDreg...
>dataSumm.aov1p$speakerIDregister = paste(dataSumm.aov1p$...
>dataSumm.aov1p
>ss.aov1p.speakerIDregister <- aov(dMFCC.SquareSum.log ~ ...
>summary(ss.aov1p.speakerIDregister)
>TukeyHSD(ss.aov1p.speakerIDregister)
---結果は省略します。
***二元配置分散分析＋交互作用効果の扱い [#e8bbb4b5]
&color(red){''この節の手順は誤りがある可能性があります。...
-お題：''dMFCC1～dMFCC12の二乗和は、母音・ソノラント子音...
--説明変数（独立変数）が２種類になりました。第１要因は「...
--データは一元配置分散分析と同じものを使います。
>data$dMFCC.SquareSum <- data$dMFCC1^2 + data$dMFCC2^2 + ...
>data$dMFCC.SquareSum.log <- log(data$dMFCC.SquareSum)
>data = transform(data, seg2 = ifelse( ( (seg == "i") | (...
>data$seg2 <- factor(data$seg2)
--交互作用効果をプロットしてみます。
>interaction.plot(data$register, data$seg2, data$dMFCC.Sq...
#ref(Rplot_21.png);
---線が平行でないので、２つの説明変数の間に交互作用があり...
-対応のない二元配置分散分析
--関数aovを使って分散分析を実行します。
> ss.aov2 <- aov(dMFCC.SquareSum.log ~ seg2 * register, d...
> summary(ss.aov2)
---「seg2 * register」は、「音素の種類、レジスタ、音素の...
--
Df Sum Sq Mean Sq F value Pr(>F)
seg2 2 1034 517.2 627.66 <2e-16 ***
register 1 61 60.7 73.68 <2e-16 ***
seg2:register 2 74 37.0 44.96 <2e-16 ***
Residuals 40132 33070 0.8
---交互作用効果は「seg2:register」と表現されます。Pr <2e-...
--交互作用効果が確認されたときは、''単純（主）効果の検定'...
---交互作用効果が有意であった場合、主効果について検討する...
---交互作用効果が有意でなければ、Tukeyの多重比較法だけでO...
---参考サイト：[[単純主効果の検定と多重比較:http://igproj...
---参考情報：'''[[R Note/統計/分散分析:http://speechresea...
--以下、作成中
***混合効果モデル：独立変数が一つかつ連続変数 [#j49fc621]
-混合効果（単）回帰分析
-参考：[[Mixed Effects Models Blog:http://mixedmodeljp.bl...
-お題：''話者の個人差''及び''音素の種類''の影響を除外して...
--[[単回帰分析:http://speechresearch.fiw-web.net/121.html...
--分析に使用するデータテーブル dataSumm.SS を作ります。各...
>data$dMFCC.SquareSum <- data$dMFCC1^2 + data$dMFCC2^2 + ...
>dataSumm.SS <- ddply(data, c('speakerID', 'seg', 'segID'...
>head(dataSumm.SS)
--
speakerID seg segID seg.length dMFCC.SquareSum
1 1 a 5169 9 59.40295
2 1 a 5213 7 24.76156
3 1 a 5311 10 27.79712
4 1 a 5325 6 22.73648
5 1 a 5331 8 19.03984
6 1 a 5458 10 46.80528
---ここで、''目的変数（従属変数）''は dMFCC.SquareSum で...
---&color(red){''重要''：ランダム要因はランダムサンプリン...
--目的・説明の両変数は対数変換します。
>dataSumm.SS$seg.length.log = log(dataSumm.SS$seg.length)
>dataSumm.SS$dMFCC.SquareSum.log = log(dataSumm.SS$dMFCC....
--lmer関数を使ってモデルを構築します。
>ss.lmer <- lmer(dMFCC.SquareSum.log ~ seg.length.log + (...
---従属変数 ~ 独立変数 + (ランダム変数) のように書く。（...
---（）の中の｜の左側の数字１は切片を示す。「(1 + seg.len...
--結果は以下の通りです。&color(red){lmer関数の実行結果は...
>summary(ss.lmer)
--
Linear mixed model fit by REML
t-tests use Satterthwaite approximations to degrees of ...
Formula: dMFCC.SquareSum.log ~ seg.length.log + (1 + seg...
Data: dataSumm.SS
REML criterion at convergence: 9438.4
Scaled residuals:
Min 1Q Median 3Q Max
-5.7769 -0.5883 0.0893 0.6904 3.0492
Random effects:
Groups Name Variance Std.Dev. Corr
seg (Intercept) 0.224867 0.47420
seg.length.log 0.038383 0.19592 -0.76
speakerID (Intercept) 0.007856 0.08864
seg.length.log 0.001212 0.03481 -0.58
Residual 0.213650 0.46222
Number of obs: 7101, groups: seg, 52; speakerID, 6
Fixed effects:
Estimate Std. Error df t value Pr(>...
(Intercept) 3.76566 0.08751 35.49000 43.033 < 2...
seg.length.log -0.15318 0.03945 37.86000 -3.883 0.00...
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1...
Correlation of Fixed Effects:
(Intr)
sg.lngth.lg -0.770
---各表示の詳細は[[Mixed Effects Models Blog - R言語を用...
---「Fixed effects:」以下が固定因子の効果量。(Intercept) ...
-文献によっては、「ランダム効果は切片を変数にする場合が多...
--このようなモデルは以下のようになります。
>ss2.lmer <- lmer(dMFCC.SquareSum.log ~ seg.length.log + ...
>summary(ss2.lmer)
--
(略)
REML criterion at convergence: 9542.9
(略)
Fixed effects:
Estimate Std. Error df t value ...
(Intercept) 3.74862 0.05653 41.00000 66.31 ...
seg.length.log -0.17120 0.01106 6872.00000 -15.48 ...
(略)
--ランダム変数が切片と傾きに与える影響を考慮したモデル（s...
>anova(ss.lmer, ss2.lmer)
--
Models:
..1: dMFCC.SquareSum.log ~ seg.length.log + (1 | speaker...
..1: seg)
object: dMFCC.SquareSum.log ~ seg.length.log + (1 + seg....
object: speakerID) + (1 + seg.length.log | seg)
Df AIC BIC logLik deviance Chisq Chi Df P...
..1 5 9541.8 9576.1 -4765.9 9531.8 ...
object 9 9447.8 9509.6 -4714.9 9429.8 101.95 4 ...
---ここで「..1」は ss2.lmer、「object」は ss.lmer です。
---まず、P値（Pr）に着目します。P値が0.05より小さいので、...
---P値が0.05より大きい（両モデルに有意な差がない）場合、A...
---対数尤度を使ったより厳密なモデル評価の方法は [[Mixed E...
-''まずは、できるだけ情報量の多いモデルを採用するべきです...
--Model failed to converge: degenerate Hessian with ○ ne...
--Model is nearly unidentifiable: very large eigenvalue
--the random-effects parameters and the residual variance...
--fixed-effect model matrix is rank deficient so dropping...
***混合効果モデル：独立変数が一つかつ離散変数 [#o1ad649c]
-混合効果（一元配置）分散分析
-参考：[[混合モデルを使って反復測定分散分析をする:http://...
--[[混合効果単回帰モデル:http://speechresearch.fiw-web.ne...
***混合効果モデル：独立変数が複数 [#f7c9443e]
-独立変数の性質によって呼び名はいろいろありますが、まとめ...
--全ての独立変数が連続：混合効果（重）回帰分析
--全ての独立変数が離散：混合効果（多元配置）分散分析
--連続・離散の独立変数が混在：混合効果共分散分析
-参考：モデルの考え方は[[混合効果モデル（変量効果モデル、...
-分析に使用するデータテーブル dataSumm.SS を作ります。
>data$dMFCC.SquareSum <- data$dMFCC1^2 + data$dMFCC2^2 + ...
>data = transform(data, seg2 = ifelse( ( (seg == "i") | (...
>data$seg2 <- factor(data$seg2)
>dataSumm.SS <- ddply(data, c('speakerID', 'seg', 'seg2',...
>dataSumm.SS$seg.length.log = log(dataSumm.SS$seg.length)
>dataSumm.SS$dMFCC.SquareSum.log = log(dataSumm.SS$dMFCC....
-お題１：''話者の個人差''の影響を除外した上で、''各音素の...
--ここで、''目的変数（従属変数）''は dMFCC.SquareSum(.log...
---seg と register は離散変数なので、混合効果（二元配置）...
--lmer関数を使ってモデルを構築します。
>ssANOVA1.lmer <- lmer(dMFCC.SquareSum.log ~ seg * regist...
---''「*」で繋げられた独立変数は交互作用を考慮''します。...
---ランダム変数は、切片・各独立変数の傾き・各独立変数の交...
---上記のモデルを実行すると、fixed-effect model matrix is...
>ssANOVA1_1.lmer <- lmer(dMFCC.SquareSum.log ~ seg * regi...
---それでも警告は消えません。続いてランダム変数の傾きを全...
>ssANOVA1_2.lmer <- lmer(dMFCC.SquareSum.log ~ seg * regi...
---それでも警告は消えません。仕方がないので、独立変数同士...
>ssANOVA1_3.lmer <- lmer(dMFCC.SquareSum.log ~ seg + regi...
--結果は以下で確認できます。詳細は省略します。
>anova(ssANOVA1_3.lmer)
>summary(ssANOVA1_3.lmer)
>step(ssANOVA1_3.lmer)
-お題２：''話者の個人差''の影響を除外した上で、''各音素の...
--''目的変数（従属変数）''は dMFCC.SquareSum(.log) で、''...
--lmer関数を使ってモデルを構築します。
>ssANOVA2.lmer <- lmer(dMFCC.SquareSum.log ~ seg2 * regis...
---今度は警告は出ませんでした。
--結果は以下で確認できます。
>anova(ssANOVA2.lmer)
--
Analysis of Variance Table of type 3 with Satterthwaite
approximation for degrees of freedom
Sum Sq Mean Sq NumDF DenDF F.value Pr...
seg2 110.023 55.011 2 5.8759 199.646 3.966...
register 3.356 3.356 1 8.1706 12.178 0.00...
seg2:register 4.192 2.096 2 7.5733 7.608 0.01...
---独立変数間の有意な交互作用と、各独立変数の有意な影響が...
--単純主効果の確認（≒多重比較）は下記のようにします。
>summary(ssANOVA2.lmer)
--
Linear mixed model fit by REML
t-tests use Satterthwaite approximations to degrees of ...
Formula: dMFCC.SquareSum.log ~ seg2 * register + (1 + se...
Data: dataSumm.SS
REML criterion at convergence: 11133.4
Scaled residuals:
Min 1Q Median 3Q Max
-5.2300 -0.6120 0.1008 0.6875 3.1298
Random effects:
Groups Name Variance Std.De...
speakerID (Intercept) 0.0091745 0.0957...
seg2sonorant 0.0070485 0.0839...
seg2vowel 0.0003944 0.0198...
registerRe-Read 0.0033398 0.0577...
seg2sonorant:registerRe-Read 0.0021252 0.0461...
seg2vowel:registerRe-Read 0.0034361 0.0586...
Residual 0.2755449 0.5249...
Number of obs: 7153, groups: speakerID, 6
Fixed effects:
Estimate Std. Error d...
(Intercept) 3.53758 0.04069 4.9380...
seg2sonorant -0.22829 0.04481 4.3800...
seg2vowel -0.22922 0.02024 8.0650...
registerRe-Read 0.15976 0.02883 5.2340...
seg2sonorant:registerRe-Read -0.17606 0.04774 11.0610...
seg2vowel:registerRe-Read -0.09066 0.03639 6.3130...
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1...
Correlation of Fixed Effects:
(Intr) sg2snr sg2vwl rgsR-R sg2s:R-R
seg2sonornt -0.679
seg2vowel -0.291 0.222
registrR-Rd -0.804 0.717 0.295
sg2snrn:R-R 0.302 -0.651 -0.092 -0.477
sg2vwl:rR-R 0.486 -0.549 -0.567 -0.741 0.353
--もしくは、lmerTestパッケージのstep関数を使います。こち...
>step(ssANOVA2.lmer)
-お題３：''話者の個人差''及び''音素の種類''の影響を除外し...
--''目的変数（従属変数）''は dMFCC.SquareSum(.log) で、''...
---seg.length.log は連続変数、register は離散変数なので、...
--lmer関数を使ってモデルを構築します。
>ssANCOVA.lmer <- lmer(dMFCC.SquareSum.log ~ seg.length.l...
--結果は以下で確認できます。
>anova(ssANCOVA.lmer)
--
Analysis of Variance Table of type 3 with Satterthwaite
approximation for degrees of freedom
Sum Sq Mean Sq NumDF DenDF F.va...
seg.length.log 3.2454 3.2454 1 38.378 15.8...
register 0.7014 0.7014 1 18.591 3.4...
seg.length.log:register 0.0262 0.0262 1 25.829 0.1...
---seg.length.logの有意な影響が見られます。レジスタの影響...
--単純主効果は以下で確認できます。
>summary(ssANCOVA.lmer)
--
(略)
Fixed effects:
Estimate Std. Error ...
(Intercept) 3.72788 0.09773 32.01...
seg.length.log -0.15572 0.04095 36.29...
registerRe-Read 0.13573 0.07338 18.59...
seg.length.log:registerRe-Read -0.01233 0.03447 25.83...
(略)
***ロジスティック回帰分析 [#p9f2d132]
***ロジスティック混合効果モデル [#k9d83231]
終了行:
*音響解析（音声コーパス）データの統計解析 [#aa6f3f70]
#contents
**準備 [#x5be6c17]
***R と R Studio のインストール [#eeb2823d]
-&pgid(,R Note);
***サンプルデータ [#ce3453b5]
#ref(CSJ_MFCC26_short.zip);
→ &color(red){''これは話者数の少ないデータなので、下記の...
//-[[宮澤幸希，白勢彩子，馬塚れい子，菊池英明，”知識獲得...
-CSJの音声データファイルの冒頭120秒をMFCC解析した結果。各...
filename CSJのファイル名。12ファイルある。
register レジスタ名。"APS"は学会講演録音。"Re-Read"は学...
speakerID 話者名。APSとRe-Readについて同じ話者がいるため...
gender 話者の性別。
seg 音素名。SclSは<cl>、SsvSは<sv>。音素ラベルのないフ...
segID 音素ID。各ファイルについて音素ごとに数値をつけて...
frameID フレームID。各ファイルについて0.01秒ごとに数値...
beg_time フレームの開始時間。
MFCC1～MFCC12 MFCC12次元の解析結果。
C0 パワーの解析結果。
dMFCC1～dMFCC12 デルタMFCC12次元の解析結果。
dC0 デルタパワーの解析結果。
--ここでの「フレーム」は[[音声解析の分析窓:http://speechr...
--解析の詳細は &pgid(,MFCC解析のツール);/[[HTK(Hidden Mar...
***パッケージインストール [#k8de8a47]
-Rを起動して以下のコマンドを実行します。初回に一回だけ実...
>install.packages("plyr") #データテーブルの集計用パッ...
>install.packages("ggplot2") #視覚化のためのパッケージ
>install.packages("reshape") #データテーブルの整形用パッ...
>install.packages("lme4") #一般化線形モデルのためのパ...
>install.packages("lmerTest") #一般化線形モデルの検定のた...
>install.packages("pwr") #検定力解析（パワーアナリシ...
**データテーブルの読み込み [#a560c213]
-ダウンロードした &ref(CSJ_MFCC26_short.zip); を解凍し、...
-Rを起動して、以下を実行((R Studioであれば、右上の「Envir...
>data <- read.table("E:/Users/m-kouki/Desktop/CSJ_MFCC26....
--パスの区切り文字は「/」または「\\」
--今回のデータの1行目はヘッダ（各列のデータの詳細を説明し...
--参考：[[R-Source 40. ファイルからデータを読み込む:http:...
-読み込まれたデータは変数 data に格納されます。以下のよう...
>data #データテーブルの中身を見る
>head(data) #データテーブルの冒頭数行を見る
>names(data) #変数名の一覧を得る
-
[1] "filename" "register" "speakerID" "gender" "se...
[8] "beg_time" "MFCC1" "MFCC2" "MFCC3" "MF...
[15] "MFCC7" "MFCC8" "MFCC9" "MFCC10" "MF...
[22] "dMFCC1" "dMFCC2" "dMFCC3" "dMFCC4" "dM...
[29] "dMFCC8" "dMFCC9" "dMFCC10" "dMFCC11" "dM...
-データの一部分を指定することができます。
>data[ , 1:4] #1列目から4列目までの全ての行
>data[ 1:4 ,] #1行目から4行目までの全ての列
**データの前処理 [#l688e65c]
***要約統計量を得る [#u5c296b9]
-以下を実行します。
>summary(data)
--[[データの尺度:http://speechresearch.fiw-web.net/114.ht...
---'''連続型変数 Area, Slope, Soil.pH, Worm.density は6つ...
--カテゴリカル変数のつもりで作成したデータが、summaryを表...
speakerID
Min. :1.000
1st Qu.:2.000
Median :4.000
Mean :3.596
3rd Qu.:5.000
Max. :6.000
--実際には、speakerIDの数値はID番号なので、平均値や中央値...
>data$speakerID <- factor(data$speakerID)
>data$segID <- factor(data$segID)
>data$frameID <- factor(data$frameID)
--もう一度summaryを確認します。
>summary(data)
-データの総数は以下で求められます。
>nrow(data)
***特定のデータを抜き出す [#lacf4d9c]
-単に特定の列だけ抜き出したい場合、例えば以下のように書け...
>data[ c("speakerID","MFCC1") ]
-条件に該当する行だけ抜き出したい場合、例えば以下のような...
>data[ C0 > 20 & gender == "F" , 1:9] # C0の値が20以上の...
-subset関数を使うやり方が便利。[[Rデータフレームから一部...
>subset(data, (C0 > 20) & (gender == "F"), select = 1:9 )
--ヘッダ名を指定して取り出すこともできる。
>subset(data, (C0 > 20) & (gender == "F"), select = c(spe...
-特定のデータのみを抜き出して新しいデータテーブルを作るこ...
>dataF <- subset(data, (gender == "F")) # 女性のみを取り...
>summary(dataF)
-欠損値を除外したい場合は na.omit 関数を使う。
>data <- na.omit(data)
--どれか一つの列に欠損値を含んでいた場合、その行ごと削除...
-複数の条件で抜き出したデータテーブルを結合する場合は、rb...
>dataV <- rbind( subset(data, seg=="i"), subset(data, se...
>summary(dataV)
--特定のデータを抜き出した後のデータテーブルをsummaryする...
seg
a :3519
o :3451
u :3304
i :1719
e :1514
aH : 0
(Other): 0
--再度factor関数にかけてやると、カウント数0のデータは消え...
>dataV$seg <- factor(dataV$seg)
>summary(dataV)
***特定のデータの計算結果を追記する [#e3780abf]
-MFCC1～MFCC12の平均を求めて、データテーブルに"MFCC.Mean"...
>data$MFCC.Mean <- (data$MFCC1 + data$MFCC2 + data$MFCC3 ...
--apply関数を使うこともできます。
>data$MFCC.Mean <- apply(data[, 9:20], 1, mean)
---各行について集計したい場合は第二引数に1、各列について...
--値のどこかにNaNを含む可能性があって、その値を除外して平...
>data$MFCC.Mean <- apply(data[, 9:20], 1, mean, na.rm=TRUE)
-dMFCC1～dMFCC12の二乗和((この値は音響的特徴の変動量の大...
>data$dMFCC.SquareSum <- data$dMFCC1^2 + data$dMFCC2^2 + ...
>summary(data)
-対数変換したデータを追記する
>data$dMFCC.SquareSum.log = log(data$dMFCC.SquareSum)
-C0を平均値が0、分散が1になるように標準化したC0.scale列を...
>data$C0.scale = scale(data$C0)
>summary(data)
--列名は C0.scale.V1 になっている。→ &color(red){''理由は...
-【発展】条件を指定してデータを追記するには、transform 関...
>data = transform(data, seg2 = ifelse( ( (seg == "i") | (...
--さらに条件を複雑にして、「短母音5種類には"sV"、長母音5...
>data = transform(data, seg3 = ifelse( ( (seg == "i") | (...
-【発展】特定のデータを含む行に隣接する行のデータを得る
--例えば、「短母音の行の直前と直後の行を探して、該当する...
--データテーブルの基本的な考え方は、「''1行の中に解析に必...
---まずは、seg列のみを取り出して新しいデータテーブルdata....
>data.seg = data[ c("seg") ]
---空欄を意味する1行1列のデータテーブルを作ります。
>seg.nothing = data.frame(seg=c("Nothing"))
---data.segを前に1列ずらしたデータテーブル data.seg.prev ...
>data.seg.prev = rbind(seg.nothing, data.seg)
>data.seg.next = rbind(data.seg, seg.nothing)
---続いて、data.seg.prev の最初の列、 data.seg.next の最...
>data.seg.prev = data.seg.prev[1:(nrow(data.seg.prev)-1) ...
>data.seg.next = data.seg.next[2:nrow(data.seg.next) , ]
---元のデータテーブル data に data.seg.prev と data.seg.n...
> data$seg.prev = data.seg.prev
> data$seg.next = data.seg.next
> head(data)
---これで、必要な情報は全て1行の中に集約されました。
--transform関数を使って該当条件のラベルを付加します。
>data = transform(data, seg.neighbor = ifelse( ( (seg.pre...
--このサンプルデータは各行が「フレーム（分析窓）」単位で...
***特定のデータの集計値を得る [#p0a9dcd9]
-ddply関数が大変便利です((浅井君、教えてくれてありがとう)...
-[[plyrパッケージのインストール:http://speechresearch.fiw...
>library(plyr)
-様々な実例は [[「plyrパッケージで君も前処理スタ☆」改め「...
-各話者・各レジスタごとに、C0の平均値と標準偏差を求めて、...
>dataSumm <- ddply(data, c('speakerID', 'register'), summ...
>dataSumm
--指定した条件ごとの集計値が出力されます。
speakerID register C0mean C0sd
1 1 APS 15.69648 2.949087
2 1 Re-Read 13.24777 2.565976
3 2 APS 15.75885 2.861640
4 2 Re-Read 13.95841 2.957040
5 3 APS 14.53344 2.584729
6 3 Re-Read 14.81330 2.511646
7 4 APS 13.91550 2.801826
8 4 Re-Read 13.18621 3.143895
9 5 APS 15.02620 2.385534
10 5 Re-Read 14.91980 3.005093
11 6 APS 15.10699 2.306960
12 6 Re-Read 14.23184 2.976182
-各話者・各レジスタごとに、C0の平均値と標準偏差を求めて、...
>data <- ddply(data, c('speakerID', 'register'), transfor...
>data
--元のデータに列が追加されています。
-各話者について、C0の値の平均値が0、標準偏差が1になるよう...
>data <- ddply(data, c('speakerID'), transform, C0.z = (C...
>head(data)
-【発展】特定条件の集計値を全行に追加する
--例えば「レジスタAPSに限定して、話者ごとにC0の平均値を求...
>#レジスタAPS以外のC0値をNAにした新しい列「C0_APS」を作り...
>data = transform(data, C0_APS = ifelse( register == "AP...
>#各話者のAPSのC0の平均値を求めて、すべての行に追記する
>data = ddply(data, c('speakerID'), transform, C0_APSmean...
-【発展】フレーム単位のデータテーブルから音素単位のデータ...
--サンプルデータ data の各行は、窓長0.01秒の時間フレーム...
--segIDについて集計値を求めます。segIDは異なる話者で重複...
> data2 <- ddply(data, c('speakerID', 'segID'), summarize...
> nrow(data2)
--
[1] 7153
---他のfactor変数を条件に加えても構いません。データの総行...
> data2 <- ddply(data, c('filename', 'register', 'speaker...
> nrow(data2)
---
[1] 7153
--よく行われる処理は、各音素について時間フレームの中心（...
> data2 <- ddply(data, c('speakerID', 'segID'), summarize...
---length(MFCC1) は各条件で抽出されたMFCC1のベクトルサイ...
---round(length(MFCC1)/2) は各音素のフレーム長の半分を四...
---MFCC1[round(length(MFCC1)/2)] は各音素の時間的中心点の...
***外れ値を除外する [#ma6808aa]
-ddply関数を使って、「話者＆レジスタごとにC0の平均値と標...
-まず、平均値＋2.0SD を上回る値は1、そうでない値は0をつけ...
>data <- ddply(data, c('speakerID', 'register'), transfor...
-続いて、平均値－2.0SD を下回る値は1、そうでない値は0をつ...
>data <- ddply(data, c('speakerID', 'register'), transfor...
-C0out1またはC0out2が1になっている行ははずれ値として削除...
>data2 = subset(data, C0out1=="0")
>data2 = subset(data2, C0out2=="0")
>nrow(data)
>nrow(data2)
***ランダムサンプリング [#x9fff2b3]
-ddply関数を使えば、指定した各条件について、該当するデー...
-話者＆レジスタ＆各音素カテゴリーごとに、C0とdC0をそれぞ...
>sampling = 10
>dataRand = ddply(data, c('speakerID', 'register', 'seg')...
>dataRand
--サンプリングしたC0とdC0の行番号は一致していない（それぞ...
--"replace = T"は、データが重複してサンプリングされること...
**データの視覚化 [#qcf6ad0a]
-データの特性を理解するために、視覚化はとても重要です。
***準備 [#add37ba8]
-ggplot2関数が大変多機能で、便利です((浅井君、教えてくれ...
-[[ggplot2パッケージのインストール:http://speechresearch....
>library(ggplot2)
>library(plyr)
>library(reshape)
-サンプルデータ &ref(CSJ_MFCC26_short.zip); をダウンロー...
-以後の説明では、短母音のみのデータテーブルを使います。
>data <- read.table("E:/Users/m-kouki/Desktop/CSJ_MFCC26....
>data$speakerID = factor(data$speakerID)
>data$segID = factor(data$segID)
>data$frameID = factor(data$frameID)
>dataV = subset(data, (seg=="i")|(seg=="e")|(seg=="o")|(s...
>summary(dataV)
***散布図 [#madd8730]
-単にplot関数に１つのデータ列を与えると、散布図になります。
>plot(dataV$C0)
#ref(Rplot_1.png);
--縦軸がC0の値、横軸は左からデータの行数です。
--ここで、分析に使う各パラメータの散布図を確認しておくと...
-plot関数に２つのデータ列を与えてみます。
>plot(dataV$MFCC2, dataV$MFCC4)
#ref(Rplot_2.png);
--ここでは、縦軸がMFCC4の値、横軸がMFCC2の値です。
--女性話者のデータに限定してみます。
>plot(subset(dataV, gender=="F")$MFCC2, subset(dataV, gen...
-ggplotを使うと、条件ごとの色分けも簡単にできます。以下は...
>ggplot(dataV) +
> geom_point(aes(x=MFCC2, y=MFCC4, color=seg)) +
> labs(title="All")
#ref(Rplot_3.png);
--女性話者のデータに限定してみます。後でグラフに機能を追...
>p = ggplot(subset(dataV, gender=="F")) +
> geom_point(aes(x=MFCC2, y=MFCC4, color=seg)) +
> labs(title="Female")
>p
--グラフの軸のサイズを指定することもできます。上記の変数 ...
>p = p + xlim(-30, 30) + ylim(-30, 30)
>p
#ref(Rplot_4.png);
-ddplyと組み合わせると色々と便利です。
--母音ごと・レジスタごとに平均値を求めて、レジスタで色分...
>dataSumm = ddply(dataV, c('seg', 'register'), summarize,...
>ggplot(dataSumm) +
> geom_text(aes(x=MFCC2mean, y=MFCC4mean, label=seg, col...
---散布図の各点を文字にしたい場合は、geom_textを使います...
#ref(Rplot_5.png);
--話者ごと・母音ごと・レジスタごとにMFCC2, MFCC4, MFCC6, ...
>dataSumm = ddply(dataV, c('speakerID', 'seg', 'register'...
>pairs(~dataSumm$MFCC2mean + dataSumm$MFCC4mean + dataSum...
#ref(Rplot_6.png);
---これはggplotを使っていません。((ggplotで同様のことはpl...
***ヒストグラム [#qcff188b]
-分析の最初のチェックとして、ヒストグラムの形状は非常に重...
--分析に用いる各要因の各水準ごとにヒストグラムをプロット...
--単峰性でも、ヒストグラムの形状がゆがんでいる場合、注意...
-全てのデータのMFCC1のヒストグラムを、通常の方法でプロット
>hist(dataV$MFCC1)
-ggplotを使用した方法
>ggplot( dataV, aes(x = MFCC1) ) + geom_histogram()
#ref(Rplot_7.png);
--全データのMFCC1は２峰性の分布になっています。各要因・各...
-ggplotとsubplotを組み合わせると、各要因・各水準のヒスト...
--話者別にヒストグラムをプロット
>ggplot(dataV) +
> geom_histogram( aes(x = MFCC1) ) +
> facet_wrap(~speakerID)
--レジスタ・話者・母音別にヒストグラムをプロット
>ggplot(dataV) +
> geom_histogram( aes(x = MFCC1) ) +
> facet_wrap(register~speakerID~seg)
--レジスタで色分けして、話者を縦グリッド・母音を横グリッ...
>ggplot(dataV) +
> geom_histogram( aes( x = MFCC1, fill=factor(register) ...
> facet_grid(speakerID~seg)
#ref(Rplot_8.png);
---/u/では全ての話者で二峰性の分布になっていることが分か...
-多次元ヒストグラム
--[[R generate 2D histogram from raw data:http://stackove...
--[[二変量データのヒストグラムと密度推定:http://www.okada...
***箱ヒゲ図（ボックスプロット）(([[統計学：Rを用いた入門...
-t検定や分散分析の前に、データの分布や有意傾向を視覚的に...
-全てのデータのMFCC1の箱ヒゲ図を、通常の方法でプロット
>boxplot(subset(dataV, select = c(MFCC1)))
--複数の列の値を同時にプロットすることもできます。
>boxplot(subset(dataV, select = c(MFCC1, MFCC2, MFCC3) ) )
#ref(Rplot_9.png);
--箱ヒゲ図の見方は以下の通りです。
(([[R+ 箱ひげ図の読み方:http://wb-nahce.info/rplus/statse...
---中央の横線：各データの全体の中央値（第2クォンタイル((...
---箱の上端の線：全体の中央値より大きいデータの中央値（第...
---箱の下端の線：全体の中央値より小さいデータの中央値（第...
---箱の中に、データ全体の約50%が含まれます。上下のヒゲは...
---ヒゲの外側にある値は、はずれ値として扱われます。
-全てのデータのMFCC1の箱ヒゲ図を、ggplotでは以下のように...
>ggplot(dataV) +
> geom_boxplot(aes(y = MFCC1, x = "MFCC1"))
-ggplotは、boxplotと特性が異なります。boxplotは異なる列の...
--レジスタで色分けして、母音ごとにMFCC1の値をプロット
>ggplot(dataV) +
> geom_boxplot(aes(y = MFCC1, x = seg, fill = register))
---ノッチを加えてプロット(('''この図にてへこんだ部分がノ...
>ggplot(dataV) +
> geom_boxplot(aes(y = MFCC1, x = seg, fill = register),...
#ref(Rplot_10.png);
--箱ヒゲ図は分布のばらつきを見るのに適していますが、''分...
>ggplot(dataV) +
> geom_violin(aes(y = MFCC1, x = seg, fill = register))
#ref(Rplot_15.png);
--レジスタで色分けして、話者ごとにMFCC1の値をプロットし、...
>ggplot(dataV) +
> geom_boxplot(aes(y = MFCC1, x = speakerID, fill = regi...
> facet_wrap( ~ seg)
#ref(Rplot_11.png);
-ggplotで、異なる列のデータを並べてプロットしたい場合は、...
>dataV.melt <- melt(subset(dataV, select = c(register, sp...
>head(dataV.melt)
--赤字のメッセージが出ますが問題ありません。dataV.meltの...
--レジスタで色分けして、MFCC1,MFCC2,MFCC3の値をプロットし...
>ggplot(dataV.melt) +
> geom_boxplot(aes(y = value, x = variable, fill = regis...
> facet_grid(speakerID ~ seg)
#ref(Rplot_12.png);
---この例では細かすぎて見づらくなっていますが、違いが出る...
***棒グラフ [#i5d11244]
-barplot関数を使った方法：[[バイオスタティスティクス 棒グ...
-ggplotを使う場合、まずはddplyを使ってデータを集計し、平...
--今回は、「各話者・各母音・各レジスター」別にMFCC2とMFCC...
--「各話者・各母音・各レジスター」別にMFCC2とMFCC4の平均...
>dataVSumm1 = ddply(dataV, c('speakerID', 'seg', 'registe...
>dataVSumm1
--reshapeパッケージのmelt関数を使って、データを並び替えま...
>dataVSumm1M = melt(dataVSumm1)
>dataVSumm1M
---%%エラーバーを表示しないなら、この時点でプロットできま...
>qplot(register, weight=value, geom="bar", data=dataVSumm...
---エラーバーを表示しないなら、その値の「各母音・各レジス...
>dataVSumm2 = ddply(dataVSumm1, c('seg', 'register'), sum...
>dataVSumm2
>dataVSumm2M = melt(dataVSumm2)
>dataVSumm2M
>
>qplot(register, weight=value, geom="bar", data=dataVSumm...
#ref(Rplot_13.png);
--エラーバーを表示する場合、その値の「各母音・各レジスタ...
>dataVSumm3 = ddply(dataVSumm1M, c('seg', 'register', 'va...
>dataVSumm3
---length(value)は、集計するデータサイズ。ここでは話者数...
--ggplot関数を使って棒グラフとエラーバーを表示します（参...
>ggplot(dataVSumm3, aes(x=variable, y=mean, fill=register...
> geom_bar(stat="identity", position="dodge") +
> geom_errorbar(aes(ymax=mean+se, ymin=mean-se), positio...
> facet_grid( ~ seg)
---"aes"の設定はgeom_barでもgeom_errorbarでも共通して使う...
#ref(Rplot_14.png);
---fillとfacetの設定を変えてプロットの形式を変更します。
--
>ggplot(dataVSumm3, aes(x=variable, y=mean, fill=seg)) +
> geom_bar(stat="identity", position="dodge") +
> geom_errorbar(aes(ymax=mean+se, ymin=mean-se), positio...
> facet_grid( ~ register)
#ref(Rplot_14a.png);
**統計解析 [#n2090f8b]
***準備 [#n9f75954]
-[[lme4パッケージとlmerTestパッケージのインストール:http:...
>library(lme4)
>library(lmerTest)
>library(plyr)
>library(reshape)
>library(pwr)
-サンプルデータ &ref(CSJ_MFCC26_short.zip); をダウンロー...
>data <- read.table("E:/Users/m-kouki/Desktop/CSJ_MFCC26....
>data$speakerID = factor(data$speakerID)
>data$segID = factor(data$segID)
>data$frameID = factor(data$frameID)
***t検定とパワーアナリシス [#l8aa50bb]
-対応のないt検定
--全てのデータのMFCC5の平均値とMFCC6の平均値の間に、有意...
--まずは箱ヒゲ図を書いてみます。
>boxplot(subset(data, select = c(MFCC5, MFCC6) ) )
#ref(Rplot_16.png);
---ほぼボックスが重複しているように見えますが、ノッチが極...
--F検定で等分散性を見てみます（Rではt.test関数で自動的に...
>var.test(data$MFCC5, data$MFCC6)
--
F test to compare two variances
data: data$MFCC5 and data$MFCC6
F = 1.017, num df = 40137, denom df = 40137, p-value = 0...
(略)
---p値は0.09なので、ぎりぎりですがMFCC5とMFCC6は「等分散...
--t検定を実施します。ここでは等分散であることを明記します。
>t.test(data$MFCC5, data$MFCC6, var.equal=TRUE)
--
Two Sample t-test
data: data$MFCC5 and data$MFCC6
t = 5.3208, df = 80274, p-value = 1.036e-07
(略)
---p値は0.001より小さいので、MFCC5とMFCC6の平均値は有意に...
--var.equalオプションなしで実行した場合、等分散であればt...
>t.test(data$MFCC4, data$MFCC9)
--
Welch Two Sample t-test
data: data$MFCC4 and data$MFCC9
t = 1.3549, df = 77047.3, p-value = 0.1754
(略)
---p値が0.05より大きいので、MFCC4とMFCC9の平均値は有意に...
-【注】t検定で比較する2群の標本数が大きく違う場合は、ウェ...
--上のデータの検定力を求めてみます。参考：[[効果量・検定...
--今回のt検定で使用した２つのデータサイズ確認
>length(data$MFCC4)
--
[1] 40138
--
>length(data$MFCC9)
--
[1] 40138
--検定力を求める
>pwr.t2n.test(n1=40138, n2=40138, sig.level=0.05, d=0.8)
---n1=40138 : データ群1のデータサイズ（サンプル数）です。
---n2=40138 : データ群2のデータサイズ（サンプル数）です。...
---d=0.8 : Cohenの効果量d には 0.8 を指定しています。これ...
---sig.level=0.05 : 有意水準0.05で統計を行う場合を想定し...
--結果は以下の通りです。
t test power calculation
n1 = 40138
n2 = 40138
d = 0.8
sig.level = 0.05
power = 1
alternative = two.sided
---power は 1 になりました。すなわち100回中統計を実施した...
---power の値は 0.8 以上であれば（100回中80回同じ結果にな...
--有意水準と検定力の関係について：[[例数設計の基礎 第８回...
-対応のあるt検定
--上でMFCC4とMFCC9の平均値に有意差がありませんでしたが、...
--そこで、話者ごとにMFCC4とMFCC9の平均値を求めて、その値...
> dataSumm1 <- ddply(data, c('speakerID'), summarize, MFC...
> dataSumm1
--
speakerID MFCC4 MFCC9
1 1 -0.36694394 -0.67979122
2 2 -0.02150908 0.28991697
3 3 0.23361956 -1.49124097
4 4 0.12591752 0.05997551
5 5 -0.97832807 0.45119489
6 6 -0.09088042 -0.02595399
--このデータは、各行が同じ話者から得られたデータなので、...
> t.test(dataSumm1$MFCC4, dataSumm1$MFCC9, paired = TRUE)
--
Paired t-test
data: dataSumm1$MFCC4 and dataSumm1$MFCC9
t = 0.1191, df = 5, p-value = 0.9098
---p値が0.05より大きいので、MFCC4とMFCC9の平均値は有意に...
--ただし、今回のデータは、''集計前よりもあまりにも数が少...
---ある数のデータセットから得られた統計の結果が信頼のでき...
--もう少し対応のあるデータを増やします。レジスタ・音素・...
> dataSumm2 <- ddply(data, c('speakerID', 'register', 'se...
> dataSumm2
--これでデータ数は454個になりました。対応のあるt検定を実...
> t.test(dataSumm2$MFCC4, dataSumm2$MFCC9, paired = TRUE)
--
Paired t-test
data: dataSumm2$MFCC4 and dataSumm2$MFCC9
t = 2.1648, df = 453, p-value = 0.03093
---p値は0.001より小さくなりました。したがって、MFCC4とMFC...
--このように、特に自然発話データの場合、どのような集計を...
***単回帰分析 [#tc65193f]
-お題：''各音素の持続時間（フレーム数）''と、''各音素のdM...
--音素が長いほど、発音に使える時間に余裕があるので、各音...
--分析に使用するデータテーブル dataSumm.SS を作ります。各...
>data$dMFCC.SquareSum <- data$dMFCC1^2 + data$dMFCC2^2 + ...
>dataSumm.SS <- ddply(data, c('speakerID', 'segID'), summ...
>head(dataSumm.SS)
--
speakerID segID seg.length dMFCC.SquareSum
1 1 4627 4 54.75966
2 1 4628 2 55.05001
3 1 4629 3 45.30628
4 1 4630 6 46.35851
5 1 4631 6 41.52738
6 1 4632 13 39.62246
-ここで、''目的変数（従属変数）''は dMFCC.SquareSum で、'...
--単回帰分析では、目的・説明の両変数は正規分布する連続変...
>par(mfrow=c(1,2))
>hist(dataSumm.SS$seg.length)
>hist(dataSumm.SS$dMFCC.SquareSum)
#ref(Rplot_17.png);
--どちらも、L字型の分布になっています。持続時間や二乗和の...
>dataSumm.SS$seg.length.log = log(dataSumm.SS$seg.length)
>dataSumm.SS$dMFCC.SquareSum.log = log(dataSumm.SS$dMFCC....
>par(mfrow=c(1,2))
>hist(dataSumm.SS$seg.length.log)
>hist(dataSumm.SS$dMFCC.SquareSum.log)
#ref(Rplot_18.png);
---これで正規分布に近づきました。
--続いて、散布図を描いてみます。
>par(mfrow=c(1,1))
>plot( dataSumm.SS$seg.length.log, dataSumm.SS$dMFCC.Squa...
#ref(Rplot_19.png);
--分布が縞々になってしまいました。これは、seg.length が整...
--相関分析にかけたところ、有意なある程度の負の相関（p-val...
>cor.test( dataSumm.SS$seg.length.log, dataSumm.SS$dMFCC....
--lm関数を使って回帰分析にかけます。
>ss.lm <- lm(dMFCC.SquareSum.log ~ seg.length.log, data =...
>summary(ss.lm)
---lm関数の引数で回帰式を指定します。
y(目的変数) = a(切片) + b(傾き)x(説明変数) → y~x
y = bx → y~-1+x
---結果は以下の通りです。
--
Call:
lm(formula = dMFCC.SquareSum.log ~ seg.length.log, data ...
Residuals:
Min 1Q Median 3Q Max
-3.00788 -0.29871 0.07117 0.36034 1.36548
---Residuals（残差）は、データと回帰モデルが予測したの値...
--
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.882726 0.014611 265.74 <2e-16 ***
seg.length.log -0.251529 0.008559 -29.39 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1...
Residual standard error: 0.512 on 6693 degrees of freedom
---(Intercept) は切片a、seg.length.log は説明変数の傾き（...
---目的変数に対して、説明変数 seg.length.log は有意な負の...
--
Multiple R-squared: 0.1143, Adjusted R-squared: 0.1142
F-statistic: 863.7 on 1 and 6693 DF, p-value: < 2.2e-16
--さきほどの図に回帰直線を追記します。
>abline(ss.lm)
--回帰診断図を描きます。
>par(mfrow=c(2,2))
>plot(ss.lm)
***一元配置分散分析＋多重比較法 [#i12035e7]
&color(red){''この節の手順は誤りがある可能性があります。...
-お題：''dMFCC1～dMFCC12の二乗和は、母音・ソノラント子音...
-３群以上の間の有意差の比較は、[[t検定:http://speechresea...
--dMFCC1～dMFCC12の二乗和を求め、[[ここ:http://speechrese...
>data$dMFCC.SquareSum <- data$dMFCC1^2 + data$dMFCC2^2 + ...
>data$dMFCC.SquareSum.log <- log(data$dMFCC.SquareSum)
>data = transform(data, seg2 = ifelse( ( (seg == "i") | (...
>data$seg2 <- factor(data$seg2)
>summary(data)
--箱ヒゲ図を表示してみます。
>data.melt.aov1 <- melt(subset(data, select = c(seg2, dMF...
>library(ggplot2)
>ggplot(data.melt.aov1) + geom_boxplot(aes(y = value, x =...
#ref(Rplot_20.png);
-対応のない一元配置分散分析
--関数aovを使って分散分析を実行します。
>ss.aov1 <- aov(dMFCC.SquareSum.log ~ seg2, data = data)
>ss.aov1
--
Call:
aov(formula = dMFCC.SquareSum.log ~ seg2, data = data)
Terms:
seg2 Residuals
Sum of Squares 1032.00 33207.52
Deg. of Freedom 2 40135
Residual standard error: 0.9096128
Estimated effects may be unbalanced
--分散分析表を見るには、summaryもしくはanova関数を使いま...
>summary(ss.aov1)
--
Df Sum Sq Mean Sq F value Pr(>F)
seg2 2 1032 516.0 623.6 <2e-16 ***
Residuals 40135 33208 0.8
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1...
---意味は以下の通りです。((分散分析表の枠は[[一元配置の分...
-------------------------------------------------------
要因 自由度 平方和 平均平方和 F-値 P-値
-------------------------------------------------------
F 2 1032 516.0 623.6 <2e-16 ***
誤差 40135 33208 0.8
-------------------------------------------------------
---P値が0.05より小さいため、説明変数（seg2）の各水準（vow...
--この結果だけでは、どの水準とどの水準に差があったのかは...
>TukeyHSD(ss.aov1)
--
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = dMFCC.SquareSum.log ~ seg2, data = da...
$seg2
diff lwr upr p...
sonorant-non-sonorant -0.1880719 -0.2295339 -0.1466100 ...
vowel-non-sonorant -0.3335195 -0.3556718 -0.3113673 ...
vowel-sonorant -0.1454476 -0.1870469 -0.1038483 ...
---diff : 平均値の差、lwr : 信頼区間の下限、upr : 信頼区...
---p値がゼロになってしまった。→ &color(red){''理由は？''};
--[[R言語で統計解析入門： くり返しのある一元配置分散分析...
>pairwise.t.test( data$dMFCC.SquareSum.log, data$seg2, p....
--
non-sonorant sonorant
sonorant <2e-16 -
vowel <2e-16 <2e-16
---各組み合わせについて有意差がある。
-対応のある一元配置分散分析
--レジスタ・話者別、また音素の種類ラベル（説明変数）別に...
> library(plyr)
> dataSumm.aov1p <- ddply(data, c('speakerID', 'register'...
> dataSumm.aov1p
--話者について対応をとり分散分析にかける場合、以下のよう...
> ss.aov1p.speaker <- aov(dMFCC.SquareSum.log ~ seg2 + sp...
> summary(ss.aov1p.speaker)
--
Df Sum Sq Mean Sq F value Pr(>F)
seg2 2 0.6355 0.3177 24.398 7.34e-07 ***
speakerID 5 0.1842 0.0368 2.829 0.0344 *
Residuals 28 0.3646 0.0130
(略)
---speakerIDの有意な影響が認められます。
---多重比較を行います。
--
> TukeyHSD(ss.aov1p.speaker)
--
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = dMFCC.SquareSum.log ~ seg2 + speakerI...
$seg2
diff lwr upr ...
sonorant-non-sonorant -0.1950974 -0.3103730 -0.07982184 ...
vowel-non-sonorant -0.3231252 -0.4384008 -0.20784963 ...
vowel-sonorant -0.1280278 -0.2433034 -0.01275223 ...
(略)
---各組み合わせについて有意差がある。
--レジスタについて対応をとり分散分析にかける場合、以下の...
> ss.aov1p.register <- aov(dMFCC.SquareSum.log ~ seg2 + s...
> summary(ss.aov1p.register)
--
Df Sum Sq Mean Sq F value Pr(>F)
seg2 2 0.6355 0.3177 19.276 3.21e-06 ***
register 1 0.0214 0.0214 1.299 0.263
Residuals 32 0.5275 0.0165
---registerの有意な影響は認められません。
---多重比較を行います。
--
>TukeyHSD(ss.aov1p.register)
--
(略)
$seg2
diff lwr up...
sonorant-non-sonorant -0.1950974 -0.3238975 -0.066297298...
vowel-non-sonorant -0.3231252 -0.4519253 -0.194325088...
vowel-sonorant -0.1280278 -0.2568279 0.000772313...
(略)
---registerの有意な影響がないにもかかわらず、レジスタにつ...
--speakerIDとregisterを結合して一つのラベル（speakerIDreg...
>dataSumm.aov1p$speakerIDregister = paste(dataSumm.aov1p$...
>dataSumm.aov1p
>ss.aov1p.speakerIDregister <- aov(dMFCC.SquareSum.log ~ ...
>summary(ss.aov1p.speakerIDregister)
>TukeyHSD(ss.aov1p.speakerIDregister)
---結果は省略します。
***二元配置分散分析＋交互作用効果の扱い [#e8bbb4b5]
&color(red){''この節の手順は誤りがある可能性があります。...
-お題：''dMFCC1～dMFCC12の二乗和は、母音・ソノラント子音...
--説明変数（独立変数）が２種類になりました。第１要因は「...
--データは一元配置分散分析と同じものを使います。
>data$dMFCC.SquareSum <- data$dMFCC1^2 + data$dMFCC2^2 + ...
>data$dMFCC.SquareSum.log <- log(data$dMFCC.SquareSum)
>data = transform(data, seg2 = ifelse( ( (seg == "i") | (...
>data$seg2 <- factor(data$seg2)
--交互作用効果をプロットしてみます。
>interaction.plot(data$register, data$seg2, data$dMFCC.Sq...
#ref(Rplot_21.png);
---線が平行でないので、２つの説明変数の間に交互作用があり...
-対応のない二元配置分散分析
--関数aovを使って分散分析を実行します。
> ss.aov2 <- aov(dMFCC.SquareSum.log ~ seg2 * register, d...
> summary(ss.aov2)
---「seg2 * register」は、「音素の種類、レジスタ、音素の...
--
Df Sum Sq Mean Sq F value Pr(>F)
seg2 2 1034 517.2 627.66 <2e-16 ***
register 1 61 60.7 73.68 <2e-16 ***
seg2:register 2 74 37.0 44.96 <2e-16 ***
Residuals 40132 33070 0.8
---交互作用効果は「seg2:register」と表現されます。Pr <2e-...
--交互作用効果が確認されたときは、''単純（主）効果の検定'...
---交互作用効果が有意であった場合、主効果について検討する...
---交互作用効果が有意でなければ、Tukeyの多重比較法だけでO...
---参考サイト：[[単純主効果の検定と多重比較:http://igproj...
---参考情報：'''[[R Note/統計/分散分析:http://speechresea...
--以下、作成中
***混合効果モデル：独立変数が一つかつ連続変数 [#j49fc621]
-混合効果（単）回帰分析
-参考：[[Mixed Effects Models Blog:http://mixedmodeljp.bl...
-お題：''話者の個人差''及び''音素の種類''の影響を除外して...
--[[単回帰分析:http://speechresearch.fiw-web.net/121.html...
--分析に使用するデータテーブル dataSumm.SS を作ります。各...
>data$dMFCC.SquareSum <- data$dMFCC1^2 + data$dMFCC2^2 + ...
>dataSumm.SS <- ddply(data, c('speakerID', 'seg', 'segID'...
>head(dataSumm.SS)
--
speakerID seg segID seg.length dMFCC.SquareSum
1 1 a 5169 9 59.40295
2 1 a 5213 7 24.76156
3 1 a 5311 10 27.79712
4 1 a 5325 6 22.73648
5 1 a 5331 8 19.03984
6 1 a 5458 10 46.80528
---ここで、''目的変数（従属変数）''は dMFCC.SquareSum で...
---&color(red){''重要''：ランダム要因はランダムサンプリン...
--目的・説明の両変数は対数変換します。
>dataSumm.SS$seg.length.log = log(dataSumm.SS$seg.length)
>dataSumm.SS$dMFCC.SquareSum.log = log(dataSumm.SS$dMFCC....
--lmer関数を使ってモデルを構築します。
>ss.lmer <- lmer(dMFCC.SquareSum.log ~ seg.length.log + (...
---従属変数 ~ 独立変数 + (ランダム変数) のように書く。（...
---（）の中の｜の左側の数字１は切片を示す。「(1 + seg.len...
--結果は以下の通りです。&color(red){lmer関数の実行結果は...
>summary(ss.lmer)
--
Linear mixed model fit by REML
t-tests use Satterthwaite approximations to degrees of ...
Formula: dMFCC.SquareSum.log ~ seg.length.log + (1 + seg...
Data: dataSumm.SS
REML criterion at convergence: 9438.4
Scaled residuals:
Min 1Q Median 3Q Max
-5.7769 -0.5883 0.0893 0.6904 3.0492
Random effects:
Groups Name Variance Std.Dev. Corr
seg (Intercept) 0.224867 0.47420
seg.length.log 0.038383 0.19592 -0.76
speakerID (Intercept) 0.007856 0.08864
seg.length.log 0.001212 0.03481 -0.58
Residual 0.213650 0.46222
Number of obs: 7101, groups: seg, 52; speakerID, 6
Fixed effects:
Estimate Std. Error df t value Pr(>...
(Intercept) 3.76566 0.08751 35.49000 43.033 < 2...
seg.length.log -0.15318 0.03945 37.86000 -3.883 0.00...
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1...
Correlation of Fixed Effects:
(Intr)
sg.lngth.lg -0.770
---各表示の詳細は[[Mixed Effects Models Blog - R言語を用...
---「Fixed effects:」以下が固定因子の効果量。(Intercept) ...
-文献によっては、「ランダム効果は切片を変数にする場合が多...
--このようなモデルは以下のようになります。
>ss2.lmer <- lmer(dMFCC.SquareSum.log ~ seg.length.log + ...
>summary(ss2.lmer)
--
(略)
REML criterion at convergence: 9542.9
(略)
Fixed effects:
Estimate Std. Error df t value ...
(Intercept) 3.74862 0.05653 41.00000 66.31 ...
seg.length.log -0.17120 0.01106 6872.00000 -15.48 ...
(略)
--ランダム変数が切片と傾きに与える影響を考慮したモデル（s...
>anova(ss.lmer, ss2.lmer)
--
Models:
..1: dMFCC.SquareSum.log ~ seg.length.log + (1 | speaker...
..1: seg)
object: dMFCC.SquareSum.log ~ seg.length.log + (1 + seg....
object: speakerID) + (1 + seg.length.log | seg)
Df AIC BIC logLik deviance Chisq Chi Df P...
..1 5 9541.8 9576.1 -4765.9 9531.8 ...
object 9 9447.8 9509.6 -4714.9 9429.8 101.95 4 ...
---ここで「..1」は ss2.lmer、「object」は ss.lmer です。
---まず、P値（Pr）に着目します。P値が0.05より小さいので、...
---P値が0.05より大きい（両モデルに有意な差がない）場合、A...
---対数尤度を使ったより厳密なモデル評価の方法は [[Mixed E...
-''まずは、できるだけ情報量の多いモデルを採用するべきです...
--Model failed to converge: degenerate Hessian with ○ ne...
--Model is nearly unidentifiable: very large eigenvalue
--the random-effects parameters and the residual variance...
--fixed-effect model matrix is rank deficient so dropping...
***混合効果モデル：独立変数が一つかつ離散変数 [#o1ad649c]
-混合効果（一元配置）分散分析
-参考：[[混合モデルを使って反復測定分散分析をする:http://...
--[[混合効果単回帰モデル:http://speechresearch.fiw-web.ne...
***混合効果モデル：独立変数が複数 [#f7c9443e]
-独立変数の性質によって呼び名はいろいろありますが、まとめ...
--全ての独立変数が連続：混合効果（重）回帰分析
--全ての独立変数が離散：混合効果（多元配置）分散分析
--連続・離散の独立変数が混在：混合効果共分散分析
-参考：モデルの考え方は[[混合効果モデル（変量効果モデル、...
-分析に使用するデータテーブル dataSumm.SS を作ります。
>data$dMFCC.SquareSum <- data$dMFCC1^2 + data$dMFCC2^2 + ...
>data = transform(data, seg2 = ifelse( ( (seg == "i") | (...
>data$seg2 <- factor(data$seg2)
>dataSumm.SS <- ddply(data, c('speakerID', 'seg', 'seg2',...
>dataSumm.SS$seg.length.log = log(dataSumm.SS$seg.length)
>dataSumm.SS$dMFCC.SquareSum.log = log(dataSumm.SS$dMFCC....
-お題１：''話者の個人差''の影響を除外した上で、''各音素の...
--ここで、''目的変数（従属変数）''は dMFCC.SquareSum(.log...
---seg と register は離散変数なので、混合効果（二元配置）...
--lmer関数を使ってモデルを構築します。
>ssANOVA1.lmer <- lmer(dMFCC.SquareSum.log ~ seg * regist...
---''「*」で繋げられた独立変数は交互作用を考慮''します。...
---ランダム変数は、切片・各独立変数の傾き・各独立変数の交...
---上記のモデルを実行すると、fixed-effect model matrix is...
>ssANOVA1_1.lmer <- lmer(dMFCC.SquareSum.log ~ seg * regi...
---それでも警告は消えません。続いてランダム変数の傾きを全...
>ssANOVA1_2.lmer <- lmer(dMFCC.SquareSum.log ~ seg * regi...
---それでも警告は消えません。仕方がないので、独立変数同士...
>ssANOVA1_3.lmer <- lmer(dMFCC.SquareSum.log ~ seg + regi...
--結果は以下で確認できます。詳細は省略します。
>anova(ssANOVA1_3.lmer)
>summary(ssANOVA1_3.lmer)
>step(ssANOVA1_3.lmer)
-お題２：''話者の個人差''の影響を除外した上で、''各音素の...
--''目的変数（従属変数）''は dMFCC.SquareSum(.log) で、''...
--lmer関数を使ってモデルを構築します。
>ssANOVA2.lmer <- lmer(dMFCC.SquareSum.log ~ seg2 * regis...
---今度は警告は出ませんでした。
--結果は以下で確認できます。
>anova(ssANOVA2.lmer)
--
Analysis of Variance Table of type 3 with Satterthwaite
approximation for degrees of freedom
Sum Sq Mean Sq NumDF DenDF F.value Pr...
seg2 110.023 55.011 2 5.8759 199.646 3.966...
register 3.356 3.356 1 8.1706 12.178 0.00...
seg2:register 4.192 2.096 2 7.5733 7.608 0.01...
---独立変数間の有意な交互作用と、各独立変数の有意な影響が...
--単純主効果の確認（≒多重比較）は下記のようにします。
>summary(ssANOVA2.lmer)
--
Linear mixed model fit by REML
t-tests use Satterthwaite approximations to degrees of ...
Formula: dMFCC.SquareSum.log ~ seg2 * register + (1 + se...
Data: dataSumm.SS
REML criterion at convergence: 11133.4
Scaled residuals:
Min 1Q Median 3Q Max
-5.2300 -0.6120 0.1008 0.6875 3.1298
Random effects:
Groups Name Variance Std.De...
speakerID (Intercept) 0.0091745 0.0957...
seg2sonorant 0.0070485 0.0839...
seg2vowel 0.0003944 0.0198...
registerRe-Read 0.0033398 0.0577...
seg2sonorant:registerRe-Read 0.0021252 0.0461...
seg2vowel:registerRe-Read 0.0034361 0.0586...
Residual 0.2755449 0.5249...
Number of obs: 7153, groups: speakerID, 6
Fixed effects:
Estimate Std. Error d...
(Intercept) 3.53758 0.04069 4.9380...
seg2sonorant -0.22829 0.04481 4.3800...
seg2vowel -0.22922 0.02024 8.0650...
registerRe-Read 0.15976 0.02883 5.2340...
seg2sonorant:registerRe-Read -0.17606 0.04774 11.0610...
seg2vowel:registerRe-Read -0.09066 0.03639 6.3130...
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1...
Correlation of Fixed Effects:
(Intr) sg2snr sg2vwl rgsR-R sg2s:R-R
seg2sonornt -0.679
seg2vowel -0.291 0.222
registrR-Rd -0.804 0.717 0.295
sg2snrn:R-R 0.302 -0.651 -0.092 -0.477
sg2vwl:rR-R 0.486 -0.549 -0.567 -0.741 0.353
--もしくは、lmerTestパッケージのstep関数を使います。こち...
>step(ssANOVA2.lmer)
-お題３：''話者の個人差''及び''音素の種類''の影響を除外し...
--''目的変数（従属変数）''は dMFCC.SquareSum(.log) で、''...
---seg.length.log は連続変数、register は離散変数なので、...
--lmer関数を使ってモデルを構築します。
>ssANCOVA.lmer <- lmer(dMFCC.SquareSum.log ~ seg.length.l...
--結果は以下で確認できます。
>anova(ssANCOVA.lmer)
--
Analysis of Variance Table of type 3 with Satterthwaite
approximation for degrees of freedom
Sum Sq Mean Sq NumDF DenDF F.va...
seg.length.log 3.2454 3.2454 1 38.378 15.8...
register 0.7014 0.7014 1 18.591 3.4...
seg.length.log:register 0.0262 0.0262 1 25.829 0.1...
---seg.length.logの有意な影響が見られます。レジスタの影響...
--単純主効果は以下で確認できます。
>summary(ssANCOVA.lmer)
--
(略)
Fixed effects:
Estimate Std. Error ...
(Intercept) 3.72788 0.09773 32.01...
seg.length.log -0.15572 0.04095 36.29...
registerRe-Read 0.13573 0.07338 18.59...
seg.length.log:registerRe-Read -0.01233 0.03447 25.83...
(略)
***ロジスティック回帰分析 [#p9f2d132]
***ロジスティック混合効果モデル [#k9d83231]
ページ名:
既存のページ名で編集する
