LMEモデルフィッティング R言語編

続き。

LMEモデル推定をR言語を用いて実践してみる。

まずRをインストールし、RGuiを開く。
今回使うパッケージはlme4とlmeTestなので、これらを読み込む。lme4のlmer()関数によってLMEモデル推定を行うが、lme4だけだと固定効果のp値が得られないため、lmeTestを追加で読み込む必要がある。

>Install.packages(“lme4”) #lme4インストール
>library(“lme4”)　#lme4読み込み
>options(repos="http://cran.ism.ac.jp/")
>install.packages("lmerTest")　#lmerTestインストール
>library(“lmerTest”) #lmerTest読み込み

次にデータを読み込んでみる。今回、4人の学生の勉強日数と点数のテストデータをCSVファイルで用意した。

> s_data <- read.csv("score_data.csv")
> s_data
   subject days    Score
1        0   10 97.88387
2        0    8 98.59650
3        0    6 97.65732
4        0    4 96.16080
5        0   10 99.59362
6        0    8 97.77845
7        0   10 98.22498
8        0   10 98.04627
9        0    8 98.25393
10       0   10 97.46261
11       0    6 98.18022
12       0    4 96.99735
13       0   13 97.90047
14       0    6 98.66602
15       0   15 98.52445
16       0   12 97.65170
17       0    8 97.65037
18       0   14 98.28197
19       0   11 97.47814
20       0   10 98.51969
21       0   11 98.76532
22       0    8 97.27131
23       0    6 98.02530
24       0   10 98.52534
25       0   14 99.38336
26       0    4 94.42770
27       0    7 98.01417
28       0   12 98.10540
29       0    4 97.87208
30       1   10 97.82623
31       1   10 96.94224
32       1   14 98.31716
33       1   13 95.45953
34       1   18 96.94513
35       1   11 93.20422
36       1    6 94.51405
37       1    8 91.38044
38       1   16 93.66502
39       1    6 93.81745
40       1   13 97.06315
41       1    8 96.63337
42       1   11 96.66075
43       1    9 96.37080
44       1   13 98.71230
45       1   18 97.57873
46       1    5 94.33944
47       1    7 94.78274
48       1   17 97.75903
49       1   11 96.82348
50       1   17 97.77848
51       1    5 91.73288
52       2   15 90.96240
53       2   12 91.43013
54       2   15 96.27122
55       2   14 93.39732
56       2   18 92.99953
57       2   11 93.39127
58       2    9 91.33970
59       2    5 87.46454
60       2   18 91.80603
61       2    4 90.71793
62       2   10 91.20598
63       2    4 90.24194
64       3    5 90.67530
65       3   17 95.64061
66       3   14 97.77640
67       3   15 96.77373
68       3   14 95.39185
69       3   14 97.48028
70       3    5 91.10883
71       3    4 88.89533
72       3   10 95.65205
73       3    5 92.51790
74       3    6 89.78949
75       3   14 94.84116
76       3   15 96.99177
77       3   13 95.21246

グラフにplotすると

>plot_data <- function(x){
for (i in 1:nrow(x)){
plot(x[i,2],x[i,3],col=x[i,1]+1,xlim=c(min(x[,2]),max(x[,2])),ylim=c(min(x[,3]),max(x[,3])),ann=F)
par(new=T)
}
}
>plot_data(s_data)

各学生の勉強時間（横軸）とテストの点数（縦軸）

さて、ここからが本題であるが、lmer関数でLMEモデル推定を行う。

第一引数にモデル関数、第二引数に解析データの行列を入れる。

 モデル関数の表記は以下のように行う。Yは従属変数ラベル、Xi(i=1,2,…)は説明変数ラベル（切片は省略するが、切片のみを記述する場合は１を入れる）、Ci(i=1,2,…)はクラスタラベルである

Y~X1+X2+X3…+(X1'+X2'+…|C1)+(X1"+X2"+…|C2)+…

例えば、上記のデータで、勉強日数に応じて点数が変化するようなモデルを考える。
勉強日数の係数（傾き）が学生共通で一定（固定効果）で、切片は学生それぞれに対して変わる（変量効果）ような場合、

Score~days+(1|subject)

となる。従って、

lmer(Score~days+(1|subject),data=s_data)

で推定が行える。結果を見るには、

> test_model1 <- lmer(Score~days+(1|subject),data=s_data)
> summary(test_model1)

とすれば見ることができる。
推定結果の見方を説明する。

Linear mixed model fit by REML t-tests use Satterthwaite approximations to
  degrees of freedom [merModLmerTest]
