4. Ancova

ANCOVA

rm(list = ls())

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(ggpubr)
library(rstatix)


Adjuntando el paquete: 'rstatix'

The following object is masked from 'package:stats':

    filter

library(broom)
library(tuneR)
library(readxl)
library(tidyverse)
library(egg)

Cargando paquete requerido: gridExtra

Adjuntando el paquete: 'gridExtra'

The following object is masked from 'package:dplyr':

    combine


Adjuntando el paquete: 'egg'

The following object is masked from 'package:ggpubr':

    ggarrange

library(tidyverse)
library(ggplot2)
library(ggpmisc)

Cargando paquete requerido: ggpp
Registered S3 methods overwritten by 'ggpp':
  method                  from   
  heightDetails.titleGrob ggplot2
  widthDetails.titleGrob  ggplot2

Adjuntando el paquete: 'ggpp'

The following objects are masked from 'package:ggpubr':

    as_npc, as_npcx, as_npcy

The following object is masked from 'package:ggplot2':

    annotate

library(broom)
library(ggplot2)
library(patchwork)
library(egg)
library(ggpubr)
library(readxl)
library(tidyverse)
library(egg)
library(tidyverse)


dat_clean=read.csv("dat_clean_modified_zscore_anchoveta.csv")


dat_clean$group <- factor(dat_clean$Subclass_n,      # Reordering group factor levels
                         levels = c("Moda 3.5 cm", "Moda 4 cm", "Moda 5 cm", "Moda 7.5 cm","Moda 10.5 cm", "Moda 11 cm", "Moda 12 cm","Moda 12.5 cm","Moda 13 cm","Moda 13.5 cm"),labels = c("3.5", "4", "5", "7.5","10.5","11","12","12.5","13","13.5"))

dat_clean$Banda <- factor(dat_clean$Banda,
  levels = c("35-45","45-90","90-170","170-260"),labels = c("35-45","45-90","90-170","170-260"))

dat=dat_clean


#dat=dat%>%
 # select(-Value)%>%
  #rename(Value="Value_linear")

anxiety <- dat %>%
  select(group,Frequency,Value)%>%
  mutate(Frecuencia=round(as.numeric(Frequency,0)),Sv=round(Value,2),Class=group)%>%
  #rename(pretest = Frecuencia, posttest = Valor, group=Class) %>%
  drop_na()%>%
  mutate(Frecuencia=as.numeric(Frecuencia))

library(dplyr)

# Especifica el tamaño del subconjunto balanceado
tamaño_subconjunto <- 1000

# Realiza el muestreo aleatorio y balanceado
subset_balanceado <- anxiety %>%
  group_by(Class) %>%
  sample_n(tamaño_subconjunto, replace = T) %>%
  ungroup()

# Verifica el tamaño del subconjunto balanceado
nrow(subset_balanceado)

[1] 10000

anxiety=subset_balanceado

# Ahora `subset_balanceado` contiene un subconjunto aleatorio y balanceado de tus datos.

Comprobar supuestos Supuesto de linealidad

Comprobar supuestos

Supuesto de linealidad

Cree un diagrama de dispersión entre la covariable (es decir, ) y la variable de resultado (es decir, frecuencia). Agregue líneas de regresión, muestre las ecuaciones correspondientes y el R2 por grupos.

ggscatter(
  anxiety, x = "Frecuencia", y = "Sv",
  color = "Class", add = "reg.line",size = 0.5, alpha=0.5
  )+
  stat_regline_equation(
    aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"), color = Class)
    )+
  theme_bw()

Warning: The dot-dot notation (`..eq.label..`) was deprecated in ggplot2 3.4.0.
ℹ Please use `after_stat(eq.label)` instead.
ℹ The deprecated feature was likely used in the ggpubr package.
  Please report the issue at <https://github.com/kassambara/ggpubr/issues>.

Hubo una relación lineal entre la puntuación de ansiedad antes y después de la prueba para cada grupo de entrenamiento, según se evaluó mediante la inspección visual de un diagrama de dispersión.

