Se presentan los datos de rendimiento de 2 variedades de arroz: - Porte alto. - Porte baja. Expuestas a 3 niveles de fertilizacion nitrogenada (0, 100 y 200 t/ha)

Carga de librerias

#install.packages("tidyverse")
#install.packages("agricolae")

library(readr)
library(tidyverse)
library(agricolae)

Carga de datos

arroz <- read.csv("data/Datos_arroz.csv", sep = ";")

Cambiar las clases de las variables

arroz$Tratamiento <- as.factor(arroz$Tratamiento)
arroz$Variedad <- as.factor(arroz$Variedad)
arroz$Nitrogeno <- as.factor(arroz$Nitrogeno)
arroz$Bloque <- as.factor(arroz$Bloque)

# muestre la estructura de la data
str(arroz)
## 'data.frame':    24 obs. of  5 variables:
##  $ Tratamiento: Factor w/ 6 levels "T1","T2","T3",..: 1 1 1 1 2 2 2 2 3 3 ...
##  $ Variedad   : Factor w/ 2 levels "VA","VB": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Nitrogeno  : Factor w/ 3 levels "N0 ","N100","N200": 1 1 1 1 2 2 2 2 3 3 ...
##  $ Bloque     : Factor w/ 4 levels "1","2","3","4": 1 2 3 4 1 2 3 4 1 2 ...
##  $ Rendimiento: num  7.5 6 7 8.5 8.5 6.8 7.3 8.4 7.6 5.9 ...

Analisis exploratorio de datos

Se lleva a cabo un analisis grafico con el objetivo de revisar el efecto de las variables sobre la respuesta (rendimiento)

arroz %>% ggplot(aes(x = Nitrogeno, y = Rendimiento, fill = Nitrogeno)) +
  geom_boxplot() +
  geom_jitter(width = 0.2) +
  labs(y = "Rendimiento (t/ha)", title = "Efecto del Nitrogeno sobre el Rto", caption = "Author: Luis Fdo Delgado")

arroz %>% ggplot(aes(x = Variedad, y = Rendimiento, fill = Variedad)) +
  geom_boxplot() +
  geom_jitter(width = 0.2) +
  labs(y = "Rendimiento (t/ha)", title = "Efecto del variedad sobre el Rto", caption = "Author: Luis Fdo Delgado")

arroz %>% ggplot(aes(x = Variedad:Nitrogeno, y = Rendimiento, fill = Variedad)) +
  geom_boxplot() +
  geom_jitter(width = 0.2) +
  labs(y = "Rendimiento (t/ha)", title = "Efecto del variedad:Nitrogeno sobre el Rto", caption = "Author: Luis Fdo Delgado")

Graficos de interaccion

interaccion <- arroz %>% group_by(Variedad, Nitrogeno) %>% 
  summarise(Rento = mean(Rendimiento), 
            Rento_sd = sd(Rendimiento))


  interaccion %>% ggplot() +
  geom_line(aes(x = Nitrogeno, y = Rento, group = Variedad, color = Variedad), size = 0.4) +
  labs(y = "Rendimiento (t/ha)", title = "Efecto del Nitrogeno interactuando con la variedad", caption = "Author: Luis Fdo Delgado")

interaccion %>% ggplot() +
  geom_line(aes(x = Variedad, y = Rento, group = Nitrogeno, color = Nitrogeno), size = 0.4) +
  labs(y = "Rendimiento (t/ha)", title = "Efecto del variedad interactuando con la Nitrogeno", caption = "Author: Luis Fdo Delgado")

Modelacion bajo un diseno de bloques completos al azar con estructura factorial

modelo1 <- aov(Rendimiento ~ Bloque + Nitrogeno*Variedad, data = arroz)
summary(modelo1) # The best model
##                    Df Sum Sq Mean Sq F value   Pr(>F)    
## Bloque              3 17.578   5.859  82.014 1.57e-09 ***
## Nitrogeno           2  3.000   1.500  20.995 4.49e-05 ***
## Variedad            1  0.042   0.042   0.583 0.456904    
## Nitrogeno:Variedad  2  2.333   1.167  16.330 0.000172 ***
## Residuals          15  1.072   0.071                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
modelo2 <- aov(Rendimiento ~ Bloque + Tratamiento, data = arroz)
summary(modelo2) 
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Bloque       3 17.578   5.859   82.01 1.57e-09 ***
## Tratamiento  5  5.375   1.075   15.05 2.19e-05 ***
## Residuals   15  1.072   0.071                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
modelo3 <- aov(Rendimiento ~ Bloque + Nitrogeno + Variedad, data = arroz)
summary(modelo3)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Bloque       3 17.578   5.859  29.254 6.15e-07 ***
## Nitrogeno    2  3.000   1.500   7.489  0.00465 ** 
## Variedad     1  0.042   0.042   0.208  0.65409    
## Residuals   17  3.405   0.200                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Comparacion de medias de la interaccion

duncan_compara <- duncan.test(modelo1, c("Nitrogeno", "Variedad"), alpha = 0.05)
duncan_compara$groups
##         Rendimiento groups
## N200:VB        8.00      a
## N100:VA        7.75     ab
## N100:VB        7.50     bc
## N0 :VA         7.25      c
## N200:VA        7.25      c
## N0 :VB         6.50      d

Grafico de barras de las medias

bar.group(duncan_compara$groups, ylim = c(0, 10),  
          xlab = "Tratamientos", ylab = "Rendimiento medio (t/ha)", col = "purple")