---
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中有相對的變異
* 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) 判斷是否要用隨機效果模型
* 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)
```
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較好
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
Sign in with Wallet
Connect another wallet