2011/07/29

RでMixed Model メモ

杉本隆さん（おすぎ）から提供いただいた、前回のRを動かした回のメモをアップします！
おすぎ、ありがとう！！

皆様、書式が少しずれてしまっているのですが、まずはアップしてしまいます。




更新：2011/06/09

【コンセプト】
□データ解析の実務者の視点からRを使用
□データを取り扱う際に遭遇する困難を明らかに
□実例プログラムは必ず入れる
□安定第一という観点からR本体で全て進める（コマンダーは参照程度）

0. プログラムを書き始めるために

Rを起動して最初に出てくるR Consoleは，結果の出力と一緒になっていて，プログラムを加工しにくい。
メニューの[ファイル]→[新しいスクリプト]でRエディタを起動する。
エディタに書いたプログラムを実行するときは，選択して[Ctrl+r]をするか，R Consoleにコピペして実行


1. データの読み込み

データの種類には色々あるが、現実場面でよく使われそうな形式として、Excelに入力したデータ（csvファイルに変換したもの）やSPSSのデータファイル（savファイル）を取り上げる。まず，csvファイルの読み込みを解説する。csvファイルを使うのは，RでExcelファイル（xlsファイル）をそのまま読もうとすると大変な理由による。

まず，読み込むときにフォルダのパス名を教えてあげなければいけないのだが，普段はこれは面倒すぎる
ので，省略するために，「作業ディレクトリ」という，データの置いてあるフォルダの指定を行う

作業ディレクトリを変更する場合は
メニューの[ファイル]→[ディレクトリの変更]でフォルダを指定。デスクトップにフォルダを作ると分かりやすいかも。

Rは，データを「オブジェクト」として認識する。基本的には，<-の矢印が向けられた先にあるものが
オブジェクトとなる。名前は自由。


そして，今度はcsvファイルをRに読み込む。これからデータ加工を始めるので，順番に作っていくことを考えて，オブジェクト名を「a」とした。
作業ディレクトリを指定しているので，"jgss_mini.csv"と，データファイル名だけでOK

a <- read.csv("jgss_mini.csv", header=T) #　header=Tは一行目が変数名という意味

ファイル名打つのが面倒ならば

a <- read.csv(file.choose(), header=T)



読み込んだデータセット確認のため


view(a)→デフォルトでは×
head(a) # 最初の6人表示
fix(a) or edit(a)

オブジェクト確認
ls()

オブジェクト削除するときは
rm(a)

全部削除は
rm(list=ls())

どうせ1つのデータセットで集中的に作業するんだから，attachでいいんじゃないか。

-----------------------------------------------

とりあえず，度数分布でも確認してみたい
table(a$DOPET)
%を出すには
prop.table(table(a$DOPET))


クロス
table(a$DOPET,a$SEXA)


library("Hmisc")
describe(a)
describe(a$DOPET)→n，欠損値，平均，度数分布まで出る。factorだろうがnon-factorだろうが関係ないよ

library("gmodels")


Frequencies and Crosstabs
http://www.statmethods.net/stats/frequencies.html

--------------------------------------------------

2. データマネジメント

(1)数量型変数を因子型変数に変換，またその逆
全て数量変数として読み込まれているので，カテゴリとして扱いたい変数を因子型に変換する。Rはこれらの区別を厳密につけたがるため

a$CASEID <- as.factor(a$CASEID)
a$SEXA <- as.factor(a$SEXA)
a$MARC <- as.factor(a$MARC)
a$XJOB1WK <- as.factor(a$XJOB1WK)
a$DOPET <- as.factor(a$DOPET)
a$DODOGE <- as.factor(a$DODOGE)
a$DODOGI <- as.factor(a$DODOGI)
a$DOCAT <- as.factor(a$DOCAT)
a$DOMAMMAL <- as.factor(a$DOMAMMAL)
a$DOBIRDA <- as.factor(a$DOBIRDA)
a$DOFISHA <- as.factor(a$DOFISHA)
a$DOREPTL <- as.factor(a$DOREPTL)
a$DOINSECT <- as.factor(a$DOINSECT)
a$DOOTHER <- as.factor(a$DOOTHER)
a$DOPNOMK <- as.factor(a$DOPNOMK)

summary(a)

で確認


以上は1つずつ変換したが，↓の3行を使えば，それだけでできる。
作業内容は，オブジェクト該当変数をfactorに変換した「af」を作り，aとmergeして，新データセット「b」を作成


af <- sapply(data.frame(a[,1:2], a[,4:6], a[,8:19]), factor)
af <- data.frame(af)
b <- merge(af, a, by=names(af), all=T) #aとafの順番を間違えないように

summary(b)

で確認


因子から数値，NAが入っている数値データはそれをやる必要性あり




(2)欠損値の定義
(3)欠損値の除去
(4)連続変数をカテゴリ変数に（例：年齢層作成）
(5)逆転項目の定義




--------------------------------

[応用]

(1)

SPSSのデータファイルを読み込んで、汎用性の高いcsvに変換して保存し，さらにそれを読み込むという作業

まず、他のデータ形式を読み込むための準備として、foreignパッケージを読み込む

