---
title: "Analysis Script of 'Simultaneously Presented Facial and Contextual Information Influences Observers’ Facial Expressions and Self-Reports of Emotion' paper"
output: html_document
---

RStudio version the script was ran with: 4.0.2  
  
```{r include=FALSE}
#Load in the packages
library(tidyverse) #version 1.3.0
library(lme4) #version 1.1-23
```

```{r}
options(decimals = 9, scipen = 9999)
```

In this script we are using generalised linear mixed-effects models (GLMM) to analyse participants facial muscle movements (i.e. corrugator, zygomatic and orbicularis oculi) and their behavioural ratings of affect (Valence and Arousal ratings using the Self-Assessment Manikin, Bradley & Lang, 1994). 

We will use the maximal random effect structure consistent with our experimental design following published recommendations (Barr, 2013; Barr et al., 2013). That is, **Subject** and **Item** (i.e. random intercepts), the **Condition x Image Category** (i.e. the by-subject random slope) and **Condition** (the by-item random slope) will be included as random main effect terms. For all measurements, the models will include main effects of **Condition** (face-context vs context-only) and **Image Category** (neutral, mutilation, threat, happy, erotic) and all the interaction combinations between these two factors. To test for main effects, we compared models with and withoutu the fixed effects of interest using the anova(test='chisq') function. To test for simple effects of interactions when the interactiosn showed statitically significant effects, we dummy coded our variables and looked at the output of the variable of interest within the GLMMs. Because corrugator, zygomatic and orbicularis (i.e. startle eyeblink) responses were measured as continuous values ranging from negative to positive, we modelled the GLMMs with Gaussian distributions. As behavioural ratings consisted of interval ratings, we modelled this data with Poisson distributions. 

## Main data analysis:

##############################################################################################################
#########################################Corrugator data analysis#############################################

###Corrugator main effects:

```{r include=FALSE}
#Read in the corrugator data
corrugator_data <- read_csv('Data.csv') %>%
  filter(Measure == 'corrugator_EMG', Outlier == 'No') #only selecting the data without the outliers i.e. Outlier == 'No'

```

```{r}
#Deviation code the image categories
corrugator_data$erotic <- scale(ifelse(corrugator_data$Image_Category == 'Erotic', 1,0),
                             center = TRUE, scale = FALSE)
corrugator_data$mutilation <- scale(ifelse(corrugator_data$Image_Category == 'Mutilation', 1,0),
                             center = TRUE, scale = FALSE)
corrugator_data$threat <- scale(ifelse(corrugator_data$Image_Category == 'Threat', 1,0),
                             center = TRUE, scale = FALSE)
corrugator_data$happy <- scale(ifelse(corrugator_data$Image_Category == 'Happy', 1,0),
                             center = TRUE, scale = FALSE)

#Deviation code the condition
corrugator_data$face_context <- scale(ifelse(corrugator_data$Condition == 'Face-context', 1,0), 
                             center = T, scale = F) 

#Declare Participant_Number and Image as factors
corrugator_data$participant <- factor(corrugator_data$Participant_Number)
corrugator_data$item <- factor(corrugator_data$Image_Set)

#Convert corrugator EMG responses to Z-scores
corrugator_data <- corrugator_data %>%
  mutate(response = scale(Response))
```

####Run the models:

```{r}
#Maximal model
corrugator.model <- lmer(response ~ face_context + mutilation + threat + erotic + happy + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy + 
                            (1 + face_context + mutilation + threat + erotic + happy + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy|participant) + 
                            (1 + face_context|item), data = corrugator_data, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 1000000)))

```

```{r}
#Effect of condition
corrugator.no_condition <- lmer(response ~ mutilation + threat + erotic + happy + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy + 
                            (1 + face_context + mutilation + threat + erotic + happy + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy|participant) + 
                            (1 + face_context|item), data = corrugator_data, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 1000000)))

corrugator_condition_effect <- anova(corrugator.model, corrugator.no_condition, type = 'chisq')
```

```{r}
#Condition results
corrugator_condition_effect
```

```{r}
#Effect of category
corrugator.no_category <- lmer(response ~ face_context + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy + 
                            (1 + face_context + mutilation + threat + erotic + happy + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy|participant) + 
                            (1 + face_context|item), data = corrugator_data, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 1000000)))

corrugator_category_effect <- anova(corrugator.model, corrugator.no_category, type = 'chisq')
```

```{r}
#Category results
corrugator_category_effect
```

```{r}
#Effect of category-condition interaction
corrugator.no_category_condition_interaction <- lmer(response ~ face_context + mutilation + threat + erotic + happy + (1 + face_context + mutilation + threat + erotic + happy + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy|participant) + 
                            (1 + face_context|item), data = corrugator_data, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 1000000)))


corrugator_category_condition_effect <- anova(corrugator.model, corrugator.no_category_condition_interaction, type = 'chisq')
```

```{r}
#Category-condition interaction results
corrugator_category_condition_effect
```

###Corrugator simple effects:

Since the **Category** main effect was significant, we are now looking at the simple effects.

####Run the models:

#####Threat category
```{r}
corrugator_data_threat_category <- corrugator_data

#Declare Participant_Number and Image as factors
corrugator_data_threat_category$participant <- factor(corrugator_data_threat_category$Participant_Number)
corrugator_data_threat_category$item <- factor(corrugator_data_threat_category$Image_Set)

# corrugator_data_threat$threat <- ifelse(corrugator_data_threat_category$Image_Category == 'Threat', 1,0) #threat is coded as the intercept
corrugator_data_threat_category$neutral <- ifelse(corrugator_data_threat_category$Image_Category == 'Neutral', 1,0)
corrugator_data_threat_category$happy <- ifelse(corrugator_data_threat_category$Image_Category == 'Happy', 1,0)
corrugator_data_threat_category$mutilation <- ifelse(corrugator_data_threat_category$Image_Category == 'Mutilation', 1,0)
corrugator_data_threat_category$erotic <- ifelse(corrugator_data_threat_category$Image_Category == 'Erotic', 1,0)

#Convert corrugator EMG responses to Z-scores
corrugator_data_threat_category <- corrugator_data_threat_category %>%
  mutate(response = scale(Response))
```

```{r}

corrugator_data_threat_category.model <- lmer(response ~ mutilation + erotic + happy + neutral + 
                            (1 + mutilation + erotic + happy + neutral|participant) + 
                            (1|item), data = corrugator_data_threat_category, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

corrugator_coefs_threat_category <- data.frame(coef(summary(corrugator_data_threat_category.model)))

corrugator_coefs_threat_category$p.z <- 2 * (1 - pnorm(abs(corrugator_coefs_threat_category$t.value)))
```

```{r}
#Threat results
corrugator_coefs_threat_category
```

#####Mutilation
```{r}
corrugator_data_mutilation_category <- corrugator_data

#Declare Participant_Number and Image as factors
corrugator_data_mutilation_category$participant <- factor(corrugator_data_mutilation_category$Participant_Number)
corrugator_data_mutilation_category$item <- factor(corrugator_data_mutilation_category$Image_Set)

corrugator_data_mutilation_category$threat <- ifelse(corrugator_data_mutilation_category$Image_Category == 'Threat', 1,0) 
corrugator_data_mutilation_category$neutral <- ifelse(corrugator_data_mutilation_category$Image_Category == 'Neutral', 1,0)
corrugator_data_mutilation_category$happy <- ifelse(corrugator_data_mutilation_category$Image_Category == 'Happy', 1,0)
# corrugator_data_mutilation_category$mutilation <- ifelse(corrugator_data_mutilation_category$Image_Category == 'Mutilation', 1,0) #mutilation is coded as the intercept
corrugator_data_mutilation_category$erotic <- ifelse(corrugator_data_mutilation_category$Image_Category == 'Erotic', 1,0)

#Convert corrugator EMG responses to Z-scores
corrugator_data_mutilation_category <- corrugator_data_mutilation_category %>%
  mutate(response = scale(Response))
```

```{r}

corrugator_data_mutilation_category.model <- lmer(response ~ threat + erotic + happy + neutral + 
                            (1 + threat + erotic + happy + neutral|participant) + 
                            (1|item), data = corrugator_data_mutilation_category, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

corrugator_coefs_mutilation_category <- data.frame(coef(summary(corrugator_data_mutilation_category.model)))

corrugator_coefs_mutilation_category$p.z <- 2 * (1 - pnorm(abs(corrugator_coefs_mutilation_category$t.value)))
```

```{r}
#Mutilation results
corrugator_coefs_mutilation_category
```

#####Happy
```{r}
corrugator_data_happy_category <- corrugator_data

#Declare Participant_Number and Image as factors
corrugator_data_happy_category$participant <- factor(corrugator_data_happy_category$Participant_Number)
corrugator_data_happy_category$item <- factor(corrugator_data_happy_category$Image_Set)

corrugator_data_happy_category$threat <- ifelse(corrugator_data_happy_category$Image_Category == 'Threat', 1,0) 
corrugator_data_happy_category$neutral <- ifelse(corrugator_data_happy_category$Image_Category == 'Neutral', 1,0)
# corrugator_data_happy_category$happy <- ifelse(corrugator_data_happy_category$Image_Category == 'Happy', 1,0) #happy is coded as the intercept
corrugator_data_happy_category$mutilation <- ifelse(corrugator_data_happy_category$Image_Category == 'Mutilation', 1,0) 
corrugator_data_happy_category$erotic <- ifelse(corrugator_data_happy_category$Image_Category == 'Erotic', 1,0)

#Convert corrugator EMG responses to Z-scores
corrugator_data_happy_category <- corrugator_data_happy_category %>%
  mutate(response = scale(Response))
```

```{r}

corrugator_data_happy_category.model <- lmer(response ~ threat + erotic + mutilation + neutral + 
                            (1 + threat + erotic + mutilation + neutral|participant) + 
                            (1|item), data = corrugator_data_happy_category, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

corrugator_coefs_happy_category <- data.frame(coef(summary(corrugator_data_happy_category.model)))

corrugator_coefs_happy_category$p.z <- 2 * (1 - pnorm(abs(corrugator_coefs_happy_category$t.value)))
```

```{r}
#Happy results
corrugator_coefs_happy_category
```

#####Erotic
```{r}
corrugator_data_erotic_category <- corrugator_data

#Declare Participant_Number and Image as factors
corrugator_data_erotic_category$participant <- factor(corrugator_data_erotic_category$Participant_Number)
corrugator_data_erotic_category$item <- factor(corrugator_data_erotic_category$Image_Set)

corrugator_data_erotic_category$threat <- ifelse(corrugator_data_erotic_category$Image_Category == 'Threat', 1,0) 
corrugator_data_erotic_category$neutral <- ifelse(corrugator_data_erotic_category$Image_Category == 'Neutral', 1,0)
corrugator_data_erotic_category$happy <- ifelse(corrugator_data_erotic_category$Image_Category == 'Happy', 1,0) 
corrugator_data_erotic_category$mutilation <- ifelse(corrugator_data_erotic_category$Image_Category == 'Mutilation', 1,0) 
# corrugator_data_erotic_category$erotic <- ifelse(corrugator_data_erotic_category$imageCategory == 'Erotic', 1,0) #Erotic is coded as the intercept

#Convert corrugator EMG responses to Z-scores
corrugator_data_erotic_category <- corrugator_data_erotic_category %>%
  mutate(response = scale(Response))
```

```{r}

corrugator_data_erotic_category.model <- lmer(response ~ threat + happy + mutilation + neutral + 
                            (1 + threat + happy + mutilation + neutral|participant) + 
                            (1|item), data = corrugator_data_erotic_category, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

corrugator_coefs_erotic_category <- data.frame(coef(summary(corrugator_data_erotic_category.model)))

corrugator_coefs_erotic_category$p.z <- 2 * (1 - pnorm(abs(corrugator_coefs_erotic_category$t.value)))
```

```{r}
#Erotic results
corrugator_coefs_erotic_category
```

