Cargar librerías necesarias Cargar la base de datos y adaptarla a dos formatos: Long para análisis de medidas repetidas y wide para análisis puntuales.

wide=read_xlsx("PEERJ_HOCKEY_DATA.xlsx",1) %>% 
  mutate(CODIGO=as.factor(as.integer(CODIGO)),
         Position_play=factor(as.integer(Position_play),levels=c(4,3,2,1),labels=c( "Forward", "Center", "Defense", "Goalie" )),
         Outcome=factor(as.integer(Outcome), levels=c(0,1),labels=c("Loss", "Win"))) %>% 
  rename(CORT_0=CORT_Awakening,CORT_1=CORT_30,CORT_2=CORT_PRE,CORT_3=CORT_POST) %>%
  mutate(CAR=CORT_1-CORT_0) %>% 
  select(starts_with("CORT_"), CAR, CODIGO,Position_play,Outcome,
         Hard_day:Salience_CAR) 
  

long=wide %>% pivot_longer(-c(CODIGO:Salience_CAR,CAR), values_to="CORT",names_to="Time") %>%
  select(CORT,Time,everything()) %>% 
  mutate(Time=factor(Time,levels=c("CORT_0","CORT_1","CORT_2","CORT_3"),
                     labels=c("C0","C1","C2","C3")))
#Exploramos las primeras lineas para ver si estamos conformes
long %>% slice(1:10) %>% knitr::kable()
CORT Time CAR CODIGO Position_play Outcome Hard_day Vigor_reach_challenge Anxiety_perception Nerviousness_perception Inspires_teammates Sacrifice_myself Compite_with_passion My_Words_motivate_teammates Be_defeated_not_failure SAM_Affective_valence SAM_AROUSAL SAM_DOMINANCE Cognitive_Axiety Somatic_Axiety Self_confidence Salience_CAR
2.30 C0 1.60 1 Defense Loss 7 8 4 5 7 8 8 7 5 6 6 7 31 25 30 69.56522
3.90 C1 1.60 1 Defense Loss 7 8 4 5 7 8 8 7 5 6 6 7 31 25 30 69.56522
4.10 C2 1.60 1 Defense Loss 7 8 4 5 7 8 8 7 5 6 6 7 31 25 30 69.56522
17.10 C3 1.60 1 Defense Loss 7 8 4 5 7 8 8 7 5 6 6 7 31 25 30 69.56522
3.81 C0 8.09 2 Center Loss 9 9 2 7 6 9 8 6 7 9 5 7 32 28 36 212.33596
11.90 C1 8.09 2 Center Loss 9 9 2 7 6 9 8 6 7 9 5 7 32 28 36 212.33596
5.52 C2 8.09 2 Center Loss 9 9 2 7 6 9 8 6 7 9 5 7 32 28 36 212.33596
9.78 C3 8.09 2 Center Loss 9 9 2 7 6 9 8 6 7 9 5 7 32 28 36 212.33596
6.50 C0 2.24 3 Forward Loss 7 8 3 5 7 9 8 7 8 8 7 7 23 19 29 34.46154
8.74 C1 2.24 3 Forward Loss 7 8 3 5 7 9 8 7 8 8 7 7 23 19 29 34.46154

1. Two-way Repeated Measures ANOVA (modelo principal del estudio)

Objetivo: analizar la evolución del cortisol a lo largo del día y las diferencias entre ganadoras y perdedoras.

Estructura del modelo:

Requisitos:

normalidad_test <- long %>%
  group_by(Time, Outcome) %>%
  shapiro_test(CORT)
normalidad_test %>% knitr::kable()
Time Outcome variable statistic p
C0 Loss CORT 0.8340646 0.0004532
C0 Win CORT 0.9460418 0.2220482
C1 Loss CORT 0.8193423 0.0002387
C1 Win CORT 0.9679101 0.6157329
C2 Loss CORT 0.9175515 0.0301867
C2 Win CORT 0.8216420 0.0006799
C3 Loss CORT 0.9792229 0.8313972
C3 Win CORT 0.9403148 0.1657810

Tenemos falta de normalidad en algunas combinaciones de TimexOutcome. Probamos a ver si una transformación logaritmica lo resuelve:

long=long %>% mutate(logCORT=log(CORT))


normalidad_test <- long %>%
  group_by(Time, Outcome) %>%
  shapiro_test(logCORT) 

normalidad_test %>% knitr::kable()
Time Outcome variable statistic p
C0 Loss logCORT 0.9502522 0.2010404
C0 Win logCORT 0.9611487 0.4619368
C1 Loss logCORT 0.9570077 0.2954854
C1 Win logCORT 0.9779618 0.8555333
C2 Loss logCORT 0.9752457 0.7254474
C2 Win logCORT 0.9806345 0.9072595
C3 Loss logCORT 0.8252102 0.0003073
C3 Win logCORT 0.9575769 0.3917023

Pues solo tenemos un grupo para el que no tenemos normalidad tras la transformación (C3-Loss). Hay que tener en cuenta que estas pruebas son resistentes a cierta desviación de la normalidad y yo creo que bastará con mencionarlo. Vamos a proceder con la ANOVA de medidas repetidas sobre los datos log-transformados.

homogeneidad_test <- long %>%
  group_by(Time) %>%
  levene_test(logCORT ~ Outcome)
homogeneidad_test %>% knitr::kable()
Time df1 df2 statistic p
C0 1 50 0.6000159 0.4422197
C1 1 50 0.7137658 0.4022211
C2 1 50 0.0001854 0.9891913
C3 1 50 0.1004199 0.7526457

No tenemos problemas por aquí. Más razón para no preocuparnos de más por la normalidad en alguna categoría.

aov_mixto <- long %>%
  # El diseño: logCORT como dependiente, CODIGO como sujeto,
  # Time como factor within, Outcome como factor between.
  anova_test(
    logCORT ~ Outcome * Time,
    wid = CODIGO,
    within = Time,
    between = Outcome
  )


get_anova_table(aov_mixto) %>% knitr::kable()
Effect DFn DFd F p p<.05 ges
Outcome 1 200 11.675 0.000767 * 0.055
Time 3 200 18.128 0.000000 * 0.214
Outcome:Time 3 200 11.053 0.000001 * 0.142

Lo importante es la interacción Outcome:Time que tiene un tamaño de efecto medio según Cohen (ges=generalized eta squared=0.142). Todo lo mayor a 0.13 se considera efecto medio, pero revisa si en deporte se usan otros valores. Algunos autores a eso lo llaman grande. es un poco arbitrario el corte. En resumen: La evolución del cortisol a lo largo del día, no es lo mismo para ganadoras que perdedoras. Supongo que harán ploff las perdedoras cuando veamos el gráfico.

