pilot1
Samah
2023-03-14
Preliminaries and data reading
# additional libraries
library("knitr")
library("janitor")
library("broom.mixed")
library("lme4")
library("emmeans")
library("tidyverse")
library("kableExtra")
library("modelr")
library("broom")
library("nlme")
library("meta")
library("metafor")
library(jtools) # Load jtools
theme_set(theme_bw())# reading the data file
pilot1_data = read_csv("252.csv")
df_shape= filter(pilot1_data, !is.na(d))
# pilot1_data = pilot1_data %>%
# select(ID, Title, d, d_var, Author)
df_shape_summary = df_shape %>%
group_by(ID, Title, Author) %>%
summarize(mean = mean(d),
mean_se = mean(d_var))Initial exploration
First, Visualizing the data to have an initial idea of how it looks.
First dividing the language group into two groups: the first one is the indo-european group which includes the English and the Spanish languages. The second group includes the rest of the languages: Japanese, Chinese, Tsimane.
df_shape$indoeuropean <- fct_relevel(as.factor(df_shape$language %in%
c("english","spanish", "german")),
"TRUE")
# creating a plot that shows the effects sizes colored per language group as well as the polynomial regression curve that fits it.
ggplot(df_shape,
aes(x = mean_age_months, y = d, color = language))+
geom_pointrange(aes(ymin = d - d_var, ymax = d + d_var),
alpha = .5) +
geom_smooth(aes(group = indoeuropean,
lty = indoeuropean),
col = "black",
method = "lm", se = FALSE,
formula = y ~ poly(x,2)) +
geom_hline(yintercept = 0, lty = 3) +
ylab("Standardized Mean Difference (d)") +
xlab("Mean age (months)") +
scale_color_discrete(name = "Language") +
scale_linetype_discrete(name = "Indo-European") +
theme(legend.position = "bottom")## Warning: Removed 8 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 8 rows containing missing values (`geom_pointrange()`).
# ggsave("first graph.png", width = 7, height = 4)Basic meta-analysis
metagen approach
using the meta-analytic function meta-gen which calculates the weights for each effects and confidence interval, pooled effect size, the heterogeneity.
m.gen <- metagen( TE= d,
seTE = d_var,
studlab = ID,
data = df_shape,
sm = "SMD",
fixed = FALSE,
random = TRUE,
method.tau = "REML",
hakn = TRUE,
title = "pilot shape bias meta-analysis"
)
# summary(m.gen)['TE']
# m.gen["TE.fixed"]
# m.gen["TE.random"]
# m.gen["w.random"]forest plot using the m-gen function object
forextobj <- forest.meta(m.gen,
sortvar = TE,
prediction = TRUE,
print.tau2 = FALSE,
leftlabs = c("Author", "g", "SE"))
forest plot from the rma model
metafor approach
using rma.mv instead of m.gen
mod <- rma.mv(yi = d,
V = d_var,
random = ~ 1 | ID,
slab = short_cite,
data = df_shape)
summary(mod)##
## Multivariate Meta-Analysis Model (k = 140; method: REML)
##
## logLik Deviance AIC BIC AICc
## -345.7749 691.5499 695.5499 701.4188 695.6381
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.1518 0.3897 25 no ID
##
## Test for Heterogeneity:
## Q(df = 139) = 959.8543, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.4669 0.0828 5.6377 <.0001 0.3046 0.6293 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Nested model.
mod_nested <- rma.mv(yi = d,
V = d_var,
random = ~ 1 | ID/exp_num,
slab = short_cite,
data = filter(df_shape, !is.na(exp_num)))
summary(mod_nested)##
## Multivariate Meta-Analysis Model (k = 133; method: REML)
##
## logLik Deviance AIC BIC AICc
## -328.8672 657.7344 663.7344 672.3828 663.9219
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.1446 0.3803 25 no ID
## sigma^2.2 0.0093 0.0964 34 no ID/exp_num
##
## Test for Heterogeneity:
## Q(df = 132) = 895.1881, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.4670 0.0838 5.5735 <.0001 0.3028 0.6312 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Forest plot from metafor.
forest(mod) 
Try this using ggplot
forest_data <-tibble(yi = mod$yi,
se = sqrt(mod$vi),
slab = mod$slab)
ggplot(forest_data,
aes(x = slab, y = yi)) +
geom_pointrange(aes(ymin = yi - se, ymax = yi + se)) +
geom_hline(yintercept = 0, lty = 2) +
coord_flip() 
#
# theme_set(theme_bw(base_size=10))
# ata.frame(ES=ROM.ma$b,SE=ROM.ma$se,Type="Summary",Study="Summary"))
# forrest_data$Study2<-factor(forrest_data$Study, levels=rev(levels(forrest_data$Study)) )
# levels(forrest_data$Study2)
# plot1<-ggplot(data=forrest_data,aes(x=Study2,y=ES,ymax=ES+(1.96*SE),ymin=ES-(1.96*SE),size=factor(Type),colour=factor(Type)))+geom_pointrange()
# plot2<-plot1+coord_flip()+geom_hline(aes(x=0), lty=2,size=1)+scale_size_manual(values=c(0.5,1))
# plot3<-plot2+xlab("Study")+ylab("log response ratio")+scale_colour_manual(values=c("grey","black"))
# plot3+theme(legend.position="none")Revised forest plot
df_shape$mean_age_months_centered36 <- df_shape$mean_age_months - 36
df_shape$log_mean_age_months <- log(df_shape$mean_age_months)
ggplot(df_shape, aes(x = short_cite, y = d,
ymin=d-sqrt(d_var)*1.96,
ymax=d+sqrt(d_var)*1.96)) +
geom_pointrange(aes(color=indoeuropean), alpha = .5, position=position_dodge2(width=.5)) +
coord_flip() +
geom_hline(yintercept = 0, lty = 2) +
aes(x=reorder(short_cite,-d, sum)) +
ylab("Standardized Mean Difference (d)") +
xlab("Citation")
#png("secondgraph.png")col.contour = c("gray75", "gray85", "gray95")
funnel(m.gen,
comb.random = TRUE,
xlim = c(-2, 4),
contour = c(0.9, 0.95, 0.99),
col.contour = col.contour)
regtest(x = d, vi = d_var,
data = df_shape)##
## Regression Test for Funnel Plot Asymmetry
##
## Model: mixed-effects meta-regression model
## Predictor: standard error
##
## Test for Funnel Plot Asymmetry: z = 6.6402, p < .0001
## Limit Estimate (as sei -> 0): b = -0.8597 (CI: -1.2814, -0.4381)# Add a legend
legend(x = 3.3, y = 0.1, cex = 0.5,
legend = c("p < 0.1", "p < 0.05", "p < 0.01"),
fill = col.contour)
#png("funnel.png")Publication bias
funnel plots using ggplot to account for moderators:
x = summary(m.gen)['TE']
y = summary(m.gen)['seTE']
m.gen["TE.fixed"]## $TE.fixed
## [1] 0.317156ter = m.gen["TE.random"]
data.gen = data.frame(x,y,ter)
ggplot(data = data.gen, mapping = aes(x=TE, y = seTE, color= )) +
geom_point() +
geom_vline(xintercept = 0.5401759)
ggplot(data = df_shape, mapping = aes(x=d, y = d_var, color= indoeuropean)) +
geom_point() +
geom_vline(xintercept = 0.5401759) +
scale_y_reverse() +
geom_smooth(method = "lm")## `geom_smooth()` using formula = 'y ~ x'
ggplot(data = df_shape, mapping = aes(x=d, y = d_var, color= language)) +
geom_point() +
geom_vline(xintercept = 0.5401759) +
scale_y_reverse() +
geom_smooth(method = "lm")## `geom_smooth()` using formula = 'y ~ x'
plotting coefficients from the rmv model: assuming that those coefficients correspond to effect sizes
Confirmatory analysis
For primary analyses, i will exclude effect sizes from clinical populations and multilingual populations.
I will investigate the hypotheses via multi-level meta-regressions using the metafor package.
In all models, I will include random effects that control for non-independence between effect sizes based on grouping by paper and grouping by experiment.
I will first fit: Shape bias ~ 1 Shape bias ~ age shape bias ~ log(age) shape bias ~ poly(age,2)
intercept:
# using the meta and metafor packages to analyze meta-analysis effect sizes
mod_intercept <- rma.mv(d ~ 1,
V = d_var,
random = ~1 | as.factor(Title) /
as.factor(exp_num),
slab = Title,
data = filter(df_shape, !is.na(exp_num)))
summary(mod_intercept)##
## Multivariate Meta-Analysis Model (k = 133; method: REML)
##
## logLik Deviance AIC BIC AICc
## -328.4943 656.9886 662.9886 671.6370 663.1761
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.1519 0.3897 26 no as.factor(Title)
## sigma^2.2 0.0027 0.0524 34 no as.factor(Title)/as.factor(exp_num)
##
## Test for Heterogeneity:
## Q(df = 132) = 895.1881, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.4776 0.0828 5.7691 <.0001 0.3154 0.6399 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1age
mod_age <- rma.mv(d ~ mean_age_months_centered36,
V = d_var,
random = ~1 | as.factor(Title) /
as.factor(exp_num),
slab = Title,
data = filter(df_shape, !is.na(exp_num)))## Warning: Rows with NAs omitted from model fitting.summary(mod_age)##
## Multivariate Meta-Analysis Model (k = 132; method: REML)
##
## logLik Deviance AIC BIC AICc
## -311.2327 622.4654 630.4654 641.9355 630.7854
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.1143 0.3380 25 no as.factor(Title)
## sigma^2.2 0.0078 0.0884 33 no as.factor(Title)/as.factor(exp_num)
##
## Test for Residual Heterogeneity:
## QE(df = 130) = 841.0606, p-val < .0001
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 27.8175, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.4438 0.0753 5.8971 <.0001 0.2963 0.5913
## mean_age_months_centered36 -0.0146 0.0028 -5.2742 <.0001 -0.0201 -0.0092
##
## intrcpt ***
## mean_age_months_centered36 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1log age
mod_log_age <- rma.mv(d ~ log_mean_age_months,
V = d_var,
random = ~1 | as.factor(Title) /
as.factor(exp_num),
slab = Title,
data = filter(df_shape, !is.na(log_mean_age_months)))
summary(mod_log_age)##
## Multivariate Meta-Analysis Model (k = 132; method: REML)
##
## logLik Deviance AIC BIC AICc
## -315.7124 631.4248 639.4248 650.8949 639.7448
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.1183 0.3439 25 no as.factor(Title)
## sigma^2.2 0.0059 0.0767 33 no as.factor(Title)/as.factor(exp_num)
##
## Test for Residual Heterogeneity:
## QE(df = 130) = 851.6190, p-val < .0001
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 18.6968, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 2.0909 0.3879 5.3904 <.0001 1.3307 2.8512 ***
## log_mean_age_months -0.4711 0.1089 -4.3240 <.0001 -0.6846 -0.2575 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Let’s look at what this means:
ggplot(df_shape,
aes(x = mean_age_months, y = d, color = language))+
geom_pointrange(aes(ymin = d - d_var, ymax = d + d_var),
alpha = .5) +
geom_smooth(aes(group = 1),
col = "black",
method = "lm", se = FALSE,
formula = y ~ log(x)) +
geom_hline(yintercept = 0, lty = 3) +
ylab("Standardized Mean Difference (d)") +
xlab("Mean age (months)") +
scale_color_discrete(name = "Language") +
scale_linetype_discrete(name = "Indo-European") +
theme(legend.position = "bottom") +
xlim(0,80)## Warning: Removed 8 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 8 rows containing missing values (`geom_pointrange()`).
polynoial age
mod_poly <- rma.mv(d ~ mean_age_months_centered36 + I(mean_age_months_centered36^2),
V = d_var,
random = ~1 | as.factor(Title) /
as.factor(exp_num),
slab = Title,
data = filter(df_shape, !is.na(log_mean_age_months)))
summary(mod_poly)##
## Multivariate Meta-Analysis Model (k = 132; method: REML)
##
## logLik Deviance AIC BIC AICc
## -303.1769 606.3538 616.3538 630.6528 616.8416
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0984 0.3138 25 no as.factor(Title)
## sigma^2.2 0.0423 0.2056 33 no as.factor(Title)/as.factor(exp_num)
##
## Test for Residual Heterogeneity:
## QE(df = 129) = 829.1803, p-val < .0001
##
## Test of Moderators (coefficients 2:3):
## QM(df = 2) = 45.6512, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb
## intrcpt 0.6197 0.0892 6.9485 <.0001 0.4449
## mean_age_months_centered36 -0.0048 0.0037 -1.2926 0.1961 -0.0121
## I(mean_age_months_centered36^2) -0.0007 0.0002 -4.1852 <.0001 -0.0011
## ci.ub
## intrcpt 0.7945 ***
## mean_age_months_centered36 0.0025
## I(mean_age_months_centered36^2) -0.0004 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model comparison and plotting AICc , what is the criteria ? cutoff
#anova(mod_log_age, mod_poly, refit = TRUE)
#anova(mod_age, mod_poly) ## the two models are not comparable, not nested
#plot_component(mod, type = "BIC")
summary(mod_intercept)$fit.stats[5,'REML']## [1] 663.1761AIC_vector <- c(summary(mod_intercept)$fit.stats[5,'REML'],
summary(mod_age)$fit.stats[5,'REML'],
summary(mod_log_age)$fit.stats[5,'REML'],
summary(mod_poly)$fit.stats[5,'REML'])
# AIC_vector <- c(AICc(mod_intercept), AICc(mod_age), AICc(mod_log_age), AICc(mod_poly))
modelnames <- c("intercept", "age", "logage", "polyage")
data.AIC = data.frame(AIC_vector, modelnames)
ggplot(data = data.AIC, mapping = aes(x = AIC_vector, y = modelnames)) +
geom_point()
Polynomial age models
Let’s start with an interaction with indoeuropean with
standard quadratic terms. This model is very interpretable.
rma.mv(d ~ mean_age_months_centered36 * indoeuropean +
I(mean_age_months_centered36^2) * indoeuropean,
V = d_var,
random = ~ 1 | ID/exp_num,
slab = short_cite,
data = filter(df_shape, !is.na(exp_num), !is.na(language)))##
## Multivariate Meta-Analysis Model (k = 129; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.1335 0.3654 23 no ID
## sigma^2.2 0.0308 0.1754 32 no ID/exp_num
##
## Test for Residual Heterogeneity:
## QE(df = 123) = 754.4436, p-val < .0001
##
## Test of Moderators (coefficients 2:6):
## QM(df = 5) = 91.5051, p-val < .0001
##
## Model Results:
##
## estimate se zval
## intrcpt 0.7034 0.0988 7.1163
## mean_age_months_centered36 0.0024 0.0042 0.5591
## indoeuropeanFALSE -0.4062 0.0960 -4.2318
## I(mean_age_months_centered36^2) -0.0010 0.0002 -4.5380
## mean_age_months_centered36:indoeuropeanFALSE -0.0338 0.0082 -4.1176
## indoeuropeanFALSE:I(mean_age_months_centered36^2) 0.0011 0.0003 3.4408
## pval ci.lb ci.ub
## intrcpt <.0001 0.5097 0.8972
## mean_age_months_centered36 0.5761 -0.0059 0.0107
## indoeuropeanFALSE <.0001 -0.5943 -0.2181
## I(mean_age_months_centered36^2) <.0001 -0.0014 -0.0005
## mean_age_months_centered36:indoeuropeanFALSE <.0001 -0.0498 -0.0177
## indoeuropeanFALSE:I(mean_age_months_centered36^2) 0.0006 0.0005 0.0017
##
## intrcpt ***
## mean_age_months_centered36
## indoeuropeanFALSE ***
## I(mean_age_months_centered36^2) ***
## mean_age_months_centered36:indoeuropeanFALSE ***
## indoeuropeanFALSE:I(mean_age_months_centered36^2) ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Next try breaking down by language. Here we can see Spanish is sparse and has a huge interaction term for some reason. Probably just overfit.
rma.mv(d ~ mean_age_months_centered36 * language +
I(mean_age_months_centered36^2) * language,
V = d_var,
random = ~ 1 | ID/exp_num,
slab = short_cite,
data = filter(df_shape, !is.na(exp_num), !is.na(language)))## Warning: Redundant predictors dropped from the model.##
## Multivariate Meta-Analysis Model (k = 129; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.2084 0.4565 23 no ID
## sigma^2.2 0.0310 0.1762 32 no ID/exp_num
##
## Test for Residual Heterogeneity:
## QE(df = 112) = 672.2426, p-val < .0001
##
## Test of Moderators (coefficients 2:17):
## QM(df = 16) = 180.4172, p-val < .0001
##
## Model Results:
##
## estimate se zval
## intrcpt 0.7275 0.5561 1.3081
## mean_age_months_centered36 -0.1196 0.1290 -0.9270
## languageenglish 0.0071 0.5691 0.0125
## languagegerman -1.0465 0.3533 -2.9618
## languagejapanese -0.1818 0.6055 -0.3002
## languagespanish 17.2964 3.2782 5.2762
## languagetsimane 1.4406 6.8350 0.2108
## I(mean_age_months_centered36^2) 0.0017 0.0052 0.3191
## mean_age_months_centered36:languageenglish 0.1178 0.1286 0.9157
## mean_age_months_centered36:languagegerman 0.1040 0.1436 0.7244
## mean_age_months_centered36:languagejapanese 0.1180 0.1295 0.9114
## mean_age_months_centered36:languagespanish 4.4113 0.8209 5.3737
## mean_age_months_centered36:languagetsimane -0.0114 0.2688 -0.0423
## languageenglish:I(mean_age_months_centered36^2) -0.0025 0.0052 -0.4757
## languagegerman:I(mean_age_months_centered36^2) -0.0011 0.0059 -0.1912
## languagejapanese:I(mean_age_months_centered36^2) -0.0022 0.0055 -0.3957
## languagespanish:I(mean_age_months_centered36^2) 0.2050 0.0407 5.0368
## pval ci.lb ci.ub
## intrcpt 0.1908 -0.3625 1.8175
## mean_age_months_centered36 0.3539 -0.3725 0.1333
## languageenglish 0.9900 -1.1083 1.1226
## languagegerman 0.0031 -1.7390 -0.3540
## languagejapanese 0.7641 -1.3686 1.0051
## languagespanish <.0001 10.8712 23.7216
## languagetsimane 0.8331 -11.9557 14.8369
## I(mean_age_months_centered36^2) 0.7496 -0.0086 0.0119
## mean_age_months_centered36:languageenglish 0.3598 -0.1343 0.3699
## mean_age_months_centered36:languagegerman 0.4688 -0.1774 0.3854
## mean_age_months_centered36:languagejapanese 0.3621 -0.1358 0.3718
## mean_age_months_centered36:languagespanish <.0001 2.8024 6.0203
## mean_age_months_centered36:languagetsimane 0.9663 -0.5382 0.5154
## languageenglish:I(mean_age_months_centered36^2) 0.6343 -0.0127 0.0078
## languagegerman:I(mean_age_months_centered36^2) 0.8483 -0.0126 0.0104
## languagejapanese:I(mean_age_months_centered36^2) 0.6923 -0.0128 0.0085
## languagespanish:I(mean_age_months_centered36^2) <.0001 0.1252 0.2848
##
## intrcpt
## mean_age_months_centered36
## languageenglish
## languagegerman **
## languagejapanese
## languagespanish ***
## languagetsimane
## I(mean_age_months_centered36^2)
## mean_age_months_centered36:languageenglish
## mean_age_months_centered36:languagegerman
## mean_age_months_centered36:languagejapanese
## mean_age_months_centered36:languagespanish ***
## mean_age_months_centered36:languagetsimane
## languageenglish:I(mean_age_months_centered36^2)
## languagegerman:I(mean_age_months_centered36^2)
## languagejapanese:I(mean_age_months_centered36^2)
## languagespanish:I(mean_age_months_centered36^2) ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1With the orthogonal polynomials, it blows up completely.
rma.mv(d ~ poly(mean_age_months_centered36,2) * language,
V = d_var,
random = ~ 1 | ID/exp_num,
slab = short_cite,
data = filter(df_shape, !is.na(exp_num), !is.na(language)))## Warning: Redundant predictors dropped from the model.##
## Multivariate Meta-Analysis Model (k = 129; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.2084 0.4565 23 no ID
## sigma^2.2 0.0310 0.1762 32 no ID/exp_num
##
## Test for Residual Heterogeneity:
## QE(df = 112) = 672.2426, p-val < .0001
##
## Test of Moderators (coefficients 2:17):
## QM(df = 16) = 180.4172, p-val < .0001
##
## Model Results:
##
## estimate se
## intrcpt 0.8310 1.3147
## poly(mean_age_months_centered36, 2)1 -15.6298 10.3513
## poly(mean_age_months_centered36, 2)2 4.0138 12.5770
## languageenglish -0.2591 1.3324
## languagegerman -1.0723 0.9990
## languagejapanese -0.3837 1.3428
## languagespanish 65.6068 12.6838
## languagetsimane 1.4194 6.3349
## poly(mean_age_months_centered36, 2)1:languageenglish 13.7848 10.3218
## poly(mean_age_months_centered36, 2)2:languageenglish -5.9723 12.5548
## poly(mean_age_months_centered36, 2)1:languagegerman 14.2181 11.5411
## poly(mean_age_months_centered36, 2)2:languagegerman -2.7003 14.1202
## poly(mean_age_months_centered36, 2)1:languagejapanese 14.4434 10.6063
## poly(mean_age_months_centered36, 2)2:languagejapanese -5.1865 13.1061
## poly(mean_age_months_centered36, 2)1:languagespanish 1086.8686 205.0889
## poly(mean_age_months_centered36, 2)2:languagespanish 492.7409 97.8288
## poly(mean_age_months_centered36, 2)1:languagetsimane -1.7898 42.3090
## zval pval
## intrcpt 0.6321 0.5273
## poly(mean_age_months_centered36, 2)1 -1.5099 0.1311
## poly(mean_age_months_centered36, 2)2 0.3191 0.7496
## languageenglish -0.1945 0.8458
## languagegerman -1.0734 0.2831
## languagejapanese -0.2857 0.7751
## languagespanish 5.1725 <.0001
## languagetsimane 0.2241 0.8227
## poly(mean_age_months_centered36, 2)1:languageenglish 1.3355 0.1817
## poly(mean_age_months_centered36, 2)2:languageenglish -0.4757 0.6343
## poly(mean_age_months_centered36, 2)1:languagegerman 1.2320 0.2180
## poly(mean_age_months_centered36, 2)2:languagegerman -0.1912 0.8483
## poly(mean_age_months_centered36, 2)1:languagejapanese 1.3618 0.1733
## poly(mean_age_months_centered36, 2)2:languagejapanese -0.3957 0.6923
## poly(mean_age_months_centered36, 2)1:languagespanish 5.2995 <.0001
## poly(mean_age_months_centered36, 2)2:languagespanish 5.0368 <.0001
## poly(mean_age_months_centered36, 2)1:languagetsimane -0.0423 0.9663
## ci.lb ci.ub
## intrcpt -1.7458 3.4078
## poly(mean_age_months_centered36, 2)1 -35.9178 4.6583
## poly(mean_age_months_centered36, 2)2 -20.6367 28.6643
## languageenglish -2.8705 2.3523
## languagegerman -3.0302 0.8857
## languagejapanese -3.0154 2.2481
## languagespanish 40.7471 90.4666 ***
## languagetsimane -10.9968 13.8355
## poly(mean_age_months_centered36, 2)1:languageenglish -6.4456 34.0152
## poly(mean_age_months_centered36, 2)2:languageenglish -30.5793 18.6347
## poly(mean_age_months_centered36, 2)1:languagegerman -8.4020 36.8382
## poly(mean_age_months_centered36, 2)2:languagegerman -30.3754 24.9748
## poly(mean_age_months_centered36, 2)1:languagejapanese -6.3445 35.2313
## poly(mean_age_months_centered36, 2)2:languagejapanese -30.8740 20.5011
## poly(mean_age_months_centered36, 2)1:languagespanish 684.9018 1488.8354 ***
## poly(mean_age_months_centered36, 2)2:languagespanish 300.9999 684.4819 ***
## poly(mean_age_months_centered36, 2)1:languagetsimane -84.7139 81.1343
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1