サンプルサイズが条件ごとに異なる一要因分散分析

今回は、一要因分散分析について考えます。典型的な教科書では、分散分析の条件ごとのサンプルサイズが揃っている場合を扱います。もしも条件ごとのサンプルサイズが大きく異なる場合はどんな問題に気をつければよいかを見てみましょう。

ちなみに似たような疑問をt検定について取り上げたことがありますので、よろしければそちらも御覧ください。

webbeginner.hatenablog.com

そもそも分散分析とは？

分散分析という名前が示す通り、２つの分布の分散を比によって比較します。主な用途は実験や調査で設定されたグループや条件の比較です。グループ間、または条件間の平均値差をグループ間のデータ変動（ある種の分散）とみなして、グループ間変動の大きさが、グループ内のデータ変動に比べて何倍大きいのかを計算します。

各グループのサンプルサイズが大きく異なる場合

統計学のテキストなどでは、各条件のサンプルサイズが揃っている場合が扱われることがほとんどです。そこで、今回はあえて、サンプルサイズの異なる３つの条件を比較してみたいと思います。

以下のデータは、R言語のrnorm関数を使って生成しました。条件Aと条件Bは５０個のサンプルがありますが、条件Cでは１０個のサンプルにしています。

上記の下半分に示されている数字（-0.283, 0.936, 0.967）が各条件の平均値です。ここで条件Aと条件 Bの差よりも条件Aと条件Cの差の方が大きいことを覚えておきましょう。

さて、このデータについてR言語のlm関数を使って一要因分散分析をした結果が以下の通りです。

R言語の出力のうち、最終行が分散分析の結果を示しています。F(2,107) = 24.49, p < .001.と出力から、論文やテキストで良く見かける表記に書き直しました。さて、この結果が示すように、分散分析の結果は、0.1％の有意水準で有意です。３つの条件のいずれかに偶然とは呼べないような差があると考えられます。実は、lm関数の出力を見ると、A条件とB条件の差、A条件とC条件の差が既に示されています。出力上半分のcoefficientsの内、groupBとgroupCという行がそれです。この出力結果は、groupAを基準とした時のgroupBとgroupCの差が示されています。一番左のEstimateの部分です。みてわかるように条件Aとの差異は条件Bより条件Cで大きくなっています。ただし、この結果をt値という統計量で見ると結果が逆転します。条件Aと条件Bの差に対応するt値は6.673、条件Aと条件Cの差に対応するt値は3.948です。t値を計算する差異には、各条件の平均の差が条件内のばらつきを基準に計算されます。例えば、サンプルサイズが小さいと平均値についての誤差は相対的に大きくなるので、t値は小さくなります。つまり、サンプルサイズが小さい条件を含む比較では、有意差が得られにくくなります。今回の数値例では、どちらの条件比較も有意差がありますが、場合によっては平均値の値の上では同じようなペアであっても有意差が出たり出なかったりすることも考えられます。これを考えるとやはり一般的なテキストに書かれているように「各条件のサンプルサイズはなるべくそろえるように」という言葉に従っておいた方が良さそうです。

R言語の一要因分散分析の関数はどんな計算をしている？

さて、条件間でサンプルサイズが違う場合の分散分析をする上で、１つ気をつけた方がよいことがあります。それは、すべてのデータの平均の計算方法です。今回のサンプルでは、３つの条件の平均値がそれぞれ-0.283, 0.936, 0.967でしたので、その３つの数字の平均を単純に計算すると0.54になります。しかしこれはデータ全体の平均値としては不適切です。各条件のサンプルサイズを考慮した平均値は0.3847235(2017/1/16 修正しました)となります。分散分析をするときに、各条件の平均値を基準としたデータのばらつきは、サンプルサイズの計算に影響を及ぼさないのですが、条件間のデータ変動を計算する差異には、全体の平均がどこになるかで結果が異なるので注意が必要です*1。

*1:R言語やMS Excelの分散分析では特に指定しなくても重み付けされた平均値が計算されるので安心してください。

R言語での変数の型確認の重要性(1)

今回は、R言語で統計解析を学び始めた人にとって分かりにくい点について書きます。

