最近になってから、Rへの情熱が復活してきたので、とりあえずSPSSで分析した修士論文のデータをRを使って検算した。
っと一言で言ってみたけれども、実際には1月くらいは七転八倒していた。

twitter上でつぶやいているといろいろな人の協力を得られるもので、@kosugitti先生にはわざわざ時間を取ってもらったり、@boz_takemura 先生 にはコードを送ってもらったり。

ありがたやありがたや。。。

受けた恩は後進に返さなきゃ(ペイフォワード!!)ということで、検算の過程で悩んだり、七転八倒した所をできるだけ素人目線でメモしておこうと思います。たくさんあるので小出しにしていきます。

まずは、ちょっと昔のエントリにも書いたけど不正なマルチバイト文字について
http://d.hatena.ne.jp/jigawa91/20110531/1306849846

前回はshift-jisモードで読み込むという方法を取ったのだけど、どうやらその方法はあまりよろしくないみたいです。本来であればUTF-8という文字形式で保存したほうがいいらしい。

まあここは「みたい」とか「らしい」でしかないのですが。。。

ただ、僕はExcelでデータを打ち込んでしまっていているし、どこが不正なマルチバイト文字なのかよくわかんねえよ。。。ということで、方法を探ってみました。

おそらくとんでもなく迂遠な方法だと思うので、もっと簡単な方法を知っている人がいれば教えてもらいたいもんです。

まずは、これがExcelで打ち込んだデータ(ダミーデータにしています)。

ここから、まずファイル→名前をつけて保存→タブ区切りテキスト

という形式で保存して、とりあえずtxtファイルを作成。

そしてそのtxtファイルをエディタで開いて(僕の場合はJeditX)、上書き保存時にUTF-8で保存する。

ここから、txtのデータをタブ区切りで読み込むようにRに命令。
law <- read.table(“blog1.txt”,head=T,sep=”,”,na.strings=”NA”)

※ちなみに、na.strings=”NA” は欠損値をNAとするという意味です。
僕は最初欠損値を ANと間違えて書いていたおかげで文字データと数字データがわけわからん事になってましたww気をつけましょう。。

うん、間違いなく間違ってるな。。たぶんタブ区切りって言うのがうまく伝わってない。。。

まあ今までRでデータ読み込む時には、csvだったのでさっきのtxtをExcelでもう一度読み込んでcsv で保存し直してみよう。このやり方は、エディタによってはタブを,で置換してcsvにするとかちょっと面倒なことも必要でしたがまあなんとか。。。

どうだーー!!
law <- read.csv(file.choose(),head=T,sep=”,”,na.strings=”NA”)

できた!!ちなみにfile.choose()は、@kosugitti 先生に教えてもらった関数で便利です。

赤線の部分に注目。しっかりとデータフレームが読み込めたことが分かる。
lawというデータフレームの中に、373のオブジェクト(ここではn)と21のvariables(ここでは変数)が入っていると確認。

Rstudioを使えば(デフォルトと配置を変えてるけど)、もとのデータを簡単に見る事もできる。
クリックすれば自動でview()が発動するので、非常に簡単です。

今日はここまで。

今日は質問紙の分析に入ろうかと思います。

まずは、質問紙の構成の確認。

このダミーデータは先行研究においては3つの因子からなる15項目の質問紙。
m_1～m_5 までが第1因子、m_6～m_10が第2因子、m_11～m_15が第3因子となっている。
まずは先行研究と同様の因子構造を使用していいかどうかを検討するためにα係数を算出しよう。

コマンドは
alpha(qp[c(“m_1″,”m_2″,”m_3″,”m_4″,”m_5”)])

これで第1因子を構成する質問項目の信頼性の検討を行う。
psychパッケージの中にあるalpha関数を使うのが便利みたい。
これまた@kosugitti先生に感謝。

ちなみにpsychパッケージの使い方よく分からんのうなんて思ったらhelp(psych)をコンソールに入れるとhelpが出てきます。英語なので怖いけど、僕のように現在完了形から見失っているような英語力でも少しは理解できる単純な英語なのでビビらず見てみよう!!そこの検索窓でalphaを検索して見てみたのがこちら。

結果の見方が書いてあるね～!!

ちなみに今回のデータの結果がこちら。

まずピンクの所をみると、law_alphaとかあるけど、なぜか小数点第一位までしか書いてないなあ。law_alphaとstd.alphaが出てるけどこれって少数点第二位まで出ないのかしらん??

