set.seed(123)
diam_geom = c(
rnorm (4, 1.8, 0.1),
rnorm (4, 2.0, 0.12),
rnorm (4, 1.9, 0.09)
)
gen = gl (3, 4, 12, paste0 ('g_', 1:3))
procedencia = gl (4,1,12, paste0('l_', 1:4))
data= data.frame(gen, procedencia, diam_geom)
head(data)
## gen procedencia diam_geom
## 1 g_1 l_1 1.743952
## 2 g_1 l_2 1.776982
## 3 g_1 l_3 1.955871
## 4 g_1 l_4 1.807051
## 5 g_2 l_1 2.015515
## 6 g_2 l_2 2.205808
library(collapsibleTree)
## Warning: package 'collapsibleTree' was built under R version 4.2.3
collapsibleTreeSummary (data,
c('procedencia',
'gen',
'diam_geom'),
collapsed = FALSE)
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.2.2
ggplot(data)+
aes (gen, diam_geom)+
geom_point (size=7,
color = 'yellow')+
facet_wrap(~procedencia)+
theme_dark()
\[H_0: \mu_{g_1}=\mu_{g2}= \mu_{g_3}\]
mod = aov (diam_geom ~ procedencia + gen, data)
summary(mod)
## Df Sum Sq Mean Sq F value Pr(>F)
## procedencia 3 0.04373 0.01458 1.334 0.3483
## gen 2 0.08908 0.04454 4.078 0.0762 .
## Residuals 6 0.06554 0.01092
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Advertencia No revisar el p valor de los bloques, es decir la procedencia. Solo interpretar el de los factores, es decir genes.
Como el p valor de los genes es mayor al 5% (0.007*100= 7%>5%), no rechaza la hipotesis nula, es decir que los fenotipos son iguales, el diametro es el mismo.
¿Valió la pena preocuparse por la procedencia de las semillas? En el análsis inferencial H= a Fvalue H = 1.334 (F-value de los bloques), cuando H>1,sugiere que si valió la pena bloquear
En ocasiones cuando el H es menor a 1, si vale la pena bloquear (artículo Carlos y profe), porque el bloqueo suele resolver problemas que no están a la vista.
Pero quitar el bloqueo porque no “sirve” y dejar el análisis sin bloqueo está mal, se deben dejar así sean menores que 1 y la conclusión al final se mantenga.
#Paso 1. Extraer residuales del modelo.
res_mod1 = mod$residuals
#2. Probando normalidad
shapiro.test(res_mod1)
##
## Shapiro-Wilk normality test
##
## data: res_mod1
## W = 0.97765, p-value = 0.9725
# Paso 3. Probando homocedasticidad (varianzas iguales)
bartlett.test(res_mod1, data$gen)
##
## Bartlett test of homogeneity of variances
##
## data: res_mod1 and data$gen
## Bartlett's K-squared = 1.709, df = 2, p-value = 0.4255
Si los patrones de los residuales son exponenciales, en forma de cono, puede ser peligroso, lo ideal es que sean constantes
En ambos casos los p-value de las pruebas fueron >5%
plot(data$diam_geom,
res_mod1,
pch=16)
Cuando los resultados de los residuales tiene patrón, se vuelve cuestionable el análisis, en este caso, no existe ese patrón, se ve algo constante.
Valió la pena bloquear, no hay diferencias en el genotipo, se cumplieron todos los supuestos.
set.seed(123)
diam_geom = c(
rnorm (20, 1.8, 0.1),
rnorm (20, 2.0, 0.12),
rnorm (20, 1.9, 0.09)
)
gen = gl (3, 20, 60, paste0 ('g_', 1:3))
procedencia = gl (4,5,60, paste0('l_', 1:4))
data= data.frame(gen, procedencia, diam_geom)
head(data)
## gen procedencia diam_geom
## 1 g_1 l_1 1.743952
## 2 g_1 l_1 1.776982
## 3 g_1 l_1 1.955871
## 4 g_1 l_1 1.807051
## 5 g_1 l_1 1.812929
## 6 g_1 l_2 1.971506
library(collapsibleTree)
collapsibleTreeSummary (data,
c('procedencia',
'gen',
'diam_geom'),
collapsed = FALSE)
library(lattice)
bwplot (diam_geom ~ gen | procedencia,
data)
mod3= aov (diam_geom ~ procedencia * gen, data)
summary(mod3)
## Df Sum Sq Mean Sq F value Pr(>F)
## procedencia 3 0.0310 0.01034 1.135 0.344
## gen 2 0.3233 0.16164 17.732 1.71e-06 ***
## procedencia:gen 6 0.0407 0.00678 0.744 0.617
## Residuals 48 0.4376 0.00912
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Usando * es porque hay repeticiones, + cuando no#
Primero ver la interacción=(procedencia:gen),da mayor al 5%, es decir que no hay interacción, si se cumple esto se puede ver el genotipo, si no, se deben ir a los gráficos (el de las cajas).