Transcript 分层线性回归模型

预备知识1：线性模型

二元回归模型
矩阵形式
最小二乘估计（ordinary least
squares,OLS）

估计量

估计量方差
其中
总平方和
fi 为预测值
 回归平方和
 残差平方和
 判定系数（coefficient of determination）R
squared


调整R squared
广义最小二乘（generalized leastsquares,GLS）

如果
阵

估计量

方差
这里
为已知协方差矩
预备知识2：固定效应模型

可加效应模型
 y ijk    i   j   ijk

i  1,  , a , j  1,  , b, k  1,  , m


2
  ijk i.i.d N 0, 

i  0,   j  0

 i
j




方差分析（ analysis of variance,
ANOVA）
假设
H01 : 1  2    a  0
H02 : 1  2    b  0
偏差平方和的分解




SS T   y ij  y..
i
j
2




  y i.  y..   y. j  y..   y ij  y i.  y. j  y..
i
j
2
i
j
2
i
j
2
 SS A  SS B  SS e
检验统计量
y i.    i   i.
y. j     j  . j

2
 ESS A  E  i   i.  .. 
 i j

2
2
 b i  a  1

i

y..    ..
同理ESS B  a    b  1
2
j
2
j



2
ESS T  E  i   j   ij  .. 
 i j

2
2
2
 b i a   j  ab  1
i
j
 ESS e  a  1b  1 2
EMSA
 SS A 
 E

 a 1 
b i2
i
a 1
a   2j
 2
 SS B 
j
2
EMSB  E




b 1
 b 1 


SS e
2


EMSe  E


 a  1b  1 
MSA
MSB
 F1 
F2 
MSe
MSe
交互效应模型
 y ijk     i   j  ij   ijk
i  1,  , a , j  1,  , b, k  1,  , m


 ijk i.i.d N 0,  2


 i  0,   j  0


i
j

ij  0, j  1,  , b


i

ij  0, i  1,  a



j



方差分析
假设
H 01 : 1  2    a  0
H 02 : 1  2    b  0
H 03 : ij  0 对一切i  1,, a j  1,, b
偏差平方和的分解

SS T   y ijk  y...
i
j

k



  y i..  y...   y. j.  y...
i
j
k

2

i
j
k


2
  y ijk  y ij.   y ij.  y i..  y. j.  y...
i
j
k
2
i
 SS A  SS B  SS e  SS AB
j
k

2
检验统计量
y i..    i   i..
y. j.     j  . j.
y ij.    i   j  ij   ij.

2
ESS A  E  i   i..  ... 
 i j k

 bm  i2  a  1 2

i
ESS B  am  2j  b  1 2
j

y...    ...



2
ESS e  E   ijk   ij.   abm  1 2
 i j k




2
ESS T  E  i   j  ij   ijk  ... 
 i j k

 bm    am   m   abm  1 2
2
i
2
ij
2
j
i
j
i
j
ESS AB  m   a  1b  1 2
2
ij
i
j
 SS A 
EMSA  E

 a 1 
 SS B 
EMSB  E

 b 1 
bm  i2
i
a 1
am  2j
j
b 1
 2
 2
m  
2
ij


SS AB
j
 
EMSAB  E
 2
 a  1b  1  a  1b  1
 SS e

2


EMSe  E


 abm  1 
MSA
 FA 
MSe
MSB
FB 
MSe
FAB
MSAB

MSe
随机效应模型
 y ijk    i   j  ij   ijk i  1, , a , j  1, , b, k  1, , m

2

i
.
i
.
d
N
0
,

ijk


i i.i.d N 0,  2

2

i
.
i
.
d
N
0
,

j



ij i.i.d N 0,  2

 诸 ijk、诸i、诸 j、诸ij 相互独立








方差分析
假设
H 01 :   0
2

H 02 :   0
2

H 03 :  2  0
偏差平方和的分解

SS T   y ijk  y...
i
j

k

2




  y i..  y...   y. j.  y...   y ijk  y ij.
i
j
k

2
i
j
  y ij.  y i..  y. j.  y...
i
j
k
 SS A  SS B  SS e  SS AB

k
2
2
i
j
k

2
检验统计量
y i..    i 
 j
j
b

 ij
j
b
  i..
   i    i.   i..
y...        ..  ...



2
 ESS A  E  i    i.  ..   i..  ... 
 i j k

 bm a  1 2  ma  1 2  a  1 2