赤線の所を見ると、Reliability if an item in dropped と書いてある。
もちろん英語にはうんざりしてるが、これはその質問項目を除いた時のα係数の変化を示している(のだろう)。
ここの数値が、上にあるα係数よりも著しく大きくなるようであればちょっといかんのだ。だから今回のデータの信頼性はまあまあオッケーだなということが分かる。

紫で囲んだ所は、それぞれの質問項目の記述統計。平均も偏差も出してくれているので明らかな偏りがあればチェックでできる。

その下のとこは何なんだろう???
質問項目間の相関係数??

結果の見方はもう少し勉強しないと行けないがなんとか、信頼性の検討はできた。

このalpha係数は本来であれば因子分析をした後にやることが多いけど、先行研究がたっぷりある場合とかはこれだけでも結構説得力はあるねえ。
しかし、α係数の小数点第二位まで出ないのかしら?

蛇足ですが、@kosugittiに教えてもらったomegaについても備忘録的に記載。
参考ページはこちら。
http://d.hatena.ne.jp/kosugitti/20091224/1261643078

omega関数は、指定した因子構造と仮定してその負荷量から計算してくれるもの。

今回のデータだとこんな風に書く。
omega(qp,nfactors=3)

なんか変なんでました。there in no packegeと言うのがでていれGPArotationというパッケージがインストールされていないのが問題みたい。Rではこういう事はよくあるので、焦らず、tools→installs packegeでインストールして再試行。

でました。ここらへんはまだまだ勉強不足。あくまでやってみただけという感じですがなんか見た目が格好いいな。
因子分析等いろいろやった後でまたチャレンジしよう。
今日はここまで。

しばらくお休みをしていましたが、また復活しました。

備忘録的にメモ。
ちょっと前にattach()をあきらめていたのだけれども、データセットをする前にattach()を使うとなんだかよくわからないことになるので、すべてのデータセットとカテゴリカルデータや順序尺度を定義づけてしまってから使うと便利になる。

その後のデータフレームを使いする時も同様だと思うので注意。

あとshift-Jisで読み込むのもやめておく。
UTF-8で保存がいいみたい。

詳細は落ち着いたら画像付きでアップする予定。

先週の日曜日にHiRoshimaRというRの勉強会があってとても勉強になった。

今日はそこで@kosugitti 先生が紹介してくれた心理学者のためのパッケージを丁寧に七転八倒することにした。

まずは、パッケージのインストールから。
Tool →Install Packeges..といくとこんな画面になる。

またしても英語でちょっとイライラするが負けられない。

二行目の空欄のとこに”psych”と入力。
多分純正のRだといろいろ選べるみたいなのだが、しばらくはこのRstudioを使うので気にせずに行こう。

ちなみに3行目は、そのパッケージの保存場所。デフォルトではユーザー名/Library/R/2.13(バージョン名)/Library となっている。
作業ディレクトリに変更してもいいとは思うんだけど、まあデフォルトのしておこう。

実行するとこんな感じ。ほんとにインストールできてるのか確認。

よしどうやら入ってるみたい。

続いてこのパッケージを読み込むためのコマンド。
ちゃんと聞いて来たもんね。

library(psych)

これで、インストールしたpsychパッケージを読み出せる。

これを前回のデータの頭に入れてみる。ちょっと面倒なので、カテゴリカルデータの定義付けの部分は無視してみてみる。
library(psych)
#インストールしたパッケージを読みこみますよ。
lowdate1<-file(“nslpns.csv”,encoding=”Shift-JIS”)
lowdate2<-read.csv(lowdate1,header=T,sep=’,’)
lowdate2$sex<-factor(lowdate2$sex,levels=0:1,labels=c(“male”,”female”))
lowdate2$job<-factor(lowdate2$job,levels=0:1,labels=c(“ns”,”lpns”))
lowdate2$jobposition<-factor(lowdate2$jobposition,levels=0:2,labels=c(“highexecutive”,”executive”,”stuff”))
lowdate2$workingform<-factor(lowdate2$workingform,levels=0:4,labels=c(“day”,”2alternation”,”3alternation”,”unusual”,”parttime”))
#ここまででカテゴリー変数の定義完了。

この次に勉強会で習ったdescribe関数をいれる。この関数は簡単に言うと、グレードアップしたsummary関数という感じ。
とりあえず、年齢を入れてみよう。

describe(lowdate2$age)

