set.seed(12345) # for reproducibility
options(knitr.kable.NA = '')

if (!require("pacman", quietly = TRUE)) install.packages("pacman")
pacman::p_load(stringr, tidyverse, parallel, afex, emmeans, cowplot, scales)
pacman::p_load_gh("RLesur/klippy")
klippy::klippy()
pacman::p_load_gh("ropenscilabs/icon")

set_sum_contrasts() # see Singmann & Kellen (2020)
## setting contr.sum globally: options(contrasts=c('contr.sum', 'contr.poly'))
emm_options(lmer.df = "satterthwaite")
options(dplyr.summarise.inform=FALSE) # 200603 regrouping 메시지 안 보이게 함.

1 Image Statistics

SHINE(spectrum, histogram, and intensity normalization and equalization; Willenbockel et al., 2010) 툴박스를 사용하여 60개 영상(뱀, 얼굴, 꽃 각 20개)에서 물체가 차지한(즉, 배경을 제외한 전경의) 영역을 분석하고, 그 결과를 A와 C에 제시하였다. 원본 영상의 시각적 차이를 줄이기 위해 모든 영상을 회색 톤(greyscale)으로 변환한 후, 히스토그램 매칭과 푸리에 진폭 매칭 절차를 진행하였다. 히스토그램 매칭 절차에서는 개별 영상의 밝기 히스토그램을 전체 영상의 밝기 히스토그램 평균에 일치시킴으로써 영상들의 밝기와 대비 차이를 제거하였다(‘histMatch’ 함수). 푸리에 진폭 매칭 절차에서는 먼저 개별 영상들을 진폭 스펙트럼과 위상 스펙트럼으로 분리하고, 진폭 스펙트럼의 회전 평균(rotational average)을 전체 평균에 일치시켰다. 그다음, 변환된 진폭 스펙트럼과 고유의 위상 스펙트럼을 재결합하여 공간주파수 에너지 차이가 최소화된 영상들을 산출하였다(‘sfMatch’ 함수). 영상 처리로 생성된(‘SHINEd’) 영상들의 저수준 시각특질을 B와 D에 제시하였다. 새 영상들을 구성하는 화소들의 밝기 분포와 공간주파수의 회전 평균이 일치하고 있음을 알 수 있다.

아래는 각 볌주 평균 영상의 luminance histogram: A는 원본(Original) 영상, B는 변환(SHINEd) 영상이다.

lum1 <- read.csv(file="data/bCFS_snake_color.lumHist.csv", header = TRUE) # color
lum2 <- read.csv(file="data/bCFS_snake_shine.lumHist.csv", header = TRUE) # shined
lum1$Category = factor(lum1$Category, levels=c(1,3,2), labels=c("Face","Snake","Flower"))
lum2$Category = factor(lum2$Category, levels=c(1,3,2), labels=c("Face","Snake","Flower"))

head(lum1)
##   Category Value
## 1     Face     0
## 2     Face     0
## 3     Face     0
## 4     Face     0
## 5     Face     0
## 6     Face     0
tail(lum1)
##       Category Value
## 97195    Snake   212
## 97196    Snake   230
## 97197    Snake   231
## 97198    Snake   233
## 97199    Snake   234
## 97200    Snake   234
head(lum2)
##   Category Value
## 1     Face     8
## 2     Face     9
## 3     Face    10
## 4     Face    11
## 5     Face    11
## 6     Face    12
tail(lum2)
##       Category Value
## 97195    Snake   235
## 97196    Snake   236
## 97197    Snake   237
## 97198    Snake   238
## 97199    Snake   239
## 97200    Snake   240

fig.lum1 <- lum1 %>% ggplot(aes(x=Value, color=Category)) +
  geom_freqpoly(binwidth = 1,  alpha = 1) + 
  labs(x = "Luminance", y = "Pixel Count") +
  coord_cartesian(ylim = c(0,1000), clip = "on") +
  scale_y_continuous(breaks = seq(0,1000,200)) +
  theme_bw(base_size = 18) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        aspect.ratio = 1,
        legend.position = c(0.75, 0.85),
        legend.title = element_blank(),
        legend.background = element_blank(),
        legend.key = element_blank())


fig.lum2 <- lum2 %>% ggplot(aes(x=Value, color=Category)) +
  geom_freqpoly(binwidth = 1,  alpha = 1) + 
  labs(x = "Luminance", y = "Pixel Count") +
  coord_cartesian(ylim = c(0,1000), clip = "on") +
  scale_y_continuous(breaks = seq(0,1000,200)) +
  theme_bw(base_size = 18) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        aspect.ratio = 1,
        legend.position = c(0.75, 0.85),
        legend.title = element_blank(),
        legend.background = element_blank(),
        legend.key = element_blank())

cowplot::plot_grid(fig.lum1, fig.lum2, axis = "b",
                   labels = c('A', 'B'), 
                   label_size = 20)

아래는 모든 영상의 푸리에 스펙트라 회전 평균: C는 원본(Original) 영상, D는 변환(SHINEd) 영상이다.

sf1 <- read.csv(file="data/bCFS_snake_color.sfplot.csv", header = TRUE) # color
sf2 <- read.csv(file="data/bCFS_snake_shine.sfplot.csv", header = TRUE) # shined

head(sf1)
##   IMname Category       Energy Frequency
## 1      1        1 474000000000         1
## 2      1        1  31304734183         2
## 3      1        1  17110979051         3
## 4      1        1   6793997665         4
## 5      1        1   3397962914         5
## 6      1        1   1464753523         6
tail(sf1)
##      IMname Category  Energy Frequency
## 5395     60        2 1286716        85
## 5396     60        2 1262494        86
## 5397     60        2 1227971        87
## 5398     60        2 1171561        88
## 5399     60        2 1354170        89
## 5400     60        2 1220981        90
head(sf2)
##   IMname Category       Energy Frequency
## 1      1        1 134000000000         1
## 2      1        1  15876622446         2
## 3      1        1   7807238423         3
## 4      1        1   3260427433         4
## 5      1        1   2247109466         5
## 6      1        1   1052493644         6
tail(sf2)
##      IMname Category   Energy Frequency
## 5395     60        2 429920.8        85
## 5396     60        2 408211.3        86
## 5397     60        2 395509.9        87
## 5398     60        2 405116.0        88
## 5399     60        2 377310.0        89
## 5400     60        2 330469.6        90

sf1$IMname = factor(sf1$IMname)
sf1$Category = factor(sf1$Category, levels=c(1,3,2), labels=c("Face","Snake","Flower"))

sf2$IMname = factor(sf2$IMname)
sf2$Category = factor(sf2$Category, levels=c(1,3,2), labels=c("Face","Snake","Flower"))


fig.sf1 <- ggplot(data = sf1, aes(x=Frequency, y=Energy, group=IMname, color=Category)) +
  geom_line(alpha = 1) +
  scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) +
  scale_y_log10(breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) +
  labs(x = "Spatial Frequency (cycles/image)", y = "Energy") +
  theme_bw(base_size = 18) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        aspect.ratio = 1,
        legend.position = c(0.75, 0.85),
        legend.title = element_blank(),
        legend.background = element_blank(),
        legend.key = element_blank())

fig.sf2 <- ggplot(data = sf2, aes(x=Frequency, y=Energy, group=IMname, color=Category)) +
  geom_line(alpha = 1) +
  scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) +
  scale_y_log10(breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) +
  labs(x = "Spatial Frequency (cycles/image)", y = "Energy") +
  theme_bw(base_size = 18) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        aspect.ratio = 1,
        legend.position = c(0.75, 0.85),
        legend.title = element_blank(),
        legend.background = element_blank(),
        legend.key = element_blank())

cowplot::plot_grid(fig.sf1, fig.sf2, axis = "b",
                   labels = c('C', 'D'), 
                   label_size = 20)



2 Experiment 1: Rating

원자극(Original)을 평정한 28명의 자료와 통제된 자극(SHINEd)을 평정한 30명의 자료를 합쳐 분석한다.

S1 <- read.csv(file="data/bCFS_snake_rating_color.csv", header = TRUE)
S2 <- read.csv(file="data/bCFS_snake_rating_shine.csv", header = TRUE)

S1$VisAttr <- factor("Original") # 실험 요인 
S2$VisAttr <- factor("SHINEd")