Rの変数オブジェクトは型宣言しないが、型がある

Rは、統計パッケージアプリケーションという捉え方と、データ処理のためのプログラム言語という捉え方があります。後者に関してRを特徴づける方法がいくつもあります。その１つが変数の定義方法です。C言語やpythonでは、変数を作成するときに、その変数の性質を決めるために型の宣言を行います。例えばxという変数が整数を扱う変数であればintという型を与え、文字列を扱う変数ならstringという型を与えます。一方、Rでは、このような型宣言なしで変数を作成できます。ただし、Rでも変数には以下のように幾つかのタイプがあります。

• NULL型：NULLつまり空のオブジェクトを扱うための型です。
• 命題：真(TRUE)または偽(FALSE)だけを扱うための型です。以下で述べるような型変換のプロセスを経ると1または0を入力として受け付けることができます。
• 文字列：文字列を扱うための型です。
• 数値：数値を扱うための型です。

型は変換が可能

オブジェクトの方は変換が可能なのですが、その方向性は次のように決まっています。

NULL型→命題型→数値型→文字列

つまり、命題型を数値型に変換できるが、文字列型は数値型には変換できません。

型確認の方法

型を確認する方法がいくつかあります。一つはclass関数です。引数として型を知りたいオブジェクトを指定すると、そのオブジェクトの型が出力されます。もう一つはstr関数です。この関数の引数としてデータフレームオブジェクトを設定すると、そのデータフレームに含まれる変数の型を一つずつ示してくれます。

Rのオブジェクトは宣言時に型を指定しない

R言語では、オブジェクトの型が、代入される値によって決まります。なので、自分がオブジェクトに何を代入するのか（したのか）を把握しないと後々、エラーにつながってきます。この投稿の続きでは、そんなエラーが現れやすい例として回帰分析を取り上げる予定です。

音の印象を尋ねるウエブページの構築

本日は、音声を聞いて、その印象を尋ねるウエブページをHTML5の標準技術で（プラグインや別のサービスなし）作ることを試みます。これができると音声を聴取する調査がぐっと簡単になりますよね。

## 概要

• 音声を提示する
• 反応を取得する

音声を提示する

音声の提示はaudioタグを使います。audioタグには再生するファイルの指定や再生方法に関して以下のようなパラメータがあります。

• controls: 再生ボタン、音量、再生位置などの表示（ trueを明示しなくて良い）
• src: 再生するファイルの指定
• loop: trueまたはfalseでループ再生するかどうかを指定
• preload: 音声をあらかじめ読み込むかどうかを指定(none, metadata, auto)
• autoplay: 音声の読み込みが完了すると同時に再生を始めるかどうかを指定

実際にはこんな感じで書きます。
サンプルのコントローラは、音楽は魔王魂さんの作成された楽曲ファイルが再生されるように設定されています。
全曲無料・フリー音楽素材/魔王魂

<div>
<audio controls src = "ファイル名" ></audio>
</div>

反応を取得する

反応の取得はいろいろな方法があります。例えば、多肢選択による反応取得だったら、ラジオボタンがいいかもしれません。

<form>
<input name = "question" type = "radio" value = "yes">好き
<input name = "question" type = "radio" value = "no">すきではない

</form>

好き
すきではない

以上のようなhtmlフォームをjavascriptでつなぎ合わせて、反応を配列に記録していけば良いでしょう。

rmarkdownファイルのテキストで改行する方法

今回は、非常に細かい話。テキスト改行のほんのちょっとしたコツです。

改行しているのに反映されない

「テキストを改行したければ、returnキーを押せばいいのでは？」と思う方もいるとは思いますが、RStudio内のrmarkdownファイルでは、ファイル内で改行していても、実際に出力されるpdf,htmlファイルでそれが反映されません。確実に改行するためには、改行コード前に空白スペースを２つ入れる必要があります。以下の例をみてください。

最初の一行。この後で改行を試みます。
改行コードの後の文。しかし実際には改行されていません。

上記の文をrmarkdownファイル内で書くと、改行されません。そこで改行の前の文末に空白スペースを２ついれると改行されます。