Vamos a hacer un análisis extra con otra metodología que a mí me gusta más aunque no sean la que te han pedido. Básicamente por que me entero mejor yo:

lmer_model <- lmer(
  logCORT ~ Outcome * Time + (1 | CODIGO),
  data = long
)





tabla_lmm_apa <- modelsummary(
  lmer_model,
  stars = TRUE,          
  statistic = "conf.int", 

  output = "tinytable",  
  fmt = fmt_decimal(digits = 3), 

  gof_map = "default",
  title = "Modelo LMM: Efectos del Resultado y el Tiempo en log(Cortisol)"
)




tabla_lmm_apa
Modelo LMM: Efectos del Resultado y el Tiempo en log(Cortisol)
(1)
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
(Intercept) 1.516***
[1.339, 1.694]
OutcomeWin 0.080
[-0.181, 0.341]
TimeC1 0.591***
[0.403, 0.778]
TimeC2 0.007
[-0.180, 0.194]
TimeC3 0.947***
[0.760, 1.134]
OutcomeWin × TimeC1 -0.283*
[-0.559, -0.008]
OutcomeWin × TimeC2 -0.005
[-0.281, 0.270]
OutcomeWin × TimeC3 -0.936***
[-1.211, -0.661]
SD (Intercept CODIGO) 0.317
SD (Observations) 0.355

El gráfico de la interacción sin transformar es este:

grafico_interaccion <- long %>%
  group_by(Time, Outcome) %>%
  summarise(
    mean_CORT = mean(CORT),
    sd_CORT = sd(CORT),
    n = n(),
    se_CORT = sd_CORT / sqrt(n)
  ) %>%
  ggplot(aes(x = Time, y = mean_CORT, group = Outcome, color = Outcome)) +
  geom_point(size = 3) +
  geom_line(linewidth = 1) +
  geom_errorbar(aes(ymin = mean_CORT - se_CORT, ymax = mean_CORT + se_CORT), width = 0.2) +
  labs(
    title = "Evolución del Cortisol (CORT) por Resultado del Partido",
    x = "Tiempo de Medición (Time)",
    y = "Cortisol Promedio (CORT)"
  ) +
  theme_minimal()
## `summarise()` has grouped output by 'Time'. You can override using the
## `.groups` argument.
print(grafico_interaccion)

Anda, lo que yo esperaba como un “ploff” de las perdedoras, es más bien lo contrario.

En realidad usamos log(CORT), y lo correcto sería hacer cálculos en log(CORT) y transformar el gráfico a variables originales. Esto es más correcto, usando medias geométricas:

ci_log_scale <- long %>%
  group_by(Time, Outcome) %>%
  get_summary_stats(logCORT, type = "mean_se") %>% 
  
  mutate(
    t_crit = qt(0.975, df = n - 1),
    
    mean_log = mean,
    ci_lower_log = mean_log - (t_crit * se),
    ci_upper_log = mean_log + (t_crit * se)
  ) %>%
  
  mutate(
    mean_CORT_geom = exp(mean_log),
    ci_lower_CORT = exp(ci_lower_log),
    ci_upper_CORT = exp(ci_upper_log)
  )

grafico_interaccion_ajustado <- ci_log_scale %>%
  ggplot(aes(x = Time, y = mean_CORT_geom, group = Outcome, color = Outcome)) +
  
  # La media geométrica y los límites transformados
  geom_point(size = 3) +
  geom_line(linewidth = 1) +
  geom_errorbar(aes(ymin = ci_lower_CORT, ymax = ci_upper_CORT), width = 0.2) +
  
  labs(
    title = "Evolución de Cortisol (Media Geométrica) por Resultado del Partido",
    caption = "Las barras de error representan el 95% CI calculado sobre log(CORT) y transformado.",
    x = "Tiempo de Medición (Time)",
    y = "Media Geométrica del Cortisol (CORT)"
  ) +
  theme_minimal()

print(grafico_interaccion_ajustado)

Algo curioso es que parece que en “C1” algo ayuda a predecir la derrota. Me desvío un poco para probar la idea:

coef_rename <- c(
  "(Intercept)" = "Constante (Intersección)",
  "I(log(CORT_1) - log(CORT_0))" = "Cambio Log-Cortisol (T1 - T0)"
)

testeandoIdea <- glm(I(Outcome=="Loss") ~ I(log(CORT_1)-log(CORT_0)), 
                 data=wide, 
                 family=binomial)


modelsummary(
  testeandoIdea,
  exponentiate = TRUE,      # ¡Importante! Muestra Odds Ratios en lugar de log-odds
  coef_map = coef_rename,   # Aplica los nombres limpios
  stars = TRUE,             # Añade los asteriscos de significancia (p < .05, etc.)
  statistic = "conf.int",   # Muestra el Intervalo de Confianza en lugar del error estándar
  title = "Odds Ratios: Predicción de Pérdida (Loss) según cambio en Cortisol",
  gof_map = c("nobs", "aic", "logLik"), # Selecciona qué estadísticos mostrar al final
  output = "gt"             # Formato de alta calidad para Word/HTML
)
Odds Ratios: Predicción de Pérdida (Loss) según cambio en Cortisol
(1)
Constante (Intersección) 0.210*
[0.050, 0.693]
Cambio Log-Cortisol (T1 - T0) 47.431**
[4.437, 996.425]
Num.Obs. 52
AIC 122.1
Log.Lik. -59.053
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 
predicciones <- predict(testeandoIdea, type = "response")

# 2. Crear el objeto ROC
# Nota: I(Outcome=="Loss") genera un vector lógico (TRUE/FALSE) que pROC entiende bien
curva_roc <- roc(wide$Outcome == "Loss", predicciones)
## Setting levels: control = FALSE, case = TRUE
## Setting direction: controls < cases
# 3. Dibujar la curva
plot(curva_roc, 
     main = "Curva ROC del Modelo de Cortisol",
     col = "#2196F3", 
     lwd = 4, 
     print.auc = TRUE,           # Imprime el valor del AUC en el gráfico
     legacy.axes = TRUE,         # Eje X como 1-Especificidad (estándar de publicación)
     xlab = "1 - Especificidad (Falsos Positivos)",
     ylab = "Sensibilidad (Verdaderos Positivos)")

# 4. Ver el AUC exacto y su intervalo de confianza
print(auc(curva_roc))
## Area under the curve: 0.7857
print(ci.auc(curva_roc))
## 95% CI: 0.655-0.9164 (DeLong)

Pues sí que se consigue predecir bantante bien las derrotas mirando el cambio entre C0 y C1. Estaría bien tener más equipos.

2. Fiabilidad de medidas repetidas

