glmmMLとlme4 ~ランダムなエフェcうと〜ってのだけど

いきなり関係ない人には全く意味わかんネタだけど

統計フリーソフトのRでGLMMを煮詰めるのに代表的なパッケージには

glmmML と lme4 てあるけど

この差はざっくりいえば

glmmML : random effectはclusterで指定 < 指定は１つのみだがIDとか数多くても計算してくれる

lme4 :説明変数に複数個指定できて、指定方法も入れ子状とかいろいろあり <便利だがIDとかの数が多い項目へのrandom effectは苦手で計算できない

何が違うのかなーてイマイチ理解が追いつかないんだけど

とりあえず比べてみると、今いじってるあるデータセットはほぼ同じ結果だった:

__________

> groupS <-read.delim(pipe("pbpaste"))
> groupS
clutch sTime small large
1 Z04 2218 sib sib
2 Z04 3281 sib sib
3 Z04 3510 sib sib
4 Z04 2547 sib sib
5 Z04 2047 sib sib
6 Z04 2114 sib sib
7 Z04 3053 sib sib
8 Z04 2730 sib sib
9 Z04 2628 sib sib
10 Z04 2603 sib sib
11 Z04 2404 sib sib
12 Z04 1684 sib sib
13 Z04 2738 sib sib
14 Z04 2514 sib sib
15 Z04 2623 sib sib
16 Z04 2708 sib sib
17 Z04 2283 sib sib
18 Z04 2074 sib sib
19 Z11 3979 sib sib
20 Z11 2723 sib sib
21 Z11 2603 sib sib
22 Z11 2552 sib sib
23 Z11 1412 sib sib
24 Z11 3858 sib sib
25 V2 2249 sib sib
26 V2 2822 sib sib
27 V2 2459 sib sib
28 V2 1980 sib sib
29 V2 3064 sib sib
30 V2 2338 sib sib
31 F1 2332 sib nonsib
32 F1 2644 sib nonsib
33 F1 2607 sib nonsib
34 F1 2943 sib nonsib
35 F1 2495 sib nonsib
36 F1 2280 sib nonsib
37 Z02 2207 sib nonsib
38 Z02 2964 sib nonsib
39 Z02 1897 sib nonsib
40 Z02 2738 sib nonsib
41 Z02 2619 sib nonsib
42 Z02 2883 sib nonsib
43 Z11 3410 sib nonsib
44 Z11 2934 sib nonsib
45 Z11 2851 sib nonsib
46 Z11 2709 sib nonsib
47 Z11 2677 sib nonsib
48 V2 2480 sib nonsib
49 V2 2364 sib nonsib
50 V2 2381 sib nonsib
51 V2 2594 sib nonsib
52 V2 2933 sib nonsib
53 V2 3296 sib nonsib
54 Z11 2950 sib nonsib
55 Z11 2190 sib nonsib
56 Z11 2420 sib nonsib
57 V2 2522 sib nonsib
58 Z11 2347 sib nonsib
59 F1 3490 nonsib sib
60 F1 1954 nonsib sib
61 F1 1465 nonsib sib
62 F1 2250 nonsib sib
63 F1 2414 nonsib sib
64 F1 2397 nonsib sib
65 Z11 2746 nonsib sib
66 Z11 2623 nonsib sib
67 Z11 707 nonsib sib
68 V2 1697 nonsib sib
69 V2 1587 nonsib sib
70 V2 1124 nonsib sib
71 V2 1036 nonsib sib
72 V1 3424 nonsib sib
73 V1 1203 nonsib sib
74 V1 2058 nonsib sib
75 V1 2601 nonsib sib
76 V1 3463 nonsib sib
77 V1 3298 nonsib sib
78 Z11 3166 nonsib sib
79 Z11 3220 nonsib sib
80 V2 3060 nonsib sib
81 V2 2377 nonsib sib
82 V2 2188 nonsib sib
83 F1 2894 nonsib nonsib
84 F1 2260 nonsib nonsib
85 F1 2229 nonsib nonsib
86 F1 3056 nonsib nonsib
87 F1 2730 nonsib nonsib
88 F1 2764 nonsib nonsib
89 Z02 1969 nonsib nonsib
90 Z02 4139 nonsib nonsib
91 Z02 1873 nonsib nonsib
92 Z02 2694 nonsib nonsib
93 Z02 2197 nonsib nonsib
94 V2 3555 nonsib nonsib
95 V2 3104 nonsib nonsib
96 V2 1959 nonsib nonsib
97 V1 3045 nonsib nonsib
98 V1 2534 nonsib nonsib
99 V1 2973 nonsib nonsib
100 V1 2674 nonsib nonsib
101 Z11 1968 nonsib nonsib
102 Z11 2009 nonsib nonsib
103 Z11 1901 nonsib nonsib
104 V1 3377 nonsib nonsib
105 V2 3321 nonsib nonsib
106 V2 3363 nonsib nonsib

