Ejercicio 1

Con la base de datos genxenv que está disponible en la biblioteca agricolae, responda lo siguiente:

¿Existen diferencias estadísticas entre los ambientes (ENV) para la variable producción (YLD)?
¿Existen diferencias estadísticas entre los genotipos (GEN) para la variable producción (YLD)?
Estime las medias de cada ambiente junto al intervalo de confianza del 95%.
Ejecute pruebas de comparaciones múltiples para la fuente de variación ambientes (ENV) e informe los resultados. Nota: independientemente del modelo que utilice, recuerde reportar el diagnóstico de los residuales.

Estadísticos Descriptivos.

library(agricolae)
library(tidyverse)

## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --

## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.2     v dplyr   1.0.7
## v tidyr   1.1.3     v stringr 1.4.0
## v readr   1.4.0     v forcats 0.5.1

## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(janitor)

## 
## Attaching package: 'janitor'

## The following objects are masked from 'package:stats':
## 
##     chisq.test, fisher.test

data(genxenv)
datos <- genxenv %>% 
  mutate(
    ENV = factor(
      ENV), GEN = factor(GEN)) 

library(skimr)
skim(datos)

Data summary
Name	datos
Number of rows	250
Number of columns	3
_______________________
Column type frequency:
factor	2
numeric	1
________________________
Group variables	None

Variable type: factor

skim_variable	n_missing	complete_rate	ordered	n_unique	top_counts
ENV	0	1	FALSE	5	1: 50, 2: 50, 3: 50, 4: 50
GEN	0	1	FALSE	50	1: 5, 2: 5, 3: 5, 4: 5

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
YLD	0	1	49.24	32.58	9.9	20.09	34.3	80.5	126.47	▇▂▂▃▁

datos %>% 
  group_by(GEN) %>% 
  summarise(promedio = mean(YLD),
            desviacion = sd(YLD),
            N = n())

## # A tibble: 50 x 4
##    GEN   promedio desviacion     N
##    <fct>    <dbl>      <dbl> <int>
##  1 1         44.7       32.6     5
##  2 2         41.3       26.3     5
##  3 3         39.7       23.4     5
##  4 4         49.2       32.3     5
##  5 5         35.2       25.3     5
##  6 6         54.5       38.7     5
##  7 7         51.2       32.5     5
##  8 8         50.8       33.8     5
##  9 9         60.8       49.5     5
## 10 10        49.1       40.5     5
## # ... with 40 more rows

datos %>% 
  group_by(ENV) %>% 
  summarise(promedio = mean(YLD),
            desviacion = sd(YLD),
            N = n())

## # A tibble: 5 x 4
##   ENV   promedio desviacion     N
##   <fct>    <dbl>      <dbl> <int>
## 1 1         24.0       7.34    50
## 2 2         16.5       3.19    50
## 3 3         78.4       9.40    50
## 4 4         94.3      16.2     50
## 5 5         32.9       4.54    50

Visualización de Descriptivos.

datos %>% 
  ggplot(mapping = aes(x = GEN, y = YLD, color = ENV)) +
  geom_point() +
  geom_line(aes(group = ENV)) + theme(axis.text.x = element_text(angle = 30, hjust = 1))

NOTAS: - A través de los estadisticos descriptivos podemos observar que los datos están estructurados por dos fuentes de variación o categorias: ambiente (ENV) y componenete genético (GEN).Que influyen sobre la producción o variable de respuesta. - En la gráfica vemos como el ambiente 3 y 4 inside en la expresión de los genes y favorece el aumento en las medias la producción. - Podríamos decir que las fuentes de variación tienen una relación aditiva.

Modelo.

\[ y_{ij} = \mu + \gamma_i + \delta_j + \epsilon_{ij}, \qquad i = 1, 2, \cdots, 6, \quad j = 1, 2, \cdots, 4 \\ \epsilon_{ij} \sim \mathcal{N}(0,\sigma^2) \quad i.i.d. \]

Hipótesis.

ENV como principal (puede ser de verificación)

\[H_0: \gamma_1 = \gamma_2 = \gamma_3 = \gamma_4 = \gamma_5 = 0 \\ H_A: \textrm{Alguna } \gamma_i \textrm{ diferente de } 0\]

Verificación GEN (puede ser principal)

\[H_0: \delta_1 = \delta_2 = \delta_3 = \delta_4 =...\delta_{50}= 0 \\ H_A: \textrm{Alguna } \delta_i \textrm{ diferente de } 0 \]

