Taller Final - Bioestadistica
Sebastian Hidalgo Jimenez y Maria Camila Chaves
2026-05-29
Abstracto
Se hizo una replicación de la metodología estadística de “Habitat selection by a threathened desert amphibian” (Hinderer et al., 2020). Con esta se obtuvo un modelo de regresión logística para identificar la preferencia de la rana leopardo Chiricaua (Lithobates chiricahuensis) por microhábitats en Nuevo Mexico, Estados Unidos. A partir del modelo se obtuvieron Odd ratios (OR) para las variables categóricas para poder observar las preferencias para cada tipo de las variables. Los resultados fueron consistentes con el modelo original donde para esta especie la presencia de cobertura arbórea, presencia de agua, cantidad de cobertura baja y su tipo, y el sustrato de lodo son preferidos, esto es consistente con las necesidades de un anfibio y la biología de esta especie en una zona tan árida como la que habita, por lo cual se debe buscar para su conservación que estas características estén presentes en su hábitat. Aunque hay ciertas diferencias en el modelo y OR debido a una falta de una descripción más explícita de los métodos estadísticos empleados que llevaron a diferencias en la metodología usada (otras fueron por decisión) y cambiando la interpretación de categorías dentro de las variables categóricas especialmente.
Abstract
A replication of the statistical methodology from ‘Habitat selection by a threatened desert amphibian’ (Hinderer et al., 2020) was conducted. Through this replication, a logistic regression model was obtained to identify the microhabitat preferences of the Chiricahua leopard frog (Lithobates chiricahuensis) in New Mexico, United States. Odds ratios (OR) for the categorical variables were derived from the model to observe preferences across each category. The results were consistent with the original model, demonstrating that tree cover, the presence of water, the amount and type of low cover, and mud substrate are preferred by this species. This aligns with the ecological requirements of an amphibian and the biology of this species in such an arid environment, highlighting that these characteristics must be maintained within their habitat for conservation purposes. However, certain differences in the model and the ORs emerged due to a lack of a more explicit description of the statistical methods originally employed, which led to variations in the methodology used (some of which were intentional decisions) and altered the interpretation of specific categories within the categorical variables.
Introducción
La degradación de hábitat y la fragmentación son una fuente de riesgo importante para las poblaciones de anfibios alrededor del mundo, por lo que frente a los cambios climáticos que con el pasar de los anos empeoran la situación para este grupo de vertebrados que debido a su piel son muy vulnerables a cambios en la temperatura que pueden llevar a la desecación y además son altamente dependientes a fuentes de agua. De este modo, Hinderer et al. (2020) evaluó los microhábitats de la rana leopardo Chiricaua (Lithobates chiricahuensis), una especie nativa del el suroeste de Nuevo Mexico y de la zona central al sureste de Arizona, siendo una especie altamente acuática que suele encontrarse en ambientes naturales como piscinas rocosas, estanques y manantiales aunque también se ha adaptado a hábitats antrópicos como tanques de agua, debido a las condiciones generalmente áridas de su área y la reducción de sus zonas acuáticas es una especie vulnerable. Teniendo en cuenta estas características Hinderer et al. (2019) utilizo el rancho privado Ladder Ranch, ubicado en Nuevo Mexico, Estados Unidos donde muestrearon 44 individuos de esta especie durante el verano de 2014 usando radiotransmisores, estas se localizaron diariamente por 8 semanas donde se tomaron características de donde las encontraban y de un hábitat aleatorio ubicado a 5m de los individuos, esto con el fin de poder identificar las características ambientales que la especie prefiere durante la temporada de lluvias y de este modo mejorar las decisiones en cuanto su conservación.
Materiales y metodos
Se utilizó la base de datos original del estudio disponible en Dryad, la cual es un CSV que contiene las siguientes columnas: DATE, TIME, FROGID, USE, UTME, UTMN, WATER, DIST, SUBSTRATE, FCOVERPCT, FCOVERTYP, OCOVERTYP, OCOVERPCT y PHOTO (las coordenadas y la foto fueron retiradas por los autores por seguridad de la población y el rancho).
Con la base de datos se trabajó en R mediante RStudio para replicar el análisis. Se utilizó un modelo de Cox estratificado por pares, que funciona como análogo a la regresión logística condicional y permite comparar directamente el microhábitat utilizado por cada rana con su punto aleatorio disponible a 5 m. Para estructurar el emparejamiento se creó la variable STRATA, que asigna un índice único a cada par rana con punto aleatorio, garantizando que la comparación sea siempre dentro del par y no entre individuos.
Previo a la creación del modelo se recodificaron las variables. Los tipos de sustrato y de cobertura baja con pocas observaciones se condensaron en una categoría “OTHER”; para la cobertura alta, todas las especies arbóreas se agruparon en “TREE SPP.” excepto Juniper y Willow, que se mantuvieron como categorías independientes por su frecuencia. La variable WATER se estandarizó a binaria (1 = presencia, 0 = ausencia) corrigiendo inconsistencias de mayúsculas en el CSV original. Adicionalmente se convirtio OCOVERPCT en binaria (1 si > 10%, 0 si ≤ 10%), pero esta variable no resultó significativa y no se incluyó en el modelo final.
Las variables del modelo final fueron WATER, FCOVERPCT, FCOVERPCT² (creada para capturar la relación no lineal con máximo alrededor del 70% de cobertura), LCOVERTYP (recodificada a partir de FCOVERTYP), OCOVERTYP y SUBSTRATE. Al igual que en el artículo, se evaluó la importancia de cada variable mediante pruebas de razón de verosimilitud (LRT), comparando el modelo completo contra versiones reducidas sin cada variable; un valor p < 0.05 indica que la variable contribuye significativamente al modelo.
A diferencia del artículo original que usa coxph, en la replicación se utilizó la función coxme, que permite incluir un efecto aleatorio por individuo (1|FROGID). Este efecto corrige la pseudorreplicación generada porque cada rana aporta múltiples observaciones, y permite que la tendencia de selección varíe entre individuos, evitando que se subestime el error estándar de los coeficientes. El efecto aleatorio se incluyó tanto en el modelo completo como en cada modelo reducido para el LRT.
Finalmente, se obtuvieron los Odds Ratios (OR) exponenciando los coeficientes del modelo (exp(coef)) y sus intervalos de confianza del 95% con confint(). Con estos se construyeron gráficas que muestran la relación entre cada variable y la probabilidad de selección de hábitat, incluyendo los IC al 95% para facilitar la interpretación de la magnitud e incertidumbre de cada efecto.
Resultados
A partir de la metodología anteriormente descrita se obtuvieron los siguientes resultados:
Modelo
A partir del resumen del modelo se observan las diferentes variables y categorías en ellas junto a su nivel de significancia. En estas se observa que las que presentan un mayor nivel de significancia son presencia de agua (WATER), el porcentaje de cobertura baja y su versión cuadrática (FCOVERPCT y FCOVERPCT^2), el substrato especialmente tierra, roca y otros (SUBSTRATE_CATSOIL, SUBSTRATECAT_ROCK, SUBSTRATECATE_OTHERS). A las anteriores variables les siguen (p < 0.01) son el tipo de cobertura superior de Willow y otras especies de árboles (SUBSTRATE_CATWILLOW y SUBSTRATE_TREESPP). Continua el sustrato de arena (SUBSTRATE_CATSAND) (p<0.05) y la cobertura baja leñosa y vegetación acuática LCOVERTYPWOODY.
Con el resumen del modelo podemos observar cómo tenemos 1036 eventos de encuentro con los individuos mientras que se tienen 2072 observaciones totales al sumarle a cada una el punto aleatorio. Posteriormente, se comprueba que el uso del efecto aleatorio por el individuo fue significativo al obtener un chi cuadrado alto, un valor p de 0 o cercano indica que la probabilidad de fuera al azar es casi nula y ver que el modelo aleatorio presenta valores más bajos de AIC y BIC que miden el ajuste del modelo aumentando de valor si son muy complejos.
Tabla 1 y 2
##
## ===== Odds Ratios del modelo final =====## OR 2.5 % 97.5 %
## WATER 2.944 1.929 4.493
## FCOVERPCT 1.087 1.064 1.111
## FCOVERPCT2 0.999 0.999 1.000
## LCOVERTYPANNUALS 1.914 0.361 10.158
## LCOVERTYPAQVEG 4.081 0.745 22.367
## LCOVERTYPOPENWATER 1.945 0.361 10.496
## LCOVERTYPOTHER 3.160 0.589 16.959
## LCOVERTYPROCK 3.654 0.712 18.745
## LCOVERTYPWOODY 7.052 1.130 44.008
## OCOVERTYP_CATJUNIPER 1.539 0.832 2.847
## OCOVERTYP_CATOTHER 1.735 0.731 4.120
## OCOVERTYP_CATTREESPP 1.919 1.234 2.984
## OCOVERTYP_CATWILLOW 2.126 1.328 3.405
## SUBSTRATE_CATOTHER 0.118 0.037 0.377
## SUBSTRATE_CATROCK 0.161 0.079 0.331
## SUBSTRATE_CATSAND 0.437 0.224 0.851
## SUBSTRATE_CATSOIL 0.257 0.131 0.504## # A tibble: 6 × 5
## Variable df chi2 p_valor Significancia
## <chr> <dbl> <dbl> <chr> <chr>
## 1 WATER 1 27 < 0.001 ***
## 2 LCOVERPCT 2 84.7 < 0.001 ***
## 3 LCOVERPCT^2 1 49.7 < 0.001 ***
## 4 LCOVERTYPE 6 32.1 < 0.001 ***
## 5 SUBSTRATE 4 49.4 < 0.001 ***
## 6 OCOVERTYPE 4 16.4 0.0025 **Siguiendo a la tabla 1, esta muestra el uso de la prueba de mayor verosimilitud (LRT) para cada variable donde se obtuvo que todas las variables utilizadas en el modelo son muy significativas, pero el tipo de cobertura superior es un poco menos significativa que las demás (p<0.01). Esto muestra que las ranas están prefiriendo ubicaciones según la presencia de agua, el porcentaje de cobertura baja (cuadrática) y su tipo, el tipo de cobertura alta/superior y el tipo de sustrato del suelo.
Luego, se obtuvieron los odd ratio (OR) para cada una de las variables. estos expresan la diferencia de probabilidad en escoger un sitio que tenga la variable presente contra un valor base cuando la variable era categórica, mientras que para variables continuas refleja la selección donde el nivel de la variable aumenta respecto a cero. A continuación, se analiza cada variable por separado.
La presencia de agua presento 2.944 (95% IC = 1.929 -4.493) más probabilidades de elegir una ubicación con agua que una sin ella.
En cuanto a la cantidad de cobertura baja se detectó que la prefieren un 1.087 mas (1.064-1.111), la relación cuadrática por el IC pareciera no ser significativa pero el hecho de estar en 1 probablemente se deba a aproximación de R. Se observa preferencia por la cobertura leñosa, roca y vegetación acuática (7.052, 3.654 y 4.081 respectivamente)., es importante destacar que los IC amplios y que pasan por 1 se deben a la separación por estratos y no una falta de efecto que fue confirmada con el LRT global de la Tabla 2, sin embargo, también indica que excepto por WOODY todas las otras no presentan significancia estadística por si solas.
Continuando con el tipo de cobertura superior o alta se observa una preferencia por los árboles Willow (sauces) y especies arbóreas en general con 2.126 y 1.919 respectivo de mayor probabilidad de seleccionar este tipo de cobertura.
En el caso del sustrato, se observa como los valores de roca, arena y tierra son bastante bajos (0.161, 0.437 y 0.257 respectivamente) lo que indica que el barro es preferido en mayor medida al ser el termino de referencia con OR = 1.Figura 1
# Figura 1: OR ~ % de cobertura baja
# Tiene Relación cuadrática. OR relativo a 0% de cobertura (referencia = 1).
b1 <- coef(modelo_final)["FCOVERPCT"]
b2 <- coef(modelo_final)["FCOVERPCT2"]
V <- vcov(modelo_final)
i1 <- which(names(coef(modelo_final)) == "FCOVERPCT")
i2 <- which(names(coef(modelo_final)) == "FCOVERPCT2")
cover_seq <- seq(0, 100, by = 1)
log_OR <- b1 * cover_seq + b2 * cover_seq^2
# Varianza del log-OR por método delta
var_logOR <- cover_seq^2 * V[i1,i1] + cover_seq^4 * V[i2,i2] +
2 * cover_seq^3 * V[i1,i2]
se_logOR <- sqrt(pmax(var_logOR, 0))
fig2_df <- data.frame(
cover = cover_seq,
OR = exp(log_OR),
OR_low = exp(log_OR - 1.96 * se_logOR),
OR_high = exp(log_OR + 1.96 * se_logOR)
)
max_cover <- cover_seq[which.max(fig2_df$OR)]
cat("\nMáximo de OR en % de cobertura:", max_cover,
"| Artículo: ~70%\n")##
## Máximo de OR en % de cobertura: 70 | Artículo: ~70%fig2 <- ggplot(fig2_df, aes(x = cover, y = OR)) +
geom_ribbon(aes(ymin = OR_low, ymax = OR_high),
fill = "grey80", color = "grey60", linetype = "dashed") +
geom_line(linewidth = 0.9) +
geom_vline(xintercept = max_cover, linetype = "dashed") +
labs(x = "% Low-Lying Cover",
y = "Odds Ratio",
caption = paste0("Selection maximized at ", max_cover, "% cover")) +
scale_x_continuous(breaks = seq(0, 100, 20)) +
theme_classic(base_size = 13)
print(fig2)
La grafica anterior (Figura 1) muestra la relación de la cobertura baja y los odds ratio de selección del hábitat. En esta se observa cómo se aumenta la probabilidad de seleccionar el hábitat hasta el punto máximo donde la baja llega al 70%, a partir de ese punto las probabilidades comienzan a disminuir, sin embargo debido a la amplitud del IC de 95% esto no puede asegurarse con certeza ya que el pico podría ser más alto (aumentando el OR a más de 30%) o más bajo.
Figura 1
# Figura 2: OR por tipo de cobertura y sustrato
# Extraer coeficientes del modelo como data.frame con IC
coef_df <- data.frame(
term = names(coef(modelo_final)),
coef = coef(modelo_final),
se = sqrt(diag(vcov(modelo_final)))
) %>%
mutate(OR = exp(coef),
OR_low = exp(coef - 1.96 * se),
OR_high = exp(coef + 1.96 * se))
# Panel a: tipo de cobertura baja
panel_a_df <- coef_df %>%
filter(grepl("^LCOVERTYP", term)) %>%
mutate(
label = recode(sub("LCOVERTYP", "", term),
ANNUALS = "Annuals",
AQVEG = "Aq Veg",
OPENWATER = "Open Water",
ROCK = "Rock",
WOODY = "Woody Cover",
OTHER = "Other"),
label = fct_reorder(label, OR, .desc = TRUE)
)
fig3a <- ggplot(panel_a_df, aes(x = label, y = OR)) +
geom_hline(yintercept = 1, linetype = "dashed", color = "grey50") +
geom_errorbar(aes(ymin = OR_low, ymax = OR_high), width = 0.2) +
geom_point(size = 3) +
annotate("text", x = 0.7, y = 1.15, label = "(No Cover)",
size = 3.2, color = "grey40", hjust = 0) +
labs(x = "Low-Lying Cover Type", y = "Odds Ratio", tag = "(a)") +
theme_classic(base_size = 12)
# Panel (b): tipo de cobertura superior
panel_b_df <- coef_df %>%
filter(grepl("^OCOVERTYP_CAT", term)) %>%
mutate(
label = recode(sub("OCOVERTYP_CAT", "", term),
WILLOW = "Willow",
JUNIPER = "Juniper",
TREESPP = "Tree spp.",
OTHER = "Other"),
label = fct_reorder(label, OR, .desc = TRUE)
)
fig3b <- ggplot(panel_b_df, aes(x = label, y = OR)) +
geom_hline(yintercept = 1, linetype = "dashed", color = "grey50") +
geom_errorbar(aes(ymin = OR_low, ymax = OR_high), width = 0.2) +
geom_point(size = 3) +
annotate("text", x = 0.7, y = 1.04, label = "(No Cover)",
size = 3.2, color = "grey40", hjust = 0) +
labs(x = "Overstory Cover Type", y = "Odds Ratio", tag = "(b)") +
theme_classic(base_size = 12)
# Panel c: tipo de sustrato
panel_c_df <- coef_df %>%
filter(grepl("^SUBSTRATE_CAT", term)) %>%
mutate(
label = recode(sub("SUBSTRATE_CAT", "", term),
SAND = "Sand",
SOIL = "Soil",
ROCK = "Rock",
OTHER = "Other"),
label = fct_reorder(label, OR, .desc = TRUE)
)
fig3c <- ggplot(panel_c_df, aes(x = label, y = OR)) +
geom_hline(yintercept = 1, linetype = "dashed", color = "grey50") +
geom_errorbar(aes(ymin = OR_low, ymax = OR_high), width = 0.2) +
geom_point(size = 3) +
annotate("text", x = 0.7, y = 1.05, label = "(Mud Substrate)",
size = 3.2, color = "grey40", hjust = 0) +
labs(x = "Substrate Type", y = "Odds Ratio", tag = "(c)") +
theme_classic(base_size = 12)
# Figura combinada y figura c por separado
fig3_combined <- fig3a / fig3b
print(fig3_combined)
Estas graficas muestran como interactúa cada tipo dentro de las variables categóricas con el OR. En primer lugar, la figura 1a muestra como el coeficiente estadístico de todos los tipos de cobertura baja se encuentra sobre la no cobertura, donde la vegetación acuática, cobertura leñosa y roca siendo más preferidas entre el resto, con la categoría OTHER siendo cercana. También es importante notar que el IC para todas excepto leñosos se encuentran por debajo del OR sin cobertura por lo que no se puede afirmar con certeza el efecto de preferencia; a esto se le suma que sus IC son bastante amplios por lo que el valor exacto de OR no puede afirmarse con certeza y puede ser tanto más alto como más bajo, donde los leñosos también se ven afectados y por lo tanto se ve una preferencia, pero no se puede afirmar con seguridad su magnitud.
La grafica 1b presenta el OR relacionado con el tipo de cobertura superior. Tanto los sauces como las otras especies de árboles presentan los OR más altos e IC no tan amplios por lo que se puede tener cierta confianza moderada de la magnitud de la preferencia. Caso contrario ocurre con Juniper y otros donde el IC queda por debajo de la línea sin cobertura, especialmente no se puede interpretar mucho de la categoría otros al tener también alto IC. Por último, es importante destacar que se ve preferencia por cualquier tipo de cobertura por sobre a ninguna, aunque en juncos y otros no se puede afirmar con certeza.
Con la última grafica de evaluación de OR y las variables esta la 1c evaluando el tipo de sustrato. En este caso se observa como todas las variables están significativamente por debajo de la línea de referencia del barro, lo que comprueba el análisis hecho antes donde el barro tiene una selección bastante más alta por la especie de rana estudiada.
Discusion
De acuerdo con los resultados previamente descritos el modelo logro describir las preferencias de Lithobates chiricahuensis por características de su hábitat. Donde las características seleccionadas suelen ser claves para evitar desecación al igual que en el artículo como lo son la cobertura baja y su tipo, el tipo de cobertura alta, la presencia de agua y sustrato de barro, que son claves en climas tan áridos como los de Nuevo Mexico. El hecho de que el agua abierta y las plantas tuvieran tan poca selección se relaciona con el hecho de que no tienen cobertura del sol y por lo tanto no ayudan mucho a termo regularse y mantenerse húmedas en esas temperaturas. Adicionalmente se vio una clara preferencia por especies arbóreas probablemente por la capacidad de sombra y los microhábitats húmedos que se generan en la hojarasca baja o entre las raíces (parte de Woody altamente seleccionado), esto a su vez explica porque la categoría otros (OTHER) en los tipos de cobertura alta fue cercana a los de OR más altos y es debido a que la hojarasca estaba incluida en ellos la cual es un hábitat bastante común para ranas en las Americas. Finalmente, en términos de conservación se concuerda con las recomendaciones del articulo original, donde para su conservación se debe impulsar que se mantengan restos leñosos, vegetación acuática, rocas, sustrato de lodo y cobertura arbórea para ayudar a que la población pueda tener un hábitat más amplio a los estanques de reproducción, estas características de hábitat normalmente se generan por si solas en la temporada de lluvias y por lo tanto se recomienda no limpiar estos restos que son en verdad el hábitat de la especie y tener más rango de movimiento en sus hábitat evita el aislamiento y mejora el flujo genético lo cual es clave en una especie amenazada como esta.
Se tienen limitaciones importantes en la replicación de estos datos principalmente debido a la descripción de la metodología del artículo. En este se omiten partes importantes del procesamiento de la base de datos y análisis estadísticos donde no se especifica como se agrupan variables categóricas (excepto el tipo de cobertura superior) lo cual dificulta la replicación y explica las diferencias con el articulo original en los valores tanto de las variables como de las categorías en estas; como en el caso de los tipos de cobertura baja, el artículo original tiene un comportamiento altamente diferente que en el actual lo cual puede deberse a la inclusión de la hojarasca en él. Adicionalmente, no se especifica que función exacta se usa para el modelado de la regresión logística condicional, que a pesar de que se utilizó coxme pudo ser coxph o clogit lo cual pudo llevar a las discrepancias encontradas, pero tambien se asume por temas de tiempo no se pudo probar todos los modelos para encontrar el adecuado al igual que tampoco se hizo la selección de variables a partir de un modelo con todas las variables de la base de datos, esto también pudo ser una fuente de variación en el modelo actual, que pudo ayudar a ajustar mas los coeficientes a los del articulo al tener más iteraciones.
A continuacion se presentan las principales diferencias con el articulo original:
Tabla 3
# Cada fila es un OR comparado entre el articulo original
# y la replicacion hecha en el Rmd del taller.
# Columnas:
# Variable -> variable o categoria comparada
# Articulo_valor -> OR reportado en Hinderer et al. 2021 (IC 95%)
# Rmd_valor -> OR obtenido en la replicacion (IC 95%)
# Estado -> Coincide / Parcial / Diferencia
# Causa_probable -> explicacion de la discrepancia
tabla <- tibble(
Variable = c(
# Presencia de agua
"WATER",
# Cobertura baja: efecto continuo y pico de la curva cuadratica
"Pico OR cobertura baja (%)", # maximo de la curva cuadratica
"FCOVERPCT (por unidad)", # efecto lineal por cada 1% de cobertura
# Categorias de tipo de cobertura baja (referencia = NONE)
"WOODY vs NONE",
"ROCK vs NONE",
"AQVEG vs NONE",
# Categorias de tipo de cobertura alta (referencia = NONE)
"Willow vs NONE",
"Tree spp. vs NONE",
"Juniper vs NONE",
# Categorias de sustrato (referencia = MUD)
"SAND vs MUD",
"SOIL vs MUD",
"ROCK vs MUD"
),
# Valores del articulo original (Hinderer et al., 2021)
# Estimado puntual con IC 95% entre parentesis
Articulo_valor = c(
"2.92 (2.78 - 3.07)",
"~70%",
"No reportado", # el articulo no desglosa el coeficiente lineal
"6.03 (2.83 - 12.85)",
"4.14 (2.01 - 8.53)",
"3.98 (1.83 - 8.65)",
"2.22 (2.09 - 2.36)",
"1.99 (1.90 - 2.09)",
"1.54 (1.38 - 1.71)",
"0.43 (0.38 - 0.49)",
"0.26 (0.23 - 0.30)",
"0.16 (0.14 - 0.18)"
),
# Valores obtenidos en la replicacion (Taller Rmd)
# El modelo usa coxme en lugar de coxph, incorporando el efecto
# aleatorio (1|FROGID) de forma explicita, lo que produce IC
# mas amplios en casi todos los OR
Rmd_valor = c(
"2.944 (1.929 - 4.493)", # mismo punto estimado, IC mas amplio
"~68 - 69%", # pico matematico exacto; articulo redondea a 70%
"1.087 (1.064 - 1.111)", # disponible en coxme, no reportado en articulo
"7.052 (1.130 - 44.008)", # IC muy amplio por efecto aleatorio
"3.654 (0.712 - 18.745)", # IC cruza 1; pierde significancia individual
"4.081 (0.745 - 22.367)", # punto estimado similar; IC mucho mas amplio
"2.126 (1.328 - 3.405)", # mismo punto estimado; IC moderadamente mas amplio
"1.919 (1.234 - 2.984)", # misma direccion; IC moderadamente mas amplio
"IC cruza 1 (NS)", # pierde significancia individual en la replicacion
"0.437 (IC similar)", # practicamente identico al articulo
"0.257 (IC similar)", # practicamente identico al articulo
"0.168 (0.082 - 0.344)" # diferencia minima; misma direccion y magnitud
),
# Estado de la comparacion:
# Coincide -> diferencia minima o nula
# Parcial -> misma direccion, magnitud o IC algo distintos
# Diferencia -> discrepancia relevante en magnitud o significancia
Estado = c(
"Diferencia", # WATER: IC mucho mas amplio
"Coincide", # pico cuadratico: practicamente igual (~1-2% de diferencia)
"Diferencia", # no comparable directamente (no reportado en articulo)
"Diferencia", # WOODY: IC enormemente mas amplio
"Diferencia", # ROCK: IC cruza 1
"Diferencia", # AQVEG: punto estimado similar pero IC muy amplio
"Parcial", # Willow: mismo punto estimado, IC algo mas amplio
"Parcial", # Tree spp.: misma direccion, IC algo mas amplio
"Diferencia", # Juniper: pierde significancia individual
"Coincide", # SAND: practicamente identico
"Coincide", # SOIL: practicamente identico
"Parcial" # ROCK: diferencia minima, misma direccion
),
# Explicacion de cada discrepancia
Causa_probable = c(
"IC mucho mas amplio; efecto aleatorio (1|FROGID) aumenta la incertidumbre de los estimados",
"Diferencia de ~1-2%; articulo redondea a 70%",
"No reportado en articulo; coeficiente disponible en coxme",
"Efecto aleatorio + estratificacion por pares genera IC muy amplios en categorias de baja frecuencia",
"IC cruza 1; efecto aleatorio redistribuye varianza reduciendo precision individual",
"Punto estimado similar; IC amplio por baja frecuencia de la categoria y efecto aleatorio",
"Mismo punto estimado; coxme produce IC mas amplios al modelar variabilidad entre ranas",
"Misma direccion e interpretacion; IC ligeramente mas amplio por efecto aleatorio",
"Pierde significancia individual; efecto aleatorio absorbe parte de la varianza de Juniper",
"Sustrato muy frecuente en los datos; robusto al cambio de modelo",
"Sustrato muy frecuente en los datos; robusto al cambio de modelo",
"Diferencia minima (~5%); misma direccion y conclusion biologica"
)
)
#Generar tabla--
kable(
tabla,
col.names = c(
"Variable / Categoria",
"Articulo — OR (IC 95%)",
"Rmd — OR (IC 95%)",
"Estado",
"Causa probable"
),
caption = "Comparacion de Odds Ratios: Hinderer et al. (2021) vs replicacion Rmd",
align = c("l", "c", "c", "c", "l")
)| Variable / Categoria | Articulo — OR (IC 95%) | Rmd — OR (IC 95%) | Estado | Causa probable |
|---|---|---|---|---|
| WATER | 2.92 (2.78 - 3.07) | 2.944 (1.929 - 4.493) | Diferencia | IC mucho mas amplio; efecto aleatorio (1|FROGID) aumenta la incertidumbre de los estimados |
| Pico OR cobertura baja (%) | ~70% | ~68 - 69% | Coincide | Diferencia de ~1-2%; articulo redondea a 70% |
| FCOVERPCT (por unidad) | No reportado | 1.087 (1.064 - 1.111) | Diferencia | No reportado en articulo; coeficiente disponible en coxme |
| WOODY vs NONE | 6.03 (2.83 - 12.85) | 7.052 (1.130 - 44.008) | Diferencia | Efecto aleatorio + estratificacion por pares genera IC muy amplios en categorias de baja frecuencia |
| ROCK vs NONE | 4.14 (2.01 - 8.53) | 3.654 (0.712 - 18.745) | Diferencia | IC cruza 1; efecto aleatorio redistribuye varianza reduciendo precision individual |
| AQVEG vs NONE | 3.98 (1.83 - 8.65) | 4.081 (0.745 - 22.367) | Diferencia | Punto estimado similar; IC amplio por baja frecuencia de la categoria y efecto aleatorio |
| Willow vs NONE | 2.22 (2.09 - 2.36) | 2.126 (1.328 - 3.405) | Parcial | Mismo punto estimado; coxme produce IC mas amplios al modelar variabilidad entre ranas |
| Tree spp. vs NONE | 1.99 (1.90 - 2.09) | 1.919 (1.234 - 2.984) | Parcial | Misma direccion e interpretacion; IC ligeramente mas amplio por efecto aleatorio |
| Juniper vs NONE | 1.54 (1.38 - 1.71) | IC cruza 1 (NS) | Diferencia | Pierde significancia individual; efecto aleatorio absorbe parte de la varianza de Juniper |
| SAND vs MUD | 0.43 (0.38 - 0.49) | 0.437 (IC similar) | Coincide | Sustrato muy frecuente en los datos; robusto al cambio de modelo |
| SOIL vs MUD | 0.26 (0.23 - 0.30) | 0.257 (IC similar) | Coincide | Sustrato muy frecuente en los datos; robusto al cambio de modelo |
| ROCK vs MUD | 0.16 (0.14 - 0.18) | 0.168 (0.082 - 0.344) | Parcial | Diferencia minima (~5%); misma direccion y conclusion biologica |
Analizando primero la prueba LRT se observa como el agua, tipo de cobertura superior, porcentaje de cobertura baja (cuadrática) y sustrato fueron los más cercanos, mientras que el tipo de cobertura baja (32.1 contra 61.11) y el porcentaje de esta misma (84.7 contra 67.78) si son bastante distintos, esto principalmente puede deberse a lo que se menciono anteriormente de la modificacion de las categorias de la variable, la estratificación por pares por esta causa y el efecto aleatorio de coxme captura varianza individual afectando el resultado de la chi-cuadrado, también se observa como el nivel de significancia de el tipo de cobertura superior se redujo ligeramente probablemente por estas mismas causas. Continuando con los OR como se observa en la tabla 3 las mayores diferencias se ven en los IC que esto puede deberse a la incertidumbre que revela el efecto aleatorio, mientras que en las categorías individuales de tipo de cobertura vegetal se observa como ROCK no presenta la significancia del artículo y su OR disminuye significativamente y OTHER tiene un OR muy alto e intervalo muy amplio, esto nuevamente a causa de la agrupación de esta variable y la separación de estratos del efecto aleatorio. Por ultimo es importante destacar como JUNIPER en el tipo de cobertura superior perdió significancia individual, lo cual afecta la interpretación comparado con el articulo y puede deberse a la ampliación de su IC debido a que parte de su varianza la toma el efecto aleatorio y sus pocas repeticiones ampliaron el rango. Aunque es importante aclarar que estos intervalos amplios muestran la variacion de una manera mas clara y son mas “reales” al reflejar la variabilidad presente en los datos del articulos de forma mas explicita
Codigo reproducible
library(survival)
library(ggplot2)
library(dplyr)
library(forcats)
library(patchwork)
library(coxme)
library(knitr)
# 1. Cargar datos
datos_raw <- read.csv("C:/Users/sebas/Documents/5to semestre/Bioestadistica/ranas/2014_habitat_data_NO_UTM.csv",stringsAsFactors = FALSE)
# 2. Preparación y recodificación de variables
datos <- datos_raw %>%
mutate(
# Variable respuesta
# USE = 1: localización donde se encontró la rana (seleccionada)
# USE = 0: localización aleatoria a 5 m (disponible, no ocupada)
USE = as.integer(USE),
# WATER: presencia de agua en el punto
# El artículo registra "y"/"Y" para presencia, "n" para ausencia por lo que se
#arregla el problema de mayusculas
WATER = ifelse(toupper(WATER) == "Y", 1L, 0L),
# FCOVERPCT: % de cobertura baja (ground-level)
FCOVERPCT = as.numeric(FCOVERPCT),
# FCOVERPCT^2: término cuadrático
# El artículo encontró relación no lineal (máximo a ~70% de cobertura)
FCOVERPCT2 = FCOVERPCT^2,
# SUBSTRATE: tipo de sustrato dominante
# El artículo condensa tipos poco frecuentes en "OTHER"
SUBSTRATE_CAT = case_when(
SUBSTRATE == "MUD" ~ "MUD",
SUBSTRATE %in% c("SAND","DRYSAND",
"WETSAND","DAMPSAND") ~ "SAND",
SUBSTRATE %in% c("SOIL","DRYSOIL","DAMPSOIL") ~ "SOIL",
SUBSTRATE %in% c("GRAVEL","COBBLE","BEDROCK") ~ "ROCK",
TRUE ~ "OTHER"
),
SUBSTRATE_CAT = relevel(factor(SUBSTRATE_CAT), ref = "MUD"),
# LCOVERTYP: tipo de cobertura baja
# Categorías del artículo: NONE , ANNUALS, OPENWATER, ROCK,
# AQVEG (vegetación acuática), WOODY, OTHER
LCOVERTYP = case_when(
FCOVERTYP == "NONE" ~ "NONE",
FCOVERTYP == "ANNUALS" ~ "ANNUALS",
FCOVERTYP == "OPENWATER" ~ "OPENWATER",
FCOVERTYP %in% c("ROCK","COBBLE") ~ "ROCK",
FCOVERTYP %in% c("TYPHA","POTA","ALGAE",
"JUNCUS","SCIRPUS") ~ "AQVEG",
FCOVERTYP %in% c("COARSEWOODYDEBRIS",
"FINEWOODYDEBRIS",
"ROOTS"
) ~ "WOODY",
TRUE ~ "OTHER"
),
LCOVERTYP = relevel(factor(LCOVERTYP), ref = "NONE"),
# OCOVERTYP_CAT: tipo de cobertura superior (overstory)
# NONE, WILLOW, JUNIPER, TREESPP (otras especies), OTHER
OCOVERTYP_CAT = case_when(
OCOVERTYP == "NONE" ~ "NONE",
OCOVERTYP == "WILLOW" ~ "WILLOW",
OCOVERTYP == "JUNIPER" ~ "JUNIPER",
OCOVERTYP %in% c("WALNUT","OAK","PINE","BEECH",
"COTTONWOOD","BIRCH","MAPLE",
"CEDAR") ~ "TREESPP",
TRUE ~ "OTHER"
),
OCOVERTYP_CAT = relevel(factor(OCOVERTYP_CAT), ref = "NONE"),
# OCOVERPRES: cobertura superior
OCOVERPRES = ifelse(OCOVERPCT > 10, 1L, 0L),
# STRATA: índice de emparejamiento
# Cada par (rana + punto aleatorio) comparte el mismo estrato.
# Las filas están organizadas en pares consecutivos: fila impar=rana,
# fila par=aleatoria.
STRATA = rep(seq_len(nrow(datos_raw) / 2), each = 2)
)
# 3. Modelo
# Regresión logística condicional (Cox estratificado por par)
# Variables incluidas según selección de autores:
# WATER, FCOVERPCT, FCOVERPCT^2, LCOVERTYP, OCOVERTYP_CAT, SUBSTRATE_CAT
modelo_final <- coxme(
Surv(rep(1, nrow(datos)), USE) ~ WATER +
FCOVERPCT +
FCOVERPCT2 +
LCOVERTYP +
OCOVERTYP_CAT +
SUBSTRATE_CAT +
strata(STRATA) +
(1|FROGID),
data = datos
)
summary(modelo_final)
coefs <- as.data.frame(coef(summary(modelo_final)))
# Agregar columna de significancia
coefs$signif <- ifelse(coefs$p < 0.001, "***",
ifelse(coefs$p < 0.01, "**",
ifelse(coefs$p < 0.05, "*",
ifelse(coefs$p < 0.1, ".",
""))))
# Tabla de resultados del modelo
colnames(coefs) <- c("coef", "exp(coef)", "se(coef)", "z", "p", "")
kable(coefs,
digits = 4,
caption = "Resultados del modelo coxme — efectos fijos")
# Crear el data frame manualmente con los datos de summary
tabla_efectos_aleatorios <- data.frame(
Modelo = c("Integrated loglik", "Penalized loglik"),
Chisq = c(414.3, 414.3),
df = c(18, 17),
p = c(0, 0),
AIC = c(378.3, 289.4),
BIC = c(380.3, 296.3)
)
# Generar la tabla
knitr::kable(tabla_efectos_aleatorios,
digits = 1,
align = c("l", "c", "c", "c", "c", "c"),
caption = "Evaluación de la necesidad del efecto aleatorio por individuo (FROGID)")
#Calculo de Odds Ratios (OR) y sus IC ;ara el modelo final
OR_tabla <- round(exp(cbind(OR = coef(modelo_final),
confint(modelo_final))), 3)
cat("\n===== Odds Ratios del modelo final =====\n")
print(OR_tabla)
# Modelos reducidos
m_sin_water <- coxme(
Surv(rep(1, nrow(datos)), USE) ~
FCOVERPCT + FCOVERPCT2 + LCOVERTYP + OCOVERTYP_CAT +
SUBSTRATE_CAT + strata(STRATA) + (1 | FROGID),
data = datos
)
m_sin_fcovpct <- coxme(
Surv(rep(1, nrow(datos)), USE) ~
WATER + LCOVERTYP + OCOVERTYP_CAT +
SUBSTRATE_CAT + strata(STRATA) + (1 | FROGID),
data = datos
)
m_sin_fcovpct2 <- coxme(
Surv(rep(1, nrow(datos)), USE) ~
WATER + FCOVERPCT + LCOVERTYP + OCOVERTYP_CAT +
SUBSTRATE_CAT + strata(STRATA) + (1 | FROGID),
data = datos
)
m_sin_lcovtyp <- coxme(
Surv(rep(1, nrow(datos)), USE) ~
WATER + FCOVERPCT + FCOVERPCT2 + OCOVERTYP_CAT +
SUBSTRATE_CAT + strata(STRATA) + (1 | FROGID),
data = datos
)
m_sin_ocovtyp <- coxme(
Surv(rep(1, nrow(datos)), USE) ~
WATER + FCOVERPCT + FCOVERPCT2 + LCOVERTYP +
SUBSTRATE_CAT + strata(STRATA) + (1 | FROGID),
data = datos
)
m_sin_sub <- coxme(
Surv(rep(1, nrow(datos)), USE) ~
WATER + FCOVERPCT + FCOVERPCT2 + LCOVERTYP + OCOVERTYP_CAT +
strata(STRATA) + (1 | FROGID),
data = datos
)
# Función auxiliar coxme
lr_test_coxme <- function(modelo_completo, modelo_reducido, nombre) {
ll_completo <- modelo_completo$loglik[2]
ll_reducido <- modelo_reducido$loglik[2]
df_completo <- modelo_completo$df[1]
df_reducido <- modelo_reducido$df[1]
chi2 <- 2 * (ll_completo - ll_reducido)
df <- round(df_completo - df_reducido)
p <- pchisq(chi2, df = df, lower.tail = FALSE)
tibble(
Variable = nombre,
df = df,
chi2 = round(chi2, 2),
p_valor = round(p, 4)
)
}
# Tabla con dplyr
tabla2 <- bind_rows(
lr_test_coxme(modelo_final, m_sin_water, "WATER"),
lr_test_coxme(modelo_final, m_sin_fcovpct, "LCOVERPCT"),
lr_test_coxme(modelo_final, m_sin_fcovpct2, "LCOVERPCT^2"),
lr_test_coxme(modelo_final, m_sin_lcovtyp, "LCOVERTYPE"),
lr_test_coxme(modelo_final, m_sin_ocovtyp, "OCOVERTYPE"),
lr_test_coxme(modelo_final, m_sin_sub, "SUBSTRATE")
) |>
mutate(
Significancia = case_when(
p_valor < 0.001 ~ "***",
p_valor < 0.01 ~ "**",
p_valor < 0.05 ~ "*",
p_valor < 0.1 ~ ".",
TRUE ~ ""
)
) |>
arrange(p_valor) |> # Primero ordenamos numéricamente
mutate(
# Pasamos el p_valor a texto: si es menor a 0.001, ponemos "< 0.001"
# Si no, mantenemos el número formateado con 4 decimales
p_valor = case_when(
p_valor < 0.001 ~ "< 0.001",
TRUE ~ sprintf("%.4f", p_valor)
)
)
print(tabla2)
# Figura 1: OR ~ % de cobertura baja
# Tiene Relación cuadrática. OR relativo a 0% de cobertura (referencia = 1).
b1 <- coef(modelo_final)["FCOVERPCT"]
b2 <- coef(modelo_final)["FCOVERPCT2"]
V <- vcov(modelo_final)
i1 <- which(names(coef(modelo_final)) == "FCOVERPCT")
i2 <- which(names(coef(modelo_final)) == "FCOVERPCT2")
cover_seq <- seq(0, 100, by = 1)
log_OR <- b1 * cover_seq + b2 * cover_seq^2
# Varianza del log-OR por método delta
var_logOR <- cover_seq^2 * V[i1,i1] + cover_seq^4 * V[i2,i2] +
2 * cover_seq^3 * V[i1,i2]
se_logOR <- sqrt(pmax(var_logOR, 0))
fig2_df <- data.frame(
cover = cover_seq,
OR = exp(log_OR),
OR_low = exp(log_OR - 1.96 * se_logOR),
OR_high = exp(log_OR + 1.96 * se_logOR)
)
max_cover <- cover_seq[which.max(fig2_df$OR)]
cat("\nMáximo de OR en % de cobertura:", max_cover,
"| Artículo: ~70%\n")
fig2 <- ggplot(fig2_df, aes(x = cover, y = OR)) +
geom_ribbon(aes(ymin = OR_low, ymax = OR_high),
fill = "grey80", color = "grey60", linetype = "dashed") +
geom_line(linewidth = 0.9) +
geom_vline(xintercept = max_cover, linetype = "dashed") +
labs(x = "% Low-Lying Cover",
y = "Odds Ratio",
caption = paste0("Selection maximized at ", max_cover, "% cover")) +
scale_x_continuous(breaks = seq(0, 100, 20)) +
theme_classic(base_size = 13)
print(fig2)
# Figura 2: OR por tipo de cobertura y sustrato
# Extraer coeficientes del modelo como data.frame con IC
coef_df <- data.frame(
term = names(coef(modelo_final)),
coef = coef(modelo_final),
se = sqrt(diag(vcov(modelo_final)))
) %>%
mutate(OR = exp(coef),
OR_low = exp(coef - 1.96 * se),
OR_high = exp(coef + 1.96 * se))
# Panel a: tipo de cobertura baja
panel_a_df <- coef_df %>%
filter(grepl("^LCOVERTYP", term)) %>%
mutate(
label = recode(sub("LCOVERTYP", "", term),
ANNUALS = "Annuals",
AQVEG = "Aq Veg",
OPENWATER = "Open Water",
ROCK = "Rock",
WOODY = "Woody Cover",
OTHER = "Other"),
label = fct_reorder(label, OR, .desc = TRUE)
)
fig3a <- ggplot(panel_a_df, aes(x = label, y = OR)) +
geom_hline(yintercept = 1, linetype = "dashed", color = "grey50") +
geom_errorbar(aes(ymin = OR_low, ymax = OR_high), width = 0.2) +
geom_point(size = 3) +
annotate("text", x = 0.7, y = 1.15, label = "(No Cover)",
size = 3.2, color = "grey40", hjust = 0) +
labs(x = "Low-Lying Cover Type", y = "Odds Ratio", tag = "(a)") +
theme_classic(base_size = 12)
# Panel (b): tipo de cobertura superior
panel_b_df <- coef_df %>%
filter(grepl("^OCOVERTYP_CAT", term)) %>%
mutate(
label = recode(sub("OCOVERTYP_CAT", "", term),
WILLOW = "Willow",
JUNIPER = "Juniper",
TREESPP = "Tree spp.",
OTHER = "Other"),
label = fct_reorder(label, OR, .desc = TRUE)
)
fig3b <- ggplot(panel_b_df, aes(x = label, y = OR)) +
geom_hline(yintercept = 1, linetype = "dashed", color = "grey50") +
geom_errorbar(aes(ymin = OR_low, ymax = OR_high), width = 0.2) +
geom_point(size = 3) +
annotate("text", x = 0.7, y = 1.04, label = "(No Cover)",
size = 3.2, color = "grey40", hjust = 0) +
labs(x = "Overstory Cover Type", y = "Odds Ratio", tag = "(b)") +
theme_classic(base_size = 12)
# Panel c: tipo de sustrato
panel_c_df <- coef_df %>%
filter(grepl("^SUBSTRATE_CAT", term)) %>%
mutate(
label = recode(sub("SUBSTRATE_CAT", "", term),
SAND = "Sand",
SOIL = "Soil",
ROCK = "Rock",
OTHER = "Other"),
label = fct_reorder(label, OR, .desc = TRUE)
)
fig3c <- ggplot(panel_c_df, aes(x = label, y = OR)) +
geom_hline(yintercept = 1, linetype = "dashed", color = "grey50") +
geom_errorbar(aes(ymin = OR_low, ymax = OR_high), width = 0.2) +
geom_point(size = 3) +
annotate("text", x = 0.7, y = 1.05, label = "(Mud Substrate)",
size = 3.2, color = "grey40", hjust = 0) +
labs(x = "Substrate Type", y = "Odds Ratio", tag = "(c)") +
theme_classic(base_size = 12)
# Figura combinada y figura c por separado
fig3_combined <- fig3a / fig3b
print(fig3_combined)
print(fig3c)
# Cada fila es un OR comparado entre el articulo original
# y la replicacion hecha en el Rmd del taller.
# Columnas:
# Variable -> variable o categoria comparada
# Articulo_valor -> OR reportado en Hinderer et al. 2021 (IC 95%)
# Rmd_valor -> OR obtenido en la replicacion (IC 95%)
# Estado -> Coincide / Parcial / Diferencia
# Causa_probable -> explicacion de la discrepancia
tabla <- tibble(
Variable = c(
# Presencia de agua
"WATER",
# Cobertura baja: efecto continuo y pico de la curva cuadratica
"Pico OR cobertura baja (%)", # maximo de la curva cuadratica
"FCOVERPCT (por unidad)", # efecto lineal por cada 1% de cobertura
# Categorias de tipo de cobertura baja (referencia = NONE)
"WOODY vs NONE",
"ROCK vs NONE",
"AQVEG vs NONE",
# Categorias de tipo de cobertura alta (referencia = NONE)
"Willow vs NONE",
"Tree spp. vs NONE",
"Juniper vs NONE",
# Categorias de sustrato (referencia = MUD)
"SAND vs MUD",
"SOIL vs MUD",
"ROCK vs MUD"
),
# Valores del articulo original (Hinderer et al., 2021)
# Estimado puntual con IC 95% entre parentesis
Articulo_valor = c(
"2.92 (2.78 - 3.07)",
"~70%",
"No reportado", # el articulo no desglosa el coeficiente lineal
"6.03 (2.83 - 12.85)",
"4.14 (2.01 - 8.53)",
"3.98 (1.83 - 8.65)",
"2.22 (2.09 - 2.36)",
"1.99 (1.90 - 2.09)",
"1.54 (1.38 - 1.71)",
"0.43 (0.38 - 0.49)",
"0.26 (0.23 - 0.30)",
"0.16 (0.14 - 0.18)"
),
# Valores obtenidos en la replicacion (Taller Rmd)
# El modelo usa coxme en lugar de coxph, incorporando el efecto
# aleatorio (1|FROGID) de forma explicita, lo que produce IC
# mas amplios en casi todos los OR
Rmd_valor = c(
"2.944 (1.929 - 4.493)", # mismo punto estimado, IC mas amplio
"~68 - 69%", # pico matematico exacto; articulo redondea a 70%
"1.087 (1.064 - 1.111)", # disponible en coxme, no reportado en articulo
"7.052 (1.130 - 44.008)", # IC muy amplio por efecto aleatorio
"3.654 (0.712 - 18.745)", # IC cruza 1; pierde significancia individual
"4.081 (0.745 - 22.367)", # punto estimado similar; IC mucho mas amplio
"2.126 (1.328 - 3.405)", # mismo punto estimado; IC moderadamente mas amplio
"1.919 (1.234 - 2.984)", # misma direccion; IC moderadamente mas amplio
"IC cruza 1 (NS)", # pierde significancia individual en la replicacion
"0.437 (IC similar)", # practicamente identico al articulo
"0.257 (IC similar)", # practicamente identico al articulo
"0.168 (0.082 - 0.344)" # diferencia minima; misma direccion y magnitud
),
# Estado de la comparacion:
# Coincide -> diferencia minima o nula
# Parcial -> misma direccion, magnitud o IC algo distintos
# Diferencia -> discrepancia relevante en magnitud o significancia
Estado = c(
"Diferencia", # WATER: IC mucho mas amplio
"Coincide", # pico cuadratico: practicamente igual (~1-2% de diferencia)
"Diferencia", # no comparable directamente (no reportado en articulo)
"Diferencia", # WOODY: IC enormemente mas amplio
"Diferencia", # ROCK: IC cruza 1
"Diferencia", # AQVEG: punto estimado similar pero IC muy amplio
"Parcial", # Willow: mismo punto estimado, IC algo mas amplio
"Parcial", # Tree spp.: misma direccion, IC algo mas amplio
"Diferencia", # Juniper: pierde significancia individual
"Coincide", # SAND: practicamente identico
"Coincide", # SOIL: practicamente identico
"Parcial" # ROCK: diferencia minima, misma direccion
),
# Explicacion de cada discrepancia
Causa_probable = c(
"IC mucho mas amplio; efecto aleatorio (1|FROGID) aumenta la incertidumbre de los estimados",
"Diferencia de ~1-2%; articulo redondea a 70%",
"No reportado en articulo; coeficiente disponible en coxme",
"Efecto aleatorio + estratificacion por pares genera IC muy amplios en categorias de baja frecuencia",
"IC cruza 1; efecto aleatorio redistribuye varianza reduciendo precision individual",
"Punto estimado similar; IC amplio por baja frecuencia de la categoria y efecto aleatorio",
"Mismo punto estimado; coxme produce IC mas amplios al modelar variabilidad entre ranas",
"Misma direccion e interpretacion; IC ligeramente mas amplio por efecto aleatorio",
"Pierde significancia individual; efecto aleatorio absorbe parte de la varianza de Juniper",
"Sustrato muy frecuente en los datos; robusto al cambio de modelo",
"Sustrato muy frecuente en los datos; robusto al cambio de modelo",
"Diferencia minima (~5%); misma direccion y conclusion biologica"
)
)
#Generar tabla
kable(
tabla,
col.names = c(
"Variable / Categoria",
"Articulo — OR (IC 95%)",
"Rmd — OR (IC 95%)",
"Estado",
"Causa probable"
),
caption = "Comparacion de Odds Ratios: Hinderer et al. (2021) vs replicacion Rmd",
align = c("l", "c", "c", "c", "l")
)Referencias
Hinderer RK, Litt AR, McCaffery M (2020). Habitat selection by a threatened desert amphibian. Ecol Evol.2021;11:536–546. https://doi.org/10.1002/ece3.7074
Referencias de Software Utilizado
Therneau T (2024). A Package for Survival Analysis in R. R package version 3.8-3, https://CRAN.R-project.org/package=survival. Terry M. Therneau, Patricia M. Grambsch (2000). Modeling Survival Data: Extending the Cox Model. Springer, New York. ISBN 0-387-98784-3.
Wickham H (2016). ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. ISBN 978-3-319-24277-4, https://ggplot2.tidyverse.org.
Wickham H, François R, Henry L, Müller K, Vaughan D (2023). dplyr: A Grammar of Data Manipulation. doi:10.32614/CRAN.package.dplyr https://doi.org/10.32614/CRAN.package.dplyr, R package version 1.1.4, https://CRAN.R-project.org/package=dplyr.
Wickham H (2023). forcats: Tools for Working with Categorical Variables (Factors). doi:10.32614/CRAN.package.forcats https://doi.org/10.32614/CRAN.package.forcats, R package version 1.0.0, https://CRAN.R-project.org/package=forcats.
Pedersen T (2025). patchwork: The Composer of Plots. doi:10.32614/CRAN.package.patchwork https://doi.org/10.32614/CRAN.package.patchwork, R package version 1.3.2, https://CRAN.R-project.org/package=patchwork.
Therneau TM (2024). coxme: Mixed Effects Cox Models. doi:10.32614/CRAN.package.coxme https://doi.org/10.32614/CRAN.package.coxme, R package version 2.2-22, https://CRAN.R-project.org/package=coxme.
Xie Y (2025). knitr: A General-Purpose Package for Dynamic Report Generation in R. R package version 1.50, https://yihui.org/knitr/. Xie Y (2015). Dynamic Documents with R and knitr, 2nd edition. Chapman and Hall/CRC, Boca Raton, Florida. ISBN 978-1498716963, https://yihui.org/knitr/. Xie Y (2014). “knitr: A Comprehensive Tool for Reproducible Research in R.” In Stodden V, Leisch F, Peng RD (eds.), Implementing Reproducible Computational Research. Chapman and Hall/CRC. ISBN 978-1466561595.