> library(glmmML)

> glmm1 <-glmmML(cbind(sTime, 4800 - sTime) ~ small*large, cluster = clutch, family= binomial, data = groupS)

> summary(glmm1)

Call: glmmML(formula = cbind(sTime, 4800 - sTime) ~ small * large, family = binomial, data = groupS, cluster = clutch)

coef se(coef) z Pr(>|z|)
(Intercept) 0.21926 0.057657 3.803 0.000143
smallsib 0.00886 0.008426 1.052 0.293000
largesib -0.34280 0.008658 -39.595 0.000000
smallsib:largesib 0.36812 0.013220 27.845 0.000000

Scale parameter in mixing distribution: 0.1404 gaussian
Std. Error: 0.04072

LR p-value for H_0: sigma = 0: 7.456e-286

Residual deviance: 31250 on 101 degrees of freedom AIC: 31260

> library(lme4)
> library(lmerTest)

> glmm2 = glmer(cbind(sTime, 4800 - sTime) ~ small*large + (1| clutch),  family= binomial, data = groupS)
> summary(glmm2)
Generalized linear mixed model fit by maximum likelihood
(Laplace Approximation) [glmerMod]
Family: binomial ( logit )
Formula:
cbind(sTime, 4800 - sTime) ~ small * large + (1 | clutch)
Data: groupS

AIC BIC logLik deviance df.resid
32202.5 32215.8 -16096.3 32192.5 101

Scaled residuals:
Min 1Q Median 3Q Max
-45.221 -10.608 0.085 8.947 45.782

Random effects:
Groups Name Variance Std.Dev.
clutch (Intercept) 0.01971 0.1404
Number of obs: 106, groups: clutch, 6

Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.219271 0.057659 3.80 0.000143
smallsib 0.008861 0.008426 1.05 0.292954
largesib -0.342805 0.008658 -39.60 < 2e-16
smallsib:largesib 0.368117 0.013220 27.85 < 2e-16

(Intercept) ***
smallsib
largesib ***
smallsib:largesib ***
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
(Intr) smllsb largsb
smallsib -0.077
largesib -0.073 0.483
smllsb:lrgs 0.029 -0.563 -0.608

____

ちなみに

傾き固定でclutchごとにランダム切片(1| clutch)から、切片固定でclutchごとにランダム傾き(0| clutch)にすると、しぬ。あまり意味がないだろうが。

____

> glmm3 = glmer(cbind(sTime, 4800 - sTime) ~ small*large + (0| clutch),  family= binomial, data = groupS)
(function (cl, name, valueClass) でエラー:
assignment of an object of class “numeric” is not valid for @‘Dim’ in an object of class “dgTMatrix”; is(value, "integer") is not TRUE

__________

はにゃ〜。。。∧( 'Θ' )∧

この２つのパッケージって、結局どんな具合でrandom effectの入りかたが変わるんでしょな？

まーcorrelationも算出してくれるから、lme4とlmeTestが便利っぽいですの。

統計に関してはコレまでも他にいーーっっぱいうーむううぅーむ...となってきてるので。

またログを整理する気になったらブログに書いておこ。忘備録として。