1 Introducción

Se presenta el análisis estadístico del presupuesto de actividad para un ejemplo de Potos flavus.

Se utiliza un diseño N-of-1 de tipo AB:

Fase A: observación.
Fase B: enriquecimiento ambiental.

Acá presento muchos análisis y gráficos, no es necesario usarlos todos. Pueden elegir los que se ajusten mejor a lo que quieran responder.

Elegí un color para cada unidad comportamental y ese mismo color se repite en todos los análisis y gráficos. En caso de agregar una unidad comportamental adicional a las que tengo, deben agregar un nuevo color también.

En la esquina inferior derecha de cada bloque encuentran el botón show para ver cada uno de los scripts.

En la esquina superior izquierda de cada bloque de código encuentran el botón para copiar y así pegar en sus scripts.

2 Instalación y carga de paquetes

paquetes <- c(
  "readxl",
  "dplyr",
  "tidyr",
  "rphylopic",
  "ggplot2",
  "purrr",
  "stringr",
  "janitor",
  "rstatix",
  "vegan",
  "writexl",
  "car",
  "qcc"
)

instalar <- paquetes[
  !paquetes %in% rownames(installed.packages())
]

if (length(instalar) > 0) {
  install.packages(
    instalar,
    dependencies = TRUE
  )
}

invisible(
  lapply(
    paquetes,
    library,
    character.only = TRUE
  )
)

2.1 Silueta de Potos flavus

Esto es solamente para agregar la figura de la especie que están trabajando en los gráficos y figuras.

Si quieren hacerlo, deben cambiar el nombre de la especie entre comillas por la de ustedes; si no, pueden saltarse este paso y eliminar todo lo que incluya add_phylopic de los gráficos.

potosimg <- get_uuid(
  name = "Potos flavus"
)

potos <- get_phylopic(potosimg)

3 Importación de los datos

Mi archivo de excel se llama Potos.xlsx. Deben reemplazar esto por el nombre del archivo de ustedes tal como aparece en la carpeta. El archivo Potos.xlsx debe encontrarse en la misma carpeta que este script.

datos <- read_excel(
  path = "Potos.xlsx",
  sheet = "Estados"
)

4 Limpieza y organización

Esto lo hago con el fin de verificar que todos los datos del excel esten correctos y no haya ningún error.

datos <- datos %>%
  clean_names() %>%
  mutate(
    dia = as.numeric(dia),

    fase = factor(
      fase,
      levels = c(
        "Observacion",
        "Enriquecimiento"
      )
    ),

    ea = as.factor(ea)
  ) %>%
  arrange(dia)

4.1 Revisión de la base de datos

4.1.1 Nombres de las variables

Son los encabezados de cada columna del excel que debe incluir: día, ea y las unidadas comportamentales que usen.

names(datos)

## [1] "dia"                       "fase"                     
## [3] "ea"                        "alimentacion"             
## [5] "territorialidad"           "locomocion"               
## [7] "descanso"                  "mantenimiento_autocuidado"
## [9] "estereotipias"

4.1.2 Primeros registros

Acá les debe mostrar las primeras filas de sus datos

head(datos)

4.1.3 Estructura de la base

Les muestra como esta constituido cada objeto de la tabla.

Verifiquen que:

Día sea un num y que incluya la totalidad de su muestreo.
Fase sea un factor de 2 niveles (Observacion y enriquecimiento).
Todas las unidades comportamentales sean num y que incluya la totalidad de su muestreo.

str(datos)

## tibble [20 × 9] (S3: tbl_df/tbl/data.frame)
##  $ dia                      : num [1:20] 1 2 3 4 5 6 7 8 9 10 ...
##  $ fase                     : Factor w/ 2 levels "Observacion",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ ea                       : Factor w/ 10 levels "Aguila","Arena",..: 7 7 7 7 7 7 7 7 7 7 ...
##  $ alimentacion             : num [1:20] 53.8 23 17.5 40 15.8 ...
##  $ territorialidad          : num [1:20] 19.5 22 18 13.3 13.6 ...
##  $ locomocion               : num [1:20] 19.3 42.8 44 39.1 53.3 ...
##  $ descanso                 : num [1:20] 0 0 0 0 0 0 0 0 0 0 ...
##  $ mantenimiento_autocuidado: num [1:20] 4.095 4.443 1.291 2.381 0.386 ...
##  $ estereotipias            : num [1:20] 3.34 7.69 19.17 5.12 16.93 ...

4.2 Resumen general

summary(datos)

##       dia                     fase           ea      alimentacion   
##  Min.   : 1.00   Observacion    :10   Ninguno :10   Min.   : 8.231  
##  1st Qu.: 5.75   Enriquecimiento:10   Recinto : 2   1st Qu.:17.103  
##  Median :10.50                        Aguila  : 1   Median :21.793  
##  Mean   :10.50                        Arena   : 1   Mean   :28.844  
##  3rd Qu.:15.25                        Ficus   : 1   3rd Qu.:47.629  
##  Max.   :20.00                        Forrajeo: 1   Max.   :56.947  
##                                       (Other) : 4                   
##  territorialidad    locomocion       descanso mantenimiento_autocuidado
##  Min.   : 3.941   Min.   :17.95   Min.   :0   Min.   :0.3861           
##  1st Qu.: 7.931   1st Qu.:34.46   1st Qu.:0   1st Qu.:0.9716           
##  Median :12.946   Median :43.89   Median :0   Median :1.1469           
##  Mean   :14.221   Mean   :46.05   Mean   :0   Mean   :1.8421           
##  3rd Qu.:18.390   3rd Qu.:53.37   3rd Qu.:0   3rd Qu.:2.1320           
##  Max.   :37.974   Max.   :83.01   Max.   :0   Max.   :6.0654           
##                                                                        
##  estereotipias   
##  Min.   : 0.000  
##  1st Qu.: 0.000  
##  Median : 2.253  
##  Mean   : 9.042  
##  3rd Qu.:12.788  
##  Max.   :49.811  
##

