練習用のデータでできそうな分析としてはχ二乗があった。

本当は残差分析という奴もあるらしいがまあ練習なのでちょっとほっとこう。

χ二乗検定は独立性の検定。

カテゴリカルデータ同士になんらかの関係があるかどうかを検定するものなので、例えば無作為に100人の外人さんを抽出したときに、人種と性別の間に一定の特徴があるかどうかみたいな検定を行うものだ。

よくある分析としては、エクセルなどの表を作成したそれぞれのセルがお互いに関係しているかどうかを検定するという考え方。

ここでいうと、ちょっと分かりにくいが性別によって人種に一定の傾向があるかどうかの検定ということになる。

まずは、Ｒで表を作成する。
これは読んで地のごとしtable というコマンドを使うみたい。

いつものように、カテゴリカルデータの定義づけも含めてコマンドを入力。

s<- read.csv(“yome.csv”,header=T,sep=’,’)
s$sex<-factor(s$sex,levels=0:1,labels=c(“male”,”female”)) s$race<-factor(s$race,levels=0:2,labels=c(“white”,”black”,’yellow’))

• table(s$sex,s$race)

これでどうだ??

イエス!!次にこのテーブルをデータフレームに設定して(ここれはvと名づけてみた)、chisq.testを使う。

なんかできたぽいが、、、よくわからんメッセージ。

警告メッセージ：
In chisq.test(v) : カイ自乗近似は不正確かもしれません

グーグル先生!!

http://www.aichi-gakuin.ac.jp/~chino/psycstat/chapter7/sec7.html

どうやらtable の中の数が少なすぎると計算そのものはしてくれるけど、出てきた数字には信頼性がないという事みたい。

ふむふむそういう場合にはどうしたらいいのかしら??

どうやらこういう場合にはフィッシャーの正確検定というのをやるらしいね。これはオッズ比というヤツともカラムみたいだけど、今んとこは無視!!

フィッシャーの正確検定をやってみる。
関数はfisher.test

できた。今日は調子がよいらしい。

ここの結果としては、人種と性別は独立しているということ。
なんとなく、データに偏りがないという認識でもいいのかな。

実際に使うときにまた考えてみることにしよっと。

今日は分散分析を行ってみようと思ったのだ。
僕がいつも分散分析はWeb Anova(http://www.hju.ac.jp/~kiriki/anova4/)を利用しているが、もちろんＲでも出来るだろう。

改めて分散分析について考えてみると、自分がよくわかってなかったことに気づく。 っということでメモ。

今回の練習用データで考えてみると、荒っぽく考えてみと分散分析は、要因の多いｔ検定。

男女によって体重の平均値を比べるときの検定はｔ検定を行うことができる。2つの平均値を比べるのでこれでよい。

ただ、黒人と白人と黄色人種によって体重が異なるかどうかを検定するためには、分散分析を使う。
えっ??だったらｔ検定を2回すればいいんじゃないの?ということになるのだが、どうやらダメらしい(第二種の過誤というのが起きるらしい)。

ということで、今回のデータで分散分析を行ってみる。

こんな感じ。関数は一元配置はoneway.test

ほほう。 anovaで見るような表ができた。
この結果だと、p-value を見るとどうやら人種と体重の相関はないみたいだ。

これは一元配置の分散分析。水準数は白、黒、黄色の3水準であるといえる。

余談だがもし、2水準で分散分析を行った場合、ｔ検定と同じ結果になるらしい。ちょっとやってみよう。

まずは、t検定。男女によって体重に差があるのか?

続いて、一元配置の分散分析を2水準で。

OK!!なんか勉強になるなあ。改めて基礎から勉強している感じ。
あれ??そういえば、下位検定はどうやって出すんだろう?
まあ、とりあえずおいておいてここまでは一元配置の分散分析。

分散分析は一元配置ではなくこの配置の数を増やしていって交互作用ということも考えることができるらしい。

ここの例で言うと、
人種ごとに体重に差があるか?(人種の主効果)
性別ごとに体重に差があるか?(体重の主効果)
人種×性別ごとに体重に差があるか?(交互作用)

ということらしい。理論的にはいくらでも増やせるみたいなんだけど、解釈が難しくなるので2元配置くらいまでが普通ということ。

ではこの2元配置の分散分析をしてみる。関数はaov

なんのことやらよーわからんので、もうちょっと勉強してみよう。summary()を追加してみる。

コマンドは、summary(aov(s$weight~s$race+s$sex))

おお!なんかそれっぽいのが出てきたぞ!!

これあってんのかな???

性別の主効果は有意、人種の主効果はなし。
交互作用は何処?

試しにWeb Anovaで計算してみると数字がちょっと違う。。。改めてやってみると入力が結構大変だったなあ。
そうか、当分散の検定か。Ｒはウェルチがデフォルトなんだ。

ううーーん。酔っ払ったので今日はここまで。
というか本を忘れたので、ちょっと尻切れトンボですな。

今日の課題(というか本を見つけて追記すること)
1) tukeyの下位検定について
2) 交互作用はどうやって出力する?
3) 当分散を仮定しない分散分析を行ってその結果とWeb anovaの結果を照会する。
4) 当分散性の検定ってどうやるんだっけ?