#####Neutral
```{r}
corrugator_data_neutral_category <- corrugator_data

#Declare Participant_Number and Image as factors
corrugator_data_neutral_category$participant <- factor(corrugator_data_neutral_category$Participant_Number)
corrugator_data_neutral_category$item <- factor(corrugator_data_neutral_category$Image_Set)

corrugator_data_neutral_category$threat <- ifelse(corrugator_data_neutral_category$Image_Category == 'Threat', 1,0) 
# corrugator_data_neutral_category$neutral <- ifelse(corrugator_data_neutral_category$Image_Category == 'Neutral', 1,0) #Neutral is coded as the intercept
corrugator_data_neutral_category$happy <- ifelse(corrugator_data_neutral_category$Image_Category == 'Happy', 1,0) 
corrugator_data_neutral_category$mutilation <- ifelse(corrugator_data_neutral_category$Image_Category == 'Mutilation', 1,0) 
corrugator_data_neutral_category$erotic <- ifelse(corrugator_data_neutral_category$Image_Category == 'Erotic', 1,0)

#Convert corrugator EMG responses to Z-scores
corrugator_data_neutral_category <- corrugator_data_neutral_category %>%
  mutate(response = scale(Response))
```

```{r}

corrugator_data_neutral_category.model <- lmer(response ~ threat + happy + mutilation + erotic + 
                            (1 + threat + happy + mutilation + erotic|participant) + 
                            (1|item), data = corrugator_data_neutral_category, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

corrugator_coefs_neutral_category <- data.frame(coef(summary(corrugator_data_neutral_category.model)))

corrugator_coefs_neutral_category$p.z <- 2 * (1 - pnorm(abs(corrugator_coefs_neutral_category$t.value)))
```

```{r}
#Neutral results
corrugator_coefs_neutral_category
```

Since the **Category x Condition** was significant, we are now looking at the simple effects.

####Run the models:

#####Threat category
```{r}
corrugator_data_threat_category_condition <- corrugator_data

#Declare Participant_Number and Image as factors
corrugator_data_threat_category_condition$participant <- factor(corrugator_data_threat_category_condition$Participant_Number)
corrugator_data_threat_category_condition$item <- factor(corrugator_data_threat_category_condition$Image_Set)

#Dummy coding for simple effects
corrugator_data_threat_category_condition$face_context <- ifelse(corrugator_data_threat_category_condition$Condition == 'Face-context', 1, 0)  

# corrugator_data_threat_category_condition$threat <- ifelse(corrugator_data_threat_category_condition$Image_Category == 'Threat', 1,0) #threat is coded as the intercept
corrugator_data_threat_category_condition$neutral <- ifelse(corrugator_data_threat_category_condition$Image_Category == 'Neutral', 1,0)
corrugator_data_threat_category_condition$happy <- ifelse(corrugator_data_threat_category_condition$Image_Category == 'Happy', 1,0)
corrugator_data_threat_category_condition$mutilation <- ifelse(corrugator_data_threat_category_condition$Image_Category == 'Mutilation', 1,0)
corrugator_data_threat_category_condition$erotic <- ifelse(corrugator_data_threat_category_condition$Image_Category == 'Erotic', 1,0)

#Z-score
corrugator_data_threat_category_condition <- corrugator_data_threat_category_condition %>%
  mutate(response = scale(Response))
```

```{r}

corrugator_data_threat_category_condition.model <- lme4::lmer(response ~ face_context + mutilation + erotic + happy + neutral + face_context:neutral + face_context:mutilation + face_context:erotic + face_context:happy + 
                            (1 + face_context + mutilation + erotic + happy + neutral + face_context:neutral + face_context:mutilation + face_context:erotic + face_context:happy|participant) + 
                            (1 + face_context|item), data = corrugator_data_threat_category_condition, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

corrugator_coefs_threat_category_condition <- data.frame(coef(summary(corrugator_data_threat_category_condition.model)))

corrugator_coefs_threat_category_condition$p.z <- 2 * (1 - pnorm(abs(corrugator_coefs_threat_category_condition$t.value)))
```

```{r}
#Theat-condition interaction results
corrugator_coefs_threat_category_condition
```
#####Mutilation
```{r}
corrugator_data_mutilation_category_condition <- corrugator_data

#Declare Participant_Number and Image as factors
corrugator_data_mutilation_category_condition$participant <- factor(corrugator_data_mutilation_category_condition$Participant_Number)
corrugator_data_mutilation_category_condition$item <- factor(corrugator_data_mutilation_category_condition$Image_Set)

#Dummy coding for simple effects
corrugator_data_mutilation_category_condition$face_context <- ifelse(corrugator_data_mutilation_category_condition$Condition == 'Face-context', 1, 0)  #intact is coded as the intercept

corrugator_data_mutilation_category_condition$threat <- ifelse(corrugator_data_mutilation_category_condition$Image_Category == 'Threat', 1,0) 
corrugator_data_mutilation_category_condition$neutral <- ifelse(corrugator_data_mutilation_category_condition$Image_Category == 'Neutral', 1,0)
corrugator_data_mutilation_category_condition$happy <- ifelse(corrugator_data_mutilation_category_condition$Image_Category == 'Happy', 1,0)
# corrugator_data_mutilation_category_condition$mutilation <- ifelse(corrugator_data_mutilation_category_condition$Image_Category == 'Mutilation', 1,0) #mutilation is coded as the intercept
corrugator_data_mutilation_category_condition$erotic <- ifelse(corrugator_data_mutilation_category_condition$Image_Category == 'Erotic', 1,0)

#Z-score
corrugator_data_mutilation_category_condition <- corrugator_data_mutilation_category_condition %>%
  mutate(response = scale(Response))
```

```{r}

corrugator_data_mutilation_category_condition.model <- lmer(response ~ face_context + threat + erotic + happy + neutral + face_context:neutral + face_context:threat + face_context:erotic + face_context:happy + 
                            (1 + face_context + threat + erotic + happy + neutral + face_context:neutral + face_context:threat + face_context:erotic + face_context:happy|participant) + 
                            (1 + face_context|item), data = corrugator_data_mutilation_category_condition, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

corrugator_coefs_mutilation_category_condition <- data.frame(coef(summary(corrugator_data_mutilation_category_condition.model)))

corrugator_coefs_mutilation_category_condition$p.z <- 2 * (1 - pnorm(abs(corrugator_coefs_mutilation_category_condition$t.value)))
```

```{r}
#Mutilation-condition interaction results
corrugator_coefs_mutilation_category_condition
```

#####Happy
```{r}
corrugator_data_happy_category_condition <- corrugator_data

#Declare Participant_Number and Image as factors
corrugator_data_happy_category_condition$participant <- factor(corrugator_data_happy_category_condition$Participant_Number)
corrugator_data_happy_category_condition$item <- factor(corrugator_data_happy_category_condition$Image_Set)

#Dummy coding for simple effects
corrugator_data_happy_category_condition$face_context <- ifelse(corrugator_data_happy_category_condition$Condition == 'Face-context', 1, 0) #Context-only is coded as the intercept

corrugator_data_happy_category_condition$threat <- ifelse(corrugator_data_happy_category_condition$Image_Category == 'Threat', 1,0) 
corrugator_data_happy_category_condition$neutral <- ifelse(corrugator_data_happy_category_condition$Image_Category == 'Neutral', 1,0)
# corrugator_data_happy_category_condition$happy <- ifelse(corrugator_data_happy_category_condition$Image_Category == 'Happy', 1,0) #happy is coded as the intercept
corrugator_data_happy_category_condition$mutilation <- ifelse(corrugator_data_happy_category_condition$Image_Category == 'Mutilation', 1,0) 
corrugator_data_happy_category_condition$erotic <- ifelse(corrugator_data_happy_category_condition$Image_Category == 'Erotic', 1,0)

#Z-score
corrugator_data_happy_category_condition <- corrugator_data_happy_category_condition %>%
  mutate(response = scale(Response))
```

```{r}

corrugator_data_happy_category_condition.model <- lmer(response ~ face_context + threat + erotic + mutilation + neutral + face_context:neutral + face_context:threat + face_context:erotic + 
                                       face_context:mutilation + 
                            (1 + face_context + threat + erotic + mutilation + neutral + face_context:neutral + face_context:threat + face_context:erotic + face_context:mutilation|participant) + 
                            (1 + face_context|item), data = corrugator_data_happy_category_condition, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 1000000)))

corrugator_coefs_happy_category_condition <- data.frame(coef(summary(corrugator_data_happy_category_condition.model)))

corrugator_coefs_happy_category_condition$p.z <- 2 * (1 - pnorm(abs(corrugator_coefs_happy_category_condition$t.value)))
```

```{r}
#Happy-condition interaction results
corrugator_coefs_happy_category_condition
```

#####Erotic
```{r}
corrugator_data_erotic_category_condition <- corrugator_data

#Declare Participant_Number and Image as factors
corrugator_data_erotic_category_condition$participant <- factor(corrugator_data_erotic_category_condition$Participant_Number)
corrugator_data_erotic_category_condition$item <- factor(corrugator_data_erotic_category_condition$Image_Set)

#Dummy coding for simple effects
corrugator_data_erotic_category_condition$face_context <- ifelse(corrugator_data_erotic_category_condition$Condition == 'Face-context', 1, 0)  #Context-only is coded as the intercept

corrugator_data_erotic_category_condition$threat <- ifelse(corrugator_data_erotic_category_condition$Image_Category == 'Threat', 1,0) 
corrugator_data_erotic_category_condition$neutral <- ifelse(corrugator_data_erotic_category_condition$Image_Category == 'Neutral', 1,0)
corrugator_data_erotic_category_condition$happy <- ifelse(corrugator_data_erotic_category_condition$Image_Category == 'Happy', 1,0) 
corrugator_data_erotic_category_condition$mutilation <- ifelse(corrugator_data_erotic_category_condition$Image_Category == 'Mutilation', 1,0) 
# corrugator_data_erotic_category_condition$erotic <- ifelse(corrugator_data_erotic_category_condition$Image_Category == 'Erotic', 1,0) #Erotic is coded as the intercept

#Z-score
corrugator_data_erotic_category_condition <- corrugator_data_erotic_category_condition %>%
  mutate(response = scale(Response))
```

```{r}

corrugator_data_erotic_category_condition.model <- lmer(response ~ face_context + threat + happy + mutilation + neutral + face_context:neutral + face_context:threat + face_context:happy + 
                                       face_context:mutilation + 
                            (1 + face_context + threat + happy + mutilation + neutral + face_context:neutral + face_context:threat + face_context:happy + face_context:mutilation|participant) + 
                            (1 + face_context|item), data = corrugator_data_erotic_category_condition, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 1000000)))

corrugator_coefs_erotic_category_condition <- data.frame(coef(summary(corrugator_data_erotic_category_condition.model)))

corrugator_coefs_erotic_category_condition$p.z <- 2 * (1 - pnorm(abs(corrugator_coefs_erotic_category_condition$t.value)))
```

```{r}
#Erotic-condition interaction results
corrugator_coefs_erotic_category_condition
```

#####Neutral
```{r}
corrugator_data_neutral_category_condition <- corrugator_data

#Declare Participant_Number and Image as factors
corrugator_data_neutral_category_condition$participant <- factor(corrugator_data_neutral_category_condition$Participant_Number)
corrugator_data_neutral_category_condition$item <- factor(corrugator_data_neutral_category_condition$Image_Set)

#Dummy coding for simple effects
corrugator_data_neutral_category_condition$face_context <- ifelse(corrugator_data_neutral_category_condition$Condition == 'Face-context', 1, 0)  #Context-only is coded as the intercept

corrugator_data_neutral_category_condition$threat <- ifelse(corrugator_data_neutral_category_condition$Image_Category == 'Threat', 1,0) 
# corrugator_data_neutral_category_condition$neutral <- ifelse(corrugator_data_neutral_category_condition$Image_Category == 'Neutral', 1,0) #Neutral is coded as the intercept
corrugator_data_neutral_category_condition$happy <- ifelse(corrugator_data_neutral_category_condition$Image_Category == 'Happy', 1,0) 
corrugator_data_neutral_category_condition$mutilation <- ifelse(corrugator_data_neutral_category_condition$Image_Category == 'Mutilation', 1,0) 
corrugator_data_neutral_category_condition$erotic <- ifelse(corrugator_data_neutral_category_condition$Image_Category == 'Erotic', 1,0)

#Z-score
corrugator_data_neutral_category_condition <- corrugator_data_neutral_category_condition %>%
  mutate(response = scale(Response))
```

```{r}

corrugator_data_neutral_category_condition.model <- lmer(response ~ face_context + threat + happy + mutilation + erotic + face_context:erotic + face_context:threat + face_context:happy + 
                                       face_context:mutilation + 
                            (1 + face_context + threat + happy + mutilation + erotic + face_context:erotic + face_context:threat + face_context:happy + face_context:mutilation|participant) + 
                            (1 + face_context|item), data = corrugator_data_neutral_category_condition, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

corrugator_coefs_neutral_category_condition <- data.frame(coef(summary(corrugator_data_neutral_category_condition.model)))

corrugator_coefs_neutral_category_condition$p.z <- 2 * (1 - pnorm(abs(corrugator_coefs_neutral_category_condition$t.value)))
```

