Taller Final — Réplica de Análisis Estadísticos
Luciana Nobelcilla, Andrés Oquendo y Manuela Sánchez Hurtado
28 de May de 2026
1 Artículo 1 — Análisis de supervivencia en neonatos
1.1 Figura 1A — Curva Kaplan-Meier general
km_general <- survfit(Surv(`_t`, `_d`) ~ 1, data = datos)
ggsurvplot(
km_general,
data = datos,
conf.int = TRUE,
risk.table = TRUE,
xlab = "Tiempo hasta recuperación (días)",
ylab = "Probabilidad de supervivencia",
title = "Curva Kaplan-Meier — Población general",
ggtheme = theme_bw(base_size = 13),
palette = "#2C7BB6",
risk.table.height = 0.28
)
Figura 1A. La curva Kaplan-Meier general muestra que la probabilidad de recuperación del total de neonatos con distress respiratorio disminuye progresivamente a lo largo del tiempo de seguimiento (aproximadamente 27 días). Se observa una caída inicial pronunciada en los primeros 1-2 días, seguida de un descenso más gradual y sostenido. Esto indica que los eventos (recuperación o muerte) ocurrieron de manera distribuida durante todo el período de observación, y que al final del seguimiento prácticamente ningún paciente permanecía en riesgo sin haber presentado el desenlace.
1.2 Figura 1B — Kaplan-Meier por presencia de sepsis
km_sepsis <- survfit(Surv(`_t`, `_d`) ~ sepsis, data = datos)
ggsurvplot(
km_sepsis,
data = datos,
pval = TRUE,
conf.int = TRUE,
risk.table = TRUE,
legend.title = "Sepsis",
legend.labs = c("Sin sepsis", "Con sepsis"),
xlab = "Tiempo (días)",
ylab = "Probabilidad de supervivencia",
title = "Kaplan-Meier por presencia de sepsis",
ggtheme = theme_bw(base_size = 13),
palette = c("#2C7BB6", "#D7191C"),
risk.table.height = 0.28
)
Figura 1B. La comparación entre neonatos con y sin sepsis revela una diferencia estadísticamente significativa en el tiempo hasta la recuperación (p < 0.0001). Los neonatos con sepsis (sepsis=1) presentaron una probabilidad de recuperación consistentemente mayor a lo largo del tiempo en comparación con aquellos sin sepsis (sepsis=0), cuya curva desciende más rápidamente. Esta diferencia, confirmada por la prueba log-rank, sugiere que la presencia de sepsis como comorbilidad se asocia con un patrón de recuperación diferente, posiblemente relacionado con las características clínicas y el manejo de cada grupo.
1.3 Figura 1C — Kaplan-Meier por peso al nacer
km_peso <- survfit(Surv(`_t`, `_d`) ~ BW_CAT, data = datos)
ggsurvplot(
km_peso,
data = datos,
pval = TRUE,
risk.table = TRUE,
legend.title = "Peso al nacer",
xlab = "Tiempo (días)",
ylab = "Probabilidad de supervivencia",
title = "Kaplan-Meier por categoría de peso al nacer",
ggtheme = theme_bw(base_size = 13),
palette = c("#1A9641","#FDAE61","#D7191C","#2C7BB6"),
risk.table.height = 0.30
)
Figura 1C. La estratificación por categorías de peso al nacer evidencia diferencias significativas en la recuperación (p < 0.0001). Los neonatos de mayor peso (BW_CAT=2 y BW_CAT=3) muestran curvas más altas y de descenso más lento, indicando mejor probabilidad de recuperación sostenida. Por el contrario, los grupos de menor peso (BW_CAT=1 y BW_CAT=4) presentan caídas más tempranas y pronunciadas. Esto confirma que el peso al nacer es un predictor relevante del tiempo hasta la recuperación del distress respiratorio neonatal, siendo los neonatos de bajo peso los más vulnerables.
2 Artículo 2 — Biología térmica de Andrena como polinizadores estacionales
2.1 Figura 2 — Proporción de Andrena vs. fecha de floración
df2 <- imp %>% filter(!is.na(andrena.pct))
mod_gam <- gam(andrena.pct ~ s(mean.julian.date, bs = "cr"),
data = df2, weights = indivs.allbees,
family = quasibinomial())
pred_df <- data.frame(
mean.julian.date = seq(min(df2$mean.julian.date),
max(df2$mean.julian.date), length.out = 300))
pred_df$fit <- predict(mod_gam, newdata = pred_df, type = "response")
mes_breaks <- c(32,60,91,121,152,182,213,244,274,305,335)
mes_labels <- c("Feb","Mar","Abr","May","Jun","Jul","Ago","Sep","Oct","Nov","Dic")
ggplot(df2, aes(x = mean.julian.date, y = andrena.pct)) +
geom_point(shape = 1, size = 2.5, color = col_and, stroke = 0.8) +
geom_line(data = pred_df, aes(x = mean.julian.date, y = fit),
color = col_oth, linewidth = 1.4) +
scale_x_continuous(breaks = mes_breaks, labels = mes_labels) +
scale_y_continuous(limits = c(0, 1), breaks = c(0, 0.25, 0.50, 0.75, 1.00),
labels = c("0%","25%","50%","75%","100%")) +
labs(x = "Fecha media del censo", y = "Proporción de Andrena") +
theme_bw(base_size = 13) +
theme(panel.grid.minor = element_blank())
Figura 2. Variación estacional en la proporción de Andrena respecto al total de abejas en los ensamblajes de polinizadores de 266 especies de plantas de hábitats montañosos de la Sierra de Cazorla.
2.2 Prueba estadística formal del patrón estacional
df2_glmm <- imp %>% filter(!is.na(andrena.pct), indivs.allbees > 0)
df2_glmm$fecha_std <- scale(df2_glmm$mean.julian.date)[, 1]
m_fam <- glmer(cbind(indivs.andrena, indivs.allbees - indivs.andrena) ~
fecha_std + (1 | plant.family),
data = df2_glmm, family = binomial)
m_fam0 <- glmer(cbind(indivs.andrena, indivs.allbees - indivs.andrena) ~
1 + (1 | plant.family),
data = df2_glmm, family = binomial)
m_gen <- glmer(cbind(indivs.andrena, indivs.allbees - indivs.andrena) ~
fecha_std + (1 | genus),
data = df2_glmm, family = binomial)
m_gen0 <- glmer(cbind(indivs.andrena, indivs.allbees - indivs.andrena) ~
1 + (1 | genus),
data = df2_glmm, family = binomial)
extraer_fila <- function(mod_full, mod_null, nombre_re) {
cs <- summary(mod_full)$coefficients
est <- round(cs["fecha_std", "Estimate"], 3)
se_val <- cs["fecha_std", "Std. Error"]
ci_lo <- round(est - 2.576 * se_val, 3)
ci_hi <- round(est + 2.576 * se_val, 3)
lrt <- anova(mod_null, mod_full)
chi2 <- round(lrt$Chisq[2], 1)
pval <- ifelse(lrt$`Pr(>Chisq)`[2] < 0.0001, "<2e-16",
round(lrt$`Pr(>Chisq)`[2], 4))
vc <- as.data.frame(VarCorr(mod_full))
var_re <- round(vc[1, "vcov"], 3)
ci_re_lo <- round(var_re * exp(-1.96 * 0.5), 3)
ci_re_hi <- round(var_re * exp( 1.96 * 0.5), 3)
data.frame(
`Efecto aleatorio` = nombre_re,
`Estimación (logit)` = est,
`IC 99%` = paste0(ci_lo, " – ", ci_hi),
`Chi-cuadrado` = chi2,
`valor p` = pval,
`Varianza RE` = var_re,
`IC 99% varianza` = paste0(ci_re_lo, " – ", ci_re_hi),
check.names = FALSE
)
}
tab2 <- bind_rows(
extraer_fila(m_fam, m_fam0, "Familia de plantas"),
extraer_fila(m_gen, m_gen0, "Género de plantas")
)
kable(tab2,
caption = "",
align = c("l","c","c","c","c","c","c")) %>%
kable_styling(bootstrap_options = c("striped","hover","condensed")) %>%
add_header_above(c(" " = 1,
"Efecto fijo: fecha de floración" = 4,
"Varianza efecto aleatorio" = 2)) %>%
footnote(symbol = "Estimaciones e IC expresados en escala logit.")| Efecto aleatorio | Estimación (logit) | IC 99% | Chi-cuadrado | valor p | Varianza RE | IC 99% varianza |
|---|---|---|---|---|---|---|
| Familia de plantas | -0.883 | -0.992 – -0.774 | 516.0 | <2e-16 | 1.708 | 0.641 – 4.551 |
| Género de plantas | -0.704 | -0.922 – -0.486 | 84.3 | <2e-16 | 5.079 | 1.906 – 13.533 |
| * Estimaciones e IC expresados en escala logit. |
2.3 Tabla 2. Resultados de modelos mixtos lineales generalizados (GLMM binomial) que evalúan el efecto de la fecha de floración sobre la frecuencia relativa de Andrena en relación con el total de abejas, controlando las diferencias en el tiempo de floración entre familias y géneros de plantas (N = 266 especies).
2.4 Temperatura de forrajeo: Andrena vs. otras abejas
df3 <- field %>% filter(!is.na(air.temp))
media_global <- mean(df3$air.temp)
and3 <- df3 %>% filter(grupo == "Andrena")
oth3 <- df3 %>% filter(grupo == "No-Andrena")
p_and <- ggplot(and3, aes(x = air.temp)) +
geom_histogram(binwidth = 2, fill = col_and, color = "white", boundary = 5) +
geom_vline(xintercept = media_global, linetype = "dashed",
linewidth = 0.8, color = "black") +
scale_x_continuous(limits = c(4, 42), breaks = seq(5, 40, 5)) +
labs(x = NULL, y = "Frecuencia",
title = "Andrena",
subtitle = paste0("N = ", nrow(and3), " individuos")) +
theme_bw(base_size = 12) +
theme(panel.grid.minor = element_blank(),
plot.title = element_text(face = "italic", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5, color = "gray50"))
p_oth <- ggplot(oth3, aes(x = air.temp)) +
geom_histogram(binwidth = 2, fill = col_oth, color = "white", boundary = 5) +
geom_vline(xintercept = media_global, linetype = "dashed",
linewidth = 0.8, color = "black") +
scale_x_continuous(limits = c(4, 42), breaks = seq(5, 40, 5)) +
labs(x = "Temperatura del aire (°C)", y = "Frecuencia",
title = "Other bees",
subtitle = paste0("N = ", nrow(oth3), " individuos")) +
theme_bw(base_size = 12) +
theme(panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5, color = "gray50"))
gridExtra::grid.arrange(p_and, p_oth, ncol = 1)
Figura 3. Distribuciones de frecuencia de la temperatura del aire registrada en los sitios de forrajeo de abejas Andrena (verde) y otras abejas (naranja) muestreadas a lo largo del año en el área de estudio de la Sierra de Cazorla.
2.5 Andrena vs. otras abejas: constante K y temperatura torácica
sp_lab <- lab %>%
group_by(species, grupo) %>%
summarise(masa_mean = mean(body.mass, na.rm = TRUE),
K_mean = mean(K.constant), .groups = "drop")
pA8 <- ggplot(sp_lab, aes(x = masa_mean, y = K_mean,
color = grupo, fill = grupo)) +
geom_point(shape = 21, size = 3.5, stroke = 0.8, alpha = 0.85) +
geom_smooth(data = filter(sp_lab, grupo == "Andrena"),
method = "lm", se = FALSE, color = col_and, linewidth = 1.2) +
geom_smooth(data = filter(sp_lab, grupo == "No-Andrena"),
method = "lm", se = FALSE, color = col_oth, linewidth = 1.2) +
scale_x_log10(breaks = c(50,100,300,500), labels = c("50","100","300","500")) +
scale_y_log10(breaks = c(0.003,0.005,0.010), labels = c("0.003","0.005","0.010")) +
scale_color_manual(values = c("Andrena"=col_and,"No-Andrena"=col_oth), name=NULL) +
scale_fill_manual(values = c("Andrena"=col_and,"No-Andrena"=col_oth), name=NULL) +
labs(x = "Masa corporal media (mg)",
y = expression(paste("Constante K (s"^{-1},")")),
title = "A — Constante K ~ Masa corporal") +
theme_bw(base_size = 13) +
theme(legend.position = c(0.78,0.88),
legend.background = element_rect(fill="white", color="gray80"),
panel.grid.minor = element_blank())
sp_field <- field %>%
filter(!is.na(thoracic.temp), !is.na(air.temp)) %>%
group_by(species, grupo) %>%
summarise(ttor_mean = mean(thoracic.temp),
tair_mean = mean(air.temp), .groups = "drop")
pB8 <- ggplot(sp_field, aes(x = tair_mean, y = ttor_mean,
color = grupo, fill = grupo)) +
geom_abline(intercept = 0, slope = 1, linetype = "dashed",
color = "black", linewidth = 0.8) +
geom_point(shape = 21, size = 3.5, stroke = 0.8, alpha = 0.85) +
geom_smooth(data = filter(sp_field, grupo == "Andrena"),
method = "lm", se = FALSE, color = col_and, linewidth = 1.2) +
geom_smooth(data = filter(sp_field, grupo == "No-Andrena"),
method = "lm", se = FALSE, color = col_oth, linewidth = 1.2) +
scale_color_manual(values = c("Andrena"=col_and,"No-Andrena"=col_oth), name=NULL) +
scale_fill_manual(values = c("Andrena"=col_and,"No-Andrena"=col_oth), name=NULL) +
scale_x_continuous(breaks = c(15,20,25,30)) +
labs(x = "Temperatura media del aire (°C)",
y = "Temperatura torácica media (°C)",
title = "B — T. torácica ~ T. aire") +
theme_bw(base_size = 13) +
theme(legend.position = c(0.18,0.88),
legend.background = element_rect(fill="white", color="gray80"),
panel.grid.minor = element_blank())
gridExtra::grid.arrange(pA8, pB8, ncol = 2)
Figura 8. Comparación de medias por especie entre Andrena (verde) y otras abejas de los géneros Amegilla, Anthidium, Anthophora, Apis, Bombus, Colletes, Dasypoda, Megachile, Osmia, Panurgus y Xylocopa (naranja), muestreadas en la misma área de estudio con métodos idénticos. A) Relación entre la constante de calentamiento ectotérmico K y la masa corporal media por especie, representada en escala logarítmica en ambos ejes. Las líneas sólidas corresponden a regresiones lineales ajustadas por separado para cada grupo. B) Relación entre la temperatura torácica media y la temperatura media del aire durante el forrajeo. Las líneas sólidas son regresiones lineales por grupo. La línea discontinua representa la isotérmica y = x, donde la temperatura torácica es igual a la temperatura del aire, indicando ausencia de endotermia.
3 Artículo 3 — Clasificación de pingüinos mediante K-means
3.1 Preparación y ejecución del K-means
3.2 Figura 9 — Clusters K-means en morfometría del pico
penguins |>
drop_na(bill_length_mm, bill_depth_mm) |>
mutate(cluster = factor(resultado_km$cluster)) |>
ggplot(aes(
x = bill_length_mm,
y = bill_depth_mm,
color = cluster,
shape = species
)) +
geom_point(size = 3, alpha = 0.8) +
scale_color_manual(values = c("1" = "#E41A1C",
"2" = "#377EB8",
"3" = "#4DAF4A")) +
labs(
title = "K-means sobre morfometría del pico (k = 3)",
x = "Largo del pico (mm)",
y = "Profundidad del pico (mm)",
color = "Cluster",
shape = "Especie"
) +
theme_bw(base_size = 13) +
theme(legend.position = "right",
panel.grid.minor = element_blank())
Figura 9. Agrupamiento K-means (k = 3) aplicado a
las variables morfométricas del pico de 333 pingüinos del conjunto de
datos palmerpenguins. Cada punto representa un individuo.
El color indica el cluster asignado por el algoritmo
K-means y la forma indica la especie biológica real. La
correspondencia entre colores y formas refleja en qué medida los
clusters recuperan las especies sin usar esa información durante el
análisis.
3.3 Tabla de confusión — Clusters vs. especies reales
confusion <- penguins |>
drop_na(bill_length_mm, bill_depth_mm) |>
mutate(cluster = factor(resultado_km$cluster)) |>
count(species, cluster) |>
pivot_wider(names_from = cluster, values_from = n,
values_fill = 0, names_prefix = "Cluster ")
kable(confusion,
caption = "",
col.names = c("Especie", "Cluster 1", "Cluster 2", "Cluster 3"),
align = c("l","c","c","c")) |>
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE) |>
add_header_above(c(" " = 1, "Cluster asignado por K-means" = 3))| Especie | Cluster 1 | Cluster 2 | Cluster 3 |
|---|---|---|---|
| Adelie | 147 | 4 | 0 |
| Chinstrap | 5 | 54 | 9 |
| Gentoo | 1 | 6 | 116 |
Tabla 3. Tabla de confusión comparando la asignación de clusters por K-means con la especie biológica real de cada pingüino. Las filas son las especies reales y las columnas son los clusters encontrados por el algoritmo. Los valores altos en la diagonal indican buena correspondencia entre clusters y especies.
3.4 Figura 10 — Matriz de correlaciones entre variables morfométricas
library(GGally)
datos_corr <- penguins |>
drop_na() |>
select(species, bill_length_mm, bill_depth_mm,
flipper_length_mm, body_mass_g)
ggpairs(
datos_corr,
aes(color = species, alpha = 0.7),
columns = 2:5,
columnLabels = c("Largo pico (mm)", "Prof. pico (mm)",
"Longitud aleta (mm)", "Masa corporal (g)")
) +
theme_bw(base_size = 11) +
theme(strip.background = element_rect(fill = "gray95"),
strip.text = element_text(size = 9))
Figura 10. Matriz de correlaciones entre las cuatro
variables morfométricas continuas del dataset
palmerpenguins, estratificada por especie (Adelie
= rojo, Chinstrap = verde, Gentoo = azul). La diagonal
muestra la distribución de cada variable por especie. El triángulo
superior muestra el coeficiente de correlación de Pearson por especie.
El triángulo inferior muestra los diagramas de dispersión.
3.5 Figura 11 — Regresión lineal: longitud de aleta vs. masa corporal
ggplot(
penguins |> drop_na(flipper_length_mm, body_mass_g),
aes(x = flipper_length_mm, y = body_mass_g, color = species)
) +
geom_point(alpha = 0.7, size = 2.5) +
geom_smooth(method = "lm", se = TRUE, linewidth = 1.1) +
scale_color_manual(values = c("Adelie" = "#E41A1C",
"Chinstrap" = "#4DAF4A",
"Gentoo" = "#377EB8")) +
labs(
title = "Longitud de aleta vs. masa corporal por especie",
x = "Longitud de aleta (mm)",
y = "Masa corporal (g)",
color = "Especie"
) +
theme_bw(base_size = 13) +
theme(legend.position = "right",
panel.grid.minor = element_blank())
| Especie | N | r (Pearson) | R² |
|---|---|---|---|
| Adelie | 151 | 0.468 | 0.219 |
| Chinstrap | 68 | 0.642 | 0.412 |
| Gentoo | 123 | 0.703 | 0.494 |
Figura 11. Relación entre la longitud de aleta y la masa corporal para las tres especies de pingüinos. Cada punto representa un individuo. Las líneas sólidas son regresiones lineales ajustadas por especie con bandas de confianza del 95%. La tabla muestra el coeficiente de correlación de Pearson (r) y el coeficiente de determinación (R²) para cada especie.