library(foreign)

そして、データファイルの置いてあるフォルダを指定してデータセット（jgss_mini.sav）を読み込む
ここでは「jgss」という名前のオブジェクトで

jgss <- read.spss("データファイルのフォルダまでのパス/jgss_mini.sav", to.data.frame=T, use.value.labels=F)

しかし，パス指定は最初はハードルが高いので，以下の作業ディレクトリの指定の方法を薦める


*[作業ディレクトリの場所指定が面倒くさい場合]

作業ディレクトリの場所にデータファイルがあれば，フォルダのパス名を省略できる。
現在の場所は

getwd()

で確認できる

プログラムでは，

setwd("ここにフォルダのパス/ファイル名")

であり，フォルダのパスは，フォルダのアドレス欄からコピペできるが，
Windowsでは，フォルダの区切り記号が「\」なため，これを「/」に変換する必要
があり，少し面倒くさい。

これで、場所指定が省略されてファイル名だけで作業できるようになる。

jgss <- read.spss("jgss_mini.sav", to.data.frame=T, use.value.labels=F)
--------------

次に，読み込んだsavファイルをcsvファイルに変換して保存する

write.table(jgss, "jgss_mini.csv", sep=",", row.names=FALSE)

あとはcsvの読み込みと同じ



20110610 mixed model multilevel model

> results1 <- lmer(mathach ~ 1 + (1| schid), data = a )
> #これはnull modelを表現している！
> ＃スクールによって切片が違うかを見ている
>
> summary(results1)
Linear mixed model fit by REML
Formula: mathach ~ 1 + (1 | schid)
Data: a
AIC BIC logLik deviance REMLdev
47123 47143 -23558 47116 47117
Random effects:
Groups Name Variance Std.Dev.
schid (Intercept) 8.614 2.9350
Residual 39.148 6.2569
Number of obs: 7185, groups: schid, 160

Fixed effects:
Estimate Std. Error t value
(Intercept) 12.6370 0.2443 51.72

> #deviance値でモデルの適合度を判断するのが通例らしい
> ＃Interclass correlation coefficient (ICC)＝8.614/8.614+39.148である。σ＾2＝39.148。

> ＃ICC＝0.18はだいぶ高いらしい。学校によって、18%の効果がある→mixed modelを使用する価値がある

> ＃学校の数はある程度（20とか？）ないと、信頼性が怪しくなる（３校とかでは怪しいのではないか・・）

>
> results2 <- lmer(mathach ~ 1 + meanses + schtype + (1 | schid), data = a)
> summary(results2)
Linear mixed model fit by REML
Formula: mathach ~ 1 + meanses + schtype + (1 | schid)
Data: a
AIC BIC logLik deviance REMLdev
46956 46991 -23473 46944 46946
Random effects:
Groups Name Variance Std.Dev.
schid (Intercept) 2.3139 1.5211
Residual 39.1614 6.2579
Number of obs: 7185, groups: schid, 160

Fixed effects:
Estimate Std. Error t value
(Intercept) 12.0963 0.1986 60.90
meanses 5.3329 0.3686 14.47
schtype 1.2253 0.3058 4.01

Correlation of Fixed Effects:
(Intr) meanss
meanses 0.245
schtype -0.698 -0.356


> #devianceが減っているので、より適切なモデルであることがわかる
> ＃ランダム切片も、8.614から2.31に減っているのでより説明できている
> ＃meansesが１ポイント上がると、5点くらい数学の点数が上がる
> ＃lme4の作者が、自由度が複雑なので、Ｐ値は出さないように設計した


> results3 <- lmer(mathach ~ meanses + schtype + ses + meanses*ses + schtype*ses + (1 + ses | schid), data = a )
> summary(results3)
Linear mixed model fit by REML
Formula: mathach ~ meanses + schtype + ses + meanses * ses + schtype * ses + (1 + ses | schid)
Data: a
AIC BIC logLik deviance REMLdev
46525 46594 -23253 46498 46505
Random effects:
Groups Name Variance Std.Dev. Corr
schid (Intercept) 2.407424 1.55159
ses 0.014619 0.12091 1.000
Residual 36.758381 6.06287
Number of obs: 7185, groups: schid, 160

Fixed effects:
Estimate Std. Error t value
(Intercept) 12.0948 0.2028 59.65
meanses 3.3202 0.3885 8.55
schtype 1.1949 0.3079 3.88
ses 2.9052 0.1483 19.59
meanses:ses 0.8462 0.2718 3.11
schtype:ses -1.5781 0.2245 -7.03

Correlation of Fixed Effects:
(Intr) meanss schtyp ses mnss:s
meanses 0.213
schtype -0.676 -0.344
ses 0.076 -0.146 -0.064
meanses:ses -0.143 0.178 -0.081 0.278
schtype:ses -0.062 -0.080 0.093 -0.679 -0.356



＃Rstudioをインストールすると便利！

投稿者 時刻: 0 件のコメント:
ラベル: ,

次回は9月2日（金）16:00-　Rの出力から結果を読む