Intraclass Correlation Coefficient (ICC) Modelo: two-way mixed, absolute agreement

Evaluar la fiabilidad intra-sujeto de C0, C1, C2 y C3.

Reportar:

ICC 95% CI Interpretación (p. ej., buena, moderada, excelente)

Entiendo que quieren hacer el análisis con 4 jueces que miden la misma cosa, o lo que es lo mismo, miden a las personas en 4 tiempos muy bien elegidos. No son tiempos cualquiera y estudian la estabilidad.

wide=wide %>% mutate(logCORT_0=log(CORT_0),
                      logCORT_1=log(CORT_1),
                      logCORT_2=log(CORT_2),
                      logCORT_3=log(CORT_3))

icc_results <- ICC(wide %>% select(starts_with("logCORT_")))

# 4. Buscar en el output: "Two-way mixed, absolute agreement"
# Normalmente corresponde al código ICC3 o ICC3k en la tabla de salida

icc_results$results %>% knitr::kable()
type ICC F df1 df2 p lower bound upper bound
Single_raters_absolute ICC1 0.2444020 2.293820 51 156 4.96e-05 0.1110287 0.4026436
Single_random_raters ICC2 0.2855229 3.298939 51 153 0.00e+00 0.1328142 0.4532458
Single_fixed_raters ICC3 0.3649724 3.298939 51 153 0.00e+00 0.2238735 0.5192460
Average_raters_absolute ICC1k 0.5640461 2.293820 51 156 4.96e-05 0.3331478 0.7294497
Average_random_raters ICC2k 0.6151626 3.298939 51 153 0.00e+00 0.3798918 0.7682987
Average_fixed_raters ICC3k 0.6968722 3.298939 51 153 0.00e+00 0.5357041 0.8120394

Ahí están todos los posibles ICC, pero no os interesan todos. Según os solicitan, quieren el ICC3, que sale un valor bajo.

De hecho dudo que os interese tampoco ese ya que hemos visto antes que las mediciones de ICC no son estables (C3 alto parece asociado a la derrota, que es algo que ocurre tras C0, C1 y C2). Lo pongo por que lo pedís.

Se puede interpretar de este modo si queréis.

La fiabilidad intra-sujeto para logCORT a través de los cuatro puntos temporales fue moderada/pobre (ICC(3,1) = 0.37, IC 95% [0.22, 0.52]). Este valor refleja la variabilidad biológica inherente al cortisol y justifica el uso de modelos de efectos mixtos para controlar la varianza residual no explicada por las diferencias entre sujetos.

O bien El 37% de la variabilidad que se observa en los niveles de cortisol se debe a diferencias reales entre los sujetos (características estables de la persona) mientras que el 63% restante se debe a fluctuaciones intra-sujeto por (y error de medición, claro).

3. Comparaciones puntuales del CAR (C1 – C0)

Comparar CAR entre ganadoras y perdedoras mediante:

t-test independiente (si log-transformado es normal) o Mann–Whitney U (si se mantiene la no normalidad) Reportar:

Media ± SD d de Cohen + 95% CI Valor p

Empecemos por la normalidad de CAR

#PAra CAR como diferencia
wide %>%  shapiro_test(CAR)
## # A tibble: 1 × 3
##   variable statistic        p
##   <chr>        <dbl>    <dbl>
## 1 CAR          0.902 0.000415

Los datos no son normales, pero tampoco es que sean ofensivamente anormales

wide %>% ggplot(aes(x=CORT_0,y=CORT_1))+
  geom_point()

Siguiendo las instrucciones vamos a seguir por la vía no paramétrica con CAR

Comparar CAR entre ganadoras y perdedoras mediante prueba de Mann-Whitney

mann_whitney_results <- wide %>% 
  wilcox_test(CAR ~ Outcome) 
mann_whitney_results %>% knitr::kable()
.y. group1 group2 n1 n2 statistic p
CAR Loss Win 28 24 532 0.000332 
add_stat_effect_size <- function(data, variable, by, ...) {
  test_result <- effectsize::rank_biserial(data[[variable]] ~ data[[by]])
  # Devolvemos el valor con 2 decimales
  tibble::tibble(effect_size = sprintf("%.2f", test_result$r_rank_biserial))
}

tabla_publicacion <- wide %>%
  select(Outcome, CAR) %>% # Selecciona tus variables
  tbl_summary(
    by = Outcome,
    statistic = list(all_continuous() ~ "{median} ({p25}, {p75})"), # Formato Mediana (IQR)
    digits = all_continuous() ~ 2,
    label = list(CAR ~ "CAR=C1-C0")
  ) %>%
  add_p(test = list(all_continuous() ~ "wilcox.test")) %>% 
  add_stat(fns = all_continuous() ~ add_stat_effect_size) %>% 
  modify_header(
    label = "**Variable**",
    effect_size = "**Size effect ($r$)**",
    p.value = "**p-valor**"
  ) %>%
  bold_labels() %>%
  italicize_levels()
## The following warnings were returned during `modify_header()`:
## ! For variable `CAR` (`Outcome`) and "estimate", "statistic", "p.value",
##   "conf.low", and "conf.high" statistics: cannot compute exact p-value with
##   ties
## ! For variable `CAR` (`Outcome`) and "estimate", "statistic", "p.value",
##   "conf.low", and "conf.high" statistics: cannot compute exact confidence
##   intervals with ties
tabla_publicacion
Variable Loss
N = 281		 Win
N = 241		 p-valor2 Size effect ( rr)
CAR=C1-C0 3.78 (2.50, 4.89) 1.84 (1.04, 2.66) <0.001 0.58
1 Median (Q1, Q3)
2 Wilcoxon rank sum test

El tamaño de efecto no suele ser d de Cohen como pides sino la r de Rosenthal/Rank-biserial, que de acuerdo a los criterios de Cohen se considera “LARGE” para valores superiores a 0.5

Si quieres añadir un IC95% es este:

effectsize::rank_biserial(CAR ~ Outcome, data = wide)
## r (rank biserial) |       95% CI
## --------------------------------
## 0.58              | [0.34, 0.76]

Lo mismo se puede mostrar de forma alternativa gráficamente, todo en uno:

# Crear el gráfico
ggbetweenstats(
  data = wide,
  x = Outcome,
  y = CAR,
  type = "nonparametric",      
  conf.level = 0.95,           
  plot.type = "boxviolin",     
  mean.plotting = TRUE,        
  centrality.plotting = TRUE,  
  title = "Differences in CAR by Outcome",
  xlab = "Outcome",
  ylab = "CAR",
  messages = FALSE,            
  package = "ggsci",           # Paleta de colores científica (Nature Publishing Group)
  palette = "nrc_npg"
) +  theme(
    plot.title = element_text(size = 14, face = "bold"),
    axis.title = element_text(size = 12)
  )

