2017-05-18

Estadística en R

  • R básico tiene funciones para las pruebas mas comunes de estadística:
    • correlación
    • regresión
    • chi-cuadrado …

Estadística en R

  • Gran cantidad de paquetes para (casi) cualquier método estadístico alternativo:
    • modelos lineales mixtos (lme4)
    • ecología de comunidades (vegan)
    • selección de modelos (AICcmodavg, MuMIn)
    • modelos lineales y estadística multivariada (MASS)

Datos de ejemplo

Juego de datos:

a <- "http://marceloarayasalas.weebly.com/uploads"
b <- "/2/5/5/2/25524573/datos_phaethornis___diastema.rds"

download.file(paste0(a, b), destfile = "datos_Phaethornis_&_diastema.RDS")

datos <- readRDS("datos_Phaethornis_&_diastema.RDS")

phaethornis <- datos[[1]]

diastema <- datos[[2]]

Datos de ejemplo

Estadística en R

Juego de datos Phaethornis longirostris

str(phaethornis)
## 'data.frame':    121 obs. of  6 variables:
##  $ coord.sing.event: int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Lek             : Factor w/ 4 levels "CCE","CCL","LOC",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ pattern         : Factor w/ 3 levels "No.ovlp","No.pattern",..: 3 3 3 2 2 1 2 3 3 3 ...
##  $ duetdur         : num  11.9 33.8 48.7 31.1 45.8 ...
##  $ sing.perf       : num  0.405 0.454 0.23 0.287 0.217 ...
##  $ distcat         : Factor w/ 2 levels "Close","Far": NA NA NA NA NA 1 NA NA NA NA ...

Estadística en R

Juego de datos Diastema spp

str(diastema)
## 'data.frame':    901 obs. of  25 variables:
##  $ species    : Factor w/ 3 levels "diastema","vocator",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ sound.files: Factor w/ 216 levels "Dd_AguaBuena_01.wav",..: 1 1 1 1 1 1 1 1 2 2 ...
##  $ selec      : chr  "2" "3" "4" "6" ...
##  $ duration   : num  0.119 0.118 0.12 0.124 0.121 ...
##  $ meanfreq   : num  3.5 3.51 3.48 3.48 3.41 ...
##  $ sd         : num  0.45 0.468 0.461 0.5 0.429 ...
##  $ median     : num  3.37 3.38 3.31 3.27 3.22 ...
##  $ Q25        : num  3.27 3.25 3.24 3.22 3.21 ...
##  $ Q75        : num  3.57 3.6 3.56 3.54 3.49 ...
##  $ IQR        : num  0.302 0.349 0.317 0.323 0.273 ...
##  $ skew       : num  5.93 6.83 7.3 8.6 10.49 ...
##  $ kurt       : num  44.4 61.9 66.8 91.6 122.2 ...
##  $ sp.ent     : num  0.775 0.798 0.775 0.76 0.694 ...
##  $ sfm        : num  0.252 0.297 0.279 0.256 0.23 ...
##  $ mode       : num  3.26 3.25 3.25 3.23 3.21 ...
##  $ centroid   : num  3.5 3.51 3.48 3.48 3.41 ...
##  $ peakf      : num  3.2 3.21 3.17 3.16 3.17 ...
##  $ meandom    : num  3.32 3.28 3.23 3.25 3.25 ...
##  $ mindom     : num  2.15 3.19 3.19 3.19 3.14 ...
##  $ maxdom     : num  4.74 3.49 3.27 3.4 3.49 ...
##  $ dfrange    : num  2.584 0.3015 0.0861 0.2153 0.3445 ...
##  $ modindx    : num  0.258 0.25 0.25 0.25 0.25 ...
##  $ startdom   : num  2.15 3.49 3.27 3.4 3.49 ...
##  $ enddom     : num  4.74 3.19 3.19 3.19 3.14 ...
##  $ dfslope    : num  21.648 -2.563 -0.716 -1.732 -2.843 ...

Cargar paquetes

Mínimo

pqts <- c("corrplot", "MASS", "multcomp", "ggplot2", 
    "vegan", "MuMIn")

A <- lapply(pqts, function(y) {
    if (!y %in% installed.packages()[, "Package"]) 
        install.packages(y)
    require(y, character.only = T)
})
## Loading required package: corrplot
## Loading required package: MASS
## Loading required package: multcomp
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## 
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
## 
##     geyser
## Loading required package: ggplot2
## Loading required package: vegan
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.4-1
## Loading required package: MuMIn

Estadística descriptiva

Mínimo

x <- phaethornis$duetdur

min(x)
## [1] 5.103

Máximo

max(x)
## [1] 438.772

Estadística descriptiva

Rango

range(x)
## [1]   5.103 438.772

Varianza

var(x)
## [1] 4449.565

Estadística descriptiva

Desviación estándar

sd(x)
## [1] 66.70506

