-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathchap9.html
637 lines (617 loc) · 87.6 KB
/
chap9.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
624
625
626
627
628
629
630
631
632
633
634
635
636
637
<!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>第 9 章 多水平模型 | 广义线性混合模型</title>
<meta name="author" content="Wang Zhen">
<meta name="description" content="9.1 设计结构类型：单水平和多水平模型定义 在第 8 章中，我们研究了构建和使用处理（或更一般地说，解释变量结构）的主要策略。在本章中，我们将注意力转向设计，或者 Fisher (1935) 所说的线性预测器的“地形”方面。简而言之，处理方面主要决定了我们如何指定 \(\symbf{X\beta}\)，而设计或地形方面则倾向于推动与使用 \(\symbf{Zb}\) 相关的指定和策略。...">
<meta name="generator" content="bookdown 0.38 with bs4_book()">
<meta property="og:title" content="第 9 章 多水平模型 | 广义线性混合模型">
<meta property="og:type" content="book">
<meta property="og:description" content="9.1 设计结构类型：单水平和多水平模型定义 在第 8 章中，我们研究了构建和使用处理（或更一般地说，解释变量结构）的主要策略。在本章中，我们将注意力转向设计，或者 Fisher (1935) 所说的线性预测器的“地形”方面。简而言之，处理方面主要决定了我们如何指定 \(\symbf{X\beta}\)，而设计或地形方面则倾向于推动与使用 \(\symbf{Zb}\) 相关的指定和策略。...">
<meta name="twitter:card" content="summary">
<meta name="twitter:title" content="第 9 章 多水平模型 | 广义线性混合模型">
<meta name="twitter:description" content="9.1 设计结构类型：单水平和多水平模型定义 在第 8 章中，我们研究了构建和使用处理（或更一般地说，解释变量结构）的主要策略。在本章中，我们将注意力转向设计，或者 Fisher (1935) 所说的线性预测器的“地形”方面。简而言之，处理方面主要决定了我们如何指定 \(\symbf{X\beta}\)，而设计或地形方面则倾向于推动与使用 \(\symbf{Zb}\) 相关的指定和策略。...">
<!-- 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="" 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="active" 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="chap9" class="section level1" number="9">
<h1>
<span class="header-section-number">第 9 章</span> 多水平模型<a class="anchor" aria-label="anchor" href="#chap9"><i class="fas fa-link"></i></a>
</h1>
<div id="sec9-1" class="section level2" number="9.1">
<h2>
<span class="header-section-number">9.1</span> 设计结构类型：单水平和多水平模型定义<a class="anchor" aria-label="anchor" href="#sec9-1"><i class="fas fa-link"></i></a>
</h2>
<p>在第 <a href="chap8.html#chap8">8</a> 章中，我们研究了构建和使用处理（或更一般地说，解释变量结构）的主要策略。在本章中，我们将注意力转向设计，或者 Fisher (1935) 所说的线性预测器的“地形”方面。简而言之，处理方面主要决定了我们如何指定 <span class="math inline">\(\symbf{X\beta}\)</span>，而设计或地形方面则倾向于推动与使用 <span class="math inline">\(\symbf{Zb}\)</span> 相关的指定和策略。</p>
<p>正如我们在第 <a href="chap2.html#chap2">2</a> 章中所看到的，广义上定义的设计或地形方面涉及<strong>重复单元</strong> (units of replication) 和<strong>观察的单元</strong> (units of observation)<a class="footnote-ref" tabindex="0" data-toggle="popover" data-content='<p>译者注：重复单元与观察的单元的定义见 <a href="chap2.html#sec2-2-3">2.2.3</a> 节。请注意勿将后者与“观察单元”混淆。若译名确实容易混淆，请谨记英文名及其含义。</p>'><sup>24</sup></a> 及其结构。这些单元是什么？我们有单元是更大单元的子集或细分吗？我们有没有有意义的单元分组——层、集群或区组 (strata, clusters or blocks)？我们如何对这些单元进行随机化？</p>
<p>本章我们主要关注的是正确识别重复单元和随后的 <span class="math inline">\(\symbf{Zb}\)</span> 的指定。我们的重复单元是否只有一种尺寸 (size)，还是有不同尺寸的重复单元？如果是后者，那么我们就有了一个<strong>多水平研究</strong> (multi-level study). 多水平研究的名称因所在领域而异。从历史上看，农业统计学家是最先认识和命名多水平研究的——他们称之为<strong>裂区实验</strong> (split-plot experiments). 在医疗中心，生物统计学家称它们为<strong>嵌套析因/集群设计</strong> (nested factorial or clustered designs). 社会科学家则称之为<strong>分层设计</strong> (hierarchical designs). 无论不同领域的术语如何，它们指的都是同一设计结构。</p>
<p>模型中表征设计的组成部分有时被称为干扰参数 (nuisance parameters). 这既不适当也不利于理解，因为“讨厌”一词常被解读为“不重要”。设计部分可能很复杂，我们如何定义线性预测器来表征它们对后续推断有重大影响。在建模练习中，当事情出现严重错误时，往往是因为在设计方面出现了问题。当存在多水平结构而未能考虑或考虑不当时，是导致事情变糟的最常见原因。</p>
<p>关于设计或地形结构的许多术语都来源于实验设计。但这并不意味着本章的讨论仅适用于设计实验。调查、观察性研究和准实验都是有结构的。它们都有观察的单元。它们都有与重复单元至少是不严格对应的实体。聚类、分层和病例对照匹配都是区组的形式。为了构建一个有用且信息丰富的模型，所有这些都需要被适当地理解和描述。在某些方面，本章对于非正式设计的实验或随机试验的研究可能更重要，原因在于参与非实验性研究的人有时可能忽视设计原则的应用。因此，为避免任何误解，每当我们进行比较性的推断时，设计原则都适用。本章使用设计术语是因为它是多水平建模的通用语言，但读者需要理解本章中的语言是广义的。</p>
<p>综上所述，让我们为多水平研究是什么以及不是什么建立一个精确的工作定义。</p>
<p>我们将一项研究定义为多水平的，当且仅当它具有至少两个随机变异源，以及构成研究架构的至少两种不同尺寸的单元。</p>
<p>或者，我们也可以根据不同的随机化方案数量来定义多水平研究。关键在于：不止一种尺寸的单元，或等价地，多个随机化方案。如果我们只有一种尺寸的单元，意味着重复单元和观察的单元是相同的，或等价地，所有处理水平或因素水平组合都使用单一的随机化方案，那么我们就不是在进行多水平研究。一种尺寸的单元：单水平。两种或更多种尺寸的单元：多水平。单一的随机化方案：单水平。两步或更多步的随机化：多水平研究。</p>
<p>最简单的多水平研究是具有多个抽样单元的完全随机设计。在这类研究中，重复单元由几个个体组成，如围栏中的动物、诊所中的患者、区 (plot) 中的植物。本章重点关注更复杂的多水平研究，因为这些研究通常会给建模师和数据分析师带来最大的困难。</p>
<p>在 <a href="chap9.html#sec9-2">9.2</a> 节中，我们将更详细地了解多水平研究是如何产生的。我们还回顾了第 <a href="chap2.html#chap2">2</a> 章中介绍的框架 ANOVA（“Fisher会做什么？”或 WWFD ?）过程，这是一种帮助我们可视化单水平和多水平研究之间区别的方法，也是一种帮助人们正确指定多水平设计组成部分的方法。</p>
<p>如果不考虑分区 (blocking)，我们就无法讨论多水平研究。在多水平研究中使用分区时，会出现单水平设计中不会出现的问题。对于单级设计，决定是将区组 (blocks) 效应视为固定效应还是随机效应通常对模型估计或推断的影响最多只是轻微的。但多水平研究则不然。具体来说，将区组效应定义为固定效应可能会让你陷入无法轻易逃脱的困境。主要的两个问题是：错误的推断统计量和虚假的不可估性。<a href="chap9.html#sec9-3">9.3</a> 节将重点讨论分区。</p>
<p>在 <a href="chap9.html#sec9-4">9.4</a> 节中，我们考虑一些示例，说明多因素处理结构和多水平设计结构如何组合在一起形成整个线性预测器。根据定义，多因素模型必须是混合模型。本章到 <a href="chap9.html#sec9-4">9.4</a> 节的内容同样适用于非高斯 GLMMs 和 LMMs.</p>
<p>在前面的章节中，特别是第 <a href="chap3.html#chap3">3</a>、<a href="chap5.html#chap5">5</a> 和 <a href="chap6.html#chap6">6</a> 章，我们讨论了边际模型和条件模型、G 侧和 R 侧协方差建模以及 GEE 型模型。在 <a href="chap9.html#sec9-5">9.5</a> 节中，我们将研究这些问题如何具体应用于多水平模型。</p>
</div>
<div id="sec9-2" class="section level2" number="9.2">
<h2>
<span class="header-section-number">9.2</span> 多水平模型的类型及其产生方式<a class="anchor" aria-label="anchor" href="#sec9-2"><i class="fas fa-link"></i></a>
</h2>
<p>重复单元和随机化方案提供了我们构建 <span class="math inline">\(\symbf{Zb}\)</span> 所需的关键信息。在本节中，我们回顾第 <a href="chap2.html#chap2">2</a> 章中的一些关键概念以及它们如何在多水平场景中体现。</p>
<div id="sec9-2-1" class="section level3" number="9.2.1">
<h3>
<span class="header-section-number">9.2.1</span> 重复单元——不仅在设计的实验中<a class="anchor" aria-label="anchor" href="#sec9-2-1"><i class="fas fa-link"></i></a>
</h3>
<p>对于设计的实验，重复单元（设计语言中的实验单元）从规划阶段到数据分析阶段都占据着重要地位。它正式定义为独立分配处理的最小实体。“独立”这个词至关重要。如果我们将动物圈养，并通过一个公共的围栏喂食器分配并应用处理，则实验单元就是围栏。动物版对着分配；我们可以单独测量它们，但它们到处理水平的分配不是独立的。我们不能在围栏内将一只动物分配给一个处理，将另一只动物分配给另一个不同的处理。实验单元思维方式深深地嵌入在实验设计的思考过程中。这大大有助于从实验设计中构造线性预测器。我们已经有了所有必要的信息。</p>
<p>对于调查、观察性研究和准实验，我们必须更加努力，但我们仍然需要识别这些要素，以便构建能够产生准确分析的线性预测器。每项调查、观察性研究等都具有类似于重复单元的概念。让我们看一个随机试验的例子，它几乎是但又不完全是一个设计的实验，来加深我们对术语的理解，然后考虑这些思想如何转化为调查和观察性研究。</p>
<p>为此，我们使用 Fisher 的术语“地形”，就像我们在第 <a href="chap2.html#chap2">2</a> 章中所做的那样。“地形”一词很有帮助，因为它不受可能妨碍讨论的设计或勘测术语的影响，但我们需要将“地形”理解为广义上的研究结构，而不是狭义上的地理或地质概念。</p>
<p>最简单的地形结构有一种单元。我们对这些单元进行观察，并将处理水平分配给这些单元。这种结构定义了一种完全随机设计 (completely randomized design, CRD)。CRD 显然是一个单水平结构。就调查而言，简单的随机样本与 CRD 具有相同的结构，因此需要考虑类似的建模因素。</p>
<p>CRD 是唯一仅具有一种单元尺寸的设计，但 CRD 并不是唯一的单水平设计。如果我们有两种尺寸的单元，比如教室和教室内的学生，并且教室接受处理分配（例如课程 A1 分配到一个教室，课程 A2 分配到另一个教室等），但我们对单个学生进行测试以查看处理的效果，那么我们就有了一个带有子抽样 (sub-sampling) 的 CRD，这是最简单的多水平结构。假设教室被分区——学校是一个明显的分区标准。现在我们的研究中有三种尺寸的单元。如果我们仍然在教室水平上接受处理（比如在任何给定的学校，一间教室接受课程 A1，另一间教室接受课程 A2），那么我们就会采用随机区组设计——如果我们仍然对单个学生进行测量，那么我们就有了子抽样。带子抽样的 CRD 和随机区组设计（带或不带子抽样）都是单因素设计，但它们仍然需要多水平模型。为什么？因为，尽管它们具有单一的处理因素
，并且需要一个将因素水平分配给重复单元的随机化方案，但它们仍然需要多个随机变量——除了观测的（基于随机效应的条件）分布之外，线性预测器中至少有一个随机效应——以充分考虑相关概率过程的变异源。</p>
<p>当我们添加第二个因素时，会出现更复杂的多水平研究。假设除了课程（A1 或 A2）作为一个因素外，我们还有一个因素 B。如果像之前一样将因素的组合分配给教室，如果我们在教室水平进行观察，那么这仍然是一个单水平实验；如果我们在单个学生水平进行观察，那么这是一个简单的多水平子抽样模型。我们有一个由因素组合定义的 a × b 个处理，但重复单元仍然是教室。然而，假设我们像以前一样将课程分配给教室，但将因素 B 的水平分配给整个学校或教室内的单个学生。现在，我们有两个不同尺寸的重复单元，并且我们需要两个随机化方案，一个用于将不同的课程（因素 A 的水平）分配给教室，另一个用于分配因素 B 的水平。这些是更复杂的多水平研究的定义标准。</p>
<p>这在观察性研究中是如何运作的？假设我们想知道暴露在环境危害中是否会增加某种健康问题的风险。一个可能的设计会涉及识别暴露于这种危害的人群和一组未暴露但在其他方面相似的人群。我们可能会进一步尝试创建匹配的配对——一个暴露和一个未暴露但具有相似特征的人。匹配在观察性研究中产生线性预测效应的方式与在设计性实验中的分区是一样的。个人显然是观察的单元。那么重复的单元是什么呢？这里的“处理”因素是危害，存在或不存在。危害是“分配”给地点的。如果我们想要将这种推断扩展到所有这类危害，重复的单元就是地点，而不是个体。如果我们只查看一个地点，我们就会相应地限制推断范围。</p>
<p>在一个调查样本中，假设我们想要找出学生的政治取向是否因高等教育机构的类型而异。例如，公立大学是否拥有更高比例的自我认同的保守派？私立文理学院是否拥有更高比例的自我认同的自由派？在这种情况下，因素是机构的类型，分配给类型的最小实体是整个学院或大学。被调查的单个学生构成子抽样。如果我们加入第二个因素，比如专业，那么特定机构中特定专业的整个项目就是关于这个因素的重复单元。现在我们有一个多水平研究：两个因素，机构类型和专业；两个重复单元，学院或大学和学院内部的项目。</p>
<p>我们在第 <a href="chap2.html#chap2">2</a> 章中介绍的框架 ANOVA (WWFD?) 过程可以是识别多水平研究并将其转化为线性预测器的特别有效的方法。</p>
</div>
<div id="sec9-2-2" class="section level3" number="9.2.2">
<h3>
<span class="header-section-number">9.2.2</span> 重温“Fisher 会怎么做？”——地形和处理部分<a class="anchor" aria-label="anchor" href="#sec9-2-2"><i class="fas fa-link"></i></a>
</h3>
<p>让我们将上述示例转换为框架 ANOVA 术语。这将有助于我们直观地了解单水平和多水平。首先，最简单的完全随机形式成为：</p>
<div class="inline-figure"><img src="figure/figure%20274.png#center" style="width:100.0%"></div>
<p>在这里，单水平结构是显而易见的—— “trt” 必须是因子，而 “unit(trt)” 是重复单元的唯一可能候选者。回想一下我们在第 <a href="chap3.html#chap3">3</a> 章中对条件和边际模型的讨论，必须将 “unit(trt)” 所在行视为“残差”，因此在高斯模型的线性预测器中不会出现，但它可能会并且通常应该作为非高斯模型线性预测器中的随机效应出现。这种残差与线性预测器效应是过度离散的一个方面，我们将在第 <a href="#chap11"><strong>??</strong></a> 章和第 <a href="#chap12"><strong>??</strong></a> 章中更深入地探讨它。</p>
<p>区组设计有一个额外特征：</p>
<div class="inline-figure"><img src="figure/figure%20275-1.png#center" style="width:100.0%"></div>
<p>现在在地形方面有两个元素。在这个有限的意义上，我们有一个多水平模型，但我们仍然只有一个因素和唯一一个可能的重复单元：区组-处理的组合。因此，从 GLMM 的角度来看，我们有一个单水平模型。同样，区组-处理效应必须视为残差，因此在高斯模型的线性预测器中不会出现，但在某些非高斯模型中可能会出现。</p>
<p>如果我们还包括子抽样，例如我们的学校/教室示例中的学生，我们就有</p>
<div class="inline-figure"><img src="figure/figure%20275-2.png#center" style="width:100.0%"></div>
<p>“unit(block)” 行下方的实线将框架 ANOVA 分为上下两部分。上部分是处理水平随机化发生的地方；而下部分是无法独立于区组-处理组合进行随机化的抽样单元。但是，它们除了区组和 unit(block) 的变异外，还有自己的变异。地形方面有三个水平：区组、单元和子单元。模型必须以某种方式考虑这三个水平。</p>
<p>对<strong>分层线性模型</strong> (hierarchical linear model, <strong>HLM</strong>) 用户特别重要的一点评论是：在 HLM 术语中，subunit(block × trt) 被称为水平 1 (level 1)，即研究中的最小单元。在之前的框架 ANOVA 中，对于没有子抽样的区组设计，区组-处理组合是代表水平 1 的项。这就是 HLM 术语可能有些不够清楚甚至令人困惑的地方。HLM 水平并不区分观察的单元和重复单元。在有子抽样的情况下，重复单元是水平 2，没有子抽样的情况下，重复单元是水平 1。这不一定是个问题，但可能会令人困惑。HLM 用户确实需要在这方面多加努力，因为 HLM 术语并不强调设计概念。</p>
<p>对于环境危害观察性研究，可能的框架 ANOVA 为：</p>
<div class="inline-figure"><img src="figure/figure%20276-1.png#center" style="width:100.0%"></div>
<p>关于正在研究的处理因素，我们仍然有一个单水平设计，但现在我们有了三层嵌套结构。这种结构极易导致<strong>伪重复</strong> (pseudo-replication)：在这种情况下，错误地将配对或个体作为重复。在这个设定中，它们显然不是。请注意，在 HLM 术语中，真正的重复，site，是模型水平 3 的组件。再次强调，HLM 水平和重复单元之间没有明确的对应关系。</p>
<p>最后，大学类型-专业的例子也是一个多水平模型。</p>
<div class="inline-figure"><img src="figure/figure%20276-2.png#center" style="width:100.0%"></div>
<p>这里我们有两个不同的重复单元：类型/机构，以及专业/项目。请注意，我们将框架 ANOVA 分为三个层次：两个水平与处理设计相关，即类型/机构和专业/项目，以及一个水平用于子抽样，即单个学生。通常，对于高斯模型，我们必须将学生效应视为残差，但对于一些非高斯模型，学生效应应在线性预测器中表示。此外，我们在处理结构方面有一些灵活性。如果不同类型的机构有不同的专业，我们可以用一个嵌套效应—— major(institution) ——来代替单独的主效应（major）和交互效应（major × institution）。</p>
<details><summary><font color="#8B2232">图 2.6</font>
</summary><img src="figure/figure%202.6.png#center" style="width:60.0%"></details><p><br>
图 2.6 中的不连通不完全设计—— 10 个区组，每个尺寸为 3，有 6 种处理，其中一组（处理 0, 1, 2）始终一起位于 5 个区组中，另一组（处理 3, 4 和 5）始终位于另五个区组中——为我们提供了多水平结构的另一种视角。框架 ANOVA</p>
<div class="inline-figure"><img src="figure/figure%20277.png#center" style="width:70.0%"></div>
<p>揭示了多水平结构的面貌。重复单元是组 (set) 的区组以及组内处理的区组内单元 (units within block for treatments within a set)。如果处理设计具有额外的结构——例如，析因设计中 “set” 对应于因素 A，每个组内的处理对应于因素 B 的水平——那么我们可以在区组水平用 A 替换 “set”，在单元水平用 B 和 A × B 替换 “trt(set)”。</p>
<p>我们可以根据特定研究的要求和限制，以无数种方式改变这些结构。Federer and King (2007) 一书中列举了裂区设计的各种变体。任何多水平研究——无论是设计的研究、观察性研究还是调查——的关键特征都是至少存在两个不同尺寸的重复单元。这意味着所有多水平模型都必须是混合模型。与最小单元相关的变异可以嵌入 <span class="math inline">\(\symbf y\mid \symbf b\)</span> 的分布中。实际上，对于具有高斯分布的观测，这是必须的。另一方面，与较大单元相关的变异必须在线性预测器中建模——也就是说，在 <span class="math inline">\(\symbf b\)</span> 中必须至少包含较大的重复单元作为随机模型效应。换句话说，如果我们有一个多水平结构，我们的线性预测器就必须包含一个 <span class="math inline">\(\symbf{Zb}\)</span> 分量。</p>
</div>
</div>
<div id="sec9-3" class="section level2" number="9.3">
<h2>
<span class="header-section-number">9.3</span> 区组在多水平模型中的作用<a class="anchor" aria-label="anchor" href="#sec9-3"><i class="fas fa-link"></i></a>
</h2>
<p>区组效应是固定的还是随机的？自从 Yates (1940) 开发了区组间信息的恢复（原始的随机区组效应方法）以来，这是持续争议的问题。这一问题对于多水平模型尤其重要，因为固定或随机的决定会影响可估性。</p>
<p>在建模中，我们有两种不可估性——真实的和虚假的。真实不可估性的例子包括仅在截距上定义的函数或 ANOVA 型模型中孤立的处理效应：<span class="math inline">\(\eta + \tau_i\)</span> 是可估的，但 <span class="math inline">\(\eta\)</span> 和 <span class="math inline">\(\tau_i\)</span> 单独是不可估的。在仅含主效应的设计中的交互效应，以及因子设计中缺少因子组合的某些交互对比，都是真实不可估性的其他常见例子。虚假不可估性的最严重的例子出现在多水平设计中。其中大多数可以直接归因于将区组当作固定效应来处理，这要么是有意为之，要么是更常见地，使用了不能识别固定效应与随机效应差异的软件，从而不适当地评估了可估性。</p>
<p>在 <a href="chap9.html#sec9-3-1">9.3.1</a> 节中，我们将更详细地讨论固定区组或随机区组的争议的各个方面。在 <a href="chap9.html#sec9-3-2">9.3.2</a> 节中，我们具体考虑了多水平模型对可估性的影响。</p>
<div id="sec9-3-1" class="section level3" number="9.3.1">
<h3>
<span class="header-section-number">9.3.1</span> 重温“固定区组效应与随机区组效应”<a class="anchor" aria-label="anchor" href="#sec9-3-1"><i class="fas fa-link"></i></a>
</h3>
<p>区组效应是固定的还是随机的？用线性预测器术语来表述这个问题，我们是把模型写成 <span class="math inline">\(\eta+\tau_i+\rho_j\)</span>，其中 <span class="math inline">\(\rho_j\)</span> 表示区组效应，是 <span class="math inline">\(\symbf{X\beta}\)</span> 的一部分；还是写成 <span class="math inline">\(\eta+\tau_i+r_j\)</span>，其中 <span class="math inline">\(r_j\)</span> 表示区组效应，有一个假定的分布，是 <span class="math inline">\(\symbf{Zb}\)</span> 的一部分？争议主要有两个方面，一个是定义性的，一个是实用性的，后者涉及两个模型的估计和推断的小样本行为。</p>
<p>让我们从回顾定义开始。根据定义，固定模型效应表示具有有限水平的因素，并且我们的设计提供了对该因素的每个水平的观察。随机模型效应来源于一个因素，其水平构成了一个总体。我们无法观察整个总体；相反，我们通过某种形式的抽样来观察总体的代表。</p>
<p>对于随机效应，我们通过模型表示；对于固定效应，我们则是详尽地观察。因此，对于固定效应的推断不会超出实际观察到的水平（使用回归模型在定量水平之间进行插值的情况除外），而对于随机效应的推断则扩展到所观察到的水平所代表的总体中。从定义的角度来看，这场争议的关键在于：我们是将观察到的区组视为整个总体，还是将它们理解为来自一个目标总体的代表？我们是否满足于将推断仅限于我们实际观察到的区组，即狭义或特定区组的推断，还是我们的目标要求进行广义或总体平均的推断？支持随机区组效应观点的人认为，根据纯粹的定义标准，假设固定区组效应是一个缺乏内在一致性的建模决策。</p>
<p>固定区组的拥护者通常会回避定义问题并关注实际问题。而实际问题是什么呢？</p>
<p>首先，如果我们将讨论限制在均衡数据（即无缺失数据的随机完全区组设计）以及根据处理差异定义的可估函数上，那么争议就没有实际意义。固定和随机区组效应的分析会产生相同的结果。</p>
<p>那么不均衡数据？在实际研究中，数据缺失常常发生。此外，在许多情况下，我们使用不完全区组设计而不是完全区组设计来提高效率 (efficiency). 我们将在第 <a href="#chap21"><strong>??</strong></a> 章中更详细地讨论后一个问题，该章涉及规划设计的精度、功效 (power) 和样本量。</p>
<p>这就提出了第二个主要的实际问题。对于不均衡区组设计，Yates (1940) 表明，恢复区组间信息并将其与区组内信息相结合可以提高效率。后来，区组间和区组内信息的综合分析被证明等价于随机区组效应混合模型分析（例如，参见 Bhattacharya, 1998）。随机区组分析更有效——前提是我们有已知的区组方差。但当我们对其未知时又该如何呢？这是一个重要的问题，因为在实际研究中，我们从来都无法确切知道。这就是真正的争议所在，这可以表述为：为了从恢复区组间信息中获益，区组方差的估计需要有多好？得到一个好的估计需要满足什么条件？</p>
<p>固定区组的主要论点基于一个假设，即在有较少区组的研究中（例如，在农业中，使用 3 - 5 个区组是常规做法），区组自由度太少，因此无法得到区块方差的良好估计，因此恢复区组间信息是适得其反的。如果该假设是真的，那么支持固定区组的论点就站得住脚。如果不是，则不然。</p>
<p>那么情况有多糟糕呢？让我们考虑一个案例，固定区组拥护者可能会称之为“随机区组模型的噩梦场景”。以下是分区安排：</p>
<div class="inline-figure"><img src="figure/figure%20279-1.png#center" style="width:30.0%"></div>
<p>我们有 3 个处理和 3 个区组，每个区组的尺寸为 2。该设计仅允许使用 2 个自由度来估计区组方差。固定区组和随机区组分析比较起来如何？为了回答这个问题，我们转向模拟。考虑两种场景：1) 所有三种处理的均值相等（例如，<span class="math inline">\(\mu_1 = \mu_2 = \mu_3 = 0\)</span>）；2) 处理之间相差 15 个单位（例如，<span class="math inline">\(\mu_1 = 0,\mu_2 = 15,\mu_3 = 30\)</span>）。场景 1 允许我们评估 I 类错误控制；场景 2 允许我们评估功效。假设数据符合高斯分布。不失一般性，设 <span class="math inline">\(\sigma^2=1\)</span>。在理论上，对于 <span class="math inline">\(\alpha = 0.05\)</span>，在场景 1 下，对于 <span class="math inline">\(H_0\colon \mu_1 = \mu_2 = \mu_3\)</span> 的拒绝率应为 5%. 在场景 2 下，随机区组分析的拒绝率取决于区组方差 <span class="math inline">\(\sigma^2_B\)</span>。对于本模拟中使用的 <span class="math inline">\(\sigma^2_B=0.25\)</span>，在理论上，固定区组和随机区组分析下的拒绝率应分别为 0.806 和 0.849. 读者可以很容易地验证，改变 <span class="math inline">\(\sigma^2_B\)</span> 不会改变基本发现。这些发现是什么？以下表格提供了 1000 次模拟实验的摘要。</p>
<ul>
<li>场景 1：所有均值相等</li>
</ul>
<div class="inline-figure"><img src="figure/figure%20279-2.png#center" style="width:80.0%"></div>
<ul>
<li>场景 2：所有均值不等</li>
</ul>
<div class="inline-figure"><img src="figure/figure%20279-3.png#center" style="width:80.0%"></div>
<p>鉴于该设计中每个处理有两个重复，且只有 1 个自由度用于残差，我们不应该对任何一种分析抱有太大期望。当然，我们不会在实际研究中采用这种设计。也就是说，如果区组太少确实会影响随机区组效应分析，我们应该在这种设计中看到这种影响被最大化。然而事实并非如此。可以肯定的是，<span class="math inline">\(\sigma^2_B\)</span> 和 <span class="math inline">\(\sigma^2\)</span> 本身的估计并不准确。然而，关于 I 类错误率和功效特性的比较显示，随机区组分析比固定区组分析更有利。</p>
<p>这是因为处理的检验统计量取决于 <span class="math inline">\(Var{\left[\symbf{K'}\left(\symbf{X'V}^{-1}\symbf{X}\right)^-\symbf{K}\right]}\)</span> 的总体估计，其中 <span class="math inline">\(\symbf K\)</span> 通过可估函数 <span class="math inline">\(\symbf{K'\beta}\)</span> 定义了 <span class="math inline">\(H_0\colon\mu_1=\mu_2=\mu_3\)</span>。重要的是这个估计的质量，而不是单独的 <span class="math inline">\(\sigma^2_B\)</span> 和 <span class="math inline">\(\sigma^2\)</span>。底线：这远不是固定区组分析的决定性胜利——如果说有什么的话，随机区组分析似乎在整体上略微优于固定区组模型——尤其是在检验处理效应方面。</p>
<p>到目前为止，我们一直在考虑单水平的情况。如果从得分来看，在定义方面，随机区组模型因对方弃权而获胜；在实际应用方面，对于区组较少的设计，两者基本上打了个平手，而对于区组较多的设计，随机区组模型略胜一筹。如果我们只关注单水平设计，这场争议就没有看上去那么严重。现在，我们将注意力转向多水平模型，这才是真正的问题所在。</p>
</div>
<div id="sec9-3-2" class="section level3" number="9.3.2">
<h3>
<span class="header-section-number">9.3.2</span> 固定区组、多水平设计与虚假的不可估性<a class="anchor" aria-label="anchor" href="#sec9-3-2"><i class="fas fa-link"></i></a>
</h3>
<p>让我们从图 2.6 中的不连通设计开始。回想，平面图为</p>
<div class="inline-figure"><img src="figure/figure%20280.png#center" style="width:30.0%"></div>
<p>我们在第 <a href="chap2.html#chap2">2</a> 章中看到，对于该设计，假定区组效应 <span class="math inline">\(\rho_j\)</span> 是固定的，则朴素线性预测器 <span class="math inline">\(\eta+\tau_i+\rho_j\)</span> 会导致边际处理均值的可估性问题，以及涉及不在共同区组中处理差异的可估性问题。例如，考虑定义为 <span class="math inline">\(\eta+\tau_i+1/10\sum_{j=1}^{10}\rho_j\)</span> 的边际均值，为了使边际均值可估，必须存在系数 <span class="math inline">\(a_{ij}\)</span>，使得 <span class="math inline">\(\sum_{i,j}a_{ij}\left(\eta+\tau_i+\rho_j\right)=\eta+\tau_i+1/10\sum_j\rho_j\)</span>。这意味着 <span class="math inline">\(\sum_{i,j}a_{ij}=1\)</span>（为获取 <span class="math inline">\(\eta\)</span>），<span class="math inline">\(\sum_ja_{ij}=1\)</span> 以及对于所有 <span class="math inline">\(i'\ne i\)</span> 有 <span class="math inline">\(\sum_ja_{i'j}=0\)</span>（为获取 <span class="math inline">\(\tau_i\)</span>），并且对于 <span class="math inline">\(j=1,2,\ldots\)</span> 有 <span class="math inline">\(\sum_ia_{ij}=1/10\)</span>。前两个条件很容易同时满足，但前两个条件排除了了第三个条件的满足。或者，我们可以证明 <span class="math inline">\(\symbf K\)</span> 不满足 <span class="math inline">\(\symbf{K'}(\symbf{X'X})^{-}\symbf{X'X}=\symbf{K'}\)</span> 准则，因此 <span class="math inline">\(\mu+\tau_i+1/10\sum_j\rho_j\)</span> 不可估。</p>
<p>从区块效应的定义角度来看，这个例子为我们首次展示了虚假的不可估性。如果我们把区块效应视为随机的，那么恢复区组间信息将使我们能够准确分析这一设计的真正内容：一个嵌套的析因设计，其中析因结构由 set 和 treatment(set) 定义。这就是我们 <a href="chap9.html#sec9-2-2">9.2.2</a> 节末尾看到的框架 ANOVA.</p>
<p>利用固定区组模型，我们可以使用不太简单的线性预测器版本，以解决可估性问题：<span class="math inline">\(\eta+\alpha_i+\beta(\alpha)_{ij}+\rho(\alpha)_{ik}\)</span>，其中 <span class="math inline">\(\alpha\)</span> 表示组效应，<span class="math inline">\(\beta(\alpha)\)</span> 表示 set 内的处理，<span class="math inline">\(\rho(\alpha)\)</span> 表示区组效应，现在视为嵌套在组内。实际上，大多数固定区组拥护者会承认，在这种设计背景下，<span class="math inline">\(\rho(\alpha)\)</span> 是一个整区效应 (whole-plot effect) ——即关于组的重复单元——因此是随机的。然而，我们应该注意到，固定效应拥护者倾向于使用普通最小二乘软件，该软件在求解估计方程、确定标准误以及（对于本次讨论而言最重要的）评估可估性时将 <span class="math inline">\(\rho(\alpha)\)</span> 视为固定的。</p>
<p>上述修复对于均衡数据“起效”（在一定程度上）。现在让我们考虑一个不均衡的多水平设计。Stroup et al. (2018) 阐释了裂区实验的分析，其中整区是不完全区组设计，其结构与上一节中所示的类似。以下显示了分区安排的基本特征。</p>
<div class="inline-figure"><img src="figure/figure%20281-1.png#center" style="width:70.0%"></div>
<p>因素 A 的水平被随机分配给贯穿每个区组的行。因素 B 的水平随机分配给区组中每行内的两个单元格，即在因素 A 的每个重复单元内部。用框架 ANOVA 的术语来说，我们有</p>
<div class="inline-figure"><img src="figure/figure%20281-2.png#center" style="width:100.0%"></div>
<p>其中 <span class="math inline">\(r\)</span> 表示区组数量：如果每种类型使用超过一个区组的话，则 <span class="math inline">\(r\)</span> 可以是 3，也可以是 6, 9 等。标准的实验设计教科书为这种设计给出的线性预测器为：</p>
<p><span class="math display">\[\eta_{ijk}=\eta+\alpha_i+\beta_j+(\alpha\beta)_{ij}+\rho_k+(\alpha\rho)_{ik}\]</span></p>
<p>如果我们把区组——即 <span class="math inline">\(\rho_k\)</span> 和 <span class="math inline">\((\alpha\rho)_{ik}\)</span> ——当作固定效应，那么所有的 A 和 A × B 边际均值都将变得不可估，所有 A 的差异也将变得不可估，以及所有涉及 A 的不同水平的 A × B 的差异——这包括给定 B 的 A 的简单效应——也将变得不可估。这提供了一个明显虚假的不可估的例子。没有人会争辩说 <span class="math inline">\((\alpha\rho)_{ij}\)</span> 应该视为固定效应，因为它是关于因子 A 的重复单元效应。</p>
<p>问题通常出现在使用为普通最小二乘法设计的软件（例如 SAS<sup>®</sup> PROC GLM）来计算此类设计的分析时。PROC GLM 无法区分固定效应和随机效应，因此不恰当地评估了边际均值、差异及所有其他可估函数，就好像 <span class="math inline">\((\alpha\rho)_{ik}\)</span> 是固定的。这就是为什么混合模型（包括所有多水平模型）的计算只能使用专门为混合模型分析设计的软件来完成。</p>
</div>
</div>
<div id="sec9-4" class="section level2" number="9.4">
<h2>
<span class="header-section-number">9.4</span> 使用多水平设计<a class="anchor" aria-label="anchor" href="#sec9-4"><i class="fas fa-link"></i></a>
</h2>
<p>在前面的章节中，我们定义了多水平结构，回顾了使用框架 ANOVA (WWFD?) 过程来识别多水平结构，并仔细研究了区组效应的适当指定对多水平分析的影响。处理完所有这些预备知识后，我们现在看看使用多水平模型的实际方面。我们首先通过多水平分析的示例开始。</p>
<div id="sec9-4-1" class="section level3" number="9.4.1">
<h3>
<span class="header-section-number">9.4.1</span> 多水平结构示例<a class="anchor" aria-label="anchor" href="#sec9-4-1"><i class="fas fa-link"></i></a>
</h3>
<p>在本节中，我们将考虑三个多水平分析的示例——指定模型、将模型转换为所需的 GLIMMIX 语句并解释相关输出。在这里，我们只关注高斯数据。在未来的章节中，我们将考虑这些相同的设计，但使用非高斯数据。</p>
<div class="rmdnote">
<div class="example">
<p><span id="exm:ex9-1" class="example"><strong>示例 9.1 (嵌套析因结构) </strong></span><br></p>
<p>此示例的数据显示在 Data Set 9.1 中。该设计遵循我们在上一节中考虑的图 2.6 的结构——朴素地称其为不连通设计，但更准确地称为嵌套析因。该模型的完整指定为：</p>
<ul>
<li><p>线性预测器：<span class="math inline">\(\eta_{ijk}=\eta+\alpha_i+\beta(\alpha)_{ij}+r(\alpha)_{ik}\)</span>，其中 <span class="math inline">\(r\left(\alpha\right)\)</span> 表示组内的区组。</p></li>
<li>
<p>分布：</p>
<ul>
<li><span class="math inline">\(y_{jk}\mid r\left(\alpha\right)_{ik}\sim N\left(\mu_{ijk},\sigma^2\right)\)</span></li>
<li><span class="math inline">\(r\left(\alpha\right)_{ik}\text{ iid }N\left(0,\sigma_R^2\right)\)</span></li>
</ul>
</li>
<li><p>连接：恒等，<span class="math inline">\(\eta_{ijk}=\mu_{ijk}\)</span></p></li>
</ul>
<p>GLIMMIX 语句为</p>
<pre class="sas"><code>proc glimmix data=ten_blk_six_trt;
class block set trt;
model y = set trt(set);
random block(set);
lsmeans set trt(set) / slice=set slicediff=set;</code></pre>
<p>相关输出：</p>
<div class="inline-figure"><img src="figure/figure%20283-1.png#center" style="width:50.0%"></div>
<p>对于区组，方差分量估计值为 <span class="math inline">\(\hat\sigma^2_R=60.55\)</span>，对于以区组效应为条件的观测，方差分量估计为 <span class="math inline">\(\hat\sigma^2=22.75\)</span>。</p>
<div class="inline-figure"><img src="figure/figure%20283-2.png#center" style="width:60.0%"></div>
<p>“trt(set)” 的 F 值检验 <span class="math inline">\(H_0\colon\)</span> 所有 <span class="math inline">\(\beta(\alpha)_{ij} = 0\)</span>，我们可以更有意义地表示为 <span class="math inline">\(H_0\colon\eta_{00}=\eta_{01}=\eta_{02}\)</span> 以及 <span class="math inline">\(\eta_{13}=\eta_{14}=\eta_{15}\)</span>。它的 <span class="math inline">\(p\)</span> 值为 0.0089，这告诉我们 1) 检验 “set” 的主效应构成了一个“无趣的假设”，2) 由于 “trt(set)” 是一个综合假设，我们需要将其分成几个组成部分，<span class="math inline">\(H_0\colon\eta_{00}=\eta_{01}=\eta_{02}\)</span> 以及 <span class="math inline">\(\eta_{13}=\eta_{14}=\eta_{15}\)</span>，使用 <code>SLICE</code> 选项以获得信息丰富的解释。</p>
<div class="inline-figure"><img src="figure/figure%20283-3.png#center" style="width:60.0%"></div>
<p>这告诉我们，我们有证据表明在组 1（<span class="math inline">\(p=0.0017\)</span>）内有处理效应，但在组 0（<span class="math inline">\(p=0.9709\)</span>）内没有。查看处理均值就能明白这一结果。</p>
<div class="inline-figure"><img src="figure/figure%20283-4.png#center" style="width:50.0%"></div>
<p>最后，SLICEDIFF 列表显示了组内处理的成对简单效应。</p>
<div class="inline-figure"><img src="figure/figure%20283-5.png#center" style="width:80.0%"></div>
<p>我们看到，组 1 中的处理效应来源于统计上显著的差异 <span class="math inline">\(\eta_{13}-\eta_{15}\)</span> 和 <span class="math inline">\(\eta_{14}-\eta_{15}\)</span>（<span class="math inline">\(p\)</span> 值分别为 0.0007 和 0.0048）。换句话说，处理 5 的均值与处理 3 和 4 的均值存在显著差异。</p>
</div>
</div>
<div class="rmdnote">
<div class="example">
<p><span id="exm:ex9-2" class="example"><strong>示例 9.2 (不完全条区 (strip-plot)) </strong></span><br></p>
<p>此示例的数据在 SAS Data and Program Library 中显示为 Data Set 9.2。处理设计为 3 × 3 析因。地形设计有九个区组，每个区组由两行两列、总共四个单元组成。因素 A 的水平随机分配给每个区组的行，但有以下限制：三个区组仅接收水平 A1 和 A2；三个区组仅接收 A1 和 A3；其他三个区组仅接收 A2 和 A3。因素 B 的水平随机分配给每个区组的列，但有以下限制：接收给定一对 A 水平的区组之一必须接收 B1 和 B2，一个必须接收 B1 和 B3，另一个必须接收 B2 和 B3。
因此，我们的多水平结构是二维不完全区组设计。</p>
<p>框架 ANOVA 过程可以帮助我们直观地看到所需的线性预测器：</p>
<div class="inline-figure"><img src="figure/figure%20284.png#center" style="width:100.0%"></div>
<p>请注意，虽然有两个随机化过程，但实际上存在三种尺寸的重复单元：用于 A 水平的行、用于 B 水平的列以及用于 A × B 组合的行列交叉。这些数据的模型为：</p>
<ul>
<li><p>线性预测器：<span class="math inline">\(\eta_{ijk}=\eta+\alpha_i+\beta_j+(\alpha\beta)_{ij}+r_k+(ar)_{ik}+(br)_{jk}\)</span></p></li>
<li>
<p>分布：</p>
<ul>
<li><span class="math inline">\(y_{ijk}\mid r_k,\left(ar\right)_{ik},\left(br\right)_{jk}\thicksim NI\left(\mu_{ijk},\sigma^2\right)\)</span></li>
<li><span class="math inline">\(r_k\mathrm{~iid~}N\left(0,\sigma_R^2\right)\)</span></li>
<li><span class="math inline">\(\left(ar\right)_k\mathrm{~iid~}N\left(0,\sigma_{AR}^2\right)\)</span></li>
<li><span class="math inline">\(\left(br\right)_k\mathrm{~iid~}N\left(0,\sigma_{BR}^2\right)\)</span></li>
</ul>
</li>
<li><p>连接：恒等，<span class="math inline">\(\eta_{ijk}=\mu_{ijk}\)</span></p></li>
</ul>
<p>该模型的 GLIMMIX 语句，包括激活 ODS 图形以使用 LSMEANS 语句 <code>MEANPLOT</code> 选项生成交互图，为：</p>
<pre class="sas"><code>ods graphics on;
proc glimmix data=ds8d2;
class block a b;
model y=a|b / ddfm=kr;
random intercept a b / subject=block;
lsmeans a*b / plot=meanplot(sliceby=a join) slicediff=a;</code></pre>
<p>与我们对不均衡混合模型的标准误和检验统计量偏差的讨论一致，我们使用 <code>DDFM=KR</code> 来调用 Satterthwaite 自由度近似和 Kenward-Roger 偏差校正。交互图如下所示：</p>
<div class="inline-figure"><img src="figure/figure%20285-1.png#center" style="width:60.0%"></div>
<p>我们没有看到 B 效应与 A3 结合的视觉证据。对于 A1 和 A2，B1 的均值似乎不同于 B2 和 B3 的均值。由 A1, A2, B2 和 B3 组成的 A-B 组合在视觉上略有不同。下面的正式统计量可以说明差异是否具有统计学意义，但主题专家需要确定这些差异是否大到足以产生影响，尽管有正式统计量。</p>
<p>方差分量估计和总体效应检验为：</p>
<div class="inline-figure"><img src="figure/figure%20285-2.png#center" style="width:50.0%"></div>
<div class="inline-figure"><img src="figure/figure%20286-1.png#center" style="width:40.0%"></div>
<p>A × B 的 <span class="math inline">\(F\)</span> 统计量提供了交互作用的有力证据（<span class="math inline">\(p =0.007\)</span>），与交互作用图中的视觉证据一致。此时，推断应关注简单效应。因此，我们看一下 LSMEANS 和简单效应的结果：</p>
<div class="inline-figure"><img src="figure/figure%20286-2.png#center" style="width:30.0%"></div>
<div class="inline-figure"><img src="figure/figure%20286-3.png#center" style="width:100.0%"></div>
<p>这些结果为上面给出的视觉特征提供了正式支持：在 A3 中，B 的差异均未达到统计显著性；均值范围在 18.47 到 21.13 之间，标准误为 2.19。在 A1 和 A2 中，B2 和 B3 均值与 B1 不同，但彼此没有差异，B1 与 B3 的差异大于 B2 与 B3 的差异。我们可以继续。例如，B1 与 A1 和 A2 组合的均值似乎与 A3 处所有 B 均值的集合相似。我们将此作为练习。我们在这里看到的是分析的基本流程。</p>
<p>请注意，对处理设计本身的推断与对任何其他析因处理设计的推断没有区别。唯一的区别在于正确定义模型的随机效应以与设计结构保持一致。我们这样做是为了确保使用正确的标准误和检验统计量。</p>
</div>
</div>
<div class="rmdnote">
<div class="example">
<p><span id="exm:ex9-3" class="example"><strong>示例 9.3 (不完全区组的响应曲面设计) </strong></span><br></p>
<p>在 <a href="chap8.html#sec8-6-1">8.6.1</a> 节中，我们为所有因素的定量处理设计引入了二阶多项式回归。回想一下，该模型的吸引力部分源于使用高效、不完全析因设计的能力。在响应面应用中，一种流行的设计是 Box-Behnken 设计 (Box and Behnken, 1960)。 此示例使用三因素 Box-Behnken 设计——请参阅响应曲面教科书，例如 Myers, Montgomery and Anderson-Cook (2016) 对此和其他响应曲面设计进行了全面描述。每个因素都有三个水平——正如响应曲面设计和分析的常见做法，我们将三个水平编码为 -1, 0 和 1。我们通过将每个因素的 水平 0 与来自其他两个因素的水平 -1 和 1 组成的 2<sup>2</sup> 析因设计结合形成处理组合来构建设计。由于我们这里有三个因素，这给我们提供了 12 个处理组合，每组包含四个处理组合。除此之外，我们添加了一种额外的处理组合，其中所有因素的水平均为 0，在设计术语中称为“中心点” (“center point”).</p>
<p>本例的数据在 SAS Data and Program Library 中显示为 Data Set 9.3. 这些数据的多水平性源于使用的不完全区组设计。与特定因素水平 0 相关联的每组四个处理组合看起来都在一起，每组两个区组，每个区组尺寸为 4. 该设计有六个这样的区组，每组两个。第七个区组由四个全部接收中心点的单元组成，完成了整个设计。请注意，存在两个随机化阶段：第一个阶段随机地将组分配给区组；第二个阶段随机地将组内的处理组合分配给区组内的单元。这种两阶段随机化使得这成为了一个多水平设计。同时请注意，这些多水平并不完全对应于设计中的单一因素（就像<a href="chap9.html#exm:ex9-1">示例 9.1</a> 和 <a href="chap9.html#exm:ex9-2">9.2</a> 那样。此处呈现这个例子正是为了表明多水平结构可以采取多种形式。建模从业者需要提高警惕。这个例子也突出了依赖食谱的弱点。学习思考这些事情比学习一个足够全面的食谱来涵盖所有建模意外情况更容易——假设你能设计出这样的食谱！</p>
<p>我们将如何使用 “WWFD?” 框架 ANOVA 过程来帮助构建线性预测器？这是：</p>
<div class="inline-figure"><img src="figure/figure%20287.png#center" style="width:80.0%"></div>
<p>函数 <span class="math inline">\(f\left(A,B,C\right)=\beta_0+\sum_{i=A,B,C}\beta_iX_i+\sum_{i=A,B,C}\beta_{ii}X_i^2+\sum_{i\neq i'=A,B,C}\beta_{ii'}X_iX_i'\)</span>，因素 A, B 和 C 的二阶响应曲面模型，其中 <span class="math inline">\(\beta\)</span>s 表示回归系数，<span class="math inline">\(X_i;i=A,B,C\)</span> 表示三个因素的定量水平。</p>
<p>这些数据的模型可写为：</p>
<ul>
<li><p>线性预测器：<span class="math inline">\(\eta_{ijkl}=\eta_{ijk}+r_{k}\)</span>。其中 <span class="math inline">\(r_l\)</span> 表示区组效应，并且 <span class="math inline">\(\eta_{ijk}=f\left(A,B,C\right)+\left(\alpha\beta\gamma\right)_{ijk}=\beta_0+\sum_{i=A,B,C}\beta_iX_i+\sum_{i=A,B,C}\beta_{ii}X_i^2+\sum_{i\neq i'=A,B,C}\beta_{ii'}X_iX_{i'}+\left(\alpha\beta\gamma\right)_{ijk}\)</span>。其中 <span class="math inline">\(\left(\alpha\beta\gamma\right)_{ijk}\)</span> 为总和项（欠拟合项）。</p></li>
<li>
<p>分布：</p>
<ul>
<li><span class="math inline">\(y_{ijkl}\mid r_l\sim N\left({\mu}_{ijkl},{\sigma}^2\right)\)</span></li>
<li><span class="math inline">\(r_l\mathrm{~iid~}N\left(0,{\sigma}_R^2\right)\)</span></li>
</ul>
</li>
<li><p>连接：恒等，<span class="math inline">\(\eta_{ijkl}=\mu_{ijkl}\)</span></p></li>
</ul>
<p>GLIMMIX 语句分三个步骤出现。在步骤 1 中，我们评估模型的欠拟合情况。假设没有欠拟合的证据，我们拟合完整的二阶响应曲面，从线性预测器中删除 <span class="math inline">\(\left(\alpha\beta\gamma\right)_{ijk}\)</span>。如果有看起来可以忽略不计的二阶项，我们可以删除它们并运行第三步以拟合更简约的模型。步骤 1 的 GLIMMIX 语句为：</p>
<pre class="sas"><code>proc glimmix;
class block aa bb cc;
model y=a|a|b|b|c|c@2 aa*bb*cc/ htype=1 ddfm=kr;
random intercept / subject=block;</code></pre>
<p>注意，我们将 AA, BB 和 CC 定义为 CLASS 变量，以形成欠拟合项 <code>AA*BB*CC</code>。我们在 <a href="chap8.html#sec8-6">8.6</a> 节介绍定量析因处理设计时使用了这个技巧。此外，我们在这里使用序贯检验策略，因为我们想在拟合 <span class="math inline">\(f(A,B,C)\)</span> 后评估欠拟合。因为我们有一个不完全区组设计，我们使用 Kenward-Roger 调整。相关输出：</p>
<div class="inline-figure"><img src="figure/figure%20288.png#center" style="width:80.0%"></div>
<p>特别是，因为此时我们正在检验欠拟合项，所以此时我们只关心 <code>AA*BB*CC</code> 的结果。<span class="math inline">\(F = 0.22\)</span> 和 <span class="math inline">\(p = 0.8776\)</span> 表明没有欠拟合的证据。我们继续执行步骤 2。</p>
<p>步骤 2 仅需要更改 CLASS 和 MODEL 语句。我们从 MODEL 中删除欠拟合项，并在 CLASS 语句中删除相应的效应。新的语句为：</p>
<pre class="sas"><code>class block;
model y=a|a|b|b|c|c@2 / htype=1 ddfm=kr;</code></pre>
<p>相关输出：</p>
<div class="inline-figure"><img src="figure/figure%20289.png#center" style="width:60.0%"></div>
<p>请注意线性-线性项 <code>A*B</code> 和二次项 <code>C*C</code>。它们的 <span class="math inline">\(F\)</span> 值均小于 1。此外，线性-线性项 <code>A*C</code> 的 <span class="math inline">\(p\)</span> 值为 0.2869。这让人怀疑这三个项是否应该保留在模型中。<code>C</code> 的线性主效应也有一个我们通常认为不显著的 <span class="math inline">\(p\)</span> 值。然而，线性-线性 <code>B*C</code> 项具有统计显著性：包含 <span class="math inline">\({\beta}_{BC}X_BX_C\)</span> 但不包含 <span class="math inline">\(\beta_C X_C\)</span> 的回归模型违反了序贯检验的精神，并且会产生因素 B 线性和二次主效应的解释问题。因此，因此，根据步骤 2 得出的合理步骤 3 简化模型将是：</p>
<p><span class="math display">\[f\left(A,B,C\right)=\beta_0+\sum_{i=A,B,C}\beta_iX_i+\sum_{i=A,B}\beta_{ii}X_i^2+\beta_{BC}X_BX_C\]</span></p>
<p>针对步骤 3 修改的 GLIMMIX 语句是：</p>
<pre class="sas"><code>proc glimmix;
class block;
model y=a|a b|b b|c/ solution htype=1 ddfm=kr;
random intercept / subject=block;</code></pre>
<p>相关输出：</p>
<div class="inline-figure"><img src="figure/figure%20290-1.png#center" style="width:60.0%"></div>
<div class="inline-figure"><img src="figure/figure%20290-2.png#center" style="width:80.0%"></div>
<p>我们可以使用这些回归系数生成响应曲面估计图，执行此操作的 SAS 程序是：</p>
<pre class="sas"><code>data rs3;
do a=-1 to 1;
do b=-1 to 1;
do c=-1 to 1;
y_hat=50.085+1.60*a+1.694*b+0.519*c-3.303*a*a-3.515*b*b-1.163*b*c;
output;
end;
end;
end;
proc sort data=rs3;
by c;
proc g3d;
by c;
plot a*b=y_hat/xticknum=3 yticknum=3;</code></pre>
<p>这将创建一个三维图，显示每个 C 水平在 A 和 B 水平上的预测响应：</p>
<div class="inline-figure"><img src="figure/figure%20291.png#center" style="width:100.0%"></div>
<p>我们可以看到，因子 C 对响应曲面的形状有轻微影响。整体图像是一个在 A 和 B 水平上略微倾斜的二维二次曲面，最大响应出现在中心点附近。本示例的要点是，准确拟合曲面需要对欠拟合以及要保留哪些效应做出适当的决定。除非我们正确地指定多水平误差结构，否则我们无法做到这一点。这个例子说明了多水平混合模型和响应曲面方法是兼容的——在这样的例子中，我们必须同时使用这两种方法。</p>
<p>此处未显示但值得一提的另一种非标准多水平结构是 Milliken and Johnson (2008) 提出的混合裂区 (mixed-up split-plot). 一个由多州组成的联合体进行了一项多因素研究。该方案要求采用裂区设计。这一部分方案内容得到了传达，但关于哪个因素被分配给整区（较大的重复单元）以及哪个因素被分配给裂区（较小的重复单元）的部分却没有传达。因此，对于不同州的数据，因素 A 和 B 的位置被颠倒了。将这种结构视为不完全区组设计——类似于不连通/嵌套析因示例——揭示了其分析实际上是一种直接但不常见的多水平建模应用。</p>
</div>
</div>
</div>
<div id="sec9-4-2" class="section level3" number="9.4.2">
<h3>
<span class="header-section-number">9.4.2</span> 多因素处理和多水平设计结构——它们如何一起拟合<a class="anchor" aria-label="anchor" href="#sec9-4-2"><i class="fas fa-link"></i></a>
</h3>
<p>上一节中的三个示例都使用了特定的处理结构。事实上，我们在将多水平模型的地形方面与处理方面相结合的方式上具有相当大的灵活性。</p>
<p>考虑上一节中的<a href="chap9.html#exm:ex9-1">示例 9.1</a>。我们使用了线性预测器 <span class="math inline">\(\eta_{ijk}=\eta+\alpha_i+\beta(\alpha)_{ij}+r(\alpha)_{ik}\)</span>。回顾我们在第 <a href="chap8.html#chap8">8</a> 章中考虑的处理结构的类型，如果处理设计的具体细节需要，我们可以将处理分量改为其他形式。例如，如果每组内的处理水平实际上都是相同的，我们的处理设计实际上是 2 × 3 析因设计，我们可以相应地修改线性预测器，用 <span class="math inline">\(\eta_{ijk}=\eta+\alpha_{i}+\beta_{j}+\left(\alpha\beta\right)_{ij}+r\left(\alpha\right)_{ik}\)</span> 替换上面的式子。如果每组中的处理水平是定量的，并且可以修正为多项式回归，我们可以使用 <span class="math inline">\(\eta_{ijk}=\beta_{0i}+\beta_{1i}X_{j}+\beta_{2i}X_{j}^{2}+r(\alpha)_{ik}\)</span>。相反，如果我们有一个更加非线性的形状，我们可能更喜欢线性化的 Hoerl，<span class="math inline">\(\eta_{ijk}=\beta_{0i}+\beta_{1i}X_{j}+\beta_{2i}\log\left(X_{j}\right)+r\left(\alpha\right)_{ik}\)</span>。注意，所有这些都是 <span class="math inline">\(\eta_{ijk}=\eta_{ij}+r\left(\alpha\right)_{ik}\)</span> 的变体。嵌套、析因和回归模型只是表示编写 <span class="math inline">\(\eta_{ij}\)</span> 以适应处理设计的不同方式。为多水平设计结构建模的组件 <span class="math inline">\(r(\alpha)_{ik}\)</span> 从不更改。请注意，它与关于 <span class="math inline">\(\eta_{ij}\)</span> 的所有变体都兼容。</p>
<p>举例说明，让我们将另一种均值模型应用于<a href="chap9.html#exm:ex9-2">示例 9.2</a>。这些数据提供了比 Data Set 9.4.1 更有趣的演示，因为增加了复杂性：示例 9.1 是基本的均衡裂区；示例 9.2 是不完全条区，具有关于因素 A 水平、因素 B 水平和 AB 处理组合的不同重复单元。然而，我们可以改变处理模型，就像我们对<a href="chap9.html#exm:ex9-1">示例 9.1</a> 所做的那样。也就是说，线性预测器的一般形式 <span class="math inline">\(\eta_{ijk}=\eta_{ij}+r_{k}+\left(ar\right)_{ik}+\left(br\right)_{jk}\)</span> 划分为处理部分 <span class="math inline">\(\eta_{ij}\)</span> 和多水平设计（或“地形”）部分 <span class="math inline">\(r_k+\left(ar\right)_{ik}+\left(br\right)_{jk}\)</span>。我们可以用我们见过的任何方式更改 <span class="math inline">\(\eta_{ij}\)</span>，同时保持 <span class="math inline">\(r_k+\left(ar\right)_{ik}+\left(br\right)_{jk}\)</span> 不变。</p>
<p>我们首先将处理方面表示为嵌套模型，即 <span class="math inline">\(\eta_{ij}=\eta+\alpha_i+\beta(\alpha)_{ij}\)</span>。实现该模型的 GLIMMIX 语句是：</p>
<pre class="sas"><code>proc glimmix data=ds8d2;
class block a b;
model y=a b(a)/ddfm=kr;
random intercept a b / subject=block;
lsmeans b(a) / slice=a slicediff=a;
lsmestimate b(a)‘b1+b2 avg for a1 and a2’ 1 1 0 1 1 0 0 0 0,
‘b1+b2 avg for a3’ 0 0 1 0 0 1 0 0 0,
‘b3 avg for a1 and a2’ 0 0 0 0 0 0 1 1 0 /
divisor=4,2,2;
estimate ‘b1&2 - b3 diff for a1&2’ b(a) 1 1 -2 1 1 -2 0 0 0,
‘b1&2 - b3 diff for a3’ b(a) 0 0 0 0 0 0 1 1 -2 /
divisor=4,2;</code></pre>
<p>我们这样做可能是为了将推断集中在给定 A 的 B 的简单效应上。请注意，使用 LSMESTIMATE 和 ESTIMATE 语句来估计我们早期分析中出现的分组的均值和差异。使用均值模型的嵌套形式可以使这些估计不易出现可估性问题，从而促成这些估计。同时，重要的是，请注意我们完全没有改变随机语句。我们改变了均值模型，将关注点重新集中在我们希望强调的处理分析方面，但我们保留了多水平设计模型。</p>
<p>相关输出（跳过了与之前重复的输出）：</p>
<div class="inline-figure"><img src="figure/figure%20293-1.png#center" style="width:60.0%"></div>
<p>B(A) 检验告诉我们，在 A 的每个水平上，B 因素的均值之间存在差异，但正如我们之前看到的，这种统计量的组合妨碍了有用的解释。最好是查看 SLICE 输出或使用 ESTIMATE 和 LSMEANSTIMATE 输出进行更仔细的观察。</p>
<div class="inline-figure"><img src="figure/figure%20293-2.png#center" style="width:60.0%"></div>
<p>这样是更好的。我们看到在 A1 和 A2 内有显著的 B 效应，但在 A3 内没有。将这些差异分开，我们有：</p>
<div class="inline-figure"><img src="figure/figure%20293-3.png#center" style="width:100.0%"></div>
<p>以及</p>
<div class="inline-figure"><img src="figure/figure%20294.png#center" style="width:100.0%"></div>
<p>LSMESTIMATE 输出分别列出 <span class="math inline">\(\begin{pmatrix}\eta_{11}+\eta_{12}+\eta_{21}+\eta_{22}\end{pmatrix}/4\text{,}\begin{pmatrix}\eta_{31}+\eta_{32}\end{pmatrix}/2\)</span> 和 <span class="math inline">\(\left(\eta_{13}+\eta_{23}\right)/2\)</span> 的估计。ESTIMATE 输出分别列出 <span class="math inline">\(\left[\left(\eta_{11}+\eta_{12}+\eta_{21}+\eta_{22}\right)/4\right]–\left[\left(\eta_{31}+\eta_{32}\right)/2\right]\)</span> 和 <span class="math inline">\(\begin{bmatrix}(\eta_{13}+\eta_{23})/2\end{bmatrix}-\eta_{33}\)</span> 的估计。我们这样做是因为<a href="chap9.html#exm:ex9-2">示例 9.2</a> 的分析表明 A1 和 A2 中 B1 和 B2 的均值相似，因此对这四个均值进行平均是有意义的。B1 和 B2 的均值在 A3 内没有差异，因此对它们进行平均也是有意义的。 B3 的均值与 B1 和 B2 不同，但 B3 与 A1 的组合与 B3 与 A2 的组合没有区别。因此就有了第三个 LSMESTIMATE 语句。ESTIMATE 语句获取这些分组之间的差异估计。当我们想知道差异的大小，而不仅仅是差异是否存在时，我们就会这样做。</p>
<p>让我们更进一步。假设 B 的水平是定量的。在交互图中，均值模式表明可能存在二次回归——当然是在 A1 和 A2 内的 B 水平上。我们可以用定性-定量模型对此进行建模：具体来说，<span class="math inline">\(\eta_{ij}=\beta_{0i}+\beta_{1i}X_j+\beta_{2i}X_j^2\)</span>，其中 <span class="math inline">\(X_j\)</span> 表示因素 B 的第 <span class="math inline">\(j\)</span> 个定量水平。GLIMMIX 语句为：</p>
<pre class="sas"><code>proc glimmix data=ds8d2;
class block a b;
model y=a xb(a) xb_sq(a)/noint s ddfm=kr;
contrast ‘quadratic equal, a1 and a2’ xb_sq(a) 1 -1 0;
contrast ‘quadratic equal, a3 v a1&2’ xb_sq(a) 1 1 -2;
contrast ‘regression equal, a1 and a2’ xb(a) 1 -1 0,
xb_sq(a) 1 -1 0;
contrast ‘regression equal, a3 v a1&2’ xb(a) 1 1 -2,
xb_sq(a) 1 1 -2;
random intercept a b / subject=block;</code></pre>
<p>变量 <code>XB=B</code> 以及 <code>XB_SQ=B*B</code>。我们必须在数据步骤中定义它们。它们允许 GLIMMIX 使用 B 作为 CLASS 变量来指定多水平设计模型，但将 <code>XB</code> 和 <code>XB_SQ</code> 分别用作 <span class="math inline">\(X_j\)</span> 和 <span class="math inline">\(X_j^2\)</span> 的直接变量。请注意，RANDOM 语句保持不变。CONTRAST 语句允许我们检验 A1 和 A2 的二次回归系数是否相等，如果相等，则检验它们是否与 A3 不同。第二组 CONTRAST 将线性和二次项合并为回归的共同检验。<code>NOINT</code> 选项为我们提供了一个满秩模型：因此，SOLUTION 向量提供了实际的回归估计，无需进行进一步的操作。相关输出：</p>
<div class="inline-figure"><img src="figure/figure%20295-1.png#center" style="width:100.0%"></div>
<p>例如，A1 的估计回归方程为 <span class="math inline">\(0.64+26.31X_j-5.71X_j^2\)</span>。这些估计值清楚地显示了 A1 和 A2 的关于 B 的回归的相似性以及 A3 的 B 的回归非常不同。我们还可以看到，A3 下的线性和二次系数在统计上并不显著——更准确地说，对于任何合理的 <span class="math inline">\(\alpha\)</span> 水平，我们无法拒绝 <span class="math inline">\(H_0\colon\beta_{13} = 0\)</span> 和 <span class="math inline">\(H_0\colon\beta_{23} = 0\)</span>。对于 A3，仅使用截距可能会产生更好的拟合效果。我们将在下面继续探讨这一点。</p>
<div class="inline-figure"><img src="figure/figure%20295-2.png#center" style="width:60.0%"></div>
<p>我们在这里展示 “type III tests of fixed effects” 更多的是因为它们没有用，而不是因为它们提供了有意义的输出（它们实际上没有）。此列表中的每一行都是过度组合的检验统计量。这三个自由度来自将截距、线性和二次回归的每个参数在因素 A 的所有水平上分别组合，并同时对其进行检验。将这些结果与上述各个回归系数的结果进行比较，你可以看到 “type III tests” 并没有提供有用的信息。</p>
<p>最后是对比结果：</p>
<div class="inline-figure"><img src="figure/figure%20295-3.png#center" style="width:80.0%"></div>
<p>这些结果证实了我们从其他角度观察到的现象。A1 和 A2 在 B 上的回归是相似的，但与 A3 的回归不同。</p>
<p>最后，我们之前注意到，各个回归系数的检验表明，仅使用截距对 A3 的响应进行建模可能更好，因为 B 似乎在 A3 中没有效应。简化模型将是</p>
<p><span class="math display">\[\eta_{ij}=\beta_{0i}+I\left(A<3\right)\beta_{1i}X_j+I\left(A<3\right)\beta_{2i}X_j^2\]</span></p>
<p>其中 <span class="math inline">\(I(A < 3)\)</span> 是指示函数，对于 A1 和 A2 水平等于 1，对于 A3 等于 0. 这消除了用于 A3 的 B 的回归。为了使用 GLIMMIX 实现这一点，我们需要在数据集中定义指示函数。在下面的程序中，对于 A1，<code>B_1=1</code>，否则为 0；对于 A2，<code>B_2=1</code>，否则为 0；对于 A3，<code>B_3=1</code>，否则为 0。GLIMMIX 语句为：</p>
<pre class="sas"><code>proc glimmix data=ds8d2;
class block a b;
model y=b_1 b_2 b_3 b_1*xb b_2*xb b_1*xb_sq b_2*xb_sq/noint s
ddfm=kr;
contrast ‘intercept equal, a1 and a2’ b_1 1 b_2 -1;
contrast ‘linear equal, a1 and a2’ b_1*xb 1 b_2*xb -1;
contrast ‘quadratic equal, a1 and a2’ b_1*xb_sq 1 b_2*xb_sq -1;
random intercept a b / subject=block;</code></pre>
<p>请注意，MODEL 语句现在根本不包含 CLASS 变量。所有三个截距均出现，并且仅显示 A1 和 A2 的线性项和二次项。注意 CONTRAST 语句的形式。每个效应只有一个系数，因为模型中的每项都是一个参数。</p>
<div class="rmdcaution">
<p>对于学生：你应该能够写出每个对比的可估函数的精确形式。</p>
</div>
<p>输出为：</p>
<div class="inline-figure"><img src="figure/figure%20296.png#center" style="width:80.0%"></div>
<p>第一个表给出了简化模型的系数估计。注意，<span class="math inline">\(\hat\beta_{03} =19.44\)</span> 估计 A3 在所有 B 水平上的预期响应。对比的结果证实，不存在统计上显著的证据来得出 A1 和 A2 的回归参数不同的结论。我们在这里没有展示它，但我们可以通过在 A1 和 A2 上合并 <span class="math inline">\(\beta_{0i},\beta_{1i}\)</span> 和 <span class="math inline">\(\beta_{2i}\)</span> 来进一步简化模型。</p>
</div>
</div>
<div id="sec9-5" class="section level2" number="9.5">
<h2>
<span class="header-section-number">9.5</span> 边际和条件多水平模型<a class="anchor" aria-label="anchor" href="#sec9-5"><i class="fas fa-link"></i></a>
</h2>
<p>我们在第 <a href="chap3.html#chap3">3</a> 章中介绍了边际模型，并在第 <a href="chap5.html#chap5">5</a> 章和第 <a href="chap6.html#chap6">6</a> 章中提供了额外的背景知识。我们已经看到，高斯仅方差分量模型可以重新表述为复合对称边际模型。对于非高斯模型，我们可以对工作的相关矩阵进行类似的处理，但与高斯情况不同，这两个参数化是不等价的。</p>
<p>对于多水平模型，条件与边际的区别在何时重要？这个问题主要在我们有负的方差分量估计时出现。为了解这是如何工作的，我们重新回顾了<a href="chap9.html#exm:ex9-2">示例 9.2</a> 中使用的设计结构，但使用了不同的数据。首先，我们考虑高斯情况。然后，我们简要评论了对于非高斯情况，我们的方法将如何不同。在未来的章节中，我们将更详细地讨论非高斯情况。</p>
<div id="sec9-5-1" class="section level3" number="9.5.1">
<h3>
<span class="header-section-number">9.5.1</span> 高斯数据<a class="anchor" aria-label="anchor" href="#sec9-5-1"><i class="fas fa-link"></i></a>
</h3>
<p>Data Set 9.4 使用与 Data Set 9.2相同的结构，但数据不同。使用我们在<a href="chap9.html#exm:ex9-2">示例 9.2</a> 中使用的相同 GLIMMIX 语句，方差分量估计为：</p>
<div class="inline-figure"><img src="figure/figure%20297-1.png#center" style="width:60.0%"></div>
<p>请注意，因素 A 的重复单元之间的方差估计 <span class="math inline">\(\hat\sigma_{AR}^2\)</span> 已设置为零。当 REML 解向量 <span class="math inline">\(\hat{\symbf\sigma}\)</span> 收敛到满足 <span class="math inline">\(\hat\sigma_{AR}^2<0\)</span> 的值时，就会出现这种情况。由于负解位于参数空间之外，通常的做法是将这些方差分量估计置零。</p>
<p>Stroup and Littell (2002) 发现置零做法会对功效特性和置信区间覆盖范围产生不利影响。建议的修复方法是允许 <span class="math inline">\(\hat\sigma_{AR}^2\)</span> 保持负值。在 GLIMMIX 中，我们通过在 PROC GLIMMIX 语句中包含 <code>NOBOUND</code> 作为选项来实现此目的。使用 <code>NOBOUND</code> 得出的方差分量估计为</p>
<div class="inline-figure"><img src="figure/figure%20297-2.png#center" style="width:60.0%"></div>
<p><code>NOBOUND</code> 如何影响对 A 和 B 处理效应的推断？</p>
<p>不带 <code>NOBOUND</code>：</p>
<div class="inline-figure"><img src="figure/figure%20298-1.png#center" style="width:60.0%"></div>
<p>带 <code>NOBOUND</code>：</p>
<div class="inline-figure"><img src="figure/figure%20298-2.png#center" style="width:60.0%"></div>
<p>我们看到了 <code>NOBOUND</code> 对这些数据的影响：主效应检验统计量因 <code>NOBOUND</code> 而增加；交互检验略有下降。<code>NOBOUND</code> 的影响是特定于设计的——主效应和交互检验不会总是受到这种影响。从全局来看，模拟结果一致表明，<code>NOBOUND</code> 结果在 I 类错误率控制和功效特性方面更准确。我们是如何评估后者的？在第 <a href="#chap21"><strong>??</strong></a> 章中，我们讨论了基于 GLMM 的功效分析。当我们说 <code>NOBOUND</code> 在功效特性方面更准确时，我们的意思是，当处理均值不等时，使用 <code>NOBOUND</code> 的拒绝率比不使用 <code>NOBOUND</code> 获得的拒绝率更符合第 <a href="#chap21"><strong>??</strong></a> 章所示方法的理论预测。</p>
<p>虽然 <code>NOBOUND</code> 产生了更准确的推断，但它给我们留下了一个尴尬的问题，即准确解释负方差的含义。例如，在这里，我们所说的 <span class="math inline">\(\hat\sigma_{AR}^2=-6.595\)</span> 是什么意思？</p>
<p>对于高斯数据，另一种方法是将受影响的方差分量重新表示为复合对称协方差。在我们的示例中，GLIMMIX 语句为：</p>
<pre class="sas"><code>proc glimmix data=ds8d2;
class block a b;
model y=a|b/ddfm=satterth;
random intercept b / subject=block;
random _residual_ / type=cs subject=block*a;</code></pre>
<p>得到的方差分量估计：</p>
<div class="inline-figure"><img src="figure/figure%20298-3.png#center" style="width:50.0%"></div>
<p>请注意，这只是 <code>NOBOUND</code> 结果的重新表达，与第 <a href="chap3.html#chap3">3</a> 章中方差分量模型与复合对称性模型的等价性类似。主要区别在于：在统计理论中，-6.595 的协方差是良定的；而 -6.595 的方差则不是。</p>
</div>
<div id="sec9-5-2" class="section level3" number="9.5.2">
<h3>
<span class="header-section-number">9.5.2</span> 非高斯模型<a class="anchor" aria-label="anchor" href="#sec9-5-2"><i class="fas fa-link"></i></a>
</h3>
<p>正如我们在第 3 章中看到的，边际复合对称模型与条件方差分量模型的等价性并不能扩展到非高斯模型。边际模型的目标是不同的参数——边际分布的均值——而条件模型的目标是 <span class="math inline">\(\symbf y\mid \symbf b\)</span> 的期望值，通常是指数族成员的可识别参数。</p>
<p>考虑到这种不等价性，当我们在非高斯 GLMM 中估计负方差分量时，<code>NOBOUND</code> 是适当的补丁。复合对称工作相关性实际上不能被视为一种可行的替代方案。在接下来的章节中，我们将回到这个主题，通过示例进行说明。例如，在第 <a href="#chap11"><strong>??</strong></a> 章中，我们将详细介绍计数数据的示例。</p>
</div>
</div>
<div id="sec9-6" class="section level2" number="9.6">
<h2>
<span class="header-section-number">9.6</span> 总结<a class="anchor" aria-label="anchor" href="#sec9-6"><i class="fas fa-link"></i></a>
</h2>
<ul>
<li>多水平研究定义为至少存在两种随机变异源的研究。</li>
<li>最简单的多水平研究是带子抽样的完全随机设计。</li>
<li>对于高斯数据，多水平研究必须在线性预测器中包含一个随机效应，以解释较大单元的变异。最小单元的变异嵌入给定随机效应的观测的条件分布。</li>
<li>对于单参数非高斯分布，与最小水平相关联的变异也可以包括在线性预测器中。不这样做是过度分散的常见原因，这一问题从第 <a href="#chap11"><strong>??</strong></a> 章开始进行更详细的讨论。</li>
<li>掌握“Fisher 会怎么做？”过程是准确构建多水平模型的有效辅助。</li>
<li>设计结构和处理结构可以混合搭配。除少数例外情况外，没有一种处理结构必然排除某种设计结构，反之亦然。</li>
<li>负方差分量估计影响功效和 I 类错误控制。通常，将负解置零会为检验统计量带来偏差，因此可能是一个问题。</li>
<li>对于高斯数据，允许方差估计保持为负值，虽然在解释方差估计本身时很难解释，但这是保持检验统计量无偏并因此保持对 I 类错误控制的最佳方法。</li>
<li>多水平模型通常可以重新表示为复合对称协方差模型。这样做可以消除负方差估计的问题，并保持检验统计量和标准误的无偏性。对于非高斯 GLMMs，将模型重新表达为复合对称协方差模型变得更加重要，因为 <code>NOBOUND</code> 选项更容易导致估计方程性能不佳。</li>
</ul>
</div>
<div id="exe9" class="section level2 unnumbered">
<h2>练习<a class="anchor" aria-label="anchor" href="#exe9"><i class="fas fa-link"></i></a>
</h2>
<ol style="list-style-type: decimal">
<li>本习题部分使用文件 Ch_8_Problem1.sas 中包含的数据。本研究比较了四种处理在暴露于降解剂后的劣化情况。数据来自七个批次，这些批次是随机选取的批次总体的代表。每个批次都有足够的材料可以分成两半——一种处理可以随机分配给半个批次，另一种处理分配给另一个半个批次。材料暴露于降解剂中；在数据集中，X 表示暴露水平，X 的值等于 1, 2, 4, 8 和 16，值越大，暴露越大。响应变量用 Y 表示；你可以假定它具有高斯分布。另外，众所周知，在本研究观察到的 X 范围内，Y 对 X 的增加呈线性响应。Y 值较低则不佳；Y 值较高则良好。各批次的处理分配如下图所示：
<ol style="list-style-type: lower-alpha">
<li>写出与本研究相关的统计模型的完整描述，与上面给出的描述和假定一致。</li>
<li>设计是连通的吗？请解释。</li>
<li>该研究的目的是确定“最佳”处理。对于研究人员来说，“最佳”意味着对增加降解物暴露带来的不良影响具有最大抵抗力的处理（Y 越大意味着抵抗力越强）。根据你在 a. 中给出的模型，使用 PROC GLIMMIX 分析数据。使用相关结果来实现研究人员的目标。哪种处理最好，为什么？处理如何排名以及为什么？</li>
</ol>
</li>
</ol>
<div class="inline-figure"><img src="figure/figure%20300.png#center" style="width:50.0%"></div>
<ol start="2" style="list-style-type: decimal">
<li>
<p>本习题使用文件 Ch_9_Problem2.sas 中包含的数据。进行了一项多州研究。处理设计为 2<sup>2</sup> 析因。数据集中的因素表示为 A 和 B，各具有两个水平，表示为 0 和 1。原始方案要求每个州使用裂区，其中因素 A 作为整区因素，因素 B 作为裂区因素。州 1、州 2 和州 3 按照原始方案正确执行；州 4 和州 5 弄错了方案：它们使用因素 B 作为全区因素，使用因素 A 作为裂区因素。响应变量是 Y，你可以假定它具有高斯分布。研究人员现在希望对整个数据集进行单一的合并分析。</p>
<ol style="list-style-type: lower-alpha">
<li>编写包含所需元素的统计模型，以实现所需的分析。</li>
<li>使用 GLIMMIX 执行与 a. 一致的分析。写一份简短的报告来描述两种处理的效应。请记住，报告处理效应“大小” (“how big”) 以及“类型” (“what kind”) 与报告效应是否具有统计显著性同样重要。特征描述至少有一部分基于数据尺度进行——因为，非统计学读者会坚持这一点。引用相关证据来支持你的描述。</li>
<li>如果你假定了随机区组效应，如果你将假定更改为固定区组，会发生什么？如果你假定了固定区组，如果你将假定更改为随机区组效应，又会发生什么？“会发生什么”包括在 GLIMMIX 语句中你需要做出哪些更改，以及输出列表中会发生哪些变化？</li>
</ol>
</li>
</ol>
<div class="rmdcaution">
<p>提示：对于 a. 利用这样一个事实：将裂区结构视为不连通不完全区组设计。在整个研究中的所有区组中，给处理组合一个一致的标识。以这种方式写出分区计划并利用“Fisher 会怎么做？”过程应该会产生一个可用的方法。</p>
</div>
<ol start="3" style="list-style-type: decimal">
<li>
<p>数据集 Ch9_Problem3a.sas 和 Ch9_Problem3b.sas 包含来自一个条区实验的数据，该实验有六个区组和 3 × 3 的处理设计—— A 有 3 个水平，B 有 3 个水平。以下是条区中一个区组的示意图。在每个区组中，A 的水平被随机分配到一行，并在整行中应用；B 的水平被随机分配到一列，并在整列中应用。每个区组中的 A 到行和 B 到列的分配都是潜在唯一的。</p>
<ol style="list-style-type: lower-alpha">
<li>显示得出此设计的框架 ANOVA 的 WWFD 过程。</li>
<li>写出 a. 得出的线性预测器。包括关于随机模型效应的假定。</li>
<li>对于以随机模型效应为条件的以下哪种响应分布，在线性预测器中包含与 <code>block*a*b</code> 对应的项是合理的？
<ol style="list-style-type: lower-roman">
<li>高斯</li>
<li>贝塔</li>
<li>二项</li>
<li>泊松</li>
<li>负二项</li>
</ol>
</li>
<li>判断题：若为你在 a. 中编写线性预测器拟合负二项模型，你隐含地假定了 <code>block*a*b</code>（单元水平）效应具有伽马分布。</li>
</ol>
</li>
</ol>
<div class="inline-figure"><img src="figure/figure%20301.png#center" style="width:30.0%"></div>
</div>
</div>
<div class="chapter-nav">
<div class="prev"><a href="chap8.html"><span class="header-section-number">8</span> 处理和解释变量结构</a></div>
<div class="next"><a href="bib.html">参考文献</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="#chap9"><span class="header-section-number">9</span> 多水平模型</a></li>
<li><a class="nav-link" href="#sec9-1"><span class="header-section-number">9.1</span> 设计结构类型：单水平和多水平模型定义</a></li>
<li>
<a class="nav-link" href="#sec9-2"><span class="header-section-number">9.2</span> 多水平模型的类型及其产生方式</a><ul class="nav navbar-nav">
<li><a class="nav-link" href="#sec9-2-1"><span class="header-section-number">9.2.1</span> 重复单元——不仅在设计的实验中</a></li>
<li><a class="nav-link" href="#sec9-2-2"><span class="header-section-number">9.2.2</span> 重温“Fisher 会怎么做？”——地形和处理部分</a></li>
</ul>
</li>
<li>
<a class="nav-link" href="#sec9-3"><span class="header-section-number">9.3</span> 区组在多水平模型中的作用</a><ul class="nav navbar-nav">
<li><a class="nav-link" href="#sec9-3-1"><span class="header-section-number">9.3.1</span> 重温“固定区组效应与随机区组效应”</a></li>
<li><a class="nav-link" href="#sec9-3-2"><span class="header-section-number">9.3.2</span> 固定区组、多水平设计与虚假的不可估性</a></li>
</ul>
</li>
<li>
<a class="nav-link" href="#sec9-4"><span class="header-section-number">9.4</span> 使用多水平设计</a><ul class="nav navbar-nav">
<li><a class="nav-link" href="#sec9-4-1"><span class="header-section-number">9.4.1</span> 多水平结构示例</a></li>
<li><a class="nav-link" href="#sec9-4-2"><span class="header-section-number">9.4.2</span> 多因素处理和多水平设计结构——它们如何一起拟合</a></li>
</ul>
</li>
<li>
<a class="nav-link" href="#sec9-5"><span class="header-section-number">9.5</span> 边际和条件多水平模型</a><ul class="nav navbar-nav">
<li><a class="nav-link" href="#sec9-5-1"><span class="header-section-number">9.5.1</span> 高斯数据</a></li>
<li><a class="nav-link" href="#sec9-5-2"><span class="header-section-number">9.5.2</span> 非高斯模型</a></li>
</ul>
</li>
<li><a class="nav-link" href="#sec9-6"><span class="header-section-number">9.6</span> 总结</a></li>
<li><a class="nav-link" href="#exe9">练习</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>