前回（7/29金）のエビ研でも、Mixed Model (Multilevel Model)の解析をRを用いて行いました。
1ヶ月あいてしまうと、データセットをRに読み込むところから躓いてしまい、解析にかかるまで時間がかかってしまいましたが、なんとか解析をすることができました。

次回は、この出力を、どのように読んで、論文等の表にするときにはどのように記述するのかをやってみたいと思います。

日　時：2011年9月2日（金）16:00～18:00
テーマ：Mixed Model (18)
場　所：医学部３号館　3階S308

前回演習で行った出力（次回に使う出力）を松長さんが送って下さいました。松長さんありがとうございます！

> a <- read.csv(file.choose(), header = T)
> attach(a)
> a$centses <- ses - meanses
> library(lme4)
要求されたパッケージ Matrix をロード中です
要求されたパッケージ lattice をロード中です

次のパッケージを付け加えます: 'Matrix'

The following object(s) are masked from 'package:base':

det


次のパッケージを付け加えます: 'lme4'

The following object(s) are masked from 'package:stats':

AIC

> results1 <- lmer(mathach ~ 1 + (1 | schid), data = a)
> summary(results1)

Linear mixed model fit by REML
Formula: mathach ~ 1 + (1 | schid)

Data: a   
AICBIClogLikdevianceREMLdev
4712347143-235584711647117

Random effects:
GroupsNameVarianceStd.Dev.
schid(Intercept)8.6142.9350
Residual39.1486.2569
Number of obs: 7185, groups: schid, 160

Fixed effects:   
EstimateStd. Errort value
(Intercept)12.63700.244351.72

> results2 <- lmer(mathach ~ 1 + meanses + schtype + (1 | schid), data = a)
> summary(results2)

Linear mixed model fit by REML
Formula: mathach ~ 1 + meanses + schtype + (1 | schid)

Data: a     
AICBIClogLikdevianceREMLdev
4695646991-234734694446946

Random effects:
GroupsNameVarianceStd.Dev.
schid(Intercept)2.31391.5211
Residual39.16146.2579
Number of obs: 7185, groups: schid, 160

Fixed effects: 
EstimateStd. Errort value
(Intercept)12.09630.198660.90
meanses5.33290.368614.47
schtype1.22530.30584.01

Correlation of Fixed Effects:
(Intr)meanss
meanses0.245 
schtype-0.698-0.356

> results3 <- lmer(mathach ~ meanses + schtype + centses + meanses*centses + schtype*centses + (1 + centses|schid), data = a)
> summary(results3)

Linear mixed model fit by REML
Formula: mathach ~ meanses + schtype + centses + meanses * centses + schtype * centses + (1 + centses | schid)
Data: a
AICBIClogLikdevianceREMLdev
4652446592-232524649646504
    
Random effects:
GroupsNameVarianceStd.Dev. Corr
schid(Intercept)2.381861.54333
centses0.101380.318410.392
Residual36.721136.05980
Number of obs: 7185, groups: schid, 160

Fixed effects:
EstimateStd. Error t value
(Intercept)12.11360.198860.93
meanses5.33910.369314.46
schtype1.21670.3064 3.97
centses 2.93880.155118.95
meanses:centses1.03890.2989 3.48
schtype:centses-1.6426 0.2398 -6.85

Correlation of Fixed Effects:
(Intr)meanssschtypcentssmnss:c
meanses0.245
schtype-0.697-0.356
centses0.0800.020-0.056
menss:cntss0.0190.079-0.0280.282 
schtyp:cnts -0.055 -0.0290.082-0.694-0.351


投稿者 時刻: 0 件のコメント:
ラベル: , ,

2011/07/28

次回7月29日（金）もRで実習

前回（6/10金）のエビ研では、実際にMixed Model (Multilevel Model)の解析をRを用いて行いました。

というわけで、次回もまた、Rを用いて実際に動かしてみたいと思います。

日　時：2011年7月29日（金）18:00～20:00
テーマ：Mixed Model (17)
場　所：医学部３号館　3階S308

前回と同じデータセットとRをご持参いただければと思います。
（前回の予告記事をご参照くださいませ。）

また、Rでデータを読み込むときの手順を土屋政雄さんのご厚意でここに掲載させていただきます。
つっちーどうもありがとうございます！
------------------------------------------------------
【データを読み込むまで】

まず，作業ディレクトリ（フォルダ）の変更
メニューの[ファイル]→[ディレクトリの変更]でフォルダを指定。デスクトップにフォルダを作ると分かりやすいかも。名前は英数字で（フォルダパスの中に日本語が入るとうまくいかない事があるため）。

getwd() #現在の作業フォルダの確認

a<-read.csv("hlm_sem.csv", header=T) # 「hlm_sem.csv」というデータファイルを，「a」という名前でRに読み込む
b<-read.csv(file.choose(), header=T) # こちらを使うと，ファイル名を入れずにすみ，ダイアログから選択できる

fix(a) #spss風な画面で読み込まれたことを確認
------------------------------------------------------

ご都合のつく方は是非いらしてくださいませ。（文責：宮本有紀）

投稿者 時刻: 0 件のコメント:
ラベル: , ,
登録: 投稿 (Atom)