Otra posibilidad

Podriamos tener una versión alternativa de CAR como CAR_Ratio=C1/C0

wide=wide %>% mutate(CAR_Ratio=CORT_1/CORT_0)
wide$CAR_Ratio %>% hist()

Esa Ratio está pidiendo un logaritmo:

wide$logCAR_Ratio=log(wide$CAR_Ratio)
wide$logCAR_Ratio %>% hist()

wide %>%  shapiro_test(logCAR_Ratio)
## # A tibble: 1 × 3
##   variable     statistic      p
##   <chr>            <dbl>  <dbl>
## 1 logCAR_Ratio     0.959 0.0726

Prueba superada. En caso de que os interese.

logCAR_Ratio = log(C1/C0)= log(C1)-log(C0)

ggbetweenstats(
  data = wide,
  x = Outcome,
  y = logCAR_Ratio,
  type = "parametric",      
  conf.level = 0.95,           
  plot.type = "boxviolin",     
  mean.plotting = TRUE,        
  centrality.plotting = TRUE,  
  title = "Differences in log(C1/C0) by Outcome",
  xlab = "Outcome",
  ylab = "log(C1/C0)",
  messages = FALSE,            
  package = "ggsci",           
  palette = "nrc_npg"
) +
  theme(
    plot.title = element_text(size = 14, face = "bold"),
    axis.title = element_text(size = 12)
  )

4. Porcentaje de cambio (pre → post)

Variable: Δ% cortisol = [(post – pre) / pre] × 100

Dado que esta variable rara vez es normal:

Usar Mann–Whitney U para comparar Win vs Loss. Reportar r o d como tamaño del efecto + CI.

wide=wide %>% mutate(Delta_CORT=((CORT_3-CORT_2)/CORT_2)*100)
hist(wide$Delta_CORT)

tabla_publicacion <- wide %>%
  select(Outcome, Delta_CORT) %>% # Selecciona tus variables
  tbl_summary(
    by = Outcome,
    statistic = list(all_continuous() ~ "{median} ({p25}, {p75})"), # Formato Mediana (IQR)
    digits = all_continuous() ~ 2,
    label = list(Delta_CORT ~ "Δ% cortisol")
  ) %>%
  add_p(test = list(all_continuous() ~ "wilcox.test")) %>% 
  add_stat(fns = all_continuous() ~ add_stat_effect_size) %>% 
  modify_header(
    label = "**Variable**",
    effect_size = "**Size effect ($r$)**",
    p.value = "**p-valor**"
  ) %>%
  bold_labels() %>%
  italicize_levels()

tabla_publicacion
Variable Loss
N = 281		 Win
N = 241		 p-valor2 Size effect ( rr)
Δ% cortisol 128.10 (95.43, 244.69) -12.60 (-32.90, 56.33) <0.001 0.74
1 Median (Q1, Q3)
2 Wilcoxon rank sum exact test

El effect size en forma de IC95%. Sale enorme el efecto en todo el IC95%.

effectsize::rank_biserial(Delta_CORT ~ Outcome, data = wide)
## r (rank biserial) |       95% CI
## --------------------------------
## 0.74              | [0.57, 0.86]

Lo mismo se puede mostrar de forma alternativa gráficamente, todo en uno:

# Crear el gráfico
ggbetweenstats(
  data = wide,
  x = Outcome,
  y = Delta_CORT,
  type = "nonparametric",      
  conf.level = 0.95,           
  plot.type = "boxviolin",     
  mean.plotting = TRUE,        
  centrality.plotting = TRUE,  
  title = "Differences in Δ% cortisol by Outcome",
  xlab = "Outcome",
  ylab = "Δ% cortisol",
  messages = FALSE,            
  package = "ggsci",           
  palette = "nrc_npg"
) +  theme(
    plot.title = element_text(size = 14, face = "bold"),
    axis.title = element_text(size = 12)
  )

5. Comparaciones por posición de juego

Prueba adecuada:

One-way ANOVA con Position (4 niveles) como factor Post-hoc Bonferroni

Las variables dependientes son C0 a C3, y también CAR=C1-C0 y posiblemente logCAR = log(C1/C0)

Reportar:

F, p η²p Comparaciones múltiples ajustadas

Según me dicen por correo la variable dependiente es cortisol, aunque tengo 4 de esas. Voy a hacerlo para las 4 mediciones de cortisol por separado. En realidad de conversaciones posteriores vimos que realmente lo que interesaba como variable dependientes era CAR (C1-C0) o bien logCAR=log(C1/C0). (Ojo que aunque lo llame logCAR no es logCar_Ratio, es el log de la ratio C1 vs C0, y solo trabajo con ella por lo normal que es.)

Para C_0

 wide %>% 
  anova_test(logCORT_0 ~ Position_play, detailed = TRUE)
## ANOVA Table (type II tests)
## 
##          Effect   SSn    SSd DFn DFd     F     p p<.05  ges
## 1 Position_play 0.642 10.128   3  48 1.015 0.394       0.06

Nada

Para C_1

 wide %>% 
  anova_test(logCORT_1 ~ Position_play, detailed = TRUE)
## ANOVA Table (type II tests)
## 
##          Effect   SSn   SSd DFn DFd     F     p p<.05   ges
## 1 Position_play 0.576 8.869   3  48 1.038 0.384       0.061

Nada

Para C_2

 wide %>% 
  anova_test(logCORT_2 ~ Position_play, detailed = TRUE)
## ANOVA Table (type II tests)
## 
##          Effect   SSn   SSd DFn DFd     F     p p<.05   ges
## 1 Position_play 2.262 13.66   3  48 2.649 0.059       0.142

Casi sale significativo. No hay pruebas múltiples que hacer.

Para C_3

 wide %>% 
  anova_test(logCORT_3 ~ Position_play, detailed = TRUE)
## ANOVA Table (type II tests)
## 
##          Effect   SSn    SSd DFn DFd     F     p p<.05   ges
## 1 Position_play 1.224 18.115   3  48 1.081 0.366       0.063

Nada

Primera variable de Cortisol dependiente de interés: CAR

wide %>% 
  anova_test(CAR ~ Position_play, detailed = TRUE)
## ANOVA Table (type II tests)
## 
##          Effect    SSn     SSd DFn DFd     F    p p<.05   ges
## 1 Position_play 15.882 270.935   3  48 0.938 0.43       0.055

Gráficamente no va a tener interés, pero lo mostramps de todos modos.

