Queremos comparar la frecuencia observada para el número de individuos presentes de siete especies de invertebrados en cinco sitios diferentes. Para ello compararemos la frecuencia observada (la real) respecto a la frecuencia esperada. Se parte de la hipótesis nula que las frecuencias de las especies no varian entre los sitios, siendo la hipótesis alternativa la diferencia entre sitios. Si nuestro test Chi cuadrado presenta un p valor inferior a 0.05, esto nos indicará la variación de frecuencias entre los sitios (es decir no todas las especies presentan las mismas frecuencias entre los diferentes sitios).
Cargamos nuestra matriz de datos, que tendrá el siguiente aspecto:
setwd(dir = "D:/R/MARKDOWN/chicuadrado")
tabla <- read.table(file = "chi.csv", sep = ";", header=TRUE, row.names=1)
tabla
## SiteA SiteB SiteC SiteD SiteE
## spA 39 25 128 13 0
## spB 13 33 13 48 23
## spC 146 98 5 21 4
## spD 32 8 2 25 2
## sdE 6 105 10 5 2
## spF 0 25 14 15 1
## spG 12 5 9 0 12
str(tabla)
## 'data.frame': 7 obs. of 5 variables:
## $ SiteA: int 39 13 146 32 6 0 12
## $ SiteB: int 25 33 98 8 105 25 5
## $ SiteC: int 128 13 5 2 10 14 9
## $ SiteD: int 13 48 21 25 5 15 0
## $ SiteE: int 0 23 4 2 2 1 12
En este caso tenemos siete especies y cinco sitios.
Pasamos nuestros datos a una tabla:
comotabla <- as.table(as.matrix(tabla))
comotabla
## SiteA SiteB SiteC SiteD SiteE
## spA 39 25 128 13 0
## spB 13 33 13 48 23
## spC 146 98 5 21 4
## spD 32 8 2 25 2
## sdE 6 105 10 5 2
## spF 0 25 14 15 1
## spG 12 5 9 0 12
Ahora representamos gráficamente la tabla por medio de círculos. Intalamos install.packages("gplots")
y lo activamos. A mayor diámetro mayor frecuencia de individuos de la especie en un sitio:
library("gplots")
## Warning: package 'gplots' was built under R version 3.4.4
##
## Attaching package: 'gplots'
## The following object is masked from 'package:stats':
##
## lowess
balloonplot(t(comotabla), main ="Especies vs Sitios", xlab ="", ylab="",
label = FALSE, show.margins = FALSE, dotsize=6)
Hacemos el test, representamos los valores observados y los esperados, y los residuos (a mayor valor más peso de la especie en ese sitio, tanto en + como en -):
test <- chisq.test(tabla)
## Warning in chisq.test(tabla): Chi-squared approximation may be incorrect
test
##
## Pearson's Chi-squared test
##
## data: tabla
## X-squared = 758.25, df = 24, p-value < 2.2e-16
test$observed
## SiteA SiteB SiteC SiteD SiteE
## spA 39 25 128 13 0
## spB 13 33 13 48 23
## spC 146 98 5 21 4
## spD 32 8 2 25 2
## sdE 6 105 10 5 2
## spF 0 25 14 15 1
## spG 12 5 9 0 12
test$expected
## SiteA SiteB SiteC SiteD SiteE
## spA 56.55172 68.18131 41.273637 28.959956 10.033370
## spB 35.86207 43.23693 26.173526 18.364850 6.362625
## spC 75.58621 91.13014 55.165740 38.707453 13.410456
## spD 19.03448 22.94883 13.892102 9.747497 3.377086
## sdE 35.31034 42.57175 25.770857 18.082314 6.264739
## spF 15.17241 18.29255 11.073415 7.769744 2.691880
## spG 10.48276 12.63849 7.650723 5.368187 1.859844
round(test$expected,1)
## SiteA SiteB SiteC SiteD SiteE
## spA 56.6 68.2 41.3 29.0 10.0
## spB 35.9 43.2 26.2 18.4 6.4
## spC 75.6 91.1 55.2 38.7 13.4
## spD 19.0 22.9 13.9 9.7 3.4
## sdE 35.3 42.6 25.8 18.1 6.3
## spF 15.2 18.3 11.1 7.8 2.7
## spG 10.5 12.6 7.7 5.4 1.9
round(test$residuals,2)
## SiteA SiteB SiteC SiteD SiteE
## spA -2.33 -5.23 13.50 -2.97 -3.17
## spB -3.82 -1.56 -2.57 6.92 6.60
## spC 8.10 0.72 -6.75 -2.85 -2.57
## spD 2.97 -3.12 -3.19 4.89 -0.75
## sdE -4.93 9.57 -3.11 -3.08 -1.70
## spF -3.90 1.57 0.88 2.59 -1.03
## spG 0.47 -2.15 0.49 -2.32 7.44
El p valor es muy bajo, nos indica diferencias en las frecuencias de las especies entre los diferentes sitios.
Representamos visualmente los residuos. Instalamos el paquete corrplot
install.packages("corrplot")
:
library(corrplot)
## Warning: package 'corrplot' was built under R version 3.4.4
## corrplot 0.84 loaded
corrplot(test$residuals, is.cor = FALSE)
Valores azules (positivos) indican asociación entre esa especie y el sitio correspondiente, mientras que valores claros indican que la especie tiende a no ocupar ese sitio.
Ahora vamos a representar el % de contribución de cada celda al total del valor del chi cuadrado, y lo representamos gráficamente:
contrib <- 100*test$residuals^2/test$statistic
round(contrib, 3)
## SiteA SiteB SiteC SiteD SiteE
## spA 0.718 3.607 24.034 1.160 1.323
## spB 1.922 0.320 0.874 6.307 5.737
## spC 8.651 0.068 6.016 1.068 0.871
## spD 1.165 1.284 1.343 3.148 0.074
## sdE 3.209 12.073 1.273 1.248 0.383
## spF 2.001 0.324 0.102 0.887 0.140
## spG 0.029 0.609 0.031 0.708 7.291
corrplot(contrib, is.cor = FALSE)
En general globos más grandes indican mayor contribución de esa especie al estadístico chi.