S1$Subject <- factor(paste('e1s', str_pad(S1$Subject, 2, pad = "0"), sep = "")) # 참가자 
S2$Subject <- factor(paste('e2s', str_pad(S2$Subject, 2, pad = "0"), sep = ""))

S1$Image <- factor(str_replace_all(S1$Image, fixed(" "), ""))  # 파일명 앞 빈칸 삭제 
S2$Image <- factor(str_replace_all(S2$Image, fixed(" "), ""))

dim(S1)
## [1] 1680    8
dim(S2)
## [1] 1800    8

SS <- rbind(S1, S2) # 두 실험 자료를 합친다. 
head(SS)
##   Subject Trial Category         Image Valence Arousal IMidx  VisAttr
## 1   e1s01     1        2 flower011.png       9       8    54 Original
## 2   e1s01     3        3  snake019.png       5       5    82 Original
## 3   e1s01     7        3  snake008.png       4       7    71 Original
## 4   e1s01     8        2 flower003.png       7       8    46 Original
## 5   e1s01    10        1   face003.png       1       9     9 Original
## 6   e1s01    11        2 flower020.png       8       8    63 Original
tail(SS)
##      Subject Trial Category        Image Valence Arousal IMidx VisAttr
## 3475   e2s31    46        3 snake015.png       1       9    55  SHINEd
## 3476   e2s31    43        3 snake016.png       1       9    56  SHINEd
## 3477   e2s31    33        3 snake017.png       1       9    57  SHINEd
## 3478   e2s31    36        3 snake018.png       1       9    58  SHINEd
## 3479   e2s31    51        3 snake019.png       1       9    59  SHINEd
## 3480   e2s31    12        3 snake020.png       1       9    60  SHINEd

SS <- SS %>% select(VisAttr, Subject, Image, Category, Valence, Arousal)
SS$Category <- factor(SS$Category, levels=c(1,3,2), labels=c("Face","Snake","Flower")) 

str(SS) 
## 'data.frame':    3480 obs. of  6 variables:
##  $ VisAttr : Factor w/ 2 levels "Original","SHINEd": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Subject : Factor w/ 58 levels "e1s01","e1s02",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Image   : Factor w/ 60 levels "face001.png",..: 31 59 48 23 3 40 21 4 11 42 ...
##  $ Category: Factor w/ 3 levels "Face","Snake",..: 3 2 2 3 1 3 3 1 1 2 ...
##  $ Valence : int  9 5 4 7 1 8 7 1 1 3 ...
##  $ Arousal : int  8 5 7 8 9 8 7 8 8 7 ...
# 1. VisAttr : 1_color, 2_shiend. Factor임을 확인할 것.
# 2. Subject ID. Factor임을 확인할 것.
# 3. Image ID. Factor임을 확인할 것.
# 4. Category : stimulus category. 1-face, 2-flower, 3-snake
# 5. Valence : 1(negative) ~ 9(positive)
# 6. Arousal : 1(weak) ~ 9(strong)

table(SS$Subject, SS$VisAttr)
##        
##         Original SHINEd
##   e1s01       60      0
##   e1s02       60      0
##   e1s03       60      0
##   e1s04       60      0
##   e1s05       60      0
##   e1s06       60      0
##   e1s07       60      0
##   e1s08       60      0
##   e1s09       60      0
##   e1s10       60      0
##   e1s11       60      0
##   e1s12       60      0
##   e1s13       60      0
##   e1s14       60      0
##   e1s15       60      0
##   e1s16       60      0
##   e1s17       60      0
##   e1s18       60      0
##   e1s19       60      0
##   e1s20       60      0
##   e1s21       60      0
##   e1s22       60      0
##   e1s23       60      0
##   e1s24       60      0
##   e1s25       60      0
##   e1s26       60      0
##   e1s27       60      0
##   e1s28       60      0
##   e2s01        0     60
##   e2s02        0     60
##   e2s03        0     60
##   e2s04        0     60
##   e2s05        0     60
##   e2s06        0     60
##   e2s07        0     60
##   e2s08        0     60
##   e2s09        0     60
##   e2s10        0     60
##   e2s11        0     60
##   e2s12        0     60
##   e2s13        0     60
##   e2s14        0     60
##   e2s15        0     60
##   e2s16        0     60
##   e2s17        0     60
##   e2s18        0     60
##   e2s19        0     60
##   e2s20        0     60
##   e2s21        0     60
##   e2s23        0     60
##   e2s24        0     60
##   e2s25        0     60
##   e2s26        0     60
##   e2s27        0     60
##   e2s28        0     60
##   e2s29        0     60
##   e2s30        0     60
##   e2s31        0     60
table(SS$Image, SS$VisAttr)
##                
##                 Original SHINEd
##   face001.png         28     30
##   face002.png         28     30
##   face003.png         28     30
##   face004.png         28     30
##   face005.png         28     30
##   face006.png         28     30
##   face007.png         28     30
##   face008.png         28     30
##   face009.png         28     30
##   face010.png         28     30
##   face011.png         28     30
##   face012.png         28     30
##   face013.png         28     30
##   face014.png         28     30
##   face015.png         28     30
##   face016.png         28     30
##   face017.png         28     30
##   face018.png         28     30
##   face019.png         28     30
##   face020.png         28     30
##   flower001.png       28     30
##   flower002.png       28     30
##   flower003.png       28     30
##   flower004.png       28     30
##   flower005.png       28     30
##   flower006.png       28     30
##   flower007.png       28     30
##   flower008.png       28     30
##   flower009.png       28     30
##   flower010.png       28     30
##   flower011.png       28     30
##   flower012.png       28     30
##   flower013.png       28     30
##   flower014.png       28     30
##   flower015.png       28     30
##   flower016.png       28     30
##   flower017.png       28     30
##   flower018.png       28     30
##   flower019.png       28     30
##   flower020.png       28     30
##   snake001.png        28     30
##   snake002.png        28     30
##   snake003.png        28     30
##   snake004.png        28     30
##   snake005.png        28     30
##   snake006.png        28     30
##   snake007.png        28     30
##   snake008.png        28     30
##   snake009.png        28     30
##   snake010.png        28     30
##   snake011.png        28     30
##   snake012.png        28     30
##   snake013.png        28     30
##   snake014.png        28     30
##   snake015.png        28     30
##   snake016.png        28     30
##   snake017.png        28     30
##   snake018.png        28     30
##   snake019.png        28     30
##   snake020.png        28     30
table(SS$Category, SS$Subject)
##         
##          e1s01 e1s02 e1s03 e1s04 e1s05 e1s06 e1s07 e1s08 e1s09 e1s10 e1s11
##   Face      20    20    20    20    20    20    20    20    20    20    20
##   Snake     20    20    20    20    20    20    20    20    20    20    20
##   Flower    20    20    20    20    20    20    20    20    20    20    20
##         
##          e1s12 e1s13 e1s14 e1s15 e1s16 e1s17 e1s18 e1s19 e1s20 e1s21 e1s22
##   Face      20    20    20    20    20    20    20    20    20    20    20
##   Snake     20    20    20    20    20    20    20    20    20    20    20
##   Flower    20    20    20    20    20    20    20    20    20    20    20
##         
##          e1s23 e1s24 e1s25 e1s26 e1s27 e1s28 e2s01 e2s02 e2s03 e2s04 e2s05
##   Face      20    20    20    20    20    20    20    20    20    20    20
##   Snake     20    20    20    20    20    20    20    20    20    20    20
##   Flower    20    20    20    20    20    20    20    20    20    20    20
##         
##          e2s06 e2s07 e2s08 e2s09 e2s10 e2s11 e2s12 e2s13 e2s14 e2s15 e2s16
##   Face      20    20    20    20    20    20    20    20    20    20    20
##   Snake     20    20    20    20    20    20    20    20    20    20    20
##   Flower    20    20    20    20    20    20    20    20    20    20    20
##         
##          e2s17 e2s18 e2s19 e2s20 e2s21 e2s23 e2s24 e2s25 e2s26 e2s27 e2s28
##   Face      20    20    20    20    20    20    20    20    20    20    20
##   Snake     20    20    20    20    20    20    20    20    20    20    20
##   Flower    20    20    20    20    20    20    20    20    20    20    20
##         
##          e2s29 e2s30 e2s31
##   Face      20    20    20
##   Snake     20    20    20
##   Flower    20    20    20

