Grooming Behavior
Jesi
2025-07-05
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:
ggplot(data, aes(x = trial, y = groomALL, 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 = "Grooming Bouts Across Trials by Treatment - Raw Data", y = "groomALL", x = "Trial", color="Hot?")
# Experiment 1
p1<-ggplot(exp1, aes(x = factor(trial), y = groomALL, 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 = "groomALL",
x = "Trial",
color = "Hot?") +
theme_minimal()
# Experiment 2
p2<-ggplot(exp2, aes(x = factor(trial), y = groomALL, 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")
#no transformations - count data
modelhot <- glmer( groomALL ~ hot + (1 | subject), data = data, family = poisson)
modelnull<-glmer( groomALL ~ 1 + (1 | subject), data = data, family = poisson)
anova(modelhot,modelnull)## Data: data
## Models:
## modelnull: groomALL ~ 1 + (1 | subject)
## modelhot: groomALL ~ hot + (1 | subject)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## modelnull 2 1132.5 1138.3 -564.24 1128.5
## modelhot 3 1119.3 1128.0 -556.63 1113.3 15.222 1 9.557e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(modelhot)## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## hot 1 15.103 15.103 15.103#significant effect of hot on groomALL
#Model with trial and experiment (model with all interactions does not converge - too complex for this dataset)
modeli<-glmer( groomALL ~ hot + trial * exp + (1 | subject), data=data, family=poisson)
modelplus <- glmer( groomALL ~ hot + trial + exp + (1 | subject), data=data, family=poisson)
modelexp<- glmer( groomALL ~ hot + trial + (1 | subject), data=data, family=poisson)
modeltrial<-glmer( groomALL ~ hot + exp + (1 | subject), data=data, family=poisson)
anova(modeli, modelplus, modelhot, modelexp, modeltrial)## Data: data
## Models:
## modelhot: groomALL ~ hot + (1 | subject)
## modeltrial: groomALL ~ hot + exp + (1 | subject)
## modelexp: groomALL ~ hot + trial + (1 | subject)
## modelplus: groomALL ~ hot + trial + exp + (1 | subject)
## modeli: groomALL ~ hot + trial * exp + (1 | subject)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## modelhot 3 1119.3 1128.0 -556.63 1113.3
## modeltrial 4 1113.6 1125.3 -552.79 1105.6 7.6684 1 0.0056196 **
## modelexp 5 1072.5 1087.2 -531.27 1062.5 43.0526 1 5.329e-11 ***
## modelplus 6 1066.9 1084.4 -527.43 1054.9 7.6684 1 0.0056196 **
## modeli 8 1052.7 1076.2 -518.36 1036.7 18.1379 2 0.0001152 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#interaction model is best fit. However, my hypothesis was simple: subjects in the hot treatment was groom more often than those in the control treatment. should i stick to the simple model? or do i need to report all of these interactions if it affects interpretability?
model <- glmer( groomALL ~ hot + trial * exp + (1 | subject), data=data, family=poisson)
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':
##
## recodeAnova(model, type="III")## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: groomALL
## Chisq Df Pr(>Chisq)
## (Intercept) 3.0991 1 0.0783356 .
## hot 19.4448 1 1.035e-05 ***
## trial 36.4974 2 1.188e-08 ***
## exp 4.2887 1 0.0383671 *
## trial:exp 18.0471 2 0.0001205 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#treatment effect persists in the interaction model. Significant interaction for trial*exp, significant effect of trial and exp
# 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.98351, p-value = 0.0952# 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 okay
Looking across trials
Trial 1 treatment model vs null
library(lmerTest)
#including exp does not improve model
model <- glmer(groomALL ~ hot +
(1 | subject),
data = T1, family=poisson)
nullmodel <- glmer(groomALL ~ 1 +
(1 | subject),
data = T1, family=poisson)
anova(model, nullmodel)## Data: T1
## Models:
## nullmodel: groomALL ~ 1 + (1 | subject)
## model: groomALL ~ hot + (1 | subject)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## nullmodel 2 290.01 293.67 -143.01 286.01
## model 3 276.85 282.34 -135.43 270.85 15.16 1 9.876e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(model, type="III")## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: groomALL
## Chisq Df Pr(>Chisq)
## (Intercept) 6.1568 1 0.01309 *
## hot 15.2094 1 9.623e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#significant improvement on the null model, significant effect of treatment in trial 1
Trial 2 treatment model vs null
#including exp does not improve model
model <- glmer( groomALL ~ hot +
(1 | subject),
data = T2, family=poisson)
nullmodel <- glmer(groomALL ~ 1 +
(1 | subject),
data = T2, family=poisson)
anova(model, nullmodel)## Data: T2
## Models:
## nullmodel: groomALL ~ 1 + (1 | subject)
## model: groomALL ~ hot + (1 | subject)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## nullmodel 2 358.26 361.92 -177.13 354.26
## model 3 360.25 365.74 -177.12 354.25 0.0106 1 0.9181summary(model)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: groomALL ~ hot + (1 | subject)
## Data: T2
##
## AIC BIC logLik -2*log(L) df.resid
## 360.2 365.7 -177.1 354.2 43
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8702 -1.3430 -0.0594 0.9798 3.0063
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 1.058 1.028
## Number of obs: 46, groups: subject, 23
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.64718 0.23507 7.007 2.43e-12 ***
## hoty 0.01058 0.10163 0.104 0.917
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## hoty -0.217#does not improve upon null model
Trial 3 treatment model vs null
#including exp does not improve model
model <- glmer( groomALL ~ hot +
(1 | subject),
data = T3, family=poisson)
nullmodel <- glmer( groomALL ~ 1 +
(1 | subject),
data = T3, family=poisson)
anova(model, nullmodel)## Data: T3
## Models:
## nullmodel: groomALL ~ 1 + (1 | subject)
## model: groomALL ~ hot + (1 | subject)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## nullmodel 2 412.72 416.37 -204.36 408.72
## model 3 405.02 410.51 -199.51 399.02 9.6951 1 0.001848 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(model, type="II")## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: groomALL
## Chisq Df Pr(>Chisq)
## hot 9.8027 1 0.001743 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#significant improvement on null model, hot grooms significantly more in trial 3
Separated by Experiment
Experiment 1 treatment model vs null
#interaction does not improve model in exp 1
model <- glmer(groomALL ~ hot +
(1 | subject),
data = exp1, family=poisson)
nullmodel <- glmer(groomALL ~ 1 +
(1 | subject),
data = exp1, family=poisson)
anova(model, nullmodel)## Data: exp1
## Models:
## nullmodel: groomALL ~ 1 + (1 | subject)
## model: groomALL ~ hot + (1 | subject)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## nullmodel 2 509.27 513.74 -252.64 505.27
## model 3 511.26 517.96 -252.63 505.26 0.0135 1 0.9075Anova(model, type="II")## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: groomALL
## Chisq Df Pr(>Chisq)
## hot 0.0135 1 0.9074#treatment not significant in experiment 1
Experiment 2
# Experiment 2 treatment model vs null
model <- glmer( groomALL ~ hot +
(1 | subject),
data = exp2, family=poisson)
nullmodel <- glmer( groomALL ~ 1 +
(1 | subject),
data = exp2, family=poisson)
anova(model, nullmodel)## Data: exp2
## Models:
## nullmodel: groomALL ~ 1 + (1 | subject)
## model: groomALL ~ hot + (1 | subject)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## nullmodel 2 471.04 475.50 -233.52 467.04
## model 3 469.45 476.15 -231.72 463.45 3.5854 1 0.05829 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(model, type="II")## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: groomALL
## Chisq Df Pr(>Chisq)
## hot 3.7896 1 0.05157 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Nearly significant effect of treatment. hot treatment=more grooming, not significant.
#combined experiments = larger sample size = more power to make a conclusion. However, patterns don’t match between experiments.
Extras——————————————————————————————
Side effect?
modelside<- glmer(groomALL~ side + (1|subject), data=data, family=poisson)
summary(modelside)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: groomALL ~ side + (1 | subject)
## Data: data
##
## AIC BIC logLik -2*log(L) df.resid
## 1134.0 1142.8 -564.0 1128.0 135
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6134 -1.4693 -0.5606 1.0459 6.7604
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 1.332 1.154
## Number of obs: 138, groups: subject, 23
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.3535 0.3442 3.933 8.4e-05 ***
## sider 0.3353 0.4934 0.679 0.497
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## sider -0.694#no effect from which side of experimental tank was nearest the window
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
#output warning based on size of quadratic variables, rescale:
data$hour2_scaled <- (data$hour_decimal^2) / 100 # or /1000
# Then refit
model_quad_scaled <- glmer(groomALL ~ hot + hour2_scaled + (1 | subject),
data = data, family = poisson)
Anova(model_quad_scaled, type="III")## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: groomALL
## Chisq Df Pr(>Chisq)
## (Intercept) 4.7266 1 0.02970 *
## hot 5.6493 1 0.01746 *
## hour2_scaled 37.5161 1 9.066e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#significant quadratic effect of time of day tested. More grooming occurs during early morning and evening hours and is reduced during midday.
#significant effect of treatment persists when time of day is taken into account.
#hot*time interaction is not significant, results in higher AIC.
#Plot
data$hot <- as.factor(data$hot)
# Center time
data$hour_c <- scale(data$hour_decimal, scale = FALSE)
data$hour2_c <- data$hour_c^2
model_interact <- glmer(groomALL ~ hour_c + hour2_c + hot + (1 | subject),
data = data, family = poisson)
# Create time sequence
hour_seq <- seq(min(data$hour_decimal), max(data$hour_decimal), length.out = 100)
hour_c_seq <- hour_seq - mean(data$hour_decimal)
hour2_c_seq <- hour_c_seq^2
# Create grid for both levels of hot
new_data <- expand.grid(
hour_c = hour_c_seq,
hot = levels(data$hot)
)
# Add hour2 and original hour for plotting
new_data$hour2_c <- new_data$hour_c^2
new_data$hour_decimal <- new_data$hour_c + mean(data$hour_decimal)
new_data$predicted <- predict(model_interact,
newdata = new_data,
type = "response",
re.form = NA)
ggplot(new_data, aes(x = hour_decimal, y = predicted, color = hot)) +
geom_line(linewidth = 1.2) +
labs(title = "Predicted Grooming Across Time of Day by Treatment",
x = "Hour of Day",
y = "Predicted Grooming Count",
color = "Treatment") +
theme_minimal()
Tank effect?
modeltank<- glmer(groomALL~ as.factor(tank) + (1|subject), data=data, family=poisson)
summary(modeltank)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: groomALL ~ as.factor(tank) + (1 | subject)
## Data: data
##
## AIC BIC logLik -2*log(L) df.resid
## 1121.2 1141.6 -553.6 1107.2 131
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6002 -1.4663 -0.5435 0.7652 6.7020
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 1.26 1.122
## Number of obs: 138, groups: subject, 23
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.99192 0.34335 5.801 6.57e-09 ***
## as.factor(tank)2 -0.31276 0.20353 -1.537 0.1244
## as.factor(tank)3 -0.04776 0.20216 -0.236 0.8132
## as.factor(tank)4 0.46409 0.26185 1.772 0.0763 .
## as.factor(tank)5 -0.97868 0.50559 -1.936 0.0529 .
## as.factor(tank)6 -0.26671 0.15540 -1.716 0.0861 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) as.()2 as.()3 as.()4 as.()5
## as.fctr(t)2 -0.312
## as.fctr(t)3 -0.317 0.585
## as.fctr(t)4 -0.298 0.447 0.450
## as.fctr(t)5 -0.677 0.216 0.218 0.200
## as.fctr(t)6 -0.373 0.762 0.766 0.585 0.256#no significant effect of home tank on grooming