```{r}
#Neutral-condition interaction results
corrugator_coefs_neutral_category_condition
```

##############################################################################################################
#########################################Zygomatic data analysis##############################################

###Zygomatic main effects:
```{r message=FALSE, warning=FALSE, include=FALSE}
#Read in the zygomatic data
zygomatic_data <- read_csv('Data.csv') %>%
  filter(Measure == 'zygomatic_EMG', Outlier == 'No')
```

```{r}
#Deviation code the image categories
zygomatic_data$erotic <- scale(ifelse(zygomatic_data$Image_Category == 'Erotic', 1,0),
                             center = TRUE, scale = FALSE)
zygomatic_data$mutilation <- scale(ifelse(zygomatic_data$Image_Category == 'Mutilation', 1,0),
                             center = TRUE, scale = FALSE)
zygomatic_data$threat <- scale(ifelse(zygomatic_data$Image_Category == 'Threat', 1,0),
                             center = TRUE, scale = FALSE)
zygomatic_data$happy <- scale(ifelse(zygomatic_data$Image_Category == 'Happy', 1,0),
                             center = TRUE, scale = FALSE)

#Deviation code the condition
zygomatic_data$face_context <- scale(ifelse(zygomatic_data$Condition == 'Face-context', 1,0), 
                             center = T, scale = F)

#Declare Participant_Number and Image as factors
zygomatic_data$participant <- factor(zygomatic_data$Participant_Number)
zygomatic_data$item <- factor(zygomatic_data$Image_Set)

#Convert zygomatic EMG responses to Z-scores
zygomatic_data <- zygomatic_data %>%
  mutate(response = scale(Response))

```

###Run the models:

####Main effects:
```{r}
#Maximal model
zygomatic.model <- lmer(response ~ face_context + mutilation + threat + erotic + happy + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy + 
                           (1 + face_context + mutilation + threat + erotic + happy + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy|participant) + 
                           (1 + face_context|item), data = zygomatic_data, REML = FALSE,
                         control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 1000000)))

```

```{r}
#Effect of condition
zygomatic.no_condition <- lmer(response ~ mutilation + threat + erotic + happy + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy + 
                           (1 + face_context + mutilation + threat + erotic + happy + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy|participant) + 
                           (1 + face_context|item), data = zygomatic_data, REML = FALSE,
                         control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 1000000)))

zygomatic_condition_effect <- anova(zygomatic.model, zygomatic.no_condition, type = 'chisq')
```

```{r}
#Condition results
zygomatic_condition_effect
```

```{r}
#Effect of category
zygomatic.no_category <- lmer(response ~ face_context + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy + 
                           (1 + face_context + mutilation + threat + erotic + happy + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy|participant) + 
                           (1 + face_context|item), data = zygomatic_data, REML = FALSE,
                         control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 1000000)))

zygomatic_category_effect <- anova(zygomatic.model, zygomatic.no_category, type = 'chisq')
```

```{r}
#Category results
zygomatic_category_effect
```

```{r}
#Effect of category-condition interaction
zygomatic_no_category_condition_interaction <- lmer(response ~ face_context + mutilation + threat + erotic + happy +  
                           (1 + face_context + mutilation + threat + erotic + happy + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy|participant) + 
                           (1 + face_context|item), data = zygomatic_data, REML = FALSE,
                         control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 1000000)))

zygomatic_category_condition_effect <- anova(zygomatic.model, zygomatic_no_category_condition_interaction, type = 'chisq')
```

```{r}
#Category-condition interaction results
zygomatic_category_condition_effect
```

###Zygomatic simple effects:

Since the main effect of **Category** was significant, we are now looking at the simple effects.

#####Threat category
```{r}
zygomatic_data_threat_category <- zygomatic_data

#Declare Participant_Number and Image as factors
zygomatic_data_threat_category$participant <- factor(zygomatic_data_threat_category$Participant_Number)
zygomatic_data_threat_category$item <- factor(zygomatic_data_threat_category$Image_Set)

# zygomatic_data_threat$threat <- ifelse(zygomatic_data_threat_category$Image_Category == 'Threat', 1,0) #threat is coded as the intercept
zygomatic_data_threat_category$neutral <- ifelse(zygomatic_data_threat_category$Image_Category == 'Neutral', 1,0)
zygomatic_data_threat_category$happy <- ifelse(zygomatic_data_threat_category$Image_Category == 'Happy', 1,0)
zygomatic_data_threat_category$mutilation <- ifelse(zygomatic_data_threat_category$Image_Category == 'Mutilation', 1,0)
zygomatic_data_threat_category$erotic <- ifelse(zygomatic_data_threat_category$Image_Category == 'Erotic', 1,0)

#Convert zygomatic EMG responses to Z-scores
zygomatic_data_threat_category <- zygomatic_data_threat_category %>%
  mutate(response = scale(Response))
```

```{r}

zygomatic_data_threat_category.model <- lmer(response ~ mutilation + erotic + happy + neutral + 
                                                      (1 + mutilation + erotic + happy + neutral|participant) + 
                                                      (1|item), data = zygomatic_data_threat_category, REML = FALSE,
                                                    control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

zygomatic_coefs_threat_category <- data.frame(coef(summary(zygomatic_data_threat_category.model)))

zygomatic_coefs_threat_category$p.z <- 2 * (1 - pnorm(abs(zygomatic_coefs_threat_category$t.value)))
```

```{r}
#Threat results
zygomatic_coefs_threat_category
```

#####Mutilation
```{r}
zygomatic_data_mutilation_category <- zygomatic_data

#declare subjID anf Item/Image as factors
zygomatic_data_mutilation_category$participant <- factor(zygomatic_data_mutilation_category$Participant_Number)
zygomatic_data_mutilation_category$item <- factor(zygomatic_data_mutilation_category$Image_Set)

zygomatic_data_mutilation_category$threat <- ifelse(zygomatic_data_mutilation_category$Image_Category == 'Threat', 1,0) 
zygomatic_data_mutilation_category$neutral <- ifelse(zygomatic_data_mutilation_category$Image_Category == 'Neutral', 1,0)
zygomatic_data_mutilation_category$happy <- ifelse(zygomatic_data_mutilation_category$Image_Category == 'Happy', 1,0)
# zygomatic_data_mutilation_category$mutilation <- ifelse(zygomatic_data_mutilation_category$Image_Category == 'Mutilation', 1,0) #mutilation is coded as the intercept
zygomatic_data_mutilation_category$erotic <- ifelse(zygomatic_data_mutilation_category$Image_Category == 'Erotic', 1,0)

#Convert zygomatic EMG responses to Z-scores
zygomatic_data_mutilation_category <- zygomatic_data_mutilation_category %>%
  mutate(response = scale(Response))
```

```{r}

zygomatic_data_mutilation_category.model <- lmer(response ~ threat + erotic + happy + neutral + 
                                                    (1 + threat + erotic + happy + neutral|participant) + 
                                                    (1|item), data = zygomatic_data_mutilation_category, REML = FALSE,
                                                  control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

zygomatic_coefs_mutilation_category <- data.frame(coef(summary(zygomatic_data_mutilation_category.model)))

zygomatic_coefs_mutilation_category$p.z <- 2 * (1 - pnorm(abs(zygomatic_coefs_mutilation_category$t.value)))
```

```{r}
#Mutilation results
zygomatic_coefs_mutilation_category
```

#####Happy
```{r}
zygomatic_data_happy_category <- zygomatic_data

#Declare Participant_Number and Image as factors
zygomatic_data_happy_category$participant <- factor(zygomatic_data_happy_category$Participant_Number)
zygomatic_data_happy_category$item <- factor(zygomatic_data_happy_category$Image_Set)

zygomatic_data_happy_category$threat <- ifelse(zygomatic_data_happy_category$Image_Category == 'Threat', 1,0) 
zygomatic_data_happy_category$neutral <- ifelse(zygomatic_data_happy_category$Image_Category == 'Neutral', 1,0)
# zygomatic_data_happy_category$happy <- ifelse(zygomatic_data_happy_category$Image_Category == 'Happy', 1,0) #happy is coded as the intercept
zygomatic_data_happy_category$mutilation <- ifelse(zygomatic_data_happy_category$Image_Category == 'Mutilation', 1,0) 
zygomatic_data_happy_category$erotic <- ifelse(zygomatic_data_happy_category$Image_Category == 'Erotic', 1,0)

#Convert zygomatic EMG responses to Z-scores
zygomatic_data_happy_category <- zygomatic_data_happy_category %>%
  mutate(response = scale(Response))
```

```{r}

zygomatic_data_happy_category.model <- lmer(response ~ threat + erotic + mutilation + neutral + 
                                                (1 + threat + erotic + mutilation + neutral|participant) + 
                                                (1|item), data = zygomatic_data_happy_category, REML = FALSE,
                                              control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

zygomatic_coefs_happy_category <- data.frame(coef(summary(zygomatic_data_happy_category.model)))

zygomatic_coefs_happy_category$p.z <- 2 * (1 - pnorm(abs(zygomatic_coefs_happy_category$t.value)))
```

```{r}
#Happy results
zygomatic_coefs_happy_category
```

#####Erotic
```{r}
zygomatic_data_erotic_category <- zygomatic_data

#Declare Participant_Number and Image as factors
zygomatic_data_erotic_category$participant <- factor(zygomatic_data_erotic_category$Participant_Number)
zygomatic_data_erotic_category$item <- factor(zygomatic_data_erotic_category$Image_Set)

zygomatic_data_erotic_category$threat <- ifelse(zygomatic_data_erotic_category$Image_Category == 'Threat', 1,0) 
zygomatic_data_erotic_category$neutral <- ifelse(zygomatic_data_erotic_category$Image_Category == 'Neutral', 1,0)
zygomatic_data_erotic_category$happy <- ifelse(zygomatic_data_erotic_category$Image_Category == 'Happy', 1,0) 
zygomatic_data_erotic_category$mutilation <- ifelse(zygomatic_data_erotic_category$Image_Category == 'Mutilation', 1,0) 
# zygomatic_data_erotic_category$erotic <- ifelse(zygomatic_data_erotic_category$Image_Category == 'Erotic', 1,0) #Erotic is coded as the intercept

#Convert zygomatic EMG responses to Z-scores
zygomatic_data_erotic_category <- zygomatic_data_erotic_category %>%
  mutate(response = scale(Response))
```

```{r}

zygomatic_data_erotic_category.model <- lmer(response ~ threat + happy + mutilation + neutral + 
                                                (1 + threat + happy + mutilation + neutral|participant) + 
                                                (1|item), data = zygomatic_data_erotic_category, REML = FALSE,
                                              control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

zygomatic_coefs_erotic_category <- data.frame(coef(summary(zygomatic_data_erotic_category.model)))

zygomatic_coefs_erotic_category$p.z <- 2 * (1 - pnorm(abs(zygomatic_coefs_erotic_category$t.value)))
```

```{r}
#Erotic results
zygomatic_coefs_erotic_category
```

#####Neutral
```{r}
zygomatic_data_neutral_category <- zygomatic_data

#Declare Participant_Number and Image as factors
zygomatic_data_neutral_category$participant <- factor(zygomatic_data_neutral_category$Participant_Number)
zygomatic_data_neutral_category$item <- factor(zygomatic_data_neutral_category$Image_Set)

zygomatic_data_neutral_category$threat <- ifelse(zygomatic_data_neutral_category$Image_Category == 'Threat', 1,0) 
# zygomatic_data_neutral_category$neutral <- ifelse(zygomatic_data_neutral_category$Image_Category == 'Neutral', 1,0) #Neutral is coded as the intercept
zygomatic_data_neutral_category$happy <- ifelse(zygomatic_data_neutral_category$Image_Category == 'Happy', 1,0) 
zygomatic_data_neutral_category$mutilation <- ifelse(zygomatic_data_neutral_category$Image_Category == 'Mutilation', 1,0) 
zygomatic_data_neutral_category$erotic <- ifelse(zygomatic_data_neutral_category$Image_Category == 'Erotic', 1,0)

#Convert zygomatic EMG responses to Z-scores
zygomatic_data_neutral_category <- zygomatic_data_neutral_category %>%
  mutate(response = scale(Response))
```