Homogeneidad de las pendientes de regresión

Esta suposición verifica que no haya interacción significativa entre la covariable y la variable de agrupación. Esto se puede evaluar de la siguiente manera:

anxiety %>% anova_test(Sv ~ Class*Frecuencia)

ANOVA Table (type II tests)

            Effect DFn  DFd       F        p p<.05   ges
1            Class   9 9980 330.260 0.00e+00     * 0.229
2       Frecuencia   1 9980 415.103 1.93e-90     * 0.040
3 Class:Frecuencia   9 9980  43.529 2.90e-77     * 0.038

Hubo homogeneidad de pendientes de regresión ya que el término de interacción no fue estadísticamente significativo, F(4, 4990) = 76.544, p = *…..en realidad fue significativa

Normalidad de los residuos

En primer lugar, debe calcular el modelo mediante . En R, puede aumentar fácilmente sus datos para agregar valores ajustados y residuales mediante la función [paquete de escobas]. Llamemos a la salida porque contiene varias métricas útiles para el diagnóstico de regresión.lm()augment(model)model.metrics

Normality of residuals You first need to compute the model using lm(). In R, you can easily augment your data to add fitted values and residuals by using the function augment(model) [broom package]. Let’s call the output model.metrics because it contains several metrics useful for regression diagnostics.

# Fit the model, the covariate goes first
model <- lm(Sv ~ Frecuencia + Class, data = anxiety)
model


Call:
lm(formula = Sv ~ Frecuencia + Class, data = anxiety)

Coefficients:
(Intercept)   Frecuencia       Class4       Class5     Class7.5    Class10.5  
  -52.31979     -0.02242     -2.24026     -3.47311     -0.41229     -1.98671  
    Class11      Class12    Class12.5      Class13    Class13.5  
   -1.03635      1.29878      4.08846      5.94212      9.23484

# Inspect the model diagnostic metrics
model.metrics <- augment(model) %>%
  select(-.hat, -.sigma, -.fitted, -.fitted) # Remove details

model.metrics

# A tibble: 10,000 × 6
      Sv Frecuencia Class  .resid    .cooksd .std.resid
   <dbl>      <dbl> <fct>   <dbl>      <dbl>      <dbl>
 1 -65.3         42 3.5   -12.1   0.000337       -1.69 
 2 -54.7        223 3.5     2.65  0.0000139       0.370
 3 -65.5        103 3.5   -10.9   0.000222       -1.51 
 4 -59.7        112 3.5    -4.91  0.0000445      -0.685
 5 -50.1        203 3.5     6.76  0.0000857       0.943
 6 -61.4        250 3.5    -3.49  0.0000264      -0.486
 7 -48.8        188 3.5     7.72  0.000109        1.08 
 8 -55.7        110 3.5    -0.884 0.00000145     -0.123
 9 -47.4        185 3.5     9.09  0.000150        1.27 
10 -45.9         69 3.5     7.97  0.000132        1.11 
# ℹ 9,990 more rows

#Assess normality of residuals using shapiro wilk test
#shapiro_test(model.metrics$.resid)

La prueba de Shapiro Wilk no fue significativa (p > 0,05), por lo que podemos suponer normalidad de los residuos

Homogeneidad de las varianzas

Homogeneidad de las varianzas ANCOVA asume que la varianza de los residuos es igual para todos los grupos. Esto se puede comprobar mediante el test de Levene:

Homogeneity of variances ANCOVA assumes that the variance of the residuals is equal for all groups. This can be checked using the Levene’s test:

model.metrics %>% levene_test(.resid ~ Class)

# A tibble: 1 × 4
    df1   df2 statistic         p
  <int> <int>     <dbl>     <dbl>
1     9  9990      83.1 1.64e-149

La prueba de Levene no fue significativa (p > 0,05), por lo que podemos suponer homogeneidad de las varianzas residuales para todos los grupos.

The Levene’s test was not significant (p > 0.05), so we can assume homogeneity of the residual variances for all groups.

