随机效应模型介绍及实例分析

一、模型定义

1.1引入

在给出模型的具体定义之前先看看下面这个案例
例：为研究家庭背景对学生成绩的影响，考虑了以下三种情形，：
1、假设学生来自A学校，为研究家庭背景对学生成绩的影响，构建了线性模型y=a+bx+u，其中y为学生成绩，a为A校学生成绩的平均水平（截距项），x为学生家庭背景（自变量），u衡量具体学生的成绩差异（随机误差）
2、假设学生来自A，B两校，构建模型y=a+cs+bx+u,s是一个二分变量，c是s的回归系数。
3、假设学生来自不同学校，学校是从全省范围内随机选取的，因为不同学校存在差别，反映在学生的平均成绩是不同的 ，这时模型需要修改为y=a+uj+bx+u（随机效应模型）,a表示学校总体的平均成绩，uj是随机变量，表示某校学生的平均成绩与全体学生的平均成绩的差值（好学校为正值，差学校为负值）
通过比较可以发现情形一、学生从A校抽取，学校是固定的，只需要考虑学生之间的差异
情形二、两个固定的学校，即考虑学生之间的差异，也要考虑学校之间的差异
情形三、学生是从全省任意选择的学校里抽取的，学校是任意的，这时我们也要考虑不同学校之间的差异，随机效应模型多了一个效应uj，uj是一个随机变量

一般来说， 究竟把一个效应看作是随机的，还是固定的，这取决于研究的目的和样品取得的方法。随机效应模型仅包括作为固定效应的截距和一组定义的随机效应

1.2模型一般形式

随机效应模型的一般形式可以写成:

$\begin{cases}\\Y=X{\beta}+Zu+{\epsilon}\\ \\E(u)=0, Var(u)={\Sigma}\\ \\{\epsilon}N(0,{\sigma}^2I_n),且与u无关\\ \end{cases}$⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧​Y=Xβ+Zu+ϵE(u)=0,Var(u)=ΣϵN(0,σ2In​),且与u无关​
其中Y是nx1维观测向量，X为固定参数β（固定效应，）的nxp设计矩阵，Z为随机参数u（随机效应）的nxq设计矩阵，ε为n维随机误差
则：
$\\E(Y)=X{\beta} \\Var(Y)=E[(Y-EY)'(Y-EY)] \\=E(Y'Y)-E(Y)'E(Y) \\=z{\Sigma}Z'+{\sigma^2I_n}=V({\theta})$E(Y)=XβVar(Y)=E[(Y−EY)′(Y−EY)]=E(Y′Y)−E(Y)′E(Y)=zΣZ′+σ2In​=V(θ)

$\\其中，θ表示Z{\Sigma}中的未知参数，它反映了随机效应u及随即及误差{\epsilon}中方差及协方差参数\\$其中，θ表示ZΣ中的未知参数，它反映了随机效应u及随即及误差ϵ中方差及协方差参数

二、模型的参数估计

2.1固定效应和随机效应的估计

2.2参数的极大似然估计

**Newton-Raphason算法

2.3参数的限制极大似然估计

三、实例分析

基于线性混合效应模型分析性别和态度对音调高低的影响
数据来自于Winter and Grawunder (2012)，该数据共分为5列（subject：主体；gender：性别；scenario：情景；attitude：态度；frequency：音调）

3.1描述性统计

根据箱线图，可以看出男性的声音比女性的声音低，但群体之间又有个体差异。为此，假设不同的随机截距项来对个体差异建模。态度不友好时音调要高于态度友好时，但男性在不同的态度时音调有更多的重叠部分。

参数估计

R语言有很多软件包可以做混合线性模型，比如lme4、nlme、lmerTest等，本文采用的lme4

参数说明：DV-因变量
fixed_factor - 固定因子（即考察的自变量）
random_factor - 随机因子
模型建立

结果

通过比较随机效应的标准差(Std.Dev)，可以发现音调在区域间的变化大于情景。 固定效应截距为256.846Hz，表明女性在态度不友好的时候平均音调为256.846 Hz； 斜率-19.721 Hz表明态度从不友好到友好音调降低19.721 Hz； 斜率-108.517 Hz表明女性比男性音调高108.516Hz.

检验结果表明，各变量都是显著的，另模型的拟合优度检验结果为R2=0.858.表明模型对数据的解释效果较好。

四、附录

4.1数据说明

研究数据来自于Winter and Grawunder (2012)，该数据共分为5列

subject：个体；
gender：性别；
scenario：情景；
attitude：态度；
frequency：音调

数据来源（http://www.bodowinter.com/uploads/1/2/9/3/129362560/politeness_data.csv）

4.2代码

#install.packages("lme4 )
#输入数据并简单观察数据结构
politeness=read.csv("E:/StudyLife/politeness_data.csv")
head(politeness);tail(politeness)
summary(politeness)
str(politeness)
colnames(politeness)
#删除缺失数据
politeness=politeness[-which(!complete.cases(politeness))]
#数据探索
par(mfrow = c(2,2))
boxplot(frequency~subject,politeness,main = '个体')
boxplot(frequency~scenario,politeness,main = '情景')
#boxplot(frequency~attitude,politeness,main = '态度')
#boxplot(frequency~gender,politeness,main = '性别')
boxplot(frequency ~ attitude*gender,col=c("white","lightgray"),politeness)

模型拟合

library(lme4)
politeness.model = lmer(frequency ~ attitude +gender + (1|subject) +(1|scenario),data=politeness)
summary(politeness.model)

**固定效应与随机效应的估计

library(lme4)
library(lmerTest)
politeness.model = lmer(frequency ~ attitude +gender + (1|subject) +(1|scenario),data=politeness)
anova(politeness.model)
ranova(politeness.model)

模型拟合优度检验

library(lme4)
library(sjstats)
#p_value(politeness.model)
politeness.model = lmer(frequency ~ attitude +gender + (1|subject) +(1|scenario),data=politeness)
r2(politeness.model)