同理 ESS B  amb  12  mb  1 2  b  1 2



2
2
ESS e  E   ijk   ij.   abm  1
 i j k


2
ESS T  E  i     j    ij  ..   ijk  ... 
 i j k

 bm a  1 2  amb  12  mab  1 2  abm  1 2


 ESS AB  a  1b  1  ma  1b  1
2

 SS A 
2
2
2
EMS

E



m


bm



A



a

1



 EMS  E SS B    2  m 2  am 2
B



 b 1 

 SS AB

EMSAB  E
   2  m 2

 a  1b  1 

 SS e 
 EMSe  E
   2

 abm  1 
2

MSA
 F1 
MSAB
MSB
F2 
MSAB
MSAB
F3 
MSe
预备知识3：三大检验
似然比检验LR
Wald检验
拉格朗日乘子检验LM
三大检验的引入
（1）模型是非线性的
 （2）约束是非线性的
 （3）扰动项分布是非正态的，


在这些情况下，F检验不再适用，通常需要
采用LR、Wald、LM其中之一来检验约束条
件是否成立。
三大检验方法共同点




这三个检验方法都是渐进等价的，他们所
用统计量的小样本分布是未知的，但大样
本下都渐进服从自由度为约束个数的卡方
分布。
三大检验方法是三种基于极大似然法的大
样本检验方法。
根据模型的特点采用不同的检验方法。
模型视为给定参数的数据生成过程的集合。
极大似然估计（ML）
（一）极大似然原理
假设对于给定样本 Y , X ， 其联合概率分布存
在
f Y , X ; 
。将该联合概率密度函数视为未知参数 
的函数，则 f Y , X ;  称为似然函数（Likelihood Function）, 即观
测到所给样本的可能性.
极大似然原理就是寻找未知参数  的估计
大，或者说寻找使得样本
Y , X 
ˆ
，使得似然函数达到最
出现的概率最大的 ˆ 。
求极大似然函数估计值的一般步骤：
（1）
（2）
（3）
（4）
写出似然函数；
对似然函数取对数，并整理；
求导数 ；
解似然方程
极大似然估计，是一种概率论在统计学的应用，
它是参数估计的方法之一。说的是已知某个随
机样本满足某种概率分布，但是其中具体的参
数不清楚，参数估计就是通过若干次试验，观
察其结果，利用结果推出参数的大概值。极大
似然估计是建立在这样的思想上：已知某个参
数能使这个样本出现的概率最大，我们当然不
会再去选择其他小概率的样本，所以干脆就把
这个参数作为估计的真实值。
极大似然估计量（MLE）的性质
（1）一致性：
ˆ
是


的一致估计量，即

lim P ˆ     )  1,  为任意给定的正数。
n
（2）渐进有效性：
ˆML
是渐进有效的且达到所有一致估计量
的Cramer-Rao下界，即是所有一致渐进正态估计量中方差最小的
（3）渐进正态性
似然比检验（LR）
检验思想：如果参数约束是有效的，那么加上
这样的约束不应该引起似然函数最大值的大幅
度降低。也就是说似然比检验的实质是在比较
有约束条件下的似然函数最大值与无约束条件
下似然函数最大值。似然比定义为有约束条件
下的似然函数最大值与无约束条件下似然函数
最大值之比。以似然比为基础可以构造一个服
从卡方分布统计量
似然比检验（LR）
1、似然比
H0 : g     C
命题：
如果约束是无效的，有约束的最大似然函数值当然不会
超过无约束的最大似然函数值，但如果约束条件“有
效”，有约束的最大值应当“接近”无约束的最大值，
这正是似然比检验的基本思路。
2
L(  ,  )