最初の一行。この後で改行を試みます。  
改行コードの後の文。今度は空白スペースを２つ挿入してあるので改行されます。

ショートカットキーとしてのタブキーを活用する

このやり方がわかって、実際に試してみると、日本語の使用者にとって不便な解決方法であることがわかります。なぜなら、空白スペースを入れるためだけに、macなら「かな」から「英数」に切り替え、スペースを入力後すぐに「英数」から「かな」に戻す操作が必要だからです*1。そこで、タブキーを「空白スペース２つ入力」のショートカットキーとして利用しましょう。実は、RStudioではタブキーを押すと「空白スペース２つ挿入」するがデフォルトの設定になっています。これはRStudioのtoolsメニューにあるGlobalOptionを選択し、codeメニューのediting内に設定があります。Tab widthの数字で変更可能です。

パラグラフを変えるには改行２回でOK

若干不思議なのですが、空白スペースをはさまなくても２回連続で改行すると、１行空きで新たな行を始めることができます。

まとめ

rmarkdownでの改行は？

• 単純な改行なら文末に空白スペース２個(タブキーで入力)＋改行コード
• 段落を変えるなら２回連続の改行コード

*1:windowsなら一旦「全角」から「半角」に切り替え、スペース入力後に再び「全角」に戻すことになります。

三平方の定理で理解する回帰分析における分散の分解

今回は、回帰分析の細かいモデルの話を中学校で習う「三平方の定理」を利用して考えてみます。

最小2乗法による推定

回帰分析では、最小２乗法という計算方法で、回帰直線をひくための２つの推定値（切片と傾き）を決めます。何を最小にするかというと実際に観測したデータとモデルによる推定値の間に生まれる誤差（残差と呼びます）を最小にしたいわけです。２乗というのは以下のような計算方法で残差を計算することによります。

• 予測したモデルにおける推定値と観測した値を引き算する。
• 引き算の答えを２乗する。この段階で、推定値が大きくなる場合でも、小さくなる場合でも、２乗した後の値の符号がプラスになります。
• ２乗した値を全て足しあわせます。（これを平方和と呼びます）

この計算過程は分散の計算方法によく似ています。分散では、データの各値の平均からのズレを２乗して全て足し合わせます。さらにこれをデータ数（または自由度）で割ると分散になります。この分散と残差の計算方法の類似を考慮に入れれば、データのばらつきが大きいときは残差（誤差）が大きいといえることが理解できるのではないでしょうか。

平方和の分解と分散説明率

さて、ここで回帰分析における平方和は次のように分解できます。

従属変数の平方和　＝　予測された値の平方和　＋　残差の平方和

そして分散説明率は以下のように評価されます。

予測された値の平方和　/ 従属変数の平方和

つまり、元々ある従属変数のばらつきの何割を予測によってカバーできたかが示されています*1。ここまでの説明は、統計学のテキストで共通に出てくるのですが、わかりそうで分かりにくい表現になっている気がします。

三平方の定理の定理で理解してみよう

三平方の定理は、三角形の斜辺の二乗はその他の２辺のそれぞれの２乗和に等しいというやつです。これが上記の平方和の分解に当てはめられるのですが、それを理解するためには、まず平方和の表現を幾何学的にしてみたいと思います。

このアイデア自体は、自分で思いついたのですが、いろいろ調べてみたところ、類似の説明はすでに以下のサイトや次の書籍で行われています *2。

朝野 煕彦 (2000). 入門 多変量解析の実際 第2版 (KS理工学専門書)

http://d.hatena.ne.jp/rikunora/20131202/p1

従属変数の平方和

まず、従属変数の各サンプルの値から従属変数全体の平均を引きます。以下の図の真ん中から丸までに引かれた線分がそれです。

これを２乗するのですが、グラフ内では先ほどの線を一辺とする正方形の面積になります。

これを足し合わせたものが従属変数の平方和です。

予測された値の平方和

従属変数の平方和と同じような計算方法で、予測された値の平方和を求められます。ただし、今度は、予測された値から従属変数の平均値を引いた線を一辺とする正方形を作り、それを足し合わせます。

