1 用交叉隨機因子模擬數據

本研究邀請受試者對臉孔表情進行分類，以反應時間作為主要因變量。實驗假設臉孔有兩種類型：來自subject的ingroup或outgroup，每張臉在刺激集中只出現一次。研究目的在於探討不同類型人臉的分類速度是否存在差異。

隨機因素1：受試者subjects（T/tau）
隨機因子2：人臉（O/omega）
固定因素1：類別（ingroup、outgroup）
- within subject：受試者將同時看到ingroup和outgroup的臉孔
- between face：ingroup或outgroup

2 需要的軟體

#loadrequiredpackages
library("lme4") #model specification / estimation
library("lmerTest") #provides p-values in the output
library("tidyverse") #data wrangling and visualisation
library("kableExtra") #organize the table

# ensure this script returns the same results on each run
set.seed(8675309)

3 建立數據生成參數

本模擬研究中，100名受試者的每一個將對所有50個刺激項目（25個ingroup和25個outgroup）作出反應，總共5000個觀察值。

3.1 數據結構

參數β₀和β₁是固定效應，我們將遇到典型刺激的典型受試者的平均RT設置為800毫秒，並假設outgroup臉孔的反應通常比ingroup臉孔慢50毫秒。

item_id：從 1 到 50 運行並索引項目編號。
category：臉孔是ingroup還是outgroup(項目編號1-25是ingroup；項目編號26-50是outgroup)。
RT：是參與者對實驗的反應時間。
基本模型： \[RT_{si}=β_0+β_1X_i+e_{si}\]
β₀：受試者的總平均RT。
β₁：ingroup和outgroup的平均RT差異。
X_i：ingroup(-0.5)；outgroup(0.5)。
e_si：受試者s與項目i的殘差。

# set all data-generating parameters
beta_0  <- 800 # intercept; i.e., the grand mean 總平均(800/ms)
beta_1  <-  50 # slope; i.e, effect of category 平均差(50/ms)

3.2 隨機效應參數

我們假設每個T_0s來自平均值為零且標準差未知的正態分佈τ0(tau_0)，平均值為零反映每個隨機效應都是偏離總平均的假設。對於逐項隨機截距，我們假設O_0i值是從平均值為零且標準差未知的正態分佈中取樣的，ω0(omega_0)。在模擬中，我們將每個受試者s的隨機截距設置為100，每個項目i的隨機截距設置為80。

模型： \[RT_{si}=β_0+T_{0s}+O_{0s}+β_1X_i+e_{si} \]
tau_0：受試者隨機截距的標準差。
omega_0：逐項隨機截距的標準差。
T_0s：每個受試者s的隨機截距。
O_0i：每個項目i的隨機截距。

omega_0 <-  80 # by-item random intercept sd
tau_0   <- 100 # by-subject random intercept sd

然而，模型與類別固定效應β1相關不足。目前，我們的模型假設每個受試者對ingroup臉孔的情緒分類比outgroup臉孔快80毫秒。但實際情況一些受試者會對ingroup/outgroup的差異更敏感，因此增加每個受試者s的隨機斜率T_1s。

模型： \[RT_{si}=β_0+T_{0s}+O_{0s}+(β_1+T_{1s})X_i+e_{si} \]
T_1s：每個受試者s的隨機斜率。
tau_1：受試者隨機斜率的標準差。
rho：受試者隨機截距和斜率之間的相關。
sigma：殘差的標準差。

tau_1   <-  40 # by-subject random slope sd
rho     <-  .2 # correlation between intercept and slope
sigma   <- 200 # residual (error) sd

最終我們所建立的模型為： \[RT_{si}=β_0+T_{0s}+O_{0s}+(β_1+T_{1s})X_i+e_{si} \] 接下來，我們會使用此數據結構來模擬抽樣。

4 模擬抽樣

Srep 1.設定受試者與刺激項目數量

本次將模擬來自100名受試者的數據對25張ingroup臉孔和25張outgroup臉孔的反應，因此將n_subj設為100，並設定n_ingroup和n_outgroup每個條件下的刺激項目數為25。

# set number of subjects and items
n_subj     <- 100 # number of subjects
n_ingroup  <-  25 # number of items in ingroup
n_outgroup <-  25 # number of items in outgroup