contrasts(SS$Category)
##        [,1] [,2]
## Face      1    0
## Snake     0    1
## Flower   -1   -1
contrasts(SS$VisAttr)
##          [,1]
## Original    1
## SHINEd     -1



2.1 Valence

정서가를 분석한다.

SS %>% group_by(VisAttr, Category, Subject) %>% # subject-level
  summarise(mn = mean(Valence), .groups="keep") %>% 
  ungroup() %>% 
  group_by(VisAttr, Category) %>% # group-level
  summarise(Mean = mean(mn),
            Std = sd(mn)) %>% 
  ungroup() %>% 
  knitr::kable()
VisAttr Category Mean Std
Original Face 2.464286 0.6595252
Original Snake 2.192857 1.1480026
Original Flower 6.973214 1.0770744
SHINEd Face 2.835000 0.9143634
SHINEd Snake 3.653333 1.1345615
SHINEd Flower 6.030000 0.9782638

Category(뱀, 얼굴, 꽃)와 VisAttr(Low-Level Visual Attributes; Original vs. SHINEd)를 고정효과로 afex::mixed()를 통해 선형혼합모형 분석을 실시하였다. 참가자(Subject)와 실험자극(Image)의 무선효과가 절편으로 포함되었다. Category는 참가자내 변인이고, VisAttr는 자극내 편인이므로 각각 기울기의 무선효과로 모형에 포함되었다.

(nc <- detectCores())
## [1] 4
cl <- makeCluster(rep("localhost", nc))
SSv.lmer <- afex::mixed(Valence ~ Category*VisAttr + (Category|Subject) + (VisAttr|Image), 
                        data = SS, method = "S", cl = cl,
                        control = lmerControl(optimizer = "bobyqa",
                                              optCtrl = list(maxfun = 1e6)))
## Fitting one lmer() model. [DONE]
## Calculating p-values. [DONE]
stopCluster(cl)
anova(SSv.lmer)
## Mixed Model Anova Table (Type 3 tests, S-method)
## 
## Model: Valence ~ Category * VisAttr + (Category | Subject) + (VisAttr | 
## Model:     Image)
## Data: SS
##                  num Df den Df        F    Pr(>F)    
## Category              2 72.461 181.9140 < 2.2e-16 ***
## VisAttr               1 64.257   4.9849   0.02906 *  
## Category:VisAttr      2 62.479  12.3349 3.055e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(SSv.lmer)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Valence ~ Category * VisAttr + (Category | Subject) + (VisAttr |  
##     Image)
##    Data: data
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+06))
## 
## REML criterion at convergence: 11155.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.0605 -0.5388 -0.0594  0.4474  6.8546 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr       
##  Image    (Intercept) 0.05584  0.2363              
##           VisAttr1    0.02111  0.1453   0.33       
##  Subject  (Intercept) 0.21383  0.4624              
##           Category1   0.39852  0.6313   -0.05      
##           Category2   0.80637  0.8980    0.27 -0.21
##  Residual             1.21951  1.1043              
## Number of obs: 3480, groups:  Image, 60; Subject, 58
## 
## Fixed effects:
##                    Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)         4.02478    0.07052 77.12863  57.076  < 2e-16 ***
## Category1          -1.37514    0.09717 78.09250 -14.152  < 2e-16 ***
## Category2          -1.10169    0.12838 69.07693  -8.581 1.72e-12 ***
## VisAttr1           -0.14800    0.06629 64.25724  -2.233   0.0291 *  
## Category1:VisAttr1 -0.03736    0.09102 64.68499  -0.410   0.6828    
## Category2:VisAttr1 -0.58224    0.12379 60.96440  -4.703 1.51e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Ctgry1 Ctgry2 VsAtt1 C1:VA1
## Category1   -0.039                            
## Category2    0.211 -0.265                     
## VisAttr1     0.070 -0.001  0.008              
## Ctgry1:VsA1 -0.001  0.072 -0.023 -0.044       
## Ctgry2:VsA1  0.008 -0.022  0.055  0.232 -0.242

범주, 실험, 범주x실험 상호작용 효과가 모두 유의미하였다.

SSv.lmer.emm <- emmeans(SSv.lmer, "Category", by = "VisAttr", lmerTest.limit = 3480)
plot(SSv.lmer.emm, horizontal = FALSE, type = "response") + facet_grid(cols = vars(VisAttr))

SSv.lmer.emm.pair <- update(pairs(SSv.lmer.emm), by = NULL, adjust = "holm") 
# 'by=NULL'은 familywise error rate correction에 필요.
SSv.lmer.emm.pair %>% summary(infer = TRUE)
##  contrast       VisAttr  estimate    SE   df lower.CL upper.CL t.ratio p.value
##  Face - Snake   Original    0.271 0.256 73.0   -0.424    0.967   1.059 0.2932 
##  Face - Flower  Original   -4.509 0.289 69.6   -5.293   -3.725 -15.624 <.0001 
##  Snake - Flower Original   -4.780 0.356 65.0   -5.750   -3.811 -13.415 <.0001 
##  Face - Snake   SHINEd     -0.818 0.240 66.1   -1.470   -0.166  -3.413 0.0022 
##  Face - Flower  SHINEd     -3.195 0.272 64.0   -3.935   -2.455 -11.755 <.0001 
##  Snake - Flower SHINEd     -2.377 0.339 61.2   -3.300   -1.453  -7.019 <.0001 
## 
## Degrees-of-freedom method: satterthwaite 
## Confidence level used: 0.95 
## Conf-level adjustment: bonferroni method for 6 estimates 
## P value adjustment: holm method for 6 tests

시각특질을 통제한 후, 정서가의 극단값이 감소하였다. 원자극을 썼을 때는 뱀과 얼굴의 정서가가 비슷했으나 통제한 후에는 뱀의 정서가가 덜 부정적인 방향으로 바뀌었다.



2.2 Arousal

각성 점수를 분석하였다.

SS %>% group_by(VisAttr, Category, Subject) %>% # subject-level
  summarise(mn = mean(Arousal)) %>% 
  ungroup() %>% 
  group_by(VisAttr, Category) %>% # group-level
  summarise(Mean = mean(mn),
            Std = sd(mn)) %>% 
  ungroup() %>% 
  knitr::kable()
VisAttr Category Mean Std
Original Face 6.819643 1.085198
Original Snake 6.983929 1.270320
Original Flower 5.941071 1.479117
SHINEd Face 5.926667 1.635898
SHINEd Snake 4.465000 1.902274
SHINEd Flower 4.511667 1.821173

각성 점수 분석을 위한 선형혼합모형은 이전과 같다.

cl <- makeCluster(rep("localhost", nc))
SSa.lmer <- afex::mixed(Arousal ~ Category*VisAttr + (Category|Subject) + (VisAttr|Image), 
                        data = SS, method = "S", cl = cl,
                        control = lmerControl(optimizer = "bobyqa",
                                              optCtrl = list(maxfun = 1e6)))
## Fitting one lmer() model. [DONE]
## Calculating p-values. [DONE]
stopCluster(cl)
anova(SSa.lmer)
## Mixed Model Anova Table (Type 3 tests, S-method)
## 
## Model: Arousal ~ Category * VisAttr + (Category | Subject) + (VisAttr | 
## Model:     Image)
## Data: SS
##                  num Df den Df       F    Pr(>F)    
## Category              2 68.796 13.3563 1.256e-05 ***
## VisAttr               1 58.444 23.2526 1.053e-05 ***
## Category:VisAttr      2 64.220  6.2886  0.003205 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(SSa.lmer)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Arousal ~ Category * VisAttr + (Category | Subject) + (VisAttr |  
##     Image)
##    Data: data
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+06))
## 
## REML criterion at convergence: 13084.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.8898 -0.4740  0.0359  0.5587  3.5413 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr       
##  Image    (Intercept) 0.05808  0.2410              
##           VisAttr1    0.03713  0.1927   0.62       
##  Subject  (Intercept) 1.55089  1.2453              
##           Category1   1.00050  1.0002   -0.28      
##           Category2   0.81083  0.9005    0.08 -0.67
##  Residual             2.11830  1.4554              
## Number of obs: 3480, groups:  Image, 60; Subject, 58
## 
## Fixed effects:
##                    Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)          5.7747     0.1684 59.8335  34.297  < 2e-16 ***
## Category1            0.5985     0.1429 66.4351   4.188 8.48e-05 ***
## Category2           -0.0502     0.1310 68.2654  -0.383 0.702685    
## VisAttr1             0.8069     0.1673 58.4437   4.822 1.05e-05 ***
## Category1:VisAttr1  -0.3604     0.1404 62.6566  -2.566 0.012693 *  
## Category2:VisAttr1   0.4526     0.1283 63.8388   3.528 0.000781 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Ctgry1 Ctgry2 VsAtt1 C1:VA1
## Category1   -0.252                            
## Category2    0.072 -0.645                     
## VisAttr1     0.051 -0.009  0.003              
## Ctgry1:VsA1 -0.009  0.080 -0.047 -0.258       
## Ctgry2:VsA1  0.003 -0.047  0.089  0.074 -0.650

