Play with R 第19期：分层线性模型/多水平

Original 荷兰心理统计联盟荷兰心理统计联盟 2023-02-03

收录于合集 #R语言 62个

Play with R 第19期：

Multilevel linear models

19.1本章将告诉我们什么？（略）

19.2分层数据

19.2.1类内相关

19.2.2分层模型的好处

19.3分层线性模型的理论

19.3.1例子

19.3.2固定系数和随机系数

19.3.2.1随机截距模型

19.3.2.2随机斜率模型

19.3.2.3随机斜率和随机截距模型

19.4多水平模型

19.4.1 评估与比较多水平模型的匹配

19.4.2 协方差结构的类型

19.5. 一些实际问题

19.5.1. 理论假设

19.5.2. 样本量和统计功效

19.5.3. 集中变量 (Centring variables)

19.6.多层模型在R中的应用步骤（以书上的例子为例）

19.7 增长模型

19.7.1 增长模型（指数）

19.7.2 一个例子（略）

19.7.3 重构数据

19.7.4 设置基本模型

19.7.5 添加时间作为固定影响

19.7.6 介绍自由截距项

19.7.7构造协方差结构

19.7.8 模型比较

19.7.9 加入高指数项

19.8 如何报告多水平模型

19.1本章将告诉我们什么？

没有实质性的内容，只是讲作者为什么要加这一章，这一章对他来说也是较新的内容，尽管如此，开始学习也不晚。

19.2分层数据

在现实生活中，数据往往是分层的，如书中图19.2所示是一个最简单的分层模式，考察了一个学校中学生的表现，那么学生是第一层级，来源的班级是第二层级。当然，还有复杂的模式，如图19.3所示，包含了三个层级：学生、班级、学校。像上述数据就称为分层数据。将那些并非主要考察的变量称为背景因素。

19.2.1类内相关

为什么分层数据是重要的呢？因为这意味着数据残差将相关。例如，我们假定抽取的学生之间的表现没有关联，但是如果他们恰巧来自同一班级，并且该班级受老师的影响较大，那么学生之间就具有了相关性。那么，在进行ANOVA时，独立性的缺失会影响统计结果。

19.2.2分层模型的好处

1.抛开回归斜率同质性假设。

如果我们要将那些非关注的变量放入协变量进行协方差分析，我们需要满足协变量与结果之间的关系在组成我们的预测变量的不同组中是相同的。但是，这种情况并不总是会发生。而分层模型不用考虑这一点。

2.对独立性假设说再见。

当使用独立方差分析时，我们必须假设数据的不同情况是独立的然而很可能事实并非如此。此外，多元回归依赖于独立的观察。但是，在某些情况下，您可能想多次测量某人（即随着时间的推移）。那么，分层模型是专门设计用来允许您对案例之间的这些关系建模的。

3.数据丢失时请大笑。

数据丢失或设计不平衡时，回归，方差分析，ANCOVA和我们介绍的其他大多数测试都会产生奇怪的结果。分层模型不需要完整的数据集，因此，如果某个时间点上的数据丢失，则无需插补，也无需删除整个案例。相反，可以使用可用数据成功估计参数，这为处理丢失的数据提供了一个相对简单的解决方案。需要强调的是，没有任何统计方法可以克服丢失的数据。应使用良好的方法，设计和研究执行力，以尽量减少缺失值，并应始终探索缺失值的原因。只是，当使用传统的统计方法进行重复测量数据时，通常需要其他方法来解决丢失的数据，这可能会带来问题，在这种情况下分层模型是更好的选择。

19.3分层线性模型的理论

两个要点：什么是多层线性模型；模型中的主要概念

如果你学过13.14章，那你已经接触过多层线性模型了，重复测量设计就是一个两层的模型（分数是第一层，嵌套在作为第二层的被试里）

19.3.1例子

用例子来具体阐述模型中的概念：整容手术对人们生活质量的影响。

数据：(Cosmetic Surgery.dat)

变量：术后生活质量（结果变量）；术前生活质量（调试变量）；手术与否（用于分组，控制组0还没做过整容手术和实验组1做了整容手术）；手术与否文本（分组更详细的说明）；手术诊所（10个）；年龄；BDI（贝克抑郁量表，）；手术原因（1.外貌改变2躯体不适原因）；手术原因文本（定义不同的组别）；性别