4.1 模擬刺激項目的抽樣

Srep 2.創建列表，列出每個項目i屬於哪個類別，以及它的隨機效應O_0i。

# simulate a sample of items
n_items <- n_ingroup + n_outgroup

items <- data.frame(
  item_id = seq_len(n_items),
  category = rep(c("ingroup", "outgroup"), c(n_ingroup, n_outgroup)),
  X_i = rep(c(-0.5, 0.5), c(n_ingroup, n_outgroup)),
  O_0i = rnorm(n = n_items, mean = 0, sd = omega_0)
)

另外，我們會使用數字預測器來表示每個刺激項目i出現在哪個類別中（對應於模型中的Xi）。由於我們預測對ingroup臉孔的反應將比outgroup臉孔更快，因此將ingroup設為-0.5，outgroup為+0.5；之後我們會把效果編碼因子乘以類別的固定效果(beta_1=50)來模擬數據，其中ingroup臉孔平均與總平均相差-25毫秒，而outgroup臉孔與總平均相差25毫秒。

# effect-code category
items$X_i <- recode(items$category, "ingroup" = -0.5, "outgroup" = +0.5)

4.2 模擬受試者的抽樣

Srep 3.創建列表，列出每個受試者及其兩個相關的隨機效應。

方法1-MASS::mvrnorm() from the MASS package.

接著我們將模擬對單個受試者的抽樣，生成一個表格，列出每個受試者及兩個相關的隨機效應。由於不能使用rnorm()從單量分佈中獨立於T_1s值對T_0s值進行採樣，而必須從雙變量正態分佈中抽取T_0s、T_1si對（每個受試者一對）。因此我們會使用mvrnorm()函數，輸入描述此分佈的三個參數——兩個variances和一個correlation來使用matrix()函數轉換為2x2方差-協方差矩陣，然後使用Sigma參數將該矩陣輸入給mvrnorm()。

# calculate random intercept / random slope covariance
covar <- rho * tau_0 * tau_1

# put values into variance-covariance matrix
cov_mx  <- matrix(
  c(tau_0^2, covar,
    covar,   tau_1^2),
  nrow = 2, byrow = TRUE)

# generate the by-subject random effects
subject_rfx <- MASS::mvrnorm(n = n_subj,
                             mu = c(T_0s = 0, T_1s = 0),
                             Sigma = cov_mx)

# combine with subject IDs
subjects <- data.frame(subj_id = seq_len(n_subj),
                       subject_rfx)

方法2：rnorm_multi() from the faux package.

另一種方法是使用faux（DeBruine，2020）中的函數rnorm_multi()，它通過指定每個變量的平均值(mu)、標準差(sd)與相關(r)，從多元正態分佈生成一個包含n個模擬值的表，它可以是單個值(應用於所有配對)、相關矩陣或相關矩陣右上三角形中的值的向量。

# sample from a multivariate random distribution 
subjects <- faux::rnorm_multi(
    n = n_subj, 
    mu = 0, # means for random effects are always 0
    sd = c(tau_0, tau_1), # set SDs
    r = rho, # set correlation, see ?rnorm_multi
    varnames = c("T_0s", "T_1s")
  )

# add subject IDs
subjects$subj_id <- seq_len(n_subj)

確認數值

data.frame(
  parameter = c("omega_0", "tau_0", "tau_1", "rho"),
  value = c(omega_0, tau_0, tau_1, rho),
  simulated = c(
    sd(items$O_0i),
    sd(subjects$T_0s),
    sd(subjects$T_1s), 
    cor(subjects$T_0s, subjects$T_1s)
  )
)

##   parameter value  simulated
## 1   omega_0  80.0 68.7759151
## 2     tau_0 100.0 91.1610878
## 3     tau_1  40.0 41.0218392
## 4       rho   0.2  0.1083145

4.3 模擬實驗

由於所有受試者對所有項目都有反應，我們可以通過使用tidyverse的函數crossing()製作一個包含受試者和項目表中所有可能組合的表格來設定模擬實驗。每個實驗都有與之相關的隨機誤差，反映了由於未知因素導致的波動，我們通過從平均值為0且SD為sigma的正態分佈中取樣進行模擬。