```{r}

zygomatic_data_neutral_category.model <- lmer(response ~ threat + happy + mutilation + erotic + 
                                                 (1 + threat + happy + mutilation + erotic|participant) + 
                                                 (1|item), data = zygomatic_data_neutral_category, REML = FALSE,
                                               control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

zygomatic_coefs_neutral_category <- data.frame(coef(summary(zygomatic_data_neutral_category.model)))

zygomatic_coefs_neutral_category$p.z <- 2 * (1 - pnorm(abs(zygomatic_coefs_neutral_category$t.value)))
```

```{r}
#Neutral results
zygomatic_coefs_neutral_category
```

Since the main effect of **Condition** was significant, we are now looking at the simple effects.

####Run the model:
```{r}
zygomatic_data_condition <- zygomatic_data

#Declare Participant_Number and Image as factors
zygomatic_data_condition$participant <- factor(zygomatic_data_condition$Participant_Number)
zygomatic_data_condition$item <- factor(zygomatic_data_condition$Image_Set)

#Dummy coding for simple effects
zygomatic_data_condition$face_context <- ifelse(zygomatic_data_condition$Condition == 'Face-context', 1, 0) 


zygomatic_data_condition <- zygomatic_data_condition %>%
  mutate(response = scale(Response))
```

```{r}
zygomatic_data_condition.model <- lmer(response ~ face_context +
                                                (1 + face_context|participant) + 
                                                (1 + face_context|item), data = zygomatic_data_condition, REML = FALSE,
                                              control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

zygomatic_coefs_condition <- data.frame(coef(summary(zygomatic_data_condition.model)))

zygomatic_coefs_condition$p.z <- 2 * (1 - pnorm(abs(zygomatic_coefs_condition$t.value)))
```

```{r}
#Condition results
zygomatic_coefs_condition
```

##############################################################################################################
###########################################Orbicularis data analysis##########################################

```{r message=FALSE, warning=FALSE, include=FALSE, paged.print=FALSE}
#Read in the orbicularis data
SBR_data <- read_csv('Data.csv') %>%
  filter(Measure == 'orbicularis_EMG', Outlier == 'No')

```

```{r}
#Deviation code the image categories
SBR_data$erotic <- scale(ifelse(SBR_data$Image_Category == 'Erotic', 1,0),
                             center = TRUE, scale = FALSE)
SBR_data$mutilation <- scale(ifelse(SBR_data$Image_Category == 'Mutilation', 1,0),
                             center = TRUE, scale = FALSE)
SBR_data$threat <- scale(ifelse(SBR_data$Image_Category == 'Threat', 1,0),
                             center = TRUE, scale = FALSE)
SBR_data$happy <- scale(ifelse(SBR_data$Image_Category == 'Happy', 1,0),
                             center = TRUE, scale = FALSE)

#Deviation code the condition
SBR_data$face_context <- scale(ifelse(SBR_data$Condition == 'Face-context', 1,0), 
                             center = T, scale = F)

#Declare Participant_Number and Image as factors
SBR_data$participant <- factor(SBR_data$Participant_Number)
SBR_data$item <- factor(SBR_data$Image_Set)

#z-score
SBR_data <- SBR_data %>%
  mutate(response = scale(Response))
```

####Run the models: 
```{r}
#Maximal model
SBR.model <- lmer(response ~ face_context + mutilation + threat + erotic + happy + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy + 
                            (1 + face_context + mutilation + threat + erotic + happy + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy|participant) + 
                            (1 + face_context|item), data = SBR_data, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))
```

```{r}
#Effect of condition
SBR.no_condition <- lmer(response ~ mutilation + threat + erotic + happy + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy + 
                            (1 + face_context + mutilation + threat + erotic + happy + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy|participant) + 
                            (1 + face_context|item), data = SBR_data, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))
```

```{r}
#Condition results
anova(SBR.model, SBR.no_condition, type = 'chisq')
```

```{r}
#Effect of category
SBR.no_category <- lmer(response ~ face_context + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy + 
                            (1 + face_context + mutilation + threat + erotic + happy + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy|participant) + 
                            (1 + face_context|item), data = SBR_data, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))
```

```{r}
#Category results
anova(SBR.model, SBR.no_category, type = 'chisq')
```

```{r}
#Effect of category-condition interaction
SBR_no_category_condition_interaction <- lmer(response ~ face_context + mutilation + threat + erotic + happy + 
                            (1 + face_context + mutilation + threat + erotic + happy + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy|participant) + 
                            (1 + face_context|item), data = SBR_data, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))
```

```{r}
#Category-condition interaction results
anova(SBR.model, SBR_no_category_condition_interaction, type = 'chisq')
```

Because the main effect of **Category** was significant, we ran the simple effects.

###Simple effects for SBR category:

#####Threat category
```{r}
SBR_data_threat <- SBR_data

#Declare Participant_Number and Image as factors
SBR_data_threat$participant <- factor(SBR_data_threat$Participant_Number)
SBR_data_threat$item <- factor(SBR_data_threat$Image_Set)

# SBR_data_threat$threat <- ifelse(SBR_data_threat$Image_Category == 'Threat', 1,0) #threat is coded as the intercept
SBR_data_threat$neutral <- ifelse(SBR_data_threat$Image_Category == 'Neutral', 1,0)
SBR_data_threat$happy <- ifelse(SBR_data_threat$Image_Category == 'Happy', 1,0)
SBR_data_threat$mutilation <- ifelse(SBR_data_threat$Image_Category == 'Mutilation', 1,0)
SBR_data_threat$erotic <- ifelse(SBR_data_threat$Image_Category == 'Erotic', 1,0)

SBR_data_threat <- SBR_data_threat %>%
  mutate(response = scale(Response))
```

```{r}
SBR_data_threat_category.model <- lmer(response ~ mutilation + erotic + happy + neutral + 
                            (1 + mutilation + erotic + happy + neutral|participant) + 
                            (1|item), data = SBR_data_threat, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

#summary(SBR_data_threat_category.model)

coefs_threat_category <- data.frame(coef(summary(SBR_data_threat_category.model)))

coefs_threat_category$p.z <- 2 * (1 - pnorm(abs(coefs_threat_category$t.value)))
```

```{r}
#Threat results
coefs_threat_category
```

#####Mutilation
```{r}
SBR_data_mutilation <- SBR_data

#Declare Participant_Number and Image as factors
SBR_data_mutilation$participant <- factor(SBR_data_mutilation$Participant_Number)
SBR_data_mutilation$item <- factor(SBR_data_mutilation$Image_Set)

SBR_data_mutilation$threat <- ifelse(SBR_data_mutilation$Image_Category == 'Threat', 1,0) 
SBR_data_mutilation$neutral <- ifelse(SBR_data_mutilation$Image_Category == 'Neutral', 1,0)
SBR_data_mutilation$happy <- ifelse(SBR_data_mutilation$Image_Category == 'Happy', 1,0)
# SBR_data_mutilation$mutilation <- ifelse(SBR_data_mutilation$Image_Category == 'Mutilation', 1,0) #mutilation is coded as the intercept
SBR_data_mutilation$erotic <- ifelse(SBR_data_mutilation$Image_Category == 'Erotic', 1,0)

SBR_data_mutilation <- SBR_data_mutilation %>%
  mutate(response = scale(Response))
```

```{r}
SBR_data_mutilation_category.model <- lmer(response ~ threat + erotic + happy + neutral + 
                            (1 + threat + erotic + happy + neutral|participant) + 
                            (1|item), data = SBR_data_mutilation, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

#summary(SBR_data_mutilation_category.model)

coefs_mutilation_category <- data.frame(coef(summary(SBR_data_mutilation_category.model)))

coefs_mutilation_category$p.z <- 2 * (1 - pnorm(abs(coefs_mutilation_category$t.value)))
```

```{r}
#Mutilation results
coefs_mutilation_category
```

#####Happy
```{r}
SBR_data_happy <- SBR_data

#Declare Participant_Number and Image as factors
SBR_data_happy$participant <- factor(SBR_data_happy$Participant_Number)
SBR_data_happy$item <- factor(SBR_data_happy$Image_Set)

SBR_data_happy$threat <- ifelse(SBR_data_happy$Image_Category == 'Threat', 1,0) 
SBR_data_happy$neutral <- ifelse(SBR_data_happy$Image_Category == 'Neutral', 1,0)
# SBR_data_happy$happy <- ifelse(SBR_data_happy$Image_Category == 'happy', 1,0) #happy is coded as the intercept
SBR_data_happy$mutilation <- ifelse(SBR_data_happy$Image_Category == 'Mutilation', 1,0) 
SBR_data_happy$erotic <- ifelse(SBR_data_happy$Image_Category == 'Erotic', 1,0)

SBR_data_happy <- SBR_data_happy %>%
  mutate(response = scale(Response))
```

```{r}
SBR_data_happy_category.model <- lmer(response ~ threat + erotic + mutilation + neutral + 
                            (1 + threat + erotic + mutilation + neutral|participant) + 
                            (1|item), data = SBR_data_happy, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

coefs_happy_category <- data.frame(coef(summary(SBR_data_happy_category.model)))

coefs_happy_category$p.z <- 2 * (1 - pnorm(abs(coefs_happy_category$t.value)))
```

```{r}
#Happy results
coefs_happy_category
```

#####Erotic
```{r}
SBR_data_erotic <- SBR_data

#Declare Participant_Number and Image as factors
SBR_data_erotic$participant <- factor(SBR_data_erotic$Participant_Number)
SBR_data_erotic$item <- factor(SBR_data_erotic$Image_Set)

SBR_data_erotic$threat <- ifelse(SBR_data_erotic$Image_Category == 'Threat', 1,0) 
SBR_data_erotic$neutral <- ifelse(SBR_data_erotic$Image_Category == 'Neutral', 1,0)
SBR_data_erotic$happy <- ifelse(SBR_data_erotic$Image_Category == 'Happy', 1,0) 
SBR_data_erotic$mutilation <- ifelse(SBR_data_erotic$Image_Category == 'Mutilation', 1,0) 
# SBR_data_erotic$erotic <- ifelse(SBR_data_erotic$Image_Category == 'Erotic', 1,0) #Erotic is coded as the intercept

SBR_data_erotic <- SBR_data_erotic %>%
  mutate(response = scale(Response))
```

```{r}

SBR_data_erotic_category.model <- lmer(response ~ threat + happy + mutilation + neutral + 
                            (1 + threat + happy + mutilation + neutral|participant) + 
                            (1|item), data = SBR_data_erotic, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

coefs_erotic_category <- data.frame(coef(summary(SBR_data_erotic_category.model)))

coefs_erotic_category$p.z <- 2 * (1 - pnorm(abs(coefs_erotic_category$t.value)))
```

```{r}
#Erotic results
coefs_erotic_category
```

#####Neutral
```{r}
SBR_data_neutral <- SBR_data

#Declare Participant_Number and Image as factors
SBR_data_neutral$participant <- factor(SBR_data_neutral$Participant_Number)
SBR_data_neutral$item <- factor(SBR_data_neutral$Image_Set)

SBR_data_neutral$threat <- ifelse(SBR_data_neutral$Image_Category == 'Threat', 1,0) 
# SBR_data_neutral$neutral <- ifelse(SBR_data_neutral$Image_Category == 'Neutral', 1,0) #Neutral is coded as the intercept
SBR_data_neutral$happy <- ifelse(SBR_data_neutral$Image_Category == 'Happy', 1,0) 
SBR_data_neutral$mutilation <- ifelse(SBR_data_neutral$Image_Category == 'Mutilation', 1,0) 
SBR_data_neutral$erotic <- ifelse(SBR_data_neutral$Image_Category == 'Erotic', 1,0)

SBR_data_neutral <- SBR_data_neutral %>%
  mutate(response = scale(Response))
```

```{r}

SBR_data_neutral_category.model <- lmer(response ~ threat + happy + mutilation + erotic + 
                            (1 + threat + happy + mutilation + erotic|participant) + 
                            (1|item), data = SBR_data_neutral, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

coefs_neutral_category <- data.frame(coef(summary(SBR_data_neutral_category.model)))

coefs_neutral_category$p.z <- 2 * (1 - pnorm(abs(coefs_neutral_category$t.value)))
```

```{r}
#Neutral results
coefs_neutral_category
```

Although we did not find a significant Condition x Category interaction, we still looked at the condition effects on the negative image categories since this is what one of our hypothesis predicted. We looked at the simple effects of the Threat and Mutilation image categories.

####Run the models:

#####Threat category
```{r}
SBR_threat_category_condition <- SBR_data

#Declare Participant_Number and Image as factors
SBR_threat_category_condition$participant <- factor(SBR_threat_category_condition$Participant_Number)
SBR_threat_category_condition$item <- factor(SBR_threat_category_condition$Image_Set)

#Dummy coding for simple effects
SBR_threat_category_condition$face_context <- ifelse(SBR_threat_category_condition$Condition == 'Face-context', 1, 0) #Context-only is coded as the intercept

# SBR_threat_category_condition$threat <- ifelse(SBR_threat_category_condition$Image_Category == 'Threat', 1,0) #threat is coded as the intercept
SBR_threat_category_condition$neutral <- ifelse(SBR_threat_category_condition$Image_Category == 'Neutral', 1,0)
SBR_threat_category_condition$happy <- ifelse(SBR_threat_category_condition$Image_Category == 'Happy', 1,0)
SBR_threat_category_condition$mutilation <- ifelse(SBR_threat_category_condition$Image_Category == 'Mutilation', 1,0)
SBR_threat_category_condition$erotic <- ifelse(SBR_threat_category_condition$Image_Category == 'Erotic', 1,0)

#Z-score
SBR_threat_category_condition <- SBR_threat_category_condition %>%
  mutate(response = scale(Response))
```

```{r}

SBR_threat_category_condition.model <- lmer(response ~ face_context + mutilation + erotic + happy + neutral + face_context:neutral + face_context:mutilation + face_context:erotic + face_context:happy + 
                            (1 + face_context + mutilation + erotic + happy + neutral + face_context:neutral + face_context:mutilation + face_context:erotic + face_context:happy|participant) + 
                            (1 + face_context|item), data = SBR_threat_category_condition, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

SBR_coefs_threat_category <- data.frame(coef(summary(SBR_threat_category_condition.model)))

SBR_coefs_threat_category$p.z <- 2 * (1 - pnorm(abs(SBR_coefs_threat_category$t.value)))
```

```{r}
#Threat results
SBR_coefs_threat_category
```
#####Mutilation
```{r}
SBR_mutilation_category_condition <- SBR_data

#Declare Participant_Number and Image as factors
SBR_mutilation_category_condition$participant <- factor(SBR_mutilation_category_condition$Participant_Number)
SBR_mutilation_category_condition$item <- factor(SBR_mutilation_category_condition$Image_Set)

#Dummy coding for simple effects
SBR_mutilation_category_condition$face_context <- ifelse(SBR_mutilation_category_condition$Condition == 'Face-context', 1, 0)  #Context-only is coded as the intercept

SBR_mutilation_category_condition$threat <- ifelse(SBR_mutilation_category_condition$Image_Category == 'Threat', 1,0) 
SBR_mutilation_category_condition$neutral <- ifelse(SBR_mutilation_category_condition$Image_Category == 'Neutral', 1,0)
SBR_mutilation_category_condition$happy <- ifelse(SBR_mutilation_category_condition$Image_Category == 'Happy', 1,0)
# SBR_mutilation_category_condition$mutilation <- ifelse(SBR_mutilation_category_condition$Image_Category == 'Mutilation', 1,0) #mutilation is coded as the intercept
SBR_mutilation_category_condition$erotic <- ifelse(SBR_mutilation_category_condition$Image_Category == 'Erotic', 1,0)

#Z-score
SBR_mutilation_category_condition <- SBR_mutilation_category_condition %>%
  mutate(response = scale(Response))
```

```{r}

SBR_mutilation_category_condition.model <- lmer(response ~ face_context + threat + erotic + happy + neutral + face_context:neutral + face_context:threat + face_context:erotic + face_context:happy + 
                            (1 + face_context + threat + erotic + happy + neutral + face_context:neutral + face_context:threat + face_context:erotic + face_context:happy|participant) + 
                            (1 + face_context|item), data = SBR_mutilation_category_condition, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

SBR_coefs_mutilation_category <- data.frame(coef(summary(SBR_mutilation_category_condition.model)))

SBR_coefs_mutilation_category$p.z <- 2 * (1 - pnorm(abs(SBR_coefs_mutilation_category$t.value)))
SBR_coefs_mutilation_category
```
##############################################################################################################
#######################################Behavioural data analysis##############################################

##########################
#########Valence data analysis
##########################

```{r}
valence_data <- read_csv('Data.csv') %>%
  filter(Measure == 'valence_rating')
```

####Main effects:
```{r}
#Declare Participant_Number and Image as factors
valence_data$participant <- factor(valence_data$Participant_Number)
valence_data$item <- factor(valence_data$Image_Set)

#Deviation code the image categories
valence_data$erotic <- scale(ifelse(valence_data$Image_Category == 'Erotic', 1,0),
                             center = TRUE, scale = FALSE)
valence_data$mutilation <- scale(ifelse(valence_data$Image_Category == 'Mutilation', 1,0),
                             center = TRUE, scale = FALSE)
valence_data$threat <- scale(ifelse(valence_data$Image_Category == 'Threat', 1,0),
                             center = TRUE, scale = FALSE)
valence_data$happy <- scale(ifelse(valence_data$Image_Category == 'Happy', 1,0),
                             center = TRUE, scale = FALSE)

#Deviation code the condition
valence_data$face_context <- scale(ifelse(valence_data$Condition == 'Face-context', 1,0),
                             center = T, scale = F)
```


```{r}
#Maximal model
valence.model <- glmer(Response ~ face_context + mutilation + threat + erotic + happy + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy +
                                  (1 + face_context + mutilation + threat + erotic + happy + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy|participant) +
                                  (1 + face_context|item), data = valence_data, family = poisson(),
                                control = glmerControl("bobyqa", calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

```

```{r}
#Condition model
valence.condition_model <- glmer(Response ~ mutilation + threat + erotic + happy + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy +
                                  (1 + face_context + mutilation + threat + erotic + happy + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy|participant) +
                                  (1 + face_context|item), data = valence_data, family = poisson(),
                                control = glmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 1000000)))

valence.condition_effect <- anova(valence.model, valence.condition_model)
```

```{r}
valence.condition_effect
```


```{r}
#Category model
valence.category_model <- glmer(Response ~ face_context + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy +
                                  (1 + face_context + mutilation + threat + erotic + happy + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy|participant) +
                                  (1 + face_context|item), data = valence_data, family = poisson(),
                                control = glmerControl("bobyqa", calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

valence.category_effect <- anova(valence.model, valence.category_model)
```

```{r}
valence.category_effect
```


```{r}
#Category-condition interaction
valence.category_condition_model <- glmer(Response ~ face_context + mutilation + threat + erotic + happy +
                                  (1 + face_context + mutilation + threat + erotic + happy + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy|participant) +
                                  (1 + face_context|item), data = valence_data, family = poisson(),
                                control = glmerControl("bobyqa", calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

valence.category_condition_effect <- anova(valence.model, valence.category_condition_model)
```

```{r}
valence.category_condition_effect
```


###Valence simple effects:

####Valence simple effects at the level of the category

#####Threat
```{r}
valence_data_threat <- valence_data

#Declare subjID anf Item/Image as factors
valence_data_threat$participant <- factor(valence_data_threat$Participant_Number)
valence_data_threat$item <- factor(valence_data_threat$Image_Set)

#Dummy code the image categories
valence_data_threat$erotic <- ifelse(valence_data_threat$Image_Category == 'Erotic', 1,0)
valence_data_threat$mutilation <- ifelse(valence_data_threat$Image_Category == 'Mutilation', 1,0)
# valence_data_threat$threat <- iflese(valence_data_threat$Image_Category == 'Threat', 1,0)
valence_data_threat$happy <- ifelse(valence_data_threat$Image_Category == 'Happy', 1,0)
valence_data_threat$neutral <- ifelse(valence_data_threat$Image_Category == 'Neutral', 1,0)
```

```{r}
#Threat model
valence.threat_model <- glmer(Response ~ mutilation + neutral + erotic + happy +
                                  (1 + mutilation + neutral + erotic + happy|participant) +
                                  (1|item), data = valence_data_threat, family = poisson(),
                                control = glmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

summary(valence.threat_model)
```

#####Mutilation
```{r}
valence_data_mutilation <- valence_data

#Declare Participant_Number and Image as factors
valence_data_mutilation$participant <- factor(valence_data_mutilation$Participant_Number)
valence_data_mutilation$item <- factor(valence_data_mutilation$Image_Set)

#Dummy code the image categories
valence_data_mutilation$erotic <- ifelse(valence_data_mutilation$Image_Category == 'Erotic', 1,0)
# valence_data_mutilation$mutilation <- ifelse(valence_data_mutilation$Image_Category == 'Mutilation', 1,0)
valence_data_mutilation$threat <- ifelse(valence_data_mutilation$Image_Category == 'Threat', 1,0)
valence_data_mutilation$happy <- ifelse(valence_data_mutilation$Image_Category == 'Happy', 1,0)
valence_data_mutilation$neutral <- ifelse(valence_data_mutilation$Image_Category == 'Neutral', 1,0)
```

```{r}
#Mutilation model
valence.mutilation_model <- glmer(Response ~ threat + neutral + erotic + happy +
                                  (1 + threat + neutral + erotic + happy|participant) +
                                  (1|item), data = valence_data_mutilation, family = poisson(),
                                control = glmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

summary(valence.mutilation_model)
```

#####Happy
```{r}
valence_data_happy <- valence_data

#Declare Participant_Number and Image as factors
valence_data_happy$participant <- factor(valence_data_happy$Participant_Number)
valence_data_happy$item <- factor(valence_data_happy$Image_Set)

#Dummy code the image categories
valence_data_happy$erotic <- ifelse(valence_data_happy$Image_Category == 'Erotic', 1,0)
valence_data_happy$mutilation <- ifelse(valence_data_happy$Image_Category == 'Mutilation', 1,0)
valence_data_happy$threat <- ifelse(valence_data_happy$Image_Category == 'Threat', 1,0)
# valence_data_happy$happy <- ifelse(valence_data_happy$Image_Category == 'happy', 1,0)
valence_data_happy$neutral <- ifelse(valence_data_happy$Image_Category == 'Neutral', 1,0)
```

```{r}
#happy model
valence.happy_model <- glmer(Response ~ threat + neutral + erotic + mutilation +
                                  (1 + threat + neutral + erotic + mutilation|participant) +
                                  (1|item), data = valence_data_happy, family = poisson(),
                                control = glmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 1000000)))

summary(valence.happy_model)
```

#####Erotic
```{r}
valence_data_erotic <- valence_data

#Declare Participant_Number and Image as factors
valence_data_erotic$participant <- factor(valence_data_erotic$Participant_Number)
valence_data_erotic$item <- factor(valence_data_erotic$Image_Set)

#Dummy code the image categories
# valence_data_erotic$erotic <- ifelse(valence_data_erotic$Image_Category == 'Erotic', 1,0)
valence_data_erotic$mutilation <- ifelse(valence_data_erotic$Image_Category == 'Mutilation', 1,0)
valence_data_erotic$threat <- ifelse(valence_data_erotic$Image_Category == 'Threat', 1,0)
valence_data_erotic$happy <- ifelse(valence_data_erotic$Image_Category == 'Happy', 1,0)
valence_data_erotic$neutral <- ifelse(valence_data_erotic$Image_Category == 'Neutral', 1,0)
```

```{r}
#Erotic model
valence.erotic_model <- glmer(Response ~ threat + neutral + happy + mutilation +
                                  (1 + threat + neutral + happy + mutilation|participant) +
                                  (1|item), data = valence_data_erotic, family = poisson(),
                                control = glmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

summary(valence.erotic_model)
```

#####Neutral
```{r}
valence_data_neutral <- valence_data

#Declare Participant_Number and Image as factors
valence_data_neutral$participant <- factor(valence_data_neutral$Participant_Number)
valence_data_neutral$item <- factor(valence_data_neutral$Image_Set)

#Dummy code the image categories
valence_data_neutral$erotic <- ifelse(valence_data_neutral$Image_Category == 'Erotic', 1,0)
valence_data_neutral$mutilation <- ifelse(valence_data_neutral$Image_Category == 'Mutilation', 1,0)
valence_data_neutral$threat <- ifelse(valence_data_neutral$Image_Category == 'Threat', 1,0)
valence_data_neutral$happy <- ifelse(valence_data_neutral$Image_Category == 'Happy', 1,0)
# valence_data_neutral$neutral <- ifelse(valence_data_neutral$Image_Category == 'Neutral', 1,0)
```