分层模型：

水平1——不同诊所的整容者。

水平2——不同的整容诊所。

分层依据：由同一个外科医生操刀做手术的整容者之间显然不是相互独立的。比起在不同诊所手术的整容者，同一诊所的整容者有更多相似的特征。

多层模型图示19.5

19.3.2固定系数和随机系数

在本书的所有回归分析中，我们假设所有回归参数都是固定参数。我们经常见到的线性模型是用两个参数来表示的：截距b0 和斜率b1 ：（Y为因变量，X为自变量，ε为误差，都随函数i变化，i代表不同的被试，即水平1不同的变量）

当我们分析一个这样的回归，我们假设b是固定的，根据数据可以计算出来。基于这个假设，我们默认这个模型适用于全部样本，因此我们可以根据同样的截距和斜率来算出样本中的每个数据。

然而，我们也可以将这些参数设置为随机参数，把他们假设为可以变化的值。

我们以11章的例子（研究人的性欲与其伴侣的性欲（图中以虚线表示）的关系，组别有安慰剂组，低剂量组和高剂量组），演示固定和随机参数搭配，图示19.6：

19.3.2.1随机截距模型

最简单的将随机参数引进模型的方式是假设在不同组中的截距是不固定的（随机截距）。在研究性欲的例子中，即假设在不同组别（安慰剂组，低剂量组，高剂量组）的被试的与伴侣的性欲关系是相同的情况下，他们的模型是在不同的起点的（即截距不同）见图19.6的第一个图。

19.3.2.2随机斜率模型

同时我们也可以假设在不同组中的斜率是不固定的（随机斜率）

在研究性欲的例子中，即假设他们的模型是在相同的起点（即截距相同），不同组别（安慰剂组，低剂量组，高剂量组）的被试的与伴侣的性欲关系是不同的（即斜率不同），见图19.6的第二幅图。这就触犯了在ANCOVA中不同情境下回归斜率一致性的假设。当这个回归斜率一致性的假设站不住脚时，用分层模型就能更清晰地估计斜率的变化。需要注意的是，现实中很少出现截距固定，斜率随机的模型，因为斜率的变化通常会引起关系整体的变化，即截距的变化，所以通常我们要假设斜率随机时，我们会同时假设截距随机。

19.3.2.3随机斜率和随机截距模型

这是现实生活中最常见的情况，模型中的斜率和截距都是变化的。

19.4多水平模型

我们已经从概念上了解有关随机截距、随机斜率、随机截距和斜率模型，接下来可以通过具体的实例来看如何操作。我们要预测人们外科手术后的生活质量QoL，可以将其表示为线性模型：QoLAfterSurgeryi= b0 +b1Surgeryi+εi

其中，模型第二层中有一个情境变量为实施手术的不同诊所，我们可以允许代表手术对生活质量影响的模型在不同的环境(诊所)中发生变化，即允许截距在诊所之间变化，或者允许斜率在诊所之间变化，或者允许两者在诊所之间变化。因此，我们在模型截距上增加一个随机变化量u0j组成总模型的截距估计(b0+ u0j),总模型方程变为Yij = (b0 + u0j ) +b1Xij+ εij

方程中的js反映了截距变化的变量水平(本例中为诊所)的水平，即第二层变量。另一种写法是把误差项去掉，它看起来像一个普通的回归方程，而截距从固定的b0变成了随机的b0j：

因此，如果我们想知道诊所7的估计截距，我们只需将第二个方程中的j替换为“诊所7”即可:

b0Clinic7 = b0 +u0Clinic7

如果我们想要包含评估手术对生活质量的影响的随机斜率，需要在整个模型的斜率上增加一个成分来测量斜率的变化程度，u1j。因此，倾斜度由b1变为(b1 + u1j)。该术语估计了整个模型的斜率(b1)，以及整个模型(u1j)在不同情境中斜率的变化程度。总模型变为

Yij = b0 + (b1 + u1j )Xij + εij