4.3 Definición de las conductas

Agreguen las unidades comportamentales que encontraron a un vector llamado conductas.

conductas <- c(
  "alimentacion",
  "territorialidad",
  "locomocion",
  "descanso",
  "mantenimiento_autocuidado",
  "estereotipias"
)

4.4 Verificación de columnas

Verifica que las columnas existan. No debe arrojar ningún aviso

columnas_faltantes <- setdiff(
  conductas,
  names(datos)
)

if (length(columnas_faltantes) > 0) {
  stop(
    paste(
      "Faltan estas columnas en la base:",
      paste(
        columnas_faltantes,
        collapse = ", "
      )
    )
  )
}

4.5 Verificación de datos faltantes

Todo debe ser igual a 0

faltantes <- datos %>%
  summarise(
    across(
      all_of(conductas),
      ~ sum(is.na(.x))
    )
  )

print(faltantes)

## # A tibble: 1 × 6
##   alimentacion territorialidad locomocion descanso mantenimiento_autocuidado
##          <int>           <int>      <int>    <int>                     <int>
## 1            0               0          0        0                         0
## # ℹ 1 more variable: estereotipias <int>

4.6 Verificación de los porcentajes

Los porcentajes de las conductas deben sumar 100 % para cada día de observación.

datos <- datos %>%
  mutate(
    suma_porcentajes = rowSums(
      across(all_of(conductas)),
      na.rm = TRUE
    )
  )

datos %>%
  dplyr::select(
    dia,
    fase,
    suma_porcentajes
  )

5 Transformación a formato largo

Organiza los datos en filas por unicadaes comportamentales para cada día.

datos_largos <- datos %>%
  pivot_longer(
    cols = all_of(conductas),
    names_to = "conducta",
    values_to = "porcentaje"
  ) %>%
  mutate(
    conducta = factor(
      conducta,
      levels = conductas,
      labels = c(
        "Alimentación",
        "Territorialidad",
        "Locomoción",
        "Descanso",
        "Mantenimiento",
        "Estereotipias"
      )
    )
  )

head(datos_largos)

6 Estadística descriptiva

Se calculan el tamaño de muestra, la media, la mediana, la desviación estándar, el error estándar, el mínimo, el máximo y el intervalo de confianza del 95 %.

resumen_descriptivo <- datos_largos %>%
  group_by(
    fase,
    conducta
  ) %>%
  summarise(
    n = sum(!is.na(porcentaje)),

    media = mean(
      porcentaje,
      na.rm = TRUE
    ),

    mediana = median(
      porcentaje,
      na.rm = TRUE
    ),

    desviacion_estandar = sd(
      porcentaje,
      na.rm = TRUE
    ),

    error_estandar =
      desviacion_estandar / sqrt(n),

    minimo = min(
      porcentaje,
      na.rm = TRUE
    ),

    maximo = max(
      porcentaje,
      na.rm = TRUE
    ),

    ic_95_inferior =
      media -
      qt(
        0.975,
        df = pmax(n - 1, 1)
      ) *
      error_estandar,

    ic_95_superior =
      media +
      qt(
        0.975,
        df = pmax(n - 1, 1)
      ) *
      error_estandar,

    .groups = "drop"
  )

print(resumen_descriptivo)

## # A tibble: 12 × 11
##    fase  conducta     n  media mediana desviacion_estandar error_estandar minimo
##    <fct> <fct>    <int>  <dbl>   <dbl>               <dbl>          <dbl>  <dbl>
##  1 Obse… Aliment…    10 23.7     20.6               13.9            4.41   8.40 
##  2 Obse… Territo…    10 13.4     13.4                5.81           1.84   4.76 
##  3 Obse… Locomoc…    10 43.2     43.9               10.5            3.31  19.3  
##  4 Obse… Descanso    10  0        0                  0              0      0    
##  5 Obse… Manteni…    10  2.38     1.67               1.87           0.592  0.386
##  6 Obse… Estereo…    10 17.3     14.2               16.1            5.09   0    
##  7 Enri… Aliment…    10 34.0     34.1               19.1            6.04   8.23 
##  8 Enri… Territo…    10 15.0     11.5               10.8            3.43   3.94 
##  9 Enri… Locomoc…    10 48.9     44.2               21.8            6.89  17.9  
## 10 Enri… Descanso    10  0        0                  0              0      0    
## 11 Enri… Manteni…    10  1.30     1.06               0.826          0.261  0.788
## 12 Enri… Estereo…    10  0.759    0                  1.03           0.326  0    
## # ℹ 3 more variables: maximo <dbl>, ic_95_inferior <dbl>, ic_95_superior <dbl>

7 Cambios entre las fases

Obtiene el cambio relativo y absoluto entre las fases de observación y enriquecimiento.

cambios <- resumen_descriptivo %>%
  dplyr::select(
    fase,
    conducta,
    media
  ) %>%
  pivot_wider(
    names_from = fase,
    values_from = media
  ) %>%
  mutate(
    cambio_absoluto =
      Enriquecimiento - Observacion,

    cambio_relativo_porcentaje = ifelse(
      Observacion == 0,
      NA_real_,
      100 * cambio_absoluto / Observacion
    )
  )

print(cambios)