ちょっとしたメモ。
1) Fix関数を使用したデータの変更と書き出し。
2) aggregateを使って一気に記述統計をまとめる(新しい練習用シートで)
3) χ二乗検定

前回、性別と人種という二つのカテゴリカルデータを定義づけたので、ちょっともう一度図にしてみる。

s<- read.csv(“yome.csv”,header=T,sep=’,’)
s$sex<-factor(s$sex,levels=0:1,labels=c(“male”,”female”))
s$race<-factor(s$race,levels=0:2,labels=c(“white”,”black”,”yellow”))
plot(s$race,s$weight)

こんな感じで散布図にしてみた。
あらかじめ、性別は、数字ではなくカテゴリカルデータとして扱っているのでplotを使ってもいい感じにしてくれるはず。

おう!なんかいい感じの図が出来てきた。
この図の形成はこれから考えるとしてと。

どうやって保存するんだろうといろいろいじってみると、

R→ファイル→別名で保存。でできるみたい。

でも。。。
サイズとかいろいろできるみたいだけどまあ今のところはいっか。

そういえば、度数分布表を作っていなかった。
練習用データの中のcutom(習慣)とInhe(遺伝)は5件法のデータなので、ヒストグラムを作る。

これは結構簡単だった。

hist(s$custom,right=FALSE)

これでよし。

こういうのをひっくるめてfurequency 関数というヤツもあるみたいなのだが、もっと根本的なことが分かってからにしよっと。

たちまち今までの復習。
Rエディタにまとめておいとく。#を使ってメモ。

とりあえず少しずつだな。

次は、分散分析にチャレンジ!!

カテゴリーデータの扱いについてべんきょーしてみる。

ということで、今使っているデータはカテゴリーデータが性別しかないので、ちょっと追加してみた。

race は人種で。0-white,1-black,2-yellowということにしてみよう。保存してみた。

Rを起動してsを読み込んでみる。

あら、そのままデータは反映されないのね。
本当ならsの次のファイルで読み込んじゃえばいいと思うのだけど。
ちょっといろいろ作りすぎてわけわからんくなってるので一度がっつり消しちゃいましょう。

編集のコンソールの消去とワークスペースの消去をしていったん奇麗にします。

これでsを入れてもでてこない。

ということで最初から読み込み直し。

s<-read.csv(“yome.csv,header=t,sep=’,’)
s
よし、だいぶ覚えてきたなあ。

さて、このデータのうち、sexとrace は性別(男、女)と人種(白人、黒人、黄色人種)だから，コイツをうまい事しないと。

ふむふむ。カテゴリーデータの定義付けというヤツが必要。

この定義付けがうまくいけば、分散分析やフループごとの散布図などが出来るらしい。

試行錯誤の結果作ったコマンドがこれ。
性別について定義付けを行う。
> s$sex<-factor(s$sex,levels=0,labels=”male”)
sのデータの変数のsex の0というのは、maleですよ。

というようなコマンド。これでsを見てみる。

よし、男性はmaleになった。続いて女子もやってみる。
といいたいところだけど、めんどうなので一気にやりたいよね。

っとかなり試行錯誤した結果がこちら

s$sex<-factor(s$sex,levels=0:1,labels=c(“male”,”female”))

ミソは、0:1 の所。これは0と1ではなく、0から1という感じで捉えるべきみたい。