Formula: Score ~ days + (1 | subject)
   Data: s_data

REML criterion at convergence: 297.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.4714 -0.6311  0.1049  0.7079  2.1640 

Random effects:
 Groups   Name        Variance Std.Dev.
 subject  (Intercept) 8.186    2.861   
 Residual             2.232    1.494   
Number of obs: 77, groups:  subject, 4

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)    
(Intercept) 91.59004    1.51029  3.61000  60.644 1.48e-06 ***
days         0.31828    0.04266 72.11000   7.461 1.53e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
     (Intr)
days -0.298

赤文字にした部分が変量効果（Random effects）、固定効果（Fixed effects）についての記述。
変量効果については、切片（Intercept）の変量効果における正規分布の分散(Variance)と標準偏差(Std.Dev)が記してある。
固定効果については、dayの係数と切片について、推定した値（Estimate）とP値（Pr）が記してある。
固定効果についてはfixef()で、それぞれのクラスタにおける変量効果についてはranef()で参照できる。

> fixef(test_model1)
(Intercept)        days 
  91.590043    0.318283 
> ranef(test_model1)
$subject
  (Intercept)
0   3.4623447
1   0.6747299
2  -3.3261432
3  -0.8109314

傾きにも変量効果を加えて再度推定してみる。

> test_model2 <- lmer(Score~days+(days|subject),data=s_data)

plotしてみる

>lmer_plot <- function(x,y){
rmat <- ranef(y)
fmat <- fixef(y)
for (i in 0:max(x[,1])+1){
func_buff <- function(x){(fmat[2]+rmat[[1]][i,2])*x+(fmat[1]+rmat[[1]][i,1])}
plot(func_buff,xlim=c(min(x[,2]),max(x[,2])),ylim=c(min(x[,3]),max(x[,3])),col=i,ann=F)
par(new=T)
}
fix_func <- function(x){fmat[2]*x+fmat[1]}
plot(fix_func,xlim=c(min(x[,2]),max(x[,2])),ylim=c(min(x[,3]),max(x[,3])),lty=2,ann=F)
par(new=T)
}

>lmer_plot(s_data,test_model2)

点線：固定効果モデル
各実線：それぞれの学生の混合効果モデル

また、切片、傾きそれぞれの変量効果の確率分布をplotすると次のようになる。

切片（定数項の係数）の確率分布

傾き（daysの係数）の確率分布

【Ubuntu12.04】Latexのパッケージ（footmisc）導入

そもそもlatexのパッケージ導入自体初めてのような。。。

Ubuntu12.04では/etc/texmf/texmf.confに設定ファイルがあるので、設定ファイルを開いてパスを確認。

TEXMFHOME = $HOME/texmf

ローカルパッケージを入れたことがなければディレクトリ自体ないと思うので、以下のコマンドでディレクトリを作成する

$ mkdir -p ~/texmf/tex/latex

作成したディレクトリにスタイルファイル（～.sty）を入れて、

$ mktexlsr ~/texmf

で導入完了。texファイルに￥usepackage{～} で使える。

---

~footmiscの場合～

以下のリンクのDownloadからzipをダウンロードして適当なところで展開する。
http://www.ctan.org/tex-archive/macros/latex/contrib/footmisc

展開したディレクトリ内で

$ platex footmisc.ins

 を実行するとfootmisc.styができるので、これを先ほどの~/texmf/tex/latexに移動してmktexlsrすればよい。

ただ、いざfootnoteを使おうとすると、

！Package footmisc Error: Can't define commands for footnote symbol.

と起こられてしまう。Hでhelpを見ると

Use e-LaTeX, or load package fixltx2e before footmisc

 どうやらfixltx2eというパッケージを入れなければならないようなので、以下のリンクのSourcesからfixltx2e.styをダウンロードする。
http://www.ctan.org/pkg/fixltx2e

同様に~/texmf/tex/latexに移動、mktexlsrする。

texファイルに￥usepackage{fixltx2e}を追記すればfootnoteが使えるようになる。

統計初心者でも2分で雰囲気理解できるLME（線形混合）モデル

この文章の主旨

• 統計初心者が統計の雰囲気を理解する
• 統計初心者がLMEモデルの雰囲気を理解する

※素人なので内容の正しさについては保証しかねます
※統計を全く触れたことない人は、〔青文字〕は読み飛ばしてください

---

Q&A

Q.統計ってなあに？
A.「おっぱいの大きい子は優しいってホント？」「二郎系ラーメンの致死量は？」そういう疑問を実際のデータから解決することです

