1 Pruebas estadísticas en R

En este curso vamos a revisar cómo hacer algunas pruebas estadísticas paramétricas y no paramétricas utilizando R.

Aquí les pongo una tabla para resumir las funciones para hacer distintas pruebas paramétricas y no paramétricas dependiendo de lo que se quiera comparar.

Variables a comparar	Prueba paramétrica	Prueba no paramétrica
Una media o mediana vs: un valor, otra media o mediana, o media o mediana de datos pareada.	t.test()	wilcox.test()
Más de dos medias o medianas	aov(), TukeyHSD()	kruskal.test()
Dos variables relación	cor()	cor()
Dos variables dependencia	lm()	-

Diferencias entre pruebas paramétricas y no paramétricas

• Las paramétricas comparan medias, las no paramétricas medianas.
• La media se ve muy afectada por valores extremos, mientras que la mediana no.
• Las pruebas paramétricas suponen una distribucióm normal de los datos, las no paramétricas, no.
• Las pruebas paramétricas suelen suponer igualdad de varianzas, las no paramétricas, no.
• Las paramétricas tienen utilidad baja con una n <= 20 - 30, mientras que las no paramétricas, tienen una utilidad alta en esos casos (tiene que ver con la varianza alta de muestras chiquitas).
• Las pruebas paramétricas tienen alto poder estadístico (son más estrictas), mientras que las no paramétricas, tienen poder medio (son más laxas).
• Ambas requieren muestras aleatorias e independientes.

Primero vamos a crear unos datos de prueba para realizar las pruebas. Supongamos que los datos representan mediciones de altura de plantas (cm) puestas bajo dos tratamientos de fertilizante (A y B). De tal manera, tenemos 20 observaciones por tratamiento.

datos <- data.frame(Tratamiento = rep(c("A","B"),each = 20),
                    Altura = c(rnorm(20,10,3),
                                    rnorm(20,20,2)))
datos

##    Tratamiento    Altura
## 1            A 12.852157
## 2            A  3.389011
## 3            A 13.363668
## 4            A 11.076183
## 5            A  8.690731
## 6            A  8.579248
## 7            A  8.781423
## 8            A 12.158028
## 9            A  4.556388
## 10           A  7.370822
## 11           A 10.432980
## 12           A 14.937462
## 13           A  9.467594
## 14           A 15.029689
## 15           A  7.353987
## 16           A  7.485911
## 17           A 11.304386
## 18           A 10.944261
## 19           A  9.823903
## 20           A  9.878039
## 21           B 16.524604
## 22           B 22.943522
## 23           B 21.102774
## 24           B 19.840130
## 25           B 20.920713
## 26           B 20.935202
## 27           B 18.032236
## 28           B 19.144688
## 29           B 22.003282
## 30           B 17.117138
## 31           B 17.867559
## 32           B 22.062799
## 33           B 21.630362
## 34           B 20.891759
## 35           B 18.515907
## 36           B 18.265348
## 37           B 20.398222
## 38           B 21.925150
## 39           B 17.732617
## 40           B 17.515467

1.1 Comparación de dos medias o medianas

1.1.1 Prueba de t (paramétrica)

Primero utilizaremos la prueba de t de student para comparar entre las medias de los dos tratamientos

t.test(datos[,2]~datos[,1])

## 
##  Welch Two Sample t-test
## 
## data:  datos[, 2] by datos[, 1]
## t = -12.23, df = 32.316, p-value = 1.145e-13
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -11.542019  -8.247341
## sample estimates:
## mean in group A mean in group B 
##        9.873793       19.768474

Esta prueba nos dice que el p-value < 0.05, por lo tanto rechazamos la hipótesis nula y concluimos que las medias son distintas entre los dos tratamientos.

1.1.2 Wilcoxon (no paramétrica)

Suponiendo que nuestros datos no se distribuyeran de manera normal podríamos utilizar la prueba de Wilcoxon para comparar entre las medianas de los dos tratamientos.