ggbetweenstats(
  data = wide,
  x = Position_play,
  y = CAR,
  type = "parametric",      
  conf.level = 0.95,           
  plot.type = "boxviolin",     
  mean.plotting = TRUE,        
  centrality.plotting = TRUE,  
  title = "Differences in CAR by Position_play",
  xlab = "Position_play",
  ylab = "CAR",
  messages = FALSE,            
  package = "ggsci",           
  palette = "nrc_npg"
) +
  theme(
    plot.title = element_text(size = 14, face = "bold"),
    axis.title = element_text(size = 12)
  )

Nada de interés aquí. No vale la pena mostrar nada. Solo decirlo de palabra. El gráfico es para que veáis los datos reales.

 wide %>% 
  anova_test(logCAR_Ratio ~ Position_play, detailed = TRUE)
## ANOVA Table (type II tests)
## 
##          Effect   SSn   SSd DFn DFd     F     p p<.05   ges
## 1 Position_play 0.193 5.123   3  48 0.604 0.616       0.036
ggbetweenstats(
  data = wide,
  x = Position_play,
  y = logCAR_Ratio,
  type = "parametric",      
  conf.level = 0.95,           
  plot.type = "boxviolin",     
  mean.plotting = TRUE,        
  centrality.plotting = TRUE,  
  title = "Differences in log(C1/C0) by Position_play",
  xlab = "Position_play",
  ylab = "log(C1/C0)",
  messages = FALSE,            
  package = "ggsci",           
  palette = "nrc_npg"
) +
  theme(
    plot.title = element_text(size = 14, face = "bold"),
    axis.title = element_text(size = 12)
  )

Variante con un ANOVA como el del primer paso pero añadiendo como variable explicativa la posición de juego

Esto no tiene interés tras comentarlo con Manuel no es una hipótesis que busquéis, pero lo dejo por que si tenéis curiosidad de ver qué diferencias hay entre posiciones y donde podéis ir cortos o largos de muestra. Lo mismo queréis retirar a los porteros o quien sabe si tratarlos como defensas… Al menos que se vean las medidas repetidas en el mismo gráfico.

Estructura del modelo:

  • Within-subject factor: Time (4 niveles: C0, C1, C2, C3)
  • Between-subject factor: Position_play
aov_mixto_pos <- long %>%
  anova_test(
    logCORT ~ Outcome * Time + Position_play,
    wid = CODIGO,
    within = Time,
    between = c(Outcome, Position_play)
  )
get_anova_table(aov_mixto_pos) %>% knitr::kable()
Effect DFn DFd F p p<.05 ges
Outcome 1 197 10.185 2e-03 * 0.049
Time 3 197 19.247 0e+00 * 0.227
Position_play 3 197 5.116 2e-03 * 0.072
Outcome:Time 3 197 11.735 4e-07 * 0.152 
lmer_model <- lmer(
  logCORT ~  Time * Position_play+(1 | CODIGO),
  data = long
)


tabla_lmm_apa <- modelsummary(
  lmer_model,
  # Configuración del modelo:
  stars = TRUE,          
  statistic = "conf.int",


  output = "tinytable", 
  fmt = fmt_decimal(digits = 3),

  gof_map = "default",
  title = "Modelo LMM: Efectos del Resultado, Posición  y el Tiempo en log(Cortisol) con interacciones"
)
tabla_lmm_apa
Modelo LMM: Efectos del Resultado, Posición y el Tiempo en log(Cortisol) con interacciones
(1)
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
(Intercept) 1.743***
[1.462, 2.025]
TimeC1 0.405*
[0.080, 0.731]
TimeC2 0.144
[-0.182, 0.469]
TimeC3 0.574***
[0.249, 0.900]
Position_playCenter -0.228
[-0.618, 0.163]
Position_playDefense -0.263
[-0.621, 0.095]
Position_playGoalie -0.292
[-0.872, 0.288]
TimeC1 × Position_playCenter 0.146
[-0.305, 0.598]
TimeC2 × Position_playCenter -0.246
[-0.698, 0.206]
TimeC3 × Position_playCenter -0.061
[-0.512, 0.391]
TimeC1 × Position_playDefense 0.045
[-0.369, 0.459]
TimeC2 × Position_playDefense -0.209
[-0.623, 0.205]
TimeC3 × Position_playDefense -0.122
[-0.536, 0.292]
TimeC1 × Position_playGoalie -0.038
[-0.709, 0.632]
TimeC2 × Position_playGoalie 0.146
[-0.525, 0.817]
TimeC3 × Position_playGoalie 0.084
[-0.586, 0.755]
SD (Intercept CODIGO) 0.296
SD (Observations) 0.421 
ci_log_scale_pos <- long %>%
  group_by(Time, Position_play) %>%
  get_summary_stats(logCORT, type = "mean_se") %>% 
  mutate(
    t_crit = qt(0.975, df = n - 1),
    mean_log = mean,
    ci_lower_log = mean_log - (t_crit * se),
    ci_upper_log = mean_log + (t_crit * se)
  ) %>%
  mutate(
    mean_CORT_geom = exp(mean_log),
    ci_lower_CORT = exp(ci_lower_log),
    ci_upper_CORT = exp(ci_upper_log)
  )

grafico_posiciones <- ci_log_scale_pos %>%
  ggplot(aes(x = Time, y = mean_CORT_geom)) +
  coord_cartesian(ylim=c(0,20))+  
  # Puntos, líneas y barras de error
  geom_point(size = 2) +
  geom_line(linewidth = 0.8) +
  geom_errorbar(aes(ymin = ci_lower_CORT, ymax = ci_upper_CORT), width = 0.2) +
  
  # CREAR LAS FACETAS (Paneles por posición)
  facet_wrap(~ Position_play, ncol = 2) + 
  
  labs(
    title = "Evolución de Cortisol por Posición de Juego",
    caption = "Las barras de error representan el 95% CI transformado de la escala logarítmica.",
    x = "Tiempo de Medición",
    y = "Media Geométrica del Cortisol",
    color = "Resultado"
  ) +
  theme_minimal() +
  theme(
    legend.position = "bottom",
    strip.background = element_rect(fill = "gray90"), 
    strip.text = element_text(face = "bold")          
  )

print(grafico_posiciones)
## Ignoring unknown labels:
## • colour : "Resultado"
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?

Variante con un ANOVA como el del primer paso pero añadiendo como variable explicativa la posición de juego y la victoria/derrota

No tiene interés, pero lo dejo para que tengáis una idea visual de medidas repetidas de cortisol y patrones que se puedan asociar con posición y predicción de derrota.

Estructura del modelo:

  • Within-subject factor: Time (4 niveles: C0, C1, C2, C3)
  • Between-subject factor: Outcome (Win vs Loss) y Position_play
