Our preregistered hypotheses:
library(tidyverse)
library(dplyr)
library(ggplot2)
library(nlme)
library(lme4)
library(ggpubr)
library(rstatix)
library(sjPlot) # for plotting lmer models
data <- read.csv("prop_look_data_20250612.csv")
data_included <- subset(data, SubjectInfo.Include == 1)
# Proplookdoor centered at chance = 0.5
data_included$propLookDoorCentered <- data_included$propLookDoor - 0.5
cat("Number of subjects:", length(unique(data_included$SubjectInfo.subjID)), "\n")
## Number of subjects: 105
cat("Number of navigators (Group 1):", length(unique(subset(data_included, Locomotor_group == 1)$SubjectInfo.subjID)), "\n")
## Number of navigators (Group 1): 37
cat("Number of pre-navigators (Group -1):", length(unique(subset(data_included, Locomotor_group == -1)$SubjectInfo.subjID)), "\n")
## Number of pre-navigators (Group -1): 68
cat("Overall proportion looking to door (mean):", mean(data_included$propLookDoor, na.rm=TRUE), "\n")
## Overall proportion looking to door (mean): 0.5195047
# Load the matched pairs data
matched_pairs <- read.csv("matched_pairs_mavew.csv")
head(matched_pairs)
## Walker_id Prewalker_id Walker_age Prewalker_age
## 1 SAX_MAVEW_2075 SAX_MAVEW_2094 366 354
## 2 SAX_MAVEW_2086 SAX_MAVEW_2030 384 371
## 3 SAX_MAVEW_2063 SAX_MAVEW_2021 397 397
## 4 SAX_MAVEW_2054 SAX_MAVEW_2015 385 373
## 5 SAX_MAVEW_2066 SAX_MAVEW_2012 404 402
## 6 SAX_MAVEW_2069 SAX_MAVEW_2097 378 365
# Calculate age difference for each pair
matched_pairs <- matched_pairs %>%
mutate(age_difference = abs(Walker_age - Prewalker_age))
# Verify all pairs meet the ≤14 days age difference criterion from preregistration
all_pairs_valid <- all(matched_pairs$age_difference <= 14)
cat("meet 14 day criterion?", all_pairs_valid, "\n")
## meet 14 day criterion? TRUE
# overall mean ages don't significantly differ (criterion from preregistration)
age_ttest <- t.test(matched_pairs$Walker_age, matched_pairs$Prewalker_age)
print(age_ttest)
##
## Welch Two Sample t-test
##
## data: matched_pairs$Walker_age and matched_pairs$Prewalker_age
## t = 0.64531, df = 63.629, p-value = 0.521
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -11.87809 23.21143
## sample estimates:
## mean of x mean of y
## 378.3636 372.6970
cat("Mean age difference between groups:",
mean(matched_pairs$Walker_age) - mean(matched_pairs$Prewalker_age), "days\n")
## Mean age difference between groups: 5.666667 days
matched_subject_ids <- c(matched_pairs$Walker_id, matched_pairs$Prewalker_id)
# matched participants only subset!
data_matched <- subset(data_included, SubjectInfo.subjID %in% matched_subject_ids)
# verify number of participants
# cat("Number of unique subjects in matched pairs:", length(unique(matched_subject_ids)), "\n")
# cat("Number of subjects in data_matched:", length(unique(data_matched$SubjectInfo.subjID)), "\n")
# pair identifier
pair_mapping <- rbind(
data.frame(SubjectInfo.subjID = matched_pairs$Walker_id, pair_id = 1:nrow(matched_pairs), Locomotor_group = 1),
data.frame(SubjectInfo.subjID = matched_pairs$Prewalker_id, pair_id = 1:nrow(matched_pairs), Locomotor_group = -1)
)
# join matched data
data_matched <- data_matched %>%
left_join(select(pair_mapping, SubjectInfo.subjID, pair_id), by = "SubjectInfo.subjID")
# look at age distribution in matched dataset
ggplot(data_matched, aes(x = SubjectInfo.Age, fill = factor(Locomotor_group))) +
geom_histogram(position = "dodge", bins = 10, alpha = 0.7) +
scale_fill_discrete(name = "Group", labels = c("-1" = "Pre-Navigators", "1" = "Navigators")) +
labs(title = "Age Distribution in Matched Sample", x = "Age (Days)", y = "Count") +
theme_minimal()
# trial average for full sample
data_trialavg <- data_included %>%
group_by(SubjectInfo.subjID, SubjectInfo.Age, SubjectInfo.Age_z, SubjectInfo.navExp, SubjectInfo.navExp_scaled, Locomotor_group) %>%
summarize(PropLookDoor = mean(propLookDoor, na.rm = TRUE), .groups = 'drop')
# trial average for matched pairs only
data_matched_trialavg <- data_matched %>%
group_by(SubjectInfo.subjID, SubjectInfo.Age, SubjectInfo.Age_z, SubjectInfo.navExp, SubjectInfo.navExp_scaled, Locomotor_group, pair_id) %>%
summarize(PropLookDoor = mean(propLookDoor, na.rm = TRUE), .groups = 'drop')
# different analysis groups
navigators <- subset(data_included, Locomotor_group == 1)
prenavigators <- subset(data_included, Locomotor_group == -1)
# trial-averaged subsets for t-tests
navigator_look <- navigators %>%
group_by(SubjectInfo.subjID) %>%
summarize(mean_propLookDoor = mean(propLookDoor, na.rm = TRUE))
prenavigator_look <- prenavigators %>%
group_by(SubjectInfo.subjID) %>%
summarize(mean_propLookDoor = mean(propLookDoor, na.rm = TRUE))
# matched-only versions
matched_navigators <- subset(data_matched, Locomotor_group == 1)
matched_prenavigators <- subset(data_matched, Locomotor_group == -1)
matched_navigator_look <- matched_navigators %>%
group_by(SubjectInfo.subjID) %>%
summarize(mean_propLookDoor = mean(propLookDoor, na.rm = TRUE))
matched_prenavigator_look <- matched_prenavigators %>%
group_by(SubjectInfo.subjID) %>%
summarize(mean_propLookDoor = mean(propLookDoor, na.rm = TRUE))
# check!!! sample sizes
# cat("FULL DATASET:\n")
# cat("Number of navigator subjects:", length(unique(navigators$SubjectInfo.subjID)), "\n")
# cat("Number of pre-navigator subjects:", length(unique(prenavigators$SubjectInfo.subjID)), "\n\n")
#
# cat("MATCHED DATASET:\n")
# cat("Number of matched navigator subjects:", length(unique(matched_navigators$SubjectInfo.subjID)), "\n")
# cat("Number of matched pre-navigator subjects:", length(unique(matched_prenavigators$SubjectInfo.subjID)), "\n")
# FULL DATASET
p_age_full <- ggplot(data_trialavg,
aes(x = SubjectInfo.Age, y = PropLookDoor, color = factor(Locomotor_group))) +
geom_point(size = 3, alpha = 0.7) +
geom_smooth(method = "lm", se = TRUE, alpha = 0.2) +
geom_hline(yintercept = 0.5, linetype = "dashed", color = "gray50") +
scale_color_manual(values = c("-1" = "red", "1" = "blue"),
labels = c("-1" = "Pre-Navigators", "1" = "Navigators"),
name = "Group") +
labs(title = "Age vs. Proportion Looking to Door",
x = "Age (Days)",
y = "Proportion Looking to Door") +
ylim(0, 1) +
theme_minimal() +
theme(
legend.position = "bottom",
axis.text = element_text(size = 8),
axis.title = element_text(size = 8),
plot.title = element_text(size = 12, hjust = 0.5)
)
# MATCHED DATASET
p_age_matched <- ggplot(data_matched_trialavg,
aes(x = SubjectInfo.Age, y = PropLookDoor, color = factor(Locomotor_group))) +
geom_point(size = 3, alpha = 0.7) +
geom_smooth(method = "lm", se = TRUE, alpha = 0.2) +
geom_hline(yintercept = 0.5, linetype = "dashed", color = "gray50") +
scale_color_manual(values = c("-1" = "red", "1" = "blue"),
labels = c("-1" = "Pre-Navigators", "1" = "Navigators"),
name = "Group") +
labs(title = "Age vs. Proportion Looking to Door (Matched)",
x = "Age (Days)",
y = "Proportion Looking to Door") +
ylim(0, 1) +
theme_minimal() +
theme(
legend.position = "bottom",
axis.text = element_text(size = 8),
axis.title = element_text(size = 8),
plot.title = element_text(size = 12, hjust = 0.5)
)
# Arrange both visualizations
ggarrange(p_age_full, p_age_matched, ncol = 2, nrow = 1, labels = "AUTO", common.legend = TRUE)
# correlation tests - Full
cor_age_all <- cor.test(
data_trialavg$SubjectInfo.Age_z,
data_trialavg$PropLookDoor
)
cor_age_navigators <- cor.test(
subset(data_trialavg, Locomotor_group == 1)$SubjectInfo.Age_z,
subset(data_trialavg, Locomotor_group == 1)$PropLookDoor
)
cor_age_prenavigators <- cor.test(
subset(data_trialavg, Locomotor_group == -1)$SubjectInfo.Age_z,
subset(data_trialavg, Locomotor_group == -1)$PropLookDoor
)
cat("Correlation between age and looking behavior (all subjects):\n")
## Correlation between age and looking behavior (all subjects):
print(cor_age_all)
##
## Pearson's product-moment correlation
##
## data: data_trialavg$SubjectInfo.Age_z and data_trialavg$PropLookDoor
## t = 2.8573, df = 103, p-value = 0.005171
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.08367993 0.43981932
## sample estimates:
## cor
## 0.2709988
cat("\nCorrelation between age and looking behavior (navigators only):\n")
##
## Correlation between age and looking behavior (navigators only):
print(cor_age_navigators)
##
## Pearson's product-moment correlation
##
## data: subset(data_trialavg, Locomotor_group == 1)$SubjectInfo.Age_z and subset(data_trialavg, Locomotor_group == 1)$PropLookDoor
## t = -0.50205, df = 35, p-value = 0.6188
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.3976809 0.2462070
## sample estimates:
## cor
## -0.08455765
cat("\nCorrelation between age and looking behavior (pre-navigators only):\n")
##
## Correlation between age and looking behavior (pre-navigators only):
print(cor_age_prenavigators)
##
## Pearson's product-moment correlation
##
## data: subset(data_trialavg, Locomotor_group == -1)$SubjectInfo.Age_z and subset(data_trialavg, Locomotor_group == -1)$PropLookDoor
## t = 2.9204, df = 66, p-value = 0.004779
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1086198 0.5336662
## sample estimates:
## cor
## 0.3382842
mod2 <- lme(propLookDoorCentered ~ Locomotor_group + SubjectInfo.navExp_scaled + SubjectInfo.Age_z +
Texture_side + Trial_num,
random = ~1|SubjectInfo.subjID,
data = data_included,
na.action = na.exclude)
summary(mod2)
## Linear mixed-effects model fit by REML
## Data: data_included
## AIC BIC logLik
## 317.9233 361.1527 -150.9617
##
## Random effects:
## Formula: ~1 | SubjectInfo.subjID
## (Intercept) Residual
## StdDev: 9.644148e-06 0.262241
##
## Fixed effects: propLookDoorCentered ~ Locomotor_group + SubjectInfo.navExp_scaled + SubjectInfo.Age_z + Texture_side + Trial_num
## Value Std.Error DF t-value p-value
## (Intercept) 0.03068222 0.01451108 1542 2.114399 0.0346
## Locomotor_group 0.00925948 0.01036933 100 0.892968 0.3740
## SubjectInfo.navExp_scaled 0.00640425 0.04915724 100 0.130281 0.8966
## SubjectInfo.Age_z 0.00956277 0.00709441 100 1.347931 0.1807
## Texture_side 0.09092640 0.00646695 1542 14.060167 0.0000
## Trial_num -0.01795647 0.02107077 1542 -0.852198 0.3942
## Correlation:
## (Intr) Lcmtr_ SbI.E_ SbI.A_ Txtr_s
## Locomotor_group 0.465
## SubjectInfo.navExp_scaled -0.503 -0.692
## SubjectInfo.Age_z -0.052 -0.280 -0.035
## Texture_side 0.033 0.002 0.001 -0.005
## Trial_num -0.726 -0.006 0.004 0.003 -0.047
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.429025073 -0.703811236 0.007309415 0.733371584 2.306628130
##
## Number of Observations: 1648
## Number of Groups: 104
plot_model(mod2, show.values = TRUE)
mod2_simple <- lm(PropLookDoor - 0.5 ~ Locomotor_group, data = data_trialavg)
# manually divide p for one-tailed
summary(mod2_simple)
##
## Call:
## lm(formula = PropLookDoor - 0.5 ~ Locomotor_group, data = data_trialavg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.150095 -0.030170 0.002598 0.030617 0.128895
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.023779 0.005199 4.574 1.34e-05 ***
## Locomotor_group 0.014034 0.005199 2.699 0.00813 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0509 on 103 degrees of freedom
## Multiple R-squared: 0.06606, Adjusted R-squared: 0.05699
## F-statistic: 7.285 on 1 and 103 DF, p-value: 0.008126
mod2_grouponly <- lm(PropLookDoor - 0.5 ~ Locomotor_group, data = data_trialavg)
mod2_group_age <- lm(PropLookDoor - 0.5 ~ Locomotor_group + SubjectInfo.Age_z, data = data_trialavg)
# mod2_group_age_navexp <- lm(PropLookDoor - 0.5 ~ Locomotor_group + SubjectInfo.Age_z + SubjectInfo.navExp_scaled, data = data_trialavg)
# compare models
anova(mod2_grouponly, mod2_group_age)
## Analysis of Variance Table
##
## Model 1: PropLookDoor - 0.5 ~ Locomotor_group
## Model 2: PropLookDoor - 0.5 ~ Locomotor_group + SubjectInfo.Age_z
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 103 0.26687
## 2 102 0.25749 1 0.0093795 3.7155 0.05669 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(summary(mod2_grouponly)$coefficients["Locomotor_group",])
## Estimate Std. Error t value Pr(>|t|)
## 0.014033530 0.005199245 2.699147560 0.008126474
print(summary(mod2_group_age)$coefficients["Locomotor_group",])
## Estimate Std. Error t value Pr(>|t|)
## 0.009558971 0.005632623 1.697072615 0.092732420
# print(summary(mod2_group_age_navexp)$coefficients["Locomotor_group",])
texture_data <- data_included %>%
group_by(Texture_side, Locomotor_group) %>%
summarize(mean_propLookDoor = mean(propLookDoor, na.rm = TRUE), .groups = 'drop')
ggplot(texture_data, aes(x = factor(Texture_side), y = mean_propLookDoor, fill = factor(Locomotor_group))) +
geom_bar(stat = "identity", position = position_dodge(), alpha = 0.7) +
geom_hline(yintercept = 0.5, linetype = "dashed", color = "gray50") +
scale_fill_manual(values = c("-1" = "red", "1" = "blue"),
labels = c("-1" = "Pre-Navigators", "1" = "Navigators"),
name = "Group") +
scale_x_discrete(labels = c("-1" = "Low SF on Door", "1" = "High SF on Door")) +
labs(title = "Effect of Texture Side on Looking Behavior",
x = "Texture Side",
y = "Mean Proportion Looking to Door") +
ylim(0, 1) +
theme_minimal() +
theme(
legend.position = "bottom",
axis.text = element_text(size = 12),
axis.title = element_text(size = 14),
plot.title = element_text(size = 16, hjust = 0.5)
)
# texture effect within each group
texture_effect_nav <- t.test(
subset(data_included, Locomotor_group == 1 & Texture_side == 1)$propLookDoor,
subset(data_included, Locomotor_group == 1 & Texture_side == -1)$propLookDoor
)
texture_effect_prenav <- t.test(
subset(data_included, Locomotor_group == -1 & Texture_side == 1)$propLookDoor,
subset(data_included, Locomotor_group == -1 & Texture_side == -1)$propLookDoor
)
cat("Texture side effect for navigators:\n")
## Texture side effect for navigators:
print(texture_effect_nav)
##
## Welch Two Sample t-test
##
## data: subset(data_included, Locomotor_group == 1 & Texture_side == 1)$propLookDoor and subset(data_included, Locomotor_group == 1 & Texture_side == -1)$propLookDoor
## t = 8.372, df = 581.5, p-value = 4.238e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.1385091 0.2234152
## sample estimates:
## mean of x mean of y
## 0.6276761 0.4467140
cat("\nTexture side effect for pre-navigators:\n")
##
## Texture side effect for pre-navigators:
print(texture_effect_prenav)
##
## Welch Two Sample t-test
##
## data: subset(data_included, Locomotor_group == -1 & Texture_side == 1)$propLookDoor and subset(data_included, Locomotor_group == -1 & Texture_side == -1)$propLookDoor
## t = 11.346, df = 1075.7, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.1496523 0.2122360
## sample estimates:
## mean of x mean of y
## 0.6000246 0.4190804
# separating models for locomotor group and navigation experience
# competition between categorical and continuous navigation predictors
# only locomotor group, no navigation experience)
mod2a <- lme(propLookDoorCentered ~ Locomotor_group + SubjectInfo.Age_z + Texture_side + Trial_num,
random = ~1|SubjectInfo.subjID,
data = data_included,
na.action = na.exclude)
# only navigation experience, no locomotor group
mod2b <- lme(propLookDoorCentered ~ SubjectInfo.navExp_scaled + SubjectInfo.Age_z + Texture_side + Trial_num,
random = ~1|SubjectInfo.subjID,
data = data_included,
na.action = na.exclude)
summary(mod2a)
## Linear mixed-effects model fit by REML
## Data: data_included
## AIC BIC logLik
## 307.2191 345.1169 -146.6096
##
## Random effects:
## Formula: ~1 | SubjectInfo.subjID
## (Intercept) Residual
## StdDev: 9.50357e-06 0.2615977
##
## Fixed effects: propLookDoorCentered ~ Locomotor_group + SubjectInfo.Age_z + Texture_side + Trial_num
## Value Std.Error DF t-value p-value
## (Intercept) 0.03064666 0.012429078 1557 2.465723 0.0138
## Locomotor_group 0.00953016 0.007363572 102 1.294231 0.1985
## SubjectInfo.Age_z 0.01000484 0.007031424 102 1.422875 0.1578
## Texture_side 0.09067458 0.006419705 1557 14.124416 0.0000
## Trial_num -0.01706587 0.020916382 1557 -0.815909 0.4147
## Correlation:
## (Intr) Lcmtr_ SbI.A_ Txtr_s
## Locomotor_group 0.177
## SubjectInfo.Age_z -0.072 -0.411
## Texture_side 0.038 0.004 -0.004
## Trial_num -0.839 -0.004 0.003 -0.046
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.42956770 -0.70529570 0.00478191 0.73681272 2.31218482
##
## Number of Observations: 1664
## Number of Groups: 105
summary(mod2b)
## Linear mixed-effects model fit by REML
## Data: data_included
## AIC BIC logLik
## 309.4208 347.2508 -147.7104
##
## Random effects:
## Formula: ~1 | SubjectInfo.subjID
## (Intercept) Residual
## StdDev: 9.608156e-06 0.2622249
##
## Fixed effects: propLookDoorCentered ~ SubjectInfo.navExp_scaled + SubjectInfo.Age_z + Texture_side + Trial_num
## Value Std.Error DF t-value p-value
## (Intercept) 0.02466102 0.01284852 1542 1.919367 0.0551
## SubjectInfo.navExp_scaled 0.03678079 0.03548366 101 1.036556 0.3024
## SubjectInfo.Age_z 0.01133360 0.00681119 101 1.663968 0.0992
## Texture_side 0.09091401 0.00646654 1542 14.059150 0.0000
## Trial_num -0.01785011 0.02106913 1542 -0.847216 0.3970
## Correlation:
## (Intr) SbI.E_ SbI.A_ Txtr_s
## SubjectInfo.navExp_scaled -0.284
## SubjectInfo.Age_z 0.092 -0.329
## Texture_side 0.036 0.003 -0.004
## Trial_num -0.817 0.001 0.002 -0.047
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.49522075 -0.70090407 0.01456427 0.73560325 2.30627875
##
## Number of Observations: 1648
## Number of Groups: 104
cat("Model 2 with both predictors:\n",
paste0("Locomotor_group: β = ", round(summary(mod2)$tTable["Locomotor_group", "Value"], 3),
", t = ", round(summary(mod2)$tTable["Locomotor_group", "t-value"], 2),
", p = ", round(summary(mod2)$tTable["Locomotor_group", "p-value"], 3)), "\n",
paste0("Navigation Experience: β = ", round(summary(mod2)$tTable["SubjectInfo.navExp_scaled", "Value"], 5),
", t = ", round(summary(mod2)$tTable["SubjectInfo.navExp_scaled", "t-value"], 2),
", p = ", round(summary(mod2)$tTable["SubjectInfo.navExp_scaled", "p-value"], 3)), "\n\n")
## Model 2 with both predictors:
## Locomotor_group: β = 0.009, t = 0.89, p = 0.374
## Navigation Experience: β = 0.0064, t = 0.13, p = 0.897
cat("Model 2 with only Locomotor_group:\n",
paste0("Locomotor_group: β = ", round(summary(mod2a)$tTable["Locomotor_group", "Value"], 3),
", t = ", round(summary(mod2a)$tTable["Locomotor_group", "t-value"], 2),
", p = ", round(summary(mod2a)$tTable["Locomotor_group", "p-value"], 3)), "\n\n")
## Model 2 with only Locomotor_group:
## Locomotor_group: β = 0.01, t = 1.29, p = 0.199
cat("Model 2 with only Navigation Experience:\n",
paste0("Navigation Experience: β = ", round(summary(mod2b)$tTable["SubjectInfo.navExp_scaled", "Value"], 5),
", t = ", round(summary(mod2b)$tTable["SubjectInfo.navExp_scaled", "t-value"], 2),
", p = ", round(summary(mod2b)$tTable["SubjectInfo.navExp_scaled", "p-value"], 3)), "\n\n")
## Model 2 with only Navigation Experience:
## Navigation Experience: β = 0.03678, t = 1.04, p = 0.302
one_tailed_p_locomotor <- summary(mod2a)$tTable["Locomotor_group", "p-value"] / 2
one_tailed_p_navexp <- summary(mod2b)$tTable["SubjectInfo.navExp_scaled", "p-value"] / 2
cat("one-tailed p-values:\n")
## one-tailed p-values:
cat("Locomotor group only model:", round(one_tailed_p_locomotor, 3), "\n")
## Locomotor group only model: 0.099
cat("Navigation experience only model:", round(one_tailed_p_navexp, 3), "\n\n")
## Navigation experience only model: 0.151