# cross subject and item IDs; add an error term
# nrow(.) is the number of rows in the table
trials <- crossing(subjects, items)  %>%
  mutate(e_si = rnorm(nrow(.), mean = 0, sd = sigma))

4.4 計算反應值

有了這個實驗結果表，再加上beta_0和beta_1，我們就有了根據上面定義的線性模型計算RT所需的全部數值： \[RT_{si}=β_0+T_{0s}+O_{0s}+(β_1+T_{1s})X_i+e_{si} \] 因此，我們通過模型計算RT，之後使用dplyr::select()來保留我們需要的列。

dat_sim <- trials %>%
  mutate(RT = beta_0 + T_0s + O_0i + (beta_1 + T_1s) * X_i + e_si) %>%
  select(subj_id, item_id, category, X_i, RT)

4.4.1 繪圖

下面對模擬數據的數值進行繪圖以方便觀察模擬結果。

ggplot(dat_sim, aes(category, RT, color = category)) +
  # predicted means
  geom_hline(yintercept = (beta_0 - 0.5*beta_1), color = "orange") +
  geom_hline(yintercept = (beta_0 + 0.5*beta_1), color = "dodgerblue") +
  # actual data
  geom_violin(alpha = 0, show.legend = FALSE) +
  stat_summary(fun = mean,geom="crossbar", show.legend = FALSE) +
  scale_color_manual(values = c("orange", "dodgerblue")) +
  ggtitle("Predicted versus simulated values")

#長條圖
par(mfrow=c(2, 1))
ingroup <- dat_sim$RT[which(dat_sim$category == "ingroup")]
hist(ingroup)
outgroup <- dat_sim$RT[which(dat_sim$category == "outgroup")]
hist(outgroup)

4.5 數據模擬功能

為了方便嘗試不同參數或生成數據集以進行Power分析，我們可以將上面的所有程式碼放入自定義函數中，默認值為本研究所使用的數值，下方的程式碼為先前幾段程式碼的濃縮版。設定完畢後，就可以使用函數my_sim_data()來產生具有默認參數的數據集，或使用my_sim_data(n_subj=500,beta_1=0)生成具有500個受試者且沒有類別影響的數據集。

# set up the custom data simulation function
my_sim_data <- function(
  n_subj      = 100,   # number of subjects
  n_ingroup  =  25,   # number of ingroup stimuli
  n_outgroup =  25,   # number of outgroup stimuli
  beta_0     = 800,   # grand mean
  beta_1     =  50,   # effect of category
  omega_0    =  80,   # by-item random intercept sd
  tau_0      = 100,   # by-subject random intercept sd
  tau_1      =  40,   # by-subject random slope sd
  rho        = 0.2,   # correlation between intercept and slope
  sigma      = 200) { # residual (standard deviation)
  
  # simulate a sample of items
  items <- data.frame(
    item_id = seq_len(n_ingroup + n_outgroup),
    category = rep(c("ingroup", "outgroup"), c(n_ingroup, n_outgroup)),
    X_i = rep(c(-0.5, 0.5), c(n_ingroup, n_outgroup)),
    O_0i = rnorm(n = n_ingroup + n_outgroup, mean = 0, sd = omega_0)
  )

  # simulate a sample of subjects
  subjects <- faux::rnorm_multi(
    n = n_subj, mu = 0, sd = c(tau_0, tau_1), r = rho, 
    varnames = c("T_0s", "T_1s")
  )
  subjects$subj_id <- 1:n_subj
  
  # cross subject and item IDs 
  crossing(subjects, items)  %>%
    mutate(
      e_si = rnorm(nrow(.), mean = 0, sd = sigma),
      RT = beta_0 + T_0s + O_0i + (beta_1 + T_1s) * X_i + e_si
    ) %>%
    select(subj_id, item_id, category, X_i, RT)
}

5 分析模擬數據

5.1 設定公式

