Analysis of Cross Sectional Dative Priming
Seamus Donnelly
library(brms)
library(readxl)
library(rmarkdown)
library(tidybayes)
library(ggplot2)
library(ggthemes)
library(tidyr)
library(readODS)
library(dplyr)
library(GGally)
library(DHARMa)
library(ggpubr)
library(tidybayes)
library(patchwork)
library(sjPlot)
XSectional <- read_ods("DOD Priming Data Sheet_Wgender_May42020 (1).ods")1 Some preliminary stuff
1.1 Missing data
Drop invalid-trials plus invalid participants.
XSectional <- subset(XSectional, drop != 1)XSectional %>%
group_by(Case_ID) %>%
slice(1) %>%
ungroup() %>%
mutate(
age_group = round(Age, 0)
) %>%
group_by(age_group) %>%
count()1.2 Descriptive Statistics
Get number of trials for each cell type
XSectional <- XSectional %>%
mutate(
Prime = ifelse(Prime_type==.5, yes="DOD", no="POD"),
Overlap = ifelse(Verb_match==.5, yes="Overlap", no="No Overlap")
)
XSectional %>%
group_by(Prime, Trial_verb) %>%
summarise( n = n()
)## `summarise()` has grouped output by 'Prime'. You can override using the
## `.groups` argument.XSectional %>%
group_by(Overlap, Trial_verb) %>%
summarise(
n = n()
)## `summarise()` has grouped output by 'Overlap'. You can override using the
## `.groups` argument.XSectional %>%
group_by(Item) %>%
summarise(
n = n()
)XSectional %>%
group_by(List) %>%
summarise(
n = n()
)Create Prime bias variable. Note we will use this same characterization for the second study, with longitudinal data as I’m worried about controlling for a post-treatment variable.
XSectional$Agec = XSectional$Age - mean(XSectional$Age)
XSectional$bias <-"DOD"
XSectional %>%
group_by(Trial_verb) %>%
summarise(mean = mean(Response, na.rm=TRUE))XSectional$bias[XSectional$Prime_verb=="passed"] = "POD"
XSectional$bias[XSectional$Prime_verb=="sent"] = "POD"
XSectional$bias[XSectional$Prime_verb=="threw"] = "POD"XSectional$Bias_POD <- ifelse(XSectional$bias=="POD", yes = .5, no=-.5 )table(XSectional$bias, XSectional$Overlap, XSectional$Prime_type )## , , = -0.5
##
##
## No Overlap Overlap
## DOD 208 206
## POD 205 210
##
## , , = 0.5
##
##
## No Overlap Overlap
## DOD 202 209
## POD 212 206XSectional$biasmatch <-"match"
XSectional$biasmatch[XSectional$bias=="POD" & XSectional$Prime == "DOD"] = "mismatch"
XSectional$biasmatch[XSectional$bias=="DOD" & XSectional$Prime == "POD"] = "mismatch"
table(XSectional$biasmatch)##
## match mismatch
## 826 832XSectional$biasmismatchc <- ifelse(XSectional$biasmatch == "mismatch", yes=.5, no=-.5)
table(XSectional$biasmatch, XSectional$Overlap, XSectional$Prime_type) ## , , = -0.5
##
##
## No Overlap Overlap
## match 205 210
## mismatch 208 206
##
## , , = 0.5
##
##
## No Overlap Overlap
## match 202 209
## mismatch 212 2061.2.1 Plots of DoDs produced in each condition.
Now let’s look at the plots of the proportion of DoDs produced for the two two-way interactions. We’ll start with prime and overlap, i.e., the lexical boost.
XSectional %>%
group_by(Prime, Overlap) %>%
summarise(mean = mean(Response, na.rm=TRUE),
sd = sd(Response, na.rm=TRUE)) %>%
ggplot(aes(y=mean, x=Prime, fill=Overlap)) +
geom_bar(stat="identity",position="dodge") +
geom_errorbar(aes(ymin=mean, ymax=mean+sd), width=.2,
position=position_dodge(.9)) ## `summarise()` has grouped output by 'Prime'. You can override using the
## `.groups` argument.
And we’ll disaggregate by verb to make sure they’re all fairly similar.
XSectional %>%
group_by(Prime, Overlap, Trial_verb) %>%
summarise(mean = mean(Response, na.rm=TRUE),
sd = sd(Response, na.rm=TRUE)) %>%
ggplot(aes(y=mean, x=Prime, fill=Overlap)) +
geom_bar(stat="identity",position="dodge") +
geom_errorbar(aes(ymin=mean, ymax=mean+sd), width=.2,
position=position_dodge(.9)) +
facet_wrap(~Trial_verb)## `summarise()` has grouped output by 'Prime', 'Overlap'. You can override using
## the `.groups` argument.
Similar pattern for lexical boost in each condition: bigger effect of
priming with lexical overlap in each condition, though not visibly
obvious in all conditions. There certainly are not some verbs showing
the opposite pattern.
Ok now time to look at prime bias
XSectional %>%
filter(Overlap=="No Overlap") %>%
group_by(Prime, biasmatch) %>%
summarise(mean = mean(Response, na.rm=TRUE),
sd = sd(Response, na.rm=TRUE)) %>%
ggplot(aes(y=mean, x=Prime, fill=biasmatch)) +
geom_bar(stat="identity",position="dodge") +
geom_errorbar(aes(ymin=mean, ymax=mean+sd), width=.2,
position=position_dodge(.9)) ## `summarise()` has grouped output by 'Prime'. You can override using the
## `.groups` argument.
XSectional %>%
filter(Overlap=="No Overlap") %>%
group_by(Prime, biasmatch, Trial_verb) %>%
summarise(mean = mean(Response, na.rm=TRUE),
sd = sd(Response, na.rm=TRUE)) %>%
ggplot(aes(y=mean, x=Prime, fill=biasmatch)) +
geom_bar(stat="identity",position="dodge") +
geom_errorbar(aes(ymin=mean, ymax=mean+sd), width=.2,
position=position_dodge(.9)) +
facet_wrap(~Trial_verb)## `summarise()` has grouped output by 'Prime', 'biasmatch'. You can override
## using the `.groups` argument.
It’s a big more difficult to figure out what’s going on here. For a couple reasons.
First: the target verb threw was always primed with eitherbrought or showed, both of which DOD biased. As a result, bias mistmatch is confounded with prime structure within the target verb threw. Likewise gave was always primed by passed or sent both of which are POD biased. As such prime bias is confounded with prime structure with verb. Neither of these structures are confounded within the whole data set – but this might make estimating an interaction between prime bias and prime structure within trial type difficult. Of the four unconfounded verbs, we see the predicted effect for three priming effects are larger for mismatching conditions than matching conditions for brought, passed and sent but not showed. If anything showed illustrates the opposite effect, but this is a bit hard to see visually, as it seems to producing a lot of DODs in all 4 conditions.
In order to remain consistent with Peter et al (2015), we’ll only consider the trials without lexical overlap. But from the graph above, it’s worthing noting that the interaction between prime bias and prime is really different across trials with and without verb overlap. Looked at one way, it looks like priming effects are larger for bias mismatch, when there is no overlap. Looked at the other way, it seems like there is a lexical boost only when the verb bias matches.
1.3 Exploring the random effect structure.
Before getting to the substantive models, I fit several variance component models to determine an appropriate random effect structure. Note that because each participant responded to each verb twice, it is in principle possible to fit a model with random effects by participant, target verb and their interaction. I suspect this will greatly increase the estimation time, and will add a bunch of difficult to interpret parameters to the model. So I’m going to fit a a variance components model with and without this random factor, and compare the fit via loo.
logit_prior = c(prior(normal(0,2), class=sd))Model with random slopes by participant, by target verb and by by participant:trial verb.
mXsec0 <- brm(Response ~ 1 + (1 | Case_ID) + (1 | Trial_verb) + (1 | Case_ID:Trial_verb), logit_prior, control=list(adapt_delta = .9), seed=12345, family="bernoulli", iter=4000, cores=4, data=XSectional, file="mXsec0")
summary(mXsec0, prob = .9)## Family: bernoulli
## Links: mu = logit
## Formula: Response ~ 1 + (1 | Case_ID) + (1 | Trial_verb) + (1 | Case_ID:Trial_verb)
## Data: XSectional (Number of observations: 1658)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Group-Level Effects:
## ~Case_ID (Number of levels: 140)
## Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 1.51 0.15 1.28 1.76 1.00 2667 4823
##
## ~Case_ID:Trial_verb (Number of levels: 838)
## Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.20 0.15 0.01 0.49 1.01 1084 2211
##
## ~Trial_verb (Number of levels: 6)
## Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 1.07 0.43 0.57 1.91 1.00 2221 3493
##
## Population-Level Effects:
## Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS
## Intercept -1.56 0.47 -2.31 -0.79 1.00 1465 1763
##
## 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).Fit a simpler model without the interaction between participant and trial.
mXsec0b <- brm(Response ~ 1 + (1 | Case_ID) + (1 | Trial_verb), logit_prior, family="bernoulli", iter=4000, cores=4, seed=12345, data=XSectional, file="mXsec0b")
summary(mXsec0b, prob=.9)## Family: bernoulli
## Links: mu = logit
## Formula: Response ~ 1 + (1 | Case_ID) + (1 | Trial_verb)
## Data: XSectional (Number of observations: 1658)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Group-Level Effects:
## ~Case_ID (Number of levels: 140)
## Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 1.50 0.14 1.28 1.75 1.00 3225 5585
##
## ~Trial_verb (Number of levels: 6)
## Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 1.04 0.43 0.56 1.85 1.00 2676 3317
##
## Population-Level Effects:
## Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS
## Intercept -1.55 0.47 -2.28 -0.79 1.00 1955 2954
##
## 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).#mXsec0<- add_criterion(mXsec0, "loo")
#mXsec0b<- add_criterion(mXsec0b, "loo")#loo(mXsec0, mXsec0b)
#model_weights(mXsec0, mXsec0b)M0b is favored over m0 and this difference was greater than 2SE, so I think that given the added complexity of the part X item random effects isn’t adding anything, then it’s not worth the extra interpretational complexity when the other covariates are added.
I’m going to make another data set for participants who produced at least 1 DoD, since a participant needs to have acquired a construction to be primed by it. This was also done in the Kumarage et al (2022) paper
df_raw_part <- XSectional %>%
group_by(Case_ID) %>%
summarise(
Observed_Mean = mean(Response, na.rm=TRUE)
)
XSectional2 <- merge(XSectional, df_raw_part, by="Case_ID")
XSectional2 %>%
mutate(
Age_dis = round(Age, 0),
Age_dis = ifelse(Age_dis == 3, yes=4, no = ifelse( # add the < 3.5 to 4 (only 1 kid)
Age_dis ==9, yes=8, no= Age_dis # add the > 8.5 to 9 (only 1 kid)
))
) %>%
group_by(Case_ID) %>%
slice(1) %>%
ungroup() %>%
group_by(Age_dis) %>%
summarise(
No_DOD = sum(Observed_Mean ==0),
Total = n(),
Prop = No_DOD/Total
)XSectional2 <- subset(XSectional2, Observed_Mean != 0)2 Testing effects of lexical overlap
I’ve coded prime using effects coding and overlap using dummy coding. With this parameterization, the effect of prime tells me the priming effect in the non-overlap condition (i.e., the abstract priming effect), rather than the main effect. The interaction then tells us how this effect is moderated by lexical overlap.
2.1 Model on Full Sample.
mXsec1 <- brm(Response ~ 1 + Prime_type*Overlap*Agec + (1 + Prime_type*Overlap| Case_ID) + (1 + Prime_type*Overlap*Agec | Trial_verb), control=list(adapt_delta=.9), logit_prior, family="bernoulli", iter=4000, cores=4, seed=12345, data=XSectional, file="mXsec1")
summary(mXsec1, prob = .9)## Family: bernoulli
## Links: mu = logit
## Formula: Response ~ 1 + Prime_type * Overlap * Agec + (1 + Prime_type * Overlap | Case_ID) + (1 + Prime_type * Overlap * Agec | Trial_verb)
## Data: XSectional (Number of observations: 1658)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Group-Level Effects:
## ~Case_ID (Number of levels: 140)
## Estimate Est.Error l-90% CI
## sd(Intercept) 1.73 0.20 1.42
## sd(Prime_type) 0.53 0.34 0.05
## sd(OverlapOverlap) 0.40 0.27 0.04
## sd(Prime_type:OverlapOverlap) 0.94 0.56 0.10
## cor(Intercept,Prime_type) -0.24 0.38 -0.79
## cor(Intercept,OverlapOverlap) 0.30 0.40 -0.46
## cor(Prime_type,OverlapOverlap) -0.12 0.44 -0.78
## cor(Intercept,Prime_type:OverlapOverlap) -0.13 0.36 -0.69
## cor(Prime_type,Prime_type:OverlapOverlap) 0.03 0.43 -0.67
## cor(OverlapOverlap,Prime_type:OverlapOverlap) -0.10 0.43 -0.78
## u-90% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 2.07 1.00 3284 5448
## sd(Prime_type) 1.14 1.00 1999 4070
## sd(OverlapOverlap) 0.91 1.00 1809 3160
## sd(Prime_type:OverlapOverlap) 1.88 1.00 1355 2888
## cor(Intercept,Prime_type) 0.47 1.00 6763 4985
## cor(Intercept,OverlapOverlap) 0.84 1.00 6762 5825
## cor(Prime_type,OverlapOverlap) 0.64 1.00 3666 4897
## cor(Intercept,Prime_type:OverlapOverlap) 0.51 1.00 6790 4075
## cor(Prime_type,Prime_type:OverlapOverlap) 0.74 1.00 2752 4264
## cor(OverlapOverlap,Prime_type:OverlapOverlap) 0.64 1.00 2797 5287
##
## ~Trial_verb (Number of levels: 6)
## Estimate
## sd(Intercept) 1.40
## sd(Prime_type) 0.62
## sd(OverlapOverlap) 0.46
## sd(Agec) 0.37
## sd(Prime_type:OverlapOverlap) 0.79
## sd(Prime_type:Agec) 0.22
## sd(OverlapOverlap:Agec) 0.67
## sd(Prime_type:OverlapOverlap:Agec) 0.41
## cor(Intercept,Prime_type) -0.19
## cor(Intercept,OverlapOverlap) -0.18
## cor(Prime_type,OverlapOverlap) 0.01
## cor(Intercept,Agec) 0.16
## cor(Prime_type,Agec) -0.02
## cor(OverlapOverlap,Agec) -0.14
## cor(Intercept,Prime_type:OverlapOverlap) -0.13
## cor(Prime_type,Prime_type:OverlapOverlap) 0.02
## cor(OverlapOverlap,Prime_type:OverlapOverlap) -0.04
## cor(Agec,Prime_type:OverlapOverlap) -0.01
## cor(Intercept,Prime_type:Agec) -0.02
## cor(Prime_type,Prime_type:Agec) -0.02
## cor(OverlapOverlap,Prime_type:Agec) 0.02
## cor(Agec,Prime_type:Agec) 0.01
## cor(Prime_type:OverlapOverlap,Prime_type:Agec) -0.02
## cor(Intercept,OverlapOverlap:Agec) -0.27
## cor(Prime_type,OverlapOverlap:Agec) 0.00
## cor(OverlapOverlap,OverlapOverlap:Agec) 0.14
## cor(Agec,OverlapOverlap:Agec) -0.21
## cor(Prime_type:OverlapOverlap,OverlapOverlap:Agec) 0.01
## cor(Prime_type:Agec,OverlapOverlap:Agec) 0.02
## cor(Intercept,Prime_type:OverlapOverlap:Agec) -0.05
## cor(Prime_type,Prime_type:OverlapOverlap:Agec) 0.01
## cor(OverlapOverlap,Prime_type:OverlapOverlap:Agec) 0.02
## cor(Agec,Prime_type:OverlapOverlap:Agec) 0.03
## cor(Prime_type:OverlapOverlap,Prime_type:OverlapOverlap:Agec) 0.01
## cor(Prime_type:Agec,Prime_type:OverlapOverlap:Agec) -0.04
## cor(OverlapOverlap:Agec,Prime_type:OverlapOverlap:Agec) -0.02
## Est.Error
## sd(Intercept) 0.49
## sd(Prime_type) 0.44
## sd(OverlapOverlap) 0.33
## sd(Agec) 0.23
## sd(Prime_type:OverlapOverlap) 0.59
## sd(Prime_type:Agec) 0.22
## sd(OverlapOverlap:Agec) 0.34
## sd(Prime_type:OverlapOverlap:Agec) 0.39
## cor(Intercept,Prime_type) 0.31
## cor(Intercept,OverlapOverlap) 0.31
## cor(Prime_type,OverlapOverlap) 0.33
## cor(Intercept,Agec) 0.30
## cor(Prime_type,Agec) 0.32
## cor(OverlapOverlap,Agec) 0.32
## cor(Intercept,Prime_type:OverlapOverlap) 0.31
## cor(Prime_type,Prime_type:OverlapOverlap) 0.33
## cor(OverlapOverlap,Prime_type:OverlapOverlap) 0.33
## cor(Agec,Prime_type:OverlapOverlap) 0.32
## cor(Intercept,Prime_type:Agec) 0.33
## cor(Prime_type,Prime_type:Agec) 0.34
## cor(OverlapOverlap,Prime_type:Agec) 0.33
## cor(Agec,Prime_type:Agec) 0.33
## cor(Prime_type:OverlapOverlap,Prime_type:Agec) 0.34
## cor(Intercept,OverlapOverlap:Agec) 0.29
## cor(Prime_type,OverlapOverlap:Agec) 0.32
## cor(OverlapOverlap,OverlapOverlap:Agec) 0.32
## cor(Agec,OverlapOverlap:Agec) 0.31
## cor(Prime_type:OverlapOverlap,OverlapOverlap:Agec) 0.32
## cor(Prime_type:Agec,OverlapOverlap:Agec) 0.34
## cor(Intercept,Prime_type:OverlapOverlap:Agec) 0.32
## cor(Prime_type,Prime_type:OverlapOverlap:Agec) 0.33
## cor(OverlapOverlap,Prime_type:OverlapOverlap:Agec) 0.33
## cor(Agec,Prime_type:OverlapOverlap:Agec) 0.33
## cor(Prime_type:OverlapOverlap,Prime_type:OverlapOverlap:Agec) 0.34
## cor(Prime_type:Agec,Prime_type:OverlapOverlap:Agec) 0.34
## cor(OverlapOverlap:Agec,Prime_type:OverlapOverlap:Agec) 0.33
## l-90% CI u-90% CI
## sd(Intercept) 0.81 2.34
## sd(Prime_type) 0.06 1.43
## sd(OverlapOverlap) 0.05 1.08
## sd(Agec) 0.09 0.79
## sd(Prime_type:OverlapOverlap) 0.07 1.88
## sd(Prime_type:Agec) 0.01 0.62
## sd(OverlapOverlap:Agec) 0.26 1.30
## sd(Prime_type:OverlapOverlap:Agec) 0.03 1.15
## cor(Intercept,Prime_type) -0.67 0.36
## cor(Intercept,OverlapOverlap) -0.66 0.37
## cor(Prime_type,OverlapOverlap) -0.54 0.55
## cor(Intercept,Agec) -0.35 0.64
## cor(Prime_type,Agec) -0.55 0.52
## cor(OverlapOverlap,Agec) -0.64 0.42
## cor(Intercept,Prime_type:OverlapOverlap) -0.63 0.41
## cor(Prime_type,Prime_type:OverlapOverlap) -0.52 0.57
## cor(OverlapOverlap,Prime_type:OverlapOverlap) -0.57 0.52
## cor(Agec,Prime_type:OverlapOverlap) -0.55 0.52
## cor(Intercept,Prime_type:Agec) -0.56 0.53
## cor(Prime_type,Prime_type:Agec) -0.57 0.54
## cor(OverlapOverlap,Prime_type:Agec) -0.53 0.58
## cor(Agec,Prime_type:Agec) -0.54 0.56
## cor(Prime_type:OverlapOverlap,Prime_type:Agec) -0.58 0.54
## cor(Intercept,OverlapOverlap:Agec) -0.71 0.25
## cor(Prime_type,OverlapOverlap:Agec) -0.52 0.53
## cor(OverlapOverlap,OverlapOverlap:Agec) -0.41 0.65
## cor(Agec,OverlapOverlap:Agec) -0.68 0.33
## cor(Prime_type:OverlapOverlap,OverlapOverlap:Agec) -0.51 0.53
## cor(Prime_type:Agec,OverlapOverlap:Agec) -0.53 0.58
## cor(Intercept,Prime_type:OverlapOverlap:Agec) -0.56 0.49
## cor(Prime_type,Prime_type:OverlapOverlap:Agec) -0.55 0.56
## cor(OverlapOverlap,Prime_type:OverlapOverlap:Agec) -0.53 0.55
## cor(Agec,Prime_type:OverlapOverlap:Agec) -0.52 0.57
## cor(Prime_type:OverlapOverlap,Prime_type:OverlapOverlap:Agec) -0.55 0.56
## cor(Prime_type:Agec,Prime_type:OverlapOverlap:Agec) -0.59 0.51
## cor(OverlapOverlap:Agec,Prime_type:OverlapOverlap:Agec) -0.57 0.52
## Rhat Bulk_ESS
## sd(Intercept) 1.00 3894
## sd(Prime_type) 1.00 2510
## sd(OverlapOverlap) 1.00 3796
## sd(Agec) 1.00 3742
## sd(Prime_type:OverlapOverlap) 1.00 3593
## sd(Prime_type:Agec) 1.00 4562
## sd(OverlapOverlap:Agec) 1.00 3822
## sd(Prime_type:OverlapOverlap:Agec) 1.00 3736
## cor(Intercept,Prime_type) 1.00 9728
## cor(Intercept,OverlapOverlap) 1.00 10736
## cor(Prime_type,OverlapOverlap) 1.00 9860
## cor(Intercept,Agec) 1.00 8984
## cor(Prime_type,Agec) 1.00 7604
## cor(OverlapOverlap,Agec) 1.00 5790
## cor(Intercept,Prime_type:OverlapOverlap) 1.00 9771
## cor(Prime_type,Prime_type:OverlapOverlap) 1.00 9882
## cor(OverlapOverlap,Prime_type:OverlapOverlap) 1.00 8900
## cor(Agec,Prime_type:OverlapOverlap) 1.00 7436
## cor(Intercept,Prime_type:Agec) 1.00 13506
## cor(Prime_type,Prime_type:Agec) 1.00 10553
## cor(OverlapOverlap,Prime_type:Agec) 1.00 8267
## cor(Agec,Prime_type:Agec) 1.00 8237
## cor(Prime_type:OverlapOverlap,Prime_type:Agec) 1.00 7142
## cor(Intercept,OverlapOverlap:Agec) 1.00 8362
## cor(Prime_type,OverlapOverlap:Agec) 1.00 7648
## cor(OverlapOverlap,OverlapOverlap:Agec) 1.00 6059
## cor(Agec,OverlapOverlap:Agec) 1.00 7953
## cor(Prime_type:OverlapOverlap,OverlapOverlap:Agec) 1.00 6341
## cor(Prime_type:Agec,OverlapOverlap:Agec) 1.00 5886
## cor(Intercept,Prime_type:OverlapOverlap:Agec) 1.00 12984
## cor(Prime_type,Prime_type:OverlapOverlap:Agec) 1.00 10445
## cor(OverlapOverlap,Prime_type:OverlapOverlap:Agec) 1.00 8802
## cor(Agec,Prime_type:OverlapOverlap:Agec) 1.00 7542
## cor(Prime_type:OverlapOverlap,Prime_type:OverlapOverlap:Agec) 1.00 6752
## cor(Prime_type:Agec,Prime_type:OverlapOverlap:Agec) 1.00 5667
## cor(OverlapOverlap:Agec,Prime_type:OverlapOverlap:Agec) 1.00 6598
## Tail_ESS
## sd(Intercept) 5297
## sd(Prime_type) 2872
## sd(OverlapOverlap) 4244
## sd(Agec) 3810
## sd(Prime_type:OverlapOverlap) 4020
## sd(Prime_type:Agec) 4651
## sd(OverlapOverlap:Agec) 3944
## sd(Prime_type:OverlapOverlap:Agec) 4549
## cor(Intercept,Prime_type) 6296
## cor(Intercept,OverlapOverlap) 6178
## cor(Prime_type,OverlapOverlap) 6611
## cor(Intercept,Agec) 5868
## cor(Prime_type,Agec) 7034
## cor(OverlapOverlap,Agec) 5726
## cor(Intercept,Prime_type:OverlapOverlap) 5602
## cor(Prime_type,Prime_type:OverlapOverlap) 6609
## cor(OverlapOverlap,Prime_type:OverlapOverlap) 7011
## cor(Agec,Prime_type:OverlapOverlap) 6513
## cor(Intercept,Prime_type:Agec) 5636
## cor(Prime_type,Prime_type:Agec) 6057
## cor(OverlapOverlap,Prime_type:Agec) 6248
## cor(Agec,Prime_type:Agec) 6306
## cor(Prime_type:OverlapOverlap,Prime_type:Agec) 6801
## cor(Intercept,OverlapOverlap:Agec) 6464
## cor(Prime_type,OverlapOverlap:Agec) 7096
## cor(OverlapOverlap,OverlapOverlap:Agec) 5876
## cor(Agec,OverlapOverlap:Agec) 7178
## cor(Prime_type:OverlapOverlap,OverlapOverlap:Agec) 6704
## cor(Prime_type:Agec,OverlapOverlap:Agec) 6745
## cor(Intercept,Prime_type:OverlapOverlap:Agec) 6051
## cor(Prime_type,Prime_type:OverlapOverlap:Agec) 6587
## cor(OverlapOverlap,Prime_type:OverlapOverlap:Agec) 7118
## cor(Agec,Prime_type:OverlapOverlap:Agec) 7007
## cor(Prime_type:OverlapOverlap,Prime_type:OverlapOverlap:Agec) 6421
## cor(Prime_type:Agec,Prime_type:OverlapOverlap:Agec) 6497
## cor(OverlapOverlap:Agec,Prime_type:OverlapOverlap:Agec) 7230
##
## Population-Level Effects:
## Estimate Est.Error l-90% CI u-90% CI Rhat
## Intercept -1.61 0.62 -2.59 -0.59 1.00
## Prime_type 1.05 0.40 0.45 1.72 1.00
## OverlapOverlap -0.66 0.35 -1.22 -0.10 1.00
## Agec 0.45 0.23 0.09 0.82 1.00
## Prime_type:OverlapOverlap 1.64 0.61 0.73 2.69 1.00
## Prime_type:Agec -0.19 0.22 -0.54 0.16 1.00
## OverlapOverlap:Agec 0.06 0.34 -0.47 0.61 1.00
## Prime_type:OverlapOverlap:Agec 0.08 0.36 -0.48 0.64 1.00
## Bulk_ESS Tail_ESS
## Intercept 2140 3566
## Prime_type 4617 3959
## OverlapOverlap 3470 4006
## Agec 3555 4526
## Prime_type:OverlapOverlap 4451 4428
## Prime_type:Agec 7011 5345
## OverlapOverlap:Agec 4065 3875
## Prime_type:OverlapOverlap:Agec 6380 4405
##
## 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).mXsec1 <- add_criterion(mXsec1, "loo")2.1.1 Some diagnostics
Trace plots (commented out for space)
#plot(mXsec1, N = 4)Compare posterior predictive distribution to data.
pp_check(mXsec1)## Using 10 posterior draws for ppc type 'dens_overlay' by default.
Take a look at some regression diagnostics (QQ plots, residuals vs fitted, over/underdispersion)
mXsec1.model.check <- createDHARMa(
simulatedResponse = t(posterior_predict(mXsec1)),
observedResponse = XSectional$Response,
fittedPredictedResponse = apply(t(posterior_epred(mXsec1)), 1, mean),
integerResponse=TRUE
)
plot(mXsec1.model.check) # standard plots
testDispersion(mXsec1.model.check)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.78235, p-value < 2.2e-16
## alternative hypothesis: two.sidedplotResiduals(mXsec1.model.check, form=XSectional$Agec)
Mild, though statistically significant patterning in residuals, and some substantial underdispersion (possibly because of the random effects). In general I don’t see anything too problematic, though; we have a lot of data, so the significance tests are likely to be positive, and the trneds aren’t huge. Also, the author of the DHARMa package says that with mixed models and large samples, small artifactual trends in residuals can be significant.
Take a look at the distribution of Pareto’s k.
XSectional %>%
mutate(
k = mXsec1$criteria$loo$diagnostics$pareto_k
) %>%
ggplot(aes(x=k)) + geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
XSectional %>%
mutate(
k = mXsec1$criteria$loo$diagnostics$pareto_k
) %>%
filter(k >= .7)One value above .7, will re-fit without this observation.
2.1.2 Inferences
Caclulate the proportion of the posterior distribution in each direction.
hypothesis(mXsec1, "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.05 0.4 0.45 1.72 118.4 0.99
## 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(mXsec1, "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.64 0.61 0.73 2.69 152.85
## 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.hypothesis(mXsec1, "Agec > 0")## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
## 1 (Agec) > 0 0.45 0.23 0.09 0.82 39 0.98 *
## ---
## '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(mXsec1, "Prime_type:Agec < 0")## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (Prime_type:Agec) < 0 -0.19 0.22 -0.54 0.16 4.85
## Post.Prob Star
## 1 0.83
## ---
## '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(mXsec1, "Prime_type:OverlapOverlap:Agec > 0")## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (Prime_type:Overl... > 0 0.08 0.36 -0.48 0.64 1.41
## Post.Prob Star
## 1 0.59
## ---
## '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.Strong evidence for priming and lexical boost, less evidence for interaction between lexical boost and age, or priming effect and age.
cbPalette <- c("#000000", "#56B4E9", "#E69F00","#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
df <- tibble( # Design matrix (minus intercept which is automatic)
Overlap = factor(rep(c(0, 0, 1, 1), 7), labels = c("No Overlap", "Overlap")),
Prime_type = rep(c(-.5, .5, -.5, .5), 7),
Prime = rep(c("POD", "DOD", "POD", "DOD"), 7),
Agec = c(-3, -3, -3, -3, -2, -2, -2, -2, -1, -1, -1, -1, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3),
Age = c(3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9)
)
m1_post <- fitted(mXsec1, newdata=df, scale="linear", re_formula=NA, summary=FALSE, ndraws=25) %>%
t() %>%
as_tibble() %>%
bind_cols(df, .) %>%
pivot_longer(starts_with("V"), names_to="sample", values_to="linpred")## Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if
## `.name_repair` is omitted as of tibble 2.0.0.
## ℹ Using compatibility `.name_repair`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.(posterior_plot1 <- m1_post %>%
mutate(
con =paste(sample, Prime)
) %>%
ggplot(aes(y=linpred, x=Age, fill=Prime, color=Prime)) + geom_line(aes(fill=Prime, color=Prime, group=con)) + facet_wrap(~Overlap) + theme_minimal() + labs(y="Linear Predictor") + scale_colour_manual(values=cbPalette) + scale_fill_manual(values=cbPalette))## Warning in geom_line(aes(fill = Prime, color = Prime, group = con)): Ignoring
## unknown aesthetics: fill
points_model1 <- XSectional %>%
group_by(Part_No, Prime, Overlap) %>%
mutate(
mean = mean(Response)
) %>%
slice(1) %>%
ungroup() %>%
dplyr::select(Age, Prime, Overlap, mean)m1_fitted_response <- fitted(mXsec1, re_formula=NA, scale="response", probs=c(.05, .95)) %>%
cbind(XSectional) %>%
dplyr::select(Age, Prime, Overlap, Estimate, Q5, Q95)
m1_fitted_plot <- ggplot() +
geom_point(data= points_model1, aes(y=mean, x=Age, colour=Prime), shape=3) +
geom_line(data= m1_fitted_response , aes(y=Estimate, x=Age, fill=Prime, color=Prime, ymin=Q5, ymax=Q95, group=Prime), size=1.15) +
geom_ribbon(data= m1_fitted_response , aes(y=Estimate, ymin=Q5, ymax=Q95, x=Age, fill=Prime, group=Prime), alpha=0.5) +
facet_wrap(~Overlap)+
theme_minimal() +
labs(y="Probability") +
theme(
strip.text.x = element_blank(),
strip.background = element_blank()
) + scale_colour_manual (values=cbPalette) + scale_fill_manual(values=cbPalette)## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning in geom_line(data = m1_fitted_response, aes(y = Estimate, x = Age, :
## Ignoring unknown aesthetics: fill, ymin, and ymaxposterior_plot1/m1_fitted_plot + plot_layout(guides = "collect") & theme(legend.position = "right") 
ggsave("plots/Fig2.png", last_plot(), dpi=500)## Saving 7 x 5 in imageThe priming effect and lexical boost are clear. It looks like there is an interaction between age and priming but the evidence for this effect is not compelling (posterior probability = .83). This is more visible on the logit scale.
2.1.3 Simplify the random effects structure
Note that in the model above, I included correlated random effects. It took forever to run, and I essentially got the prior back for most of those correlations (modes of 0 with really large errors). I’m going to fit the model again with uncorrelated random effects, to make sure I get comparable results–in other words I want to make sure that all the uncertainty I’m adding to the mode isn’t affecting the posterior distributions for other parameters. UPDATE: We get roughly the same results. I’ve suppressed the output from these models in this HTML file, but the code is in the R Markdown file for those interested.
It’s also likely that I’ve overparameterized the model by including random effects for age and its interaction. I’ll try removing those.
Interesting to see it’s still not less underdispersed.
Let’s look at the results.
90% of the posterior distribution is in the predicted direction. Just to contextualize these results, I’ll see which model–with or without random effects by age–fits better. Given the high pareto’s k value in model 1, I’ll use reloo.
#loo1 <- loo(mXsec1)
#loo2 <- reloo(mXsec1, loo1) #
#loo3 <- loo(mXsec1bb)#loo_compare(loo2, loo3)Modest preference for model 1 (within a standard error. )
I also re-fit the model without observation with problematic pareto’s k. Check R markdown file for code/results
2.2 Only participants who produced at least 1 DOD.
Now I’ll re-fit the model on the participants who produced at least 1 DOD.
XSectional2 = filter(XSectional2, Observed_Mean > 0)
mXsec2 <- brm(Response ~ 1 + Prime_type*Overlap*Agec + (1 + Prime_type*Overlap| Case_ID) + (1 + Prime_type*Overlap*Agec | Trial_verb), logit_prior, family="bernoulli", iter=4000, cores=4, data=XSectional2, seed=12345, file="mXsec2")
summary(mXsec2, prob = .9)## Family: bernoulli
## Links: mu = logit
## Formula: Response ~ 1 + Prime_type * Overlap * Agec + (1 + Prime_type * Overlap | Case_ID) + (1 + Prime_type * Overlap * Agec | Trial_verb)
## Data: XSectional2 (Number of observations: 1348)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Group-Level Effects:
## ~Case_ID (Number of levels: 114)
## Estimate Est.Error l-90% CI
## sd(Intercept) 1.29 0.17 1.03
## sd(Prime_type) 0.61 0.36 0.08
## sd(OverlapOverlap) 0.44 0.28 0.04
## sd(Prime_type:OverlapOverlap) 0.96 0.57 0.10
## cor(Intercept,Prime_type) -0.37 0.35 -0.83
## cor(Intercept,OverlapOverlap) 0.31 0.37 -0.40
## cor(Prime_type,OverlapOverlap) -0.19 0.42 -0.80
## cor(Intercept,Prime_type:OverlapOverlap) -0.22 0.36 -0.74
## cor(Prime_type,Prime_type:OverlapOverlap) 0.07 0.43 -0.63
## cor(OverlapOverlap,Prime_type:OverlapOverlap) -0.14 0.43 -0.79
## u-90% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 1.58 1.00 3792 5572
## sd(Prime_type) 1.23 1.00 1939 3093
## sd(OverlapOverlap) 0.92 1.00 1762 3284
## sd(Prime_type:OverlapOverlap) 1.94 1.00 1368 2842
## cor(Intercept,Prime_type) 0.32 1.00 7628 5708
## cor(Intercept,OverlapOverlap) 0.84 1.00 7158 5376
## cor(Prime_type,OverlapOverlap) 0.58 1.00 4502 5778
## cor(Intercept,Prime_type:OverlapOverlap) 0.46 1.00 7557 4419
## cor(Prime_type,Prime_type:OverlapOverlap) 0.76 1.00 3744 5586
## cor(OverlapOverlap,Prime_type:OverlapOverlap) 0.63 1.00 3501 5203
##
## ~Trial_verb (Number of levels: 6)
## Estimate
## sd(Intercept) 1.38
## sd(Prime_type) 0.67
## sd(OverlapOverlap) 0.45
## sd(Agec) 0.36
## sd(Prime_type:OverlapOverlap) 0.76
## sd(Prime_type:Agec) 0.23
## sd(OverlapOverlap:Agec) 0.69
## sd(Prime_type:OverlapOverlap:Agec) 0.43
## cor(Intercept,Prime_type) -0.19
## cor(Intercept,OverlapOverlap) -0.18
## cor(Prime_type,OverlapOverlap) -0.01
## cor(Intercept,Agec) 0.17
## cor(Prime_type,Agec) -0.02
## cor(OverlapOverlap,Agec) -0.13
## cor(Intercept,Prime_type:OverlapOverlap) -0.12
## cor(Prime_type,Prime_type:OverlapOverlap) 0.02
## cor(OverlapOverlap,Prime_type:OverlapOverlap) -0.05
## cor(Agec,Prime_type:OverlapOverlap) -0.01
## cor(Intercept,Prime_type:Agec) -0.02
## cor(Prime_type,Prime_type:Agec) -0.02
## cor(OverlapOverlap,Prime_type:Agec) 0.02
## cor(Agec,Prime_type:Agec) 0.00
## cor(Prime_type:OverlapOverlap,Prime_type:Agec) -0.02
## cor(Intercept,OverlapOverlap:Agec) -0.27
## cor(Prime_type,OverlapOverlap:Agec) -0.01
## cor(OverlapOverlap,OverlapOverlap:Agec) 0.13
## cor(Agec,OverlapOverlap:Agec) -0.21
## cor(Prime_type:OverlapOverlap,OverlapOverlap:Agec) 0.01
## cor(Prime_type:Agec,OverlapOverlap:Agec) 0.05
## cor(Intercept,Prime_type:OverlapOverlap:Agec) -0.04
## cor(Prime_type,Prime_type:OverlapOverlap:Agec) -0.00
## cor(OverlapOverlap,Prime_type:OverlapOverlap:Agec) 0.02
## cor(Agec,Prime_type:OverlapOverlap:Agec) 0.03
## cor(Prime_type:OverlapOverlap,Prime_type:OverlapOverlap:Agec) 0.00
## cor(Prime_type:Agec,Prime_type:OverlapOverlap:Agec) -0.03
## cor(OverlapOverlap:Agec,Prime_type:OverlapOverlap:Agec) 0.00
## Est.Error
## sd(Intercept) 0.50
## sd(Prime_type) 0.45
## sd(OverlapOverlap) 0.33
## sd(Agec) 0.23
## sd(Prime_type:OverlapOverlap) 0.57
## sd(Prime_type:Agec) 0.22
## sd(OverlapOverlap:Agec) 0.34
## sd(Prime_type:OverlapOverlap:Agec) 0.39
## cor(Intercept,Prime_type) 0.31
## cor(Intercept,OverlapOverlap) 0.31
## cor(Prime_type,OverlapOverlap) 0.33
## cor(Intercept,Agec) 0.30
## cor(Prime_type,Agec) 0.32
## cor(OverlapOverlap,Agec) 0.33
## cor(Intercept,Prime_type:OverlapOverlap) 0.31
## cor(Prime_type,Prime_type:OverlapOverlap) 0.33
## cor(OverlapOverlap,Prime_type:OverlapOverlap) 0.33
## cor(Agec,Prime_type:OverlapOverlap) 0.33
## cor(Intercept,Prime_type:Agec) 0.33
## cor(Prime_type,Prime_type:Agec) 0.34
## cor(OverlapOverlap,Prime_type:Agec) 0.33
## cor(Agec,Prime_type:Agec) 0.33
## cor(Prime_type:OverlapOverlap,Prime_type:Agec) 0.33
## cor(Intercept,OverlapOverlap:Agec) 0.29
## cor(Prime_type,OverlapOverlap:Agec) 0.31
## cor(OverlapOverlap,OverlapOverlap:Agec) 0.32
## cor(Agec,OverlapOverlap:Agec) 0.31
## cor(Prime_type:OverlapOverlap,OverlapOverlap:Agec) 0.32
## cor(Prime_type:Agec,OverlapOverlap:Agec) 0.33
## cor(Intercept,Prime_type:OverlapOverlap:Agec) 0.32
## cor(Prime_type,Prime_type:OverlapOverlap:Agec) 0.34
## cor(OverlapOverlap,Prime_type:OverlapOverlap:Agec) 0.33
## cor(Agec,Prime_type:OverlapOverlap:Agec) 0.34
## cor(Prime_type:OverlapOverlap,Prime_type:OverlapOverlap:Agec) 0.34
## cor(Prime_type:Agec,Prime_type:OverlapOverlap:Agec) 0.34
## cor(OverlapOverlap:Agec,Prime_type:OverlapOverlap:Agec) 0.33
## l-90% CI u-90% CI
## sd(Intercept) 0.79 2.33
## sd(Prime_type) 0.09 1.51
## sd(OverlapOverlap) 0.05 1.06
## sd(Agec) 0.08 0.80
## sd(Prime_type:OverlapOverlap) 0.07 1.86
## sd(Prime_type:Agec) 0.02 0.62
## sd(OverlapOverlap:Agec) 0.28 1.34
## sd(Prime_type:OverlapOverlap:Agec) 0.03 1.20
## cor(Intercept,Prime_type) -0.66 0.35
## cor(Intercept,OverlapOverlap) -0.66 0.36
## cor(Prime_type,OverlapOverlap) -0.55 0.54
## cor(Intercept,Agec) -0.35 0.64
## cor(Prime_type,Agec) -0.54 0.51
## cor(OverlapOverlap,Agec) -0.64 0.44
## cor(Intercept,Prime_type:OverlapOverlap) -0.62 0.42
## cor(Prime_type,Prime_type:OverlapOverlap) -0.52 0.56
## cor(OverlapOverlap,Prime_type:OverlapOverlap) -0.58 0.50
## cor(Agec,Prime_type:OverlapOverlap) -0.55 0.54
## cor(Intercept,Prime_type:Agec) -0.56 0.51
## cor(Prime_type,Prime_type:Agec) -0.56 0.54
## cor(OverlapOverlap,Prime_type:Agec) -0.53 0.57
## cor(Agec,Prime_type:Agec) -0.54 0.56
## cor(Prime_type:OverlapOverlap,Prime_type:Agec) -0.57 0.53
## cor(Intercept,OverlapOverlap:Agec) -0.70 0.24
## cor(Prime_type,OverlapOverlap:Agec) -0.53 0.51
## cor(OverlapOverlap,OverlapOverlap:Agec) -0.42 0.63
## cor(Agec,OverlapOverlap:Agec) -0.68 0.34
## cor(Prime_type:OverlapOverlap,OverlapOverlap:Agec) -0.52 0.53
## cor(Prime_type:Agec,OverlapOverlap:Agec) -0.51 0.59
## cor(Intercept,Prime_type:OverlapOverlap:Agec) -0.57 0.50
## cor(Prime_type,Prime_type:OverlapOverlap:Agec) -0.55 0.55
## cor(OverlapOverlap,Prime_type:OverlapOverlap:Agec) -0.52 0.56
## cor(Agec,Prime_type:OverlapOverlap:Agec) -0.52 0.60
## cor(Prime_type:OverlapOverlap,Prime_type:OverlapOverlap:Agec) -0.55 0.55
## cor(Prime_type:Agec,Prime_type:OverlapOverlap:Agec) -0.58 0.52
## cor(OverlapOverlap:Agec,Prime_type:OverlapOverlap:Agec) -0.54 0.55
## Rhat Bulk_ESS
## sd(Intercept) 1.00 4103
## sd(Prime_type) 1.00 3506
## sd(OverlapOverlap) 1.00 4222
## sd(Agec) 1.00 3377
## sd(Prime_type:OverlapOverlap) 1.00 4354
## sd(Prime_type:Agec) 1.00 5257
## sd(OverlapOverlap:Agec) 1.00 4418
## sd(Prime_type:OverlapOverlap:Agec) 1.00 4378
## cor(Intercept,Prime_type) 1.00 9401
## cor(Intercept,OverlapOverlap) 1.00 11022
## cor(Prime_type,OverlapOverlap) 1.00 10975
## cor(Intercept,Agec) 1.00 11116
## cor(Prime_type,Agec) 1.00 9385
## cor(OverlapOverlap,Agec) 1.00 7622
## cor(Intercept,Prime_type:OverlapOverlap) 1.00 12403
## cor(Prime_type,Prime_type:OverlapOverlap) 1.00 11282
## cor(OverlapOverlap,Prime_type:OverlapOverlap) 1.00 8047
## cor(Agec,Prime_type:OverlapOverlap) 1.00 7806
## cor(Intercept,Prime_type:Agec) 1.00 15347
## cor(Prime_type,Prime_type:Agec) 1.00 12307
## cor(OverlapOverlap,Prime_type:Agec) 1.00 9037
## cor(Agec,Prime_type:Agec) 1.00 8420
## cor(Prime_type:OverlapOverlap,Prime_type:Agec) 1.00 7303
## cor(Intercept,OverlapOverlap:Agec) 1.00 10398
## cor(Prime_type,OverlapOverlap:Agec) 1.00 8440
## cor(OverlapOverlap,OverlapOverlap:Agec) 1.00 7017
## cor(Agec,OverlapOverlap:Agec) 1.00 7922
## cor(Prime_type:OverlapOverlap,OverlapOverlap:Agec) 1.00 6949
## cor(Prime_type:Agec,OverlapOverlap:Agec) 1.00 5911
## cor(Intercept,Prime_type:OverlapOverlap:Agec) 1.00 12569
## cor(Prime_type,Prime_type:OverlapOverlap:Agec) 1.00 11826
## cor(OverlapOverlap,Prime_type:OverlapOverlap:Agec) 1.00 8238
## cor(Agec,Prime_type:OverlapOverlap:Agec) 1.00 7946
## cor(Prime_type:OverlapOverlap,Prime_type:OverlapOverlap:Agec) 1.00 6912
## cor(Prime_type:Agec,Prime_type:OverlapOverlap:Agec) 1.00 5780
## cor(OverlapOverlap:Agec,Prime_type:OverlapOverlap:Agec) 1.00 5877
## Tail_ESS
## sd(Intercept) 5221
## sd(Prime_type) 3903
## sd(OverlapOverlap) 4247
## sd(Agec) 2728
## sd(Prime_type:OverlapOverlap) 4584
## sd(Prime_type:Agec) 5442
## sd(OverlapOverlap:Agec) 4259
## sd(Prime_type:OverlapOverlap:Agec) 4604
## cor(Intercept,Prime_type) 6198
## cor(Intercept,OverlapOverlap) 6371
## cor(Prime_type,OverlapOverlap) 6269
## cor(Intercept,Agec) 6385
## cor(Prime_type,Agec) 6523
## cor(OverlapOverlap,Agec) 6715
## cor(Intercept,Prime_type:OverlapOverlap) 5702
## cor(Prime_type,Prime_type:OverlapOverlap) 5519
## cor(OverlapOverlap,Prime_type:OverlapOverlap) 6197
## cor(Agec,Prime_type:OverlapOverlap) 7099
## cor(Intercept,Prime_type:Agec) 5719
## cor(Prime_type,Prime_type:Agec) 5973
## cor(OverlapOverlap,Prime_type:Agec) 6385
## cor(Agec,Prime_type:Agec) 6424
## cor(Prime_type:OverlapOverlap,Prime_type:Agec) 6004
## cor(Intercept,OverlapOverlap:Agec) 6721
## cor(Prime_type,OverlapOverlap:Agec) 7208
## cor(OverlapOverlap,OverlapOverlap:Agec) 6373
## cor(Agec,OverlapOverlap:Agec) 6811
## cor(Prime_type:OverlapOverlap,OverlapOverlap:Agec) 6921
## cor(Prime_type:Agec,OverlapOverlap:Agec) 6594
## cor(Intercept,Prime_type:OverlapOverlap:Agec) 6231
## cor(Prime_type,Prime_type:OverlapOverlap:Agec) 6001
## cor(OverlapOverlap,Prime_type:OverlapOverlap:Agec) 6558
## cor(Agec,Prime_type:OverlapOverlap:Agec) 6593
## cor(Prime_type:OverlapOverlap,Prime_type:OverlapOverlap:Agec) 6446
## cor(Prime_type:Agec,Prime_type:OverlapOverlap:Agec) 6388
## cor(OverlapOverlap:Agec,Prime_type:OverlapOverlap:Agec) 6991
##
## Population-Level Effects:
## Estimate Est.Error l-90% CI u-90% CI Rhat
## Intercept -1.06 0.60 -2.00 -0.07 1.00
## Prime_type 1.05 0.41 0.42 1.72 1.00
## OverlapOverlap -0.65 0.32 -1.17 -0.14 1.00
## Agec 0.35 0.22 0.01 0.68 1.00
## Prime_type:OverlapOverlap 1.68 0.60 0.79 2.69 1.00
## Prime_type:Agec -0.21 0.23 -0.57 0.15 1.00
## OverlapOverlap:Agec 0.05 0.34 -0.49 0.60 1.00
## Prime_type:OverlapOverlap:Agec 0.05 0.36 -0.53 0.62 1.00
## Bulk_ESS Tail_ESS
## Intercept 3041 3631
## Prime_type 5499 4709
## OverlapOverlap 4978 5311
## Agec 4897 4705
## Prime_type:OverlapOverlap 5755 5296
## Prime_type:Agec 8160 4768
## OverlapOverlap:Agec 5551 4996
## Prime_type:OverlapOverlap:Agec 7864 5514
##
## 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).mXsec2 <- add_criterion(mXsec2, "loo")#plot(mXsec2, N = 4)Check residuals
mXsec2.model.check <- createDHARMa(
simulatedResponse = t(posterior_predict(mXsec2)),
observedResponse = XSectional2$Response,
fittedPredictedResponse = apply(t(posterior_epred(mXsec2)), 1, mean),
integerResponse=TRUE
)
plot(mXsec2.model.check) # standard plots
testDispersion(mXsec2.model.check)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.79347, p-value < 2.2e-16
## alternative hypothesis: two.sidedplotResiduals(mXsec2.model.check, form=XSectional2$Agec)
Same results are before (underdispersion and small trends in the
residuals vs fitted), which does not vary by age.
hypothesis(mXsec2, "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.05 0.41 0.42 1.72 89.91 0.99
## 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(mXsec2, "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.68 0.6 0.79 2.69 241.42
## 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.hypothesis(mXsec2, "Agec > 0")## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
## 1 (Agec) > 0 0.35 0.22 0.01 0.68 20.39 0.95 *
## ---
## '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(mXsec2, "Prime_type:Agec < 0")## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (Prime_type:Agec) < 0 -0.21 0.23 -0.57 0.15 5.33
## Post.Prob Star
## 1 0.84
## ---
## '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(mXsec2, "Prime_type:OverlapOverlap:Agec > 0")## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (Prime_type:Overl... > 0 0.05 0.36 -0.53 0.62 1.22
## 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.We get pretty much the same pattern: clear lexical boost and clear priming, but no effect of age on the priming effect or the lexical boost.
I’m going to make similar plots as before. I’ve not included most of the code so as to maintain readability of this file
df <- tibble( # Design matrix (minus intercept which is automatic)
Overlap = factor(rep(c(0, 0, 1, 1), 7), labels = c("No Overlap", "Overlap")),
Prime_type = rep(c(-.5, .5, -.5, .5), 7),
Prime = rep(c("POD", "DOD", "POD", "DOD"), 7),
Agec = c(-3, -3, -3, -3, -2, -2, -2, -2, -1, -1, -1, -1, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3),
Age = c(3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9)
)
m2_post <- fitted(mXsec2, newdata=df, scale="linear", re_formula=NA, summary=FALSE, ndraws=25) %>%
t() %>%
as_tibble() %>%
bind_cols(df, .) %>%
pivot_longer(starts_with("V"), names_to="sample", values_to="linpred")
(posterior_plot2 <- m2_post %>%
mutate(
con =paste(sample, Prime)
) %>%
ggplot(aes(y=linpred, x=Age, fill=Prime, color=Prime)) + geom_line(aes(fill=Prime, color=Prime, group=con)) + facet_wrap(~Overlap) + theme_minimal() + labs(y="Linear Predictor") + scale_colour_manual (values=cbPalette) + scale_fill_manual(values=cbPalette))## Warning in geom_line(aes(fill = Prime, color = Prime, group = con)): Ignoring
## unknown aesthetics: fill
points_model2 <- XSectional %>%
group_by(Part_No, Prime, Overlap) %>%
mutate(
mean = mean(Response)
) %>%
slice(1) %>%
ungroup() %>%
dplyr::select(Age, Prime, Overlap, mean)m2_fitted_response <- fitted(mXsec2, re_formula=NA, scale="response", probs=c(.05, .95)) %>%
cbind(XSectional2) %>%
dplyr::select(Age, Prime, Overlap, Estimate, Q5, Q95)
m2_fitted_plot <- ggplot() +
geom_point(data= points_model2, aes(y=mean, x=Age, colour=Prime), shape=3) +
geom_line(data= m2_fitted_response , aes(y=Estimate, x=Age, fill=Prime, color=Prime, ymin=Q5, ymax=Q95, group=Prime), size=1.1) +
geom_ribbon(data= m2_fitted_response , aes(y=Estimate, ymin=Q5, ymax=Q95, x=Age, fill=Prime, group=Prime), alpha=0.5) +
facet_wrap(~Overlap)+
theme_minimal() +
labs(y="Probability") +
theme(
strip.text.x = element_blank(),
strip.background = element_blank()
) + scale_colour_manual (values=cbPalette) + scale_fill_manual(values=cbPalette)## Warning in geom_line(data = m2_fitted_response, aes(y = Estimate, x = Age, :
## Ignoring unknown aesthetics: fill, ymin, and ymaxposterior_plot2/m2_fitted_plot + plot_layout(guides = "collect") & theme(legend.position = "right") 
Check model with simpler random effects.
#loo1 <- loo(mXsec2)
#loo2 <- reloo(mXsec2, loo1) #
#loo3 <- loo(mXsec2bb)
#loo_compare(loo2, loo3)2.3 Table with models.
APA Style table with parameters of the model.
tab_model(mXsec1, mXsec2, 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",
"Age.c",
"Prime * Overlap",
"Prime * Age.c",
"Overlap * Age.c",
"Prime * Overlap * Age.c"
))| Model 1 | Model 2 | |
|---|---|---|
| Predictors | Log-Odds | Log-Odds |
| Intercept |
-1.62 (-2.82 : -0.32) |
-1.06 (-2.24 : 0.17) |
| Prime (DOD) |
1.04 (0.29 : 1.92) |
1.04 (0.25 : 1.89) |
| Overlap |
-0.66 (-1.34 : 0.02) |
-0.65 (-1.30 : 0.00) |
| Age.c |
0.45 (0.00 : 0.91) |
0.35 (-0.07 : 0.77) |
| Prime * Overlap |
1.61 (0.51 : 2.95) |
1.64 (0.60 : 2.96) |
| Prime * Age.c |
-0.19 (-0.62 : 0.24) |
-0.20 (-0.66 : 0.22) |
| Overlap * Age.c |
0.05 (-0.64 : 0.76) |
0.04 (-0.62 : 0.76) |
| Prime * Overlap * Age.c |
0.07 (-0.62 : 0.79) |
0.04 (-0.67 : 0.76) |
| N | 140 Case_ID | 114 Case_ID |
| 6 Trial_verb | 6 Trial_verb | |
| Observations | 1658 | 1348 |
3 Prime Bias.
Ok now I’ll look at the effect of primebias match. Peter et al (2015) looked at this in the verb match condition only. For consistency’s sake, I’ll do that here.
3.1 Full Sample
XSectional3 <- filter(XSectional, Overlap=="No Overlap")
table(XSectional3$biasmatch, XSectional3$Prime_type) ##
## -0.5 0.5
## match 205 202
## mismatch 208 212mXsec3 <- brm(Response ~ 1 +Prime_type*biasmatch*Agec+ (1 +Prime_type*biasmatch | Case_ID) + (1 +Prime_type*biasmatch*Agec | Trial_verb), logit_prior, family="bernoulli", iter=4000, cores=4, seed=12345,file="mXsec3", data=XSectional3)
summary(mXsec3 , prob = .9)## Family: bernoulli
## Links: mu = logit
## Formula: Response ~ 1 + Prime_type * biasmatch * Agec + (1 + Prime_type * biasmatch | Case_ID) + (1 + Prime_type * biasmatch * Agec | Trial_verb)
## Data: XSectional3 (Number of observations: 827)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Group-Level Effects:
## ~Case_ID (Number of levels: 140)
## Estimate Est.Error l-90% CI
## sd(Intercept) 1.76 0.27 1.34
## sd(Prime_type) 0.56 0.42 0.05
## sd(biasmatchmismatch) 0.60 0.41 0.05
## sd(Prime_type:biasmatchmismatch) 0.97 0.70 0.09
## cor(Intercept,Prime_type) -0.05 0.42 -0.72
## cor(Intercept,biasmatchmismatch) 0.12 0.41 -0.59
## cor(Prime_type,biasmatchmismatch) -0.03 0.44 -0.74
## cor(Intercept,Prime_type:biasmatchmismatch) 0.03 0.41 -0.65
## cor(Prime_type,Prime_type:biasmatchmismatch) -0.13 0.45 -0.81
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) 0.01 0.44 -0.72
## u-90% CI Rhat Bulk_ESS
## sd(Intercept) 2.21 1.00 2909
## sd(Prime_type) 1.34 1.00 2824
## sd(biasmatchmismatch) 1.33 1.00 1446
## sd(Prime_type:biasmatchmismatch) 2.29 1.00 2074
## cor(Intercept,Prime_type) 0.65 1.00 9105
## cor(Intercept,biasmatchmismatch) 0.76 1.00 7147
## cor(Prime_type,biasmatchmismatch) 0.70 1.00 2666
## cor(Intercept,Prime_type:biasmatchmismatch) 0.70 1.00 8527
## cor(Prime_type,Prime_type:biasmatchmismatch) 0.66 1.00 4156
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) 0.72 1.00 4616
## Tail_ESS
## sd(Intercept) 4332
## sd(Prime_type) 4036
## sd(biasmatchmismatch) 3142
## sd(Prime_type:biasmatchmismatch) 3634
## cor(Intercept,Prime_type) 6268
## cor(Intercept,biasmatchmismatch) 5528
## cor(Prime_type,biasmatchmismatch) 4504
## cor(Intercept,Prime_type:biasmatchmismatch) 5843
## cor(Prime_type,Prime_type:biasmatchmismatch) 4952
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) 5767
##
## ~Trial_verb (Number of levels: 6)
## Estimate
## sd(Intercept) 1.33
## sd(Prime_type) 0.57
## sd(biasmatchmismatch) 0.74
## sd(Agec) 0.40
## sd(Prime_type:biasmatchmismatch) 0.97
## sd(Prime_type:Agec) 0.43
## sd(biasmatchmismatch:Agec) 0.34
## sd(Prime_type:biasmatchmismatch:Agec) 0.71
## cor(Intercept,Prime_type) -0.11
## cor(Intercept,biasmatchmismatch) 0.09
## cor(Prime_type,biasmatchmismatch) -0.05
## cor(Intercept,Agec) 0.22
## cor(Prime_type,Agec) -0.05
## cor(biasmatchmismatch,Agec) 0.06
## cor(Intercept,Prime_type:biasmatchmismatch) -0.13
## cor(Prime_type,Prime_type:biasmatchmismatch) -0.02
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) -0.09
## cor(Agec,Prime_type:biasmatchmismatch) -0.05
## cor(Intercept,Prime_type:Agec) -0.03
## cor(Prime_type,Prime_type:Agec) 0.00
## cor(biasmatchmismatch,Prime_type:Agec) -0.02
## cor(Agec,Prime_type:Agec) -0.03
## cor(Prime_type:biasmatchmismatch,Prime_type:Agec) 0.03
## cor(Intercept,biasmatchmismatch:Agec) 0.01
## cor(Prime_type,biasmatchmismatch:Agec) 0.02
## cor(biasmatchmismatch,biasmatchmismatch:Agec) -0.01
## cor(Agec,biasmatchmismatch:Agec) -0.04
## cor(Prime_type:biasmatchmismatch,biasmatchmismatch:Agec) 0.01
## cor(Prime_type:Agec,biasmatchmismatch:Agec) -0.03
## cor(Intercept,Prime_type:biasmatchmismatch:Agec) 0.14
## cor(Prime_type,Prime_type:biasmatchmismatch:Agec) -0.04
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 0.08
## cor(Agec,Prime_type:biasmatchmismatch:Agec) 0.08
## cor(Prime_type:biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) -0.05
## cor(Prime_type:Agec,Prime_type:biasmatchmismatch:Agec) -0.09
## cor(biasmatchmismatch:Agec,Prime_type:biasmatchmismatch:Agec) -0.00
## Est.Error
## sd(Intercept) 0.54
## sd(Prime_type) 0.47
## sd(biasmatchmismatch) 0.55
## sd(Agec) 0.27
## sd(Prime_type:biasmatchmismatch) 0.72
## sd(Prime_type:Agec) 0.38
## sd(biasmatchmismatch:Agec) 0.32
## sd(Prime_type:biasmatchmismatch:Agec) 0.57
## cor(Intercept,Prime_type) 0.33
## cor(Intercept,biasmatchmismatch) 0.32
## cor(Prime_type,biasmatchmismatch) 0.33
## cor(Intercept,Agec) 0.31
## cor(Prime_type,Agec) 0.34
## cor(biasmatchmismatch,Agec) 0.33
## cor(Intercept,Prime_type:biasmatchmismatch) 0.32
## cor(Prime_type,Prime_type:biasmatchmismatch) 0.33
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) 0.33
## cor(Agec,Prime_type:biasmatchmismatch) 0.33
## cor(Intercept,Prime_type:Agec) 0.33
## cor(Prime_type,Prime_type:Agec) 0.34
## cor(biasmatchmismatch,Prime_type:Agec) 0.33
## cor(Agec,Prime_type:Agec) 0.34
## cor(Prime_type:biasmatchmismatch,Prime_type:Agec) 0.33
## cor(Intercept,biasmatchmismatch:Agec) 0.33
## cor(Prime_type,biasmatchmismatch:Agec) 0.34
## cor(biasmatchmismatch,biasmatchmismatch:Agec) 0.34
## cor(Agec,biasmatchmismatch:Agec) 0.34
## cor(Prime_type:biasmatchmismatch,biasmatchmismatch:Agec) 0.34
## cor(Prime_type:Agec,biasmatchmismatch:Agec) 0.34
## cor(Intercept,Prime_type:biasmatchmismatch:Agec) 0.32
## cor(Prime_type,Prime_type:biasmatchmismatch:Agec) 0.34
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 0.33
## cor(Agec,Prime_type:biasmatchmismatch:Agec) 0.33
## cor(Prime_type:biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 0.33
## cor(Prime_type:Agec,Prime_type:biasmatchmismatch:Agec) 0.34
## cor(biasmatchmismatch:Agec,Prime_type:biasmatchmismatch:Agec) 0.34
## l-90% CI
## sd(Intercept) 0.68
## sd(Prime_type) 0.05
## sd(biasmatchmismatch) 0.07
## sd(Agec) 0.06
## sd(Prime_type:biasmatchmismatch) 0.08
## sd(Prime_type:Agec) 0.03
## sd(biasmatchmismatch:Agec) 0.02
## sd(Prime_type:biasmatchmismatch:Agec) 0.05
## cor(Intercept,Prime_type) -0.64
## cor(Intercept,biasmatchmismatch) -0.45
## cor(Prime_type,biasmatchmismatch) -0.59
## cor(Intercept,Agec) -0.34
## cor(Prime_type,Agec) -0.59
## cor(biasmatchmismatch,Agec) -0.51
## cor(Intercept,Prime_type:biasmatchmismatch) -0.63
## cor(Prime_type,Prime_type:biasmatchmismatch) -0.57
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) -0.62
## cor(Agec,Prime_type:biasmatchmismatch) -0.59
## cor(Intercept,Prime_type:Agec) -0.57
## cor(Prime_type,Prime_type:Agec) -0.56
## cor(biasmatchmismatch,Prime_type:Agec) -0.57
## cor(Agec,Prime_type:Agec) -0.59
## cor(Prime_type:biasmatchmismatch,Prime_type:Agec) -0.53
## cor(Intercept,biasmatchmismatch:Agec) -0.55
## cor(Prime_type,biasmatchmismatch:Agec) -0.54
## cor(biasmatchmismatch,biasmatchmismatch:Agec) -0.56
## cor(Agec,biasmatchmismatch:Agec) -0.59
## cor(Prime_type:biasmatchmismatch,biasmatchmismatch:Agec) -0.55
## cor(Prime_type:Agec,biasmatchmismatch:Agec) -0.58
## cor(Intercept,Prime_type:biasmatchmismatch:Agec) -0.41
## cor(Prime_type,Prime_type:biasmatchmismatch:Agec) -0.58
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) -0.48
## cor(Agec,Prime_type:biasmatchmismatch:Agec) -0.48
## cor(Prime_type:biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) -0.60
## cor(Prime_type:Agec,Prime_type:biasmatchmismatch:Agec) -0.64
## cor(biasmatchmismatch:Agec,Prime_type:biasmatchmismatch:Agec) -0.56
## u-90% CI
## sd(Intercept) 2.33
## sd(Prime_type) 1.47
## sd(biasmatchmismatch) 1.77
## sd(Agec) 0.90
## sd(Prime_type:biasmatchmismatch) 2.34
## sd(Prime_type:Agec) 1.14
## sd(biasmatchmismatch:Agec) 0.96
## sd(Prime_type:biasmatchmismatch:Agec) 1.78
## cor(Intercept,Prime_type) 0.44
## cor(Intercept,biasmatchmismatch) 0.60
## cor(Prime_type,biasmatchmismatch) 0.51
## cor(Intercept,Agec) 0.70
## cor(Prime_type,Agec) 0.50
## cor(biasmatchmismatch,Agec) 0.59
## cor(Intercept,Prime_type:biasmatchmismatch) 0.42
## cor(Prime_type,Prime_type:biasmatchmismatch) 0.54
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) 0.47
## cor(Agec,Prime_type:biasmatchmismatch) 0.50
## cor(Intercept,Prime_type:Agec) 0.51
## cor(Prime_type,Prime_type:Agec) 0.56
## cor(biasmatchmismatch,Prime_type:Agec) 0.53
## cor(Agec,Prime_type:Agec) 0.53
## cor(Prime_type:biasmatchmismatch,Prime_type:Agec) 0.57
## cor(Intercept,biasmatchmismatch:Agec) 0.56
## cor(Prime_type,biasmatchmismatch:Agec) 0.56
## cor(biasmatchmismatch,biasmatchmismatch:Agec) 0.55
## cor(Agec,biasmatchmismatch:Agec) 0.52
## cor(Prime_type:biasmatchmismatch,biasmatchmismatch:Agec) 0.56
## cor(Prime_type:Agec,biasmatchmismatch:Agec) 0.54
## cor(Intercept,Prime_type:biasmatchmismatch:Agec) 0.65
## cor(Prime_type,Prime_type:biasmatchmismatch:Agec) 0.53
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 0.60
## cor(Agec,Prime_type:biasmatchmismatch:Agec) 0.62
## cor(Prime_type:biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 0.49
## cor(Prime_type:Agec,Prime_type:biasmatchmismatch:Agec) 0.48
## cor(biasmatchmismatch:Agec,Prime_type:biasmatchmismatch:Agec) 0.55
## Rhat
## sd(Intercept) 1.00
## sd(Prime_type) 1.00
## sd(biasmatchmismatch) 1.00
## sd(Agec) 1.00
## sd(Prime_type:biasmatchmismatch) 1.00
## sd(Prime_type:Agec) 1.00
## sd(biasmatchmismatch:Agec) 1.00
## sd(Prime_type:biasmatchmismatch:Agec) 1.00
## cor(Intercept,Prime_type) 1.00
## cor(Intercept,biasmatchmismatch) 1.00
## cor(Prime_type,biasmatchmismatch) 1.00
## cor(Intercept,Agec) 1.00
## cor(Prime_type,Agec) 1.00
## cor(biasmatchmismatch,Agec) 1.00
## cor(Intercept,Prime_type:biasmatchmismatch) 1.00
## cor(Prime_type,Prime_type:biasmatchmismatch) 1.00
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) 1.00
## cor(Agec,Prime_type:biasmatchmismatch) 1.00
## cor(Intercept,Prime_type:Agec) 1.00
## cor(Prime_type,Prime_type:Agec) 1.00
## cor(biasmatchmismatch,Prime_type:Agec) 1.00
## cor(Agec,Prime_type:Agec) 1.00
## cor(Prime_type:biasmatchmismatch,Prime_type:Agec) 1.00
## cor(Intercept,biasmatchmismatch:Agec) 1.00
## cor(Prime_type,biasmatchmismatch:Agec) 1.00
## cor(biasmatchmismatch,biasmatchmismatch:Agec) 1.00
## cor(Agec,biasmatchmismatch:Agec) 1.00
## cor(Prime_type:biasmatchmismatch,biasmatchmismatch:Agec) 1.00
## cor(Prime_type:Agec,biasmatchmismatch:Agec) 1.00
## cor(Intercept,Prime_type:biasmatchmismatch:Agec) 1.00
## cor(Prime_type,Prime_type:biasmatchmismatch:Agec) 1.00
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 1.00
## cor(Agec,Prime_type:biasmatchmismatch:Agec) 1.00
## cor(Prime_type:biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 1.00
## cor(Prime_type:Agec,Prime_type:biasmatchmismatch:Agec) 1.00
## cor(biasmatchmismatch:Agec,Prime_type:biasmatchmismatch:Agec) 1.00
## Bulk_ESS
## sd(Intercept) 4968
## sd(Prime_type) 4929
## sd(biasmatchmismatch) 3744
## sd(Agec) 3363
## sd(Prime_type:biasmatchmismatch) 4002
## sd(Prime_type:Agec) 4563
## sd(biasmatchmismatch:Agec) 4248
## sd(Prime_type:biasmatchmismatch:Agec) 3795
## cor(Intercept,Prime_type) 12215
## cor(Intercept,biasmatchmismatch) 11084
## cor(Prime_type,biasmatchmismatch) 8277
## cor(Intercept,Agec) 9530
## cor(Prime_type,Agec) 8029
## cor(biasmatchmismatch,Agec) 7645
## cor(Intercept,Prime_type:biasmatchmismatch) 12113
## cor(Prime_type,Prime_type:biasmatchmismatch) 8874
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) 8044
## cor(Agec,Prime_type:biasmatchmismatch) 7949
## cor(Intercept,Prime_type:Agec) 11753
## cor(Prime_type,Prime_type:Agec) 9907
## cor(biasmatchmismatch,Prime_type:Agec) 8856
## cor(Agec,Prime_type:Agec) 7834
## cor(Prime_type:biasmatchmismatch,Prime_type:Agec) 6715
## cor(Intercept,biasmatchmismatch:Agec) 13420
## cor(Prime_type,biasmatchmismatch:Agec) 10293
## cor(biasmatchmismatch,biasmatchmismatch:Agec) 7577
## cor(Agec,biasmatchmismatch:Agec) 7731
## cor(Prime_type:biasmatchmismatch,biasmatchmismatch:Agec) 6975
## cor(Prime_type:Agec,biasmatchmismatch:Agec) 5347
## cor(Intercept,Prime_type:biasmatchmismatch:Agec) 10977
## cor(Prime_type,Prime_type:biasmatchmismatch:Agec) 9306
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 8347
## cor(Agec,Prime_type:biasmatchmismatch:Agec) 6997
## cor(Prime_type:biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 6800
## cor(Prime_type:Agec,Prime_type:biasmatchmismatch:Agec) 5473
## cor(biasmatchmismatch:Agec,Prime_type:biasmatchmismatch:Agec) 5163
## Tail_ESS
## sd(Intercept) 5708
## sd(Prime_type) 4883
## sd(biasmatchmismatch) 4107
## sd(Agec) 3036
## sd(Prime_type:biasmatchmismatch) 3552
## sd(Prime_type:Agec) 4654
## sd(biasmatchmismatch:Agec) 4729
## sd(Prime_type:biasmatchmismatch:Agec) 5009
## cor(Intercept,Prime_type) 6530
## cor(Intercept,biasmatchmismatch) 6149
## cor(Prime_type,biasmatchmismatch) 5279
## cor(Intercept,Agec) 6231
## cor(Prime_type,Agec) 6071
## cor(biasmatchmismatch,Agec) 6427
## cor(Intercept,Prime_type:biasmatchmismatch) 6046
## cor(Prime_type,Prime_type:biasmatchmismatch) 5944
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) 6774
## cor(Agec,Prime_type:biasmatchmismatch) 6585
## cor(Intercept,Prime_type:Agec) 5293
## cor(Prime_type,Prime_type:Agec) 5325
## cor(biasmatchmismatch,Prime_type:Agec) 6397
## cor(Agec,Prime_type:Agec) 6255
## cor(Prime_type:biasmatchmismatch,Prime_type:Agec) 6579
## cor(Intercept,biasmatchmismatch:Agec) 5870
## cor(Prime_type,biasmatchmismatch:Agec) 5981
## cor(biasmatchmismatch,biasmatchmismatch:Agec) 6248
## cor(Agec,biasmatchmismatch:Agec) 6719
## cor(Prime_type:biasmatchmismatch,biasmatchmismatch:Agec) 6630
## cor(Prime_type:Agec,biasmatchmismatch:Agec) 6345
## cor(Intercept,Prime_type:biasmatchmismatch:Agec) 6286
## cor(Prime_type,Prime_type:biasmatchmismatch:Agec) 6155
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 5896
## cor(Agec,Prime_type:biasmatchmismatch:Agec) 6281
## cor(Prime_type:biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 6726
## cor(Prime_type:Agec,Prime_type:biasmatchmismatch:Agec) 6878
## cor(biasmatchmismatch:Agec,Prime_type:biasmatchmismatch:Agec) 6374
##
## Population-Level Effects:
## Estimate Est.Error l-90% CI u-90% CI Rhat
## Intercept -1.96 0.62 -2.95 -0.94 1.00
## Prime_type 0.48 0.53 -0.35 1.34 1.00
## biasmatchmismatch 0.45 0.49 -0.36 1.21 1.00
## Agec 0.54 0.28 0.11 1.00 1.00
## Prime_type:biasmatchmismatch 0.97 0.83 -0.30 2.36 1.00
## Prime_type:Agec 0.01 0.38 -0.59 0.60 1.00
## biasmatchmismatch:Agec -0.14 0.29 -0.59 0.33 1.00
## Prime_type:biasmatchmismatch:Agec -0.38 0.58 -1.29 0.54 1.00
## Bulk_ESS Tail_ESS
## Intercept 3566 3753
## Prime_type 5658 5426
## biasmatchmismatch 4764 4094
## Agec 4832 4560
## Prime_type:biasmatchmismatch 5606 4108
## Prime_type:Agec 6400 5299
## biasmatchmismatch:Agec 5775 4599
## Prime_type:biasmatchmismatch:Agec 6155 4991
##
## 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).mXsec3 <- add_criterion(mXsec3, "loo")#plot(mXsec3, N=4)pp_check(mXsec3)## Using 10 posterior draws for ppc type 'dens_overlay' by default.
Residuals vs fitted
mXsec3.model.check <- createDHARMa(
simulatedResponse = t(posterior_predict(mXsec3)),
observedResponse = XSectional3$Response,
fittedPredictedResponse = apply(t(posterior_epred(mXsec3)), 1, mean),
integerResponse=TRUE
)
plot(mXsec3.model.check) # standard plots
testDispersion(mXsec3.model.check)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.67471, p-value < 2.2e-16
## alternative hypothesis: two.sidedplotResiduals(mXsec3.model.check, form=XSectional3$Agec)
hypothesis(mXsec3, "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.48 0.53 -0.35 1.34 5.18 0.84
## 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(mXsec3, "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.97 0.83 -0.3 2.36 8.47
## Post.Prob Star
## 1 0.89
## ---
## '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(mXsec3, "Prime_type:Agec> 0")## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (Prime_type:Agec) > 0 0.01 0.38 -0.59 0.6 1.04
## Post.Prob Star
## 1 0.51
## ---
## '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(mXsec3, "Prime_type:biasmatchmismatch:Agec < 0") # Reduction in bias with age## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (Prime_type:biasm... < 0 -0.38 0.58 -1.29 0.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.XSectional3 %>%
mutate(
k = mXsec3$criteria$loo$diagnostics$pareto_k
) %>%
ggplot(aes(x=k)) + geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
XSectional3 %>%
mutate(
k = mXsec3$criteria$loo$diagnostics$pareto_k
) %>%
filter(k >= .7)labels = list("match" = "Match",
"mismatch" = "Mis-Match")
labeller <- function(variable,value){
return(labels[value])
}
df <- tibble( # Design matrix (minus intercept which is automatic)
biasmatch = factor(rep(c(0, 0, 1, 1), 7), labels = c("match", "mismatch")),
Prime_type = rep(c(-.5, .5, -.5, .5), 7),
Prime = rep(c("POD", "DOD", "POD", "DOD"), 7),
Agec = c(-3, -3, -3, -3, -2, -2, -2, -2, -1, -1, -1, -1, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3),
Age = c(3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9)
)
m3_post <- fitted(mXsec3, newdata=df, scale="linear", re_formula=NA, summary=FALSE, ndraws=25) %>%
t() %>%
as_tibble() %>%
bind_cols(df, .) %>%
pivot_longer(starts_with("V"), names_to="sample", values_to="linpred")
(posterior_plot3 <- m3_post %>%
mutate(
con =paste(sample, Prime)
) %>%
ggplot(aes(y=linpred, x=Age, fill=Prime, color=Prime)) + geom_line(aes(fill=Prime, color=Prime, group=con)) + facet_wrap(~biasmatch, labeller=labeller) + theme_minimal() + labs(y="Linear Predictor") + scale_colour_colorblind() + scale_fill_colorblind() + scale_colour_manual (values=cbPalette) + scale_fill_manual(values=cbPalette))## Warning in geom_line(aes(fill = Prime, color = Prime, group = con)): Ignoring
## unknown aesthetics: fill## Warning: The `labeller` API has been updated. Labellers taking `variable` and `value`
## arguments are now deprecated.
## ℹ See labellers documentation.## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
points_model3 <- XSectional %>%
filter(Overlap=="No Overlap") %>%
group_by(Part_No, Prime, biasmatch) %>%
mutate(
mean = mean(Response)
) %>%
slice(1) %>%
ungroup() %>%
dplyr::select(Age, Prime, biasmatch, mean)m3_fitted_response <- fitted(mXsec3, re_formula=NA, scale="response", probs=c(.05, .95)) %>%
cbind(XSectional3) %>%
dplyr::select(Age, Prime, biasmatch, Estimate, Q5, Q95)
m3_fitted_plot <- ggplot() +
geom_point(data= points_model3, aes(y=mean, x=Age, colour=Prime), shape=3) +
geom_line(data= m3_fitted_response , aes(y=Estimate, x=Age, fill=Prime, color=Prime, group=Prime), size=1.1) +
geom_ribbon(data= m3_fitted_response , aes(y=Estimate, ymin=Q5, ymax=Q95, x=Age, fill=Prime, group=Prime), alpha=0.5) +
facet_wrap(~biasmatch) +
theme_minimal() +
labs(y="Probability") +
theme(
strip.text.x = element_blank(),
strip.background = element_blank()
) + scale_colour_colorblind() + scale_fill_colorblind() + scale_colour_manual (values=cbPalette) + scale_fill_manual(values=cbPalette)## Warning in geom_line(data = m3_fitted_response, aes(y = Estimate, x = Age, :
## Ignoring unknown aesthetics: fill## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.posterior_plot3/m3_fitted_plot + plot_layout(guides = "collect") & theme(legend.position = "right") 
ggsave("plots/Fig3.png", last_plot(), dpi=500)## Saving 7 x 5 in image3.2 Only Participants Who Produced 1 DoD
XSectional4 <- filter(XSectional2, Overlap=="No Overlap")
mXsec4 <-brm(Response ~ 1 +Prime_type*biasmatch*Agec + (1 +Prime_type*biasmatch | Case_ID) + (1 +Prime_type*biasmatch*Agec | Trial_verb), logit_prior, family="bernoulli", iter=4000, cores=4, control=list(adapt_delta=.95), seed=12345,file="mXsec4", data=XSectional4)
mXsec4 <- add_criterion(mXsec4, "loo")#plot(mXsec4, N=4)mXsec4.model.check <- createDHARMa(
simulatedResponse = t(posterior_predict(mXsec4)),
observedResponse = XSectional4$Response,
fittedPredictedResponse = apply(t(posterior_epred(mXsec4)), 1, mean),
integerResponse=TRUE
)
plot(mXsec4.model.check) # standard plots
testDispersion(mXsec4.model.check)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.70578, p-value < 2.2e-16
## alternative hypothesis: two.sidedplotResiduals(mXsec4.model.check, form=XSectional4$Agec)
summary(mXsec4, prob = .9)## Family: bernoulli
## Links: mu = logit
## Formula: Response ~ 1 + Prime_type * biasmatch * Agec + (1 + Prime_type * biasmatch | Case_ID) + (1 + Prime_type * biasmatch * Agec | Trial_verb)
## Data: XSectional4 (Number of observations: 673)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Group-Level Effects:
## ~Case_ID (Number of levels: 114)
## Estimate Est.Error l-90% CI
## sd(Intercept) 1.35 0.24 0.98
## sd(Prime_type) 0.59 0.44 0.05
## sd(biasmatchmismatch) 0.53 0.37 0.05
## sd(Prime_type:biasmatchmismatch) 1.06 0.75 0.09
## cor(Intercept,Prime_type) -0.08 0.42 -0.74
## cor(Intercept,biasmatchmismatch) 0.04 0.41 -0.64
## cor(Prime_type,biasmatchmismatch) -0.03 0.44 -0.74
## cor(Intercept,Prime_type:biasmatchmismatch) -0.05 0.40 -0.69
## cor(Prime_type,Prime_type:biasmatchmismatch) -0.15 0.45 -0.81
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) 0.00 0.44 -0.72
## u-90% CI Rhat Bulk_ESS
## sd(Intercept) 1.76 1.00 3380
## sd(Prime_type) 1.42 1.00 2208
## sd(biasmatchmismatch) 1.23 1.00 1160
## sd(Prime_type:biasmatchmismatch) 2.48 1.01 1387
## cor(Intercept,Prime_type) 0.65 1.00 7433
## cor(Intercept,biasmatchmismatch) 0.71 1.00 6054
## cor(Prime_type,biasmatchmismatch) 0.70 1.00 2689
## cor(Intercept,Prime_type:biasmatchmismatch) 0.64 1.00 6343
## cor(Prime_type,Prime_type:biasmatchmismatch) 0.66 1.00 2913
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) 0.72 1.00 3853
## Tail_ESS
## sd(Intercept) 5245
## sd(Prime_type) 3726
## sd(biasmatchmismatch) 2928
## sd(Prime_type:biasmatchmismatch) 2888
## cor(Intercept,Prime_type) 5832
## cor(Intercept,biasmatchmismatch) 5669
## cor(Prime_type,biasmatchmismatch) 4910
## cor(Intercept,Prime_type:biasmatchmismatch) 4999
## cor(Prime_type,Prime_type:biasmatchmismatch) 4447
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) 6289
##
## ~Trial_verb (Number of levels: 6)
## Estimate
## sd(Intercept) 1.32
## sd(Prime_type) 0.55
## sd(biasmatchmismatch) 0.65
## sd(Agec) 0.39
## sd(Prime_type:biasmatchmismatch) 1.17
## sd(Prime_type:Agec) 0.48
## sd(biasmatchmismatch:Agec) 0.36
## sd(Prime_type:biasmatchmismatch:Agec) 0.91
## cor(Intercept,Prime_type) -0.10
## cor(Intercept,biasmatchmismatch) 0.08
## cor(Prime_type,biasmatchmismatch) -0.04
## cor(Intercept,Agec) 0.21
## cor(Prime_type,Agec) -0.05
## cor(biasmatchmismatch,Agec) 0.07
## cor(Intercept,Prime_type:biasmatchmismatch) -0.16
## cor(Prime_type,Prime_type:biasmatchmismatch) -0.03
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) -0.10
## cor(Agec,Prime_type:biasmatchmismatch) -0.06
## cor(Intercept,Prime_type:Agec) -0.05
## cor(Prime_type,Prime_type:Agec) -0.01
## cor(biasmatchmismatch,Prime_type:Agec) -0.02
## cor(Agec,Prime_type:Agec) -0.05
## cor(Prime_type:biasmatchmismatch,Prime_type:Agec) 0.06
## cor(Intercept,biasmatchmismatch:Agec) 0.01
## cor(Prime_type,biasmatchmismatch:Agec) 0.02
## cor(biasmatchmismatch,biasmatchmismatch:Agec) -0.01
## cor(Agec,biasmatchmismatch:Agec) -0.04
## cor(Prime_type:biasmatchmismatch,biasmatchmismatch:Agec) 0.01
## cor(Prime_type:Agec,biasmatchmismatch:Agec) -0.03
## cor(Intercept,Prime_type:biasmatchmismatch:Agec) 0.17
## cor(Prime_type,Prime_type:biasmatchmismatch:Agec) -0.05
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 0.10
## cor(Agec,Prime_type:biasmatchmismatch:Agec) 0.10
## cor(Prime_type:biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) -0.10
## cor(Prime_type:Agec,Prime_type:biasmatchmismatch:Agec) -0.12
## cor(biasmatchmismatch:Agec,Prime_type:biasmatchmismatch:Agec) 0.01
## Est.Error
## sd(Intercept) 0.51
## sd(Prime_type) 0.47
## sd(biasmatchmismatch) 0.51
## sd(Agec) 0.27
## sd(Prime_type:biasmatchmismatch) 0.75
## sd(Prime_type:Agec) 0.41
## sd(biasmatchmismatch:Agec) 0.32
## sd(Prime_type:biasmatchmismatch:Agec) 0.66
## cor(Intercept,Prime_type) 0.33
## cor(Intercept,biasmatchmismatch) 0.32
## cor(Prime_type,biasmatchmismatch) 0.33
## cor(Intercept,Agec) 0.31
## cor(Prime_type,Agec) 0.33
## cor(biasmatchmismatch,Agec) 0.33
## cor(Intercept,Prime_type:biasmatchmismatch) 0.31
## cor(Prime_type,Prime_type:biasmatchmismatch) 0.33
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) 0.33
## cor(Agec,Prime_type:biasmatchmismatch) 0.33
## cor(Intercept,Prime_type:Agec) 0.34
## cor(Prime_type,Prime_type:Agec) 0.34
## cor(biasmatchmismatch,Prime_type:Agec) 0.33
## cor(Agec,Prime_type:Agec) 0.34
## cor(Prime_type:biasmatchmismatch,Prime_type:Agec) 0.33
## cor(Intercept,biasmatchmismatch:Agec) 0.34
## cor(Prime_type,biasmatchmismatch:Agec) 0.34
## cor(biasmatchmismatch,biasmatchmismatch:Agec) 0.33
## cor(Agec,biasmatchmismatch:Agec) 0.33
## cor(Prime_type:biasmatchmismatch,biasmatchmismatch:Agec) 0.33
## cor(Prime_type:Agec,biasmatchmismatch:Agec) 0.34
## cor(Intercept,Prime_type:biasmatchmismatch:Agec) 0.32
## cor(Prime_type,Prime_type:biasmatchmismatch:Agec) 0.33
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 0.34
## cor(Agec,Prime_type:biasmatchmismatch:Agec) 0.33
## cor(Prime_type:biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 0.33
## cor(Prime_type:Agec,Prime_type:biasmatchmismatch:Agec) 0.34
## cor(biasmatchmismatch:Agec,Prime_type:biasmatchmismatch:Agec) 0.33
## l-90% CI
## sd(Intercept) 0.68
## sd(Prime_type) 0.04
## sd(biasmatchmismatch) 0.06
## sd(Agec) 0.05
## sd(Prime_type:biasmatchmismatch) 0.16
## sd(Prime_type:Agec) 0.03
## sd(biasmatchmismatch:Agec) 0.02
## sd(Prime_type:biasmatchmismatch:Agec) 0.09
## cor(Intercept,Prime_type) -0.63
## cor(Intercept,biasmatchmismatch) -0.46
## cor(Prime_type,biasmatchmismatch) -0.58
## cor(Intercept,Agec) -0.35
## cor(Prime_type,Agec) -0.59
## cor(biasmatchmismatch,Agec) -0.49
## cor(Intercept,Prime_type:biasmatchmismatch) -0.64
## cor(Prime_type,Prime_type:biasmatchmismatch) -0.57
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) -0.62
## cor(Agec,Prime_type:biasmatchmismatch) -0.59
## cor(Intercept,Prime_type:Agec) -0.60
## cor(Prime_type,Prime_type:Agec) -0.56
## cor(biasmatchmismatch,Prime_type:Agec) -0.57
## cor(Agec,Prime_type:Agec) -0.59
## cor(Prime_type:biasmatchmismatch,Prime_type:Agec) -0.50
## cor(Intercept,biasmatchmismatch:Agec) -0.55
## cor(Prime_type,biasmatchmismatch:Agec) -0.54
## cor(biasmatchmismatch,biasmatchmismatch:Agec) -0.55
## cor(Agec,biasmatchmismatch:Agec) -0.59
## cor(Prime_type:biasmatchmismatch,biasmatchmismatch:Agec) -0.54
## cor(Prime_type:Agec,biasmatchmismatch:Agec) -0.58
## cor(Intercept,Prime_type:biasmatchmismatch:Agec) -0.39
## cor(Prime_type,Prime_type:biasmatchmismatch:Agec) -0.58
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) -0.48
## cor(Agec,Prime_type:biasmatchmismatch:Agec) -0.46
## cor(Prime_type:biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) -0.62
## cor(Prime_type:Agec,Prime_type:biasmatchmismatch:Agec) -0.65
## cor(biasmatchmismatch:Agec,Prime_type:biasmatchmismatch:Agec) -0.54
## u-90% CI
## sd(Intercept) 2.27
## sd(Prime_type) 1.45
## sd(biasmatchmismatch) 1.62
## sd(Agec) 0.87
## sd(Prime_type:biasmatchmismatch) 2.56
## sd(Prime_type:Agec) 1.27
## sd(biasmatchmismatch:Agec) 0.98
## sd(Prime_type:biasmatchmismatch:Agec) 2.12
## cor(Intercept,Prime_type) 0.46
## cor(Intercept,biasmatchmismatch) 0.60
## cor(Prime_type,biasmatchmismatch) 0.52
## cor(Intercept,Agec) 0.70
## cor(Prime_type,Agec) 0.51
## cor(biasmatchmismatch,Agec) 0.59
## cor(Intercept,Prime_type:biasmatchmismatch) 0.38
## cor(Prime_type,Prime_type:biasmatchmismatch) 0.52
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) 0.46
## cor(Agec,Prime_type:biasmatchmismatch) 0.49
## cor(Intercept,Prime_type:Agec) 0.51
## cor(Prime_type,Prime_type:Agec) 0.56
## cor(biasmatchmismatch,Prime_type:Agec) 0.53
## cor(Agec,Prime_type:Agec) 0.51
## cor(Prime_type:biasmatchmismatch,Prime_type:Agec) 0.59
## cor(Intercept,biasmatchmismatch:Agec) 0.56
## cor(Prime_type,biasmatchmismatch:Agec) 0.57
## cor(biasmatchmismatch,biasmatchmismatch:Agec) 0.55
## cor(Agec,biasmatchmismatch:Agec) 0.52
## cor(Prime_type:biasmatchmismatch,biasmatchmismatch:Agec) 0.56
## cor(Prime_type:Agec,biasmatchmismatch:Agec) 0.53
## cor(Intercept,Prime_type:biasmatchmismatch:Agec) 0.66
## cor(Prime_type,Prime_type:biasmatchmismatch:Agec) 0.51
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 0.63
## cor(Agec,Prime_type:biasmatchmismatch:Agec) 0.62
## cor(Prime_type:biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 0.47
## cor(Prime_type:Agec,Prime_type:biasmatchmismatch:Agec) 0.47
## cor(biasmatchmismatch:Agec,Prime_type:biasmatchmismatch:Agec) 0.56
## Rhat
## sd(Intercept) 1.00
## sd(Prime_type) 1.00
## sd(biasmatchmismatch) 1.00
## sd(Agec) 1.00
## sd(Prime_type:biasmatchmismatch) 1.00
## sd(Prime_type:Agec) 1.00
## sd(biasmatchmismatch:Agec) 1.00
## sd(Prime_type:biasmatchmismatch:Agec) 1.00
## cor(Intercept,Prime_type) 1.00
## cor(Intercept,biasmatchmismatch) 1.00
## cor(Prime_type,biasmatchmismatch) 1.00
## cor(Intercept,Agec) 1.00
## cor(Prime_type,Agec) 1.00
## cor(biasmatchmismatch,Agec) 1.00
## cor(Intercept,Prime_type:biasmatchmismatch) 1.00
## cor(Prime_type,Prime_type:biasmatchmismatch) 1.00
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) 1.00
## cor(Agec,Prime_type:biasmatchmismatch) 1.00
## cor(Intercept,Prime_type:Agec) 1.00
## cor(Prime_type,Prime_type:Agec) 1.00
## cor(biasmatchmismatch,Prime_type:Agec) 1.00
## cor(Agec,Prime_type:Agec) 1.00
## cor(Prime_type:biasmatchmismatch,Prime_type:Agec) 1.00
## cor(Intercept,biasmatchmismatch:Agec) 1.00
## cor(Prime_type,biasmatchmismatch:Agec) 1.00
## cor(biasmatchmismatch,biasmatchmismatch:Agec) 1.00
## cor(Agec,biasmatchmismatch:Agec) 1.00
## cor(Prime_type:biasmatchmismatch,biasmatchmismatch:Agec) 1.00
## cor(Prime_type:Agec,biasmatchmismatch:Agec) 1.00
## cor(Intercept,Prime_type:biasmatchmismatch:Agec) 1.00
## cor(Prime_type,Prime_type:biasmatchmismatch:Agec) 1.00
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 1.00
## cor(Agec,Prime_type:biasmatchmismatch:Agec) 1.00
## cor(Prime_type:biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 1.00
## cor(Prime_type:Agec,Prime_type:biasmatchmismatch:Agec) 1.00
## cor(biasmatchmismatch:Agec,Prime_type:biasmatchmismatch:Agec) 1.00
## Bulk_ESS
## sd(Intercept) 3785
## sd(Prime_type) 4334
## sd(biasmatchmismatch) 3079
## sd(Agec) 2922
## sd(Prime_type:biasmatchmismatch) 3199
## sd(Prime_type:Agec) 3088
## sd(biasmatchmismatch:Agec) 3990
## sd(Prime_type:biasmatchmismatch:Agec) 2863
## cor(Intercept,Prime_type) 9639
## cor(Intercept,biasmatchmismatch) 9102
## cor(Prime_type,biasmatchmismatch) 7549
## cor(Intercept,Agec) 7930
## cor(Prime_type,Agec) 7068
## cor(biasmatchmismatch,Agec) 6867
## cor(Intercept,Prime_type:biasmatchmismatch) 10182
## cor(Prime_type,Prime_type:biasmatchmismatch) 7240
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) 6639
## cor(Agec,Prime_type:biasmatchmismatch) 6927
## cor(Intercept,Prime_type:Agec) 9588
## cor(Prime_type,Prime_type:Agec) 8990
## cor(biasmatchmismatch,Prime_type:Agec) 6700
## cor(Agec,Prime_type:Agec) 6009
## cor(Prime_type:biasmatchmismatch,Prime_type:Agec) 5940
## cor(Intercept,biasmatchmismatch:Agec) 9891
## cor(Prime_type,biasmatchmismatch:Agec) 8662
## cor(biasmatchmismatch,biasmatchmismatch:Agec) 7870
## cor(Agec,biasmatchmismatch:Agec) 7754
## cor(Prime_type:biasmatchmismatch,biasmatchmismatch:Agec) 6462
## cor(Prime_type:Agec,biasmatchmismatch:Agec) 5586
## cor(Intercept,Prime_type:biasmatchmismatch:Agec) 8070
## cor(Prime_type,Prime_type:biasmatchmismatch:Agec) 7668
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 6450
## cor(Agec,Prime_type:biasmatchmismatch:Agec) 6648
## cor(Prime_type:biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 5704
## cor(Prime_type:Agec,Prime_type:biasmatchmismatch:Agec) 4921
## cor(biasmatchmismatch:Agec,Prime_type:biasmatchmismatch:Agec) 5118
## Tail_ESS
## sd(Intercept) 5680
## sd(Prime_type) 4051
## sd(biasmatchmismatch) 3808
## sd(Agec) 2853
## sd(Prime_type:biasmatchmismatch) 3175
## sd(Prime_type:Agec) 3562
## sd(biasmatchmismatch:Agec) 3901
## sd(Prime_type:biasmatchmismatch:Agec) 3397
## cor(Intercept,Prime_type) 6464
## cor(Intercept,biasmatchmismatch) 6171
## cor(Prime_type,biasmatchmismatch) 6337
## cor(Intercept,Agec) 5586
## cor(Prime_type,Agec) 6564
## cor(biasmatchmismatch,Agec) 6404
## cor(Intercept,Prime_type:biasmatchmismatch) 6179
## cor(Prime_type,Prime_type:biasmatchmismatch) 6135
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) 5805
## cor(Agec,Prime_type:biasmatchmismatch) 6169
## cor(Intercept,Prime_type:Agec) 5618
## cor(Prime_type,Prime_type:Agec) 6237
## cor(biasmatchmismatch,Prime_type:Agec) 5833
## cor(Agec,Prime_type:Agec) 6717
## cor(Prime_type:biasmatchmismatch,Prime_type:Agec) 6307
## cor(Intercept,biasmatchmismatch:Agec) 6067
## cor(Prime_type,biasmatchmismatch:Agec) 6343
## cor(biasmatchmismatch,biasmatchmismatch:Agec) 6203
## cor(Agec,biasmatchmismatch:Agec) 6640
## cor(Prime_type:biasmatchmismatch,biasmatchmismatch:Agec) 6427
## cor(Prime_type:Agec,biasmatchmismatch:Agec) 5958
## cor(Intercept,Prime_type:biasmatchmismatch:Agec) 5306
## cor(Prime_type,Prime_type:biasmatchmismatch:Agec) 6233
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 6214
## cor(Agec,Prime_type:biasmatchmismatch:Agec) 6879
## cor(Prime_type:biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 5952
## cor(Prime_type:Agec,Prime_type:biasmatchmismatch:Agec) 6388
## cor(biasmatchmismatch:Agec,Prime_type:biasmatchmismatch:Agec) 6041
##
## Population-Level Effects:
## Estimate Est.Error l-90% CI u-90% CI Rhat
## Intercept -1.42 0.59 -2.40 -0.46 1.00
## Prime_type 0.49 0.51 -0.31 1.31 1.00
## biasmatchmismatch 0.53 0.43 -0.18 1.20 1.00
## Agec 0.44 0.26 0.02 0.86 1.00
## Prime_type:biasmatchmismatch 1.03 0.88 -0.33 2.50 1.00
## Prime_type:Agec -0.01 0.40 -0.64 0.61 1.00
## biasmatchmismatch:Agec -0.18 0.28 -0.62 0.28 1.00
## Prime_type:biasmatchmismatch:Agec -0.36 0.66 -1.38 0.67 1.00
## Bulk_ESS Tail_ESS
## Intercept 2506 3662
## Prime_type 5087 4663
## biasmatchmismatch 4151 4092
## Agec 4517 4529
## Prime_type:biasmatchmismatch 4399 4273
## Prime_type:Agec 4760 5052
## biasmatchmismatch:Agec 5402 4445
## Prime_type:biasmatchmismatch:Agec 4212 4052
##
## 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).pp_check(mXsec4)## Using 10 posterior draws for ppc type 'dens_overlay' by default.
hypothesis(mXsec4, "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.49 0.51 -0.31 1.31 5.39 0.84
## 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(mXsec4, "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 1.03 0.88 -0.33 2.5 8.56
## Post.Prob Star
## 1 0.9
## ---
## '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(mXsec4, "Prime_type:Agec> 0")## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (Prime_type:Agec) > 0 -0.01 0.4 -0.64 0.61 0.91
## Post.Prob Star
## 1 0.48
## ---
## '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(mXsec4, "Prime_type:biasmatchmismatch:Agec < 0") # Reduction in bias with age## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (Prime_type:biasm... < 0 -0.36 0.66 -1.38 0.67 2.84
## Post.Prob Star
## 1 0.74
## ---
## '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.labels = list("match" = "Match",
"mismatch" = "Mis-Match")
labeller <- function(variable,value){
return(labels[value])
}
df <- tibble( # Design matrix (minus intercept which is automatic)
biasmatch = factor(rep(c(0, 0, 1, 1), 7), labels = c("match", "mismatch")),
Prime_type = rep(c(-.5, .5, -.5, .5), 7),
Prime = rep(c("POD", "DOD", "POD", "DOD"), 7),
Agec = c(-3, -3, -3, -3, -2, -2, -2, -2, -1, -1, -1, -1, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3),
Age = c(3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9)
)
m4_post <- fitted(mXsec4, newdata=df, scale="linear", re_formula=NA, summary=FALSE, ndraws=25) %>%
t() %>%
as_tibble() %>%
bind_cols(df, .) %>%
pivot_longer(starts_with("V"), names_to="sample", values_to="linpred")
(posterior_plot4 <- m4_post %>%
mutate(
con =paste(sample, Prime)
) %>%
ggplot(aes(y=linpred, x=Age, fill=Prime, color=Prime)) + geom_line(aes(fill=Prime, color=Prime, group=con)) + facet_wrap(~biasmatch, labeller=labeller) + theme_minimal() + labs(y="Linear Predictor") +scale_colour_manual(values=cbPalette) +scale_fill_manual(values=cbPalette))## Warning in geom_line(aes(fill = Prime, color = Prime, group = con)): Ignoring
## unknown aesthetics: fill## Warning: The `labeller` API has been updated. Labellers taking `variable` and `value`
## arguments are now deprecated.
## ℹ See labellers documentation.
points_model4 <- XSectional2 %>%
filter(Overlap=="No Overlap") %>%
group_by(Part_No, Prime, biasmatch) %>%
mutate(
mean = mean(Response)
) %>%
slice(1) %>%
ungroup() %>%
dplyr::select(Age, Prime, biasmatch, mean)m4_fitted_response <- fitted(mXsec4, re_formula=NA, scale="response", probs=c(.05, .95)) %>%
cbind(XSectional4) %>%
dplyr::select(Age, Prime, biasmatch, Estimate, Q5, Q95)
m4_fitted_plot <- ggplot() +
geom_point(data= points_model3, aes(y=mean, x=Age, colour=Prime), shape=3) +
geom_line(data= m3_fitted_response , aes(y=Estimate, x=Age, fill=Prime, color=Prime, group=Prime), size=1.1) +
geom_ribbon(data= m3_fitted_response , aes(y=Estimate, ymin=Q5, ymax=Q95, x=Age, fill=Prime, group=Prime), alpha=0.5) +
facet_wrap(~biasmatch) +
theme_minimal() +
labs(y="Probability") +
theme(
strip.text.x = element_blank(),
strip.background = element_blank()
) +scale_colour_manual(values=cbPalette) +scale_fill_manual(values=cbPalette)## Warning in geom_line(data = m3_fitted_response, aes(y = Estimate, x = Age, :
## Ignoring unknown aesthetics: fillposterior_plot4/m4_fitted_plot + plot_layout(guides = "collect") & theme(legend.position = "right") 
3.3 Table
tab_model(mXsec1, mXsec2, mXsec3, mXsec4, show.ci=.90, show.re.var=FALSE, show.r2=FALSE, show.icc=FALSE, transform=NULL, collapse.ci=TRUE, show.p=TRUE, ci.hyphen = " : ", dv.labels=c("Model 1", "Model 2", "Model 3", "Model 4"), string.stat="", file="Table1", bpe="mean",
pred.labels=c(
"Intercept",
"Prime (DOD)",
"Overlap",
"Age.c",
"Prime * Overlap",
"Prime * Age.c",
"Overlap * Age.c",
"Prime * Overlap * Age.c",
"Bias Mismatch (POD)",
"Prime * Bias Mismatch ",
"Bias Mismatch * Age.c",
"Prime * Bias Mismatch * Age.c"
))| Model 1 | Model 2 | Model 3 | Model 4 | |
|---|---|---|---|---|
| Predictors | Log-Odds | Log-Odds | Log-Odds | Log-Odds |
| Intercept |
-1.61 (-2.59 : -0.59) |
-1.06 (-2.00 : -0.07) |
-1.96 (-2.95 : -0.94) |
-1.42 (-2.40 : -0.46) |
| Prime (DOD) |
1.05 (0.45 : 1.72) |
1.05 (0.42 : 1.72) |
0.48 (-0.35 : 1.34) |
0.49 (-0.31 : 1.31) |
| Overlap |
-0.66 (-1.22 : -0.10) |
-0.65 (-1.17 : -0.14) |
||
| Age.c |
0.45 (0.09 : 0.82) |
0.35 (0.01 : 0.68) |
0.54 (0.11 : 1.00) |
0.44 (0.02 : 0.86) |
| Prime * Overlap |
1.64 (0.73 : 2.69) |
1.68 (0.79 : 2.69) |
||
| Prime * Age.c |
-0.19 (-0.54 : 0.16) |
-0.21 (-0.57 : 0.15) |
0.01 (-0.59 : 0.60) |
-0.01 (-0.64 : 0.61) |
| Overlap * Age.c |
0.06 (-0.47 : 0.61) |
0.05 (-0.49 : 0.60) |
||
| Prime * Overlap * Age.c |
0.08 (-0.48 : 0.64) |
0.05 (-0.53 : 0.62) |
||
| Bias Mismatch (POD) |
0.45 (-0.36 : 1.21) |
0.53 (-0.18 : 1.20) |
||
| Prime * Bias Mismatch |
0.97 (-0.30 : 2.36) |
1.03 (-0.33 : 2.50) |
||
| Bias Mismatch * Age.c |
-0.14 (-0.59 : 0.33) |
-0.18 (-0.62 : 0.28) |
||
| Prime * Bias Mismatch * Age.c |
-0.38 (-1.29 : 0.54) |
-0.36 (-1.38 : 0.67) |
||
| N | 140 Case_ID | 114 Case_ID | 140 Case_ID | 114 Case_ID |
| 6 Trial_verb | 6 Trial_verb | 6 Trial_verb | 6 Trial_verb | |
| Observations | 1658 | 1348 | 827 | 673 |
4 Exploratory analysis of changes in verb bias across time
mXsec5 <- brm(Response ~ Agec + (1 | Case_ID) + (1 + Agec | Trial_verb), prior=logit_prior, family="bernoulli", iter=4000, cores=4, seed=12345, data=XSectional, file="mXsec5")
mXsec5 <- add_criterion(mXsec5, "loo")summary(mXsec5)## Family: bernoulli
## Links: mu = logit
## Formula: Response ~ Agec + (1 | Case_ID) + (1 + Agec | Trial_verb)
## Data: XSectional (Number of observations: 1658)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Group-Level Effects:
## ~Case_ID (Number of levels: 140)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 1.44 0.14 1.17 1.74 1.00 3175 5103
##
## ~Trial_verb (Number of levels: 6)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 1.12 0.46 0.54 2.26 1.00 3287 4681
## sd(Agec) 0.15 0.13 0.01 0.48 1.00 2641 4050
## cor(Intercept,Agec) 0.03 0.52 -0.89 0.91 1.00 10737 4965
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -1.55 0.51 -2.53 -0.49 1.00 2482 4244
## Agec 0.38 0.14 0.11 0.65 1.00 4004 4792
##
## 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).This code produces model predicted values for each verb conditional on time (for the figure below), using the fixed and (by item) random effects in the model above. It takes a while to run through the loop so I have saved the data and commented most of this out
#Agec=c(-3, -2, -1, 0, 1, 2, 3)
#Age = c(3, 4, 5, 6, 7, 8, 9)
#Trial_verb = c(rep("brought", 7), rep("gave", 7), rep("passed", 7), rep("sent", 7), rep("showed", 7), rep("threw", 7))
#newdata <- tibble(
# Agec = rep(Agec, 6),
# Age = rep(Age, 6),
# Trial_verb = c(rep("brought", 7), rep("gave", 7), rep("passed", 7), rep("sent", 7), rep("showed", 7), rep("threw", 7))
#)
#posterior <- posterior_samples(mXsec5) %>%
# dplyr::select("b_Intercept", "b_Agec", starts_with("r_Trial_verb")) %>%
# dplyr::rename("brought_Intercept" = "r_Trial_verb[brought,Intercept]", "gave_Intercept" = "r_Trial_verb[gave,Intercept]", "passed_Intercept" = "r_Trial_verb[passed,Intercept]", "sent_Intercept" = "r_Trial_verb[sent,Intercept]", "showed_Intercept" = "r_Trial_verb[showed,Intercept]", "threw_Intercept" = "r_Trial_verb[threw,Intercept]","brought_Agec" = "r_Trial_verb[brought,Agec]", "gave_Agec" = "r_Trial_verb[gave,Agec]", "passed_Agec" = "r_Trial_verb[passed,Agec]", "sent_Agec" = "r_Trial_verb[sent,Agec]", "showed_Agec" = "r_Trial_verb[showed,Agec]", "threw_Agec" = "r_Trial_verb[threw,Agec]") %>%
# mutate(
# sample = paste("Samp", row_number(), sep="_")
# ) %>%
# relocate(sample) %>%
# pivot_longer(-"sample", names_to = "condition", values_to="value") %>%
# separate(condition, c("condition", "effect"), "_")
#uni <- unique(posterior$sample)
#for(i in 1:8000){
# pos <- filter(posterior, sample == uni[i])
# df<- tibble(
# Fixed = filter(pos, condition=="b" & effect=="Intercept")$value + filter(pos, condition=="b" & effect=="Agec")$value * Agec,
# Brought = Fixed + filter(pos, condition=="brought" & effect=="Intercept")$value + filter(pos, condition=="brought" & effect=="Agec")$value* Agec,
# Gave = Fixed + filter(pos, condition=="gave" & effect=="Intercept")$value + filter(pos, condition=="gave" & effect=="Agec")$value * Agec,
# Passed = Fixed + filter(pos, condition=="passed" & effect=="Intercept")$value + filter(pos, condition=="passed" & effect=="Agec")$value * Agec,
# Sent = Fixed + filter(pos, condition=="sent" & effect=="Intercept")$value + filter(pos, condition=="sent" & effect=="Agec")$value * Agec,
# Showed = Fixed + filter(pos, condition=="showed" & effect=="Intercept")$value + filter(pos, condition=="showed" & effect=="Agec")$value * Agec,
# Threw = Fixed + filter(pos, condition=="threw" & effect=="Intercept")$value + filter(pos, condition=="threw" & effect=="Agec")$value * Agec
# )
# out <- tibble(X = c(df$Brought, df$Gave, df$Passed, df$Sent, df$Showed, df$Threw))
# names(out)[1] <- pos$sample[1]
# newdata = tibble(newdata, out)
#}
#saveRDS(newdata, "plots/by_word.rds")
newdata <- readRDS("plots/by_word.rds")points <- XSectional %>%
group_by(Part_No, Trial_verb) %>%
mutate(
mean= mean(Response)
) %>%
slice(1) %>%
ungroup() %>%
dplyr::select(Age, Trial_verb, mean)
summ <- newdata %>%
pivot_longer(-c("Agec", "Age", "Trial_verb"), names_to = "sample", values_to="value") %>%
mutate(
prob = inv_logit_scaled(value)
) %>%
group_by(Age, Trial_verb) %>%
summarise(
median = median(prob),
Q95 = quantile(prob, .95),
Q5 = quantile(prob, .05)
) %>%
ungroup()## `summarise()` has grouped output by 'Age'. You can override using the `.groups`
## argument.ggplot() +
geom_line(data=summ, aes(x=Age, y = median, ymin=Q5, ymax=Q95, fill=Trial_verb, color=Trial_verb, linetype=Trial_verb, group=Trial_verb), size=1.25) +
geom_hline(yintercept=.5, linetype="dotted") +
theme_minimal() +
#geom_point(data=points, aes(x=Age, y=mean, fill=Trial_verb, color=Trial_verb), shape=3) +
ylab("Proportion DODs Produced") +
labs(color="Verb", fill="Verb", linetype="Verb") +
scale_x_continuous(breaks = round(seq(min(summ$Age), max(summ$Age), by = 1),1)) + scale_colour_colorblind() + scale_fill_colorblind()## Warning in geom_line(data = summ, aes(x = Age, y = median, ymin = Q5, ymax =
## Q95, : Ignoring unknown aesthetics: ymin, ymax, and fill
ggsave("plots/Fig5.png", last_plot(), dpi=500)## Saving 7 x 5 in image