범주, 실험, 범주x실험 상호작용 효과가 모두 유의미하였다.

SSa.lmer.emm <- emmeans(SSa.lmer, "Category", by = "VisAttr", lmerTest.limit = 3480)
plot(SSa.lmer.emm, horizontal = FALSE, type = "response") + facet_grid(cols = vars(VisAttr))

SSa.lmer.emm.pair <- update(pairs(SSa.lmer.emm), by = NULL, adjust = "holm")
SSa.lmer.emm.pair %>% summary(infer = TRUE)
##  contrast       VisAttr  estimate    SE   df lower.CL upper.CL t.ratio p.value
##  Face - Snake   Original  -0.1643 0.362 69.1   -1.147    0.819 -0.454  1.0000 
##  Face - Flower  Original   0.8786 0.329 71.7   -0.013    1.770  2.674  0.0278 
##  Snake - Flower Original   1.0429 0.296 74.8    0.240    1.845  3.523  0.0029 
##  Face - Snake   SHINEd     1.4617 0.334 59.4    0.550    2.374  4.374  0.0003 
##  Face - Flower  SHINEd     1.4150 0.300 60.0    0.595    2.235  4.711  0.0001 
##  Snake - Flower SHINEd    -0.0467 0.267 60.8   -0.775    0.681 -0.175  1.0000 
## 
## Degrees-of-freedom method: satterthwaite 
## Confidence level used: 0.95 
## Conf-level adjustment: bonferroni method for 6 estimates 
## P value adjustment: holm method for 6 tests

뱀의 각성 점수는 시각특질을 통제하지 않으면 얼굴과 비슷하고, 통제하면 꽃과 비슷했다.





3 Experiment 2: Breaking CFS

원자극을 사용한 실험과 통제 자극을 사용한 실험 결과를 합쳐 분석한다. 각각 20명씩 총 40명의 자료이다.

E1 <- read.csv(file="data/bCFS_snake_1_color.csv", header = TRUE)
E2 <- read.csv(file="data/bCFS_snake_2_shine.csv", header = TRUE)

E1$VisAttr <- factor("Original") # 실험 요인
E2$VisAttr <- factor("SHINEd")

E1$Subject <- factor(paste('e1s', str_pad(E1$Subject, 2, pad = "0"), sep = "")) # 참가자 
E2$Subject <- factor(paste('e2s', str_pad(E2$Subject, 2, pad = "0"), sep = ""))

E1$Image <- factor(str_replace_all(E1$Image, fixed(" "), ""))  # 파일명 앞 빈칸 삭제 
E2$Image <- factor(str_replace_all(E2$Image, fixed(" "), ""))

E1$LocCorr <- as.numeric(E1$LocCorr==1) # binarizing
E1$CatCorr <- as.numeric(E1$CatCorr==1) 
E2$LocCorr <- as.numeric(E2$LocCorr==1)
E2$CatCorr <- as.numeric(E2$CatCorr==1) 

E1$RT <- E1$RT*1000
E2$RT <- E2$RT*1000

EE <- rbind(E1, E2) # 두 실험 자료를 합친다. 
head(EE)
##   Subject Trial Block BlockTrial Category Location         Image ImgContrast
## 1   e1s01     1     1          1        2        1 flower020.png           1
## 2   e1s01     2     1          2        1        2   face001.png           1
## 3   e1s01     3     1          3        1        1   face002.png           1
## 4   e1s01     4     1          4        2        2 flower012.png           1
## 5   e1s01     5     1          5        2        2 flower016.png           1
## 6   e1s01     6     1          6        1        2   face007.png           1
##       RT LocResp LocCorr CatResp CatCorr tindex  VisAttr
## 1 2072.4       1       1       2       1    116 Original
## 2 1838.7       2       1       1       1      4 Original
## 3 1470.6       1       1       1       1      7 Original
## 4 4109.9       2       1       2       1     71 Original
## 5 4477.1       2       1       2       1     95 Original
## 6 2306.5       2       1       1       1     40 Original
tail(EE)
##       Subject Trial Block BlockTrial Category Location         Image
## 14395   e2s20   355     3        115        2        2 flower006.png
## 14396   e2s20   356     3        116        1        2   face004.png
## 14397   e2s20   357     3        117        2        1 flower007.png
## 14398   e2s20   358     3        118        3        1  snake003.png
## 14399   e2s20   359     3        119        3        1  snake010.png
## 14400   e2s20   360     3        120        2        1 flower013.png
##       ImgContrast     RT LocResp LocCorr CatResp CatCorr tindex VisAttr
## 14395           1 3492.1       1       0       7       0    275  SHINEd
## 14396           1 6264.0       1       0       7       0    262  SHINEd
## 14397           1 5579.2       1       1       3       0    278  SHINEd
## 14398           1 1137.6       1       1       7       0    255  SHINEd
## 14399           1 6162.7       1       1       3       1    297  SHINEd
## 14400           1 2856.5       1       1       3       0    314  SHINEd

EE <- EE %>% select(VisAttr, Subject, Image, Category, RT, LocCorr, CatCorr)
EE$Category <- factor(EE$Category, levels=c(1,3,2), labels=c("Face","Snake","Flower")) 

str(EE) # 실험, 참가자, 이미지 factor 여부와 level 확인. 
## 'data.frame':    14400 obs. of  7 variables:
##  $ VisAttr : Factor w/ 2 levels "Original","SHINEd": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Subject : Factor w/ 40 levels "e1s01","e1s02",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Image   : Factor w/ 60 levels "face001.png",..: 40 1 2 32 36 7 55 21 57 20 ...
##  $ Category: Factor w/ 3 levels "Face","Snake",..: 3 1 1 3 3 1 2 3 2 1 ...
##  $ RT      : num  2072 1839 1471 4110 4477 ...
##  $ LocCorr : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ CatCorr : num  1 1 1 1 1 1 1 1 1 1 ...

table(EE$Subject, EE$VisAttr)
##        
##         Original SHINEd
##   e1s01      360      0
##   e1s02      360      0
##   e1s03      360      0
##   e1s04      360      0
##   e1s05      360      0
##   e1s06      360      0
##   e1s07      360      0
##   e1s08      360      0
##   e1s09      360      0
##   e1s10      360      0
##   e1s11      360      0
##   e1s12      360      0
##   e1s13      360      0
##   e1s14      360      0
##   e1s15      360      0
##   e1s16      360      0
##   e1s17      360      0
##   e1s18      360      0
##   e1s19      360      0
##   e1s20      360      0
##   e2s01        0    360
##   e2s02        0    360
##   e2s03        0    360
##   e2s04        0    360
##   e2s05        0    360
##   e2s06        0    360
##   e2s07        0    360
##   e2s08        0    360
##   e2s09        0    360
##   e2s10        0    360
##   e2s11        0    360
##   e2s12        0    360
##   e2s13        0    360
##   e2s14        0    360
##   e2s15        0    360
##   e2s16        0    360
##   e2s17        0    360
##   e2s18        0    360
##   e2s19        0    360
##   e2s20        0    360
table(EE$Image, EE$VisAttr)
##                
##                 Original SHINEd
##   face001.png        120    120
##   face002.png        120    120
##   face003.png        120    120
##   face004.png        120    120
##   face005.png        120    120
##   face006.png        120    120
##   face007.png        120    120
##   face008.png        120    120
##   face009.png        120    120
##   face010.png        120    120
##   face011.png        120    120
##   face012.png        120    120
##   face013.png        120    120
##   face014.png        120    120
##   face015.png        120    120
##   face016.png        120    120
##   face017.png        120    120
##   face018.png        120    120
##   face019.png        120    126
##   face020.png        120    114
##   flower001.png      120    120
##   flower002.png      120    120
##   flower003.png      120    120
##   flower004.png      120    120
##   flower005.png      120    120
##   flower006.png      120    120
##   flower007.png      120    120
##   flower008.png      120    120
##   flower009.png      120    120
##   flower010.png      120    120
##   flower011.png      120    120
##   flower012.png      120    120
##   flower013.png      120    120
##   flower014.png      120    120
##   flower015.png      120    120
##   flower016.png      120    120
##   flower017.png      120    120
##   flower018.png      120    120
##   flower019.png      120    120
##   flower020.png      120    120
##   snake001.png       120    120
##   snake002.png       120    120
##   snake003.png       120    120
##   snake004.png       120    120
##   snake005.png       120    120
##   snake006.png       120    120
##   snake007.png       120    120
##   snake008.png       120    120
##   snake009.png       120    120
##   snake010.png       120    120
##   snake011.png       120    120
##   snake012.png       120    120
##   snake013.png       120    120
##   snake014.png       120    120
##   snake015.png       120    120
##   snake016.png       120    120
##   snake017.png       120    120
##   snake018.png       120    120
##   snake019.png       120    120
##   snake020.png       120    120