残差の平方和

残差の平方和は、各サンプルの値から、それに対応する予測値を引きます。この線を一辺とする正方形を作り、その面積を足し合わせます。

３種類の平方和を三平方の定理に当てはめる

ここまでで３種類の平方和を計算しました。これを三平方の定理に当てはめると斜辺は従属変数の平方和、残りの二辺は、予測された値の平方和と残差の平方和となります。図にすると以下のようになります。

そしてこの図は、以下のような内容の理解に役立ちます。

予測値と残差の相関は０

予測値の平方和と残差の平方和は、直角になるような位置に置かれています。これは、それぞれの値が互いに影響を及ぼさない（つまり相関が0）ことを覚えるのに役立ちます。

従属変数と予測値の指す方向の類似性がモデルの説明率理解に役立つ

予測値の平方和と従属変数の平方和の位置関係は、それに対応する２辺の角度で表せます。そして、この角度は、予測値の平方和が大きくなればなるほど小さくなります。言い換えれば、２辺が指し示す方向がよく似ているということです。これは、予測値の平方和が大きいほど、モデルとしてうまく実データを説明できていることを覚えるのに役立つでしょう。

逆に、モデルとして実データがうまく説明できていないような場合は、予測値の平方和と従属変数の平方和の関係は直角に近くなります。

まとめ

今回は、回帰分析の各変数の散らばりを三平方の定理で考えてみました。まとめとして、回帰分析において、従属変数のデータのばらつきは、予測値のデータのばらつきと残差のばらつきに分解できることがわかりました。できるだけ予測値のデータのばらつきを大きくすることを目指して、いろいろな現象を説明、予測してみましょう。

*1:なので分散説明率は0から1の値をとります

*2:他の人も使っている説明方法という意味で安心して主張できます

相関と回帰は何が違うか？

今回は、相関分析と回帰分析は何が違うのかを考えます。

相関と回帰の共通点、類似点

相関分析と回帰分析はどちらも、２つの連続量変数（数量型データ）の関係を調べるために行います。ここでの「連続量変数」とは、単に数字で表したデータという意味ではなく、条件が決められています。まず、この変数の尺度水準は間隔尺度です。例えば、１と２の間の１と２と３の間の１は等距離であるということです。また、この変数は確率分布として正規分布を仮定します。
以上のような共通点に加えて、二つの分析によって求めたいことがとても似ています。おおざっぱにいってしまえば、どちらも一方の変数の値の取り方の大小が、他方の値の取り方の大小と結びついているかを調べる分析です。

因果関係の有無が違う

次に、相関分析と回帰分析の違いを考えましょう。決定的に違うのは、因果関係の想定の有無です。相関分析では、変数間の因果関係を想定しません。どちらが原因、または結果でもよくて、相互に影響し合っているのでも構いません。

一方、回帰分析では、想定としての因果を明確にします*1。原因となる側の変数は独立変数、予測変数などと呼ばれます。これに対して結果となる側の変数は、従属変数、応答変数などと呼ばれます。回帰分析は、その予測を式のモデルで表すことができます。

Y = a + bX + e

ここでYは従属変数（つまり予測される側）、Xは独立変数（つまり予測する側）です。aは切片、bは係数または傾き、eは残差と呼ばれています。eを除いた部分は、中学校の一次関数で出てきています。一次関数、つまり、直線を引くための式です。ここで、先ほどの切片と傾きという言葉の意味がわかると思います。切片は横軸の値（独立変数）が0のときの縦軸の値（従属変数）であり、傾きは、直線がどのような角度になるのかを表しています。回帰分析の詳細については、別の投稿でまた書きます。

R言語を利用するときの注意

相関分析では、相関係数を計算します。相関係数は-1から1までの値をとる指標で、-1に近いほど、「一方の変数の値が大きくなるほど、他方の変数の値が小さくなる」ような関係を示しています。この関係のことを負の相関関係と呼びます。また、1に近いほど「一方の変数の値が大きくなるほど、他方の変数の値も大きくなる」関係を示します。こちらは正の相関関係と呼びます。-1と1の中間、つまり0になったときは、2つの変数間に大小関係の結びつきがないことを示しています。このような場合を無相関と呼んだり、２変数が独立の関係であると言ったりします。R言語では、相関係数を求めるためにcor関数を用います。例えばvarAとvarBの相関係数を求めるなら

