Anova no parametrico
kruskal.test(dietas~trats)
##
## Kruskal-Wallis rank sum test
##
## data: dietas by trats
## Kruskal-Wallis chi-squared = 2.34, df = 2, p-value = 0.3104
data(InsectSprays)
attach(InsectSprays)
str(InsectSprays) #Permite ver las características de las variables
## 'data.frame': 72 obs. of 2 variables:
## $ count: num 10 7 20 14 14 12 10 23 17 20 ...
## $ spray: Factor w/ 6 levels "A","B","C","D",..: 1 1 1 1 1 1 1 1 1 1 ...
tapply(count, spray, mean) #Permite ver la media de los factores (Tratamientos)
## A B C D E F
## 14.500000 15.333333 2.083333 4.916667 3.500000 16.666667
tapply(count, spray, sd) #Permite ver la desviación estándar de los factores.
## A B C D E F
## 4.719399 4.271115 1.975225 2.503028 1.732051 6.213378
tapply(count, spray, length) #Permite ver la cantidad de datos.
## A B C D E F
## 12 12 12 12 12 12
andeva2<-aov(count~spray, data=InsectSprays)
summary(andeva2)#no normales
## Df Sum Sq Mean Sq F value Pr(>F)
## spray 5 2669 533.8 34.7 <2e-16 ***
## Residuals 66 1015 15.4
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
shapiro.test(residuals(andeva2))
##
## Shapiro-Wilk normality test
##
## data: residuals(andeva2)
## W = 0.96006, p-value = 0.02226
#Homogeneidad de varianzas
library(car)
leveneTest(andeva2)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 3.8214 0.004223 **
## 66
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
kruskal.test(count~spray, data=InsectSprays)
##
## Kruskal-Wallis rank sum test
##
## data: count by spray
## Kruskal-Wallis chi-squared = 54.691, df = 5, p-value = 1.511e-10
#Post-hoc (No paramétricas)
#LSD test (Least Significant Difference) (Comparación Planeada)
pairwise.wilcox.test(count,spray, p.adj = "none", exact=F)
##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: count and spray
##
## A B C D E
## B 0.5812 - - - -
## C 3.8e-05 3.8e-05 - - -
## D 7.8e-05 6.9e-05 0.0027 - -
## E 3.4e-05 3.4e-05 0.0526 0.1744 -
## F 0.4342 0.9078 3.4e-05 7.0e-05 3.4e-05
##
## P value adjustment method: none
#Bonferroni (recomendado para tamaño de muestras iguales, medias entre trat. diferentes)
pairwise.wilcox.test(count,spray, p.adj = "bonf",exact=F)
##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: count and spray
##
## A B C D E
## B 1.00000 - - - -
## C 0.00058 0.00058 - - -
## D 0.00117 0.00104 0.03977 - -
## E 0.00051 0.00051 0.78860 1.00000 -
## F 1.00000 1.00000 0.00052 0.00105 0.00052
##
## P value adjustment method: bonferroni
ta <- c(40, 48, 54, 42, 55, 38, 42, 49, 41, 50)
tb <- c(22, 24, 23, 29, 29, 27, 27, 29, 33, 30)
tc <- c(24, 23, 24, 36, 27, 22, 31, 33, 34, 29)
tctl <- c(39, 34, 39, 40, 38, 36, 36, 40, 39, 37)
crecimiento <- c(ta,tb,tc,tctl) #permite agrupar el crecimiento en un solo vector
crecimiento
## [1] 40 48 54 42 55 38 42 49 41 50 22 24 23 29 29 27 27 29 33 30 24 23 24 36 27
## [26] 22 31 33 34 29 39 34 39 40 38 36 36 40 39 37
trats <- gl(4,10,label=c("ta","tb","tc", "tctl")) #asocia los valores de crecimiento al los tratamientos
length(trats)==length(crecimiento) #confirmamos que tanto las variables independientes como los factores tengan la misma longitud, en caso no de tener la misma longitud, tendremos que utilizar NA, para el valor faltante.
## [1] TRUE
is.factor(trats) #Confirmamos que el vector es factor
## [1] TRUE
#Ejecutamos el ANOVA
anova3 <- aov(crecimiento~trats)
#Resumen
summary(anova3)
## Df Sum Sq Mean Sq F value Pr(>F)
## trats 3 2307.1 769.0 39.51 1.76e-11 ***
## Residuals 36 700.7 19.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Pueba Post Hoc test
TukeyHSD(anova3)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = crecimiento ~ trats)
##
## $trats
## diff lwr upr p adj
## tb-ta -18.6 -23.91377 -13.28623 0.0000000
## tc-ta -17.6 -22.91377 -12.28623 0.0000000
## tctl-ta -8.1 -13.41377 -2.78623 0.0012165
## tc-tb 1.0 -4.31377 6.31377 0.9569349
## tctl-tb 10.5 5.18623 15.81377 0.0000322
## tctl-tc 9.5 4.18623 14.81377 0.0001498
ANOVA FACTORIAL
data(ToothGrowth)
str(ToothGrowth) #Tenemos dos variables numericas y un factor
## 'data.frame': 60 obs. of 3 variables:
## $ len : num 4.2 11.5 7.3 5.8 6.4 10 11.2 11.2 5.2 7 ...
## $ supp: Factor w/ 2 levels "OJ","VC": 2 2 2 2 2 2 2 2 2 2 ...
## $ dose: num 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ...
head(ToothGrowth)
## len supp dose
## 1 4.2 VC 0.5
## 2 11.5 VC 0.5
## 3 7.3 VC 0.5
## 4 5.8 VC 0.5
## 5 6.4 VC 0.5
## 6 10.0 VC 0.5
andeva3<-aov(len ~ supp + dose, data=ToothGrowth)
summary(andeva3)
## Df Sum Sq Mean Sq F value Pr(>F)
## supp 1 205.4 205.4 11.45 0.0013 **
## dose 1 2224.3 2224.3 123.99 6.31e-16 ***
## Residuals 57 1022.6 17.9
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
andeva4<-aov(len ~ supp * dose, data=ToothGrowth) #Observar que se utiliza el "*", para encontrar la interaccion entre "supp"" y "dose".
summary(andeva4)
## Df Sum Sq Mean Sq F value Pr(>F)
## supp 1 205.4 205.4 12.317 0.000894 ***
## dose 1 2224.3 2224.3 133.415 < 2e-16 ***
## supp:dose 1 88.9 88.9 5.333 0.024631 *
## Residuals 56 933.6 16.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
etchrate <- read.table("https://raw.githubusercontent.com/osoramirez/bases_datos/master/etchrate.csv", sep="", header=TRUE)
etchrate
## RF run rate
## 1 160 1 575
## 2 160 2 542
## 3 160 3 530
## 4 160 4 539
## 5 160 5 570
## 6 180 1 565
## 7 180 2 593
## 8 180 3 590
## 9 180 4 579
## 10 180 5 610
## 11 200 1 600
## 12 200 2 651
## 13 200 3 610
## 14 200 4 637
## 15 200 5 629
## 16 220 1 725
## 17 220 2 700
## 18 220 3 715
## 19 220 4 685
## 20 220 5 710
head(etchrate)#muestra el encabezado
## RF run rate
## 1 160 1 575
## 2 160 2 542
## 3 160 3 530
## 4 160 4 539
## 5 160 5 570
## 6 180 1 565
str(etchrate) #observamos la caracteristica de cada variable. Note que no hay variables categóricas o factores
## 'data.frame': 20 obs. of 3 variables:
## $ RF : int 160 160 160 160 160 180 180 180 180 180 ...
## $ run : int 1 2 3 4 5 1 2 3 4 5 ...
## $ rate: int 575 542 530 539 570 565 593 590 579 610 ...
attach(etchrate)
grp.means<-tapply(rate,RF,mean)
grp.means
## 160 180 200 220
## 551.2 587.4 625.4 707.0
grp.sd<-tapply(rate,RF,sd)
grp.sd
## 160 180 200 220
## 20.01749 16.74216 20.52559 15.24795
grp.var<-tapply(rate,RF,var)
grp.var
## 160 180 200 220
## 400.7 280.3 421.3 232.5
# Convertimos cada variable a factor
etchrate$RF <- as.factor(etchrate$RF)
etchrate$run <- as.factor(etchrate$run)
str(etchrate) #Note que ahora nuestras variables de "RF" y "run" son factores
## 'data.frame': 20 obs. of 3 variables:
## $ RF : Factor w/ 4 levels "160","180","200",..: 1 1 1 1 1 2 2 2 2 2 ...
## $ run : Factor w/ 5 levels "1","2","3","4",..: 1 2 3 4 5 1 2 3 4 5 ...
## $ rate: int 575 542 530 539 570 565 593 590 579 610 ...
# Corremos modelo del anova
etch.rate.aov <- aov(rate~RF,data=etchrate)
summary(etch.rate.aov)
## Df Sum Sq Mean Sq F value Pr(>F)
## RF 3 66871 22290 66.8 2.88e-09 ***
## Residuals 16 5339 334
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Exploramos visualmente
boxplot(rate~RF, col=456,xlab="RF Power (W)", ylab=expression(paste("Observed Etch Rate (",ring(A),"/min)")))
# Media global
(erate.mean <- mean(etchrate$rate))
## [1] 617.75
# desviación estándar global
(erate.sd <- sd(etchrate$rate))
## [1] 61.6483
#Diseño unifactorial de multiples niveles
bacterias<-read.table("https://wpd.ugr.es/~bioestad/wp-content/uploads/supuesto2.txt", header = TRUE)
bacterias
## Calidad Fabricante
## 1 120 1
## 2 240 2
## 3 240 3
## 4 300 4
## 5 300 5
## 6 240 1
## 7 360 2
## 8 270 3
## 9 240 4
## 10 360 5
## 11 300 1
## 12 180 2
## 13 300 3
## 14 300 4
## 15 240 5
## 16 360 1
## 17 180 2
## 18 360 3
## 19 360 4
## 20 360 5
## 21 240 1
## 22 300 2
## 23 360 3
## 24 360 4
## 25 360 5
## 26 180 1
## 27 240 2
## 28 300 3
## 29 360 4
## 30 360 5
## 31 144 1
## 32 360 2
## 33 360 3
## 34 360 4
## 35 360 5
## 36 300 1
## 37 360 2
## 38 360 3
## 39 360 4
## 40 300 5
## 41 240 1
## 42 360 2
## 43 300 3
## 44 300 4
## 45 360 5
#Debemos transformar la variable referente a los niveles del factor fijo como factor para poder hacer los cálculos de forma adecuada
bacterias$Fabricante<- factor(bacterias$Fabricante)
bacterias$Fabricante
## [1] 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3
## [39] 4 5 1 2 3 4 5
## Levels: 1 2 3 4 5
#Para calcular la tabla ANOVA primero hacemos uso de la función “aov” de la siguiente forma:
mod <- aov(Calidad ~ Fabricante, data = bacterias)
summary(mod) # TABLA ANOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## Fabricante 4 57363 14341 3.976 0.00827 **
## Residuals 40 144272 3607
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Varianza dentro de los factores (Varianza residual) #3607
#Varianza entre los factores # = 1192.667
#Varianza Total= 3607+1192.667=4799.667
## Diseño en Bloques Completos Aleatorizados
abetos <- read.table("https://wpd.ugr.es/~bioestad/wp-content/uploads/supuesto3.txt", header=TRUE)
abetos
## y Tratamiento Abeto
## 1 7 1 1
## 2 9 2 1
## 3 10 3 1
## 4 8 1 2
## 5 9 2 2
## 6 10 3 2
## 7 9 1 3
## 8 9 2 3
## 9 12 3 3
## 10 10 1 4
## 11 9 2 4
## 12 12 3 4
## 13 11 1 5
## 14 12 2 5
## 15 14 3 5
## 16 8 1 6
## 17 10 2 6
## 18 9 3 6
## 19 7 1 7
## 20 8 2 7
## 21 7 3 7
## 22 8 1 8
## 23 8 2 8
## 24 7 3 8
## 25 7 1 9
## 26 9 2 9
## 27 10 3 9
## 28 8 1 10
## 29 9 2 10
## 30 10 3 10
#al no ser directamente motivo de estudio, los abetos es un factor secundario que recibe el nombre de bloque.
#Variable respuesta: Número de semillas
#Factor: Tratamiento que tiene tres niveles. Es un factor de efectos fijos ya que viene decidido qué niveles concretos se van a utilizar.
#Bloque: Abeto que tiene diez niveles. Es un factor de efectos fijos ya que viene decidido qué niveles concretos se van a utilizar.
#Modelo completo: Los tres tratamientos se prueban en cada bloque exactamente una vez.
#Tamaño del experimento: Número total de observaciones (30).
# las observaciones en una sola columna y a continuación especificado su tratamiento y su bloque correspondiente.
abetos$Tratamiento = factor(abetos$Tratamiento)
abetos$Tratamiento
## [1] 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3
## Levels: 1 2 3
abetos$Abeto = factor(abetos$Abeto)
abetos$Abeto
## [1] 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 6 6 6 7 7 7 8 8 8 9
## [26] 9 9 10 10 10
## Levels: 1 2 3 4 5 6 7 8 9 10
mod = aov(y ~ Tratamiento + Abeto, data = abetos)
summary(mod) # TABLA ANOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## Tratamiento 2 16.2 8.100 9.228 0.00174 **
## Abeto 9 54.8 6.089 6.937 0.00026 ***
## Residuals 18 15.8 0.878
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# La eficacia de este diseño depende de los efectos de los bloques, Un valor grande de F de los bloques (6.937) implica que el factor bloque tiene un efecto grande
#¿Cómo saber cuándo se puede prescindir de los bloques?
#La respuesta la tenemos en el valor de la F experimental de los bloques,
#se ha comprobado que si dicho valor es mayor que 3, no conviene prescindir de los bloques para efectuar los contrastes.
#La interacción entre el factor bloque y los tratamientos vamos a estudiarla analíticamente mediante el Test de Interacción de un grado de Tukey
#Aditividad entre los bloques y tratamientos
library(daewr)
## Registered S3 method overwritten by 'DoE.base':
## method from
## factorize.factor conf.design
Tukey1df(abetos)
## Source df SS MS F Pr>F
## A 2 16.2 8.1
## B 9 54.8 6.0889
## Error 18 15.8 71.1
## NonAdditivity 1 3.5573 3.5573 4.94 0.0401
## Residual 17 12.2427 0.7202
#Puesto que el p-valor (Pr > F) es 1 no rechazamos la hipótesis nula de no interacción,
#es decir, no hay interacción entre los tratamientos aplicados y los abetos
#Normalidad
shapiro.test(mod$residuals) #0.3935 normales
##
## Shapiro-Wilk normality test
##
## data: mod$residuals
## W = 0.96415, p-value = 0.3935
qqnorm(mod$residuals)
#Homogeneidad de Varianzas
bartlett.test(abetos$y, abetos$Tratamiento)
##
## Bartlett test of homogeneity of variances
##
## data: abetos$y and abetos$Tratamiento
## Bartlett's K-squared = 4.1729, df = 2, p-value = 0.1241
bartlett.test(abetos$y, abetos$Abeto)
##
## Bartlett test of homogeneity of variances
##
## data: abetos$y and abetos$Abeto
## Bartlett's K-squared = 4.0723, df = 9, p-value = 0.9066
#Independencia
layout(matrix(c(1,2,3,4),2,2))
plot(mod)
#Post hoc
library(agricolae)
duncan=duncan.test(mod, "Tratamiento", group = T)
duncan
## $statistics
## MSerror Df Mean CV
## 0.8777778 18 9.2 10.18367
##
## $parameters
## test name.t ntr alpha
## Duncan Tratamiento 3 0.05
##
## $duncan
## Table CriticalRange
## 2 2.971152 0.8802727
## 3 3.117384 0.9235973
##
## $means
## y std r Min Max Q25 Q50 Q75
## 1 8.3 1.337494 10 7 11 7.25 8 8.75
## 2 9.2 1.135292 10 8 12 9.00 9 9.00
## 3 10.1 2.183270 10 7 14 9.25 10 11.50
##
## $comparison
## NULL
##
## $groups
## y groups
## 3 10.1 a
## 2 9.2 b
## 1 8.3 c
##
## attr(,"class")
## [1] "group"
#En el apartado “$groups” concluimos que los tres tratamientos difieren significativamente entre sí.
#Se observa que la concentración media del número de semillas es mayor con el Tratamiento3 (10.1) y menor con el Tratamiento1 (8.3).
duncan=duncan.test(mod, "Abeto", group = T)
duncan
## $statistics
## MSerror Df Mean CV
## 0.8777778 18 9.2 10.18367
##
## $parameters
## test name.t ntr alpha
## Duncan Abeto 10 0.05
##
## $duncan
## Table CriticalRange
## 2 2.971152 1.607151
## 3 3.117384 1.686250
## 4 3.209655 1.736161
## 5 3.273593 1.770746
## 6 3.320327 1.796026
## 7 3.355651 1.815133
## 8 3.382941 1.829895
## 9 3.404326 1.841462
## 10 3.421226 1.850604
##
## $means
## y std r Min Max Q25 Q50 Q75
## 1 8.666667 1.5275252 3 7 10 8.0 9 9.5
## 10 9.000000 1.0000000 3 8 10 8.5 9 9.5
## 2 9.000000 1.0000000 3 8 10 8.5 9 9.5
## 3 10.000000 1.7320508 3 9 12 9.0 9 10.5
## 4 10.333333 1.5275252 3 9 12 9.5 10 11.0
## 5 12.333333 1.5275252 3 11 14 11.5 12 13.0
## 6 9.000000 1.0000000 3 8 10 8.5 9 9.5
## 7 7.333333 0.5773503 3 7 8 7.0 7 7.5
## 8 7.666667 0.5773503 3 7 8 7.5 8 8.0
## 9 8.666667 1.5275252 3 7 10 8.0 9 9.5
##
## $comparison
## NULL
##
## $groups
## y groups
## 5 12.333333 a
## 4 10.333333 b
## 3 10.000000 b
## 10 9.000000 bc
## 2 9.000000 bc
## 6 9.000000 bc
## 1 8.666667 bc
## 9 8.666667 bc
## 8 7.666667 c
## 7 7.333333 c
##
## attr(,"class")
## [1] "group"
## Diseño en Cuadrados Latinos
#Hemos estudiado en el apartado anterior que los diseños en bloques completos aleatorizados
#utilizan un factor de control o variable de bloque con objeto de eliminar su influencia en
#la variable respuesta y así reducir el error experimental. Los diseños en cuadrados latinos
#utilizan dos variables de bloque para reducir el error experimental.
#Características:
#Se controlan tres fuentes de variabilidad, un factor principal y dos factores de bloque.
#Cada uno de los factores tiene el mismo número de niveles, K .
#Cada nivel del factor principal aparece una vez en cada fila y una vez en cada columna.
#No hay interacción entre los factores.
#Variable respuesta: Rendimiento
#Factor: Tiempo de reposo que tiene seis niveles. Es un factor de efectos fijos ya que viene decidido que niveles concretos se van a utilizar.
#Bloques: Lotes y Concentraciones, ambos con seis niveles y ambos son factores de efectos fijos.
#Tamaño del experimento: Número total de observaciones (36).
latino <- read.table("https://wpd.ugr.es/~bioestad/wp-content/uploads/supuesto5.txt", header = TRUE)
latino
## Observaciones Lote Concentraciones Tiempo_de_reposo
## 1 12 Lote1 1 A
## 2 24 Lote1 2 B
## 3 10 Lote1 3 C
## 4 18 Lote1 4 D
## 5 21 Lote1 5 E
## 6 18 Lote1 6 F
## 7 21 Lote2 1 B
## 8 26 Lote2 2 C
## 9 24 Lote2 3 D
## 10 16 Lote2 4 E
## 11 20 Lote2 5 F
## 12 21 Lote2 6 A
## 13 20 Lote3 1 C
## 14 16 Lote3 2 D
## 15 19 Lote3 3 E
## 16 18 Lote3 4 F
## 17 16 Lote3 5 A
## 18 19 Lote3 6 B
## 19 22 Lote4 1 D
## 20 15 Lote4 2 E
## 21 14 Lote4 3 F
## 22 19 Lote4 4 A
## 23 27 Lote4 5 B
## 24 17 Lote4 6 C
## 25 15 Lote5 1 E
## 26 13 Lote5 2 F
## 27 17 Lote5 3 A
## 28 25 Lote5 4 B
## 29 21 Lote5 5 C
## 30 22 Lote5 6 D
## 31 17 Lote6 1 F
## 32 11 Lote6 2 A
## 33 12 Lote6 3 B
## 34 22 Lote6 4 C
## 35 14 Lote6 5 D
## 36 20 Lote6 6 E
latino$Lote <- factor(latino$Lote)
latino$Concentraciones <- factor(latino$Concentraciones)
latino$Tiempo_de_reposo <- factor(latino$Tiempo_de_reposo)
mod1<-aov(formula = Observaciones ~ Lote + Concentraciones + Tiempo_de_reposo,
data = latino)
summary(mod1)
## Df Sum Sq Mean Sq F value Pr(>F)
## Lote 5 99.6 19.91 1.149 0.368
## Concentraciones 5 70.6 14.11 0.814 0.553
## Tiempo_de_reposo 5 117.9 23.58 1.361 0.281
## Residuals 20 346.6 17.33
#Observando los valores de los p-valores, 0.281, 0.368 y 0.553; mayores respectivamente que el nivel de significación del 5%, deducimos que ningún efecto es significativo.
## Diseños Factoriales
#considerar dos o más factores y estudiar el efecto conjunto
#que dichos factores producen sobre la variable respuesta.
#Variable respuesta: Tiempo de supervivencia
#Factor: Tipo de veneno que tiene tres niveles. Es un factor de efectos fijos ya que viene decidido qué niveles concretos se van a utilizar.
#Factor: Tipo de antídoto que tiene cuatro niveles. Es un factor de efectos fijos ya que viene decidido qué niveles concretos se van a utilizar.
#Tamaño del experimento: Número total de observaciones (12).
factorial <- read.table("https://wpd.ugr.es/~bioestad/wp-content/uploads/supuesto8.txt", header=TRUE)
factorial
## Tiempo Veneno Antídoto
## 1 4.5 1 1
## 2 2.9 2 1
## 3 2.1 3 1
## 4 11.0 1 2
## 5 6.1 2 2
## 6 3.7 3 2
## 7 4.5 1 3
## 8 3.5 2 3
## 9 2.5 3 3
## 10 7.1 1 4
## 11 10.2 2 4
## 12 3.6 3 4
factorial$Antídoto<-factor(factorial$Antídoto)
factorial$Veneno<-factor(factorial$Veneno)
mod <- aov(Tiempo~ Veneno + Antídoto , data = factorial )
mod
## Call:
## aov(formula = Tiempo ~ Veneno + Antídoto, data = factorial)
##
## Terms:
## Veneno Antídoto Residuals
## Sum of Squares 30.58667 39.40917 23.89333
## Deg. of Freedom 2 3 6
##
## Residual standard error: 1.995551
## Estimated effects may be unbalanced
summary(mod)
## Df Sum Sq Mean Sq F value Pr(>F)
## Veneno 2 30.59 15.293 3.840 0.0844 .
## Antídoto 3 39.41 13.136 3.299 0.0995 .
## Residuals 6 23.89 3.982
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#El modelo con replicación
factorial <- read.table("https://wpd.ugr.es/~bioestad/wp-content/uploads/supuesto9.txt", header=TRUE)
factorial
## Tiempo Veneno Antidoto
## 1 4.5 1 1
## 2 4.3 1 1
## 3 2.9 2 1
## 4 2.3 2 1
## 5 2.1 3 1
## 6 2.3 3 1
## 7 11.0 1 2
## 8 7.2 1 2
## 9 6.1 2 2
## 10 12.4 2 2
## 11 3.7 3 2
## 12 2.9 3 2
## 13 4.5 1 3
## 14 7.6 1 3
## 15 3.5 2 3
## 16 4.0 2 3
## 17 2.5 3 3
## 18 2.2 3 3
## 19 7.1 1 4
## 20 6.2 1 4
## 21 10.2 2 4
## 22 3.8 2 4
## 23 3.6 3 4
## 24 3.3 3 4
factorial$Veneno <- factor(factorial$Veneno)
factorial$Antidoto<-factor(factorial$Antidoto)
mod <- aov(Tiempo~ Veneno * Antidoto , data = factorial )
mod
## Call:
## aov(formula = Tiempo ~ Veneno * Antidoto, data = factorial)
##
## Terms:
## Veneno Antidoto Veneno:Antidoto Residuals
## Sum of Squares 60.44333 60.26167 20.36333 53.51000
## Deg. of Freedom 2 3 6 12
##
## Residual standard error: 2.111674
## Estimated effects may be unbalanced
summary(mod)
## Df Sum Sq Mean Sq F value Pr(>F)
## Veneno 2 60.44 30.222 6.777 0.0107 *
## Antidoto 3 60.26 20.087 4.505 0.0245 *
## Veneno:Antidoto 6 20.36 3.394 0.761 0.6138
## Residuals 12 53.51 4.459
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Por lo tanto la interacción entre ambos factores no es significativa y debemos eliminarla del modelo.
#Construimos de nuevo la Tabla ANOVA en la que sólo figurarán los efectos principales
mod <- aov(Tiempo~ Veneno + Antidoto , data = factorial )
mod
## Call:
## aov(formula = Tiempo ~ Veneno + Antidoto, data = factorial)
##
## Terms:
## Veneno Antidoto Residuals
## Sum of Squares 60.44333 60.26167 73.87333
## Deg. of Freedom 2 3 18
##
## Residual standard error: 2.025851
## Estimated effects may be unbalanced
summary(mod)
## Df Sum Sq Mean Sq F value Pr(>F)
## Veneno 2 60.44 30.222 7.364 0.0046 **
## Antidoto 3 60.26 20.087 4.894 0.0117 *
## Residuals 18 73.87 4.104
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
factorial <- read.table("https://wpd.ugr.es/~bioestad/wp-content/uploads/supuesto10.txt", header = TRUE)
factorial
## Altura Carbono Presion Velocidad
## 1 10 10 25 200
## 2 11 12 25 200
## 3 2 14 25 200
## 4 3 10 25 250
## 5 2 12 25 250
## 6 4 14 25 250
## 7 5 10 30 200
## 8 5 12 30 200
## 9 -3 14 30 200
## 10 -1 10 30 250
## 11 -3 12 30 250
## 12 1 14 30 250
factorial$Carbono <- factor(factorial$Carbono)
factorial$Velocidad <- factor(factorial$Velocidad)
factorial$Presion <- factor(factorial$Presion)
mod <- aov(Altura~ Carbono + Presion + Velocidad + Carbono*Presion + Carbono*Velocidad + Presion*Velocidad , data = factorial )
#Altura: Nombre de la columna de las observaciones
#Carbono: Nombre de la columna en la que está representado el primer factor
#Presion: Nombre de la columna en la que está representado el segundo factor
#Velocidad: Nombre de la columna en la que está representado el tercer factor
#Carbono*Presion, Carbono*Velocidad y Presion*Velocidad hace referencia a las distintas interacciones.
#data= data.frame en el que están guardados los datos
summary(mod)
## Df Sum Sq Mean Sq F value Pr(>F)
## Carbono 2 24.50 12.25 147 0.00676 **
## Presion 1 65.33 65.33 784 0.00127 **
## Velocidad 1 48.00 48.00 576 0.00173 **
## Carbono:Presion 2 1.17 0.58 7 0.12500
## Carbono:Velocidad 2 75.50 37.75 453 0.00220 **
## Presion:Velocidad 1 1.33 1.33 16 0.05719 .
## Residuals 2 0.17 0.08
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod <- aov(Altura~ Carbono + Presion + Velocidad + Carbono*Velocidad, data = factorial )
summary(mod)
## Df Sum Sq Mean Sq F value Pr(>F)
## Carbono 2 24.50 12.25 22.97 0.003019 **
## Presion 1 65.33 65.33 122.50 0.000105 ***
## Velocidad 1 48.00 48.00 90.00 0.000220 ***
## Carbono:Velocidad 2 75.50 37.75 70.78 0.000215 ***
## Residuals 5 2.67 0.53
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Diseño factorial de tres factores con replicación
factorial <- read.table("https://wpd.ugr.es/~bioestad/wp-content/uploads/supuesto11.txt", header = TRUE)
factorial
## Altura Carbono Presion Velocidad
## 1 10 10 25 200
## 2 20 10 25 200
## 3 11 12 25 200
## 4 9 12 25 200
## 5 2 14 25 200
## 6 -1 14 25 200
## 7 3 10 25 250
## 8 5 10 25 250
## 9 2 12 25 250
## 10 5 12 25 250
## 11 4 14 25 250
## 12 7 14 25 250
## 13 5 10 30 200
## 14 9 10 30 200
## 15 5 12 30 200
## 16 4 12 30 200
## 17 -3 14 30 200
## 18 -2 14 30 200
## 19 -1 10 30 250
## 20 -3 10 30 250
## 21 -3 12 30 250
## 22 2 12 30 250
## 23 1 14 30 250
## 24 3 14 30 250
factorial$Carbono <- factor(factorial$Carbono)
factorial$Velocidad <- factor(factorial$Velocidad)
factorial$Presion <- factor(factorial$Presion)
mod <- aov(Altura~ Carbono + Presion + Velocidad + Carbono*Presion + Carbono*Velocidad + Presion*Velocidad + Carbono*Velocidad*Presion, data = factorial )
#Altura: Nombre de la columna de las observaciones
#Carbono: Nombre de la columna en la que está representado el primer factor
#Presion: Nombre de la columna en la que está representado el segundo factor
#Velocidad: Nombre de la columna en la que está representado el tercer factor
#Carbono*Presion, Carbono*Velocidad, Presion*Velocidad y Carbono*Velocidad*Presion hace referencia a las distintas interacciones.
summary(mod)
## Df Sum Sq Mean Sq F value Pr(>F)
## Carbono 2 88.08 44.04 5.683 0.018350 *
## Presion 1 150.00 150.00 19.355 0.000866 ***
## Velocidad 1 80.67 80.67 10.409 0.007270 **
## Carbono:Presion 2 14.25 7.12 0.919 0.425122
## Carbono:Velocidad 2 230.58 115.29 14.876 0.000564 ***
## Presion:Velocidad 1 1.50 1.50 0.194 0.667799
## Carbono:Presion:Velocidad 2 1.75 0.88 0.113 0.894175
## Residuals 12 93.00 7.75
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod <- aov(Altura~ Carbono + Presion + Velocidad + Carbono*Presion + Carbono*Velocidad + Presion*Velocidad, data = factorial )
summary(mod)
## Df Sum Sq Mean Sq F value Pr(>F)
## Carbono 2 88.08 44.04 6.507 0.010038 *
## Presion 1 150.00 150.00 22.164 0.000336 ***
## Velocidad 1 80.67 80.67 11.919 0.003886 **
## Carbono:Presion 2 14.25 7.12 1.053 0.375033
## Carbono:Velocidad 2 230.58 115.29 17.035 0.000178 ***
## Presion:Velocidad 1 1.50 1.50 0.222 0.645047
## Residuals 14 94.75 6.77
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod <- aov(Altura~ Carbono + Presion + Velocidad + Carbono*Velocidad, data = factorial )
summary(mod)
## Df Sum Sq Mean Sq F value Pr(>F)
## Carbono 2 88.08 44.04 6.776 0.006856 **
## Presion 1 150.00 150.00 23.077 0.000166 ***
## Velocidad 1 80.67 80.67 12.410 0.002612 **
## Carbono:Velocidad 2 230.58 115.29 17.737 6.91e-05 ***
## Residuals 17 110.50 6.50
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Ejercicios
#Se realiza un estudio del contenido de azufre en cinco yacimientos de carbón.
#Se toman muestras aleatoriamente de cada uno de los yacimientos y se analizan.
#Los datos del porcentaje de azufre por muestra se indican en la tabla adjunta.
#diseño unifactorial totalmente aleatorizado de efectos fijos no-equilibrado.
porcentaje <- read.table("https://wpd.ugr.es/~bioestad/wp-content/uploads/guiado1-1.txt", header = TRUE)
porcentaje
## Azufre Yacimiento
## 1 151 1
## 2 192 1
## 3 108 1
## 4 204 1
## 5 214 1
## 6 176 1
## 7 117 1
## 8 169 2
## 9 64 2
## 10 90 2
## 11 141 2
## 12 101 2
## 13 128 2
## 14 159 2
## 15 156 2
## 16 122 3
## 17 132 3
## 18 139 3
## 19 133 3
## 20 154 3
## 21 104 3
## 22 225 3
## 23 149 3
## 24 130 3
## 25 75 4
## 26 126 4
## 27 69 4
## 28 62 4
## 29 90 4
## 30 120 4
## 31 32 4
## 32 73 4
## 33 80 5
## 34 90 5
## 35 124 5
## 36 82 5
## 37 72 5
## 38 57 5
## 39 118 5
## 40 54 5
## 41 130 5
#Variable respuesta: Contenido de Azufre
#Factor: Tipo de yacimiento con cinco niveles. Es un factor de Efectos fijos ya que viene decidido qué niveles concretos se van a utilizar
porcentaje$Yacimiento<-factor(porcentaje$Yacimiento)
porcentaje$Yacimiento
## [1] 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 5 5 5 5 5 5
## [39] 5 5 5
## Levels: 1 2 3 4 5
mod <- aov(Azufre ~ Yacimiento, data = porcentaje)
summary(mod)
## Df Sum Sq Mean Sq F value Pr(>F)
## Yacimiento 4 40433 10108 8.534 5.97e-05 ***
## Residuals 36 42640 1184
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#2. Si se rechaza la hipótesis nula que las medias de porcentaje de azufre en los cinco yacimientos es la misma, determinar que medias difieren entre sí utilizando el método de comparaciones múltiples de Tukey.
mod.tukey <- TukeyHSD(mod, ordered = TRUE)
mod.tukey
## Tukey multiple comparisons of means
## 95% family-wise confidence level
## factor levels have been ordered
##
## Fit: aov(formula = Azufre ~ Yacimiento, data = porcentaje)
##
## $Yacimiento
## diff lwr upr p adj
## 5-4 8.791667 -39.217257 56.80059 0.9841364
## 2-4 45.125000 -4.275775 94.52577 0.0873389
## 3-4 62.236111 14.227188 110.24503 0.0057086
## 1-4 85.125000 33.990340 136.25966 0.0002709
## 2-5 36.333333 -11.675590 84.34226 0.2131394
## 3-5 53.444444 6.868947 100.01994 0.0177365
## 1-5 76.333333 26.542032 126.12463 0.0008288
## 3-2 17.111111 -30.897812 65.12003 0.8429902
## 1-2 40.000000 -11.134660 91.13466 0.1865081
## 1-3 22.888889 -26.902412 72.68019 0.6809794
#3. Estudiar las hipótesis de modelo: Homocedasticidad (Homogeneidad de las varianzas por grupo), Independencia y Normalidad.
bartlett.test(porcentaje$Azufre, porcentaje$Yacimiento)
##
## Bartlett test of homogeneity of variances
##
## data: porcentaje$Azufre and porcentaje$Yacimiento
## Bartlett's K-squared = 1.2655, df = 4, p-value = 0.8672
#independencia
layout(matrix(c(1,2,3,4),2,2))
plot(mod)
shapiro.test(mod$residuals)
##
## Shapiro-Wilk normality test
##
## data: mod$residuals
## W = 0.97902, p-value = 0.6384
# Ejercicio 2
#Se realiza un estudio sobre el efecto del fotoperiodo y del genotipo en el periodo latente de infección del moho de cebada aislado AB3.
#Se obtienen cincuenta hojas de cuatro genotipos distintos.
#Cada grupo es infectado y posteriormente expuesto a diferente fotoperiodo.
#Los distintos fotoperiodos se trataron como bloques y se obtuvieron los siguientes datos de los totales para los bloques y tratamientos. La respuesta anotada es el número de días hasta la aparición de síntomas visibles.
#diseño en bloques completos aleatorizados.
#Factor: Genotipo que tiene cuatro niveles. Es un factor de efectos fijos ya que viene decidido qué niveles concretos se van a utilizar.
#Bloque: Fotoperiodo que tiene cinco niveles. Es un factor de efectos fijos ya que viene decidido qué niveles concretos se van a utilizar.
dias <- read.table("https://wpd.ugr.es/~bioestad/wp-content/uploads/guiado2txt.txt", header = TRUE)
dias
## Días Fotoperiodo Genotipo
## 1 630 1 1
## 2 640 1 2
## 3 640 1 3
## 4 660 1 4
## 5 610 2 1
## 6 630 2 2
## 7 630 2 3
## 8 660 2 4
## 9 560 3 1
## 10 600 3 2
## 11 650 3 3
## 12 620 3 4
## 13 570 4 1
## 14 620 4 2
## 15 620 4 3
## 16 610 4 4
## 17 590 5 1
## 18 620 5 2
## 19 580 5 3
## 20 630 5 4
dias$Fotoperiodo = factor(dias$Fotoperiodo)
dias$Fotoperiodo
## [1] 1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4 5 5 5 5
## Levels: 1 2 3 4 5
dias$Genotipo = factor(dias$Genotipo)
dias$Genotipo
## [1] 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4
## Levels: 1 2 3 4
mod = aov(Días ~ Genotipo + Fotoperiodo, data = dias)
summary(mod)
## Df Sum Sq Mean Sq F value Pr(>F)
## Genotipo 3 5255 1751.7 5.041 0.0173 *
## Fotoperiodo 4 5030 1257.5 3.619 0.0371 *
## Residuals 12 4170 347.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#2. En caso de que influyan significativamente alguno de los dos factores, extraer conclusiones utilizando el método de Duncan.
library(agricolae)
duncan <-duncan.test(mod, "Genotipo" )
duncan
## $statistics
## MSerror Df Mean CV
## 347.5 12 618.5 3.013962
##
## $parameters
## test name.t ntr alpha
## Duncan Genotipo 4 0.05
##
## $duncan
## Table CriticalRange
## 2 3.081307 25.68782
## 3 3.225244 26.88778
## 4 3.312453 27.61481
##
## $means
## Días std r Min Max Q25 Q50 Q75
## 1 592 28.63564 5 560 630 570 590 610
## 2 622 14.83240 5 600 640 620 620 630
## 3 624 27.01851 5 580 650 620 630 640
## 4 636 23.02173 5 610 660 620 630 660
##
## $comparison
## NULL
##
## $groups
## Días groups
## 4 636 a
## 3 624 a
## 2 622 a
## 1 592 b
##
## attr(,"class")
## [1] "group"
duncan <-duncan.test(mod, "Fotoperiodo" )
duncan
## $statistics
## MSerror Df Mean CV
## 347.5 12 618.5 3.013962
##
## $parameters
## test name.t ntr alpha
## Duncan Fotoperiodo 5 0.05
##
## $duncan
## Table CriticalRange
## 2 3.081307 28.71986
## 3 3.225244 30.06145
## 4 3.312453 30.87430
## 5 3.370172 31.41228
##
## $means
## Días std r Min Max Q25 Q50 Q75
## 1 642.5 12.58306 4 630 660 637.5 640 645.0
## 2 632.5 20.61553 4 610 660 625.0 630 637.5
## 3 607.5 37.74917 4 560 650 590.0 610 627.5
## 4 605.0 23.80476 4 570 620 600.0 615 620.0
## 5 605.0 23.80476 4 580 630 587.5 605 622.5
##
## $comparison
## NULL
##
## $groups
## Días groups
## 1 642.5 a
## 2 632.5 ab
## 3 607.5 b
## 4 605.0 b
## 5 605.0 b
##
## attr(,"class")
## [1] "group"
bartlett.test(dias$Días, dias$Fotoperiodo)
##
## Bartlett test of homogeneity of variances
##
## data: dias$Días and dias$Fotoperiodo
## Bartlett's K-squared = 3.0629, df = 4, p-value = 0.5474
bartlett.test(dias$Días, dias$Genotipo)
##
## Bartlett test of homogeneity of variances
##
## data: dias$Días and dias$Genotipo
## Bartlett's K-squared = 1.6252, df = 3, p-value = 0.6537
plot(mod)
shapiro.test(mod$residuals)
##
## Shapiro-Wilk normality test
##
## data: mod$residuals
## W = 0.95316, p-value = 0.4176
#Ejercicio 3
library(readxl)
plantasejemplo3<-read_excel("Data/plantasejemplo3.xlsx")
#¿Se puede afirmar que los distintos niveles de agua influyen en la longitud del tallo de los guisantes? ¿Y el tipo de planta?
#Variable respuesta: Longitud del tallo
#Factor A: Nivel del agua con tres niveles. Es un factor de Efectos fijos.
#Factor B: Tipo de planta con dos niveles. Es un factor de Efectos fijos.
attach(plantasejemplo3)
mod<-aov(long_tallo ~ nivel_agua*Tipo_de_planta)
summary(mod) #da igual si cambio las variables categoricas por letras o numeros, tampoco importa el orden en el que coloque las columnas
## Df Sum Sq Mean Sq F value Pr(>F)
## nivel_agua 2 10774 5387 572.56 < 2e-16 ***
## Tipo_de_planta 1 1246 1246 132.44 1.58e-12 ***
## nivel_agua:Tipo_de_planta 2 514 257 27.32 1.75e-07 ***
## Residuals 30 282 9
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# PCA
#Construccion
library(vegan)
## Loading required package: permute
## Loading required package: lattice
## Registered S3 methods overwritten by 'vegan':
## method from
## plot.rda klaR
## predict.rda klaR
## print.rda klaR
## This is vegan 2.5-7
##
## Attaching package: 'vegan'
## The following object is masked from 'package:outliers':
##
## scores
library(ade4)#Algunas opciones requieren este paquete!
library(PerformanceAnalytics)
## Loading required package: xts
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
##
## Attaching package: 'PerformanceAnalytics'
## The following objects are masked from 'package:agricolae':
##
## kurtosis, skewness
## The following object is masked from 'package:graphics':
##
## legend
#Para el ejemplo vamos a utilizar el dataset “doubs”, disponible en Vegan:
data("doubs")
#Dado que solo vamos a trabajar con los datos ambientales, debemos seleccionar el elemento “env”, en el cual están los datos ambientales:
env<-doubs$env
head(env)
## dfs alt slo flo pH har pho nit amm oxy bdo
## 1 3 934 6.176 84 79 45 1 20 0 122 27
## 2 22 932 3.434 100 80 40 2 20 10 103 19
## 3 102 914 3.638 180 83 52 5 22 5 105 35
## 4 185 854 3.497 253 80 72 10 21 0 110 13
## 5 215 849 3.178 264 81 84 38 52 20 80 62
## 6 324 846 3.497 286 79 60 20 15 0 102 53
#En ocasiones podemos necesitar eliminar elementos desde R, sin afectar el archivo de origen del cual cargamos nuestros datos. En este caso vamos a eliminar la fila #8 de nuestro set de datos:
env<-env[-8, ] #elimina fila (fila, columna)
#Note que R sobreescribe sobre el mismo vector!!
chart.Correlation(env)
cor(env)
## dfs alt slo flo pH har
## dfs 1.00000000 -0.93893906 -0.7555877 0.94734109 0.01522514 0.73030747
## alt -0.93893906 1.00000000 0.7658280 -0.86280579 -0.05027785 -0.78549505
## slo -0.75558775 0.76582805 1.0000000 -0.71656547 -0.27713925 -0.66664648
## flo 0.94734109 -0.86280579 -0.7165655 1.00000000 0.03312637 0.73662526
## pH 0.01522514 -0.05027785 -0.2771392 0.03312637 1.00000000 0.08451511
## har 0.73030747 -0.78549505 -0.6666465 0.73662526 0.08451511 1.00000000
## pho 0.47322419 -0.43698762 -0.3997436 0.37868776 -0.07940250 0.37318607
## nit 0.73874252 -0.75263077 -0.6070823 0.59315083 -0.04046292 0.53495392
## amm 0.40795909 -0.38107709 -0.3492304 0.29252515 -0.12200180 0.29613600
## oxy -0.57051149 0.42467525 0.4937437 -0.42109447 0.19239804 -0.37363740
## bdo 0.43542349 -0.38238374 -0.3346889 0.29518914 -0.16170564 0.33697473
## pho nit amm oxy bdo
## dfs 0.4732242 0.73874252 0.4079591 -0.5705115 0.4354235
## alt -0.4369876 -0.75263077 -0.3810771 0.4246753 -0.3823837
## slo -0.3997436 -0.60708231 -0.3492304 0.4937437 -0.3346889
## flo 0.3786878 0.59315083 0.2925252 -0.4210945 0.2951891
## pH -0.0794025 -0.04046292 -0.1220018 0.1923980 -0.1617056
## har 0.3731861 0.53495392 0.2961360 -0.3736374 0.3369747
## pho 1.0000000 0.80093149 0.9699345 -0.7575015 0.9091698
## nit 0.8009315 1.00000000 0.8022323 -0.6867146 0.6832300
## amm 0.9699345 0.80223230 1.0000000 -0.7462700 0.9028247
## oxy -0.7575015 -0.68671462 -0.7462700 1.0000000 -0.8398165
## bdo 0.9091698 0.68322998 0.9028247 -0.8398165 1.0000000
#A continuación podemos elaborar el PCA a partir del dataset:
acp<-rda(env,scale=T) #Scale realiza estandarización variables; utiliza la matriz de correlaciones.
#La opción summary despliega la información obtenida en el PCA:
summary(acp)
##
## Call:
## rda(X = env, scale = T)
##
## Partitioning of correlations:
## Inertia Proportion
## Total 11 1
## Unconstrained 11 1
##
## Eigenvalues, and their contribution to the correlations
##
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Eigenvalue 6.446 2.2252 0.99825 0.39831 0.36224 0.25361 0.16106
## Proportion Explained 0.586 0.2023 0.09075 0.03621 0.03293 0.02306 0.01464
## Cumulative Proportion 0.586 0.7883 0.87903 0.91524 0.94817 0.97123 0.98587
## PC8 PC9 PC10 PC11
## Eigenvalue 0.11062 0.022949 0.017476 0.004378
## Proportion Explained 0.01006 0.002086 0.001589 0.000398
## Cumulative Proportion 0.99593 0.998013 0.999602 1.000000
##
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores: 4.189264
##
##
## Species scores
##
## PC1 PC2 PC3 PC4 PC5 PC6
## dfs 1.10786 0.4901 -0.20230 0.045871 -0.187500 -0.20053
## alt -1.06604 -0.5607 0.15135 0.193929 0.096582 -0.05541
## slo -0.95746 -0.5295 -0.25933 -0.352365 0.066943 -0.39476
## flo 0.98732 0.6262 -0.23731 0.009137 -0.105599 -0.31209
## pH -0.02679 0.4991 1.14319 -0.014864 -0.005398 -0.17594
## har 0.91493 0.5514 -0.09538 -0.128923 0.639664 0.07145
## pho 1.02238 -0.6481 0.20016 -0.178923 0.040772 -0.04615
## nit 1.14057 -0.1628 0.05081 -0.283125 -0.272940 0.20434
## amm 0.96776 -0.7382 0.18740 -0.198493 -0.021140 0.04733
## oxy -0.99282 0.4993 0.04744 -0.543295 -0.012373 0.06953
## bdo 0.96051 -0.7274 0.09630 0.089149 0.178489 -0.14521
##
##
## Site scores (weighted sums of species scores)
##
## PC1 PC2 PC3 PC4 PC5 PC6
## 1 -1.35309 -1.0614 -0.626184 -1.14841 -1.0494 -1.821e+00
## 2 -1.05499 -0.7841 0.195705 0.90757 -1.7294 2.655e-01
## 3 -0.97457 -0.4890 1.340010 0.61152 -0.8592 -7.288e-01
## 4 -0.90281 -0.3118 0.000372 0.17712 0.2052 5.288e-01
## 5 -0.45666 -0.6973 0.550276 1.16964 1.2500 1.273e-01
## 6 -0.81070 -0.7590 -0.318284 0.77249 -0.2537 1.263e-01
## 7 -0.85277 -0.1858 0.231151 -0.37159 1.3564 -3.616e-01
## 9 -0.27926 -0.4610 0.061754 1.60909 1.2127 8.870e-01
## 10 -0.59145 -0.5554 -1.595293 -0.35281 0.6397 -2.559e-01
## 11 -0.34078 0.3167 0.005014 -1.23777 0.7883 -2.345e-01
## 12 -0.44165 0.3209 -0.697214 -0.46903 0.3928 5.963e-01
## 13 -0.39855 0.6314 0.003511 -0.90415 1.0778 5.002e-05
## 14 -0.22649 0.7350 0.946755 -0.87823 0.8730 3.814e-02
## 15 -0.21927 1.0432 2.269502 -0.19576 -0.1024 -2.360e-01
## 16 -0.16778 0.2507 -0.340084 -0.54136 -0.1054 6.195e-01
## 17 0.14914 0.3628 -0.171537 -0.14337 -0.1372 1.435e+00
## 18 0.08633 0.3674 -0.238531 -0.44014 -0.2901 1.015e+00
## 19 0.10967 0.4821 0.232435 -0.28363 -0.7316 9.414e-01
## 20 0.18575 0.3732 -0.268777 -0.69333 -0.9729 1.131e+00
## 21 0.16766 0.3112 -0.834583 0.21603 -0.7147 5.100e-01
## 22 0.13041 0.4842 -0.108135 0.18812 -0.2895 -5.716e-01
## 23 1.28519 -1.3164 0.715652 -0.57204 0.6499 -1.021e+00
## 24 1.01542 -0.4735 0.019519 1.43289 0.5981 1.440e-01
## 25 2.10059 -2.1406 0.361157 -1.21146 -0.1786 4.424e-01
## 26 0.89379 -0.1213 -0.671652 0.86581 -0.3046 -8.871e-02
## 27 0.61092 0.3178 -0.139402 0.31511 -0.8178 -9.684e-01
## 28 0.82353 0.8569 0.802337 -0.04239 -0.8747 5.103e-02
## 29 0.67793 1.0652 -1.729365 0.28387 0.2944 -1.028e+00
## 30 0.83450 1.4380 0.003890 0.93621 0.0730 -1.543e+00
summary(acp,scaling=1)
##
## Call:
## rda(X = env, scale = T)
##
## Partitioning of correlations:
## Inertia Proportion
## Total 11 1
## Unconstrained 11 1
##
## Eigenvalues, and their contribution to the correlations
##
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Eigenvalue 6.446 2.2252 0.99825 0.39831 0.36224 0.25361 0.16106
## Proportion Explained 0.586 0.2023 0.09075 0.03621 0.03293 0.02306 0.01464
## Cumulative Proportion 0.586 0.7883 0.87903 0.91524 0.94817 0.97123 0.98587
## PC8 PC9 PC10 PC11
## Eigenvalue 0.11062 0.022949 0.017476 0.004378
## Proportion Explained 0.01006 0.002086 0.001589 0.000398
## Cumulative Proportion 0.99593 0.998013 0.999602 1.000000
##
## Scaling 1 for species and site scores
## * Sites are scaled proportional to eigenvalues
## * Species are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores: 4.189264
##
##
## Species scores
##
## PC1 PC2 PC3 PC4 PC5 PC6
## dfs 1.447 1.0896 -0.6715 0.24106 -1.03324 -1.3207
## alt -1.393 -1.2467 0.5024 1.01912 0.53222 -0.3649
## slo -1.251 -1.1772 -0.8608 -1.85173 0.36890 -2.5998
## flo 1.290 1.3924 -0.7878 0.04801 -0.58192 -2.0554
## pH -0.035 1.1098 3.7948 -0.07811 -0.02975 -1.1587
## har 1.195 1.2260 -0.3166 -0.67751 3.52493 0.4706
## pho 1.336 -1.4410 0.6644 -0.94026 0.22468 -0.3040
## nit 1.490 -0.3619 0.1687 -1.48786 -1.50406 1.3457
## amm 1.264 -1.6413 0.6221 -1.04311 -0.11649 0.3117
## oxy -1.297 1.1100 0.1575 -2.85509 -0.06818 0.4579
## bdo 1.255 -1.6174 0.3197 0.46849 0.98358 -0.9563
##
##
## Site scores (weighted sums of species scores)
##
## PC1 PC2 PC3 PC4 PC5 PC6
## 1 -1.03580 -0.47737 -0.1886364 -0.218531 -0.19043 -2.765e-01
## 2 -0.80759 -0.35267 0.0589557 0.172702 -0.31383 4.031e-02
## 3 -0.74604 -0.21992 0.4036746 0.116366 -0.15591 -1.107e-01
## 4 -0.69110 -0.14022 0.0001121 0.033704 0.03725 8.029e-02
## 5 -0.34957 -0.31363 0.1657693 0.222572 0.22683 1.933e-02
## 6 -0.62059 -0.34137 -0.0958822 0.146997 -0.04605 1.918e-02
## 7 -0.65279 -0.08358 0.0696338 -0.070711 0.24615 -5.490e-02
## 9 -0.21377 -0.20732 0.0186032 0.306194 0.22006 1.347e-01
## 10 -0.45276 -0.24979 -0.4805779 -0.067136 0.11609 -3.886e-02
## 11 -0.26087 0.14244 0.0015104 -0.235536 0.14305 -3.560e-02
## 12 -0.33808 0.14434 -0.2100341 -0.089252 0.07127 9.055e-02
## 13 -0.30509 0.28397 0.0010576 -0.172050 0.19558 7.596e-06
## 14 -0.17338 0.33057 0.2852075 -0.167119 0.15842 5.792e-03
## 15 -0.16785 0.46919 0.6836817 -0.037251 -0.01858 -3.583e-02
## 16 -0.12843 0.11274 -0.1024494 -0.103016 -0.01913 9.406e-02
## 17 0.11417 0.16318 -0.0516750 -0.027282 -0.02491 2.178e-01
## 18 0.06609 0.16523 -0.0718567 -0.083754 -0.05265 1.542e-01
## 19 0.08395 0.21681 0.0700204 -0.053972 -0.13277 1.429e-01
## 20 0.14219 0.16784 -0.0809685 -0.131933 -0.17656 1.717e-01
## 21 0.12834 0.13995 -0.2514160 0.041108 -0.12969 7.744e-02
## 22 0.09983 0.21776 -0.0325754 0.035797 -0.05254 -8.679e-02
## 23 0.98382 -0.59208 0.2155883 -0.108853 0.11794 -1.550e-01
## 24 0.77731 -0.21296 0.0058801 0.272664 0.10854 2.186e-02
## 25 1.60801 -0.96276 0.1087977 -0.230529 -0.03241 6.717e-02
## 26 0.68419 -0.05457 -0.2023335 0.164755 -0.05528 -1.347e-02
## 27 0.46766 0.14295 -0.0419945 0.059962 -0.14840 -1.470e-01
## 28 0.63041 0.38542 0.2417020 -0.008066 -0.15873 7.749e-03
## 29 0.51895 0.47910 -0.5209667 0.054018 0.05342 -1.561e-01
## 30 0.63881 0.64676 0.0011718 0.178152 0.01325 -2.343e-01
#Scaling 1: Gráfico biplot de distancias. “Las distancias entre los objetos
#en el gráfico biplot son las aproximaciones de las distancias Euclideanas
#in un espacio multidimensional”,
#es decir, los ángulos entre los vectores NO reflejan sus correlaciones.
#Scaling 2: Gráfico de correlaciones.
#“Las distancias entre los objetos en el gráfico biplot NO son las aproximaciones de las
#distancias Euclideanas in un espacio multidimensional”, es decir, los ángulos entre los descriptores reflejan las correlaciones.
# Graficas
par(mfrow = c(1, 2)) #divide la pantalla para obtener los dos gráficos
biplot(acp, scaling = 1, main = "PCA - scaling: 1")
biplot(acp, main = "PCA - scaling: 2")
prin_comp <- prcomp(env, scale. = T)
summary(prin_comp)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 2.539 1.4917 0.99912 0.63112 0.60186 0.50360 0.40132
## Proportion of Variance 0.586 0.2023 0.09075 0.03621 0.03293 0.02306 0.01464
## Cumulative Proportion 0.586 0.7883 0.87903 0.91524 0.94817 0.97123 0.98587
## PC8 PC9 PC10 PC11
## Standard deviation 0.33260 0.15149 0.13220 0.06617
## Proportion of Variance 0.01006 0.00209 0.00159 0.00040
## Cumulative Proportion 0.99593 0.99801 0.99960 1.00000