contrasts(EE$Category)
##        [,1] [,2]
## Face      1    0
## Snake     0    1
## Flower   -1   -1
contrasts(EE$VisAttr)
##          [,1]
## Original    1
## SHINEd     -1



3.1 Localization Accuracy

표적의 위치 판단의 정확률을 분석하였다.

EE %>% group_by(VisAttr, Category, Subject) %>% # subject-level
  summarise(mn = mean(LocCorr)) %>% 
  ungroup() %>% 
  group_by(VisAttr, Category) %>% # group-level
  summarise(Mean = mean(mn),
            Std = sd(mn)) %>% 
  ungroup() %>% 
  knitr::kable()
VisAttr Category Mean Std
Original Face 0.9541667 0.0390400
Original Snake 0.9545833 0.0534261
Original Flower 0.9070833 0.0872448
SHINEd Face 0.9620833 0.0560959
SHINEd Snake 0.9537500 0.0754361
SHINEd Flower 0.9350000 0.0993620 
cl <- makeCluster(rep("localhost", nc))
EEloc.glmer <- afex::mixed(LocCorr ~ Category*VisAttr + (Category|Subject) + (VisAttr|Image), 
                           data = EE, method = "LRT", cl = cl,
                           family = binomial(link="logit"), 
                           control = glmerControl(optimizer = "bobyqa",
                                                  optCtrl = list(maxfun = 1e6)))
## Fitting 4 (g)lmer() models.
stopCluster(cl)

Category의 효과만 유의하였다. Original 자극을 썼을 때 뱀과 꽃의 차이만 있었다. 아래 사후분석 결과도 같다.

anova(EEloc.glmer) # EEloc.glamer$anova_table을 프린트한다.
## Mixed Model Anova Table (Type 3 tests, LRT-method)
## 
## Model: LocCorr ~ Category * VisAttr + (Category | Subject) + (VisAttr | 
## Model:     Image)
## Data: EE
## Df full model: 15
##                  Df  Chisq Chi Df Pr(>Chisq)  
## Category         13 7.8053      2    0.02019 *
## VisAttr          14 2.0555      1    0.15166  
## Category:VisAttr 13 1.7362      2    0.41975  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(EEloc.glmer)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: LocCorr ~ Category * VisAttr + (Category | Subject) + (VisAttr |  
##     Image)
##    Data: data
## Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+06))
## 
##      AIC      BIC   logLik deviance df.resid 
##   5371.6   5485.3  -2670.8   5341.6    14385 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -10.5779   0.1181   0.1677   0.2337   1.2128 
## 
## Random effects:
##  Groups  Name        Variance Std.Dev. Corr       
##  Image   (Intercept) 0.08665  0.2944              
##          VisAttr1    0.10727  0.3275   0.28       
##  Subject (Intercept) 0.98424  0.9921              
##          Category1   0.32752  0.5723   -0.45      
##          Category2   0.03812  0.1952    0.64 -0.62
## Number of obs: 14400, groups:  Image, 60; Subject, 40
## 
## Fixed effects:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         3.47315    0.17232  20.156   <2e-16 ***
## Category1           0.10676    0.13264   0.805   0.4209    
## Category2           0.21376    0.10279   2.080   0.0376 *  
## VisAttr1           -0.24794    0.17181  -1.443   0.1490    
## Category1:VisAttr1  0.02913    0.12943   0.225   0.8220    
## Category2:VisAttr1  0.09598    0.09323   1.030   0.3032    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Ctgry1 Ctgry2 VsAtt1 C1:VA1
## Category1   -0.284                            
## Category2    0.213 -0.511                     
## VisAttr1    -0.019  0.003  0.009              
## Ctgry1:VsA1  0.004 -0.037  0.003 -0.292       
## Ctgry2:VsA1  0.005  0.004 -0.002  0.217 -0.496
EEloc.glmer.emm <- emmeans(EEloc.glmer, "Category", by = "VisAttr", 
                           lmerTest.limit = 3480, type = "response")
plot(EEloc.glmer.emm, horizontal = FALSE) + facet_grid(cols = vars(VisAttr))

EEloc.glmer.emm.pair <- update(pairs(EEloc.glmer.emm), by = NULL, adjust = "holm")
EEloc.glmer.emm.pair %>% summary(infer = TRUE)
##  contrast       VisAttr  odds.ratio    SE  df asymp.LCL asymp.UCL z.ratio
##  Face / Snake   Original      0.840 0.235 Inf     0.402      1.76 -0.622 
##  Face / Flower  Original      1.789 0.566 Inf     0.777      4.12  1.839 
##  Snake / Flower Original      2.128 0.515 Inf     1.125      4.03  3.125 
##  Face / Snake   SHINEd        0.961 0.273 Inf     0.453      2.04 -0.141 
##  Face / Flower  SHINEd        1.314 0.435 Inf     0.549      3.15  0.825 
##  Snake / Flower SHINEd        1.368 0.337 Inf     0.714      2.62  1.272 
##  p.value
##  1.0000 
##  0.3296 
##  0.0107 
##  1.0000 
##  1.0000 
##  0.8132 
## 
## Confidence level used: 0.95 
## Conf-level adjustment: bonferroni method for 6 estimates 
## Intervals are back-transformed from the log odds ratio scale 
## P value adjustment: holm method for 6 tests 
## Tests are performed on the log odds ratio scale



3.2 Breaking (= Localization) Time

표적 위치를 보고하는 데 걸린 시간을 분석하였다.

numCorr <- EE %>% filter(LocCorr==1) %>% nrow() # correct RTs
numAntici <- EE %>% filter(LocCorr==1 & RT < 200) %>% nrow() # anticipatory correct RTs

EEcorr <- EE %>% filter(LocCorr == 1 & RT > 200)
EEcorr %>% group_by(VisAttr, Category, Subject) %>% # subject-level
  summarise(mn = mean(RT)) %>% 
  ungroup() %>% 
  group_by(VisAttr, Category) %>% # group-level
  summarise(Mean = mean(mn),
            Std = sd(mn)) %>% 
  ungroup() %>% 
  knitr::kable()
VisAttr Category Mean Std
Original Face 1523.138 415.5439
Original Snake 2112.743 717.7735
Original Flower 2673.397 763.0071
SHINEd Face 1767.412 489.4127
SHINEd Snake 2490.043 677.2309
SHINEd Flower 2585.240 786.3497

전체 13600개의 정반응 중에서 0.2초보다 빠른 반응은 2개였다(0.0147059%).
Breaking Time이 정규 분포를 이루지 않으므로 log 값을 변환한다.

rt.dens <- ggplot(EEcorr, aes(x=RT)) + 
  geom_density() + 
  theme_bw(base_size = 18) +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank())

rt.log.dens <- EEcorr %>% 
  mutate(logRT = log(RT)) %>% 
  ggplot(aes(x=logRT)) + 
  geom_density() + 
  theme_bw(base_size = 18) + 
  labs(x = "Log(RT)") + 
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank()) 

plot_grid(rt.dens, rt.log.dens, ncol = 2, label_size = 20)

cl <- makeCluster(rep("localhost", nc))
EErt.log <- afex::mixed(log(RT) ~ Category*VisAttr + (Category|Subject) + (VisAttr|Image), 
                           data = EEcorr, method = "S", cl = cl,
                           control = lmerControl(optimizer = "bobyqa",
                                                 optCtrl = list(maxfun = 1e8)))