wilcox.test(datos[,2]~datos[,1])

## 
##  Wilcoxon rank sum test
## 
## data:  datos[, 2] by datos[, 1]
## W = 0, p-value = 1.451e-11
## alternative hypothesis: true location shift is not equal to 0

Esta prueba nos dice que el p-value < 0.05, por lo tanto rechazamos la hipótesis nula y concluimos que las medias son distintas entre los dos tratamientos.

Ahora, para gráficas los resultados podríamos hacer lo siguiente

A<-datos[which(datos[,1]=="A"),]
B<-datos[which(datos[,1]=="B"),]

avg_A<-mean(A[,2])
avg_A

## [1] 9.873793

avg_B<-mean(B[,2])
avg_B

## [1] 19.76847

std.error<-function (x) sd(x)/sqrt(sum(!is.na(x)))

se_A<-std.error(A[,2])
se_A

## [1] 0.681566

se_B<-std.error(B[,2])
se_B

## [1] 0.435915

y.means<-c(avg_A,avg_B)
y.se<-c(se_A,se_B)

#Para hacer una gráfica
barx<-barplot(c(avg_A,avg_B),names.arg = c("A","B"),
        main ="Crecimiento de plantas (cm) bajo dos tratamientos de fertilizante",
        col = c("white","white"),
        ylim=c(0,25))

#Barras de error se pueden poner pero es más complicado
error.bar <- function(x, y, upper, lower=upper, length=0.1,...){
  if(length(x) != length(y) | length(y) !=length(lower) | length(lower) != length(upper))
    stop("vectors must be same length")
  arrows(x,y+ upper, x, y-lower, angle=90, code=3, length=length, ...)
}

error.bar(barx,y.means, 1.96*y.se)

1.1.3 Ejercicio t student

Hacer una prueba de t de student para otra base de datos y concluir si las medias son distintas o no.

Los nuevos datos serían los siguientes

datos <- data.frame(Tratamiento = rep(c("C","D"),each = 20),
                    Altura = c(rnorm(20,30,3),
                                    rnorm(20,5,3)))

1.2 Pruebas sobre supuestos de algunas pruebas paramétricas: Normalidad y homocedasticidad

1.2.1 Levene, prueba de homogeneidad de varianzas

Esta prueba sirve para probar si la varianza es homogenea entre los niveles de un factor. Por ejemplo, en el ejemplo pasado podríamos ver si la varianza es homogenea entre los dos niveles del factor tratamiento, es decir, C y D. Algunas pruebas paramétricas tienen el supuesto de la homogeneidad de varianzas (p.ej., ANOVA), por lo cual es necesario probar si se cumple con este requisito o no. Esta prueba podría servir para decidir si utilizar una prueba paramétrica (p.ej., ANOVA) o no (p.ej., Kruskal-Wallis).

Para hacer esta prueba hay que cargar el paquete “car”.

Vamos a utilizar los mismos datos que en el ejercicio anterior para hacer estas pruebas.

En esta prueba la hipótesis nula es que las varianzas son homogeneas.

library(car)

## Loading required package: carData

leveneTest(datos$Altura~datos$Tratamiento)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.3286 0.5699
##       38

Esto nos dice que los datos son homoscedasticos (i.e., los dos grupos tienen la misma varianza). Debido a que la hipótesis nula es que los datos son homoscedásticos, el valor de significancia (0.9) > 0.5, por lo cual no se rechaza la hipótesis nula.

1.2.2 Shapiro-Wilk, prueba de normalidad

Esta prueba sirve para probar si los datos se distribuyen de manera normal. Este es otro supuesto para las pruebas paramétricas, por lo cual, esta prueba también puede servir para elegir la prueba para analizar los datos (prueba paramétrica o no). En esta prueba la hipótesis nula es que los datos se distribuyen de manera normal.

do.call("rbind",with(datos,tapply(datos$Altura,datos$Tratamiento,function(x) unlist(shapiro.test(x)[c("statistic", "p.value")])))) 