接下來，將對模擬數據進行分析。 lmer() 的第一個參數是定義線性模型結構的模型公式，公式將計算RT。 - 公式： \[RT \sim 1+X\_i+(1 \mid item\_id+(1+X\_i \mid subj\_id)\]

RT：反應時間。
1：大截距(beta_0)。
X_i：是項目i的ingroup/outgroup操作的預測變量。
(1|item_id)：指定一個by-subject隨機截取(O_0i)，允許截距(1)在項目(item_id)的隨機因子上變化，用於估計O_0i基礎參數，即omega_0。
(1+X_i|subj_id)：指定特定受試者的隨機截距(T_0s)加上類別的特定受試者隨機斜率(T_1s)，允許類別的截距和效果（由X_i編碼）隨受試者的隨機因素(subj_id)而變化，為估計T_0s和T_1s值下的三個參數的指令，即tau_0、tau_1和rho。

5.2 解讀lmer摘要

lme4 函數的其他參數是在其中找到值的數據框的名稱(dat_sim)，之後載入lmerTest，將p-values使用Satterthwaite近似得出，lmer()中的默認估計方法「限制似然估計(REML=TRUE)」是最合適的(Luke，2017)，並使用summary()函數查看結果。

第一段為模型擬合一般資訊，包括convergence、最小值與最大值、Q1與Q3以及中位數等；第二段為按受試者隨機效應的隨機效應參數的估計值，Groups的值為subj_id，tau_0和tau_1的估計值分別為93.81和51.41，受試者隨機截距和斜率之間的估計相關性為-0.05，逐項隨機效應的參數估計值，omega_0的估計值為63.98，殘差估計值為203.13；第三段為固定效應的參數估計值，即β₀和β₁，估計值約為818.03和44.64，並提供標準誤、估計的自由度(Satterthwaite方法)、t-values以及p-values。

※補充：原論文模擬之隨機效應參數的估計值，tau_0和tau_1的估計值分別為91.74和57.43，受試者隨機截距和斜率之間的估計相關性為0.12，omega_0的估計值為63.81，殘差估計值為203.18；固定效應的參數估計值，β₀和β₁估計值約為807.72 a和39.47

# fit a linear mixed-effects model to data
mod_sim <- lmer(RT ~ 1 + X_i + (1 | item_id) + (1 + X_i | subj_id),
                data = dat_sim)

summary(mod_sim, corr = FALSE)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: RT ~ 1 + X_i + (1 | item_id) + (1 + X_i | subj_id)
##    Data: dat_sim
## 
## REML criterion at convergence: 67732.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8300 -0.6755  0.0014  0.6799  3.6306 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  subj_id  (Intercept)  8801     93.81        
##           X_i          2643     51.41   -0.05
##  item_id  (Intercept)  4094     63.98        
##  Residual             41260    203.13        
## Number of obs: 5000, groups:  subj_id, 100; item_id, 50
## 
## Fixed effects:
##             Estimate Std. Error     df t value Pr(>|t|)    
## (Intercept)   818.03      13.35 120.73  61.290   <2e-16 ***
## X_i            44.64      19.67  54.57   2.269   0.0272 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

另外，也可以使用broom.mixed::tidy()輸出經過整理的固定跟隨機效果表格。

其中，「effect」為表示是固定效應參數(fixed)還是隨機效應參數(ran_pars)；「group」按隨機效應參數的隨機因子列名稱(或殘差)指定組；「term」是指固定效應的預測項，也是隨機效應的參數。同時，我們在broom.mixed::tidy() 的標準輸出加入了sim，broom.mixed::tidy(mod_sim) ，因此可以將先前設定的模擬參數與該數據集中的估計參數進行比較；最後也列出了標準誤、t-values、估計的自由度和固定效應的p-values。

# get a tidy table of results
broom.mixed::tidy(mod_sim) %>% 
  mutate_if(is.numeric, round, 3) %>%
  mutate(
    parameter = c("beta_0", "beta_1", "omega_0", "tau_0", "tau_1", "rho", "sigma"),
    value = c(beta_0, beta_1, omega_0, tau_0, tau_1, rho, sigma),
  ) %>%
  select(1:3, 9, 10, 4:8) %>%
  knitr::kable()%>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = F)