## # A tibble: 6 × 5
##   conducta    Observacion Enriquecimiento cambio_absoluto cambio_relativo_porc…¹
##   <fct>             <dbl>           <dbl>           <dbl>                  <dbl>
## 1 Alimentaci…       23.7           34.0             10.4                    43.9
## 2 Territoria…       13.4           15.0              1.60                   11.9
## 3 Locomoción        43.2           48.9              5.67                   13.1
## 4 Descanso           0              0                0                      NA  
## 5 Mantenimie…        2.38           1.30            -1.07                  -45.2
## 6 Estereotip…       17.3            0.759          -16.6                   -95.6
## # ℹ abbreviated name: ¹cambio_relativo_porcentaje

8 Series temporales

8.1 Serie temporal separada por conducta

Acá obtienen gráficos sobre el comportamiento de cada unidad comportamental a lo largo de todo el muestreo.

NOTA Mi gráfico de descanso está en 0 porque decidí no tenerlo en cuenta para mi análisis final, pero a ustedes les va a salir un gráfico para cada unidad comportamental que tengan en el excel.

La linea punteada en mi caso diferencia la linea base comportamental y la fase de enriquecimientos. Si su muestreo es de más de 10 días para cada fase, deben moverla editando el valor de xintercept en geom_vline

grafico_temporal <- ggplot(
  datos_largos,
  aes(
    x = dia,
    y = porcentaje,
    group = 1
  )
) +
  geom_line(
    linewidth = 0.7
  ) +
  geom_point(
    size = 2
  ) +

  geom_vline(
    xintercept = 10.5,
    linetype = "dashed",
    linewidth = 0.8
  ) +

  facet_wrap(
    ~ conducta,
    scales = "free_y",
    ncol = 2
  ) +

  scale_x_continuous(
    breaks = 1:20
  ) +

  labs(
    title = "Presupuesto de actividad de Potos flavus",
    subtitle = paste(
      "La línea discontinua indica el inicio",
      "del enriquecimiento ambiental"
    ),
    x = "Día de observación",
    y = "Porcentaje del tiempo (%)"
  ) +

  theme_classic(
    base_size = 12
  ) +

  theme(
    strip.text = element_text(
      face = "bold"
    ),

    axis.text.x = element_text(
      angle = 45,
      hjust = 1
    )
  )

grafico_temporal

8.2 Serie temporal con todas las conductas

Combina en un gráfico el resumen del comportamiento de todas las unidades comportamentales a lo largo del tiempo

Con ggplot pueden editar todo lo que quieran del gráfico: Títulos, colores, tamaños, ubicaciones, etc.

grafico_temporal_total <- ggplot(
  datos_largos,
  aes(
    x = dia,
    y = porcentaje,
    color = conducta,
    group = conducta
  )
) +

  geom_line(
    linewidth = 1
  ) +

  geom_point(
    size = 2
  ) +

  geom_vline(
    xintercept = 10.5,
    linetype = "dashed",
    linewidth = 0.8,
    color = "black"
  ) +

  scale_x_continuous(
    breaks = 1:20
  ) +

  scale_color_manual(
    values = c(
      "Alimentación" = "#b6d7a8",
      "Territorialidad" = "#a2c4c9",
      "Locomoción" = "#b4a7d6",
      "Descanso" = "#d5a6bd",
      "Mantenimiento" = "#9fc5e8",
      "Estereotipias" = "#f9cb9c"
    )
  ) +

  labs(
    title = "Presupuesto de actividad de Potos flavus",
    subtitle = paste(
      "La línea discontinua indica el inicio",
      "del enriquecimiento ambiental"
    ),
    x = "Día de observación",
    y = "Porcentaje del tiempo (%)",
    color = "Unidad comportamental"
  ) +

  theme_classic(
    base_size = 13
  ) +

  add_phylopic(
    img = potos,
    x = 19,
    y = 70,
    height = 20
  ) +

  theme(
    legend.position = "right",

    legend.text = element_text(
      size = 14
    ),

    legend.title = element_text(
      face = "bold"
    ),

    axis.text.x = element_text(
      angle = 45,
      hjust = 1
    )
  )

grafico_temporal_total

9 Medias e intervalos de confianza

Para este gráfico se excluye la conducta de descanso mediante conducta != "Descanso", Si quieren utilizar todas las conductas deben borrar esa parte

Con ggplot pueden editar todo lo que quieran del gráfico: Títulos, colores, tamaños, ubicaciones, etc.

resumen_descriptivo2 <- subset(
  resumen_descriptivo,
  conducta != "Descanso"
)

grafico_medias <- ggplot(
  resumen_descriptivo2,
  aes(
    x = conducta,
    y = media,
    fill = fase
  )
) +

  geom_col(
    position = position_dodge(
      width = 0.8
    ),
    width = 0.7
  ) +

  geom_errorbar(
    aes(
      ymin = pmax(
        ic_95_inferior,
        0
      ),
      ymax = ic_95_superior
    ),
    position = position_dodge(
      width = 0.8
    ),
    width = 0.2
  ) +

  geom_text(
    aes(
      label = paste0(
        round(media, 1),
        "%"
      )
    ),
    position = position_dodge(
      width = 0.8
    ),
    vjust = -0.5,
    size = 3.4
  ) +

  labs(
    title = "Presupuesto de actividad por fase",
    x = NULL,
    y = "Media del tiempo observado (%)",
    fill = NULL
  ) +

  theme_classic(
    base_size = 12
  ) +

  theme(
    axis.title.x = element_text(
      size = 16,
      face = "bold"
    ),

    axis.title.y = element_text(
      size = 16,
      face = "bold"
    ),

    axis.text.x = element_text(
      size = 14,
      angle = 0,
      hjust = 0.5
    ),

    axis.text.y = element_text(
      size = 14
    ),

    legend.text = element_text(
      size = 14
    ),

    legend.position = "top"
  )