aov_mixto_pos <- long %>%
  anova_test(
    logCORT ~ Outcome * Time + Position_play,
    wid = CODIGO,
    within = Time,
    between = c(Outcome, Position_play)
  )
get_anova_table(aov_mixto_pos) %>% knitr::kable()
Effect DFn DFd F p p<.05 ges
Outcome 1 197 10.185 2e-03 * 0.049
Time 3 197 19.247 0e+00 * 0.227
Position_play 3 197 5.116 2e-03 * 0.072
Outcome:Time 3 197 11.735 4e-07 * 0.152

Sale todo interesante. Veamoslo en una tabla de tipo publicación:

lmer_model <- lmer(
  logCORT ~ Outcome * Time + Position_play+(1 | CODIGO),
  data = long
)


tabla_lmm_apa <- modelsummary(
  lmer_model,
  # Configuración del modelo:
  stars = TRUE,          
  statistic = "conf.int",


  output = "tinytable", 
  fmt = fmt_decimal(digits = 3),

  gof_map = "default",
  title = "Modelo LMM: Efectos del Resultado, Posición  y el Tiempo en log(Cortisol) con interacciones"
)
tabla_lmm_apa
Modelo LMM: Efectos del Resultado, Posición y el Tiempo en log(Cortisol) con interacciones
(1)
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
(Intercept) 1.717***
[1.482, 1.953]
OutcomeWin 0.100
[-0.156, 0.357]
TimeC1 0.591***
[0.403, 0.778]
TimeC2 0.007
[-0.180, 0.194]
TimeC3 0.947***
[0.760, 1.134]
Position_playCenter -0.244+
[-0.511, 0.023]
Position_playDefense -0.316*
[-0.560, -0.071]
Position_playGoalie -0.220
[-0.617, 0.176]
OutcomeWin × TimeC1 -0.283*
[-0.559, -0.008]
OutcomeWin × TimeC2 -0.005
[-0.281, 0.270]
OutcomeWin × TimeC3 -0.936***
[-1.211, -0.661]
SD (Intercept CODIGO) 0.303
SD (Observations) 0.355

En general los delanteros tienen valores más altos de cortisol que el resto, diferenciandose significativamente de los defensas, y casi significativamente de los defensas. Los porteros no van muy atrás, pero deben ser pocos y el intervalo de confianza es ancho.

De forma gráfica sería esto… Tal vez algo complicado:

ci_log_scale_pos <- long %>%
  group_by(Time, Outcome, Position_play) %>%
  get_summary_stats(logCORT, type = "mean_se") %>% 
  mutate(
    t_crit = qt(0.975, df = n - 1),
    mean_log = mean,
    ci_lower_log = mean_log - (t_crit * se),
    ci_upper_log = mean_log + (t_crit * se)
  ) %>%
  mutate(
    mean_CORT_geom = exp(mean_log),
    ci_lower_CORT = exp(ci_lower_log),
    ci_upper_CORT = exp(ci_upper_log)
  )

grafico_posiciones <- ci_log_scale_pos %>%
  ggplot(aes(x = Time, y = mean_CORT_geom, group = Outcome, color = Outcome)) +
  coord_cartesian(ylim=c(0,20))+  
  # Puntos, líneas y barras de error
  geom_point(size = 2) +
  geom_line(linewidth = 0.8) +
  geom_errorbar(aes(ymin = ci_lower_CORT, ymax = ci_upper_CORT), width = 0.2) +
  
  # CREAR LAS FACETAS (Paneles por posición)
  facet_wrap(~ Position_play, ncol = 2) + 
  
  labs(
    title = "Evolución de Cortisol por Posición de Juego",
    subtitle = "Comparación entre Resultado (Outcome) y Tiempo (Time)",
    caption = "Las barras de error representan el 95% CI transformado de la escala logarítmica.",
    x = "Tiempo de Medición",
    y = "Media Geométrica del Cortisol",
    color = "Resultado"
  ) +
  theme_minimal() +
  theme(
    legend.position = "bottom",
    strip.background = element_rect(fill = "gray90"), 
    strip.text = element_text(face = "bold")          
  )

print(grafico_posiciones)

Tendería a pensar que el CAR=C1-C0, es un predictor de la derrota cuando se mide en delanteros y centrales. A defensas y porteros no les afecta.

¿Es CAR un predictor de la derrota cuando nos limitamos a centrales y delanteros?

No tiene interés, pero lo miro solo por curiosear.

wideLimitado=wide %>% filter(Position_play %in% c("Forward","Center"))

coef_rename <- c(
  "(Intercept)" = "Constante (Intersección)",
  "I(log(CORT_1) - log(CORT_0))" = "Cambio Log(C1 / C0)"
)

testeandoIdea <- glm(I(Outcome=="Loss") ~ I(log(CORT_1)-log(CORT_0)), 
                 data=wideLimitado, 
                 family=binomial)


modelsummary(
  testeandoIdea,
  exponentiate = TRUE,      # ¡Importante! Muestra Odds Ratios en lugar de log-odds
  coef_map = coef_rename,   # Aplica los nombres limpios
  stars = TRUE,             # Añade los asteriscos de significancia (p < .05, etc.)
  statistic = "conf.int",   # Muestra el Intervalo de Confianza en lugar del error estándar
  title = "Odds Ratios: Predicción de Pérdida (Loss) según cambio en Cortisol",
  gof_map = c("nobs", "aic", "logLik"), # Selecciona qué estadísticos mostrar al final
  output = "gt"             # Formato de alta calidad para Word/HTML
)
Odds Ratios: Predicción de Pérdida (Loss) según cambio en Cortisol
(1)
Constante (Intersección) 0.069*
[0.004, 0.540]
Cambio Log(C1 / C0) 567.279*
[8.998, 316190.640]
Num.Obs. 27
AIC 99.1
Log.Lik. -47.565
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

Impresionante la OR.

predicciones <- predict(testeandoIdea, type = "response")

curva_roc <- roc(wideLimitado$Outcome == "Loss", predicciones)
## Setting levels: control = FALSE, case = TRUE
## Setting direction: controls < cases
# 3. Dibujar la curva
plot(curva_roc, 
     main = "Curva ROC del Modelo de Cortisol sobre Medios y Delanteros",
     col = "#2196F3", 
     lwd = 4, 
     print.auc = TRUE,           # Imprime el valor del AUC en el gráfico
     legacy.axes = TRUE,         # Eje X como 1-Especificidad (estándar de publicación)
     xlab = "1 - Especificidad (Falsos Positivos)",
     ylab = "Sensibilidad (Verdaderos Positivos)")

