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:

Within-subject factor: Time (4 niveles: C0, C1, C2, C3)
Between-subject factor: Outcome (Win vs Loss)

Requisitos:

Verificar normalidad (Kolmogorov–Smirnov)
Verificar homogeneidad (Levene)
Verificar esfericidad (Mauchly)
Aplicar corrección Greenhouse–Geisser si procede
Reportar:
F, p
η²p -95% CI de los efectos
Efecto principal de Time, de Outcome y la interacción Time × Outcome

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 = 28¹	Win N = 24¹	p-valor²	Size effect ( $r$ )
CAR=C1-C0	3.78 (2.50, 4.89)	1.84 (1.04, 2.66)	<0.001	0.58
¹ Median (Q1, Q3)
² 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 = 28¹	Win N = 24¹	p-valor²	Size effect ( $r$ )
Δ% cortisol	128.10 (95.43, 244.69)	-12.60 (-32.90, 56.33)	<0.001	0.74
¹ Median (Q1, Q3)
² 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 (NO TENGO CLARA LA VARIABLE DEPENDIENTE!!)

Prueba adecuada:

One-way ANOVA con Position (4 niveles) como factor Post-hoc Bonferroni 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

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

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

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?

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-Cortisol (T1 - T0)"
)

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-Cortisol (T1 - T0)	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?

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

Dado que la mayoría de variables no son normales:

Spearman’s rho entre: cortisol variables ansiedad (CSAI-2R) valencia, arousal, dominio (SAM) variables motivacionales Reportar:

rho p 95% CI si es posible

Las variables de que dispongo son estas:

CORTISOL:

“CORT_0” “CORT_1” “CORT_2” “CORT_3”

“CAR” “CAR_Ratio”Delta_CORT”

Me falta saber donde encajan las que restan en esas categorías. Hago un intento pero me tenéis que decir si está bien y corregir lo que esté mal:

ansiedad (CSAI-2R)

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

valencia, arousal, dominio (SAM)

“Inspires_teammates” “Sacrifice_myself”
“Cognitive_Axiety” “Somatic_Axiety” “Self_confidence” “Salience_CAR”
“SAM_Affective_valence” “SAM_AROUSAL” “SAM_DOMINANCE”

variables motivacionales

“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?

7. Transformaciones necesarias

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

hockey

2025-12-10

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

2. Fiabilidad de medidas repetidas

3. Comparaciones puntuales del CAR (C1 – C0)

Empecemos por la normalidad de CAR

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

Otra posibilidad

4. Porcentaje de cambio (pre → post)

5. Comparaciones por posición de juego (NO TENGO CLARA LA VARIABLE DEPENDIENTE!!)

Para C_0

Para C_1

Para C_2

Para C_3

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

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

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

6. Correlaciones psicofisiológicas

CORTISOL:

ansiedad (CSAI-2R)

valencia, arousal, dominio (SAM)

variables motivacionales

7. Transformaciones necesarias