同样，我们可以把误差项放到一个单独的方程中，使它与我们熟悉的线性模型的联系起来。它现在看起来像一个普通的回归方程，而斜率从一个固定的b1变成了一个随机的b1j：

综上，我们通过具体的随机斜率和截距项来表示整个模型：

Yij = (b0 + u0j ) + (b1 + u1j )Xij + εij

同样地，我们可以把整个模型拆分成几个具体的线性模型：

在此基础上，我们需要考虑增加另一个预测变量，例如术前生活质量，于是模型方程变为QoLAfter Surgeryi = b0 +b1Surgeryi+b2QoLBefore Surgeryi +εi

在此模型中请记住，i表示第一层变量，在本例中是我们的被试。因此，我们可以预测一个人在手术后的生活质量，用他们的名字代替i:

QoL AfterSam = b0 +b1SurgerySam+b2QoLBeforeSam +εSam

现在，如果我们想让手术对术后生活质量的影响在不同的情况下有所不同那么我们就简单地用b0j代替b0。如果我们想要让手术对术后生活质量影响的斜率在不同的情况下变化那么我们就用b1j代替b1:

请记住，方程中的j与第二册情境变量(在本例中为诊所)有关。所以，如果我们想要预测一个人的分数，不仅要从其名字来预测，还要从其所去的诊所来预测。假设Sam在7号诊所做手术，然后我们可以替换如下的is和js：

QoLAfter SurgerySam,Clinic 7

=b0Clinic7 +b1Clinic7SurgerySam, Clinic7 +b2QoLBefore SurgerySam, Clinic7 +εSam, Clinic7

综上，我们在多水平模型中所做的实际上是一种特殊的回归，即允许截距或者斜率，或者两者，在不同的情境中变化。其实质改变的是对于每一个允许随机变化的参数，我们得到了一个参数的可变性的估计值和参数本身。

19.4.1 评估与比较多水平模型的匹配

正如在logistic回归(第8章)中，多水平模型的整体拟合采用卡方似然比检验(见第18.4.3节)和R报告−2log-likelihood (-2LL，见第8.3.1节)。本质上，log-likelihood的值越小越好。R还产生了两个调整版本的对数似然值，已在第7.6. 3节中简要描述。这两种情况都可以用与log-likelihood相同的方式来解释，但是它们在很多方面都被已修正。

Akaike’s information criterion (AIC): 这基本上是对拟合优度的一个度量，它根据模型的复杂性进行修正，其已考虑了那些已经估计的参数的数量。

Schwarz’s Bayesian criterion (BIC): 这个统计数据与AIC是可比较的，尽管它稍微保守一些(它对被估计参数的数量进行了更严格的校正)。当样本量较大且参数数量较少时，应采用该方法。

一般情况下，较小的AIC和BIC参数值意味着模型的适合度较高。许多作者建议从基本模型开始建构一个所有参数都为固定的，并适当加入随机系数和探索confounding variables 的多水平模型 (Raudenbush &Bryk, 2002;Twisk,2006) 。这样做的一个好处是可以在使参数随机化或添加变量时有效比较模型的适合度。为了比较模型，我们只需从原模型的值中减去新模型的log-likelihood。

19.4.2 协方差结构的类型

如果在多水平模型中存在随机效应或重复测量，就必须确定数据的协方差结构。如果你有随机效应和重复测量，那么你可以为其指定不同的协方差结构。协方差结构简单地说明了方差-协方差矩阵的形式(矩阵中对角线元素是方差，非对角线元素是协方差)。这个矩阵有多种形式，我们需要告诉R我们所认为它是什么形式。当然,我们可能并不知道它的准确形式，因此采用不同定义的协方差结构和使用拟合优度指数(AIC和BIC)来看是否改变模型的协方差结构可以提高了匹配度。

R会用协方差结构作为估计模型参数的起点。因此，根据我们选择的协方差结构将会得到不同的结果。如果你指定一个协方差结构太简单,那么你很有可能得到TypeI error (发现参数是显著但实际上并不是)，但如果你指定一个太复杂的模型也会得到一个TypeII error (发现参数不显著但事实上却相反)。我们将介绍四种在R中最为常见的协方差结构:

1. Variance components 基本方差结构（独立模型）

这个协方差结构非常简单，其假设所有随机效应都是独立的(即矩阵中所有的协方差都是0)。随机效应的方差被认为是相同的(即在矩阵中均为1)而求结果变量的方差和。这种协方差结构有时被称为独立模型。

2. Diagonal对角线结构

方差结构与基本方差结构相似，不同之处在于其方差被认为是异质的(即矩阵的对角线由不同的方差项组成)。这个结构再次假设方差是独立的，因此，所有的协方差都是0。

3. AR

这种类型代表一阶自回归结构。其说明方差之间的关系以一种系统化方式变化。如果你把矩阵的行和列想象成时间点，那么它假设重复测量的相关性在相邻的时间点是最高的。因此,在第一列,时间点1和2之间的相关性是ρ，假设值为0.3。当我们搬到时间点3,时间点1和3之间的相关性是ρ2值为0.09。它的值减少了，因为时间点1的分数与时间2的分数的相关性比其与时间3的分数的关系更大。时间点4时,相关性为ρ3值为0.027。因此,相邻时间点之间的相关性是假定为ρ,两个间隔点上的时间点1和2的相关性为ρ2,时间点1与间隔3的相关性则为ρ3。因此，随着时间点的推移，分数之间的相关性越来越小。方差被假定为同质的（即在矩阵中均为1），但是在这个协方差结构的另一个版本中可以是异质的。这种结构通常用于重复测量数据(特别是长时间测量，例如增长模型)。

4. Unstructured 非结构化

这种协方差结构是普遍的，协方差被认为是完全不可预测的，因为其不符合一个系统化的结构模式。

19.5. 一些实际问题

19.5.1. 理论假设

多级线性模型是回归的扩展，因此所有的回归假设都适用于多级模型(参见7.7.2节)。但是，需要注意的是，对于独立性和独立误差的假设有时可以通过一个多水平模型来解决，因为这个模型的目的是考虑由更高水平变量引起的情况之间的相关性。因此，如果一个有两个水平或3个水平的变量导致了独立性的缺乏，那么多层模型应该可以解决这个问题(尽管不总是这样)。在多层模型中还有两个与随机系数相关的附加假设。假设这些系数在整个模型中呈正态分布。因此，所有不同的截距和斜率都围绕这个模型呈正态分布。

19.5.2. 样本量和统计功效

对于测定固定效应或随机效应系数所需要的样本和统计量的问题比较复杂。从本质上说，关键信息是数据越多越好。随着模型中引入更多的水平，需要估计更多的参数，样本量也需要越大。Kreft和de Leeuw得出结论，如果您正在寻找跨水平间的交互作用，那么对于多水平变量的样本量至少需要20个区组，并且每个区组所包含的人数不能太少。

19.5.3. 集中变量 (Centring variables)

“Centring”指的是将一个变量转换成围绕一个定点数的偏差的过程。这个不动点可以是你选择的任何值，但通常用的是总平均值。例如Z分数的计算。

在多层模型中有两种典型的集中形式:总体平均集中 (grand mean centring)和组平均集中 (group mean centring)。总平均集中意味着对于给定的变量，我们用每个观测分数减去所有观测分数的平均值。组均值集中是指对于给定的变量，我们用每个观测分数减去组内所有观测分数的平均值。通常是只有一个水平的预测因子这两种集中形式都可以实现。而如果使用总平均集中通常是指只有一个水平的变量围绕只有两个水平变量。

在在多元回归中，截距表示所有预测因子取值为0时的结果值。这一数据处理过程的用处是当预测因子的值为0而没有意义时，如果预测因子以其平均值为中心，那么截距的值就是所预测因子的平均值。

然而在多层模型中，centring variables的作用稍显复杂。本质上，如果你用原始分数拟合一个多水平模型，然后再用以总体均值为中心的预测因子拟合同一个模型，那么得到的模型是等价的。因此，总体均值中心不会改变模型，但它会改变你对参数的解释(你不能把它们当作原始分数来解释)。当使用组均值时，要复杂得多。在这种情况下，原始分数模型在固定部分和随机部分都不等同于中心模型。一个例外是，只有截距是随机的，组均值作为二级变量重新引入模型(Kreft &de Leeuw, 1998)。