```{r}
#Neutral model
valence.neutral_model <- glmer(Response ~ threat + erotic + happy + mutilation +
                                  (1 + threat + erotic + happy + mutilation|participant) +
                                  (1|item), data = valence_data_neutral, family = poisson(),
                                control = glmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

summary(valence.neutral_model)
```

###Valence simple effects at the level of interaction effects:

#####Threat
```{r}
valence_threat_category_condition <- valence_data

#Declare Participant_Number and Image as factors
valence_threat_category_condition$participant <- factor(valence_threat_category_condition$Participant_Number)
valence_threat_category_condition$item <- factor(valence_threat_category_condition$Image_Set)

#Dummy code the image categories
valence_threat_category_condition$erotic <- ifelse(valence_threat_category_condition$Image_Category == 'Erotic', 1,0)
valence_threat_category_condition$mutilation <- ifelse(valence_threat_category_condition$Image_Category == 'Mutilation', 1,0)
# valence_threat_category_condition$threat <- ifelse(valence_threat_category_condition$Image_Category == 'Threat', 1,0)
valence_threat_category_condition$happy <- ifelse(valence_threat_category_condition$Image_Category == 'Happy', 1,0)
valence_threat_category_condition$neutral <- ifelse(valence_threat_category_condition$Image_Category == 'Neutral', 1,0)

#Deviation code the condition
valence_threat_category_condition$face_context <- scale(ifelse(valence_threat_category_condition$Condition == 'Face-context', 1,0), center = T, scale = F) #Context-only is coded as the intercept
```

```{r}
#Threat model
valence.threat_category_condition_model <- glmer(Response ~ face_context + mutilation + neutral + erotic + happy + face_context:mutilation + face_context:neutral + face_context:erotic + face_context:happy +
                                  (1 + face_context + mutilation + neutral + erotic + happy + face_context:mutilation + face_context:neutral + face_context:erotic + face_context:happy|participant) +
                                  (1 + face_context|item), data = valence_threat_category_condition, family = poisson(),
                                control = glmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

summary(valence.threat_category_condition_model)
```

#####Mutilation
```{r}
valence_mutilation_category_condition <- valence_data

#Declare Participant_Number and Image as factors
valence_mutilation_category_condition$participant <- factor(valence_mutilation_category_condition$Participant_Number)
valence_mutilation_category_condition$item <- factor(valence_mutilation_category_condition$Image_Set)

#Dummy code the image categories
valence_mutilation_category_condition$erotic <- ifelse(valence_mutilation_category_condition$Image_Category == 'Erotic', 1,0)
# valence_mutilation_category_condition$mutilation <- ifelse(valence_mutilation_category_condition$Image_Category == 'Mutilation', 1,0)
valence_mutilation_category_condition$threat <- ifelse(valence_mutilation_category_condition$Image_Category == 'Threat', 1,0)
valence_mutilation_category_condition$happy <- ifelse(valence_mutilation_category_condition$Image_Category == 'Happy', 1,0)
valence_mutilation_category_condition$neutral <- ifelse(valence_mutilation_category_condition$Image_Category == 'Neutral', 1,0)

#Deviation code the condition
valence_mutilation_category_condition$face_context <- scale(ifelse(valence_mutilation_category_condition$Condition == 'Face-context', 1,0), center = T, scale = F)
```

```{r}
#Mutilation model
valence.mutilation_category_condition_model <- glmer(Response ~ face_context + threat + neutral + erotic + happy + face_context:threat + face_context:neutral + face_context:erotic + face_context:happy +
                                  (1 + face_context + threat + neutral + erotic + happy + face_context:threat + face_context:neutral + face_context:erotic + face_context:happy|participant) +
                                  (1 + face_context|item), data = valence_mutilation_category_condition, family = poisson(),
                                control = glmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

summary(valence.mutilation_category_condition_model)
```

#####Happy
```{r}
valence_happy_category_condition <- valence_data

#Declare Participant_Number and Image as factors
valence_happy_category_condition$participant <- factor(valence_happy_category_condition$Participant_Number)
valence_happy_category_condition$item <- factor(valence_happy_category_condition$Image_Set)

#Dummy code the image categories
valence_happy_category_condition$erotic <- ifelse(valence_happy_category_condition$Image_Category == 'Erotic', 1,0)
valence_happy_category_condition$mutilation <- ifelse(valence_happy_category_condition$Image_Category == 'Mutilation', 1,0)
valence_happy_category_condition$threat <- ifelse(valence_happy_category_condition$Image_Category == 'Threat', 1,0)
# valence_happy_category_condition$happy <- ifelse(valence_happy_category_condition$Image_Category == 'happy', 1,0)
valence_happy_category_condition$neutral <- ifelse(valence_happy_category_condition$Image_Category == 'Neutral', 1,0)

#Deviation code the condition
valence_happy_category_condition$face_context <- scale(ifelse(valence_happy_category_condition$Condition == 'Face-context', 1,0), center = T, scale = F)
```

```{r}
#happy model
valence.happy_category_condition_model <- glmer(Response ~ face_context + threat + neutral + erotic + mutilation + face_context:threat + face_context:neutral + face_context:erotic + face_context:mutilation +
                                  (1 + face_context + threat + neutral + erotic + mutilation + face_context:threat + face_context:neutral + face_context:erotic + face_context:mutilation|participant) +
                                  (1 + face_context|item), data = valence_happy_category_condition, family = poisson(),
                                control = glmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 1000000)))

summary(valence.happy_category_condition_model)
```

#####Erotic
```{r}
valence_erotic_category_condition <- valence_data

#Declare Participant_Number and Image as factors
valence_erotic_category_condition$participant <- factor(valence_erotic_category_condition$Participant_Number)
valence_erotic_category_condition$item <- factor(valence_erotic_category_condition$Image_Set)

#Dummy code the image categories
# valence_erotic_category_condition$erotic <- ifelse(valence_erotic_category_condition$Image_Category == 'Erotic', 1,0)
valence_erotic_category_condition$mutilation <- ifelse(valence_erotic_category_condition$Image_Category == 'Mutilation', 1,0)
valence_erotic_category_condition$threat <- ifelse(valence_erotic_category_condition$Image_Category == 'Threat', 1,0)
valence_erotic_category_condition$happy <- ifelse(valence_erotic_category_condition$Image_Category == 'Happy', 1,0)
valence_erotic_category_condition$neutral <- ifelse(valence_erotic_category_condition$Image_Category == 'Neutral', 1,0)

#Deviation code the condition
valence_erotic_category_condition$face_context <- scale(ifelse(valence_erotic_category_condition$Condition == 'Face-context', 1,0), center = T, scale = F)
```

```{r}
#Erotic model
valence.erotic_category_condition_model <- glmer(Response ~ face_context + threat + neutral + happy + mutilation + face_context:threat + face_context:neutral + face_context:happy + face_context:mutilation +
                                  (1 + face_context + threat + neutral + happy + mutilation + face_context:threat + face_context:neutral + face_context:happy + face_context:mutilation|participant) +
                                  (1 + face_context|item), data = valence_erotic_category_condition, family = poisson(),
                                control = glmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

summary(valence.erotic_category_condition_model)
```

#####Neutral
```{r}
valence_neutral_category_condition <- valence_data

#Declare Participant_Number and Image as factors
valence_neutral_category_condition$participant <- factor(valence_neutral_category_condition$Participant_Number)
valence_neutral_category_condition$item <- factor(valence_neutral_category_condition$Image_Set)

#Dummy code the image categories
valence_neutral_category_condition$erotic <- ifelse(valence_neutral_category_condition$Image_Category == 'Erotic', 1,0)
valence_neutral_category_condition$mutilation <- ifelse(valence_neutral_category_condition$Image_Category == 'Mutilation', 1,0)
valence_neutral_category_condition$threat <- ifelse(valence_neutral_category_condition$Image_Category == 'Threat', 1,0)
valence_neutral_category_condition$happy <- ifelse(valence_neutral_category_condition$Image_Category == 'Happy', 1,0)
# valence_neutral_category_condition$neutral <- ifelse(valence_neutral_category_condition$Image_Category == 'Neutral', 1,0)

#Deviation code the condition
valence_neutral_category_condition$face_context <- scale(ifelse(valence_neutral_category_condition$Condition == 'Face-context', 1,0), center = T, scale = F)
```

```{r}
#Neutral model
valence.neutral_category_condition_model <- glmer(Response ~ face_context + threat + erotic + happy + mutilation + face_context:threat + face_context:erotic + face_context:happy + face_context:mutilation +
                                  (1 + face_context + threat + erotic + happy + mutilation + face_context:threat + face_context:erotic + face_context:happy + face_context:mutilation|participant) +
                                  (1 + face_context|item), data = valence_neutral_category_condition, family = poisson(),
                                control = glmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 1000000)))

summary(valence.neutral_category_condition_model)
```


##########################
#########Arousal data analysis
##########################

```{r}
arousal_data <- read_csv('Data.csv') %>%
  filter(Measure == 'arousal_rating')
```

####Main effects:
```{r}
#Declare Participant_Number and Image as factors
arousal_data$participant <- factor(arousal_data$Participant_Number)
arousal_data$item <- factor(arousal_data$Image_Set)

#Deviation code the image categories
arousal_data$erotic <- scale(ifelse(arousal_data$Image_Category == 'Erotic', 1,0),
                             center = TRUE, scale = FALSE)
arousal_data$mutilation <- scale(ifelse(arousal_data$Image_Category == 'Mutilation', 1,0),
                             center = TRUE, scale = FALSE)
arousal_data$threat <- scale(ifelse(arousal_data$Image_Category == 'Threat', 1,0),
                             center = TRUE, scale = FALSE)
arousal_data$happy <- scale(ifelse(arousal_data$Image_Category == 'Happy', 1,0),
                             center = TRUE, scale = FALSE)

#Deviation code the condition
arousal_data$face_context <- scale(ifelse(arousal_data$Condition == 'Face-context', 1,0),
                             center = T, scale = F)
```

```{r}
#Maximal model
arousal.model <- glmer(Response ~ face_context + mutilation + threat + erotic + happy + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy +
                                  (1 + face_context + mutilation + threat + erotic + happy + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy|participant) +
                                  (1 + face_context|item), data = arousal_data, family = poisson(),
                                control = glmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

```

```{r}
#Condition model
arousal.condition_model <- glmer(Response ~ mutilation + threat + erotic + happy + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy +
                                  (1 + face_context + mutilation + threat + erotic + happy + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy|participant) +
                                  (1 + face_context|item), data = arousal_data, family = poisson(),
                                control = glmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

arousal.condition_effect <- anova(arousal.model, arousal.condition_model)
```

```{r}
arousal.condition_effect
```


```{r}
#Category model
arousal.category_model <- glmer(Response ~ face_context + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy +
                                  (1 + face_context + mutilation + threat + erotic + happy + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy|participant) +
                                  (1 + face_context|item), data = arousal_data, family = poisson(),
                                control = glmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

arousal.category_effect <- anova(arousal.model, arousal.category_model)
```

```{r}
arousal.category_effect
```


```{r}
#Category-condition interaction
arousal.category_condition_model <- glmer(Response ~ face_context + mutilation + threat + erotic + happy +
                                  (1 + face_context + mutilation + threat + erotic + happy + face_context:mutilation + face_context:threat + face_context:erotic + face_context:happy|participant) +
                                  (1 + face_context|item), data = arousal_data, family = poisson(),
                                control = glmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

arousal.category_condition_effect <- anova(arousal.model, arousal.category_condition_model)
```

```{r}
arousal.category_condition_effect
```


###Arousal simple effects at the level of the category:

#####Threat
```{r}
arousal_data_threat <- arousal_data

#Declare Participant_Number and Image as factors
arousal_data_threat$participant <- factor(arousal_data_threat$Participant_Number)
arousal_data_threat$item <- factor(arousal_data_threat$Image_Set)

#Dummy code the image categories
arousal_data_threat$erotic <- ifelse(arousal_data_threat$Image_Category == 'Erotic', 1,0)
arousal_data_threat$mutilation <- ifelse(arousal_data_threat$Image_Category == 'Mutilation', 1,0)
# arousal_data_threat$threat <- iflese(arousal_data_threat$Image_Category == 'Threat', 1,0)
arousal_data_threat$happy <- ifelse(arousal_data_threat$Image_Category == 'Happy', 1,0)
arousal_data_threat$neutral <- ifelse(arousal_data_threat$Image_Category == 'Neutral', 1,0)
```

