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)
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
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.
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 |
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
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.
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 |