Analysis of Dative Priming Data at 42 months
Seamus Donnelly
1 Opening Files
library(brms)
library(readxl)
library(writexl)
library(DHARMa)
library(rmarkdown)
library(tidybayes)
library(ggplot2)
library(tidyr)
library(dplyr)
library(GGally)
library(ggpubr)
library(tidybayes)
library(patchwork)
library(sjPlot)
library(lubridate)
library(zoo)
library(readODS)DOD_42mo <- read_ods("DOD_42mo_blind.ods")
N_total <- nrow(DOD_42mo) # sanity checkCreate mismatch variables. I’m also creating some character versions of many factors to make plotting a bit easier later.
DOD_42mo <- DOD_42mo %>%
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_42mo$bias <-"DOD" # set as DOD bias
DOD_42mo$bias[DOD_42mo$Prime.verb=="passed"] = "POD" # Change POD-biased prime verbs.
DOD_42mo$bias[DOD_42mo$Prime.verb=="sent"] = "POD"
DOD_42mo$bias[DOD_42mo$Prime.verb=="threw"] = "POD"
DOD_42mo$biasmatch <-"match" # set to match
DOD_42mo$biasmatch[DOD_42mo$bias=="POD" & DOD_42mo$Prime == "DOD"] = "mismatch" # Change to mismatch.
DOD_42mo$biasmatch[DOD_42mo$bias=="DOD" & DOD_42mo$Prime == "POD"] = "mismatch"
table(DOD_42mo$biasmatch)##
## match mismatch
## 541 539DOD_42mo$biasmismatchc <- ifelse(DOD_42mo$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_42mo <- DOD_42mo %>%
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_42mo, Trial_Problem==1)) == (nrow(filter(DOD_42mo, List==6))/12) + (nrow(filter(DOD_42mo, List==8))/12)## [1] TRUEN_miss <- DOD_42mo %>%
filter(is.na(RESPONSE.CODE)) %>%
nrow()
N_Prob = nrow(filter(DOD_42mo, Trial_Problem==1)) See how many missing values we have by participant, and make indicator variables for (a) participants missing more than 10 trials and (b) participants missing more than 6 trials.
DOD_42mo %>%
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: 12 × 2
## # Groups: N_Missing [12]
## N_Missing n
## <int> <int>
## 1 0 25
## 2 1 19
## 3 2 15
## 4 3 10
## 5 4 3
## 6 5 1
## 7 6 6
## 8 8 4
## 9 9 1
## 10 10 2
## 11 11 1
## 12 12 3DOD_42mo <- DOD_42mo %>%
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), # Participants missing 10
MightDrop = ifelse(N_Missing > 6, 1, 0)
) %>%
full_join(DOD_42mo, .)## Joining with `by = join_by(Blind_ID)`N_drop = DOD_42mo %>%
filter(Drop == 1 & !is.na(RESPONSE.CODE)) %>% # save number of additional missing trials.
count()And by item.
DOD_42mo %>%
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: 11 × 2
## # Groups: Target.verb [11]
## Target.verb N_Missing
## <chr> <int>
## 1 brought 41
## 2 gave 38
## 3 passed 33
## 4 push 2
## 5 pushed 1
## 6 putting 1
## 7 rolled 1
## 8 sent 31
## 9 showed 46
## 10 threw 45
## 11 vroomed 1All non-target verbs (putting, rolled etc) 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_42mo %>%
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.254 0.436
## 2 DOD Overlap 0.219 0.414
## 3 POD No Overlap 0.253 0.435
## 4 POD Overlap 0.168 0.374DOD_42mo %>%
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.252 0.435
## 2 DOD mismatch 0.219 0.415
## 3 POD match 0.199 0.400
## 4 POD mismatch 0.219 0.414The differences in missingness between condition don’t look drastic. To sum up things, we have a fair amount of missing data, including 6 participants who are missing 10 or more trials. There are another 11 missing more than six trials. I’m going to exclude the 6 participants missing 10 or more, and I’m going to include the others, but test whether their inclusion affects the results in any meaningful way. The patterns of missingness aren’t drastically different across the experimental conditions.
DOD_42mo <- filter(DOD_42mo, !is.na(RESPONSE.CODE) & Drop==0) # 835
nrow(DOD_42mo) == N_total - N_miss - N_drop## n
## [1,] TRUEDescriptive Statistics
DOD_42mo %>%
summarise(
min = min(Session.Age),
max = max(Session.Age)
)## min max
## 1 41 : 27 43 : 53 Plot Data
(p1 <- DOD_42mo %>%
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="Prop DOD") +
theme(legend.position = "none")
) +
theme(text = element_text(size=20)) ## `summarise()` has grouped output by 'Prime'. You can override using the
## `.groups` argument.
ggsave("plots/pres1.png", last_plot(), dpi=500)## Saving 7 x 5 in imageDOD_42mo %>%
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.
Looks good–the difference between DODs in the DOD vs POD condition is bigger when there is lexical overlap for all verbs.
And let’s look at the pattern of responses when we only consider participants who produced at least 1 DOD.
DOD_42mo <- DOD_42mo %>%
group_by(Blind_ID) %>%
mutate(
Prod_DOD = sum(RESPONSE.CODE)
) %>%
dplyr::select(Blind_ID, Prod_DOD) %>%
slice(1) %>%
ungroup() %>%
full_join(DOD_42mo, ., by="Blind_ID")
writexl::write_xlsx(DOD_42mo, "DOD_42mo_prepped.xlsx")DOD_42mo2 <- filter(DOD_42mo, Prod_DOD > 0)
(p1b <- DOD_42mo2 %>%
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.
Look at it by verb, in case anything weird is going on.
DOD_42mo2 %>%
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.
Maybe a bit different for brought but we’re working with less data here and the overlap condition appears to be at floor.
Plot proportion of DOD responses by each cell in two-way interaction between biasmatch and prime type.
(p1c <- DOD_42mo %>%
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(.9)) +
ylim(0, .6) +
theme_minimal() +
labs(y="Raw Mean")
)## `summarise()` has grouped output by 'Prime'. You can override using the
## `.groups` argument.
DOD_42mo %>%
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)) %>%
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(.9)) +
facet_wrap(~Target.verb) ## `summarise()` has grouped output by 'Prime', 'biasmatch'. You can override
## using the `.groups` argument.
Notice that threw only had DODs in prime bias match and PODs in prime bias mismatch. This is because the two verbs that primed threw were both DOD-biased verbs. The opposite should be true for gave; however, on some trials participants produced gave instead of the target verb provided by the experimenter. We can see the numbers of those below:
DOD_42mo %>%
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 14
## 2 brought DOD mismatch 19
## 3 brought POD match 15
## 4 brought POD mismatch 15
## 5 gave DOD match 1
## 6 gave DOD mismatch 35
## 7 gave POD match 34
## 8 gave POD mismatch 2
## 9 passed DOD match 9
## 10 passed DOD mismatch 21
## # ℹ 12 more rowsAnd we can see more information about these 3 specific trials.
DOD_42mo %>%
filter(Overlap == "No Overlap") %>%
filter(Target.verb=="gave" & Prime=="DOD" & biasmatch == "match")## Blind_ID Gender Session.Date Session.Age Sess.Age.(Days) List Item
## 1 aOusW0Av 1 2018-03-02 00:00:00 UTC 42 : 20 1300 4 2r
## Prime.verb Target.verb Prime.Type Verb.match RESPONSE.CODE Prime Overlap
## 1 brought gave 0.5 -0.5 1 DOD No Overlap
## bias biasmatch biasmismatchc Trial_Drop Trial_Problem N_Missing Drop
## 1 DOD match -0.5 0 0 5 0
## MightDrop Prod_DOD
## 1 0 2DOD_42mo %>%
filter(Overlap == "No Overlap") %>%
filter(Target.verb=="gave" & Prime=="POD" & biasmatch == "mismatch")## Blind_ID Gender Session.Date Session.Age Sess.Age.(Days) List Item
## 1 AUD17Yu2 0 2018-05-17 00:00:00 UTC 42 : 21 1299 2 2r
## 2 LkVZQCTn 1 2019-04-15 00:00:00 UTC 42 : 23 1300 5 5
## Prime.verb Target.verb Prime.Type Verb.match RESPONSE.CODE Prime Overlap
## 1 brought gave -0.5 -0.5 0 POD No Overlap
## 2 brought gave -0.5 -0.5 1 POD No Overlap
## bias biasmatch biasmismatchc Trial_Drop Trial_Problem N_Missing Drop
## 1 DOD mismatch 0.5 0 0 0 0
## 2 DOD mismatch 0.5 0 0 3 0
## MightDrop Prod_DOD
## 1 0 0
## 2 0 34 Lexical Boost
4.1 All Participants
BRMS’s default priors on random effects are a bit wide (the 95% CI on the prior extends out to 10). Given that these are going to be difficult to estimate in the current design (very few words), I’ve decided to use slightly more informative priors. These will pull random effects away from really implausible values. These still aren’t especially informative – the 95% CI on the prior should extend out to a bit past 4, which is still a large value for a random effect sd.
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 complicate the model. So I’ll fit two variance component models with the two random effects specifications and compare them to one another.
m42mo0<- brm(RESPONSE.CODE ~ 1 + (1 | Blind_ID) + (1 | Target.verb), family="bernoulli", control=list(adapt_delta=.9), seed = 12345, prior=logit_prior, iter=4000, cores=4, data=DOD_42mo, file="42mo/output/m42mo0.rds")m42mo0b <- 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_42mo, file="42mo/output/m42mo0b.rds")Add loo for loo_compare
Simpler model fits numerically better.
Here’s the model of the lexical boost.
m42mo1 <- 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, iter = 8000, control=list(adapt_delta=.9), data=DOD_42mo, file="42mo/output/m42mo1.rds")
m42mo1<- add_criterion(m42mo1, "loo")Trace plots look fine (commented out for space)
#plot(m42mo1, N = 4)summary(m42mo1, 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_42mo (Number of observations: 835)
## Draws: 4 chains, each with iter = 8000; warmup = 4000; thin = 1;
## total post-warmup draws = 16000
##
## Group-Level Effects:
## ~Blind_ID (Number of levels: 84)
## Estimate Est.Error l-90% CI
## sd(Intercept) 1.97 0.36 1.43
## sd(Prime.Type) 1.01 0.57 0.12
## sd(OverlapOverlap) 0.67 0.48 0.06
## sd(Prime.Type:OverlapOverlap) 1.24 0.79 0.13
## cor(Intercept,Prime.Type) 0.17 0.38 -0.50
## cor(Intercept,OverlapOverlap) 0.30 0.41 -0.49
## cor(Prime.Type,OverlapOverlap) -0.02 0.43 -0.71
## cor(Intercept,Prime.Type:OverlapOverlap) 0.30 0.39 -0.44
## cor(Prime.Type,Prime.Type:OverlapOverlap) 0.08 0.43 -0.63
## cor(OverlapOverlap,Prime.Type:OverlapOverlap) 0.08 0.44 -0.67
## u-90% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 2.59 1.00 5509 10366
## sd(Prime.Type) 1.99 1.00 3650 4778
## sd(OverlapOverlap) 1.56 1.00 3793 6675
## sd(Prime.Type:OverlapOverlap) 2.65 1.00 4183 6509
## cor(Intercept,Prime.Type) 0.74 1.00 11305 11088
## cor(Intercept,OverlapOverlap) 0.86 1.00 12058 11767
## cor(Prime.Type,OverlapOverlap) 0.69 1.00 9652 10993
## cor(Intercept,Prime.Type:OverlapOverlap) 0.83 1.00 13278 10883
## cor(Prime.Type,Prime.Type:OverlapOverlap) 0.76 1.00 9979 11589
## cor(OverlapOverlap,Prime.Type:OverlapOverlap) 0.76 1.00 9206 10994
##
## ~Target.verb (Number of levels: 6)
## Estimate Est.Error l-90% CI
## sd(Intercept) 1.92 0.67 1.05
## sd(Prime.Type) 0.69 0.54 0.06
## sd(OverlapOverlap) 1.40 0.79 0.30
## sd(Prime.Type:OverlapOverlap) 1.08 0.77 0.10
## cor(Intercept,Prime.Type) -0.22 0.43 -0.84
## cor(Intercept,OverlapOverlap) 0.10 0.37 -0.53
## cor(Prime.Type,OverlapOverlap) -0.02 0.44 -0.73
## cor(Intercept,Prime.Type:OverlapOverlap) -0.11 0.42 -0.76
## cor(Prime.Type,Prime.Type:OverlapOverlap) -0.01 0.45 -0.74
## cor(OverlapOverlap,Prime.Type:OverlapOverlap) 0.08 0.44 -0.65
## u-90% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 3.16 1.00 8794 10773
## sd(Prime.Type) 1.75 1.00 8526 8453
## sd(OverlapOverlap) 2.87 1.00 6528 5039
## sd(Prime.Type:OverlapOverlap) 2.55 1.00 8017 6655
## cor(Intercept,Prime.Type) 0.55 1.00 19703 12768
## cor(Intercept,OverlapOverlap) 0.70 1.00 14563 12081
## cor(Prime.Type,OverlapOverlap) 0.70 1.00 8616 11531
## cor(Intercept,Prime.Type:OverlapOverlap) 0.61 1.00 19914 11715
## cor(Prime.Type,Prime.Type:OverlapOverlap) 0.73 1.00 13474 12413
## cor(OverlapOverlap,Prime.Type:OverlapOverlap) 0.77 1.00 14492 13115
##
## Population-Level Effects:
## Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS
## Intercept -2.87 0.91 -4.34 -1.38 1.00 5493
## Prime.Type 0.40 0.64 -0.60 1.46 1.00 10383
## OverlapOverlap -0.99 0.86 -2.49 0.29 1.00 6499
## Prime.Type:OverlapOverlap 1.29 1.01 -0.33 2.96 1.00 10068
## Tail_ESS
## Intercept 8356
## Prime.Type 10281
## OverlapOverlap 8575
## Prime.Type:OverlapOverlap 10476
##
## 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(m42mo1, "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.4 0.64 -0.6 1.46 2.96 0.75
## 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(m42mo1, "Prime.Type:OverlapOverlap > 0")## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (Prime.Type:Overl... > 0 1.29 1.01 -0.33 2.96 10.01
## Post.Prob Star
## 1 0.91
## ---
## '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(m42mo1, "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.7 1.02 0.1 3.41 22.63
## 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.m42mo1.model.check <- createDHARMa(
simulatedResponse = t(posterior_predict(m42mo1)),
observedResponse = DOD_42mo$RESPONSE.CODE,
fittedPredictedResponse = apply(t(posterior_epred(m42mo1)), 1, mean),
integerResponse=TRUE
)
plot(m42mo1.model.check) # standard plots
testDispersion(m42mo1.model.check)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.68528, p-value = 0.000625
## alternative hypothesis: two.sidedMild underdispersion and mild trend in residuals vs fitted.
Drop problem trials and re-fit. Drop participants with intermediate missing data and re-fit. Drop trials with high Pareto’s k and re-fit. None of these make a huge difference, so I’ve excluded this output from this html file, but see R Markdown file for full details.
4.2 DOD-Producers Only
Now we’ll consider only the kids who produced at least one DOD. The following output is really similar to that from the previous model, so I won’t add as much annotation.
m42mo2 <- 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, iter = 8000, data=DOD_42mo2, file="42mo/output/m42mo2.rds")
m42mo2<- add_criterion(m42mo2, "loo")Chains look fine
#plot(m42mo2, N=4)hypothesis(m42mo2, "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.53 0.59 -0.37 1.51 5.02 0.83
## 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(m42mo2, "Prime.Type:OverlapOverlap > 0")## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (Prime.Type:Overl... > 0 1.71 0.92 0.27 3.26 36.83
## Post.Prob Star
## 1 0.97 *
## ---
## '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(m42mo2, "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 2.24 0.95 0.81 3.85 110.89
## Post.Prob Star
## 1 0.99 *
## ---
## '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.m42mo2.model.check <- createDHARMa(
simulatedResponse = t(posterior_predict(m42mo2)),
observedResponse = DOD_42mo2$RESPONSE.CODE,
fittedPredictedResponse = apply(t(posterior_epred(m42mo2)), 1, mean),
integerResponse=TRUE
)
plot(m42mo2.model.check) # standard plots
testDispersion(m42mo2.model.check)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.73824, p-value = 0.00025
## alternative hypothesis: two.sidedMild underdispersion and mild trend in residuals vs fitted.
tab_model(m42mo1, m42mo2, 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 |
-2.87 (-4.68 : -1.04) |
-1.31 (-2.86 : 0.36) |
| Prime (DOD) |
0.38 (-0.82 : 1.71) |
0.50 (-0.55 : 1.78) |
| Overlap |
-0.93 (-2.89 : 0.62) |
-0.83 (-2.58 : 0.69) |
| Prime * Overlap |
1.28 (-0.72 : 3.35) |
1.67 (-0.03 : 3.64) |
| N | 84 Blind_ID | 46 Blind_ID |
| 6 Target.verb | 6 Target.verb | |
| Observations | 835 | 474 |
5 Bias Mismatch
Next consider how the biasmatch variable. For this, we’ll only consider the trials without lexical overlap.
DOD_42mo3 <- filter(DOD_42mo, Overlap=="No Overlap")
m42mo3 <- brm(RESPONSE.CODE ~ 1 +Prime.Type*biasmatch +
(1 + Prime.Type*biasmatch| Blind_ID) +
(1 + Prime.Type*biasmatch| Target.verb),
logit_prior,
family="bernoulli",
iter=8000,
seed = 12345,
cores=4,
control=list(adapt_delta=.90),
file="42mo/output/m42mo3.rds",
data=DOD_42mo3)
m42mo3 <- add_criterion(m42mo3, "loo")
summary(m42mo3, 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_42mo3 (Number of observations: 385)
## Draws: 4 chains, each with iter = 8000; warmup = 4000; thin = 1;
## total post-warmup draws = 16000
##
## Group-Level Effects:
## ~Blind_ID (Number of levels: 84)
## Estimate Est.Error l-90% CI
## sd(Intercept) 1.59 0.43 0.94
## sd(Prime.Type) 0.84 0.61 0.07
## sd(biasmatchmismatch) 0.73 0.54 0.06
## sd(Prime.Type:biasmatchmismatch) 0.86 0.67 0.07
## cor(Intercept,Prime.Type) 0.08 0.43 -0.64
## cor(Intercept,biasmatchmismatch) -0.00 0.43 -0.70
## cor(Prime.Type,biasmatchmismatch) -0.08 0.45 -0.78
## cor(Intercept,Prime.Type:biasmatchmismatch) 0.00 0.44 -0.71
## cor(Prime.Type,Prime.Type:biasmatchmismatch) -0.11 0.45 -0.80
## cor(biasmatchmismatch,Prime.Type:biasmatchmismatch) -0.04 0.45 -0.76
## u-90% CI Rhat Bulk_ESS
## sd(Intercept) 2.34 1.00 6204
## sd(Prime.Type) 1.98 1.00 7717
## sd(biasmatchmismatch) 1.75 1.00 5001
## sd(Prime.Type:biasmatchmismatch) 2.16 1.00 9748
## cor(Intercept,Prime.Type) 0.75 1.00 21820
## cor(Intercept,biasmatchmismatch) 0.71 1.00 20313
## cor(Prime.Type,biasmatchmismatch) 0.67 1.00 10017
## cor(Intercept,Prime.Type:biasmatchmismatch) 0.72 1.00 28098
## cor(Prime.Type,Prime.Type:biasmatchmismatch) 0.66 1.00 17476
## cor(biasmatchmismatch,Prime.Type:biasmatchmismatch) 0.70 1.00 14833
## Tail_ESS
## sd(Intercept) 8648
## sd(Prime.Type) 9801
## sd(biasmatchmismatch) 7893
## sd(Prime.Type:biasmatchmismatch) 9696
## cor(Intercept,Prime.Type) 12541
## cor(Intercept,biasmatchmismatch) 11062
## cor(Prime.Type,biasmatchmismatch) 12772
## cor(Intercept,Prime.Type:biasmatchmismatch) 11944
## cor(Prime.Type,Prime.Type:biasmatchmismatch) 14305
## cor(biasmatchmismatch,Prime.Type:biasmatchmismatch) 14713
##
## ~Target.verb (Number of levels: 6)
## Estimate Est.Error l-90% CI
## sd(Intercept) 1.90 0.69 1.00
## sd(Prime.Type) 0.87 0.69 0.07
## sd(biasmatchmismatch) 0.61 0.55 0.04
## sd(Prime.Type:biasmatchmismatch) 1.07 0.82 0.08
## cor(Intercept,Prime.Type) -0.12 0.43 -0.78
## cor(Intercept,biasmatchmismatch) -0.01 0.45 -0.74
## cor(Prime.Type,biasmatchmismatch) -0.00 0.45 -0.73
## cor(Intercept,Prime.Type:biasmatchmismatch) -0.16 0.44 -0.81
## cor(Prime.Type,Prime.Type:biasmatchmismatch) -0.02 0.45 -0.74
## cor(biasmatchmismatch,Prime.Type:biasmatchmismatch) -0.01 0.45 -0.74
## u-90% CI Rhat Bulk_ESS
## sd(Intercept) 3.22 1.00 12185
## sd(Prime.Type) 2.22 1.00 9903
## sd(biasmatchmismatch) 1.67 1.00 11659
## sd(Prime.Type:biasmatchmismatch) 2.67 1.00 12247
## cor(Intercept,Prime.Type) 0.64 1.00 26245
## cor(Intercept,biasmatchmismatch) 0.72 1.00 30724
## cor(Prime.Type,biasmatchmismatch) 0.75 1.00 17335
## cor(Intercept,Prime.Type:biasmatchmismatch) 0.61 1.00 28688
## cor(Prime.Type,Prime.Type:biasmatchmismatch) 0.71 1.00 20275
## cor(biasmatchmismatch,Prime.Type:biasmatchmismatch) 0.73 1.00 15062
## Tail_ESS
## sd(Intercept) 12432
## sd(Prime.Type) 11221
## sd(biasmatchmismatch) 10063
## sd(Prime.Type:biasmatchmismatch) 10860
## cor(Intercept,Prime.Type) 12343
## cor(Intercept,biasmatchmismatch) 12098
## cor(Prime.Type,biasmatchmismatch) 13476
## cor(Intercept,Prime.Type:biasmatchmismatch) 12286
## cor(Prime.Type,Prime.Type:biasmatchmismatch) 13617
## cor(biasmatchmismatch,Prime.Type:biasmatchmismatch) 13882
##
## Population-Level Effects:
## Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS
## Intercept -2.40 0.91 -3.89 -0.93 1.00 9445
## Prime.Type 0.92 0.85 -0.43 2.33 1.00 15549
## biasmatchmismatch -0.66 0.66 -1.75 0.38 1.00 15817
## Prime.Type:biasmatchmismatch -0.57 1.24 -2.55 1.51 1.00 15214
## Tail_ESS
## Intercept 10923
## Prime.Type 11824
## biasmatchmismatch 11919
## Prime.Type:biasmatchmismatch 11467
##
## 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).#plot(m42mo3, N=4)hypothesis(m42mo3, "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.92 0.85 -0.43 2.33 6.8 0.87
## 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(m42mo3, "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.57 1.24 -2.55 1.51 0.44
## Post.Prob Star
## 1 0.31
## ---
## '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.m42mo3.model.check <- createDHARMa(
simulatedResponse = t(posterior_predict(m42mo3)),
observedResponse = DOD_42mo3$RESPONSE.CODE,
fittedPredictedResponse = apply(t(posterior_epred(m42mo3)), 1, mean),
integerResponse=TRUE
)
plot(m42mo3.model.check) # standard plots
testDispersion(m42mo3.model.check)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.65452, p-value = 0.001375
## alternative hypothesis: two.sidedMild underdispersion and mild trend in residuals vs fitted.
p3 <- fitted(m42mo3, re_formula=NA, scale="response") %>%
cbind(DOD_42mo3) %>%
dplyr::select(Prime, biasmatch, Estimate, Est.Error) %>%
ggplot(aes(y=Estimate, x=biasmatch, fill=Prime )) +
geom_bar(stat="identity",position="dodge") +
geom_errorbar(aes(ymin=Estimate, ymax=Estimate+Est.Error), width=.2,
position=position_dodge(.9)) +
theme_minimal() +
ylim(0, .6) +
labs(y="Predicted Probability") +
theme(legend.position = "none")
p1c | p3 + plot_layout(guides = "collect") & theme(legend.position = "right") 
One thing stands out about that graph: we see the opposite pattern of
results in the raw data and in the model-implied values. Recall that
when we plotted the effect of prime surprisal by verb, we saw mostly the
opposite effect by verb as what we saw in the full sample. This means
the discrpenacy above could reflect Simpson’s Paradox: an effect on the
conditional distribution disappearing in the marginal distribution.
If the issue is Simpson’s Paradox, and not some sort of coding mistake, when we get rid of random effects by verb, the coefficients should change direction. I’ve fit the model without random effects by target word, and indeed the coefficient for the prime surprisal effect becomes positive. I’ve excluded it here for space reasons, but it’s in the R Markdown file. See longer explanation in 54 month script.
m42mo3e <- brm(RESPONSE.CODE ~ 1 +Prime.Type*biasmatch + (1 + Prime.Type*biasmatch| Blind_ID), seed = 12345, logit_prior, family="bernoulli", iter=4000, cores=4, file="42mo/output/m42mo3e.rds", control=list(adapt_delta=.9), data=DOD_42mo3)
summary(m42mo3e, prob=.9)5.1 DOD-Producers Only
DOD_42mo4 <-
DOD_42mo3 %>%
filter(Prod_DOD > 0)m42mo4 <- brm(RESPONSE.CODE ~ 1 +Prime.Type*biasmatch + (1 + Prime.Type*biasmatch| Blind_ID) + (1 + Prime.Type*biasmatch| Target.verb), logit_prior, family="bernoulli", seed = 12345, iter=8000, cores=4, file="42mo/output/m4.rds", data=DOD_42mo4)
summary(m42mo4, 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_42mo4 (Number of observations: 218)
## Draws: 4 chains, each with iter = 8000; warmup = 4000; thin = 1;
## total post-warmup draws = 16000
##
## Group-Level Effects:
## ~Blind_ID (Number of levels: 46)
## Estimate Est.Error l-90% CI
## sd(Intercept) 0.48 0.34 0.04
## sd(Prime.Type) 0.92 0.64 0.08
## sd(biasmatchmismatch) 0.68 0.50 0.06
## sd(Prime.Type:biasmatchmismatch) 0.89 0.69 0.07
## cor(Intercept,Prime.Type) 0.03 0.44 -0.69
## cor(Intercept,biasmatchmismatch) -0.06 0.44 -0.76
## cor(Prime.Type,biasmatchmismatch) -0.10 0.45 -0.79
## cor(Intercept,Prime.Type:biasmatchmismatch) -0.01 0.44 -0.73
## cor(Prime.Type,Prime.Type:biasmatchmismatch) -0.12 0.45 -0.81
## cor(biasmatchmismatch,Prime.Type:biasmatchmismatch) -0.04 0.45 -0.76
## u-90% CI Rhat Bulk_ESS
## sd(Intercept) 1.11 1.00 4555
## sd(Prime.Type) 2.11 1.00 4651
## sd(biasmatchmismatch) 1.61 1.00 4581
## sd(Prime.Type:biasmatchmismatch) 2.20 1.00 6634
## cor(Intercept,Prime.Type) 0.73 1.00 9340
## cor(Intercept,biasmatchmismatch) 0.68 1.00 10202
## cor(Prime.Type,biasmatchmismatch) 0.67 1.00 8120
## cor(Intercept,Prime.Type:biasmatchmismatch) 0.72 1.00 13959
## cor(Prime.Type,Prime.Type:biasmatchmismatch) 0.66 1.00 12418
## cor(biasmatchmismatch,Prime.Type:biasmatchmismatch) 0.71 1.00 12526
## Tail_ESS
## sd(Intercept) 5913
## sd(Prime.Type) 6143
## sd(biasmatchmismatch) 6970
## sd(Prime.Type:biasmatchmismatch) 7363
## cor(Intercept,Prime.Type) 10236
## cor(Intercept,biasmatchmismatch) 11021
## cor(Prime.Type,biasmatchmismatch) 10877
## cor(Intercept,Prime.Type:biasmatchmismatch) 11419
## cor(Prime.Type,Prime.Type:biasmatchmismatch) 11917
## cor(biasmatchmismatch,Prime.Type:biasmatchmismatch) 13725
##
## ~Target.verb (Number of levels: 6)
## Estimate Est.Error l-90% CI
## sd(Intercept) 1.92 0.67 1.02
## sd(Prime.Type) 0.78 0.64 0.06
## sd(biasmatchmismatch) 0.59 0.51 0.04
## sd(Prime.Type:biasmatchmismatch) 0.98 0.78 0.07
## cor(Intercept,Prime.Type) -0.09 0.44 -0.76
## cor(Intercept,biasmatchmismatch) -0.01 0.45 -0.73
## cor(Prime.Type,biasmatchmismatch) -0.02 0.45 -0.75
## cor(Intercept,Prime.Type:biasmatchmismatch) -0.15 0.44 -0.81
## cor(Prime.Type,Prime.Type:biasmatchmismatch) -0.04 0.45 -0.77
## cor(biasmatchmismatch,Prime.Type:biasmatchmismatch) -0.03 0.45 -0.76
## u-90% CI Rhat Bulk_ESS
## sd(Intercept) 3.17 1.00 7059
## sd(Prime.Type) 2.02 1.00 7312
## sd(biasmatchmismatch) 1.58 1.00 6675
## sd(Prime.Type:biasmatchmismatch) 2.51 1.00 9354
## cor(Intercept,Prime.Type) 0.66 1.00 17928
## cor(Intercept,biasmatchmismatch) 0.72 1.00 15754
## cor(Prime.Type,biasmatchmismatch) 0.72 1.00 12621
## cor(Intercept,Prime.Type:biasmatchmismatch) 0.63 1.00 17186
## cor(Prime.Type,Prime.Type:biasmatchmismatch) 0.71 1.00 14574
## cor(biasmatchmismatch,Prime.Type:biasmatchmismatch) 0.72 1.00 12167
## Tail_ESS
## sd(Intercept) 9822
## sd(Prime.Type) 7394
## sd(biasmatchmismatch) 7488
## sd(Prime.Type:biasmatchmismatch) 7011
## cor(Intercept,Prime.Type) 11664
## cor(Intercept,biasmatchmismatch) 11705
## cor(Prime.Type,biasmatchmismatch) 12742
## cor(Intercept,Prime.Type:biasmatchmismatch) 11271
## cor(Prime.Type,Prime.Type:biasmatchmismatch) 11510
## cor(biasmatchmismatch,Prime.Type:biasmatchmismatch) 12817
##
## Population-Level Effects:
## Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS
## Intercept -1.04 0.84 -2.41 0.33 1.00 4010
## Prime.Type 1.18 0.79 -0.05 2.51 1.00 9993
## biasmatchmismatch -0.66 0.59 -1.63 0.26 1.00 10641
## Prime.Type:biasmatchmismatch -0.85 1.16 -2.71 1.06 1.00 10108
## Tail_ESS
## Intercept 6876
## Prime.Type 9384
## biasmatchmismatch 9432
## Prime.Type:biasmatchmismatch 9533
##
## 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).m42mo4 <- add_criterion(m42mo4, "loo")#plot(m42mo4, N=4)hypothesis(m42mo4, "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.18 0.79 -0.05 2.51 16.78 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(m42mo4, "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.85 1.16 -2.71 1.06 0.29
## Post.Prob Star
## 1 0.22
## ---
## '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.m42mo4.model.check <- createDHARMa(
simulatedResponse = t(posterior_predict(m42mo4)),
observedResponse = DOD_42mo4$RESPONSE.CODE,
fittedPredictedResponse = apply(t(posterior_epred(m42mo4)), 1, mean),
integerResponse=TRUE
)
plot(m42mo4.model.check) # standard plots
testDispersion(m42mo4.model.check)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.79162, p-value = 0.02875
## alternative hypothesis: two.sidedMild underdispersion and mild trend in residuals vs fitted.
tab_model(m42mo1, m42mo2, m42mo3, m42mo4, 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="Table42", 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 |
-2.87 (-4.34 : -1.38) |
-1.30 (-2.59 : 0.02) |
-2.40 (-3.89 : -0.93) |
-1.04 (-2.41 : 0.33) |
| Prime (DOD) |
0.40 (-0.60 : 1.46) |
0.53 (-0.37 : 1.51) |
0.92 (-0.43 : 2.33) |
1.18 (-0.05 : 2.51) |
| Overlap |
-0.99 (-2.49 : 0.29) |
-0.86 (-2.22 : 0.37) |
||
| Prime * Overlap |
1.29 (-0.33 : 2.96) |
1.71 (0.27 : 3.26) |
||
| Bias Mismatch (POD) |
-0.66 (-1.75 : 0.38) |
-0.66 (-1.63 : 0.26) |
||
| Prime * Bias Mismatch |
-0.57 (-2.55 : 1.51) |
-0.85 (-2.71 : 1.06) |
||
| N | 84 Blind_ID | 46 Blind_ID | 84 Blind_ID | 46 Blind_ID |
| 6 Target.verb | 6 Target.verb | 6 Target.verb | 6 Target.verb | |
| Observations | 835 | 474 | 385 | 218 |