よしできた。ちなみにパッケージを入れていないRでdescribe()って打ち込んでも「関数が見つからない」と返ってきます。パッケージには他にもいろいろあるみたいなのだが、とりあえずこのpsychのパッケージには心理系の分析に便利な関数がいろいろあるみたいなので、しばらくこのままのデータで練習してみようと思います。!

次回やること
1) describe.by
2) corr.test
3) multi.hist

ここまで追ってみよっと。
因子分析はもっと簡単なデータで改めてやることにしよっと。

分散分析までなんとかかんとかできるようになったので、本格的に自分の研究のデータを分析しようと思ったのだ。

本当は時間がないのでSPSSを使おうとしたのだけど、N先生から
“逃げ場を無くせ!!”との激。

ちなみにN先生は、ヒゲ付き棒グラフを美しく仕上げるために七転八倒している。

よーーしわかったよ。やるぜ!!

ということで、データを読み込むところから始める。

データの概要は、

だいたい300人分のデータで看護師さんと准看護師さんのデータ。
前半部は、年齢、性別、勤務形態などのデモグラフィックデータで後半部は3つの質問紙の質問項目のスコア。
ファイルは、nslpnsとしてcsv方式で保存(もちろん作業ディレクトリに)

これを読み込んだら、subset関数などを使ってソートを掛けながら少しずづ分析に向けて進めていくのだ。

エクセルで変数名を分かりやすく形成して、よし読み込むのみ。

Rでやれば、あとからどんなふうにデータを形成したのかがよくわかるのだ。

まずはRstudioを開いて、file→new→Rscript でエディタを開く。

今回は、看護師と准看護師さんのデータなので、nslpnsと名前をつけとこう。

エディタの左上のフロッピーマークをクリック。

こんな感じのが出てくるのでファイル名を付けて保存。
この時にwhereの部分をRの作業ディレクトリに変えておくのを忘れちゃだめ。

フフ。コマンドをエディタにいれてみる。
データの読取はもう出来るもんね。

s<- read.csv(“nslpns.csv”,header=T,sep=’,’)

!!

なんだよ。不正なマルチバイト文字って。

この後、約3時間ほどマルチバイト文字と戦う。

txt方式で読み込もうとしたり、csvデータをいじってみたり。

よくわかんないが、マルチバイト文字という奴は文字形式の問題らしいが、元データはきっちりしてるつもりなんだけどな。。。
UTFとshift JISの問題だとかなんとかようわからんので誰か教えてください。

いろいろ四苦八苦して、新たなコードをネット上で見つける。一応これで読み込めた。

lowdate1<-file(“nslpns.csv”,encoding=”Shift-JIS”)
lowdate2<-read.csv(lowdate1,,header=T,sep=’,’)

よくわかってないけど、たぶん一行目はcsvを読み込むときにShiftJISで読み込んでねということみたい。

膨大な量のデータなので、どっかで間違えて邪悪な文字を使ってしまってだんだろう。

ふう。燃え尽きた。続きはまた今度にしよう。

そうそう。今日は二元配置の分散分析にチャレンジ。
おっと対応はありませんね。

Rstdioを使ったことでずいぶん直感的な使い方ができるようになった。

二元配置の分散分析のコマンドは、summary(aov(s$weight~s$race*s$sex))

従属変数が最初で~でつないで*で交互にみたいな感じかな。
この*は:でもいいらしいがまあこのままで。

Rstudioのスクリーンショットを見てみる。

右上のエディタで作った二元配置の分散分析.Rを左上のコンソールに読み込んだ形式になる(RunAll)。

読み込むと左下ワークスペース(黒い円の部分)にデータフレームがセットされてる。

ここではコマンド1行目のs(いわゆるローデータ)とコマンド7目のs2(分散分析の結果)がデータセットされていると考えたらいいのだろう。

この読み込んだデータフレームは、保存をするとずっとそのままになる。
コンソールの消去は、とりあえず机の上を片付けるようなものだけど、ワークスペースの消去は頭の中をリセットする様なものかな。

ところで、結果だが、p値を見ると、このデータではsexのみが有意。

だから、体重に関しては性の主効果は有意、人種の主効果、及び交互作用はなしという見方になる。

酔っ払っているので今日はここまで。

ちょっとanobakun http://www11.atpages.jp/~riseki/pukiwikiplus/index.php?ANOVA%B7%AF

をやってみたのだけど、これは詳しい人が作成したプログラムを読み込むタイプなのでまた別カテゴリーでやってみます。