Nivel de Significancia.

usaremos un nivel de significancia del 5% (0.05)

Modelo.

modeloENVGEN <- aov(YLD ~ ENV + GEN,
               data = datos)
summary(modeloENVGEN)

##              Df Sum Sq Mean Sq F value   Pr(>F)    
## ENV           4 242922   60731 876.774  < 2e-16 ***
## GEN          49   7830     160   2.307 2.83e-05 ***
## Residuals   196  13576      69                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Nota:

Como el valor p (2e-16) para la fuente de variación “ENV” es menor (pero por mucho) que el nivel de significancia (0.05), existe evidencia para rechazar la hipótesis nula, es decir, que al menos un par de “ambientes” poseen diferencia estadística para la variable producción (YLD).
Como el valor p (2.83e-05) para la fuente de variación “GEN” es menor que el nivel de significancia (0.05), existe evidencia para rechazar la hipótesis nula, es decir, que al menos un par de “Genes” poseen diferencia estadística para la variable producción (YLD).

Diagnóstico del Modelo.

Antes de estimar las medias realizaremos diagnóstico del Modelo:

residualesENVGEN <- residuals(modeloENVGEN)

Normalidad

shapiro.test(residualesENVGEN)

## 
##  Shapiro-Wilk normality test
## 
## data:  residualesENVGEN
## W = 0.97896, p-value = 0.0009167

library(ggpubr)
ggqqplot(residualesENVGEN)

Notas: - Como el valor p (0.0009167) es menor que el nivel de significancia (0,05) podemos decir que existe evidencia para rechazar la hipótesis nula, Es decir, los residuales del modelo no se distribuyen de forma normal. - En la gráfica observamos los desvios de algunos datos que influyen en la no normalidad de los residuales del modelo.

Homocedasticidad.

library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:dplyr':
## 
##     recode

## The following object is masked from 'package:purrr':
## 
##     some

leveneTest(residualesENVGEN ~ datos$ENV)

## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value    Pr(>F)    
## group   4  8.4299 2.183e-06 ***
##       245                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

leveneTest(residualesENVGEN ~ datos$GEN)

## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group  49  0.6269 0.9727
##       200

NOTA: - para la fuente de varianza “ENV” cuyo p(2.183e-06) es menor que la hipÓtesis nula (0,05) no existe homocedasticidad, es decir, la variabilidad de los diferentes ambiente es discrepante.

para la fuente de varianza “GEN” cuyo p(0.9727) es mayor que la hipótesis nula (0,05) existe homocedasticidad, es decir, la variabilidad de los diferentes genes no es discrepante.

Estmaciones promedio e IC 95%

library(ggeffects)
ggeffect(model = modeloENVGEN)

## $ENV
## # Predicted values of YLD
## 
## ENV | Predicted |         95% CI
## --------------------------------
## 1   |     24.05 | [21.73, 26.37]
## 2   |     16.54 | [14.22, 18.86]
## 3   |     78.43 | [76.11, 80.75]
## 4   |     94.34 | [92.02, 96.66]
## 5   |     32.86 | [30.54, 35.18]
## 
## $GEN
## # Predicted values of YLD
## 
## GEN | Predicted |         95% CI
## --------------------------------
## 1   |     44.74 | [37.40, 52.08]
## 7   |     51.21 | [43.87, 58.55]
## 13  |     57.65 | [50.31, 64.99]
## 19  |     52.12 | [44.78, 59.46]
## 26  |     44.67 | [37.33, 52.01]
## 32  |     50.19 | [42.85, 57.53]
## 38  |     53.12 | [45.78, 60.46]
## 50  |     46.11 | [38.77, 53.45]
## 
## attr(,"class")
## [1] "ggalleffects" "list"        
## attr(,"model.name")
## [1] "modeloENVGEN"

NOTA: - Como ya hemos observado en el gráfico de estadísticos predictivos, los ambientes 3 y 4 favorecen la expresióm de los genes que inciden en una mayor producción con medias de 78.43 y 94.34 respectivamente. - El intervalo de confianza con el 95% de seguridad para la media del ambiente 3 va desde 76.11 hasta 80.75. Para el ambiente 4 va de 92.02 hasta 96.66.

Comparaciones Múltiples:

library(broom)
TukeyHSD(modeloENVGEN, which = "ENV") %>% tidy()

