1 Opening Files

library(brms)
library(readxl)
library(writexl)
library(rmarkdown)
library(tidybayes)
library(DHARMa)
library(ggplot2)
library(tidyr)
library(dplyr)
library(GGally)
library(ggpubr)
library(tidybayes)
library(patchwork)
library(sjPlot)
library(lubridate)
library(zoo)
library(readODS)
DOD_54mo <- read_ods("DOD_54mo_blind.ods") 
N_total <- nrow(DOD_54mo) # sanity check
DOD_42mo <- read_ods("DOD_42mo_prepped.ods") %>%
  group_by(Blind_ID) %>%
  slice(1) %>%
  dplyr::select("Blind_ID", "Prod_DOD") %>%
  rename("DOD_42" = "Prod_DOD") %>%
  ungroup()
DOD_48mo <- read_ods("DOD_48mo_prepped.ods") %>%
  group_by(Blind_ID) %>%
  slice(1) %>%
  dplyr::select("Blind_ID", "Prod_DOD") %>%
  rename("DOD_48" = "Prod_DOD") %>%
  ungroup()

Create mismatch variables. I’m also creating some character versions of many factors to make plotting a bit easier later.

DOD_54mo <- DOD_54mo %>% 
 mutate(
    Prime = ifelse(Prime.Type==.5, yes="DOD", no="POD"), # Prime type, character version
    Overlap = ifelse(Verb.match==.5, yes="Overlap", no="No Overlap") # overlap
    )

DOD_54mo$bias <-"DOD" # set as DOD bias
DOD_54mo$bias[DOD_54mo$Prime.verb=="passed"] = "POD" # Change POD-biased prime verbs. 
DOD_54mo$bias[DOD_54mo$Prime.verb=="sent"] = "POD"
DOD_54mo$bias[DOD_54mo$Prime.verb=="threw"] = "POD"

DOD_54mo$biasmatch <-"match" # set to match

DOD_54mo$biasmatch[DOD_54mo$bias=="POD" & DOD_54mo$Prime == "DOD"] = "mismatch" # Change to mismatch. 
DOD_54mo$biasmatch[DOD_54mo$bias=="DOD" & DOD_54mo$Prime == "POD"] = "mismatch"

table(DOD_54mo$biasmatch)
## 
##    match mismatch 
##      529      527
DOD_54mo$biasmismatchc <- ifelse(DOD_54mo$biasmatch == "mismatch", yes=.5, no=-.5)

2 Missing data

There is one type of trial coded as an admin error that I don’t want to remove at the time being – trials 12r on lists 6 and 8.

DOD_54mo <- DOD_54mo %>%
  mutate(
    Trial_Drop = ifelse(is.na(RESPONSE.CODE), yes=1, no=0), 
    Trial_Problem = ifelse(List==6 & Item == "12r", yes = 1, no =ifelse(
      List==8 & Item == "12r", yes = 1, no =0
    ))
  )

nrow(filter(DOD_54mo, Trial_Problem==1)) == (nrow(filter(DOD_54mo, List==6))/12) + (nrow(filter(DOD_54mo, List==8))/12)
## [1] TRUE
N_miss <- DOD_54mo %>%
  filter(is.na(RESPONSE.CODE)) %>%
  nrow()

See how many missing values we have by participant

DOD_54mo %>%
  group_by(Blind_ID) %>%
  mutate(
    N_Missing = sum(is.na(RESPONSE.CODE))
  ) %>%
  slice(1) %>%
  select("Blind_ID", "N_Missing") %>%
  group_by(N_Missing) %>%
  count()
## # A tibble: 5 × 2
## # Groups:   N_Missing [5]
##   N_Missing     n
##       <int> <int>
## 1         0    67
## 2         1    12
## 3         2     3
## 4         3     5
## 5        10     1
DOD_54mo <- DOD_54mo %>%
  group_by(Blind_ID) %>%
  mutate(
    N_Missing = sum(is.na(RESPONSE.CODE))
  ) %>%
  slice(1) %>%
  select("Blind_ID", "N_Missing") %>%
  mutate(
    Drop = ifelse(N_Missing >= 10, 1, 0), 
    MightDrop = ifelse(N_Missing > 6, 1, 0)
  ) %>%
  full_join(DOD_54mo, .)
## Joining with `by = join_by(Blind_ID)`
N_drop = DOD_54mo %>%
  filter(Drop == 1 & !is.na(RESPONSE.CODE)) %>%
  count()

One participant is missing 10 trials.

And by item.

DOD_54mo %>%
  group_by(Target.verb) %>%
  mutate(
    N_Missing = sum(is.na(RESPONSE.CODE))
  ) %>%
  slice(1) %>%
  select("Target.verb", "N_Missing") %>%
  group_by(Target.verb) 
## # A tibble: 7 × 2
## # Groups:   Target.verb [7]
##   Target.verb N_Missing
##   <chr>           <int>
## 1 bought              2
## 2 brought             9
## 3 gave                9
## 4 passed              1
## 5 sent                3
## 6 showed              6
## 7 threw              13

All non-target verbs have been removed. It doesn’t make sense to estimate random effects for a level of a random factor with only a single observation.

And let’s just check the proportion of missing trials across the various experimental conditions.

DOD_54mo %>%
  mutate(
    Missing = ifelse(is.na(RESPONSE.CODE), yes=1, no = 0)
  ) %>%
  group_by(Prime, Overlap) %>%
  summarise(
    mean = mean(Missing), 
    sd = sd(Missing)
  )
## `summarise()` has grouped output by 'Prime'. You can override using the
## `.groups` argument.
## # A tibble: 4 × 4
## # Groups:   Prime [2]
##   Prime Overlap      mean    sd
##   <chr> <chr>       <dbl> <dbl>
## 1 DOD   No Overlap 0.0478 0.214
## 2 DOD   Overlap    0.0288 0.167
## 3 POD   No Overlap 0.0586 0.235
## 4 POD   Overlap    0.0295 0.170
DOD_54mo %>%
  mutate(
    Missing = ifelse(is.na(RESPONSE.CODE), yes=1, no = 0)
  ) %>%
  group_by(Prime, biasmatch) %>%
  summarise(
    mean = mean(Missing), 
    sd = sd(Missing)
  )
## `summarise()` has grouped output by 'Prime'. You can override using the
## `.groups` argument.
## # A tibble: 4 × 4
## # Groups:   Prime [2]
##   Prime biasmatch   mean    sd
##   <chr> <chr>      <dbl> <dbl>
## 1 DOD   match     0.0377 0.191
## 2 DOD   mismatch  0.0379 0.191
## 3 POD   match     0.0379 0.191
## 4 POD   mismatch  0.0494 0.217

The differences in missingness between condition don’t look drastic. To sum up things, we have a fair amount of missing data, including 1 participants who are missing 10 or more trials. I’m goig to remove the participant who is missing 10 trials.

DOD_54mo <- filter(DOD_54mo, !is.na(RESPONSE.CODE) & Drop==0) # 1011

nrow(DOD_54mo) == N_total - N_miss - N_drop
##         n
## [1,] TRUE
DOD_54mo %>%
  summarise(
    min = min(Session.Age),
    max = max(Session.Age)
  )
##       min    max
## 1 54 : 10 55 : 5
(p1 <- DOD_54mo %>% 
  group_by(Prime, Overlap) %>%
  summarise(mean = mean(RESPONSE.CODE, na.rm=TRUE), 
            sd = sd(RESPONSE.CODE, na.rm=TRUE)) %>%
  
ggplot(aes(y=mean, x=Overlap, fill=Prime)) +
  geom_bar(stat="identity",position="dodge") +
 geom_errorbar(aes(ymin=mean, ymax=mean+sd), width=.2,
                 position=position_dodge(.9)) +
  ylim(0, .6) + 
  theme_minimal() + 
  labs(y="Raw Mean") + 
  theme(legend.position = "none")
)
## `summarise()` has grouped output by 'Prime'. You can override using the
## `.groups` argument.

