forked from wangzhen89/GLMM
-
Notifications
You must be signed in to change notification settings - Fork 0
/
chap7.html
623 lines (600 loc) · 82 KB
/
chap7.html
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
<!DOCTYPE html>
<html lang="en">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<meta charset="utf-8">
<meta name="viewport" content="width=device-width, initial-scale=1, shrink-to-fit=no">
<title>第 7 章 推断（二） | 广义线性混合模型</title>
<meta name="author" content="Wang Zhen">
<meta name="description" content="7.1 介绍 在第 6 章中，我们重点讨论了可估和可预测函数的推断。在本章中，我们将注意力转向 \(\symbf\sigma\) 的元素，即方差和协方差分量。使用第 4 章开头的约定，当我们提到“协方差分量”时，我们指的是 \(\symbf\sigma\) 的所有部分：方差、协方差，在某些情况下还包括相关参数。 我们可能对协方差分量的推断感兴趣，因为我们的目标针对的是...">
<meta name="generator" content="bookdown 0.38 with bs4_book()">
<meta property="og:title" content="第 7 章 推断（二） | 广义线性混合模型">
<meta property="og:type" content="book">
<meta property="og:description" content="7.1 介绍 在第 6 章中，我们重点讨论了可估和可预测函数的推断。在本章中，我们将注意力转向 \(\symbf\sigma\) 的元素，即方差和协方差分量。使用第 4 章开头的约定，当我们提到“协方差分量”时，我们指的是 \(\symbf\sigma\) 的所有部分：方差、协方差，在某些情况下还包括相关参数。 我们可能对协方差分量的推断感兴趣，因为我们的目标针对的是...">
<meta name="twitter:card" content="summary">
<meta name="twitter:title" content="第 7 章 推断（二） | 广义线性混合模型">
<meta name="twitter:description" content="7.1 介绍 在第 6 章中，我们重点讨论了可估和可预测函数的推断。在本章中，我们将注意力转向 \(\symbf\sigma\) 的元素，即方差和协方差分量。使用第 4 章开头的约定，当我们提到“协方差分量”时，我们指的是 \(\symbf\sigma\) 的所有部分：方差、协方差，在某些情况下还包括相关参数。 我们可能对协方差分量的推断感兴趣，因为我们的目标针对的是...">
<!-- JS --><script src="https://cdnjs.cloudflare.com/ajax/libs/clipboard.js/2.0.6/clipboard.min.js" integrity="sha256-inc5kl9MA1hkeYUt+EC3BhlIgyp/2jDIyBLS6k3UxPI=" crossorigin="anonymous"></script><script src="https://cdnjs.cloudflare.com/ajax/libs/fuse.js/6.4.6/fuse.js" integrity="sha512-zv6Ywkjyktsohkbp9bb45V6tEMoWhzFzXis+LrMehmJZZSys19Yxf1dopHx7WzIKxr5tK2dVcYmaCk2uqdjF4A==" crossorigin="anonymous"></script><script src="https://kit.fontawesome.com/6ecbd6c532.js" crossorigin="anonymous"></script><script src="libs/jquery-3.6.0/jquery-3.6.0.min.js"></script><meta name="viewport" content="width=device-width, initial-scale=1, shrink-to-fit=no">
<link href="libs/bootstrap-4.6.0/bootstrap.min.css" rel="stylesheet">
<script src="libs/bootstrap-4.6.0/bootstrap.bundle.min.js"></script><script src="libs/bs3compat-0.7.0/transition.js"></script><script src="libs/bs3compat-0.7.0/tabs.js"></script><script src="libs/bs3compat-0.7.0/bs3compat.js"></script><link href="libs/bs4_book-1.0.0/bs4_book.css" rel="stylesheet">
<script src="libs/bs4_book-1.0.0/bs4_book.js"></script><script type="text/x-mathjax-config">
MathJax.Hub.Config({
"HTML-CSS": {
fonts: ["STIX-Web"]
},
SVG: {
font: "STIX-Web"
},
TeX: {Augment: {
Definitions: {macros: {symbf: 'Symbf'}},
Parse: {prototype: {
csMathchar0mi: function (name, mchar) {
var MML = MathJax.ElementJax.mml;
var def = {};
if (Array.isArray(mchar)) {def = mchar[1]; mchar = mchar[0]}
this.Push(this.mmlToken(MML.mi(MML.entity("#x"+mchar)).With(def)));
},
Symbf: function (name) {
var MML = MathJax.ElementJax.mml;
var math = this.ParseArg(name);
this.Push(MML.mstyle(math).With({mathvariant: "bold"}));
}
}}
}}
});
</script><script src="https://cdnjs.cloudflare.com/ajax/libs/autocomplete.js/0.38.0/autocomplete.jquery.min.js" integrity="sha512-GU9ayf+66Xx2TmpxqJpliWbT5PiGYxpaG8rfnBEk1LL8l1KGkRShhngwdXK1UgqhAzWpZHSiYPc09/NwDQIGyg==" crossorigin="anonymous"></script><script src="https://cdnjs.cloudflare.com/ajax/libs/mark.js/8.11.1/mark.min.js" integrity="sha512-5CYOlHXGh6QpOFA/TeTylKLWfB3ftPsde7AnmhuitiTX4K5SqCLBeKro6sPS8ilsz1Q4NRx3v8Ko2IBiszzdww==" crossorigin="anonymous"></script><!-- CSS --><style type="text/css">
div.hanging-indent{margin-left: 1.5em; text-indent: -1.5em;}
</style>
<link rel="stylesheet" href="style.css">
</head>
<body data-spy="scroll" data-target="#toc">
<div class="container-fluid">
<div class="row">
<header class="col-sm-12 col-lg-3 sidebar sidebar-book"><a class="sr-only sr-only-focusable" href="#content">Skip to main content</a>
<div class="d-flex align-items-start justify-content-between">
<h1>
<a href="index.html" title="现代概念、方法和应用">广义线性混合模型</a>:
<small class="text-muted">现代概念、方法和应用</small>
</h1>
<button class="btn btn-outline-primary d-lg-none ml-2 mt-1" type="button" data-toggle="collapse" data-target="#main-nav" aria-expanded="true" aria-controls="main-nav"><i class="fas fa-bars"></i><span class="sr-only">Show table of contents</span></button>
</div>
<div id="main-nav" class="collapse-lg">
<form role="search">
<input id="search" class="form-control" type="search" placeholder="Search" aria-label="Search">
</form>
<nav aria-label="Table of contents"><h2>Table of contents</h2>
<ul class="book-toc list-unstyled">
<li><a class="" href="index.html">译者序</a></li>
<li><a class="" href="%E6%89%89%E9%A1%B5.html">扉页</a></li>
<li><a class="" href="%E7%9B%AE%E5%BD%95.html">目录</a></li>
<li><a class="" href="secpre.html">前言</a></li>
<li class="book-part">第一篇：基本背景</li>
<li><a class="" href="chap1.html"><span class="header-section-number">1</span> 建模基础</a></li>
<li><a class="" href="chap2.html"><span class="header-section-number">2</span> 设计要务</a></li>
<li><a class="" href="chap3.html"><span class="header-section-number">3</span> 搭建舞台</a></li>
<li><a class="" href="%E6%90%AD%E5%BB%BA%E8%88%9E%E5%8F%B0.html">►搭建舞台</a></li>
<li class="book-part">第二篇：估计和推断理论</li>
<li><a class="" href="chap4.html"><span class="header-section-number">4</span> GLMM 之前的估计和推断基础知识</a></li>
<li><a class="" href="chap5.html"><span class="header-section-number">5</span> GLMM 估计</a></li>
<li><a class="" href="chap6.html"><span class="header-section-number">6</span> 推断（一）</a></li>
<li><a class="active" href="chap7.html"><span class="header-section-number">7</span> 推断（二）</a></li>
<li class="book-part">第三篇：应用</li>
<li><a class="" href="chap8.html"><span class="header-section-number">8</span> 处理和解释变量结构</a></li>
<li><a class="" href="chap9.html"><span class="header-section-number">9</span> 多水平模型</a></li>
<li class="book-part">—</li>
<li><a class="" href="bib.html">参考文献</a></li>
</ul>
<div class="book-extra">
</div>
</nav>
</div>
</header><main class="col-sm-12 col-md-9 col-lg-7" id="content"><div id="chap7" class="section level1" number="7">
<h1>
<span class="header-section-number">第 7 章</span> 推断（二）<a class="anchor" aria-label="anchor" href="#chap7"><i class="fas fa-link"></i></a>
</h1>
<div id="sec7-1" class="section level2" number="7.1">
<h2>
<span class="header-section-number">7.1</span> 介绍<a class="anchor" aria-label="anchor" href="#sec7-1"><i class="fas fa-link"></i></a>
</h2>
<p>在第 <a href="chap6.html#chap6">6</a> 章中，我们重点讨论了可估和可预测函数的推断。在本章中，我们将注意力转向 <span class="math inline">\(\symbf\sigma\)</span> 的元素，即方差和协方差分量。使用第 <a href="chap4.html#chap4">4</a> 章开头的约定，当我们提到“协方差分量”时，我们指的是 <span class="math inline">\(\symbf\sigma\)</span> 的所有部分：方差、协方差，在某些情况下还包括相关参数。</p>
<p>我们可能对协方差分量的推断感兴趣，因为我们的目标针对的是 <span class="math inline">\(\symbf\sigma\)</span> 本身的元素，或者我们可能只需要评估相互竞争的协方差模型作为一种手段来达到目的——将“正确的”协方差模型作为对可估或可预测函数进行推断的前奏。直接关注协方差分量的常见场景包括遗传学，其中协方差分量度量关注性状的遗传力，以及质量改进，其中协方差分量确定问题区域。“达到目的的手段”的应用包括检验关于协方差的假定（例如，我们能否证明在单向方差分析中可以假定方差相等？），或者在重复测量或空间数据的相互竞争的相关误差模型中进行选择。</p>
<p>我们可将协方差分量的推断分为三大类：1) 正式检验；2) 拟合统计量和 3) 区间估计。当我们有明显的嵌套层次结构时，我们只能对协方差分量进行正式检验。我们在 <a href="chap7.html#sec7-2">7.2</a> 节中介绍了正式的协方差分量检验。正式检验并不总是可能的，甚至是不可取的。我们经常想要比较非嵌套的协方差模型。在这种情况下，我们必须使用拟合统计量。我们考虑非嵌套比较的例子，并在 <a href="chap7.html#sec7-3">7.3</a> 节中介绍拟合统计量。正式检验和拟合统计量要求我们有良定的似然。当我们对非高斯 GLMM 使用伪似然估计时，情况并非如此。因此，非高斯 GLMM 为协方差分量的推断带来了其他线性模型前所未有的挑战。我们在 <a href="chap7.html#sec7-2-3">7.2.3</a> 节中考虑了用于检验的这些问题，并在 <a href="chap7.html#sec7-3">7.3</a> 节中回顾了用于拟合统计量的这些问题。最后，<a href="chap7.html#sec7-4">7.4</a> 节重点讨论协方差分量的置信区间构建。</p>
</div>
<div id="sec7-2" class="section level2" number="7.2">
<h2>
<span class="header-section-number">7.2</span> 协方差分量的正式检验<a class="anchor" aria-label="anchor" href="#sec7-2"><i class="fas fa-link"></i></a>
</h2>
<p>虽然检验协方差分量的方法有很多，但这里我们将考虑三种方法，并主要关注一种方法——似然比检验。似然比检验是协方差分量的一种正式检验，足够灵活，可以与 GLMM 及其所有特殊情况（LMM, GLM 和 LM）一起使用。我们还简单地考虑基于 ANOVA 的检验，因为它们在线性模型历史中很重要。然而，与基于 ANOVA 的方差分量估计一样，基于 ANOVA 的检验只能用于仅方差分量的 LMMs，因此它们在当代线性建模中的价值有限。最后，我们考虑基于 Wald 的检验。我们这样做主要是因为他们的小样本表现非常糟糕，从而向学习线性模型的学生发出警告。</p>
<div id="sec7-2-1" class="section level3" number="7.2.1">
<h3>
<span class="header-section-number">7.2.1</span> 仅方差分量 LMMs 的基于 ANOVA 的检验<a class="anchor" aria-label="anchor" href="#sec7-2-1"><i class="fas fa-link"></i></a>
</h3>
<p>顾名思义，基于 ANOVA 的方差分量检验涉及将 ANOVA 表中的 F 检验应用于代表随机模型效应的变异源。我们可以通过例子来说明基于 ANOVA 的检验。考虑第 <a href="chap3.html#chap3">3</a> 章 Data Set 3.2 使用的四处理、八地点数据。这些数据的 ANOVA 为</p>
<div class="inline-figure"><img src="figure/figure%20210.png#center" style="width:80.0%"></div>
<p>通过检查期望均方，我们发现地点的 <span class="math inline">\(F\)</span> 值可用作 <span class="math inline">\(H_0\colon{\sigma}_L^2=0\)</span> 的检验统计量。给定 <span class="math inline">\(F = 3.41\)</span> 且 <span class="math inline">\(p = 0.0135\)</span>，对于任何大于 <span class="math inline">\(p\)</span> 值的 <span class="math inline">\(\alpha\)</span>，有足够的证据得出 <span class="math inline">\(H_0\colon{\sigma}_L^2>0\)</span>。正如 Self and Liang (1987) 指出的，我们应该将 ANOVA 的 <span class="math inline">\(p\)</span> 值除以 2，因为根据定义，ANOVA <span class="math inline">\(F\)</span> 检验是双侧的，而方差分量的检验必须是单侧检验。即，<span class="math inline">\(H_0\colon{\sigma}_L^2\ne 0\)</span> 不能作为备择假设，因为 <span class="math inline">\(\sigma_L^2\)</span> 不能为负。此检验的正确 <span class="math inline">\(p\)</span> 值为 <span class="math inline">\(0.0135/2\cong0.0068\)</span>。</p>
<p>所有 ANOVA 检验都以这种方式进行。我们通过检验期望均方来构造 <span class="math inline">\(F\)</span> 检验。对于更复杂的仅方差分量模型，分母可以是均方的线性组合。例如，假设我们有一个裂区实验，其中有两个固定效应因子 A 和 C 以及一个随机效应因子 B. A 和 B 的因子组合随机分配给整区单元，C 的水平随机分配给裂区单元。我们可以将所得 LMM 写为</p>
<ul>
<li>线性预测器：<span class="math inline">\(\eta_{ijkl}=\eta+\alpha_i+b_j+(ab)_{ij}+\gamma_k+(\alpha\gamma)_{ik}+(bc)_{ij}+(abc)_{ijk}+w_{ijk}\)</span>，其中 <span class="math inline">\(\alpha\)</span> 表示 A 效应，<span class="math inline">\(b\)</span> 表示 B 效应，<span class="math inline">\(\gamma\)</span> 表示 C 效应，<span class="math inline">\(w\)</span> 表示整区效应</li>
<li>分布：<span class="math inline">\(y_{ijkl}\mid b_j,(ab)_{ij},(ac)_{jk},(abc)_{ijk},w_{ijl}\mathrm{~iid~}N\left(\mu_{ijkl},\sigma_S^2\right)\)</span>，其中 <span class="math inline">\(\sigma^2_S\)</span> 表示裂区单元之间的方差。<span class="math inline">\(b_j\text{ iid }N\left(0,\sigma_B^2\right)\text{,}\left(ab\right)_{ij}\text{ iid }N\left(0,\sigma_{AB}^2\right)\text{,}\left(bc\right)_{jk}\text{ iid }N\left(0,\sigma_{BC}^2\right)\text{,}\left(abc\right)_{ijk}\text{ iid }N\left(0,\sigma_{ABC}^2\right)\)</span> 以及 <span class="math inline">\(w_{ijk}\mathrm{~iid~}N\big(0,\sigma_W^2\big)\)</span>。</li>
<li>连接：<span class="math inline">\(\eta_{ijkl}=\mu_{ijkl}\)</span>
</li>
</ul>
<p>ANOVA 表为</p>
<div class="inline-figure"><img src="figure/figure%20211.png#center" style="width:70.0%"></div>
<p>代表随机处理效应的变异源以粗体显示。用于检验 <span class="math inline">\(\sigma_{BC}^2\)</span> 和 <span class="math inline">\(\sigma_{ABC}^2\)</span> 的 <span class="math inline">\(F\)</span> 值是简单的 MS 之比。另一方面，分离 <span class="math inline">\(\sigma_{AB}^2\)</span> 和 <span class="math inline">\(\sigma_{B}^2\)</span> 所需的比值需要更复杂的项。具体而言，<span class="math inline">\(\sigma_{AB}^2\)</span> 的分母项必须估计 <span class="math inline">\(\sigma_{S}^{2}+4\sigma_{ABC}^{2}+2\sigma_{W}^{2}\)</span>，即 EMS(AB) 中除 <span class="math inline">\(\sigma_{AB}^2\)</span> 之外的元素。所得检验统计量为</p>
<p><span class="math display">\[F_{AB}=\frac{MS(AB)}{MS(ABC)+MSW-MSR}\]</span></p>
<p>其中 MSW 和 MSR 分别表示全区和残差（裂区）MS。对于 <span class="math inline">\(\sigma^2_B\)</span>，</p>
<p><span class="math display">\[F_B=\frac{MS(B)}{MS(AB)+MS(BC)-MS(ABC)}\]</span></p>
<p>给出所需的检验统计量。这两个检验需要对分母自由度进行 Satterthwaite 近似才能确定 <span class="math inline">\(p\)</span> 值。</p>
<p>这是我们使用 ANOVA EMS 法所能达到的极限。对于仅方差分量 LMM 之外的模型，我们需要寻找其他方法来检验协方差分量。</p>
</div>
<div id="sec7-2-2" class="section level3" number="7.2.2">
<h3>
<span class="header-section-number">7.2.2</span> 用于协方差分量检验的 Wald 统计量以及为什么不应使用<a class="anchor" aria-label="anchor" href="#sec7-2-2"><i class="fas fa-link"></i></a>
</h3>
<p>协方差分量检验的 Wald 统计量使用 <span class="math inline">\(\symbf{V}_A\left(\hat{\symbf{\sigma}}\right)\)</span> 的对角元、<span class="math inline">\(\symbf{\sigma}\)</span> 的渐近协方差阵、协方差分量估计，在第 5 章 <a href="chap5.html#sec5-4-2">5.4.2</a> 节中给出。对于第 <span class="math inline">\(i\)</span> 个协方差分量，渐近方差为</p>
<p><span class="math display" id="eq:7-1">\[\begin{align}
V_{A,ii}=2\times \operatorname{trace}\Bigg[\symbf{P}\Bigg(\frac{\partial \symbf V(\sigma)}{\partial\sigma_i}\Bigg)\symbf{P}\Bigg(\frac{\partial \symbf V(\sigma)}{\partial\sigma_i}\Bigg)\Bigg]^{-1}
\tag{7.1}
\end{align}\]</span></p>
<p>其中 <span class="math inline">\(\symbf{P}=\left[\symbf{V}\left(\hat{\symbf{\sigma}}\right)\right]^{-1}-\left[\symbf{V}\left(\hat{\symbf{\sigma}}\right)\right]^{-1}\symbf{X}\left(\symbf{X}^{\prime}\left[\symbf{V}\left(\hat{\symbf{\sigma}}\right)\right]^{-1}\symbf{X}\right)^{-}\symbf{X}^{\prime}\left[\symbf{V}\left(\hat{\symbf{\sigma}}\right)\right]^{-1}\)</span>。</p>
<p>用于检验 <span class="math inline">\(H_0\colon\sigma_i=0\)</span> 的 Wald 统计量为</p>
<p><span class="math display">\[Wald_{\sigma_i}=\frac{\hat{\sigma}_i}{\sqrt{V_{A,ii}}}\]</span></p>
<p><span class="math inline">\(Wald_{\sigma_i}\)</span> 使用标准高斯 <span class="math inline">\(N (0, 1)\)</span> 分布进行计算。例如，使用 4 处理、8 地点的 Data Set 3.2，我们有以下结果</p>
<ul>
<li>方差估计：<span class="math inline">\(\hat{\symbf{\sigma}}=\begin{bmatrix}\hat{{\sigma}}_L^2\\\hat{{\sigma}}^2\end{bmatrix}=\begin{bmatrix}1.735\\2.877\end{bmatrix}\)</span>
</li>
<li>渐近方差：<span class="math inline">\(\symbf{V}_A\left(\hat{\symbf{\sigma}}\right)=\begin{bmatrix}1.771&-0.197\\-0.197&0.788\end{bmatrix}\)</span>
</li>
</ul>
<p>用于检验 <span class="math inline">\(H_0\colon\sigma_L=0\)</span> 的 Wald 统计量为</p>
<p><span class="math display">\[Wald\left(\sigma_L^2\right)=\frac{1.735}{\sqrt{1.771}}=1.30\]</span></p>
<p><span class="math inline">\(p\)</span> 值为 0.1922。回想，对于基于 ANOVA 的检验，<span class="math inline">\(F = 3.41\)</span> 且 <span class="math inline">\(p = 0.0135\)</span>。</p>
<p>这是典型的协方差分量 Wald 检验：无意义的结果——通常是非常保守的结果。发生这种情况的原因源于方差分量估计的分布。对于仅方差分量的 LMM，<span class="math inline">\(\hat{\symbf\sigma}\)</span> 的元素取决于 <span class="math inline">\(\chi^2\)</span> 随机变量的线性组合，即强右偏的分布。中心极限定理对于 Wald 统计量确实成立，但对于方差估计，收敛到正态的速度非常慢。因此，除非样本量达到数千（最好是数万或数十万），否则 Wald 统计量一无是处。对于较小的样本量，协方差分量估计的抽样分布明显是非高斯的。</p>
<div class="rmdcaution">
<p><strong>底线</strong>：任何受过教育的数据分析师不应使用 Wald 协方差分量检验，但具有数万或更多观测的超大型研究可能除外。</p>
</div>
</div>
<div id="sec7-2-3" class="section level3" number="7.2.3">
<h3>
<span class="header-section-number">7.2.3</span> 协方差分量的似然比检验<a class="anchor" aria-label="anchor" href="#sec7-2-3"><i class="fas fa-link"></i></a>
</h3>
<p>协方差分量的似然比检验类似于第 <a href="chap6.html#chap6">6</a> 章中 <a href="chap6.html#eq:6-20">(6.20)</a> 定义的可估函数的似然比检验。保持线性预测器的固定效应部分为常数，我们在 <span class="math inline">\(H_0\cdot \symbf \sigma=\symbf\sigma_0\)</span> 下拟合模型，确定似然，然后拟合全模型。所得似然比统计量为</p>
<p><span class="math display" id="eq:7-2">\[\begin{align}
LR\left(\symbf\sigma_0\right)=2\ell\left(\hat{\symbf\sigma}\right)-2\ell\left(\hat{\symbf\sigma}_0\right)
\tag{7.2}
\end{align}\]</span></p>
<p>注意到</p>
<p><span class="math display">\[LR\left(\symbf\sigma_0\right)=-2\log\left[\frac{f\left(\hat{\symbf\sigma}_0\right)}{f\left(\hat{\symbf\sigma}\right)}\right]\]</span></p>
<p>其中 <span class="math inline">\(f(\hat{\symbf \sigma})\)</span> 表示在协方差向量估计 <span class="math inline">\(\hat{\symbf \sigma}\)</span> 处计算的似然。如果我们用高斯数据拟合 LMM 模型，我们应使用 REML 估计 <span class="math inline">\(\symbf\sigma\)</span>，从而使用 REML 似然产生似然比统计量。</p>
<p><span class="math display">\[-2\log\left[\frac{f_R\left(\hat{\symbf\sigma}_0\right)}{f_R\left(\hat{\symbf\sigma}\right)}\right]\]</span></p>
<p>其中 <span class="math inline">\(f_R(\hat{\symbf \sigma})\)</span> 表示 <a href="chap5.html#eq:5-27">(5.27)</a> 所示的 REML 似然。无论我们根据全似然还是 REML 似然计算似然比，<span class="math inline">\(LR(\sigma_0)\)</span> 具有近似 <span class="math inline">\(\chi^2_\nu\)</span> 分布的这一结果都成立。自由度 <span class="math inline">\(\nu\)</span> 等于 <span class="math inline">\(\symbf\sigma\)</span> 中协方差参数的数量与 <span class="math inline">\(\symbf\sigma_0\)</span> 中协方差参数的数量之差。由于是卡方近似，似然比统计量的另一种符号是 <span class="math inline">\(\chi^2_{LR}\)</span>。</p>
<p>举例说明，让我们回到 4 处理、8 地点的 Data Set 3.2 并检验 <span class="math inline">\(H_0\colon\sigma^2_L=0\)</span>。在全模型下</p>
<p><span class="math display">\[\symbf\sigma=\begin{bmatrix}\sigma_L^2\\\sigma^2\end{bmatrix}\]</span></p>
<p>在 <span class="math inline">\(H_0\)</span> 下，<span class="math inline">\(\symbf\sigma_0 = \begin{bmatrix}\sigma^2 \end{bmatrix}\)</span> ，即在原假设下，通过移除随机地点效应来简化模型。在全模型下，REML 似然为 <span class="math inline">\(4.4477× 10^{-28}\)</span>，因此 −2 倍 REML 对数似然为 125.95。在简化模型下，REML 似然为 <span class="math inline">\(4.4148×10^{-29}\)</span>，−2 倍 REML 对数似然为 130.58. 所得似然比统计量为 <span class="math inline">\(LR(\sigma_L^2=0)= 4.62\)</span>。使用 <span class="math inline">\(\chi^2_1\)</span> 评估 <span class="math inline">\(LR\)</span> ——因为全模型有两个方差分量，而简化模型有一个，因此 <span class="math inline">\(\nu = 1\)</span> ——我们得到 <span class="math inline">\(p = 0.0158\)</span>。这与基于 ANOVA 的检验一致。</p>
<div id="sec7-2-3-1" class="section level4" number="7.2.3.1">
<h4>
<span class="header-section-number">7.2.3.1</span> 单向方差分析：方差齐性检验<a class="anchor" aria-label="anchor" href="#sec7-2-3-1"><i class="fas fa-link"></i></a>
</h4>
<p>让我们考虑另一个例子。Data Set 7.1 包含来自四种处理的完全随机设计数据。假定数据是高斯分布的，我们可以考虑以下模型</p>
<ul>
<li>线性预测器：<span class="math inline">\(\eta_i=\eta+\tau_i\)</span>
</li>
<li>分布：通常的等方差假定：<span class="math inline">\(y_{ij}\mathrm{~iid~}N\left(\mu_i,\sigma^2\right)\)</span>；不等方差假定：<span class="math inline">\(y_{ij}\sim NI\left(\mu_i,\sigma^2_i\right)\)</span>
</li>
<li>连接：<span class="math inline">\(\eta_i=\mu_i\)</span>
</li>
</ul>
<p>在继续推断处理效应 <span class="math inline">\(\tau_i\)</span> 之前，我们需要解决方差问题：所有处理的方差相等还是不相等？</p>
<p>我们不能通过基于 ANOVA 的检验来解决这个问题，但我们可以使用似然比检验。这里，检验具有三个自由度，因为不等方差模型具有四个方差参数，<span class="math inline">\(\symbf\sigma'=\begin{bmatrix}\sigma_1^2&\sigma_2^2&\sigma_3^2&\sigma_4^2\end{bmatrix}\)</span>，而等方差模型具有一个方差参数，<span class="math inline">\(\symbf\sigma=\sigma_0^2\)</span>。由此，<span class="math inline">\(\nu=4-1=3\)</span>。在原假设 <span class="math inline">\(H_0\colon\)</span> 所有 <span class="math inline">\(\sigma_i^2=\sigma^2,i=1,2,3,4\)</span> 下，共同方差估计为 <span class="math inline">\(\hat\sigma^2=0.90\)</span> 以及 <span class="math inline">\(-2\ell\left(\hat{\symbf\sigma}_0\right)=38.25\)</span>。在全模型下，<span class="math inline">\(-2\ell\left(\hat{\symbf\sigma}\right)=35.49\)</span>，这给出了似然比统计量 <span class="math inline">\(\chi_{LR}^2=2.86,p=0.4133\)</span>。似然比检验没有提供拒绝等方差模型的证据。</p>
<div class="rmdcaution">
<p><strong>旁注</strong>：在“一般”线性模型时代，当强调 i.i.d. 误差范式的优越性时，如果我们确实拒绝了 <span class="math inline">\(H_0\)</span>，它将引发一场小型危机，以寻找方差稳定变换的形式。而在当代的建模中，我们则直接利用不等方差模型继续进行可估函数的推断工作。</p>
</div>
</div>
<div id="sec7-2-3-2" class="section level4" number="7.2.3.2">
<h4>
<span class="header-section-number">7.2.3.2</span> 重复测量示例：选择简约协方差模型<a class="anchor" aria-label="anchor" href="#sec7-2-3-2"><i class="fas fa-link"></i></a>
</h4>
<p>我们对重复测量数据（更一般地，相关误差模型）的主要讨论在第 <a href="#chap17"><strong>??</strong></a> 章和第 <a href="#chap18"><strong>??</strong></a> 章中。然而，协方差分量检验最重要的应用之一，就是在完整的协方差模型中进行选择。考虑到这一点，让我们看一个例子，这个例子是第 <a href="#chap17"><strong>??</strong></a> 章和第 <a href="#chap18"><strong>??</strong></a> 章示例的预告，并说明如何使用似然比检验来选择协方差模型。</p>
<p>Data Set 7.2 包含对总共 12 名个体进行六次测量的重复测量数据，其中 6 名个体随机分配到两种处理之一，另 6 名个体随机分配到另一种处理。“重复测量数据”是指来自任何研究的数据，其中在研究过程中按计划时间观察到与处理相关的重复单元。这些观察可能是每小时、每天、每月进行一次，也可能是不等间隔的，有策略地计划在研究背景下被认为重要的时间进行。因此，描述 Data Set 7.2 的另一种方式是完全随机设计，其中有两个处理和随机分配给每个处理的六个重复。对每个个体进行六次预先计划的观察；例如在处理开始后的 1, 2, 4, 8, 16 和 24 周。</p>
<p>假定高斯数据，我们可以将模型的一般形式写为：</p>
<ul>
<li>线性预测器：<span class="math inline">\(\eta_{ijk}=\alpha_{ij}+s_{ik}\)</span>，其中 <span class="math inline">\(\alpha_{ij}\)</span> 是处理 <span class="math inline">\(i\)</span> 和时间 <span class="math inline">\(j\)</span> 对线性预测变量的组合固定效应，<span class="math inline">\(s_{ik}\)</span> 是分配给第 <span class="math inline">\(i\)</span> 个处理的第 <span class="math inline">\(k\)</span> 个个体的随机效应。我们可以根据研究目标以多种方式划分 <span class="math inline">\(\alpha_{ij}\)</span>。最常见的两种是因子效应划分，<span class="math inline">\(\alpha_{ij}=\alpha+\tau_i+\gamma_j+(\tau\gamma)_{ij}\)</span>，其中 <span class="math inline">\(\alpha\)</span> 表示截距，<span class="math inline">\(\tau_i\)</span> 是处理主效应，<span class="math inline">\(\gamma_j\)</span> 是时间主效应，<span class="math inline">\((\tau\gamma)_{ij}\)</span> 是时间与治疗的交互作用；以及多项式回归划分</li>
</ul>
<p><span class="math display">\[\alpha_{ij}=\beta_{0i}+\sum_{m=1}^5\beta_{mi}T_j^m\]</span></p>
<p>其中，<span class="math inline">\(\beta_{0i}\)</span> 表示第 <span class="math inline">\(i\)</span> 个处理的截距，<span class="math inline">\(\beta_{mi}\)</span> 表示第 <span class="math inline">\(i\)</span> 个处理的 <span class="math inline">\(m\)</span> 阶回归系数，<span class="math inline">\(T_{jm}\)</span> 表示第 <span class="math inline">\(j\)</span> 次观测的 <span class="math inline">\(m\)</span> 次方。因此，<span class="math inline">\(m = 1\)</span> 表示线性回归项，<span class="math inline">\(m = 2\)</span> 表示二次项，等等。</p>
<ul>
<li><p>分布：<span class="math inline">\(s_{ik}\text{ iid }N\left(0,\sigma_S^2\right)\)</span>。令 <span class="math inline">\(\symbf{y'}_{ik}\mid s_{ik}=\begin{pmatrix}\begin{bmatrix}y_{i1k}&y_{i2k}&y_{i3k}&y_{i4k}&y_{i5k}&y_{i6k}\end{bmatrix}|s_{ik}\end{pmatrix}\)</span>。<span class="math inline">\(\left(\symbf{y}\mid\symbf{s}\right)\sim N\left(\symbf{\mu}\mid\symbf{s},\boldsymbol\Sigma\right)\)</span> 其中 <span class="math inline">\(\symbf{\mu}=\symbf{1}_{ts}\otimes E\Big(\symbf{y}_{ik}\mid s_{ik}\Big)\)</span>，<span class="math inline">\(\boldsymbol{\Sigma}=\symbf{I}_{ts}\otimes Var\left(\symbf{y}_{ik}\mid s_{ik}\right)\)</span> 以及 <span class="math inline">\(ts\)</span> 等于处理数 <span class="math inline">\(t\)</span> × 个体数 <span class="math inline">\(s\)</span>。</p></li>
<li><p>连接：恒等，<span class="math inline">\(\symbf\eta=\symbf\mu\)</span></p></li>
</ul>
<p>方便起见，将第 <span class="math inline">\(ik\)</span> 个个体对应的协方差阵 <span class="math inline">\(\boldsymbol\Sigma\)</span> 的分块对角元表示为 <span class="math inline">\(\boldsymbol\Sigma_{ik}\)</span>。为什么要对 <span class="math inline">\(y_{ijk}\mid s_{ik}\)</span> 的方差进行所有阐述？为什么它们不像前面例子中的其他高斯线性模型那样 <span class="math inline">\(\operatorname{iid}N\left(\mu_{ijk},\sigma^2\right)\)</span>？因为当我们对同一个体进行一段时间的观测时，我们必须考虑到这些观测具有相关性的可能，也就是说，在给定的个体内，<span class="math inline">\(Cov\left(y_{ijk},y_{ij'k}\mid s_{ik}\right)\)</span> 对于时间 <span class="math inline">\(j \ne j\)</span>，不是必然为零。</p>
<p>我们将协方差阵 <span class="math inline">\(\boldsymbol\Sigma_{ik}\)</span> 的一般形式写为</p>
<p><span class="math display" id="eq:7-3">\[\begin{align}
\boldsymbol\Sigma_{ik}=\begin{bmatrix}\sigma_1^2&\sigma_{12}&\sigma_{13}&\sigma_{14}&\sigma_{15}&\sigma_{16}\\&\sigma_2^2&\sigma_{23}&\sigma_{24}&\sigma_{25}&\sigma_{26}\\&&\sigma_3^2&\sigma_{34}&\sigma_{35}&\sigma_{36}\\&&&\sigma_4^2&\sigma_{45}&\sigma_{46}\\&&&&\sigma_5^2&\sigma_{56}\\&&&&&\sigma_6^2\end{bmatrix}
\tag{7.3}
\end{align}\]</span></p>
<p>其中在时间 <span class="math inline">\(j\)</span> 时的观测之间的方差 <span class="math inline">\(\sigma_j^2=Var\left(y_{ijk}\mid s_{ik}\right)\)</span>，<span class="math inline">\(\sigma_{ij}=Cov\left(y_{ijk},y_{ij'k}\mid S_{ik}\right)\)</span>。这称为<strong>非结构化协方差阵</strong> (unstructured covariance matrix). 潜在地，每次观测都具有独特的方差，并且同一个体每两次的观测都具有独特的协方差。对于 <span class="math inline">\(J\)</span> 次观测，非结构化协方差阵中有 <span class="math inline">\(J+[J(J-1)/2]\)</span> 个协方差参数组成 <span class="math inline">\(\symbf\sigma\)</span>。不管怎么想，这都不是一个简约的模型。</p>
<p>在另一个极端，如果我们假定所有 <span class="math inline">\(y_{ijk}\mid s_{ik}\sim N\left(\mu_{ijk},\sigma^2\right)\)</span>，因此 <span class="math inline">\(\boldsymbol\Sigma_{ik} = \symbf I\sigma^2\)</span>，那么我们就有了一个<strong>时间裂区协方差模型</strong> (split-plot-in-time covariance model)，之所以这么称呼是因为它与裂区实验是相同的模型（我们将在第 <a href="chap9.html#chap9">9</a> 章中更详细地讨论裂区和其他多水平混合模型）。时间裂区模型也称为<strong>独立模型</strong> (independence model). 在第 <a href="chap3.html#chap3">3</a> 章中，我们看到时间裂区或独立模型也可表示为复合对称模型。</p>
<p>非结构化模型和独立/复合对称模型都不适合大多数重复测量分析。就检验假设的统计功效而言，使用非结构化模型通常是浪费且成本高昂的。另一方面，独立/复合对称模型无法解释重复测量之间的重要相关性。当确实存在不可忽略的相关性时，这会导致 I 类错误率过高。我们通常可以找到折衷 (middle ground)，即一种充分考虑相关性的协方差模型，但比非结构化模型更简约。这样做可以让我们完全控制 I 类错误率，而不会不必要地牺牲功效。</p>
<p>我们将在第 <a href="#chap17"><strong>??</strong></a> 章和第 <a href="#chap18"><strong>??</strong></a> 章更详细地讨论协方差模型。目前，两个常见的折衷模型是<strong>一阶事前相关模型</strong> (first-order ante-dependence model) ——通常简称为 <strong>ANTE(1)</strong> ——以及<strong>一阶自回归模型</strong> (first-order auto-regressive model) ——简写为 <strong>AR(1)</strong>. ANTE(1) 的协方差阵为</p>
<p><span class="math display" id="eq:7-4">\[\begin{align}
\boldsymbol{\Sigma}_{{_{tk}}}=\begin{bmatrix}\sigma_{1}^{2}&\sigma_{1}\sigma_{2}\rho_{1}&\sigma_{1}\sigma_{3}\rho_{1}\rho_{2}&\sigma_{1}\sigma_{4}\rho_{1}\rho_{2}\rho_{3}&\sigma_{1}\sigma_{5}\rho_{1}\rho_{2}\rho_{3}\rho_{4}&\sigma_{1}\sigma_{6}\rho_{1}\rho_{2}\rho_{3}\rho_{4}\rho_{5}\\&\sigma_{2}^{2}&\sigma_{2}\sigma_{3}\rho_{2}&\sigma_{2}\sigma_{4}\rho_{2}\rho_{3}&\sigma_{2}\sigma_{5}\rho_{2}\rho_{3}\rho_{4}&\sigma_{2}\sigma_{6}\rho_{2}\rho_{3}\rho_{4}\rho_{5}\\&&\sigma_{3}^{2}&\sigma_{3}\sigma_{4}\rho_{3}&\sigma_{3}\sigma_{5}\rho_{3}\rho_{4}&\sigma_{3}\sigma_{6}\rho_{3}\rho_{4}\rho_{5}\\&&&\sigma_{4}^{2}&\sigma_{4}\sigma_{5}\rho_{4}&\sigma_{4}\sigma_{6}\rho_{4}\rho_{5}\\&&&&\sigma_{5}^{2}&\sigma_{5}\sigma_{6}\rho_{5}\\&&&&&\sigma_{6}^{2}\end{bmatrix}
\tag{7.4}
\end{align}\]</span></p>
<p>其中，<span class="math inline">\(\rho_j\)</span> 表示 <span class="math inline">\(y_{ijk}\mid s_{ik}\)</span> 和 <span class="math inline">\(y_{i,j+1,k}\mid s_{ik}\)</span> 之间的相关性，<span class="math inline">\(j=1,2,\ldots,5\)</span>。同一个体非相邻观测之间的相关性由间隔中包含的 <span class="math inline">\(\rho_j\)</span> 的乘积来建模。在这个例子中，对于六次重复测量，ANTE(1) 允许我们将相对于非结构化模型的协方差参数的数量从 21 个减少到 11 个。</p>
<p>如果我们进一步假定所有方差和所有相邻观测的相关性都相等，即所有 <span class="math inline">\(\sigma_j^2 =\sigma^2\)</span> 且所有 <span class="math inline">\(\rho_j=\rho\)</span>，那么我们就有一个 AR(1) 协方差模型。AR(1) 模型只有 2 个协方差参数。在许多应用中，AR(1) 提供了个体内相关性的适当模型，从而在不牺牲 I 类错误控制的情况下将功效最大化。</p>
<p>对于 Data Set 7.2，竞争模型的 REML -2 倍对数似然值为</p>
<div class="inline-figure"><img src="figure/figure%20217.png#center" style="width:80.0%"></div>
<p>似然比检验用于比较当前模型与下一行更简单的模型。例如，似然比 12.17 检验原假设，即除了 ANTE(1) 模型之外，需要指定非结构化模型所需的额外 10 个协方差参数是否均为零。<span class="math inline">\(p\)</span> 值为 0.2765 意味着我们无法拒绝这一假设——相对于 ANTE(1) 模型，非结构化模型并未增加任何价值，而是浪费且降低了功效。</p>
<p>从表中，我们看到我们未能拒绝所有相邻观测相关性均相等且所有方差均相等的原假设——也就是说，ANTE(1) 模型的额外复杂性不会为这些数据增加任何价值——但我们确实拒绝了 AR(1) 相关系数为零的假设。因此，对于这些数据，我们可以得出结论，个体内存在非平凡的相关性，并且，如果仅限于此处所示的协方差模型选择，则 AR(1) 模型提供了对该相关性进行充分建模的最简约的方法。后续分析应使用 AR(1) 模型。</p>
<p>只要我们满足两个条件，似然比方法就可以很好地工作。首先，我们必须有一个合理的似然：REML 或全似然。我们不能使用伪似然进行似然比检验。其次，比较的模型必须是嵌套的。此处，独立模型是 AR(1) + 随机个体效应模型（<span class="math inline">\(\rho=0\)</span>）的子集，而 AR(1) 模型又是 ANTE(1) 模型（具有如上所示的方差和相关性公式）的子集，ANTE(1) 模型又是非结构化模型的子集。如果模型不是嵌套的，我们就不能使用似然比方法。例如，在重复测量分析中经常使用的无额外 <span class="math inline">\(s_{ik}\)</span> 效应的 AR(1) 模型不能与使用似然比检验的复合对称模型进行比较，因为这两个模型都有 2 个协方差参数。在下一节中，我们讨论伪似然问题。在 <a href="chap7.html#sec7-3">7.3</a> 节中，我们讨论非嵌套模型的比较问题。</p>
</div>
</div>
<div id="sec7-2-4" class="section level3" number="7.2.4">
<h3>
<span class="header-section-number">7.2.4</span> GLMMs 的 PL 与积分近似的结果<a class="anchor" aria-label="anchor" href="#sec7-2-4"><i class="fas fa-link"></i></a>
</h3>
<p>当我们将相关误差模型拟合到高斯数据时，个体内相关性自然成为 <span class="math inline">\(\symbf R\)</span> 矩阵的一部分，回想 <span class="math inline">\(\symbf R = Var(\symbf y\mid \symbf b)\)</span>，给定随机效应的观测的条件方差。对于上一节描述的所有重复测量模型，<span class="math inline">\(\symbf R=\symbf I \otimes \boldsymbol \Sigma_{ik}\)</span>。</p>
<p>当我们将重复测量模型拟合到非高斯数据时，这是如何工作的？对于指数族的非高斯成员和所有拟似然，<span class="math inline">\(Var(\symbf y\mid \symbf b)\)</span> 至少部分依赖于 <span class="math inline">\(E(\symbf y\mid\symbf b )\)</span>。这意味着我们不能像处理高斯数据那样只编写具有相关结构的 <span class="math inline">\(\symbf R\)</span> 矩阵，而忽略方差函数。</p>
<p>解决这个问题有两种方法：“LMM 类比”和“费舍尔会怎么做？”前者得到工作相关矩阵和边际模型，也称为 <strong>R 侧建模</strong> (R-side modeling) 或 <strong>GEE 型建模</strong> (GEE-type modeling). 后者导致通过线性预测器而不是 <span class="math inline">\(Var(\symbf y\mid \symbf b)\)</span> 来对相关性建模。这称为 <strong>G 侧建模</strong> (G-side modeling). 上一节中的 AR(1) + 随机个体效应模型能最好地说明这种差异。</p>
<div id="sec7-2-4-1" class="section level4" number="7.2.4.1">
<h4>
<span class="header-section-number">7.2.4.1</span> R 侧或工作相关模型<a class="anchor" aria-label="anchor" href="#sec7-2-4-1"><i class="fas fa-link"></i></a>
</h4>
<p>对于非高斯 GLMMs，回想，在矩阵形式中，<span class="math inline">\(Var\left(\symbf{y}\mid\symbf{b}\right)=\symbf{V}_\mu^{1/2}\symbf{A}\symbf{V}_\mu^{1/2}\)</span>，其中 <span class="math inline">\(\symbf V_\mu\)</span> 表示方差函数矩阵，<span class="math inline">\(\symbf A\)</span> 表示尺度参数矩阵。我们在第 <a href="chap5.html#chap5">5</a> 章中看到，对 <span class="math inline">\(\symbf y\mid\symbf b\)</span> 条件分布中的协方差建模的一种方法是，用<strong>工作相关性</strong> (working correlation)，又称<strong>工作协方差阵</strong> (working covariance matrix)，替换尺度矩阵。对于 AR(1) 重复测量模型，我们将工作协方差阵写为</p>
<p><span class="math display" id="eq:7-5">\[\begin{align}
\symbf{A}_W=\phi\begin{bmatrix}1&\rho&\rho^2&\rho^3&\rho^4&\rho^5\\&1&\rho&\rho^2&\rho^3&\rho^4\\&&1&\rho&\rho^2&\rho^3\\&&&1&\rho&\rho^2\\&&&&1&\rho\\&&&&&1\end{bmatrix}
\tag{7.5}
\end{align}\]</span></p>
<p>例如，假设我们有二项数据。Data Set 7.3 给出了这样一个数据集，其设计结构与 Data Set 7.2 相同，但具有二项而不是高斯响应变量。用于该数据的 <strong>R 侧 AR(1) GLMM</strong> 为：</p>
<ul>
<li>线性预测器：与高斯模型相同，即 <span class="math inline">\(\eta_{ijk}=\alpha_{ij}+s_{ik}\)</span>，其中 <span class="math inline">\(\alpha_{ij}\)</span> 表示处理-时间分量，可按照高斯模型的描述进行划分。</li>
<li>分布：
<ul>
<li>
<span class="math inline">\(s_{ik}\mathrm{~iid~}N\left(0,\sigma_S^2\right)\)</span>。这称为 G 侧，因为它描述了线性预测器中随机效应的分布，并且一般地，<span class="math inline">\(Var\left(\symbf{b}\right)=\symbf{G}\)</span>。</li>
<li>
<span class="math inline">\(y_{ijk}\mid s_{ik}\sim\mathrm{Binomial}\left(n_{ijk},\pi_{ijk}\right)\)</span>，但对方差进行拟似然校正：<span class="math inline">\(Var\left(\symbf{y}_{ik}\mid s_{ik}\right)=\symbf{V}_\mu^{1/2}\symbf{A}_W\symbf{V}_\mu^{1/2}\)</span>，其中 <span class="math inline">\(\symbf y_{ik}\)</span> 表示第 <span class="math inline">\(ik\)</span> 个个体的 6 次观测的向量（如上面针对高斯模型所述），以及 <span class="math inline">\(\symbf A_W\)</span> 为 <a href="chap7.html#eq:7-5">(7.5)</a> 中给出的工作协方差。因此：
<span class="math display">\[Var\left(\symbf{y}\mid\symbf{s}\right)=\symbf{I}_{ts}\otimes\left(\symbf{V}_\mu^{1/2}\symbf{A}_W\symbf{V}_\mu^{1/2}\right)\]</span>
</li>
</ul>
</li>
<li>连接：<span class="math inline">\(\mathrm{logit}\left(\pi_{ijk}\right)\)</span>
</li>
</ul>
<p>GLMM 从业者通常将此称为“GEE 型”模型。由于该模型在线性预测器中包含随机个体效应，因此它不是真正的 GEE 模型。真正的 GEE 模型是完全边际的：线性预测器仅包含固定效应，即 <span class="math inline">\(\symbf\eta=\symbf{X\beta}\)</span> 并且所有的方差结构都包含于 <span class="math inline">\(\symbf{V}_\mu^{1/2}\symbf{A}_W\symbf{V}_\mu^{1/2}\)</span>。GEE 型模型具有工作协方差阵，但线性预测器中可能包含其他随机效应。</p>
<p>GEE 型模型的估计必须使用某种形式的线性化：第 <a href="chap5.html#chap5">5</a> 章描述的伪似然、Breslow and Clayton (1993) 描述的惩罚拟似然或明确地使用广义估计方程，如 Liang and Zeger (1986) 所述。在本次讨论中，我们使用伪似然法。</p>
<p>我们现在描述 G 侧的替代方法，然后比较并对比这两种方法的模型选择问题。</p>
</div>
<div id="sec7-2-4-2" class="section level4" number="7.2.4.2">
<h4>
<span class="header-section-number">7.2.4.2</span> Fisher 会怎么做？G 侧方法<a class="anchor" aria-label="anchor" href="#sec7-2-4-2"><i class="fas fa-link"></i></a>
</h4>
<p>回想第 <a href="chap2.html#chap2">2</a> 章，模型构建的“费舍尔会怎么做？”过程，包括编写“地形”（或设计）方面和处理方面的框架 ANOVA，然后将它们组合起来。对这些数据执行此操作会得到：</p>
<div class="inline-figure"><img src="figure/figure%20219.png#center" style="width:100.0%"></div>
<p>正如我们在第 <a href="chap3.html#chap3">3</a> 章中看到的，当我们使用组合框架 ANOVA 来构建高斯数据的模型时，我们必须将最后一行解释为“残差”，以解释高斯分布的尺度参数 <span class="math inline">\(\sigma^2\)</span>，或者在更复杂的情况下，如重复测量，<span class="math inline">\(\symbf y\mid\symbf b\)</span> 的协方差阵。然而，我们不需要对指数族的单参数成员进行此操作，例如本例中的二项。</p>
<p>这就提出了一个问题：二项数据重复测量之间的相关性出现在何处？从概念上讲，在二项 GLMM 中，我们假定 <span class="math inline">\(\pi_{ijk}\)</span> 因处理和时间效应而系统地变异，并且它会受到个体间和个体内变异的随机扰动。同一个体重复测量之间的相关性必须是后者的一部分：个体内 <span class="math inline">\(\pi_{ijk}\)</span> 的随机扰动。换言之：在应用连接之前对随机扰动 <span class="math inline">\(\pi_{ijk}\)</span> 中的相关性进行建模更有意义，还是在应用连接之后通过工作相关性进行建模更有意义。前者要求个体内变异是线性预测器的一部分，而且给予我们一个真实的似然。后者给了我们一个拟似然和一个不对应于任何可行的概率机制的过程。</p>
<p>“Fisher 会怎么做”给了我们如下的 GLMM：</p>
<ul>
<li><p>线性预测器：<span class="math inline">\(\eta_{ijk}=\alpha_{ij}+s_{ik}+w_{ijk}\)</span></p></li>
<li><p>分布：<span class="math inline">\(y_{ijk}\mid s_{ik},w_{ijk}\sim\mathrm{Binomial}\left(n_{ijk},\pi_{ijk}\right)\)</span>，<span class="math inline">\(s_{ik}\mathrm{~iid~}N\left(0,\sigma_s^2\right)\)</span>，<span class="math inline">\(\mathbf{w}_{ik}\sim N\left(\mathbf{0},\boldsymbol\Sigma_{ik}\right)\)</span>，其中 <span class="math inline">\(\symbf{w}'_{ik}=\begin{bmatrix}w_{i1k}&w_{i2k}&w_{i3k}&w_{i4k}&w_{i5k}&w_{i6k}\end{bmatrix}\)</span> 以及
<span class="math display" id="eq:7-6">\[\begin{align}
\boldsymbol\Sigma_{ik}=\sigma_W^2\begin{bmatrix}1&\rho&\rho^2&\rho^3&\rho^4&\rho^5\\&1&\rho&\rho^2&\rho^3&\rho^4\\&&1&\rho&\rho^2&\rho^3\\&&&1&\rho&\rho^2\\&&&&1&\rho\\&&&&&1\end{bmatrix}
\tag{7.6}
\end{align}\]</span></p></li>
<li><p>连接：<span class="math inline">\(\eta_{ijk}=\operatorname{logit}\left(\pi_{ijk}\right)\)</span></p></li>
</ul>
<p>这是一个 G 侧模型，因为除了 <span class="math inline">\(y_{ijk}\mid s_{ik},w_{ijk}\)</span> 的条件分布之外，所有的随机变异性都出现在线性预测器中。请注意，该模型中的所有分布都是真实的概率分布。这意味着可以使用伪似然或积分近似来进行估计和推断。这对于 GLMMs 的协方差模型选择具有重要意义。</p>
</div>
<div id="sec7-2-4-3" class="section level4" number="7.2.4.3">
<h4>
<span class="header-section-number">7.2.4.3</span> R 侧与 G 侧：协方差模型选择的后果<a class="anchor" aria-label="anchor" href="#sec7-2-4-3"><i class="fas fa-link"></i></a>
</h4>
<p>对于具有工作相关结构的非高斯模型，我们必须使用伪似然估计。对于仅 G 侧的非高斯模型，我们可以使用伪似然或积分近似。为什么这很重要？我们可以通过检查似然比统计量来回答这个问题。</p>
<p>我们知道，在实践中，我们将检验协方差分量的似然比统计量计算为 <span class="math inline">\(\chi_{LR}^2=2\ell\left(\hat{\symbf\sigma}\right)-2\ell\left(\hat{\symbf\sigma}_0\right)\)</span>，其中 <span class="math inline">\(\ell\left(\hat{\symbf\sigma}\right)\)</span> 表示在 <span class="math inline">\(\hat{\symbf\sigma}\)</span> 处计算的对数似然。使用积分近似，作为估计过程的一部分我们可获得 <span class="math inline">\(\ell\left(\hat{\symbf\sigma}\right)\)</span> 。但是使用伪似然时，我们从不计算 <span class="math inline">\(\ell\left(\hat{\symbf\sigma}\right)\)</span>。我们构建了一个伪似然，因此我们最接近 <span class="math inline">\(\ell\left(\hat{\symbf\sigma}\right)\)</span> 的是 <span class="math inline">\(p\ell\left(\hat{\symbf\sigma}\right)\)</span>，即在 <span class="math inline">\(\hat{\symbf\sigma}\)</span> 处计算的伪似然。使用此方法进行似然比检验的问题在于，全模型伪似然所基于的伪变量 <span class="math inline">\(\symbf y^*\)</span> 与简化模型所基于的伪变量不同，因此伪似然也不同。虽然我们可以计算一个伪似然比 <span class="math inline">\(2p\ell\left(\hat{\sigma}\right)-2p\ell\left(\hat{\sigma}_0\right)\)</span>，但由于全模型下的伪变量 <span class="math inline">\(\symbf y^*_0\)</span> 和原假设（即简化模型）下的伪变量 <span class="math inline">\(\symbf y_0\)</span> 不具有可比性，因此它的解释存疑
。</p>
<p>对于 Data Set 7.3，下表显示了 R 侧模型的伪对数似然以及使用拉普拉斯法计算的、具有我们在高斯示例（Data Set 7.2）中考虑的协方差结构的 G 侧模型的对数似然。请注意，R 侧复合对称 (CS) 和独立模型并不等价（正如我们在第 <a href="chap3.html#chap3">3</a> 章中看到的），因此每个模型都有不同的伪似然。尽管表中没有明确说明，但所有对数似然和伪对数似然统计量均已乘以 -2.</p>
<div class="inline-figure"><img src="figure/figure%20221.png#center" style="width:90.0%"></div>
<p>如果我们使用伪似然统计量，检验非结构化协方差模型与 ANTE(1) 协方差模型的 <span class="math inline">\(p\)</span> 值为 0.0008，这（虚假地）表明我们必须使用非结构化工作协方差模型来充分解释个体内相关。G 侧方法给出了更准确的刻画：带有额外随机个体效应 <span class="math inline">\(s_{ik}\)</span> 的 AR(1) 模型是足够的。</p>
<div class="rmdcaution">
<p>我们如何知道 R 侧结果是虚假的？其实 Data Set 7.3 是模拟的。该数据是使用 AR(1) + 随机个体模型生成的。这是一个典型的模拟数据集：R 测伪似然比检验未能可靠地选择 AR(1) + 随机个体模型；而 G 侧似然比检验做到了。</p>
</div>
<p>总之，对于协方差分量的正式检验，似然比检验是我们可使用的最通用 (versatile) 且最普遍适用 (generally applicable) 的工具。对于高斯数据，使用 REML 对数似然来构造似然比统计量。对于非高斯数据，我们必须使用一种允许计算真实似然的估计方法，例如通过拉普拉斯或高斯-厄尔米特求积的积分近似。这将非高斯 GLMM 的正式协方差分量检验限制为仅 G 侧模型，即没有工作相关结构的模型。对于许多 GLMMs 来说，我们可以在 G 侧模型和 R 侧、GEE 型模型之间进行选择。如果优先考虑正式的协方差检验，请尽可能使用 G 侧模型。</p>
<p>上一段不应被解读为对 R 侧、GEE 型建模的全面谴责。有许多情况是 GEE 型模型是唯一可行的方法，而在其他情况下，GEE 是首选方法。我们将在第 <a href="#chap11"><strong>??</strong></a>、<a href="#chap12"><strong>??</strong></a> 章（过度分散）和第 <a href="#chap17"><strong>??</strong></a> 章（非高斯重复测量）中进一步讨论这个问题。</p>
<p>请注意，我们讨论的所有检验都要求固定效应 <span class="math inline">\(\symbf{X\beta}\)</span> 保持不变。在原假设下，全模型和简化模型之间唯一允许变化的是模型的随机部分，对于 LMM 为 <span class="math inline">\(\symbf {Zb},\symbf G\)</span> 和 <span class="math inline">\(\symbf R\)</span>，对于 GLMMs 为 <span class="math inline">\(\symbf {Zb},\symbf G\)</span> 或 <span class="math inline">\(\symbf A_W\)</span>。</p>
<p>现在，我们将注意力转向协方差模型的比较，对于这种比较，正式的检验要么不可能，要么根本不合适。</p>
</div>
</div>
</div>
<div id="sec7-3" class="section level2" number="7.3">
<h2>
<span class="header-section-number">7.3</span> 用拟和统计量比较协方差模型<a class="anchor" aria-label="anchor" href="#sec7-3"><i class="fas fa-link"></i></a>
</h2>
<p>假设我们想了解复合对称模型的拟合与没有额外个体随机效应的 AR(1) 模型的拟合相比如何。两个模型都有 2 个协方差参数—— <span class="math inline">\(\sigma^2\)</span> 和 <span class="math inline">\(\rho\)</span>，尽管我们在这两个模型中对这些参数的解释不同。两者都不是另一个的子集。我们无法构建正式的检验。</p>
<p>即使可能进行正式检验，许多人认为，对于许多应用来说，假设检验并不是识别适当、简约模型的关键。一些人认为检验是不合适的；另一些人只是表示，虽然检验没有错，但它也不能提供信息。有关这些问题的更广泛讨论，见 Burnham and Anderson (2002) 和 Konishi and Kitagawa( 2008).</p>
<p>信息准则提供了正式检验的替代方案。目前应用最广泛的两种信息准则是 AIC-Akaike(1974) 信息准则和 Schwarz (1978) 提出的 BIC-Bayes 信息准则。Burnham and Anderson (1998) 建议对 Akaike 准则进行小样本校正——校正的信息准则缩写为 AICC. 本节简要概述了信息准则背后的理由。然后，我们使用 Data Set 7.2 和 7.3 更新示例以说明其应用。</p>
<div id="sec7-3-1" class="section level3" number="7.3.1">
<h3>
<span class="header-section-number">7.3.1</span> AIC 和 AICC<a class="anchor" aria-label="anchor" href="#sec7-3-1"><i class="fas fa-link"></i></a>
</h3>
<p>AIC 的根源可追溯到 Kullbach and Leibler (1951)，他们提出了两个模型之间差异的度量，他们用函数 <span class="math inline">\(f\)</span> 表示“真相” (“truth”)，用函数 <span class="math inline">\(g\)</span> 表示“近似模型” (“approximating model). 按照他们的符号，“真相”只取决于数据，由此得到函数 <span class="math inline">\(f(\symbf y)\)</span>，而近似模型取决于模型和参数，因此表示为 <span class="math inline">\(g(\symbf y\mid\symbf\theta)\)</span>。对于连续数据，Kullback and Leiber 首先使用了一个被 Konishi and Kitagawa 称为 <strong>K-L 信息</strong>的量，记为 <span class="math inline">\(I(f,g)\)</span>，定义为：</p>
<p><span class="math display">\[E_f\left[\log\!\left(\frac{f\left(\symbf{y}\right)}{g\left(\symbf{y}\mid\symbf{\theta}\right)}\right)\right]\]</span></p>
<p>请注意，K-L 信息在概念上类似于似然比统计量。此外，完美的近似模型将产生 <span class="math inline">\(\frac{f\left(\symbf y\right)}{g\left(\symbf y \mid \symbf\theta\right)}=1\)</span> 的比值，因此对数比等于 0。我们可将 <span class="math inline">\(I(f,g)\)</span> 重写为 <span class="math inline">\(E_f\left\{\log\Bigl[f\left(\symbf{y}\right)\Bigr]\right\}-E_f\left\{\log\Bigl[g\bigl(\symbf{y}\left|\symbf{\theta}\right)\Bigr]\right\}\)</span>。令 <span class="math inline">\(E_f\left\{\log\!\left[f\left(\symbf{y}\right)\right]\right\}=C\)</span>，其中 <span class="math inline">\(C\)</span> 为满足 <span class="math inline">\(-E_f\left\{\log\Bigl[g\bigl(\symbf{y}\mid\symbf{\theta}\bigr)\Bigr]\right\}=I\bigl(f,g\bigr)-C\)</span> 的常数。式 <span class="math inline">\(I(f,g)-C\)</span> 称为 <strong>K-L 距离</strong>。</p>
<p>Akaike 表明 <span class="math inline">\(-2\ell\left(\hat{\symbf{\theta}}\right)+2p\)</span>（其中 <span class="math inline">\(p\)</span> 表示参数向量 <span class="math inline">\(\symbf\theta\)</span> 中的参数数量）提供了 K-L 距离的良好近似。<strong>AIC</strong> (Akaike information criterion) 定义为</p>
<p><span class="math display" id="eq:7-7">\[\begin{align}
AIC=-2\ell\left(\hat{\symbf{\theta}}\right)+2p
\tag{7.7}
\end{align}\]</span></p>
<p>注意，我们可将 AIC 视为 -2 乘以对数似然值减去一个等于模型参数数量的惩罚项。我们知道，如果我们向模型中添加额外的参数，我们总是可以增加对数似然值，直到 <span class="math inline">\(\ell(\symbf y)\)</span> 给出最大的可能值。这就是我们在第 <a href="chap6.html#chap6">6</a> 章中讨论的偏差统计量的基本思想。然而，我们会遇到一个边际收益递减的点，并且拟合模型的整体目的是为数据的变异或表现提供一个简约的解释。AIC 对向模型中添加新参数进行惩罚；除非对数似然值的增加超过了添加到模型中的参数数量，否则我们认为添加这些参数是适得其反的。</p>
<p>Burnham and Anderson (1998) 建议进行小样本校正。令</p>
<p><span class="math display">\[n^*=\begin{cases}N\;\text{ 若 }\;\ell(\mathbf{\theta})\;\text{ 为全似然}\\N-\operatorname{rank}(\symbf{X})\;\text{ 若 }\ell(\symbf{\theta})\text{ 为}\;\text{REML 似然}\end{cases}\]</span></p>
<p>其中 <span class="math inline">\(N\)</span> 表示观测总数。小样本校正为</p>
<p><span class="math display">\[\frac{2p(p+1)}{n^*-p-1}\]</span></p>
<p>为 AIC 加上校正项得到小样本调整，或“校正” AIC：</p>
<p><span class="math display" id="eq:7-8">\[\begin{align}
AICC=AIC+\frac{2p\left(p+1\right)}{n^*-p-1}=-2\ell\left(\symbf \theta\right)+\frac{2pn^*}{n^*-p-1}
\tag{7.8}
\end{align}\]</span></p>
</div>
<div id="sec7-3-2" class="section level3" number="7.3.2">
<h3>
<span class="header-section-number">7.3.2</span> BIC<a class="anchor" aria-label="anchor" href="#sec7-3-2"><i class="fas fa-link"></i></a>
</h3>
<p>Schwarz (1978) 提出了一种使用贝叶斯观点的替代方法。我们有许多相互竞争的模型，其一被认为是“正确的”模式。对于每个相互竞争的模型，你都有一个它是“正确的”模型的先验概率。假定所有竞争模型是等可能的，Schwarz 推导出了<strong>贝叶斯信息准则</strong> (Bayesian information
criterion, BIC)：</p>
<p><span class="math display">\[BIC=-2\ell\left(\hat{\symbf{\theta}}\right)+p\log s\]</span></p>
<p>其中 <span class="math inline">\(s\)</span> 表示个体数。</p>
</div>
<div id="sec7-3-3" class="section level3" number="7.3.3">
<h3>
<span class="header-section-number">7.3.3</span> 协方差模型比较的应用<a class="anchor" aria-label="anchor" href="#sec7-3-3"><i class="fas fa-link"></i></a>
</h3>
<p>对于使用 REML 估计的 LMMs，由于固定效应参数不出现在 REML 似然中，所以对数似然的参数向量为 <span class="math inline">\(\symbf\sigma\)</span>。参数数量等于协方差参数的数量。保持设计尺寸值 <span class="math inline">\(n^*\)</span> 和 <span class="math inline">\(s\)</span> 不变。为了使用 REML 评估 LMM 协方差模型，三个信息准则可写为</p>
<p><span class="math display" id="eq:7-9">\[\begin{align}
&\text{AIC:}-2\ell_R\left(\hat{\symbf\sigma}\right)+2p\notag\\
&\text{AICC: AIC}+\frac{2p\left(p+1\right)}{n^*-p-1}=-2\ell_R\left(\hat{\symbf{\sigma}}\right)+\frac{2pn^*}{n^*-p-1}
\tag{7.9}
\\
&\text{BIC: }-2\ell_R\left(\hat{\symbf{\sigma}}\right)+p\log\left(s\right)\notag
\end{align}\]</span></p>
<p>其中 <span class="math inline">\(\ell_R\left(\hat{\symbf{\sigma}}\right)\)</span> 表示在 <span class="math inline">\(\hat{\symbf\sigma}\)</span> 处计算的 REML 对数似然。</p>
<p>对于使用积分近似的 GLMMs，参数向量包括协方差参数和固定效应。如果我们令 <span class="math inline">\(p\)</span> 等于协方差分量的数量——以保持与 <a href="chap7.html#eq:7-9">(7.9)</a> 的一致性——并回想固定效应参数的数量等于 <span class="math inline">\(\operatorname{rank} (\symbf X)\)</span>，我们可使用积分近似将 GLMMs 的信息准则写为</p>
<p><span class="math display" id="eq:7-10">\[\begin{align}
&\text{AIC: }-2\ell\left(\hat{\symbf\sigma},\hat{\symbf\beta}\right)+2\left[p+\operatorname{rank}\left(\symbf{X}\right)\right]\notag\\
&\text{AICC: }-2\ell\left(\hat{\symbf{\sigma}},\hat{\symbf{\beta}}\right)+\frac{2\left[p+\operatorname{rank}\left(\symbf{X}\right)\right]n^*}{n^*-\left[p+\operatorname{rank}\left(\symbf{X}\right)\right]-1}
\tag{7.10}
\\&\text{BIC: }-2\ell\Big(\hat{\symbf\sigma},\hat{\symbf{\beta}}\Big)+\Big[p+\operatorname{rank}(\symbf{X})\Big]\text{log}(s)\notag
\end{align}\]</span></p>
<p>其中 <span class="math inline">\(\ell\Big(\hat{\symbf\sigma},\hat{\symbf{\beta}}\Big)\)</span> 表示在参数估计处计算的全对数似然。同样，如果我们保持模型的固定效应部分不变，那么竞争模型的信息准则之间的差异仅由协方差参数导致，从而使我们能够比较各种协方差模型。</p>
<p>在 <a href="chap7.html#sec7-2">7.2</a> 节中，我们讨论了使用似然比检验来比较协方差模型的两个例子：使用 REML 的高斯数据的 Data Set 7.2 和通过拉普拉斯法使用积分近似的非高斯数据 Data Set 7.3. 继续这些示例，让我们看看信息准则的应用。</p>
<p>首先，我们考虑高斯数据。基于 REML 的信息准则为：</p>
<div class="inline-figure"><img src="figure/figure%20225-1.png#center" style="width:70.0%"></div>
<p>回想 AIC 的起源，它是对 K-L 距离（即近似模型与“真相”之间的差异）的近似—— AIC 值越小越好。我们也以类似的方式解释其他信息准则。根据所有三个准则，“获胜”模型是 AR(1). 虽然 ANTE(1) 和非结构化模型都具有更大的对数似然值，但该改进并未抵消协方差参数数量的增加。我们可以看到信息准则的吸引力——我们可以快速扫描所有竞争模型的信息并做出选择。请注意，随着协方差参数数量的增加，AICC 相对于 AIC 的校正值会大幅增加。Burnham and Anderson 建议，AICC 是这三个信息准则中最可靠的，可作为这些数据的协方差模型选择的指南。</p>
<p>现在，二项 G 侧重复测量 GLMM 的基于全似然的信息准则为：</p>
<div class="inline-figure"><img src="figure/figure%20225-2.png#center" style="width:70.0%"></div>
<p>同样，结果与我们使用似然比检验发现的结果一致。AR(1) 模型提供了最简约的拟合——我们可将非结构化和 ANTE(1) 模型膨胀的信息准则解释为它们对协方差结构“过度建模” (“over-model”) 的证据，以及将复合对称/独立模型膨胀的信息准则解释为对协方差结构“建模不足” (“under-model”) 的证据。过度建模会损害功效，而建模不足则无法控制 I 类错误。</p>
</div>
</div>
<div id="sec7-4" class="section level2" number="7.4">
<h2>
<span class="header-section-number">7.4</span> 区间估计<a class="anchor" aria-label="anchor" href="#sec7-4"><i class="fas fa-link"></i></a>
</h2>
<p>当研究的目标针对协方差分量本身时，我们可能对这些分量的区间估计更感兴趣，而不是检验关于它们的假设。我们简要地讨论一下这个问题。有兴趣的读者可以在专门讨论方差分量的出版物中找到对此主题的更深入讨论，例如 Searle et al. (1992). 这里我们考虑三种构建置信区间的方法。第一个基于标准高斯 <span class="math inline">\(N(0,1)\)</span> 分布，我们提到它主要是因为许多人出于无知而使用它（但不应使用），学习线性模型的学生应该知道这一点。第二个是基于 <span class="math inline">\(\chi^2\)</span> 近似的 Wald 程序。第三个基于轮廓似然 (profile likelihood).</p>
<div id="首先不该做什么" class="section level3 unnumbered">
<h3>首先：不该做什么<a class="anchor" aria-label="anchor" href="#%E9%A6%96%E5%85%88%E4%B8%8D%E8%AF%A5%E5%81%9A%E4%BB%80%E4%B9%88"><i class="fas fa-link"></i></a>
</h3>
<p>之前，我们看到一种永远不应该用于协方差分量的检验统计量是使用正态近似的 Wald 统计量，即对于检验 <span class="math inline">\(H_0\colon\sigma_i=0\)</span>，</p>
<p><span class="math display">\[Wald_{\sigma_i}=\frac{\hat{\sigma}_i}{\sqrt{V_{A_{_{ii}}}}}\]</span></p>
<p>其中 <span class="math inline">\(A_{_{ii}}\)</span> 为 <span class="math inline">\(\hat\sigma_i\)</span> 的渐近方差。它基于假定 <span class="math display">\[Z=\frac{\hat{{\sigma}}_i-{\sigma}_i}{\sqrt{V_{A,ii}}}\]</span> 具有近似 <span class="math inline">\(N(0,1)\)</span> 分布。我们直言不讳：“任何受过教育的数据分析师不应使用 Wald 协方差分量检验。”除非样本量很大，否则正态近似非常差。</p>
<p>与此相关的是根据正态近似导出的置信区间</p>
<p><span class="math display">\[\hat{\sigma}_i\pm{Z}_{\alpha_2}\sqrt{V_{A,ii}}\]</span></p>
<p>其中 <span class="math inline">\((1− \alpha)\)</span> 表示置信水平。同样，对于小于数万或数十万的样本量，正态近似非常差。我们强烈反对这种构建置信区间的方法。重复一遍：任何受过教育的数据分析师都不应该以这种方式构建置信区间。</p>
</div>
<div id="sec7-4-1" class="section level3" number="7.4.1">
<h3>
<span class="header-section-number">7.4.1</span> 基于 <span class="math inline">\(\chi^2\)</span> 的 Wald 法<a class="anchor" aria-label="anchor" href="#sec7-4-1"><i class="fas fa-link"></i></a>
</h3>
<p>这种方法只是对在入门统计课程中教授的关于方差的标准区间估计的一个详细说明。令 <span class="math inline">\(\sigma_i\)</span> 为感兴趣的协方差分量。其估计 <span class="math inline">\(\hat \sigma_i\)</span> 具有近似分布</p>
<p><span class="math display">\[\frac{\sigma_i\chi_\nu^2}\nu\]</span></p>
<p>经计算得到</p>
<p><span class="math display">\[\frac{\nu\times\hat{\sigma}_i}{\chi_{\nu,1-\alpha/2}^2}<\sigma_i<\frac{\nu\times\hat{\sigma}_i}{\chi_{\nu,\alpha/2}^2}\]</span></p>
<p>我们用它来构建置信区间，其中 <span class="math inline">\((1− \alpha)\)</span> 表示置信水平。</p>
</div>
<div id="sec7-4-2" class="section level3" number="7.4.2">
<h3>
<span class="header-section-number">7.4.2</span> 基于似然的方法<a class="anchor" aria-label="anchor" href="#sec7-4-2"><i class="fas fa-link"></i></a>
</h3>
<p>此程序基于似然比统计量 <span class="math inline">\(\chi_{LR}^2=2\ell\left(\hat{\symbf\sigma}\right)-2\ell\left(\hat{\symbf\sigma}_0\right)\)</span>。令 <span class="math inline">\(\sigma_{0,i}\)</span> 为 <span class="math inline">\(H_0\)</span> 下 <span class="math inline">\(\sigma_i\)</span> 置信区间的感兴趣参数。我们可将 <span class="math inline">\(\ell(\symbf\sigma_0)\)</span> 重写为 <span class="math inline">\(\ell(\symbf\sigma_{0,i},\hat{\symbf\sigma}_{i'})\)</span>，其中 <span class="math inline">\(\symbf\sigma_{i'}\)</span> 定义为移除 <span class="math inline">\(\sigma_i\)</span> 的协方差参数向量，并且 <span class="math inline">\(\hat{\symbf\sigma}_{i'}\)</span> 是其最大似然或 REML 估计，具体取决于我们用于估计协方差分量的程序。遵循 Patiwan (2001)，我们将此函数称为<strong>轮廓似然</strong> (profile likelihood).</p>
<p>定义</p>
<p><span class="math display">\[\ell\left(\sigma_{i^{\prime}}\mid\tilde{\symbf\sigma}\right)=\sup_{\sigma_{i'}}\left[\ell\left(\tilde{\sigma}\mid\symbf\sigma_{i^{\prime}}\right)\right]\]</span></p>
<p>其中我们将 <span class="math inline">\(\sigma_i\)</span> 保持为常数 <span class="math inline">\(\tilde\sigma\)</span>。<span class="math inline">\(\sigma_i\)</span> 的轮廓似然置信区间是所有满足 <span class="math inline">\(2\bigg[\ell\big(\hat{\symbf\sigma}\big)-\ell\big(\symbf\sigma\big|\tilde{\sigma}\big)\bigg]<\chi_{1,1-\alpha}^2\)</span> 的 <span class="math inline">\(\sigma_i\)</span> 的集合，其中 <span class="math inline">\((1− \alpha)\)</span> 表示置信水平。</p>
<p>SAS<sup>®</sup> PROC GLIMMIX 提供了计算协方差分量的基于 <span class="math inline">\(\chi^2\)</span> 的 Wald 置信区间和轮廓似然置信区间的选项，如果需要，还提供了 GEE 型模型中工作协方差参数的选项。对于这些模型，轮廓似然使用伪对数似然。有关更多详细信息，请参阅 SAS/STAT/PROC GLIMMIX Version 9.22 On-Line Documentation (2023).</p>
</div>
</div>
<div id="sec7-5" class="section level2" number="7.5">
<h2>
<span class="header-section-number">7.5</span> 总结<a class="anchor" aria-label="anchor" href="#sec7-5"><i class="fas fa-link"></i></a>
</h2>
<p>术语“协方差分量”统指 GLMM 的 <span class="math inline">\(\symbf G\)</span> 和 <span class="math inline">\(\symbf R\)</span> 矩阵中的方差、协方差和相关参数。出于符号目的，它们包含第 <a href="chap6.html#chap6">6</a> 章和第 <a href="chap7.html#chap7">7</a> 章中经常提到的向量 <span class="math inline">\(\symbf\sigma\)</span>。</p>
<p>对协方差分量的推断有三种基本形式：</p>
<ul>
<li>关于 <span class="math inline">\(\symbf\sigma\)</span> 的一个或多个分量的正式推断
<ul>
<li>GLMM 主要工具：似然比检验</li>
<li>可能的原因包括
<ul>
<li>验证假定，例如方差齐性</li>
<li>协方差分量本身是内在兴趣</li>
<li>在竞争的协方差模型中进行选择，例如重复测量分析</li>
<li>后者<a class="footnote-ref" tabindex="0" data-toggle="popover" data-content="<p>译者注：这里的“后者”是指在“竞争的协方差模型中进行选择”。</p>"><sup>22</sup></a>仅在模型嵌套时才可能进行</li>
</ul>
</li>
</ul>
</li>
<li>非正式推断
<ul>
<li>使用拟合统计量，例如 AIC, AICC, BIC</li>
<li>拟合统计量 = 对数似然 - 惩罚</li>
<li>用于在竞争的协方差结构中进行选择</li>
<li>当结构没有嵌套时特别有用</li>
</ul>
</li>
<li>区间估计
<ul>
<li>Wald（简单但不推荐）</li>
<li>基于似然的——轮廓似然</li>
</ul>
</li>
</ul>
<p>对协方差分量的推断需要保持模型的固定效应部分 <span class="math inline">\(\symbf{X\beta}\)</span> 不变。否则，检验会混淆 <span class="math inline">\(\symbf\beta\)</span> 和 <span class="math inline">\(\symbf\sigma\)</span>，因此无效。</p>
<p>对于高斯数据：</p>
<ul>
<li>使用受限似然</li>
<li>似然比检验使用 REML 似然</li>
<li>拟合统计量是惩罚 REML 似然</li>
<li>区间估计为轮廓 REML 似然</li>
</ul>
<p>对于非高斯数据：</p>
<ul>
<li>推断需要条件模型和积分近似（拉普拉斯或求积）。基于伪似然的协方差推断是未定义且无效的。</li>
<li>根据定义，边际模型是一种拟似然模型。在 GLIMMIX 中，这意味着仅通过伪似然进行估计。对于伪似然，拟合统计量未定义。检验需要计算 <span class="math inline">\(H_0\)</span> 和 <span class="math inline">\(H_A\)</span> 下的伪似然。但 <span class="math inline">\(H_0\)</span> 和 <span class="math inline">\(H_A\)</span> 各自隐含着不同的伪变量，因此产生的伪似然是不可比的。</li>
<li>同样的问题也适用于使用条件模型的伪似然估计。</li>
<li>这意味着刚接触 GLMM 的学生和数据分析师要小心：伪似然是默认的 GLIMMIX 算法——你必须使用 <code>METHOD=LAPLACE</code> 或 <code>METHOD=QUADRATURE</code> 覆盖它，才能对非高斯数据的协方差分量进行推断。</li>
</ul>
</div>
<div id="exe7" class="section level2 unnumbered">
<h2>练习<a class="anchor" aria-label="anchor" href="#exe7"><i class="fas fa-link"></i></a>
</h2>
<ol style="list-style-type: decimal">
<li>本习题涉及 SAS Data and Program Library 中的文件 ch7_prob1.sas，其中给出了数据以及两个 GLIMMIX 运行。数据来自一项涉及三个因素的研究：A、B 和存储时间。因素 A 和 B 分别是制造工艺和精加工过程，各有两个水平。由每种 A × B 组合制造的产品批次被置于存储单元中——每种 A × B 组合有六个存储单元。在存储 1 周、2 周、4 周、8 周和 16 周后，从每批产品中抽取样本以检查产品的强度。由于抽样是破坏性的，所以可以假定在不同周数的测量是相互独立的。同时，也可以假定“强度”具有高斯分布。有两个 GLIMMIX 程序，分别称为 “/*Run1*/” 和 “/*Run2*/”.
<ol style="list-style-type: lower-alpha">
<li>给出这两个运行的模型的完整描述。</li>
<li>Run 1 用于估计方差分量的默认方法的名称是什么？</li>
<li>就混合模型理论而言，这两次运行的主要区别是什么？</li>
<li>首选哪种方法（即你应该报告哪种方法）？简单解释一下。</li>
<li>正式地，根据你在 a. 中给出的模型参数，“storage_time*a*b” 在这些运行中检验的是什么？</li>
<li>对于该项目研究团队的非统计学家，你会如何解释 “storage_time*a*b” 检验及其结果？提示：在 “Run 3” 中提供了一个交互图来帮助你回答这个问题。</li>
</ol>
</li>
<li>本习题涉及文件 ch7_prob2.sas 中给出的数据和两个 GLIMMIX 运行。数据来自不完全区组设计，具有七种处理和七个区组。假定数据是高斯的。
<ol style="list-style-type: lower-alpha">
<li>使用 GLIMMIX <code>COVTEST</code> 选项来检验在给定区组效应的条件下，观测条件分布的方差齐性假设（即处理间方差相等）。</li>
<li>根据 a. 中的结果，为这些数据编写适当的 LMM.</li>
<li>使用 b.，你对处理均值响应之间的差异得出了什么结论？假定 trt 5 为“对照”或参考处理。举出适当的证据。</li>
</ol>
</li>
<li>本习题涉及文件 ch7_prob3.sas 中给出的数据和数据步语句。数据来自具有三种处理和六个区组的随机完全区组设计。响应变量 <span class="math inline">\(Y\)</span> 是计数。
<ol style="list-style-type: lower-alpha">
<li>运行标准随机区组效应 LMM 分析。记下方差分量估计、处理的 <span class="math inline">\(F\)</span> 检验以及处理均值的估计（和标准误）。</li>
<li>使用对数计数（一种标准的 GLMM 之前的方差稳定变换）重新运行分析。请注意，“log” 是指自然对数。对于零计数，你可以在获取对数之前向计数添加一个小增量 <span class="math inline">\(\varepsilon > 0\)</span> ——例如 SAS 语句 <code>if log=0 then log_y=log(y+0.1)</code>。获取方差分量估计、处理的 <span class="math inline">\(F\)</span> 值以及数据尺度上的估计值（和标准误）或处理均值。</li>
<li>假设给定区组的计数具有泊松分布，以 GLMM 形式重新运行分析。使用与 a. 和 b. 相同的线性预测器。获取方差分量估计、处理的 <span class="math inline">\(F\)</span> 值以及数据尺度上的估计值（和标准误）或处理均值。</li>
<li>根据 c. 获取 GLMM 的 Pearson 卡方/DF 拟合统计量。</li>
<li>重新运行 c. 中的 GLMM，但向线性预测器添加随机 block*trt 效应，并获得 <span class="math inline">\(H_0\colon\sigma_{BT}^2 = 0\)</span> 的似然比检验。d. 中的统计量与 e. 的检验具有本质上相同的解释。简单地解释一下。</li>
<li>描述 c. 和 e. 之间的线性预测器的变化对 <span class="math inline">\(F\)</span> 检验和数据尺度上处理均值（和标准误）的估计的影响。</li>
<li>选择最适合报告的 GLMM（c. 或 e.），并对 GLMM、对数变换 LMM 和未变换 LMM 的方差分量、<span class="math inline">\(F\)</span> 以及数据尺度的估计结果进行比较并对比。</li>
</ol>
</li>
<li>本习题涉及文件 ch7_prob4.sas 中给出的数据和数据步语句。这是一组重复测量的数据。两种治疗（在数据集中编码为 “0” 和 “1”）被随机分配给每个个体（在数据集中表示为 ID。注意 trt “0” 的 ID=1 与 trt 的 ID=1 不同 “1” ——它们是嵌套的，而不是交叉的）。从 WEEK=0 开始每周观察每个受试者。WEEK=0 时的观测不是基线测量。假定这是在对个体进行治疗之后进行的第一次测量。还假定响应变量是高斯的。
<ol style="list-style-type: lower-alpha">
<li>完成下表以获取指定的合理的协方差模型。假定 week 是一个 CLASS 变量。根据此表中包含的信息回答 b. 和 c.
<img src="figure/figure%20230.png#center" style="width:80.0%">
</li>
<li>协方差模型复杂性的增加对治疗与周的交互作用的检验有何影响？Kenward-Roger 调整对此有何影响？</li>
<li>验证（在这些模型中）ANTE(1) 协方差模型是否提供最佳拟合。为什么？</li>
<li>对于 ANTE(1) 模型，使用三明治估计（而非模型基于的统计量）确定治疗与周交互作用检验的 <span class="math inline">\(F\)</span> 值和 <span class="math inline">\(p\)</span> 值。使用经典三明治估计和 Morel 校正的三明治估计分别完成。</li>
<li>将 ANTE(1) 模型的三明治估计结果与 AR(1) 和 UN 模型的三明治估计结果进行比较。相比于协方差模型对基于模型的估计的影响，协方差模型的选择如何影响 <span class="math inline">\(F\)</span> 值和 <span class="math inline">\(p\)</span> 值？</li>
</ol>
</li>
<li>本习题涉及文件 Ch7_prob5.sas. 数据来自二项响应的重复测量研究。有 12 个个体（通过变量 ID 编码）；随机分配 6 名个体至每种治疗。每种治疗在五个等距时间点上观测（两种治疗观测的时间点相同）。SAS 程序显示了研究人员试图确定用于分析的正确协方差结构。总的来说，研究人员首先假定了一个非结构化的工作协方差模型。然后进行了检验：以查看是否存在序列相关的证据 (“DiagR” test)，若存在，那么是否有理由使用非结构化协方差模型，或者是否可以将结构简化为 Toeplitz 模型。研究人员还从 ANTE(1) 开始，经过一系列简化，最终简化为 AR(1). 假定简化后得到了 Toeplitz 和 AR(1) 工作协方差模型，那么根据提供的拟合统计量（例如，<span class="math inline">\(-2\log \text{PL}\)</span> 和广义 <span class="math inline">\(\chi^2/DF\)</span>）比较这两个模型，以选择最佳模型。
<ol style="list-style-type: lower-alpha">
<li>现在暂时不要怀疑研究人员的方法。根据研究者的逻辑，你会“排除”和“纳入”以下哪种协方差结构？
<img src="figure/figure%20231.png#center" style="width:80.0%">
</li>
<li>研究人员的方法存在致命缺陷，导致 a. 的回答无法使用。该缺陷是什么（你应该能用最多一两句话来回答）。</li>
<li>作为一般原则，应该如何重新指定模型才能真正合理地进行协方差模型选择？</li>
<li>根据你在 c. 中的回答重新进行模型选择。在此基础上你选择什么协方差模型？记录你的步骤和支持证据。</li>
<li>无需详细说明，根据你在 d. 中选择的模型，你对是否存在治疗和/或时间效应会得出什么结论。证据是什么？</li>
</ol>
</li>
</ol>
</div>
</div>
<div class="chapter-nav">
<div class="prev"><a href="chap6.html"><span class="header-section-number">6</span> 推断（一）</a></div>
<div class="next"><a href="chap8.html"><span class="header-section-number">8</span> 处理和解释变量结构</a></div>
</div></main><div class="col-md-3 col-lg-2 d-none d-md-block sidebar sidebar-chapter">
<nav id="toc" data-toggle="toc" aria-label="On this page"><h2>On this page</h2>
<ul class="nav navbar-nav">
<li><a class="nav-link" href="#chap7"><span class="header-section-number">7</span> 推断（二）</a></li>
<li><a class="nav-link" href="#sec7-1"><span class="header-section-number">7.1</span> 介绍</a></li>
<li>
<a class="nav-link" href="#sec7-2"><span class="header-section-number">7.2</span> 协方差分量的正式检验</a><ul class="nav navbar-nav">
<li><a class="nav-link" href="#sec7-2-1"><span class="header-section-number">7.2.1</span> 仅方差分量 LMMs 的基于 ANOVA 的检验</a></li>
<li><a class="nav-link" href="#sec7-2-2"><span class="header-section-number">7.2.2</span> 用于协方差分量检验的 Wald 统计量以及为什么不应使用</a></li>
<li><a class="nav-link" href="#sec7-2-3"><span class="header-section-number">7.2.3</span> 协方差分量的似然比检验</a></li>
<li><a class="nav-link" href="#sec7-2-4"><span class="header-section-number">7.2.4</span> GLMMs 的 PL 与积分近似的结果</a></li>
</ul>
</li>
<li>
<a class="nav-link" href="#sec7-3"><span class="header-section-number">7.3</span> 用拟和统计量比较协方差模型</a><ul class="nav navbar-nav">
<li><a class="nav-link" href="#sec7-3-1"><span class="header-section-number">7.3.1</span> AIC 和 AICC</a></li>
<li><a class="nav-link" href="#sec7-3-2"><span class="header-section-number">7.3.2</span> BIC</a></li>
<li><a class="nav-link" href="#sec7-3-3"><span class="header-section-number">7.3.3</span> 协方差模型比较的应用</a></li>
</ul>
</li>
<li>
<a class="nav-link" href="#sec7-4"><span class="header-section-number">7.4</span> 区间估计</a><ul class="nav navbar-nav">
<li><a class="nav-link" href="#%E9%A6%96%E5%85%88%E4%B8%8D%E8%AF%A5%E5%81%9A%E4%BB%80%E4%B9%88">首先：不该做什么</a></li>
<li><a class="nav-link" href="#sec7-4-1"><span class="header-section-number">7.4.1</span> 基于 \(\chi^2\) 的 Wald 法</a></li>
<li><a class="nav-link" href="#sec7-4-2"><span class="header-section-number">7.4.2</span> 基于似然的方法</a></li>
</ul>
</li>
<li><a class="nav-link" href="#sec7-5"><span class="header-section-number">7.5</span> 总结</a></li>
<li><a class="nav-link" href="#exe7">练习</a></li>
</ul>
<div class="book-extra">
<ul class="list-unstyled">
</ul>
</div>
</nav>
</div>
</div>
</div> <!-- .container -->
<footer class="bg-primary text-light mt-5"><div class="container"><div class="row">
<div class="col-12 col-md-6 mt-3">
<p>"<strong>广义线性混合模型</strong>: 现代概念、方法和应用" was written by Wang Zhen. It was last built on 2024-05-19.</p>
</div>
<div class="col-12 col-md-6 mt-3">
<p>This book was built by the <a class="text-light" href="https://bookdown.org">bookdown</a> R package.</p>
</div>
</div></div>
</footer><!-- dynamically load mathjax for compatibility with self-contained --><script>
(function () {
var script = document.createElement("script");
script.type = "text/javascript";
var src = "true";
if (src === "" || src === "true") src = "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.9/latest.js?config=TeX-MML-AM_CHTML";
if (location.protocol !== "file:")
if (/^https?:/.test(src))
src = src.replace(/^https?:/, '');
script.src = src;
document.getElementsByTagName("head")[0].appendChild(script);
})();
</script><script type="text/x-mathjax-config">const popovers = document.querySelectorAll('a.footnote-ref[data-toggle="popover"]');
for (let popover of popovers) {
const div = document.createElement('div');
div.setAttribute('style', 'position: absolute; top: 0, left:0; width:0, height:0, overflow: hidden; visibility: hidden;');
div.innerHTML = popover.getAttribute('data-content');
var has_math = div.querySelector("span.math");
if (has_math) {
document.body.appendChild(div);
MathJax.Hub.Queue(["Typeset", MathJax.Hub, div]);
MathJax.Hub.Queue(function() {
popover.setAttribute('data-content', div.innerHTML);
document.body.removeChild(div);
})
}
}
</script>
</body>
</html>