## Fitting one lmer() model. [DONE]
## Calculating p-values. [DONE]
stopCluster(cl)
anova(EErt.log)
## Mixed Model Anova Table (Type 3 tests, S-method)
## 
## Model: log(RT) ~ Category * VisAttr + (Category | Subject) + (VisAttr | 
## Model:     Image)
## Data: EEcorr
##                  num Df den Df       F    Pr(>F)    
## Category              2 72.695 64.3754 < 2.2e-16 ***
## VisAttr               1 41.318  1.7239  0.196435    
## Category:VisAttr      2 72.613  7.3259  0.001264 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(EErt.log)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(RT) ~ Category * VisAttr + (Category | Subject) + (VisAttr |  
##     Image)
##    Data: data
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+08))
## 
## REML criterion at convergence: 14798.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.2267 -0.6731 -0.1047  0.5645  4.7671 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr       
##  Image    (Intercept) 0.009171 0.09577             
##           VisAttr1    0.004558 0.06751  0.59       
##  Subject  (Intercept) 0.068522 0.26177             
##           Category1   0.005029 0.07092  -0.54      
##           Category2   0.001155 0.03399   0.48 -0.56
##  Residual             0.166772 0.40838             
## Number of obs: 13598, groups:  Image, 60; Subject, 40
## 
## Fixed effects:
##                    Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)         7.53875    0.04334 44.72497 173.949  < 2e-16 ***
## Category1          -0.23428    0.02135 82.61387 -10.972  < 2e-16 ***
## Category2           0.05440    0.01895 64.59626   2.871  0.00553 ** 
## VisAttr1           -0.05573    0.04244 41.31836  -1.313  0.19644    
## Category1:VisAttr1 -0.01974    0.01738 79.71480  -1.136  0.25950    
## Category2:VisAttr1 -0.03772    0.01433 67.14245  -2.633  0.01051 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Ctgry1 Ctgry2 VsAtt1 C1:VA1
## Category1   -0.271                            
## Category2    0.129 -0.491                     
## VisAttr1     0.035  0.000  0.000              
## Ctgry1:VsA1  0.000  0.343 -0.193 -0.340       
## Ctgry2:VsA1  0.000 -0.208  0.469  0.174 -0.489

Category와 Category x VisAttr 효과가 유의하였다.

EErt.log.emm <- emmeans(EErt.log, "Category", by = "VisAttr", 
                        lmerTest.limit = 13598, type = "response")
plot(EErt.log.emm, horizontal = FALSE) + facet_grid(cols = vars(VisAttr))

EErt.log.emm.pair <- update(pairs(EErt.log.emm), by = NULL, adjust = "holm")
EErt.log.emm.pair %>% summary(infer = TRUE)
##  contrast       VisAttr  ratio     SE   df lower.CL upper.CL t.ratio p.value
##  Face / Snake   Original 0.763 0.0398 73.5    0.662    0.879 -5.183  <.0001 
##  Face / Flower  Original 0.612 0.0340 82.1    0.527    0.711 -8.854  <.0001 
##  Snake / Flower Original 0.802 0.0401 65.3    0.700    0.919 -4.415  0.0001 
##  Face / Snake   SHINEd   0.736 0.0255 74.2    0.670    0.808 -8.854  <.0001 
##  Face / Flower  SHINEd   0.714 0.0281 71.8    0.642    0.794 -8.562  <.0001 
##  Snake / Flower SHINEd   0.970 0.0302 68.9    0.892    1.056 -0.974  0.3337 
## 
## Degrees-of-freedom method: satterthwaite 
## Confidence level used: 0.95 
## Conf-level adjustment: bonferroni method for 6 estimates 
## Intervals are back-transformed from the log scale 
## P value adjustment: holm method for 6 tests 
## Tests are performed on the log scale

Low-level visual attributes를 통제하기 전에는 세 범주들 간의 차이가 모두 유의하였다. 통제후에는 뱀과 꽃의 차이가 사라졌다.



3.3 Categorization Accuracy

표적 범주 보고의 정확율을 분석하였다.

EE %>% group_by(VisAttr, Category, Subject) %>% # subject-level
  summarise(mn = mean(CatCorr)) %>% 
  ungroup() %>% 
  group_by(VisAttr, Category) %>% # group-level
  summarise(Mean = mean(mn),
            Std = sd(mn)) %>% 
  ungroup() %>% 
  knitr::kable()
VisAttr Category Mean Std
Original Face 0.9475000 0.0400201
Original Snake 0.9245833 0.0762074
Original Flower 0.8700000 0.0993914
SHINEd Face 0.9391667 0.1372532
SHINEd Snake 0.9016667 0.1882437
SHINEd Flower 0.8733333 0.1959868 
cl <- makeCluster(rep("localhost", nc))
EEcat.glmer <- afex::mixed(CatCorr ~ Category*VisAttr + (Category|Subject) + (VisAttr|Image), 
                           data = EE, method = "LRT", cl = cl,
                           family = binomial(link="logit"), 
                           control = glmerControl(optimizer = "bobyqa",
                                                  optCtrl = list(maxfun = 1e6)))
## Fitting 4 (g)lmer() models.
stopCluster(cl)
anova(EEcat.glmer)
## Mixed Model Anova Table (Type 3 tests, LRT-method)
## 
## Model: CatCorr ~ Category * VisAttr + (Category | Subject) + (VisAttr | 
## Model:     Image)
## Data: EE
## Df full model: 15
##                  Df   Chisq Chi Df Pr(>Chisq)    
## Category         13 21.0362      2  2.704e-05 ***
## VisAttr          14  0.8436      1     0.3584    
## Category:VisAttr 13  0.3621      2     0.8344    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(EEcat.glmer)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: CatCorr ~ Category * VisAttr + (Category | Subject) + (VisAttr |  
##     Image)
##    Data: data
## Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+06))
## 
##      AIC      BIC   logLik deviance df.resid 
##   6879.2   6992.8  -3424.6   6849.2    14385 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -10.1036   0.1309   0.1941   0.2838   2.6724 
## 
## Random effects:
##  Groups  Name        Variance Std.Dev. Corr       
##  Image   (Intercept) 0.08643  0.2940              
##          VisAttr1    0.06150  0.2480   0.07       
##  Subject (Intercept) 1.42733  1.1947              
##          Category1   0.30445  0.5518   -0.48      
##          Category2   0.15552  0.3944    0.57 -0.76
## Number of obs: 14400, groups:  Image, 60; Subject, 40
## 
## Fixed effects:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         3.01180    0.19940  15.105  < 2e-16 ***
## Category1           0.37601    0.12455   3.019  0.00254 ** 
## Category2           0.14879    0.10853   1.371  0.17039    
## VisAttr1           -0.18193    0.19755  -0.921  0.35708    
## Category1:VisAttr1 -0.01719    0.11740  -0.146  0.88356    
## Category2:VisAttr1  0.05275    0.09715   0.543  0.58714    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Ctgry1 Ctgry2 VsAtt1 C1:VA1
## Category1   -0.308                            
## Category2    0.334 -0.648                     
## VisAttr1    -0.016 -0.003 -0.001              
## Ctgry1:VsA1 -0.003 -0.068  0.043 -0.330       
## Ctgry2:VsA1 -0.003  0.045 -0.069  0.363 -0.661
EEcat.glmer.emm <- emmeans(EEcat.glmer, "Category", by = "VisAttr", 
                           lmerTest.limit = 3480, type = "response")
plot(EEcat.glmer.emm, horizontal = FALSE) + facet_grid(cols = vars(VisAttr))

EEcat.glmer.emm.pair <- update(pairs(EEcat.glmer.emm), by = NULL, adjust = "holm")
EEcat.glmer.emm.pair %>% summary(infer = TRUE)
##  contrast       VisAttr  odds.ratio    SE  df asymp.LCL asymp.UCL z.ratio
##  Face / Snake   Original       1.17 0.326 Inf     0.561      2.44 0.565  
##  Face / Flower  Original       2.51 0.654 Inf     1.259      4.99 3.521  
##  Snake / Flower Original       2.14 0.457 Inf     1.221      3.76 3.573  
##  Face / Snake   SHINEd         1.35 0.401 Inf     0.613      2.95 0.998  
##  Face / Flower  SHINEd         2.42 0.676 Inf     1.155      5.05 3.155  
##  Snake / Flower SHINEd         1.80 0.411 Inf     0.981      3.29 2.556  
##  p.value
##  0.6369 
##  0.0021 
##  0.0021 
##  0.6369 
##  0.0064 
##  0.0318 
## 
## Confidence level used: 0.95 
## Conf-level adjustment: bonferroni method for 6 estimates 
## Intervals are back-transformed from the log odds ratio scale 
## P value adjustment: holm method for 6 tests 
## Tests are performed on the log odds ratio scale

