前回（12/16）は共分散構造分析の第一回でした。
米倉佑貴さんから共分散構造分析の基礎から教えていただいています。
わかりやすい！です。
また記事を書きたいと思いますが、とりあえず次回日程のみメモ。
2012年1月13日（金）
共分散構造分析(2)
です。
2011/12/25
2011/11/24
次回は12月16日（金）18:00- 共分散構造分析
ブログアップがなかなか追いつかず、気付くと広告（1ヶ月以上更新がないと出てきちゃう）が表示されてしまっています。。
エビ研で長らく取り上げてきたMixed modelあるいはマルチレベルモデルが前回で一区切りとなりました。
前回は、土屋政雄さんがご自身の資料を用いて線形混合モデルについてお話ししてくださり、みなで入れ子構造についてなどなど話し合いました。
また、土屋さんはこの資料（線形混合モデルの入門スライド）を一般に公開してくださっています！
つっちー、どうもありがとう！！！
また、データを皆で共有したいときの話について、米倉佑貴さんがInter-University Consortium for Political and Social Research (ICPSR)のGuide to Social Science Data Preparation and Archiving （社会科学データの準備と保存の手引き）の紹介をしてくださいました。
原版はICPSRのウェブサイトの中から、日本語版は東京大学社会科学研究所の中の社会調査・データアーカイブ研究センター内 (http://ssjda.iss.u-tokyo.ac.jp/text2.html) からダウンロードすることができます。
さて、次回からは共分散構造分析に入ります。よねちゃん、どうぞよろしくお願いいたします。
日 時：2011年12月16日（金）18:00～20:00
テーマ：共分散構造分析 (1)
場 所：医学部３号館 3階S308
ご都合のつく方は是非いらしてくださいませ。（文責：宮本有紀）
エビ研で長らく取り上げてきたMixed modelあるいはマルチレベルモデルが前回で一区切りとなりました。
前回は、土屋政雄さんがご自身の資料を用いて線形混合モデルについてお話ししてくださり、みなで入れ子構造についてなどなど話し合いました。
また、土屋さんはこの資料（線形混合モデルの入門スライド）を一般に公開してくださっています！
つっちー、どうもありがとう！！！
また、データを皆で共有したいときの話について、米倉佑貴さんがInter-University Consortium for Political and Social Research (ICPSR)のGuide to Social Science Data Preparation and Archiving （社会科学データの準備と保存の手引き）の紹介をしてくださいました。
原版はICPSRのウェブサイトの中から、日本語版は東京大学社会科学研究所の中の社会調査・データアーカイブ研究センター内 (http://ssjda.iss.u-tokyo.ac.jp/text2.html) からダウンロードすることができます。
さて、次回からは共分散構造分析に入ります。よねちゃん、どうぞよろしくお願いいたします。
日 時：2011年12月16日（金）18:00～20:00
テーマ：共分散構造分析 (1)
場 所：医学部３号館 3階S308
ご都合のつく方は是非いらしてくださいませ。（文責：宮本有紀）
2011/10/04
次回は10月28日（金）18:00- Rの出力から結果を読むのと
前回やったことも書きたいのですが、とりあえず次回のエビ研日程だけ。
次回も、Mixed Model (Multilevel Model)の解析の出力を、どのように読んで、論文等の表にするときにはどのように記述するのかをやってみたいと思います。
日 時：2011年10月28日（金）18:00～20:00
テーマ：Mixed Model (19)
場 所：医学部３号館 3階S308
ご都合のつく方は是非いらしてくださいませ。（文責：宮本有紀）
次回も、Mixed Model (Multilevel Model)の解析の出力を、どのように読んで、論文等の表にするときにはどのように記述するのかをやってみたいと思います。
日 時：2011年10月28日（金）18:00～20:00
テーマ：Mixed Model (19)
場 所：医学部３号館 3階S308
ご都合のつく方は是非いらしてくださいませ。（文責：宮本有紀）
2011/09/01
次回は9月2日（金）16:00- Rの出力から結果を読む（再掲）
次回エビ研情報を再掲します。
次回は、Mixed Model (Multilevel Model)の解析の出力を、どのように読んで、論文等の表にするときにはどのように記述するのかをやってみたいと思います。
日 時：2011年9月2日（金）16:00～18:00
テーマ：Mixed Model (18)
場 所：医学部３号館 3階S308
ご都合のつく方は是非いらしてくださいませ。（文責：宮本有紀）
次回は、Mixed Model (Multilevel Model)の解析の出力を、どのように読んで、論文等の表にするときにはどのように記述するのかをやってみたいと思います。
日 時：2011年9月2日（金）16:00～18:00
テーマ：Mixed Model (18)
場 所：医学部３号館 3階S308
ご都合のつく方は是非いらしてくださいませ。（文責：宮本有紀）
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をインストールすると便利！
おすぎ、ありがとう！！
皆様、書式が少しずれてしまっているのですが、まずはアップしてしまいます。
更新：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をインストールすると便利！
次回は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
Random effects:
Number of obs: 7185, groups: schid, 160
Fixed effects:
> 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
Random effects:
Number of obs: 7185, groups: schid, 160
Fixed effects:
Correlation of Fixed Effects:
> 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
Random effects:
Number of obs: 7185, groups: schid, 160
Fixed effects:
Correlation of Fixed Effects:
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
| 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 |
Fixed effects:
| Estimate | Std. Error | t value | |
| (Intercept) | 12.6370 | 0.2443 | 51.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
| 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 |
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 |
> 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
| AIC | BIC | logLik | deviance | REMLdev |
| 46524 | 46592 | -23252 | 46496 | 46504 |
Random effects:
| Groups | Name | Variance | Std.Dev. | Corr |
| schid | (Intercept) | 2.38186 | 1.54333 | |
| centses | 0.10138 | 0.31841 | 0.392 | |
| Residual | 36.72113 | 6.05980 |
Fixed effects:
| Estimate | Std. Error | t value | |
| (Intercept) | 12.1136 | 0.1988 | 60.93 |
| meanses | 5.3391 | 0.3693 | 14.46 |
| schtype | 1.2167 | 0.3064 | 3.97 |
| centses | 2.9388 | 0.1551 | 18.95 |
| meanses:centses | 1.0389 | 0.2989 | 3.48 |
| schtype:centses | -1.6426 | 0.2398 | -6.85 |
Correlation of Fixed Effects:
| (Intr) | meanss | schtyp | centss | mnss:c | |
| meanses | 0.245 | ||||
| schtype | -0.697 | -0.356 | |||
| centses | 0.080 | 0.020 | -0.056 | ||
| menss:cntss | 0.019 | 0.079 | -0.028 | 0.282 | |
| schtyp:cnts | -0.055 | -0.029 | 0.082 | -0.694 | -0.351 |
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風な画面で読み込まれたことを確認
------------------------------------------------------
ご都合のつく方は是非いらしてくださいませ。（文責：宮本有紀）
というわけで、次回もまた、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風な画面で読み込まれたことを確認
------------------------------------------------------
ご都合のつく方は是非いらしてくださいませ。（文責：宮本有紀）
2011/06/02
次回6月10日（金）はRで実習
次回のエビ研では、実際にMixed Model (Multilevel Model)の解析を
皆でしてみようというご提案をいただき、是非したい！ということになりました。
参考文献として下記を使います。
Estimating Multilevel Models using SPSS, Stata, SAS, and R
JeremyJ.Albright and Dani M. Marinova
http://www.indiana.edu/~statmath/stat/all/hlm/hlm.pdf
できれば一緒に行いたいので、解析はフリーソフトのRを使いたいと思います。
データセットは，このなかの↓This is the data set.のdata setを用います。
http://www.ats.ucla.edu/stat/hlm/seminars/hlm6/mlm_hlm6_seminar.htm
このデータセットはSPSSのファイルなので，可能ならばCSVに変換しておいてください。
ご自分でRとパッケージのlme4もインストールしておいていただければ幸いです。
一連のやり方が分からなかったら，「r パッケージ インストール」
と検索してみれば，たいてい画像付きの解説が出てくるようです。
SPSSのデータをCSV形式に変換できなかったり
上記の手順のどこかができなかった場合には宮本までご連絡くださいませ。
当日は、可能な限りPC持参で、文献のプリントアウトが必要な方は各自でご準備いただければ幸いです。
日 時：2011年6月10日（金）18:00～20:00
テーマ：Mixed Model (16)
場 所：医学部３号館 3階S308
ご都合のつく方は是非いらしてくださいませ。（文責：宮本有紀）
皆でしてみようというご提案をいただき、是非したい！ということになりました。
参考文献として下記を使います。
Estimating Multilevel Models using SPSS, Stata, SAS, and R
JeremyJ.Albright and Dani M. Marinova
http://www.indiana.edu/~statmath/stat/all/hlm/hlm.pdf
できれば一緒に行いたいので、解析はフリーソフトのRを使いたいと思います。
データセットは，このなかの↓This is the data set.のdata setを用います。
http://www.ats.ucla.edu/stat/hlm/seminars/hlm6/mlm_hlm6_seminar.htm
このデータセットはSPSSのファイルなので，可能ならばCSVに変換しておいてください。
ご自分でRとパッケージのlme4もインストールしておいていただければ幸いです。
一連のやり方が分からなかったら，「r パッケージ インストール」
と検索してみれば，たいてい画像付きの解説が出てくるようです。
SPSSのデータをCSV形式に変換できなかったり
上記の手順のどこかができなかった場合には宮本までご連絡くださいませ。
当日は、可能な限りPC持参で、文献のプリントアウトが必要な方は各自でご準備いただければ幸いです。
日 時：2011年6月10日（金）18:00～20:00
テーマ：Mixed Model (16)
場 所：医学部３号館 3階S308
ご都合のつく方は是非いらしてくださいませ。（文責：宮本有紀）
2011/05/31
次回Mixed Model続きの日程は6月10日（金）です
みなさまこんにちは。前回3/23からだいぶあいてしまいました。
次回のエビ研日程をお知らせします。
日 時：2011年6月10日（金）18:00～20:00
テーマ：Mixed Model (16)
場 所：医学部３号館 3階S308
内容は、Mixed Modelです。
方法がわかりやすく書いてある論文をお持ちの方は是非宮本までご一報ください。
当日持ち込みも歓迎です。
また、現在計画中の研究や、取ったデータで解析方法について考えるのも歓迎です。
ご都合のつく方は是非いらしてくださいませ。（文責：宮本有紀）
次回のエビ研日程をお知らせします。
日 時：2011年6月10日（金）18:00～20:00
テーマ：Mixed Model (16)
場 所：医学部３号館 3階S308
内容は、Mixed Modelです。
方法がわかりやすく書いてある論文をお持ちの方は是非宮本までご一報ください。
当日持ち込みも歓迎です。
また、現在計画中の研究や、取ったデータで解析方法について考えるのも歓迎です。
ご都合のつく方は是非いらしてくださいませ。（文責：宮本有紀）
2011/03/25
次回日程は4月に決めたいと思います
前回（3月23日（水）13:30～）の勉強会にいらしてくださった皆様、どうもありがとうございました。前回はマルチレベル解析を行った論文を読み、さらに、メンバーが現在執筆中の論文の表を見せていただきました。それぞれ、いろいろなご活動をなさりながら、寒い中ご参加いただきありがとうございました。
次回開催日程は、新年度のみなさまの予定と今後の電力の供給状況等を考慮し、4月以降に決定したいと思います。
今のところ、4月～5月のあいだでまずは一度開催したいと考えております。
ご要望やお問い合わせがありましたら、宮本までご連絡くださいませ。
（文責：宮本有紀）
次回開催日程は、新年度のみなさまの予定と今後の電力の供給状況等を考慮し、4月以降に決定したいと思います。
今のところ、4月～5月のあいだでまずは一度開催したいと考えております。
ご要望やお問い合わせがありましたら、宮本までご連絡くださいませ。
（文責：宮本有紀）
2011/03/22
時間変更
みなさまたびたび申し訳ありません。
明日のエビ研は、計画停電等の可能性を考慮し、18時から予定していた勉強会を13時30分からに変更させていただきます。
夜の時間で予定していてくださっていた方には申し訳ありませんが、ご了承いただければ幸いです。
日 時：2011年3月23日（水）13:30～15:30
テーマ：Mixed Model (15)
場 所：医学部３号館 3階S308
内容は、
1. メンバーの研究の解析について
2. 残りの時間は下記の論文を読む予定
O'Campo P, Caughy MO, Nettles SM. Partner abuse or violence, parenting and neighborhood influences on children's behavioral problems. Soc Sci Med. 2010 May;70(9):1404-15. Epub 2010 Feb 12.
PMID: 20163906
です。
ご都合のつく方は是非いらしてくださいませ。（文責：宮本有紀）
明日のエビ研は、計画停電等の可能性を考慮し、18時から予定していた勉強会を13時30分からに変更させていただきます。
夜の時間で予定していてくださっていた方には申し訳ありませんが、ご了承いただければ幸いです。
日 時：2011年3月23日（水）13:30～15:30
テーマ：Mixed Model (15)
場 所：医学部３号館 3階S308
内容は、
1. メンバーの研究の解析について
2. 残りの時間は下記の論文を読む予定
O'Campo P, Caughy MO, Nettles SM. Partner abuse or violence, parenting and neighborhood influences on children's behavioral problems. Soc Sci Med. 2010 May;70(9):1404-15. Epub 2010 Feb 12.
PMID: 20163906
です。
ご都合のつく方は是非いらしてくださいませ。（文責：宮本有紀）
2011/03/17
次回Mixed Model続きの日程は3月23日（水）です
みなさま、先日、次回日程を2月23日と予告していましたが、諸般の事情により、一ヶ月後の3月23日に変更させていただきました。予定を空けてくださっていた方、大変申し訳ありませんでした。
また、次回は、Mixed Modelの論文を読みつつ、この春論文を提出された皆さんやご卒業される皆さんをねぎらい祝うことも同時に行いたいと思っています。
次回みんなで読んでみる予定にしているのは下記の論文です。
O'Campo P, Caughy MO, Nettles SM. Partner abuse or violence, parenting and neighborhood influences on children's behavioral problems. Soc Sci Med. 2010 May;70(9):1404-15. Epub 2010 Feb 12.
PMID: 20163906
日 時：2011年3月23日（水）18:00～20:00（→13:30～に時間変更）
テーマ：Mixed Model (15)
場 所：医学部３号館 3階S308
--
Mixed Modelを用いた研究で解析方法について丁寧に記述してある論文をお持ちの方は、是非ご紹介ください。
ご関心のある方は次回も是非いらしてください！（文責：宮本有紀）
また、次回は、Mixed Modelの論文を読みつつ、この春論文を提出された皆さんやご卒業される皆さんをねぎらい祝うことも同時に行いたいと思っています。
次回みんなで読んでみる予定にしているのは下記の論文です。
O'Campo P, Caughy MO, Nettles SM. Partner abuse or violence, parenting and neighborhood influences on children's behavioral problems. Soc Sci Med. 2010 May;70(9):1404-15. Epub 2010 Feb 12.
PMID: 20163906
日 時：2011年3月23日（水）
テーマ：Mixed Model (15)
場 所：医学部３号館 3階S308
--
Mixed Modelを用いた研究で解析方法について丁寧に記述してある論文をお持ちの方は、是非ご紹介ください。
ご関心のある方は次回も是非いらしてください！（文責：宮本有紀）
2011/02/09
次回Mixed Model続きの日程は調整中
次回Mixed Model続きは、諸事情により、当初予定していた開催日程を変更させていただきます。日程は決定し次第アップします。よろしくお願いいたします。
2011/01/28
次回Mixed Model続きの日程は調整中
以下の記事を「次回Mixed Model続きの日程は2月23日（水）」というタイトルで書いておりましたが、諸事情により、開催日程を延期させていただきたく御願い申し上げます。
また日程が決定し次第アップします。申し訳ありません。
ご不明な点は精神看護学分野の宮本<yyuki-tkyアットマークumin.ac.jp>までご連絡くださいませ。
前回ご参加の皆様、ありがとうございました！
とりあえず、次回もMixed Model続きです。
日 時：2011年2月23日（水）18:30～20:00 日程調整中
テーマ：Mixed Model (15)
場 所：医学部３号館 3階S308
次回みんなで読んでみる予定なのは下記の論文です。
O'Campo P, Caughy MO, Nettles SM. Partner abuse or violence, parenting and neighborhood influences on children's behavioral problems. Soc Sci Med. 2010 May;70(9):1404-15. Epub 2010 Feb 12.
PMID: 20163906
Growth Modelは継続的にデータを取っていく（時間を連続して）ときなんかによくて、Mixed Modelは時間をカテゴリみたいに入れたりするときに使ったり。ということや、
SPSSのマニュアルとか、読んでみると良いこととか、
ICCが0.05-0.2の間だとMulti-levelは適しているようだ、とか、
Multi-level分析の結果は、標準化された値が出るわけではないので、変数間の影響力の程度の違いを見るのであれば変数の値を標準化させてから投入すると良いけれども、どちらかというと、この変数の値がこれくらいあがると従属変数にこのような影響が、というのを見ることが多いのではないか、とか
メモがいっぱいで、整理していると時間がかかっちゃうので、ばらばらしたまま書いておきます。
前回読んだのは下記です。
Cho SH, Ketefian S, Barkauskas VH, Smith DG. The effects of nurse staffing on adverse events, morbidity, mortality, and medical costs. Nurs Res. 2003 Mar-Apr;52(2):71-9.
PMID: 12657982
ほんとにばらばらとしたメモですが。。。
Mixed Modelを用いた研究で解析方法について丁寧に記述してある論文をお持ちの方は、是非ご紹介ください。
ご関心のある方は次回も是非いらしてください！（文責：宮本有紀）
また日程が決定し次第アップします。申し訳ありません。
ご不明な点は精神看護学分野の宮本<yyuki-tkyアットマークumin.ac.jp>までご連絡くださいませ。
前回ご参加の皆様、ありがとうございました！
とりあえず、次回もMixed Model続きです。
日 時：2011年
テーマ：Mixed Model (15)
場 所：医学部３号館 3階S308
次回みんなで読んでみる予定なのは下記の論文です。
O'Campo P, Caughy MO, Nettles SM. Partner abuse or violence, parenting and neighborhood influences on children's behavioral problems. Soc Sci Med. 2010 May;70(9):1404-15. Epub 2010 Feb 12.
PMID: 20163906
Growth Modelは継続的にデータを取っていく（時間を連続して）ときなんかによくて、Mixed Modelは時間をカテゴリみたいに入れたりするときに使ったり。ということや、
SPSSのマニュアルとか、読んでみると良いこととか、
ICCが0.05-0.2の間だとMulti-levelは適しているようだ、とか、
Multi-level分析の結果は、標準化された値が出るわけではないので、変数間の影響力の程度の違いを見るのであれば変数の値を標準化させてから投入すると良いけれども、どちらかというと、この変数の値がこれくらいあがると従属変数にこのような影響が、というのを見ることが多いのではないか、とか
メモがいっぱいで、整理していると時間がかかっちゃうので、ばらばらしたまま書いておきます。
前回読んだのは下記です。
Cho SH, Ketefian S, Barkauskas VH, Smith DG. The effects of nurse staffing on adverse events, morbidity, mortality, and medical costs. Nurs Res. 2003 Mar-Apr;52(2):71-9.
PMID: 12657982
ほんとにばらばらとしたメモですが。。。
Mixed Modelを用いた研究で解析方法について丁寧に記述してある論文をお持ちの方は、是非ご紹介ください。
ご関心のある方は次回も是非いらしてください！（文責：宮本有紀）
登録:
投稿 (Atom)