Q.LMEモデルとは？
A.統計の本などにもあまり載っていなく、統計モデルの中でも2番目くらいに扱いが難しいと言われているモデルです〔1番は階層ベイズモデル〕。ただ、雰囲気だけ（一変数一次）ならなんとか理解できるはず。。。。

---

簡単のため、一変数一次の場合について考えてみる。
学生の勉強時間とテストの点数を例にとる。

あるAさんについて、複数のテストを行わせ、その勉強時間と点数に関してプロットしたところ、次のようになった。

このとき、〔線形モデルで〕

点数=a×勉強時間+b+ε

aは傾き　bは切片　εは誤差

と表現でき、「Aさんは勉強すればするほど高い点数がとれる」ということが分かる。
また、ここから色々頑張ると、この結果が偶然起こり得るかどうかを確認することもできる。

このように、
①モデル（どのように表現するか）
②有意かどうか（偶然起こったのではないか）
を調べたい。

次に、同様に14名の学生に関して、勉強時間と点数をプロットすると、次のようになった。

このとき、学生によって傾き（点数の上がり方）や、切片（勉強をしなくても取れる点数）は異なるものの、全ての学生が勉強するほど点数が上がっており、「すべての学生が勉強時間に比例して点数が上がる」ことは言えそうである。

このまま学生という括りを無視し、全てのデータを統合してしまうと、次の図のようになり、
勉強時間と点数の間に、関係は確認されなくなってしまう。〔回帰式は勿論、相関係数も有意にはならない〕

これでは「すべての学生が勉強時間に比例して点数が上がる」ことは言えない。
学生によって傾きや切片が異なる要因が明確で、測れるものであれば、その要因を変数にとることで、モデル化できる〔重回帰分析等〕が、この場合これは難しい。
単に頭の回転が速いのかもしれないし、効率の良い勉強方法を知っているのかもしれないし、親の教え方が上手いのかもしれない。分からない。

また、この例と逆に、統合したデータでは傾向が見られても、実状全く意味を持たない結果であることもある。
〔有意な回帰式や相関係数が見られても、グループ内で有意でなければ意味がない。勿論、安直に正規化したり、相関係数の平均をとったりしてはいけない〕

このように、「学生Aさん」等の、グループ毎では全てある傾向を示しているものの、違うグループでは傾きや切片といったパラメータが異なり、かつパラメータが異なる要因が明確でない場合、LMEモデルを用いるのが適切である。

先の勉強時間と点数の例を、一変数一次のLMEモデルで表現すると、次のようになる。

点数=（a_f + a_r）×勉強時間+ (b_f + b_r)+ε

a_f,a_rは傾き　b_f,b_rは切片　εは誤差

a_fおよびb_fを固定効果（fixed effect）、a_rおよびb_rを変量効果（random effect）と呼ぶ。
固定効果はグループ間で共通しているパラメータ、変量効果はグループ内で共通しているパラメータである。
勉強時間と点数の例でいえば、変量効果が「個人差」を表している。

LMEモデルの利点は、「変量効果」を取り入れることで、個人差等の、実状がよくわからないような「ランダムな要素」も組み込んで表現できる点である。

LMEモデル解析はRのlmerを使うと簡単にできます。
いつか詳しく説明しますが、とりあえずサンプルコードを張っておきます。

---

実際にRを用いて混合モデル解析をやってみる。（サンプルデータの内容はあまり気にしないで）

$ install.packages("lme4")　#パッケージのインストール
$ library("lme4")  #インポート 

$ data1 <- read.csv("test.csv")　#サンプルデータ（CSVファイル）の読み込み
$ data1
   subject difficulty pupil

1        1        0.1  70.0
2        1        0.2  71.0
3        1        0.3  72.0
4        1        0.4  73.0
5        1        0.5  73.5
6        1        0.6  75.0
7        1        0.7  74.0
8        1        0.8  76.0
9        1        0.9  79.0
10       2        0.1  72.0
11       2        0.2  74.0
12       2        0.3  76.0
13       2        0.4  78.0
14       2        0.5  79.5
15       2        0.6  82.0
16       2        0.7  82.0
17       2        0.8  85.0
18       2        0.9  89.0
19       3        0.1  73.0
20       3        0.2  74.5
21       3        0.3  76.0
22       3        0.4  77.5
23       3        0.5  78.5
24       3        0.6  80.5
25       3        0.7  80.0
26       3        0.8  82.5
27       3        0.9  86.0

