1. Setup

Our preregistered hypotheses:

  1. Infants with active navigational walking experience (navigators) will attend more to parts of a scene that afford navigation (open doorway) compared to parts that do not (paintings).
  2. Attention to navigational affordances will increase with individual experience with active navigation.
  3. If navigational affordances can be learned from passive experience, looking behavior would relate to chronological age even after accounting for navigational experience.
library(tidyverse)
library(dplyr)
library(ggplot2)
library(nlme)
library(lme4)
library(ggpubr)
library(rstatix)
library(sjPlot) # for plotting lmer models

2. Data Preparation

2.1 Load Data From Standard Pipeline

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

2.2 Integrate Navigation Experience Data

nav_exp_data <- read.csv("MAVEW_participants_with_exp.csv")
head(nav_exp_data)
##               ID      Group Age WalkExp CrawlExp
## 1 SAX_MAVEW_2001     Walker 394      10      154
## 2 SAX_MAVEW_2002 Pre-walker 401      NA      114
## 3 SAX_MAVEW_2003     Walker 363      44      163
## 4 SAX_MAVEW_2004 Pre-walker 277      NA       49
## 5 SAX_MAVEW_2005 Pre-walker 382      NA      142
## 6 SAX_MAVEW_2006 Pre-walker 299      NA      102
# left join to keep participants from prop looking data
data_with_exp <- data_included %>%
  left_join(nav_exp_data, by = c("SubjectInfo.subjID" = "ID"))

data_with_exp <- data_with_exp %>%
  mutate(
    # walking experience is 0 for pre-walkers (instead of 'NA')
    WalkExp = ifelse(Locomotor_group == -1, 0, WalkExp),
    CrawlExp = ifelse(is.na(CrawlExp), 0, CrawlExp), 
    # flag subjects with missing experience data
    exp_data_missing = is.na(WalkExp) & Locomotor_group == 1
  )

# check!!
matched_count <- sum(!is.na(data_with_exp$WalkExp) | data_with_exp$Locomotor_group == -1)
total_count <- nrow(data_with_exp)
cat(matched_count, "out of", total_count, "\n")
## 1669 out of 1685
# upgrade main dataset to have navigation experience info
data_included <- data_with_exp %>%
  rename(
    SubjectInfo.navExp = WalkExp,
    SubjectInfo.crawlExp = CrawlExp
  ) %>%
  mutate(
    # Age from nav_exp_data if available, otherwise use existing age data (they are the same)
    SubjectInfo.Age = ifelse(is.na(Age), SubjectInfo.testAge, Age)
  )

data_included <- data_included %>%
  mutate(
    # z-score age, mean zero unit variance
    SubjectInfo.Age_z = as.vector(scale(SubjectInfo.Age)),
    # scale navigation experience to 0-1, 0 for pre-walkers
    SubjectInfo.navExp_scaled = ifelse(Locomotor_group == -1, 0,
                                       ifelse(is.na(SubjectInfo.navExp), NA,
                                              (SubjectInfo.navExp - min(SubjectInfo.navExp[!is.na(SubjectInfo.navExp)], na.rm = TRUE)) / 
                                              (max(SubjectInfo.navExp, na.rm = TRUE) - min(SubjectInfo.navExp[!is.na(SubjectInfo.navExp)], na.rm=TRUE))))
  )

navigator_exp_summary <- data_included %>%
  filter(Locomotor_group == 1) %>%
  summarize(
    avg_nav_exp = mean(SubjectInfo.navExp, na.rm = TRUE),
    min_nav_exp = min(SubjectInfo.navExp, na.rm = TRUE),
    max_nav_exp = max(SubjectInfo.navExp, na.rm = TRUE)
  )

cat("Navigator experience (days): Mean =", navigator_exp_summary$avg_nav_exp, 
    ", Min =", navigator_exp_summary$min_nav_exp, 
    ", Max =", navigator_exp_summary$max_nav_exp, "\n")
## Navigator experience (days): Mean = 46.25997 , Min = 3 , Max = 156
# write.csv(data_with_exp, file = "experiment_1_mavew.csv", row.names = FALSE)

2.3 Integrate Age-Matched Pairs Data

# 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()

2.4 Create Derived Datasets for Analysis

# 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")

3. Hypothesis 1: Navigators vs. Pre-navigators Group Differences

3.1 Group Comparison Tests