grafico_medias

10 Boxplots

Recuerden que en mi análisis no tengo en cuenta descanso, por eso está en 0. Pero a ustedes si les salen correctamente todos los boxplots

Con ggplot pueden editar todo lo que quieran del gráfico: Títulos, colores, tamaños, ubicaciones, etc.

grafico_boxplot <- ggplot(
  datos_largos,
  aes(
    x = fase,
    y = porcentaje,
    fill = fase
  )
) +

  geom_boxplot(
    width = 0.6,
    outlier.shape = NA,
    alpha = 0.75
  ) +

  geom_jitter(
    width = 0.08,
    size = 2
  ) +

  facet_wrap(
    ~ conducta,
    scales = "free_y",
    ncol = 2
  ) +

  labs(
    title = "Distribución de las conductas por fase",
    x = NULL,
    y = "Porcentaje del tiempo (%)"
  ) +

  theme_classic(
    base_size = 12
  ) +

  theme(
    legend.position = "none",

    strip.text = element_text(
      face = "bold"
    )
  )

grafico_boxplot

11 Normalidad por conducta y fase

Se utiliza la prueba de Shapiro-Wilk cuando existen entre 3 y x observaciones, son independientes y medibles cuantitativamente.

El script ya está diseñado para que les arroje el resultado en la columna 5 evaluando automáticamente el valor de p : Sin evidencia contra normalidad, No normal o No evaluable

normalidad <- datos_largos %>%
  group_by(
    conducta,
    fase
  ) %>%
  summarise(
    n = sum(!is.na(porcentaje)),

    p_shapiro = if (
      n >= 3 &&
      n <= 5000 &&
      sd(
        porcentaje,
        na.rm = TRUE
      ) > 0
    ) {
      shapiro.test(
        porcentaje
      )$p.value
    } else {
      NA_real_
    },

    resultado = case_when(
      is.na(p_shapiro) ~
        "No evaluable",

      p_shapiro < 0.05 ~
        "No normal",

      TRUE ~
        "Sin evidencia contra normalidad"
    ),

    .groups = "drop"
  )

print(normalidad)

## # A tibble: 12 × 5
##    conducta        fase                n  p_shapiro resultado                   
##    <fct>           <fct>           <int>      <dbl> <chr>                       
##  1 Alimentación    Observacion        10  0.155     Sin evidencia contra normal…
##  2 Alimentación    Enriquecimiento    10  0.0427    No normal                   
##  3 Territorialidad Observacion        10  0.690     Sin evidencia contra normal…
##  4 Territorialidad Enriquecimiento    10  0.164     Sin evidencia contra normal…
##  5 Locomoción      Observacion        10  0.110     Sin evidencia contra normal…
##  6 Locomoción      Enriquecimiento    10  0.603     Sin evidencia contra normal…
##  7 Descanso        Observacion        10 NA         No evaluable                
##  8 Descanso        Enriquecimiento    10 NA         No evaluable                
##  9 Mantenimiento   Observacion        10  0.129     Sin evidencia contra normal…
## 10 Mantenimiento   Enriquecimiento    10  0.0000738 No normal                   
## 11 Estereotipias   Observacion        10  0.140     Sin evidencia contra normal…
## 12 Estereotipias   Enriquecimiento    10  0.00126   No normal

12 Homogeneidad de varianzas

Se aplica la prueba de Levene para cada conducta suponiendo que son independientes y medibles cuantitativamente.

El script ya está diseñado para que les arroje el resultado en la columna 3 evaluando automáticamente el valor de p : Varianzas no homogéneas o Sin evidencia de heterogeneidad

homogeneidad <- datos_largos %>%
  group_by(conducta) %>%
  group_modify(
    ~ {
      prueba <- car::leveneTest(
        porcentaje ~ fase,
        data = .x
      )

      tibble(
        p_levene = prueba$`Pr(>F)`[1]
      )
    }
  ) %>%
  mutate(
    resultado = ifelse(
      p_levene < 0.05,
      "Varianzas no homogéneas",
      "Sin evidencia de heterogeneidad"
    )
  )

print(homogeneidad)

## # A tibble: 6 × 3
## # Groups:   conducta [6]
##   conducta         p_levene resultado                      
##   <fct>               <dbl> <chr>                          
## 1 Alimentación      0.0319  Varianzas no homogéneas        
## 2 Territorialidad   0.234   Sin evidencia de heterogeneidad
## 3 Locomoción        0.0320  Varianzas no homogéneas        
## 4 Descanso        NaN       <NA>                           
## 5 Mantenimiento     0.0505  Sin evidencia de heterogeneidad
## 6 Estereotipias     0.00447 Varianzas no homogéneas

13 Prueba de Wilcoxon

La prueba compara las fases de observación y enriquecimiento para cada conducta. Se usa si no hay normalidad en los datos

Utilicé la corrección de Holm method = "holm" para controlar falsos positivos por múltiples comparaciones.

Al utilizar paired = TRUE, ambas fases deben tener el mismo número de observaciones y los registros deben estar en el orden correcto para que cada dato de la fase A corresponda con un dato de la fase B. Si no son datos pareados como en mi caso, paired = FALSE

wilcoxon_resultados <- datos_largos %>%
  group_by(conducta) %>%
  wilcox_test(
    porcentaje ~ fase,
    paired = FALSE,
    exact = FALSE
  ) %>%
  adjust_pvalue(
    method = "holm"
  ) %>%
  add_significance(
    "p.adj"
  )

print(wilcoxon_resultados)