##   statistic.W    p.value
## C   0.8829228 0.01997661
## D   0.9554031 0.45656139

Esto nos dice que los dos niveles del factor tienen distribución normal, ya que su valor de significancia (p-value) > 0.5.

1.3 Comparación de más de dos medias o medianas

En esta sección veremos las pruebas para comparar más de dos medias o medianas. Primero, crearemos nuestra base de datos. En este caso supongamos que tenemos los datos de área basal de tres comunidades distintas de plantas. En cada comunidad se muestrearon 40 parcelas para obtener estos datos. El área basal de cada planta se calcula como el área de un círculo cuyo diámetro es igual al diámetro del tallo. Para calcular el área basal de la comunidad se suman todos los valores y usualmente se estandarizan a una ha. De tal manera, los datos de área basal se encuentran en m^2 / ha.

datos <- data.frame(Comunidad = rep(c("Com1","Com2","Com3"),each = 40),
                    AB = c(rnorm(40,35,4),
                               rnorm(40,42,8),
                               rnorm(40,25,10)))

1.3.1 ANOVA (paramétrica)

Primero cargamos los paquetes que vamos a utilizar

#Cargar los paquetes necesarios
library("stats")# aov y shapiro y kruskal-wallis
library("car")#levene
library("pgirmess")#tukey no paramétrico (kruskalmc)

Recordar que la hipótesis nula del ANOVA es que las medias son iguales.

#Anova
anova_com<-aov(datos$AB~datos$Comunidad)
summary(anova_com)

##                  Df Sum Sq Mean Sq F value   Pr(>F)    
## datos$Comunidad   2   3734  1867.1   28.55 7.98e-11 ***
## Residuals       117   7651    65.4                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Debido a que el valor de significancia es <= 0.5 rechazamos la hipótesis nula y concluimos que no todas las medias son iguales (es decir, que por lo menos hay una distinta). Si el ANOVA sale significativo, después quisieramos ver la media de cuáles grupos difiere entre sí y cuáles no. Para eso vamos a utilizar la prueba de Tukey.

1.3.2 TukeyHSD

Esta prueba nos sirve para hacer comparaciones pareadas entre las combinaciones de niveles de un factor. La hipótesis nula es que las medias comparadas son iguales.

TukeyHSD(anova_com)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = datos$AB ~ datos$Comunidad)
## 
## $`datos$Comunidad`
##                 diff         lwr       upr     p adj
## Com2-Com1   4.792443   0.4998027  9.085083 0.0246536
## Com3-Com1  -8.685531 -12.9781715 -4.392891 0.0000139
## Com3-Com2 -13.477974 -17.7706145 -9.185334 0.0000000

Las tres comparaciones tienen un valor de significancia <= 0.5 por lo cual, rechazamos la hipótesis nula y concluimos que las tres medias son distintas entre sí.

1.3.3 Kruskal-Wallis (no paramétrica)

El equivalente no paramétrico del ANOVA es la prueba de Kruskal-Wallis. La hipótesis nula de esta prueba es que las medianas son iguales

kruskal_pais<-kruskal.test(datos$AB~datos$Comunidad)
kruskal_pais

## 
##  Kruskal-Wallis rank sum test
## 
## data:  datos$AB by datos$Comunidad
## Kruskal-Wallis chi-squared = 34.645, df = 2, p-value = 2.998e-08

La prueba nos arroja un valor de significancia <= 0.5 por lo cual rechazamos la hipótesis nula y concluimos que al menos una mediana es distinta a las demás. Para saber cuáles medianas son distintas entre sí realizamos la siguiente prueba de comparación multiple después de realizar una prueba de Kruskal-Wallis.

1.3.4 TukeyHSD no paramétrico

# Para leer la ayuda de la función
??kruskalmc

## starting httpd help server ... done

tukey.no.par<-kruskalmc(datos$AB~datos$Comunidad)
tukey.no.par

