
---
title: "整合分析（Meta-Analysis）"
author: "許清芳，鄭中平"
date: 2024-04-27
output: slidy_presentation
---

## 整合分析是什麼？

整合分析是分析多個研究結果，可以用來綜合研究結果，瞭解研究異質性（研究結果變異），或是找到造成變異的可能原因。由於分析單位不是人，而是分析，是「分析的分析」，所以稱為「meta」-analysis。

## 分析人與分析研究

### 分析人，我們假設要分析同學的身高：
- 收集同學的身高，換成相同單位
- 算平均以及標準差
- 看看性別是否是造成身高有變異的原因之一
- 擔心一下有否同學缺席，影響結論

### 分析研究，我們分析研究的效果大小：
- 收集研究的效果量，換成相同的效果量指標
- 算平均以及標準差
- 看看是否有某些研究特性，造成研究效果有變異
- 擔心一下有沒有研究缺席，影響結論

## 整合分析的材料是什麼？可以回答什麼問題？

### 整合分析的材料
- 整合分析的材料是一些研究（通常是某個特定議題）
- 紀錄每個研究的效果量，與/或可能影響效果的因素（描述研究而非個人的變項）

### 整合分析可以回答什麼問題
- 某特定議題研究的平均效果，對分歧的研究下綜合結論
- 某特定議題研究的變異，瞭解研究結果可以變化多大
- 哪些研究特性（例如，實驗操控方法、測量方法、研究對象等等）可能影響研究結果，解釋過去的分歧，文獻間的不一致，導引未來研究


## 看一個例子

### 恐懼管理理論
恐懼管理理論認為人類會因為意識到自己終將死亡而心生恐懼，透過認同國家、民族、文化價值或世界觀等，我們得以超愈自己終將死亡的恐懼。恐懼管理理論的研究常以實驗法進行，操弄死亡顯著性高與低的情境，看看參與者是否表現出較多的文化價值認同。

### 資料來源
Burke、Martens與Faucher（2010）以164篇文章中的277個研究，探索實驗的死亡顯著性效果大小，並考慮參與者特性、抽樣區域等研究屬性的可能調節效果，底下顯示分析抽樣區域的可能調節效果。

## 看一下資料
```{r pkgs, echo=TRUE, cache=FALSE, message=FALSE, warning=FALSE}
#整體設定，含載入套件
source("https://raw.githubusercontent.com/ChungPingCheng/R4BS2/main/R4BS_setup.R")
```


```{r }
#讀進資料
dta <- read.table("https://raw.githubusercontent.com/ChungPingCheng/R4BS2/main/Data/TMTmeta.txt", header=T, stringsAsFactors = TRUE)
```


```{r }
#顯示前六筆資料
#程式報表11.1
dim(dta)
head(dta)
```

## 描述性統計

```{r eval=FALSE}
#依不同研究群，看看效果量等的描述性統計
#程式報表11.2
dta |>
  select(-Study) |>
  gtsummary::tbl_summary(by=Group,
                statistic = list(all_continuous() ~ "{mean} ({sd})"))
```

```{r out.width = "35%",echo = FALSE}
dta |>
  select(-Study) |>
  gtsummary::tbl_summary(by=Group,
                statistic = list(all_continuous() ~ "{mean} ({sd})"))|>
            as_gt() |>
            gt::gtsave(file = "tbl11_2.png")

knitr::include_graphics("tbl11_2.png")
```

```{r eval=FALSE}
#依不同區域，看看效果量等的描述性統計
#程式報表11.3
dta |>
  select(-Study) |>
  gtsummary::tbl_summary(by=Region,
                  statistic = list(all_continuous() ~ "{mean} ({sd})"))
```

```{r out.width = "50%",echo = FALSE}
dta |>
  select(-Study) |>
  gtsummary::tbl_summary(by=Region,
                  statistic = list(all_continuous() ~ "{mean} ({sd})")) |>
            as_gt() |>
            gt::gtsave(file = "tbl11_3.png")

knitr::include_graphics("tbl11_3.png")
```

## 資料圖形化

```{r}
#效果量直方圖
#圖11.2
ggplot(data=dta, 
       aes(ES, after_stat(density)))+
  geom_histogram(binwidth = function(x) 2 * IQR(x) / (length(x)^(1/3)),
                 col=1, fill=8)+
  labs(x="研究所得效果（相關係數量）",
       y="密度")
```

## 整合分析 

```{r }
#利用相關與樣本數計算相關的抽樣變異，並放進原先資料
dta <- metafor::escalc(measure="COR", ri=ES, ni=N, data=dta)
head(dta)
```

```{r }
#計算所有研究的平均效果
#程式報表11.4
summary(res <- metafor::rma(yi=yi, vi=vi, data=dta))
```

- 平均效果是 r = .373（ z = 29.70, p < .001），顯著不為0，信賴區間則是.348到.397。綜合這277個研究結果，死亡顯著性操弄有效果。
-  $\tau^2$ 是效果量變異的估計；$I^2$是這些變異與總變異（研究間變異加上抽樣造成變異）間的比值，介於0-1間；$H^2$是總變異除以抽樣造成變異
-  結果顯示，$\tau$ = .184，對於介於正負1間的相關來說不算太小，$I^2$與$H^2$則顯示研究間變異占總變異的 90.46\%、總變異是抽樣變異的10.48倍（亦即研究間變異是抽樣變異的9.48倍），也顯示有不小的研究間變異。Q檢定拒絕研究間同質的虛無假設（Q(df= 276) = 12088.78, p < .001），支持以上觀察。
- 以上顯示研究間存在異質性，使用預設模型（隨機效果模型）來估計與檢定應屬適當。