続いてrace(人種)
s$race<-factor(s$race,levels=0:2,labels=c(“white”,”black”,”yellow”))

sのデータのraceという変数に右辺の式を代入しますよ。
sのデータのraceという変数は、0～2まであって、そのラベルは
white black yellow ですよ。という感じかな。

よーしいってみよう!

きたーーーーー!!!
ちなみに一口に試行錯誤とは書いてありますが、だいたい2時間くらいはかかってますけどね。

今日のまとめのRエディタ。

(メモ)
この定義付けは最初にエディタで一気にしたほうがいいと思う。なにがなんやらわからんくなる。

次回はこれを使って、散布図を作成してみよっと。

今までのデータを保存して終了。
そうか、でもyome.Rを持ってれば後からそれをコピペすればいいんだ。

思いつきでt検定をしてみた。

t.test(s$weight~s$sex)

これは、sというデータの中で体重が性別ごとに異なるかt検定をやりなさいというコマンド(だと思う)。

でたでた。なんかいろいろ出てますがな。蕁麻疹が。

Welch Two Sample t-test

多分ウェルチを使った2つの標本のt検定を行ったという感じ。

data: s$weight by s$sex

データはsの体重を性別で。

t = 3.7311, df = 18.305, p-value = 0.001493

自由度(df) P値は(p-value)

alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:

くーーー蕁麻疹。

yahoo翻訳したれ。

対立仮説: 手段の真の違いは0 95パーセントの信頼区間と等しくはありません:

???

8.21321 29.32525
sample estimates:

mean in group 0 mean in group 1
69.69231 50.92308

ここら辺は気にせずに行こう。
なんか、5%水準で有意ということでいいっぽい気がするが。。。

ちょっと本(青木先生のヤツ)をごそごそ。
いろいろ載ってんなあ。まだまだぜんぜん理解できん。

どうやらRではウェルチがデフォルトらしい。
ウェルチは、等分散を想定しない場合の検定。

もし、等分散を想定して検定するならば

var.equal=TRUEを追加するらしい。

ということでやってみた。

OK!! できましてん。
じっくりもぞもぞですな。

今回の疑問
1 5%水準じゃない場合(1%とか10%)を出力できないか?
（追加編集）
daihikoさんからのコメントで解決。
対立仮説: 平均値の真の差異は、0と等しくない
95%信頼区間

という意味。「差異は、0と等しくない」ってことは、差があるってことだよね。
それに対して、帰無仮説は「差異は、0と等しい」。

有意水準はp-valueのところを見ればいい。ここが.05未満だったら5%水準、
.01未満だったら1%水準、というように見る。

そもそもの和訳が間違っていたということらしい。
今回のデータでは，ｐ値が0.001493なので1%水準でも有意ということ。

そっかもっともっと統計のことを本質的に勉強しないと。
2 等分散かどうかの検定(ExcelでいうF検定)はどうするの

ちなみに昨日Rの勉強会だったので、このブログをN先生に見せたのだが、なんと身長はtallではなく、heightが正しいらしい。

めちゃくちゃ恥ずかしいが、今更治すのも面倒なのでこのまま押し切る事にしよう。

僕のししょー(生理、実験系)は、
「データをGetしたらまず散布図を作りなさい。そしてその散布図を眺めなさい」

と僕に教えてくれた。

ということでまずは散布図を作成する。
一番簡単そうなtall(正しくはheightね)とweightでチャレンジ。

plot(tall,weight)
これで散布図は出てくるはず。

あれ?オブジェクトがありませんときたもんだ。

「!!」

そうだattach(s)がない。
このattach(s)は、コマンドの最初に入れた時に以降のコマンドは全部sのデータを利用するというような意味で、エディタで一気にコマンドを入れてしまうような場合に便利(と理解している)。

今回は、コンソールでいろいろやるのでいちいちコマンドにsのデータを使うというメッセージを入れてみる。

s$と入れるのが正しいらしい。
plot(s$tall,s$weight)

これでsのデータの中から、tallとweight の散布図を作れ というコマンドになる。

よし。できた!!美しい散布図だ。このファイルは、形成できるのかな??論文で使えるのかしら?保存は?まあまだまだわかんないことばかりだが、放置して次に行こう。

