Statistical Analysis for MSc Project - Auditory Enrichment Effects on a Solitary Captive African Elephant- Mondula
Caitlin Murray
2024-08-02
Loading Libraries
# ggplot2 for data visualization
library(ggplot2)
# dplyr for data manipulation (part of tidyverse, but included separately for clarity)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union# tidyr for data tidying (part of tidyverse, but included separately for clarity)
library(tidyr)
# knitr for creating dynamic reports
library(knitr)
# lme4 for fitting linear and generalized linear mixed-effects models
library(lme4)## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack# multcomp for multiple comparison procedures
library(multcomp) ## Loading required package: mvtnorm## Loading required package: survival## Loading required package: TH.data## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select##
## Attaching package: 'TH.data'## The following object is masked from 'package:MASS':
##
## geyser# reshape2 for transforming data between wide and long formats
library(reshape2)##
## Attaching package: 'reshape2'## The following object is masked from 'package:tidyr':
##
## smiths# broom for tidying and organising model output
library(broom)BEHAVIOURAL ANALYSIS
Total Time Spent Per Behaviour Across Treatment Phases
Data Preparation
# Reading in the data set into dataframe
total_per_obs <- read.csv("~/Desktop/Mondy/Total per obs .csv")
# Viewing the data set
head(total_per_obs)Calculating Summary Statistics
# Calculate summary statistics
summary_stats <- total_per_obs %>%
summarise(across(Foraging:Play, list(mean = ~mean(.),
median = ~median(.),
sd = ~sd(.),
min = ~min(.),
max = ~max(.))))
# View the summary statistics
print(summary_stats)## Foraging_mean Foraging_median Foraging_sd Foraging_min Foraging_max
## 1 2344.18 2319.5 675.8739 1030 3600
## Locomotion_mean Locomotion_median Locomotion_sd Locomotion_min Locomotion_max
## 1 721.44 680 376.255 0 1777
## Attentive_mean Attentive_median Attentive_sd Attentive_min Attentive_max
## 1 184.48 93.5 267.3975 0 1236
## Inactive_mean Inactive_median Inactive_sd Inactive_min Inactive_max
## 1 209.32 67 289.3201 0 997
## Self.Maintenance_mean Self.Maintenance_median Self.Maintenance_sd
## 1 128.4 69.5 175.1053
## Self.Maintenance_min Self.Maintenance_max Vocalisation_mean
## 1 0 844 1.7
## Vocalisation_median Vocalisation_sd Vocalisation_min Vocalisation_max
## 1 0 6.145348 0 33
## Stereotyping_mean Stereotyping_median Stereotyping_sd Stereotyping_min
## 1 0.24 0 1.697056 0
## Stereotyping_max Play_mean Play_median Play_sd Play_min Play_max
## 1 12 8.58 0 36.12907 0 236# View the summary statistics in a regular table format
summary_stats_table <- summary_stats %>%
pivot_longer(cols = everything(),
names_to = c("Behaviour", ".value"),
names_sep = "_") %>%
arrange(Behaviour)
# Display the table
kable(summary_stats_table, caption = "Summary Statistics for Each Behavioural Category")| Behaviour | mean | median | sd | min | max |
|---|---|---|---|---|---|
| Attentive | 184.48 | 93.5 | 267.397452 | 0 | 1236 |
| Foraging | 2344.18 | 2319.5 | 675.873907 | 1030 | 3600 |
| Inactive | 209.32 | 67.0 | 289.320064 | 0 | 997 |
| Locomotion | 721.44 | 680.0 | 376.254970 | 0 | 1777 |
| Play | 8.58 | 0.0 | 36.129071 | 0 | 236 |
| Self.Maintenance | 128.40 | 69.5 | 175.105333 | 0 | 844 |
| Stereotyping | 0.24 | 0.0 | 1.697056 | 0 | 12 |
| Vocalisation | 1.70 | 0.0 | 6.145348 | 0 | 33 |
# Function to calculate summary statistics
summary_stats <- function(total_per_obs, treatment) {
total_per_obs %>%
filter(Treatment == treatment) %>%
summarise(
Foraging_mean = mean(Foraging, na.rm = TRUE),
Foraging_median = median(Foraging, na.rm = TRUE),
Foraging_sd = sd(Foraging, na.rm = TRUE),
Foraging_min = min(Foraging, na.rm = TRUE),
Foraging_max = max(Foraging, na.rm = TRUE),
Locomotion_mean = mean(Locomotion, na.rm = TRUE),
Locomotion_median = median(Locomotion, na.rm = TRUE),
Locomotion_sd = sd(Locomotion, na.rm = TRUE),
Locomotion_min = min(Locomotion, na.rm = TRUE),
Locomotion_max = max(Locomotion, na.rm = TRUE),
Attentive_mean = mean(Attentive, na.rm = TRUE),
Attentive_median = median(Attentive, na.rm = TRUE),
Attentive_sd = sd(Attentive, na.rm = TRUE),
Attentive_min = min(Attentive, na.rm = TRUE),
Attentive_max = max(Attentive, na.rm = TRUE),
Inactive_mean = mean(Inactive, na.rm = TRUE),
Inactive_median = median(Inactive, na.rm = TRUE),
Inactive_sd = sd(Inactive, na.rm = TRUE),
Inactive_min = min(Inactive, na.rm = TRUE),
Inactive_max = max(Inactive, na.rm = TRUE),
Self.Maintenance_mean = mean(Self.Maintenance, na.rm = TRUE),
Self.Maintenance_median = median(Self.Maintenance, na.rm = TRUE),
Self.Maintenance_sd = sd(Self.Maintenance, na.rm = TRUE),
Self.Maintenance_min = min(Self.Maintenance, na.rm = TRUE),
Self.Maintenance_max = max(Self.Maintenance, na.rm = TRUE),
Vocalisation_mean = mean(Vocalisation, na.rm = TRUE),
Vocalisation_median = median(Vocalisation, na.rm = TRUE),
Vocalisation_sd = sd(Vocalisation, na.rm = TRUE),
Vocalisation_min = min(Vocalisation, na.rm = TRUE),
Vocalisation_max = max(Vocalisation, na.rm = TRUE),
Stereotyping_mean = mean(Stereotyping, na.rm = TRUE),
Stereotyping_median = median(Stereotyping, na.rm = TRUE),
Stereotyping_sd = sd(Stereotyping, na.rm = TRUE),
Stereotyping_min = min(Stereotyping, na.rm = TRUE),
Stereotyping_max = max(Stereotyping, na.rm = TRUE),
Play_mean = mean(Play, na.rm = TRUE),
Play_median = median(Play, na.rm = TRUE),
Play_sd = sd(Play, na.rm = TRUE),
Play_min = min(Play, na.rm = TRUE),
Play_max = max(Play, na.rm = TRUE)
)
}
# Calculate summary statistics for each treatment
treatments <- unique(total_per_obs$Treatment)
summary_results <- bind_rows(lapply(treatments, function(treatment) {
summary_stats(total_per_obs, treatment) %>%
mutate(Treatment = treatment)
}))
# Display the summary statistics
print(summary_results)## Foraging_mean Foraging_median Foraging_sd Foraging_min Foraging_max
## 1 2519.7 2543.5 639.7190 1686 3467
## 2 2106.1 1993.5 730.9854 1072 3301
## 3 2689.1 2575.0 548.9294 2044 3600
## 4 2182.1 2300.0 603.0101 1353 3095
## 5 2223.9 2072.0 777.0875 1030 3141
## Locomotion_mean Locomotion_median Locomotion_sd Locomotion_min Locomotion_max
## 1 685.5 625.0 459.8198 25 1468
## 2 878.8 806.5 507.1396 111 1777
## 3 587.3 680.0 305.4199 0 929
## 4 680.9 619.5 200.6968 392 1071
## 5 774.7 720.0 338.6047 291 1413
## Attentive_mean Attentive_median Attentive_sd Attentive_min Attentive_max
## 1 65.4 38.5 95.84154 0 318
## 2 227.1 131.5 287.47694 0 905
## 3 44.1 10.5 63.69624 0 187
## 4 257.0 148.5 358.40883 0 1236
## 5 328.8 297.0 310.30086 0 821
## Inactive_mean Inactive_median Inactive_sd Inactive_min Inactive_max
## 1 201.8 124.0 279.7212 0 935
## 2 225.2 86.0 342.6176 0 997
## 3 106.4 53.0 142.3612 0 417
## 4 338.9 184.5 357.0174 12 925
## 5 174.3 0.0 283.2710 0 675
## Self.Maintenance_mean Self.Maintenance_median Self.Maintenance_sd
## 1 127.0 89.5 147.14166
## 2 122.6 62.0 142.50786
## 3 170.4 88.5 219.30354
## 4 137.3 52.5 252.54441
## 5 84.7 69.0 94.31631
## Self.Maintenance_min Self.Maintenance_max Vocalisation_mean
## 1 0 479 0.6
## 2 17 486 4.2
## 3 0 716 2.7
## 4 0 844 0.0
## 5 0 289 1.0
## Vocalisation_median Vocalisation_sd Vocalisation_min Vocalisation_max
## 1 0 1.897367 0 6
## 2 0 10.293472 0 33
## 3 0 8.538150 0 27
## 4 0 0.000000 0 0
## 5 0 3.162278 0 10
## Stereotyping_mean Stereotyping_median Stereotyping_sd Stereotyping_min
## 1 0.0 0 0.000000 0
## 2 0.0 0 0.000000 0
## 3 0.0 0 0.000000 0
## 4 0.0 0 0.000000 0
## 5 1.2 0 3.794733 0
## Stereotyping_max Play_mean Play_median Play_sd Play_min Play_max Treatment
## 1 0 0.0 0 0.00000 0 0 Pre
## 2 0 32.0 0 76.38499 0 236 CM
## 3 0 0.0 0 0.00000 0 0 NS
## 4 0 3.8 0 12.01666 0 38 WN
## 5 12 7.1 0 18.07669 0 57 Post# Display the table
kable(summary_results, caption = "Summary Statistics for Each Behavioural Category per Treatment")| Foraging_mean | Foraging_median | Foraging_sd | Foraging_min | Foraging_max | Locomotion_mean | Locomotion_median | Locomotion_sd | Locomotion_min | Locomotion_max | Attentive_mean | Attentive_median | Attentive_sd | Attentive_min | Attentive_max | Inactive_mean | Inactive_median | Inactive_sd | Inactive_min | Inactive_max | Self.Maintenance_mean | Self.Maintenance_median | Self.Maintenance_sd | Self.Maintenance_min | Self.Maintenance_max | Vocalisation_mean | Vocalisation_median | Vocalisation_sd | Vocalisation_min | Vocalisation_max | Stereotyping_mean | Stereotyping_median | Stereotyping_sd | Stereotyping_min | Stereotyping_max | Play_mean | Play_median | Play_sd | Play_min | Play_max | Treatment |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 2519.7 | 2543.5 | 639.7190 | 1686 | 3467 | 685.5 | 625.0 | 459.8198 | 25 | 1468 | 65.4 | 38.5 | 95.84154 | 0 | 318 | 201.8 | 124.0 | 279.7212 | 0 | 935 | 127.0 | 89.5 | 147.14166 | 0 | 479 | 0.6 | 0 | 1.897367 | 0 | 6 | 0.0 | 0 | 0.000000 | 0 | 0 | 0.0 | 0 | 0.00000 | 0 | 0 | Pre |
| 2106.1 | 1993.5 | 730.9854 | 1072 | 3301 | 878.8 | 806.5 | 507.1396 | 111 | 1777 | 227.1 | 131.5 | 287.47694 | 0 | 905 | 225.2 | 86.0 | 342.6176 | 0 | 997 | 122.6 | 62.0 | 142.50786 | 17 | 486 | 4.2 | 0 | 10.293472 | 0 | 33 | 0.0 | 0 | 0.000000 | 0 | 0 | 32.0 | 0 | 76.38499 | 0 | 236 | CM |
| 2689.1 | 2575.0 | 548.9294 | 2044 | 3600 | 587.3 | 680.0 | 305.4199 | 0 | 929 | 44.1 | 10.5 | 63.69624 | 0 | 187 | 106.4 | 53.0 | 142.3612 | 0 | 417 | 170.4 | 88.5 | 219.30354 | 0 | 716 | 2.7 | 0 | 8.538150 | 0 | 27 | 0.0 | 0 | 0.000000 | 0 | 0 | 0.0 | 0 | 0.00000 | 0 | 0 | NS |
| 2182.1 | 2300.0 | 603.0101 | 1353 | 3095 | 680.9 | 619.5 | 200.6968 | 392 | 1071 | 257.0 | 148.5 | 358.40883 | 0 | 1236 | 338.9 | 184.5 | 357.0174 | 12 | 925 | 137.3 | 52.5 | 252.54441 | 0 | 844 | 0.0 | 0 | 0.000000 | 0 | 0 | 0.0 | 0 | 0.000000 | 0 | 0 | 3.8 | 0 | 12.01666 | 0 | 38 | WN |
| 2223.9 | 2072.0 | 777.0875 | 1030 | 3141 | 774.7 | 720.0 | 338.6047 | 291 | 1413 | 328.8 | 297.0 | 310.30086 | 0 | 821 | 174.3 | 0.0 | 283.2710 | 0 | 675 | 84.7 | 69.0 | 94.31631 | 0 | 289 | 1.0 | 0 | 3.162278 | 0 | 10 | 1.2 | 0 | 3.794733 | 0 | 12 | 7.1 | 0 | 18.07669 | 0 | 57 | Post |
# Check the structure of the data to ensure columns are correctly named
str(total_per_obs)## 'data.frame': 50 obs. of 10 variables:
## $ Session.ID : chr "Pre 1" "Pre 2" "Pre 3" "Pre 4" ...
## $ Treatment : chr "Pre" "Pre" "Pre" "Pre" ...
## $ Foraging : int 2809 3357 3467 1959 1954 2278 2893 1926 1686 2868 ...
## $ Locomotion : int 379 218 25 1371 912 652 580 648 1468 602 ...
## $ Attentive : int 96 0 12 46 318 31 0 91 60 0 ...
## $ Inactive : int 45 0 19 210 283 160 88 935 278 0 ...
## $ Self.Maintenance: int 271 25 77 14 133 479 39 0 102 130 ...
## $ Vocalisation : int 0 0 0 0 0 0 0 0 6 0 ...
## $ Stereotyping : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Play : int 0 0 0 0 0 0 0 0 0 0 ...Calculating the Average Time Spent doing each Behaviour per Treatment Phase
# Calculate the average time spent exhibiting each behaviour in each treatment phase
average_duration <- total_per_obs %>%
group_by(Treatment) %>%
summarise(across(Foraging:Play, \(x) mean(x, na.rm = TRUE)))
print(average_duration)## # A tibble: 5 × 9
## Treatment Foraging Locomotion Attentive Inactive Self.Maintenance Vocalisation
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 CM 2106. 879. 227. 225. 123. 4.2
## 2 NS 2689. 587. 44.1 106. 170. 2.7
## 3 Post 2224. 775. 329. 174. 84.7 1
## 4 Pre 2520. 686. 65.4 202. 127 0.6
## 5 WN 2182. 681. 257 339. 137. 0
## # ℹ 2 more variables: Stereotyping <dbl>, Play <dbl># Reshape data to long format for easier analysis
data_long <- total_per_obs %>%
pivot_longer(cols = Foraging:Play, names_to = "Behaviour", values_to = "TimeSpent")
# Convert variables to factors for GLM analysis
data_long$Treatment <- factor(data_long$Treatment)
data_long$Behaviour <- factor(data_long$Behaviour)Post- Hoc Analyses comparing time spent in behaviour
# Subset the data for each treatment
data_Pre <-filter(data_long, Treatment == "Pre")
data_CM <- filter(data_long, Treatment == "CM")
data_NS <- filter(data_long, Treatment == "NS")
data_Post <- filter(data_long, Treatment == "Post")
data_WN <- filter(data_long, Treatment == "WN")
# Fit GLMs for each treatment group
model_CM <- glm(TimeSpent ~ Behaviour, data = data_CM)
model_NS <- glm(TimeSpent ~ Behaviour, data = data_NS)
model_Post <- glm(TimeSpent ~ Behaviour, data = data_Post)
model_Pre <- glm(TimeSpent ~ Behaviour, data = data_Pre)
model_WN <- glm(TimeSpent ~ Behaviour, data = data_WN)
# Perform Tukey's HSD test for post-hoc analysis within each treatment
post_hoc_CM <- glht(model_CM, linfct = mcp(Behaviour = "Tukey"))
post_hoc_NS <- glht(model_NS, linfct = mcp(Behaviour = "Tukey"))
post_hoc_Post <- glht(model_Post, linfct = mcp(Behaviour = "Tukey"))
post_hoc_Pre <- glht(model_Pre, linfct = mcp(Behaviour = "Tukey"))
post_hoc_WN <- glht(model_WN, linfct = mcp(Behaviour = "Tukey"))
# Output summaries of the post-hoc tests
summary(post_hoc_CM)##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: glm(formula = TimeSpent ~ Behaviour, data = data_CM)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Foraging - Attentive == 0 1879.0 159.5 11.779 < 0.001 ***
## Inactive - Attentive == 0 -1.9 159.5 -0.012 1.00000
## Locomotion - Attentive == 0 651.7 159.5 4.085 0.00118 **
## Play - Attentive == 0 -195.1 159.5 -1.223 0.92538
## Self.Maintenance - Attentive == 0 -104.5 159.5 -0.655 0.99804
## Stereotyping - Attentive == 0 -227.1 159.5 -1.424 0.84649
## Vocalisation - Attentive == 0 -222.9 159.5 -1.397 0.85890
## Inactive - Foraging == 0 -1880.9 159.5 -11.791 < 0.001 ***
## Locomotion - Foraging == 0 -1227.3 159.5 -7.694 < 0.001 ***
## Play - Foraging == 0 -2074.1 159.5 -13.002 < 0.001 ***
## Self.Maintenance - Foraging == 0 -1983.5 159.5 -12.435 < 0.001 ***
## Stereotyping - Foraging == 0 -2106.1 159.5 -13.203 < 0.001 ***
## Vocalisation - Foraging == 0 -2101.9 159.5 -13.177 < 0.001 ***
## Locomotion - Inactive == 0 653.6 159.5 4.097 0.00111 **
## Play - Inactive == 0 -193.2 159.5 -1.211 0.92897
## Self.Maintenance - Inactive == 0 -102.6 159.5 -0.643 0.99826
## Stereotyping - Inactive == 0 -225.2 159.5 -1.412 0.85219
## Vocalisation - Inactive == 0 -221.0 159.5 -1.385 0.86432
## Play - Locomotion == 0 -846.8 159.5 -5.309 < 0.001 ***
## Self.Maintenance - Locomotion == 0 -756.2 159.5 -4.741 < 0.001 ***
## Stereotyping - Locomotion == 0 -878.8 159.5 -5.509 < 0.001 ***
## Vocalisation - Locomotion == 0 -874.6 159.5 -5.483 < 0.001 ***
## Self.Maintenance - Play == 0 90.6 159.5 0.568 0.99922
## Stereotyping - Play == 0 -32.0 159.5 -0.201 1.00000
## Vocalisation - Play == 0 -27.8 159.5 -0.174 1.00000
## Stereotyping - Self.Maintenance == 0 -122.6 159.5 -0.769 0.99466
## Vocalisation - Self.Maintenance == 0 -118.4 159.5 -0.742 0.99570
## Vocalisation - Stereotyping == 0 4.2 159.5 0.026 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)summary(post_hoc_NS)##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: glm(formula = TimeSpent ~ Behaviour, data = data_NS)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Foraging - Attentive == 0 2.645e+03 1.081e+02 24.477 <0.001 ***
## Inactive - Attentive == 0 6.230e+01 1.081e+02 0.577 0.9991
## Locomotion - Attentive == 0 5.432e+02 1.081e+02 5.027 <0.001 ***
## Play - Attentive == 0 -4.410e+01 1.081e+02 -0.408 0.9999
## Self.Maintenance - Attentive == 0 1.263e+02 1.081e+02 1.169 0.9408
## Stereotyping - Attentive == 0 -4.410e+01 1.081e+02 -0.408 0.9999
## Vocalisation - Attentive == 0 -4.140e+01 1.081e+02 -0.383 0.9999
## Inactive - Foraging == 0 -2.583e+03 1.081e+02 -23.900 <0.001 ***
## Locomotion - Foraging == 0 -2.102e+03 1.081e+02 -19.450 <0.001 ***
## Play - Foraging == 0 -2.689e+03 1.081e+02 -24.885 <0.001 ***
## Self.Maintenance - Foraging == 0 -2.519e+03 1.081e+02 -23.308 <0.001 ***
## Stereotyping - Foraging == 0 -2.689e+03 1.081e+02 -24.885 <0.001 ***
## Vocalisation - Foraging == 0 -2.686e+03 1.081e+02 -24.860 <0.001 ***
## Locomotion - Inactive == 0 4.809e+02 1.081e+02 4.450 <0.001 ***
## Play - Inactive == 0 -1.064e+02 1.081e+02 -0.985 0.9767
## Self.Maintenance - Inactive == 0 6.400e+01 1.081e+02 0.592 0.9990
## Stereotyping - Inactive == 0 -1.064e+02 1.081e+02 -0.985 0.9767
## Vocalisation - Inactive == 0 -1.037e+02 1.081e+02 -0.960 0.9799
## Play - Locomotion == 0 -5.873e+02 1.081e+02 -5.435 <0.001 ***
## Self.Maintenance - Locomotion == 0 -4.169e+02 1.081e+02 -3.858 0.0029 **
## Stereotyping - Locomotion == 0 -5.873e+02 1.081e+02 -5.435 <0.001 ***
## Vocalisation - Locomotion == 0 -5.846e+02 1.081e+02 -5.410 <0.001 ***
## Self.Maintenance - Play == 0 1.704e+02 1.081e+02 1.577 0.7642
## Stereotyping - Play == 0 -1.776e-13 1.081e+02 0.000 1.0000
## Vocalisation - Play == 0 2.700e+00 1.081e+02 0.025 1.0000
## Stereotyping - Self.Maintenance == 0 -1.704e+02 1.081e+02 -1.577 0.7642
## Vocalisation - Self.Maintenance == 0 -1.677e+02 1.081e+02 -1.552 0.7790
## Vocalisation - Stereotyping == 0 2.700e+00 1.081e+02 0.025 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)summary(post_hoc_Post)##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: glm(formula = TimeSpent ~ Behaviour, data = data_Post)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Foraging - Attentive == 0 1895.1 150.4 12.604 < 0.001 ***
## Inactive - Attentive == 0 -154.5 150.4 -1.028 0.97040
## Locomotion - Attentive == 0 445.9 150.4 2.966 0.06039 .
## Play - Attentive == 0 -321.7 150.4 -2.140 0.38923
## Self.Maintenance - Attentive == 0 -244.1 150.4 -1.623 0.73615
## Stereotyping - Attentive == 0 -327.6 150.4 -2.179 0.36445
## Vocalisation - Attentive == 0 -327.8 150.4 -2.180 0.36384
## Inactive - Foraging == 0 -2049.6 150.4 -13.632 < 0.001 ***
## Locomotion - Foraging == 0 -1449.2 150.4 -9.638 < 0.001 ***
## Play - Foraging == 0 -2216.8 150.4 -14.744 < 0.001 ***
## Self.Maintenance - Foraging == 0 -2139.2 150.4 -14.227 < 0.001 ***
## Stereotyping - Foraging == 0 -2222.7 150.4 -14.783 < 0.001 ***
## Vocalisation - Foraging == 0 -2222.9 150.4 -14.784 < 0.001 ***
## Locomotion - Inactive == 0 600.4 150.4 3.993 0.00154 **
## Play - Inactive == 0 -167.2 150.4 -1.112 0.95444
## Self.Maintenance - Inactive == 0 -89.6 150.4 -0.596 0.99894
## Stereotyping - Inactive == 0 -173.1 150.4 -1.151 0.94527
## Vocalisation - Inactive == 0 -173.3 150.4 -1.153 0.94497
## Play - Locomotion == 0 -767.6 150.4 -5.105 < 0.001 ***
## Self.Maintenance - Locomotion == 0 -690.0 150.4 -4.589 < 0.001 ***
## Stereotyping - Locomotion == 0 -773.5 150.4 -5.144 < 0.001 ***
## Vocalisation - Locomotion == 0 -773.7 150.4 -5.146 < 0.001 ***
## Self.Maintenance - Play == 0 77.6 150.4 0.516 0.99959
## Stereotyping - Play == 0 -5.9 150.4 -0.039 1.00000
## Vocalisation - Play == 0 -6.1 150.4 -0.041 1.00000
## Stereotyping - Self.Maintenance == 0 -83.5 150.4 -0.555 0.99933
## Vocalisation - Self.Maintenance == 0 -83.7 150.4 -0.557 0.99932
## Vocalisation - Stereotyping == 0 -0.2 150.4 -0.001 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)summary(post_hoc_Pre)##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: glm(formula = TimeSpent ~ Behaviour, data = data_Pre)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Foraging - Attentive == 0 2.454e+03 1.351e+02 18.171 < 0.001 ***
## Inactive - Attentive == 0 1.364e+02 1.351e+02 1.010 0.97314
## Locomotion - Attentive == 0 6.201e+02 1.351e+02 4.591 < 0.001 ***
## Play - Attentive == 0 -6.540e+01 1.351e+02 -0.484 0.99973
## Self.Maintenance - Attentive == 0 6.160e+01 1.351e+02 0.456 0.99982
## Stereotyping - Attentive == 0 -6.540e+01 1.351e+02 -0.484 0.99973
## Vocalisation - Attentive == 0 -6.480e+01 1.351e+02 -0.480 0.99975
## Inactive - Foraging == 0 -2.318e+03 1.351e+02 -17.161 < 0.001 ***
## Locomotion - Foraging == 0 -1.834e+03 1.351e+02 -13.580 < 0.001 ***
## Play - Foraging == 0 -2.520e+03 1.351e+02 -18.655 < 0.001 ***
## Self.Maintenance - Foraging == 0 -2.393e+03 1.351e+02 -17.714 < 0.001 ***
## Stereotyping - Foraging == 0 -2.520e+03 1.351e+02 -18.655 < 0.001 ***
## Vocalisation - Foraging == 0 -2.519e+03 1.351e+02 -18.650 < 0.001 ***
## Locomotion - Inactive == 0 4.837e+02 1.351e+02 3.581 0.00836 **
## Play - Inactive == 0 -2.018e+02 1.351e+02 -1.494 0.81105
## Self.Maintenance - Inactive == 0 -7.480e+01 1.351e+02 -0.554 0.99934
## Stereotyping - Inactive == 0 -2.018e+02 1.351e+02 -1.494 0.81108
## Vocalisation - Inactive == 0 -2.012e+02 1.351e+02 -1.490 0.81350
## Play - Locomotion == 0 -6.855e+02 1.351e+02 -5.075 < 0.001 ***
## Self.Maintenance - Locomotion == 0 -5.585e+02 1.351e+02 -4.135 < 0.001 ***
## Stereotyping - Locomotion == 0 -6.855e+02 1.351e+02 -5.075 < 0.001 ***
## Vocalisation - Locomotion == 0 -6.849e+02 1.351e+02 -5.071 < 0.001 ***
## Self.Maintenance - Play == 0 1.270e+02 1.351e+02 0.940 0.98211
## Stereotyping - Play == 0 -4.263e-14 1.351e+02 0.000 1.00000
## Vocalisation - Play == 0 6.000e-01 1.351e+02 0.004 1.00000
## Stereotyping - Self.Maintenance == 0 -1.270e+02 1.351e+02 -0.940 0.98212
## Vocalisation - Self.Maintenance == 0 -1.264e+02 1.351e+02 -0.936 0.98260
## Vocalisation - Stereotyping == 0 6.000e-01 1.351e+02 0.004 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)summary(post_hoc_WN)## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: glm(formula = TimeSpent ~ Behaviour, data = data_WN)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Foraging - Attentive == 0 1.925e+03 1.345e+02 14.312 <0.01 ***
## Inactive - Attentive == 0 8.190e+01 1.345e+02 0.609 0.9988
## Locomotion - Attentive == 0 4.239e+02 1.345e+02 3.151 0.0349 *
## Play - Attentive == 0 -2.532e+02 1.345e+02 -1.882 0.5631
## Self.Maintenance - Attentive == 0 -1.197e+02 1.345e+02 -0.890 0.9870
## Stereotyping - Attentive == 0 -2.570e+02 1.345e+02 -1.911 0.5437
## Vocalisation - Attentive == 0 -2.570e+02 1.345e+02 -1.911 0.5438
## Inactive - Foraging == 0 -1.843e+03 1.345e+02 -13.703 <0.01 ***
## Locomotion - Foraging == 0 -1.501e+03 1.345e+02 -11.160 <0.01 ***
## Play - Foraging == 0 -2.178e+03 1.345e+02 -16.194 <0.01 ***
## Self.Maintenance - Foraging == 0 -2.045e+03 1.345e+02 -15.202 <0.01 ***
## Stereotyping - Foraging == 0 -2.182e+03 1.345e+02 -16.222 <0.01 ***
## Vocalisation - Foraging == 0 -2.182e+03 1.345e+02 -16.222 <0.01 ***
## Locomotion - Inactive == 0 3.420e+02 1.345e+02 2.543 0.1776
## Play - Inactive == 0 -3.351e+02 1.345e+02 -2.491 0.1997
## Self.Maintenance - Inactive == 0 -2.016e+02 1.345e+02 -1.499 0.8085
## Stereotyping - Inactive == 0 -3.389e+02 1.345e+02 -2.519 0.1874
## Vocalisation - Inactive == 0 -3.389e+02 1.345e+02 -2.519 0.1867
## Play - Locomotion == 0 -6.771e+02 1.345e+02 -5.034 <0.01 ***
## Self.Maintenance - Locomotion == 0 -5.436e+02 1.345e+02 -4.041 <0.01 **
## Stereotyping - Locomotion == 0 -6.809e+02 1.345e+02 -5.062 <0.01 ***
## Vocalisation - Locomotion == 0 -6.809e+02 1.345e+02 -5.062 <0.01 ***
## Self.Maintenance - Play == 0 1.335e+02 1.345e+02 0.992 0.9756
## Stereotyping - Play == 0 -3.800e+00 1.345e+02 -0.028 1.0000
## Vocalisation - Play == 0 -3.800e+00 1.345e+02 -0.028 1.0000
## Stereotyping - Self.Maintenance == 0 -1.373e+02 1.345e+02 -1.021 0.9715
## Vocalisation - Self.Maintenance == 0 -1.373e+02 1.345e+02 -1.021 0.9715
## Vocalisation - Stereotyping == 0 5.684e-14 1.345e+02 0.000 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)# Purpose: To compare the means of different behavioral categories (e.g., Foraging, Attentive, Inactive, etc.) in relation to the time spent on these behaviors, using Tukey's HSD test to identify which pairs of behavioral means are significantly different from one another.
# Interpretation: Overall, Foraging consistently shows a significantly higher time spent compared to Attentive across all data sets. Locomotion also shows a significant increase in time compared to Attentive in most cases. Other comparisons either do not show significant differences or are only marginally significant in some data sets. Fitting a Simple GLM
# Fit a Generalised Linear Model (GLM) to analyse the effect of treatment and behaviour on time spent carrying out each behaviour
# Convert Treatment to a factor if it isn't already
data_long$Treatment <- factor(data_long$Treatment)
# Set 'Pre' as the reference level
data_long$Treatment <- relevel(data_long$Treatment, ref = "Pre")
# Running model with both main effects of Treatment & Behaviour, and their interaction
model <- glm(TimeSpent ~ Treatment * Behaviour, data = data_long)
summary(model)##
## Call:
## glm(formula = TimeSpent ~ Treatment * Behaviour, data = data_long)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1193.90 -88.57 -1.45 3.40 1194.90
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 65.40 98.01 0.667 0.50504
## TreatmentCM 161.70 138.61 1.167 0.24416
## TreatmentNS -21.30 138.61 -0.154 0.87796
## TreatmentPost 263.40 138.61 1.900 0.05820
## TreatmentWN 191.60 138.61 1.382 0.16775
## BehaviourForaging 2454.30 138.61 17.706 < 2e-16
## BehaviourInactive 136.40 138.61 0.984 0.32576
## BehaviourLocomotion 620.10 138.61 4.474 1.03e-05
## BehaviourPlay -65.40 138.61 -0.472 0.63734
## BehaviourSelf.Maintenance 61.60 138.61 0.444 0.65702
## BehaviourStereotyping -65.40 138.61 -0.472 0.63734
## BehaviourVocalisation -64.80 138.61 -0.467 0.64043
## TreatmentCM:BehaviourForaging -575.30 196.03 -2.935 0.00355
## TreatmentNS:BehaviourForaging 190.70 196.03 0.973 0.33130
## TreatmentPost:BehaviourForaging -559.20 196.03 -2.853 0.00459
## TreatmentWN:BehaviourForaging -529.20 196.03 -2.700 0.00727
## TreatmentCM:BehaviourInactive -138.30 196.03 -0.706 0.48095
## TreatmentNS:BehaviourInactive -74.10 196.03 -0.378 0.70565
## TreatmentPost:BehaviourInactive -290.90 196.03 -1.484 0.13869
## TreatmentWN:BehaviourInactive -54.50 196.03 -0.278 0.78116
## TreatmentCM:BehaviourLocomotion 31.60 196.03 0.161 0.87203
## TreatmentNS:BehaviourLocomotion -76.90 196.03 -0.392 0.69508
## TreatmentPost:BehaviourLocomotion -174.20 196.03 -0.889 0.37479
## TreatmentWN:BehaviourLocomotion -196.20 196.03 -1.001 0.31756
## TreatmentCM:BehaviourPlay -129.70 196.03 -0.662 0.50863
## TreatmentNS:BehaviourPlay 21.30 196.03 0.109 0.91353
## TreatmentPost:BehaviourPlay -256.30 196.03 -1.307 0.19189
## TreatmentWN:BehaviourPlay -187.80 196.03 -0.958 0.33869
## TreatmentCM:BehaviourSelf.Maintenance -166.10 196.03 -0.847 0.39738
## TreatmentNS:BehaviourSelf.Maintenance 64.70 196.03 0.330 0.74155
## TreatmentPost:BehaviourSelf.Maintenance -305.70 196.03 -1.559 0.11976
## TreatmentWN:BehaviourSelf.Maintenance -181.30 196.03 -0.925 0.35565
## TreatmentCM:BehaviourStereotyping -161.70 196.03 -0.825 0.40998
## TreatmentNS:BehaviourStereotyping 21.30 196.03 0.109 0.91353
## TreatmentPost:BehaviourStereotyping -262.20 196.03 -1.338 0.18188
## TreatmentWN:BehaviourStereotyping -191.60 196.03 -0.977 0.32902
## TreatmentCM:BehaviourVocalisation -158.10 196.03 -0.807 0.42048
## TreatmentNS:BehaviourVocalisation 23.40 196.03 0.119 0.90505
## TreatmentPost:BehaviourVocalisation -263.00 196.03 -1.342 0.18056
## TreatmentWN:BehaviourVocalisation -192.20 196.03 -0.980 0.32751
##
## (Intercept)
## TreatmentCM
## TreatmentNS
## TreatmentPost .
## TreatmentWN
## BehaviourForaging ***
## BehaviourInactive
## BehaviourLocomotion ***
## BehaviourPlay
## BehaviourSelf.Maintenance
## BehaviourStereotyping
## BehaviourVocalisation
## TreatmentCM:BehaviourForaging **
## TreatmentNS:BehaviourForaging
## TreatmentPost:BehaviourForaging **
## TreatmentWN:BehaviourForaging **
## TreatmentCM:BehaviourInactive
## TreatmentNS:BehaviourInactive
## TreatmentPost:BehaviourInactive
## TreatmentWN:BehaviourInactive
## TreatmentCM:BehaviourLocomotion
## TreatmentNS:BehaviourLocomotion
## TreatmentPost:BehaviourLocomotion
## TreatmentWN:BehaviourLocomotion
## TreatmentCM:BehaviourPlay
## TreatmentNS:BehaviourPlay
## TreatmentPost:BehaviourPlay
## TreatmentWN:BehaviourPlay
## TreatmentCM:BehaviourSelf.Maintenance
## TreatmentNS:BehaviourSelf.Maintenance
## TreatmentPost:BehaviourSelf.Maintenance
## TreatmentWN:BehaviourSelf.Maintenance
## TreatmentCM:BehaviourStereotyping
## TreatmentNS:BehaviourStereotyping
## TreatmentPost:BehaviourStereotyping
## TreatmentWN:BehaviourStereotyping
## TreatmentCM:BehaviourVocalisation
## TreatmentNS:BehaviourVocalisation
## TreatmentPost:BehaviourVocalisation
## TreatmentWN:BehaviourVocalisation
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 96067.48)
##
## Null deviance: 263071848 on 399 degrees of freedom
## Residual deviance: 34584292 on 360 degrees of freedom
## AIC: 5764.1
##
## Number of Fisher Scoring iterations: 2# Running model includes only the main effects of Treatment & Behaviour.
model_2 <- glm(TimeSpent ~ Treatment + Behaviour, data = data_long)
# Testing which model fits data better
anova(model_2,model)# Calculating p-value
1-pchisq(3014.8, 28)## [1] 0# p value = 0
# Improvement in model fit by including the interaction term (in Model 1) is statistically significant
# Model with interaction (1) improves the model fit.
# Foraging is found to be significantly influenced by treatment, specifically it is decreased significantly during the CM, WN and Post phases. # Saving the output into a table to include in paper
# Tidy the model output
model_summary <- tidy(model)
# Saving the summary to a CSV file
write.csv(model_summary, file = "model_summary.csv", row.names = FALSE)Residual Analysis
# Residual plots for model with interaction
par(mfrow = c(2, 2))
plot(residuals(model), fitted(model), main = "Residuals vs Fitted")
abline(h = 0, col = "red")
# Q-Q plot of residuals
qqnorm(residuals(model))
qqline(residuals(model), col = "red")
# Scale-Location plot
plot(fitted(model), sqrt(abs(residuals(model))), main = "Scale-Location")
abline(h = 0, col = "red")
# Residuals vs Leverage plot
plot(hatvalues(model), residuals(model), main = "Residuals vs Leverage")
abline(h = 0, col = "red")
# Reset graphics layout
par(mfrow = c(1, 1))
# Residual plots look slightly off, but this is believed to be due to the amount of 0's within the data set.
# We continue as they are adequate. Analysis of Treatment Effects on Time Spent for Each Behavior - GLM
# List of behaviors to analyze
behaviors <- c("Foraging", "Attentive", "Locomotion", "Inactive",
"Vocalisation", "Self.Maintenance", "Play", "Stereotypy")
# Loop through each behavior and perform the analysis
for (behavior in behaviors) {
# Subset the data for the current behavior
data_behavior <- subset(data_long, Behaviour == behavior)
# Check if the Treatment variable has more than one level
if (length(unique(data_behavior$Treatment)) > 1) {
# Fit the model with Treatment as a predictor
model_with_treatment <- glm(TimeSpent ~ Treatment, data = data_behavior, family = gaussian())
# Fit the null model (no predictors, only intercept)
null_model <- glm(TimeSpent ~ 1, data = data_behavior, family = gaussian())
# Perform LRT to compare the two models
lrt_result <- anova(null_model, model_with_treatment, test = "Chisq")
# Print the results for the current behavior
cat("LRT Results for", behavior, ":\n")
print(lrt_result)
cat("\n---------------------------------------\n")
} else {
# If there's only one level in Treatment, skip this behavior
cat("Skipping", behavior, "- Treatment has only one level.\n")
cat("\n---------------------------------------\n")
}
}## LRT Results for Foraging :
## Analysis of Deviance Table
##
## Model 1: TimeSpent ~ 1
## Model 2: TimeSpent ~ Treatment
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 49 22383471
## 2 45 19911508 4 2471964 0.2322
##
## ---------------------------------------
## LRT Results for Attentive :
## Analysis of Deviance Table
##
## Model 1: TimeSpent ~ 1
## Model 2: TimeSpent ~ Treatment
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 49 3503568
## 2 45 2885664 4 617905 0.04703 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ---------------------------------------
## LRT Results for Locomotion :
## Analysis of Deviance Table
##
## Model 1: TimeSpent ~ 1
## Model 2: TimeSpent ~ Treatment
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 49 6936822
## 2 45 6451547 4 485275 0.4956
##
## ---------------------------------------
## LRT Results for Inactive :
## Analysis of Deviance Table
##
## Model 1: TimeSpent ~ 1
## Model 2: TimeSpent ~ Treatment
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 49 4101599
## 2 45 3812413 4 289186 0.4912
##
## ---------------------------------------
## LRT Results for Vocalisation :
## Analysis of Deviance Table
##
## Model 1: TimeSpent ~ 1
## Model 2: TimeSpent ~ Treatment
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 49 1850.5
## 2 45 1732.1 4 118.4 0.5452
##
## ---------------------------------------
## LRT Results for Self.Maintenance :
## Analysis of Deviance Table
##
## Model 1: TimeSpent ~ 1
## Model 2: TimeSpent ~ Treatment
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 49 1502432
## 2 45 1464547 4 37885 0.884
##
## ---------------------------------------
## LRT Results for Play :
## Analysis of Deviance Table
##
## Model 1: TimeSpent ~ 1
## Model 2: TimeSpent ~ Treatment
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 49 63960
## 2 45 56752 4 7207.7 0.2215
##
## ---------------------------------------
## Skipping Stereotypy - Treatment has only one level.
##
## ---------------------------------------# In general, this output suggests that, except for a possible marginal effect on Attentive behavior, Treatment does not significantly influence the time spent on different behaviors in your data set. Residual Analysis
# Residual plots for the current behavior model
par(mfrow = c(2, 2))
# Plot residuals vs fitted values
plot(residuals(model_with_treatment), fitted(model_with_treatment),
main = paste("Residuals vs Fitted for", behavior))
abline(h = 0, col = "red")
# Q-Q plot of residuals
qqnorm(residuals(model_with_treatment))
qqline(residuals(model_with_treatment), col = "red")
# Scale-Location plot
plot(fitted(model_with_treatment), sqrt(abs(residuals(model_with_treatment))),
main = paste("Scale-Location for", behavior))
abline(h = 0, col = "red")
# Residuals vs Leverage plot
plot(hatvalues(model_with_treatment), residuals(model_with_treatment),
main = paste("Residuals vs Leverage for", behavior))## Warning in plot.window(...): axis(1, *): range of values ( 0) is small wrt |M|
## = 0.1 --> not pretty() abline(h = 0, col = "red")
# Reset graphics layout
par(mfrow = c(1, 1))
# Residuals are slightly better here than previously.Pairwise t-Tests for Treatment Levels in Attentive Behavior
# Subset the data for the Attentive behavior as value is close to 0.05
data_attentive <- subset(data_long, Behaviour == "Attentive")
# Pairwise t-tests between Treatment levels for Attentive behavior
pairwise_attentive <- pairwise.t.test(data_attentive$TimeSpent, data_attentive$Treatment, p.adjust.method = "bonferroni")
# Print the pairwise comparisons results
print(pairwise_attentive)##
## Pairwise comparisons using t tests with pooled SD
##
## data: data_attentive$TimeSpent and data_attentive$Treatment
##
## Pre CM NS Post
## CM 1.00 - - -
## NS 1.00 1.00 - -
## Post 0.25 1.00 0.16 -
## WN 0.98 1.00 0.67 1.00
##
## P value adjustment method: bonferroni# Pairwise comparisons did not find any significant differences between treatment groups for behaviour 'Attentive'.Visualising Main Effects of Treatment on Behaviour
# Set the order of the treatments
data_long$Treatment <- factor(data_long$Treatment, levels = c("Pre", "CM", "NS", "WN", "Post"))
# Calculate IQR and define outliers
data_cleaned <- data_long %>%
group_by(Treatment, Behaviour) %>%
mutate(Q1 = quantile(TimeSpent, 0.25, na.rm = TRUE),
Q3 = quantile(TimeSpent, 0.75, na.rm = TRUE),
IQR = Q3 - Q1,
lower_bound = Q1 - 1.5 * IQR,
upper_bound = Q3 + 1.5 * IQR) %>%
filter(TimeSpent >= lower_bound & TimeSpent <= upper_bound)
# Assigning colour to each behaviour
behaviour_colors <- c("Attentive" = "#4E79A7", # Blue
"Foraging" = "#F28E2B", # Orange
"Inactive" = "#E15759", # Red
"Locomotion" = "#76B7B2", # Teal
"Play" = "#59A14F", # Green
"Self.Maintenance" = "#EDC948", # Yellow
"Stereotyping" = "#AF7AA1", # Purple
"Vocalisation" = "#FF9DA7") # Pink
# Plotting boxplot of time spent exhibiting each behaviour during each treatment phase without outliers
ggplot(data_cleaned, aes(x = Behaviour, y = TimeSpent, fill = Behaviour)) +
geom_boxplot() +
facet_wrap(~ Treatment) +
scale_fill_manual(values = behaviour_colors) +
labs(y = "Total time (secs) spent exhibiting each behaviour",
x = "Behaviour") +
coord_cartesian(ylim = c(0, max(data_cleaned$TimeSpent, na.rm = TRUE))) + # Adjust y-axis limits if needed
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
text = element_text(family = "Arial", color = "#333333"),
plot.title = element_text(size = 14, face = "bold", margin = margin(b = 10)),
axis.title = element_text(size = 10),
axis.title.x = element_text(size = 10, margin = margin(t = 10)), # Move x-axis title down
legend.title = element_text(size = 12),
legend.text = element_text(size = 10))
ENCLOSURE USE ANALYSIS
Data Preparation
# Create the dataset with cleaned column names
enclosure_use_data <- data.frame(
Treatment = c("Pre_Enrichment", "Classical_Music", "Natural_Sound", "White_Noise", "Post_Enrichment"),
A = c(4382, 8491, 3083, 10356, 7863),
B = c(664, 1008, 1858, 1292, 1304),
C = c(2949, 2774, 1923, 3743, 2223),
D = c(2072, 2431, 2318, 4324, 612),
Inner_Circle = c(10048, 7393, 12682, 3167, 3010),
Outer_Circle = c(13333, 7655, 10850, 8502, 15773)
)
# Purpose: Define a dataset where each treatment's effect is recorded across different locations (A, B, C, D, Inner Circle, Outer Circle).
# Reshape the data to long format
data_long_location <- enclosure_use_data %>%
pivot_longer(cols = A:Outer_Circle, names_to = "Location", values_to = "Time")
# Convert columns to factors
data_long_location$Treatment <- as.factor(data_long_location$Treatment)
data_long_location$Location <- as.factor(data_long_location$Location)
# Set 'Pre' as the reference level
data_long_location$Treatment <- relevel(data_long_location$Treatment, ref = "Pre_Enrichment")Fitting a Simple Linear Model to analyse effect of Treatment on Time spent in each Location
# Fitting a Linear model
model_lm <- lm(Time ~ Location + Treatment, data = data_long_location)
summary(model_lm)##
## Call:
## lm(formula = Time ~ Location + Treatment, data = data_long_location)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4111.4 -953.3 99.6 1139.8 5239.1
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7140.2 1643.3 4.345 0.000314 ***
## LocationB -5609.8 1800.1 -3.116 0.005439 **
## LocationC -4112.6 1800.1 -2.285 0.033392 *
## LocationD -4483.6 1800.1 -2.491 0.021657 *
## LocationInner_Circle 425.0 1800.1 0.236 0.815758
## LocationOuter_Circle 4387.6 1800.1 2.437 0.024256 *
## TreatmentClassical_Music -616.0 1643.3 -0.375 0.711708
## TreatmentNatural_Sound -122.3 1643.3 -0.074 0.941396
## TreatmentPost_Enrichment -443.8 1643.3 -0.270 0.789857
## TreatmentWhite_Noise -344.0 1643.3 -0.209 0.836301
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2846 on 20 degrees of freedom
## Multiple R-squared: 0.694, Adjusted R-squared: 0.5564
## F-statistic: 5.041 on 9 and 20 DF, p-value: 0.001264# Purpose: Fit a linear model to examine the effects of Location and Treatment on Time.
# Interpretation: Subject spends significantly more time in the Outer Circle Overall but spends most time in the A section when indoors (same amount of time as inner circle).Conducting an LRT to test for significance of treatment on time spent in each location
# Fit a simpler model, e.g., without the Treatment factor
model_simple <- lm(Time ~ Location, data = data_long_location)
# Perform LRT
anova(model_simple, model_lm)# The high p-value (0.9958) in the ANOVA table indicates that the Treatment variable does not significantly improve the model's fit compared to the model that only includes Location. In other words, the effect of Treatment on Time is not statistically significant when considering the variability explained by Location.Plotting Residuals vs Fitted Values
# Residuals vs Fitted Values Plot
plot(model_lm$residuals ~ model_lm$fitted.values, main = "Residuals vs Fitted",
xlab = "Fitted Values", ylab = "Residuals")
abline(h = 0, col = "red")
# Q-Q Plot
qqnorm(model_lm$residuals, main = "Normal Q-Q Plot")
qqline(model_lm$residuals, col = "red")
# Visualising the Enclosure Use of the Subject
# Your data
treatments <- c("Pre-Enrichment", "Classical Music", "Natural Sound", "White Noise", "Post Enrichment")
A <- c(4382, 8491, 3083, 10356, 7863)
B <- c(664, 1008, 1858, 1292, 1304)
C <- c(2949, 2774, 1923, 3743, 2223)
D <- c(2072, 2431, 2318, 4324, 612)
inner_circle <- c(10048, 7393, 12682, 3167, 3010)
outer_circle <- c(13333, 7655, 10850, 8502, 15773)
# Create a data frame
circle_data <- data.frame(
Treatment = treatments,
A = A,
B = B,
C = C,
D = D,
Inner_Circle = inner_circle,
Outer_Circle = outer_circle
)
# Reshape the data to long format for ggplot
circle_data_long <- circle_data %>%
pivot_longer(cols = A:Outer_Circle, names_to = "Location", values_to = "Time")
# Convert 'Treatment' to a factor to ensure correct ordering in the plot
circle_data_long$Treatment <- factor(circle_data_long$Treatment, levels = treatments)
# Rename the locations to remove underscores and keep A-D as they are
circle_data_long$Location <- factor(circle_data_long$Location,
labels = c("A", "B", "C", "D", "Inner Circle", "Outer Circle"))
# Create a professional color palette
professional_colors <- c("A" = "#4E79A7", "B" = "#F28E2B", "C" = "#E15759",
"D" = "#76B7B2", "Inner Circle" = "#59A14F", "Outer Circle" = "#EDC948")
# Create the bar plot showing time spent in each section across treatment phase
ggplot(circle_data_long, aes(x = Treatment, y = Time, fill = Location)) +
geom_bar(stat = "identity", position = "dodge") +
scale_fill_manual(values = professional_colors) +
labs(
x = "Treatment Phase",
y = "Time Spent (seconds)",
fill = "Location") +
theme_minimal() +
theme(
text = element_text(family = "Arial", color = "#333333"),
plot.title = element_text(size = 14, face = "bold", margin = margin(b = 10)), # Adjusted margin to increase spacing
axis.title = element_text(size = 12),
legend.title = element_text(size = 12),
legend.text = element_text(size = 10)
)
# Creating a Heat Map to visualise the Indoor Enclosure Use
# Manually creating the data frame
data <- data.frame(
Treatment = c("Pre-Enrichment", "Classical Music", "Natural Sound", "White Noise", "Post Enrichment"),
A1 = c(436, 5202, 96, 5742, 2547),
A2 = c(1828, 1153, 2360, 3661, 3338),
A3 = c(1574, 1376, 273, 672, 884),
A4 = c(544, 760, 354, 281, 1094),
B1 = c(312, 260, 911, 797, 145),
B2 = c(326, 721, 947, 495, 1092),
B3 = c(0, 0, 0, 0, 67),
B4 = c(26, 27, 0, 0, 0),
C1 = c(934, 1077, 45, 2275, 367),
C2 = c(457, 484, 1279, 31, 1281),
C3 = c(1066, 881, 553, 1269, 418),
C4 = c(492, 332, 46, 168, 157),
D1 = c(464, 1077, 883, 1009, 24),
D2 = c(53, 835, 899, 454, 423),
D3 = c(94, 32, 99, 770, 95),
D4 = c(1461, 487, 437, 2091, 70)
)
# Convert the data frame to long format
data_long_location <- melt(data, id.vars = "Treatment", variable.name = "Location", value.name = "Usage")
# Extract rows and columns from Location
data_long_location$Row <- factor(substr(data_long_location$Location, 2, 2), levels = c("1", "2", "3", "4"))
data_long_location$Column <- factor(substr(data_long_location$Location, 1, 1), levels = c("A", "B", "C", "D"))
# Create the heatmap with facet_wrap showing enclosure use in indoor house
ggplot(data_long_location, aes(x = Column, y = Row, fill = Usage)) +
geom_tile(color = "white") +
scale_fill_gradient(low = "white", high = "steelblue", name = "Time (s)") +
labs(x = "Column",
y = "Row") +
facet_wrap(~ Treatment, ncol = 3) + # Adjust number of columns in the grid
theme_minimal() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
axis.text.y = element_text(size = 10),
strip.text = element_text(size = 12),
legend.position = "bottom",
legend.title = element_text(size = 12),
legend.text = element_text(size = 10),
legend.key.width = unit(2, "cm"), # Adjust the width of the legend key
plot.title = element_text(size = 14, face = "bold"),
panel.grid = element_blank()
)