꽃의 범주 판단이 visual attributes 통제 후 조금 정확해졌지만, 주목할 만한 변화는 아니었다.





4 Plots

지금까지의 결과를 그림으로 정리하였다.

pd <- position_dodge(.3)

# Valence
SSv.P <- as_tibble(SSv.lmer.emm) %>%
  ggplot(aes(x=Category, y=emmean, ymin=lower.CL, ymax=upper.CL, 
             color=VisAttr, shape=VisAttr)) +
  geom_point(size = 4, position = pd) +
  geom_errorbar(width = .2, position = pd) +
  geom_line(aes(group = VisAttr), position = pd) + 
  scale_color_manual(values = c("red", "black")) +
  coord_cartesian(ylim = c(1,9), clip = "on") +
  scale_y_continuous(breaks = seq(1,9)) +
  labs(x = "Category", y = "Valence") +
  theme_bw(base_size = 18) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        aspect.ratio = 1,
        legend.position = c(0.35, 0.85),
        legend.title = element_blank(),
        legend.background = element_blank(),
        legend.key = element_blank())

# Arousal
SSa.P <- as_tibble(SSa.lmer.emm) %>%
  ggplot(aes(x=Category, y=emmean, ymin=lower.CL, ymax=upper.CL, 
             color=VisAttr, shape=VisAttr)) +
  geom_point(size = 4, position = pd) +
  geom_errorbar(width = .2, position = pd) +
  geom_line(aes(group = VisAttr), position = pd) + 
  scale_color_manual(values = c("red", "black")) +
  coord_cartesian(ylim = c(1,9), clip = "on") +
  scale_y_continuous(breaks = seq(1,9)) +
  labs(x = "Category", y = "Arousal") +
  theme_bw(base_size = 18) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        aspect.ratio = 1,
        legend.position = c(0.6, 0.2),
        legend.title = element_blank(),
        legend.background = element_blank(),
        legend.key = element_blank())

plot_grid(SSv.P, SSa.P, ncol = 2, labels = c('E', 'F'), label_size = 20)

# Localization Accuracy
EEloc.P <- as_tibble(EEloc.glmer.emm) %>%
  ggplot(aes(x=Category, y=prob, ymin=asymp.LCL, ymax=asymp.UCL, 
             color=VisAttr, shape=VisAttr)) +
  geom_point(size = 4, position = pd) +
  geom_errorbar(width = .2, position = pd) +
  geom_line(aes(group = VisAttr), position = pd) + 
  scale_color_manual(values = c("red", "black")) +
  coord_cartesian(ylim = c(0.5,1), clip = "on") +
  scale_y_continuous(breaks = seq(0.5,1,0.1)) +
  labs(x = "Category", y = "Localization Accuracy") +
  theme_bw(base_size = 18) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        aspect.ratio = 1,
        legend.position = c(0.65, 0.4),
        legend.title = element_blank(),
        legend.background = element_blank(),
        legend.key = element_blank())

# Categorization Accuracy
EEcat.P <- as_tibble(EEcat.glmer.emm) %>%
  ggplot(aes(x=Category, y=prob, ymin=asymp.LCL, ymax=asymp.UCL, 
             color=VisAttr, shape=VisAttr)) +
  geom_point(size = 4, position = pd) +
  geom_errorbar(width = .2, position = pd) +
  geom_line(aes(group = VisAttr), position = pd) + 
  scale_color_manual(values = c("red", "black")) +
  coord_cartesian(ylim = c(0.3,1), clip = "on") +
  scale_y_continuous(breaks = seq(0.3,1,0.1)) +
  labs(x = "Category", y = "Categorization Accuracy") +
  theme_bw(base_size = 18) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        aspect.ratio = 1,
        legend.position = c(0.65, 0.4),
        legend.title = element_blank(),
        legend.background = element_blank(),
        legend.key = element_blank())

plot_grid(EEloc.P, EEcat.P, ncol = 2, labels = c('G', 'H'), label_size = 20)

# Localization RT = Breaking Time
EErt.P <- as_tibble(EErt.log.emm) %>%
  mutate(BT = response/1000,
         LCL = lower.CL/1000,
         UCL = upper.CL/1000) %>%
  ggplot(aes(x=Category, y=BT, ymin=LCL, ymax=UCL, 
             color=VisAttr, shape=VisAttr)) +
  geom_point(size = 4, position = pd) +
  geom_errorbar(width = .2, position = pd) +
  geom_line(aes(group = VisAttr), position = pd) + 
  scale_color_manual(values = c("red", "black")) +
  coord_cartesian(ylim = c(0,3), clip = "on") +
  scale_y_continuous(breaks = seq(0,3,0.5)) +
  labs(x = "Category", y = "Breaking Time (s)") +
  theme_bw(base_size = 18) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        aspect.ratio = 1,
        legend.position = c(0.6, 0.2),
        legend.title = element_blank(),
        legend.background = element_blank(),
        legend.key = element_blank())

plot_grid(EErt.P, ncol = 2, align = "hv", axis = "b",
          labels = c('I'), label_size = 20)





5 BLUPS: Arousal and Breaking Time

Best Linear Unbiased Predictions (BLUPs)은 선형혼합모형에서 추정한 무선 효과를 말한다. BLUPs을 예를 들어 설명하자면 다음과 같다. 만약 참가자들이 꽃보다 뱀을 봤을 때 더 각성되는지 알고 싶다면, 각성수준의 평균 변화량은 고정효과(Category)의 slope로 추정하고 개인차는 무선 기울기(random slope; (Category|Subject))로 추정할 수 있다. 꽃보다 뱀을 보고 많이 각성된 참가자의 기울기는 절대값이 큰 양의 값으로 추정될 것이고, 뱀보다 꽃에 더 각성되는 참가자의 기울기는 음의 값으로 추정될 것이다. 여기서 무선 기울기는 개별 참가자의 특성을 반영하는 추정치라고 볼 수 있다. 참가자가 20명이라면, 선형혼합모형에서 이러한 기울기 값을 20개 얻을 수 있다. 이를 BLUPs이라고 한다. 만약 우리가 참가자들의 fMRI 자료를 가지로 있다면, BLUPs을 regressor로 활용할 수 있다(사례). 본 분석에서는 참가자들이 아니라, 자극의 BLUPs을 활용할 것이다.

가설은 다음과 같다: 각성을 크게 유발한 뱀 사진일수록 연속점멸 억제를 빨리 벗어날 것이다. 연구 1의 각성 점수에서 뱀 자극은 얼굴과 꽃에 비해 저수준 시각특질 통제 여부에 따라 극적인 변화를 보였다(그림 B). 뱀의 각성점수는 평정 전에는 얼굴 수준으로 높았지만 평정 후에는 꽃 수준으로 낮아졌다. 연구 2에서도 뱀 자극은 시각특질 통제의 영향을 두드러지게 받았다(그림 E). 뱀과 꽃에 대한 반응시간 차이가 통제 후에 사라졌다. 이에 비해, 정서가(valence) 점수와 얼굴 및 꽃 자극은 시각특질 통제의 영향을 덜 받았으므로 위 가설을 검증할 조건으로 맞지 않다.

연구 1과 연구 2의 참가자는 달랐지만 자극은 같았기 때문에 자극의 BLUPs을 활용하여 위 가설을 검증할 수 있다. (1) 연구 1의 각성 점수를 위한 선형혼합모형에서 (VisAttr|Image)으로 추정한 BLUPs은 자극특질 통제 전후로 각성 점수가 변화한 정도를 의미한다. 어떤 자극의 BLUP이 음수라면, 통제 후 그 자극의 각성점수가 작아졌다는 것이다. (2) 연구 2의 breaking time을 위한 선형혼합모형에서 (VisAttr|Image)으로 추정한 BLUPs은 자극특질 통제 전후로 탐지 반응시간이 변화한 정도를 의미한다. 어떤 자극의 BLUP이 양수라면, 통제 후 그 자극이 CFS에서 더 늦게 탐지되었다는 것이다. 위 가설에 따르면, (1)과 (2)의 BLUPs은 역상관을 보일 것이다(단측 가설).