Outliers

An outlier is a point that has an extreme outcome variable value. The presence of outliers may affect the interpretation of the model.

Outliers can be identified by examining the standardized residual (or studentized residual), which is the residual divided by its estimated standard error. Standardized residuals can be interpreted as the number of standard errors away from the regression line.

Las observaciones cuyos residuos estandarizados son mayores que 3 en valor absoluto son posibles valores atípicos

model.metrics %>% 
  filter(abs(.std.resid) > 3) %>%
  as.data.frame()

       Sv Frecuencia Class    .resid      .cooksd .std.resid
1  -76.86         59   3.5 -23.21750 0.0011646303  -3.238789
2  -75.90         53   3.5 -22.39201 0.0011088290  -3.123680
3  -30.85        181   7.5  25.93990 0.0012122347   3.618192
4  -30.48         46   7.5  23.28334 0.0012350603   3.248082
5  -35.72        204   7.5  21.58554 0.0008749987   3.010894
6  -30.82         35   7.5  22.69674 0.0012295393   3.166346
7  -33.01        193   7.5  24.04893 0.0010616724   3.354464
8  -35.21        204   7.5  22.09554 0.0009168342   3.082032
9  -28.99        160   7.5  27.32910 0.0013233925   3.811931
10 -31.66         46   7.5  22.10334 0.0011130468   3.083469
11 -36.01        252   7.5  22.37164 0.0010936781   3.120816
12 -30.85        181   7.5  25.93990 0.0012122347   3.618192
13 -34.20        212   7.5  23.28489 0.0010384280   3.247965
14 -80.33         67  10.5 -24.52143 0.0012676396  -3.420635
15 -77.25         54  10.5 -21.73288 0.0010462048  -3.031735
16 -78.60         72  10.5 -22.67934 0.0010652061  -3.163637
17 -78.13         76  10.5 -22.11966 0.0009994862  -3.085541
18 -81.40         51  10.5 -25.95014 0.0015096457  -3.620069
19 -79.75         40  10.5 -24.54674 0.0014139145  -3.424395
20 -77.25         54  10.5 -21.73288 0.0010462048  -3.031735
21 -82.45         55  10.5 -26.91046 0.0015977425  -3.753997
22 -79.66         55    12 -27.40595 0.0016743448  -3.823143
23 -77.07         81    12 -24.23306 0.0011900653  -3.380330
24 -75.79         81    12 -22.95306 0.0010676661  -3.201780
25 -31.31         92    12  21.77355 0.0009284796   3.037188
26 -77.54         68    12 -24.99451 0.0013247553  -3.486639
27 -74.57         36    12 -22.74191 0.0012487804  -3.172673
28 -24.69        169  13.5  22.18375 0.0008993250   3.094292
29 -23.10         92  13.5  22.04749 0.0008955581   3.075300
30 -23.64        170  13.5  23.25617 0.0009901023   3.243881

COMPUTATION

El orden de las variables es importante a la hora de calcular ANCOVA. Primero desea eliminar el efecto de la covariable, es decir, desea controlarla, antes de ingresar su variable o interés principal.

¡La covariable va primero (y no hay interacción)! Si no lo haces en orden, obtendrás resultados diferentes.

#COMPUTATION
res.aov <- anxiety %>%
anova_test(Sv ~ Frecuencia + Class)
get_anova_table(res.aov)

ANOVA Table (type II tests)

      Effect DFn  DFd       F        p p<.05   ges
1 Frecuencia   1 9989 399.784 3.08e-87     * 0.038
2      Class   9 9989 318.072 0.00e+00     * 0.223

Después del ajuste por la puntuación de ansiedad previa a la prueba, hubo una diferencia estadísticamente significativa en la puntuación de ansiedad posterior a la prueba entre los grupos, F(4, 4994) = 580.579, p < 0,0001 *.

POST HOC TEST