DOD_54mo %>% 
  group_by(Prime, Overlap, Target.verb) %>%
  summarise(mean = mean(RESPONSE.CODE, na.rm=TRUE), 
            sd = sd(RESPONSE.CODE, na.rm=TRUE)) %>%
  
ggplot(aes(y=mean, x=Overlap, fill=Prime)) +
  geom_bar(stat="identity",position="dodge") +
 geom_errorbar(aes(ymin=mean, ymax=mean+sd), width=.2,
                 position=position_dodge(.9)) +
  facet_wrap(~Target.verb) 
## `summarise()` has grouped output by 'Prime', 'Overlap'. You can override using
## the `.groups` argument.

This is totally different than the other time points – it looks like there is a lexical boost for some verbs (showed, brought and threw, but not others (gave, passed or sent).

And let’s look at the pattern of responses when we only consider participants who produced at least 1 DOD.

DOD_54mo <- DOD_54mo %>%
  group_by(Blind_ID) %>%
  mutate(
    Prod_DOD = sum(RESPONSE.CODE)
  ) %>%
  dplyr::select(Blind_ID, Prod_DOD) %>%
  slice(1) %>%
  ungroup() %>%
  full_join(DOD_54mo, ., by="Blind_ID") %>%
  left_join(., DOD_42mo, by="Blind_ID") %>%
  left_join(., DOD_48mo, by="Blind_ID")

write_xlsx(DOD_54mo, "DOD_54mo_prepped.xlsx")
DOD_54mo2 <- filter(DOD_54mo, !(Prod_DOD==0 & DOD_48==0 & DOD_42==0))

DOD_54mo2 %>%
  group_by(Blind_ID) %>%
  slice(1) %>%
  ungroup() %>%
  count()
## # A tibble: 1 × 1
##       n
##   <int>
## 1    75
(p1b <- DOD_54mo2 %>% 
  group_by(Prime, Overlap) %>%
  summarise(mean = mean(RESPONSE.CODE, na.rm=TRUE), 
            sd = sd(RESPONSE.CODE, na.rm=TRUE)) %>%
  
ggplot(aes(y=mean, x=Overlap, fill=Prime)) +
  geom_bar(stat="identity",position="dodge") +
 geom_errorbar(aes(ymin=mean, ymax=mean+sd), width=.2,
                 position=position_dodge(.9)) +
  ylim(0, .9) + 
  theme_minimal() + 
  labs(y="Raw Mean") + 
  theme(legend.position = "none")
)
## `summarise()` has grouped output by 'Prime'. You can override using the
## `.groups` argument.

DOD_54mo2 %>% 
  group_by(Prime, Overlap, Target.verb) %>%
  summarise(mean = mean(RESPONSE.CODE, na.rm=TRUE), 
            sd = sd(RESPONSE.CODE, na.rm=TRUE)) %>%
  
ggplot(aes(y=mean, x=Overlap, fill=Prime)) +
  geom_bar(stat="identity",position="dodge") +
 geom_errorbar(aes(ymin=mean, ymax=mean+sd), width=.2,
                 position=position_dodge(.9)) +
  facet_wrap(~Target.verb) 
## `summarise()` has grouped output by 'Prime', 'Overlap'. You can override using
## the `.groups` argument.

Plot proprortion of DOD responses by each cell in two-way interaction between biasmatch and prime type.

(p1c <- DOD_54mo %>% 
  filter(Overlap == "No Overlap") %>%
  group_by(Prime, biasmatch) %>%
  summarise(mean = mean(RESPONSE.CODE, na.rm=TRUE), 
            sd = sd(RESPONSE.CODE, na.rm=TRUE)) %>%
  
ggplot(aes(y=mean, x=biasmatch, fill=Prime)) +
  geom_bar(stat="identity",position="dodge") +
 geom_errorbar(aes(ymin=mean, ymax=mean+sd), width=.2,
                 position=position_dodge(.75)) +
  ylim(0, .9) + 
  theme_minimal() + 
  labs(y="Raw Mean")
)
## `summarise()` has grouped output by 'Prime'. You can override using the
## `.groups` argument.

cbPalette <- c("#000000",  "#56B4E9", "#E69F00","#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
DOD_54mo %>% 
  filter(Overlap == "No Overlap") %>%
  group_by(Prime, biasmatch, Target.verb) %>%
  summarise(mean = mean(RESPONSE.CODE, na.rm=TRUE), 
            sd = sd(RESPONSE.CODE, na.rm=TRUE),
            total=n()) %>%

ggplot(aes(y=mean, x=biasmatch, fill=Prime, label=total)) +
  geom_bar(stat="identity",position="dodge") +
 geom_errorbar(aes(ymin=mean, ymax=mean+sd), width=.2,
                 position=position_dodge(.9)) + ylab("Prop DOD") +
  facet_wrap(~Target.verb) + scale_fill_manual(values=cbPalette)+ theme_minimal()
## `summarise()` has grouped output by 'Prime', 'biasmatch'. You can override
## using the `.groups` argument.

ggsave("plots/D1.png", last_plot(), dpi=500)
## Saving 7 x 5 in image
DOD_54mo %>% 
  filter(Overlap == "No Overlap") %>%
  group_by(Target.verb, Prime, biasmatch) %>%
  count()
## # A tibble: 22 × 4
## # Groups:   Target.verb, Prime, biasmatch [22]
##    Target.verb Prime biasmatch     n
##    <chr>       <chr> <chr>     <int>
##  1 brought     DOD   match        19
##  2 brought     DOD   mismatch     18
##  3 brought     POD   match        23
##  4 brought     POD   mismatch     18
##  5 gave        DOD   match         2
##  6 gave        DOD   mismatch     46
##  7 gave        POD   match        41
##  8 gave        POD   mismatch      1
##  9 passed      DOD   match        11
## 10 passed      DOD   mismatch     20
## # ℹ 12 more rows

Again, it looks like the pattern is different on the levels of aggregate data and indiviaul verbs, which was the case with the 42 month verbs.

DOD_54mo %>%
   filter(Overlap == "No Overlap") %>% 
   filter(Target.verb=="gave" & Prime=="DOD" & biasmatch == "match")
##   Blind_ID Zoom Gender Session.Age List Item Prime.verb Target.verb Prime.Type
## 1 9hKLS9my    1      1     54 : 22    3    2     showed        gave        0.5
## 2 LkVZQCTn    1      1     54 : 15    7    2     showed        gave        0.5
##   Verb.match RESPONSE.CODE COVID.Zoom.session Prime    Overlap bias biasmatch
## 1       -0.5             1                  1   DOD No Overlap  DOD     match
## 2       -0.5             1                 NA   DOD No Overlap  DOD     match
##   biasmismatchc Trial_Drop Trial_Problem N_Missing Drop MightDrop Prod_DOD
## 1          -0.5          0             0         3    0         0        1
## 2          -0.5          0             0         1    0         0        1
##   DOD_42 DOD_48
## 1      3      0
## 2      3      0
DOD_54mo %>%
   filter(Overlap == "No Overlap") %>% 
   filter(Target.verb=="gave" & Prime=="POD" & biasmatch == "mismatch")
##   Blind_ID Zoom Gender Session.Age List Item Prime.verb Target.verb Prime.Type
## 1 XncBwsN2   NA      1     54 : 22    8  12r       gave        gave       -0.5
##   Verb.match RESPONSE.CODE COVID.Zoom.session Prime    Overlap bias biasmatch
## 1       -0.5             0                 NA   POD No Overlap  DOD  mismatch
##   biasmismatchc Trial_Drop Trial_Problem N_Missing Drop MightDrop Prod_DOD
## 1           0.5          0             1         0    0         0        2
##   DOD_42 DOD_48
## 1     NA      0

3 Lexical Boost

3.1 All Participants

logit_prior = prior(normal(0,2), class=sd)

First lets consider whether we should include random effects by the interaction between ID and Item. In principle this is possible, but I think it will greatly complicated the model. So I’ll fit two variance component models with the two random effects specifications and compare them to one another.

m54mo0<- brm(RESPONSE.CODE ~ 1 + (1 | Blind_ID) + (1 | Target.verb),  family="bernoulli", prior=logit_prior, seed=12345, iter=4000, cores=4, data=DOD_54mo, file="54mo/output/m54mo0.rds")
m54mo0b <- brm(RESPONSE.CODE ~ 1 + (1 | Blind_ID) + (1 | Target.verb) + (1 | Blind_ID:Target.verb),  seed=12345,  family="bernoulli", prior=logit_prior, iter=4000, cores=4, data=DOD_54mo, file="54mo/output/m54mo0b.rds")

Add loo for loo_compare

#loo1 <- loo(m54mo0b)
#loo2 <- reloo(m54mo0b, loo1)
#loo3 <- loo(m54mo0)
#loo_compare(loo2, loo3)

Here’s the model of the lexical boost.

m54mo1 <- brm(RESPONSE.CODE ~ 1 + Prime.Type*Overlap  +  (1 +  Prime.Type*Overlap | Blind_ID) + (1 +  Prime.Type*Overlap   | Target.verb), logit_prior,family="bernoulli", cores=4, seed=12345, control=list(adapt_delta=.9),  data=DOD_54mo, file="54mo/output/m54mo1.rds")

m54mo1 <- add_criterion(m54mo1, "loo")

Chains look okay.

#plot(m54mo1, N=4)
summary(m54mo1, prob=.9)
##  Family: bernoulli 
##   Links: mu = logit 
## Formula: RESPONSE.CODE ~ 1 + Prime.Type * Overlap + (1 + Prime.Type * Overlap | Blind_ID) + (1 + Prime.Type * Overlap | Target.verb) 
##    Data: DOD_54mo (Number of observations: 1011) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Group-Level Effects: 
## ~Blind_ID (Number of levels: 87) 
##                                               Estimate Est.Error l-90% CI
## sd(Intercept)                                     1.69      0.32     1.20
## sd(Prime.Type)                                    0.71      0.48     0.06
## sd(OverlapOverlap)                                0.62      0.43     0.06
## sd(Prime.Type:OverlapOverlap)                     1.04      0.69     0.12
## cor(Intercept,Prime.Type)                         0.15      0.40    -0.56
## cor(Intercept,OverlapOverlap)                     0.28      0.41    -0.48
## cor(Prime.Type,OverlapOverlap)                    0.10      0.44    -0.65
## cor(Intercept,Prime.Type:OverlapOverlap)          0.14      0.41    -0.58
## cor(Prime.Type,Prime.Type:OverlapOverlap)        -0.09      0.44    -0.77
## cor(OverlapOverlap,Prime.Type:OverlapOverlap)     0.04      0.45    -0.71
##                                               u-90% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)                                     2.23 1.00     1723     2802
## sd(Prime.Type)                                    1.57 1.00     1366     1999
## sd(OverlapOverlap)                                1.42 1.00     1111     2123
## sd(Prime.Type:OverlapOverlap)                     2.29 1.00     1362     1950
## cor(Intercept,Prime.Type)                         0.76 1.00     5057     3135
## cor(Intercept,OverlapOverlap)                     0.85 1.00     4252     3166
## cor(Prime.Type,OverlapOverlap)                    0.78 1.00     2112     2659
## cor(Intercept,Prime.Type:OverlapOverlap)          0.76 1.00     5092     2712
## cor(Prime.Type,Prime.Type:OverlapOverlap)         0.66 1.00     2978     3276
## cor(OverlapOverlap,Prime.Type:OverlapOverlap)     0.76 1.00     2561     3161
## 
## ~Target.verb (Number of levels: 6) 
##                                               Estimate Est.Error l-90% CI
## sd(Intercept)                                     1.81      0.62     1.02
## sd(Prime.Type)                                    0.58      0.51     0.04
## sd(OverlapOverlap)                                0.88      0.55     0.13
## sd(Prime.Type:OverlapOverlap)                     1.47      1.02     0.14
## cor(Intercept,Prime.Type)                         0.00      0.44    -0.71
## cor(Intercept,OverlapOverlap)                    -0.36      0.38    -0.88
## cor(Prime.Type,OverlapOverlap)                    0.01      0.44    -0.71
## cor(Intercept,Prime.Type:OverlapOverlap)         -0.09      0.39    -0.69
## cor(Prime.Type,Prime.Type:OverlapOverlap)        -0.06      0.45    -0.76
## cor(OverlapOverlap,Prime.Type:OverlapOverlap)     0.06      0.43    -0.66
##                                               u-90% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)                                     2.99 1.00     2038     2561
## sd(Prime.Type)                                    1.54 1.00     2662     2029
## sd(OverlapOverlap)                                1.88 1.00     1927     1998
## sd(Prime.Type:OverlapOverlap)                     3.35 1.00     1871     2606
## cor(Intercept,Prime.Type)                         0.72 1.00     6343     2947
## cor(Intercept,OverlapOverlap)                     0.35 1.00     4546     3046
## cor(Prime.Type,OverlapOverlap)                    0.74 1.00     3684     3281
## cor(Intercept,Prime.Type:OverlapOverlap)          0.60 1.00     5935     3263
## cor(Prime.Type,Prime.Type:OverlapOverlap)         0.69 1.00     2913     2794
## cor(OverlapOverlap,Prime.Type:OverlapOverlap)     0.74 1.00     3473     3580
## 
## Population-Level Effects: 
##                           Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS
## Intercept                    -3.19      0.84    -4.51    -1.78 1.00     1739
## Prime.Type                    0.84      0.57    -0.06     1.76 1.00     4098
## OverlapOverlap               -0.56      0.66    -1.62     0.50 1.00     2148
## Prime.Type:OverlapOverlap     0.72      1.07    -0.97     2.54 1.00     3476
##                           Tail_ESS
## Intercept                     2536
## Prime.Type                    3202
## OverlapOverlap                2331
## Prime.Type:OverlapOverlap     2751
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
hypothesis(m54mo1, "Prime.Type > 0")
## Hypothesis Tests for class b:
##         Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
## 1 (Prime.Type) > 0     0.84      0.57    -0.06     1.76      14.62      0.94
##   Star
## 1     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
hypothesis(m54mo1, "Prime.Type:OverlapOverlap > 0")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (Prime.Type:Overl... > 0     0.72      1.07    -0.97     2.54       3.28
##   Post.Prob Star
## 1      0.77     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
hypothesis(m54mo1, "Intercept + OverlapOverlap + .5*Prime.Type + .5*Prime.Type:OverlapOverlap > 
           Intercept + OverlapOverlap - .5*Prime.Type - .5*Prime.Type:OverlapOverlap ")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (Intercept+Overla... > 0     1.56      1.05    -0.08     3.32      15.39
##   Post.Prob Star
## 1      0.94     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
m54mo1.model.check <- createDHARMa(
  simulatedResponse = t(posterior_predict(m54mo1)), 
  observedResponse = DOD_54mo$RESPONSE.CODE, 
  fittedPredictedResponse = apply(t(posterior_epred(m54mo1)), 1, mean), 
  integerResponse=TRUE
)

plot(m54mo1.model.check) # standard plots

testDispersion(m54mo1.model.check)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.73983, p-value = 0.0045
## alternative hypothesis: two.sided

Sensitivity analyses (without zoom sessions, without problem trials and without trials with high Pareto’s k) in R Markdown file.

3.2 Grammatical Only

m54mo2 <- brm(RESPONSE.CODE ~ 1 + Prime.Type*Overlap  +  (1 +  Prime.Type*Overlap | Blind_ID) + (1 +  Prime.Type*Overlap   | Target.verb), logit_prior,family="bernoulli", cores=4, seed=12345, control=list(adapt_delta=.9), data=DOD_54mo2, file="54mo/output/m54mo2.rds")
m54mo2 <- add_criterion(m54mo2, "loo")

summary(m54mo2, prob=.9)
##  Family: bernoulli 
##   Links: mu = logit 
## Formula: RESPONSE.CODE ~ 1 + Prime.Type * Overlap + (1 + Prime.Type * Overlap | Blind_ID) + (1 + Prime.Type * Overlap | Target.verb) 
##    Data: DOD_54mo2 (Number of observations: 872) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Group-Level Effects: 
## ~Blind_ID (Number of levels: 75) 
##                                               Estimate Est.Error l-90% CI
## sd(Intercept)                                     1.51      0.31     1.03
## sd(Prime.Type)                                    0.72      0.49     0.07
## sd(OverlapOverlap)                                0.61      0.42     0.06
## sd(Prime.Type:OverlapOverlap)                     1.00      0.68     0.09
## cor(Intercept,Prime.Type)                         0.14      0.41    -0.57
## cor(Intercept,OverlapOverlap)                     0.29      0.40    -0.44
## cor(Prime.Type,OverlapOverlap)                    0.11      0.44    -0.66
## cor(Intercept,Prime.Type:OverlapOverlap)          0.12      0.41    -0.60
## cor(Prime.Type,Prime.Type:OverlapOverlap)        -0.10      0.46    -0.79
## cor(OverlapOverlap,Prime.Type:OverlapOverlap)     0.03      0.44    -0.69
##                                               u-90% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)                                     2.06 1.00     1735     2613
## sd(Prime.Type)                                    1.62 1.00     1366     1845
## sd(OverlapOverlap)                                1.36 1.00     1050     2128
## sd(Prime.Type:OverlapOverlap)                     2.26 1.00     1388     2091
## cor(Intercept,Prime.Type)                         0.77 1.00     4482     2596
## cor(Intercept,OverlapOverlap)                     0.85 1.00     3763     3080
## cor(Prime.Type,OverlapOverlap)                    0.78 1.00     1632     2191
## cor(Intercept,Prime.Type:OverlapOverlap)          0.76 1.00     5009     2700
## cor(Prime.Type,Prime.Type:OverlapOverlap)         0.70 1.00     2872     3136
## cor(OverlapOverlap,Prime.Type:OverlapOverlap)     0.75 1.00     2851     3106
## 
## ~Target.verb (Number of levels: 6) 
##                                               Estimate Est.Error l-90% CI
## sd(Intercept)                                     1.82      0.62     1.02
## sd(Prime.Type)                                    0.57      0.49     0.04
## sd(OverlapOverlap)                                0.89      0.56     0.14
## sd(Prime.Type:OverlapOverlap)                     1.41      1.03     0.12
## cor(Intercept,Prime.Type)                         0.00      0.44    -0.72
## cor(Intercept,OverlapOverlap)                    -0.36      0.39    -0.88
## cor(Prime.Type,OverlapOverlap)                    0.02      0.45    -0.71
## cor(Intercept,Prime.Type:OverlapOverlap)         -0.09      0.41    -0.73
## cor(Prime.Type,Prime.Type:OverlapOverlap)        -0.05      0.45    -0.77
## cor(OverlapOverlap,Prime.Type:OverlapOverlap)     0.05      0.43    -0.67
##                                               u-90% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)                                     3.01 1.00     2678     3006
## sd(Prime.Type)                                    1.50 1.00     2628     2497
## sd(OverlapOverlap)                                1.91 1.00     1576     1410
## sd(Prime.Type:OverlapOverlap)                     3.41 1.00     1473     2318
## cor(Intercept,Prime.Type)                         0.72 1.00     6932     3133
## cor(Intercept,OverlapOverlap)                     0.36 1.00     4358     3070
## cor(Prime.Type,OverlapOverlap)                    0.74 1.00     3174     3233
## cor(Intercept,Prime.Type:OverlapOverlap)          0.61 1.00     5257     2815
## cor(Prime.Type,Prime.Type:OverlapOverlap)         0.69 1.00     2683     3104
## cor(OverlapOverlap,Prime.Type:OverlapOverlap)     0.74 1.00     3221     3459
## 
## Population-Level Effects: 
##                           Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS
## Intercept                    -2.80      0.83    -4.16    -1.43 1.00     1456
## Prime.Type                    0.84      0.55    -0.01     1.72 1.00     3608
## OverlapOverlap               -0.53      0.63    -1.54     0.52 1.00     2501
## Prime.Type:OverlapOverlap     0.80      1.05    -0.76     2.57 1.00     2816
##                           Tail_ESS
## Intercept                     1959
## Prime.Type                    2843
## OverlapOverlap                2744
## Prime.Type:OverlapOverlap     2150
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
hypothesis(m54mo2, "Prime.Type > 0", seed=2)
## Hypothesis Tests for class b:
##         Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
## 1 (Prime.Type) > 0     0.84      0.55    -0.01     1.72       17.6      0.95
##   Star
## 1     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
hypothesis(m54mo2, "Prime.Type:OverlapOverlap > 0")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (Prime.Type:Overl... > 0      0.8      1.05    -0.76     2.57       4.12
##   Post.Prob Star
## 1       0.8     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
hypothesis(m54mo2, "Intercept + OverlapOverlap + .5*Prime.Type + .5*Prime.Type:OverlapOverlap > 
           Intercept + OverlapOverlap - .5*Prime.Type - .5*Prime.Type:OverlapOverlap ")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (Intercept+Overla... > 0     1.64         1     0.11      3.4      23.24
##   Post.Prob Star
## 1      0.96    *
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Looks like the posterior probability for the priming effect is rounded up to 95. Get the exact value

posterior_samples(m54mo2) %>%
  dplyr::select(b_Prime.Type) %>%
  filter(b_Prime.Type > 0) %>%
  count()/4000
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.
##         n
## 1 0.94625

Chains look good.

#plot(m54mo2, N=4)
m54mo2.model.check <- createDHARMa(
  simulatedResponse = t(posterior_predict(m54mo2)), 
  observedResponse = DOD_54mo2$RESPONSE.CODE, 
  fittedPredictedResponse = apply(t(posterior_epred(m54mo2)), 1, mean), 
  integerResponse=TRUE
)

plot(m54mo2.model.check) # standard plots

testDispersion(m54mo2.model.check)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.74686, p-value = 0.006
## alternative hypothesis: two.sided
tab_model(m54mo1, m54mo2, show.re.var=FALSE, show.r2=FALSE, show.icc=FALSE, transform=NULL, collapse.ci=TRUE, p.style="scientific", ci.hyphen = " : ", dv.labels=c("Model 1",  "Model 2"), string.stat="",
          pred.labels=c(
            "Intercept", 
            "Prime (DOD)", 
            "Overlap", 
            "Prime * Overlap"
          ))
  Model 1 Model 2
Predictors Log-Odds Log-Odds
Intercept -3.21
(-4.80 : -1.42)
-2.81
(-4.45 : -1.09)
Prime (DOD) 0.83
(-0.27 : 1.99)
0.83
(-0.23 : 1.94)
Overlap -0.57
(-1.91 : 0.73)
-0.55
(-1.80 : 0.73)
Prime * Overlap 0.70
(-1.35 : 3.04)
0.75
(-1.17 : 3.14)
N 87 Blind_ID 75 Blind_ID
6 Target.verb 6 Target.verb
Observations 1011 872

4 Bias Mismatch

4.1 Full Data

DOD_54mo3 <- filter(DOD_54mo, Overlap=="No Overlap")
m54mo3 <- brm(RESPONSE.CODE ~ 1 +Prime.Type*biasmatch + (1 + Prime.Type*biasmatch| Blind_ID) + (1 + Prime.Type*biasmatch| Target.verb), seed=12345, control=list(adapt_delta=.9),   logit_prior, family="bernoulli", iter=4000, cores=4, file="54mo/output/m54mo3.rds", data=DOD_54mo3)
summary(m54mo3, prob=.9)
##  Family: bernoulli 
##   Links: mu = logit 
## Formula: RESPONSE.CODE ~ 1 + Prime.Type * biasmatch + (1 + Prime.Type * biasmatch | Blind_ID) + (1 + Prime.Type * biasmatch | Target.verb) 
##    Data: DOD_54mo3 (Number of observations: 479) 
##   Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup draws = 8000
## 
## Group-Level Effects: 
## ~Blind_ID (Number of levels: 87) 
##                                                     Estimate Est.Error l-90% CI
## sd(Intercept)                                           1.32      0.49     0.53
## sd(Prime.Type)                                          1.31      0.86     0.12
## sd(biasmatchmismatch)                                   1.05      0.69     0.11
## sd(Prime.Type:biasmatchmismatch)                        1.21      0.90     0.10
## cor(Intercept,Prime.Type)                               0.00      0.42    -0.69
## cor(Intercept,biasmatchmismatch)                        0.13      0.40    -0.56
## cor(Prime.Type,biasmatchmismatch)                       0.05      0.43    -0.67
## cor(Intercept,Prime.Type:biasmatchmismatch)             0.12      0.44    -0.64
## cor(Prime.Type,Prime.Type:biasmatchmismatch)           -0.13      0.45    -0.80
## cor(biasmatchmismatch,Prime.Type:biasmatchmismatch)     0.05      0.45    -0.70
##                                                     u-90% CI Rhat Bulk_ESS
## sd(Intercept)                                           2.16 1.00     1693
## sd(Prime.Type)                                          2.88 1.00     1834
## sd(biasmatchmismatch)                                   2.31 1.00     1635
## sd(Prime.Type:biasmatchmismatch)                        2.95 1.00     2880
## cor(Intercept,Prime.Type)                               0.68 1.00     4830
## cor(Intercept,biasmatchmismatch)                        0.76 1.00     4750
## cor(Prime.Type,biasmatchmismatch)                       0.73 1.00     3122
## cor(Intercept,Prime.Type:biasmatchmismatch)             0.77 1.00     7231
## cor(Prime.Type,Prime.Type:biasmatchmismatch)            0.64 1.00     5336
## cor(biasmatchmismatch,Prime.Type:biasmatchmismatch)     0.75 1.00     6242
##                                                     Tail_ESS
## sd(Intercept)                                           1537
## sd(Prime.Type)                                          3691
## sd(biasmatchmismatch)                                   2889
## sd(Prime.Type:biasmatchmismatch)                        3572
## cor(Intercept,Prime.Type)                               4828
## cor(Intercept,biasmatchmismatch)                        5632
## cor(Prime.Type,biasmatchmismatch)                       4850
## cor(Intercept,Prime.Type:biasmatchmismatch)             5769
## cor(Prime.Type,Prime.Type:biasmatchmismatch)            5873
## cor(biasmatchmismatch,Prime.Type:biasmatchmismatch)     6400
## 
## ~Target.verb (Number of levels: 6) 
##                                                     Estimate Est.Error l-90% CI
## sd(Intercept)                                           1.99      0.75     1.00
## sd(Prime.Type)                                          1.50      0.94     0.18
## sd(biasmatchmismatch)                                   0.86      0.72     0.06
## sd(Prime.Type:biasmatchmismatch)                        1.18      0.92     0.08
## cor(Intercept,Prime.Type)                               0.11      0.41    -0.60
## cor(Intercept,biasmatchmismatch)                        0.04      0.45    -0.69
## cor(Prime.Type,biasmatchmismatch)                      -0.05      0.44    -0.76
## cor(Intercept,Prime.Type:biasmatchmismatch)            -0.02      0.44    -0.73
## cor(Prime.Type,Prime.Type:biasmatchmismatch)            0.02      0.44    -0.71
## cor(biasmatchmismatch,Prime.Type:biasmatchmismatch)    -0.06      0.44    -0.76
##                                                     u-90% CI Rhat Bulk_ESS
## sd(Intercept)                                           3.39 1.00     4709
## sd(Prime.Type)                                          3.24 1.00     3254
## sd(biasmatchmismatch)                                   2.24 1.00     4062
## sd(Prime.Type:biasmatchmismatch)                        3.01 1.00     5048
## cor(Intercept,Prime.Type)                               0.73 1.00     8203
## cor(Intercept,biasmatchmismatch)                        0.75 1.00     7804
## cor(Prime.Type,biasmatchmismatch)                       0.69 1.00     7655
## cor(Intercept,Prime.Type:biasmatchmismatch)             0.71 1.00    10772
## cor(Prime.Type,Prime.Type:biasmatchmismatch)            0.74 1.00     8030
## cor(biasmatchmismatch,Prime.Type:biasmatchmismatch)     0.69 1.00     6905
##                                                     Tail_ESS
## sd(Intercept)                                           5032
## sd(Prime.Type)                                          3340
## sd(biasmatchmismatch)                                   4017
## sd(Prime.Type:biasmatchmismatch)                        4229
## cor(Intercept,Prime.Type)                               6101
## cor(Intercept,biasmatchmismatch)                        4789
## cor(Prime.Type,biasmatchmismatch)                       6239
## cor(Intercept,Prime.Type:biasmatchmismatch)             5760
## cor(Prime.Type,Prime.Type:biasmatchmismatch)            6046
## cor(biasmatchmismatch,Prime.Type:biasmatchmismatch)     6942
## 
## Population-Level Effects: 
##                              Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS
## Intercept                       -3.03      0.98    -4.66    -1.48 1.00     2717
## Prime.Type                       0.87      1.11    -0.92     2.73 1.00     4092
## biasmatchmismatch               -0.99      0.92    -2.60     0.33 1.00     3753
## Prime.Type:biasmatchmismatch     0.55      1.51    -1.79     3.16 1.00     4918
##                              Tail_ESS
## Intercept                        4108
## Prime.Type                       4534
## biasmatchmismatch                3870
## Prime.Type:biasmatchmismatch     4869
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
m54mo3 <- add_criterion(m54mo3, "loo")

Chains look fine.

#plot(m54mo3)
hypothesis(m54mo3, "Prime.Type > 0")
## Hypothesis Tests for class b:
##         Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
## 1 (Prime.Type) > 0     0.87      1.11    -0.92     2.73       4.17      0.81
##   Star
## 1     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
hypothesis(m54mo3, "Prime.Type:biasmatchmismatch > 0") # Bias mismatch
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (Prime.Type:biasm... > 0     0.55      1.51    -1.79     3.16       1.76
##   Post.Prob Star
## 1      0.64     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
m54mo3.model.check <- createDHARMa(
  simulatedResponse = t(posterior_predict(m54mo3)), 
  observedResponse = DOD_54mo3$RESPONSE.CODE, 
  fittedPredictedResponse = apply(t(posterior_epred(m54mo3)), 1, mean), 
  integerResponse=TRUE
)

plot(m54mo3.model.check) # standard plots

testDispersion(m54mo3.model.check)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.61694, p-value = 0.002
## alternative hypothesis: two.sided

4.2 Quick check of aggregation bias.

When looking at the data aggregated across verb, it looks like there is the predicted interaction between prime structure and bias mismatch (for which these is no evidence in the model above). However, I think the apparent bias effect reflects aggregation bias (see prime surprisal effects by verb above). I’ll fit models (a) without random effects by verbs, (b) with just random intercepts by verbs, and (c) without random effects but without gave and threw.

m54mo_aggreg1<- brm(RESPONSE.CODE ~ 1 +Prime.Type*biasmatch + (1 + Prime.Type*biasmatch| Blind_ID), seed=12345,  logit_prior, family="bernoulli", iter=4000, cores=4, control=list(adapt_delta=.9),   file="54mo/output/m54mo_aggreg1.rds", data=DOD_54mo3)

summary(m54mo_aggreg1, prob=.9)
##  Family: bernoulli 
##   Links: mu = logit 
## Formula: RESPONSE.CODE ~ 1 + Prime.Type * biasmatch + (1 + Prime.Type * biasmatch | Blind_ID) 
##    Data: DOD_54mo3 (Number of observations: 479) 
##   Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup draws = 8000
## 
## Group-Level Effects: 
## ~Blind_ID (Number of levels: 87) 
##                                                     Estimate Est.Error l-90% CI
## sd(Intercept)                                           0.85      0.40     0.15
## sd(Prime.Type)                                          0.69      0.51     0.05
## sd(biasmatchmismatch)                                   0.78      0.53     0.07
## sd(Prime.Type:biasmatchmismatch)                        0.82      0.64     0.06
## cor(Intercept,Prime.Type)                               0.02      0.44    -0.71
## cor(Intercept,biasmatchmismatch)                        0.04      0.43    -0.65
## cor(Prime.Type,biasmatchmismatch)                      -0.01      0.44    -0.73
## cor(Intercept,Prime.Type:biasmatchmismatch)             0.02      0.44    -0.71
## cor(Prime.Type,Prime.Type:biasmatchmismatch)           -0.15      0.46    -0.82
## cor(biasmatchmismatch,Prime.Type:biasmatchmismatch)    -0.07      0.45    -0.76
##                                                     u-90% CI Rhat Bulk_ESS
## sd(Intercept)                                           1.50 1.00     1295
## sd(Prime.Type)                                          1.65 1.00     2049
## sd(biasmatchmismatch)                                   1.74 1.00     1154
## sd(Prime.Type:biasmatchmismatch)                        2.05 1.00     2317
## cor(Intercept,Prime.Type)                               0.73 1.00     4871
## cor(Intercept,biasmatchmismatch)                        0.74 1.00     2997
## cor(Prime.Type,biasmatchmismatch)                       0.71 1.00     2530
## cor(Intercept,Prime.Type:biasmatchmismatch)             0.74 1.00     6038
## cor(Prime.Type,Prime.Type:biasmatchmismatch)            0.66 1.00     4491
## cor(biasmatchmismatch,Prime.Type:biasmatchmismatch)     0.69 1.00     4919
##                                                     Tail_ESS
## sd(Intercept)                                           1546
## sd(Prime.Type)                                          3078
## sd(biasmatchmismatch)                                   2540
## sd(Prime.Type:biasmatchmismatch)                        3038
## cor(Intercept,Prime.Type)                               4365
## cor(Intercept,biasmatchmismatch)                        3711
## cor(Prime.Type,biasmatchmismatch)                       4227
## cor(Intercept,Prime.Type:biasmatchmismatch)             5143
## cor(Prime.Type,Prime.Type:biasmatchmismatch)            6185
## cor(biasmatchmismatch,Prime.Type:biasmatchmismatch)     5949
## 
## Population-Level Effects: 
##                              Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS
## Intercept                       -2.53      0.32    -3.10    -2.05 1.00     4183
## Prime.Type                      -0.35      0.52    -1.24     0.49 1.00     5423
## biasmatchmismatch               -0.26      0.50    -1.13     0.50 1.00     2561
## Prime.Type:biasmatchmismatch     2.43      0.84     1.11     3.84 1.00     3418
##                              Tail_ESS
## Intercept                        4181
## Prime.Type                       5229
## biasmatchmismatch                3301
## Prime.Type:biasmatchmismatch     3441
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
hypothesis(m54mo_aggreg1, "Prime.Type > 0")
## Hypothesis Tests for class b:
##         Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
## 1 (Prime.Type) > 0    -0.35      0.52    -1.24     0.49       0.32      0.24
##   Star
## 1     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
hypothesis(m54mo_aggreg1, "Prime.Type:biasmatchmismatch > 0")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (Prime.Type:biasm... > 0     2.43      0.84     1.11     3.84     887.89
##   Post.Prob Star
## 1         1    *
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
m54mo_aggreg2<- brm(RESPONSE.CODE ~ 1 +Prime.Type*biasmatch + (1 + Prime.Type*biasmatch| Blind_ID) + (1 | Target.verb), seed=12345,  logit_prior, family="bernoulli", iter=4000, control=list(adapt_delta=.9), cores=4, file="54mo/output/m54mo_aggreg2.rds", data=DOD_54mo3)

summary(m54mo_aggreg2, prob=.9)
##  Family: bernoulli 
##   Links: mu = logit 
## Formula: RESPONSE.CODE ~ 1 + Prime.Type * biasmatch + (1 + Prime.Type * biasmatch | Blind_ID) + (1 | Target.verb) 
##    Data: DOD_54mo3 (Number of observations: 479) 
##   Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup draws = 8000
## 
## Group-Level Effects: 
## ~Blind_ID (Number of levels: 87) 
##                                                     Estimate Est.Error l-90% CI
## sd(Intercept)                                           1.40      0.51     0.57
## sd(Prime.Type)                                          1.24      0.84     0.11
## sd(biasmatchmismatch)                                   0.92      0.65     0.09
## sd(Prime.Type:biasmatchmismatch)                        1.10      0.84     0.09
## cor(Intercept,Prime.Type)                              -0.06      0.42    -0.73
## cor(Intercept,biasmatchmismatch)                        0.08      0.42    -0.61
## cor(Prime.Type,biasmatchmismatch)                       0.01      0.44    -0.73
## cor(Intercept,Prime.Type:biasmatchmismatch)             0.09      0.44    -0.65
## cor(Prime.Type,Prime.Type:biasmatchmismatch)           -0.15      0.45    -0.82
## cor(biasmatchmismatch,Prime.Type:biasmatchmismatch)    -0.00      0.44    -0.71
##                                                     u-90% CI Rhat Bulk_ESS
## sd(Intercept)                                           2.26 1.00     1837
## sd(Prime.Type)                                          2.79 1.00     1889
## sd(biasmatchmismatch)                                   2.14 1.00     1355
## sd(Prime.Type:biasmatchmismatch)                        2.74 1.00     3047
## cor(Intercept,Prime.Type)                               0.67 1.00     5462
## cor(Intercept,biasmatchmismatch)                        0.74 1.00     5452
## cor(Prime.Type,biasmatchmismatch)                       0.73 1.00     3381
## cor(Intercept,Prime.Type:biasmatchmismatch)             0.77 1.00     7549
## cor(Prime.Type,Prime.Type:biasmatchmismatch)            0.65 1.00     6580
## cor(biasmatchmismatch,Prime.Type:biasmatchmismatch)     0.71 1.00     6826
##                                                     Tail_ESS
## sd(Intercept)                                           1829
## sd(Prime.Type)                                          3548
## sd(biasmatchmismatch)                                   2663
## sd(Prime.Type:biasmatchmismatch)                        4183
## cor(Intercept,Prime.Type)                               5380
## cor(Intercept,biasmatchmismatch)                        5823
## cor(Prime.Type,biasmatchmismatch)                       4802
## cor(Intercept,Prime.Type:biasmatchmismatch)             5844
## cor(Prime.Type,Prime.Type:biasmatchmismatch)            5836
## cor(biasmatchmismatch,Prime.Type:biasmatchmismatch)     6623
## 
## ~Target.verb (Number of levels: 6) 
##               Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     1.89      0.69     1.01     3.19 1.00     4358     6141
## 
## Population-Level Effects: 
##                              Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS
## Intercept                       -3.15      0.95    -4.71    -1.62 1.00     3159
## Prime.Type                       1.00      0.74    -0.15     2.23 1.00     4843
## biasmatchmismatch               -0.37      0.60    -1.42     0.53 1.00     4488
## Prime.Type:biasmatchmismatch     0.17      1.20    -1.75     2.13 1.00     4258
##                              Tail_ESS
## Intercept                        3904
## Prime.Type                       3699
## biasmatchmismatch                3769
## Prime.Type:biasmatchmismatch     4260
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
hypothesis(m54mo_aggreg2, "Prime.Type > 0")
## Hypothesis Tests for class b:
##         Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
## 1 (Prime.Type) > 0        1      0.74    -0.15     2.23      12.36      0.93
##   Star
## 1     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
hypothesis(m54mo_aggreg2, "Prime.Type:biasmatchmismatch > 0")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (Prime.Type:biasm... > 0     0.17       1.2    -1.75     2.13       1.24
##   Post.Prob Star
## 1      0.55     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
m54mo_aggreg3<- DOD_54mo3 %>%
  filter(Target.verb != "gave" & Target.verb != "threw") %>%
 brm(RESPONSE.CODE ~ 1 +Prime.Type*biasmatch + (1 + Prime.Type*biasmatch| Blind_ID), seed=12345,  logit_prior, family="bernoulli", iter=4000, cores=4, control=list(adapt_delta=.9), file="54mo/output/m54mo_aggreg3.rds", data=.)

summary(m54mo_aggreg3, prob=.9)
##  Family: bernoulli 
##   Links: mu = logit 
## Formula: RESPONSE.CODE ~ 1 + Prime.Type * biasmatch + (1 + Prime.Type * biasmatch | Blind_ID) 
##    Data: . (Number of observations: 307) 
##   Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup draws = 8000
## 
## Group-Level Effects: 
## ~Blind_ID (Number of levels: 87) 
##                                                     Estimate Est.Error l-90% CI
## sd(Intercept)                                           1.14      0.63     0.17
## sd(Prime.Type)                                          1.29      0.91     0.10
## sd(biasmatchmismatch)                                   1.23      0.87     0.10
## sd(Prime.Type:biasmatchmismatch)                        1.80      1.27     0.14
## cor(Intercept,Prime.Type)                               0.10      0.44    -0.65
## cor(Intercept,biasmatchmismatch)                        0.07      0.43    -0.65
## cor(Prime.Type,biasmatchmismatch)                       0.08      0.44    -0.68
## cor(Intercept,Prime.Type:biasmatchmismatch)             0.02      0.43    -0.70
## cor(Prime.Type,Prime.Type:biasmatchmismatch)           -0.07      0.45    -0.76
## cor(biasmatchmismatch,Prime.Type:biasmatchmismatch)    -0.01      0.45    -0.74
##                                                     u-90% CI Rhat Bulk_ESS
## sd(Intercept)                                           2.24 1.00     1095
## sd(Prime.Type)                                          2.98 1.00     1871
## sd(biasmatchmismatch)                                   2.86 1.00     1296
## sd(Prime.Type:biasmatchmismatch)                        4.17 1.00     1627
## cor(Intercept,Prime.Type)                               0.77 1.00     3083
## cor(Intercept,biasmatchmismatch)                        0.76 1.00     3445
## cor(Prime.Type,biasmatchmismatch)                       0.77 1.00     2867
## cor(Intercept,Prime.Type:biasmatchmismatch)             0.73 1.00     4115
## cor(Prime.Type,Prime.Type:biasmatchmismatch)            0.69 1.00     3689
## cor(biasmatchmismatch,Prime.Type:biasmatchmismatch)     0.71 1.00     4110
##                                                     Tail_ESS
## sd(Intercept)                                           1426
## sd(Prime.Type)                                          3127
## sd(biasmatchmismatch)                                   2377
## sd(Prime.Type:biasmatchmismatch)                        2363
## cor(Intercept,Prime.Type)                               4144
## cor(Intercept,biasmatchmismatch)                        4342
## cor(Prime.Type,biasmatchmismatch)                       4562
## cor(Intercept,Prime.Type:biasmatchmismatch)             5042
## cor(Prime.Type,Prime.Type:biasmatchmismatch)            5116
## cor(biasmatchmismatch,Prime.Type:biasmatchmismatch)     5741
## 
## Population-Level Effects: 
##                              Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS
## Intercept                       -3.08      0.60    -4.20    -2.26 1.00     2408
## Prime.Type                       0.34      0.82    -1.01     1.66 1.00     3383
## biasmatchmismatch               -0.89      0.99    -2.76     0.43 1.00     2062
## Prime.Type:biasmatchmismatch     0.03      1.37    -2.17     2.30 1.00     3116
##                              Tail_ESS
## Intercept                        3419
## Prime.Type                       3942
## biasmatchmismatch                2100
## Prime.Type:biasmatchmismatch     3534
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
hypothesis(m54mo_aggreg3, "Prime.Type > 0")
## Hypothesis Tests for class b:
##         Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
## 1 (Prime.Type) > 0     0.34      0.82    -1.01     1.66       2.15      0.68
##   Star
## 1     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
hypothesis(m54mo_aggreg3, "Prime.Type:biasmatchmismatch > 0")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (Prime.Type:biasm... > 0     0.03      1.37    -2.17      2.3       1.01
##   Post.Prob Star
## 1       0.5     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

4.3 Grammatical Only

DOD_54mo4 <-
  DOD_54mo2 %>%
  filter(Overlap=="No Overlap")
m54mo4 <- brm(RESPONSE.CODE ~ 1 +Prime.Type*biasmatch + (1 + Prime.Type*biasmatch| Blind_ID) + (1 + Prime.Type*biasmatch| Target.verb), seed=12345,  logit_prior, family="bernoulli", iter=4000,  control=list(adapt_delta=.9),cores=4, file="54mo/output/m54mo4.rds", data=DOD_54mo4)
summary(m54mo4, prob=.9)
##  Family: bernoulli 
##   Links: mu = logit 
## Formula: RESPONSE.CODE ~ 1 + Prime.Type * biasmatch + (1 + Prime.Type * biasmatch | Blind_ID) + (1 + Prime.Type * biasmatch | Target.verb) 
##    Data: DOD_54mo4 (Number of observations: 412) 
##   Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup draws = 8000
## 
## Group-Level Effects: 
## ~Blind_ID (Number of levels: 75) 
##                                                     Estimate Est.Error l-90% CI
## sd(Intercept)                                           1.13      0.52     0.25
## sd(Prime.Type)                                          1.36      0.87     0.15
## sd(biasmatchmismatch)                                   1.12      0.71     0.12
## sd(Prime.Type:biasmatchmismatch)                        1.24      0.92     0.10
## cor(Intercept,Prime.Type)                              -0.01      0.42    -0.70
## cor(Intercept,biasmatchmismatch)                        0.11      0.41    -0.58
## cor(Prime.Type,biasmatchmismatch)                       0.06      0.43    -0.67
## cor(Intercept,Prime.Type:biasmatchmismatch)             0.09      0.44    -0.66
## cor(Prime.Type,Prime.Type:biasmatchmismatch)           -0.14      0.45    -0.81
## cor(biasmatchmismatch,Prime.Type:biasmatchmismatch)     0.03      0.45    -0.71
##                                                     u-90% CI Rhat Bulk_ESS
## sd(Intercept)                                           2.01 1.00     1249
## sd(Prime.Type)                                          2.93 1.00     2075
## sd(biasmatchmismatch)                                   2.42 1.00     1559
## sd(Prime.Type:biasmatchmismatch)                        2.98 1.00     3267
## cor(Intercept,Prime.Type)                               0.69 1.00     5718
## cor(Intercept,biasmatchmismatch)                        0.76 1.00     4226
## cor(Prime.Type,biasmatchmismatch)                       0.75 1.00     3402
## cor(Intercept,Prime.Type:biasmatchmismatch)             0.77 1.00     8717
## cor(Prime.Type,Prime.Type:biasmatchmismatch)            0.64 1.00     6976
## cor(biasmatchmismatch,Prime.Type:biasmatchmismatch)     0.75 1.00     6522
##                                                     Tail_ESS
## sd(Intercept)                                           1275
## sd(Prime.Type)                                          4145
## sd(biasmatchmismatch)                                   3757
## sd(Prime.Type:biasmatchmismatch)                        4841
## cor(Intercept,Prime.Type)                               5451
## cor(Intercept,biasmatchmismatch)                        4763
## cor(Prime.Type,biasmatchmismatch)                       5016
## cor(Intercept,Prime.Type:biasmatchmismatch)             5686
## cor(Prime.Type,Prime.Type:biasmatchmismatch)            6468
## cor(biasmatchmismatch,Prime.Type:biasmatchmismatch)     5967
## 
## ~Target.verb (Number of levels: 6) 
##                                                     Estimate Est.Error l-90% CI
## sd(Intercept)                                           1.98      0.74     1.00
## sd(Prime.Type)                                          1.31      0.91     0.13
## sd(biasmatchmismatch)                                   0.84      0.71     0.06
## sd(Prime.Type:biasmatchmismatch)                        1.17      0.91     0.09
## cor(Intercept,Prime.Type)                               0.10      0.41    -0.60
## cor(Intercept,biasmatchmismatch)                        0.04      0.44    -0.69
## cor(Prime.Type,biasmatchmismatch)                      -0.05      0.45    -0.75
## cor(Intercept,Prime.Type:biasmatchmismatch)            -0.02      0.44    -0.73
## cor(Prime.Type,Prime.Type:biasmatchmismatch)            0.01      0.45    -0.72
## cor(biasmatchmismatch,Prime.Type:biasmatchmismatch)    -0.06      0.45    -0.77
##                                                     u-90% CI Rhat Bulk_ESS
## sd(Intercept)                                           3.37 1.00     5486
## sd(Prime.Type)                                          3.03 1.00     3755
## sd(biasmatchmismatch)                                   2.24 1.00     5548
## sd(Prime.Type:biasmatchmismatch)                        2.93 1.00     6732
## cor(Intercept,Prime.Type)                               0.75 1.00    11966
## cor(Intercept,biasmatchmismatch)                        0.75 1.00    15059
## cor(Prime.Type,biasmatchmismatch)                       0.70 1.00     8959
## cor(Intercept,Prime.Type:biasmatchmismatch)             0.71 1.00    14491
## cor(Prime.Type,Prime.Type:biasmatchmismatch)            0.74 1.00    10468
## cor(biasmatchmismatch,Prime.Type:biasmatchmismatch)     0.70 1.00     7465
##                                                     Tail_ESS
## sd(Intercept)                                           6007
## sd(Prime.Type)                                          4250
## sd(biasmatchmismatch)                                   4627
## sd(Prime.Type:biasmatchmismatch)                        5155
## cor(Intercept,Prime.Type)                               6006
## cor(Intercept,biasmatchmismatch)                        5696
## cor(Prime.Type,biasmatchmismatch)                       6420
## cor(Intercept,Prime.Type:biasmatchmismatch)             5512
## cor(Prime.Type,Prime.Type:biasmatchmismatch)            6561
## cor(biasmatchmismatch,Prime.Type:biasmatchmismatch)     6717
## 
## Population-Level Effects: 
##                              Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS
## Intercept                       -2.72      0.98    -4.33    -1.16 1.00     3574
## Prime.Type                       0.87      1.05    -0.81     2.61 1.00     6421
## biasmatchmismatch               -0.99      0.90    -2.56     0.36 1.00     4954
## Prime.Type:biasmatchmismatch     0.63      1.50    -1.65     3.16 1.00     6980
##                              Tail_ESS
## Intercept                        5102
## Prime.Type                       5307
## biasmatchmismatch                4619
## Prime.Type:biasmatchmismatch     4703
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
m4 <- add_criterion(m54mo4, "loo")
hypothesis(m54mo4, "Prime.Type > 0")
## Hypothesis Tests for class b:
##         Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
## 1 (Prime.Type) > 0     0.87      1.05    -0.81     2.61       4.38      0.81
##   Star
## 1     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
hypothesis(m54mo4, "Prime.Type:biasmatchmismatch > 0") # Bias mismatch
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (Prime.Type:biasm... > 0     0.63       1.5    -1.65     3.16       1.86
##   Post.Prob Star
## 1      0.65     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
m54mo4.model.check <- createDHARMa(
  simulatedResponse = t(posterior_predict(m54mo4)), 
  observedResponse = DOD_54mo4$RESPONSE.CODE, 
  fittedPredictedResponse = apply(t(posterior_epred(m54mo4)), 1, mean), 
  integerResponse=TRUE
)

plot(m54mo4.model.check) # standard plots

testDispersion(m54mo4.model.check)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.62811, p-value = 0.00175
## alternative hypothesis: two.sided
tab_model(m54mo1, m54mo2, m54mo3, m54mo4,  show.ci=.90, show.re.var=FALSE, show.r2=FALSE, show.icc=FALSE, transform=NULL, collapse.ci=TRUE, p.style="scientific", ci.hyphen = " : ", dv.labels=c("Model 1",  "Model 2", "Model 3", "Model 4"), string.stat="", file="Table54", bpe="mean",
           pred.labels=c(
            "Intercept", 
            "Prime (DOD)", 
            "Overlap", 
            "Prime * Overlap",
            "Bias Mismatch (POD)", 
            "Prime * Bias Mismatch "
          ))
  Model 1 Model 2 Model 3 Model 4
Predictors Log-Odds Log-Odds Log-Odds Log-Odds
Intercept -3.19
(-4.51 : -1.78)
-2.80
(-4.16 : -1.43)
-3.03
(-4.66 : -1.48)
-2.72
(-4.33 : -1.16)
Prime (DOD) 0.84
(-0.06 : 1.76)
0.84
(-0.01 : 1.72)
0.87
(-0.92 : 2.73)
0.87
(-0.81 : 2.61)
Overlap -0.56
(-1.62 : 0.50)
-0.53
(-1.54 : 0.52)
Prime * Overlap 0.72
(-0.97 : 2.54)
0.80
(-0.76 : 2.57)
Bias Mismatch (POD) -0.99
(-2.60 : 0.33)
-0.99
(-2.56 : 0.36)
Prime * Bias Mismatch 0.55
(-1.79 : 3.16)
0.63
(-1.65 : 3.16)
N 87 Blind_ID 75 Blind_ID 87 Blind_ID 75 Blind_ID
6 Target.verb 6 Target.verb 6 Target.verb 6 Target.verb
Observations 1011 872 479 412