似然比：
L( ˆ , ˆ 2 )
2
ˆ
ˆ )
L
(

,

无约束模型似然函数值：
2
有约束模型似然函数值： L(  ,  )
显然 0    1。如果原假设是真，则  趋近于1；
如果  太小，则约束无效，拒绝原假设。
可以证明，对大样本来说，检验统计量为，
LR  2ln   2 ln L( ˆ , ˆ 2 )  ln L(  ,  2 )  ~  2 (q)
拒绝域， LR  
2
1
(q )
似然比检验另一种表达，
LR  2 ln   n(ln e*' e*  ln ee) ~  2 (q)
'
e e 有约束模型残差平方和；
**
ee无约束模型残差平方和；
Wald检验
检验思想：如果约束是有效的，那么在没有约束情
况下估计出来的估计量应该渐进地满足约束条件，
因为MLE是一致的。
以无约束估计量为基础可以构造一个Wald统计量，
这个统计量也服从卡方分布
H0 : g     C
Wald检验
如果约束条件为真，则 g   MLE   C  0 不应该显著异于
零，其中  MLE 是无约束极大似然估计值。当 g   MLE   C
显著异于零时，约束条件无效，拒绝原假设。检验统计
量。 Wald检验实际基于g（ β ）和C之间的距离。
 


1
 
a
W  ( g ˆ  C) Var g (ˆ )  C  ( g ˆ  C) ~  2 (q)


Wald只需要估计无约束模型，但需要计算渐进协方差矩
阵。
在线性约束条件下， Wald检验
H0 : R  r
1
a
W  ( Rˆ  r )  Rˆ ( X X ) R ( Rˆ  r ) ~  2 (q)
2
拒绝域，
1
W  2 (q)
Wald统计量另一种表达形式，
n(e*' e*  ee)
W 
~  2 (q)
ee
'
e e 有约束模型残差平方和；
* *
ee无约束模型残差平方和；
拉格朗日乘子检验（LM）
检验思想：在约束条件下，可以用拉格朗日方法构
造目标函数。如果约束有效，则最大化拉格朗日函
数所得估计量应位于最大化无约束所得参数估计值
附近。
这里也是构造一个LM统计量该统计量服从卡方分布。
拉格朗日乘子检验（LM）
拉格朗日乘子检验（LM），又称为Score检验。该检验基
于约束模型，无需估计无约束模型。
假设约束条件为 H0 : g    C ，在约束条件下最大化对
数似然函数，另 
格朗日函数为
表示拉格朗日乘子向量，此时，拉
LnL ( )  LnL( )   g( )  C 
约束条件下最大化问题就是求解下式根，
LnL ( ) LnL( )

 g   0


LnL ( )
 g(  ) 
 g ( )  C  0
其中，g 是矩阵g= 
的转置


  

如果约束成立，对数似然函数值不会有显著变化。这就

意味着在一阶条件下，第二项应该很小，特别是
应该很小。
因此，约束条件是否成立检验转化成检验 H0： =0 ，这就
是拉格朗日乘子检验的思想。
但是直接检验 H0： =0
比较困难，有一个等价而简单
的方法。如果约束条件成立，在约束估计值处计算对数似然
函数的导数应该近似为零，如果该值显著异于零，则约束条
件不成立，拒绝原假设。
对数似然函数的导数就是得分向量，因此，LM检验就是检验
约束条件下参数估计值的得分向量值是否显著异于零，因而，
LM检验又称为得分检验。
在最大似然估计过程中，通过解似然方程 S (ˆ)  0
，
可以求出无约束估计量 ˆ ；如果计算有约束估计量

在此处得分，则 S ( ) 一般不为零，但是如果约束有
效，则 S ( ) 趋近于零。
在原假设成立条件下，
a
LM  S ( ) I ( ) S ( ) ~  (q)
1
2
对于线性约束
将有关量代入上式得，
ne*' X ( X ' X ) 1 X ' e*
2
2
LM=
=
nR
~

(q)
'
e*e*
'
e e 有约束模型残差平方和；
**
R 2是e*对X 回归的拟合优度；
拒绝域，
LM  nR   (q)
2
2
LM统计量另一种表达形式，
n(e e  ee)
W
~  2 (q)
ee
'
* *
'
* *
'
e e 有约束模型残差平方和；
**
ee无约束模型残差平方和；
LR、 Wald 、LM关系(一般情况下成立)：
Wald  LR  LM
对于似然比检验，既需要估计有约束的模型，也需要估计
无约束的模型；对于Wald检验，只需要估计无约束模型；
对于LM检验，只需要估计有约束的模型。一般情况下，由
于估计有约束模型相对更复杂，所有Wald检验最为常用。
对于小样本而言，似然比检验的渐进性最好，LM检验也较
好，Wald检验有时会拒绝原假设，其小样本性质不尽如人
意。
多层线性 模型
hierarchical linear model （HLM）
概念
分层线性模型（hierarchical linear model HLM
）又名多层线性模型 （Multilevel Linear Model
MLM）、层次线性模型（Hierarch Linear
Mode1）、多层分析（Multilevel Analysis/Model
）。
 HLM又被通俗的称为“回归的回归”。
 一般线性回归和多重线性回归都是发生在单一
层面，HLM相对于更适用于嵌套数据（nest
data）。”

假设




由于个体行为不仅受个体自身特征的影响，也受到
其所处环境（群体/层次）的影响。
相对于不同层次的数据，传统的线性模型在进行变
异分解时，对群组效应分离不出，而增大模型的误
差项。
而且不同群体的变异来源也可能分布不同，可能满
足不了传统回归的方差齐性假设。在模型应用方面，
不同群体（层次）的数据，也不能应用同一模型。
鉴于传统方法的局限性，分层技术则解决了这些生
态谬误（Ecological Fallacy）。
两个层面的假设：
个体层面：这个与普通的回归分析相同，
只考虑自变量X对因变量Y的影响。
 群组层面：群组因素W分别对个体层面中
回归系数和截距的影响。

数学模型：

个体层面：

群组层面：

涉及到多个群组层次的时候原理与之类似，可以把较低级层次的群组，
如不同的乡镇层面与不同的县市层面，可以这样理解，乡镇即是一个
个体，群组即是不同的县市。
更多层次的可以这样理解，一直是下一层对上一层回归系数和截距的
回归。
与普通的“回归的回归”不同的是，整个计算过程通过迭代过程完成。



合并模型


固定（非随机项）
为固定参数

随机项

随机参数
和
。
模型的GLS表达方式

令

则
固定j不独立

假设

这里

则
独立。

GLS 迭代估计固定系数 。
不同j 相互
矩阵形式
这里
X内部变量（within），W外部变量
（between），WX交互作用(cross-level )

这里
独立
并相互
大矩阵形式
这里

方差记为
记
 GLS 估计

应用




多层线性模型的适用范围非常广，凡是具有嵌套和分层
的数据均可使用多层线性模型进行分析。
此外，多层线性模型还可以用于纵向研究。采用多层分
析的方法处理重复测量数据与时间变量之间的关系。
在多层结构中可以对非平衡测量数据得到参数的有效估
计。因此用多层分析法处理重复测量的数据，不要求所
有的观测个体有相同的观测次数。
在纵向调查研究中，由于各种各样的原因，被试个体观
测值部分缺失的情况时有发生，因此多层分析法处理缺
失数据而不影响参数估计精度的这一特征，使得多层分
析法处理在处理纵向观测数据时，比传统多元重复测量
方法有很大的优势。
处理多元重复测量数据，多层分析
法优点：

多层分析法通过考虑测量水平和个体水平
不同的差异，明确表示出个体在水平1（不
同测量点）的变化情况，因而对于数据的
解释（个体随时间的增长趋势）是在个体
与重复测量交互作用基础上的解释，即不
仅包含了不同测量点的差异，而且包含了
个体之间存在的差异。

多层分析法对数据资料较传统多元重复测
量方法有较低的要求，对于重复测量的次
数和重复测量之间的时间跨度都没有严格
的限制。不同个体可以有不同的测量次数，
测量与测量之间的时间跨度也可以不同。

多层分析模型可以定义重复观测变量之间复杂
的协方差结构，并且对所定义的不同的协方差
结构进行显著性检验。在多层分析模型中，通
过定义第一水平和第二水平的随机变异来解释
个体随时间的复杂变化情况，当数据满足传统
多变量重复测量模型对数据的要求和假设时，
层次分析法得到与传统固定效应多元重复测量
模型相同的参数估计和假设检验结果。用多层
分析模型可以考虑更高一层的变量，如不同地
区儿童对个体增长的影响。
缺点
用于多层分析模型的参数估计方法较传统
估计参数的方法要复杂得多。
 不能处理变量之间间接的影响关系和处理
复杂的观测变量和潜变量之间的关系。

中学数据（High school, hg）
主要变量：
 1 数学成绩 （math achievement, mathach）
 2社会经济状况(socio-economic status, ses)
 1977个学生， 40个学校：21个公立学校
（public school），19个天主教学校
(catholic school)
 2个sector: public 和 catholic

1线性回归：fit<-lm(mathach~ses, hs)
2固定效应模型
summary(lm( mathach ~ Sector/ses-1, hs))
 主效应Sector，交互效应Sector:ses
 summary(lm( mathach ~ ses * Sector-1, hs))
 主效应Sector,ses，交互效应Sector:ses
 lm1<-lm(matach~factor(school)/ses-1,hs)
 主效应school,交互效应 school:ses
 Summary(lm(matach~factor(school)*ses-1,hs))
 主效应school,ses,交互效应 school:ses

Wald test












ddu <- up( hs, ~ factor(school))
dim(ddu)
ind <- ddu$Sector == "Catholic"
L <- rbind( "Catholic" = ind,"Public" = 1-ind)
L <- L/apply(L,1, sum)
L <- cbind( rbind( L, 0,0), rbind( 0,0,L))
rownames( L ) <- c("Cath Int", "Pub Int", "Cath Slope",
"Pub Slope")
L%*%lml$coefficients
wald (lml, L)
diffmat <- rbind( "Int" = c( -1, 1, 0, 0), Slope = c( 0 , 0, -1, 1))
diffmat %*% L%*%lml$coefficients
wald (lml, diffmat %*% L)#different of sector
3多元方差分析（MANOVA）
library(nlme)
 lcoefs <- coef( lmList( mathach ~ ses |
factor(school), hs))
 lm.mult <- lm( as.matrix(lcoefs) ~Sector,
up( hs, ~ factor(school)))
 summary(lm.mult)

4学校之间模型
fit.eco <- lm( mathach ~ ses,up( hs, ~
factor(school), all = T))
 summary( fit.eco )

5分层线性模型
library(nlme)
 fit <- lme( mathach ~ ses * Sector, hs, random
= ~ 1 + ses | school)
 summary(fit)
 wald(fit, 'ses')
# overall test for 'ses'
 wald(fit, 'Sector')

6加入背景变量（ contextual
variables ）的HLM








fitc <- lme( mathach ~ ses * Sector + cvar(ses,school),
hs, random = ~ 1 + ses | school)
summary( fitc )
wald( fitc )
wald( fitc, -1) # overall test of FE model
wald( fitc, 'ses') # overall test of all 'ses' effects
wald( fitc, 'Sector') # overall test of all 'Sector' effects
wald( fitc, ':') # overall test of all interaction effects
cvar and cvars are intended to create contextual
variables in model formulas. If 'x' is numerical, cvar is
equivalent to capply(x,id,mean) and cvars is equivalent
to capply(x,id,sum).
Wald检验
L <- list( 'Effect of ses' = rbind(
"Within-school" = c( 0,1,0,0,0),
"Contextual" = c( 0,0,0,1,0), "Compositional" =
c( 0,1,0,1,0)))
 wald ( fitc , L )
 dvar is equivalent to x - cvar(x,by) and creates what is
commonly known as a version of 'x' that is 'centered
within groups' (CWG). It creates the correct matrix
for a factor so that the between group interpretation
of the effect of cvar(x,by) is that of the 'between group'
or 'compositional' effect of the factor.

Using CWG instead of raw effect








fitcd <- update( fitc, . ~ dvar(ses,school)*Sector +
cvar(ses,school))
fitcd <- lme( mathach ~ dvar(ses,school)*Sector +
cvar(ses,school), hs, random = ~ 1 + ses | school）
summary( fitcd )
summary( fitc )
fitca <- update( fitc, random = ~ 1 + dvar( ses, school )
| school)
summary( fitca )
summary( fitc )
anova( fitca, fitc )
Diagnostics with Level 1 residuals
qqnorm( fitc )
 qqnorm( fitc , abline = c(0,1), id = .01 )
 plot( fitc ) # not as generous as in lm, need to
make your own: # This is a plot of residuals (zscore) versus fitted value
 diag.df <- data.frame( resid = resid(fitc), fitted
= fitted(fitc))
 summary( lm( resid ~ fitted, diag.df))

Diagnostics with Level 2 residuals













hs$ses.m <- with( hs, cvar( ses, school))
fitc <- lme( mathach ~ ses * Sector + ses.m, hs, random = ~ 1 + ses | school )
some( coef ( fitc ) ) # BLUP in each cluster
some( ranef ( fitc ) ) # RE in each cluster # Note: coef( fitc ) = fixed.effect +
random.effect
re <- ranef( fitc, aug = T) # creates data frame with up( dd, ~ school) variables
some( re )
plot( re ) # special plot method
qqnorm( fitc, ~ ranef(.) | Sector, id = .05)
library(p3d)
Init3d( family='serif', cex = 1.2)
Plot3d( `(Intercept)` ~ ses + ses.m | Sector, re) # note funny names in backquotes
Ell3d()
Id3d( pad=1 )