Analysis of Cross Sectional Dative Priming
Longitudinal Analyses: 42 months, 48 months, 54 months, Longitudinal
library(brms)
library(readxl)
library(rmarkdown)
library(tidybayes)
library(ggplot2)
library(tidyr)
library(dplyr)
library(GGally)
library(DHARMa)
library(ggpubr)
library(tidybayes)
library(patchwork)
library(sjPlot)
XSectional <- read_excel("DOD Priming Data Sheet_Wgender_May42020 (1).xlsx")1 Some preliminary stuff
1.1 Missing data
Invalid trials
XSectional %>%
filter(trial_drop==1) Total invalid trial types. Total number of invalid trials per participant.
XSectional %>%
filter(trial_drop==1) %>%
summarise(
A_Error = sum(Admin.Error, na.rm = TRUE),
V_Error = sum(Verb.error, na.rm = TRUE),
Non_allow = sum(Non.allowable.prompt==TRUE, na.rm = TRUE),
N_Targ = sum(Non.Target.Response == TRUE, na.rm = TRUE),
Attn = sum(Attentiveness == TRUE, na.rm = TRUE)
)XSectional %>%
group_by(Part_No) %>%
mutate(
n_miss = sum(trial_drop, na.rm = TRUE)
) %>%
slice(1) %>%
ungroup(Part_No) %>%
filter(n_miss > 0) %>%
dplyr::select(Case_ID, List, trial_drop, n_miss)Drop invalid-trials plus invalid participants.
XSectional <- subset(XSectional, drop != 1)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).Add loo for loo_compare
mXsec0<- add_criterion(mXsec0, "loo")
mXsec0b<- add_criterion(mXsec0b, "loo")loo(mXsec0, mXsec0b)## Output of model 'mXsec0':
##
## Computed from 8000 by 1658 log-likelihood matrix
##
## Estimate SE
## elpd_loo -737.2 22.9
## p_loo 112.0 4.6
## looic 1474.3 45.7
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 1656 99.9% 2358
## (0.5, 0.7] (ok) 2 0.1% 2975
## (0.7, 1] (bad) 0 0.0% <NA>
## (1, Inf) (very bad) 0 0.0% <NA>
##
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'mXsec0b':
##
## Computed from 8000 by 1658 log-likelihood matrix
##
## Estimate SE
## elpd_loo -735.2 22.8
## p_loo 101.3 4.2
## looic 1470.5 45.7
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 1657 99.9% 3163
## (0.5, 0.7] (ok) 1 0.1% 3893
## (0.7, 1] (bad) 0 0.0% <NA>
## (1, Inf) (very bad) 0 0.0% <NA>
##
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.
##
## Model comparisons:
## elpd_diff se_diff
## mXsec0b 0.0 0.0
## mXsec0 -1.9 0.4model_weights(mXsec0, mXsec0b)## Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.
## Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.## mXsec0 mXsec0b
## 6.685358e-07 9.999993e-01M0b 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 plot the random intercepts and the average number of double object datives produced by participant.
df <- mXsec0b %>%
spread_draws(b_Intercept, r_Case_ID[Case_ID,Intercept])## Warning: `gather_()` was deprecated in tidyr 1.2.0.
## i Please use `gather()` instead.
## i The deprecated feature was likely used in the tidybayes package.
## Please report the issue at <https://github.com/mjskay/tidybayes/issues/new>.df <- df %>%
mutate(
Part_Mean = inv_logit_scaled(b_Intercept + r_Case_ID)
)
df_part <- df %>%
group_by(Case_ID) %>%
summarise(
Predict_Mean = mean(Part_Mean),
Predict_SD = sd(Part_Mean)
)
df_raw_part <- XSectional %>%
group_by(Case_ID) %>%
summarise(
Observed_Mean = mean(Response, na.rm=TRUE)
)
df_combined <- merge(df_part, df_raw_part, by="Case_ID")
hist(df_combined$Observed_Mean)
hist(df_combined$Predict_Mean )
There are a lot of 0s because some children did not produce a single DoD in this task.
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
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.
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.
## i Using compatibility `.name_repair`.(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"))## Warning: 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.5) +
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()
) ## Warning: Ignoring unknown aesthetics: fill, ymin, ymaxposterior_plot1/m1_fitted_plot + plot_layout(guides = "collect") & theme(legend.position = "right") 
The 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.
summary(mXsec1b , 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 u-90% CI Rhat
## sd(Intercept) 1.76 0.17 1.49 2.05 1.00
## sd(Prime_type) 0.49 0.31 0.04 1.03 1.00
## sd(OverlapOverlap) 0.37 0.25 0.03 0.83 1.00
## sd(Prime_type:OverlapOverlap) 0.94 0.50 0.12 1.77 1.01
## Bulk_ESS Tail_ESS
## sd(Intercept) 2974 5136
## sd(Prime_type) 1237 2238
## sd(OverlapOverlap) 1341 3046
## sd(Prime_type:OverlapOverlap) 1089 1799
##
## ~Trial_verb (Number of levels: 6)
## Estimate Est.Error l-90% CI u-90% CI Rhat
## sd(Intercept) 1.39 0.51 0.77 2.38 1.00
## sd(Prime_type) 0.53 0.38 0.06 1.26 1.00
## sd(OverlapOverlap) 0.46 0.33 0.06 1.06 1.00
## sd(Agec) 0.34 0.23 0.06 0.75 1.00
## sd(Prime_type:OverlapOverlap) 0.75 0.55 0.06 1.80 1.00
## sd(Prime_type:Agec) 0.21 0.20 0.01 0.58 1.00
## sd(OverlapOverlap:Agec) 0.64 0.35 0.24 1.28 1.00
## sd(Prime_type:OverlapOverlap:Agec) 0.39 0.35 0.02 1.06 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 3104 5026
## sd(Prime_type) 2157 1995
## sd(OverlapOverlap) 2614 2846
## sd(Agec) 1642 1921
## sd(Prime_type:OverlapOverlap) 2357 3010
## sd(Prime_type:Agec) 3741 3778
## sd(OverlapOverlap:Agec) 2032 2909
## sd(Prime_type:OverlapOverlap:Agec) 2713 2921
##
## Population-Level Effects:
## Estimate Est.Error l-90% CI u-90% CI Rhat
## Intercept -1.62 0.60 -2.57 -0.63 1.00
## Prime_type 1.00 0.35 0.48 1.57 1.00
## OverlapOverlap -0.55 0.30 -1.01 -0.06 1.00
## Agec 0.46 0.22 0.11 0.81 1.00
## Prime_type:OverlapOverlap 1.48 0.53 0.68 2.38 1.00
## Prime_type:Agec -0.18 0.21 -0.51 0.16 1.00
## OverlapOverlap:Agec 0.06 0.33 -0.43 0.58 1.00
## Prime_type:OverlapOverlap:Agec 0.08 0.35 -0.48 0.64 1.00
## Bulk_ESS Tail_ESS
## Intercept 2061 3050
## Prime_type 4386 4479
## OverlapOverlap 3783 3430
## Agec 2430 3566
## Prime_type:OverlapOverlap 4360 3301
## Prime_type:Agec 6507 5375
## OverlapOverlap:Agec 2546 3359
## Prime_type:OverlapOverlap:Agec 5077 4466
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).hypothesis(mXsec1b, "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.35 -0.48 0.64 1.48
## Post.Prob Star
## 1 0.6
## ---
## '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(mXsec1b, "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.18 0.21 -0.51 0.16 4.33
## Post.Prob Star
## 1 0.81
## ---
## '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(mXsec1b, "Agec > 0")## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
## 1 (Agec) > 0 0.46 0.22 0.11 0.81 42.24 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(mXsec1b, "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.48 0.53 0.68 2.38 180.82
## 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(mXsec1b, "Prime_type > 0")## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
## 1 (Prime_type) > 0 1 0.35 0.48 1.57 199 1
## 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.It’s also likely that I’ve overparameterized the model by including random effects for age and its interaction. I’ll try removing those.
mXsec1bb <- brm(Response ~ 1 + Prime_type*Overlap*Agec + (1 + Prime_type*Overlap| Case_ID) + (1 + Prime_type*Overlap | Trial_verb), control = list(adapt_delta=.98), logit_prior, family="bernoulli", iter=4000, cores=4, seed=12345, data=XSectional, file="mXsec1bb")
m1bb <- add_criterion(mXsec1bb, "loo")mXsec1bb.model.check <- createDHARMa(
simulatedResponse = t(posterior_predict(mXsec1bb)),
observedResponse = XSectional$Response,
fittedPredictedResponse = apply(t(posterior_epred(mXsec1bb)), 1, mean),
integerResponse=TRUE
)
plot(mXsec1bb.model.check) # standard plots
testDispersion(mXsec1bb.model.check)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.79468, p-value < 2.2e-16
## alternative hypothesis: two.sidedplotResiduals(mXsec1bb.model.check, form=XSectional$Agec)
Interesting to see it’s still not less underdispersed.
Let’s look at the results.
summary(mXsec1bb, prob = .9)## Family: bernoulli
## Links: mu = logit
## Formula: Response ~ 1 + Prime_type * Overlap * Agec + (1 + Prime_type * Overlap | Case_ID) + (1 + Prime_type * Overlap | 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.67 0.19 1.38
## sd(Prime_type) 0.51 0.33 0.05
## sd(OverlapOverlap) 0.40 0.27 0.04
## sd(Prime_type:OverlapOverlap) 0.95 0.54 0.11
## cor(Intercept,Prime_type) -0.24 0.39 -0.79
## cor(Intercept,OverlapOverlap) 0.29 0.40 -0.47
## cor(Prime_type,OverlapOverlap) -0.11 0.43 -0.77
## cor(Intercept,Prime_type:OverlapOverlap) -0.15 0.35 -0.69
## cor(Prime_type,Prime_type:OverlapOverlap) 0.04 0.43 -0.66
## cor(OverlapOverlap,Prime_type:OverlapOverlap) -0.10 0.44 -0.77
## u-90% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 1.99 1.00 2857 4821
## sd(Prime_type) 1.10 1.00 1905 3571
## sd(OverlapOverlap) 0.88 1.00 1460 2687
## sd(Prime_type:OverlapOverlap) 1.89 1.00 1351 2591
## cor(Intercept,Prime_type) 0.47 1.00 5803 4779
## cor(Intercept,OverlapOverlap) 0.84 1.00 5219 5433
## cor(Prime_type,OverlapOverlap) 0.63 1.00 3100 5135
## cor(Intercept,Prime_type:OverlapOverlap) 0.48 1.00 6409 5029
## cor(Prime_type,Prime_type:OverlapOverlap) 0.73 1.00 2521 4769
## cor(OverlapOverlap,Prime_type:OverlapOverlap) 0.67 1.00 2416 4532
##
## ~Trial_verb (Number of levels: 6)
## Estimate Est.Error l-90% CI
## sd(Intercept) 1.39 0.50 0.80
## sd(Prime_type) 0.58 0.42 0.06
## sd(OverlapOverlap) 0.51 0.32 0.08
## sd(Prime_type:OverlapOverlap) 0.74 0.56 0.07
## cor(Intercept,Prime_type) -0.30 0.40 -0.85
## cor(Intercept,OverlapOverlap) -0.41 0.37 -0.89
## cor(Prime_type,OverlapOverlap) 0.04 0.43 -0.67
## cor(Intercept,Prime_type:OverlapOverlap) -0.16 0.41 -0.77
## cor(Prime_type,Prime_type:OverlapOverlap) 0.04 0.44 -0.69
## cor(OverlapOverlap,Prime_type:OverlapOverlap) -0.08 0.43 -0.75
## u-90% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 2.34 1.00 2664 4894
## sd(Prime_type) 1.35 1.00 2764 3093
## sd(OverlapOverlap) 1.10 1.00 2366 2302
## sd(Prime_type:OverlapOverlap) 1.80 1.00 2876 3125
## cor(Intercept,Prime_type) 0.44 1.00 6438 5577
## cor(Intercept,OverlapOverlap) 0.31 1.00 5653 6014
## cor(Prime_type,OverlapOverlap) 0.73 1.00 6241 6034
## cor(Intercept,Prime_type:OverlapOverlap) 0.57 1.00 7643 5820
## cor(Prime_type,Prime_type:OverlapOverlap) 0.74 1.00 5721 6585
## cor(OverlapOverlap,Prime_type:OverlapOverlap) 0.65 1.00 5951 6337
##
## Population-Level Effects:
## Estimate Est.Error l-90% CI u-90% CI Rhat
## Intercept -1.61 0.62 -2.60 -0.61 1.00
## Prime_type 1.04 0.39 0.47 1.68 1.00
## OverlapOverlap -0.59 0.34 -1.15 -0.03 1.00
## Agec 0.46 0.14 0.23 0.69 1.00
## Prime_type:OverlapOverlap 1.55 0.59 0.68 2.54 1.00
## Prime_type:Agec -0.22 0.17 -0.51 0.06 1.00
## OverlapOverlap:Agec -0.01 0.13 -0.23 0.21 1.00
## Prime_type:OverlapOverlap:Agec 0.14 0.26 -0.29 0.57 1.00
## Bulk_ESS Tail_ESS
## Intercept 2148 3012
## Prime_type 3517 2806
## OverlapOverlap 3266 3825
## Agec 1869 4170
## Prime_type:OverlapOverlap 3894 4532
## Prime_type:Agec 6641 6234
## OverlapOverlap:Agec 8058 5656
## Prime_type:OverlapOverlap:Agec 6928 6045
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).hypothesis(mXsec1bb, "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.14 0.26 -0.29 0.57 2.4
## Post.Prob Star
## 1 0.71
## ---
## '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(mXsec1bb, "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.22 0.17 -0.51 0.06 8.95
## 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(mXsec1bb, "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.55 0.59 0.68 2.54 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(mXsec1bb, "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.04 0.39 0.47 1.68 128.03 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.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"))## Warning: 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.5) +
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()
) ## Warning: Ignoring unknown aesthetics: fill, ymin, ymaxposterior_plot2/m2_fitted_plot + plot_layout(guides = "collect") & theme(legend.position = "right") 
ggsave("plots/Fig3.png", last_plot(), dpi=500)## Saving 7 x 5 in image# One of the pareto k values was exactly .7. I've re-fit the model without that observation. Nothing important changed, but here is the code for estimating this model.
XSectional2 %>%
mutate(
k = mXsec2$criteria$loo$diagnostics$pareto_k
) %>%
ggplot(aes(x=k)) + geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
XSectional2 %>%
mutate(
k = mXsec2$criteria$loo$diagnostics$pareto_k
) %>%
filter(k < .7) %>%
brm(Response ~ 1 + Prime_type*Overlap*Agec + (1 + Prime_type*Overlap| Case_ID) + (1 + Prime_type*Overlap*Agec | Trial_verb), prior=logit_prior, family="bernoulli", iter=4000, cores=4, control=list(adapt_delta=.95), seed=12345, data=., file="m2_no_low_k") %>%
summary(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: . (Number of observations: 1347)
## 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.55 0.34 0.06
## sd(OverlapOverlap) 0.43 0.27 0.05
## sd(Prime_type:OverlapOverlap) 1.02 0.57 0.14
## cor(Intercept,Prime_type) -0.34 0.36 -0.82
## cor(Intercept,OverlapOverlap) 0.33 0.37 -0.37
## cor(Prime_type,OverlapOverlap) -0.16 0.43 -0.80
## cor(Intercept,Prime_type:OverlapOverlap) -0.25 0.34 -0.75
## cor(Prime_type,Prime_type:OverlapOverlap) 0.05 0.42 -0.65
## cor(OverlapOverlap,Prime_type:OverlapOverlap) -0.14 0.43 -0.79
## u-90% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 1.57 1.00 3243 5385
## sd(Prime_type) 1.14 1.00 1562 2518
## sd(OverlapOverlap) 0.92 1.00 1458 3108
## sd(Prime_type:OverlapOverlap) 1.98 1.00 1107 2844
## cor(Intercept,Prime_type) 0.36 1.00 4690 4469
## cor(Intercept,OverlapOverlap) 0.85 1.00 4928 4957
## cor(Prime_type,OverlapOverlap) 0.59 1.00 3148 5373
## cor(Intercept,Prime_type:OverlapOverlap) 0.38 1.00 5025 4761
## cor(Prime_type,Prime_type:OverlapOverlap) 0.74 1.00 2551 4116
## cor(OverlapOverlap,Prime_type:OverlapOverlap) 0.62 1.00 2518 4893
##
## ~Trial_verb (Number of levels: 6)
## Estimate
## sd(Intercept) 1.40
## sd(Prime_type) 0.67
## sd(OverlapOverlap) 0.50
## sd(Agec) 0.36
## sd(Prime_type:OverlapOverlap) 0.77
## sd(Prime_type:Agec) 0.23
## sd(OverlapOverlap:Agec) 0.66
## sd(Prime_type:OverlapOverlap:Agec) 0.45
## cor(Intercept,Prime_type) -0.18
## cor(Intercept,OverlapOverlap) -0.18
## cor(Prime_type,OverlapOverlap) -0.02
## cor(Intercept,Agec) 0.16
## cor(Prime_type,Agec) -0.02
## cor(OverlapOverlap,Agec) -0.14
## cor(Intercept,Prime_type:OverlapOverlap) -0.11
## cor(Prime_type,Prime_type:OverlapOverlap) 0.02
## cor(OverlapOverlap,Prime_type:OverlapOverlap) -0.06
## cor(Agec,Prime_type:OverlapOverlap) 0.00
## 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.26
## cor(Prime_type,OverlapOverlap:Agec) -0.00
## cor(OverlapOverlap,OverlapOverlap:Agec) 0.13
## cor(Agec,OverlapOverlap:Agec) -0.21
## cor(Prime_type:OverlapOverlap,OverlapOverlap:Agec) 0.00
## cor(Prime_type:Agec,OverlapOverlap:Agec) 0.04
## cor(Intercept,Prime_type:OverlapOverlap:Agec) -0.06
## cor(Prime_type,Prime_type:OverlapOverlap:Agec) -0.01
## cor(OverlapOverlap,Prime_type:OverlapOverlap:Agec) 0.04
## cor(Agec,Prime_type:OverlapOverlap:Agec) 0.02
## cor(Prime_type:OverlapOverlap,Prime_type:OverlapOverlap:Agec) 0.01
## cor(Prime_type:Agec,Prime_type:OverlapOverlap:Agec) -0.03
## cor(OverlapOverlap:Agec,Prime_type:OverlapOverlap:Agec) 0.01
## Est.Error
## sd(Intercept) 0.50
## sd(Prime_type) 0.45
## sd(OverlapOverlap) 0.34
## sd(Agec) 0.24
## sd(Prime_type:OverlapOverlap) 0.59
## sd(Prime_type:Agec) 0.22
## sd(OverlapOverlap:Agec) 0.35
## sd(Prime_type:OverlapOverlap:Agec) 0.41
## cor(Intercept,Prime_type) 0.31
## cor(Intercept,OverlapOverlap) 0.31
## cor(Prime_type,OverlapOverlap) 0.33
## cor(Intercept,Agec) 0.31
## cor(Prime_type,Agec) 0.32
## cor(OverlapOverlap,Agec) 0.33
## cor(Intercept,Prime_type:OverlapOverlap) 0.32
## 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.33
## 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.33
## 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.33
## 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.80 2.36
## sd(Prime_type) 0.09 1.51
## sd(OverlapOverlap) 0.06 1.15
## sd(Agec) 0.08 0.79
## sd(Prime_type:OverlapOverlap) 0.06 1.90
## sd(Prime_type:Agec) 0.02 0.64
## sd(OverlapOverlap:Agec) 0.25 1.33
## sd(Prime_type:OverlapOverlap:Agec) 0.03 1.26
## cor(Intercept,Prime_type) -0.65 0.36
## cor(Intercept,OverlapOverlap) -0.67 0.36
## cor(Prime_type,OverlapOverlap) -0.55 0.53
## cor(Intercept,Agec) -0.38 0.65
## cor(Prime_type,Agec) -0.54 0.50
## cor(OverlapOverlap,Agec) -0.64 0.41
## cor(Intercept,Prime_type:OverlapOverlap) -0.62 0.45
## cor(Prime_type,Prime_type:OverlapOverlap) -0.52 0.56
## cor(OverlapOverlap,Prime_type:OverlapOverlap) -0.60 0.49
## cor(Agec,Prime_type:OverlapOverlap) -0.52 0.53
## cor(Intercept,Prime_type:Agec) -0.55 0.53
## cor(Prime_type,Prime_type:Agec) -0.57 0.53
## cor(OverlapOverlap,Prime_type:Agec) -0.53 0.56
## cor(Agec,Prime_type:Agec) -0.55 0.55
## cor(Prime_type:OverlapOverlap,Prime_type:Agec) -0.57 0.55
## cor(Intercept,OverlapOverlap:Agec) -0.70 0.25
## cor(Prime_type,OverlapOverlap:Agec) -0.53 0.52
## cor(OverlapOverlap,OverlapOverlap:Agec) -0.41 0.63
## cor(Agec,OverlapOverlap:Agec) -0.69 0.33
## cor(Prime_type:OverlapOverlap,OverlapOverlap:Agec) -0.53 0.53
## cor(Prime_type:Agec,OverlapOverlap:Agec) -0.51 0.58
## cor(Intercept,Prime_type:OverlapOverlap:Agec) -0.58 0.48
## cor(Prime_type,Prime_type:OverlapOverlap:Agec) -0.56 0.54
## cor(OverlapOverlap,Prime_type:OverlapOverlap:Agec) -0.52 0.59
## cor(Agec,Prime_type:OverlapOverlap:Agec) -0.53 0.55
## cor(Prime_type:OverlapOverlap,Prime_type:OverlapOverlap:Agec) -0.54 0.55
## cor(Prime_type:Agec,Prime_type:OverlapOverlap:Agec) -0.58 0.54
## cor(OverlapOverlap:Agec,Prime_type:OverlapOverlap:Agec) -0.54 0.54
## Rhat Bulk_ESS
## sd(Intercept) 1.00 2466
## sd(Prime_type) 1.00 2557
## sd(OverlapOverlap) 1.00 2990
## sd(Agec) 1.00 2389
## sd(Prime_type:OverlapOverlap) 1.00 2941
## sd(Prime_type:Agec) 1.00 4553
## sd(OverlapOverlap:Agec) 1.00 3456
## sd(Prime_type:OverlapOverlap:Agec) 1.00 3436
## cor(Intercept,Prime_type) 1.00 6519
## cor(Intercept,OverlapOverlap) 1.00 7827
## cor(Prime_type,OverlapOverlap) 1.00 8048
## cor(Intercept,Agec) 1.00 6898
## cor(Prime_type,Agec) 1.00 6069
## cor(OverlapOverlap,Agec) 1.00 5808
## cor(Intercept,Prime_type:OverlapOverlap) 1.00 7683
## cor(Prime_type,Prime_type:OverlapOverlap) 1.00 8035
## cor(OverlapOverlap,Prime_type:OverlapOverlap) 1.00 6081
## cor(Agec,Prime_type:OverlapOverlap) 1.00 6868
## cor(Intercept,Prime_type:Agec) 1.00 10278
## cor(Prime_type,Prime_type:Agec) 1.00 10031
## cor(OverlapOverlap,Prime_type:Agec) 1.00 8005
## cor(Agec,Prime_type:Agec) 1.00 7232
## cor(Prime_type:OverlapOverlap,Prime_type:Agec) 1.00 5610
## cor(Intercept,OverlapOverlap:Agec) 1.00 6644
## cor(Prime_type,OverlapOverlap:Agec) 1.00 5824
## cor(OverlapOverlap,OverlapOverlap:Agec) 1.00 5579
## cor(Agec,OverlapOverlap:Agec) 1.00 6215
## cor(Prime_type:OverlapOverlap,OverlapOverlap:Agec) 1.00 6105
## cor(Prime_type:Agec,OverlapOverlap:Agec) 1.00 4511
## cor(Intercept,Prime_type:OverlapOverlap:Agec) 1.00 10720
## cor(Prime_type,Prime_type:OverlapOverlap:Agec) 1.00 8369
## cor(OverlapOverlap,Prime_type:OverlapOverlap:Agec) 1.00 8166
## cor(Agec,Prime_type:OverlapOverlap:Agec) 1.00 7360
## cor(Prime_type:OverlapOverlap,Prime_type:OverlapOverlap:Agec) 1.00 6458
## cor(Prime_type:Agec,Prime_type:OverlapOverlap:Agec) 1.00 5293
## cor(OverlapOverlap:Agec,Prime_type:OverlapOverlap:Agec) 1.00 6293
## Tail_ESS
## sd(Intercept) 3867
## sd(Prime_type) 2594
## sd(OverlapOverlap) 3339
## sd(Agec) 2086
## sd(Prime_type:OverlapOverlap) 3535
## sd(Prime_type:Agec) 3709
## sd(OverlapOverlap:Agec) 3953
## sd(Prime_type:OverlapOverlap:Agec) 3146
## cor(Intercept,Prime_type) 5807
## cor(Intercept,OverlapOverlap) 6203
## cor(Prime_type,OverlapOverlap) 6164
## cor(Intercept,Agec) 6127
## cor(Prime_type,Agec) 6198
## cor(OverlapOverlap,Agec) 5743
## cor(Intercept,Prime_type:OverlapOverlap) 5661
## cor(Prime_type,Prime_type:OverlapOverlap) 6144
## cor(OverlapOverlap,Prime_type:OverlapOverlap) 6223
## cor(Agec,Prime_type:OverlapOverlap) 6179
## cor(Intercept,Prime_type:Agec) 5895
## cor(Prime_type,Prime_type:Agec) 5948
## cor(OverlapOverlap,Prime_type:Agec) 6309
## cor(Agec,Prime_type:Agec) 5992
## cor(Prime_type:OverlapOverlap,Prime_type:Agec) 6664
## cor(Intercept,OverlapOverlap:Agec) 5650
## cor(Prime_type,OverlapOverlap:Agec) 5786
## cor(OverlapOverlap,OverlapOverlap:Agec) 6026
## cor(Agec,OverlapOverlap:Agec) 6207
## cor(Prime_type:OverlapOverlap,OverlapOverlap:Agec) 6572
## cor(Prime_type:Agec,OverlapOverlap:Agec) 6841
## cor(Intercept,Prime_type:OverlapOverlap:Agec) 6101
## cor(Prime_type,Prime_type:OverlapOverlap:Agec) 6369
## cor(OverlapOverlap,Prime_type:OverlapOverlap:Agec) 6500
## cor(Agec,Prime_type:OverlapOverlap:Agec) 6866
## cor(Prime_type:OverlapOverlap,Prime_type:OverlapOverlap:Agec) 6229
## cor(Prime_type:Agec,Prime_type:OverlapOverlap:Agec) 6256
## cor(OverlapOverlap:Agec,Prime_type:OverlapOverlap:Agec) 7199
##
## Population-Level Effects:
## Estimate Est.Error l-90% CI u-90% CI Rhat
## Intercept -1.03 0.62 -1.99 -0.03 1.00
## Prime_type 1.03 0.40 0.42 1.67 1.00
## OverlapOverlap -0.68 0.34 -1.22 -0.14 1.00
## Agec 0.35 0.22 0.01 0.70 1.00
## Prime_type:OverlapOverlap 1.72 0.60 0.82 2.72 1.00
## Prime_type:Agec -0.20 0.22 -0.56 0.15 1.00
## OverlapOverlap:Agec 0.08 0.34 -0.45 0.62 1.00
## Prime_type:OverlapOverlap:Agec -0.01 0.38 -0.62 0.59 1.00
## Bulk_ESS Tail_ESS
## Intercept 1854 2210
## Prime_type 3439 3584
## OverlapOverlap 2672 3592
## Agec 2961 3316
## Prime_type:OverlapOverlap 3580 4321
## Prime_type:Agec 6127 5623
## OverlapOverlap:Agec 3740 4133
## Prime_type:OverlapOverlap:Agec 5008 4711
##
## 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).Check model with simpler random effects.
mXsec2bb <- brm(Response ~ 1 + Prime_type*Overlap*Agec + (1 + Prime_type*Overlap| Case_ID) + (1 + Prime_type*Overlap | Trial_verb), logit_prior, family="bernoulli", iter=4000, cores=4, control=list(adapt_delta=.95), data=XSectional2, file="mXsec2bb")summary(mXsec2bb, prob = .9)## Family: bernoulli
## Links: mu = logit
## Formula: Response ~ 1 + Prime_type * Overlap * Agec + (1 + Prime_type * Overlap | Case_ID) + (1 + Prime_type * Overlap | 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.25 0.16 0.99
## sd(Prime_type) 0.57 0.34 0.07
## sd(OverlapOverlap) 0.42 0.27 0.04
## sd(Prime_type:OverlapOverlap) 0.97 0.56 0.11
## cor(Intercept,Prime_type) -0.35 0.35 -0.83
## cor(Intercept,OverlapOverlap) 0.31 0.38 -0.40
## cor(Prime_type,OverlapOverlap) -0.16 0.42 -0.79
## cor(Intercept,Prime_type:OverlapOverlap) -0.24 0.35 -0.75
## cor(Prime_type,Prime_type:OverlapOverlap) 0.06 0.42 -0.64
## cor(OverlapOverlap,Prime_type:OverlapOverlap) -0.13 0.44 -0.78
## u-90% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 1.52 1.00 3030 5250
## sd(Prime_type) 1.17 1.00 1758 3085
## sd(OverlapOverlap) 0.89 1.00 1428 2652
## sd(Prime_type:OverlapOverlap) 1.93 1.00 1257 2360
## cor(Intercept,Prime_type) 0.33 1.00 4958 5444
## cor(Intercept,OverlapOverlap) 0.84 1.00 4658 4682
## cor(Prime_type,OverlapOverlap) 0.58 1.00 3373 5037
## cor(Intercept,Prime_type:OverlapOverlap) 0.42 1.00 5219 4303
## cor(Prime_type,Prime_type:OverlapOverlap) 0.75 1.00 2802 4271
## cor(OverlapOverlap,Prime_type:OverlapOverlap) 0.64 1.00 2478 4335
##
## ~Trial_verb (Number of levels: 6)
## Estimate Est.Error l-90% CI
## sd(Intercept) 1.38 0.48 0.81
## sd(Prime_type) 0.60 0.43 0.08
## sd(OverlapOverlap) 0.51 0.33 0.08
## sd(Prime_type:OverlapOverlap) 0.71 0.55 0.06
## cor(Intercept,Prime_type) -0.29 0.39 -0.84
## cor(Intercept,OverlapOverlap) -0.39 0.37 -0.89
## cor(Prime_type,OverlapOverlap) 0.03 0.42 -0.66
## cor(Intercept,Prime_type:OverlapOverlap) -0.15 0.41 -0.77
## cor(Prime_type,Prime_type:OverlapOverlap) 0.02 0.44 -0.70
## cor(OverlapOverlap,Prime_type:OverlapOverlap) -0.09 0.43 -0.76
## u-90% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 2.30 1.00 2979 5188
## sd(Prime_type) 1.41 1.00 2659 3113
## sd(OverlapOverlap) 1.12 1.00 2768 2902
## sd(Prime_type:OverlapOverlap) 1.75 1.00 3061 3720
## cor(Intercept,Prime_type) 0.43 1.00 6608 5449
## cor(Intercept,OverlapOverlap) 0.30 1.00 5894 4843
## cor(Prime_type,OverlapOverlap) 0.72 1.00 6721 5885
## cor(Intercept,Prime_type:OverlapOverlap) 0.58 1.00 9100 5927
## cor(Prime_type,Prime_type:OverlapOverlap) 0.74 1.00 7204 6657
## cor(OverlapOverlap,Prime_type:OverlapOverlap) 0.65 1.00 5950 6608
##
## Population-Level Effects:
## Estimate Est.Error l-90% CI u-90% CI Rhat
## Intercept -1.03 0.60 -1.99 -0.02 1.00
## Prime_type 1.03 0.39 0.44 1.67 1.00
## OverlapOverlap -0.59 0.32 -1.11 -0.08 1.00
## Agec 0.35 0.13 0.15 0.57 1.00
## Prime_type:OverlapOverlap 1.56 0.55 0.72 2.48 1.00
## Prime_type:Agec -0.23 0.18 -0.52 0.06 1.00
## OverlapOverlap:Agec -0.01 0.13 -0.23 0.21 1.00
## Prime_type:OverlapOverlap:Agec 0.13 0.27 -0.31 0.57 1.00
## Bulk_ESS Tail_ESS
## Intercept 1899 3133
## Prime_type 3730 3774
## OverlapOverlap 3287 4750
## Agec 3542 5393
## Prime_type:OverlapOverlap 3819 4503
## Prime_type:Agec 7014 6104
## OverlapOverlap:Agec 8450 6184
## Prime_type:OverlapOverlap:Agec 6531 5776
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).hypothesis(mXsec2bb, "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.03 0.39 0.44 1.67 99 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(mXsec2bb, "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.56 0.55 0.72 2.48 306.69
## 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(mXsec2bb, "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.13 0.15 0.57 274.86 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(mXsec2bb, "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.23 0.18 -0.52 0.06 9.51
## 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(mXsec2bb, "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.13 0.27 -0.31 0.57 2.23
## Post.Prob Star
## 1 0.69
## ---
## '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.#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")
mXsec3 <- 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)# Some larger Pareto's k vaues .7. I've re-fit the model without those observations. Nothing important changed, but here is the code for estimating this model.
mXsec3_no_low_k <- XSectional3 %>%
mutate(
k = mXsec3$criteria$loo$diagnostics$pareto_k
) %>%
filter(k < .7) %>%
brm(Response ~1 + Prime_type*biasmatch*Agec+ (1 +Prime_type*biasmatch | Case_ID) + (1 +Prime_type*biasmatch*Agec | Trial_verb), prior=logit_prior, family="bernoulli", iter=4000, cores=4, control=list(adapt_delta=.95), seed=12345, data=., file="mXsec3_no_low_k")
summary(mXsec3_no_low_k, 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: . (Number of observations: 824)
## 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.87 0.28 1.44
## sd(Prime_type) 0.59 0.45 0.05
## sd(biasmatchmismatch) 0.52 0.38 0.05
## sd(Prime_type:biasmatchmismatch) 1.06 0.73 0.10
## cor(Intercept,Prime_type) -0.10 0.42 -0.76
## cor(Intercept,biasmatchmismatch) 0.10 0.43 -0.64
## cor(Prime_type,biasmatchmismatch) -0.06 0.45 -0.77
## cor(Intercept,Prime_type:biasmatchmismatch) 0.08 0.40 -0.62
## cor(Prime_type,Prime_type:biasmatchmismatch) -0.18 0.45 -0.84
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) -0.01 0.44 -0.72
## u-90% CI Rhat Bulk_ESS
## sd(Intercept) 2.35 1.00 3083
## sd(Prime_type) 1.45 1.00 2303
## sd(biasmatchmismatch) 1.23 1.00 1273
## sd(Prime_type:biasmatchmismatch) 2.39 1.00 1276
## cor(Intercept,Prime_type) 0.64 1.00 8665
## cor(Intercept,biasmatchmismatch) 0.77 1.00 6362
## cor(Prime_type,biasmatchmismatch) 0.68 1.00 3097
## cor(Intercept,Prime_type:biasmatchmismatch) 0.71 1.00 7792
## cor(Prime_type,Prime_type:biasmatchmismatch) 0.62 1.00 3078
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) 0.72 1.00 3832
## Tail_ESS
## sd(Intercept) 5218
## sd(Prime_type) 3522
## sd(biasmatchmismatch) 2716
## sd(Prime_type:biasmatchmismatch) 3002
## cor(Intercept,Prime_type) 5916
## cor(Intercept,biasmatchmismatch) 5464
## cor(Prime_type,biasmatchmismatch) 4622
## cor(Intercept,Prime_type:biasmatchmismatch) 5249
## cor(Prime_type,Prime_type:biasmatchmismatch) 4627
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) 5221
##
## ~Trial_verb (Number of levels: 6)
## Estimate
## sd(Intercept) 1.37
## sd(Prime_type) 0.61
## sd(biasmatchmismatch) 0.71
## sd(Agec) 0.43
## sd(Prime_type:biasmatchmismatch) 1.02
## sd(Prime_type:Agec) 0.41
## sd(biasmatchmismatch:Agec) 0.36
## sd(Prime_type:biasmatchmismatch:Agec) 0.67
## cor(Intercept,Prime_type) -0.13
## cor(Intercept,biasmatchmismatch) 0.09
## cor(Prime_type,biasmatchmismatch) -0.06
## cor(Intercept,Agec) 0.21
## cor(Prime_type,Agec) -0.05
## cor(biasmatchmismatch,Agec) 0.06
## cor(Intercept,Prime_type:biasmatchmismatch) -0.14
## cor(Prime_type,Prime_type:biasmatchmismatch) -0.01
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) -0.08
## cor(Agec,Prime_type:biasmatchmismatch) -0.04
## cor(Intercept,Prime_type:Agec) -0.04
## cor(Prime_type,Prime_type:Agec) 0.00
## cor(biasmatchmismatch,Prime_type:Agec) -0.01
## 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.05
## cor(Prime_type:biasmatchmismatch,biasmatchmismatch:Agec) 0.01
## cor(Prime_type:Agec,biasmatchmismatch:Agec) -0.02
## cor(Intercept,Prime_type:biasmatchmismatch:Agec) 0.13
## cor(Prime_type,Prime_type:biasmatchmismatch:Agec) -0.04
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 0.06
## 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.49
## sd(biasmatchmismatch) 0.53
## sd(Agec) 0.29
## sd(Prime_type:biasmatchmismatch) 0.73
## sd(Prime_type:Agec) 0.36
## sd(biasmatchmismatch:Agec) 0.32
## sd(Prime_type:biasmatchmismatch:Agec) 0.55
## 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.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.32
## cor(Prime_type,Prime_type:Agec) 0.33
## 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.33
## cor(Prime_type:Agec,biasmatchmismatch:Agec) 0.34
## cor(Intercept,Prime_type:biasmatchmismatch:Agec) 0.33
## 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.33
## l-90% CI
## sd(Intercept) 0.69
## sd(Prime_type) 0.05
## sd(biasmatchmismatch) 0.06
## sd(Agec) 0.08
## sd(Prime_type:biasmatchmismatch) 0.10
## sd(Prime_type:Agec) 0.03
## sd(biasmatchmismatch:Agec) 0.03
## sd(Prime_type:biasmatchmismatch:Agec) 0.05
## cor(Intercept,Prime_type) -0.65
## cor(Intercept,biasmatchmismatch) -0.45
## cor(Prime_type,biasmatchmismatch) -0.59
## cor(Intercept,Agec) -0.33
## cor(Prime_type,Agec) -0.60
## cor(biasmatchmismatch,Agec) -0.50
## cor(Intercept,Prime_type:biasmatchmismatch) -0.65
## cor(Prime_type,Prime_type:biasmatchmismatch) -0.55
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) -0.61
## cor(Agec,Prime_type:biasmatchmismatch) -0.57
## cor(Intercept,Prime_type:Agec) -0.57
## cor(Prime_type,Prime_type:Agec) -0.54
## cor(biasmatchmismatch,Prime_type:Agec) -0.55
## cor(Agec,Prime_type:Agec) -0.58
## cor(Prime_type:biasmatchmismatch,Prime_type:Agec) -0.52
## cor(Intercept,biasmatchmismatch:Agec) -0.56
## cor(Prime_type,biasmatchmismatch:Agec) -0.54
## cor(biasmatchmismatch,biasmatchmismatch:Agec) -0.57
## cor(Agec,biasmatchmismatch:Agec) -0.59
## cor(Prime_type:biasmatchmismatch,biasmatchmismatch:Agec) -0.54
## cor(Prime_type:Agec,biasmatchmismatch:Agec) -0.57
## cor(Intercept,Prime_type:biasmatchmismatch:Agec) -0.45
## cor(Prime_type,Prime_type:biasmatchmismatch:Agec) -0.59
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) -0.50
## cor(Agec,Prime_type:biasmatchmismatch:Agec) -0.47
## cor(Prime_type:biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) -0.58
## cor(Prime_type:Agec,Prime_type:biasmatchmismatch:Agec) -0.63
## cor(biasmatchmismatch:Agec,Prime_type:biasmatchmismatch:Agec) -0.56
## u-90% CI
## sd(Intercept) 2.37
## sd(Prime_type) 1.54
## sd(biasmatchmismatch) 1.73
## sd(Agec) 0.96
## sd(Prime_type:biasmatchmismatch) 2.38
## sd(Prime_type:Agec) 1.12
## sd(biasmatchmismatch:Agec) 0.96
## sd(Prime_type:biasmatchmismatch:Agec) 1.74
## cor(Intercept,Prime_type) 0.44
## cor(Intercept,biasmatchmismatch) 0.60
## cor(Prime_type,biasmatchmismatch) 0.49
## cor(Intercept,Agec) 0.69
## cor(Prime_type,Agec) 0.51
## cor(biasmatchmismatch,Agec) 0.58
## cor(Intercept,Prime_type:biasmatchmismatch) 0.42
## cor(Prime_type,Prime_type:biasmatchmismatch) 0.54
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) 0.48
## cor(Agec,Prime_type:biasmatchmismatch) 0.51
## cor(Intercept,Prime_type:Agec) 0.51
## cor(Prime_type,Prime_type:Agec) 0.55
## cor(biasmatchmismatch,Prime_type:Agec) 0.54
## cor(Agec,Prime_type:Agec) 0.53
## cor(Prime_type:biasmatchmismatch,Prime_type:Agec) 0.58
## cor(Intercept,biasmatchmismatch:Agec) 0.53
## cor(Prime_type,biasmatchmismatch:Agec) 0.58
## cor(biasmatchmismatch,biasmatchmismatch:Agec) 0.54
## cor(Agec,biasmatchmismatch:Agec) 0.51
## cor(Prime_type:biasmatchmismatch,biasmatchmismatch:Agec) 0.56
## cor(Prime_type:Agec,biasmatchmismatch:Agec) 0.54
## cor(Intercept,Prime_type:biasmatchmismatch:Agec) 0.64
## cor(Prime_type,Prime_type:biasmatchmismatch:Agec) 0.52
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 0.60
## cor(Agec,Prime_type:biasmatchmismatch:Agec) 0.61
## cor(Prime_type:biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 0.50
## cor(Prime_type:Agec,Prime_type:biasmatchmismatch:Agec) 0.48
## cor(biasmatchmismatch:Agec,Prime_type:biasmatchmismatch:Agec) 0.54
## 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) 3892
## sd(Prime_type) 3837
## sd(biasmatchmismatch) 3350
## sd(Agec) 2906
## sd(Prime_type:biasmatchmismatch) 3990
## sd(Prime_type:Agec) 4452
## sd(biasmatchmismatch:Agec) 4272
## sd(Prime_type:biasmatchmismatch:Agec) 3509
## cor(Intercept,Prime_type) 9665
## cor(Intercept,biasmatchmismatch) 10778
## cor(Prime_type,biasmatchmismatch) 7010
## cor(Intercept,Agec) 7461
## cor(Prime_type,Agec) 7234
## cor(biasmatchmismatch,Agec) 6216
## cor(Intercept,Prime_type:biasmatchmismatch) 11609
## cor(Prime_type,Prime_type:biasmatchmismatch) 8803
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) 7076
## cor(Agec,Prime_type:biasmatchmismatch) 7852
## cor(Intercept,Prime_type:Agec) 11336
## cor(Prime_type,Prime_type:Agec) 10349
## cor(biasmatchmismatch,Prime_type:Agec) 8094
## cor(Agec,Prime_type:Agec) 7997
## cor(Prime_type:biasmatchmismatch,Prime_type:Agec) 6668
## cor(Intercept,biasmatchmismatch:Agec) 11230
## cor(Prime_type,biasmatchmismatch:Agec) 10076
## cor(biasmatchmismatch,biasmatchmismatch:Agec) 8071
## cor(Agec,biasmatchmismatch:Agec) 7843
## cor(Prime_type:biasmatchmismatch,biasmatchmismatch:Agec) 6289
## cor(Prime_type:Agec,biasmatchmismatch:Agec) 5832
## cor(Intercept,Prime_type:biasmatchmismatch:Agec) 10513
## cor(Prime_type,Prime_type:biasmatchmismatch:Agec) 9068
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 8179
## cor(Agec,Prime_type:biasmatchmismatch:Agec) 7577
## cor(Prime_type:biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 6427
## cor(Prime_type:Agec,Prime_type:biasmatchmismatch:Agec) 4650
## cor(biasmatchmismatch:Agec,Prime_type:biasmatchmismatch:Agec) 5189
## Tail_ESS
## sd(Intercept) 5695
## sd(Prime_type) 4561
## sd(biasmatchmismatch) 3220
## sd(Agec) 3034
## sd(Prime_type:biasmatchmismatch) 3977
## sd(Prime_type:Agec) 4393
## sd(biasmatchmismatch:Agec) 4795
## sd(Prime_type:biasmatchmismatch:Agec) 3940
## cor(Intercept,Prime_type) 6153
## cor(Intercept,biasmatchmismatch) 6263
## cor(Prime_type,biasmatchmismatch) 6356
## cor(Intercept,Agec) 6486
## cor(Prime_type,Agec) 6120
## cor(biasmatchmismatch,Agec) 6110
## cor(Intercept,Prime_type:biasmatchmismatch) 6606
## cor(Prime_type,Prime_type:biasmatchmismatch) 6394
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) 6220
## cor(Agec,Prime_type:biasmatchmismatch) 6373
## cor(Intercept,Prime_type:Agec) 5847
## cor(Prime_type,Prime_type:Agec) 6123
## cor(biasmatchmismatch,Prime_type:Agec) 5967
## cor(Agec,Prime_type:Agec) 6460
## cor(Prime_type:biasmatchmismatch,Prime_type:Agec) 6771
## cor(Intercept,biasmatchmismatch:Agec) 6074
## cor(Prime_type,biasmatchmismatch:Agec) 6449
## cor(biasmatchmismatch,biasmatchmismatch:Agec) 6622
## cor(Agec,biasmatchmismatch:Agec) 6245
## cor(Prime_type:biasmatchmismatch,biasmatchmismatch:Agec) 6479
## cor(Prime_type:Agec,biasmatchmismatch:Agec) 6610
## cor(Intercept,Prime_type:biasmatchmismatch:Agec) 5267
## cor(Prime_type,Prime_type:biasmatchmismatch:Agec) 6218
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 6276
## cor(Agec,Prime_type:biasmatchmismatch:Agec) 6443
## cor(Prime_type:biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 6693
## cor(Prime_type:Agec,Prime_type:biasmatchmismatch:Agec) 6029
## cor(biasmatchmismatch:Agec,Prime_type:biasmatchmismatch:Agec) 6382
##
## Population-Level Effects:
## Estimate Est.Error l-90% CI u-90% CI Rhat
## Intercept -2.09 0.65 -3.13 -1.04 1.00
## Prime_type 0.61 0.55 -0.25 1.51 1.00
## biasmatchmismatch 0.54 0.49 -0.26 1.28 1.00
## Agec 0.63 0.30 0.18 1.10 1.00
## Prime_type:biasmatchmismatch 0.83 0.85 -0.47 2.30 1.00
## Prime_type:Agec -0.09 0.39 -0.71 0.54 1.00
## biasmatchmismatch:Agec -0.22 0.30 -0.68 0.25 1.00
## Prime_type:biasmatchmismatch:Agec -0.31 0.58 -1.23 0.61 1.00
## Bulk_ESS Tail_ESS
## Intercept 2573 3691
## Prime_type 4784 4712
## biasmatchmismatch 4097 4205
## Agec 3946 4537
## Prime_type:biasmatchmismatch 4774 4842
## Prime_type:Agec 5425 4924
## biasmatchmismatch:Agec 5285 4537
## Prime_type:biasmatchmismatch:Agec 4558 4136
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).hypothesis(mXsec3_no_low_k, "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.61 0.55 -0.25 1.51 7.3 0.88
## 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_no_low_k, "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.83 0.85 -0.47 2.3 5.84
## Post.Prob Star
## 1 0.85
## ---
## '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_no_low_k, "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.09 0.39 -0.71 0.54 0.67
## Post.Prob Star
## 1 0.4
## ---
## '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_no_low_k, "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.31 0.58 -1.23 0.61 2.57
## Post.Prob Star
## 1 0.72
## ---
## '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)
)
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"))## Warning: Ignoring unknown aesthetics: fill## Warning: The labeller API has been updated. Labellers taking `variable` and
## `value` arguments are now deprecated. See labellers documentation.
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.5) +
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()
) ## Warning: Ignoring unknown aesthetics: fillposterior_plot3/m3_fitted_plot + plot_layout(guides = "collect") & theme(legend.position = "right") 
ggsave("plots/Fig4.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.XSectional4 %>%
mutate(
k = mXsec4$criteria$loo$diagnostics$pareto_k
) %>%
ggplot(aes(x=k)) + geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
XSectional4 %>%
mutate(
k = mXsec4$criteria$loo$diagnostics$pareto_k
) %>%
filter(k < .7) %>%
brm(Response ~1 + Prime_type*biasmatch*Agec+ (1 +Prime_type*biasmatch | Case_ID) + (1 +Prime_type*biasmatch*Agec | Trial_verb), prior=logit_prior, family="bernoulli", iter=4000, cores=4, seed=12345, data=., file="mXsec4_no_low_k") %>%
summary(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: . (Number of observations: 673)
## Draws: 4 chains, each with iter = 8000; warmup = 4000; thin = 1;
## total post-warmup draws = 16000
##
## Group-Level Effects:
## ~Case_ID (Number of levels: 114)
## Estimate Est.Error l-90% CI
## sd(Intercept) 1.35 0.24 0.97
## sd(Prime_type) 0.58 0.42 0.05
## sd(biasmatchmismatch) 0.54 0.39 0.04
## sd(Prime_type:biasmatchmismatch) 1.06 0.74 0.09
## cor(Intercept,Prime_type) -0.08 0.42 -0.75
## cor(Intercept,biasmatchmismatch) 0.04 0.41 -0.64
## cor(Prime_type,biasmatchmismatch) -0.01 0.44 -0.73
## cor(Intercept,Prime_type:biasmatchmismatch) -0.05 0.40 -0.69
## cor(Prime_type,Prime_type:biasmatchmismatch) -0.14 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 6828
## sd(Prime_type) 1.39 1.00 5867
## sd(biasmatchmismatch) 1.27 1.00 2755
## sd(Prime_type:biasmatchmismatch) 2.46 1.00 3826
## cor(Intercept,Prime_type) 0.65 1.00 18148
## cor(Intercept,biasmatchmismatch) 0.72 1.00 16189
## cor(Prime_type,biasmatchmismatch) 0.72 1.00 7462
## cor(Intercept,Prime_type:biasmatchmismatch) 0.62 1.00 15903
## cor(Prime_type,Prime_type:biasmatchmismatch) 0.65 1.00 8496
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) 0.72 1.00 9141
## Tail_ESS
## sd(Intercept) 10184
## sd(Prime_type) 8873
## sd(biasmatchmismatch) 5670
## sd(Prime_type:biasmatchmismatch) 6358
## cor(Intercept,Prime_type) 10868
## cor(Intercept,biasmatchmismatch) 10835
## cor(Prime_type,biasmatchmismatch) 11383
## cor(Intercept,Prime_type:biasmatchmismatch) 12064
## cor(Prime_type,Prime_type:biasmatchmismatch) 11176
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) 12332
##
## ~Trial_verb (Number of levels: 6)
## Estimate
## sd(Intercept) 1.32
## sd(Prime_type) 0.56
## sd(biasmatchmismatch) 0.67
## sd(Agec) 0.39
## sd(Prime_type:biasmatchmismatch) 1.16
## sd(Prime_type:Agec) 0.48
## sd(biasmatchmismatch:Agec) 0.36
## sd(Prime_type:biasmatchmismatch:Agec) 0.92
## 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.02
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) -0.09
## cor(Agec,Prime_type:biasmatchmismatch) -0.06
## cor(Intercept,Prime_type:Agec) -0.06
## cor(Prime_type,Prime_type:Agec) -0.00
## cor(biasmatchmismatch,Prime_type:Agec) -0.03
## cor(Agec,Prime_type:Agec) -0.06
## cor(Prime_type:biasmatchmismatch,Prime_type:Agec) 0.06
## cor(Intercept,biasmatchmismatch:Agec) 0.00
## 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.18
## cor(Prime_type,Prime_type:biasmatchmismatch:Agec) -0.04
## 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.13
## cor(biasmatchmismatch:Agec,Prime_type:biasmatchmismatch:Agec) 0.01
## Est.Error
## sd(Intercept) 0.52
## sd(Prime_type) 0.48
## sd(biasmatchmismatch) 0.52
## sd(Agec) 0.28
## sd(Prime_type:biasmatchmismatch) 0.76
## 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.33
## cor(Prime_type,Prime_type:Agec) 0.33
## cor(biasmatchmismatch,Prime_type:Agec) 0.33
## cor(Agec,Prime_type:Agec) 0.33
## 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.33
## cor(Prime_type:biasmatchmismatch,biasmatchmismatch:Agec) 0.33
## cor(Prime_type:Agec,biasmatchmismatch:Agec) 0.33
## 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.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.13
## sd(Prime_type:Agec) 0.04
## sd(biasmatchmismatch:Agec) 0.03
## sd(Prime_type:biasmatchmismatch:Agec) 0.09
## cor(Intercept,Prime_type) -0.62
## cor(Intercept,biasmatchmismatch) -0.45
## cor(Prime_type,biasmatchmismatch) -0.58
## cor(Intercept,Agec) -0.34
## cor(Prime_type,Agec) -0.58
## cor(biasmatchmismatch,Agec) -0.48
## cor(Intercept,Prime_type:biasmatchmismatch) -0.65
## cor(Prime_type,Prime_type:biasmatchmismatch) -0.56
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) -0.62
## cor(Agec,Prime_type:biasmatchmismatch) -0.59
## cor(Intercept,Prime_type:Agec) -0.59
## cor(Prime_type,Prime_type:Agec) -0.55
## cor(biasmatchmismatch,Prime_type:Agec) -0.56
## 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.53
## cor(biasmatchmismatch,biasmatchmismatch:Agec) -0.56
## cor(Agec,biasmatchmismatch:Agec) -0.58
## cor(Prime_type:biasmatchmismatch,biasmatchmismatch:Agec) -0.54
## cor(Prime_type:Agec,biasmatchmismatch:Agec) -0.58
## cor(Intercept,Prime_type:biasmatchmismatch:Agec) -0.38
## cor(Prime_type,Prime_type:biasmatchmismatch:Agec) -0.59
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) -0.46
## cor(Agec,Prime_type:biasmatchmismatch:Agec) -0.46
## cor(Prime_type:biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) -0.63
## 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.31
## sd(Prime_type) 1.47
## sd(biasmatchmismatch) 1.65
## sd(Agec) 0.91
## sd(Prime_type:biasmatchmismatch) 2.57
## sd(Prime_type:Agec) 1.27
## sd(biasmatchmismatch:Agec) 0.97
## sd(Prime_type:biasmatchmismatch:Agec) 2.17
## cor(Intercept,Prime_type) 0.46
## cor(Intercept,biasmatchmismatch) 0.59
## cor(Prime_type,biasmatchmismatch) 0.52
## cor(Intercept,Agec) 0.69
## cor(Prime_type,Agec) 0.51
## cor(biasmatchmismatch,Agec) 0.60
## cor(Intercept,Prime_type:biasmatchmismatch) 0.38
## cor(Prime_type,Prime_type:biasmatchmismatch) 0.53
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) 0.47
## cor(Agec,Prime_type:biasmatchmismatch) 0.48
## cor(Intercept,Prime_type:Agec) 0.49
## cor(Prime_type,Prime_type:Agec) 0.55
## cor(biasmatchmismatch,Prime_type:Agec) 0.52
## cor(Agec,Prime_type:Agec) 0.49
## cor(Prime_type:biasmatchmismatch,Prime_type:Agec) 0.60
## cor(Intercept,biasmatchmismatch:Agec) 0.55
## 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.52
## cor(Intercept,Prime_type:biasmatchmismatch:Agec) 0.67
## cor(Prime_type,Prime_type:biasmatchmismatch:Agec) 0.52
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 0.63
## cor(Agec,Prime_type:biasmatchmismatch:Agec) 0.63
## cor(Prime_type:biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 0.46
## cor(Prime_type:Agec,Prime_type:biasmatchmismatch:Agec) 0.45
## 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) 10591
## sd(Prime_type) 9092
## sd(biasmatchmismatch) 8230
## sd(Agec) 7201
## sd(Prime_type:biasmatchmismatch) 8226
## sd(Prime_type:Agec) 7444
## sd(biasmatchmismatch:Agec) 8914
## sd(Prime_type:biasmatchmismatch:Agec) 7868
## cor(Intercept,Prime_type) 22788
## cor(Intercept,biasmatchmismatch) 21885
## cor(Prime_type,biasmatchmismatch) 17635
## cor(Intercept,Agec) 17993
## cor(Prime_type,Agec) 16841
## cor(biasmatchmismatch,Agec) 14233
## cor(Intercept,Prime_type:biasmatchmismatch) 22615
## cor(Prime_type,Prime_type:biasmatchmismatch) 17299
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) 15660
## cor(Agec,Prime_type:biasmatchmismatch) 14417
## cor(Intercept,Prime_type:Agec) 22300
## cor(Prime_type,Prime_type:Agec) 18797
## cor(biasmatchmismatch,Prime_type:Agec) 16551
## cor(Agec,Prime_type:Agec) 14918
## cor(Prime_type:biasmatchmismatch,Prime_type:Agec) 11864
## cor(Intercept,biasmatchmismatch:Agec) 23123
## cor(Prime_type,biasmatchmismatch:Agec) 18913
## cor(biasmatchmismatch,biasmatchmismatch:Agec) 15872
## cor(Agec,biasmatchmismatch:Agec) 14724
## cor(Prime_type:biasmatchmismatch,biasmatchmismatch:Agec) 12773
## cor(Prime_type:Agec,biasmatchmismatch:Agec) 11711
## cor(Intercept,Prime_type:biasmatchmismatch:Agec) 19767
## cor(Prime_type,Prime_type:biasmatchmismatch:Agec) 18010
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 15006
## cor(Agec,Prime_type:biasmatchmismatch:Agec) 15299
## cor(Prime_type:biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 12982
## cor(Prime_type:Agec,Prime_type:biasmatchmismatch:Agec) 11377
## cor(biasmatchmismatch:Agec,Prime_type:biasmatchmismatch:Agec) 9914
## Tail_ESS
## sd(Intercept) 12027
## sd(Prime_type) 10058
## sd(biasmatchmismatch) 8920
## sd(Agec) 7811
## sd(Prime_type:biasmatchmismatch) 7655
## sd(Prime_type:Agec) 10148
## sd(biasmatchmismatch:Agec) 9548
## sd(Prime_type:biasmatchmismatch:Agec) 9414
## cor(Intercept,Prime_type) 12961
## cor(Intercept,biasmatchmismatch) 12533
## cor(Prime_type,biasmatchmismatch) 12304
## cor(Intercept,Agec) 12222
## cor(Prime_type,Agec) 13321
## cor(biasmatchmismatch,Agec) 12636
## cor(Intercept,Prime_type:biasmatchmismatch) 12539
## cor(Prime_type,Prime_type:biasmatchmismatch) 12327
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch) 13375
## cor(Agec,Prime_type:biasmatchmismatch) 13220
## cor(Intercept,Prime_type:Agec) 12191
## cor(Prime_type,Prime_type:Agec) 13029
## cor(biasmatchmismatch,Prime_type:Agec) 12659
## cor(Agec,Prime_type:Agec) 13072
## cor(Prime_type:biasmatchmismatch,Prime_type:Agec) 13325
## cor(Intercept,biasmatchmismatch:Agec) 11681
## cor(Prime_type,biasmatchmismatch:Agec) 13151
## cor(biasmatchmismatch,biasmatchmismatch:Agec) 12710
## cor(Agec,biasmatchmismatch:Agec) 13419
## cor(Prime_type:biasmatchmismatch,biasmatchmismatch:Agec) 13761
## cor(Prime_type:Agec,biasmatchmismatch:Agec) 12976
## cor(Intercept,Prime_type:biasmatchmismatch:Agec) 12024
## cor(Prime_type,Prime_type:biasmatchmismatch:Agec) 13017
## cor(biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 12948
## cor(Agec,Prime_type:biasmatchmismatch:Agec) 13067
## cor(Prime_type:biasmatchmismatch,Prime_type:biasmatchmismatch:Agec) 13631
## cor(Prime_type:Agec,Prime_type:biasmatchmismatch:Agec) 13567
## cor(biasmatchmismatch:Agec,Prime_type:biasmatchmismatch:Agec) 13094
##
## Population-Level Effects:
## Estimate Est.Error l-90% CI u-90% CI Rhat
## Intercept -1.41 0.60 -2.37 -0.43 1.00
## Prime_type 0.48 0.52 -0.34 1.32 1.00
## biasmatchmismatch 0.53 0.46 -0.21 1.24 1.00
## Agec 0.44 0.27 0.02 0.86 1.00
## Prime_type:biasmatchmismatch 1.04 0.87 -0.30 2.54 1.00
## Prime_type:Agec -0.01 0.40 -0.65 0.63 1.00
## biasmatchmismatch:Agec -0.18 0.29 -0.64 0.30 1.00
## Prime_type:biasmatchmismatch:Agec -0.38 0.65 -1.44 0.65 1.00
## Bulk_ESS Tail_ESS
## Intercept 7478 8949
## Prime_type 11905 10043
## biasmatchmismatch 10742 8587
## Agec 10814 9486
## Prime_type:biasmatchmismatch 11794 10319
## Prime_type:Agec 12481 9797
## biasmatchmismatch:Agec 13273 10325
## Prime_type:biasmatchmismatch:Agec 11300 9325
##
## 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).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"))## Warning: 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.5) +
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()
) ## Warning: Ignoring unknown aesthetics: fillposterior_plot4/m4_fitted_plot + plot_layout(guides = "collect") & theme(legend.position = "right") 
ggsave("plots/Fig5.png", last_plot(), dpi=500)## Saving 7 x 5 in image3.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)## Warning: There were 1 divergent transitions after warmup. Increasing adapt_delta
## above 0.8 may help. See http://mc-stan.org/misc/warnings.html#divergent-
## transitions-after-warmup## 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.18 1.73 1.00 2528 4227
##
## ~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.45 0.55 2.28 1.00 2182 3655
## sd(Agec) 0.14 0.13 0.00 0.47 1.00 2212 2875
## cor(Intercept,Agec) 0.03 0.52 -0.89 0.93 1.00 5336 4190
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -1.53 0.49 -2.51 -0.55 1.00 1446 2530
## Agec 0.38 0.13 0.11 0.64 1.00 2248 3412
##
## 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, group=Trial_verb)) +
geom_hline(yintercept=.5, linetype="dotted") +
geom_ribbon(data=summ, aes(x=Age, y=median, ymin=Q5, ymax=Q95, group=Trial_verb, fill=Trial_verb), alpha=0.3) +
theme_minimal() +
geom_point(data=points, aes(x=Age, y=mean, fill=Trial_verb, color=Trial_verb), shape=3) +
ylab("Proportion") +
labs(color="Verb", fill="Verb") +
scale_x_continuous(breaks = round(seq(min(summ$Age), max(summ$Age), by = 1),1)) ## Warning: Ignoring unknown aesthetics: ymin, ymax, fill
ggsave("plots/Fig6.png", last_plot(), dpi=500)## Saving 7 x 5 in imagesessionInfo()## R version 4.1.0 (2021-05-18)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19044)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=Dutch_Netherlands.1252 LC_CTYPE=Dutch_Netherlands.1252
## [3] LC_MONETARY=Dutch_Netherlands.1252 LC_NUMERIC=C
## [5] LC_TIME=Dutch_Netherlands.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] sjPlot_2.8.10 patchwork_1.1.1 ggpubr_0.4.0 DHARMa_0.4.6
## [5] GGally_2.1.2 dplyr_1.0.7 tidyr_1.2.0 ggplot2_3.3.5
## [9] tidybayes_3.0.2 rmarkdown_2.13 readxl_1.3.1 brms_2.17.0
## [13] Rcpp_1.0.8.3
##
## loaded via a namespace (and not attached):
## [1] backports_1.4.1 systemfonts_1.0.4 plyr_1.8.6
## [4] igraph_1.2.11 splines_4.1.0 svUnit_1.0.6
## [7] crosstalk_1.2.0 gap.datasets_0.0.5 TH.data_1.1-0
## [10] rstantools_2.1.1 inline_0.3.19 digest_0.6.29
## [13] foreach_1.5.2 htmltools_0.5.2 fansi_1.0.2
## [16] magrittr_2.0.1 checkmate_2.0.0 doParallel_1.0.17
## [19] modelr_0.1.8 RcppParallel_5.1.5 matrixStats_0.61.0
## [22] xts_0.12.1 sandwich_3.0-1 prettyunits_1.1.1
## [25] colorspace_2.0-3 ggdist_3.1.1 textshaping_0.3.6
## [28] xfun_0.30 callr_3.7.0 crayon_1.5.0
## [31] jsonlite_1.8.0 lme4_1.1-28 iterators_1.0.14
## [34] survival_3.3-1 zoo_1.8-9 glue_1.6.2
## [37] gtable_0.3.1 emmeans_1.7.2 sjstats_0.18.1
## [40] sjmisc_2.8.9 distributional_0.3.0 car_3.0-12
## [43] pkgbuild_1.3.1 rstan_2.21.3 abind_1.4-5
## [46] scales_1.2.1 mvtnorm_1.1-3 DBI_1.1.2
## [49] ggeffects_1.1.1 rstatix_0.7.0 miniUI_0.1.1.1
## [52] performance_0.8.0 xtable_1.8-4 diffobj_0.3.5
## [55] stats4_4.1.0 StanHeaders_2.21.0-7 DT_0.21
## [58] datawizard_0.3.0 htmlwidgets_1.5.4 threejs_0.3.3
## [61] arrayhelpers_1.1-0 RColorBrewer_1.1-3 posterior_1.2.1
## [64] ellipsis_0.3.2 pkgconfig_2.0.3 reshape_0.8.8
## [67] loo_2.5.0 farver_2.1.0 qgam_1.3.4
## [70] sass_0.4.0 utf8_1.2.2 labeling_0.4.2
## [73] effectsize_0.6.0.1 tidyselect_1.1.2 rlang_1.0.6
## [76] reshape2_1.4.4 later_1.3.0 munsell_0.5.0
## [79] cellranger_1.1.0 tools_4.1.0 cli_3.3.0
## [82] generics_0.1.3 sjlabelled_1.1.8 broom_0.7.12
## [85] ggridges_0.5.3 evaluate_0.15 stringr_1.4.0
## [88] fastmap_1.1.0 ragg_1.2.2 yaml_2.3.5
## [91] processx_3.5.2 knitr_1.37 purrr_0.3.4
## [94] nlme_3.1-155 mime_0.12 projpred_2.0.2
## [97] gap_1.3-1 compiler_4.1.0 bayesplot_1.9.0
## [100] shinythemes_1.2.0 rstudioapi_0.13 gamm4_0.2-6
## [103] ggsignif_0.6.3 tibble_3.1.6 bslib_0.3.1
## [106] stringi_1.7.6 highr_0.9 parameters_0.17.0
## [109] ps_1.6.0 Brobdingnag_1.2-8 lattice_0.20-45
## [112] Matrix_1.4-0 nloptr_2.0.0 markdown_1.1
## [115] shinyjs_2.1.0 tensorA_0.36.2 vctrs_0.3.8
## [118] pillar_1.8.1 lifecycle_1.0.3 jquerylib_0.1.4
## [121] bridgesampling_1.1-2 estimability_1.3 insight_0.16.0
## [124] httpuv_1.6.5 R6_2.5.1 promises_1.2.0.1
## [127] gridExtra_2.3 codetools_0.2-18 boot_1.3-28
## [130] colourpicker_1.1.1 MASS_7.3-55 gtools_3.9.2
## [133] assertthat_0.2.1 withr_2.5.0 shinystan_2.6.0
## [136] multcomp_1.4-18 mgcv_1.8-39 bayestestR_0.11.5
## [139] parallel_4.1.0 grid_4.1.0 coda_0.19-4
## [142] minqa_1.2.4 carData_3.0-5 lubridate_1.8.0
## [145] shiny_1.7.1 base64enc_0.1-3 dygraphs_1.1.1.6