## Multiple comparison test after Kruskal-Wallis 
## p.value: 0.05 
## Comparisons
##           obs.dif critical.dif difference
## Com1-Com2  16.750     18.62079      FALSE
## Com1-Com3  28.525     18.62079       TRUE
## Com2-Com3  45.275     18.62079       TRUE

En este caso nos muestra las diferencias entre los grupos, así como la diferencia que la prueba toma como crítica. En esta caso las medianas de la comunidad 1 y 3, así como la 2 y 3 sí difieren significativamente, mientras que la comunidad 1 y 2 no.

Al comparar los resultados con la prueba paramétrica y no paramétrica nos puede llamar la atención el por qué con la prueba paramétrica vemos diferencias significativas entre todos lod grupos, mientras que en la no paramétrica, no. Esto tiene que ver con el hecho de que las pruebas paramétricas tienen mayor poder estadístico. Por lo tanto, es recomendable que si los datos cumplen con los requisitos de las pruebas paramétricas (comentados al principio) se utilicen estas pruebas.

Si graficamos directamente la prueba del ANOVA vemos los siguientes gráficos, sobre la distribución de los residuos; sin embargo, estas no suelen ser las gráficas utilizadas para reportar un ANOVA, por lo cual, se hará lo siguiente.

plot(anova_com)

1.3.5 ANOVA gráfico

C1<-datos$AB[which(datos[,1]=="Com1")]
C2<-datos$AB[which(datos[,1]=="Com2")]
C3<-datos$AB[which(datos[,1]=="Com3")]

C_todos<-cbind(C1,C2,C3)

medias_C_todos<-apply(C_todos,2,mean)

std.error<-function (x) sd(x)/sqrt(sum(!is.na(x)))

se_C_todos<-apply(C_todos,2,std.error)

barx<-barplot(medias_C_todos,names.arg = c("Com1","Com2","Com2"),
        main ="Área basal (m^2 / ha) de tres comunidades vegetales",
        col = rep("white",3),
        ylim=c(0,50))

error.bar <- function(x, y, upper, lower=upper, length=0.1,...){
  if(length(x) != length(y) | length(y) !=length(lower) | length(lower) != length(upper))
    stop("vectors must be same length")
  arrows(x,y+ upper, x, y-lower, angle=90, code=3, length=length, ...)
}

error.bar(barx,medias_C_todos, 1.96*se_C_todos)

1.3.6 Ejercicio ANOVA

• Hacer una misma prueba con los siguientes datos

datos <- data.frame(Comunidad = rep(c("Com1","Com2","Com3"),each = 40),
                    AB = c(rnorm(40,35,7),
                               rnorm(40,30,8),
                               rnorm(40,33,10)))

1.4 Relación entre dos variables continuas

En este caso vamos a ver las pruebas paramétricas y no paramétricas para analizar la relación entre dos variables continuas.

En este caso vamos a utilizar la base de datos que se uso en el curso básico de R: https://github.com/JonathanVSV/JonathanVSV.github.io/blob/master/Curso_basico_R/BD.csv

Al entrar a la página presionar la tecla Alt y dar click en el botón de “Raw” para descargarlo a su disco duro.

1.4.1 Correlación

En este caso vamos a analizar la relación entre la producción de cereales y la población total. La hipótesis nula de esta prueba es que no existe relación alguna entre las variables. Para leer más información de los datos: https://data.worldbank.org/.

datos<-read.csv("BD.csv")

cor_data<-datos[,c(7,12)]

colnames(cor_data)<-c("Prod_Cereales","Poblacion_total")

cor(cor_data)

##                 Prod_Cereales Poblacion_total
## Prod_Cereales       1.0000000       0.9582025
## Poblacion_total     0.9582025       1.0000000

En este caso vemos que existe una relación fuerte entre las dos variables analizadas. Sin embargo, para calcula la significancia de la prueba realizamos lo siguiente:

cor.test(cor_data[,1],cor_data[,2])

