Analysis - personality of fostered frogs
Olivia Feagles, Yusan Yang
2026-04-03
Preliminaries
Load libraries
library(knitr) #for R Markdown
library(MASS) #for negative binomial glm (glm.nb)
library(lme4) #for mixed models
library(emmeans) #for posthoc
library(car) #for Anova
library(survival) #for survival analysis
library(coxme)
library(rptR) #for repeatability analysis
library(MuMIn) #for model selection (dredge)
library(ggfortify) #for plotting survival analysis
library(ggsignif) #for labeling significance in ggplots
library(GGally) #for correlation matrix
library(DescTools) #for Kendall's W test
library(ggpubr)
library(tidyverse) #for data processing, put last to avoid function maskingData processing
data_raw <- read.csv("Jacinta_prelim_040326.csv",
fileEncoding="UTF-8-BOM", # to avoid reading errors
na.strings = c("","N/A","#VALUE!","unknown")) # set things that represents NA
dataJR <- data_raw %>%
mutate_at(c("ID","Sex","Morph","Treatment"), as.factor) %>%
mutate(Treatment = fct_relevel(Treatment, "pure-bio","cross-foster","solitary"))
#generate data averaged by ID
data.avg <- dataJR %>% group_by(ID, Sex, Morph, Treatment) %>%
summarize(across(c(Emerge.sec, Contact.sec, Grids.prop), mean, na.rm = TRUE),
.groups = "drop"
)Figures
figcolorJR = c("darkolivegreen3", "dodgerblue", "darkorange")
figcolorJR2 = c("darkolivegreen", "dodgerblue4", "darkorange4")Fig_emerge <-
ggplot(data.avg,
aes(x= Treatment, y= Emerge.sec, fill = Treatment,color = Treatment)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(width = 0.1,
alpha = 0.3,
size = 3) +
scale_fill_manual(values = figcolorJR) +
scale_color_manual(values = figcolorJR2) +
theme_classic(base_family = "sans", base_size = 14) +
theme(legend.position = "none") +
labs(x = "Rearing Treatment",
y = "Emergence Latency (sec)")
Fig_contact <-
ggplot(data.avg,
aes(x= Treatment, y= Contact.sec, fill = Treatment, color = Treatment)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(width = 0.1,
alpha = 0.3,
size = 3) +
scale_fill_manual(values = figcolorJR) +
scale_color_manual(values = figcolorJR2) +
theme_classic(base_family = "sans", base_size = 14) +
theme(legend.position = "none") +
labs(x = "Rearing Treatment",
y = "Contact Time (sec)")
Fig_explore <- ggplot(data.avg, aes(x= Treatment, y= Grids.prop,
fill = Treatment,
color = Treatment)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(width = 0.1,
alpha = 0.3,
size = 3) +
scale_fill_manual(values = figcolorJR) +
scale_color_manual(values = figcolorJR2) +
theme_classic(base_family = "sans", base_size = 14) +
theme(legend.title = element_text(size = 12)) +
labs(x = "Rearing Treatment",
y = "Proportion of Area Explored")ggpubr::ggarrange(Fig_emerge, Fig_contact, Fig_explore,
ncol = 3, nrow = 1,
common.legend = TRUE, legend = "right",
labels = c("A", "B", "C"))
Analysis: Trt and Sex
model_emerge <- lmer(Emerge.sec ~ Treatment + Sex + (1 | ID), data = dataJR)
model_contact <- lmer(Contact.sec ~ Treatment + Sex + (1 | ID), data = dataJR)
model_grids <- lmer(Grids.prop ~ Treatment + Sex + (1 | ID), data = dataJR)Latency to emerge
summary(model_emerge)## Linear mixed model fit by REML ['lmerMod']
## Formula: Emerge.sec ~ Treatment + Sex + (1 | ID)
## Data: dataJR
##
## REML criterion at convergence: 377.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.1442 -0.7172 -0.2419 0.2884 2.5235
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0 0.0
## Residual 283837 532.8
## Number of obs: 28, groups: ID, 14
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 607.1 172.9 3.512
## Treatmentcross-foster 729.8 291.0 2.508
## Treatmentsolitary 359.5 236.7 1.519
## SexM -420.2 228.7 -1.838
##
## Correlation of Fixed Effects:
## (Intr) Trtmn- Trtmnt
## Trtmntcrss- -0.149
## Trtmntsltry -0.456 0.325
## SexM -0.567 -0.449 -0.069
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')Anova(model_emerge)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Emerge.sec
## Chisq Df Pr(>Chisq)
## Treatment 6.8436 2 0.03265 *
## Sex 3.3767 1 0.06612 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1posthoc:
emmeans(model_emerge, ~ Treatment) %>% pairs()## contrast estimate SE df t.ratio p.value
## (pure-bio) - (cross-foster) -730 291 10 -2.508 0.0731
## (pure-bio) - solitary -360 237 10 -1.519 0.3234
## (cross-foster) - solitary 370 310 10 1.196 0.4819
##
## Results are averaged over the levels of: Sex
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimatesLatency to contact
summary(model_contact)## Linear mixed model fit by REML ['lmerMod']
## Formula: Contact.sec ~ Treatment + Sex + (1 | ID)
## Data: dataJR
##
## REML criterion at convergence: 234.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.1384 -0.7232 -0.2520 0.4246 2.5568
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.0 0.00
## Residual 717.8 26.79
## Number of obs: 28, groups: ID, 14
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 30.500 8.692 3.509
## Treatmentcross-foster -6.500 14.631 -0.444
## Treatmentsolitary -22.500 11.902 -1.890
## SexM -3.500 11.499 -0.304
##
## Correlation of Fixed Effects:
## (Intr) Trtmn- Trtmnt
## Trtmntcrss- -0.149
## Trtmntsltry -0.456 0.325
## SexM -0.567 -0.449 -0.069
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')Anova(model_contact)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Contact.sec
## Chisq Df Pr(>Chisq)
## Treatment 3.6062 2 0.1648
## Sex 0.0926 1 0.7608Proportion explored
summary(model_grids)## Linear mixed model fit by REML ['lmerMod']
## Formula: Grids.prop ~ Treatment + Sex + (1 | ID)
## Data: dataJR
##
## REML criterion at convergence: -31
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6167 -0.9104 0.1679 0.4327 1.8624
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.00000 0.0000
## Residual 0.01142 0.1069
## Number of obs: 28, groups: ID, 14
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.21765 0.03467 6.278
## Treatmentcross-foster 0.09345 0.05836 1.601
## Treatmentsolitary -0.06731 0.04747 -1.418
## SexM 0.01529 0.04586 0.333
##
## Correlation of Fixed Effects:
## (Intr) Trtmn- Trtmnt
## Trtmntcrss- -0.149
## Trtmntsltry -0.456 0.325
## SexM -0.567 -0.449 -0.069
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')Anova(model_grids)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Grids.prop
## Chisq Df Pr(>Chisq)
## Treatment 6.7685 2 0.0339 *
## Sex 0.1111 1 0.7389
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1posthoc:
emmeans(model_grids, ~ Treatment) %>% pairs()## contrast estimate SE df t.ratio p.value
## (pure-bio) - (cross-foster) -0.0934 0.0584 10 -1.601 0.2895
## (pure-bio) - solitary 0.0673 0.0475 10 1.418 0.3688
## (cross-foster) - solitary 0.1608 0.0621 10 2.589 0.0641
##
## Results are averaged over the levels of: Sex
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimatesAnalysis: repeatability
Preview - scatter plots
data_wide <- dataJR %>%
pivot_wider(
id_cols = c(ID, Sex, Morph, Treatment),
names_from = Replicate,
values_from = c(Emerge.sec, Contact.sec, Grids.prop)
)Cor_Emerge <-
ggplot(data_wide, aes(x = Emerge.sec_1, y = Emerge.sec_2, color = Treatment)) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "gray") +
geom_point(size = 3, alpha = 0.7) +
scale_color_manual(values = figcolorJR) +
theme_classic(base_size = 15) +
coord_fixed(ratio = 1, xlim = c(0, max(dataJR$Emerge.sec)), ylim = c(0, max(dataJR$Emerge.sec))) +
labs(title = "Emergence Latency (sec)", x = "Trial 1", y = "Trial 2")
Cor_Contact <-
ggplot(data_wide, aes(x = Contact.sec_1, y = Contact.sec_2, color = Treatment)) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "gray") +
geom_point(size = 3, alpha = 0.7) +
scale_color_manual(values = figcolorJR) +
theme_classic(base_size = 15) +
coord_fixed(ratio = 1, xlim = c(0, max(dataJR$Contact.sec)), ylim = c(0, max(dataJR$Contact.sec))) +
labs(title = "Contact Latency (sec)", x = "Trial 1", y = "Trial 2")
Cor_Grids <-
ggplot(data_wide, aes(x = Grids.prop_1, y = Grids.prop_2, color = Treatment)) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "gray") +
geom_point(size = 3, alpha = 0.7) +
scale_color_manual(values = figcolorJR) +
theme_classic(base_size = 15) +
coord_fixed(ratio = 1, xlim = c(0, max(dataJR$Grids.prop)), ylim = c(0, max(dataJR$Grids.prop)))+
labs(title = "Proportion Explored", x = "Trial 1", y = "Trial 2")ggarrange(Cor_Emerge, Cor_Contact, Cor_Grids,
ncol = 3, nrow = 1,
common.legend = TRUE,
legend = "right",
labels = "AUTO")
LMM
might be problematic because only 2 measurement per ID, and lots of singular models (i.e. no random effect)
rep_trial_Emerge_adj <- rpt(Emerge.sec ~ Treatment + (1 | ID),
grname = "ID",
data = dataJR,
datatype = "Gaussian",
nboot = 1000,
npermut = 1000)## Bootstrap Progress:
## Permutation Progress for ID :rep_trial_Emerge_adj##
##
## Repeatability estimation using the lmm method
##
## Repeatability for ID
## R = 0.125
## SE = 0.163
## CI = [0, 0.514]
## P = 0.411 [LRT]
## 0.261 [Permutation]rep_trial_Contact_adj <- rpt(Contact.sec ~ Treatment + (1 | ID),
grname = "ID",
data = dataJR,
datatype = "Gaussian",
nboot = 1000,
npermut = 1000)## Bootstrap Progress:
## Permutation Progress for ID :rep_trial_Contact_adj##
##
## Repeatability estimation using the lmm method
##
## Repeatability for ID
## R = 0
## SE = 0.121
## CI = [0, 0.401]
## P = 1 [LRT]
## 1 [Permutation]rep_trial_Grids_adj <- rpt(Grids.prop ~ Treatment + (1 | ID),
grname = "ID",
data = dataJR,
datatype = "Gaussian",
nboot = 1000,
npermut = 1000)## Bootstrap Progress:
## Permutation Progress for ID :rep_trial_Grids_adj##
##
## Repeatability estimation using the lmm method
##
## Repeatability for ID
## R = 0
## SE = 0.13
## CI = [0, 0.437]
## P = 1 [LRT]
## 1 [Permutation]Non-parametric
data_wide_Emerge <- dataJR %>%
select(ID, Replicate, Emerge.sec) %>%
pivot_wider(names_from = Replicate, values_from = Emerge.sec) %>%
select(-ID) %>%
as.matrix()
KendallW(data_wide_Emerge , correct = TRUE, test = TRUE)##
## Kendall's coefficient of concordance Wt
##
## data: data_wide_Emerge
## Kendall chi-squared = 24.831, df = 20, subjects = 21, raters = 2,
## p-value = 0.208
## alternative hypothesis: Wt is greater 0
## sample estimates:
## Wt
## 0.6207792data_wide_Contact <- dataJR %>%
select(ID, Replicate, Contact.sec) %>%
pivot_wider(names_from = Replicate, values_from = Contact.sec) %>%
select(-ID) %>%
as.matrix()
KendallW(data_wide_Contact , correct = TRUE, test = TRUE)##
## Kendall's coefficient of concordance Wt
##
## data: data_wide_Contact
## Kendall chi-squared = 21.071, df = 20, subjects = 21, raters = 2,
## p-value = 0.3929
## alternative hypothesis: Wt is greater 0
## sample estimates:
## Wt
## 0.5267826data_wide_Grids <- dataJR %>%
select(ID, Replicate, Grids.prop) %>%
pivot_wider(names_from = Replicate, values_from = Grids.prop) %>%
select(-ID) %>%
as.matrix()
KendallW(data_wide_Grids , correct = TRUE, test = TRUE)##
## Kendall's coefficient of concordance Wt
##
## data: data_wide_Grids
## Kendall chi-squared = 24.016, df = 20, subjects = 21, raters = 2,
## p-value = 0.2417
## alternative hypothesis: Wt is greater 0
## sample estimates:
## Wt
## 0.6003899