## # A tibble: 6 × 10
##   conducta        .y.    group1 group2    n1    n2 statistic         p     p.adj
##   <fct>           <chr>  <chr>  <chr>  <int> <int>     <dbl>     <dbl>     <dbl>
## 1 Alimentación    porce… Obser… Enriq…    10    10        38   0.385     1      
## 2 Territorialidad porce… Obser… Enriq…    10    10        52   0.910     1      
## 3 Locomoción      porce… Obser… Enriq…    10    10        48   0.910     1      
## 4 Descanso        porce… Obser… Enriq…    10    10        50 NaN       NaN      
## 5 Mantenimiento   porce… Obser… Enriq…    10    10        71   0.121     0.485  
## 6 Estereotipias   porce… Obser… Enriq…    10    10        93   0.00103   0.00514
## # ℹ 1 more variable: p.adj.signif <chr>

14 NAP

El índice NAP (Non-overlap of All Pairs) representa la proporción de comparaciones entre las fases A y B que presentan un cambio en la dirección esperada.

14.1 Función para calcular NAP

Acá NO deben editar nada

calcular_nap <- function(
    fase_a,
    fase_b,
    direccion = "aumento"
) {

  fase_a <- fase_a[
    !is.na(fase_a)
  ]

  fase_b <- fase_b[
    !is.na(fase_b)
  ]

  if (
    length(fase_a) == 0 ||
    length(fase_b) == 0
  ) {
    return(NA_real_)
  }

  comparaciones <- outer(
    fase_b,
    fase_a,
    FUN = "-"
  )

  if (direccion == "aumento") {
    favorables <- sum(
      comparaciones > 0
    )

    empates <- sum(
      comparaciones == 0
    )
  } else if (direccion == "disminucion") {
    favorables <- sum(
      comparaciones < 0
    )

    empates <- sum(
      comparaciones == 0
    )
  } else {
    stop(
      paste(
        "La dirección debe ser",
        "'aumento' o 'disminucion'."
      )
    )
  }

  total <- length(fase_a) *
    length(fase_b)

  nap <- (
    favorables +
      0.5 * empates
  ) / total

  return(nap)
}

14.2 NAP para todas las conductas

Se debe ajustar la dirección aumento o disminucion de acuerdo a lo que esperaban con cada enriquecimiento para cada unidad comportamental.

El script ya está diseñado para que les arroje la interpretación en la columna 3 evaluando automáticamente el valor de NAP : Débil, Moderado, Fuerte o Contrario al esperado

direcciones <- tibble(
  conducta = factor(
    levels(datos_largos$conducta),
    levels = levels(datos_largos$conducta)
  ),

  direccion = c(
    "aumento",     # Alimentación
    "aumento",     # Territorialidad
    "aumento",     # Locomoción
    "aumento",     # Descanso
    "aumento",     # Mantenimiento
    "disminucion"  # Estereotipias
  )
)

resultados_nap <- datos_largos %>%
  left_join(
    direcciones,
    by = "conducta"
  ) %>%
  group_by(
    conducta,
    direccion
  ) %>%
  summarise(
    nap = calcular_nap(
      fase_a = porcentaje[
        fase == "Observacion"
      ],

      fase_b = porcentaje[
        fase == "Enriquecimiento"
      ],

      direccion = first(direccion)
    ),

    .groups = "drop"
  ) %>%
  mutate(
    interpretacion = case_when(
      is.na(nap) ~
        "No evaluable",

      nap < 0.50 ~
        "Cambio contrario al esperado",

      nap < 0.66 ~
        "Efecto débil",

      nap < 0.92 ~
        "Efecto moderado",

      TRUE ~
        "Efecto fuerte"
    )
  )

print(resultados_nap)

## # A tibble: 6 × 4
##   conducta        direccion     nap interpretacion              
##   <fct>           <chr>       <dbl> <chr>                       
## 1 Alimentación    aumento      0.62 Efecto débil                
## 2 Territorialidad aumento      0.48 Cambio contrario al esperado
## 3 Locomoción      aumento      0.52 Efecto débil                
## 4 Descanso        aumento      0.5  Efecto débil                
## 5 Mantenimiento   aumento      0.29 Cambio contrario al esperado
## 6 Estereotipias   disminucion  0.93 Efecto fuerte

14.3 Gráfico de NAP

Con ggplot pueden editar todo lo que quieran del gráfico: Títulos, colores, tamaños, ubicaciones, etc.

colores_nap <- c(
  "Alimentación" = "#b6d7a8",
  "Descanso" = "#d5a6bd",
  "Estereotipias" = "#f9cb9c",
  "Locomoción" = "#b4a7d6",
  "Mantenimiento" = "#9fc5e8",
  "Territorialidad" = "#a2c4c9"
)

grafico_nap <- ggplot(
  resultados_nap,
  aes(
    x = reorder(conducta, nap),
    y = nap,
    fill = conducta
  )
) +

  geom_col(
    width = 0.7
  ) +

  scale_fill_manual(
    values = colores_nap
  ) +

  geom_hline(
    yintercept = 0.50,
    linetype = "dashed"
  ) +

  geom_hline(
    yintercept = 0.75,
    linetype = "dotted"
  ) +

  geom_hline(
    yintercept = 1,
    linetype = "dotted"
  ) +

  geom_text(
    aes(
      label = round(nap, 2)
    ),
    hjust = -0.2
  ) +

  coord_flip() +

  scale_y_continuous(
    limits = c(0, 1.08),
    breaks = seq(
      0,
      1,
      0.2
    )
  ) +

  add_phylopic(
    img = potos,
    x = 1.5,
    y = 0.9,
    height = 1.3
  ) +

  labs(
    title = "Efecto del enriquecimiento ambiental",
    subtitle = "Non-overlap of All Pairs (NAP)",
    x = NULL,
    y = "NAP"
  ) +

  theme_classic(
    base_size = 12
  ) +

  theme(
    legend.position = "none"
  )