Se pueden realizar comparaciones por pares para identificar qué grupos son diferentes. Se aplica la corrección de pruebas múltiples de Bonferroni. Esto se puede hacer fácilmente usando la función [paquete rstatix], un envoltorio alrededor del paquete, que debe instalarse. Emmeans significa medias marginales estimadas (también conocidas como medias mínimas cuadráticas o medias ajustadas).emmeans_test()emmeans

#POST HOC TEST
# Pairwise comparisons
# Pairwise comparisons
library(emmeans)

Welcome to emmeans.
Caution: You lose important information if you filter this package's results.
See '? untidy'

pwc <- anxiety %>% 
  emmeans_test(
    Sv ~ Class, covariate = Frecuencia,
    p.adjust.method = "bonferroni"
    )
pwc

# A tibble: 45 × 9
   term     .y.   group1 group2    df statistic         p     p.adj p.adj.signif
 * <chr>    <chr> <chr>  <chr>  <dbl>     <dbl>     <dbl>     <dbl> <chr>       
 1 Frecuen… Sv    3.5    4       9989      6.98 3.08e- 12 1.39e- 10 ****        
 2 Frecuen… Sv    3.5    5       9989     10.8  3.63e- 27 1.64e- 25 ****        
 3 Frecuen… Sv    3.5    7.5     9989      1.29 1.99e-  1 1   e+  0 ns          
 4 Frecuen… Sv    3.5    10.5    9989      6.19 6.13e- 10 2.76e-  8 ****        
 5 Frecuen… Sv    3.5    11      9989      3.23 1.24e-  3 5.58e-  2 ns          
 6 Frecuen… Sv    3.5    12      9989     -4.05 5.20e-  5 2.34e-  3 **          
 7 Frecuen… Sv    3.5    12.5    9989    -12.7  6.46e- 37 2.91e- 35 ****        
 8 Frecuen… Sv    3.5    13      9989    -18.5  2.50e- 75 1.12e- 73 ****        
 9 Frecuen… Sv    3.5    13.5    9989    -28.7  2.97e-174 1.34e-172 ****        
10 Frecuen… Sv    4      5       9989      3.84 1.23e-  4 5.52e-  3 **          
# ℹ 35 more rows

# Display the adjusted means of each group
# Also called as the estimated marginal means (emmeans)
emm=get_emmeans(pwc)
emm

# A tibble: 10 × 8
   Frecuencia Class emmean    se    df conf.low conf.high method      
        <dbl> <fct>  <dbl> <dbl> <dbl>    <dbl>     <dbl> <chr>       
 1       152. 3.5    -55.7 0.227  9989    -56.2     -55.3 Emmeans test
 2       152. 4      -58.0 0.227  9989    -58.4     -57.5 Emmeans test
 3       152. 5      -59.2 0.227  9989    -59.7     -58.8 Emmeans test
 4       152. 7.5    -56.2 0.227  9989    -56.6     -55.7 Emmeans test
 5       152. 10.5   -57.7 0.227  9989    -58.2     -57.3 Emmeans test
 6       152. 11     -56.8 0.227  9989    -57.2     -56.3 Emmeans test
 7       152. 12     -54.4 0.227  9989    -54.9     -54.0 Emmeans test
 8       152. 12.5   -51.6 0.227  9989    -52.1     -51.2 Emmeans test
 9       152. 13     -49.8 0.227  9989    -50.2     -49.4 Emmeans test
10       152. 13.5   -46.5 0.228  9989    -47.0     -46.1 Emmeans test

emm2=emm %>%
  mutate(emmean_abs=abs(emmean), sv=10*log10(emmean_abs))
  

emm2