```{r}
#Threat model
arousal.threat_model <- glmer(Response ~ mutilation + neutral + erotic + happy +
                                  (1 + mutilation + neutral + erotic + happy|participant) +
                                  (1|item), data = arousal_data_threat, family = poisson(),
                                control = glmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))
```

```{r}
summary(arousal.threat_model)
```

#####Mutilation
```{r}
arousal_data_mutilation <- arousal_data

#Declare Participant_Number and Image as factors
arousal_data_mutilation$participant <- factor(arousal_data_mutilation$Participant_Number)
arousal_data_mutilation$item <- factor(arousal_data_mutilation$Image_Set)

#Dummy code the image categories
arousal_data_mutilation$erotic <- ifelse(arousal_data_mutilation$Image_Category == 'Erotic', 1,0)
# arousal_data_mutilation$mutilation <- ifelse(arousal_data_mutilation$Image_Category == 'Mutilation', 1,0)
arousal_data_mutilation$threat <- ifelse(arousal_data_mutilation$Image_Category == 'Threat', 1,0)
arousal_data_mutilation$happy <- ifelse(arousal_data_mutilation$Image_Category == 'Happy', 1,0)
arousal_data_mutilation$neutral <- ifelse(arousal_data_mutilation$Image_Category == 'Neutral', 1,0)
```

```{r}
#Mutilation model
arousal.mutilation_model <- glmer(Response ~ threat + neutral + erotic + happy +
                                  (1 + threat + neutral + erotic + happy|participant) +
                                  (1|item), data = arousal_data_mutilation, family = poisson(),
                                control = glmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))
```

```{r}
summary(arousal.mutilation_model)
```

#####Happy
```{r}
arousal_data_happy <- arousal_data

#Declare subjID anf Item/Image as factors
arousal_data_happy$participant <- factor(arousal_data_happy$Participant_Number)
arousal_data_happy$item <- factor(arousal_data_happy$Image_Set)

#Dummy code the image categories
arousal_data_happy$erotic <- ifelse(arousal_data_happy$Image_Category == 'Erotic', 1,0)
arousal_data_happy$mutilation <- ifelse(arousal_data_happy$Image_Category == 'Mutilation', 1,0)
arousal_data_happy$threat <- ifelse(arousal_data_happy$Image_Category == 'Threat', 1,0)
# arousal_data_happy$happy <- ifelse(arousal_data_happy$Image_Category == 'happy', 1,0)
arousal_data_happy$neutral <- ifelse(arousal_data_happy$Image_Category == 'Neutral', 1,0)
```

```{r}
#happy model
arousal.happy_model <- glmer(Response ~ threat + neutral + erotic + mutilation +
                                  (1 + threat + neutral + erotic + mutilation|participant) +
                                  (1|item), data = arousal_data_happy, family = poisson(),
                                control = glmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 1000000)))
```

```{r}
summary(arousal.happy_model)
```

#####Erotic
```{r}
arousal_data_erotic <- arousal_data

#Declare subjID anf Item/Image as factors
arousal_data_erotic$participant <- factor(arousal_data_erotic$Participant_Number)
arousal_data_erotic$item <- factor(arousal_data_erotic$Image_Set)

#Dummy code the image categories
# arousal_data_erotic$erotic <- ifelse(arousal_data_erotic$Image_Category == 'Erotic', 1,0)
arousal_data_erotic$mutilation <- ifelse(arousal_data_erotic$Image_Category == 'Mutilation', 1,0)
arousal_data_erotic$threat <- ifelse(arousal_data_erotic$Image_Category == 'Threat', 1,0)
arousal_data_erotic$happy <- ifelse(arousal_data_erotic$Image_Category == 'Happy', 1,0)
arousal_data_erotic$neutral <- ifelse(arousal_data_erotic$Image_Category == 'Neutral', 1,0)
```

```{r}
#Erotic model
arousal.erotic_model <- glmer(Response ~ threat + neutral + happy + mutilation +
                                  (1 + threat + neutral + happy + mutilation|participant) +
                                  (1|item), data = arousal_data_erotic, family = poisson(),
                                control = glmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))
```

```{r}
summary(arousal.erotic_model)
```

#####Neutral
```{r}
arousal_data_neutral <- arousal_data

#Declare subjID anf Item/Image as factors
arousal_data_neutral$participant <- factor(arousal_data_neutral$Participant_Number)
arousal_data_neutral$item <- factor(arousal_data_neutral$Image_Set)

#Dummy code the image categories
arousal_data_neutral$erotic <- ifelse(arousal_data_neutral$Image_Category == 'Erotic', 1,0)
arousal_data_neutral$mutilation <- ifelse(arousal_data_neutral$Image_Category == 'Mutilation', 1,0)
arousal_data_neutral$threat <- ifelse(arousal_data_neutral$Image_Category == 'Threat', 1,0)
arousal_data_neutral$happy <- ifelse(arousal_data_neutral$Image_Category == 'Happy', 1,0)
# arousal_data_neutral$neutral <- ifelse(arousal_data_neutral$Image_Category == 'Neutral', 1,0)
```

```{r}
#Neutral model
arousal.neutral_model <- glmer(Response ~ threat + erotic + happy + mutilation +
                                  (1 + threat + erotic + happy + mutilation|participant) +
                                  (1|item), data = arousal_data_neutral, family = poisson(),
                                control = glmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))
```

```{r}
summary(arousal.neutral_model)
```

###Arousal simple effects at the level of the condition:
```{r}
arousal_condition_simple_effects <- arousal_data

#Declare subjID anf Item/Image as factors
arousal_condition_simple_effects$participant <- factor(arousal_condition_simple_effects$Participant_Number)
arousal_condition_simple_effects$item <- factor(arousal_condition_simple_effects$Image_Set)

#Dummy code the image categories
arousal_condition_simple_effects$condition <- ifelse(arousal_condition_simple_effects$Condition == 'Face-context', 1,0)
```

```{r}
arousal_condition_simple_effects.model <- glmer(Response ~ condition + (1 + condition|participant) + (1|item), data = arousal_condition_simple_effects,
                                                family = poisson(),
                                                control = glmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))
```

```{r}
summary(arousal_condition_simple_effects.model)
```


## Supplementary analyses

Section B. In this section of the supplemental materials file we test whether including the outliers has a drastic effect on our facial electromyography results. We perform the same analysis to the main analysis (everything above from this point) and removed the condition "Outliers == 'No'" when loading and filtering the data for the corrugator, zygomatic and orbicularis data.

Section C.In this section of the supplemental materials file we assess whether the Condition effects observed in the Valence and Arousal ratings and the Corrugator and Zygomatic facial muscles are driven solely by the facial information available in the Face-context condition. We run a model with the data collected solely from the Context-only condition.

####Valence analysis
```{r}
valence_data_context_only <- read_csv('Data.csv') %>%
  filter(Measure == 'valence_rating', Condition == 'Context-only')
```

####Main effects:
```{r}
#Declare Participant_Number and Image as factors
valence_data_context_only$participant <- factor(valence_data_context_only$Participant_Number)
valence_data_context_only$item <- factor(valence_data_context_only$Image_Set)

#Deviation code the image categories
valence_data_context_only$erotic <- scale(ifelse(valence_data_context_only$Image_Category == 'Erotic', 1,0),
                             center = TRUE, scale = FALSE)
valence_data_context_only$mutilation <- scale(ifelse(valence_data_context_only$Image_Category == 'Mutilation', 1,0),
                             center = TRUE, scale = FALSE)
valence_data_context_only$threat <- scale(ifelse(valence_data_context_only$Image_Category == 'Threat', 1,0),
                             center = TRUE, scale = FALSE)
valence_data_context_only$happy <- scale(ifelse(valence_data_context_only$Image_Category == 'Happy', 1,0),
                             center = TRUE, scale = FALSE)
```


```{r}
#Maximal model
valence_data_context_only.model <- glmer(Response ~ mutilation + threat + erotic + happy + 
                                  (1 + mutilation + threat + erotic + happy|participant) +
                                  (1|item), data = valence_data_context_only, family = poisson(),
                                control = glmerControl("bobyqa", calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

```

```{r}
#Effect of category
valence_data_context_only_no_category.model <- glmer(Response ~ 
                                  (1 + mutilation + threat + erotic + happy|participant) +
                                  (1|item), data = valence_data_context_only, family = poisson(),
                                control = glmerControl("bobyqa", calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

valence_category_context_only_effect <- anova(valence_data_context_only.model, valence_data_context_only_no_category.model)
```

```{r}
#Category results
valence_category_context_only_effect
```

####Arousal analysis
```{r}
arousal_data_context_only <- read_csv('Data.csv') %>%
  filter(Measure == 'arousal_rating', Condition == 'Context-only')
```

####Main effects:
```{r}
#Declare Participant_Number and Image as factors
arousal_data_context_only$participant <- factor(arousal_data_context_only$Participant_Number)
arousal_data_context_only$item <- factor(arousal_data_context_only$Image_Set)

#Deviation code the image categories
arousal_data_context_only$erotic <- scale(ifelse(arousal_data_context_only$Image_Category == 'Erotic', 1,0),
                             center = TRUE, scale = FALSE)
arousal_data_context_only$mutilation <- scale(ifelse(arousal_data_context_only$Image_Category == 'Mutilation', 1,0),
                             center = TRUE, scale = FALSE)
arousal_data_context_only$threat <- scale(ifelse(arousal_data_context_only$Image_Category == 'Threat', 1,0),
                             center = TRUE, scale = FALSE)
arousal_data_context_only$happy <- scale(ifelse(arousal_data_context_only$Image_Category == 'Happy', 1,0),
                             center = TRUE, scale = FALSE)

```

```{r}
#Maximal model
arousal_data_context_only.model <- glmer(Response ~ mutilation + threat + erotic + happy + 
                                  (1 + mutilation + threat + erotic + happy|participant) +
                                  (1|item), data = arousal_data_context_only, family = poisson(),
                                control = glmerControl("bobyqa", calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

```

```{r}
#Effect of category
arousal_data_context_only_no_category.model <- glmer(Response ~ 
                                  (1 + mutilation + threat + erotic + happy|participant) +
                                  (1|item), data = arousal_data_context_only, family = poisson(),
                                control = glmerControl("bobyqa", calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

arousal_category_context_only_effect <- anova(arousal_data_context_only.model, arousal_data_context_only_no_category.model)
```

```{r}
#Category results
arousal_category_context_only_effect
```

####Corrugator analysis
```{r include=FALSE}
#Read in the corrugator data
corrugator_data_context_only <- read_csv('Data.csv') %>%
  filter(Measure == 'corrugator_EMG', Outlier == 'No', Condition == 'Context-only') #only selecting the data without the outliers i.e. Outlier == 'No'
```

```{r}
#Deviation code the image categories
corrugator_data_context_only$erotic <- scale(ifelse(corrugator_data_context_only$Image_Category == 'Erotic', 1,0),
                             center = TRUE, scale = FALSE)
corrugator_data_context_only$mutilation <- scale(ifelse(corrugator_data_context_only$Image_Category == 'Mutilation', 1,0),
                             center = TRUE, scale = FALSE)
corrugator_data_context_only$threat <- scale(ifelse(corrugator_data_context_only$Image_Category == 'Threat', 1,0),
                             center = TRUE, scale = FALSE)
corrugator_data_context_only$happy <- scale(ifelse(corrugator_data_context_only$Image_Category == 'Happy', 1,0),
                             center = TRUE, scale = FALSE)

#Declare Participant_Number and Image as factors
corrugator_data_context_only$participant <- factor(corrugator_data_context_only$Participant_Number)
corrugator_data_context_only$item <- factor(corrugator_data_context_only$Image_Set)

#Convert corrugator EMG responses to Z-scores
corrugator_data_context_only <- corrugator_data_context_only %>%
  mutate(response = scale(Response))
```

####Run the models:

```{r}
#Maximal model
corrugator_data_context_only.model <- lmer(response ~ mutilation + threat + erotic + happy +
                            (1 + mutilation + threat + erotic + happy|participant) +
                            (1|item), data = corrugator_data_context_only, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 1000000)))
```

```{r}
#Effect of category
corrugator_data_context_only_no_category.model <- lmer(response ~ (1 + mutilation + threat + erotic + happy|participant) +
                            (1|item), data = corrugator_data_context_only, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 1000000)))

corrugator_category_context_only_effect <- anova(corrugator_data_context_only.model, corrugator_data_context_only_no_category.model, type = 'chisq')
```