grafico_nap

15 Función para Tau-U entre fases

Acá NO deben editar nada

calcular_tau_u <- function(
    fase_a,
    fase_b,
    direccion = "aumento"
) {

  fase_a <- fase_a[
    !is.na(fase_a)
  ]

  fase_b <- fase_b[
    !is.na(fase_b)
  ]

  if (
    length(fase_a) == 0 ||
    length(fase_b) == 0
  ) {
    return(NA_real_)
  }

  diferencias <- outer(
    fase_b,
    fase_a,
    FUN = "-"
  )

  if (direccion == "aumento") {
    favorables <- sum(
      diferencias > 0
    )

    desfavorables <- sum(
      diferencias < 0
    )
  } else if (direccion == "disminucion") {
    favorables <- sum(
      diferencias < 0
    )

    desfavorables <- sum(
      diferencias > 0
    )
  } else {
    stop(
      paste(
        "La dirección debe ser",
        "'aumento' o 'disminucion'."
      )
    )
  }

  total <- length(fase_a) *
    length(fase_b)

  tau_u <- (
    favorables -
      desfavorables
  ) / total

  return(tau_u)
}

15.1 Tau-U para todas las conductas

El script ya está diseñado para que les arroje la interpretación en la columna 3 evaluando automáticamente el valor de tau-u : Pequeño, Moderado, Grande, Muy grande o Contrario al esperado

resultados_tau <- datos_largos %>%
  left_join(
    direcciones,
    by = "conducta"
  ) %>%
  group_by(
    conducta,
    direccion
  ) %>%
  summarise(
    tau_u = calcular_tau_u(
      fase_a = porcentaje[
        fase == "Observacion"
      ],

      fase_b = porcentaje[
        fase == "Enriquecimiento"
      ],

      direccion = first(direccion)
    ),

    .groups = "drop"
  ) %>%
  mutate(
    interpretacion = case_when(
      is.na(tau_u) ~
        "No evaluable",

      tau_u < 0 ~
        "Cambio contrario al esperado",

      tau_u < 0.20 ~
        "Efecto pequeño",

      tau_u < 0.60 ~
        "Efecto moderado",

      tau_u < 0.80 ~
        "Efecto grande",

      TRUE ~
        "Efecto muy grande"
    )
  )

print(resultados_tau)

## # A tibble: 6 × 4
##   conducta        direccion   tau_u interpretacion              
##   <fct>           <chr>       <dbl> <chr>                       
## 1 Alimentación    aumento      0.24 Efecto moderado             
## 2 Territorialidad aumento     -0.04 Cambio contrario al esperado
## 3 Locomoción      aumento      0.04 Efecto pequeño              
## 4 Descanso        aumento      0    Efecto pequeño              
## 5 Mantenimiento   aumento     -0.42 Cambio contrario al esperado
## 6 Estereotipias   disminucion  0.86 Efecto muy grande

16 PERMANOVA

Se usa si no hay normalidad multivariante.

La PERMANOVA evalúa si la composición general del presupuesto de actividad difiere entre las fases.

matriz_conductas <- datos %>%
  dplyr::select(
    all_of(conductas)
  )

distancia_bray <- vegdist(
  matriz_conductas,
  method = "bray"
)

set.seed(123)

permanova <- adonis2(
  distancia_bray ~ fase,
  data = datos,
  permutations = 9999
)

print(permanova)

## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 9999
## 
## adonis2(formula = distancia_bray ~ fase, data = datos, permutations = 9999)
##          Df SumOfSqs      R2      F Pr(>F)  
## Model     1  0.16707 0.14866 3.1431  0.046 *
## Residual 18  0.95680 0.85134                
## Total    19  1.12387 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

16.1 Homogeneidad de dispersiones multivariadas

Esta prueba permite comprobar el supuesto de homogeneidad de dispersiones para interpretar la PERMANOVA.

dispersion <- betadisper(
  distancia_bray,
  group = datos$fase
)

anova_dispersion <- anova(
  dispersion
)

print(anova_dispersion)

## Analysis of Variance Table
## 
## Response: Distances
##           Df   Sum Sq   Mean Sq F value Pr(>F)
## Groups     1 0.010222 0.0102218  1.0499 0.3191
## Residuals 18 0.175253 0.0097363

16.2 Prueba de permutaciones

set.seed(123)

permutest_dispersion <- permutest(
  dispersion,
  permutations = 9999
)

print(permutest_dispersion)

## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 9999
## 
## Response: Distances
##           Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
## Groups     1 0.010222 0.0102218 1.0499   9999 0.3246
## Residuals 18 0.175253 0.0097363

Si la prueba de dispersión presenta un valor de \(p > 0.05\), no existe evidencia fuerte de diferencias en la variabilidad interna de las fases.

Si presenta un valor de \(p < 0.05\), la significancia de la PERMANOVA podría estar relacionada parcialmente con diferencias de dispersión y no únicamente con diferencias entre los centroides.

17 PCoA con distancia de Bray-Curtis

pcoa <- cmdscale(
  distancia_bray,
  eig = TRUE,
  k = 2
)

porcentaje_explicado <- round(
  100 *
    pcoa$eig[1:2] /
    sum(
      pcoa$eig[
        pcoa$eig > 0
      ]
    ),
  1
)

datos_pcoa <- data.frame(
  dia = datos$dia,
  fase = datos$fase,
  eje_1 = pcoa$points[, 1],
  eje_2 = pcoa$points[, 2]
)