Error estándar

sd(x)/sqrt(length(x))
## [1] 6.064097

Estadística descriptiva

Quartiles

quantile(x)
##      0%     25%     50%     75%    100% 
##   5.103  23.133  45.970  90.619 438.772
quantile(x, probs = c(0.025, 0.975))
##    2.5%   97.5% 
##   5.801 236.978

Estadística descriptiva

Resumen

summary(x)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   5.103  23.130  45.970  68.700  90.620 438.800

Estadística descriptiva

Resumen

summary(phaethornis)
##  coord.sing.event  Lek           pattern      duetdur       
##  Min.   :  1      CCE:28   No.ovlp   :29   Min.   :  5.103  
##  1st Qu.: 31      CCL:22   No.pattern:65   1st Qu.: 23.133  
##  Median : 61      LOC:41   Ovlp      :27   Median : 45.970  
##  Mean   : 61      SUR:30                   Mean   : 68.699  
##  3rd Qu.: 91                               3rd Qu.: 90.619  
##  Max.   :121                               Max.   :438.772  
##    sing.perf         distcat  
##  Min.   :-1.00000   Close:18  
##  1st Qu.:-0.30568   Far  :26  
##  Median : 0.06376   NA's :77  
##  Mean   :-0.04186             
##  3rd Qu.: 0.25334             
##  Max.   : 1.06368

Estadística descriptiva

Frecuencias (cuentas)

table(phaethornis$Lek)
## 
## CCE CCL LOC SUR 
##  28  22  41  30

Revisión de supuestos

Normalidad

shapiro.test(x)
## 
##  Shapiro-Wilk normality test
## 
## data:  x
## W = 0.77516, p-value = 2.482e-12

Revisión de supuestos

Normalidad

mi.paleta <- colorRampPalette(c("#E37222", "#66B9BF"))

cols <- adjustcolor(mi.paleta(10), alpha.f = 0.4)

hist(x, breaks = 8, col = sample(cols, 1))

Revisión de supuestos

Normalidad

Revisión de supuestos

Normalidad (ggplot)

library(ggplot2)

ggplot(phaethornis, aes(x = duetdur)) + 
geom_histogram(color = sample(cols, 1), fill = sample(cols, 
    1)) + 
theme_classic()

Revisión de supuestos

Normalidad

qqnorm(x, col = sample(cols, 1), pch = 20, cex = 2)
qqline(x, col = 2)

Revisión de supuestos

Normalidad (ggplot)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Revisión de supuestos

Normalidad y transformación de variables

shapiro.test(log(x))
## 
##  Shapiro-Wilk normality test
## 
## data:  log(x)
## W = 0.98882, p-value = 0.4284

Revisión de supuestos

Normalidad y transformación de variables

hist(log(x), breaks = 15, col = sample(cols, 1))

Revisión de supuestos

Normalidad y transformación de variables

qqnorm(log(x), col = sample(cols, 1), pch = 20, cex = 2)
qqline(log(x), col = 2)

Revisión de supuestos

Igualdad de las varianzas

phaet2 <- phaethornis[phaethornis$Lek %in% c("SUR", 
    "CCL"), ]

var.test(phaet2$sing.perf ~ phaet2$Lek)
## 
##  F test to compare two variances
## 
## data:  phaet2$sing.perf by phaet2$Lek
## F = 1.3138, num df = 21, denom df = 29, p-value = 0.4887
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.5981496 3.0434723
## sample estimates:
## ratio of variances 
##           1.313809

Prueba de hipótesis

  • Prueba T de Student (de 2 muestras)
    • Comparar dos variables numéricas continuas
  • Supuestos
    • Observaciones (o pares) son independientes
    • Las variables siguen una distribución normal

Prueba de hipótesis

Prueba T de Student (de 2 muestras):

data(sleep)

t.test(extra ~ group, data = sleep)
## 
##  Welch Two Sample t-test
## 
## data:  extra by group
## t = -1.8608, df = 17.776, p-value = 0.07939
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -3.3654832  0.2054832
## sample estimates:
## mean in group 1 mean in group 2 
##            0.75            2.33

Prueba de hipótesis

Prueba T de Student (de 2 muestras):

x <- sleep$extra[sleep$group == 1]
y <- sleep$extra[sleep$group == 2]

t.test(x, y)
## 
##  Welch Two Sample t-test
## 
## data:  x and y
## t = -1.8608, df = 17.776, p-value = 0.07939
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -3.3654832  0.2054832
## sample estimates:
## mean of x mean of y 
##      0.75      2.33

Prueba de hipótesis

T pareada

t.test(extra ~ group, data = sleep, paired = T)
## 
##  Paired t-test
## 
## data:  extra by group
## t = -4.0621, df = 9, p-value = 0.002833
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.4598858 -0.7001142
## sample estimates:
## mean of the differences 
##                   -1.58