# 4. Ver el AUC exacto y su intervalo de confianza
print(auc(curva_roc))
## Area under the curve: 0.8722
print(ci.auc(curva_roc))
## 95% CI: 0.7285-1 (DeLong)

Bastante impresionante. Si quieres saber cómo va a terminar un partido fíjate nada más que el CAR de centrales y delanteros y sabrás si van a perder.

Se nota en C3 más la derrota en unas posiciones que otras?

Sin interés, pero lo miramos.

aov_mixto_pos <- long %>%
  anova_test(
    logCORT ~ Outcome * Time + Outcome*Position_play,
    wid = CODIGO,
    within = Time,
    between = c(Outcome, Position_play)
  )
get_anova_table(aov_mixto_pos) %>% knitr::kable()
Effect DFn DFd F p p<.05 ges
Outcome 1 194 10.445 1.0e-03 * 0.051
Time 3 194 19.739 0.0e+00 * 0.234
Position_play 3 194 5.246 2.0e-03 * 0.075
Outcome:Time 3 194 12.035 3.0e-07 * 0.157
Outcome:Position_play 3 194 2.678 4.8e-02 * 0.040

Bah, significativo pero por muy poco. Hipótesis muy complicada como para defenderla con esta muestra. Ya íbamos justillos con la anterior.

6. Correlaciones psicofisiológicas

Ansiedad competitiva – CSAI-2R (estado precompetitivo)

C.A. → Cognitive Anxiety S.A. → Somatic Anxiety S.C. → Self-Confidence 2. Ansiedad percibida / activación subjetiva (ad hoc, relacionada con CSAI-2R)

Anxiety → Ansiedad somática percibida Nervousness → Activación fisiológica subjetiva 3. Afecto básico – SAM (Self-Assessment Manikin)

Affective Valence (SAM) Arousal (SAM) Dominance (SAM) 4. Estado psicobiosocial / disposición competitiva (estado)

I am vigorous Today is a hard day 5. Motivación, compromiso y liderazgo (psicosocial)

Inspire my teammates Sacrificing myself Compete with passion My words help the team Defeat is failure

CSAI-2R: “Cognitive_Axiety” “Somatic_Axiety” “Self_confidence”

SAM: valencia, arousal, dominio (SAM)

DOMINIO (SAM): “SAM_Affective_valence” “SAM_AROUSAL” “SAM_DOMINANCE”

         "Salience_CAR"

variables motivacionales

“Hard_day” “Vigor_reach_challenge” “Anxiety_perception” “Nerviousness_perception”

“Inspires_teammates” “Sacrifice_myself”

“Compite_with_passion” “My_Words_motivate_teammates”

“Be_defeated_not_failure”

TTe lo he puesto por areas, que yo creo que se ve bien diferenciado. Adjunto una matriz de correlaciones que se hizo en sus momento para que le eches un vistazo también. Se me ocurre que podríamos hacer un PCA para todas estas variables y determinar cuales o que grupo de ellas son las que explican mayor variabilidad e la varianza de los datos. Dime que te parece y si te encaja??

variables motivacionales

“Hard_day” “Vigor_reach_challenge” “Anxiety_perception” “Nerviousness_perception” “Inspires_teammates” “Sacrifice_myself”

“Compite_with_passion” “My_Words_motivate_teammates” “Be_defeated_not_failure”

Supongo que se quiere una tabla de las variables que se porten bien del cortisol con cada una de las variables de cada categoría. En total 3 tablas. ¿Es así? ¿Solo una grande de todas contra todas?

Análisis PCA

Tenemos múltibles variables y cada una medida en una escala diferente. Lo normal es tipificar cada una de ellas y si nos interesa saber cómo explicar la variabilidad existente, hacemos un Análisis de Componentes Principales (PCA). En realidad creo que queréis hacer otra cosa, pero empezamos por aquí.

La matriz de variables de interés estará formada por las variables de la matriz de correlaciones que me has enviado y las medidas de cortisol en los 4 tiempos (C0, C1, C2, C3).

7. Transformaciones necesarias

Todas las medidas de cortisol deben analizarse sobre valores log-transformados, pero presentarse en gráficos con valores sin transformar

dfZ=wide %>%
  select(
    CORT_0, CORT_1, CORT_2, CORT_3,
    Cognitive_Axiety, Somatic_Axiety, Self_confidence,
    SAM_Affective_valence, SAM_AROUSAL, SAM_DOMINANCE,
    Hard_day, Vigor_reach_challenge,
    Anxiety_perception, Nerviousness_perception,
    Inspires_teammates, Sacrifice_myself,
    Compite_with_passion, My_Words_motivate_teammates,
    Be_defeated_not_failure,
    SAM_Affective_valence, SAM_AROUSAL,SAM_DOMINANCE  
  ) %>%
  mutate(across(everything(), scale)) %>%
  as.data.frame()

Veamos una tabla de correlaciones de Spearman entre ellas:

correlation_matrix <- cor(dfZ, method = "spearman", use = "pairwise.complete.obs")
# Mostrarla en forma de publicación redondeando a dos decimales

round(correlation_matrix,2) %>% 
  as.data.frame() %>%
  rownames_to_column(var = "Variable") %>%
  knitr::kable()