grafico_pcoa <- ggplot(
  datos_pcoa,
  aes(
    x = eje_1,
    y = eje_2,
    shape = fase
  )
) +

  geom_point(
    size = 3
  ) +

  geom_text(
    aes(
      label = dia
    ),
    nudge_y = 0.01,
    size = 3
  ) +

  stat_ellipse(
    aes(
      group = fase
    ),
    type = "norm",
    linewidth = 0.7,
    level = 0.80
  ) +

  labs(
    title = "Ordenación del presupuesto de actividad",
    subtitle = "PCoA basada en distancia de Bray-Curtis",

    x = paste0(
      "Eje 1 (",
      porcentaje_explicado[1],
      "%)"
    ),

    y = paste0(
      "Eje 2 (",
      porcentaje_explicado[2],
      "%)"
    ),

    shape = "Fase"
  ) +

  theme_classic(
    base_size = 12
  )

grafico_pcoa

18 PCA

##Transformación CLR

La transformación CLR se emplea para analizar datos composicionales. Los valores iguales o inferiores a cero se reemplazan mediante un pseudoconteo.

transformar_clr <- function(
    x,
    pseudoconteo = 0.01
) {

  x <- as.matrix(x)

  # Reemplazar ceros
  x[x <= 0] <- pseudoconteo

  # Convertir a proporciones por fila
  x <- x / rowSums(x)

  # Calcular logaritmos
  log_x <- log(x)

  # Restar la media logarítmica de cada fila
  clr <- log_x - rowMeans(log_x)

  return(clr)
}

matriz_clr <- transformar_clr(
  datos %>%
    dplyr::select(
      all_of(conductas)
    ),
  pseudoconteo = 0.01
)

colnames(matriz_clr) <- conductas

18.1 PCA de los datos composicionales

pca_clr <- prcomp(
  matriz_clr,
  center = TRUE,
  scale. = FALSE
)

summary(pca_clr)

## Importance of components:
##                           PC1     PC2    PC3     PC4     PC5       PC6
## Standard deviation     3.2304 0.66884 0.4905 0.45328 0.13068 5.597e-16
## Proportion of Variance 0.9198 0.03943 0.0212 0.01811 0.00151 0.000e+00
## Cumulative Proportion  0.9198 0.95918 0.9804 0.99849 1.00000 1.000e+00

18.1.1 Varianza explicada

varianza_pca <- round(
  100 * (
    pca_clr$sdev^2 /
      sum(pca_clr$sdev^2)
  ),
  1
)

datos_pca <- data.frame(
  dia = datos$dia,
  fase = datos$fase,
  CP1 = pca_clr$x[, 1],
  CP2 = pca_clr$x[, 2]
)

18.2 Gráfico del PCA

grafico_pca <- ggplot(
  datos_pca,
  aes(
    x = CP1,
    y = CP2,
    shape = fase
  )
) +

  geom_point(
    size = 3
  ) +

  geom_text(
    aes(
      label = dia
    ),
    nudge_y = 0.15,
    size = 3
  ) +

  stat_ellipse(
    aes(
      group = fase
    ),
    level = 0.80,
    type = "norm"
  ) +

  labs(
    title = "PCA del presupuesto de actividad",
    subtitle = "Datos transformados mediante CLR",

    x = paste0(
      "Componente 1 (",
      varianza_pca[1],
      "%)"
    ),

    y = paste0(
      "Componente 2 (",
      varianza_pca[2],
      "%)"
    ),

    shape = "Fase"
  ) +

  theme_classic(
    base_size = 12
  )

grafico_pca

18.3 Cargas de los componentes principales

Las cargas permiten identificar cuáles conductas contribuyen principalmente a cada componente.

cargas_pca <- as.data.frame(
  pca_clr$rotation[, 1:2]
) %>%
  tibble::rownames_to_column(
    "conducta"
  )

print(cargas_pca)

##                    conducta         PC1         PC2
## 1              alimentacion -0.26720554 -0.04527623
## 2           territorialidad -0.20951401 -0.24654041
## 3                locomocion -0.09582792  0.69534330
## 4                  descanso -0.13821140  0.28444013
## 5 mantenimiento_autocuidado -0.19408877 -0.60485717
## 6             estereotipias  0.90484764 -0.08310962

19 Cartas de individuos por conducta

Acá pueden editar o agregar condutas de acuerdo a sus unidades comportamentales.

crear_carta_individuos <- function(
    variable,
    nombre_conducta
) {

  qcc(
    data = variable,
    type = "xbar.one",
    labels = datos$dia,
    title = paste(
      "Carta de individuos:",
      nombre_conducta
    ),
    xlab = "Día",
    ylab = "Porcentaje del tiempo"
  )
}

19.1 Alimentación

crear_carta_individuos(
  datos$alimentacion,
  "Alimentación"
)

## List of 11
##  $ call      : language qcc(data = variable, type = "xbar.one", labels = datos$dia, title = paste("Carta de individuos:",      nombre_con| __truncated__
##  $ type      : chr "xbar.one"
##  $ data.name : chr "variable"
##  $ data      : num [1:20, 1] 53.8 23 17.5 40 15.8 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ statistics: Named num [1:20] 53.8 23 17.5 40 15.8 ...
##   ..- attr(*, "names")= chr [1:20] "1" "2" "3" "4" ...
##  $ sizes     : int [1:20] 1 1 1 1 1 1 1 1 1 1 ...
##  $ center    : num 28.8
##  $ std.dev   : num 12.6
##  $ nsigmas   : num 3
##  $ limits    : num [1, 1:2] -8.84 66.53
##   ..- attr(*, "dimnames")=List of 2
##  $ violations:List of 2
##  - attr(*, "class")= chr "qcc"

19.2 Territorialidad

crear_carta_individuos(
  datos$territorialidad,
  "Territorialidad"
)