Prueba de hipótesis

  • Prueba de Wilcoxon
    • Comparar dos variables numéricas continuas
    • Equivalente no paramétrico de la T de Student
  • Supuestos
    • Observaciones (o pares) son independientes

Prueba de hipótesis

  • Prueba de Wilcoxon
wilcox.test(x = x, y = y)
## Warning in wilcox.test.default(x = x, y = y): cannot compute exact p-value
## with ties
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  x and y
## W = 25.5, p-value = 0.06933
## alternative hypothesis: true location shift is not equal to 0

Prueba de hipótesis

Prueba de Wilcoxon pareada

wilcox.test(x = x, y = y, paired = TRUE)
## Warning in wilcox.test.default(x = x, y = y, paired = TRUE): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(x = x, y = y, paired = TRUE): cannot compute
## exact p-value with zeroes
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  x and y
## V = 0, p-value = 0.009091
## alternative hypothesis: true location shift is not equal to 0

Prueba de hipótesis

  • Correlación de Pearson
    • Mide la dependencia lineal entre 2 variables
  • Supuestos
    • Normalidad
    • Relación es lineal
    • Homocedasticidad

Prueba de hipótesis

Prueba de hipótesis

Correlación de Pearson

x2 <- diastema$meanfreq
y2 <- diastema$duration

cor(x = x2, y = y2)
## [1] -0.5760356
cor.test(x = x2, y = y2)
## 
##  Pearson's product-moment correlation
## 
## data:  x2 and y2
## t = -21.129, df = 899, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.6180935 -0.5306894
## sample estimates:
##        cor 
## -0.5760356

Prueba de hipótesis

  • Correlación de Spearman (no paramétrica)
    • Alternativa no paramétrica a pearson

Prueba de hipótesis

cor.test(x = x2, y = y2, method = "spearman")
## Warning in cor.test.default(x = x2, y = y2, method = "spearman"): Cannot
## compute exact p-value with ties
## 
##  Spearman's rank correlation rho
## 
## data:  x2 and y2
## S = 185300000, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.5200349

Matrices de correlación

Evaluar colinealidad de variables

mat.cor <- cor(diastema[, sapply(diastema, is.numeric)])

head(mat.cor)
##            duration   meanfreq          sd     median        Q25
## duration  1.0000000 -0.5760356  0.56732529 -0.5592035 -0.5771093
## meanfreq -0.5760356  1.0000000 -0.18767105  0.9936121  0.9921899
## sd        0.5673253 -0.1876710  1.00000000 -0.1491603 -0.2092945
## median   -0.5592035  0.9936121 -0.14916026  1.0000000  0.9909145
## Q25      -0.5771093  0.9921899 -0.20929449  0.9909145  1.0000000
## Q75      -0.4985858  0.9880212 -0.07466188  0.9905400  0.9775198
##                  Q75        IQR       skew       kurt     sp.ent       sfm
## duration -0.49858583 0.08962441  0.5730427  0.4960304 0.08517835 0.3707118
## meanfreq  0.98802120 0.42731043 -0.8678404 -0.8075071 0.60924625 0.1867889
## sd       -0.07466188 0.50424029  0.1482723  0.1041929 0.47988756 0.8222156
## median    0.99054000 0.44359825 -0.8836339 -0.8265777 0.62590284 0.2226624
## Q25       0.97751975 0.34944645 -0.8602681 -0.8032396 0.57542156 0.1594589
## Q75       1.00000000 0.53914231 -0.8769634 -0.8272779 0.69983409 0.3041923
##                mode   centroid      peakf    meandom     mindom     maxdom
## duration -0.5801976 -0.5760356 -0.3442749 -0.5520184 -0.4395534 -0.4990648
## meanfreq  0.9876984  1.0000000  0.3746895  0.9716538  0.7623498  0.9346644
## sd       -0.1694371 -0.1876710 -0.3876864 -0.1818978 -0.2652301 -0.1128540
## median    0.9932076  0.9936121  0.3483108  0.9723942  0.7490638  0.9371465
## Q25       0.9874968  0.9921899  0.3781663  0.9700893  0.7595913  0.9296068
## Q75       0.9788687  0.9880212  0.3167142  0.9603656  0.7343151  0.9314592
##             dfrange   modindx   startdom     enddom    dfslope
## duration -0.1826537 0.1169091 -0.2745944 -0.5801700 -0.6637003
## meanfreq  0.4188517 0.1699135  0.6900515  0.9092519  0.5756101
## sd        0.1678878 0.1198796 -0.1650753 -0.1855753 -0.2832770
## median    0.4392773 0.1901571  0.6769922  0.9110402  0.5704379
## Q25       0.4148611 0.1716737  0.6830710  0.9107326  0.5754538
## Q75       0.4494820 0.2101372  0.6766135  0.8957976  0.5375374

Matrices de correlación