# Full sample group comparison
full_group_diff <- t.test(
  subset(data_trialavg, Locomotor_group == 1)$PropLookDoor,
  subset(data_trialavg, Locomotor_group == -1)$PropLookDoor,
  alternative = "greater"
)

# Matched sample group comparison
matched_group_diff <- t.test(
  subset(data_matched_trialavg, Locomotor_group == 1)$PropLookDoor,
  subset(data_matched_trialavg, Locomotor_group == -1)$PropLookDoor,
  alternative = "greater"
)

cat("Group difference in full dataset (Navigators vs Pre-navigators):\n")
## Group difference in full dataset (Navigators vs Pre-navigators):
print(full_group_diff)
## 
##  Welch Two Sample t-test
## 
## data:  subset(data_trialavg, Locomotor_group == 1)$PropLookDoor and subset(data_trialavg, Locomotor_group == -1)$PropLookDoor
## t = 2.5892, df = 65.711, p-value = 0.005917
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
##  0.009981952         Inf
## sample estimates:
## mean of x mean of y 
## 0.5378124 0.5097454
cat("\nGroup difference in matched dataset (Navigators vs Pre-navigators):\n")
## 
## Group difference in matched dataset (Navigators vs Pre-navigators):
print(matched_group_diff)
## 
##  Welch Two Sample t-test
## 
## data:  subset(data_matched_trialavg, Locomotor_group == 1)$PropLookDoor and subset(data_matched_trialavg, Locomotor_group == -1)$PropLookDoor
## t = 2.1638, df = 63.768, p-value = 0.01711
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
##  0.006372771         Inf
## sample estimates:
## mean of x mean of y 
## 0.5415431 0.5136708

3.2 Group Differences Visualization

# visualize the full dataset?
ggplot(data_trialavg %>% filter(Locomotor_group %in% c(-1, 1)), 
       aes(x = factor(Locomotor_group), y = PropLookDoor, fill = factor(Locomotor_group))) +
  geom_boxplot(width = 0.5, alpha = 0.7) +
  geom_jitter(width = 0.2, alpha = 0.5) +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "gray50") +
  scale_x_discrete(labels = c("-1" = "Pre-Navigators", "1" = "Navigators")) +
  scale_fill_manual(values = c("-1" = "red", "1" = "blue"), 
                    labels = c("-1" = "Pre-Navigators", "1" = "Navigators"),
                    name = "Group") +
  labs(title = "Looking to Door Side by Group",
       x = "Group",
       y = "Proportion Looking to Door Side") +
  ylim(0, 1) +
  theme_minimal() +
  theme(
    legend.position = "none",
    axis.text = element_text(size = 12),
    axis.title = element_text(size = 14),
    plot.title = element_text(size = 16, hjust = 0.5)
  )

ggplot(data_matched_trialavg, 
       aes(x = factor(Locomotor_group), y = PropLookDoor, fill = factor(Locomotor_group))) +
  geom_boxplot(width = 0.5, alpha = 0.7) +
  geom_jitter(width = 0.2, alpha = 0.5) +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "gray50") +
  scale_x_discrete(labels = c("-1" = "Pre-Navigators", "1" = "Navigators")) +
  scale_fill_manual(values = c("-1" = "red", "1" = "blue"), 
                    labels = c("-1" = "Pre-Navigators", "1" = "Navigators"),
                    name = "Group") +
  labs(title = "Looking to Door Side by Group (Matched Pairs Only)",
       x = "Group",
       y = "Proportion Looking to Door Side") +
  ylim(0, 1) +
  theme_minimal() +
  theme(
    legend.position = "none",
    axis.text = element_text(size = 12),
    axis.title = element_text(size = 14),
    plot.title = element_text(size = 16, hjust = 0.5)
  )

3.3 Comparison to Chance Level

# T-tests for the full dataset
navigators_test <- t.test(navigator_look$mean_propLookDoor, mu = 0.5, alternative = "greater")
prenavigators_test <- t.test(prenavigator_look$mean_propLookDoor, mu = 0.5, alternative = "greater")
all_subjects_test <- t.test(data_trialavg$PropLookDoor, mu = 0.5, alternative = "greater")

# T-tests for the matched pairs only
matched_navigators_test <- t.test(matched_navigator_look$mean_propLookDoor, mu = 0.5, alternative = "greater")
matched_prenavigators_test <- t.test(matched_prenavigator_look$mean_propLookDoor, mu = 0.5, alternative = "greater")
all_matched_test <- t.test(data_matched_trialavg$PropLookDoor, mu = 0.5, alternative = "greater")

