# CARGA DE LIBRERÍAS
library(readxl)
library(glmmTMB)
library(DHARMa)
## This is DHARMa 0.4.7. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')
library(performance)
library(gtools)
library(emmeans)
## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'
library(glmnet)
## S'està carregant el paquet requerit: Matrix
## Loaded glmnet 4.1-8
##
## S'està adjuntant el paquet: 'glmnet'
## L'objecte següent està emmascarat per 'package:gtools':
##
## na.replace
library(ggplot2)
library(dplyr)
##
## S'està adjuntant el paquet: 'dplyr'
## Els següents objectes estan emmascarats des de 'package:stats':
##
## filter, lag
## Els següents objectes estan emmascarats des de 'package:base':
##
## intersect, setdiff, setequal, union
library(tidyr)
##
## S'està adjuntant el paquet: 'tidyr'
## Els següents objectes estan emmascarats des de 'package:Matrix':
##
## expand, pack, unpack
# LECTURA DE DATOS
xdata <- read_excel("C:/Users/prof/Desktop/TFG/Resultats_Arruis_SOCIAL.xlsx")
# INSPECCIÓN DE LOS DATOS
head(xdata)
## # A tibble: 6 × 12
## subject correct sex age.months trial.nr condition corr.side select Species
## <chr> <dbl> <chr> <dbl> <dbl> <chr> <chr> <dbl> <chr>
## 1 Blue 1 M 11 1 Touching.c… R 1 Ammotr…
## 2 Blue 1 M 11 2 Control R 1 Ammotr…
## 3 Blue 1 M 11 3 Pointing.c… L 1 Ammotr…
## 4 Blue 0 M 11 4 Pointing.a… R 1 Ammotr…
## 5 Blue 1 M 11 5 Touching.p… R 1 Ammotr…
## 6 Blue 0 M 11 6 Gazing.asy… R 1 Ammotr…
## # ℹ 3 more variables: ID <chr>, Session <dbl>, Trial <dbl>
str(xdata)
## tibble [576 × 12] (S3: tbl_df/tbl/data.frame)
## $ subject : chr [1:576] "Blue" "Blue" "Blue" "Blue" ...
## $ correct : num [1:576] 1 1 1 0 1 0 1 1 1 1 ...
## $ sex : chr [1:576] "M" "M" "M" "M" ...
## $ age.months: num [1:576] 11 11 11 11 11 11 11 11 11 11 ...
## $ trial.nr : num [1:576] 1 2 3 4 5 6 7 8 9 10 ...
## $ condition : chr [1:576] "Touching.crossed" "Control" "Pointing.crossed" "Pointing.asymmetrical" ...
## $ corr.side : chr [1:576] "R" "R" "L" "R" ...
## $ select : num [1:576] 1 1 1 1 1 1 1 1 1 1 ...
## $ Species : chr [1:576] "Ammotragus lervia" "Ammotragus lervia" "Ammotragus lervia" "Ammotragus lervia" ...
## $ ID : chr [1:576] "ID 14428" "ID 14428" "ID 14428" "ID 14428" ...
## $ Session : num [1:576] 1 1 1 1 1 1 1 1 1 1 ...
## $ Trial : num [1:576] 1 2 3 4 5 6 7 8 9 10 ...
# ESCALADO DE VARIABLES
xdata$z.age = as.vector(scale(xdata$age.months))
xdata$z.trial = as.vector(scale(xdata$trial.nr))
# MODELO COMPLETO
full = glmmTMB(correct ~ sex + z.age + condition + z.trial + corr.side,
data = xdata, family = binomial)
# VERIFICACIONES
any(is.na(xdata))
## [1] TRUE
colSums(is.na(xdata))
## subject correct sex age.months trial.nr condition corr.side
## 0 0 0 0 0 0 0
## select Species ID Session Trial z.age z.trial
## 288 0 0 0 0 0 0
check_collinearity(full)
## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
## sex 1.93 [1.73, 2.19] 1.39 0.52 [0.46, 0.58]
## z.age 1.92 [1.72, 2.18] 1.39 0.52 [0.46, 0.58]
## condition 1.03 [1.00, 1.69] 1.01 0.98 [0.59, 1.00]
## z.trial 1.01 [1.00, 735.17] 1.00 0.99 [0.00, 1.00]
## corr.side 1.02 [1.00, 2.10] 1.01 0.98 [0.48, 1.00]
testDispersion(full)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.0003, p-value = 1
## alternative hypothesis: two.sided
# COMPARACIÓN CON MODELO NULO
null = glmmTMB(correct ~ z.trial + corr.side,
data = xdata, family = binomial)
anova(full, null, test = "Chisq")
## Data: xdata
## Models:
## null: correct ~ z.trial + corr.side, zi=~0, disp=~1
## full: correct ~ sex + z.age + condition + z.trial + corr.side, zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## null 3 735.70 748.76 -364.85 729.70
## full 13 687.15 743.78 -330.58 661.15 68.54 10 8.476e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# TEST INDIVIDUALES DE EFECTOS
drop1(full, test = "Chisq")
## Single term deletions
##
## Model:
## correct ~ sex + z.age + condition + z.trial + corr.side
## Df AIC LRT Pr(>Chi)
## <none> 687.15
## sex 1 689.74 4.582 0.03231 *
## z.age 1 685.89 0.731 0.39250
## condition 8 734.86 63.710 8.678e-11 ***
## z.trial 1 689.79 4.640 0.03124 *
## corr.side 1 711.20 26.043 3.340e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# RESUMEN DEL MODELO
summary(full)
## Family: binomial ( logit )
## Formula: correct ~ sex + z.age + condition + z.trial + corr.side
## Data: xdata
##
## AIC BIC logLik -2*log(L) df.resid
## 687.2 743.8 -330.6 661.2 563
##
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.31327 0.29393 1.066 0.28652
## sexM 0.58446 0.27324 2.139 0.03243 *
## z.age 0.11180 0.13077 0.855 0.39257
## conditionGazing.asymmetrical -0.04374 0.36743 -0.119 0.90524
## conditionGazing.central 0.49120 0.37281 1.318 0.18765
## conditionPointing.asymmetrical -0.40028 0.36907 -1.085 0.27813
## conditionPointing.crossed 0.16544 0.36797 0.450 0.65301
## conditionPointing.proximal 0.35815 0.37072 0.966 0.33400
## conditionTouching.asymmetrical 1.73735 0.43813 3.965 7.33e-05 ***
## conditionTouching.crossed 1.50624 0.41953 3.590 0.00033 ***
## conditionTouching.proximal 1.73047 0.43791 3.952 7.76e-05 ***
## z.trial -0.20344 0.09490 -2.144 0.03205 *
## corr.sideR -0.95624 0.19108 -5.004 5.60e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
round(confint(full), digits = 2)
## 2.5 % 97.5 % Estimate
## (Intercept) -0.26 0.89 0.31
## sexM 0.05 1.12 0.58
## z.age -0.14 0.37 0.11
## conditionGazing.asymmetrical -0.76 0.68 -0.04
## conditionGazing.central -0.24 1.22 0.49
## conditionPointing.asymmetrical -1.12 0.32 -0.40
## conditionPointing.crossed -0.56 0.89 0.17
## conditionPointing.proximal -0.37 1.08 0.36
## conditionTouching.asymmetrical 0.88 2.60 1.74
## conditionTouching.crossed 0.68 2.33 1.51
## conditionTouching.proximal 0.87 2.59 1.73
## z.trial -0.39 -0.02 -0.20
## corr.sideR -1.33 -0.58 -0.96
# CONTRASTES ENTRE NIVELES DE LA CONDICIÓN
emmeans(full, pairwise ~ condition, type = "response")
## $emmeans
## condition prob SE df asymp.LCL asymp.UCL
## Control 0.532 0.0652 Inf 0.405 0.655
## Gazing.asymmetrical 0.521 0.0654 Inf 0.394 0.645
## Gazing.central 0.650 0.0615 Inf 0.522 0.759
## Pointing.asymmetrical 0.432 0.0648 Inf 0.312 0.561
## Pointing.crossed 0.573 0.0644 Inf 0.445 0.692
## Pointing.proximal 0.619 0.0630 Inf 0.490 0.733
## Touching.asymmetrical 0.866 0.0413 Inf 0.763 0.928
## Touching.crossed 0.837 0.0454 Inf 0.727 0.908
## Touching.proximal 0.865 0.0415 Inf 0.762 0.928
##
## Results are averaged over the levels of: sex, corr.side
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale
##
## $contrasts
## contrast odds.ratio SE df null
## Control / Gazing.asymmetrical 1.045 0.3840 Inf 1
## Control / Gazing.central 0.612 0.2280 Inf 1
## Control / Pointing.asymmetrical 1.492 0.5510 Inf 1
## Control / Pointing.crossed 0.848 0.3120 Inf 1
## Control / Pointing.proximal 0.699 0.2590 Inf 1
## Control / Touching.asymmetrical 0.176 0.0771 Inf 1
## Control / Touching.crossed 0.222 0.0930 Inf 1
## Control / Touching.proximal 0.177 0.0776 Inf 1
## Gazing.asymmetrical / Gazing.central 0.586 0.2180 Inf 1
## Gazing.asymmetrical / Pointing.asymmetrical 1.428 0.5280 Inf 1
## Gazing.asymmetrical / Pointing.crossed 0.811 0.2980 Inf 1
## Gazing.asymmetrical / Pointing.proximal 0.669 0.2480 Inf 1
## Gazing.asymmetrical / Touching.asymmetrical 0.168 0.0738 Inf 1
## Gazing.asymmetrical / Touching.crossed 0.212 0.0890 Inf 1
## Gazing.asymmetrical / Touching.proximal 0.170 0.0743 Inf 1
## Gazing.central / Pointing.asymmetrical 2.439 0.9150 Inf 1
## Gazing.central / Pointing.crossed 1.385 0.5170 Inf 1
## Gazing.central / Pointing.proximal 1.142 0.4290 Inf 1
## Gazing.central / Touching.asymmetrical 0.288 0.1270 Inf 1
## Gazing.central / Touching.crossed 0.362 0.1530 Inf 1
## Gazing.central / Touching.proximal 0.290 0.1280 Inf 1
## Pointing.asymmetrical / Pointing.crossed 0.568 0.2100 Inf 1
## Pointing.asymmetrical / Pointing.proximal 0.468 0.1750 Inf 1
## Pointing.asymmetrical / Touching.asymmetrical 0.118 0.0519 Inf 1
## Pointing.asymmetrical / Touching.crossed 0.149 0.0627 Inf 1
## Pointing.asymmetrical / Touching.proximal 0.119 0.0523 Inf 1
## Pointing.crossed / Pointing.proximal 0.825 0.3060 Inf 1
## Pointing.crossed / Touching.asymmetrical 0.208 0.0909 Inf 1
## Pointing.crossed / Touching.crossed 0.262 0.1100 Inf 1
## Pointing.crossed / Touching.proximal 0.209 0.0915 Inf 1
## Pointing.proximal / Touching.asymmetrical 0.252 0.1110 Inf 1
## Pointing.proximal / Touching.crossed 0.317 0.1340 Inf 1
## Pointing.proximal / Touching.proximal 0.254 0.1120 Inf 1
## Touching.asymmetrical / Touching.crossed 1.260 0.6050 Inf 1
## Touching.asymmetrical / Touching.proximal 1.007 0.5000 Inf 1
## Touching.crossed / Touching.proximal 0.799 0.3840 Inf 1
## z.ratio p.value
## 0.119 1.0000
## -1.318 0.9264
## 1.085 0.9766
## -0.450 1.0000
## -0.966 0.9888
## -3.965 0.0024
## -3.590 0.0100
## -3.952 0.0025
## -1.434 0.8847
## 0.965 0.9889
## -0.569 0.9997
## -1.084 0.9766
## -4.066 0.0016
## -3.695 0.0068
## -4.052 0.0017
## 2.377 0.2966
## 0.873 0.9943
## 0.354 1.0000
## -2.821 0.1092
## -2.397 0.2853
## -2.806 0.1134
## -1.529 0.8423
## -2.034 0.5194
## -4.853 <.0001
## -4.519 0.0002
## -4.840 <.0001
## -0.519 0.9999
## -3.589 0.0100
## -3.197 0.0375
## -3.574 0.0106
## -3.134 0.0454
## -2.723 0.1397
## -3.119 0.0475
## 0.481 0.9999
## 0.014 1.0000
## -0.467 0.9999
##
## Results are averaged over the levels of: sex, corr.side
## P value adjustment: tukey method for comparing a family of 9 estimates
## Tests are performed on the log odds ratio scale
# Calcular % correcto por sujeto y condición
prop <- xdata %>%
group_by(subject, condition) %>%
summarise(
prop_correct = mean(correct == 1) * 100,
.groups = "drop"
)
# Extraer las condiciones disponibles
condiciones <- unique(prop$condition)
condiciones_experimentales <- setdiff(condiciones, "Control")
# Valores de Control (todos los sujetos)
control_vals <- prop %>%
filter(condition == "Control") %>%
pull(prop_correct)
# Inicialitzar una nova taula de resultats per al test bilateral
wilcox_vs_50_2 <- data.frame(
Condition = character(),
W = numeric(),
p_value = numeric(),
Significant = logical(),
stringsAsFactors = FALSE
)
cat("### Comparació de cada condició vs 50% (Wilcoxon signed-rank test - TWO.SIDED) ###\n")
## ### Comparació de cada condició vs 50% (Wilcoxon signed-rank test - TWO.SIDED) ###
for (cond in unique(prop$condition)) {
cat("\nCondició:", cond, "\n")
# Valors de la condició
cond_vals <- prop %>%
filter(condition == cond) %>%
pull(prop_correct)
# Convertir a proporció (0.5 = 50%)
cond_vals_prop <- cond_vals / 100
# Test de Wilcoxon contra 0.5 (bilateral)
test <- suppressWarnings(
wilcox.test(cond_vals_prop,
mu = 0.5,
alternative = "two.sided",
exact = NULL,
correct = TRUE)
)
print(test)
# Guardar resultats
wilcox_vs_50_2 <- rbind(wilcox_vs_50_2, data.frame(
Condition = cond,
W = test$statistic,
p_value = test$p.value,
Significant = test$p.value < 0.05
))
}
##
## Condició: Control
##
## Wilcoxon signed rank test with continuity correction
##
## data: cond_vals_prop
## V = 6, p-value = 0.8501
## alternative hypothesis: true location is not equal to 0.5
##
##
## Condició: Gazing.asymmetrical
##
## Wilcoxon signed rank test with continuity correction
##
## data: cond_vals_prop
## V = 13.5, p-value = 1
## alternative hypothesis: true location is not equal to 0.5
##
##
## Condició: Gazing.central
##
## Wilcoxon signed rank test with continuity correction
##
## data: cond_vals_prop
## V = 15, p-value = 0.05334
## alternative hypothesis: true location is not equal to 0.5
##
##
## Condició: Pointing.asymmetrical
##
## Wilcoxon signed rank test with continuity correction
##
## data: cond_vals_prop
## V = 4, p-value = 0.1983
## alternative hypothesis: true location is not equal to 0.5
##
##
## Condició: Pointing.crossed
##
## Wilcoxon signed rank test with continuity correction
##
## data: cond_vals_prop
## V = 24, p-value = 0.407
## alternative hypothesis: true location is not equal to 0.5
##
##
## Condició: Pointing.proximal
##
## Wilcoxon signed rank test with continuity correction
##
## data: cond_vals_prop
## V = 21, p-value = 0.265
## alternative hypothesis: true location is not equal to 0.5
##
##
## Condició: Touching.asymmetrical
##
## Wilcoxon signed rank test with continuity correction
##
## data: cond_vals_prop
## V = 36, p-value = 0.01333
## alternative hypothesis: true location is not equal to 0.5
##
##
## Condició: Touching.crossed
##
## Wilcoxon signed rank test with continuity correction
##
## data: cond_vals_prop
## V = 28, p-value = 0.02107
## alternative hypothesis: true location is not equal to 0.5
##
##
## Condició: Touching.proximal
##
## Wilcoxon signed rank test with continuity correction
##
## data: cond_vals_prop
## V = 36, p-value = 0.01356
## alternative hypothesis: true location is not equal to 0.5
# Mostrar la taula resum
cat("\n### Resum: Test bilateral contra 50% per condició ###\n")
##
## ### Resum: Test bilateral contra 50% per condició ###
print(wilcox_vs_50_2)
## Condition W p_value Significant
## V Control 6.0 0.85010674 FALSE
## V1 Gazing.asymmetrical 13.5 1.00000000 FALSE
## V2 Gazing.central 15.0 0.05333685 FALSE
## V3 Pointing.asymmetrical 4.0 0.19825640 FALSE
## V4 Pointing.crossed 24.0 0.40701603 FALSE
## V5 Pointing.proximal 21.0 0.26496151 FALSE
## V6 Touching.asymmetrical 36.0 0.01332833 TRUE
## V7 Touching.crossed 28.0 0.02106789 TRUE
## V8 Touching.proximal 36.0 0.01356012 TRUE
# Calcular media y SEM para cada condición
summary_stats <- prop %>%
group_by(condition) %>%
summarise(
mean_prop = mean(prop_correct),
sem_prop = sd(prop_correct) / sqrt(n()),
.groups = "drop"
)
# Añadir columna de significancia a summary_stats
summary_stats <- summary_stats %>%
left_join(wilcox_vs_50_2, by = c("condition" = "Condition"))
# Función para crear etiquetas de significancia
summary_stats <- summary_stats %>%
mutate(sig_label = ifelse(Significant, "*", ""))
# Gráfico
ggplot(summary_stats, aes(x = condition, y = mean_prop, fill = condition)) +
geom_bar(stat = "identity", color = "black", width = 0.7) +
geom_errorbar(aes(ymin = mean_prop - sem_prop, ymax = mean_prop + sem_prop), width = 0.2) +
geom_text(aes(label = sig_label, y = mean_prop + sem_prop + 3), size = 8, color = "red") +
geom_hline(yintercept = 50, linetype = "dashed") +
labs(
title = "% de aciertos por condición",
subtitle = "Estrellas (*) indican significancia frente a Control (Wilcoxon, p < 0.05)",
x = "Condición",
y = "% correcto",
fill = "Condición"
) +
theme_minimal() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "none"
)
ggplot(summary_stats, aes(x = condition, y = mean_prop, fill = condition)) +
geom_bar(stat = "identity", color = "black", width = 0.7) +
geom_errorbar(aes(ymin = mean_prop - sem_prop, ymax = mean_prop + sem_prop), width = 0.2) +
geom_text(aes(label = sig_label, y = mean_prop + sem_prop + 3), size = 8, color = "red") +
geom_hline(yintercept = 50, linetype = "dashed") +
labs(
title = "% de encerts per condició",
subtitle = "Estrelles (*) idiquen significancia davant el 50% d'encert (Wilcoxon, p < 0.05)",
x = "",
y = "% correcte",
fill = "Condició"
) +
scale_fill_discrete(labels = c("Control",
"Gazing asymmetrical",
"Gazing central",
"Pointing asymmetrical",
"Pointing crossed",
"Pointing proximal",
"Touching asymmetrical",
"Touching crossed",
"Touching proximal"
)) +
theme_minimal() +
theme(
legend.position = "right",
axis.text.x = element_blank()
)
# GRÁFICO GLOBAL POR CONDICIÓN
ggplot(prop, aes(x = condition, y = prop_correct, fill = factor(condition))) +
geom_bar(stat = "identity", position = "dodge") +
labs(
title = "Rendimient individual per condició",
x = "",
y = "% intents correctes",
fill = "Condició"
) +
scale_fill_discrete(labels = c("Control",
"Gazing asymmetrical",
"Gazing central",
"Pointing asymmetrical",
"Pointing crossed",
"Pointing proximal",
"Touching asymmetrical",
"Touching crossed",
"Touching proximal"
)) +
geom_hline(yintercept = 50, linetype = "dotted", color = "black", size = 1) +
theme_bw() +
theme(
legend.position = "right",
axis.text.x = element_blank()
) +
facet_wrap(~subject)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
ggplot(xdata, aes(x = condition, y = correct)) +
geom_violin(aes(fill = factor(condition)), trim = FALSE) +
labs(x = "", y = "1 = Correcte, 0 = Incorrecte", fill = "Condició") +
scale_fill_discrete(labels = c("Control",
"Gazing asymmetrical",
"Gazing central",
"Pointing asymmetrical",
"Pointing crossed",
"Pointing proximal",
"Touching asymmetrical",
"Touching crossed",
"Touching proximal"
)) +
theme_bw() +
theme(legend.position = "right", axis.text.x = element_blank() )
ggplot(prop, aes(x = condition, y = prop_correct, fill = condition)) +
geom_boxplot(alpha = 0.7, outlier.shape = 16, outlier.size = 2) +
labs(
title = "Dispersió dels percentatges d'encert entre individus",
x = "",
y = "Percentatge d'encerts (%)",
fill = "Condició"
) +
scale_fill_brewer(
palette = "Set2",
labels = c("Control",
"Gazing asymmetrical",
"Gazing central",
"Pointing asymmetrical",
"Pointing crossed",
"Pointing proximal",
"Touching asymmetrical",
"Touching crossed",
"Touching proximal")
) +
geom_hline(yintercept = 50, linetype = "dashed", color = "black") +
theme_minimal() +
theme(
legend.position = "right",
axis.text.x = element_blank()
)
## Warning in RColorBrewer::brewer.pal(n, pal): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