自分でできないことだけど、理解しつつやろうと思う今日この頃。
奥が深いな。

ということで怒涛の更新。

前回の設定でもう一度開いてみる。

ほうほう。作業ディレクトリの中にあるファイルは全て右下で管理してくれている。左下は正しくコンソールなんだな。

えっとこのめんどくさいブツブツ書いてあるのを消すのはどうだったけ??

edit→Clearconsole でOK。

配置に関して、多分それぞれいろいろやり方があるんだろうけど、僕の好みにちょっと変えてみよう。おおおサイズも結構変わるんだね。

変更は前回のRstdio→Prefarances→Panelayoutで。

僕は今まで流れに沿って、左にコンソール。右にエディタを配置。ワークスペースブラウザとFileはもともと使っていなかったのでちょっと小さめに配置してみた。

さて、これでちょっと昔ののデータを読み込んでみよう。

おおすげー。右上にタブが出来てる。Spacesに近い感じなんだろう。ということは同時に何個もプラグラムファイルを開けるんだな。何個か開いてみる。

うん。これは間違いなく使いやすい。お金が取れるレベルだ。

では次にちょとRを走らせてみよう。前回やったt検定と一元配置の分散分析にチャレンジ。

エディタ右上のRun Allをクリック。ちなみにRun Lineは読んで字のごとく一行のみコマンドを実行するみたいだ。

左のコンソールに source(‘~/) とかなんとか書いてある。
これはエディタの内容を読み込んだという状態のことを指すのだろう。

!!

左下になんかある!!そうかそういうことか、このワークスペースは、コンソールに読み込んでいる内容を分かりやすく表示してくれているんだな!これは便利。

試しにコンソールにコマンドを打ち込んでみる。

おみごと。前と同じデータが出てきました。パチパチ。

このまま保存すれば、次回はこの状態のままワークスペースは復活する。もし最初からしたければ、
WorkSpace →Clear all でOK。

いやあRstudio恐るべし。

多分これでRstudioはほぼ大丈夫かな?

(追加)
ちょっと追加
Rstudio → Preferences → Editing でHighlight selected line
にすれば、エディタの実行中にコマンドにラインがはいるので便利。

本当は二元配置の分散分析を仕様と思ったのだが、話題のRstudioをインストールしてみる。

ダウンロード先はこちら。

http://www.rstudio.org/

こんな画面があって。

こんな画面があるので、デスクトップの方をダウンロード。

続いてまたしても英語でちょっと緊張するが、リンゴのマークを信じましょう。

さてダウンロードしたファイルをインストールしてみる。アプリのマークはRをちょっと可愛らしくしたような感じ。

このアプリそのものはRに乗っかるような感じでできているらしく、もともとのRは自動的に読み取られるみたいだ。

早速起動。

うん。なんか出てきた。結構直感的なインターフェースだな。。。

でも全部英語だよ。こんちくしょう。

たちまち環境設定をみてみよう。
Rstudio→Preferencesの画像。

まず一行目のDefault CRAN mirrorがわからん。なんなんじゃろこれ?

2行目のInitial working directoryは多分作業ディレクトリの変更。これは前と同じようにデスクトップの中のRフォルダにしてみる。

その下のRestore .RData into workspace at startup は最初からなので、多分設定の保存だろう。

Restore .RData into workspace at startup

は保存の設定。なにも考えずaskのままでいいと思う。

他のとこもいろいろ見てみたが正直よくわからんなあ。とくにAppearance のとこはちんぷんかんぷん。

ただ多分、このRstudioのキモは Pane Layout になるんだろうな。

左上にプログラムソースが開き、その下がコンソール。
右上がワークスペースブラウザ(未だにいまいち意味が分かっていないのだが)でその下にちょっとしたグラフか。

使いこなしたら結構楽しそうだな。2元配置の分散分析はちょっとおやすみしようかな。。。

t検定について考えていたらまたわからんくなったので、改めてやって見ることにした。

最初の頃の気持ちを思い出して丁寧にやっていく。だって結果がでて喜んでるだけでは結局使いにくいSPSSと同じだもん。ちゃんと勉強しましょう。

まずは、等分散性の検定について検討(F検定と呼ばれるもの)。

t検定とは違うところは、複数あることか。ちょっとおんなじ感じでやってみっと。3種類ある人種でトライ。

むむむエラーだ。なんでだろう。ちょっと本を読んでみよう。

わかった(多分ね)!!