## 
##  Pearson's product-moment correlation
## 
## data:  cor_data[, 1] and cor_data[, 2]
## t = 42.366, df = 160, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9433939 0.9691984
## sample estimates:
##       cor 
## 0.9582025

En este caso la significancia es <= 0.5 por lo que rechazamos la hipótesis nula y concluimos que sí existe una relación entre las dos variables.

1.4.2 Gráfica correlación

1.4.3 Otros índices de correlación

Por default, la función cor y cor.test utilizan la correlación de Pearson. A pesar de ello, se puede indicar otra función para calculat la correlación entre variables (p.ej., Spearman o Kendall). Ver la ayuda de la función cor para ver cómo hacerlo.

help(cor)

¿Qué dice?

1.4.4 Ejercicio correlación

• Hacer un análisis de correlación con otras dos variables de la hoja de datos (GDP y producción de ganado).
• Contestar si las variables están correlacionadas significativamente o no

1.5 Causalidad entre dos variables continuas

Para analizar la causalidad entre dos variables continuas se suele utilizar el análisis de regresión. Para ello, vamos a utilizar los mismos datos que en el ejercicio pasado, sólo que ahora usaremos las columnas de año y emisiones de CO2 y sólo usaremos los datos de México.

1.5.1 Regresión lineal

En el análisis de regresión lineal hay dos hipótesis nulas: 1) que el intercepto de la recta es = 0 y 2) que la pendiente de la recta = 0.

datos<-read.csv("BD.csv")

lm_data<-datos[,c(1,2,3)]
colnames(lm_data)<-c("Country","Year","C02")
lm_data<-subset(lm_data,Country == "Mexico")

modelo_lineal<-lm(lm_data$C02~lm_data$Year)

summary(modelo_lineal)

## 
## Call:
## lm(formula = lm_data$C02 ~ lm_data$Year)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -48464 -17398  -2655  14815  61524 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -1.743e+07  3.909e+05  -44.60   <2e-16 ***
## lm_data$Year  8.919e+03  1.967e+02   45.35   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 22530 on 52 degrees of freedom
## Multiple R-squared:  0.9753, Adjusted R-squared:  0.9749 
## F-statistic:  2056 on 1 and 52 DF,  p-value: < 2.2e-16

En este caso vemos que se rechaza la hipótesis nula tanto del intercepto como de la pendiente. Por ello, concluimos que 1) el intercepto de la recta es distinto de 0 y 2) que la pendiente de la recta es distinta de 0. Este último implica que el efecto de la variable independiente (año) sobre la variable dependiente (C02) es significativa.

1.5.2 Gráfica regresión lineal

plot(lm_data$Year,lm_data$C02,type="p",ylab="C02",xlab = "Year")
abline(modelo_lineal$coefficients,col="red")

1.5.3 Regresión lineal

• En r también puedes hacer glm, similar al lm, pero con la función glm

1.5.4 Ejercicio regresión lineal

• Hacer un análisis de regresión lineal para ver si GDP depende de la producción agrícola (crop).
• Concluir cuál es su R^2, si es significativa la relación, si el intercepto es distinto de 0.

2 Estadística multivariada

Ahora veremos rápidamente cómo realizar algunas de las pruebas más frecuentes utilizadas en estadística multivariada.

Algunos de las funciones para realizar análisis de estadística multivariada más comunes incluyen:

• Análisis RTQ
  • cor()
  • dist()
• Cluster
  • dist()
  • hclust()
  • cutree()
• Ordenación
  • pca()
  • cca()
  • cmdscale()
  • metaMDS()
  • decorana()
  • procrustes()
  • metaMDS()
  • decorana()
  • procrustes()
• RDA
  • rda()
• RDA parcial
  • var.part()

2.1 Análisis de Cluster

Esta análisis de agrupamiento nos va a decir cuáles unidades muestreales se parecen más entre sí de acuerdo a un conjunto de variables. En este caso vamos a utilizar los datos de varespec que vienen en el paquete vegan. Para leer sobre los datos utilicen ??varespec. Los datos contienen la información de composición de especies de 24 sitios. En este caso, vamos a calcular índices de disimilitud entre los sitios (de acuerdo a su composición) y después se realiza el análisis de cluster o aglomeramiento.