### 森林圖，圖示結果
```{r}
#圖11.3
metaviz::viz_forest(res)
```
- 圖中的每一個線段，對應一個研究，其中黑點是效果量，線段則是95\%信賴區間
- 圖中虛線標示效果量為0，如果信賴區間壓到這條虛線，代表效果可能為0，不顯著

## 檢查出版偏誤

### 漏斗圖 
```{r }
#檢查出版偏誤，先繪製 funnelplot
#圖11.4
metaviz::viz_funnel(res)
```

- 橫軸是效果量（愈向右愈大），縱軸是效果量標準誤（即vi開根號，愈下方愈大）- 黑色實線三角稱做信賴輪廓（confidence contours），是以資料的平均效果量推估之效果量與標準誤分布範圍；如果不存在出版偏誤，點（觀察效果量）應該散布在黑色三角
- 藍色區域的三角稱做顯著輪廓（significance contours），藍色區域裡邊的黑線，是5\%顯著範圍，在此三角形之外是顯著結果（p<.05），藍色區域外邊的黑線，是1\%顯著範圍，在此三角形之外則是顯著結果（p<.01）。


```{r }
#程式報表11.5前半
ranktest(res)
```

```{r }
#程式報表11.5後半
regtest(res)
```

- 若沒有出版偏誤，效果量與標準誤（或變異數）應該無關
- rank correlation test 計算效果量與抽樣變異間相關
- 迴歸方式檢驗效果量與標準誤（或變異數）之斜率

## 整合調節分析

```{r }
#後續要分析區域與研究群的調節效果
#先排除區域不明與第四組研究群體資料
dta <- dta |>
  dplyr::filter(Region != 99 & Group < 4) |> 
  dplyr::mutate( Region = factor(Region),
                 Group = factor(Group))
```

```{r }
#抽取277筆中的13%，依區域與 ES 排序，繪製森林圖，標示區域
set.seed(201408)
dta_13 <- dta |> 
  dplyr::sample_frac(.13) |> 
  arrange(Region, ES)
```


```{r }
#圖11.5
with(dta_13, metafor::forest(yi, vi, slab = Region))
```


```{r }
#檢查區域效果。排除區域不明與第四組研究群體資料
#程式報表11.6
dta <- 
  dta |> arrange(Region, ES)
res_rgn <- rma(yi=yi, vi=vi, mods = ~ Region, data=dta)
res_rgn |> summary()
```

```{r }
#繪製 forest，標示區域（不只抽出一部份）
#圖11.6
p1 <- metaviz::viz_forest(res_rgn, group=dta[,'Region'],
              xlab = "r", x_limit=c(-.80,1.1), type="study_only")
p2 <- metaviz::viz_forest(res_rgn, group=dta[,'Region'],
              xlab = "r", x_limit=c(-.80,1.1),type="summary_only")
p1/p2+plot_layout(heights = c(7, 1))
```

```{r }
#檢查出版偏誤，繪製 funnelplot，設定點的形狀
#圖11.7
metaviz::viz_funnel(res_rgn, group=dta[,"Region"],group_legend=T)+
    scale_shape_manual(values = c(15:17,0,1)) 
```

```{r }
#每組抽取10筆，依 Group 與 ES 排序，繪製森林圖，以 Group 標示
set.seed(2014)
dta1 <- dta[dta$Group==1,][sample(145,10),]
dta2 <- dta[dta$Group==2,][sample(36,10),]
dta3 <- dta[dta$Group==3,][sample(93,10),]
dta30 <- rbind(dta1,dta2,dta3)

dta30 <- dta30 |> 
  arrange(Group, ES)

#圖11.8
with(dta30, metafor::forest(yi, vi, slab = Group))
```

## 整合分析前的準備

- 選定一個主題，收集相關研究
- 把各研究中的效果量，與/或研究特性摘出來
- 整理資料，效果量要同一個指標，研究特性可能要分組

### 整合分析前的準備(細部版)

- 論文中最好說明收集策略（哪個資料庫？那幾個關鍵字？）
- 如果有排除，要說明排除標準
- 常常會變成一張流程圖，說明排除標準以及留下來的研究數量
- 摘出研究特性並適當分組其實很累，常常來來回回
- R 或網路上有一些程式可以協助效果量指標轉換


## 小結
- 相較於一般的文獻回顧，整合分析是量化地總結過去的文獻，有時有一錘定音之效
- 一般希望納入的分析不要太少，但也看過八篇文章做整合分析的，這與選擇主題是否有很多研究有關。如果要做整合調節分析，納入的分析就要更多了。
- 有些人認為整合分析不用收資料，做起來比較容易，但可能需要看比較多的文獻
- 綜合效果的平均數可能意義不大，找到研究異質性的原因可能貢獻較大。

## 習題

### 去個人化與反社會行為

當人們感到個人責任時，通常會遵守社會規範。然而，在匿名環境中，這些規範經常被忽略。許多實驗通過改變匿名群體的大小或特性，來觀察違規行為的變化。

### 資料
Postmes 與 Spears(1998) 以整合分析的方式瞭解去個人化與反社會行為的關係，記錄了七十個研究的資料，包括四個變項，依序是研究編號、效果量（以r作為效果量）、樣本數、匿名團體的成員屬性（1、2、3分別標示學生、小孩與其他）與參與者是否在群體中（1、2 分別標示在與不在）。

### 小題

(1)請用下列方式讀入資料，並繪製森林圖。
dta <- read.csv("https://raw.githubusercontent.com/ChungPingCheng/R4BS2/main/Data/Deindividuation.csv", header=T) 

(2)請計算平均效果量，並檢驗是否為0。

(3)請繪製漏斗圖，並檢驗是否有出版偏誤。

(4)請檢驗匿名團體的參與者是否在群體中，是否會影響研究中效果的大小。