# A tibble: 10 × 10
   Frecuencia Class emmean    se    df conf.low conf.high method      emmean_abs
        <dbl> <fct>  <dbl> <dbl> <dbl>    <dbl>     <dbl> <chr>            <dbl>
 1       152. 3.5    -55.7 0.227  9989    -56.2     -55.3 Emmeans te…       55.7
 2       152. 4      -58.0 0.227  9989    -58.4     -57.5 Emmeans te…       58.0
 3       152. 5      -59.2 0.227  9989    -59.7     -58.8 Emmeans te…       59.2
 4       152. 7.5    -56.2 0.227  9989    -56.6     -55.7 Emmeans te…       56.2
 5       152. 10.5   -57.7 0.227  9989    -58.2     -57.3 Emmeans te…       57.7
 6       152. 11     -56.8 0.227  9989    -57.2     -56.3 Emmeans te…       56.8
 7       152. 12     -54.4 0.227  9989    -54.9     -54.0 Emmeans te…       54.4
 8       152. 12.5   -51.6 0.227  9989    -52.1     -51.2 Emmeans te…       51.6
 9       152. 13     -49.8 0.227  9989    -50.2     -49.4 Emmeans te…       49.8
10       152. 13.5   -46.5 0.228  9989    -47.0     -46.1 Emmeans te…       46.5
# ℹ 1 more variable: sv <dbl>

Los datos se ajustan a la media +/- error estándar. La puntuación media de ansiedad fue estadísticamente significativa mayor en el grp1 (16,4 +/- 0,15) en comparación con el grp2 (15,8 +/- 0,12) y el grp3 (13,5 +/_ 0,11), p < 0,001.

Informe

Se realizó un ANCOVA para determinar el efecto de los ejercicios en la puntuación de ansiedad después de controlar la puntuación de ansiedad basal de los participantes.

El análisis post hoc se realizó con un ajuste de Bonferroni. La puntuación media de ansiedad fue estadísticamente significativa mayor en el grp1 (16,4 +/- 0,15) en comparación con el grp2 (15,8 +/- 0,12) y el grp3 (13,5 +/_ 0,11), p < 0,001.

# Visualization: line plots with p-values
pwc <- pwc %>% add_xy_position(x = "Class", fun = "mean_se")

ggline(get_emmeans(pwc), x = "Class", y = "emmean") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, color=Class), width = 1) + 
  stat_pvalue_manual(pwc, hide.ns = TRUE, tip.length = FALSE) +
  labs(
    subtitle = get_test_label(res.aov, detailed = TRUE),
    caption = get_pwc_label(pwc)
  )+
  theme_presentation(base_size = 11)+
  
  
    scale_color_manual(name="Longitud (cm)",values =c("#5f5f5f","#0000ff","#000080","#00bf00","#008000","#ffff00","#ff8000","#ff00bf","#ff0000","#a6533c"))+
  
    #scale_color_viridis_d(option = "C")+
  # Modificación para rotar las etiquetas del eje x
theme(axis.text.x = element_text(angle = 0, vjust = 0.5, hjust=1),legend.position = "none")+
  labs(x="Longitud (cm)",y="Emmean (Sv,dB)")

# Visualization: line plots with p-values
pwc <- pwc %>% add_xy_position(x = "Class", fun = "mean_se")

ggline(get_emmeans(pwc), x = "Class", y = "emmean") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, color=Class), width = 1) + 
  #stat_pvalue_manual(pwc, hide.ns = T, tip.length = FALSE) +
  labs(
    subtitle = get_test_label(res.aov, detailed = TRUE),
    caption = get_pwc_label(pwc)
  )+
  theme_presentation(base_size = 11)+
  
  
    scale_color_manual(name="Longitud (cm)",values =c("#5f5f5f","#0000ff","#000080","#00bf00","#008000","#ffff00","#ff8000","#ff00bf","#ff0000","#a6533c"))+
  
    #scale_color_viridis_d(option = "C")+
  # Modificación para rotar las etiquetas del eje x
theme(axis.text.x = element_text(angle = 0, vjust = 0.5, hjust=1),legend.position = "none")+
  labs(x="Longitud (cm)",y="Emmean (Sv,dB)")

Referencias para buscar en la web: ANCOVA

1-s2.0-S0022191021001335-mmc1.docx (live.com)

Acclimation, duration and intensity of cold exposure determine the rate of cold stress accumulation and mortality in Drosophila suzukii - ScienceDirect