## # A tibble: 10 x 7
##    term  contrast null.value estimate conf.low conf.high adj.p.value
##    <chr> <chr>         <dbl>    <dbl>    <dbl>     <dbl>       <dbl>
##  1 ENV   2-1               0    -7.51   -12.1      -2.92    1.09e- 4
##  2 ENV   3-1               0    54.4     49.8      59.0     7.57e-14
##  3 ENV   4-1               0    70.3     65.7      74.9     7.57e-14
##  4 ENV   5-1               0     8.81     4.23     13.4     3.15e- 6
##  5 ENV   3-2               0    61.9     57.3      66.5     7.57e-14
##  6 ENV   4-2               0    77.8     73.2      82.4     7.57e-14
##  7 ENV   5-2               0    16.3     11.7      20.9     1.22e-13
##  8 ENV   4-3               0    15.9     11.3      20.5     1.23e-13
##  9 ENV   5-3               0   -45.6    -50.2     -41.0     7.57e-14
## 10 ENV   5-4               0   -61.5    -66.1     -56.9     7.57e-14

NOTA: - Como el valor “p” es menor que el nivel de significancia (0.05), existe evidencia para rechazar la hipótesis nula. - Existen diferencias estadisticamente significativas entre todos los niveles de la fuente de variación “ambiente”

TukeyHSD(modeloENVGEN, which = "ENV") %>%
  tidy() %>% 
  ggplot(mapping = aes(x = contrast, y = estimate)) +
  geom_point() +
  geom_errorbar(mapping = aes(ymin = conf.low, ymax = conf.high),
                width = 0.2) +
  geom_hline(yintercept = 0, lty = 2, color = "orange") +
  coord_flip()

NOTA: ObservaciOnes de la comparaciones entre medias puntuales, intervalos y sus diferencias:

No se espera que exista igualdad entre, almenos, dos niveles del factor “ambiente” en futuras replicas del experimento.
Los ambientes 2-1 y 5-1 presentan las menores diferencias entre sus medias puntuales y las medias bajas y altas de sus intervalos.
El ambiente 4 saca gran diferencia a la media puntual (e intervalos) de sus pares: 1 (70.29547) y 2 (77.80087).

Ejercicio 2.

Datos: Areas cultivadas y producción Agricola en Antioquia desde 1990. Usaremos los datos registrados en el año 2020 del Rubro “Aguacate”.

library(tidyverse)
library(janitor)
datos_aguacate <- read_csv("Areas_cultivadas_y_produccion_agr_cola_en_Antioquia_desde_1990.csv") %>%  
  clean_names() %>% 
  filter(rubro == "Aguacate", 
         ano == "2020")

## 
## -- Column specification --------------------------------------------------------
## cols(
##   Tipo = col_character(),
##   Rubro = col_character(),
##   Subregion = col_character(),
##   Año = col_double(),
##   Municipio = col_character(),
##   `Área Total` = col_double(),
##   `Área Producción` = col_double(),
##   `Volumen Producción` = col_double()
## )

datos_aguacate

## # A tibble: 33 x 8
##    tipo      rubro   subregion   ano municipio        area_total area_produccion
##    <chr>     <chr>   <chr>     <dbl> <chr>                 <dbl>           <dbl>
##  1 Permanen~ Aguaca~ Suroeste   2020 Amagá                   7               3  
##  2 Permanen~ Aguaca~ Nordeste   2020 Amalfi                 63               0  
##  3 Permanen~ Aguaca~ Occidente  2020 Anzá                   14.8             9.8
##  4 Permanen~ Aguaca~ Urabá      2020 Apartadó              105             100  
##  5 Permanen~ Aguaca~ Oriente    2020 Argelia                56              44  
##  6 Permanen~ Aguaca~ Norte      2020 Briceño                54              26  
##  7 Permanen~ Aguaca~ Urabá      2020 Carepa                 27              10  
##  8 Permanen~ Aguaca~ Occidente  2020 Dabeiba                39.5            17.5
##  9 Permanen~ Aguaca~ Oriente    2020 El Carmen De Vi~       42              42  
## 10 Permanen~ Aguaca~ Occidente  2020 Frontino               15.2             9.3
## # ... with 23 more rows, and 1 more variable: volumen_produccion <dbl>

Objetivo y Modelo.

Objetivo: Predecir el volúmen de producción en función del área de producción.
Modelo: \(y_i = \beta_0 + \beta_1 \times X + \epsilon_i\), donde \(X:\) área de producción.