Variable CORT_0 CORT_1 CORT_2 CORT_3 Cognitive_Axiety Somatic_Axiety Self_confidence SAM_Affective_valence SAM_AROUSAL SAM_DOMINANCE Hard_day Vigor_reach_challenge Anxiety_perception Nerviousness_perception Inspires_teammates Sacrifice_myself Compite_with_passion My_Words_motivate_teammates Be_defeated_not_failure
CORT_0 1.00 0.75 0.49 0.19 -0.12 -0.08 -0.05 0.03 -0.06 -0.06 0.17 0.01 0.09 0.07 0.05 -0.16 -0.13 0.03 -0.12
CORT_1 0.75 1.00 0.45 0.40 0.21 0.12 -0.20 0.02 0.08 -0.13 0.14 0.02 0.12 0.27 -0.09 -0.27 -0.24 -0.05 0.04
CORT_2 0.49 0.45 1.00 0.25 0.08 0.09 0.05 0.19 -0.09 0.02 0.15 0.25 0.15 0.08 0.27 0.00 0.05 0.14 -0.02
CORT_3 0.19 0.40 0.25 1.00 0.27 0.34 -0.01 0.08 0.46 0.16 0.08 0.06 0.14 0.34 0.28 0.03 -0.04 0.17 0.20
Cognitive_Axiety -0.12 0.21 0.08 0.27 1.00 0.37 -0.21 -0.09 0.28 -0.15 0.20 0.12 0.24 0.31 -0.27 -0.05 0.14 0.11 0.03
Somatic_Axiety -0.08 0.12 0.09 0.34 0.37 1.00 0.14 0.32 0.57 0.34 0.24 0.16 0.20 0.73 0.15 0.26 -0.12 0.22 -0.04
Self_confidence -0.05 -0.20 0.05 -0.01 -0.21 0.14 1.00 0.47 -0.03 0.45 0.16 0.39 -0.27 0.04 0.12 0.42 0.12 0.01 -0.28
SAM_Affective_valence 0.03 0.02 0.19 0.08 -0.09 0.32 0.47 1.00 0.06 0.41 0.21 0.35 -0.26 0.32 0.05 0.27 0.04 0.08 -0.27
SAM_AROUSAL -0.06 0.08 -0.09 0.46 0.28 0.57 -0.03 0.06 1.00 0.13 0.13 0.02 0.11 0.49 0.27 0.27 -0.07 0.24 -0.07
SAM_DOMINANCE -0.06 -0.13 0.02 0.16 -0.15 0.34 0.45 0.41 0.13 1.00 0.14 0.33 -0.14 0.17 0.29 0.39 0.32 0.43 -0.19
Hard_day 0.17 0.14 0.15 0.08 0.20 0.24 0.16 0.21 0.13 0.14 1.00 0.10 0.15 0.24 -0.13 0.41 0.16 0.28 0.02
Vigor_reach_challenge 0.01 0.02 0.25 0.06 0.12 0.16 0.39 0.35 0.02 0.33 0.10 1.00 -0.11 -0.08 0.04 0.32 0.35 0.22 -0.34
Anxiety_perception 0.09 0.12 0.15 0.14 0.24 0.20 -0.27 -0.26 0.11 -0.14 0.15 -0.11 1.00 0.22 0.08 -0.07 -0.09 0.12 0.12
Nerviousness_perception 0.07 0.27 0.08 0.34 0.31 0.73 0.04 0.32 0.49 0.17 0.24 -0.08 0.22 1.00 0.07 0.20 -0.14 0.17 0.03
Inspires_teammates 0.05 -0.09 0.27 0.28 -0.27 0.15 0.12 0.05 0.27 0.29 -0.13 0.04 0.08 0.07 1.00 0.21 -0.03 0.39 0.20
Sacrifice_myself -0.16 -0.27 0.00 0.03 -0.05 0.26 0.42 0.27 0.27 0.39 0.41 0.32 -0.07 0.20 0.21 1.00 0.38 0.42 -0.10
Compite_with_passion -0.13 -0.24 0.05 -0.04 0.14 -0.12 0.12 0.04 -0.07 0.32 0.16 0.35 -0.09 -0.14 -0.03 0.38 1.00 0.20 -0.05
My_Words_motivate_teammates 0.03 -0.05 0.14 0.17 0.11 0.22 0.01 0.08 0.24 0.43 0.28 0.22 0.12 0.17 0.39 0.42 0.20 1.00 -0.02
Be_defeated_not_failure -0.12 0.04 -0.02 0.20 0.03 -0.04 -0.28 -0.27 -0.07 -0.19 0.02 -0.34 0.12 0.03 0.20 -0.10 -0.05 -0.02 1.00

He comprobado unos cuantos números al azar y parece coincidir exactamente con la vuestra. No veo que haya que cambiar nada. Por si queréis tener una versión equivalente más gráfica podríamos hacer algo de este estilo:

corrplot::corrplot(correlation_matrix, method = "color", type = "upper",
         tl.col = "black", tl.srt = 45, 
         addCoef.col = "black", number.cex = 0.6,
         diag = FALSE)

Os pondría las etiquetas bonitas vuestras si os interesa esta otra versión.

Reorganizamos el orden de las variablas basadas en cluster jerarquico de acuerdo al valor de las correlaciones:

corrplot::corrplot(correlation_matrix, method = "color", type = "upper",
         tl.col = "black", tl.srt = 45,
         addCoef.col = "black", number.cex = 0.6,
         diag = FALSE,
         order = "hclust",
         hclust.method = "complete")

Lo que estáis viendo consiste en poner cerca las variables que presentan más correlaciones entre sí, de forma que podáis ver si las agrupaciones de variables tienen más sentido. Veréis que las correlaciones ahora van más en forma de bloques. Es una ayuda visual solo.

Análisis PCA

Empezamos viendo un scree plot

d_pca_resultado <- prcomp(dfZ,center=TRUE,scale=FALSE)

  ncp=ncol(dfZ)
  factoextra::fviz_screeplot(d_pca_resultado, ncp=ncp,addlabels = TRUE,choice="variance") 
## Warning in geom_bar(stat = "identity", fill = barfill, color = barcolor, :
## Ignoring empty aesthetic: `width`.

Prece que podríais elegir unas 3 dimensiones que acumularían un 19.4%+16.7%+11.5% = 47.6% más o menos de la variabilidad, y luego aparece un codo que te invita a cortar. Es arbitrario y visual, tampoco hay que tomárselo literalmente.

Análisis Factorial exploratorio

Dado que el PCA sugiere que hay 3 factores que explican casi la mitad de la variabilidad, nos preguntamos quienes podrían ser. Vamos a intentar una rotación varimax.

graficoCargas <-function(fit1,Variables){
  
  
  cargas_tidy = fit1 %>% tidy() %>%
    select(-uniqueness) %>% 
    rename(Variable=variable) %>% 
    mutate(Variable=factor(Variable,levels=rev(Variables))) %>% 
    pivot_longer(cols = -Variable, names_to = "Factor", values_to = "carga") %>%
    arrange(Variable,Factor)
  
  ggplot(cargas_tidy, aes(x = Variable, y = carga, fill = Factor)) +
    geom_bar(stat = "identity") +
    facet_wrap(~ Factor,ncol = 3)+
    coord_flip()+
    theme_minimal()+
    theme(axis.text.x = element_text(size = 6),   # Tamaño del texto del eje X
          axis.text.y = element_text(size = 6))+
    theme(legend.position = "none")
  
}

d_fit0 <- factanal(dfZ, factors = 3, scores = "Bartlett",rotation="varimax")

graficoCargas(d_fit0,dfZ %>% names())

Os pongo las cargas en forma de gráficos en lugar de en forma de números por que yo creo que se entiende mejor. El tercer factor parece relacionado con Cortisol (azul). El primero (rojo) es un potaje de muchas cosasy el segundo (verde) es un potaje de otras.

Lo mismo vosotros que conocéis las variables os sugiere algo.Tal vez queráis limitar el número de variables de la matriz de correlaciones y reintentarlo con un subconjunto de variables más pequeño. Normalmente queréis que cada item cargue (tenga una barra grande) en solo uno de cada factor. No hacen gracia los items que reflejan valores altos en dos factores.