cat("Navigators looking to door vs. chance (0.5):\n")
## Navigators looking to door vs. chance (0.5):
print(navigators_test)
## 
##  One Sample t-test
## 
## data:  navigator_look$mean_propLookDoor
## t = 4.1436, df = 36, p-value = 9.893e-05
## alternative hypothesis: true mean is greater than 0.5
## 95 percent confidence interval:
##  0.5224058       Inf
## sample estimates:
## mean of x 
## 0.5378124
cat("\nPre-navigators looking to door vs. chance (0.5):\n")
## 
## Pre-navigators looking to door vs. chance (0.5):
print(prenavigators_test)
## 
##  One Sample t-test
## 
## data:  prenavigator_look$mean_propLookDoor
## t = 1.6657, df = 67, p-value = 0.05022
## alternative hypothesis: true mean is greater than 0.5
## 95 percent confidence interval:
##  0.4999872       Inf
## sample estimates:
## mean of x 
## 0.5097454
cat("\nAll subjects looking to door vs. chance (0.5):\n")
## 
## All subjects looking to door vs. chance (0.5):
print(all_subjects_test)
## 
##  One Sample t-test
## 
## data:  data_trialavg$PropLookDoor
## t = 3.8385, df = 104, p-value = 0.0001064
## alternative hypothesis: true mean is greater than 0.5
## 95 percent confidence interval:
##  0.511146      Inf
## sample estimates:
## mean of x 
## 0.5196357
cat("Matched navigators looking to door vs. chance (0.5):\n")
## Matched navigators looking to door vs. chance (0.5):
print(matched_navigators_test)
## 
##  One Sample t-test
## 
## data:  matched_navigator_look$mean_propLookDoor
## t = 4.4295, df = 32, p-value = 5.168e-05
## alternative hypothesis: true mean is greater than 0.5
## 95 percent confidence interval:
##  0.5256567       Inf
## sample estimates:
## mean of x 
## 0.5415431
cat("\nMatched pre-navigators looking to door vs. chance (0.5):\n")
## 
## Matched pre-navigators looking to door vs. chance (0.5):
print(matched_prenavigators_test)
## 
##  One Sample t-test
## 
## data:  matched_prenavigator_look$mean_propLookDoor
## t = 1.5483, df = 32, p-value = 0.06569
## alternative hypothesis: true mean is greater than 0.5
## 95 percent confidence interval:
##  0.4987147       Inf
## sample estimates:
## mean of x 
## 0.5136708
cat("\nAll matched subjects looking to door vs. chance (0.5):\n")
## 
## All matched subjects looking to door vs. chance (0.5):
print(all_matched_test)
## 
##  One Sample t-test
## 
## data:  data_matched_trialavg$PropLookDoor
## t = 4.17, df = 65, p-value = 4.602e-05
## alternative hypothesis: true mean is greater than 0.5
## 95 percent confidence interval:
##  0.5165599       Inf
## sample estimates:
## mean of x 
## 0.5276069

3.4 Preregistered Model 1: Age-matched Groups

# Age matched groups, testing categorical effect of locomotor status
# uses ONLY the matched participants!!
mod1 <- lme(propLookDoorCentered ~ Locomotor_group + SubjectInfo.Age_z + Texture_side + Trial_num, 
            random = ~1|SubjectInfo.subjID, 
            data = data_matched, 
            na.action = na.exclude)
summary(mod1)
## Linear mixed-effects model fit by REML
##   Data: data_matched 
##        AIC      BIC    logLik
##   234.9489 269.5844 -110.4744
## 
## Random effects:
##  Formula: ~1 | SubjectInfo.subjID
##          (Intercept) Residual
## StdDev: 1.003309e-05 0.264963
## 
## Fixed effects:  propLookDoorCentered ~ Locomotor_group + SubjectInfo.Age_z +      Texture_side + Trial_num 
##                         Value   Std.Error  DF   t-value p-value
## (Intercept)        0.02945582 0.016161809 978  1.822557  0.0687
## Locomotor_group    0.01300359 0.008219527  63  1.582036  0.1186
## SubjectInfo.Age_z  0.00670922 0.010208549  63  0.657216  0.5134
## Texture_side       0.08864451 0.008208915 978 10.798565  0.0000
## Trial_num         -0.01008598 0.026733886 978 -0.377273  0.7061
##  Correlation: 
##                   (Intr) Lcmtr_ SbI.A_ Txtr_s
## Locomotor_group    0.022                     
## SubjectInfo.Age_z -0.251 -0.081              
## Texture_side       0.052  0.003 -0.009       
## Trial_num         -0.826 -0.003  0.007 -0.062
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.38999625 -0.68679486  0.02036072  0.71513657  2.24178008 
## 
## Number of Observations: 1046
## Number of Groups: 66
# I don't know how to best plot this
plot_model(mod1, show.values = TRUE)