Gráficas

datos_aguacate %>% 
  ggplot(mapping = aes(x = area_produccion, y = volumen_produccion)) +
  geom_point() +
  geom_smooth(se = FALSE) +
  geom_smooth(method = "lm", se = FALSE, color = "red")

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

## `geom_smooth()` using formula 'y ~ x'

datos_aguacate %>% 
  ggplot(mapping = aes(x = area_produccion)) +
  geom_density()

datos_aguacate %>% 
  ggplot(mapping = aes(x = volumen_produccion)) +
  geom_density()

NOTA:

La relación entre las variables predictora(x) y de predicción(y), área_producción y Volúmen_producción, respectivamente, parece ser lineal con dirección positiva. Con pocos elementos que se desligan de esa tendencia.
Las gráficas de densidad muestran quela mayoría de las áreas de producción, para el cultivo de aguacate, van desde las 0 hasta las 90 áreas, aproximadamente. y los volúmenes de producción se ubican entre los 0 hasta los 1200 unidades de volúmen.

Hipótesis

\[H_0: \rho = 0 \\ H_A: \rho \neq 0 \]

Nivel de significancia

Nivil de significancia del 5% (0.05)

Normalidad

Area de producción

shapiro.test(datos_aguacate$area_produccion)

## 
##  Shapiro-Wilk normality test
## 
## data:  datos_aguacate$area_produccion
## W = 0.48999, p-value = 1.403e-09

NOTA: Como el valor de p (1.403e-09) es menor que el nivel de significancia, existe evidencia para rechazar la hipótesis nula, es decir, la variable “área_de_producción” no posee distribución normal.

volúmen de producción.

shapiro.test(datos_aguacate$volumen_produccion)

## 
##  Shapiro-Wilk normality test
## 
## data:  datos_aguacate$volumen_produccion
## W = 0.37889, p-value = 1.059e-10

NOTA: Como el valor de p (1.059e-10) es menor que el nivel de significancia, existe evidencia para rechazar la hipótesis nula, es decir, la variable “Volúmen_de_producción” no posee distribución normal.

Correlación

Para determinar la correlación usaremos el método No paramétrico de spearman.

cor.test(x = datos_aguacate$area_produccion, y = datos_aguacate$volumen_produccion, method = "spearman",
         conf.level = 0.95)

## Warning in cor.test.default(x = datos_aguacate$area_produccion, y =
## datos_aguacate$volumen_produccion, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  datos_aguacate$area_produccion and datos_aguacate$volumen_produccion
## S = 474.7, p-value = 3.319e-14
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.9206721

NOTA: - Como el valor p (3.319e-14) es menor que el nivel de significancia, existe evidencia para rechazar la hipótesis nula, es decir,la correlacion entre las variables no es igual a cero. - El valor de \(\rho= 0.9206721\) indica una correlación lineal positiva fuerte. - Por su significancia estadística, existe una probabilidad alta de volver a encontrar estos parámetros en otros análisis experimentales.

Modelo.

modelo_regresion <- lm(volumen_produccion ~ area_produccion, data = datos_aguacate)
summary(modelo_regresion)

## 
## Call:
## lm(formula = volumen_produccion ~ area_produccion, data = datos_aguacate)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3388.7   -46.7    36.2   126.6  2838.3 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -150.551    166.004  -0.907    0.371    
## area_produccion   13.865      1.101  12.590    1e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 855.7 on 31 degrees of freedom
## Multiple R-squared:  0.8364, Adjusted R-squared:  0.8311 
## F-statistic: 158.5 on 1 and 31 DF,  p-value: 1.004e-13

NOTA: - Como el valor de p (1e-13) es menor al valor de significancia (0,05), rechazamos la hipótesis nula, es decir, la pendiente es diferente de cero. (estadísticamente significativo) - El estimativo para la pendiente nos indica que por cada unidad de área productiva que se aumente, la producción aumentará, en promedio, a 13.865 unidades de volúmen. - Del 100% de la variabilidad en el volúmen de la producción, aproximadamente, el 83% se debe al cambio en el área de producción.

Residuales.

residuales <- residuals(modelo_regresion)

Normalidad.

shapiro.test(residuales)

## 
##  Shapiro-Wilk normality test
## 
## data:  residuales
## W = 0.62223, p-value = 5.238e-08

NOTA:

Como el valor p (5.238e-08) es menor que el nivel de significancia (0,05) podemos argumentar que existe evidencia para rechazar la hipótesis nula, Es decir, los residuales del modelo no se distribuyen de forma normal.

Homocedasticidad.

library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

bptest(modelo_regresion)

## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_regresion
## BP = 26.351, df = 1, p-value = 2.847e-07

NOTA:

Como el valor p (2.847e-07) es menor que el nivel de significancia (0,05) podemos argumentar que existe evidencia para rechazar la hipótesis nula, Es decir, nuestro modelo no es homocedástico.
Necesitamos ajustar el modelo.

Predicción

Puntual:

predict(object = modelo_regresion,
        newdata = data.frame(area_produccion = c(150.0, 300)))

##        1        2 
## 1929.209 4008.970

NOTA:

Para un área_produccion de \(150.0\), nuestro modelo predice una producción, en promedio, de \(1929.209\) unidades de volúmen

Intervalo:

predict(object = modelo_regresion,
        newdata = data.frame(area_produccion = 150.0), interval = "confidence")

##        fit      lwr      upr
## 1 1929.209 1572.206 2286.213

NOTA:

Para un área_produccion de \(150.0\), nuestro modelo predice una producción, en promedio, de \(1929.209\) unidades de volúmen. Pero, ademas, podemos obtener un intervalo de confianza en promedios que van desde \(1572.206\) hasta \(2286.213\) unidades de volúmen de producción.

Ejercicio 3.

Suponga que va a planificar un experimento en un tema de alta relevancia para usted, puede ser cualquier caso hipotético (cultivos, animales, agroindustria, etc.), informe lo siguiente:

¿Cuál o cuáles son las fuentes de variación que va a evaluar en su experimento?.
¿Cuáles niveles tienen? Reporte cuál es su variable respuesta y tentativamente en qué unidades estaría medida.
¿Quién es su unidad experimental?
¿Cómo va a realizar las mediciones?
¿Cuántas repeticiones tendría para cada tratamiento?
¿Existen factores de ruido que puedan ser bloqueados?

Pregunta a responder.

¿Cómo afecta a la producción, el tipo de clón de cacao, el tipo de suelo y la altitud?

Fuentes de Variación.

Clón de Cacao, Tipo de suelo y altitud

Niveles:

Clon Cacao: L46-H57, L46-H75, L46-H88 y CCN-51.
Tipo de Suelo: 1, 2 y 3
Altitud: 100, 600, 1200 msnm

Variable respuesta.

Variable “y” o de respuesta: Producción en kg

Unidad experimental.

39 arbustos de cada uno de los clones en fase productiva.En un área 625 m2
156 arbustos * 3 tipos de suelo = 468 arbustos.
468 arbustos* 3 altitudes = 1404 arbustos

Mediciones.

Se entandarizarán:

Estado fisiológico de los arbustos desde Plantación, en área de establecimiento final, hasta etapa productiva.
prácticas de fertilización y fitosanitarias.
Labores culturales: podas, cosecha, post cosecha.

Luego se medirá:

producción total por clon, tipo de suelo y altitud.

Repeticiones.

Serán 36 repeticiones.

Bloqueo

Elementos del clima: Precipitación, temperatura, viento. (Relacionado con la altitud)

NOTA:El modelo del diseño experimental será uno factorial de mas de dos vías o, en su defecto, un Diseño de Bloque Completamente Aleatorizado, según la relación que presenten las fuentes de variación en el análisis de los datos.

Diseno_Experimento_Regresion_Lineal

Estefani_Cardenas_Victor_Palacio

2021-11-05

Ejercicio 1

Estadísticos Descriptivos.

Visualización de Descriptivos.

Modelo.

Hipótesis.

ENV como principal (puede ser de verificación)

Verificación GEN (puede ser principal)

Nivel de Significancia.

Modelo.

Diagnóstico del Modelo.

Normalidad

Homocedasticidad.

Estmaciones promedio e IC 95%

Comparaciones Múltiples:

Ejercicio 2.

Objetivo y Modelo.

Gráficas

Hipótesis

Nivel de significancia

Normalidad

Correlación

Modelo.

Residuales.

Normalidad.

Homocedasticidad.

Predicción

Puntual:

Intervalo:

Ejercicio 3.

Pregunta a responder.

Fuentes de Variación.

Niveles:

Variable respuesta.

Unidad experimental.

Mediciones.

Repeticiones.

Bloqueo