Ejercicio con datos de reflectancia espectral
Yeimy Maryury Montilla Montilla; Carlos Eduardo Caceres Gonzalez
02/05/2021
Cargue de los datos
library(readxl)
Reflectance <- read_excel("F:/UNAL_MAESTRIA_UNAL/METODOS_MULTIVARIADOS/Trabajo2/datos.xlsx",
sheet = "Spectral_Reflectance")
data=data.frame(Reflectance);data## X560_nm X720_nm Species
## 1 9.33 19.14 SS
## 2 8.74 19.55 SS
## 3 9.31 19.24 SS
## 4 8.27 16.37 SS
## 5 10.22 25.00 SS
## 6 10.13 25.32 SS
## 7 10.42 27.12 SS
## 8 10.62 26.28 SS
## 9 15.25 38.89 SS
## 10 16.22 36.67 SS
## 11 17.24 40.74 SS
## 12 12.77 67.50 SS
## 13 12.07 33.03 JL
## 14 11.03 32.37 JL
## 15 12.48 31.31 JL
## 16 12.12 33.33 JL
## 17 15.38 40.00 JL
## 18 14.21 40.48 JL
## 19 9.69 33.90 JL
## 20 14.35 40.15 JL
## 21 38.71 77.14 JL
## 22 44.74 78.57 JL
## 23 36.67 71.43 JL
## 24 37.21 45.00 JL
## 25 8.73 23.27 LP
## 26 7.94 20.87 LP
## 27 8.37 22.16 LP
## 28 7.86 21.78 LP
## 29 8.45 26.32 LP
## 30 6.79 22.73 LP
## 31 8.34 26.67 LP
## 32 7.54 24.87 LP
## 33 14.04 44.44 LP
## 34 13.51 37.93 LP
## 35 13.33 37.93 LP
## 36 12.77 60.87 LPPunto 1: MANOVA
library(ggplot2)
ggplot(data=data, aes(X560_nm,X720_nm, color=Species)) + geom_point() +geom_smooth(method = "lm", se = FALSE) + ggtitle("Reflectancia espectral 560 nm vs 720 nm")## `geom_smooth()` using formula 'y ~ x'
Mediante el grafico se puede suponer que la reflectancia espectral en 560 nm y 720nm es diferente para cada una de las especies. pero quizás es afectada por posibles valores atípicos.
Se considera como hipótesis que las tres especies son iguales.
\[\begin{bmatrix} \mu_{560 nm}\\ \mu_{720 nm} \end{bmatrix}_{Specie = JL} = \begin{bmatrix} \mu_{560 nm}\\ \mu_{720 nm} \end{bmatrix}_{Specie = LP} = \begin{bmatrix} \mu_{560 nm}\\ \mu_{720 nm} \end{bmatrix}_{Specie = SS} \]
Manova = manova(cbind(data$X560_nm,data$X720_nm) ~ data$Species)
summary(Manova, test = 'Wilks')## Df Wilks approx F num Df den Df Pr(>F)
## data$Species 2 0.67704 3.4452 4 64 0.013 *
## Residuals 33
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1El valor de Pr(>F) está indicando que las especies JL, LP y SS son diferentes considerando los dos datos de reflectancia espectral.
Punto 2: ANOVA por respuesta espectral
Respuesta espectral de 560 nm
\[H_0: \mu_{LS}=\mu_{LP}=\mu_{SS}\\ H_a: H_0~ \text{es falsa}\]
Anova_560=aov(data$X560_nm ~ data$Species)
summary(Anova_560)## Df Sum Sq Mean Sq F value Pr(>F)
## data$Species 2 965.2 482.6 7.415 0.00219 **
## Residuals 33 2147.7 65.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Es aproximadamente site (7) veces más grade la varianza entre las especies que entre las medidas de reflectancia de cada especie. Por tanto, se puede decir que H0 es falsa y las especies no son iguales considerando solo la respuesta espectral de 720 nm.
boxplot(data$X560_nm ~ data$Species)
Al hacer una gráfica de boxplot se evidencia que la variabilidad de las tres especies es muy diferente y se confirma con el cálculo de varianza. Especialmente en la especie JL.
tapply(data$X560_nm, data$Species, var)## JL LP SS
## 178.387827 7.414354 9.444588Respuesta espectral de 720 nm
\[H_0: \mu_{LS}=\mu_{LP}=\mu_{SS}\\ H_a: H_0~ \text{es falsa}\]
Anova_720=aov(data$X720_nm ~ data$Species)
summary(Anova_720)## Df Sum Sq Mean Sq F value Pr(>F)
## data$Species 2 2027 1013.4 4.437 0.0196 *
## Residuals 33 7537 228.4
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Es aproximadamente cuatro (4) veces más grade la varianza entre las especies que entre las medidas de reflectancia de cada especie. Por tanto, se puede decir que H0 es falsa y las especies no son iguales considerando solo la respuesta espectral de 720 nm.
boxplot(data$X720_nm ~ data$Species)
Al hacer una gráfica de boxplot se evidencia que la variabilidad de la especies JL y LP quizás es Similar y el cálculo de varianza no permite suponer una diferencia marcada.
tapply(data$X720_nm, data$Species, var)## JL LP SS
## 332.1606 149.1506 203.8704Punto 3: Test de correlación
ggplot(data,aes(X560_nm,X720_nm,fill = Species))+
facet_grid(Species~., scales = 'free')+
geom_point(size=2, shape=23)+
geom_smooth(method="lm")+
theme_bw() +scale_y_continuous(limits = c(0, 100))## `geom_smooth()` using formula 'y ~ x'
Conforme a la gráfica se puede suponer que existe correlación entre las refelectacias por cada especie. ### Correlación General
cor(data$X560_nm,data$X720_nm, method = "pearson",)## [1] 0.8130816cor(data$X560_nm,data$X720_nm, method = "spearman")## [1] 0.8809371El resultado indica que existe correlación entre las dos variables. Se emplea el método “Pearson” suponiendo que es una relación lineal. Igualmente se aplicó el método “spearman” donde también se obtuvo un valor significativo.
Se realizó el ejercicio para cada una de las especies y se obtienen valores muy similares. Por tanto, se puede decir que entre los datos de reflectancia de 720_nm y 560_nm están correlacionado. En este sentido es más confiable el Manova.
Correlación en la especie JL
data_JL=subset(data,data$Species=="JL")
cor(data_JL$X560_nm,data_JL$X720_nm, method = "pearson",)## [1] 0.911518cor(data_JL$X560_nm,data_JL$X720_nm, method = "spearman")## [1] 0.8531469Correlación en la especie LP
data_LP=subset(data,data$Species=="LP")
cor(data_LP$X560_nm,data_JL$X720_nm, method = "pearson",)## [1] 0.8794062cor(data_LP$X560_nm,data_JL$X720_nm, method = "spearman")## [1] 0.5384615Correlación en la especie SS
data_SS=subset(data,data$Species=="SS")
cor(data_SS$X560_nm,data_JL$X720_nm, method = "pearson",)## [1] 0.9532662cor(data_SS$X560_nm,data_JL$X720_nm, method = "spearman")## [1] 0.8811189Punto 4: Normalidad por respuesta espectral
\[H_0: \text{La distribución es normal}\\ H_a: \text{La distribución NO es normal}\]
Respuesta espectral de 560 nm
X560_escal=scale(data$X560_nm)
X720_escal=scale(data$X720_nm)
data_scal=data.frame(X560_escal,X720_escal)ggplot(data,aes(x=X560_nm))+geom_histogram(aes(x=X560_nm, y= ..density..),alpha=0.4)+ geom_density()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(data_scal,aes(x=X560_escal))+geom_histogram(aes(x=X560_escal, y= ..density..),alpha=0.4)+ geom_density()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(data, aes(x=X560_nm,fill = Species))+facet_grid(Species~., scales = 'free')+geom_histogram(aes(x=X560_nm, y= ..density..))+ geom_density(alpha=0.4)+ scale_x_continuous(limits = c(-10, 100))+scale_y_continuous(limits = c(0, 0.2))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 6 rows containing missing values (geom_bar).
shapiro.test(data$X560_nm)##
## Shapiro-Wilk normality test
##
## data: data$X560_nm
## W = 0.64673, p-value = 4.538e-08Las gráficas anteriores permiten suponer que los datos de reflectancia de 560 nm no siguen una distribución normal, incluso se puede ver una distribución bimodal. Esto se confirma con la prueba de normalidad donde el “Shapiro-Wilk” donde se obtiene un P-Value muy pequeño. ### Respuesta espectral de 720_nm nm
ggplot(data, aes(x=X720_nm))+geom_histogram(aes(x=X720_nm, y= ..density..),alpha=0.4)+ geom_density()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(data_scal,aes(x=X720_escal))+geom_histogram(aes(x=X720_escal, y= ..density..),alpha=0.4)+ geom_density()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(data, aes(x=X720_nm,fill = Species))+facet_grid(Species~., scales = 'free')+geom_histogram(aes(x=X720_nm, y= ..density..))+ geom_density(alpha=0.4)+
scale_x_continuous(limits = c(-10, 100))+scale_y_continuous(limits = c(0, 0.2))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 6 rows containing missing values (geom_bar).
shapiro.test(data$X720_nm)##
## Shapiro-Wilk normality test
##
## data: data$X720_nm
## W = 0.84251, p-value = 0.0001286Las gráficas anteriores permiten suponer que los datos de reflectancia de 720 nm no siguen una distribución normal, incluso se puede ver una distribución bimodal. Esto se confirma con la prueba de normalidad donde el “Shapiro-Wilk” donde se obtiene un P-Value muy pequeño.
Punto 5: Normalidad Bivariada
library(royston)## The Royston's H test has been moved to 'MVN' package. 'royston'
## package will no longer be supported. Please use 'MVN' package for
## further analysis.royston.test(data[,-3])
## $H
## [1] 39.74407
##
## $p.value
## [1] 1.576579e-09
##
## $distribution
## [1] "Data analyzed have a non-normal distribution."Conforme a la prueba de normalidad multivariante de Royston se establece que los datos de refelctancia no tienen una distribución normal.
Punto 6: Igualdad de varianzas univariante
Respuesta espectral de 560 nm
ggplot(data, aes(x=X560_nm,fill = Species))+ geom_density(alpha=0.4)+
scale_x_continuous(limits = c(-10, 100))+scale_y_continuous(limits = c(0, 0.2))
library(car)## Loading required package: carDatalibrary(carData)
leveneTest(y = data$X560_nm, group = data$Species, center = "median")## Warning in leveneTest.default(y = data$X560_nm, group = data$Species, center =
## "median"): data$Species coerced to factor.## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 2 4.6452 0.01669 *
## 33
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1A partir de la prueba de leveneTest, el cual es menos susceptible a la No normalidad de los datos, se puede decir que al menos una especie presenta varianzas diferentes considerando la reflectancia de 560 nm. ### Respuesta espectral de 720 nm
ggplot(data, aes(x=X720_nm,fill = Species))+ geom_density(alpha=0.4) +
scale_x_continuous(limits = c(-10, 100))+scale_y_continuous(limits = c(0, 0.2))
bartlett.test(data$X720_nm ~ data$Species)##
## Bartlett test of homogeneity of variances
##
## data: data$X720_nm by data$Species
## Bartlett's K-squared = 1.7462, df = 2, p-value = 0.4177leveneTest(data$X720_nm,data$Species, "median")## Warning in leveneTest.default(data$X720_nm, data$Species, "median"):
## data$Species coerced to factor.## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 2 0.3817 0.6857
## 33A partir de la prueba de leveneTest, la cual es menos susceptible a la No normalidad de los datos, se puede decir que no se evidencia una diferencia clara entre la varianza de las especies.
Punto 7: Igualdad de matrices de varianzas y covarianzas de los tres tratamientos
Y = cbind(data$X560_nm, data$X720_nm)
# Matrix de varianzas y covarianza
df_Species = split(data[,1:2], data$Species)
S = lapply(df_Species, cov); S## $JL
## X560_nm X720_nm
## X560_nm 178.3878 221.8819
## X720_nm 221.8819 332.1606
##
## $LP
## X560_nm X720_nm
## X560_nm 7.414354 27.71932
## X720_nm 27.719318 149.15058
##
## $SS
## X560_nm X720_nm
## X560_nm 9.444588 28.47334
## X720_nm 28.473339 203.87040biotools::boxM(Y, data$Species)##
## Box's M-test for Homogeneity of Covariance Matrices
##
## data: Y
## Chi-Sq (approx.) = 41.969, df = 6, p-value = 1.865e-07De acuerdo con la prueba Box’s se puede determinar que las matrices de varianza generada por cada especie son diferentes.
Punto 8: outiliers univariantes
Respuesta espectral de 560 nm
library(outliers)
# Datos atipicos univariante
grubbs.test(data$X560_nm)##
## Grubbs test for one outlier
##
## data: data$X560_nm
## G = 3.22758, U = 0.69386, p-value = 0.008372
## alternative hypothesis: highest value 44.74 is an outlierRespuesta espectral de 720 nm
library(outliers)
# Datos atipicos univariante
grubbs.test(data$X720_nm)##
## Grubbs test for one outlier
##
## data: data$X720_nm
## G = 2.58808, U = 0.80316, p-value = 0.121
## alternative hypothesis: highest value 78.57 is an outlierlibrary(outliers)
# Datos atipicos univariante
grubbs.test(data$X560_nm)##
## Grubbs test for one outlier
##
## data: data$X560_nm
## G = 3.22758, U = 0.69386, p-value = 0.008372
## alternative hypothesis: highest value 44.74 is an outliergrubbs.test(data$X720_nm)##
## Grubbs test for one outlier
##
## data: data$X720_nm
## G = 2.58808, U = 0.80316, p-value = 0.121
## alternative hypothesis: highest value 78.57 is an outlierLa prueba grubbs indica que en la respuesta espectral de 560 nm los datos superiores a 44.74 correponden a valores atipicos, mientras que en la respuesta espectral de 720 nm son los datos superiores a 78.57.
Punto 9: outliers bivariantes con T2
Medias_R=apply(data[,1:2],2,mean)
T2=c()
for(j in 1:dim(data)[1])# desde 1 hasta el numero de filas
{
T2[j]=c((t(t(data[j,1:2])-(Medias_R)))%*%solve(var(data[,1:2]))%*%as.matrix(t(data[j,1:2])-(Medias_R)))
j
Medias_R
data[j,1:2]
}p=2
n=36
ftabla=qf(p = 0.05,df1 = 2,df2 = 34,lower.tail = F)
TT2=((p*(n-1))/(n-p))*ftabla;
Out=ifelse(T2>TT2,"Outliers","No Outliers")
Color=ifelse(T2>TT2,"Red","Green")data_Outliers=data.frame(data$X560_nm,data$X720_nm,data$Species,T2,TT2,Out,Color)
plot(T2,col=Color,pch=19,cex=0.85,xlab="Observación")
grid(20,20,col="lightblue")
abline(h=TT2)
ggplot(data=data_Outliers, aes(data.X560_nm,data.X720_nm, color=Out))+ geom_point()
Punto 10:comparaciones de medias por cada respuesta
library(ICSNP)## Loading required package: mvtnorm## Loading required package: ICSlibrary(doBy)
data_JL_LP=subset(data,data$Species!="SS")
summaryBy(X560_nm + X720_nm ~ Species, data=data_JL_LP, FUN=mean)## Species X560_nm.mean X720_nm.mean
## 1 JL 21.555000 46.3925
## 2 LP 9.805833 30.8200HotellingsT2(as.matrix(data_JL_LP[,1:2])~as.matrix(data_JL_LP[,3]), mu = c(0,0))##
## Hotelling's two sample T2-test
##
## data: as.matrix(data_JL_LP[, 1:2]) by as.matrix(data_JL_LP[, 3])
## T.2 = 4.2569, df1 = 2, df2 = 21, p-value = 0.02806
## alternative hypothesis: true location difference is not equal to c(0,0)Las medias entre las especies JL y LP No son iguales.
data_SS_LP=subset(data,data$Species!="JL")
summaryBy(X560_nm + X720_nm ~ Species, data=data_SS_LP, FUN=mean)## Species X560_nm.mean X720_nm.mean
## 1 LP 9.805833 30.82000
## 2 SS 11.543333 30.15167df_Species = split(data_SS_LP[,1:2], data_SS_LP$Species)
HotellingsT2(as.matrix(data_SS_LP[,1:2])~as.matrix(data_SS_LP[,3]), mu = c(0,0))##
## Hotelling's two sample T2-test
##
## data: as.matrix(data_SS_LP[, 1:2]) by as.matrix(data_SS_LP[, 3])
## T.2 = 2.4676, df1 = 2, df2 = 21, p-value = 0.109
## alternative hypothesis: true location difference is not equal to c(0,0)Se puede asumir que el vector de medias de las especies JL y LP son iguales
data_SS_JL=subset(data,data$Species!="LP")
summaryBy(X560_nm + X720_nm ~ Species, data=data_SS_JL, FUN=mean)## Species X560_nm.mean X720_nm.mean
## 1 JL 21.55500 46.39250
## 2 SS 11.54333 30.15167df_Species = split(data_SS_JL[,1:2], data_SS_JL$Species)
HotellingsT2(as.matrix(data_SS_JL[,1:2])~as.matrix(data_SS_JL[,3]), mu = c(0,0))##
## Hotelling's two sample T2-test
##
## data: as.matrix(data_SS_JL[, 1:2]) by as.matrix(data_SS_JL[, 3])
## T.2 = 3.2938, df1 = 2, df2 = 21, p-value = 0.05699
## alternative hypothesis: true location difference is not equal to c(0,0)Las medias entre las especies JL y SS es posiblemente igual considerando el p-valor.
Punto 11: MANOVA - intervalos simultáneos
Siguiendo la ecuación:
\[ \require{cancel} T_{ki} -T_{li}=\bar X_{ki} -\bar X_{li} \pm t_{n-g} (\frac{\alpha}{pg(g-1)})\sqrt {\frac{w_{ii}}{n-g}(\frac{1}{n_{k}}+\frac{1}{n_{l}})} \]
alpha= 0.05
p=2
g=3
n=36
gl=n-g
prob= alpha/(p*g*(g-1))
tc=qt(prob,gl, lower.tail = FALSE);tc## [1] 2.806761factor= ((1/12)+(1/12))
#Calculo de W
part=split(data[,1:2],data$Species)
S1= cov(part$JL)
S2= cov(part$LP)
S3= cov(part$SS)
W=(11*S1)+(11*S2)+(11*S3)
W11=W[1,1]
W22=W[2,2]intervalos para la refelctancia de 560 nm
medias_560= tapply(data$X560_nm,data$Species,mean);medias_560## JL LP SS
## 21.555000 9.805833 11.543333Especies JL y LP
medias_1= medias_560[1]
medias_2= medias_560[2]
dif_medias=medias_1-medias_2
ls=dif_medias+(tc*sqrt((W11/(n-g))*factor));ls## JL
## 20.99319li=dif_medias-(tc*sqrt((W11/(n-g))*factor));li## JL
## 2.505142El intervalo se encuentra entre 2.50 y 20.99. por tanto las especies JL y LP son diferentes en la reflectancia espectral 570 nm.
Especies LP y SS
medias_1= medias_560[2]
medias_2= medias_560[3]
dif_medias=medias_1-medias_2;
ls=dif_medias+(tc*sqrt((W11/(n-g))*factor));ls## LP
## 7.506525li=dif_medias-(tc*sqrt((W11/(n-g))*factor));li## LP
## -10.98152El intervalo se encuentra entre -10.98 y7.50 por tanto las especies LP y SS No son diferentes en la reflectancia espectral 570 nm.
Especies SS y JL
medias_1= medias_560[3]
medias_2= medias_560[1]
dif_medias=medias_1-medias_2;
ls=dif_medias+(tc*sqrt((W11/(n-g))*factor));ls## SS
## -0.7676418li=dif_medias-(tc*sqrt((W11/(n-g))*factor));li## SS
## -19.25569El intervalo se encuentra entre -19.26 y -0.77 por tanto las especies SS y JL son diferentes en la reflectancia espectral 570 nm.
intervalos para la refelctancia de 720 nm
medias_720= tapply(data$X720_nm,data$Species,mean);medias_720## JL LP SS
## 46.39250 30.82000 30.15167Especies JL y LP
medias_1= medias_720[1]
medias_2= medias_720[2]
dif_medias=medias_1-medias_2;
ls=dif_medias+(tc*sqrt((W22/(n-g))*factor));ls## JL
## 32.88947li=dif_medias-(tc*sqrt((W22/(n-g))*factor));li## JL
## -1.744474El intervalo se encuentra entre -1.74 y 32.89 por tanto las especies JL y LP No son diferentes en la reflectancia espectral 570 nm.
Especies LP y SS
medias_1= medias_720[2]
medias_2= medias_720[3]
dif_medias=medias_1-medias_2;
ls=dif_medias+(tc*sqrt((W22/(n-g))*factor));ls## LP
## 17.98531li=dif_medias-(tc*sqrt((W22/(n-g))*factor));li## LP
## -16.64864El intervalo se encuentra entre -16.65 y 17.98 por tanto las especies LP y SS No son diferentes en la reflectancia espectral 570 nm.
Especies SS y JL
medias_1= medias_720[3]
medias_2= medias_720[1]
dif_medias=medias_1-medias_2;
ls=dif_medias+(tc*sqrt((W22/(n-g))*factor));ls## SS
## 1.076141li=dif_medias-(tc*sqrt((W22/(n-g))*factor));li## SS
## -33.55781El intervalo se encuentra entre -33-56 y 1.08 por tanto las especies SS y JL No son diferentes en la reflectancia espectral 570 nm
Punto 12: MANOVA - variable de bloqueo
Se incoroporó una variable de bloqueo(dos localidades;L1 y L2) para la mitad de los datos de cada especie y se ejecutó el MANOVA.
loc=c(rep("L1",6),rep("L2",6))
data_loc=cbind(data,loc);data_loc## X560_nm X720_nm Species loc
## 1 9.33 19.14 SS L1
## 2 8.74 19.55 SS L1
## 3 9.31 19.24 SS L1
## 4 8.27 16.37 SS L1
## 5 10.22 25.00 SS L1
## 6 10.13 25.32 SS L1
## 7 10.42 27.12 SS L2
## 8 10.62 26.28 SS L2
## 9 15.25 38.89 SS L2
## 10 16.22 36.67 SS L2
## 11 17.24 40.74 SS L2
## 12 12.77 67.50 SS L2
## 13 12.07 33.03 JL L1
## 14 11.03 32.37 JL L1
## 15 12.48 31.31 JL L1
## 16 12.12 33.33 JL L1
## 17 15.38 40.00 JL L1
## 18 14.21 40.48 JL L1
## 19 9.69 33.90 JL L2
## 20 14.35 40.15 JL L2
## 21 38.71 77.14 JL L2
## 22 44.74 78.57 JL L2
## 23 36.67 71.43 JL L2
## 24 37.21 45.00 JL L2
## 25 8.73 23.27 LP L1
## 26 7.94 20.87 LP L1
## 27 8.37 22.16 LP L1
## 28 7.86 21.78 LP L1
## 29 8.45 26.32 LP L1
## 30 6.79 22.73 LP L1
## 31 8.34 26.67 LP L2
## 32 7.54 24.87 LP L2
## 33 14.04 44.44 LP L2
## 34 13.51 37.93 LP L2
## 35 13.33 37.93 LP L2
## 36 12.77 60.87 LP L2Reflect_join=cbind(data$X560_nm,data$X720_nm)
modM2=manova(Reflect_join~data_loc$Species*data_loc$loc)
summary(modM2, test="Wilks")## Df Wilks approx F num Df den Df Pr(>F)
## data_loc$Species 2 0.52677 5.4783 4 58 0.0008230 ***
## data_loc$loc 1 0.54684 12.0162 2 29 0.0001581 ***
## data_loc$Species:data_loc$loc 2 0.69750 2.8618 4 58 0.0311414 *
## Residuals 30
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1El resultado indica que hay diferencias entre los grupos.