Pause
Jesi
2025-07-13
data<-Behavioral.Trials.All.Data#Separate trials
T1<- data[data$trial=="1",]
T2<- data[data$trial=="2",]
T3<- data[data$trial=="3",]
data23<- rbind(T2,T3)
#separate experiments
exp1<-(data[data$exp=="1",])
exp2<-(data[data$exp=="2",])
#separate groups
hot<-(data[data$hot=="y",])
control<-(data[data$hot=="n",])# Reorder subjects based on the 'hot' group
data <- data %>%
mutate(subject = factor(subject, levels = unique(subject[order(hot)])))Raw data plot:
#pauseCPA is time spent paused for at least 3 seconds during avoidance trial
ggplot(data, aes(x = trial, y = pauseCPA, color = hot)) +
stat_summary(fun = mean, geom = "line", aes(group = hot)) +
stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) +
stat_summary(fun = mean, geom = "point", size = 3) +
scale_color_manual(values = c("n" = "grey50", "y" = "red")) +
labs(title = "Pause Across Trials by Treatment - Raw Data", y = "pauseCPA", x = "Trial", color="Hot?")
# Experiment 1
p1<-ggplot(exp1, aes(x = factor(trial), y = pauseCPA, color = hot, group = hot)) +
geom_jitter(width = 0.2, alpha = 0.5, size = 1.8) + # raw data
stat_summary(fun = mean, geom = "line", aes(group = hot), linewidth = 1) + # mean line
stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) + # error bars
stat_summary(fun = mean, geom = "point", size = 3, shape = 18) + # mean points
scale_color_manual(values = c("n" = "grey50", "y" = "red")) +
labs(title = "Exp 1 - Data + Means and SE" ,
y = "pauseCPA",
x = "Trial",
color = "Hot?") +
theme_minimal()
# Experiment 2
p2<-ggplot(exp2, aes(x = factor(trial), y = pauseCPA, color = hot, group = hot)) +
geom_jitter(width = 0.2, alpha = 0.5, size = 1.8) +
stat_summary(fun = mean, geom = "line", aes(group = hot), linewidth = 1) +
stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) +
stat_summary(fun = mean, geom = "point", size = 3, shape = 18) +
scale_color_manual(values = c("n" = "grey50", "y" = "red")) +
labs(title = "Exp 2 - Data + Means and SE",
y = "",
x = "Trial",
color = "Hot?") +
theme_minimal()
library(patchwork)## Warning: package 'patchwork' was built under R version 4.3.3# Combine plots and collect shared legend
(p1 + p2) +
plot_layout(ncol = 2, guides = "collect") &
theme(legend.position = "bottom")
#weird pattern
library(car)## Warning: package 'car' was built under R version 4.3.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.3.3##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode# Step 1: Shift data to ensure all values are positive
min_val <- min(data$pauseCPA, na.rm = TRUE)
data$pauseCPA_shift <- data$pauseCPA + abs(min_val) + 1
# Step 2: Estimate the Box-Cox transformation lambda
lambda_result <- powerTransform(pauseCPA_shift ~ 1, data = data)
summary(lambda_result)## bcPower Transformation to Normality
## Est Power Rounded Pwr Wald Lwr Bnd Wald Upr Bnd
## Y1 0.7637 0.76 0.5964 0.9309
##
## Likelihood ratio test that transformation parameter is equal to 0
## (log transformation)
## LRT df pval
## LR test, lambda = (0) 139.5154 1 < 2.22e-16
##
## Likelihood ratio test that no transformation is needed
## LRT df pval
## LR test, lambda = (1) 6.688076 1 0.009706#log transformation would be best
# Apply the estimated lambda to transform the data
data$pauseCPA_trans <- log(data$pauseCPA_shift)
hist(data$pauseCPA)
hist(data$pauseCPA_trans)
#untranformed looks more normal
modelhot <- lmer( pauseCPA ~ hot + (1 | subject), data = data)
anova(modelhot)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## hot 82491 82491 1 114 7.2019 0.008367 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#significant effect of hot on pauseCPA. Hot treatment subjects significantly less likely to spend time not moving. Escape behavior?
#Model with trial and experiment
#tested full interaction models as well, did not improve fit
model<-lmer( pauseCPA ~ hot + exp + (1 | subject), data = data)
modeli <- lmer( pauseCPA ~ hot * exp + (1 | subject), data = data)
modeltrial <- lmer( pauseCPA ~ hot + trial + (1 | subject), data = data)
anova(model, modeli, modeltrial, modelhot)## refitting model(s) with ML (instead of REML)## Data: data
## Models:
## modelhot: pauseCPA ~ hot + (1 | subject)
## model: pauseCPA ~ hot + exp + (1 | subject)
## modeli: pauseCPA ~ hot * exp + (1 | subject)
## modeltrial: pauseCPA ~ hot + trial + (1 | subject)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## modelhot 4 1752.2 1763.9 -872.10 1744.2
## model 5 1748.6 1763.2 -869.30 1738.6 5.6005 1 0.01796 *
## modeli 6 1749.7 1767.3 -868.87 1737.7 0.8632 1 0.35283
## modeltrial 6 1755.6 1773.2 -871.80 1743.6 0.0000 0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(model)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## hot 88668 88668 1 113 8.0562 0.005378 **
## exp 62067 62067 1 113 5.6393 0.019248 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#simple model with exp is best fit - interactions are not significant
#treatment and experiment significantly affect pause time
# Residuals vs Fitted
plot(model, which = 1) # same as plot(fitted(model), resid(model))
# QQ Plot for normality
qqnorm(resid(model))
qqline(resid(model), col = "red")
# Histogram
hist(resid(model), breaks = 20, main = "Histogram of Residuals", xlab = "Residuals")
shapiro.test(resid(model)) # Note: sensitive to sample size##
## Shapiro-Wilk normality test
##
## data: resid(model)
## W = 0.98601, p-value = 0.1738# Plot absolute residuals vs fitted values
plot(fitted(model), abs(resid(model)),
xlab = "Fitted values", ylab = "Absolute Residuals",
main = "Homoscedasticity Check")
abline(lm(abs(resid(model)) ~ fitted(model)), col = "red")
#residuals look good for raw variable
By experiment
modelexp1 <- lmer( pauseCPA ~ hot + (1 | subject), data = exp1)
anova(modelexp1)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## hot 223.02 223.02 1 21 0.0353 0.8527modelexp2 <- lmer( pauseCPA ~ hot + (1 | subject), data = exp2)
anova(modelexp2)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## hot 25359 25359 1 21 2.7611 0.1114#effect not significant in either experiemnt separately
Include pausePOST - all pausing (post-treatment + CPA test)
#pausePOST is time spent paused for at least 3 seconds during the 5 minutes post-treatment
data$pause<- data$pauseCPA + data$pausePOSTggplot(data, aes(x = trial, y = pause, color = hot)) +
stat_summary(fun = mean, geom = "line", aes(group = hot)) +
stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) +
stat_summary(fun = mean, geom = "point", size = 3) +
scale_color_manual(values = c("n" = "grey50", "y" = "red")) +
labs(title = "Pause Across Trials by Treatment - Raw Data", y = "pause", x = "Trial", color="Hot?")
# Experiment 1
p1<-ggplot(exp1, aes(x = factor(trial), y = pause, color = hot, group = hot)) +
geom_jitter(width = 0.2, alpha = 0.5, size = 1.8) + # raw data
stat_summary(fun = mean, geom = "line", aes(group = hot), linewidth = 1) + # mean line
stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) + # error bars
stat_summary(fun = mean, geom = "point", size = 3, shape = 18) + # mean points
scale_color_manual(values = c("n" = "grey50", "y" = "red")) +
labs(title = "Exp 1 - Data + Means and SE" ,
y = "pause",
x = "Trial",
color = "Hot?") +
theme_minimal()
# Experiment 2
p2<-ggplot(exp2, aes(x = factor(trial), y = pause, color = hot, group = hot)) +
geom_jitter(width = 0.2, alpha = 0.5, size = 1.8) +
stat_summary(fun = mean, geom = "line", aes(group = hot), linewidth = 1) +
stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) +
stat_summary(fun = mean, geom = "point", size = 3, shape = 18) +
scale_color_manual(values = c("n" = "grey50", "y" = "red")) +
labs(title = "Exp 2 - Data + Means and SE",
y = "",
x = "Trial",
color = "Hot?") +
theme_minimal()
library(patchwork)
# Combine plots and collect shared legend
(p1 + p2) +
plot_layout(ncol = 2, guides = "collect") &
theme(legend.position = "bottom")
#same pattern when pausePOST is included, but more dramatic
library(car)
# Step 1: Shift data to ensure all values are positive
min_val <- min(data$pause, na.rm = TRUE)
data$pause_shift <- data$pause + abs(min_val) + 1
# Step 2: Estimate the Box-Cox transformation lambda
lambda_result <- powerTransform(pause_shift ~ 1, data = data)
summary(lambda_result)## bcPower Transformation to Normality
## Est Power Rounded Pwr Wald Lwr Bnd Wald Upr Bnd
## Y1 0.8935 1 0.7045 1.0825
##
## Likelihood ratio test that transformation parameter is equal to 0
## (log transformation)
## LRT df pval
## LR test, lambda = (0) 175.5975 1 < 2.22e-16
##
## Likelihood ratio test that no transformation is needed
## LRT df pval
## LR test, lambda = (1) 1.145237 1 0.28455#log transformation would be best
# Apply the estimated lambda to transform the data
data$pause_trans <- log(data$pause_shift)
hist(data$pause)
hist(data$pause_trans)
#untranformed looks more normal
modelhot <- lmer( pause ~ hot + (1 | subject), data = data)
anova(modelhot)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## hot 364603 364603 1 114 21.687 8.746e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#significant effect of hot on pause. Hot treatment subjects significantly less likely to spend time not moving. Escape behavior?
#Model with trial and experiment
#I tested full models with interactions as well but did not improve fit.
model<-lmer( pause ~ hot + exp + (1 | subject), data = data)
modeli <- lmer( pause ~ hot * exp + (1 | subject), data = data)
modeltrial <- lmer( pause ~ hot + trial + (1 | subject), data = data)
anova(model, modeli, modeltrial, modelhot)## refitting model(s) with ML (instead of REML)## Data: data
## Models:
## modelhot: pause ~ hot + (1 | subject)
## model: pause ~ hot + exp + (1 | subject)
## modeli: pause ~ hot * exp + (1 | subject)
## modeltrial: pause ~ hot + trial + (1 | subject)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## modelhot 4 1810.8 1822.5 -901.40 1802.8
## model 5 1809.4 1824.0 -899.70 1799.4 3.3918 1 0.06552 .
## modeli 6 1810.5 1828.0 -899.24 1798.5 0.9327 1 0.33417
## modeltrial 6 1812.2 1829.8 -900.09 1800.2 0.0000 0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(model)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## hot 376399 376399 1 113 22.8563 5.301e-06 ***
## exp 55702 55702 1 113 3.3824 0.06852 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#simple model is best fit - interactions are not significant
# Residuals vs Fitted
plot(modelhot, which = 1) # same as plot(fitted(model), resid(model))
# QQ Plot for normality
qqnorm(resid(modelhot))
qqline(resid(modelhot), col = "red")
# Histogram
hist(resid(modelhot), breaks = 20, main = "Histogram of Residuals", xlab = "Residuals")
shapiro.test(resid(modelhot)) # Note: sensitive to sample size##
## Shapiro-Wilk normality test
##
## data: resid(modelhot)
## W = 0.98821, p-value = 0.2902# Plot absolute residuals vs fitted values
plot(fitted(modelhot), abs(resid(modelhot)),
xlab = "Fitted values", ylab = "Absolute Residuals",
main = "Homoscedasticity Check")
abline(lm(abs(resid(modelhot)) ~ fitted(modelhot)), col = "red")
#residuals look good for raw variable
By experiment
modelexp1 <- lmer( pause ~ hot + (1 | subject), data = exp1)
anova(modelexp1)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## hot 41.514 41.514 1 21 0.0054 0.9422modelexp2 <- lmer( pause ~ hot + (1 | subject), data = exp2)
anova(modelexp2)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## hot 54124 54124 1 21 4.5077 0.04579 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#effect is significant in experiment 2, but not exp 1.
# Fit the model with hot as fixed effect, subject as random
model_pause <- lmer(pause ~ hot + (1 | subject), data = data)
# Get estimated marginal means for hot vs. control
emm_pause <- emmeans(model_pause, ~ hot)
# View the summary
summary(emm_pause)## hot emmean SE df lower.CL upper.CL
## n 596 52.3 24.1 488 703
## y 493 52.3 24.1 385 601
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95Extras —————————————————————————————–
Side effect?
modelside<- lmer(pause~hot + side + (1|subject), data=data)
summary(modelside)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: pause ~ hot + side + (1 | subject)
## Data: data
##
## REML criterion at convergence: 1770.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5260 -0.6986 0.1036 0.7021 2.8400
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 49981 223.6
## Residual 16812 129.7
## Number of obs: 138, groups: subject, 23
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 687.55 67.23 22.18 10.226 7.37e-10 ***
## hoty -102.80 22.08 114.00 -4.657 8.75e-06 ***
## sider -192.18 95.90 21.00 -2.004 0.0581 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) hoty
## hoty -0.164
## sider -0.682 0.000#no significant effect of side, trend to locomote more if treatment was given on side away from window (but closer to overhead light)
Time effect
library(hms)## Warning: package 'hms' was built under R version 4.3.3library(lubridate)## Warning: package 'lubridate' was built under R version 4.3.3##
## Attaching package: 'lubridate'## The following object is masked from 'package:hms':
##
## hms## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union# Step 1: Add ":00" to make it HH:MM:SS
data$time_clean <- paste0(data$time, ":00")
# Step 2: Parse as hms time
data$time_parsed <- hms::as_hms(data$time_clean)
# Step 3: Extract decimal hour
data$hour_decimal <- hour(data$time_parsed) + minute(data$time_parsed)/60data$hour2 <- data$hour_decimal^2
model_time_quad <- lmer(pause ~ hour_decimal + hour2 + (1 | subject), data = data)
summary(model_time_quad)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: pause ~ hour_decimal + hour2 + (1 | subject)
## Data: data
##
## REML criterion at convergence: 1767.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.14047 -0.54267 -0.01535 0.69490 2.68194
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 45844 214.1
## Residual 15241 123.5
## Number of obs: 138, groups: subject, 23
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 539.220 364.249 121.655 1.480 0.141
## hour_decimal 38.031 57.475 119.009 0.662 0.509
## hour2 -2.869 2.220 119.596 -1.293 0.199
##
## Correlation of Fixed Effects:
## (Intr) hr_dcm
## hour_deciml -0.987
## hour2 0.973 -0.995#no significant effect of time of day on pause behavior