しかし、いちいちこれをすべての組み合わせでやるのが面倒くさいよね。

ということで総当たりの散布図

plot(s)

ししょー!!
たしかにわかりやすいですたい!!散布図万歳!!

ほほう。やはり身長と体重の散布図を見ると相関はありそうですね。ということで関数corを使ってみる。$を忘れないように。

cor(s$tall,s$weight)

さらに総当たり戦もやってみる。

cor(s)

よし。なんかいい感じ。
noはIDにしとけばよかったかな。

相関表を見てみるといろいろ気になる事が出てきた。

1．性別は、男性0,女性1にしてるけどこれはカテゴリカルデータだよね。数字データのままで大丈夫?
2．custom(習慣),Inhe(遺伝) は5段階の順序尺度だけどこれもこのままでいいのか?
3．相関の有意は出ないのかな?

しばらくは検討を重ねる。

今作っている妄想体重データで統計的な検定ができそうなこと。

1.男女の体重差のt検定
2.分散分析(ちょっと勉強しなおさなきゃな)
3.体重とその他の変数との相関検定
4.体重を従属変数とした重回帰分析(カテゴリカルデータ含む)

他になにかあるかな。。。

前回の続き。
ちょっと行き詰まっていたが、なんとか回復。
なんでうまく行かないかが分からないままなのでここでは触れずにもう一度Rを起動した。

今までのおさらいという意味を含めて一気にやってみる。

まずは、Rエディタを新規で開いてこんなコマンドを書いてみる。

これは、Rコンソールの中にコマンドが読み込まれて待機している状態。

そこでweight_sdを呼び出してみよう。

でた。体重の標準偏差。

こういう計算をするにはExcelを使った方がラクかもしれないが、実際データをいじくろうとすると、表の状態でさがすより思いついた計算をどんどん呼び出せるので便利かもしれない。

ここで復習。
最近知ったのだが、コマンドの前に#をつけるとそのコマンドは実行されないらしい。
ということでさきほどのエディタに解説を加えて実行。

よし!!できた。

聞いた話によると、プログラムを作成するにあたってアンダーバーは滅多に使われないそうです。その他の記号はコマンドの一部と認識されてしまうことがあるのであまり使わないほうがよさそう。これはしっかり癖にしちゃおっと。

あとこれからエディタを作るときは、基本的に#をつけて解説を入れるようにしよう。慣れるまでは日本語で。

というか、いつも簡単に出来ているように見えるかもしれないが、何度も何度もミスを繰り返してなんとか出来ているというのが実際なのだ。ここまではなんとか出来たので、次からは本格的に統計をやってみよっと。

今日は、ビールを飲まなかったので夜分遅くではあるのだが少しRをいじってみよう。

Rを起動し、sを入力てみると、しっかりと前回のデータが読み込まれている。

このまま関数を入れてみる。

max(age)

あれ??年齢の一番大きいヤツを呼び出したいんだけどな。

このオブジェクトがありませんは、うまくいってませんよというメッセージ。age という変数が分かってないということみたい。

いろいろ調べるとコイツがミソだ。

attach(s)

これは、以下のコマンドは(s)のデータの中にありますよという宣言のようなもの。

ということで
>attach(s)

を打ち込むと
コンソールは待機状態(sを入れて次のコマンド待ってるぜ!という感じかな)

続いて
>max(age)

きたー!!

ということで年齢の最大値はできるようになった。

勢いに乗ってmini(最小値)、mean(平均値)、sd(標準偏差)も入れてみよう。

よし。ここまでは完璧だ!!

ちなみに僕が使っている本。
勉強会を始める前に、ちょっとやってみようと思って挫折した本がこちら。

最近みてみると言葉は分かりやすいけど、うまく行かなくなるととたんに難しくなるという感じの本でした。統計の解説が分かりやすいのではじめての人にはおすすめ。
お値段2700円なーーり。

続いてこちら。これは勉強会がはじまってから、N先生のおすすめで購入。
プログラムがたくさん書いてあって、自分一人では絶対読めない漢字だけど、最近ほんの少し分かってきた。著者の青木先生はネット上でもたくさん情報を公開されているのでかなりおすすめの本だと思う。

しかし、こういう本を書ける人の頭のなかっていったいどうなってるんだろう。

すげーな。