1 Cómo usar este documento

Este es un archivo R Markdown (.Rmd): mezcla texto explicativo con bloques de código de R (las zonas grises que empiezan con ```{r}).

Tienes dos formas de usarlo en clase:

  1. Copiar y pegar: abre RStudio, crea un script nuevo y copia cada bloque de código para ejecutarlo línea por línea (Ctrl + Enter). Ideal para principiantes.
  2. Tejer (Knit): abre este .Rmd en RStudio y presiona el botón Knit. R ejecutará todo el código y generará un informe en HTML con los resultados y las gráficas incrustadas.

Todo el código de este documento fue probado y funciona. Si algo falla, casi siempre es porque falta instalar un paquete (ver la siguiente sección).

2 Instalación y carga de paquetes

En R, un paquete es una caja de herramientas que alguien ya programó. Solo se instala una vez (con install.packages) y se carga en cada sesión (con library).

# EJECUTA ESTO UNA SOLA VEZ (quita el # y córrelo en la consola).
# Requiere conexión a internet. Puede tardar unos minutos.
install.packages(c("agricolae",   # datos y estadística clásica del agro
                   "agridat",     # +390 conjuntos de datos agrícolas reales
                   "ggplot2",     # el paquete estrella para gráficos
                   "dplyr"))      # para resumir y manipular tablas
# ESTO SÍ se ejecuta cada vez que abres R:
library(agricolae)
library(agridat)
library(ggplot2)
library(dplyr)

¿Qué contienen?

  • agricolae: paquete peruano (UNAM La Molina) hecho específicamente para investigación agrícola. Trae datos como sweetpotato (camote), cotton (algodón) y soil (análisis de suelos).
  • agridat: una enciclopedia con 392 conjuntos de datos reales de experimentos agrícolas (maíz, soya, trigo, plagas, suelos, etc.).
  • ggplot2: para construir gráficos “por capas”, sumando elementos con el símbolo +.
  • dplyr: para calcular promedios por grupo de forma muy legible (usa el operador %>%).

3 Explorar un conjunto de datos

Regla de oro: antes de graficar o calcular nada, mira cómo son tus datos.

Usaremos sweetpotato de agricolae: un experimento donde se midió el rendimiento (yield) de camote bajo 4 tratamientos virales (virus).

data(sweetpotato)   # carga el conjunto de datos en memoria

head(sweetpotato)   # muestra las primeras filas
##   virus yield
## 1    cc  28.5
## 2    cc  21.7
## 3    cc  23.0
## 4    fc  14.9
## 5    fc  10.6
## 6    fc  13.1
str(sweetpotato)    # estructura: qué variables hay y de qué tipo
## 'data.frame':    12 obs. of  2 variables:
##  $ virus: Factor w/ 4 levels "cc","fc","ff",..: 1 1 1 2 2 2 3 3 3 4 ...
##  $ yield: num  28.5 21.7 23 14.9 10.6 13.1 41.8 39.2 28 38.2 ...
dim(sweetpotato)    # número de filas y columnas
## [1] 12  2

Interpretación de str:

  • virus es un factor (variable categórica) con 4 niveles: cc, fc, ff, oo.
  • yield es numérica (el rendimiento medido).

4 Estadística descriptiva simple

4.1 Sobre toda la muestra

summary(sweetpotato$yield)     # mínimo, máximo, media, mediana, cuartiles
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   10.60   20.00   28.25   27.62   38.45   41.80
mean(sweetpotato$yield)        # media (promedio)
## [1] 27.625
median(sweetpotato$yield)      # mediana
## [1] 28.25
sd(sweetpotato$yield)          # desviación estándar (qué tan dispersos están los datos)
## [1] 11.07874

4.2 Resumen por grupo (lo más útil en el agro)

Casi siempre queremos comparar tratamientos. Con dplyr calculamos, para cada virus, el número de datos, la media, la desviación estándar y el error estándar (la media dividida por la raíz de n):

resumen <- sweetpotato %>%
  group_by(virus) %>%
  summarise(
    n     = n(),
    media = mean(yield),
    de    = sd(yield),
    ee    = sd(yield) / sqrt(n())   # error estándar de la media
  )

resumen
## # A tibble: 4 × 5
##   virus     n media    de    ee
##   <fct> <int> <dbl> <dbl> <dbl>
## 1 cc        3  24.4  3.61  2.08
## 2 fc        3  12.9  2.16  1.25
## 3 ff        3  36.3  7.33  4.23
## 4 oo        3  36.9  4.3   2.48

Ya se ve que los tratamientos ff y oo rinden más que cc y fc. Vamos a visualizarlo.

5 Gráficos con ggplot2

La lógica de ggplot2 siempre es la misma:

ggplot(datos, aes(x = ..., y = ...)) +   # 1) qué variables van en cada eje
  geom_XXX() +                            # 2) qué tipo de gráfico (puntos, cajas, barras...)
  labs(...) +                             # 3) títulos y nombres de ejes
  theme_minimal()                         # 4) apariencia

5.1 Diagrama de cajas (boxplot)

Ideal para comparar la distribución de una variable numérica entre varios grupos. La caja muestra el 50 % central de los datos; la línea del medio es la mediana.

ggplot(sweetpotato, aes(x = virus, y = yield, fill = virus)) +
  geom_boxplot(alpha = 0.8) +
  geom_jitter(width = 0.12, alpha = 0.5) +   # agrega los puntos reales, desplazados
  labs(title = "Rendimiento de camote por tratamiento viral",
       x = "Tratamiento (virus)",
       y = "Rendimiento (kg/parcela)") +
  theme_minimal() +
  theme(legend.position = "none")            # ocultamos la leyenda (es redundante)

Cómo leerlo: cajas más altas = mayor rendimiento. Cajas más “largas” = datos más dispersos. Los puntos por fuera de los bigotes serían valores atípicos.

5.2 Diagrama de barras usando medias (± error)

Para “barras de promedios” el orden correcto es: 1) calcular las medias, 2) graficarlas. Aquí usamos un experimento de control de plagas (beall.webworms de agridat): se contó el número de gusanos por parcela con y sin insecticida (spray).

data(beall.webworms)

# 1) calculamos la media y el error estándar por grupo
res_plaga <- beall.webworms %>%
  group_by(spray) %>%
  summarise(media = mean(y),
            ee    = sd(y) / sqrt(n()))

res_plaga
## # A tibble: 2 × 3
##   spray media     ee
##   <fct> <dbl>  <dbl>
## 1 N     1.13  0.0527
## 2 Y     0.458 0.0292
# 2) graficamos las medias como barras, con barras de error
ggplot(res_plaga, aes(x = spray, y = media, fill = spray)) +
  geom_col(width = 0.6) +                                  # barras a partir de valores ya calculados
  geom_errorbar(aes(ymin = media - ee, ymax = media + ee), # barras de error (± 1 EE)
                width = 0.2) +
  scale_x_discrete(labels = c("N" = "Sin insecticida",
                              "Y" = "Con insecticida")) +
  labs(title = "Número medio de gusanos por parcela",
       x = NULL,
       y = "Gusanos por parcela (media ± EE)") +
  theme_minimal() +
  theme(legend.position = "none")

Interpretación agronómica: el insecticida reduce a menos de la mitad la población media de gusanos. Las barras de error, cortas y sin traslape, sugieren que la diferencia es consistente.

Nota importante: geom_col() dibuja la altura que TÚ le das (las medias que calculaste). No lo confundas con geom_bar(), que por defecto cuenta filas.

5.3 Diagrama de dispersión y correlación

Sirve para ver la relación entre dos variables numéricas. Usaremos australia.soybean (agridat): contenido de proteína vs. aceite en granos de soya.

Primero medimos la correlación con un número, r, entre -1 y 1:

data(australia.soybean)

# Coeficiente de correlación de Pearson
r <- cor(australia.soybean$protein, australia.soybean$oil)
round(r, 3)
## [1] -0.758
# Prueba estadística: ¿la correlación es significativa?
cor.test(australia.soybean$protein, australia.soybean$oil)
## 
##  Pearson's product-moment correlation
## 
## data:  australia.soybean$protein and australia.soybean$oil
## t = -24.956, df = 462, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.7939767 -0.7160695
## sample estimates:
##      cor 
## -0.75771

r es cercano a -0.76: correlación negativa y fuerte. Ahora el gráfico, con una línea de tendencia (geom_smooth):

ggplot(australia.soybean, aes(x = protein, y = oil)) +
  geom_point(alpha = 0.5, color = "darkgreen") +
  geom_smooth(method = "lm", se = TRUE, color = "black") +  # recta de regresión
  labs(title = paste0("Proteína vs. aceite en soya (r = ", round(r, 2), ")"),
       x = "Proteína (%)",
       y = "Aceite (%)") +
  theme_minimal()

Interpretación: a mayor proteína, menor aceite. Es una relación conocida en mejoramiento de soya: seleccionar por más proteína tiende a bajar el aceite. La recta resume esa tendencia.

5.4 (Extra) Mapa de calor de correlaciones

Cuando tienes muchas variables numéricas (por ejemplo un análisis de suelos completo), un mapa de calor muestra todas las correlaciones de un vistazo. Usamos soil de agricolae.

data(soil)

# Elegimos algunas propiedades del suelo
variables <- c("pH", "MO", "CIC", "P", "K", "clay", "sand", "Ca", "Mg")

# Matriz de correlaciones -> tabla larga para ggplot
matriz  <- round(cor(soil[, variables]), 2)
tabla   <- as.data.frame(as.table(matriz))
names(tabla) <- c("Var1", "Var2", "cor")

ggplot(tabla, aes(x = Var1, y = Var2, fill = cor)) +
  geom_tile(color = "white") +
  geom_text(aes(label = cor), size = 3) +
  scale_fill_gradient2(low = "#b2182b", mid = "white", high = "#2166ac",
                       midpoint = 0, limits = c(-1, 1)) +
  labs(title = "Correlación entre propiedades del suelo",
       x = NULL, y = NULL, fill = "r") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Cómo leerlo: azul = correlación positiva, rojo = negativa. Por ejemplo, clay (arcilla) y sand (arena) están fuertemente correlacionadas de forma negativa: donde hay más de una, hay menos de la otra.

6 Bonus: comparar tratamientos con ANOVA

Esto es un poco más avanzado, pero es el pan de cada día en la experimentación agropecuaria. El ANOVA responde: ¿las diferencias entre tratamientos son reales o pura casualidad? Y la prueba de Tukey (HSD) de agricolae nos dice cuáles tratamientos difieren, agrupándolos con letras.

# ¿El rendimiento del camote depende del virus?
modelo <- aov(yield ~ virus, data = sweetpotato)
summary(modelo)   # si Pr(>F) < 0.05, hay diferencias significativas
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## virus        3 1170.2   390.1   17.34 0.000733 ***
## Residuals    8  179.9    22.5                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Comparación de medias con letras (a, b, c...)
comparacion <- HSD.test(modelo, "virus", group = TRUE)
comparacion$groups
##       yield groups
## oo 36.90000      a
## ff 36.33333     ab
## cc 24.40000     bc
## fc 12.86667      c

Cómo leer las letras: tratamientos que comparten una letra no difieren significativamente; los que no comparten ninguna letra sí difieren. Así se reportan resultados en artículos y tesis del agro.

7 Ejercicios propuestos

Para practicar con los mismos paquetes (todos los datos ya están instalados):

  1. Descriptiva: con agricolae::cotton (rendimiento de algodón por sitio), calcula la media y la desviación estándar de yield para cada site. Pista: cambia virus por site y sweetpotato por cotton en el código de la Sección 4.
  2. Boxplot: grafica yield según site en cotton. ¿Qué sitio rinde más?
  3. Barras de medias: en beall.webworms, compara el número medio de gusanos según trt (tratamiento) en lugar de spray.
  4. Correlación: en agridat::lasrosas.corn (maíz de precisión en Argentina), grafica la dispersión entre nitrógeno aplicado (nitro) y rendimiento (yield), y calcula su correlación.
  5. Reto: en lasrosas.corn, haz un boxplot del rendimiento (yield) según el nivel de nitrógeno (nf).

8 Apéndice: nota sobre los paquetes de datos

  • agricolae y agridat son fuentes confiables y muy usadas en docencia agropecuaria. Con data(package = "agridat") puedes ver la lista completa de sus 392 conjuntos de datos.
  • Si buscas específicamente datos de soya, agridat incluye varios: australia.soybean, nass.soybean, davidian.soybean, gauch.soy, pacheco.soybean, entre otros.
  • Sobre soymap: no corresponde a un paquete estándar de datos en CRAN que se pueda verificar, así que en este material se reemplazó por los datos de soya de agridat (ya probados). Si tienes una fuente concreta para soymap, compártela y adaptamos los ejemplos.

Para ver la ficha de ayuda de cualquier conjunto de datos y entender sus variables, usa ?:

?sweetpotato
?australia.soybean
?lasrosas.corn