```{r}
#Category results
corrugator_category_context_only_effect
```

###Corrugator simple effects:

Since the **Category** main effect was significant, we are now looking at the simple effects.

####Run the models:

#####Threat category
```{r}
corrugator_data_context_only_threat_category <- corrugator_data_context_only

#Declare Participant_Number and Image as factors
corrugator_data_context_only_threat_category$participant <- factor(corrugator_data_context_only_threat_category$Participant_Number)
corrugator_data_context_only_threat_category$item <- factor(corrugator_data_context_only_threat_category$Image_Set)

# corrugator_data_context_only_threat_category$threat <- ifelse(corrugator_data_context_only_threat_category$Image_Category == 'Threat', 1,0) #threat is coded as the intercept
corrugator_data_context_only_threat_category$neutral <- ifelse(corrugator_data_context_only_threat_category$Image_Category == 'Neutral', 1,0)
corrugator_data_context_only_threat_category$happy <- ifelse(corrugator_data_context_only_threat_category$Image_Category == 'Happy', 1,0)
corrugator_data_context_only_threat_category$mutilation <- ifelse(corrugator_data_context_only_threat_category$Image_Category == 'Mutilation', 1,0)
corrugator_data_context_only_threat_category$erotic <- ifelse(corrugator_data_context_only_threat_category$Image_Category == 'Erotic', 1,0)

#Convert corrugator EMG responses to Z-scores
corrugator_data_context_only_threat_category <- corrugator_data_context_only_threat_category %>%
  mutate(response = scale(Response))
```

```{r}

corrugator_data_context_only_threat_category.model <- lmer(response ~ mutilation + erotic + happy + neutral +
                            (1 + mutilation + erotic + happy + neutral|participant) +
                            (1|item), data = corrugator_data_context_only_threat_category, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

corrugator_coefs_context_only_threat_category <- data.frame(coef(summary(corrugator_data_context_only_threat_category.model)))

corrugator_coefs_context_only_threat_category$p.z <- 2 * (1 - pnorm(abs(corrugator_coefs_context_only_threat_category$t.value)))
```

```{r}
#Threat results
corrugator_coefs_context_only_threat_category
```

#####Mutilation
```{r}
corrugator_data__context_only_mutilation_category <- corrugator_data_context_only

#Declare Participant_Number and Image as factors
corrugator_data__context_only_mutilation_category$participant <- factor(corrugator_data__context_only_mutilation_category$Participant_Number)
corrugator_data__context_only_mutilation_category$item <- factor(corrugator_data__context_only_mutilation_category$Image_Set)

corrugator_data__context_only_mutilation_category$threat <- ifelse(corrugator_data__context_only_mutilation_category$Image_Category == 'Threat', 1,0)
corrugator_data__context_only_mutilation_category$neutral <- ifelse(corrugator_data__context_only_mutilation_category$Image_Category == 'Neutral', 1,0)
corrugator_data__context_only_mutilation_category$happy <- ifelse(corrugator_data__context_only_mutilation_category$Image_Category == 'Happy', 1,0)
# corrugator_data__context_only_mutilation_category$mutilation <- ifelse(corrugator_data__context_only_mutilation_category$Image_Category == 'Mutilation', 1,0) #mutilation is coded as the intercept
corrugator_data__context_only_mutilation_category$erotic <- ifelse(corrugator_data__context_only_mutilation_category$Image_Category == 'Erotic', 1,0)

#Convert corrugator EMG responses to Z-scores
corrugator_data__context_only_mutilation_category <- corrugator_data__context_only_mutilation_category %>%
  mutate(response = scale(Response))
```

```{r}

corrugator_data__context_only_mutilation_category.model <- lmer(response ~ threat + erotic + happy + neutral +
                            (1 + threat + erotic + happy + neutral|participant) +
                            (1|item), data = corrugator_data_mutilation_category, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

corrugator_coefs_context_only_mutilation_category <- data.frame(coef(summary(corrugator_data__context_only_mutilation_category.model)))

corrugator_coefs_context_only_mutilation_category$p.z <- 2 * (1 - pnorm(abs(corrugator_coefs_context_only_mutilation_category$t.value)))
```

```{r}
#Mutilation results
corrugator_coefs_context_only_mutilation_category
```

#####Happy
```{r}
corrugator_data_context_only_happy_category <- corrugator_data_context_only

#Declare Participant_Number and Image as factors
corrugator_data_context_only_happy_category$participant <- factor(corrugator_data_context_only_happy_category$Participant_Number)
corrugator_data_context_only_happy_category$item <- factor(corrugator_data_context_only_happy_category$Image_Set)

corrugator_data_context_only_happy_category$threat <- ifelse(corrugator_data_context_only_happy_category$Image_Category == 'Threat', 1,0)
corrugator_data_context_only_happy_category$neutral <- ifelse(corrugator_data_context_only_happy_category$Image_Category == 'Neutral', 1,0)
# corrugator_data_context_only_happy_category$happy <- ifelse(corrugator_data_context_only_happy_category$Image_Category == 'Happy', 1,0) #happy is coded as the intercept
corrugator_data_context_only_happy_category$mutilation <- ifelse(corrugator_data_context_only_happy_category$Image_Category == 'Mutilation', 1,0)
corrugator_data_context_only_happy_category$erotic <- ifelse(corrugator_data_context_only_happy_category$Image_Category == 'Erotic', 1,0)

#Convert corrugator EMG responses to Z-scores
corrugator_data_context_only_happy_category <- corrugator_data_context_only_happy_category %>%
  mutate(response = scale(Response))
```

```{r}

corrugator_data_context_only_happy_category.model <- lmer(response ~ threat + erotic + mutilation + neutral +
                            (1 + threat + erotic + mutilation + neutral|participant) +
                            (1|item), data = corrugator_data_context_only_happy_category, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

corrugator_coefs_context_only_happy_category <- data.frame(coef(summary(corrugator_data_context_only_happy_category.model)))

corrugator_coefs_context_only_happy_category$p.z <- 2 * (1 - pnorm(abs(corrugator_coefs_context_only_happy_category$t.value)))
```

```{r}
#Happy results
corrugator_coefs_context_only_happy_category
```

#####Erotic
```{r}
corrugator_data_context_only_erotic_category <- corrugator_data_context_only

#Declare Participant_Number and Image as factors
corrugator_data_context_only_erotic_category$participant <- factor(corrugator_data_context_only_erotic_category$Participant_Number)
corrugator_data_context_only_erotic_category$item <- factor(corrugator_data_context_only_erotic_category$Image_Set)

corrugator_data_context_only_erotic_category$threat <- ifelse(corrugator_data_context_only_erotic_category$Image_Category == 'Threat', 1,0)
corrugator_data_context_only_erotic_category$neutral <- ifelse(corrugator_data_context_only_erotic_category$Image_Category == 'Neutral', 1,0)
corrugator_data_context_only_erotic_category$happy <- ifelse(corrugator_data_context_only_erotic_category$Image_Category == 'Happy', 1,0)
corrugator_data_context_only_erotic_category$mutilation <- ifelse(corrugator_data_context_only_erotic_category$Image_Category == 'Mutilation', 1,0)
# corrugator_data_context_only_erotic_category$erotic <- ifelse(corrugator_data_context_only_erotic_category$imageCategory == 'Erotic', 1,0) #Erotic is coded as the intercept

#Convert corrugator EMG responses to Z-scores
corrugator_data_context_only_erotic_category <- corrugator_data_context_only_erotic_category %>%
  mutate(response = scale(Response))
```

```{r}

corrugator_data_context_only_erotic_category.model <- lmer(response ~ threat + happy + mutilation + neutral +
                            (1 + threat + happy + mutilation + neutral|participant) +
                            (1|item), data = corrugator_data_context_only_erotic_category, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

corrugator_coefs_context_only_erotic_category <- data.frame(coef(summary(corrugator_data_context_only_erotic_category.model)))

corrugator_coefs_context_only_erotic_category$p.z <- 2 * (1 - pnorm(abs(corrugator_coefs_context_only_erotic_category$t.value)))
```

```{r}
#Erotic results
corrugator_coefs_context_only_erotic_category
```

#####Neutral
```{r}
corrugator_data_context_only_neutral_category <- corrugator_data_context_only

#Declare Participant_Number and Image as factors
corrugator_data_context_only_neutral_category$participant <- factor(corrugator_data_context_only_neutral_category$Participant_Number)
corrugator_data_context_only_neutral_category$item <- factor(corrugator_data_context_only_neutral_category$Image_Set)

corrugator_data_context_only_neutral_category$threat <- ifelse(corrugator_data_context_only_neutral_category$Image_Category == 'Threat', 1,0)
# corrugator_data_context_only_neutral_category$neutral <- ifelse(corrugator_data_context_only_neutral_category$Image_Category == 'Neutral', 1,0) #Neutral is coded as the intercept
corrugator_data_context_only_neutral_category$happy <- ifelse(corrugator_data_context_only_neutral_category$Image_Category == 'Happy', 1,0)
corrugator_data_context_only_neutral_category$mutilation <- ifelse(corrugator_data_context_only_neutral_category$Image_Category == 'Mutilation', 1,0)
corrugator_data_context_only_neutral_category$erotic <- ifelse(corrugator_data_context_only_neutral_category$Image_Category == 'Erotic', 1,0)

#Convert corrugator EMG responses to Z-scores
corrugator_data_context_only_neutral_category <- corrugator_data_context_only_neutral_category %>%
  mutate(response = scale(Response))
```

```{r}

corrugator_data_context_only_neutral_category.model <- lmer(response ~ threat + happy + mutilation + erotic +
                            (1 + threat + happy + mutilation + erotic|participant) +
                            (1|item), data = corrugator_data_context_only_neutral_category, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 100000)))

corrugator_coefs_context_only_neutral_category <- data.frame(coef(summary(corrugator_data_context_only_neutral_category.model)))

corrugator_coefs_context_only_neutral_category$p.z <- 2 * (1 - pnorm(abs(corrugator_coefs_context_only_neutral_category$t.value)))
```

```{r}
#Neutral results
corrugator_coefs_context_only_neutral_category
```


###Zygomatic analysis
```{r include=FALSE}
#Read in the zygomatic data
zygomatic_data_context_only <- read_csv('Data.csv') %>%
  filter(Measure == 'zygomatic_EMG', Outlier == 'No', Condition == 'Context-only') #only selecting the data without the outliers i.e. Outlier == 'No'

```

```{r}
#Deviation code the image categories
zygomatic_data_context_only$erotic <- scale(ifelse(zygomatic_data_context_only$Image_Category == 'Erotic', 1,0),
                             center = TRUE, scale = FALSE)
zygomatic_data_context_only$mutilation <- scale(ifelse(zygomatic_data_context_only$Image_Category == 'Mutilation', 1,0),
                             center = TRUE, scale = FALSE)
zygomatic_data_context_only$threat <- scale(ifelse(zygomatic_data_context_only$Image_Category == 'Threat', 1,0),
                             center = TRUE, scale = FALSE)
zygomatic_data_context_only$happy <- scale(ifelse(zygomatic_data_context_only$Image_Category == 'Happy', 1,0),
                             center = TRUE, scale = FALSE)

#Declare Participant_Number and Image as factors
zygomatic_data_context_only$participant <- factor(zygomatic_data_context_only$Participant_Number)
zygomatic_data_context_only$item <- factor(zygomatic_data_context_only$Image_Set)

#Convert corrugator EMG responses to Z-scores
zygomatic_data_context_only <- zygomatic_data_context_only %>%
  mutate(response = scale(Response))
```

####Run the models:

```{r}
#Maximal model
zygomatic_data_context_only.model <- lmer(response ~ mutilation + threat + erotic + happy +
                            (1 + mutilation + threat + erotic + happy|participant) +
                            (1|item), data = zygomatic_data_context_only, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 1000000)))

```

```{r}
#Effect of category
zygomatic_data_context_only_no_category.model <- lmer(response ~ (1 + mutilation + threat + erotic + happy|participant) +
                            (1|item), data = zygomatic_data_context_only, REML = FALSE,
                          control = lmerControl('bobyqa', calc.derivs = FALSE, optCtrl = list(maxfun = 1000000)))

zygomatic_category_context_only_effect <- anova(zygomatic_data_context_only.model, zygomatic_data_context_only_no_category.model, type = 'chisq')
```

```{r}
#Category results
zygomatic_category_context_only_effect
```