# delta Arousal
AA <- ranef(SSa.lmer$full_model)$Image
colnames(AA) <- c("intercept", "Arousal")
AA <- AA %>% rownames_to_column() %>% 
  filter(str_detect(rowname, 'snake')) %>% 
  select(rowname, Arousal) %>% 
  column_to_rownames(var = "rowname")

# deltal Valence
VV <- ranef(SSv.lmer$full_model)$Image
colnames(VV) <- c("intercept", "Valence")
VV <- VV %>% rownames_to_column() %>% 
  filter(str_detect(rowname, 'snake')) %>% 
  select(rowname, Valence) %>% 
  column_to_rownames(var = "rowname")

# delta BT 
RR <- ranef(EErt.log$full_model)$Image
colnames(RR) <- c("intercept", "BreakingTime")
RR <- RR %>% rownames_to_column() %>% 
  filter(str_detect(rowname, 'snake')) %>% 
  select(rowname, BreakingTime) %>% 
  column_to_rownames(var = "rowname")

BLUPS <- cbind(AA, VV, RR)
cor.test(formula = ~ Arousal + Valence, data = BLUPS,
         method = "pearson", alternative = "less")
## 
##  Pearson's product-moment correlation
## 
## data:  Arousal and Valence
## t = 0.41589, df = 18, p-value = 0.6588
## alternative hypothesis: true correlation is less than 0
## 95 percent confidence interval:
##  -1.0000000  0.4596011
## sample estimates:
##        cor 
## 0.09755861
cor.test(formula = ~ Arousal + BreakingTime, data = BLUPS,
         method = "pearson", alternative = "less")
## 
##  Pearson's product-moment correlation
## 
## data:  Arousal and BreakingTime
## t = -2.009, df = 18, p-value = 0.02989
## alternative hypothesis: true correlation is less than 0
## 95 percent confidence interval:
##  -1.00000000 -0.05839684
## sample estimates:
##        cor 
## -0.4279619
B.P <- ggplot(BLUPS, aes(x=Arousal, y=Valence)) +
  geom_point(size = 4) +
  geom_smooth(method=lm) +
  labs(x = expression(Delta*" Arousal"), 
       y = expression(Delta*" Valence")) +
  theme_bw(base_size = 18) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        aspect.ratio = 1)

R.P <- ggplot(BLUPS, aes(x=Arousal, y=BreakingTime)) +
  geom_point(size = 4) +
  geom_smooth(method=lm) +
  labs(x = expression(Delta*" Arousal"), 
       y = expression(Delta*" Breaking Time")) +
  theme_bw(base_size = 18) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        aspect.ratio = 1)

plot_grid(B.P, R.P, ncol = 2, align = "hv", axis = "b",
          labels = c('J', 'K'), label_size = 20)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

분석 결과는 가설에 부합하였다. 저수준 시각특질을 통제한 후에 각성점수가 낮아진 뱀 사진일수록 CFS를 벗어나는 데 걸린 시간도 증가하였다(그림 F). 이 결과는 (1) 각성수준이 CFS 지속시간을 결정한다는 점, (2) CFS 지속시간을 결정하는 요인에 의식적인 접근(평정)이 가능하다는 점, (2) 저수준 시각특질이 뱀 탐지에 영향을 끼친다는 점을 시사한다.

추가 질문은 다음과 같다. (1) Valence는 어떤가? 상관이 유의하지 않았다. (2) 이 결과를 다른 자극 범주에도 확장할 수 있는가? 뱀과 얼굴을 합쳐서 분석하면 상관이 사라졌다. (3) 얼굴은 저수준 시각특질 통제의 영향을 별로 받지 않는데, 그럼 얼굴은 무엇인가?

한편, 시각특질 통제 전후의 반응시간을 비교했을 때, 뱀 자극의 반응시간 변화량이 통계적으로 유의하지 않았다.

EErt.log.emm2 <- emmeans(EErt.log, "VisAttr", by = "Category", 
                         lmerTest.limit = 13598, type = "response")
plot(EErt.log.emm2, horizontal = FALSE, type) + facet_grid(cols = vars(Category))

EErt.log.emm2.pair <- update(pairs(EErt.log.emm2), by = NULL, adjust = "holm")
EErt.log.emm2.pair %>% summary(infer = TRUE) 
##  contrast          Category ratio     SE   df lower.CL upper.CL t.ratio p.value
##  Original / SHINEd Face      0.86 0.0688 50.3    0.705     1.05 -1.885  0.1595 
##  Original / SHINEd Snake     0.83 0.0781 46.6    0.656     1.05 -1.984  0.1595 
##  Original / SHINEd Flower    1.00 0.0974 46.0    0.788     1.28  0.036  0.9716 
## 
## Degrees-of-freedom method: satterthwaite 
## Confidence level used: 0.95 
## Conf-level adjustment: bonferroni method for 3 estimates 
## Intervals are back-transformed from the log scale 
## P value adjustment: holm method for 3 tests 
## Tests are performed on the log scale





6 Session Info

sessionInfo()
## R version 4.0.2 (2020-06-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] icon_0.1.0        klippy_0.0.0.9500 scales_1.1.1      cowplot_1.0.0    
##  [5] emmeans_1.5.0     afex_0.27-2       lme4_1.1-23       Matrix_1.2-18    
##  [9] forcats_0.5.0     dplyr_1.0.2       purrr_0.3.4       readr_1.3.1      
## [13] tidyr_1.1.2       tibble_3.0.3      ggplot2_3.3.2     tidyverse_1.3.0  
## [17] stringr_1.4.0     pacman_0.5.1     
## 
## loaded via a namespace (and not attached):
##  [1] TH.data_1.0-10      minqa_1.2.4         colorspace_1.4-1   
##  [4] ellipsis_0.3.1      rio_0.5.16          estimability_1.3   
##  [7] fs_1.5.0            rstudioapi_0.11     farver_2.0.3       
## [10] fansi_0.4.1         mvtnorm_1.1-1       lubridate_1.7.9    
## [13] xml2_1.3.2          codetools_0.2-16    splines_4.0.2      
## [16] knitr_1.29          jsonlite_1.7.1      nloptr_1.2.2.2     
## [19] broom_0.7.0.9001    dbplyr_1.4.4        compiler_4.0.2     
## [22] httr_1.4.2          backports_1.1.9     assertthat_0.2.1   
## [25] cli_2.0.2           htmltools_0.5.0     tools_4.0.2        
## [28] lmerTest_3.1-2      coda_0.19-3         gtable_0.3.0       
## [31] glue_1.4.2          reshape2_1.4.4      Rcpp_1.0.5.2       
## [34] carData_3.0-4       cellranger_1.1.0    vctrs_0.3.4        
## [37] nlme_3.1-149        xfun_0.16           openxlsx_4.1.5     
## [40] rvest_0.3.6         lifecycle_0.2.0     statmod_1.4.34     
## [43] MASS_7.3-52         zoo_1.8-8           hms_0.5.3          
## [46] sandwich_2.5-1      yaml_2.2.1          curl_4.3           
## [49] stringi_1.4.6       highr_0.8           boot_1.3-25        
## [52] zip_2.1.1           rlang_0.4.7         pkgconfig_2.0.3    
## [55] evaluate_0.14       lattice_0.20-41     labeling_0.3       
## [58] tidyselect_1.1.0    plyr_1.8.6          magrittr_1.5       
## [61] R6_2.4.1            generics_0.0.2      multcomp_1.4-13    
## [64] DBI_1.1.0           pillar_1.4.6        haven_2.3.1        
## [67] foreign_0.8-80      withr_2.2.0         mgcv_1.8-33        
## [70] survival_3.2-3      abind_1.4-5         modelr_0.1.8       
## [73] crayon_1.3.4        car_3.0-9           rmarkdown_2.3      
## [76] grid_4.0.2          readxl_1.3.1        data.table_1.13.0  
## [79] blob_1.2.1          reprex_0.3.0        digest_0.6.25      
## [82] xtable_1.8-4        numDeriv_2016.8-1.1 munsell_0.5.0

View on Github