effect	group	term	parameter	value	estimate	std.error	statistic	df	p.value
fixed	NA	(Intercept)	beta_0	800.0	818.030	13.347	61.290	120.731	0.000
fixed	NA	X_i	beta_1	50.0	44.639	19.671	2.269	54.573	0.027
ran_pars	subj_id	sd__(Intercept)	omega_0	80.0	93.813	NA	NA	NA	NA
ran_pars	subj_id	cor__(Intercept).X_i	tau_0	100.0	-0.048	NA	NA	NA	NA
ran_pars	subj_id	sd__X_i	tau_1	40.0	51.408	NA	NA	NA	NA
ran_pars	item_id	sd__(Intercept)	rho	0.2	63.983	NA	NA	NA	NA
ran_pars	Residual	sd__Observation	sigma	200.0	203.125	NA	NA	NA	NA

5.3 參數設定

透過前面的介紹，我們可以知道如何進行參數設定並產生模擬數據，並估計beta_0、beta_1、tau_0、tau_1、rho、omega_0和sigma的值進行分析，因此我們可以按照研究需求更改參數設定。

6 計算Power

首先，我們創建一個函數single_run()，使用my_sim_data()生成數據，對其進行分析，並回傳一個包含參數估計值的表格，以此對應一次的Power模擬。

# set up the power function
single_run <- function(filename = NULL, ...) {
  # ... is a shortcut that forwards any additional arguments to my_sim_data()
  dat_sim <- my_sim_data(...)
  mod_sim <- lmer(RT ~ X_i + (1 | item_id) + (1 + X_i | subj_id),
                dat_sim)
  
  sim_results <- broom.mixed::tidy(mod_sim)
  
  # append the results to a file if filename is set
  if (!is.null(filename)) {
    append <- file.exists(filename) # append if the file exists
    write_csv(sim_results, filename, append = append)
  }
  
  # return the tidy table
  sim_results
}

更改參數並產生模擬示例：將每組中的項目數量增加到50並將類別的影響減少到20毫秒。

# run one model with default parameters
single_run()

## # A tibble: 7 × 8
##   effect   group    term                 estim…¹ std.e…² stati…³    df   p.value
##   <chr>    <chr>    <chr>                  <dbl>   <dbl>   <dbl> <dbl>     <dbl>
## 1 fixed    <NA>     (Intercept)          812.       14.0   57.8   92.1  3.63e-74
## 2 fixed    <NA>     X_i                   47.4      23.4    2.02  50.7  4.85e- 2
## 3 ran_pars subj_id  sd__(Intercept)       80.0      NA     NA     NA   NA       
## 4 ran_pars subj_id  cor__(Intercept).X_i   0.162    NA     NA     NA   NA       
## 5 ran_pars subj_id  sd__X_i               39.9      NA     NA     NA   NA       
## 6 ran_pars item_id  sd__(Intercept)       79.1      NA     NA     NA   NA       
## 7 ran_pars Residual sd__Observation      201.       NA     NA     NA   NA       
## # … with abbreviated variable names ¹estimate, ²std.error, ³statistic

# run one model with new parameters
single_run(n_ingroup = 50, n_outgroup = 45, beta_1 = 20)

## # A tibble: 7 × 8
##   effect   group    term                estim…¹ std.e…² stati…³    df    p.value
##   <chr>    <chr>    <chr>                 <dbl>   <dbl>   <dbl> <dbl>      <dbl>
## 1 fixed    <NA>     (Intercept)         7.92e+2    13.0   61.0   178.  2.29e-121
## 2 fixed    <NA>     X_i                 1.88e+1    17.1    1.10  101.  2.75e-  1
## 3 ran_pars subj_id  sd__(Intercept)     9.93e+1    NA     NA      NA  NA        
## 4 ran_pars subj_id  cor__(Intercept).X… 7.88e-2    NA     NA      NA  NA        
## 5 ran_pars subj_id  sd__X_i             3.55e+1    NA     NA      NA  NA        
## 6 ran_pars item_id  sd__(Intercept)     7.89e+1    NA     NA      NA  NA        
## 7 ran_pars Residual sd__Observation     1.99e+2    NA     NA      NA  NA        
## # … with abbreviated variable names ¹estimate, ²std.error, ³statistic

另外，可以使用purrr::map_df函數重複運行模擬並將結果保存到數據表中，模擬完成後，我們可以回傳數據查看固定效應的估計值。

filename <- "sims/sims.csv" # change for new analyses
if (!file.exists(filename)) {
  # run simulations and save to a file
  reps <- 100
  sims <- purrr::map_df(1:reps, ~single_run(filename))
}

