picture-perception.qmd

---
title: "Picture perception analysis"
format: html
---

```{r}
install.packages("emmeans")
install.packages("brms")
install.packages("tidybayes")

library(here)
library(tidyverse)
library(lme4)
library(emmeans)
library(brms)
library(tidybayes)
```

```{r}
d <- read_csv(here("data","kisumu vocab data.csv"))
```
# visualization

```{r}
ms <- d |>
  group_by(child, condition, age) |>
  summarise(accuracy = mean(accuracy))
```

Age distribution. 

```{r}
ggplot(ms, aes(x = age)) + 
  geom_histogram(binwidth = .5)
```



```{r}
ggplot(ms, aes(x = age, y = accuracy, col = condition)) + 
  geom_jitter(alpha = .5, height = .02, width = 0) + 
  geom_smooth(method = "lm") + 
  ylab("Accuracy") +
  xlab("Age (months)")
```

# Preregistered analyses

Let's set objects to the reference level. We'll also center age for interpretability. 

These models fail to converge pretty substantially with the maximal preregistered random effects structure. 

Also, good practice to center age to make the intercept more interpretable. (Helps with convergence too often).

```{r}
d$condition <- relevel(factor(d$condition), ref = "obj")
d$age_centered <- scale(d$age, scale = FALSE, center = TRUE)[,1]
```

First model. 

```{r}
mod <- glmer(accuracy ~ condition + (condition | child) + (condition | target_word),
             data = d, family = "binomial")
summary(mod)
```

We remove random effects following our general lab protocol. This results in crossed intercepts but no slopes. 

```{r}
mod_pruned <- glmer(accuracy ~ condition + (1 | child) + (1 | target_word), 
             data = d, family = "binomial")
summary(mod_pruned)
```

This shows pretty clearly that age is worse.

We fit the same preregistered model with age. 

```{r}
mod_age <- glmer(accuracy ~ condition * age_centered + 
                   (condition | child) + 
                   (condition * age_centered | target_word), 
              data = d, family = "binomial")
summary(mod_age)
```

Again we prune - random slopes are singular, but the intercepts converge. 

```{r}
mod_age_pruned <- glmer(accuracy ~ condition * age_centered + 
                   (1 | child) + 
                   (1 | target_word), 
              data = d, family = "binomial")
summary(mod_age_pruned)
```

You also get a pretty clear age interaction: BW is worse and grows more slowly. 

The other conditions are definitely better and don't appear to be distinguishable from object. 

# More viz

Because this is a within-subjects analysis, we can look at the effect for each kid by basic subtraction.

```{r}
ms <- ms |>
  group_by(child) |>
  mutate(effect = accuracy - accuracy[condition == "obj"])

ggplot(filter(ms, condition != "obj"),
       aes(x = age, y = effect)) + 
  geom_jitter(alpha = .3, height = .02, width = 0) + 
  geom_smooth(method = "lm") + 
  geom_hline(yintercept = 0, lty = 2) + 
  facet_wrap(~condition) + 
  ylab("Accuracy compared to object") +
  xlab("Age (months)")
```
This confirms our analysis above. 

Now running exploratory Bayesian analyses with full models - first in pre-registered model without age
```{r}
bayes_mod <- brm(formula = accuracy ~ condition + (condition | child) + (condition | target_word),
                 family=bernoulli(link="logit"), data = d, seed = 1)

joint_tests(bayes_mod)

#no effect of condition in model without age

```
Get posterior estimates + 95% credible intervals

```{r}
bayes.mod %>%
  as_draws_df() %>%
  select(starts_with("b_")) %>% 
  mutate(obj = b_Intercept, 
         bw = b_Intercept + b_conditionbw, 
     	   cartoon = b_Intercept + b_conditioncartoon,
         photo = b_Intercept + b_conditionphoto) %>%
  select(obj, photo, cartoon, bw) %>%
  mean_hdi() %>%
  pivot_longer(cols = -c(.width:.interval),
               names_to = "index",
               values_to = "value") %>% 
  select(index, value) %>% 
  mutate(index = ifelse(str_detect(index, fixed(".")), index, str_c(index, ".mean"))) %>% 
  separate(index, into = c("parameter", "type"), sep = "\\.") %>% 
  pivot_wider(names_from = type, 
              values_from = value) %>%
  mutate_at(vars(2:4), plogis)
```

Second, running exploratory Bayesian analyses with full model with age
```{r}
bayes_mod_age <- brm(formula = accuracy ~ condition * age_centered + 
                   (condition | child) + 
                   (condition * age_centered | target_word),
                 family=bernoulli(link="logit"), data = d, seed = 1)

joint_tests(bayes_mod_age)

#main effect of age, no effect of condition

#model term             df1 df2 F.ratio  Chisq p.value
#condition                3 Inf   1.536  4.608  0.2027
#age_centered             1 Inf  15.569 15.569  0.0001
#condition:age_centered   3 Inf   1.594  4.782  0.1886
```

Get posterior estimates + 95% credible intervals
```{r}
bayes_mod_age %>%
  as_draws_df() %>%
  select(starts_with("b_")) %>% 
  mutate(obj = b_Intercept, 
         bw = b_Intercept + b_conditionbw, 
     	   cartoon = b_Intercept + b_conditioncartoon,
         photo = b_Intercept + b_conditionphoto) %>%
  select(obj, photo, cartoon, bw) %>%
  mean_hdi() %>%
  pivot_longer(cols = -c(.width:.interval),
               names_to = "index",
               values_to = "value") %>% 
  select(index, value) %>% 
  mutate(index = ifelse(str_detect(index, fixed(".")), index, str_c(index, ".mean"))) %>% 
  separate(index, into = c("parameter", "type"), sep = "\\.") %>% 
  pivot_wider(names_from = type, 
              values_from = value) %>%
  mutate_at(vars(2:4), plogis)
```