3.5 Model 1 Subject-Level Sanity Check

mod1_subject <- lm(PropLookDoor - 0.5 ~ Locomotor_group + SubjectInfo.Age_z, 
                   data = data_matched_trialavg)
summary(mod1_subject)
## 
## Call:
## lm(formula = PropLookDoor - 0.5 ~ Locomotor_group + SubjectInfo.Age_z, 
##     data = data_matched_trialavg)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.133233 -0.033025 -0.001341  0.029243  0.125593 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)   
## (Intercept)       0.024463   0.007182   3.406  0.00115 **
## Locomotor_group   0.013422   0.006462   2.077  0.04189 * 
## SubjectInfo.Age_z 0.007988   0.008069   0.990  0.32598   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05233 on 63 degrees of freedom
## Multiple R-squared:  0.08245,    Adjusted R-squared:  0.05332 
## F-statistic:  2.83 on 2 and 63 DF,  p-value: 0.06651
one_tailed_p_subject <- summary(mod1_subject)$coefficients["Locomotor_group", "Pr(>|t|)"] / 2
cat("One-tailed p-value (subject-level model):", round(one_tailed_p_subject, 3), "\n\n")
## One-tailed p-value (subject-level model): 0.021
plot_model(mod1_subject, title = "mod1 subject-level", show.values = TRUE)

3.6 Model 1 Without Age (Exploratory)

# Model 1 without age in age-matched sample
# redundant age control after matching?

mod1_no_age <- lme(propLookDoorCentered ~ Locomotor_group + Texture_side + Trial_num, 
                  random = ~1|SubjectInfo.subjID, 
                  data = data_matched, 
                  na.action = na.exclude)

summary(mod1_no_age)
## Linear mixed-effects model fit by REML
##   Data: data_matched 
##        AIC      BIC    logLik
##   226.0495 255.7429 -107.0247
## 
## Random effects:
##  Formula: ~1 | SubjectInfo.subjID
##          (Intercept)  Residual
## StdDev: 9.906563e-06 0.2648908
## 
## Fixed effects:  propLookDoorCentered ~ Locomotor_group + Texture_side + Trial_num 
##                       Value   Std.Error  DF   t-value p-value
## (Intercept)      0.03211938 0.015641146 978  2.053519  0.0403
## Locomotor_group  0.01344049 0.008190366  64  1.641012  0.1057
## Texture_side     0.08869375 0.008206336 978 10.807960  0.0000
## Trial_num       -0.01020255 0.026726010 978 -0.381746  0.7027
##  Correlation: 
##                 (Intr) Lcmtr_ Txtr_s
## Locomotor_group  0.002              
## Texture_side     0.051  0.002       
## Trial_num       -0.852 -0.003 -0.062
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.3918288 -0.6821729  0.0204076  0.7085152  2.1904022 
## 
## Number of Observations: 1046
## Number of Groups: 66
cat("Model 1 with age:", 
    paste0("β = ", round(summary(mod1)$tTable["Locomotor_group", "Value"], 3),
           ", t = ", round(summary(mod1)$tTable["Locomotor_group", "t-value"], 2),
           ", p = ", round(summary(mod1)$tTable["Locomotor_group", "p-value"], 3)), "\n")
## Model 1 with age: β = 0.013, t = 1.58, p = 0.119
cat("Model 1 without age:", 
    paste0("β = ", round(summary(mod1_no_age)$tTable["Locomotor_group", "Value"], 3),
           ", t = ", round(summary(mod1_no_age)$tTable["Locomotor_group", "t-value"], 2),
           ", p = ", round(summary(mod1_no_age)$tTable["Locomotor_group", "p-value"], 3)), "\n")
## Model 1 without age: β = 0.013, t = 1.64, p = 0.106
# add one-tailedp
one_tailed_p <- summary(mod1_no_age)$tTable["Locomotor_group", "p-value"] / 2
cat("one-tailed p-value:", round(one_tailed_p, 3), "\n\n")
## one-tailed p-value: 0.053
par(mfrow = c(1, 2))
plot_model(mod1, title = "Model 1 with age", show.values = TRUE)