3群以上の時は、bartlettを使うみたい。

bartlett.test(s$weight,s$race)

さあ出てこい!!

いえーーす!!これもp-value = 0.04765なので等分散性が棄却されます。サンプル数が少ないもんね。

ということで、等分散を規定しない分散分析。

ひつこいけどこれがデフォルトなんだよね。oneway.test

っとこれをみると、 p-value =0.3214 なので体重に人種の影響はないことがわかった。

一応、等分散を規定した上での分散分析も練習のためにしてみる。
t検定の時と同じでvar.equal=TRUEを付ける。

うん。p-value = 0.138で結果はかわらない。でも思ったより数字に変化があるんだな。

ところで、Ｒで分散分析には他にもいろんな種類があるみたいだ。本で読んだ限りでは、
aov関数(summaryも追加して結構一般的に使える、2元配置もできる)
anova関数(結構高度らしい)
anovakun関数(ツイッターで教えてもらったが、まだちょっと僕には早そうだ)

ここらへんを一応やってみよう。
ちなみに今気づいたのだが、summaryという関数はなんかざっくり見やすくしてくれるコマンドみたい。
試しにちょっとsをぶち込んでみる。

ほほう。なんかいい感じの表ができてる。たぶん、customとInheは順序尺度であることを一度どこかで定義付けせんとあかんみたいやな(なんで大阪弁だよ)。これはまた今後にしよっと。

ではまず、aovから。比較のためoneway.testの等分散とそうじゃないやつを上に並べてみた。

!! あれ?aovはデフォルトが等分散を規定しているのか??なぞだなあ。

ちなみにanovaだとこうなる。

こっちもそうだ。僕の中の認識では、デフォルトが等分散を規定していなくて、等分散よりも統計的なハードルが高いので気にしなくても大丈夫という感じだったんだけどちょっと違うのかな。。。
これは調べても見つからない。誰かに聞くべきか。ということで悩んでいたらツイッター上でmentionが。そもそもanovaは等分散を想定しているものだから、等分散以外はやっちゃいけないとのこと。

そういう場合は、oneway.testでやるのが普通らしい。
でもだったら2元配置の分散分析をするときはどうなんだろう?twoway.testみたいなのがあるのかな?
ということで次は二元配置の分散分析にチャレンジ。そして、anovakunもやってみよっと。

ところで、分散分析ではTukeyの下位検定が出てくるものが多い。
Rではこうなる。

TukeyHSD(aov(s$weight~s$race))

ok!だいぶいい感じですね。

ちょっといまいちわかってないのでもう一度最初からやってみる。

まずは、等分散性の検定について検討(F検定と呼ばれるもの)。
僕が習った方法としては、t検定では(分散分析でも)まず実際の分析にかける前に、このF検定を行って調べたい両者が同じ分散を持っているか、否かを判断する必要があるという。

もちろん、Ｒにもその関数はあった。

var.test というものだ。

今回のデータでは、体重のばらつきが、男性でも女性でも一緒なのか?という意味になる。バラつきが一緒な場合と、ばらつきが違う場合では検定の方法が違いますよということらしい。

(これは青木先生のページを見ると、必ずしもやらなくても最初から等分散を規定しない方法でやればいいとのことだが、難しいことはよくわからないので一応原則通りやってみる)

早速やってみよう。

ふむ。なんか出てきた。だいぶ分かってきた感じもするなあ。

これの見方は、えっとp-value = 0.03822 だから。。。

これは当分散性が棄却されたということに違いない(間違っていたら誰か教えてくれ)。

ということで、この場合には体重のばらつきは男性と女性で同じではないので、そのための分析をしなさいね。ということになるのだ。

Ｒではデフォルトでウェルチなので(等分散を規定しない)

t.test(s$weight~s$sex) となりました。

これは、 p-value = 0.001493 なので有意。男女間で体重には有意な隔たりがあるということだな。

※こういったデータの場合は普通は当分散になるはずだけど(実際に僕の少ない経験ではここが棄却されることはほとんどない。もし、棄却されたらなんか間違ってるかと思って見直す感じ)たぶん、この練習データはあまりにも人数が少ないので棄却されたんだろうと思う。

ということで、もし等分散性が棄却されなかったらどうするか?

その場合には、var.equal=TRUE を追加すればいい。

確かにあんまり数字としては変わらないね。。。

状況に応じて使い分けるという感じかな。。

だれかわかり安く教えてくれる人がいればお願いします。