library(vegan)
library(cluster)

data(varespec)

datos_clust<-varespec

# Cálculo de índices de disimilitud
ssa<-vegdist(datos_clust,method="euclidean",binary=FALSE,diag=FALSE)
# Análisis de aglomeramiento
Tsa<-agnes(ssa,diss=F,metric="euclidean",method="average")

2.1.1 Gráfica cluster

Para obtener el gráfico, utilizamos los siguientes comandos. En este caso, el eje de las Y muestra la disimilitud; por lo tanto, los sitios que se encuentran aglomerados en valores bajos de disimilitud corresponden a los sitios más parecidos entre sí.

2.1.2 Cortar dendrograma

A partir del dendrograma se puede definir un número de grupos, que se podría interpretar como tipos de comunidades o asociaciones entre especies. De tal manera, podemos asignarle un grupo a cada sitio. Esto equivale a tomar un valor de dismilitud de corte que arroje una división de los sitios en 4 grupos. En este caso, se podría interpretar que se tomó un valor de corte de aprox. 110.

grupos<-cutree(Tsa,4)
agrupamiento<-data.frame(Sitio=row.names(datos_clust),Grupo=grupos)
agrupamiento

##    Sitio Grupo
## 1     18     1
## 2     15     2
## 3     24     2
## 4     27     3
## 5     23     2
## 6     19     2
## 7     22     2
## 8     16     2
## 9     28     3
## 10    13     1
## 11    14     2
## 12    20     2
## 13    25     2
## 14     7     1
## 15     5     1
## 16     6     1
## 17     3     4
## 18     4     1
## 19     2     4
## 20     9     4
## 21    12     4
## 22    10     4
## 23    11     1
## 24    21     2

2.2 Ordenaciones

• Consta de una serie de análisis multivariados que buscan analizar la redundancia de una serie de variables explicativas.
• Busca explicar la mayor proporción de la variación observada en función de ejes principales.
• A continuación permite evaluar cuáles variables son las más correlacionadas entre sí y en comparación con los ejes principales.

2.2.1 Tipos de ordenación

• Basados en distancia
  • directos
  • indirectos
• Basados en eigenvectores
  • directos
    • explicativos
  • indirectos
    • exploratorios

2.2.2 RDA

Similar a PCA pero es un método directo o con restriciones (constrained).

Vamos a utilizar los datos del 2014 de la base de datos del “world bank”.

library(vegan)

datos<-read.csv("BD.csv")
datos<-subset(datos,datos$Year==2014)

#Scale =T para cuando tienes variables en diferentes unidades
pca.atributos<-rda(datos[,c(3:12)]~datos$Country,scale=T)
summary(pca.atributos)