plot_model(mod1_no_age, title = "Model 1 without age", show.values = TRUE)

4. Hypothesis 2: Effect of Navigation Experience

4.1 Preregistered Model 3: Navigators Only

data_included_navOnly <- subset(data_included, Locomotor_group == 1)
mod3 <- lme(propLookDoorCentered ~ SubjectInfo.navExp_scaled + SubjectInfo.Age_z + Texture_side + Trial_num, 
            random = ~1|SubjectInfo.subjID, 
            data = data_included_navOnly, 
            na.action = na.exclude)
summary(mod3)
## Linear mixed-effects model fit by REML
##   Data: data_included_navOnly 
##        AIC      BIC    logLik
##   137.6901 168.0479 -61.84505
## 
## Random effects:
##  Formula: ~1 | SubjectInfo.subjID
##          (Intercept) Residual
## StdDev: 8.575235e-06 0.263822
## 
## Fixed effects:  propLookDoorCentered ~ SubjectInfo.navExp_scaled + SubjectInfo.Age_z +      Texture_side + Trial_num 
##                                 Value  Std.Error  DF   t-value p-value
## (Intercept)                0.04398068 0.02663262 532  1.651384  0.0993
## SubjectInfo.navExp_scaled  0.01070637 0.04952899  33  0.216164  0.8302
## SubjectInfo.Age_z         -0.00787903 0.01333937  33 -0.590660  0.5588
## Texture_side               0.09131804 0.01107320 532  8.246765  0.0000
## Trial_num                 -0.00829637 0.03601516 532 -0.230358  0.8179
##  Correlation: 
##                           (Intr) SbI.E_ SbI.A_ Txtr_s
## SubjectInfo.navExp_scaled -0.537                     
## SubjectInfo.Age_z         -0.253 -0.065              
## Texture_side               0.043  0.001  0.000       
## Trial_num                 -0.678  0.008 -0.003 -0.064
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.41416031 -0.62330218  0.01651915  0.68639934  2.11700037 
## 
## Number of Observations: 570
## Number of Groups: 36
plot_model(mod3, show.values = TRUE)

4.2 Model 3 Subject-Level Sanity Check

data_trialavg_navOnly <- subset(data_trialavg, Locomotor_group == 1)
mod3_subject <- lm(PropLookDoor - 0.5 ~ SubjectInfo.navExp_scaled + SubjectInfo.Age_z, 
                   data = data_trialavg_navOnly)
summary(mod3_subject)
## 
## Call:
## lm(formula = PropLookDoor - 0.5 ~ SubjectInfo.navExp_scaled + 
##     SubjectInfo.Age_z, data = data_trialavg_navOnly)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.128548 -0.038849 -0.001745  0.028161  0.123316 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)  
## (Intercept)                0.040238   0.016896   2.382   0.0232 *
## SubjectInfo.navExp_scaled  0.011462   0.042675   0.269   0.7899  
## SubjectInfo.Age_z         -0.007833   0.011495  -0.681   0.5003  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05699 on 33 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.01538,    Adjusted R-squared:  -0.04429 
## F-statistic: 0.2578 on 2 and 33 DF,  p-value: 0.7743
one_tailed_p_navexp_mod3 <- summary(mod3_subject)$coefficients["SubjectInfo.navExp_scaled", "Pr(>|t|)"] / 2
one_tailed_p_age_mod3 <- summary(mod3_subject)$coefficients["SubjectInfo.Age_z", "Pr(>|t|)"] / 2

cat("SubjectInfo.navExp_scaled:", round(one_tailed_p_navexp_mod3, 3), "\n")
## SubjectInfo.navExp_scaled: 0.395
cat("SubjectInfo.Age_z:", round(one_tailed_p_age_mod3, 3), "\n\n")
## SubjectInfo.Age_z: 0.25
plot_model(mod3_subject, title = "mod3 subject-level (navigators only)", show.values = TRUE)

5. Hypothesis 3: Age Effects and Combined Model

5.1 Age and Looking Behavior Visualization

# 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)

5.2 Age Correlation Analysis

# 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

5.3 Preregistered Model 2: Full Sample Analysis

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)

5.4 Model 2 Subject-Level Sanity Checks

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",])

6. Additional Factors

6.1 Texture Side Effects

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

6.2 Competing Predictors Analysis

# 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