---
title: "Imperative Benefit Evaluation"
author: "Anna ChromÃ¡"
date: "`r Sys.Date()`"
output: html_document
---
***
```{r echo=F, message=F, warning=F}
library(tidyverse)
library(lsr)
library(lmerTest)
library(ggeffects)
```

```{r echo=F, message=F, warning=F}
# import data
benefiteval <- read_csv("benefiteval_publish.csv")
```

```{r echo=F, message=F, warning=F}
# colors for figures
#colors <- c("#C00000","#29A1C7")
colorsgrey <- c("grey60","black")
```

### Overview over the Sample

```{r echo=F, message=F, warning=F}
benefiteval |> distinct(Participant, Gender) |> count(Gender)
benefiteval |> summarise(n_distinct(Participant), mean(Age), sd(Age))
```

```{r echo=F, message=F, warning=F}
## relevel the critical factors
benefiteval <- benefiteval |>
  mutate(Number = fct_relevel(Number, "sg", "pl"),
         Person.of.modal = fct_relevel(Person.of.modal, "1P modal", "2P modal"))
```

```{r echo=T, message=F, warning=F}
# mean benefit rating per item
mean.rating.per.item <- benefiteval |>
  group_by(Person.of.modal, Number, Item, ItemNo) |>
  summarise(mean.rated.benefit = round(mean(Rated.benefit), 1),
            lconf = round(ciMean(Rated.benefit, conf=0.95)[,1], 1),
            hconf = round(ciMean(Rated.benefit, conf=0.95)[,2], 1)) |>
  ungroup()

#library(openxlsx)
#write.xlsx(mean.rating.per.item, "benefit_rating_individual_items.xlsx")

# mean benefit rating per participant
mean.rating.per.part <- benefiteval |>
  group_by(Person.of.modal, Number, Participant) |>
  summarise(mean.rated.benefit = mean(Rated.benefit)) |>
  ungroup()

# mean benefit rating per group
mean.rating.per.group <- benefiteval |>
  group_by(Person.of.modal, Number) |>
  summarise(mean.rated.benefit = mean(Rated.benefit)) |>
  ungroup()
```

### Figure 2
Mean rating of the IMPERATIVE BENEFIT according to PERSON OF THE MODAL VERB and NUMBER. Small dots are mean ratings for individual items, big dots are group means.

```{r fig.align='center', out.width="70%", echo=F, message=F, warning=F}
#Fig2_mean.rating.per.group.item <-
mean.rating.per.group |>
  ggplot(aes(Number, mean.rated.benefit, color = Person.of.modal)) +
  scale_color_manual(values = colorsgrey) +
  theme_bw() +
  theme(legend.position = "top",
        axis.text = element_text(size = 14),
        axis.title = element_text(size = 14)) +
  coord_cartesian(ylim = c(0,100)) +
  geom_point(data = mean.rating.per.item,
             position = position_dodge(width = 0.2), alpha = 0.2) +
  geom_point(position = position_dodge(width = 0.2), size = 3) +
  geom_hline(yintercept = 49.5, linetype = "dashed") +
  annotate("text", x = 0, y = 0, color = "grey50", size = 4, hjust = 0, vjust = 0,
           label = "speaker") +
  annotate("text", x = 0, y = 100, color = "grey50", size = 4, hjust = 0, vjust = 1,
           label = "addressee") +
  labs(x = NULL,
       y = "Mean Rating of Imperative Benefit",
       color = NULL)

#ggsave("Fig2_mean.rating.per.group.item.png",
#       Fig2_mean.rating.per.group.item, device='png', units = "px", width = 5000, height = 5000, dpi=1200)
```

### Figure 3
Mean rating of the IMPERATIVE BENEFIT according to PERSON OF THE MODAL VERB and NUMBER. Data points are mean ratings from individual participants.

```{r fig.align='center', out.width="70%", echo=F, message=F, warning=F}
#Fig3_mean.rating.per.part <-
mean.rating.per.part %>%
  ggplot(aes(Number, mean.rated.benefit, color = Person.of.modal)) +
  scale_color_manual(values = colorsgrey) +
  theme_bw() +
  theme(legend.position = "top",
        axis.text = element_text(size = 14),
        axis.title = element_text(size = 14)) +
  coord_cartesian(ylim = c(0,100)) +
  geom_point(position = position_dodge(width = 0.2), alpha = 0.2) +
  geom_hline(yintercept = 49.5, linetype = "dashed") +
  annotate("text", x = 0, y = 0, color = "grey50", size = 4, hjust = 0, vjust = 0,
           label = "speaker") +
  annotate("text", x = 0, y = 100, color = "grey50", size = 4, hjust = 0, vjust = 1,
           label = "addressee") +
  labs(x = NULL,
       y = "Mean Rating of Imperative Benefit",
       color = NULL)

#ggsave("Fig3_mean.rating.per.part.png",
#       Fig3_mean.rating.per.part, device='png', units = "px", width = 5000, height = 5000, dpi=1200)
```

### Mixed-effect model 1
The primary model of benefit rating.

CONTRASTS
```{r echo=F, message=F, warning=F}
contrasts(benefiteval$Person.of.modal)
contrasts(benefiteval$Number)
```

```{r echo=T, message=F, warning=F}
# linear model with random intercepts and by-item slopes
# without main effect of Number
m1 <- lmer(Rated.benefit ~ Person.of.modal + Person.of.modal:Number + (1|Participant) +                                                           (1+Person.of.modal+Number|Item),
           data = benefiteval, control = lmerControl(optimizer="bobyqa"))
```

```{r echo=F, message=F, warning=F}
summary(m1)
```

