Variación espacial y temporal en la talla de Nerita scabricosta en áreas protegidas del Pacífico norte de Costa Rica
Bioestadística General I · Proyecto de investigación · Los Burgaditos
Valeria Castro-Madriz, Delany Medina-Ureña, Shariel Vargas-Soto, Alberto González-Soto e Itzel Espinoza-Quirós
12 de noviembre de 2026
Profesora: Karol Ulate-Naranjo
Debajo de cada resultado, el enlace Ver comandos y salida de R despliega los comandos con los que se hizo cada prueba estadística y su salida. En el Anexo B están los scripts completos.
Resumen
Los burgados (Nerita scabricosta) son caracoles del intermareal rocoso alto que se recolectan como alimento y carnada, por lo que su tamaño puede reflejar la presión de extracción. Este estudio evaluó si el largo y el ancho de la concha difieren entre cuatro áreas silvestres protegidas del Área de Conservación Tempisque (Costa Rica), entre años y épocas climáticas, y en relación con la temperatura, la precipitación y el semáforo del Programa Nacional de Monitoreo Ecológico (PRONAMEC). Se midieron 1 019 individuos en 36 muestreos entre 2022 y 2026. El largo medio fue de 2.08 cm. La Reserva Natural Absoluta Cabo Blanco, donde no se permite la extracción, presentó los burgados más grandes, aunque el tamaño del efecto fue pequeño (ε² = 0.01) y la diferencia no se mantuvo al usar el muestreo como unidad de análisis (p = 0.740). El tamaño varió más entre años que entre áreas: dentro de la Reserva, el año explicó casi la mitad de la variación del largo (ε² = 0.48). El ANOVA de dos vías para 2025 y 2026 mostró una interacción significativa entre área y año (η² parcial = 0.16), pero el año coincidió con la época de muestreo, por lo que ambos factores no pudieron separarse. El largo medio no se asoció con la temperatura ni con la precipitación de los 90 días previos, y el ancho siguió un patrón similar. Según el semáforo del PRONAMEC, 26 de 34 muestreos fueron «regular», 4 «bueno», 4 quedaron sin categoría y ninguno fue «malo». En conjunto, el tamaño de los burgados dependió más de las características del sitio y del momento del muestreo que del tipo de manejo, por lo que se recomienda un monitoreo balanceado entre sitios, años y épocas, y calibrar los intervalos del protocolo con datos locales.
Palabras clave: gasterópodos, áreas silvestres protegidas, talla de burgados, intermareal rocoso, PRONAMEC.
Abstract
Burgados (Nerita scabricosta) are snails found in the upper rocky intertidal zone that are collected for food and bait; therefore, their size may reflect harvesting pressure. This study evaluated whether shell length and width differ across four protected wilderness areas within the Tempisque Conservation Area (Costa Rica), across years and seasons, and in relation to temperature, precipitation, and the traffic-light system of the National Ecological Monitoring Program (PRONAMEC). A total of 1,019 individuals were measured during 36 sampling events between 2022 and 2026. The mean length was 2.08 cm. The Cabo Blanco Absolute Nature Reserve, where harvesting is prohibited, had the largest snails, although the effect size was small (ε² = 0.01) and the difference was not significant when using the sampling event as the unit of analysis (p = 0.740). Size varied more between years than between areas: within the Reserve, the year accounted for nearly half of the variation in length (ε² = 0.48). The two-way ANOVA for 2025 and 2026 showed a significant interaction between area and year (partial η² = 0.16), but the year coincided with the sampling season, so the two factors could not be separated. The average length was not associated with temperature or precipitation during the preceding 90 days, and the width followed a similar pattern. According to the PRONAMEC traffic-light system, 26 of 34 samples were rated “fair,” 4 were “good,” 4 were unrated, and none were “poor.” Overall, the size of burgados depended more on site characteristics and the timing of the sampling than on management type; therefore, it is recommended to conduct balanced monitoring across sites, years, and seasons, and to calibrate the protocol intervals using local data.
Keywords: gastropods, protected areas, burgado size, rocky intertidal zone, PRONAMEC.
Introducción
Cuando la marea baja, el estrato rocoso queda al descubierto, convirtiéndose en uno de los ambientes más desafiantes de la costa, puesto que se presenta una mayor exposición a la radiación solar, olas con mayor fuerza de impacto, así como el cambiante nivel del agua durante el día (Kunze et al., 2021; Leeuwis & Gamperl, 2022). Según Baliwe et al. (2022), a pesar de dichas condiciones, en estos hábitats se pueden encontrar caracoles, mejillones, algas, entre otros organismos, que más allá de sostener la vida del litoral, representan un recurso importante para las comunidades costeras (Faria et al., 2024; Gutiérrez-Rubio & Arzola-González, 2024).
Es por ello que el tamaño del cuerpo de los gasterópodos es utilizado como una señal del estado de las poblaciones explotadas. En un área protegida de Sudáfrica, las zonas donde está prohibida la extracción de caracoles de diversas especies de interés comercial, presentaron una conservación en cuanto a las densidades y las estructuras de sus tallas naturales (Baliwe et al., 2022). En las islas Azores, los moluscos de la especie con mayor interés fueron más grandes dentro de un área sin acceso humano, en contraparte, la especie menos valorada no mostró diferencias significativas en su tallaje (Faria et al., 2024). Ambos ejemplos sugieren que estas respuestas dependen de la especie y de las características del lugar de muestreo, a su vez, que el medir la forma de las conchas pueden ayudar a determinar los efectos de la exposición de estos organismos, a la intervención humana prolongada a lo largo de los años (Doyle et al., 2022).
El burgado Nerita scabricosta, es un caracol del intermareal alto, el cual posee una amplia distribución en el Pacífico oriental, donde recientemente se ha documentado una expansión de su rango geográfico asociada a procesos de tropicalización (Zarzyczny et al., 2026). Debido a su tamaño, peso y abundancia es considerado como un recurso pesquero potencial, extrayéndose como alimento y como carnada (Chacón-Monge et al., 2019; Gutiérrez-Rubio & Arzola-González, 2024).
Sin embargo, la talla de estos gasterópodos no depende únicamente de la presión por su recolección. Factores físicos como el oleaje influyen en la estructura de tallas en costas rocosas (Rodríguez-Grimón et al., 2024), al igual que la temperatura, dado que en estos ambientes las rocas pueden superar los 50 °C y elevar la temperatura corporal de los caracoles cerca de ese umbral (Marshall et al., 2023). Para evitar la mortalidad, estos organismos buscan microhábitats más frescos y ajustan su tolerancia térmica de acuerdo con las condiciones locales de cada sitio (Hui et al., 2022; Hui et al., 2025). En adición, las variaciones climáticas pueden alterar el reclutamiento de juveniles; por ejemplo, en Lottia gigantea, se ha observado que las anomalías cálidas incrementan la llegada de nuevos individuos, y con ello, la proporción de individuos pequeños (Nielsen et al., 2025).
Por otra parte, la época del año también puede influir en la talla en direcciones opuestas. Durante la época lluviosa, la mayor disponibilidad de agua reduce el estrés por desecación en el intermareal alto, condiciones que pueden favorecer la actividad de forrajeo y el crecimiento de los gasterópodos, ya que la desecación se asocia con pérdidas de peso y menor supervivencia en caracoles litorínidos (Maia & Troncoso, 2024). No obstante, debido a que N. scabricosta se reproduce exclusivamente en ese período (Collin et al., 2017), la llegada de juveniles podría disminuir la talla media de la población. Por estas razones, al evaluar las variaciones en la talla se deben considerar no solo el nivel de protección del área de estudio, sino también las características locales del sitio, la variabilidad interanual y las condiciones climáticas y estacionales.
En Costa Rica, el Sistema Nacional de Áreas de Conservación (SINAC) cuenta con el Programa Nacional de Monitoreo Ecológico (PRONAMEC), creado con el fin de generar y difundir información científica confiable que respalde las decisiones de manejo de la biodiversidad (Fallas-Flores et al., 2024). Uno de sus indicadores para el litoral rocoso es el largo medio de la concha de los burgados, que se clasifica mediante un sistema de semáforo como bueno, regular o malo según el grado de extracción que sugiere (SINAC, 2016). El Área de Conservación Tempisque (ACT), en el Pacífico Norte, alberga áreas silvestres protegidas con distinto régimen de manejo, desde una reserva natural absoluta en la que no se permite ninguna extracción hasta refugios de vida silvestre y un área marina de manejo con usos regulados, lo que permite comparar poblaciones bajo distintos niveles de protección. Con base en ello, el presente estudio tuvo como objetivo evaluar si el tamaño de la concha de N. scabricosta difiere entre cuatro áreas silvestres protegidas del ACT, entre años y entre épocas, si se asocia con la temperatura y la precipitación previas al muestreo, y cómo se clasifican sus poblaciones según el indicador del PRONAMEC. Se plantearon las hipótesis de que la Reserva Natural Absoluta Cabo Blanco, donde no se permite la extracción, presentaría individuos de mayor tamaño que en las áreas con uso regulado; que el tamaño sería mayor en periodos con mayor precipitación y menor temperatura, debido a un menor estrés por desecación, y que además diferirían entre épocas, sin una dirección definida.
Metodología
Área de estudio: El estudio se realizó en el litoral rocoso de la costa pacífica de la península de Nicoya (Guanacaste y Puntarenas, Costa Rica), dentro del Área de Conservación Tempisque (ACT). Se muestrearon 12 sitios distribuidos en cuatro áreas silvestres protegidas con distinto régimen de manejo: (a) la Reserva Natural Absoluta Cabo Blanco (en adelante, Reserva Cabo Blanco), el área con mayor protección, con los sitios San Miguel, Punta Flor y Playa Cabo Blanco, donde se excluye toda actividad extractiva (Decreto Ley N.º 3763, 1966; Alvarado et al., 2012); (b) el Área Marina de Manejo Cabo Blanco (en adelante, Área Marina Cabo Blanco), con los sitios Cabuya, Tambor y Playa Cuestas, donde la pesca y el turismo están regulados (Decreto Ejecutivo N.º 40442-MINAE, 2017); (c) el Refugio Nacional de Vida Silvestre Camaronal, con los sitios Boca Ora, Camaronal Este y Carrillo (Decreto Ejecutivo N.º 23640-MIRENEM, 1994); y (d) el Refugio Nacional de Vida Silvestre Ostional, con los sitios La India, Mojón 69 y San Juanillo (Ley N.º 9348, 2016). En la Reserva Cabo Blanco la pesca no está permitida (Arias et al., 2015) y su área marina está excluida del Área Marina de Manejo (Decreto Ejecutivo N.º 40442-MINAE, 2017, art. 2); en el Área Marina Cabo Blanco la pesca está permitida con licencia (Decreto Ejecutivo N.º 40442-MINAE, 2017, arts. 7 y 9), y en los refugios de Camaronal y Ostional también está permitida (Arias et al., 2015). En el Área Marina Cabo Blanco habitan comunidades de pesca artesanal (Cabuya, Tambor y Manzanillo), y en el primer kilómetro desde la pleamar ordinaria solo se permite la pesca con cuerda de mano (SINAC, 2017, como se citó en Solís Rivera & Fonseca Borrás, 2018). Esta descripción del régimen de uso no constituye una medición de la intensidad de vigilancia ni de la presión de extracción.
Colecta de datos biológicos: Los muestreos se realizaron en 2022, 2023, 2025 y 2026, durante las épocas seca (diciembre a abril) y lluviosa (mayo a noviembre), con una periodicidad irregular. La colecta siguió el Indicador 2 (tallas de burgados) del Protocolo para el Monitoreo Ecológico de las Playas Rocosas (PRONAMEC; SINAC, 2016). Cada muestreo se realizó durante bajamares de sicigia, entre una hora antes y una hora después del pico de marea baja. Según el protocolo, en cada sitio y fecha se seleccionan al azar 30 individuos de Nerita scabricosta en el intermareal rocoso alto, tomados de distintas agregaciones para no sobrerrepresentar un solo grupo. A cada individuo se le midió el largo y el ancho de la concha con un calibrador Vernier, con una precisión de 0.01 cm. La especie se identificó según el protocolo PRONAMEC, que define este indicador para N. scabricosta, la especie de Nerita dominante y sujeta a extracción en el intermareal alto del Pacífico costarricense (Azofeifa-Solano et al., 2018).
Depuración de la base de datos: La base inicial contenía 1 596 registros. Se excluyeron 328 registros copiados dos veces, 60 registros digitados en milímetros que duplicaban registros ya presentes en centímetros, 64 registros de otras especies o sin medición, cinco registros con largo menor a 0.5 cm, incoherentes con el resto de su muestreo, y 120 registros de cuatro muestreos de 2022 en los que el ancho superó al largo de forma sistemática (razón ancho/largo mediana > 1), lo que indica un criterio de medición distinto. La base final incluyó 1 019 individuos. Además, se anularon 34 valores de ancho y se conservó el largo de esos individuos: 30 de un muestreo en el que el ancho no correspondía a una medición (era exactamente la décima parte del largo) y cuatro en los que el ancho era menor que 0.4 veces el largo, una proporción imposible para la forma de la concha. También se mantuvieron seis individuos aislados con ancho mayor que el largo. Estos 40 registros (30 + 4 + 6) se marcaron como dudosos para un análisis de sensibilidad. Como los datos fueron recolectados por terceros y no se contó con las hojas de campo, no fue posible verificarlos con los registros originales. La época de cada muestreo se asignó según el mes de colecta. Todo el proceso se realizó con un script de R.
El número de individuos por muestreo varió entre 1 y 30: 32 de los 36 muestreos tuvieron 30 individuos, dos tuvieron 28 y 25, y dos tuvieron menos de 10 (Mojón 69 en marzo de 2025, con un individuo, y Tambor en febrero de 2026, con cinco). Por la misma razón, no fue posible determinar por qué esos muestreos quedaron incompletos. Todos los individuos se incluyeron en las comparaciones entre áreas, años y épocas, en las que el individuo es la unidad de análisis. Los dos muestreos con menos de 10 individuos se excluyeron solo de los análisis basados en el largo medio de cada muestreo (la clasificación del PRONAMEC y la correlación con el clima), porque su promedio no se consideró representativo.
Obtención y acoplamiento de datos ambientales: El Instituto Meteorológico Nacional suministró registros horarios del periodo 2022–2025 de las estaciones automáticas San José-Pinilla (72-1493), Paquera (72-1594), Santa Rita (72-1895) y ASADA Artola (74-0816) (Instituto Meteorológico Nacional \[IMN\], 2026). Los registros se resumieron en valores diarios y un día se consideró válido si contaba con al menos 20 horas con datos. La estación Santa Rita se excluyó porque su serie presentó un desfase de fechas creciente respecto a las demás estaciones, detectado mediante correlación cruzada. Cada área se asignó a la estación utilizable más cercana en línea recta: los sitios de Cabo Blanco a Paquera (aproximadamente 30 km) y los de Ostional a San José-Pinilla (25–31 km). Los sitios de Camaronal no se incluyeron en el análisis climático por encontrarse a más de 50 km de cualquier estación utilizable. Para cada muestreo se calcularon la temperatura media del aire y la precipitación acumulada de los 90 días previos (día −90 a día −1), con al menos 80 % de días válidos; la precipitación acumulada se estimó como el promedio diario de los días válidos multiplicado por 90. Se eligió una ventana larga porque el crecimiento de la concha es acumulativo (Frank, 1969) y no existen estimaciones publicadas de la tasa de crecimiento de N. scabricosta (Chacón-Monge et al., 2019). Los registros suministrados abarcaron hasta el 31 de diciembre de 2025, por lo que los muestreos de 2026 no contaron con datos climáticos y no se incluyeron en este análisis. Estas variables representan condiciones climáticas regionales y no el microclima de los sitios.
Análisis descriptivo: Para el largo y el ancho de la concha se calcularon el tamaño de muestra (n), la media, la desviación estándar (DE), el error estándar, la mediana, los cuartiles, los valores mínimo y máximo, el coeficiente de variación (CV), la asimetría y la curtosis. Los cálculos se hicieron por área, año y época, y para las combinaciones de área y año, y de área, año y época. La curtosis se reportó sin centrar, de modo que el valor de referencia de una distribución normal es 3. Para describir la forma de las distribuciones se construyeron tablas de frecuencias del largo. El número de clases se definió con la regla de Sturges (Daniel & Cross, 2013), lo que dio 11 clases de 0.4 cm. Estos mismos intervalos se aplicaron a todas las áreas para que las distribuciones fueran comparables entre sí, y se representaron como histogramas de frecuencia relativa. La variación del largo entre áreas y años se representó con diagramas de caja. Estas medidas se examinaron antes de las pruebas de hipótesis para evaluar la forma de las distribuciones y orientar la elección entre pruebas paramétricas y no paramétricas.
Cada muestreo se clasificó según el ámbito de variación permisible del Indicador 2 del PRONAMEC (SINAC, 2016), que evalúa el largo medio de los individuos medidos en cada sitio y fecha mediante un sistema de semáforo con tres categorías: verde o bueno (> 2.5–4 cm), que según el protocolo indica condiciones aptas y ausencia de extracción; amarillo o regular (1.5–2.5 cm), que sugiere un grado medio de extracción; y rojo o malo (< 1 cm), que sugiere un exceso de extracción. El protocolo no define una categoría para los promedios entre 1 y 1.5 cm, por lo que estos muestreos se registraron como “sin categoría”. Los dos muestreos con menos de 10 individuos no se clasificaron.
Análisis inferencial: El plan de análisis y los criterios para elegir cada prueba de comparación (entre áreas, años y épocas), de asociación con el clima y con el semáforo del PRONAMEC se registraron antes de realizar las pruebas de hipótesis. El ANOVA de dos vías y el análisis con el muestreo como unidad se agregaron después, como análisis complementarios. Todas las pruebas fueron de dos colas y se evaluaron con un nivel de significancia α = 0.05. Para comparar tres o más grupos se ajustó primero un ANOVA de una vía y se evaluaron la normalidad de los residuos (Shapiro-Wilk) y la homogeneidad de varianzas (Levene). Si se cumplían ambos supuestos se usó el ANOVA con comparaciones de Tukey; si los residuos eran normales pero las varianzas diferían, el ANOVA de Welch con pruebas t por pares sin varianza común; y si no había normalidad, la prueba de Kruskal-Wallis con comparaciones por pares de Wilcoxon. En ambos casos las comparaciones por pares se ajustaron con Bonferroni. Para comparar dos grupos se evaluaron la normalidad de cada grupo y la igualdad de varianzas (prueba F); se aplicó la prueba t de Student o de Welch si los datos eran normales, y la de Wilcoxon-Mann-Whitney para muestras independientes si no lo eran. Esta última supone observaciones independientes y aleatorias, una variable continua y distribuciones de forma similar entre grupos; cuando las varianzas difirieron, la diferencia se interpretó como una tendencia de un grupo a presentar valores mayores y no como un desplazamiento de la mediana. El tamaño del efecto se expresó como η² en el ANOVA y como ε² en la prueba de Kruskal-Wallis.
Con este criterio se comparó el largo de la concha entre las cuatro áreas para todo el periodo 2022–2026 y entre años dentro de cada área. Las épocas seca y lluviosa se compararon solo en las combinaciones de área y año muestreadas en ambas épocas: la Reserva Cabo Blanco y el Área Marina Cabo Blanco en 2022, y Ostional en 2025. En las comparaciones entre años (cuatro pruebas) y entre épocas (tres pruebas), los valores de p se ajustaron con Bonferroni. La asociación entre el largo medio de cada muestreo y la temperatura media y la precipitación acumulada de los 90 días previos se evaluó mediante correlación, con el muestreo como unidad de análisis, porque todos los individuos de un muestreo comparten los mismos valores climáticos (n = 18 muestreos con estación asignada y al menos 10 individuos); se aplicó la correlación de Pearson si ambas variables seguían una distribución normal y la de Spearman en caso contrario. La asociación entre el área y la categoría del PRONAMEC se evaluó con la prueba exacta de Fisher (n = 34 muestreos), porque varias frecuencias esperadas eran menores de cinco.
Como análisis de sensibilidad, las comparaciones entre áreas, años y épocas se repitieron sin los 40 registros marcados como dudosos, y la correlación con el clima se repitió sin el sitio Punta Flor, que presentó tallas pequeñas en todos sus muestreos.Como análisis complementario no incluido en el plan metodológico inicial, se ajustó un ANOVA de dos vías con el área, el año y su interacción como factores. Se usaron solo los años 2025 y 2026, los únicos en que se muestrearon las cuatro áreas, de modo que todas las combinaciones de área y año tuvieran datos (n = 691 individuos). Como el diseño fue desbalanceado, se usaron sumas de cuadrados tipo II, y el tamaño del efecto se expresó como η² parcial. Se evaluaron la normalidad de los residuos (Shapiro-Wilk) y la homogeneidad de varianzas entre las ocho combinaciones de área y año (Levene); como las varianzas difirieron, la tabla se repitió con errores estándar robustos a la heterocedasticidad (HC3). La interacción se representó con las medias de cada combinación y su intervalo de confianza de 95 %. Además, para evaluar el efecto de la pseudorreplicación, la comparación entre áreas se repitió con la prueba de Kruskal-Wallis usando el largo medio de cada muestreo como unidad (n = 34 muestreos).
Finalmente, el ancho se analizó como complemento con las mismas pruebas, incluidos la correlación con el clima y el ANOVA de dos vías. Los análisis se realizaron en R versión 4.6.1 (R Core Team, 2026), con los paquetes readxl (Wickham & Bryan, 2026) para leer los datos, writexl (Ooms & Denney, 2026) para exportar los resultados, moments (Komsta & Novomestky, 2022) para la asimetría y la curtosis, y car (Fox & Weisberg, 2019) para la prueba de Levene y el ANOVA de dos vías.
Resultados
En total se analizaron 1 019 individuos de Nerita scabricosta en 36 muestreos realizados entre 2022 y 2026. La Reserva Cabo Blanco presentó el mayor largo medio y la mayor variabilidad, y fue el área con más individuos mayores de 3 cm (29.4 %, frente a 10.9 % en el Área Marina Cabo Blanco y menos de 1 % en Camaronal y Ostional). Camaronal fue el área con menor variabilidad. En el ancho, la Reserva Cabo Blanco también presentó el valor medio y la variabilidad más altos (Tabla 1).
Tabla 1. Estadística descriptiva del largo y ancho de la concha (cm) de Nerita scabricosta en cuatro áreas silvestres protegidas del Área de Conservación Tempisque, Costa Rica, 2022–2026.
| Variable | Área | n | Media ± DE | Mediana | Mín–Máx | CV (%) | Asimetría | Curtosis |
|---|---|---|---|---|---|---|---|---|
| Largo (cm) | RNA | 238 | 2.25 ± 0.97 | 2.31 | 0.60–4.20 | 43.1 | 0.04 | 1.66 |
| Largo (cm) | AMM | 275 | 2.08 ± 0.64 | 2.00 | 0.67–4.78 | 30.9 | 0.88 | 4.14 |
| Largo (cm) | Camaronal | 270 | 2.02 ± 0.49 | 2.10 | 0.77–3.12 | 24.2 | -0.59 | 2.74 |
| Largo (cm) | Ostional | 236 | 1.98 ± 0.55 | 2.00 | 0.83–3.13 | 27.7 | -0.04 | 1.92 |
| Largo (cm) | Total | 1019 | 2.08 ± 0.69 | 2.09 | 0.60–4.78 | 33.1 | 0.40 | 2.98 |
| Ancho (cm) | RNA | 208 | 1.96 ± 0.73 | 2.10 | 0.59–3.41 | 37.0 | -0.27 | 1.66 |
| Ancho (cm) | AMM | 275 | 1.59 ± 0.48 | 1.50 | 0.43–2.99 | 30.2 | 0.71 | 3.45 |
| Ancho (cm) | Camaronal | 269 | 1.64 ± 0.41 | 1.70 | 0.56–2.51 | 24.7 | -0.72 | 2.78 |
| Ancho (cm) | Ostional | 233 | 1.57 ± 0.46 | 1.60 | 0.63–2.60 | 29.5 | 0.02 | 2.04 |
| Ancho (cm) | Total | 985 | 1.68 ± 0.54 | 1.66 | 0.43–3.41 | 32.3 | 0.31 | 2.58 |
Nota: RNA: Reserva Natural Absoluta Cabo Blanco; AMM: Área Marina de Manejo Cabo Blanco; Camaronal: Refugio Nacional de Vida Silvestre Camaronal; Ostional: Refugio Nacional de Vida Silvestre Ostional. n: número de individuos; DE: desviación estándar; CV: coeficiente de variación. Curtosis no centrada (valor de referencia de la distribución normal = 3). El n del ancho es menor porque se excluyeron anchos no medidos o con valores imposibles.
Las distribuciones del largo difirieron en forma entre áreas (Fig. 1; Tabla 1). En la Reserva Cabo Blanco la distribución fue aplanada, con un grupo de individuos entre 1 y 1.4 cm y otro entre 2.6 y 3.8 cm. La mayoría de los individuos pequeños de esta área provino del sitio Punta Flor (54 de 68). El Área Marina Cabo Blanco y Camaronal mostraron un solo pico, entre 1.8 y 2.6 cm, con asimetría positiva en el AMM y negativa en Camaronal. Ostional también presentó una distribución aplanada, con individuos repartidos de forma similar entre 1 y 3 cm. En todas las áreas fueron escasos los individuos menores de 1 cm (3.5 % del total).
Fig. 1. Distribución de frecuencias relativas del largo de la concha de Nerita scabricosta en cuatro áreas silvestres protegidas del Área de Conservación Tempisque, Costa Rica, 2022–2026. Intervalos de clase de 0.4 cm (regla de Sturges). RNA: Reserva Natural Absoluta Cabo Blanco; AMM: Área Marina de Manejo Cabo Blanco.
El largo de la concha difirió entre áreas (Kruskal-Wallis, H = 10.34, gl = 3, p = 0.016), aunque el tamaño del efecto fue pequeño (ε² = 0.01; Fig. 2). Se aplicó esta prueba, con un nivel de significancia α = 0.05nza de 95 %, porque los residuos no siguieron una distribución normal (Shapiro-Wilk, W = 0.994, p < 0.001) y las varianzas difirieron entre áreas (Levene, F = 79.57, gl = 3 y 1 015, p < 0.001). La Reserva Cabo Blanco difirió de Camaronal (p = 0.023) y de Ostional (p = 0.020), pero no del Área Marina Cabo Blanco (p = 0.470); esta última, Camaronal y Ostional no difirieron entre sí (p = 1.000). Al excluir los registros dudosos, la diferencia entre áreas se mantuvo y fue mayor (H = 40.36, gl = 3, p < 0.001, ε² = 0.04), y la Reserva Cabo Blanco difirió de las otras tres áreas (p < 0.001).
Ver comandos y salida de R
# Comparación del largo entre áreas (2022–2026)
m_area <- aov(Largo ~ Area, data = bio)
shapiro.test(residuals(m_area)) # normalidad de residuos#>
#> Shapiro-Wilk normality test
#>
#> data: residuals(m_area)
#> W = 0.99361, p-value = 0.0002346leveneTest(Largo ~ Area, data = bio, center = mean) # homogeneidad de varianzas#> Levene's Test for Homogeneity of Variance (center = mean)
#> Df F value Pr(>F)
#> group 3 79.569 < 2.2e-16 ***
#> 1015
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1kw_area <- kruskal.test(Largo ~ Area, data = bio); kw_area#>
#> Kruskal-Wallis rank sum test
#>
#> data: Largo by Area
#> Kruskal-Wallis chi-squared = 10.341, df = 3, p-value = 0.01588unname(kw_area$statistic / (nrow(bio) - 1)) # épsilon cuadrado#> [1] 0.01015839pairwise.wilcox.test(bio$Largo, bio$Area, p.adjust.method = "bonferroni", exact = FALSE)#>
#> Pairwise comparisons using Wilcoxon rank sum test with continuity correction
#>
#> data: bio$Largo and bio$Area
#>
#> RNA Cabo Blanco AMM Cabo Blanco RNVS Camaronal
#> AMM Cabo Blanco 0.470 - -
#> RNVS Camaronal 0.023 1.000 -
#> RNVS Ostional 0.020 1.000 1.000
#>
#> P value adjustment method: bonferroni# Sensibilidad: sin los 40 registros marcados como dudosos
kruskal.test(Largo ~ Area, data = subset(bio, !Sensibilidad))#>
#> Kruskal-Wallis rank sum test
#>
#> data: Largo by Area
#> Kruskal-Wallis chi-squared = 40.356, df = 3, p-value = 8.958e-09
Fig. 2. Largo de la concha de Nerita scabricosta por área silvestre protegida, Área de Conservación Tempisque, Costa Rica, 2022–2026. La línea central indica la mediana, la caja el rango intercuartílico, los bigotes 1.5 veces ese rango y los círculos los valores atípicos. n: número de individuos. RNA: Reserva Natural Absoluta Cabo Blanco; AMM: Área Marina de Manejo Cabo Blanco.
El largo varió entre años dentro de cada área, sobre todo en la Reserva y el Área Marina Cabo Blanco (Fig. 3). En la Reserva Cabo Blanco, la mediana pasó de 1.12 cm en 2023 a 3.15 cm en 2025. con un nivel de significancia α = 0.05nza de 95 %, las diferencias entre años fueron significativas en la Reserva Cabo Blanco, el Área Marina Cabo Blanco y Camaronal, pero no en Ostional. En la Reserva Cabo Blanco y en Camaronal, los residuos no siguieron una distribución normal (Shapiro-Wilk, W = 0.986, p = 0.024, y W = 0.970, p < 0.001) y las varianzas difirieron entre años (Levene, F = 38.88, gl = 3 y 234, p < 0.001, y F = 5.13, gl = 2 y 267, p = 0.007), por lo que se aplicó la prueba de Kruskal-Wallis (Reserva Cabo Blanco: H = 113.47, gl = 3, p < 0.001, ε² = 0.48; Camaronal: H = 12.99, gl = 2, p = 0.006, ε² = 0.05).
En el Área Marina Cabo Blanco no se detectaron desviaciones de la normalidad de los residuos (W = 0.994, p = 0.303), pero las varianzas difirieron (F = 29.43, gl = 3 y 271, p < 0.001), por lo que se usó el ANOVA de Welch, que no requiere varianzas iguales (F = 22.41, gl = 3 y 120.6, p < 0.001, η² = 0.26). En Ostional, ninguno de los dos años siguió una distribución normal (2025: W = 0.955; 2026: W = 0.927; ambos p < 0.001), aunque las varianzas fueron similares (F = 1.32, gl = 145 y 89, p = 0.151); la prueba de Wilcoxon-Mann-Whitney no detectó diferencias (W = 7 512, p = 0.258). En la Reserva Cabo Blanco, todos los años difirieron entre sí excepto 2022 y 2026. En el Área Marina Cabo Blanco, el largo de 2025 fue mayor que el de los demás años (p < 0.001). En Camaronal, el largo de 2026 superó al de 2025 (p < 0.001).
Ver comandos y salida de R
# Largo entre años dentro de cada área (4 pruebas, Bonferroni x 4)
por_area <- split(bio, bio$Area)
for (a in names(por_area)) {
d <- por_area[[a]]; d$Anio <- factor(d$Anio)
cat("\n=====", a, "=====\n")
if (nlevels(d$Anio) > 2) {
print(shapiro.test(residuals(aov(Largo ~ Anio, data = d))))
print(leveneTest(Largo ~ Anio, data = d, center = mean))
} else {
print(tapply(d$Largo, d$Anio, function(x) shapiro.test(x)$p.value))
print(var.test(Largo ~ Anio, data = d))
}
}#>
#> ===== RNA Cabo Blanco =====
#>
#> Shapiro-Wilk normality test
#>
#> data: residuals(aov(Largo ~ Anio, data = d))
#> W = 0.98647, p-value = 0.0237
#>
#> Levene's Test for Homogeneity of Variance (center = mean)
#> Df F value Pr(>F)
#> group 3 38.879 < 2.2e-16 ***
#> 234
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> ===== AMM Cabo Blanco =====
#>
#> Shapiro-Wilk normality test
#>
#> data: residuals(aov(Largo ~ Anio, data = d))
#> W = 0.99368, p-value = 0.3025
#>
#> Levene's Test for Homogeneity of Variance (center = mean)
#> Df F value Pr(>F)
#> group 3 29.426 < 2.2e-16 ***
#> 271
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> ===== RNVS Camaronal =====
#>
#> Shapiro-Wilk normality test
#>
#> data: residuals(aov(Largo ~ Anio, data = d))
#> W = 0.96974, p-value = 1.776e-05
#>
#> Levene's Test for Homogeneity of Variance (center = mean)
#> Df F value Pr(>F)
#> group 2 5.127 0.006533 **
#> 267
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> ===== RNVS Ostional =====
#> 2025 2026
#> 1.207055e-04 8.239215e-05
#>
#> F test to compare two variances
#>
#> data: Largo by Anio
#> F = 1.3237, num df = 145, denom df = 89, p-value = 0.1511
#> alternative hypothesis: true ratio of variances is not equal to 1
#> 95 percent confidence interval:
#> 0.9017949 1.9100931
#> sample estimates:
#> ratio of variances
#> 1.323681p_anios <- c(
RNA = kruskal.test(Largo ~ factor(Anio), data = por_area[["RNA Cabo Blanco"]])$p.value,
AMM = oneway.test(Largo ~ factor(Anio), data = por_area[["AMM Cabo Blanco"]], var.equal = FALSE)$p.value,
Camaronal = kruskal.test(Largo ~ factor(Anio), data = por_area[["RNVS Camaronal"]])$p.value,
Ostional = wilcox.test(Largo ~ factor(Anio), data = por_area[["RNVS Ostional"]], exact = FALSE)$p.value)
oneway.test(Largo ~ factor(Anio), data = por_area[["AMM Cabo Blanco"]], var.equal = FALSE)#>
#> One-way analysis of means (not assuming equal variances)
#>
#> data: Largo and factor(Anio)
#> F = 22.413, num df = 3.00, denom df = 120.61, p-value = 1.331e-11p.adjust(p_anios, method = "bonferroni")#> RNA AMM Camaronal Ostional
#> 7.865311e-24 5.324126e-11 6.050473e-03 2.582155e-01
Fig. 3. Largo de la concha de Nerita scabricosta por área y año, Área de Conservación Tempisque, Costa Rica, 2022–2026. La línea central indica la mediana, la caja el rango intercuartílico, los bigotes 1.5 veces ese rango y los círculos los valores atípicos. Las posiciones vacías corresponden a años sin muestreo en esa área. RNA: Reserva Natural Absoluta Cabo Blanco; AMM: Área Marina de Manejo Cabo Blanco.
La diferencia entre épocas no fue consistente entre áreas. En los tres casos evaluados, al menos una de las épocas se desvió de la normalidad (Shapiro-Wilk: Reserva Cabo Blanco 2022, época seca, W = 0.915, p < 0.001; Área Marina Cabo Blanco 2022, época seca, W = 0.925, p = 0.001; Ostional 2025, época lluviosa, W = 0.926, p < 0.001), por lo que se aplicó la prueba de Wilcoxon-Mann-Whitney para muestras independientes, con un nivel de significancia α = 0.05nza de 95 %. Las varianzas fueron similares entre épocas en el Área Marina Cabo Blanco (F = 0.58, gl = 59 y 29, p = 0.075) y en Ostional (F = 0.62, gl = 55 y 89, p = 0.060), pero difirieron en la Reserva Cabo Blanco (F = 3.04, gl = 57 y 29, p = 0.002).
En la Reserva Cabo Blanco en 2022, el largo fue mayor en la época lluviosa que en la seca (W = 187, p < 0.001; diferencia de 1.08 cm, IC 95 %: 0.69–1.43); como las varianzas difirieron, este resultado indica que los individuos de la época lluviosa tendieron a ser más grandes, y no estrictamente un desplazamiento de la mediana. En el Área Marina Cabo Blanco en 2022 ocurrió lo contrario (W = 1 450, p < 0.001; diferencia de 0.48 cm, IC 95 %: 0.32–0.63). En Ostional en 2025 no hubo diferencia entre épocas (W = 2 450, p = 1.000). Los resultados de años y épocas no cambiaron al excluir los registros dudosos.
Ver comandos y salida de R
# Largo entre épocas (3 combinaciones área-año, Bonferroni x 3)
casos <- list("RNA 2022" = subset(bio, Area == "RNA Cabo Blanco" & Anio == 2022),
"AMM 2022" = subset(bio, Area == "AMM Cabo Blanco" & Anio == 2022),
"Ostional 2025" = subset(bio, Area == "RNVS Ostional" & Anio == 2025))
p_epocas <- sapply(names(casos), function(k) {
d <- casos[[k]]; d$Epoca <- factor(d$Epoca, levels = c("Seca", "Lluviosa"))
cat("\n=====", k, "=====\n")
print(tapply(d$Largo, d$Epoca, function(x) shapiro.test(x)$p.value))
print(var.test(Largo ~ Epoca, data = d))
w <- wilcox.test(Largo ~ Epoca, data = d, conf.int = TRUE, exact = FALSE); print(w)
w$p.value
})#>
#> ===== RNA 2022 =====
#> Seca Lluviosa
#> 0.0005995446 0.3350223431
#>
#> F test to compare two variances
#>
#> data: Largo by Epoca
#> F = 3.0428, num df = 57, denom df = 29, p-value = 0.001693
#> alternative hypothesis: true ratio of variances is not equal to 1
#> 95 percent confidence interval:
#> 1.547280 5.592256
#> sample estimates:
#> ratio of variances
#> 3.042758
#>
#>
#> Wilcoxon rank sum test with continuity correction
#>
#> data: Largo by Epoca
#> W = 187, p-value = 1.872e-09
#> alternative hypothesis: true location shift is not equal to 0
#> 95 percent confidence interval:
#> -1.4299859 -0.6900353
#> sample estimates:
#> difference in location
#> -1.080044
#>
#>
#> ===== AMM 2022 =====
#> Seca Lluviosa
#> 0.001247608 0.578786486
#>
#> F test to compare two variances
#>
#> data: Largo by Epoca
#> F = 0.5777, num df = 59, denom df = 29, p-value = 0.07486
#> alternative hypothesis: true ratio of variances is not equal to 1
#> 95 percent confidence interval:
#> 0.2945203 1.0565591
#> sample estimates:
#> ratio of variances
#> 0.5777037
#>
#>
#> Wilcoxon rank sum test with continuity correction
#>
#> data: Largo by Epoca
#> W = 1450, p-value = 2.547e-06
#> alternative hypothesis: true location shift is not equal to 0
#> 95 percent confidence interval:
#> 0.3200539 0.6300577
#> sample estimates:
#> difference in location
#> 0.479974
#>
#>
#> ===== Ostional 2025 =====
#> Seca Lluviosa
#> 0.4289879323 0.0000733006
#>
#> F test to compare two variances
#>
#> data: Largo by Epoca
#> F = 0.6236, num df = 55, denom df = 89, p-value = 0.06042
#> alternative hypothesis: true ratio of variances is not equal to 1
#> 95 percent confidence interval:
#> 0.3913667 1.0212781
#> sample estimates:
#> ratio of variances
#> 0.6235991
#>
#>
#> Wilcoxon rank sum test with continuity correction
#>
#> data: Largo by Epoca
#> W = 2450.5, p-value = 0.7812
#> alternative hypothesis: true location shift is not equal to 0
#> 95 percent confidence interval:
#> -0.2169777 0.1879537
#> sample estimates:
#> difference in location
#> -0.02382919p.adjust(p_epocas, method = "bonferroni")#> RNA 2022 AMM 2022 Ostional 2025
#> 5.616677e-09 7.639524e-06 1.000000e+00Con un nivel de significancia α = 0.05, el largo medio de los muestreos no se asoció con las condiciones climáticas de los 90 días previos (n = 18 muestreos; Fig. A1). Se aplicó la correlación de Spearman porque la temperatura (Shapiro-Wilk, W = 0.864, p = 0.014) y la precipitación (W = 0.848, p = 0.008) no siguieron una distribución normal, aunque el largo medio sí lo hizo (W = 0.956, p = 0.535). La correlación no fue significativa ni con la temperatura (ρ = 0.35, p = 0.150) ni con la precipitación (ρ = 0.34, p = 0.169), y la correlación de Pearson dio resultados equivalentes (temperatura: r = 0.36, IC 95 %: −0.13 a 0.71, p = 0.142; precipitación: r = 0.15, IC 95 %: −0.34 a 0.58, p = 0.555). El resultado se mantuvo al excluir el sitio Punta Flor (temperatura: ρ = 0.26, p = 0.324; precipitación: ρ = 0.20, p = 0.469). En estos muestreos, la temperatura media varió entre 25.9 y 28.3 °C, y la precipitación acumulada, entre 4 y 1 090 mm.
Ver comandos y salida de R
# Correlación entre el largo medio por muestreo y el clima de 90 días (n = 18)
p4 <- as.data.frame(read_excel("inferencial/Inferencial_Nerita.xlsx", sheet = "Datos_P4_muestreos"))
sapply(p4[c("Largo_medio", "T_90d", "P_90d")], function(x) shapiro.test(x)$p.value)#> Largo_medio T_90d P_90d
#> 0.535446017 0.014300224 0.007758464cor.test(p4$Largo_medio, p4$T_90d, method = "spearman", exact = FALSE)#>
#> Spearman's rank correlation rho
#>
#> data: p4$Largo_medio and p4$T_90d
#> S = 626.29, p-value = 0.1499
#> alternative hypothesis: true rho is not equal to 0
#> sample estimates:
#> rho
#> 0.3536719cor.test(p4$Largo_medio, p4$P_90d, method = "spearman", exact = FALSE)#>
#> Spearman's rank correlation rho
#>
#> data: p4$Largo_medio and p4$P_90d
#> S = 640.45, p-value = 0.1687
#> alternative hypothesis: true rho is not equal to 0
#> sample estimates:
#> rho
#> 0.3390639cor.test(p4$Largo_medio, p4$T_90d, method = "pearson")#>
#> Pearson's product-moment correlation
#>
#> data: p4$Largo_medio and p4$T_90d
#> t = 1.5466, df = 16, p-value = 0.1415
#> alternative hypothesis: true correlation is not equal to 0
#> 95 percent confidence interval:
#> -0.1277385 0.7082585
#> sample estimates:
#> cor
#> 0.3606391cor.test(p4$Largo_medio, p4$P_90d, method = "pearson")#>
#> Pearson's product-moment correlation
#>
#> data: p4$Largo_medio and p4$P_90d
#> t = 0.60288, df = 16, p-value = 0.555
#> alternative hypothesis: true correlation is not equal to 0
#> 95 percent confidence interval:
#> -0.3416027 0.5758392
#> sample estimates:
#> cor
#> 0.1490363# Sensibilidad: sin Punta Flor
sin_pf <- subset(p4, Sitio != "Punta Flor")
cor.test(sin_pf$Largo_medio, sin_pf$T_90d, method = "spearman", exact = FALSE)#>
#> Spearman's rank correlation rho
#>
#> data: sin_pf$Largo_medio and sin_pf$T_90d
#> S = 500.74, p-value = 0.3239
#> alternative hypothesis: true rho is not equal to 0
#> sample estimates:
#> rho
#> 0.2636233cor.test(sin_pf$Largo_medio, sin_pf$P_90d, method = "spearman", exact = FALSE)#>
#> Spearman's rank correlation rho
#>
#> data: sin_pf$Largo_medio and sin_pf$P_90d
#> S = 547.22, p-value = 0.4686
#> alternative hypothesis: true rho is not equal to 0
#> sample estimates:
#> rho
#> 0.1952697La categoría del PRONAMEC se asoció con el área (prueba exacta de Fisher, p = 0.015; n = 34 muestreos). Se usó esta prueba en lugar de χ² porque varias frecuencias esperadas fueron menores de cinco. Tres de los cuatro muestreos en la categoría verde (bueno) correspondieron a la Reserva Cabo Blanco, y los muestreos sin categoría se presentaron en la Reserva Cabo Blanco (Punta Flor) y en Ostional (Mojón 69). En conjunto, la mayoría de los muestreos se ubicó en la categoría amarilla (regular) y ninguno en la roja (malo) (Tabla 2).
Ver comandos y salida de R
# Clasificación PRONAMEC por muestreo (n >= 10) y prueba exacta de Fisher
bio$Muestreo <- paste(bio$Sitio, as.Date(bio$Fecha))
mue <- do.call(data.frame, aggregate(Largo ~ Muestreo + Area, data = bio,
FUN = function(x) c(n = length(x), media = mean(x))))
names(mue)[3:4] <- c("n", "Largo_medio")
mue <- subset(mue, n >= 10)
mue$Categoria <- cut(mue$Largo_medio, c(-Inf, 1, 1.5, 2.5, 4, Inf), right = FALSE,
labels = c("Malo", "Sin categoría", "Regular", "Bueno", "Sobre el ámbito"))
tab_p5 <- table(mue$Area, droplevels(mue$Categoria)); tab_p5#>
#> Sin categoría Regular Bueno
#> RNA Cabo Blanco 2 3 3
#> AMM Cabo Blanco 0 8 1
#> RNVS Camaronal 0 9 0
#> RNVS Ostional 2 6 0fisher.test(tab_p5)#>
#> Fisher's Exact Test for Count Data
#>
#> data: tab_p5
#> p-value = 0.01494
#> alternative hypothesis: two.sidedTabla 2. Número de muestreos de Nerita scabricosta por área silvestre protegida según el ámbito de variación permisible del largo de la concha del PRONAMEC, Área de Conservación Tempisque, Costa Rica, 2022–2026.
| Área | Bueno | Regular | Sin categoría | Malo | Total |
|---|---|---|---|---|---|
| RNA | 3 | 3 | 2 | 0 | 8 |
| AMM | 1 | 8 | 0 | 0 | 9 |
| Camaronal | 0 | 9 | 0 | 0 | 9 |
| Ostional | 0 | 6 | 2 | 0 | 8 |
| Total | 4 | 26 | 4 | 0 | 34 |
Nota. Cada muestreo corresponde a un sitio en una fecha y se clasificó según el largo medio de sus individuos (SINAC, 2016, Indicador 2): bueno, > 2.5–4 cm; regular, 1.5–2.5 cm; malo, < 1 cm. Sin categoría: largo medio entre 1 y 1.5 cm, intervalo no definido en el protocolo. Se excluyeron dos muestreos con menos de 10 individuos. RNA: Reserva Natural Absoluta Cabo Blanco; AMM: Área Marina de Manejo Cabo Blanco; Camaronal: Refugio Nacional de Vida Silvestre Camaronal; Ostional: Refugio Nacional de Vida Silvestre Ostional.
En el ANOVA de dos vías para 2025 y 2026 (n = 691), la interacción entre el área y el año fue significativa (F = 42.46, gl = 3 y 683, p < 0.001, η² parcial = 0.16), lo que indica que el cambio en el largo entre años no fue el mismo en todas las áreas (Fig. 4). De 2025 a 2026, el largo medio disminuyó en la Reserva Cabo Blanco (de 3.27 a 2.13 cm) y en el Área Marina Cabo Blanco (de 2.54 a 1.85 cm), aumentó en Camaronal (de 1.87 a 2.18 cm) y varió poco en Ostional (de 2.03 a 1.89 cm). El área (F = 47.27, gl = 3 y 683, p < 0.001, η² parcial = 0.17) y el año (F = 49.44, gl = 1 y 683, p < 0.001, η² parcial = 0.07) también fueron significativos, pero, por la interacción, sus efectos no deben interpretarse por separado. Los residuos se desviaron levemente de la normalidad (Shapiro-Wilk, W = 0.994, p = 0.010) y las varianzas difirieron entre combinaciones (Levene, F = 14.45, gl = 7 y 683, p < 0.001); con errores robustos a la heterocedasticidad, los tres términos se mantuvieron significativos (interacción: F = 41.11, gl = 3 y 683, p < 0.001). Por último, con el largo medio de cada muestreo como unidad, la diferencia entre áreas no fue significativa (Kruskal-Wallis, H = 1.26, gl = 3, p = 0.740; n = 34 muestreos).
Ver comandos y salida de R
# ANOVA de dos vías Área x Año (2025-2026, análisis no planificado)
d2 <- subset(bio, Anio %in% c(2025, 2026)); d2$Anio <- factor(d2$Anio)
table(d2$Area, d2$Anio)#>
#> 2025 2026
#> RNA Cabo Blanco 60 60
#> AMM Cabo Blanco 90 65
#> RNVS Camaronal 90 90
#> RNVS Ostional 146 90modelo2 <- aov(Largo ~ Area * Anio, data = d2)
Anova(modelo2, type = 2) # sumas de cuadrados tipo II#> Anova Table (Type II tests)
#>
#> Response: Largo
#> Sum Sq Df F value Pr(>F)
#> Area 49.836 3 47.274 < 2.2e-16 ***
#> Anio 17.373 1 49.441 4.971e-12 ***
#> Area:Anio 44.759 3 42.458 < 2.2e-16 ***
#> Residuals 240.006 683
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1shapiro.test(residuals(modelo2))#>
#> Shapiro-Wilk normality test
#>
#> data: residuals(modelo2)
#> W = 0.99423, p-value = 0.009925leveneTest(Largo ~ Area * Anio, data = d2)#> Levene's Test for Homogeneity of Variance (center = median)
#> Df F value Pr(>F)
#> group 7 14.454 < 2.2e-16 ***
#> 683
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(lm(Largo ~ Area * Anio, data = d2), type = 2, white.adjust = "hc3") # errores robustos#> Analysis of Deviance Table (Type II tests)
#>
#> Response: Largo
#> Df F Pr(>F)
#> Area 3 101.751 < 2.2e-16 ***
#> Anio 1 19.060 1.463e-05 ***
#> Area:Anio 3 41.106 < 2.2e-16 ***
#> Residuals 683
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1interaction.plot(d2$Anio, d2$Area, d2$Largo, xlab = "Año", ylab = "Largo medio (cm)",
trace.label = "Área")
# Kruskal-Wallis con el muestreo como unidad (pseudorreplicación)
kruskal.test(Largo_medio ~ Area, data = mue)#>
#> Kruskal-Wallis rank sum test
#>
#> data: Largo_medio by Area
#> Kruskal-Wallis chi-squared = 1.2556, df = 3, p-value = 0.7397
Fig. 4. Largo medio de la concha de Nerita scabricosta (± IC 95 %) por área silvestre protegida y año, Área de Conservación Tempisque, Costa Rica, 2025–2026 (n = 691 individuos). Las líneas unen las medias de cada área para mostrar la interacción. RNA: Reserva Natural Absoluta Cabo Blanco; AMM: Área Marina de Manejo Cabo Blanco.
El ancho mostró un patrón similar al del largo. Difirió entre áreas (Kruskal-Wallis, H = 52.98, gl = 3, p < 0.001, ε² = 0.05; n = 985); a diferencia del largo, la Reserva Cabo Blanco difirió de las otras tres áreas (p < 0.001), y el Área Marina Cabo Blanco difirió de Camaronal (p = 0.007). Entre años, el ancho varió en la Reserva Cabo Blanco (Kruskal-Wallis, H = 76.61, gl = 2, p < 0.001, ε² = 0.37), el Área Marina Cabo Blanco (ANOVA de Welch, F = 27.71, gl = 3 y 116, p < 0.001, η² = 0.30) y Camaronal (Kruskal-Wallis, H = 10.89, gl = 2, p = 0.017, ε² = 0.04), pero no en Ostional (Wilcoxon-Mann-Whitney, W = 6 622, p = 1.000). Entre épocas, fue mayor en la época lluviosa en la Reserva Cabo Blanco en 2022 (W = 194, p < 0.001), mayor en la seca en el Área Marina Cabo Blanco en 2022 (t de Student, t = 5.27, gl = 88, p < 0.001) y no difirió en Ostional en 2025 (W = 2 450, p = 1.000); los valores de p entre años y épocas están ajustados con Bonferroni. El ancho medio de los muestreos no se correlacionó con la temperatura (Spearman, ρ = 0.26, p = 0.308) ni con la precipitación (ρ = 0.27, p = 0.300; n = 17 muestreos, porque uno no tuvo valores de ancho válidos).
En el ANOVA de dos vías para 2025 y 2026 (n = 687), la interacción entre el área y el año también fue significativa (F = 46.72, gl = 3 y 679, p < 0.001, η² parcial = 0.17), y se mantuvo con errores robustos a la heterocedasticidad (F = 44.57, p < 0.001); los residuos no se desviaron de la normalidad (W = 0.996, p = 0.051), pero las varianzas difirieron (Levene, F = 17.27, gl = 7 y 679, p < 0.001). Con el ancho medio de cada muestreo como unidad, tampoco hubo diferencia entre áreas (Kruskal-Wallis, H = 3.80, gl = 3, p = 0.284; n = 33 muestreos). En resumen, el ancho llevó a las mismas conclusiones que el largo en casi todas las comparaciones. La diferencia estuvo en las comparaciones entre áreas: la Reserva Cabo Blanco se distinguió de las otras tres áreas, incluida el Área Marina Cabo Blanco, y el Área Marina Cabo Blanco difirió de Camaronal.
Ver comandos y salida de R
# Resultados del ancho: generados con las mismas funciones del script 03 y 04
res <- as.data.frame(read_excel("inferencial/Inferencial_Nerita.xlsx", sheet = "Resumen_pruebas"))
kable(subset(res, Variable == "Ancho"), caption = "Pruebas de hipótesis para el ancho (script 03).")| Pregunta | Grupo | Variable | Factor | n | Niveles | Prueba | Estadistico | gl | p | Efecto | p_Bonferroni | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 3 | P1 Área (ancho, complemento) | Todos | Ancho | Area | 985 | RNA, AMM, Camaronal, Ostional | Kruskal-Wallis | H = 52.98 | 3 | 0.000000 | ε² = 0.054 | NA |
| 12 | P2 Año (ancho, complemento) | RNA | Ancho | Anio | 208 | 2022, 2025, 2026 | Kruskal-Wallis | H = 76.61 | 2 | 0.000000 | ε² = 0.370 | 0.00e+00 |
| 13 | P2 Año (ancho, complemento) | AMM | Ancho | Anio | 275 | 2022, 2023, 2025, 2026 | ANOVA de Welch | F = 27.71 | 3, 116 | 0.000000 | η² = 0.296 | 0.00e+00 |
| 14 | P2 Año (ancho, complemento) | Camaronal | Ancho | Anio | 269 | 2022, 2025, 2026 | Kruskal-Wallis | H = 10.89 | 2 | 0.004310 | ε² = 0.041 | 1.72e-02 |
| 15 | P2 Año (ancho, complemento) | Ostional | Ancho | Anio | 233 | 2025, 2026 | Wilcoxon (Mann-Whitney) | W = 6622 | NA | 0.669000 | Dif. de ubicación (Hodges-Lehmann) = 0.04 cm (IC 95 %: -0.12 a 0.18) | 1.00e+00 |
| 22 | P3 Época (ancho, complemento) | RNA 2022 | Ancho | Epoca | 88 | Seca, Lluviosa | Wilcoxon (Mann-Whitney) | W = 194 | NA | 0.000000 | Dif. de ubicación (Hodges-Lehmann) = -0.89 cm (IC 95 %: -1.10 a -0.61) | 0.00e+00 |
| 23 | P3 Época (ancho, complemento) | AMM 2022 | Ancho | Epoca | 90 | Seca, Lluviosa | t de Student | t = 5.27 | 88 | 0.000001 | Dif. de medias = 0.32 cm (IC 95 %: 0.20 a 0.44) | 2.90e-06 |
| 24 | P3 Época (ancho, complemento) | Ostional 2025 | Ancho | Epoca | 144 | Seca, Lluviosa | Wilcoxon (Mann-Whitney) | W = 2450 | NA | 0.938000 | Dif. de ubicación (Hodges-Lehmann) = 0.00 cm (IC 95 %: -0.14 a 0.17) | 1.00e+00 |
np <- "inferencial/No_planificados_Nerita.xlsx"
kable(read_excel(np, sheet = "A3_ancho_clima"), caption = "Correlación entre el ancho medio y el clima (script 04).")| Variable | n | rho | p_spearman | r | IC95 | p_pearson |
|---|---|---|---|---|---|---|
| T_90d | 17 | 0.26 | 0.308 | 0.30 | -0.21 a 0.68 | 0.235 |
| P_90d | 17 | 0.27 | 0.300 | 0.09 | -0.41 a 0.55 | 0.727 |
kable(read_excel(np, sheet = "A3_ancho_ANOVA"), caption = "ANOVA de dos vías para el ancho (script 04).")| Termino | gl | F | p | eta2_parcial | F_HC3 | p_HC3 |
|---|---|---|---|---|---|---|
| Area | 3 | 47.66 | 0 | 0.174 | 144.08 | 0.0e+00 |
| Anio | 1 | 44.29 | 0 | 0.061 | 21.56 | 4.1e-06 |
| Area:Anio | 3 | 46.72 | 0 | 0.171 | 44.57 | 0.0e+00 |
| Residuals | 679 | NA | NA | NA | NA | NA |
Discusión y conclusiones
La protección absoluta de la Reserva Cabo Blanco se asocia con burgados más grandes, pero esa señal es débil frente a la variación natural entre sitios. La reserva, el área con mayor protección, presentó el mayor largo medio y la mayor proporción de individuos mayores de 3 cm, un patrón esperable donde no se permite la extracción: en zonas sin extracción de otras costas rocosas se han registrado individuos más grandes y abundantes de las especies explotadas, aunque ese efecto dependió del tipo de sustrato y de la región (Baliwe et al., 2022). Nerita scabricosta se extrae en la pesca de pequeña escala en Costa Rica (Azofeifa-Solano et al., 2018) y como alimento y carnada sin regulación en México (Gutiérrez-Rubio & Arzola-González, 2024), por lo que la extracción es una explicación plausible. Sin embargo, la evidencia de este estudio es débil: la reserva difirió de Camaronal y de Ostional, pero no del Área Marina Cabo Blanco, y el área explicó una proporción muy pequeña de la variación. En el ancho, la reserva sí difirió de las otras tres áreas, aunque el tamaño del efecto también fue pequeño. Además, la categoría de manejo es solo una aproximación a la presión de extracción, que no se midió (Nielsen et al., 2025), y las diferencias de talla atribuidas a la extracción pueden confundirse con diferencias ambientales entre sitios, como la exposición al oleaje (Doyle et al., 2022). Por otra parte, el área sin extracción estuvo representada únicamente por la Reserva Cabo Blanco, por lo que esta comparación no permite separar el efecto de la protección de las características particulares de esa área.
Por lo tanto, las hipótesis planteadas recibieron un respaldo parcial: si bien la Reserva Cabo Blanco presentó individuos ligeramente mayores, este efecto fue débil y quedó superado por la marcada variación entre sitios y fechas de muestreo; además, la talla de N. scabricosta no siguió un patrón estacional consistente ni se asoció con la precipitación o la temperatura regional. Este patrón coincide con lo observado en el Golfo de California, donde la talla de la especie difirió entre islas cercanas según las condiciones de cada sitio (Gutiérrez-Rubio & Arzola-González, 2024), y en un plano distinto, el genómico, con lo reportado en Lottia gigantea, donde la variación climática dejó una huella selectiva mayor que la cosecha, aunque las áreas marinas protegidas sí mostraron una tendencia creciente en la abundancia de individuos grandes (Nielsen et al., 2025). En conjunto, estos resultados indican que, en el periodo evaluado, la identidad del sitio y el momento del muestreo influyen más en la estructura de tallas que la categoría de manejo del área.
La variación entre sitios de una misma área ayuda a explicar estos resultados. La distribución aplanada de la Reserva Cabo Blanco refleja dos grupos de individuos de diferente tamaño: la mayoría de los individuos pequeños provino de Punta Flor (54 de 68), mientras que San Miguel y Playa Cabo Blanco concentraron los individuos grandes. Punta Flor presentó tallas pequeñas tanto en 2022 como en 2023, lo que sugiere una característica estable del sitio y no un evento puntual. En el Golfo de California, el ancho medio de N. scabricosta varió de 12.2 a 26.7 mm entre islas cercanas, y los autores lo relacionaron con la pendiente, la disponibilidad de grietas y macroalgas, el oleaje y la temperatura (Gutiérrez-Rubio & Arzola-González, 2024). La exposición al oleaje también se ha asociado con la talla de otros gasterópodos (Rodríguez-Grimón et al., 2024) y con la distribución de N. scabricosta (Azofeifa-Solano et al., 2018). Como en este estudio no se midieron estas condiciones, no es posible separar su efecto del efecto del área.
De la misma manera, la composición de sitios explica parte de la variación entre años. En la Reserva Cabo Blanco, el muestreo de 2023 incluyó únicamente Punta Flor, mientras que en 2025 este sitio no se muestreó, de modo que el contraste entre esos años mezcla el efecto del año con el del sitio. La variación restante podría relacionarse con que N. scabricosta se reproduce solo durante la época lluviosa, sin sincronía con la marea ni con la fase lunar (Collin et al., 2017), y con que el tiempo de eclosión de sus cápsulas varía entre pozas de marea (Collin et al., 2016), lo que podría generar diferencias en el momento del reclutamiento entre sitios y años. En N. japonica, el reclutamiento modificó la densidad poblacional sin cambiar la talla media de los adultos (Takada, 2008), mientras que en Lottia gigantea una anomalía cálida aumentó el reclutamiento y la proporción de individuos pequeños (Nielsen et al., 2025). Estos procesos podrían contribuir a la variación interanual observada, aunque no pueden confirmarse con los datos disponibles.
La ausencia de un patrón estacional consistente apunta en la misma dirección. En la Reserva y el Área Marina Cabo Blanco en 2022, las dos épocas no incluyeron los mismos sitios: en ambos casos, la época seca incluyó un sitio que no se muestreó en la lluviosa. En cambio, en Ostional en 2025, el único caso en que se muestrearon los mismos tres sitios en ambas épocas, no se detectó diferencia. Esto sugiere que las diferencias entre épocas reflejan principalmente qué sitios se muestrearon y no un efecto de la estacionalidad sobre la talla.
El ANOVA de dos vías, aunque no estaba incluido en el plan inicial, agrega un matiz a esta interpretación. La interacción entre el área y el año indica que el largo cambió en direcciones opuestas entre 2025 y 2026 según el área. Sin embargo, en esos años el año coincidió con la época: en la Reserva y el Área Marina Cabo Blanco, los muestreos de 2025 se hicieron en la época lluviosa y los de 2026 en la seca, mientras que en Camaronal ocurrió lo contrario, y en las tres áreas el largo fue mayor en el año muestreado en la época lluviosa. Por ello, la interacción podría reflejar en parte un efecto de la época que no se observó de forma consistente en 2022, y este diseño no permite separar ambos factores.
La falta de asociación con el clima puede deberse tanto a las limitaciones de los datos como a la biología de la especie. Los 18 muestreos compartieron solo 11 valores distintos de temperatura, porque los sitios muestreados en una misma fecha recibieron el mismo valor, lo que reduce la capacidad de las pruebas para detectar una asociación. Además, la temperatura del aire registrada a 25–31 km no representa la que experimenta el caracol. En costas rocosas tropicales, la roca alcanza más de 50 °C y la temperatura corporal de especies de Nerita puede acercarse a los 50 °C (Marshall et al., 2023), muy por encima del intervalo de 25.9 a 28.3 °C registrado por las estaciones. Las especies de Nerita también seleccionan microhábitats entre 7 y 14 °C más frescos que su límite térmico (Hui et al., 2022), y su tolerancia térmica responde a la temperatura local de cada sitio (Hui et al., 2025). Por ello, el clima regional difícilmente captura las condiciones térmicas que influyen en el crecimiento.
Según el sistema de semáforo del PRONAMEC, la mayoría de los muestreos se ubicó en la categoría amarilla (regular) y ninguno en la roja (malo). De acuerdo con el protocolo, la categoría regular sugiere un grado medio de extracción y la roja, un exceso de extracción (SINAC, 2016), por lo que en ninguna de las áreas se detectaron señales de sobreexplotación severa. Sin embargo, la categoría regular también fue frecuente en la Reserva Cabo Blanco (tres de ocho muestreos), donde la extracción está prohibida, y abarcó todos los muestreos de Camaronal, lo que sugiere que el intervalo de 1.5 a 2.5 cm refleja en buena medida la talla habitual de la especie en estos sitios y no necesariamente extracción. Esto es coherente con que los umbrales del protocolo se basan en estudios de hace casi cuatro décadas (Ortega, 1986, 1987, como se citó en SINAC, 2016) y no se han calibrado para cada sitio. La asociación entre la categoría y el área (prueba exacta de Fisher, p = 0.015) se debió principalmente a que la Reserva Cabo Blanco concentró tres de los cuatro muestreos en la categoría verde (bueno), lo cual es consistente con la interpretación del protocolo, según la cual esta categoría indica ausencia de extracción (SINAC, 2016). Por otra parte, los muestreos sin categoría, dos de Punta Flor (Reserva Cabo Blanco) y dos de Mojón 69 (Ostional), se deben a que sus promedios se ubicaron entre 1.0 y 1.5 cm, intervalo para el cual el protocolo no define una categoría; estas observaciones no necesariamente representan un deterioro de las poblaciones, sino una limitación del protocolo para clasificar muestreos con predominio de individuos pequeños. Además, al evaluar solo el promedio, el indicador no distingue una población homogénea de individuos medianos de otra con una mezcla de tallas, como la de la Reserva Cabo Blanco. Definir este intervalo, calibrar los umbrales con datos locales y complementar el promedio con la distribución de frecuencias de tallas mejoraría la utilidad del indicador para el manejo costero.
Los individuos menores de 1 cm representaron solo el 3.5 % del total. Esta escasez probablemente refleja un sesgo de muestreo más que una ausencia real de juveniles: N. scabricosta se refugia en grietas y hendiduras para evitar la desecación y la depredación, y forma agregaciones con fines de protección (Gutiérrez-Rubio & Arzola-González, 2024). Si durante el muestreo no se revisan las grietas ni el interior de las agregaciones, los juveniles refugiados quedan subrepresentados en las mediciones.
Estos resultados deben interpretarse considerando varias limitaciones. Las pruebas entre áreas, años y épocas trataron a cada individuo como una observación independiente, aunque los individuos de un mismo muestreo comparten las condiciones del sitio; esta pseudorreplicación tiende a subestimar los valores de p (Hurlbert, 1984). En un análisis complementario, con el largo medio de cada muestreo como unidad (n = 34 muestreos), la diferencia entre áreas no fue significativa (p = 0.740). Además, el resultado entre áreas dependió de los sitios incluidos: al excluir los registros dudosos, que correspondían principalmente a Punta Flor en 2023, un sitio con ejemplares pequeños, la diferencia de la Reserva Cabo Blanco aumentó. Por último, los datos fueron recolectados por terceros y no se contó con las hojas de campo, lo que impidió verificar los registros dudosos, y el diseño fue desbalanceado, porque no todas las áreas se muestrearon en todos los años ni en ambas épocas; en particular, la época lluviosa de 2026 estaba en curso mientras se realizó este estudio, por lo que de ese periodo solo se contó con los muestreos de Camaronal de mayo de 2026.
Para futuras investigaciones sobre N. scabricosta en el litoral rocoso, se recomienda mantener los mismos sitios en todas las campañas y en ambas épocas, ya que la elección del diseño de monitoreo condiciona los parámetros que se pueden estimar (Wright et al., 2022). También se recomienda registrar cuantitativamente la presión de extracción, la afluencia turística y la exposición al oleaje; muestrear activamente dentro de grietas y agregaciones para estimar mejor la abundancia de juveniles menores de 1 cm; instalar, en la medida de lo posible, sensores de temperatura sobre las rocas intermareales para capturar el microclima que experimenta la especie, y analizar los datos con modelos lineales mixtos que incluyan el sitio de muestreo como efecto aleatorio.
En conclusión, la estructura de tallas de N. scabricosta en las áreas protegidas del Pacífico norte de Costa Rica depende principalmente de las características de cada sitio y del momento del muestreo, más que de la época o de las condiciones climáticas regionales. La Reserva Cabo Blanco presentó tallas mayores, pero esa diferencia es pequeña y no se sostiene al usar el muestreo como unidad, por lo que no es posible atribuirla con seguridad a la prohibición de extracción. Estos resultados destacan la necesidad de diseños de monitoreo balanceados y de calibrar los intervalos del protocolo PRONAMEC con datos locales.
Bibliografía
Alvarado, J. J., Cortés, J., Esquivel, M. F., & Salas, E. (2012). Costa Rica’s marine protected areas: Status and perspectives. Revista de Biología Tropical, 60(1), 129–142. https://doi.org/10.15517/rbt.v60i1.2753
Arias, A., Cinner, J. E., Jones, R. E., & Pressey, R. L. (2015). Levels and drivers of fishers’ compliance with marine protected areas. Ecology and Society, 20(4), 19. https://doi.org/10.5751/ES-07999-200419
Azofeifa-Solano, J. C., Salas-Moya, C., Cubero-Campos, Y., & Sibaja-Cordero, J. A. (2018). Influencia de la vigilancia en una zona rocosa dentro de un área marina protegida, Pacífico Central de Costa Rica. Revista de Biología Tropical, 66(3), 984–995. https://doi.org/10.15517/rbt.v66i3.30487
Baliwe, N. G., Pfaff, M. C., & Branch, G. M. (2022). Assessing the effects of no-take zones in a marine protected area spanning two ecoregions and rock substrate types. Frontiers in Marine Science, 9, 893260. https://doi.org/10.3389/fmars.2022.893260
Chacón-Monge, J. L., Cubero-Campos, Y., & Sibaja-Cordero, J. A. (2019). Efecto de un área marina protegida sobre la talla de la concha de dos moluscos del intermareal en el Pacífico Central de Costa Rica. Intropica, 14(2), 138–147. https://doi.org/10.21676/23897864.3365
Collin, R., Kerr, K., Contolini, G., & Ochoa, I. (2017). Reproductive cycles in tropical intertidal gastropods are timed around tidal amplitude cycles. Ecology and Evolution, 7(15), 5977–5991. https://doi.org/10.1002/ece3.3166
Collin, R., Roof, K. E., & Spangler, A. (2016). Hatching plasticity in the tropical gastropod Nerita scabricosta. Invertebrate Biology, 135(2), 87–96. https://doi.org/10.1111/ivb.12119
Daniel, W. W., & Cross, C. L. (2013). Biostatistics: A foundation for analysis in the health sciences (10.ª ed.). John Wiley & Sons.
Decreto Ejecutivo N.º 23640-MIRENEM. (1994, 28 de septiembre). Declara Refugio Nacional de Vida Silvestre Camaronal. La Gaceta N.º 184. San José, Costa Rica.
Decreto Ejecutivo N.º 40442-MINAE. (2017, 9 de agosto). Creación del Área Marina de Manejo Cabo Blanco. La Gaceta N.º 151. San José, Costa Rica.
Decreto Ley N.º 3763. (1966, 19 de octubre). Creación de la Reserva Natural Absoluta Cabo Blanco. Asamblea Legislativa de la República de Costa Rica. San José, Costa Rica.
Doyle, D., Frias, J., Gammell, M. P., Lynch, M., & Nash, R. (2022). Assessing the morphological impacts of long-term harvesting in intertidal gastropods using historical data and morphometric tools. Journal of Molluscan Studies, 88(3), eyac019. https://doi.org/10.1093/mollus/eyac019
Fallas-Flores, S., Abarca-Guerrero, L., & Ulate-Naranjo, K. (2024). Propuesta de mejoras para la implementación del monitoreo de residuos sólidos mediante el PRONAMEC en las Áreas Marinas Protegidas de Costa Rica. Tecnología en Marcha, 37(4), 31–47. https://doi.org/10.18845/tm.v37i4.6882
Faria, J., Vale, M., Ribeiro, P., Hawkins, S. J., & Martins, G. M. (2024). The ecological value of fully enforced, no-entry, marine protected areas: A case study of harvested limpets. Aquatic Conservation: Marine and Freshwater Ecosystems, 34(2). https://doi.org/10.1002/aqc.4086
Fox, J., & Weisberg, S. (2019). An R companion to applied regression (3.ª ed.). Sage. https://www.john-fox.ca/Companion/
Frank, P. W. (1969). Growth rates and longevity of some gastropod mollusks on the coral reef at Heron Island. Oecologia, 2, 232–250. https://doi.org/10.1007/BF00379161
Gutiérrez-Rubio, Y., & Arzola-González, J. F. (2024). Densidad y estructura de tallas de Nerita scabricosta (Mollusca: Neritidae) en la zona intermareal de las islas del noroeste del Golfo de California, México. Revista Científica FCV-LUZ, 34, e34377. https://doi.org/10.52973/rcfcv-e34377
Hui, T. Y., Crickenberger, S., Chan, S. H. M., Todd, P. A., & Williams, G. A. (2025). How does local temperature shape thermal tolerance? A test using congeneric snails on tropical rocky shores. The Biological Bulletin, 249(1), 13–29. https://doi.org/10.1086/741537
Hui, T. Y., Crickenberger, S., Lau, J. W. T., & Williams, G. A. (2022). Why are ‘suboptimal’ temperatures preferred in a tropical intertidal ectotherm? Journal of Animal Ecology, 91(7), 1400–1415. https://doi.org/10.1111/1365-2656.13690
Hurlbert, S. H. (1984). Pseudoreplication and the design of ecological field experiments. Ecological Monographs, 54(2), 187–211. https://doi.org/10.2307/1942661
Instituto Meteorológico Nacional. (2026). Oficio IMN-DIM-CM-192-2026: Remisión de información meteorológica de las estaciones San José-Pinilla, Paquera, Santa Rita y ASADA Artola, periodo 2022–2025. Departamento de Información Meteorológica, San José, Costa Rica.
Komsta, L., & Novomestky, F. (2022). moments: Moments, cumulants, skewness, kurtosis and related tests (R package version 0.14.1). https://doi.org/10.32614/CRAN.package.moments
Kunze, C., Wölfelschneider, M., & Rölfer, L. (2021). Multiple driver impacts on rocky intertidal systems: The need for an integrated approach. Frontiers in Marine Science, 8, 667168. https://doi.org/10.3389/fmars.2021.667168
Leeuwis, R. H., & Gamperl, A. K. (2022). Adaptations and plastic phenotypic responses of marine animals to the environmental challenges of the high intertidal zone. En S. J. Hawkins et al. (Eds.), Oceanography and Marine Biology: An Annual Review (Vol. 60, pp. 625–680). CRC Press. https://doi.org/10.1201/9781003288602-13
Ley N.º 9348. (2016, 2 de marzo). Ley reguladora del Refugio Nacional de Vida Silvestre Ostional. La Gaceta N.º 43. San José, Costa Rica.
Maia, R. C., & Troncoso, J. S. (2024). The effect of extreme climatic events on littorinid snails in two estuarine environments, temperate (NW Spain) and tropical (NE Brazil). Ocean and Coastal Research, 72(suppl. 1), e24017. https://doi.org/10.1590/2675-2824072.23060
Marshall, D. J., Mustapha, N., & Monaco, C. J. (2023). Conservation of thermal physiology in tropical intertidal snails following an evolutionary transition to a cooler ecosystem: Climate change implications. Conservation Physiology, 11(1), coad056. https://doi.org/10.1093/conphys/coad056
Nielsen, E. S., Walkes, S., Sones, J. L., Fenberg, P. B., Paz-García, D. A., Grosberg, R. K., Sanford, E., & Bay, R. A. (2025). The effects of latitudinal gradients, climatic anomalies, and size-selective harvesting on the adaptive potential of an intertidal gastropod. Evolutionary Applications, 18(9), e70159. https://doi.org/10.1111/eva.70159
Ooms, J., & Denney, B. (2026). writexl: Export data frames to Excel 'xlsx' format (R package version 2.0.1). https://doi.org/10.32614/CRAN.package.writexl
R Core Team. (2026). R: A language and environment for statistical computing. R Foundation for Statistical Computing. https://www.R-project.org/
Rodríguez-Grimón, R. O., Campos, N. H., & Castro, Í. B. (2024). Efecto de la exposición al oleaje en el tamaño y el índice de condición de los neogastrópodos en Santa Marta, Colombia. Acta Biológica Colombiana, 29(1), 24–31. https://doi.org/10.15446/abc.v29n1.89120
Sistema Nacional de Áreas de Conservación. (2016). Protocolo PRONAMEC: Protocolo para el monitoreo ecológico de las playas rocosas. Proyecto Consolidación de las Áreas Marinas Protegidas, Programa de las Naciones Unidas para el Desarrollo / Fondo para el Medio Ambiente Mundial. San José, Costa Rica.
Solís Rivera, V., & Fonseca Borrás, M. (2018). El área marina de Cabo Blanco y su relación con las actividades de pesca de pequeña escala: Plan general de manejo. CoopeSoliDar R.L.
Takada, Y. (2008). Contrasting characteristics in increasing and decreasing phases of the Nerita japonica (Mollusca: Gastropoda) population on a boulder shore. Population Ecology, 50, 391–403. https://doi.org/10.1007/s10144-008-0091-7
Wickham, H., & Bryan, J. (2026). readxl: Read Excel files (R package version 1.5.0). https://doi.org/10.32614/CRAN.package.readxl
Wright, A. D., Grant, E. H. C., & Zipkin, E. F. (2022). A comparison of monitoring designs to assess wildlife community parameters across spatial scales. Ecological Applications, 32(6), e2621. https://doi.org/10.1002/eap.2621
Zarzyczny, K. M., Paz-García, D. A., Williams, S. T., Rius, M., MacLean, M., & Fenberg, P. B. (2026). Tropicalisation of rocky shore gastropod communities along the Baja California Peninsula. Hydrobiologia. https://doi.org/10.1007/s10750-026-06429-z
Anexos
Anexo A. Figura complementaria
Fig. A1. Relación entre el largo medio de la concha de Nerita scabricosta por muestreo y (A) la temperatura media y (B) la precipitación acumulada de los 90 días previos, Área de Conservación Tempisque, Costa Rica, 2022–2025 (n = 18 muestreos). Cada punto representa un muestreo. RNA: Reserva Natural Absoluta Cabo Blanco; AMM: Área Marina de Manejo Cabo Blanco.
Anexo B. Scripts de R
Scripts completos con los que se depuró la base de datos y se generaron todas las pruebas estadísticas, tablas y figuras de este informe.
B1. 01_limpieza_datos_Burgaditos.R
Depuración de la base de datos (requiere los archivos crudos; no se ejecuta al tejer este informe).
# =====================================================================
# 01_limpieza_datos.R
# Proyecto: Talla de Nerita scabricosta y régimen de protección
# Propósito: reproducir la depuración de la base biológica y de los
# datos meteorológicos a partir de los archivos crudos.
# Entradas : BDNeritas_Bioestadistica.xlsx y datos_meteorologicos_integrados.xlsx
# (se seleccionan a mano en una ventana al correr el script)
# Salidas : carpeta "datos_limpios" junto al archivo de la base biológica:
# datos_limpios/bio_limpia.csv
# datos_limpios/clima_diario.csv
# datos_limpios/analisisB_muestreos.csv
# datos_limpios/bitacora_limpieza.csv
# datos_limpios/BD_Nerita_limpia.xlsx (las cuatro tablas, una por hoja)
# Paquetes : readxl, dplyr, writexl
# (instalar una vez: install.packages(c("readxl", "dplyr", "writexl")))
# Los archivos crudos NUNCA se modifican.
# =====================================================================
library(readxl)
library(dplyr)
# Seleccionar en la ventana: 1) la base biológica cruda (BDNeritas...xlsx)
ruta_bio <- file.choose()
# 2) la base meteorológica cruda (datos_meteorologicos_integrados.xlsx)
ruta_met <- file.choose()
# Los resultados se guardan en "datos_limpios", junto a la base biológica
salida <- file.path(dirname(ruta_bio), "datos_limpios")
dir.create(salida, showWarnings = FALSE)
cat("Los archivos limpios se guardarán en:", salida, "\n")
# ---------------------------------------------------------------------
# 1. Base biológica: lectura y estructura
# ---------------------------------------------------------------------
# na = "NA": la fila "Ausente" (ID 359) trae el texto "NA" en vez de una celda vacía
bio_cruda <- read_excel(ruta_bio,
na = c("", "NA"), .name_repair = "unique_quiet")
names(bio_cruda) <- trimws(names(bio_cruda)) # "Codigo " trae un espacio
bio <- bio_cruda %>%
select(ID, Fecha, Epoca_original = Epoca, Area = `Area Conservacion`,
Sitio, Codigo, Norte_CRTM05 = Latitud, Este_CRTM05 = Longitud,
Especie, Largo = `Largo (cm)`, Ancho = `Ancho (cm)`) %>% # se descartan 3 columnas vacías/basura
mutate(Fecha = as.Date(Fecha),
Anio = as.integer(format(Fecha, "%Y")),
Mes = as.integer(format(Fecha, "%m")),
Incluir = TRUE,
Motivo = "",
Sensibilidad = FALSE)
str(bio)
bitacora <- data.frame(regla = character(), filas = integer(), accion = character())
registrar <- function(regla, filas, accion) {
bitacora[nrow(bitacora) + 1, ] <<- list(regla, filas, accion)
}
excluir <- function(datos, condicion, motivo) {
condicion[is.na(condicion)] <- FALSE
datos$Motivo[condicion & datos$Incluir] <- motivo
datos$Incluir[condicion] <- FALSE
datos
}
# ---------------------------------------------------------------------
# 2. Reglas de depuración (biológica)
# ---------------------------------------------------------------------
# R1. Solo Nerita: se excluyen Vasula melones y registros "Ausente"
r1 <- bio$Especie != "Nerita sp"
bio <- excluir(bio, r1, "R1: otra especie o sin organismo")
registrar("R1 Otra especie / Ausente", sum(r1), "Excluir")
# R2. Largo igual a 0 o vacío (valor imposible)
r2 <- bio$Incluir & (is.na(bio$Largo) | bio$Largo <= 0)
bio <- excluir(bio, r2, "R2: largo = 0 o vacío")
registrar("R2 Largo = 0 o vacío", sum(r2), "Excluir")
# R3a. Muestreos copiados dos veces en la base: si la segunda mitad (por ID)
# de un muestreo repite exactamente la primera, se excluye la copia.
# (En 2022 los IDs 1116-1474 son copia exacta de los IDs 1-359.)
copia <- rep(FALSE, nrow(bio))
for (m in unique(paste(bio$Sitio, bio$Fecha))) {
idx <- which(paste(bio$Sitio, bio$Fecha) == m & bio$Incluir)
idx <- idx[order(bio$ID[idx])]
k <- length(idx)
if (k >= 10 && k %% 2 == 0) {
a <- idx[1:(k / 2)]; b <- idx[(k / 2 + 1):k]
if (all(bio$Largo[a] == bio$Largo[b]) && all(bio$Ancho[a] == bio$Ancho[b])) {
copia[b] <- TRUE
}
}
}
bio <- excluir(bio, copia, "R3a: copia duplicada de un muestreo")
registrar("R3a Muestreos copiados dos veces", sum(copia), "Excluir la copia (IDs mayores)")
cat("IDs de la copia:", range(bio$ID[copia]), "\n")
# R3. Registros en mm (largo > 5 cm es imposible para la especie).
# Se verifica que sean duplicados de filas en cm del mismo sitio y fecha.
en_mm <- bio$Incluir & bio$Largo > 5
dup_verificados <- 0
for (i in which(en_mm)) {
candidatos <- bio %>%
filter(Sitio == bio$Sitio[i], Fecha == bio$Fecha[i], Largo <= 5, Incluir)
coincide <- any(abs(candidatos$Largo - bio$Largo[i] / 10) <= 0.011 &
abs(candidatos$Ancho - bio$Ancho[i] / 10) <= 0.011)
dup_verificados <- dup_verificados + coincide
}
cat("Filas en mm:", sum(en_mm), "| con duplicado exacto en cm:", dup_verificados, "\n")
bio <- excluir(bio, en_mm, "R3: registro en mm duplicado de uno en cm")
registrar("R3 Registros en mm duplicados", sum(en_mm), "Excluir")
# R4. Valores 10 veces menores que el resto (0 < largo < 0,5 cm),
# incoherentes y no verificables sin hojas de campo
r4 <- bio$Incluir & bio$Largo > 0 & bio$Largo < 0.5
bio <- excluir(bio, r4, "R4: valor incoherente (< 0,5 cm)")
registrar("R4 Largo < 0,5 cm", sum(r4), "Excluir")
# R5. Muestreos con ancho > largo de forma sistemática
# (razón ancho/largo mediana del muestreo > 1): el "largo" no es comparable
razon_muestreo <- bio %>%
filter(Incluir) %>%
group_by(Sitio, Fecha) %>%
summarise(razon_mediana = median(Ancho / Largo), .groups = "drop")
print(arrange(razon_muestreo, desc(razon_mediana)), n = 6)
muestreos_invertidos <- razon_muestreo %>% filter(razon_mediana > 1)
r5 <- bio$Incluir &
paste(bio$Sitio, bio$Fecha) %in% paste(muestreos_invertidos$Sitio, muestreos_invertidos$Fecha)
bio <- excluir(bio, r5, "R5: muestreo con ancho > largo sistemático")
registrar("R5 Muestreos ancho > largo sistemático", sum(r5), "Excluir")
# R6. Ancho calculado como largo/10 (fórmula, no medición): se anula el ancho
r6 <- bio$Incluir & !is.na(bio$Ancho) & abs(bio$Ancho / bio$Largo - 0.1) < 1e-6
bio$Ancho[r6] <- NA
bio$Sensibilidad[r6] <- TRUE
registrar("R6 Ancho = largo/10", sum(r6), "Ancho = NA; largo se usa, marcado para sensibilidad")
# R7. Casos aislados con ancho > largo: se mantienen, marcados
r7 <- bio$Incluir & !is.na(bio$Ancho) & bio$Ancho > bio$Largo
bio$Sensibilidad[r7] <- TRUE
registrar("R7 Ancho > largo aislado", sum(r7), "Se mantiene, marcado para sensibilidad")
# R8. Ancho anormalmente pequeño (ancho < 0.4 x largo): forma imposible para
# una concha globosa (mediana de la razón ancho/largo = 0.79). Probable
# error de digitación del ancho: se anula el ancho, se conserva el largo.
r8 <- bio$Incluir & !is.na(bio$Ancho) & bio$Ancho / bio$Largo < 0.4
cat("IDs con ancho < 0.4 x largo:", bio$ID[r8], "\n")
bio$Ancho[r8] <- NA
bio$Sensibilidad[r8] <- TRUE
registrar("R8 Ancho < 0.4 x largo", sum(r8), "Ancho = NA; largo se usa, marcado para sensibilidad")
# ---------------------------------------------------------------------
# 3. Variables para el análisis
# ---------------------------------------------------------------------
bio <- bio %>%
mutate(
Especie = "Nerita scabricosta",
# Época según la definición de la metodología: seca dic-abr, lluviosa may-nov
Epoca = factor(ifelse(Mes %in% c(12, 1:4), "Seca", "Lluviosa"),
levels = c("Seca", "Lluviosa")),
# Régimen de protección (análisis principal), según si se permite la pesca:
# restrictivo = RNA Cabo Blanco (sin pesca; Arias et al., 2015; Decreto 40442-MINAE Art. 2)
# uso regulado = AMM Cabo Blanco (pesca con licencia, Decreto 40442-MINAE Art. 7 y 9)
# y RNVS Camaronal y Ostional (pesca permitida; Arias et al., 2015)
Proteccion = factor(ifelse(Area == "RNA Cabo Blanco", "Restrictivo", "Uso regulado"),
levels = c("Restrictivo", "Uso regulado")),
# Agrupación secundaria (la original del pre-registro): complejo Cabo Blanco vs refugios
Complejo = factor(ifelse(grepl("^RNVS", Area), "Refugios", "Complejo Cabo Blanco"),
levels = c("Complejo Cabo Blanco", "Refugios")),
Area = factor(Area),
Razon_AL = round(Ancho / Largo, 3),
Analisis_A = Incluir & Anio >= 2025,
Analisis_B = Incluir & Anio %in% c(2022, 2023, 2025)
)
registrar("Época reasignada", sum(bio$Epoca != bio$Epoca_original), "Según definición seca dic-abr")
bio_limpia <- bio %>% filter(Incluir) %>%
select(ID, Fecha, Anio, Mes, Epoca, Area, Proteccion, Complejo, Sitio, Codigo,
Especie, Largo, Ancho, Razon_AL, Analisis_A, Analisis_B, Sensibilidad)
cat("\nFilas incluidas:", nrow(bio_limpia), "de", nrow(bio), "\n")
print(table(bio_limpia$Proteccion[bio_limpia$Analisis_A]))
print(table(bio_limpia$Complejo[bio_limpia$Analisis_A]))
# ---------------------------------------------------------------------
# 4. Datos meteorológicos: serie diaria
# ---------------------------------------------------------------------
met <- read_excel(ruta_met)
names(met) <- c("Estacion", "Anio", "Mes", "Dia", "Hora", "Precip", "Temp",
"Dir_viento", "Vel_viento", "Humedad", "Radiacion")
# Santa Rita se excluye: sus fechas están desfasadas (desfase creciente,
# ~10 d en 2022 a ~60 d en 2025) y la serie termina el 01/11/2025.
met <- met %>% filter(Estacion != "SANTA_RITA_NANDAYURE_GUANACASTE")
# Control de rangos (se reporta; no hubo valores fuera de rango)
print(summary(met[, c("Precip", "Temp", "Humedad")]))
# Un día es válido si tiene al menos 20 horas con dato
clima_diario <- met %>%
mutate(Fecha = as.Date(sprintf("%d-%02d-%02d", Anio, Mes, Dia))) %>%
group_by(Estacion, Fecha) %>%
summarise(horas_T = sum(!is.na(Temp)),
T_media = ifelse(horas_T >= 20, mean(Temp, na.rm = TRUE), NA),
T_max = ifelse(horas_T >= 20, max(Temp, na.rm = TRUE), NA),
T_min = ifelse(horas_T >= 20, min(Temp, na.rm = TRUE), NA),
horas_P = sum(!is.na(Precip)),
Precip_mm = ifelse(horas_P >= 20, sum(Precip, na.rm = TRUE), NA),
.groups = "drop")
clima_diario %>% group_by(Estacion) %>%
summarise(dias = n(), dias_T_validos = sum(!is.na(T_media))) %>% print()
# ---------------------------------------------------------------------
# 5. Análisis B: acople muestreo - temperatura de 90 días previos
# ---------------------------------------------------------------------
# Asignación por distancia (hoja "Distancias"): Cabo Blanco -> Paquera (~30 km),
# Ostional -> Pinilla (~25-31 km). Camaronal queda a > 50 km de toda
# estación utilizable y se excluye.
estacion_area <- c("RNA Cabo Blanco" = "PAQUERA",
"AMM Cabo Blanco" = "PAQUERA",
"RNVS Ostional" = "SAN JOSE_PINILLA",
"RNVS Camaronal" = NA)
muestreos <- bio_limpia %>%
filter(Analisis_B) %>%
group_by(Fecha, Area, Sitio) %>%
summarise(n = n(),
Largo_medio = mean(Largo),
Largo_mediana = median(Largo),
Largo_de = sd(Largo),
filas_sensibilidad = sum(Sensibilidad),
.groups = "drop") %>%
mutate(Estacion = unname(estacion_area[as.character(Area)]))
temp_90d <- function(estacion, fecha) {
if (is.na(estacion)) return(c(NA, NA))
ventana <- clima_diario %>%
filter(Estacion == estacion, Fecha >= fecha - 90, Fecha <= fecha - 1)
c(mean(ventana$T_media, na.rm = TRUE), sum(!is.na(ventana$T_media)) / 90)
}
res <- t(mapply(temp_90d, muestreos$Estacion, muestreos$Fecha))
muestreos$T_media_90d <- res[, 1]
muestreos$Cobertura <- res[, 2]
muestreos <- muestreos %>%
mutate(Estado = case_when(
is.na(Estacion) ~ "Sin estacion",
n < 10 ~ "n muy bajo",
Cobertura < 0.8 ~ "Cobertura < 80%",
TRUE ~ "Listo"))
print(muestreos %>% select(Fecha, Sitio, n, Largo_medio, Estacion, T_media_90d, Estado),
n = 40)
cat("Muestreos listos:", sum(muestreos$Estado == "Listo"),
"| valores de temperatura distintos:",
length(unique(round(muestreos$T_media_90d[muestreos$Estado == "Listo"], 1))), "\n")
# ---------------------------------------------------------------------
# 6. Guardar
# ---------------------------------------------------------------------
write.csv(bio_limpia, file.path(salida, "bio_limpia.csv"), row.names = FALSE)
write.csv(clima_diario, file.path(salida, "clima_diario.csv"), row.names = FALSE)
write.csv(muestreos, file.path(salida, "analisisB_muestreos.csv"), row.names = FALSE)
write.csv(bitacora, file.path(salida, "bitacora_limpieza.csv"), row.names = FALSE)
# Versión en Excel de la base actualizada (ítem 28 de la rúbrica):
# un solo archivo con una hoja por tabla
library(writexl)
write_xlsx(list(bio_limpia = bio_limpia,
analisisB_muestreos = muestreos,
clima_diario = clima_diario,
bitacora_limpieza = bitacora),
file.path(salida, "BD_Nerita_limpia.xlsx"))
print(bitacora)B2. 02_descriptiva_Burgaditos.R
Estadística descriptiva, Tabla 1, Tabla 2 y figuras 1 y 3.
# =====================================================================
# 02_descriptiva.R · Los Burgaditos
# Estadística descriptiva de la talla de Nerita scabricosta (2022-2026)
# y del clima de los 90 días previos a cada muestreo.
#
# Datos: BD_Nerita_limpia.xlsx (generado por 01_limpieza_datos.R).
# Al correr el script se abre una ventana para seleccionar ese archivo.
# Resultados: carpeta "descriptiva" junto al Excel, con
# - Descriptiva_Nerita.xlsx (todas las tablas, una por hoja)
# - figuras PNG a 330 dpi, en escala de grises (formato RBT)
#
# Paquetes: readxl (obligatorio); writexl y moments (opcionales).
# install.packages(c("readxl", "writexl", "moments"))
# =====================================================================
library(readxl)
# ---------------------------------------------------------------------
# 0. Preparación
# ---------------------------------------------------------------------
# Asimetría y curtosis: paquete moments si está; si no, mismas fórmulas.
if (requireNamespace("moments", quietly = TRUE)) {
skewness <- moments::skewness
kurtosis <- moments::kurtosis
} else {
message("Paquete 'moments' no disponible: se usan fórmulas equivalentes.")
skewness <- function(x) { x <- x[!is.na(x)]; d <- x - mean(x)
sqrt(length(x)) * sum(d^3) / sum(d^2)^1.5 }
kurtosis <- function(x) { x <- x[!is.na(x)]; d <- x - mean(x) # normal = 3
length(x) * sum(d^4) / sum(d^2)^2 }
}
# Seleccionar BD_Nerita_limpia.xlsx en la ventana
ruta_bd <- if (exists("RUTA_BD")) RUTA_BD else file.choose()
bio <- as.data.frame(read_excel(ruta_bd, sheet = "bio_limpia"))
clima <- as.data.frame(read_excel(ruta_bd, sheet = "clima_diario"))
bio$Fecha <- as.Date(bio$Fecha)
clima$Fecha <- as.Date(clima$Fecha)
# Etiquetas cortas y orden fijo de los factores
etiquetas <- c("RNA Cabo Blanco" = "RNA", "AMM Cabo Blanco" = "AMM",
"RNVS Camaronal" = "Camaronal", "RNVS Ostional" = "Ostional")
bio$Area <- factor(etiquetas[bio$Area], levels = etiquetas)
bio$Anio <- factor(bio$Anio)
bio$Epoca <- factor(bio$Epoca, levels = c("Seca", "Lluviosa"))
bio$Muestreo <- paste(bio$Sitio, bio$Fecha)
carpeta <- file.path(dirname(ruta_bd), "descriptiva")
dir.create(carpeta, showWarnings = FALSE)
guardar_figura <- function(archivo, dibujo, ancho_cm = 16.5, alto_cm = 12) {
png(file.path(carpeta, archivo), width = ancho_cm / 2.54, height = alto_cm / 2.54,
units = "in", res = 330, pointsize = 9)
dibujo()
dev.off()
cat("Figura guardada:", archivo, "\n")
}
tablas <- list() # aquí se juntan todas las tablas para exportar al final
# ---------------------------------------------------------------------
# 1. Diseño de muestreo: cuántos individuos y muestreos hay en cada celda
# (muestra qué combinaciones de factores existen y cuáles faltan)
# ---------------------------------------------------------------------
n_ind <- as.data.frame.matrix(table(bio$Area, interaction(bio$Anio, bio$Epoca, sep = " ")))
muestreos <- unique(bio[, c("Muestreo", "Area", "Anio", "Epoca")])
n_mue <- as.data.frame.matrix(table(muestreos$Area,
interaction(muestreos$Anio, muestreos$Epoca, sep = " ")))
cat("\n=== Individuos por área, año y época ===\n"); print(n_ind)
cat("\n=== Muestreos (sitio-fecha) por área, año y época ===\n"); print(n_mue)
tablas$Diseno_individuos <- cbind(Area = rownames(n_ind), n_ind)
tablas$Diseno_muestreos <- cbind(Area = rownames(n_mue), n_mue)
# ---------------------------------------------------------------------
# 2. Tendencia central, dispersión y forma (largo y ancho)
# ---------------------------------------------------------------------
descriptiva <- function(x) {
x <- x[!is.na(x)]
if (length(x) < 2) { # grupo sin datos (p. ej. ancho de Punta Flor 2023)
return(c(n = length(x), Media = NA, DE = NA, EE = NA, Mediana = NA, Q1 = NA, Q3 = NA,
Min = NA, Max = NA, CV = NA, Asimetria = NA, Curtosis = NA))
}
c(n = length(x), Media = mean(x), DE = sd(x), EE = sd(x) / sqrt(length(x)),
Mediana = median(x), Q1 = unname(quantile(x, 0.25)), Q3 = unname(quantile(x, 0.75)),
Min = min(x), Max = max(x), CV = sd(x) / mean(x) * 100,
Asimetria = skewness(x), Curtosis = kurtosis(x))
}
# Aplica descriptiva() a cada combinación de los factores indicados
por_grupo <- function(variable, factores) {
grupos <- split(bio[[variable]], bio[factores], drop = TRUE, sep = " | ")
res <- t(sapply(grupos, descriptiva))
data.frame(Grupo = rownames(res), round(res, 3), row.names = NULL)
}
for (v in c("Largo", "Ancho")) {
cat("\n\n############", v, "(cm) ############\n")
cat("\n--- Total ---\n"); print(round(descriptiva(bio[[v]]), 3))
for (f in list("Area", "Anio", "Epoca", c("Area", "Anio"), c("Area", "Anio", "Epoca"))) {
nombre <- paste(v, paste(f, collapse = "_"), sep = "_")
tablas[[nombre]] <- por_grupo(v, f)
cat("\n--- Por", paste(f, collapse = " x "), "---\n"); print(tablas[[nombre]])
}
}
cat("\nNota: el ancho tiene", sum(is.na(bio$Ancho)),
"valores faltantes (reglas R6 y R8 de la limpieza: ancho calculado como largo/10",
"o ancho < 0.4 x largo).\n")
# ---------------------------------------------------------------------
# 3. Distribución de frecuencias del largo (intervalos de Sturges,
# comunes para todas las áreas)
# ---------------------------------------------------------------------
k_sturges <- nclass.Sturges(bio$Largo) # 1 + log2(n)
k_scott <- nclass.scott(bio$Largo)
cat("\nNúmero de clases - Sturges:", k_sturges, "| Scott:", k_scott, "\n")
amplitud <- ceiling(diff(range(bio$Largo)) / k_sturges * 20) / 20 # redondeo a 0.05 cm
cortes <- seq(floor(min(bio$Largo) * 20) / 20, by = amplitud, length.out = k_sturges + 2)
cortes <- cortes[seq_len(which(cortes >= max(bio$Largo))[1])]
cat("Amplitud de clase:", amplitud, "cm | límites:", cortes, "\n")
tabla_frecuencias <- function(x) {
f <- table(cut(x, cortes, include.lowest = TRUE, right = FALSE, dig.lab = 3))
data.frame(Clase = names(f),
Marca = round(head(cortes, -1) + amplitud / 2, 3),
f = as.vector(f),
fr = round(as.vector(f) / length(x), 3),
F_acum = cumsum(as.vector(f)),
Fr_acum = round(cumsum(as.vector(f)) / length(x), 3))
}
tablas$Frecuencias_Total <- tabla_frecuencias(bio$Largo)
for (a in levels(bio$Area)) {
tablas[[paste0("Frecuencias_", a)]] <- tabla_frecuencias(bio$Largo[bio$Area == a])
}
print(tablas$Frecuencias_Total)
# ---------------------------------------------------------------------
# 4. Figuras (sin título interno: en la revista el texto va en el pie)
# ---------------------------------------------------------------------
# Fig. D1. Histogramas de frecuencia relativa por área, mismos intervalos y escala
fr_max <- max(sapply(levels(bio$Area), function(a)
max(hist(bio$Largo[bio$Area == a], breaks = cortes, plot = FALSE, right = FALSE)$counts /
sum(bio$Area == a))))
guardar_figura("Fig_D1_histogramas_area.png", alto_cm = 14, function() {
par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))
for (a in levels(bio$Area)) {
x <- bio$Largo[bio$Area == a]
h <- hist(x, breaks = cortes, plot = FALSE, right = FALSE)
h$counts <- h$counts / length(x) * 100
plot(h, main = "", xlab = "Largo de la concha (cm)", ylab = "Frecuencia relativa (%)",
ylim = c(0, fr_max * 100), col = "grey75", border = "white")
mtext(paste0(a, " (n = ", length(x), ")"), side = 3, adj = 0, line = 0.3, font = 2)
}
})
# Fig. D2. Polígonos de frecuencia relativa superpuestos (comparación de formas)
guardar_figura("Fig_D2_poligonos_area.png", function() {
par(mar = c(4, 4, 1, 1))
marcas <- head(cortes, -1) + amplitud / 2
plot(NA, xlim = range(cortes), ylim = c(0, fr_max * 100),
xlab = "Largo de la concha (cm)", ylab = "Frecuencia relativa (%)")
tipos <- 1:4; puntos <- c(19, 1, 17, 2)
for (i in seq_along(levels(bio$Area))) {
a <- levels(bio$Area)[i]; x <- bio$Largo[bio$Area == a]
fr <- hist(x, breaks = cortes, plot = FALSE, right = FALSE)$counts / length(x) * 100
lines(marcas, fr, lty = tipos[i]); points(marcas, fr, pch = puntos[i])
}
legend("topright", legend = levels(bio$Area), lty = tipos, pch = puntos, bty = "n")
})
# Fig. D3. Diagramas de caja por área y año (las celdas sin datos quedan vacías)
guardar_figura("Fig_D3_cajas_area_anio.png", function() {
par(mar = c(5, 4, 1, 1))
grupo <- interaction(bio$Anio, bio$Area, sep = "\n", lex.order = FALSE)
grupo <- factor(grupo, levels = as.vector(t(outer(levels(bio$Area), levels(bio$Anio),
function(a, y) paste(y, a, sep = "\n")))))
boxplot(bio$Largo ~ grupo, drop = FALSE, las = 1, xaxt = "n", xlab = "",
ylab = "Largo de la concha (cm)", col = "grey85", outcex = 0.6)
niveles <- levels(grupo)
axis(1, at = seq_along(niveles), labels = sub("\n.*", "", niveles), cex.axis = 0.8)
pos <- tapply(seq_along(niveles), sub(".*\n", "", niveles), mean)
mtext(names(pos), side = 1, at = pos, line = 2.5, font = 2)
abline(v = seq(length(levels(bio$Anio)) + 0.5, by = length(levels(bio$Anio)),
length.out = 3), lty = 3, col = "grey50")
})
# Fig. D4. Diagramas de caja por época en los únicos casos con ambas épocas
# en la misma área y año
estratos <- with(bio, tapply(Epoca, paste(Area, Anio), function(e) length(unique(e)) == 2))
estratos <- names(estratos)[estratos]
cat("\nÁrea-año con ambas épocas:", estratos, "\n")
guardar_figura("Fig_D4_cajas_epoca.png", ancho_cm = 12, function() {
par(mar = c(5, 4, 1, 1))
sub_e <- bio[paste(bio$Area, bio$Anio) %in% estratos, ]
sub_e <- sub_e[order(sub_e$Area, sub_e$Anio, sub_e$Epoca), ]
grupo <- paste(sub_e$Area, sub_e$Anio, sub_e$Epoca, sep = "\n")
grupo <- factor(grupo, levels = unique(grupo))
boxplot(sub_e$Largo ~ grupo, las = 1, xlab = "", xaxt = "n",
ylab = "Largo de la concha (cm)", outcex = 0.6,
col = ifelse(grepl("Seca", levels(grupo)), "white", "grey70"))
axis(1, at = seq_along(levels(grupo)), labels = levels(grupo), padj = 1,
cex.axis = 0.8) # blanco = seca, gris = lluviosa
})
# ---------------------------------------------------------------------
# 5. Ámbito de variación permisible PRONAMEC (SINAC, 2016, Indicador 2)
# El protocolo evalúa el PROMEDIO de los individuos de cada muestreo:
# Bueno: > 2.5 a 4 cm de largo (> 2.5 a 3.5 cm de ancho)
# Regular: 1.5 a 2.5 cm · Malo: < 1 cm
# 1 a < 1.5 cm no tiene categoría en el protocolo (vacío del protocolo)
# ---------------------------------------------------------------------
categoria_pronamec <- function(x, tope_bueno) {
cut(x, breaks = c(-Inf, 1, 1.5, 2.5, tope_bueno, Inf), right = FALSE,
labels = c("Malo", "Sin categoría", "Regular", "Bueno", "Sobre el ámbito"))
}
# Nota: con right = FALSE, 2.5 exacto queda como "Bueno"; el protocolo
# dice "sobre 2.5" para Bueno y "1.5 a 2.5" para Regular. Revisen si algún
# promedio cae exactamente en 2.5 (se imprime abajo).
resumen_mue <- aggregate(cbind(Largo, Ancho) ~ Muestreo + Area + Anio + Epoca + Sitio + Fecha,
data = bio, FUN = mean, na.action = na.pass, na.rm = TRUE)
names(resumen_mue)[7:8] <- c("Largo_medio", "Ancho_medio")
resumen_mue$n <- as.vector(table(bio$Muestreo)[resumen_mue$Muestreo])
resumen_mue$Largo_max <- as.vector(tapply(bio$Largo, bio$Muestreo, max)[resumen_mue$Muestreo])
resumen_mue$Cat_largo <- categoria_pronamec(resumen_mue$Largo_medio, 4)
resumen_mue$Cat_ancho <- categoria_pronamec(resumen_mue$Ancho_medio, 3.5)
resumen_mue <- resumen_mue[order(resumen_mue$Area, resumen_mue$Fecha), ]
print(resumen_mue[, c("Area", "Sitio", "Fecha", "n", "Largo_medio", "Ancho_medio",
"Cat_largo", "Cat_ancho")], digits = 3)
cat("Promedios exactamente en un límite (1, 1.5 o 2.5 cm):",
sum(round(resumen_mue$Largo_medio, 3) %in% c(1, 1.5, 2.5)), "\n")
cat("Muestreos con n < 10 (su promedio es poco confiable):",
sum(resumen_mue$n < 10), "\n")
tablas$PRONAMEC_por_muestreo <- data.frame(
resumen_mue[, c("Area", "Sitio", "Fecha", "Anio", "Epoca", "n")],
round(resumen_mue[, c("Largo_medio", "Ancho_medio", "Largo_max")], 3),
resumen_mue[, c("Cat_largo", "Cat_ancho")])
tab_cat <- table(resumen_mue$Area, resumen_mue$Cat_largo)
cat("\n=== Muestreos por área y categoría PRONAMEC (largo medio) ===\n"); print(tab_cat)
tablas$PRONAMEC_area_categoria <- cbind(Area = rownames(tab_cat), as.data.frame.matrix(tab_cat))
tab_cat_anio <- table(resumen_mue$Anio, resumen_mue$Cat_largo)
tablas$PRONAMEC_anio_categoria <- cbind(Anio = rownames(tab_cat_anio),
as.data.frame.matrix(tab_cat_anio))
# ---------------------------------------------------------------------
# 6. Clima de los 90 días previos a cada muestreo
# Misma regla que el pre-registro: días -90 a -1, al menos 80 % de días
# válidos. Camaronal no tiene estación utilizable; no hay datos del IMN
# para 2026.
# ---------------------------------------------------------------------
estacion_area <- c(RNA = "PAQUERA", AMM = "PAQUERA",
Ostional = "SAN JOSE_PINILLA", Camaronal = NA)
clima_90d <- function(estacion, fecha) {
if (is.na(estacion)) return(c(NA, NA, NA, NA))
v <- clima[clima$Estacion == estacion & clima$Fecha >= fecha - 90 &
clima$Fecha <= fecha - 1, ]
cobT <- sum(!is.na(v$T_media)) / 90; cobP <- sum(!is.na(v$Precip_mm)) / 90
c(T_media_90d = if (cobT >= 0.8) mean(v$T_media, na.rm = TRUE) else NA,
# precipitación acumulada estimada = promedio diario x 90 días
P_acum_90d = if (cobP >= 0.8) mean(v$Precip_mm, na.rm = TRUE) * 90 else NA,
Cobertura_T = cobT, Cobertura_P = cobP)
}
resumen_mue$Estacion <- unname(estacion_area[as.character(resumen_mue$Area)])
cl <- t(mapply(clima_90d, resumen_mue$Estacion, resumen_mue$Fecha))
resumen_mue <- cbind(resumen_mue, round(cl, 3))
con_clima <- subset(resumen_mue, !is.na(T_media_90d) & !is.na(P_acum_90d) & n >= 10)
cat("\nMuestreos con clima completo y n >= 10:", nrow(con_clima), "\n")
tablas$Clima_por_muestreo <- resumen_mue[, c("Area", "Sitio", "Fecha", "n", "Estacion",
"Largo_medio", "T_media_90d", "P_acum_90d",
"Cobertura_T", "Cobertura_P")]
clima_desc <- do.call(rbind, lapply(split(con_clima, paste(con_clima$Estacion, con_clima$Anio)),
function(d) data.frame(Estacion = d$Estacion[1], Anio = d$Anio[1], Muestreos = nrow(d),
T_media = round(mean(d$T_media_90d), 2), T_DE = round(sd(d$T_media_90d), 2),
P_media = round(mean(d$P_acum_90d), 1), P_DE = round(sd(d$P_acum_90d), 1))))
print(clima_desc, row.names = FALSE)
tablas$Clima_resumen <- clima_desc
cat("Valores distintos de temperatura:", length(unique(round(con_clima$T_media_90d, 1))),
"| de precipitación:", length(unique(round(con_clima$P_acum_90d))), "\n")
cat("Correlación (Pearson) entre temperatura y precipitación de 90 d:",
round(cor(con_clima$T_media_90d, con_clima$P_acum_90d), 2),
"(si es alta, las dos variables dicen casi lo mismo en la regresión múltiple)\n")
# Fig. D5. Talla media del muestreo frente a cada variable climática
guardar_figura("Fig_D5_clima_talla.png", alto_cm = 8, function() {
par(mfrow = c(1, 2), mar = c(4, 4, 1, 1))
simbolo <- c(RNA = 19, AMM = 1, Ostional = 2)[as.character(con_clima$Area)]
plot(Largo_medio ~ T_media_90d, data = con_clima, pch = simbolo,
xlab = "Temperatura media de 90 d (\u00b0C)", ylab = "Largo medio (cm)")
plot(Largo_medio ~ P_acum_90d, data = con_clima, pch = simbolo,
xlab = "Precipitaci\u00f3n acumulada de 90 d (mm)", ylab = "Largo medio (cm)")
legend("topright", legend = c("RNA", "AMM", "Ostional"), pch = c(19, 1, 2), bty = "n")
})
# ---------------------------------------------------------------------
# 7. Tablas para el informe, ya en formato RBT (punto decimal, media ± DE)
# ---------------------------------------------------------------------
# Tabla 1. Descriptiva del largo y el ancho por área (y total)
fila_tabla1 <- function(x) {
d <- descriptiva(x)
c(n = d[["n"]],
"Media ± DE" = sprintf("%.2f ± %.2f", d[["Media"]], d[["DE"]]),
Mediana = sprintf("%.2f", d[["Mediana"]]),
"Mín–Máx" = sprintf("%.2f–%.2f", d[["Min"]], d[["Max"]]),
"CV (%)" = sprintf("%.1f", d[["CV"]]),
"Asimetría" = sprintf("%.2f", d[["Asimetria"]]),
Curtosis = sprintf("%.2f", d[["Curtosis"]]))
}
tabla1 <- do.call(rbind, lapply(c("Largo", "Ancho"), function(v) {
filas <- rbind(t(sapply(levels(bio$Area), function(a) fila_tabla1(bio[[v]][bio$Area == a]))),
Total = fila_tabla1(bio[[v]]))
data.frame(Variable = paste(v, "(cm)"), "Área" = rownames(filas), filas,
check.names = FALSE, row.names = NULL)
}))
cat("\n=== Tabla 1 ===\n"); print(tabla1, row.names = FALSE)
tablas$Tabla1_informe <- tabla1
# Tabla 2. Muestreos por área y categoría PRONAMEC (largo medio)
# Decisión del grupo (pre-registro v3): se excluyen los muestreos con menos de
# 10 individuos, porque su promedio no es representativo.
n_min <- 10
mue_t2 <- subset(resumen_mue, n >= n_min)
cat("\nMuestreos excluidos de la Tabla 2 por n <", n_min, ":",
nrow(resumen_mue) - nrow(mue_t2), "\n")
categorias <- c("Bueno", "Regular", "Sin categoría", "Malo", "Sobre el ámbito")
t2 <- table(mue_t2$Area, factor(mue_t2$Cat_largo, levels = categorias))
t2 <- t2[, colSums(t2) > 0 | colnames(t2) %in% c("Bueno", "Regular", "Malo"), drop = FALSE]
tabla2 <- data.frame("Área" = rownames(t2), as.data.frame.matrix(t2),
Total = rowSums(t2), check.names = FALSE, row.names = NULL)
tabla2 <- rbind(tabla2, data.frame("Área" = "Total", t(colSums(tabla2[, -1])),
check.names = FALSE))
cat("\n=== Tabla 2 ===\n"); print(tabla2, row.names = FALSE)
tablas$Tabla2_informe <- tabla2
# ---------------------------------------------------------------------
# 8. Exportar todas las tablas a un solo Excel
# ---------------------------------------------------------------------
# Las tablas del informe van primero en el Excel
tablas <- tablas[c("Tabla1_informe", "Tabla2_informe",
setdiff(names(tablas), c("Tabla1_informe", "Tabla2_informe")))]
names(tablas) <- substr(names(tablas), 1, 31) # Excel acepta máx. 31 caracteres
if (requireNamespace("writexl", quietly = TRUE)) {
writexl::write_xlsx(tablas, file.path(carpeta, "Descriptiva_Nerita.xlsx"))
cat("\nTablas guardadas en:", file.path(carpeta, "Descriptiva_Nerita.xlsx"), "\n")
} else {
for (t in names(tablas)) write.csv(tablas[[t]], file.path(carpeta, paste0(t, ".csv")),
row.names = FALSE)
cat("\n'writexl' no está instalado: tablas guardadas como CSV en", carpeta, "\n")
}
# La hoja Clima_por_muestreo la usará el script inferencial (regresión múltiple).B3. 03_inferencial_Burgaditos.R
Estadística inferencial (regla de decisión del pre-registro), Fig. 2 y Fig. A1.
# =====================================================================
# 03_inferencial.R · Los Burgaditos
# Pruebas de hipótesis según el pre-registro v3 (26/10/2026, hoja
# "Plan v3 (vigente)"). Correr DESPUÉS de 02_descriptiva.R.
#
# P1 Área (4 niveles), todo el periodo 2022-2026
# P2 Año dentro de cada área
# P3 Época dentro de área-año (RNA 2022, AMM 2022, Ostional 2025)
# P4 Correlación (Pearson o Spearman): largo medio del muestreo y temperatura /
# precipitación de 90 d (cambio pedido por la docente el 02/11/2026)
# P5 Prueba exacta de Fisher: categoría PRONAMEC x área (muestreos con n >= 10)
#
# Regla de decisión (fijada antes de correr):
# - 3 o más grupos: residuos normales (Shapiro-Wilk) y varianzas homogéneas
# (Levene si está el paquete car; si no, Fligner-Killeen) -> ANOVA + Tukey.
# Normal pero varianzas distintas -> ANOVA de Welch + pairwise.t.test (Bonferroni).
# No normal -> Kruskal-Wallis + pairwise.wilcox.test (Bonferroni).
# - 2 grupos: normalidad por grupo + var.test -> t de Student o de Welch;
# no normal -> Wilcoxon (Mann-Whitney) con estimador de Hodges-Lehmann.
#
# Datos: BD_Nerita_limpia.xlsx (se selecciona en la ventana que se abre).
# Resultados: carpeta "inferencial" junto al Excel (figuras + Inferencial_Nerita.xlsx)
# Paquetes: readxl (obligatorio); car y writexl (opcionales).
# =====================================================================
library(readxl)
alpha <- 0.05
tiene_car <- requireNamespace("car", quietly = TRUE)
if (!tiene_car) message("Paquete 'car' no disponible: la homogeneidad de varianzas ",
"se evalúa con fligner.test() (Fligner-Killeen).")
# ---------------------------------------------------------------------
# 0. Datos
# ---------------------------------------------------------------------
ruta_bd <- if (exists("RUTA_BD")) RUTA_BD else file.choose()
bio <- as.data.frame(read_excel(ruta_bd, sheet = "bio_limpia"))
clima <- as.data.frame(read_excel(ruta_bd, sheet = "clima_diario"))
bio$Fecha <- as.Date(bio$Fecha)
clima$Fecha <- as.Date(clima$Fecha)
etiquetas <- c("RNA Cabo Blanco" = "RNA", "AMM Cabo Blanco" = "AMM",
"RNVS Camaronal" = "Camaronal", "RNVS Ostional" = "Ostional")
bio$Area <- factor(etiquetas[bio$Area], levels = etiquetas)
bio$Anio <- factor(bio$Anio)
bio$Epoca <- factor(bio$Epoca, levels = c("Seca", "Lluviosa"))
bio$Muestreo <- paste(bio$Sitio, bio$Fecha)
carpeta <- file.path(dirname(ruta_bd), "inferencial")
dir.create(carpeta, showWarnings = FALSE)
guardar_figura <- function(archivo, dibujo, ancho_cm = 16.5, alto_cm = 12) {
png(file.path(carpeta, archivo), width = ancho_cm / 2.54, height = alto_cm / 2.54,
units = "in", res = 330, pointsize = 9)
dibujo(); dev.off()
cat("Figura guardada:", archivo, "\n")
}
resumen <- data.frame() # una fila por prueba
posthoc <- data.frame() # comparaciones por pares
anotar <- function(..., p_Bonferroni = NA)
resumen <<- rbind(resumen, data.frame(..., p_Bonferroni = p_Bonferroni, check.names = FALSE))
# Convierte la matriz de p de pairwise.*.test a formato largo
pares_a_tabla <- function(pp, pregunta, grupo) {
m <- pp$p.value
idx <- which(!is.na(m), arr.ind = TRUE)
data.frame(Pregunta = pregunta, Grupo = grupo, Metodo = pp$method,
Comparacion = paste(rownames(m)[idx[, 1]], "-", colnames(m)[idx[, 2]]),
p_ajustado = signif(m[idx], 3), row.names = NULL)
}
# ---------------------------------------------------------------------
# Función general: compara la variable entre los niveles de un factor
# aplicando la regla de decisión del pre-registro
# ---------------------------------------------------------------------
comparar <- function(datos, variable, factor, pregunta, grupo = "Todos") {
d <- datos[!is.na(datos[[variable]]), ]
d[[factor]] <- droplevels(factor(d[[factor]]))
y <- d[[variable]]; g <- d[[factor]]; k <- nlevels(g)
cat("\n================================================================\n")
cat(pregunta, "|", grupo, "|", variable, "~", factor, "| n =", nrow(d), "\n")
print(tapply(y, g, function(x) round(c(n = length(x), media = mean(x),
DE = sd(x), mediana = median(x)), 3)))
if (k == 2) {
sh <- tapply(y, g, function(x) if (length(x) >= 3) shapiro.test(x)$p.value else NA)
vt <- var.test(y ~ g)$p.value
cat("Shapiro-Wilk por grupo (p):", round(sh, 4), "| var.test (p):", signif(vt, 3), "\n")
if (all(sh > alpha, na.rm = TRUE)) {
prueba <- t.test(y ~ g, var.equal = vt > alpha)
nombre <- if (vt > alpha) "t de Student" else "t de Welch"
efecto <- sprintf("Dif. de medias = %.2f cm (IC 95 %%: %.2f a %.2f)",
-diff(prueba$estimate), prueba$conf.int[1], prueba$conf.int[2])
est <- sprintf("t = %.2f", prueba$statistic); gl <- round(prueba$parameter, 1)
} else {
prueba <- wilcox.test(y ~ g, conf.int = TRUE, exact = FALSE)
nombre <- "Wilcoxon (Mann-Whitney)"
efecto <- sprintf("Dif. de ubicación (Hodges-Lehmann) = %.2f cm (IC 95 %%: %.2f a %.2f)",
prueba$estimate, prueba$conf.int[1], prueba$conf.int[2])
est <- sprintf("W = %.0f", prueba$statistic); gl <- NA
}
print(prueba)
anotar(Pregunta = pregunta, Grupo = grupo, Variable = variable, Factor = factor,
n = nrow(d), Niveles = paste(levels(g), collapse = ", "), Prueba = nombre,
Estadistico = est, gl = as.character(gl), p = signif(prueba$p.value, 3),
Efecto = efecto)
return(invisible(prueba$p.value))
}
# 3 o más grupos
modelo <- aov(y ~ g)
res <- residuals(modelo)
p_norm <- shapiro.test(if (length(res) > 5000) sample(res, 5000) else res)$p.value
if (tiene_car) {
p_var <- car::leveneTest(y ~ g, center = mean)[1, "Pr(>F)"]; nom_var <- "Levene"
} else {
p_var <- fligner.test(y ~ g)$p.value; nom_var <- "Fligner-Killeen"
}
cat("Shapiro-Wilk de residuos (p):", signif(p_norm, 3), "|", nom_var, "(p):",
signif(p_var, 3), "\n")
if (p_norm > alpha && p_var > alpha) {
tab <- summary(modelo)[[1]]; print(tab)
nombre <- "ANOVA de una vía"
est <- sprintf("F = %.2f", tab[1, "F value"])
gl <- paste(tab[1, "Df"], tab[2, "Df"], sep = ", ")
p <- tab[1, "Pr(>F)"]
efecto <- sprintf("η² = %.3f", tab[1, "Sum Sq"] / sum(tab[, "Sum Sq"]))
tk <- TukeyHSD(modelo)$g; print(round(tk, 4))
ph <- data.frame(Pregunta = pregunta, Grupo = grupo, Metodo = "Tukey HSD",
Comparacion = rownames(tk), p_ajustado = signif(tk[, "p adj"], 3),
row.names = NULL)
} else if (p_norm > alpha) {
w <- oneway.test(y ~ g, var.equal = FALSE); print(w)
nombre <- "ANOVA de Welch"
est <- sprintf("F = %.2f", w$statistic)
gl <- paste(round(w$parameter, 1), collapse = ", "); p <- w$p.value
tab <- summary(modelo)[[1]]
efecto <- sprintf("η² = %.3f", tab[1, "Sum Sq"] / sum(tab[, "Sum Sq"]))
pp <- pairwise.t.test(y, g, pool.sd = FALSE, p.adjust.method = "bonferroni"); print(pp)
ph <- pares_a_tabla(pp, pregunta, grupo)
} else {
kw <- kruskal.test(y ~ g); print(kw)
nombre <- "Kruskal-Wallis"
est <- sprintf("H = %.2f", kw$statistic); gl <- as.character(kw$parameter)
p <- kw$p.value
efecto <- sprintf("ε² = %.3f", kw$statistic / (length(y) - 1))
pp <- pairwise.wilcox.test(y, g, p.adjust.method = "bonferroni", exact = FALSE); print(pp)
ph <- pares_a_tabla(pp, pregunta, grupo)
}
posthoc <<- rbind(posthoc, cbind(ph[, 1:2], Variable = variable, ph[, -(1:2)]))
anotar(Pregunta = pregunta, Grupo = grupo, Variable = variable, Factor = factor,
n = nrow(d), Niveles = paste(levels(g), collapse = ", "), Prueba = nombre,
Estadistico = est, gl = gl, p = signif(p, 3), Efecto = efecto)
invisible(p)
}
# ---------------------------------------------------------------------
# P1. Área, todo el periodo
# ---------------------------------------------------------------------
comparar(bio, "Largo", "Area", "P1 Área")
comparar(subset(bio, !Sensibilidad), "Largo", "Area", "P1 Área (sensibilidad)")
comparar(bio, "Ancho", "Area", "P1 Área (ancho, complemento)")
guardar_figura("Fig_I1_cajas_area.png", ancho_cm = 12, function() {
par(mar = c(4, 4, 1, 1))
boxplot(Largo ~ Area, data = bio, xlab = "", ylab = "Largo de la concha (cm)",
col = "grey85", outcex = 0.6, las = 1,
names = paste0(levels(bio$Area), "\n(n = ", table(bio$Area), ")"), padj = 0.5)
})
# ---------------------------------------------------------------------
# P2. Año dentro de cada área (Bonferroni entre las 4 áreas)
# ---------------------------------------------------------------------
correr_familia <- function(subconjuntos, variable, factor, pregunta) {
inicio <- nrow(resumen) + 1
for (nombre in names(subconjuntos)) comparar(subconjuntos[[nombre]], variable, factor,
pregunta, nombre)
filas <- inicio:nrow(resumen)
resumen$p_Bonferroni[filas] <<- signif(pmin(1, resumen$p[filas] * length(filas)), 3)
}
por_area <- split(bio, bio$Area)
correr_familia(por_area, "Largo", "Anio", "P2 Año")
correr_familia(lapply(por_area, subset, !Sensibilidad), "Largo", "Anio", "P2 Año (sensibilidad)")
correr_familia(por_area, "Ancho", "Anio", "P2 Año (ancho, complemento)")
# ---------------------------------------------------------------------
# P3. Época dentro de área-año (Bonferroni entre los 3 estratos)
# ---------------------------------------------------------------------
estratos <- list("RNA 2022" = subset(bio, Area == "RNA" & Anio == "2022"),
"AMM 2022" = subset(bio, Area == "AMM" & Anio == "2022"),
"Ostional 2025" = subset(bio, Area == "Ostional" & Anio == "2025"))
cat("\nSitios por época en cada estrato (limitación declarada en el pre-registro):\n")
for (e in names(estratos)) { cat(e, "\n"); print(table(estratos[[e]]$Sitio, estratos[[e]]$Epoca)) }
correr_familia(estratos, "Largo", "Epoca", "P3 Época")
correr_familia(lapply(estratos, subset, !Sensibilidad), "Largo", "Epoca", "P3 Época (sensibilidad)")
correr_familia(estratos, "Ancho", "Epoca", "P3 Época (ancho, complemento)")
# ---------------------------------------------------------------------
# P4. Correlación entre el largo medio del muestreo y el clima de 90 d
# (cambio solicitado por la docente el 02/11/2026: correlación en vez
# de regresión, porque no se busca predecir la talla a partir del clima)
# Unidad = muestreo (todos los caracoles comparten el mismo clima)
# Regla: si ambas variables son normales (Shapiro-Wilk) -> Pearson;
# si alguna no lo es -> Spearman. La otra se reporta como complemento.
# ---------------------------------------------------------------------
mue <- aggregate(Largo ~ Muestreo + Area + Sitio + Fecha, data = bio, FUN = mean)
names(mue)[names(mue) == "Largo"] <- "Largo_medio"
mue$n <- as.vector(table(bio$Muestreo)[mue$Muestreo])
estacion_area <- c(RNA = "PAQUERA", AMM = "PAQUERA", Ostional = "SAN JOSE_PINILLA", Camaronal = NA)
clima_90d <- function(estacion, fecha) {
if (is.na(estacion)) return(c(NA, NA))
v <- clima[clima$Estacion == estacion & clima$Fecha >= fecha - 90 & clima$Fecha <= fecha - 1, ]
c(if (sum(!is.na(v$T_media)) / 90 >= 0.8) mean(v$T_media, na.rm = TRUE) else NA,
if (sum(!is.na(v$Precip_mm)) / 90 >= 0.8) mean(v$Precip_mm, na.rm = TRUE) * 90 else NA)
}
cl <- t(mapply(clima_90d, unname(estacion_area[as.character(mue$Area)]), mue$Fecha))
mue$T_90d <- cl[, 1]; mue$P_90d <- cl[, 2]
datos_P4 <- subset(mue, !is.na(T_90d) & !is.na(P_90d) & n >= 10)
cat("\n================================================================\n")
cat("P4 Clima | muestreos:", nrow(datos_P4), "\n")
print(datos_P4[order(datos_P4$Fecha), c("Area", "Sitio", "Fecha", "n", "Largo_medio", "T_90d", "P_90d")],
row.names = FALSE, digits = 4)
# Supuesto de normalidad de cada variable (Pearson requiere normalidad bivariada)
normal_P4 <- sapply(c("Largo_medio", "T_90d", "P_90d"), function(v) {
s <- shapiro.test(datos_P4[[v]])
cat(sprintf("Shapiro-Wilk %s: W = %.3f, p = %.3f\n", v, s$statistic, s$p.value))
anotar(Pregunta = "P4 Clima (supuesto)", Grupo = v, Variable = v, Factor = "", n = nrow(datos_P4),
Niveles = "", Prueba = "Shapiro-Wilk", Estadistico = sprintf("W = %.3f", s$statistic),
gl = NA, p = signif(s$p.value, 3), Efecto = "")
s$p.value > alpha
})
cat("Valores distintos: temperatura", length(unique(round(datos_P4$T_90d, 2))),
"| precipitación", length(unique(round(datos_P4$P_90d, 1))),
"(muestreos del mismo día comparten clima: hay empates)\n")
correlacion <- function(datos, etiqueta) {
for (x in c("T_90d", "P_90d")) {
principal <- if (normal_P4[["Largo_medio"]] && normal_P4[[x]]) "pearson" else "spearman"
for (m in c("pearson", "spearman")) {
ct <- cor.test(datos$Largo_medio, datos[[x]], method = m, exact = FALSE)
ic <- if (m == "pearson") sprintf("IC 95 %%: %.2f a %.2f", ct$conf.int[1], ct$conf.int[2]) else ""
anotar(Pregunta = etiqueta, Grupo = x, Variable = "Largo_medio", Factor = x,
n = nrow(datos), Niveles = ifelse(m == principal, "Prueba principal", "Complemento"),
Prueba = ifelse(m == "pearson", "Pearson", "Spearman"),
Estadistico = sprintf("%s = %.2f", ifelse(m == "pearson", "r", "ρ"), ct$estimate),
gl = ifelse(m == "pearson", as.character(ct$parameter), NA),
p = signif(ct$p.value, 3), Efecto = ic)
print(ct)
}
}
}
correlacion(datos_P4, "P4 Clima (correlación)")
correlacion(subset(datos_P4, Sitio != "Punta Flor"), "P4 Clima (sensibilidad sin Punta Flor)")
guardar_figura("Fig_I2_correlacion_clima.png", alto_cm = 8, function() {
par(mfrow = c(1, 2), mar = c(4, 4, 1, 1))
simbolo <- c(RNA = 19, AMM = 1, Ostional = 2)[as.character(datos_P4$Area)]
plot(Largo_medio ~ T_90d, data = datos_P4, pch = simbolo,
xlab = "Temperatura media de 90 d (\u00b0C)", ylab = "Largo medio (cm)")
mtext("A", side = 3, adj = 0, line = 0, font = 2)
plot(Largo_medio ~ P_90d, data = datos_P4, pch = simbolo,
xlab = "Precipitaci\u00f3n acumulada de 90 d (mm)", ylab = "Largo medio (cm)")
mtext("B", side = 3, adj = 0, line = 0, font = 2)
legend("topright", legend = c("RNA", "AMM", "Ostional"), pch = c(19, 1, 2), bty = "n")
})
# ---------------------------------------------------------------------
# P5. Prueba exacta de Fisher: categoría PRONAMEC x área (n >= 10)
# ---------------------------------------------------------------------
mue_P5 <- subset(mue, n >= 10)
mue_P5$Categoria <- cut(mue_P5$Largo_medio, c(-Inf, 1, 1.5, 2.5, 4, Inf), right = FALSE,
labels = c("Malo", "Sin categoría", "Regular", "Bueno", "Sobre el ámbito"))
tab_P5 <- table(mue_P5$Area, droplevels(mue_P5$Categoria))
cat("\n================================================================\n")
cat("P5 PRONAMEC | muestreos:", nrow(mue_P5), "\n"); print(tab_P5)
fis <- fisher.test(tab_P5); print(fis)
anotar(Pregunta = "P5 PRONAMEC", Grupo = "Todos", Variable = "Categoría (largo medio)",
Factor = "Area", n = nrow(mue_P5), Niveles = paste(colnames(tab_P5), collapse = ", "),
Prueba = "Exacta de Fisher", Estadistico = "", gl = NA, p = signif(fis$p.value, 3),
Efecto = "", p_Bonferroni = NA)
# ---------------------------------------------------------------------
# Resumen final y exportación
# ---------------------------------------------------------------------
cat("\n\n#################### RESUMEN DE PRUEBAS ####################\n")
print(resumen[, c("Pregunta", "Grupo", "Variable", "n", "Prueba", "Estadistico", "gl",
"p", "p_Bonferroni", "Efecto")], row.names = FALSE)
cat("\n#################### COMPARACIONES POR PARES ####################\n")
print(posthoc, row.names = FALSE)
salida <- list(Resumen_pruebas = resumen, Post_hoc = posthoc,
Datos_P4_muestreos = datos_P4, Tabla_P5 = cbind(Area = rownames(tab_P5),
as.data.frame.matrix(tab_P5)))
if (requireNamespace("writexl", quietly = TRUE)) {
writexl::write_xlsx(salida, file.path(carpeta, "Inferencial_Nerita.xlsx"))
cat("\nResultados guardados en:", file.path(carpeta, "Inferencial_Nerita.xlsx"), "\n")
} else {
for (t in names(salida)) write.csv(salida[[t]], file.path(carpeta, paste0(t, ".csv")),
row.names = FALSE)
cat("\nResultados guardados como CSV en:", carpeta, "\n")
}B4. 04_anova_dos_vias_Burgaditos.R
Análisis no planificados: ANOVA de dos vías, Kruskal-Wallis por muestreo y ancho.
# =====================================================================
# 04_anova_dos_vias.R · Los Burgaditos
# Análisis NO incluidos en el pre-registro v3 (registrados en la bitácora
# como no planificados):
#
# A1 ANOVA de dos vías Área x Año con interacción (pedido de la docente
# para la entrega final). Se usan solo 2025 y 2026, los únicos años
# en que se muestrearon las cuatro áreas: así todas las celdas
# área-año tienen datos y la interacción se puede estimar completa.
# Diseño desbalanceado -> sumas de cuadrados tipo II (car::Anova).
# Supuestos: Shapiro-Wilk de residuos y Levene entre celdas.
# Gráfico de interacción con medias ± IC 95 %.
# A2 Kruskal-Wallis con el largo medio de cada muestreo como unidad
# (n = 34 muestreos con n >= 10), para evaluar la pseudorreplicación.
#
# Correr DESPUÉS de 03_inferencial.R. Datos: BD_Nerita_limpia.xlsx.
# Paquetes: readxl, car (writexl opcional).
# =====================================================================
library(readxl)
library(car)
alpha <- 0.05
ruta_bd <- if (exists("RUTA_BD")) RUTA_BD else file.choose()
bio <- as.data.frame(read_excel(ruta_bd, sheet = "bio_limpia"))
bio$Fecha <- as.Date(bio$Fecha)
etiquetas <- c("RNA Cabo Blanco" = "RNA", "AMM Cabo Blanco" = "AMM",
"RNVS Camaronal" = "Camaronal", "RNVS Ostional" = "Ostional")
bio$Area <- factor(etiquetas[bio$Area], levels = etiquetas)
bio$Muestreo <- paste(bio$Sitio, bio$Fecha)
carpeta <- file.path(dirname(ruta_bd), "inferencial")
dir.create(carpeta, showWarnings = FALSE)
# ---------------------------------------------------------------------
# A1. ANOVA de dos vías (2025-2026)
# ---------------------------------------------------------------------
d <- subset(bio, Anio %in% c(2025, 2026))
d$Anio <- factor(d$Anio)
cat("\nCeldas área x año (n de individuos):\n"); print(table(d$Area, d$Anio))
modelo <- aov(Largo ~ Area * Anio, data = d)
tabla <- Anova(modelo, type = 2)
print(tabla)
sc <- tabla[["Sum Sq"]]; names(sc) <- rownames(tabla)
eta_p <- sc[1:3] / (sc[1:3] + sc["Residuals"]) # η² parcial
print(round(eta_p, 3))
sw <- shapiro.test(residuals(modelo))
lev <- leveneTest(Largo ~ Area * Anio, data = d)
print(sw); print(lev)
# Como las varianzas difieren entre celdas, se repite la tabla con errores
# robustos a heterocedasticidad (HC3) para confirmar las conclusiones.
tabla_hc3 <- Anova(lm(Largo ~ Area * Anio, data = d), type = 2, white.adjust = "hc3")
print(tabla_hc3)
# Medias por celda ± IC 95 %
medias <- aggregate(Largo ~ Area + Anio, data = d, FUN = function(x)
c(n = length(x), media = mean(x), ic = qt(0.975, length(x) - 1) * sd(x) / sqrt(length(x))))
medias <- do.call(data.frame, medias); names(medias)[3:5] <- c("n", "media", "ic95")
print(medias, digits = 3)
png(file.path(carpeta, "Fig_I3_interaccion_area_anio.png"),
width = 16.5 / 2.54, height = 11 / 2.54, units = "in", res = 330, pointsize = 9)
par(mar = c(4, 4.2, 1, 1))
x <- c(1, 2); col_area <- c("black", "grey40", "grey60", "grey20")
pch_area <- c(19, 1, 17, 2); lty_area <- c(1, 2, 3, 4)
plot(NA, xlim = c(0.8, 2.6), ylim = range(medias$media - medias$ic95, medias$media + medias$ic95),
xaxt = "n", xlab = "Año", ylab = "Largo medio de la concha (cm)")
axis(1, at = x, labels = levels(d$Anio))
for (i in seq_along(levels(d$Area))) {
m <- medias[medias$Area == levels(d$Area)[i], ]
xi <- x + (i - 2.5) * 0.03
arrows(xi, m$media - m$ic95, xi, m$media + m$ic95, angle = 90, code = 3, length = 0.03, col = col_area[i])
lines(xi, m$media, lty = lty_area[i], col = col_area[i])
points(xi, m$media, pch = pch_area[i], col = col_area[i], bg = "white")
}
legend("topright", legend = levels(d$Area), pch = pch_area, lty = lty_area, col = col_area, bty = "n")
dev.off()
res_A1 <- data.frame(Termino = rownames(tabla), SC = round(tabla[["Sum Sq"]], 2), gl = tabla$Df,
F = round(tabla[["F value"]], 2), p = signif(tabla[["Pr(>F)"]], 3),
eta2_parcial = c(round(eta_p, 3), NA))
res_sup <- data.frame(Supuesto = c("Shapiro-Wilk residuos", "Levene celdas"),
Estadistico = c(sprintf("W = %.3f", sw$statistic), sprintf("F = %.2f", lev[1, "F value"])),
gl = c(NA, paste(lev$Df, collapse = " y ")),
p = signif(c(sw$p.value, lev[1, "Pr(>F)"]), 3))
res_hc3 <- data.frame(Termino = rownames(tabla_hc3), gl = tabla_hc3$Df,
F = round(tabla_hc3$F, 2), p = signif(tabla_hc3[["Pr(>F)"]], 3))
# ---------------------------------------------------------------------
# A2. Kruskal-Wallis con el muestreo como unidad
# ---------------------------------------------------------------------
mue <- aggregate(Largo ~ Muestreo + Area, data = bio, FUN = function(x) c(n = length(x), media = mean(x)))
mue <- do.call(data.frame, mue); names(mue)[3:4] <- c("n", "Largo_medio")
mue <- subset(mue, n >= 10)
kw <- kruskal.test(Largo_medio ~ Area, data = mue)
print(kw)
res_A2 <- data.frame(Prueba = "Kruskal-Wallis (muestreo como unidad)", n = nrow(mue),
H = round(kw$statistic, 2), gl = kw$parameter, p = signif(kw$p.value, 3))
if (requireNamespace("writexl", quietly = TRUE))
writexl::write_xlsx(list(A1_ANOVA_dos_vias = res_A1, A1_supuestos = res_sup, A1_robusto_HC3 = res_hc3,
A1_medias = medias, A2_KW_muestreo = res_A2),
file.path(carpeta, "No_planificados_Nerita.xlsx"))
# ---------------------------------------------------------------------
# A3. Ancho (complemento): correlación con el clima y ANOVA de dos vías
# (solicitado por el equipo el 11/11/2026, misma regla que el largo)
# ---------------------------------------------------------------------
p4 <- as.data.frame(read_excel(file.path(carpeta, "Inferencial_Nerita.xlsx"), sheet = "Datos_P4_muestreos"))
p4$Fecha <- as.Date(p4$Fecha)
an <- aggregate(Ancho ~ Sitio + Fecha, data = bio, FUN = function(x) c(nA = length(x), Ancho_medio = mean(x)))
an <- do.call(data.frame, an); names(an)[3:4] <- c("n_ancho", "Ancho_medio")
p4a <- merge(p4, an, by = c("Sitio", "Fecha"))
p4a <- subset(p4a, n_ancho >= 10)
cat("\nA3 clima-ancho | muestreos:", nrow(p4a), "\n")
print(sapply(c("Ancho_medio", "T_90d", "P_90d"), function(v) shapiro.test(p4a[[v]])$p.value))
cor_res <- do.call(rbind, lapply(c("T_90d", "P_90d"), function(v) {
s <- suppressWarnings(cor.test(p4a$Ancho_medio, p4a[[v]], method = "spearman"))
p <- cor.test(p4a$Ancho_medio, p4a[[v]], method = "pearson")
data.frame(Variable = v, n = nrow(p4a), rho = round(s$estimate, 2), p_spearman = signif(s$p.value, 3),
r = round(p$estimate, 2), IC95 = paste(round(p$conf.int, 2), collapse = " a "),
p_pearson = signif(p$p.value, 3))
}))
print(cor_res)
da <- subset(d, !is.na(Ancho))
mod_a <- aov(Ancho ~ Area * Anio, data = da)
tab_a <- Anova(mod_a, type = 2); print(tab_a)
sca <- tab_a[["Sum Sq"]]; eta_a <- sca[1:3] / (sca[1:3] + sca[4])
sw_a <- shapiro.test(residuals(mod_a)); lev_a <- leveneTest(Ancho ~ Area * Anio, data = da)
hc3_a <- Anova(lm(Ancho ~ Area * Anio, data = da), type = 2, white.adjust = "hc3")
print(sw_a); print(lev_a); print(hc3_a)
res_A3 <- data.frame(Termino = rownames(tab_a), gl = tab_a$Df, F = round(tab_a[["F value"]], 2),
p = signif(tab_a[["Pr(>F)"]], 3), eta2_parcial = c(round(eta_a, 3), NA),
F_HC3 = round(hc3_a$F, 2), p_HC3 = signif(hc3_a[["Pr(>F)"]], 3))
if (requireNamespace("writexl", quietly = TRUE))
writexl::write_xlsx(list(A1_ANOVA_dos_vias = res_A1, A1_supuestos = res_sup, A1_robusto_HC3 = res_hc3,
A1_medias = medias, A2_KW_muestreo = res_A2,
A3_ancho_clima = cor_res, A3_ancho_ANOVA = res_A3),
file.path(carpeta, "No_planificados_Nerita.xlsx"))Anexo C. Información de la sesión de R
sessionInfo()#> R version 4.6.1 (2026-06-24 ucrt)
#> Platform: x86_64-w64-mingw32/x64
#> Running under: Windows 11 x64 (build 26200)
#>
#> Matrix products: default
#> LAPACK version 3.12.1
#>
#> locale:
#> [1] LC_COLLATE=Spanish_Costa Rica.utf8 LC_CTYPE=Spanish_Costa Rica.utf8
#> [3] LC_MONETARY=Spanish_Costa Rica.utf8 LC_NUMERIC=C
#> [5] LC_TIME=Spanish_Costa Rica.utf8
#>
#> time zone: America/Costa_Rica
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] knitr_1.52 car_3.1-5 carData_3.0-6 readxl_1.6.0
#>
#> loaded via a namespace (and not attached):
#> [1] vctrs_0.7.3 cli_3.6.6 rlang_1.3.0 xfun_0.62
#> [5] Formula_1.2-6 otel_0.2.0 jsonlite_2.0.0 glue_1.8.1
#> [9] htmltools_0.5.9 sass_0.4.10 writexl_2.0.1 rmarkdown_2.32
#> [13] cellranger_1.1.0 evaluate_1.0.5 jquerylib_0.1.4 tibble_3.3.1
#> [17] abind_1.4-8 fastmap_1.2.0 yaml_2.3.12 lifecycle_1.0.5
#> [21] compiler_4.6.1 pkgconfig_2.0.3 rstudioapi_0.19.0 digest_0.6.39
#> [25] R6_2.6.1 pillar_1.11.1 moments_0.14.1 magrittr_2.0.5
#> [29] bslib_0.12.0 tools_4.6.1 cachem_1.1.0