/
anovakun_482.txt
3036 lines (2761 loc) · 160 KB
/
anovakun_482.txt
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# 【ANOVA君:要因計画のタイプと水準数を入力することにより,分散分析を行う関数】
# 1)フリーの統計ソフトウェア「R」で動作する関数
# 2)被験者間要因(独立測度),被験者内要因(反復測度)のいずれか,または,両方を含む各タイプの分散分析を扱う
# 3)引数としては,最初にデータフレーム名,次に計画のタイプ(""で囲むこと)を入力し,その後,各要因の水準数を順に入力する
# (作成:井関龍太)
#
# 【使用法】
# anovakunに読み込むためのデータフレームは,以下のような形式で作っておく
# 1)被験者間要因はタテ,被験者内要因はヨコに並べる
# 2)被験者内要因を表す行はデータとして読み込まない(下の例では,点線より上の行は読み込まず,点線より下のデータのみ読み込む)
# 3)被験者間要因を表す列はデータとして読み込む(下の例では,a1などの文字を含む列;被験者間要因数が増えるたびにラベル用の列を増やす)
# 4)被験者間要因を表すラベルは,水準ごとに別の文字または数字を用いる(同じラベル=同じ水準と見なされる)
# 5)被験者内,被験者間要因ともに,前の方の要因から順に各水準のデータを入れ子状に整理して並べること(例を参照)
# 6)下の被験者1~6の列は例の説明のためにつけたものなので,実際のデータフレームには必要ない
#
# [AsBC計画の例](被験者間要因でデータ数が同じ)
# b1 b1 b2 b2
# c1 c2 c1 c2
# ------------------------------------
# a1 12 9 14 13 ---被験者1
# a1 13 10 14 12 ---被験者2
# a1 11 10 13 15 ---被験者3
# a2 18 12 16 15 ---被験者4
# a2 17 14 15 14 ---被験者5
# a2 15 13 18 15 ---被験者6
#
# データフレームを代入した変数名をxとすると,
#
# > anovakun(x, "AsBC", 2, 2, 2)
#
# のようにして関数を呼び出す
#
# [ABsC計画の例](被験者間要因でデータ数がふぞろい)
# c1 c2
# -------------------------------
# a1 b1 3.5 4.2 ---被験者1
# a1 b1 2.7 3.2 ---被験者2
# a1 b2 2.5 3.8 ---被験者3
# a1 b3 4.0 3.9 ---被験者4
# a2 b1 3.3 4.0 ---被験者5
# a2 b1 1.4 2.5 ---被験者6
# a2 b2 3.7 4.2 ---被験者7
# a2 b2 2.2 4.2 ---被験者8
# a2 b3 1.3 2.1 ---被験者9
# a2 b3 3.4 3.9 ---被験者10
#
# データフレームを代入した変数名をxとすると,
#
# > anovakun(x, "ABsC", 2, 3, 2)
#
# のようにして関数を呼び出す
#
# 【オプション】
# 1)long……long = Tとすると,ロング形式のデータを読み込んで処理する
# 2)type2……type2 = Tとすると,平方和の計算方法をタイプⅡに切り替える
# 3)nopost……nopost = Tとすると,下位検定を実行しない
# 4)tech……テクニカルアウトプット;tech = Tとすると,データフレームをリストでつないだ形式で結果を出力する
# 5)data.frame……data.frame = Tとすると,計算に使用したデータフレームを出力する(関数中でdatと表現されているデータフレーム)
# 6)copy……copy = Tとすると,出力結果をクリップボードにもコピーする(クリップボードの内容は上書きされる)
# 7)holm……holm = Tとすると,多重比較の方法がHolmの方法になる
# 8)hc……hc = Tとすると,多重比較の方法がHolland-Copenhaverの方法になる
# 9)s2r……s2r = Tとすると,Shafferの方法のための仮説数の計算法を具体的な棄却のパターンに基づく方法に変更する(Rasmussenのアルゴリズムに基づく)
# 10)s2d……s2d = Tとすると,Shafferの方法のための仮説数の計算法を具体的な棄却のパターンに基づく方法に変更する(Donoghueのアルゴリズムに基づく)
# 11)fs1……fs1 = Tとすると,ステップ1の基準をステップ2の基準に置き換えた方法でShafferの方法を行う
# 12)fs2r……fs2r = Tとすると,Shaffer2の方法とF-Shaffer1の方法を組み合わせた方法でShafferの方法を行う(Rasmussenのアルゴリズムに基づく)
# 13)fs2d……fs2d = Tとすると,Shaffer2の方法とF-Shaffer1の方法を組み合わせた方法でShafferの方法を行う(Donoghueのアルゴリズムに基づく)
# 14)hc, s2r……hc = Tかつs2r = T(または,fs2r = T,s2d = T,fs2d = T)とすると,Shaffer2の基準を用いてHolland-Copenhaverの方法を行う
# 15)holm, hc……holm = Tかつhc = Tとすると,Holmの調整基準にSidakの不等式を用いた多重比較(Holm-Sidak法)を行う
# 16)welch……welch = Tとすると,多重比較の際にKeselman-Keselman-Shafferの統計量とWelch-Satterthwaiteの近似自由度を用いる
# 17)criteria……criteria = Tとすると,多重比較の出力において,調整済みp値の代わりに調整済みの有意水準を表示する
# 18)lb……lb = Tとすると,すべての被験者内効果に対してイプシロンの下限値(Lower Bound)による自由度の調整を行う(Geisser-Greenhouseの保守的検定)
# 19)gg……gg = Tとすると,すべての被験者内効果に対してGreehnouse-Geisserのイプシロンによる自由度の調整を行う
# 20)hf……hf = Tとすると,すべての被験者内効果に対してHuynh-Feldtのイプシロンによる自由度の調整を行う
# 21)cm……cm = Tとすると,すべての被験者内効果に対してChi-Mullerのイプシロンによる自由度の調整を行う
# 22)auto……auto = Tとすると,球面性検定が有意であった被験者内効果に対してGreehnouse-Geisserのイプシロンによる自由度の調整を行う
# 23)mau……mau = Tとすると,球面性検定の方法がMauchlyの球面性検定になる
# 24)har……har = Tとすると,球面性検定の方法がHarrisの多標本球面性検定になる
# 25)iga……iga = Tとすると,改良版一般近似検定を行う;イプシロンの代わりに各種の推定値を算出し,分散分析の際に適用する
# 26)ciga……ciga = Tとすると,修正改良版一般近似検定を行う;イプシロンの代わりに各種の推定値を算出し,分散分析の際に適用する
# 27)eta……eta = Tとすると,分散分析表にイータ二乗を追加する
# 28)peta……peta = Tとすると,分散分析表に偏イータ二乗を追加する
# 29)geta……geta = Tとすると,分散分析表に一般化イータ二乗を追加する;geta = "要因ラベル(A, B, C...)"とすると,指定した要因を測定要因と見なして一般化イータ二乗を計算する;
# 複数の要因を測定要因として指定するには,例えば,geta = c("A", "C")のようにする(計画に含まれない要因の指定は無効になる)
# 30)eps……eps = Tとすると,分散分析表にイプシロン二乗を追加する
# 31)peps……peps = Tとすると,分散分析表に偏イプシロン二乗を追加する
# 32)geps……geps = Tとすると,分散分析表に一般化イプシロン二乗を追加する;geps = "要因ラベル(A, B, C...)"とすると,指定した要因を測定要因と見なして一般化イプシロン二乗を計算する;
# 複数の要因を測定要因として指定するには,例えば,geps = c("A", "C")のようにする(計画に含まれない要因の指定は無効になる)
# 33)omega……omega = Tとすると,分散分析表にオメガ二乗(加算モデル)を追加する
# 34)omegana……omegana = Tとすると,分散分析表にオメガ二乗(非加算モデル)を追加する
# 35)pomega……pomega = Tとすると,分散分析表に偏オメガ二乗を追加する
# 36)gomega……gomega = Tとすると,分散分析表に一般化オメガ二乗(加算モデル)を追加する;gomega = "要因ラベル(A, B, C...)"とすると,指定した要因を測定要因と見なして一般化オメガ二乗を計算する;
# 複数の要因を測定要因として指定するには,例えば,gomega = c("A", "C")のようにする(計画に含まれない要因の指定は無効になる)
# 37)gomegana……gomegana = Tとすると,分散分析表に一般化オメガ二乗(非加算モデル)を追加する;gomegana = "要因ラベル(A, B, C...)"とすると,指定した要因を測定要因と見なして一般化オメガ二乗を計算する;
# 複数の要因を測定要因として指定するには,例えば,gomegana = c("A", "C")のようにする(計画に含まれない要因の指定は無効になる)
# 38)prep……prep = Tとすると,分散分析表にp_repを追加する
# 39)nesci……nesci = Tとすると,出力を指定している効果量について非心F分布に基づく信頼区間を算出する;このオプション単独では機能しないので注意;反復測定要因を含む効果と非加算モデルに基づく効果量については信頼区間を計算しない
# 40)besci……besci = Tとすると,出力を指定している効果量についてブートストラップに基づく信頼区間を算出する;このオプション単独では機能しないので注意
# 41)cilmd……cilmd = Tとすると,記述統計量の表に差分調整型のLoftus-Massonの信頼区間を追加する
# 42)cilm……cilm = Tとすると,記述統計量の表にLoftus-Massonの信頼区間を追加する
# 43)cind……cind = Tとすると,記述統計量の表に差分調整型の正規化に基づく信頼区間を追加する
# 44)cin……cin = Tとすると,記述統計量の表に正規化に基づく信頼区間を追加する
# 45)ciml……ciml = Tとすると,記述統計量の表にマルチレベルモデルに基づく信頼区間を追加する(lmerTestパッケージが必要)
# 46)cipaird……cipaird = Tとすると,差分調整型のペアワイズ信頼区間を出力する
# 47)cipair……cipair = Tとすると,ペアワイズ信頼区間を出力する
# 48)bgraph……bgraph = "信頼区間のオプション名"とすると,信頼区間つきの棒グラフを出力する;信頼区間は内側・外側の順に2つまで指定できる;
# 例えば,bgraph = c("cind", "ciml")とすると,差分調整型の正規化に基づく信頼区間とマルチレベルモデルに基づく信頼区間を描画する;
# なお,このオプションは3要因までの計画についてのみ機能する
#
# [オプション使用の例](テクニカルアウトプットによる出力とHolmの方法による多重比較を指定)
#
# > anovakun(x, "AsB", 2, 2, tech = T, holm = T)
#
# 【技術情報】
# 1)anovakunを構成する関数は,仕様上は,最大で26要因までの計画に対応できる;この上限は,要因を表すラベルとしてアルファベット26文字(LETTERSとletters)を使用していることによる
# 2)ci.calc関数によるマルチレベルモデルに基づく信頼区間の計算には,lmerTestパッケージを使用している;自由度の推定にはlsmeans関数を用いている(Kenward-Roger法)
# 3)ss.calc関数は,デフォルトでは,タイプⅢ平方和の計算法に基づいて分散分析を行う
# 4)epsilon.calc関数は,デフォルトでは,被験者内要因を含むデータに対してMendozaの多標本球面性検定を行う(近似カイ二乗による)
# 5)epsilon.calc関数は,オプション指定により,被験者内要因を含むデータに対してMauchlyの球面性検定を行う(近似カイ二乗による)
# 6)epsilon.calc関数は,オプション指定により,被験者内要因を含むデータに対してHarrisの多標本球面性検定を行う(近似カイ二乗による)
# 7)epsilon.calc関数によるHuynh-Feldtのイプシロンの計算法は,Lecoutre(1991)の修正に基づく
# 8)epsilon.calc関数によるIGAとCIGAは,非加重平均を想定した計算法に基づく(Algina, 1997; Algina & Oshima, 1995)
# 9)anova.modeler関数は,IGAとCIGAを用いた際には,multiplier(b_hatとc_hat)によって調整した後のF値をF値の列に表示する
# 10)anova.modeler関数による加算モデルに基づくオメガ二乗,一般化オメガ二乗の計算式は,被験者と各要因の間のすべての交互作用が存在しないことを仮定するモデルによる(Dodd & Schultz, 1973を参照)
# 11)mod.Bon関数は,デフォルトでは,Shafferの方法による多重比較を行う(任意の棄却パターンにおける可能な真の帰無仮説の最大数に基づく方法)
# 12mod.Bon関数におけるShafferの方法,Holland-Copenhaverの方法の有意水準の計算は,Rasmussen(1993),Donoghue(2004)のアルゴリズムに基づく(オプションによる)
# 13)mod.Bon関数による多重比較では,p値の低い対から順に基準値よりも値が低いか否かを判定し,いったん有意でない対が見つかったら,
# 以降の対についてはすべて判断を保留する(p値が基準値を下回っても*マークを表示しない);調整済みp値の表示の際には,調整済みp値が
# 既定の有意水準(5%)を下回った対にのみ*マークを表示する
# 14)pro.fraction関数は,単純主効果の検定において誤差項をプールしない(水準別誤差項を使用;サブセットに分散分析を再適用するのと同じ)
#
# 【このファイルの含む関数】
# hmean……調和平均を計算する
# read.clip……クリップボードの情報を読み込む
# anovakun……プロセス全体の制御を行う
# uni.long……データフレームをロング形式に変形する
# ci.calc……平均の信頼区間を計算する
# ci.bars……信頼区間つきの棒グラフを作る
# elematch……文字列中のすべての要素を含む文字列をマッチングする
# expand.gmatrix……ベクトルの組み合わせを作る
# sig.sign……有意水準に合わせて記号を表示する
# epsilon.calc……Greenhouse-GeisserとHuynh-Feldtのイプシロンを計算する
# ss.calc……平方和を計算する
# qlambda.ncf……非心F分布のパラメータの信頼限界を算出する
# anova.modeler……分散分析表を作る
# mpginv……ムーア・ペンローズの逆行列を計算する
# wj.calc……Welch-Jamesアプローチに基づく近似検定を行う
# mod.Bon……修正Bonferroniの方法による多重比較を行う
# post.analyses……下位検定を行う
# pro.fraction……効果のタイプに適した下位検定を割り当てる
# boot.esci……ブートストラップ法に基づいて効果量の信頼区間を計算する
# boot.anova……主分析と単純主効果の検定の効果量を一括して計算する
# boot.inter……単純主効果の検定の効果量のみを計算する
# anova.output……出力の種類ごとに設定を割り当てる
# table.out……データフレームに書式を与えて出力する
# mod.Bon.out……修正Bonferroniの方法による多重比較の結果を出力する
# simeffects.out……単純主効果の検定の結果を出力する
# anovatan……指定した要因についてデータを分割して単純効果の検定を行う
#
# 【バージョン情報】
# 1)anovakun version 1.0.0(R 2.5.1作成;2007/09/03公開)
# ・分散分析,単純主効果の検定,多重比較
# 2)anovakun version 2.0.0(R 2.5.1作成;2007/10/01公開)
# ・球面性検定とイプシロン調整;epsilon.calc関数の追加とそれに伴う変更
# ・平方和の計算方法を修正;ss.calc関数の修正とelematch関数の追加,それに伴う変更
# ・多重比較を行う際に,他の要因の効果を考慮した上でのMSeを用いてt統計量を計算するように修正;mod.Bon関数と関連する部分の改修
# 3)anovakun version 2.1.0(R 2.5.1作成;2007/11/01公開)
# ・平方和の計算方法を変更し,高速化を試みる;ss.calc関数の変更;その他,最適化
# 4)anovakun version 2.2.0(R 2.5.1,R 2.6.0作成;2007/12/03公開)
# ・QR分解の方法をLAPACKに変更し,高速化を試みる;その他の修正
# 5)anovakun version 3.0.0(R 2.5.1,R 2.6.0作成;2008/01/04公開)
# ・タイプⅢ平方和を計算する機能を追加(タイプⅢをデフォルトに設定);ss.calc関数の変更とそれに伴う変更
# ・Shafferの方法のための可能な真の帰無仮説の数の計算アルゴリズムを追加;mod.Bon関数の変更とshaffer2,fshaffer1,fshaffer2オプションの追加
# ・多重比較の出力において調整済みp値を表示する機能を追加(デフォルトに設定);mod.Bon関数の変更とそれに伴う変更
# 6)anovakun version 3.1.0(R 2.5.1,R 2.6.0作成;2008/02/01公開)
# ・aov関数のアルゴリズムを利用して誤差平方和の計算の高速化を試みる;ss.calc関数の変更;その他の修正
# 7)anovakun version 3.2.0(R 2.5.1,R 2.6.0作成;2008/04/01公開)
# ・Shafferの方法のための別バージョンのアルゴリズムを追加;mod.Bon関数の変更とs2d,fs2dオプションの追加
# ・オプション名の変更;“shaffer2,fshaffer1,fshaffer2”を“s2r,fs1,fs2r”に変更
# ・効果量とp_repを計算する機能を追加;anova.modeler関数,anovatab.out関数ほかの変更
# 8)anovakun version 4.0.0(R 2.5.1,R 2.6.0,R 2.7.0作成;2008/06/02公開)
# ・Mendozaの多標本球面性検定とHarrisの多標本球面性検定を行う機能を追加(Mendozaをデフォルトに設定);epsilon.calc関数と関連する部分の変更
# ・Huynhの改良版一般近似検定とAlgina-Lecoutreの修正改良版一般近似検定を行う機能を追加;epsilon.calc関数及び関連する箇所の変更
# ・取り扱い可能な要因の数を拡張;anovakun関数,anova.modeler関数,epsilon.calc関数,pro.fraction関数ほかの変更
# ・複数の効果量をオプション指定したときに,すべての指標が同時に出力されるように変更;anova.modeler関数,pro.fraction関数ほかの変更
# 9)anovakun version 4.1.0(R 2.9.0,R 2.9.1,R 2.9.2作成;2009/09/01公開)
# ・欠損値があるケース(数値以外のデータを含む行)を除外して分析するように変更
# ・混合要因計画で被験者間要因のみを取り出して下位検定を行う際に警告メッセージが表示される点を修正
# 10)anovakun version 4.1.1(R 2.9.2,R 2.10.0,R 2.10.1作成;2010/01/04公開)
# ・Rのバージョン情報の表示を修正
# ・出力関数の表示エラーを修正;anovatab.out関数,simeffects.out関数の修正
# 11)anovakun version 4.2.0(R 2.10.0,R 2.10.1,R 2.11.1作成;2010/07/01公開)
# ・効果量の指標にイータ二乗,一般化イータ二乗,オメガ二乗,一般化オメガ二乗を追加;anova.modeler関数,anovatab.out関数ほかの変更
# ・データフレームの作成方法を変更;被験者間ラベルに文字列を使用しても指定した水準とデータがずれないようにする
# ・Lecoutre(1991)に基づいてHyunh-Feldtのイプシロンの計算式を変更;epsilon.calc関数の修正
# 12)anovakun version 4.3.0(R 2.14.1,R 2.14.2,R 2.15.0作成;2012/07/02公開)
# ・効果量の指標にイプシロン二乗,偏イプシロン二乗,一般化イプシロン二乗を追加;anova.modeler関数,pro.fraction関数ほかの変更
# ・オメガ二乗,偏オメガ二乗,一般化オメガ二乗の計算式を修正;anova.modeler関数,pro.fraction関数ほかの修正
# ・非加算モデルに基づくオメガ二乗,一般化オメガ二乗のオプションを追加;anova.modeler関数,pro.fraction関数ほかの修正
# ・結果をクリップボードにも出力するオプションcopyを追加(Linuxでは無効);anovakun関数の変更
# 13)anovakun version 4.3.1(R 2.15.0,R 2.15.1作成;2012/09/03公開)
# ・copy機能がMacで働かないエラーを修正;anovakun関数の修正
# 14)anovakun version 4.3.2(R 2.15.2作成;2013/01/04公開)
# ・read.clip関数を追加
# ・nopostオプションを追加;anovakun関数,post.analysis関数の変更
# ・copy機能がLinuxで働くように修正(xclipをインストールしている場合のみ機能);anovakun関数の変更
# ・出力の桁数に合わせて表の長さを調節するように変更;bstat.out関数,anovatab.out関数,mod.Bon.out関数,simeffects.out関数,post.out関数,each.out関数の変更
# 15)anovakun version 4.3.3(R 2.15.3,R 3.0.0作成;2013/05/07公開)
# ・R 3.0.0で動作するように修正;ss.calc関数の修正
# ・3要因以上の分析で一次の交互作用が有意だった際の多重比較の独立・非独立測度への割り当てが誤っていたのを修正;pro.fraction関数の修正
# ・出力桁数の調整;anovatab.out関数ほかの修正
# 16)anovakun version 4.4.0(R 3.0.1,R 3.0.2作成;2013/11/01公開)
# ・データフレームの変形手順を変更し,ロング形式のデータも扱えるように変更;anovakun関数の変更
# ・出力時の要因名,水準名を任意に指定できるように変更;anovakun関数,anova.modeler関数,post.analyses関数,pro.fraction関数ほかの変更
# ・出力表示の調整;bstat.out関数,anovatab.out関数,mod.Bon.out関数,simeffects.out関数,post.out関数,each.out関数の変更
# ・一般化効果量における複数の測定要因を指定する方法を変更(cまたはlist形式を使用する);anova.modeler関数,pro.fraction関数の変更
# 17)anovakun version 4.5.0(R 3.0.2作成;2014/02/03公開)
# ・信頼区間を計算する機能を追加;ci.calc関数の追加,anovakun関数,bstat.out関数の変更
# ・データフレーム変形部分をuni.long関数として分離;uni.long関数の追加,anovakun関数の変更
# ・3要因以上の計画について単純効果を扱うための関数を追加;anovatan関数の追加
# 18)anovakun version 4.5.1(R 3.0.2作成;2014/03/03公開)
# ・2つ以上の被験者間要因を含む計画での信頼区間算出の誤りを修正;ci.calc関数の修正
# ・信頼区間のオプションを追加;ci.calc関数の変更
# 19)anovakun version 4.6.0(R 3.0.2,R 3.1.0作成;2014/06/02公開)
# ・球面性検定の有意確率の算出を漸近展開を求める方式に変更;epsilon.calc関数の修正
# ・Chi-Mullerのイプシロンを計算し,自由度調整を行う機能を追加;epsilon.calc関数,anova.modeler関数ほかの変更
# ・サンプルサイズの多い反復測定要因を扱った場合の計算を高速化;ss.calc関数,epsilon.calc関数,anova.modeler関数ほかの変更
# 20)anovakun version 4.6.1(R 3.1.0作成;2014/07/01公開)
# ・4.6.0における変更に伴ってIGA,CIGA,anovatanを適用できなくなっていた点を修正;anova.modeler関数,anovatan関数の修正
# 21)anovakun version 4.6.2(R 3.1.0,R 3.1.1作成;2014/09/01公開)
# ・マルチレベルモデルに基づく信頼区間を算出するためのパッケージをlmerTestに変更;ci.calc関数の変更
# ・信頼区間つきの棒グラフを出力する機能の追加;ci.bars関数の追加
# ・多重比較の検定統計量としてWelch方式の統計量を計算する機能を追加;mpginv関数,wj.calc関数の追加,mod.Bon関数ほかの修正
# 22)anovakun version 4.7.0(R 3.1.0,R 3.1.1,R 3.1.2作成;2015/01/05公開)
# ・非心F分布に基づいて効果量の信頼区間を計算する機能を追加;qlamdba.ncf関数の追加;anova.modeler関数ほかの変更
# ・ブートストラップ法に基づいて効果量の信頼区間を計算する機能を追加;boot.esci関数,boot.anova関数,boot.inter関数の追加;その他関数の変更
# ・出力関数の変更;anova.output関数,table.out関数の追加;その他の関数の変更と削除
# 23)anovakun version 4.7.1(R 3.1.2,R 3.1.3作成;2015/04/01公開)
# ・3つ以上の被験者内要因を含むデザインでの自由度の割り当てのエラーを修正;episolon.calc関数の修正
# ・type2オプションが機能しなくなっていた点を修正;ss.calc関数の修正
# 24)anvoakun version 4.7.2(R 3.2.1,R 3.2.2作成;2015/10/01公開)
# ・単純主効果の検定の出力割り当てエラーほかの修正;pro.fraction関数の修正
# 25)anovakun version 4.8.0(R 3.2.2,R 3.2.3作成;2016/02/01公開)
# ・Huynhの改良版一般近似検定とAlgina-Lecoutreの修正改良版一般近似検定の修正;epsilon.calc関数ほかの修正
# 26)anovakun version 4.8.1(R 3.4.0,R 3.4.1作成;2017/08/01公開)
# ・long形式のデータを入力したときに起こるソートのエラーを修正;uni.long関数の修正
# 27)anovakun version 4.8.2(R 3.4.2作成;2018/0104公開)
# ・techオプションを指定したときにブートストラップ効果量が出力されないエラーを修正;anovakun関数の修正
#
# 【この関数の使用に関して】
# 1)anovakunとこれを構成する関数(コード群)は,自由に使用,改変,再配布していただいて結構です。
# 2)ただし,改変を加えたものを公開する際には,改変を加えたことを明記し,メインの関数の名前をanovakun以外の名前に変更してください。
# 3)anovakunとこれを構成する関数(コード群)の使用によって生じるいかなる結果に関しても作成者は責任を負いかねますのでご了承ください。
# 調和平均を計算する関数
hmean <- function(datvector){
return(length(datvector)/sum(1/datvector))
}
# クリップボードの情報を読み込む関数
# read.tableのラッパー関数;read.tableの出力先以外のオプションをすべて指定できる
# OSにかかわらず同じ書式で機能するようにしてある;ただし,LinuxはUbuntuでのみテストしており,xclipをインストールしていることを前提とする
read.clip <- function(...){
# OSごとのクリップボードを出力先に指定
plat.info <- .Platform
if(sum(grep("windows", plat.info)) != 0){# Windowsの場合
outboard <- "clipboard"
}else if(sum(grep("mac", plat.info)) != 0){# Macの場合
outboard <- pipe("pbpaste")
}else if(sum(grep("linux", R.version$system)) != 0){# Linuxの場合(xclipをインストールしている必要がある)
system("xclip -o | xclip -sel primary")
outboard <- "clipboard"
}else{# いずれのOSでもない場合
stop(message = "Unknown operating system!!")
}
# read.table関数を実行
read.table(file = outboard, ...)
}
# ANOVA君本体:プロセス全体の制御を行う関数
anovakun <- function(dataset, design, ..., long = FALSE, type2 = FALSE, nopost = FALSE, tech = FALSE, data.frame = FALSE, copy = FALSE,
holm = FALSE, hc = FALSE, s2r = FALSE, s2d = FALSE, fs1 = FALSE, fs2r = FALSE, fs2d = FALSE, welch = FALSE, criteria = FALSE,
lb = FALSE, gg = FALSE, hf = FALSE, cm = FALSE, auto = FALSE, mau = FALSE, har = FALSE, iga = FALSE, ciga = FALSE,
eta = FALSE, peta = FALSE, geta = NA, eps = FALSE, peps = FALSE, geps = NA, omega = FALSE, omegana = FALSE, pomega = FALSE,
gomega = NA, gomegana = NA, prep = FALSE, nesci = FALSE, besci = FALSE, cilmd = FALSE, cilm = FALSE, cind = FALSE, cin = FALSE,
ciml = FALSE, cipaird = FALSE, cipair = FALSE, bgraph = c(NA, NA)){
maxfact <- nchar(design) - 1# 実験計画全体における要因数
# データフレームの変形
datform <- uni.long(dataset = dataset, design = design, ... = ..., long = long)
dat <- datform$dat
factnames <- datform$factnames
flev <- datform$flev
miscase <- datform$miscase
# 記述統計量を計算する
if(sum(is.na(bgraph)) < 2){# 棒グラフに指定された信頼区間のオプションはオンにする
eval(parse(text = paste0(bgraph, " <- TRUE")))
}
baseresults <- ci.calc(dat = dat, design = design, factnames = factnames,
cilmd = cilmd, cilm = cilm, cind = cind, cin = cin, ciml = ciml, cipaird = cipaird, cipair = cipair)
# anova.modelerにデータフレームを送り,分散分析の結果を得る
mainresults <- anova.modeler(dat = dat, design = design, factnames = factnames, type2 = type2, lb = lb, gg = gg, hf = hf,
cm = cm, auto = auto, mau = mau, har = har, iga = iga, ciga = ciga, eta = eta, peta = peta, geta = geta,
eps = eps, peps = peps, geps = geps, omega = omega, omegana = omegana, pomega = pomega, gomega = gomega,
gomegana = gomegana, prep = prep, nesci = nesci)
# post.analysesにデータフレームと分散分析の結果を送り,下位検定の結果を得る
postresults <- post.analyses(dat = dat, design = design, factnames = factnames, mainresults = mainresults, type2 = type2,
nopost = nopost, holm = holm, hc = hc, s2r = s2r, s2d = s2d, fs1 = fs1, fs2r = fs2r, fs2d = fs2d, welch = welch,
criteria = criteria, lb = lb, gg = gg, hf = hf, cm = cm, auto = auto, mau = mau, har = har, iga = iga, ciga = ciga,
eta = eta, peta = peta, geta = geta, eps = eps, peps = peps, geps = geps, omega = omega, omegana = omegana, pomega = pomega,
gomega = gomega, gomegana = gomegana, prep = prep, nesci = nesci)
# 指定のあった効果量についてブートストラップ信頼区間を計算する
if(besci){
bootes <- boot.esci(dat, design, factnames = factnames, type2 = type2, nopost = nopost, mainresults = mainresults,
postresults = postresults, lb = lb, gg = gg, hf = hf, cm = cm, auto = auto, mau = mau, har = har, iga = iga, ciga = ciga,
eta = eta, peta = peta, geta = geta, eps = eps, peps = peps, geps = geps, omega = omega, omegana = omegana,
pomega = pomega, gomega = gomega, gomegana = gomegana, prep = prep, B = 2000, conf.level = 0.95)
mainresults$besci.info <- bootes$besci.info
mainresults$bescitab <- bootes$bescitab[[1]]
if(length(bootes$bescitab) > 1){# 単純主効果の結果も報告する場合
postnames <- names(postresults)
intnames <- postnames[nchar(postnames) == 3]# 単純主効果の検定を行ったリストのラベル
for(i in 1:length(intnames)){
postresults[[intnames[i]]]$bescitab <- bootes$bescitab[[i + 1]]
}
}
}
# 基本情報の取得
info1 <- paste0("[ ", design, "-Type Design ]")# 要因計画の型
info2 <- paste0("This output was generated by anovakun 4.8.2 under ", strsplit(R.version$version.string, " \\(")[[1]][1], ".")# バージョン情報など
info3 <- paste0("It was executed on ", date(), ".")# 実行日時
exe.info <- c(info1, info2, info3)
# Unbalancedデザイン(データ数ふぞろい)の場合,プロンプトを追加
if(length(unique(baseresults$bstatist$n)) != 1){
if(type2 == TRUE) mainresults$ano.info1 <- append(mainresults$ano.info1, c("== This data is UNBALANCED!! ==", "== Type II SS is applied. =="))
else mainresults$ano.info1 <- append(mainresults$ano.info1, c("== This data is UNBALANCED!! ==", "== Type III SS is applied. =="))
}
# 除外したケースの報告
if(miscase != 0){
baseresults$bstat.info1 <- append(baseresults$bstat.info1, paste0("== The number of removed case is ", miscase, ". =="))
}
# 結果を表示する
if(copy){# 指定があった場合,出力をクリップボードにコピー
plat.info <- .Platform
if(sum(grep("windows", plat.info)) != 0){# Windowsの場合
sink("clipboard", split = TRUE)
}else if(sum(grep("mac", plat.info)) != 0){# Macの場合
tclip <- pipe("pbcopy", "w")
sink(tclip, split = TRUE)
}else if(sum(grep("linux", R.version$system)) != 0){# Linxの場合(xclipをインストールしている必要がある)
tclip <- pipe("xclip -selection clipboard")
sink(tclip, split = TRUE)
}
}
if(tech){# データフレーム形式での出力の場合
postnames <- names(postresults)
intnames <- postnames[nchar(postnames) == 3]# 単純主効果の検定を行ったリストのラベル
if(length(intnames) > 0){
for(i in intnames){
postresults[[i]] <- postresults[[i]][-(7:9)]# sim.dmat,sim.cellN,sim.flevをカット
}
}
if(is.null(mainresults$besci.info)){
retlist <-list("INFORMATION" = rbind(info1, info2, info3),
"DESCRIPTIVE STATISTICS" = baseresults,
"SPHERICITY INDICES" = list(mainresults$epsi.info1, mainresults$epsitab),
"ANOVA TABLE" = list(mainresults$ano.info1, mainresults$anovatab, mainresults$nescitab),
"POST ANALYSES" = postresults)
}else{
retlist <-list("INFORMATION" = rbind(info1, info2, info3),
"DESCRIPTIVE STATISTICS" = baseresults,
"SPHERICITY INDICES" = list(mainresults$epsi.info1, mainresults$epsitab),
"ANOVA TABLE" = list(mainresults$ano.info1, mainresults$anovatab, mainresults$nescitab),
"EFFECT SIZE INFORMATION" = list(mainresults$besci.info, mainresults$bescitab),
"POST ANALYSES" = postresults)
}
if(data.frame == TRUE){
names(dat) <- c("s", factnames, "y")
retlist <- c(retlist, list("DATA.FRAME" = dat))# 計算に使用したデータフレームを付加
}
return(retlist)
}else{# 表形式での出力の場合
if(data.frame == TRUE){
names(dat) <- c("s", factnames, "y")
anova.output(maxfact = maxfact, exe.info = exe.info, baseresults = baseresults,
mainresults = mainresults, postresults = postresults)
return(list("DATA.FRAME" = dat))# 計算に使用したデータフレームを付加
}else{
anova.output(maxfact = maxfact, exe.info = exe.info, baseresults = baseresults,
mainresults = mainresults, postresults = postresults)
}
}
if(copy){
sink()
if(plat.info$OS.type != "windows"){# Mac,Linuxの場合
close(tclip)
}
}
# 指定があった場合には,棒グラフを出力
if(sum(is.na(bgraph)) < 2 & maxfact <= 3){
ci.bars(dat, design, factnames = factnames, inn.tier = bgraph[1], out.tier = bgraph[2])
}
}
# データフレームをロング形式に変形する関数
uni.long <- function(dataset, design, ..., long = FALSE){
maxfact <- nchar(design) - 1# 実験計画全体における要因数
spos <- match("s", strsplit(design, "")[[1]], nomatch = 0)
betlen <- spos - 1# 被験者間要因の数
withlen <- maxfact - spos + 1# 被験者内要因の数
levlist <- list(...)
# データフレームの変形
if(long){# ロング形式のデータフレームのヘッダをANOVA君に適合したものに変更
# 指定されたデザインとデータフレームの列数が合致しない場合にはメッセージを表示して終了
if((maxfact + 2) != ncol(dataset)){
stop(message = "\"anovakun\" has stopped working...\nThe entered design does not match the data.")
}
# ヘッダの差し替え
misv <- suppressWarnings(as.numeric(sapply(dataset[, ncol(dataset)], function(x) as.vector(x))))# 数値以外の値をNAに強制変換
dat <- cbind(data.frame(lapply(dataset[, 1:(maxfact+1), drop = FALSE], function(x) factor(x, levels = unique(x)))), misv)# ラベルをfactor型に変換
if(setequal(names(dataset), paste0("V", 1:ncol(dataset))) == TRUE | is.null(names(dataset)) == TRUE){# 入力時に指定されていない場合
factnames <- LETTERS[1:maxfact]# 要因名のラベルとしてアルファベットを使用
levlist <- lapply(1:maxfact, function(x) paste0(letters[x], 1:length(unique(dataset[, x+1]))))
for(i in 1:maxfact){# 被験者間要因の水準名を小文字アルファベットにする
levels(dat[, i+1]) <- levlist[[i]]# 水準名を上書き
}
}else{# 入力時に指定されている場合
factnames <- names(dataset)[-c(1, maxfact+2)]# もとのヘッダを保存;要因名のラベルとして使用
}
dat <- dat[do.call(order, dat[(ncol(dat)-1):1]), ]
names(dat) <- c("s", LETTERS[1:maxfact], "y")
flev <- sapply(dat, nlevels)[-c(1, maxfact+2)]# 各要因の水準数
attributes(flev) <- NULL# attributesを消去
# 欠損値の除去
if(!anyNA(misv)){# 欠損がなかった場合
miscase <- 0# 欠損ケースの数
}else{# 欠損があった場合
misid <- dataset[is.na(misv), 1]# 欠損のあったID
dat <- dat[rowSums(sapply(misid, function(x) x == dat[, 1])) == 0, ]# datからNAを含むIDの行(ケース)を除く
miscase <- length(misid)# 欠損ケースの数
}
}else{# ワイド形式のデータフレームをロング形式に変換
# 各要因の水準数と要因ラベルの設定
if(is.null(names(levlist))){# 入力時に指定されていない場合
flev <- unlist(levlist)# 各要因の水準数
factnames <- LETTERS[1:maxfact]# 要因名のラベル
levlist <- lapply(1:length(flev), function(x) paste0(letters[x], 1:flev[x]))
}else{# 入力時に指定されている場合
flev <- sapply(levlist, length)# 各要因の水準数
attributes(flev) <- NULL# attributesを消去
factnames <- names(levlist)# 要因名のラベル
}
# 水準数1が指定されていたらメッセージを表示して終了;存在しないオプション名が指定されると水準数1と解釈される
if(min(flev) == 1){
stop(message = "\"anovakun\" has stopped working...\nSome factor specifies only one level.\nOr maybe non-existent option-names are requested.")
}
# 指定されたデザインとデータフレームの列数が合致しない場合にはメッセージを表示して終了
if((betlen + ifelse(is.na(prod(flev[(betlen+1):length(flev)])), 1, prod(flev[(betlen+1):length(flev)]))) != ncol(dataset)){
stop(message = "\"anovakun\" has stopped working...\nThe entered design does not match the data.")
}
# 欠損値の除去
misv <- suppressWarnings(as.numeric(sapply(dataset[, (betlen+1):ncol(dataset)], function(x) as.vector(x))))# 数値以外の値をNAに強制変換
cdata <- array(misv, c(nrow(dataset), ncol(dataset) - betlen))# NAを含む行列
compcase <- complete.cases(cdata)# 完全ケース
dataset <- dataset[compcase, ]# datasetからNAを含む行(ケース)を除く
miscase <- sum(!compcase)# 欠損ケースの数
depv <- as.vector(cdata[compcase, ])# 従属変数をベクトル化
slab <- rep(paste0("s", 1:nrow(dataset)), ncol(dataset) - betlen)# 被験者のラベル
dat <- data.frame(factor(slab, levels = unique(slab)))# 被験者のラベルをデータフレームにする
# 被験者間要因を含む場合
if(betlen > 0){
betlab <- data.frame(lapply(dataset[, 1:betlen, drop = FALSE], function(x) factor(x, levels = unique(x))))# 被験者間要因のラベルをfactor型に変換;水準名をもとの順序にそろえる
for(i in 1:betlen){# 被験者間要因の水準名を小文字アルファベットにする
levels(betlab[, i]) <- levlist[[i]]# 水準名を上書き
}
betlab <- do.call("rbind", replicate(ncol(dataset) - betlen, betlab, simplify = FALSE))# くりかえし数のぶんだけ増やす
dat <- cbind(dat, betlab)# データフレームに追加
}
# 被験者内要因を含む場合
if(withlen > 0){
eachcue <- c(sapply(1:withlen, function(x) prod(flev[(betlen+x):maxfact])), 1)[-1]
timescue <- prod(flev[(betlen+1):maxfact]) / (eachcue * flev[(betlen+1):maxfact])
withlab <- data.frame(lapply(1:withlen, function(x) factor(rep(levlist[[betlen+x]], each = nrow(dataset) * eachcue[x],
times = timescue[x]), levels = levlist[[betlen+x]])))# 各被験者内要因の各水準を表すデータフレーム
dat <- cbind(dat, withlab)# データフレームに追加
}
dat <- cbind(dat, depv)# 従属変数をデータフレームに追加
names(dat) <- c("s", LETTERS[1:maxfact], "y")
}
return(list("dat" = dat, "factnames" = factnames, "flev" = flev, "miscase" = miscase))
}
# 平均の信頼区間を計算する関数
ci.calc <- function(dat, design, factnames = NA, conf.level = 0.95, cilmd = FALSE, cilm = FALSE, cind = FALSE, cin = FALSE, ciml = FALSE, cipaird = FALSE, cipair = FALSE){
maxfact <- nchar(design) - 1
spos <- match("s", strsplit(design, "")[[1]], nomatch = 0)
betlen <- spos - 1# 被験者間要因の数
withlen <- maxfact - spos + 1# 被験者内要因の数
cellN <- length(unique(dat$s))# 被験者間要因をつぶしての全被験者の数
flev <- sapply(2:(maxfact+1), function(x) nlevels(dat[, x]))# 各要因の水準数
comlev <- prod(flev[(betlen+1):maxfact])# すべての被験者内要因の組み合わせ水準数
bstat.info1 <- NULL
bstat.info2 <- NULL
# factnamesが省略されているときは名前をつける
if(anyNA(factnames)) factnames <- LETTERS[1:maxfact]
# 各条件ごとの平均と標準偏差を計算する
tabbase <- tapply(dat$y, dat[, (maxfact+1):2], function(x) x)
sncol <- sapply(tabbase, function(x) length(x))# セルごとのデータ数を計算
mncol <- sapply(tabbase, function(x) mean(x, na.rm = TRUE))# セルごとの平均を計算
sdcol <- sapply(tabbase, function(x) sd(x, na.rm = TRUE))# セルごとの標準偏差を計算
# 記述統計量の表において各要因の各水準を表すためのラベル(数字列)を作成
maincols <- expand.grid(lapply((maxfact+1):2, function(x) levels(dat[, x])))
maincols <- maincols[, order(maxfact:1)]# アルファベット順に列を並べ替え
# 記述統計量をデータフレームにまとめる
bstatist <- data.frame(maincols, sncol, mncol, sdcol)
names(bstatist) <- c(factnames, "n", "Mean", "S.D.")# 要因のラベルほかを列名に設定
# Loftus-Massonの信頼区間の計算
if(cilmd || cilm){
# プールした平方和と自由度の計算
pSS <- sum(dat$y^2) - sum(tapply(dat$y, interaction(dat[, 2:(maxfact+1)]), function(x) sum(x)^2/length(x)))
pDF <- nrow(dat) - prod(flev)
# 被験者内計画の場合は,被験者の誤差項の影響を除く
if(betlen == 0){
pSS <- pSS - (comlev * sum(tapply(dat$y, list(dat$s), function(x) mean(x)^2)) - sum(dat$y)^2/nrow(dat))
pDF <- pDF - (cellN - 1)
}
ci.tval <- qt((1 - conf.level)/2, pDF, lower.tail = FALSE) * sqrt(pSS / pDF * mean(1/sncol))
# 指定のあった指標を表に追加
if(cilmd){
bstat.info1 <- append(bstat.info1, "== Loftus-Masson's Difference-Adjusted Pooled Confidence Intervals ==")
bstatist <- cbind(bstatist, "CILMD-L" = mncol - ci.tval / sqrt(2), "CILMD-U" = mncol + ci.tval / sqrt(2))
}
if(cilm){
bstat.info1 <- append(bstat.info1, "== Loftus-Masson's Pooled Confidence Intervals ==")
bstatist <- cbind(bstatist, "CILM-L" = mncol - ci.tval, "CILM-U" = mncol + ci.tval)
}
}
# 正規化に基づく信頼区間の計算
if(cind || cin){
# 計画に合った計算法を適用
if(withlen == 0){# 被験者間計画の場合
nd.info <- "== Difference-Adjusted Confidence Intervals for Independent Means =="
n.info <- "== Confidence Intervals for Independent Means =="
vcom <- as.vector(tapply(dat$y, dat[, (betlen+1):2], function(x) sd(x)/sqrt(length(x))))
ci.tval <- qt((1 - conf.level)/2, sncol - 1, lower.tail = FALSE) * vcom
}else{# 反復測定要因を含む場合
nd.info <- "== Cousineau-Morey-Baguley's Difference-Adjusted Normalized Confidence Intervals =="
n.info <- "== Cousineau-Morey's Normalized Confidence Intervals =="
if(betlen == 0){# 被験者間要因がないとき(被験者内計画の場合)
caldat <- list(dat)
}else{# 1つ以上の被験者間要因があるとき(混合要因計画の場合)
caldat <- split(dat, dat[, (maxfact+1-withlen):2])
}
wdat <- lapply(caldat, function(x) data.frame(matrix(x$y[order(eval(parse(text = paste0("x[, ", c(1, (maxfact+1):(betlen+2)), "]"))))], length(x$y)/comlev, comlev)))# データをワイド形式のデータフレームに変換;tapplyによるアクセス順に対応させるための並べ替え
normdat <- lapply(wdat, function(x) x - matrix(rowMeans(x, na.rm = TRUE), nrow(x), 1) + mean(colMeans(x, na.rm = TRUE), na.rm = TRUE))# 個人ごとのデータによって正規化したデータ
nsd <- unlist(lapply(normdat, function(x) sapply(x, function(y) sd(y, na.rm = TRUE))), use.names = FALSE)# 正規化したデータを使って標準偏差を計算
ci.tval <- qt((1 - conf.level)/2, sncol - 1, lower.tail = FALSE) * sqrt(comlev / (comlev - 1)) * nsd/sqrt(sncol)
}
# 指定のあった指標を表に追加
if(cind){
bstat.info1 <- append(bstat.info1, nd.info)
bstatist <- cbind(bstatist, "CIND-L" = mncol - ci.tval / sqrt(2), "CIND-U" = mncol + ci.tval / sqrt(2))
}
if(cin){
bstat.info1 <- append(bstat.info1, n.info)
bstatist <- cbind(bstatist, "CIN-L" = mncol - ci.tval, "CIN-U" = mncol + ci.tval)
}
}
# マルチレベルモデルに基づく信頼区間の計算
if(ciml){
loadstate <- search()# この時点でのパッケージ呼び出し状態を記録
pkchk <- library(lmerTest, logical.return = TRUE)# パッケージの呼び出し
if(pkchk == FALSE) stop("Please install lmerTest package to use ciml option!!")# パッケージがなければストップ
def.contr <- suppressWarnings(lmerControl()$checkControl)# lmerControlのデフォルト設定を保存
options(lmerControl = list(check.nobs.vs.nlev = "ignore", check.nlev.gtr.1 = "ignore", check.nobs.vs.nRE = "ignore"))# 各種警告メッセージをオフ
# モデル式の作成
modeleq <- paste0("y ~ ", gsub(", ", " * ", toString(LETTERS[maxfact:1])))# 固定効果の指定
if(withlen > 0){# 変量効果の指定;ランダム切片
if(withlen == 1){
randeq <- " + (1 | s)"
}else{
randeq <- paste0(" + (1 | s)", paste0(" + (1 | ", LETTERS[(betlen+1):maxfact], ":s)", collapse = ""))
}
modeleq <- paste0(modeleq, randeq)
}
# 推定と信頼区間の算出
ml.model <- lmer(formula(modeleq), na.action = na.omit, REML = TRUE, data = dat)
lscis <- lsmeans(ml.model, test.effs = paste0(LETTERS[maxfact:1], collapse = ":"), ddf = "Kenward-Roger")
options(lmerControl = def.contr)# lmerControlの設定をもどす
newstate <- search()# 現時点でのパッケージ呼び出し状態
detlist <- newstate[!is.element(newstate, loadstate)]# 新たに呼び出したパッケージだけ選択
detdummy <- sapply(detlist, function(x) detach(x, character.only = TRUE))# パッケージの解除
# 指定のあった指標を表に追加
bstat.info1 <- append(bstat.info1, "== Blouin-Riopelle's Multilevel-Based Confidence Intervals ==")
bstatist <- cbind(bstatist, "CIML-L" = lscis$lsmeans.table$"Lower CI", "CIML-U" = lscis$lsmeans.table$"Upper CI")
}
# 記述統計量の表に信頼区間が追加されている場合は情報を追記
if(ncol(bstatist) != (maxfact + 3)) bstat.info1 <- append(bstat.info1, paste0("== ", 100 * conf.level, "% confidence intervals are calculated. =="))
ciset <- list("bstat.info1" = bstat.info1, "bstatist" = bstatist)
# ペアワイズ信頼区間の計算
if(cipaird || cipair){
if(withlen == 0){# 被験者内要因がないとき(被験者間計画の場合)
bstat.info2 <- "*** CAUTION! Pairwise confidence intervals are not suitable for between-subject designs. ***"
ciset <- append(ciset, list("bstat.info2" = bstat.info2))
}else{# 被験者内要因があるとき
if(betlen == 0){# 被験者間要因がないとき(被験者内計画の場合)
caldat <- list(dat)
betlabels <- data.frame(row.names = 1:(comlev * (comlev - 1) / 2))
}else{# 1つ以上の被験者間要因があるとき(混合要因計画の場合)
caldat <- split(dat, dat[, (maxfact+1-withlen):2])
betlabels <- expand.grid(lapply((maxfact+1-withlen):2, function(x) levels(dat[, x])))
betlabels <- betlabels[, order(betlen:1), drop = FALSE]# アルファベット順に列を並べ替え
betlabels <- betlabels[rep(1:nrow(betlabels), each = comlev * (comlev - 1) / 2), , drop = FALSE]
row.names(betlabels) <- NULL
names(betlabels) <- factnames[1:betlen]
}
wdat <- lapply(caldat, function(x) data.frame(matrix(x$y[order(eval(parse(text = paste0("x[, ", c(1, (maxfact+1):(betlen+2)), "]"))))], length(x$y)/comlev, comlev)))# データをワイド形式のデータフレームに変換;tapplyによるアクセス順に対応させるための並べ替え
ddat <- lapply(wdat, function(x) combn(comlev, 2, function(y) x[, y[1]] - x[, y[2]]))# ペアごとの差分得点
pairbase <- lapply(ddat, function(x) data.frame("Diff" = colMeans(x), "n" = nrow(x), "S.E." = apply(x, 2, function(y) sd(y)/sqrt(length(y)))))# ペアワイズの差の平均,サンプルサイズ,標準誤差
pairtab <- do.call(rbind, pairbase)
rownames(pairtab) <- NULL
crosslevels <- expand.grid(lapply((maxfact+1):(betlen+2), function(x) levels(dat[, x])), stringsAsFactors = FALSE)
comblabels <- apply(crosslevels, 1, function(x) gsub(", ", ".", toString(x[withlen:1])))
withlabels <- rep(combn(comlev, 2, function(x) gsub(", ", "-", toString(comblabels[x]))), times = prod(flev[pmin(1:betlen, betlen)]))
pairtab <- cbind(betlabels, "Pairs" = withlabels, pairtab)
ci.tval <- qt((1 - conf.level)/2, pairtab$"n" - 1, lower.tail = FALSE) * pairtab$"S.E."
# 指定のあった指標を表に追加
if(cipaird){
bstat.info2 <- append(bstat.info2, "== Franz-Loftus's Difference-Adjusted Pairwise Confidence Intervals ==")
pairtab <- cbind(pairtab, "CIPRD-L" = pairtab$Diff - ci.tval / sqrt(2), "CIPRD-U" = pairtab$Diff + ci.tval / sqrt(2))
}
if(cipair){
bstat.info2 <- append(bstat.info2, "== Franz-Loftus's Pairwise Confidence Intervals ==")
pairtab <- cbind(pairtab, "CIPR-L" = pairtab$Diff - ci.tval, "CIPR-U" = pairtab$Diff + ci.tval)
}
bstat.info2 <- append(bstat.info2, paste0("== ", 100 * conf.level, "% confidence intervals are calculated. =="))
ciset <- append(ciset, list("bstat.info2" = bstat.info2, "pairtab" = pairtab))
}
}
return(ciset)
}
# 信頼区間つきの棒グラフを作る関数
# 3要因までの計画に対応
# barmin:棒グラフy軸の下限値
ci.bars <- function(dat, design, factnames = NA, conf.level = 0.95, inn.tier = "cind", out.tier = NA, main = NULL, ylab = NULL, barmin = 0){
maxfact <- nchar(design) - 1
flev <- sapply(2:(maxfact + 1), function(x) nlevels(dat[, x]))# 各要因の水準数
# factnamesが省略されているときは名前をつける
if(anyNA(factnames)) factnames <- LETTERS[1:maxfact]
# 信頼区間の計算
eval(parse(text = paste0(c("cilmd", "cilm", "cind", "cin", "ciml"), "v <- FALSE")))
eval(parse(text = paste0(inn.tier, "v <- TRUE")))
eval(parse(text = paste0(out.tier, "v <- TRUE")))
baseresults <- ci.calc(dat, design, factnames = factnames, conf.level = conf.level, cilmd = cilmdv, cilm = cilmv, cind = cindv, cin = cinv, ciml = cimlv)
bstatist <- baseresults$bstatist
# グラフィックデバイスを開く
if(sum(grep("mac", .Platform)) != 0){# Macの場合
quartz(title = paste0("Bar Graph for ", design, "-Type Design"))
}else{# その他のOSの場合
x11(title = paste0("Bar Graph for ", design, "-Type Design"))
}
# グラフパラメータの設定
preset <- par(no.readonly = TRUE)
par(lwd = 1.5, cex = 1.8, plt = c(0.2, 0.7, 0.2, 0.87), mgp = c(2.5, 0.5, 0), xpd = TRUE)
if(is.na(out.tier) == TRUE){# 外側エラーバーがない場合
mMax <- max(bstatist[, paste0(toupper(inn.tier), "-U")])# 最大値
mMin <- min(bstatist[, paste0(toupper(inn.tier), "-L")])# 最小値
}else{# 外側エラーバーがある場合
mMax <- max(bstatist[, paste0(toupper(out.tier), "-U")])# 最大値
mMin <- min(bstatist[, paste0(toupper(out.tier), "-L")])# 最小値
}
mcoef <- 10^(nchar(abs(round(mMax)))-2)
betbars <- rep(c(1, rep(0, flev[maxfact] - 1)), prod(flev)/flev[maxfact])# バーの間隔
if(length(flev) == 3){
gdivs <- (1:prod(flev))[(1:prod(flev) %% (prod(flev)/flev[1])) == 0]
gdivs <- gdivs[-length(gdivs)]
betbars[gdivs + 1] <- 2
}
barcolor <- gray(0.5/flev[maxfact] * flev[maxfact]:1 + 0.45)
# グラフの描画
barplot(matrix(bstatist$Mean, prod(flev)/flev[maxfact], flev[maxfact]),
las = 2, space = betbars, tcl = 0.25, ylab = ylab,
ylim = c(pmin(barmin, floor(mMin/mcoef) * mcoef), ceiling(mMax/mcoef) * mcoef),
xlim = c(0.5, length(betbars) + (prod(flev)/flev[maxfact])^(1 * (sum(betbars) != 1)) + 1 * sum(betbars == 2)),
col = barcolor,
xpd = FALSE, beside = TRUE)
box(bty = "l")# x軸の線を引く
title(main, line = 0.8)# タイトルを書く
# x軸のラベル
if(maxfact == 1){# 1要因のときは水準名の代わりに要因名を表示する
xlabpos <- flev/2 + 1
text(xlabpos, 0, pos = 1, offset = 1, labels = factnames)
}else{
for(i in pmax(maxfact - 1, 1):1){# x軸に水準名を表示する
xlabframe <- unique(bstatist[, 1:i, drop = FALSE])
xtf <- betbars >= (maxfact - i)
xtf[1] <- TRUE
xcue <- (1:length(betbars))[xtf] + mean(1:(prod(flev)/prod(flev[1:i]))) - 1
xlabpos <- cumsum(c(0.5, betbars[-1]))[xcue] + xcue# x軸ラベルの位置
text(xlabpos, pmin(barmin, floor(mMin/mcoef) * mcoef), pos = 1, offset = 1.2 * pmax(maxfact - i, 1) - 0.4, labels = xlabframe[, i])
}
}
# 凡例
legend(length(betbars) + sum(betbars) + 0.5, ceiling(mMax/mcoef) * mcoef,
legend = levels(bstatist[, maxfact]), cex = 0.9, x.intersp = 0.6, y.intersp = 1,
pch = 22, pt.bg = barcolor, pt.cex = 1.6, # 条件の別を表すマークの設定
bty = "n")
# エラーバー
if(is.na(out.tier) == FALSE){# 外側エラーバーが指定されている場合
arrows(cumsum(betbars) + 1:prod(flev) - 0.5, bstatist[, paste0(toupper(out.tier), "-L")], # 外側エラーバー
cumsum(betbars) + 1:prod(flev) - 0.5, bstatist[, paste0(toupper(out.tier), "-U")],
angle = 90, length = 0.04, code = 3, lwd = 2)
}
arrows(cumsum(betbars) + 1:prod(flev) - 0.5, bstatist[, paste0(toupper(inn.tier), "-L")], # 内側エラーバー
cumsum(betbars) + 1:prod(flev) - 0.5, bstatist[, paste0(toupper(inn.tier), "-U")],
angle = 90, length = 0.06, code = 3, lwd = 2)
par(preset)# グラフパラメータの設定を元に戻す
}
# 文字列中のすべての要素を含む文字列をマッチングする関数
# grepとの違いは“A:C”などの文字列を照合パターンとした場合に“A:B:C”のように間に別の文字を挟んだ文字列もマッチと判定する点
# 照合パターンが1文字の場合はgrepと同じ結果を返す
elematch <- function(Mstrings, stex){
# マッチングする文字列を分解して,それぞれgrep関数を適用
matchlist <- lapply(strsplit(Mstrings, "")[[1]], function(x) grep(x, stex))
# 文字列の各要素とマッチした値の共通部分のみ取り出す
buffer <- matchlist[[1]]
# 文字列が1文字のときはgrepの結果をそのまま返す
if(length(matchlist) != 1){
for(i in 2:length(matchlist)){
buffer <- buffer[is.element(buffer, matchlist[[i]])]
}
}
return(buffer)
}
# ベクトルの組み合わせを作る関数
# expand.gridと同様の組み合わせを作る;結果をmatrix型で返す
expand.gmatrix <- function(...){
elem <- list(...)
if(is.list(elem[[1]])) elem <- unlist(elem, recursive = FALSE)
elemsize <- sapply(elem, length)
rcue <- c(1, cumprod(elemsize)[-length(elemsize)])
rlev <- cumprod(elemsize)
levmax <- prod(elemsize)
expmat <- mapply(function(w, x, y, z) rep.int(rep.int(w, rep.int(y, x)), prod(elemsize)/z), elem, elemsize, rcue, rlev)
return(expmat)
}
# 有意水準に合わせて記号を表示する関数
sig.sign <- function(pvalue){
ifelse(is.na(pvalue), "",
ifelse(pvalue < 0.001, "***",
ifelse(pvalue < 0.01, "**",
ifelse(pvalue < 0.05, "*",
ifelse(pvalue < 0.10, "+", "ns")))))
}
# Greenhouse-GeisserとHuynh-Feldtのイプシロンを計算する関数
# 被験者内要因を含まない計画を投入すると適切に動作しないので注意
epsilon.calc <- function(dat, design, mau = FALSE, har = FALSE, iga = FALSE, ciga = FALSE, lb = FALSE, gg = FALSE, hf = FALSE,
cm = FALSE, autov = NULL, flev = NULL, cellN = NULL, esboot = FALSE){
maxfact <- nchar(design) - 1# 実験計画全体における要因数
spos <- match("s", strsplit(design, "")[[1]], nomatch = 0)
betlen <- spos - 1# 被験者間要因の数
withlen <- maxfact - spos + 1# 被験者内要因の数
if(is.null(flev)) flev <- sapply(2:(maxfact+1), function(x) nlevels(dat[, x]))# 各要因の水準数
if(is.null(cellN)) cellN <- length(unique(dat$s))# 被験者間要因をつぶしての全被験者の数を取得
# 被験者内要因の特定
replabel <- (maxfact - withlen + 2):(maxfact + 1)
repnum <- sort(replabel, decreasing = TRUE)# 後の方の要因のラベルを先に並べる
repmat <- flev[repnum - 1]# 各要因の水準数をベクトル化
rl <- prod(repmat)# 全被験者内要因の組み合わせ水準数を取得
# データフレームを分割し,共分散行列を作る
if(betlen == 0){# 被験者間要因がないときはデータフレームを分割しない
othlabel <- 1
ol <- 1
othN <- cellN# 被験者間要因がないときはcellNと同じ値を代入
covmatrices <- list(cov(do.call(cbind, split(dat$y, dat[, repnum]))))
}else{# 被験者間要因の組み合わせ水準ごとにデータフレームを分割
# 被験者間要因の特定
othlabel <- 1:betlen + 1
othnum <- sort(othlabel, decreasing = TRUE)# 後の方の要因のラベルを先に並べる
othmat <- flev[othnum - 1]# 各要因の水準数をベクトル化
ol <- prod(othmat)# 全被験者間要因の組み合わせ水準数
sdat <- split(dat$y, dat[, othnum])
othN <- sapply(sdat, length) / rl# 被験者間要因の各組み合わせにおける被験者数をベクトル化
covmatrices <- lapply(sdat, function(x) cov(matrix(x, ncol = rl)))
}
# 複数の共分散行列をプール
tm <- Reduce(f = "+", x = lapply(1:ol, function(w) (othN[w] - 1) * covmatrices[[w]])) / (cellN - ol)
# 正規直交対比行列を作る;被験者内要因の数に合わせて異なるパターンを得る
combmat <- expand.gmatrix(lapply(repmat - 1, function(x) 0:x))
cuemat <- cbind(1:nrow(combmat), colSums(t(combmat != 0) * 10^((withlen-1):0)), rowSums(combmat != 0))
ortho.ord <- unlist(lapply(0:withlen, function(x) cuemat[cuemat[, 3] == x, 1][sort.list(cuemat[cuemat[, 3] == x, 2])]))
ortho.helm <- Reduce(f = kronecker, x = lapply(flev[replabel - 1], function(w) cbind(1, contr.helmert(w))))
ortho.helm <- ortho.helm[, ortho.ord]# 解釈しやすい順に並べ替え
ortho.coef <- t(ortho.helm[, 2:rl, drop = FALSE])# 直交対比のパターンのみを取り出す;行が1のときにベクトルに変換されないようにdrop = FALSEを使う
ortho.denomi <- rowSums(ortho.coef^2)^(1/2)
# パターンを直交対比行列にする
orthoM <- ortho.coef / ortho.denomi
# 被験者内要因の数によって処理を変更する
if(withlen == 1){# 被験者内要因が1つのとき
matdivider <- flev[replabel - 1] - 1# 正規直交行列を分割する際の行数
effect.name <- paste0(names(dat)[replabel], collapse = ":")# 効果のラベルを作る
ss.name <- paste0(paste0(c("s", LETTERS[othlabel - 1]), collapse = ":"), ":", effect.name)
# 正規直交対比行列関連の処理
seportM <- list(orthoM)
totalM <- seportM
gmd <- matdivider
otoM <- orthoM %*% tm %*% t(orthoM)# 共分散行列と正規直交対比行列をかける
ss.er <- (cellN - ol) * sum(diag(otoM))
df.er <- (cellN - ol) * nrow(otoM)
otoM <- list(otoM)
}else{# 被験者内要因が複数あるとき
ctlframe <- unlist(lapply(1:withlen, function(y) combn(replabel, y, function(x) x, simplify = FALSE)), recursive = FALSE)
matdivider <- unlist(lapply(ctlframe, function(x) prod(flev[x - 1] - 1)))# 正規直交行列を分割する際の行数
effect.name <- unlist(lapply(ctlframe, function(x) paste0(names(dat)[x], collapse = ":")))# 効果のラベルを作る
ss.name <- paste0(paste0(c("s", LETTERS[othlabel - 1]), collapse = ":"), ":", effect.name)
# ラベルの追加
effect.name <- c("Global", effect.name)
# 正規直交対比行列を被験者内要因の水準数によって分割
divpoint <- mapply(function(x, y) seq(x, y), cumsum(matdivider) - (matdivider - 1), cumsum(matdivider), SIMPLIFY = FALSE)
seportM <- lapply(divpoint, function(x) orthoM[x, , drop = FALSE])
totalM <- c(list(orthoM), seportM)
gmd <- c(rl - 1, matdivider)
otoM <- lapply(totalM, function(x) x %*% tm %*% t(x))# 共分散行列と正規直交対比行列をかける
ss.er <- sapply(otoM, function(x) (cellN - ol) * sum(diag(x)))[-1]
df.er <- sapply(otoM, function(x) (cellN - ol) * nrow(x))[-1]
}
ss.errs <- matrix(c(ss.er, df.er), ncol = 2, dimnames = list(ss.name, c("ss.col", "df.col")))
# イプシロンを計算する
LB.ep <- 1 / gmd
GG.ep <- sapply(otoM, function(x) sum(diag(x))^2 / (nrow(x) * sum(x^2)))
HF.ep <- ((cellN - ol + 1) * gmd * GG.ep - 2) / (gmd * (cellN - ol - gmd * GG.ep))# Lecoutre(1991)の修正
va <- (cellN - ol - 1) + (cellN - ol) * (cellN - ol - 1) / 2
CM.ep <- pmax(LB.ep, HF.ep * (va - 2) * (va - 4) / va^2)
# 球面性検定の実施
if(esboot){# ブートストラップ計算の場合はスキップ
if(lb){
sph.ep <- LB.ep[min(withlen, 2):length(gmd)]
}else if(gg){
sph.ep <- GG.ep[min(withlen, 2):length(gmd)]
}else if(hf){
sph.ep <- HF.ep[min(withlen, 2):length(gmd)]
}else if(cm){
sph.ep <- CM.ep[min(withlen, 2):length(gmd)]
}else if(!is.null(autov)){
sph.ep <- GG.ep[min(withlen, 2):length(gmd)]
sph.ep[autov] <- 1
}else{
sph.ep <- rep(1, nrow(ss.errs))
}
return(list("ss.errs" = ss.errs, "sph.ep" = sph.ep))
}else{
if(mau){# Mauchlyの球面性検定
# プロンプトの準備
epsi.info1 <- paste0("== Mauchly's Sphericity Test and Epsilons ==")
lamlab <- "W"
eps.Lambda <- sapply(otoM, function(x) det(x) / (sum(diag(x)) / nrow(x))^nrow(x))
eps.m <- 1 - (2 * gmd^2 + gmd + 2) / (6 * (cellN - ol) * gmd)
epsChi <- -(cellN - ol) * eps.m * log(eps.Lambda)
if(any(min(othN) < gmd)){# 各群の被験者数が正規直交対比行列の行数を下回るときは妥当なカイ二乗値を計算できない
epsChi[min(othN) < gmd] <- NA
eps.Lambda[min(othN) < gmd] <- NA
epsi.info1 <- paste0(epsi.info1, "\n",
"*** CAUTION! The test of SPHERICITY is INVALID because of small sample size. ***", "\n",
"*** The minimum sample size for valid computation is N = ", max(gmd) + 1, " at each group. ***")
}
eps.df <- gmd * (gmd + 1)/ 2 - 1
eps.p1 <- pchisq(epsChi, ifelse(eps.df == 0, NA, eps.df), lower.tail = FALSE)
eps.p2 <- pchisq(epsChi, eps.df + 4, lower.tail = FALSE)
eps.w2 <- (gmd + 2) * (gmd - 1) * (gmd - 2) * (2 * gmd^3 + 6 * gmd^2 + 3 * gmd + 2) / (288 * gmd^2 * (cellN - ol)^2 * eps.m^2)
eps.p <- eps.p1 + eps.w2 * (eps.p2 - eps.p1)
}else if(har){# Harrisの多標本球面性検定
# プロンプトの準備
epsi.info1 <- paste0("== Harris's Multisample Sphericity Test and Epsilons ==")
lamlab <- "h_hat"
proA <- lapply(totalM, function(x) lapply(covmatrices, function(y) x %*% y %*% t(x)))
epsTr <- lapply(proA, function(y) sapply(y, function(x) sum(diag(x))))
epsSq <- lapply(proA, function(y) sapply(y, function(x) sum(diag(x %*% x))))
eps.Lambda <- sapply(1:length(proA), function(x) sum((othN - 1) * epsTr[[x]])^2 / sum((othN - 1) * epsSq[[x]]))
epsChi <- pmax(0, ((cellN - ol) * gmd / 2) * ((cellN - ol) * gmd / eps.Lambda - 1))# 負の値は0にそろえる
if(any(min(othN) < gmd)){# 各群の被験者数が正規直交対比行列の行数を下回るときは妥当なカイ二乗値を計算できない
epsChi[min(othN) < gmd] <- NA
eps.Lambda[min(othN) < gmd] <- NA
epsi.info1 <- paste0(epsi.info1, "\n",
"*** CAUTION! The test of SPHERICITY is INVALID because of small sample size. ***", "\n",
"*** The minimum sample size for valid computation is N = ", max(gmd) + 1, " at each group. ***")
}
eps.df <- ((ol * gmd * (gmd + 1)) / 2) - 1
eps.p0 <- pchisq(epsChi, ifelse(eps.df == 0, NA, eps.df), lower.tail = FALSE)
eps.p6 <- pchisq(epsChi, eps.df + 6, lower.tail = FALSE)
eps.p4 <- pchisq(epsChi, eps.df + 4, lower.tail = FALSE)
eps.p2 <- pchisq(epsChi, eps.df + 2, lower.tail = FALSE)
eps.p <- pmax(0, eps.p0 + ((gmd^3 + 3 * gmd^2 - 8 * gmd - 12 - 200/gmd) * eps.p6 / 12
+ (-2 * gmd^3 - 5 * gmd^2 + 7 * gmd + 12 + 420/gmd) * eps.p4 / 8
+ (gmd^3 + 2 * gmd^2 - gmd - 2 - 216/gmd) * eps.p2 / 4
+ (-2 * gmd^3 - 3 * gmd^2 + gmd + 436/gmd) * eps.p0 / 24
) / (cellN - ol))
}else{# Mendozaの多標本球面性検定
# プロンプトの準備
epsi.info1 <- paste0("== Mendoza's Multisample Sphericity Test and Epsilons ==")
lamlab <- "Lambda"
proA <- lapply(totalM, function(x) lapply(1:ol, function(y) x %*% (othN[y] * covmatrices[[y]]) %*% t(x)))
eps.m <- 1 - ((((cellN-ol) * gmd^2 * (gmd + 1) * (2 * gmd + 1) - (2*(cellN-ol) * gmd^2)) *
sum(1/(othN-1)) - 4) / (6 * (cellN-ol) * gmd * (ol * gmd * (gmd + 1) - 2)))
eps.m[is.nan(eps.m)] <- 0# NaNが出たところには0を代入
menL1 <- log(cellN-ol) * ((cellN-ol)*(gmd)/2) - sapply(gmd, function(x) sum(log(othN-1) * ((othN-1) * x / 2)))
menL2 <- sapply(proA, function(y) sum(sapply(y, function(x) determinant(x, logarithm = TRUE)$modulus[1]) * (othN-1)/2))
menL3 <- sapply(proA, function(y) sum(diag(Reduce(f = "+", x = y)/nrow(y[[1]]))))
menL3 <- log(ifelse(menL3 < 0, NA, menL3)) * ((cellN-ol) * gmd / 2)# トレースが負のときは対数に変換できないのでNAに置き換え
menL <- menL1 + menL2 - menL3
epsChi <- - 2 * eps.m * menL
eps.Lambda <- exp(menL)
if(any(min(othN) < gmd)){# 各群の被験者数が正規直交対比行列の行数を下回るときは妥当なカイ二乗値を計算できない
epsChi[min(othN) < gmd] <- NA
eps.Lambda[min(othN) < gmd] <- NA
epsi.info1 <- paste0(epsi.info1, "\n",
"*** CAUTION! The test of SPHERICITY is INVALID because of small sample size. ***", "\n",
"*** The minimum sample size for valid computation is N = ", max(gmd) + 1, " at each group. ***")
}