### Mixed-effect model 2
Model of the magnitude of the preference for each of the conditions.

```{r echo=T, message=F, warning=F}
# transform the scale into two scales for speaker and addressee benefit separately
benefiteval.two.scales <-
  benefiteval |>
  mutate(spe.likelihood = abs(Rated.benefit-99),
         add.likelihood = Rated.benefit) |>
  select(-Rated.benefit) |>
  pivot_longer(cols = c(spe.likelihood, add.likelihood),
               names_to = "benefit",
               values_to = "likelihood") |>
  mutate(benefit = fct_recode(benefit, Speaker = "spe.likelihood", Addressee = "add.likelihood"),
         benefit = fct_relevel(benefit, c("Speaker", "Addressee"))) |>
  filter((benefit == "Speaker"   & Person.of.modal == "1P modal") |
         (benefit == "Addressee" & Person.of.modal == "2P modal"))
```

```{r echo=T, message=F, warning=F}
m2 <- lmer(likelihood ~ benefit + Number:benefit + (1|Participant) + (1+benefit|Item),
                     data = benefiteval.two.scales, control = lmerControl(optimizer="bobyqa"))
```

```{r echo=F, message=F, warning=F}
summary(m2)
```

### Figure 4
Predictions from model 2<br>
Model-predicted values of the two-scale transformed IMPERATIVE BENEFIT with 95% CI for benefiting speaker in 1P modal vs addressee in 2P modal condition with regard to NUMBER.

```{r fig.align='center', out.width="70%", echo=F, message=F, warning=F}
#Fig4_predicted.two.scales <- 
plot(
  ggpredict(m2, terms = c("Number","benefit"), type = "fixed", ci.level = 0.95)
) +
  theme_bw() +
  theme(legend.position = "top",
        axis.text = element_text(size = 14),
        axis.title = element_text(size = 14)) +
  coord_cartesian(ylim = c(0,100)) +
  geom_hline(yintercept = 50, linetype = "dashed") +
  scale_color_manual(values = colorsgrey) +
  annotate("text", x = 0, y = 0, color = "grey50", size = 4, hjust = 0, vjust = 0,
           label = "speaker") +
  annotate("text", x = 0, y = 100, color = "grey50", size = 4, hjust = 0, vjust = 1,
           label = "addressee") +
  labs(title = NULL,
       x = NULL,
       y = "Predicted Magnitude of Rating",
       color = NULL)

#ggsave("Fig4_predicted.two.scales.png",
#       Fig4_predicted.two.scales, device='png', units = "px", width = 5000, height = 5000, dpi=1200)
```

## Individual Items in Detail

### Figure 5
Mean ratings with 95% CI of IMPERATIVE BENEFIT in individual test items according to PERSON OF THE MODAL VERB and NUMBER.

```{r fig.align='center', fig.height=10, out.width="100%", echo=F, message=F, warning=F}
#Fig5_mean.rating.per.item <-
mean.rating.per.item |>
  ggplot(aes(mean.rated.benefit, Item, color = Person.of.modal)) +
  facet_wrap(. ~ Number) +
  scale_color_manual(values = colorsgrey) +
  geom_vline(xintercept = 50, linetype = "dashed") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text = element_text(size = 14),
        axis.title = element_text(size = 14),
        legend.text = element_text(size = 14)) +
  geom_point(position = position_dodge(width = 0.3)) +
  geom_errorbar(aes(xmin = lconf, xmax = hconf), width = 0.3, position = position_dodge(width = 0.3)) +
  annotate("text", x = 0, y = 1, color = "grey50", angle = 90, size = 5, hjust = 0, vjust = 1,
           label = "speaker benefit") +
  annotate("text", x = 100, y = 1, color = "grey50", angle = 90, size = 5, hjust = 0, vjust = 0,
           label = "addressee benefit") +
  labs(y = NULL,
       x = "Mean Rating of Imperative Benefit",
       color = NULL)
  
#ggsave("Fig5_mean.rating.per.item.png",
#       Fig5_mean.rating.per.item, device='png', units = "px", width = 10000, height = 7500, dpi=1200)
```

```{r echo=F, message=F, warning=F}
# calculate mean rating difference per item (rating for 2P modal minus 1P modal)
mean.diff.per.item <- mean.rating.per.item %>%
  select(-lconf, -hconf) %>%
  pivot_wider(names_from = Person.of.modal, values_from = mean.rated.benefit) %>%
  mutate(diff = `2P modal`-`1P modal`)
```

### Figure 6
Difference between 2P modal and 1P modal rating according to test item and NUMBER.

```{r fig.align='center', fig.height=10, out.width="100%", echo=F, message=F, warning=F}
#Fig6_mean.diff.per.item <-
mean.diff.per.item %>%
  ggplot(aes(diff,reorder(Item, diff))) +
  facet_grid(cols = vars(Number)) +
  theme_bw() +
  theme(axis.text = element_text(size = 14),
        axis.title = element_text(size = 14)) +
  geom_point() +
  annotate("text", x = 0, y = "22 nech ... tÃ¡c", color = "grey50", angle = 90, size = 5, hjust = 1, vjust = 1,
           label = "no difference") +
  annotate("text", x = 100, y = "22 nech ... tÃ¡c", color = "grey50", angle = 90, size = 5, hjust = 1, vjust = 0,
           label = "maximum expected difference") +
  labs(x = "Difference in Mean Rating",
       y = NULL)

#ggsave("Fig6_mean.diff.per.item.png",
#       Fig6_mean.diff.per.item, device='png', units = "px", width = 10000, height = 7500, dpi=1200)
```