# read saved simulation data
sims <- read_csv(filename, col_types = cols(
# makes sure plots display in this order
group = col_factor(ordered = TRUE),
term = col_factor(ordered = TRUE)
))
sims %>%
filter(effect == "fixed") %>%
select(term, estimate, p.value)

## # A tibble: 200 × 3
##    term        estimate  p.value
##    <fct>          <dbl>    <dbl>
##  1 (Intercept)    793.  5.25e-86
##  2 X_i             71.0 1.74e- 3
##  3 (Intercept)    820.  3.03e-76
##  4 X_i             45.9 6.48e- 2
##  5 (Intercept)    793.  2.68e-86
##  6 X_i             28.2 1.95e- 1
##  7 (Intercept)    811.  2.83e-83
##  8 X_i             66.3 5.08e- 3
##  9 (Intercept)    811.  1.55e-72
## 10 X_i             56.7 2.97e- 2
## # … with 190 more rows

表中的每一行都是固定效應參數的估計值和來自單次Monte-Carlo模擬計算的相關p-values((Intercept)=β₀和X_i= β₁)。我們需要計算顯著運行的比例。首先，我們計算p-value < alpha，其中alpha是誤報率(如.05)，在效果顯著為TRUE(T=1)，在效果不顯著為FALSE(F=0)，並計算比例。

# calculate mean estimates and power for specified alpha
alpha <- 0.05

sims %>% 
  filter(effect == "fixed") %>%
  group_by(term) %>%
  summarise(
    mean_estimate = mean(estimate),
    mean_se = mean(std.error),
    power = mean(p.value < alpha),
    .groups = "drop"
  )

## # A tibble: 2 × 4
##   term        mean_estimate mean_se power
##   <fct>               <dbl>   <dbl> <dbl>
## 1 (Intercept)         803.     15.3  1   
## 2 X_i                  45.9    23.6  0.53

上表為Power分析結果，第二行的Power為0.53是拒絕β₁原假設的概率，是模型中與X_i相關的係數(H₀:β₁=0)。如果我們想了解Power如何隨不同的參數設置而變化，就需要使用single_run()以不同的參數值重新進行模擬。

7 ANOVA檢定

許多研究會計算平均每個受試者在ingroup和outgroup刺激的反應時間，並使用配對樣本T檢定或ANOVA進行比較。本次研究將使用afex::aov_ez進行分析。

# aggregate by subject and analyze with ANOVA
dat_subj <- dat_sim %>%
  group_by(subj_id, category, X_i) %>%
  summarise(RT = mean(RT), .groups = "drop")

afex::aov_ez(
  id = "subj_id",
  dv = "RT",
  within = "category",
  data = dat_subj
)

## Anova Table (Type 3 tests)
## 
## Response: RT
##     Effect    df     MSE         F  ges p.value
## 1 category 1, 99 2971.82 33.52 *** .043   <.001
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

或者按項目，計算每個項目的所有受試者得分的平均值。

# aggregate by item and analyze with ANOVA
dat_item <-  dat_sim %>%
  group_by(item_id, category, X_i) %>%
  summarise(RT = mean(RT), .groups = "drop")

afex::aov_ez(
  id = "item_id",
  dv = "RT",
  between = "category",
  data = dat_item
)

## Converting to factor: category

## Contrasts set to contr.sum for the following variables: category

## Anova Table (Type 3 tests)
## 
## Response: RT
##     Effect    df     MSE      F  ges p.value
## 1 category 1, 48 4506.53 5.53 * .103    .023
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

另外，可以創建一個Power分析函數，函數使用我們從my_sim_data()生成數據的過程來模擬數據，創建兩個數據集，並使用ANOVA對其進行分析，並且只需回傳類別效應的p-values，將Power計算為在模擬當中拒絕虛無假設的百分比。