$ model1 <- lmer(pupil ~ difficulty + (difficulty|subject),data=datas)　#LMEモデル解析にかけて結果をmodel1へ
$ summary(model1)　#解析結果を見る（t値しか出ないので注意）

Linear mixed model fit by REML ['lmerMod']
Formula: pupil ~ difficulty + (difficulty | subject)
   Data: datas
REML criterion at convergence: 79.8
Scaled residuals:
     Min       1Q   Median       3Q      Max
-1.92414 -0.30170  0.05265  0.42106  1.92576
Random effects:
 Groups   Name        Variance Std.Dev. Corr
 subject  (Intercept)  1.1758  1.0843    
          difficulty  23.7130  4.8696   0.60
 Residual              0.7722  0.8788    
Number of obs: 27, groups:  subject, 3
Fixed effects:
            Estimate Std. Error t value
(Intercept)  70.1389     0.7265   96.55
difficulty   14.5000     2.8868    5.02
Correlation of Fixed Effects:
           (Intr)
difficulty 0.397 

---

参考
http://oku.edu.mie-u.ac.jp/~okumura/stat/meta.html
http://ja.wikipedia.org/wiki/%E6%B7%B7%E5%90%88%E3%83%A2%E3%83%87%E3%83%AB
http://www.lowtem.hokudai.ac.jp/plantecol/akihiro/obenkyou/GLMMexample.pdf　著者（隅田）
http://www.slideshare.net/yutamura1/ss-42303827?ref=http%3A%2F%2Fsssslide.com%2Fwww.slideshare.net%2Fyutamura1%2Fss-42303827

【windows7】パーティション縮小のためのデフラグと不要な移動負荷ファイルの削除

HDDからSSDに換装しようと思ったが、SSDに付いて来たSamsung Data Migrationがデフラグ中にエラーで中断するので、EASEUS TODO BACKUPで頑張ってみる

まず、HDDの状態確認のために以下のリンクからMydeflagをダウンロード＆インストール
http://www.mydefrag.com

「Analyze only」でHDDの状態を確認する。おそらく、論理クラスタ番号が大きい部分（可視化画面で上の方）に移動不可ファイルがたまっている（可視化画面のピンク線）。マウスポインタを合わせるとパスが確認できるので、移動不可ファイルの内容を確認しておく。

まずやるべきことをリストアップする
①デフラグソフトで移動不可で不要なファイル（system volume informationとかCSCとか）を削除する
②HDDの縮小作業のためにデフラグする
③元HHDのパーティションをターゲットSSDの容量より縮小させる

---

①デフラグソフトで移動可能なファイルを削除する
system volume informationは、自動復元ポイントです。削除すると、復元出来なくなるので注意。

• 以下を参考に不要な自動復元ポイントを削除する。

http://www.office-qa.com/win/win141.htm

• 確認のため、以下を参考にsystem volume informationに管理者のアクセス権限を追加

http://www.nzs.info/tip/Windows/Tips_Windows_Defrag_SystemVolumeInfo_20121012.html

• コマンドラインを管理者権限で開き、

$ cd C:¥"system volume information" 

復元ポイントファイル（{英数字の羅列}）削除を確認

次に、CSCを削除する。CSCはネットワーク共有フォルダのオフラインキャッシュである。

• コントロールパネル▶︎同期センター▶︎オフラインファイルの管理▶︎オフラインファイルを無効にする▶︎再起動　（これをしないとオフラインキャッシュにアクセスできない）
• 管理者権限でコマンドプロンプトを開き、
  $ takeown /r /f C:\Windows\CSC
  $ rd /s C:\Windows\CSC

その他、移動不可ファイルで不必要なファイルを随時削除する

---

②HDDの縮小作業のためにデフラグする
Mydeflagで「Consolidate free space」を実行する。
実行後も移動不可ファイルが上部に溜まっているようなら、システムファイルを移動出来る有料のデフラグソフトを使った方がいいかも。

---

③HDDのパーティションを縮小させる
スタート▶︎検索窓に「ハードディスク」▶︎ハードディスクパーティションの作成とフォーマット
Cドライブを選択し、右クリック▶︎パーティションの縮小
ターゲットSSDよりも容量を小さくする

---

これで大体900GBまでしか縮小できなかったパーティションが400弱まで縮小できた。
しかし本題のクローン作業でEASEUS TODO BACKUPでのクローンでは「セクターの読み込みに失敗しました」とエラーが。。。
どうやらいじっているうちに不整合が起こったみたいなので、CHKDSKを実行してみる（死ぬほど長い）

本題のクローン作業については次回