cor(varA,varB)

と書きます。さらに、この相関係数が0であるという帰無仮説を立てて、相関係数の有意性を検定するにはcor.test関数をつかいます。

cor.test(varA,varB)

ここで引数のとり方を逆にしても同じ値が返されます。どちらが原因とも結果とも決めていないので、強弁する強さのみが示されていると考えられます。

さて、２変数の関係を調べる場合、相関係数のみから結論を導くと誤った結論になることがあります。例えば、変数の中に外れ値と呼べるような極端に大きい値や小さい値が含まれていると、相関係数の絶対値は大きくなってしまいます。なので、相関関係を調べるときには、必ず散布図も書きましょう。散布図を書くにはplot関数をつかいます。

plot(varA,varB)

この例の場合、最初の引数であるvarAが横軸、２番目の引数であるvarBが縦軸の値になります。相関係数の場合と同じで、どちらを先にしても良いのですが、グラフ上のどちらの軸になっているかを覚えておきましょう。これは後で回帰分析のときに大事なポイントになります。

実際に散布図を描いてみると、正の相関が強く見られる場合では、データポイントが右肩上がりになるように分布し、負の相関関係が強く見られる場合では、データポイントの分布が右肩下がりになります。無相関の場合は、全体にまんべんなく分布します。

次に、R言語を利用して回帰分析をするには、lm関数をつかいます。

res <- lm(varB ~ varA)
summary(res)

上記の例では、resというオブジェクトに回帰分析の結果を保存しています。これをprint関数の引数として表示すると切片と傾きの推定値しか出力されません。summary関数を使うと、独立変数の傾きの有意性やモデル全体の説明率が表示されます。
また、２つの変数のうちどちらを従属変数に設定するかで、切片や傾きの推定値は異なるので注意が必要です。

相関分析と同じように、回帰分析でも散布図が便利です。以下のようなコードで散布図に回帰直線を加えられます。

plot(varA,varB)
abline(lm(varB ~ varA))

ここで重要なのは、回帰分析ではどちらの変数を横軸、縦軸にとってもよかったのですが、回帰分析では、独立変数を横軸、従属変数を縦軸にとらなければならないことです。これを誤ると、回帰直線が散布図上のデータに全く沿わない形で描かれてしまいます。特に、plot関数では最初の引数が横軸になりますが、lm関数の式は従属変数から書き始めるので注意しましょう。

まとめ

以上をまとめると次のようになります。

• 2変数に関して、特に因果関係を想定せずに、共変動する関係にあるかどうかだけを知りたいなら相関分析（cor関数、cor.test関数）
• 2変数に関して、因果関係を想定し、一方が他方に対して持つ影響力の有無や、その予測を知りたい場合は回帰分析（lm関数）

これらを意識しながら、適切な分析方法を選ぶようにしましょう。

*1:想定した因果関係が実在するとは限らないので注意が必要です。むしろそのような因果関係があるかを検証する方法の１つとして回帰分析があります。

アクセス数の合計が200,000を超えました！ありがとうございます！

本日は、ご報告とお礼です。

本ブログのアクセス数の合計が200,000を超えました。前回のご報告でアクセス数100,000突破が今年２０１６年２月のことでしたので、１年かからずに100,000アクセスをいただいたことになります。

最近のアクセス傾向

これまでのアクセス傾向とそう変わりなく、単純に1日あたりのアクセス数がじわじわと伸びた結果、意外に早い200,000アクセスとなりました。前にもご報告しましたが、このブログは平日のアクセス数が休日のおよそ２倍あります。そして、今年の夏に気がついたのですが、夏休みになってもほとんどアクセス数が減りませんでした。どうやら学生よりも社会人の方にたくさんアクセスしていただいたようです。実務に十分な情報を提供できるよう今後も分かりやすい表現にこだわっていきたいと思います。