关于进行centring variables的决定是非常复杂的，需要在给定的分析中自己做出决定。Centring是克服预测变量间多重共线性的有效方法。当预测因子的零点没有意义时，它也很有帮助。最后，带有centring的预测因子的多级模型趋向于更稳定，并且来自这些模型的估计可以或多或少被视为相互独立，这可能是可取的。

这就产生了一个问题，究竟是总体均值集中还是组均值集中更好。做统计的人通常会把注意力集中在他们做事情的最好方法上，但是最好的方法往往取决于你真正想做的是什么。有些人根据一些统计标准来决定是使用群体集中还是大均值集中;然而，在不集中、组均值集中和总体均值集中之间没有统计上正确的选择。

19.6.多水平模型在R中的应用步骤（以书上的例子为例）

第一步：下载并导入需要的安装包。

install.packages("car"); install.packages("ggplot2"); install.packages("nlme"); install.packages("reshape")

然后加载这些安装包进R:

library(car);ibrary(ggplot2);library(nlme);library(reshape)

第二步：导入数据。

surgeryData=read.delim("Cosmetic Surgery.dat", header = TRUE)

第三步：画图

这一步可以利用ggplot包进行。这一步主要是为了观察数据的分布。以书上的例子为例，我们可以简单地观察两种条件下，也就是整容手术和正在等待手术的被试的基线生活质量和术后生活质量之间的关系。

PS：书上还演示了如果忽略数据的结构并进行ANOVA和ANCOVA分析，数据结果会如何。从而说明多水平模型其实是并不可怕，只是在以前的基础上进行的简单扩展。详见书本19.6.4和 19.6.5节。

第四步：评估进行对水平模型的条件

在诸如此类的多级分析中，第一步是首先评估是否需要这样做。首先，我们需要固定一个基线模型，其中我们只包括截距。第二，拟合一个模型，该模型允许截距在不同条件下变化。第三，我们对这两个模型进行比较，看看是否由于允许截距的变化而提高了模型的拟合度。如果是的话，便需要进行多水平模型分析。

在本例中，我们首先需要让R建立一个只包含截距的基线模型，使用gls()函数(广义最小二乘)来实现。