## 
## Call:
## rda(formula = datos[, c(3:12)] ~ datos$Country, scale = T) 
## 
## Partitioning of correlations:
##               Inertia Proportion
## Total              10          1
## Constrained        10          1
## Unconstrained       0          0
## 
## Eigenvalues, and their contribution to the correlations 
## 
## Importance of components:
##                         RDA1   RDA2
## Eigenvalue            8.8711 1.1289
## Proportion Explained  0.8871 0.1129
## Cumulative Proportion 0.8871 1.0000
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                         RDA1   RDA2
## Eigenvalue            8.8711 1.1289
## Proportion Explained  0.8871 0.1129
## Cumulative Proportion 0.8871 1.0000
## 
## 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:  2.114743 
## 
## 
## Species scores
## 
##                                                 RDA1      RDA2
## CO2.emissions..kt.                            0.6106 -0.272743
## Livestock.production.index..2004.2006...100. -0.6687 -0.001619
## Food.production.index..2004.2006...100.      -0.5868 -0.320824
## Crop.production.index..2004.2006...100.      -0.6219 -0.245941
## Cereal.production..metric.tons.               0.6079 -0.278685
## GDP.per.capita..current.US..                  0.6216 -0.246693
## Life.expectancy.at.birth..total..years.       0.6549  0.135545
## Urban.population                              0.6622 -0.092946
## Urban.population....of.total.                 0.5953  0.304754
## Population..total                             0.6622 -0.093022
## 
## 
## Site scores (weighted sums of species scores)
## 
##         RDA1    RDA2
## 54  -0.03996  1.7262
## 108 -1.47497 -0.8977
## 162  1.51493 -0.8285
## 
## 
## Site constraints (linear combinations of constraining variables)
## 
##         RDA1    RDA2
## 54  -0.03996  1.7262
## 108 -1.47497 -0.8977
## 162  1.51493 -0.8285
## 
## 
## Biplot scores for constraining variables
## 
##                         RDA1    RDA2
## datos$CountryMexico -0.02314  0.9997
## datos$CountryUSA     0.87736 -0.4798
## 
## 
## Centroids for factor constraints
## 
##                            RDA1    RDA2
## datos$CountryGuatemala -1.47497 -0.8977
## datos$CountryMexico    -0.03996  1.7262
## datos$CountryUSA        1.51493 -0.8285

Aquí podemos ver la proporción de la varianza explicada por el modelo: 0.88, así como los eigenvalores y los valores de cada variable y país. Para entenderlo mejor, vamos a graficarlo.

2.2.2.1 Gráfica RDA

Aquí podemos ver, por ejemplo, que los indices de producción de comida, cultivos agrícolas y ganadería son mayores en Guatemala que en México y EUA. Por el contrario, la población urbana es mayor es EUA que en México y Guatemala.

plot(pca.atributos)

2.2.2.2 Ejercicio RDA

Utilizar los datos de varespec

• Ver de que se trata la hoja de datos
• Hacer un pca de sitios en función de las especies presentes

data(varespec)
help(varespec)

pca_varespec<-rda(varespec)
plot(pca_varespec)

Aquí podemos ver que los sitios 9, 10, 12 y 2 son los que tienen mayor presencia de la especies Cladstel, mientras que los sitios 27, 28, dePleusch.

2.2.2.3 Ejercicio RDA Scaling 1

Aquí vamos a utilizar la opción de scaling = 1, para interpretar las relaciones entre los sitios. En este ejemplo, podemos interpretar que los sitios 9, 10, 12, 2 se parecen más entre sí; mientras que 27 y 28 se parecen también entre sí. Por otro lado, los primeros se encuentran en el extremo positivo del eje principal 1 y en el lado negativo del eje principal 2.

pca_varespec<-rda(varespec)
#scaling=1: Plano de distancias. Relación de distancias se mantiene, sirve para 
#interpretar objetos (renglones), i.e., sitios.
biplot(pca_varespec,scaling=1)

2.2.2.4 Ejercicio RDA Scaling 2

Aquí vamos a utilizar la opción de scaling = 2, para interpretar las relaciones entre las especies. En esta caso podemos ver que las especies Cladrang y Clarbu están muy relacionadas con el eje principal 2, mientras que Cladstel con el eje 1. Por otro lado, las especies Cladstel, Cladarang, Cladarbu y Pleuschr son las más importantes para explicar la variación presente en los datos.

pca_varespec<-rda(varespec)
#scaling=2: Plano de correlaciones. Relación de correlaciones se mantiene, sirve 
#para interpretar atributos (columnas), i.e., especies.
biplot(pca_varespec,scaling=2)

3 Recuerden

• Guarden todos sus códigos y reúsenlos para crear códigos cada vez más complejos.
• En sus códigos utilicen nombres de variables que faciliten el reuso de las rutinas.
• Practiquen frecuentemente la programación en R.
• Utilicen la ayuda (help) o el sitio stackoverflow.com.
• Este fue un curso introductorio, muchos temas y trucos les quedan pendientes por aprender.