---
# System prepended metadata

title: HLM Model (Mixed Effect Model)
tags: [Fixed Effect Model, Econometrics, 計量經濟, HLM, Random Effect Model, Mixed Effect]

---

---
GA: UA-159972578-2
---

###### tags: `計量經濟` `Econometrics` `HLM` `Mixed Effect` `Random Effect Model` `Fixed Effect Model` 

# HLM Model (Mixed Effect Model)

### 前情提要
* Panel：排除了group effect後，自變數和目標變數有什麼影響
* Mixed Effect：想知道group如何影響自變數和目標變數的關係

### 在回歸中使用變數解釋分群變異
* 分群變數(群間)：班級
* 群組變數(群內)：班級的士氣
* 一般變數
* 係數效果反映在：
    * 截距：分群變數（Random的分佈，有一個平均值和變異數）
    * 斜率
    * 係數
    * 效果

## 1. HLM (Hierarchical Linear Model) 

### (1) 資料結構
又稱LMM(Linear Mixed Model)及(Multilevel Model)
* 當你的資料符合下面兩個狀況時，則選擇用Multilevel Analysis而不使用OLS Regression
    1. Dependence as a nuisance：殘差獨立性假設未過
    2. Dependence as an interesting phenomenon：在不同level中有相對的變異

![](https://i.imgur.com/GvhrkNX.png =80%x)

* Level 1變數（micro level）：巢形於群組下的變數e.g.學生
* Level 2變數（macro level）：group level，位於較高層次的變數e.g.學校
* Y (應變數)是在Level 1


### (2) HLM 在做什麼？
1. 計算組間變異(between variance)、組內變異(within variance)
2. 固定效果(fixed effect)、隨機效果(random effect)
3. 估計兩個參數：截距(intercept) & 斜率(slope)
4. 計算ICC & rwg 兩個重要指標


### (3) HLM指標

<b> ◉ 可靠度(Reliability) </b>
* 當我們想用一個平均值當族群變數的代表時，要去確認reliability
* 意義：測量的結果是robust的，大家差異不大
* 指標：ICC2


<b> ◉ 組內變異：rwg (reliability within group) </b>
* rwg = 1- Observed Group Variance/Expected Random Variance
* 利用一個維度項目估計「組內一致性」
* 用rwg（reliability within group）去確認每個組內是否有共識（對某一自變數）
    * 若>0.7代表可以用組平均去代表
    * <0.7代表他們之間有太大的變異，不能直接使用
    * 0.7是通則，但不是一定的，可以用rwg.sim去計算該值
* ranvar的A為該自變數的scale（例如：1-5分，就是(5^2-1)/12=2）


```{r}
RWG.HRS<-rwg(bhr2000$HRS,bhr2000$GRP,ranvar=10.00) #ranvar=(11^2-1)/12=10
mean(RWG.HRS[,2]) # 0.7353417

# example 10 groups & 5 response options ,repeat 10000 times
RWG.OUT<-rwg.sim(gsize=10, nresp=5, nrep=10000)
summary(RWG.OUT) #0.95 upper value=0.53

# 看在不同信心水準底下的rwg值
quantile(RWG.OUT,c(.90,.95,.99))
```

<b>◉ 組間變異：ICC1, ICC2 (intraclass correlation coefficient)</b>

* rwg是對x做；ICC是對y做
* ICC1(即R2)：群體類別可以解釋目標變數的變異的比率，y會隨著x而變
* ICC2：群組層次的信度，y在群組內夠不夠穩定（reliablility）、有可預測性，用來檢視群體平均數的信度值
* 用ANOVA去判斷要用哪一個模型
* ICC大的話提供做Mixed Effect的理由

```
hrs.mod<-aov(HRS~as.factor(GRP),data=bhr2000)
summary(hrs.mod) # 查看intercept及residual
```
```
##                  Df Sum Sq Mean Sq F value Pr(>F)    
## as.factor(GRP)   98   3371   34.40    12.5 <2e-16 ***
## Residuals      5301  14591    2.75                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

```
ICC1(hrs.mod) # 0.1741008
ICC2(hrs.mod) # 0.9199889
```



## 2. HLM for Multilevel Data

### (1) 資料是否需要用到HLM模型
* 群組間平均值(截距)有無顯著差異
* 群組間自變數(斜率)有無顯著差異


### (2) 判斷是否要用隨機效果模型

![](https://i.imgur.com/49x67k5.png =30%x)

* Yij（i: individual, j: group）
* 𝛾00：固定效果的截距項
* u0j：隨著group做變動的variance
* rij：隨著個體、群體變異
```{r}
N1<-lme(WBEING~1,random=~1|GRP,data=bh1996,control=list(opt="optim"))
VarCorr(N1) # 查看intercept及residual的variance
```
```
## GRP = pdLogChol(1) 
##             Variance   StdDev   
## (Intercept) 0.03580079 0.1892110
## Residual    0.78949727 0.8885366
```
<b>◉ 固定效果模型</b>
+ 目標變數：WBEING(well being)
+ ~1：代表只放截距項

<b>◉ 隨機效果模型</b>
* 該截距項帶有一個random的效果
* random by |後面的東西，即GRP
```{r}
N0<-gls(WBEING~1,data=bh1996,control=list(opt="optim")) # without Random Effect
```
可以用likelyhood公式(給定模型a,b的observe data的機率)看兩個模型之間有無顯著差異
或是用ANOVA直接幫你算：
```{r}
anova(N0, N1)
```
```
##    Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## N0     1  2 19540.17 19553.98 -9768.084                        
## N1     2  3 19353.34 19374.06 -9673.669 1 vs 2 188.8303  <.0001
```

* 若不顯著：代表兩者差不多，就可以選簡單的模型
* 若顯著：則要使用比較複雜的Random模型（AIC較小的）


## 3. 改善模型
<b>◉ 固定效果模型</b>
```{r}
M1<-lme(WBEING~HRS,random=~1|GRP,data=bh1996,method="ML")
summary(M1)
M2<-lme(WBEING~HRS+G.HRS,random=~1|GRP,data=bh1996,method="ML")
summary(M2)  # HRS & G.HRS are signifacant
anova(M1,M2) # p < 0.0001
```
```
##    Model df      AIC      BIC    logLik   Test L.Ratio p-value
## M1     1  4 19234.89 19262.52 -9613.445                       
## M2     2  5 19200.92 19235.46 -9595.461 1 vs 2 35.9672  <.0001
```
使用M2較好
```{r}
v = sapply(list(N1,M1,M2),VarCorr)[1:2,]
colnames(v) = c("N1","M1","M2"); rownames(v) = c("u0j","rij"); v
```
```
##     N1 (NULL)    M1 (HRS)     M2(G.HRS)         
## u0j "0.03580079" "0.02426467" "0.0130006"  # 組間差異
## rij "0.78949727" "0.77996426" "0.7800093"  # 組內差異
```
* intercept：組間變異
* residual：組內變異
 
上值為變異，看變異是否有下降（被變數解釋）
HRS和G.HRS都幫整個模型（N1）解釋掉了不少組間變異，但組內變異比較難
```{r}
M3<-lme(WBEING~HRS+LEAD+G.HRS,random=~LEAD|GRP,data=bh1996,method="ML")
M3a<-update(M3,random=~1|GRP)
```
![](https://i.imgur.com/4hMNXbl.png =80%x)

N1~M2在組間變異下降幅度較大；M3a則是在組內變異下降較多
加上Leader的影響，M3a好很多
<br>

<b>◉ 支線（隨機效果模型增加一個LEAD）</b>

M3在隨機效果模型上加上領導風格
```{r}
M4<-lme(WBEING~HRS+LEAD+G.HRS+LEAD:G.HRS,random=~LEAD|GRP,data=bh1996,method="ML")
```
M4再加上交互作用
```{r}
anova(M3,M4)
```
```
##    Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## M3     1  8 17810.35 17865.61 -8897.177                        
## M4     2  9 17809.61 17871.77 -8895.806 1 vs 2 2.741624  0.0978
```
M3較好
![](https://i.imgur.com/Alnkxpt.png =55%x)
u2j不是我們關心的，隨著隨機效果變數增加，就會多一個uij；但我們想看的是組間/內差異
從數字去比較，看到M3和M4在解釋組內變異效果上相較之前M2有顯著的提升（與M3a有一樣的效果）
```
anova(M3a,M3)
```
```
##     Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## M3a     1  6 17834.15 17875.59 -8911.076                        
## M3      2  8 17810.35 17865.61 -8897.177 1 vs 2 27.79872  <.0001

```
一次可以比較一個變化，雖然不能直接比較M2跟M3，但從跟M3a的比較來看M3較好


## 4. 選擇模型
* 指標：變異解釋量
* Variance Explained = 1 – (Var with Predictor/Var without Predictor)
```{r}
explained = function (x) 100*(1- as.numeric(x[2])/as.numeric(x[1]))
apply(v, 1, explained)       # %explained N1 to M1
apply(v[,2:3], 1, explained) # %explained M1 to M2
apply(v[,3:4], 1, explained) # %explained M2 to M3a
apply(v[,c(1,4)], 1, explained) # %explained N1 to M3a
```
N1 to M1
```
##       u0j       rij 
## 32.223088  1.207479
```
M1 to M2
```
##          u0j          rij 
## 46.421690466 -0.005774624
```
M2 to M3a
```
##      u0j      rij 
## -5.03746 17.06781
```
N1 to M3a
```
##      u0j      rij 
## 61.85699 18.06447
```
將method都統一的話，就可以一次比較N1~M4了
```{r}
N1<-lme(WBEING~1,random=~1|GRP,data=bh1996,method = "ML")
anova(N1,M1,M2,M3a,M3,M4)
```

```
##     Model df      AIC      BIC    logLik   Test   L.Ratio p-value
## N1      1  3 19347.56 19368.28 -9670.778                         
## M1      2  4 19234.89 19262.52 -9613.445 1 vs 2  114.6660  <.0001
## M2      3  5 19200.92 19235.46 -9595.461 2 vs 3   35.9672  <.0001
## M3a     4  6 17834.15 17875.59 -8911.076 3 vs 4 1368.7704  <.0001
## M3      5  8 17810.35 17865.61 -8897.177 4 vs 5   27.7987  <.0001
## M4      6  9 17809.61 17871.77 -8895.806 5 vs 6    2.7416  0.0978
```
可以看到M3a較M2好；M3又比M3a好；M4效果下降了，所以應該使用M3