# power function for ANOVA
my_anova_power <- function(...) {
  dat_sim <- my_sim_data(...)
  
  dat_subj <-  dat_sim %>%
    group_by(subj_id, category, X_i) %>%
    summarise(RT = mean(RT), .groups = "drop")
  
  dat_item <-  dat_sim %>%
    group_by(item_id, category, X_i) %>%
    summarise(RT = mean(RT), .groups = "drop")

  a_subj <- afex::aov_ez(id = "subj_id",
                         dv = "RT",
                         within = "category",
                         data = dat_subj)
  suppressMessages(
    # check contrasts message is annoying
    a_item <- afex::aov_ez(
      id = "item_id",
      dv = "RT",
      between = "category",
      data = dat_item
    )
  )
  
  list(
    "subj" = a_subj$anova_table$`Pr(>F)`,
    "item" = a_item$anova_table$`Pr(>F)`
  )
}

使用默認參數運行此函數以確定每個分析檢測 50 ms 類別效應的能力。

# run simulations and calculate power 
reps <- 100
anova_sims <- purrr::map_df(1:reps, ~my_anova_power())

alpha <- 0.05
power_subj <- mean(anova_sims$subj < alpha)
power_item <- mean(anova_sims$item < alpha)

by-subjects ANOVA具有 0.92 的功效，而by-items ANOVA具有0.54的功效。這不僅是設計內部與設計之間或受試者數量與項目數量的結果，而是某些匯總分析的誤報率過高的結果。將類別的影響設置為0以計算誤報率，這是當沒有實際影響時卻結論為有影響的概率(Type I error)。

# run simulations and calculate the false positive rate 
reps <- 100
anova_fp <- purrr::map_df(1:reps, ~my_anova_power(beta_1 = 0))

false_pos_subj <- mean(anova_fp$subj < alpha)
false_pos_item <- mean(anova_fp$item < alpha)

理想情況下，誤報率將等於alpha，我們在此處將其設置為“r alpha”。按受試者集合分析的誤報率大幅膨脹為“r false_pos_subj”，而按項目集合分析的誤報率接近標準值“r false_pos_item”，這是平均項目和分析項目間因素的結果，也是使用混合效應模型分析交叉因子數據的優勢之一。

8 Sensitivity分析

8.1 設定Power計算函數

# set up the power function
single_run <- function(filename = NULL, ...) {
  # ... is a shortcut that forwards any arguments to my_sim_data()
  dat_sim <- my_sim_data(...)

  # run lmer and capture any warnings
  ww <- ""
  suppressMessages(suppressWarnings(
    mod_sim <- withCallingHandlers({
      lmer(RT ~ X_i + (1 | item_id) + (1 + X_i | subj_id),
                  dat_sim, REML = FALSE)},
      warning = function(w) { ww <<- w$message }
    )
  ))
  
  # get results table and add rep number and any warnings
  sim_results <- broom.mixed::tidy(mod_sim) %>%
    mutate(warnings = ww)
  
  # add columns for the specified parameters
  params <- list(...)
  for (name in names(params)) {
    sim_results[name] <- params[name]
  }
  
  # append the results to a file if filename is set
  if (!is.null(filename)) {
    append <- file.exists(filename) # append if the file exists
    write_csv(sim_results, filename, append = append)
  }

  sim_results
}

8.1.1 範例一：Effect size

範例一將設定一個數據表，包含將測試的所有參數組合，如：從0到100毫秒的類別效應並重複100次，間隔為 10，其他參數為默認值。

filename1 <- "sims/sens1.csv"
nreps <- 100 # number of replications per parameter combo

params <- crossing(
  rep        = 1:nreps, # repeats each combo nreps times
  n_subj     = 100, # number of subjects
  n_ingroup  =  25, # number of ingroup stimuli
  n_outgroup =  25, # number of outgroup stimuli
  beta_0     = 800, # grand mean
  beta_1     = seq(0, 100, by = 10), # effect of category
  omega_0    = 100, # by-item random intercept sd
  tau_0      =  80, # by-subject random intercept sd
  tau_1      =  40, # by-subject random slope sd
  rho        = 0.2, # correlation between intercept and slope
  sigma      = 200  # residual (standard deviation)
) %>%
  select(-rep) # remove rep column

產生的表格有1100行，因此會運行1100次模擬，並利用下面的程式碼將結果保存到命名文件中。

if (!file.exists(filename1)) {
  # run a simulation for each row of params
  # and save to a file on each rep
  sims1 <- purrr::pmap_df(params, single_run, filename = filename1)
}