## List of 11
##  $ call      : language qcc(data = variable, type = "xbar.one", labels = datos$dia, title = paste("Carta de individuos:",      nombre_con| __truncated__
##  $ type      : chr "xbar.one"
##  $ data.name : chr "variable"
##  $ data      : num [1:20, 1] 19.5 22 18 13.3 13.6 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ statistics: Named num [1:20] 19.5 22 18 13.3 13.6 ...
##   ..- attr(*, "names")= chr [1:20] "1" "2" "3" "4" ...
##  $ sizes     : int [1:20] 1 1 1 1 1 1 1 1 1 1 ...
##  $ center    : num 14.2
##  $ std.dev   : num 7.66
##  $ nsigmas   : num 3
##  $ limits    : num [1, 1:2] -8.77 37.21
##   ..- attr(*, "dimnames")=List of 2
##  $ violations:List of 2
##  - attr(*, "class")= chr "qcc"

19.3 Locomoción

crear_carta_individuos(
  datos$locomocion,
  "Locomoción"
)

## List of 11
##  $ call      : language qcc(data = variable, type = "xbar.one", labels = datos$dia, title = paste("Carta de individuos:",      nombre_con| __truncated__
##  $ type      : chr "xbar.one"
##  $ data.name : chr "variable"
##  $ data      : num [1:20, 1] 19.3 42.8 44 39.1 53.3 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ statistics: Named num [1:20] 19.3 42.8 44 39.1 53.3 ...
##   ..- attr(*, "names")= chr [1:20] "1" "2" "3" "4" ...
##  $ sizes     : int [1:20] 1 1 1 1 1 1 1 1 1 1 ...
##  $ center    : num 46.1
##  $ std.dev   : num 9.71
##  $ nsigmas   : num 3
##  $ limits    : num [1, 1:2] 16.9 75.2
##   ..- attr(*, "dimnames")=List of 2
##  $ violations:List of 2
##  - attr(*, "class")= chr "qcc"

19.4 Mantenimiento

crear_carta_individuos(
  datos$mantenimiento_autocuidado,
  "Mantenimiento"
)

## List of 11
##  $ call      : language qcc(data = variable, type = "xbar.one", labels = datos$dia, title = paste("Carta de individuos:",      nombre_con| __truncated__
##  $ type      : chr "xbar.one"
##  $ data.name : chr "variable"
##  $ data      : num [1:20, 1] 4.095 4.443 1.291 2.381 0.386 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ statistics: Named num [1:20] 4.095 4.443 1.291 2.381 0.386 ...
##   ..- attr(*, "names")= chr [1:20] "1" "2" "3" "4" ...
##  $ sizes     : int [1:20] 1 1 1 1 1 1 1 1 1 1 ...
##  $ center    : num 1.84
##  $ std.dev   : num 1.21
##  $ nsigmas   : num 3
##  $ limits    : num [1, 1:2] -1.8 5.49
##   ..- attr(*, "dimnames")=List of 2
##  $ violations:List of 2
##  - attr(*, "class")= chr "qcc"

19.5 Estereotipias

crear_carta_individuos(
  datos$estereotipias,
  "Estereotipias"
)

## List of 11
##  $ call      : language qcc(data = variable, type = "xbar.one", labels = datos$dia, title = paste("Carta de individuos:",      nombre_con| __truncated__
##  $ type      : chr "xbar.one"
##  $ data.name : chr "variable"
##  $ data      : num [1:20, 1] 3.34 7.69 19.17 5.12 16.93 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ statistics: Named num [1:20] 3.34 7.69 19.17 5.12 16.93 ...
##   ..- attr(*, "names")= chr [1:20] "1" "2" "3" "4" ...
##  $ sizes     : int [1:20] 1 1 1 1 1 1 1 1 1 1 ...
##  $ center    : num 9.04
##  $ std.dev   : num 6.96
##  $ nsigmas   : num 3
##  $ limits    : num [1, 1:2] -11.8 29.9
##   ..- attr(*, "dimnames")=List of 2
##  $ violations:List of 2
##  - attr(*, "class")= chr "qcc"

Análisis del presupuesto de actividad URRAS

Ejemplo de Potos flavus 25552

Santiago Ballesteros Rocha

26/08/2026

1 Introducción

2 Instalación y carga de paquetes

2.1 Silueta de Potos flavus

3 Importación de los datos

4 Limpieza y organización

4.1 Revisión de la base de datos

4.1.1 Nombres de las variables

4.1.2 Primeros registros

4.1.3 Estructura de la base

4.2 Resumen general

4.3 Definición de las conductas

4.4 Verificación de columnas

4.5 Verificación de datos faltantes

4.6 Verificación de los porcentajes

5 Transformación a formato largo

6 Estadística descriptiva

7 Cambios entre las fases

8 Series temporales

8.1 Serie temporal separada por conducta

8.2 Serie temporal con todas las conductas

9 Medias e intervalos de confianza

10 Boxplots

11 Normalidad por conducta y fase

12 Homogeneidad de varianzas

13 Prueba de Wilcoxon

14 NAP

14.1 Función para calcular NAP

14.2 NAP para todas las conductas

14.3 Gráfico de NAP

15 Función para Tau-U entre fases

15.1 Tau-U para todas las conductas

16 PERMANOVA

16.1 Homogeneidad de dispersiones multivariadas

16.2 Prueba de permutaciones

17 PCoA con distancia de Bray-Curtis

18 PCA

18.1 PCA de los datos composicionales

18.1.1 Varianza explicada

18.2 Gráfico del PCA

18.3 Cargas de los componentes principales

19 Cartas de individuos por conducta

19.1 Alimentación

19.2 Territorialidad

19.3 Locomoción

19.4 Mantenimiento

19.5 Estereotipias