interceptOnly<gls(Post_QoL~1,data=surgeryData,method=“ML")summary(interceptOnly)

接下来，我们需要建立相同的模型，但这一次允许截距在不同条件下变化，在本例中，我们希望它们在不同的诊所中变化。我们使用lme(linear mixed effect)函数来实现这一点。

randomInterceptOnly<-lme(Post_QoL~1,data=surgeryData, random = ~1|Clinic, method = "ML") summary(randomInterceptOnly)

然后我们需要比较这两个模型。一种比较简单的方法是使用如下函数进行。

anova(interceptOnly,randomInterceptOnly

结果如下：

-2LL的变化服从卡方分布，因此得到的结果是χ2(1) = 107.65, p < .0001。我们可以得到截距在不同诊所下变化差异显著，需要进行多水平模型分析。

第五步：添加固定效应。

现在有了一个带有截距的基线模型，我们可以通过添加预测因子来开始构建最终的模型。让我们首先添加一个外科手术变量，定义为一个人是否接受了手术，或者是否在等待名单上。可以通过lme()函数实现。

randomInterceptSurgery<-lme(Post_QoL~ Surgery, data =surgeryData, random = ~1|Clinic, method = "ML") summary(randomInterceptSurgery)

结果发现手术似乎不是一个重要的预测因素。在假设中还包括了基线生活质量这一假设预测因素，所以我们也应该加上这个固定的影响。

randomInterceptSurgeryQoL<-lme(Post_QoL~ Surgery + Base_QoL, data =surgeryData,random=~1|Clinic,method="ML") summary(randomInterceptSurgeryQoL)

结果发现BIC和AIC均较前一模型有所下降，例如，BIC从1924.62降至1865.59。这意味着这个模型拟合的更好。同样，使用anova()函数比较模型。

anova(randomInterceptOnly,randomInterceptSurger,randomInterceptSurgeryQoL)

第六步：引入随机斜率

我们已经看到，包含随机截距对于这个模型非常重要，它显著地改变了log-likelihood。图2显示不同的诊所有不同的斜率。因此，我们现在可以看看增加一个随机斜率是否对建立模型有利。R语言如下：

addRandomSlope<lme(Post_QoL~Surgery+Base_QoL,data=surgeryData,random=~Surgery|Clinic,method="ML")

summary(addRandomSlope) anova(randomInterceptSurgeryQoL,addRandomSlope)

第七步：向模型中添加交互项

我们现在可以通过添加另一个变量来构建模型。我们测量的一个变量是这个人做整容手术的原因（Surgery:Reason）是为了解决身体问题还是纯粹为了虚荣？我们可以将这个变量添加到模型中，并观察它是否与手术相互作用来预测生活质量。我们的模型将简单地扩展以合并这些新项，每个项都有一个回归系数。本例中重新建立了两个模型，所用R语言如下：

addReason<lme(Post_QoL~Surgery+Base_QoL+Reason,data=surgeryData,random=~Surgery|Clinic, method = "ML")

finalModel<-lme(Post_QoL~Surgery+ Base_QoL +Reason+ Surgery:Reason, data =surgeryData, random = ~Surgery|Clinic, method = "ML")

同样使用anova()函数来比较不同模型之间的拟合度。

anova(addRandomSlope,addReason,finalModel)

如果想获得置信区间这些参数，可以使用intervals()函数来实现。其中，我们简单地指定我们想要的置信区间的模型，以及置信区间的级别作为比例(即， 0.99产生99%的置信区间)。

例如：

intervals(finalModel, 0.90)

intervals(finalModel, 0.95)

intervals(finalModel, 0.99)

最后，我们可以对这个finalModel进行解释。考虑到手术的原因是身体方面或外貌会对因变量的影响存在交互作用。为了打破这种相互作用，我们可以分别对两个原因组重新分析。可以通过以下两句来实现：

physical Subset<-surgeryData$Reason==1

cosmeticSubset<-surgeryData$Reason==0

接下来，我们再围绕这两项原因建立新的模型。如下：

physicalModel<-lme(Post_QoL~Surgery+Base_QoL, data =surgeryData, random = ~Surgery|Clinic,subset=physicalSubset, method = "ML")

cosmeticModel<-lme(Post_QoL~Surgery+Base_QoL,data=surgeryData,random=~Surgery|Clinic,subset=cosmeticSubset, method = "ML")

最后，我们获取这两个模型的总的数据信息。

summary(physicalModel)

summary(cosmeticModel)

19.7 增长模型

应用增长模型旨在发现某一变量在时间维度上的变化率。

19.7.1 增长模型（指数）

图19.10给出了可能存在的增长模型。下图展示了三条指数曲线，随着指数增大，曲线会越来越陡峭，时间维度上的变化会越来越快。

关于数据模型适配问题：1、你可以适配指数模型的指数直到比你所拥有的时间点少一个。2、指数模型被定义为一个简单有效的函数。

19.7.2 一个例子（略）

19.7.3 重构数据

应用多水平模型研究我们在时间维度上的测量数据面临的问题是需要一个与我们所用的不同的数据形式。对于重测实验设计我们通常将一个被试设置为一行：在这个例子中，重测时间变量被设置为四列（见3.9.4）。在第三章中我们称这种数据行为为宽数据形式。如果我们想要进行传统的重测方差分析，这种数据输出会很好；然而，对于多水平模型我们需要时间变量列展示，我们称之为长数据类型。因此我们需要重构数据。

19.7.4 设置基本模型

现在我们设置好了数据形式，我们可以分析数据了。基本上，我们按照之前的例子的方法进行分析。这里只有一点不同：因为我们在分析时间序列数据，我们必须对协方差结构进行建模。做这个最常用的办法是推算出一次回归协方差结构；要提醒你的是，这意味着与时间紧密联系的数据点被认为是比时间上数据点的距离更相关。在其他方面我们像先前例子一样设置模型。

首先我们适配仅包含截距的基线模型，像先前的例子一样，用gls() function:

intercept <-gls(Life_Satisfaction~ 1, data =restructuredData, method = "ML", na.action = na.exclude)

要记得我们要让R构造一个命名为intercept的变量，我们已经命名过Life-Satisfaction 是输出变量，它仅被截距项预测（函数中的~1）。函数剩下的部分命名数据（data=restructured-Data）和怎样估计模型（method=“ML”）。这是一个重要的和先前例子不同之处。我们已经包含了一个新选项 na.action=na.exclude。这是因为现在的数据文件（不像先前的例子）包含缺失数据，这些缺失值在数据文件中被命名为“NA”。这个选项告诉R当它遇到“NA”时该怎么做，我们在既往的例子中设置过（见R’s Souls‘Tip 19‘2）。没有这个选项模型会返回错误提示。

接下来我们需要适配同样的模型，但是这一次允许截距项varyacrosscontexts。像先前的例子一样，我们用lme() function:

randomIntercept<-lme(Life_Satisfaction~ 1, data =restructuredData, random = ~1|Person, method = "ML", na.action = na.exclude, control = list(opt="optim"))

这条命令的形式和先前例子相同。我们构造了一个模型（命名为randomIntercept）仅仅用截距预测生活满意度（Life-Satisfaction~1），但是也允许截距项vary across people（random=~1|Persom）。记住数据集中Person变量是一个用来预测数据数据是否来自同一个人的数值变量。我们再一次用最大似然估计（method=“ML”），包含缺失数据点（na.action=na.exclude）。然后要注意的是，这里有一个我们之前没有遇到过的新选项：control=list（opt=“option”），这个选项可以使得R更好的估计模型。通常情况下这个选项是设置好的默认值，但是对于用默认值无法计算的数据，我们必须要设置一下使得这些数据能被计算（见R’s Souls’Tip 19.4）。

19.7.5 添加时间作为固定影响

在一个增长曲线分析中，我们主要对“时间”这个固定的影响变量感兴趣。在我们的数据集中时间作为因变量，生活满意度分数在基线0、6个月、12个月和18个月时被记录下来。在先前的例子中，我们单独用lme()function构造了我们的模型；然而对于这个例子我们有很多选项(method = “ML”, na.action = na.exclude,control= list(opt=“optim”))，在我们每一个新模型中要分清这些选项。我们可以通过以下命令很快地把时间作为一预测变量加入到先前模型中：

timeRI<-update(randomIntercept, .~. + Time)

19.7.6 介绍自由截距项

我们可以很简单地用update()往模型中添加一个自由截距项重新命名模型中的自由部分。我们用update() function构造一个新模型（称为timeRS），这对于先前的模型（timeRI）非常重要，但是把模型的自由部分更新成了random = ~Time|Person:

timeRS<-update(timeRI,random= ~Time|Person)

19.7.7构造协方差结构

现在我们有了一个基本自由影响模型，我们接下来接下来介绍协方差结构或误差的建模主题。我们通过用correlation=选项，x是指先前定义的几个协方差结构之一。我们可以用update()function往先前模型中添加一个协方差结构来构造一个新模型（命名为ARModel），这个模型和先前的模型很相像，但是加入了一个first-order autoregressive covariance structure：

ARModel<update(timeRS,correlation=corAR1(0, form = ~Time|Person))

19.7.8 模型比较

截止到现在我们已经构造了五个模型，我们通过以下方法检查模型与数据的匹配度，我们可以通过以下命令比较五个模型：

我们通过以下命令快速检查最后的模型和进行置信区间的参数估计：Summary(ARModel); intervals(ARModel)

19.7.9 加入高指数项

我们可以用update()往模型中添加二次项来构造新模型（timeQuardratic）：

timeQuadratic<-update(ARModel, ~.+ I(Time^2))

同样的我们可以往模型中添加三次项，然后通过方差分析与一次项模型、二次项模型相比较并获得最后的结果：

我们的模型已经将时间作为预测项，我们需要消除这个预测项，用多项式进行替换。这个操作可以通过以下命令进行：

polyModel<-update(ARModel, .~ poly(Time,3))

19.8 如何报告多水平模型

对于最后的模型报告你可以像这样展示：