# read saved simulation data
# NB: col_types is set for warnings in case 
#     the first 1000 rows don't have any
ct <- cols(warnings = col_character(),
           # makes sure plots display in this order
           group = col_factor(ordered = TRUE),
           term = col_factor(ordered = TRUE))
sims1 <- read_csv(filename1, col_types = ct)

下面的程式碼將計算每組的平均估計值和Power。

# calculate mean estimates and power for specified alpha
alpha <- 0.05

power1 <- sims1 %>% 
  filter(effect == "fixed", term == "X_i") %>%
  group_by(term, beta_1) %>%
  summarise(
    mean_estimate = mean(estimate),
    mean_se = mean(std.error),
    power = mean(p.value < alpha),
    .groups = "drop"
  ) 

power1 %>%
  ggplot(aes(beta_1, power)) +
  geom_point() +
  geom_smooth(se = FALSE) +
  ylim(0, 1) +
  scale_x_continuous(name = "Effect of category in ms (beta_1)",
                     breaks = seq(0, 100, 10)) +
  ggtitle("Power for designs varying in effect size")

8.1.2 範例二：受試者與刺激的數量

範例二將模擬10到50個受試者(間隔為10)和10到25個刺激(間隔為5)的20種組合，每一種組合重複50次。

filename2 <- "sims/sens2.csv"
nreps <- 50 # number of replications per parameter combo

params <- crossing(
  rep       = 1:nreps, 
  n_subj    = seq(10, 50, by = 10),
  n_ingroup = seq(10, 25, by = 5),
  beta_0    = 800, # grand mean
  beta_1    = 100, # effect of category
  omega_0   = 100, # by-item random intercept sd
  tau_0     =  80, # by-subject random intercept sd
  tau_1     =  40, # by-subject random slope sd
  rho       = 0.2, # correlation between intercept and slope
  sigma     = 200  # residual (standard deviation)
) %>%
  mutate(n_outgroup = n_ingroup) %>%
  select(-rep) # remove rep column

保存表格到命名文件當中。

if (!file.exists(filename2)) {
  # run a simulation for each row of params
  # and save to a file on each rep
  sims2 <- purrr::pmap_df(params, single_run, filename = filename2)
}

# read saved simulation data
sims2 <- read_csv(filename2, col_types = ct)

計算平均估計值與Power。

# calculate mean estimates and power for specified alpha
alpha <- 0.05

power2 <- sims2 %>% 
  filter(effect == "fixed", term == "X_i") %>%
  group_by(term, n_subj, n_ingroup) %>%
  summarise(
    mean_estimate = mean(estimate),
    mean_se = mean(std.error),
    power = mean(p.value < alpha),
    .groups = "drop"
  ) 

power2 %>%
  ggplot(aes(n_subj, n_ingroup, fill = power)) +
  geom_tile(show.legend = FALSE) +
  geom_text(aes(label = round(power, 2)), 
            color = "black", size = 6) +
  scale_x_continuous(name = "Number of subjects",
                     breaks = seq(10, 50, 10)) +
  scale_y_continuous(name = "Number of items/group",
                     breaks = seq(5, 25, 5)) +
  scale_fill_viridis_c(limits = c(0, 1)) +
  ggtitle("Power for designs varying in number of subjects and items")

8.1.3 範例三：隨機截距SDs

範例三將模擬具有50個受試者、10個刺激以及按刺激和按受試者隨機截取SD的設計，並重複50次，範圍從20到100(間隔20)。

filename3 <- "sims/sens3.csv"
nreps <- 50 # number of replications per parameter combo

params <- crossing(
  rep        = 1:nreps,
  n_subj     =  50, # number of subjects
  n_ingroup  =  10, # number of ingroup items
  n_outgroup =  10, # number of outgroup items
  beta_0     = 800, # grand mean
  beta_1     =  50, # effect of category
  omega_0    = seq(20, 100, by = 20), # by-item random intercept sd
  tau_0      = seq(20, 100, by = 20), # by-subject random intercept sd
  tau_1      =  40, # by-subject random slope sd
  rho        = 0.2, # correlation between intercept and slope
  sigma      = 200  # residual (standard deviation)
) %>%
  select(-rep) # remove rep column