2019-11-10

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:

diastema <- read.csv("diastema.csv")

Datos de ejemplo

Juego de datos Diastema spp

str(diastema)
## 'data.frame':    901 obs. of  25 variables:
##  $ species    : Factor w/ 3 levels "diastema","hylaeformis",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ sound.files: Factor w/ 168 levels "10y12pm.wav",..: 17 17 17 17 17 17 17 17 18 18 ...
##  $ selec      : int  2 3 4 6 7 8 9 10 1 2 ...
##  $ 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
## corrplot 0.84 loaded
## 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.5-6
## Loading required package: MuMIn

Estadística descriptiva

Mínimo

x <- diastema$meanfreq

min(x)
## [1] 2.64822

Máximo

max(x)
## [1] 5.393333

Estadística descriptiva

Rango

range(x)
## [1] 2.648220 5.393333

Varianza

var(x)
## [1] 0.4473477

Estadística descriptiva

Desviación estándar

sd(x)
## [1] 0.6688406

Error estándar

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

Estadística descriptiva

Quartiles

quantile(x)
##       0%      25%      50%      75%     100% 
## 2.648220 3.405402 3.784998 4.298062 5.393333
quantile(x, probs = c(0.025, 0.975))
##     2.5%    97.5% 
## 2.724807 5.207910

Estadística descriptiva

Resumen

summary(x)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.648   3.405   3.785   3.860   4.298   5.393

Estadística descriptiva

Resumen

summary(diastema)
##         species                        sound.files      selec        
##  diastema   :563   D.vocator_AguaBuena_01.wav:  8   Min.   :   1.00  
##  hylaeformis:135   D.vocator_AguaBuena_03.wav:  8   1st Qu.:   4.00  
##  vocator    :203   D.vocator_AguaBuena_04.wav:  8   Median :   8.00  
##                    D.vocator_AguaBuena_07.wav:  8   Mean   :  20.66  
##                    D.vocator_GOLFITO_01.wav  :  8   3rd Qu.:  14.00  
##                    D.vocator_GOLFITO_02.wav  :  8   Max.   :1026.00  
##                    (Other)                   :853                    
##     duration          meanfreq           sd             median     
##  Min.   :0.02601   Min.   :2.648   Min.   :0.1182   Min.   :2.565  
##  1st Qu.:0.06048   1st Qu.:3.405   1st Qu.:0.2885   1st Qu.:3.255  
##  Median :0.10213   Median :3.785   Median :0.3502   Median :3.804  
##  Mean   :0.09298   Mean   :3.860   Mean   :0.3810   Mean   :3.831  
##  3rd Qu.:0.11810   3rd Qu.:4.298   3rd Qu.:0.4642   3rd Qu.:4.331  
##  Max.   :0.18032   Max.   :5.393   Max.   :0.9598   Max.   :5.477  
##                                                                    
##       Q25             Q75             IQR                skew        
##  Min.   :2.455   Min.   :2.607   Min.   :0.008685   Min.   : 0.9212  
##  1st Qu.:3.215   1st Qu.:3.496   1st Qu.:0.186969   1st Qu.: 2.6701  
##  Median :3.667   Median :4.076   Median :0.288732   Median : 5.1661  
##  Mean   :3.724   Mean   :4.013   Mean   :0.289331   Mean   : 5.7503  
##  3rd Qu.:4.213   3rd Qu.:4.488   3rd Qu.:0.381443   3rd Qu.: 8.6154  
##  Max.   :5.171   Max.   :5.646   Max.   :1.463277   Max.   :13.7573  
##                                                                      
##       kurt             sp.ent            sfm                mode      
##  Min.   :  2.425   Min.   :0.4050   Min.   :0.008168   Min.   :2.455  
##  1st Qu.: 10.662   1st Qu.:0.6815   1st Qu.:0.120140   1st Qu.:3.226  
##  Median : 36.287   Median :0.7558   Median :0.203346   Median :3.738  
##  Mean   : 52.533   Mean   :0.7345   Mean   :0.225889   Mean   :3.790  
##  3rd Qu.: 86.779   3rd Qu.:0.8142   3rd Qu.:0.317360   3rd Qu.:4.281  
##  Max.   :200.870   Max.   :0.9501   Max.   :0.760965   Max.   :5.667  
##                                                                       
##     centroid         peakf            meandom          mindom     
##  Min.   :2.648   Min.   :0.01406   Min.   :2.110   Min.   :1.895  
##  1st Qu.:3.405   1st Qu.:2.67695   1st Qu.:3.221   1st Qu.:2.929  
##  Median :3.785   Median :3.35073   Median :3.646   Median :3.488  
##  Mean   :3.860   Mean   :2.98061   Mean   :3.751   Mean   :3.441  
##  3rd Qu.:4.298   3rd Qu.:4.03547   3rd Qu.:4.264   3rd Qu.:4.005  
##  Max.   :5.393   Max.   :5.21253   Max.   :5.398   Max.   :4.823  
##                                                                   
##      maxdom         dfrange          modindx          startdom    
##  Min.   :2.584   Min.   :0.0000   Min.   :0.0000   Min.   :1.895  
##  1st Qu.:3.359   1st Qu.:0.1723   1st Qu.:0.2000   1st Qu.:3.144  
##  Median :3.962   Median :0.3445   Median :0.2857   Median :3.618  
##  Mean   :4.009   Mean   :0.5685   Mean   :0.3223   Mean   :3.614  
##  3rd Qu.:4.651   3rd Qu.:0.8613   3rd Qu.:0.3750   3rd Qu.:4.264  
##  Max.   :6.029   Max.   :2.7563   Max.   :1.0000   Max.   :5.685  
##                                                                   
##      enddom         dfslope       
##  Min.   :1.895   Min.   :-25.580  
##  1st Qu.:3.187   1st Qu.: -1.761  
##  Median :3.747   Median :  0.000  
##  Mean   :3.805   Mean   :  4.839  
##  3rd Qu.:4.307   3rd Qu.: 11.539  
##  Max.   :6.029   Max.   : 49.834  
## 

Estadística descriptiva

Frecuencias (cuentas)

table(diastema$species)
## 
##    diastema hylaeformis     vocator 
##         563         135         203

Revisión de supuestos

Normalidad

shapiro.test(x)
## 
##  Shapiro-Wilk normality test
## 
## data:  x
## W = 0.97673, p-value = 8.306e-11

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(diastema, aes(x = meanfreq)) + 
  
  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.98142, p-value = 2.695e-09

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

diast2 <- diastema[diastema$species %in% c("diastema", "vocator"), ]

var.test(diast2$duration ~ diast2$species)
## 
##  F test to compare two variances
## 
## data:  diast2$duration by diast2$species
## F = 7.6358, num df = 562, denom df = 202, p-value < 2.2e-16
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  6.042163 9.531010
## sample estimates:
## ratio of variances 
##           7.635828

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 = 185300310, 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)
##                selec   duration   meanfreq          sd     median
## selec     1.00000000 -0.1146174  0.1092018 -0.06607895  0.1083299
## duration -0.11461742  1.0000000 -0.5760356  0.56732529 -0.5592035
## meanfreq  0.10920176 -0.5760356  1.0000000 -0.18767105  0.9936121
## sd       -0.06607895  0.5673253 -0.1876710  1.00000000 -0.1491603
## median    0.10832987 -0.5592035  0.9936121 -0.14916026  1.0000000
## Q25       0.11202386 -0.5771093  0.9921899 -0.20929449  0.9909145
##                 Q25         Q75         IQR       skew       kurt
## selec     0.1120239  0.10137154 0.003000185 -0.1269901 -0.1051796
## duration -0.5771093 -0.49858583 0.089624414  0.5730427  0.4960304
## meanfreq  0.9921899  0.98802120 0.427310432 -0.8678404 -0.8075071
## sd       -0.2092945 -0.07466188 0.504240288  0.1482723  0.1041929
## median    0.9909145  0.99054000 0.443598249 -0.8836339 -0.8265777
## Q25       1.0000000  0.97751975 0.349446453 -0.8602681 -0.8032396
##              sp.ent         sfm       mode   centroid      peakf
## selec    0.05677291 -0.03036253  0.1166506  0.1092018 -0.1281589
## duration 0.08517835  0.37071177 -0.5801976 -0.5760356 -0.3442749
## meanfreq 0.60924625  0.18678887  0.9876984  1.0000000  0.3746895
## sd       0.47988756  0.82221565 -0.1694371 -0.1876710 -0.3876864
## median   0.62590284  0.22266244  0.9932076  0.9936121  0.3483108
## Q25      0.57542156  0.15945890  0.9874968  0.9921899  0.3781663
##             meandom      mindom     maxdom     dfrange     modindx
## selec     0.1225399  0.09395709  0.1332595  0.07830479 -0.03061407
## duration -0.5520184 -0.43955344 -0.4990648 -0.18265368  0.11690907
## meanfreq  0.9716538  0.76234976  0.9346644  0.41885165  0.16991349
## sd       -0.1818978 -0.26523014 -0.1128540  0.16788779  0.11987961
## median    0.9723942  0.74906384  0.9371465  0.43927729  0.19015714
## Q25       0.9700893  0.75959129  0.9296068  0.41486113  0.17167366
##            startdom     enddom       dfslope
## selec     0.1271745  0.1261979  0.0004864947
## duration -0.2745944 -0.5801700 -0.6637003273
## meanfreq  0.6900515  0.9092519  0.5756101288
## sd       -0.1650753 -0.1855753 -0.2832770413
## median    0.6769922  0.9110402  0.5704378915
## Q25       0.6830710  0.9107326  0.5754537628

Matrices de correlación

Matrices de correlación

Evaluar colinealidad de variables (paquete “corrplot”)

corrplot(corr = mat.cor, method = "ellipse", 
         col = mi.paleta(200), tl.col = "black")

Matrices de correlación

Matrices de correlación

Evaluar colinealidad de variables (paquete “corrplot”)

corrplot.mixed(corr = mat.cor[1:6, 1:6], 
               upper = "number", lower = "ellipse")

Matrices de correlación

Prueba de hipótesis

  • Regresión lineal
    • Modela la relación (lineal) entre una variable continua (respuesta) y una o mas variables predictoras
  • Supuestos
    • Homocedasticidad
    • Independencia de los errores (residuales)
    • Predictores no son colineares

Prueba de hipótesis

Regresión lineal

reg <- lm(diastema$duration ~ diastema$duration)
## Warning in model.matrix.default(mt, mf, contrasts): the response appeared
## on the right-hand side and was dropped
## Warning in model.matrix.default(mt, mf, contrasts): problem with term 1 in
## model.matrix: no columns are assigned
summary(reg)
## 
## Call:
## lm(formula = diastema$duration ~ diastema$duration)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.066970 -0.032502  0.009153  0.025117  0.087339 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.092979   0.001171    79.4   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03515 on 900 degrees of freedom

Prueba de hipótesis

Regresión lineal

coef(reg)
## (Intercept) 
##  0.09297865
confint(reg)
##                  2.5 %     97.5 %
## (Intercept) 0.09068052 0.09527679

Prueba de hipótesis

Regresión lineal

plot(diastema$meanfreq, diastema$duration * 1000, pch = 20, 
  col = adjustcolor(sample(cols, 1), alpha.f = 0.3), cex = 2, 
  ylab = "Duración (ms)", 
  xlab = "Frecuencia promedio (kHz)")
abline(reg)

Prueba de hipótesis

Regresión lineal

ggplot(diastema, aes(x = meanfreq, y = duration)) + 
  geom_point(fill = sample(cols, 1), size = 2, shape = 21) + 
  ylab("Duración (ms)") +  xlab("Frecuencia promedio (kHz)") +
  geom_smooth(method = "lm", col = "red2")

Prueba de hipótesis

  • ANOVA
    • Compara promedios entre grupos
  • Supuestos
    • Homocedasticidad
    • Normalidad de las variables
    • Normalidad de los residuales
    • Independencia de las observaciones
reg3 <- lm(diast2$duration ~ diast2$species)

Prueba de hipótesis

Regresión lineal (ANOVA)

anova(reg3)
## Analysis of Variance Table
## 
## Response: diast2$duration
##                 Df  Sum Sq Mean Sq F value    Pr(>F)    
## diast2$species   1 0.80098 0.80098  2568.7 < 2.2e-16 ***
## Residuals      764 0.23823 0.00031                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Prueba de hipótesis

ANOVA

reg3 <- lm(diast2$duration ~ diast2$species)
anova(reg3)
## Analysis of Variance Table
## 
## Response: diast2$duration
##                 Df  Sum Sq Mean Sq F value    Pr(>F)    
## diast2$species   1 0.80098 0.80098  2568.7 < 2.2e-16 ***
## Residuals      764 0.23823 0.00031                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Prueba de hipótesis

ANOVA: normalidad de residuales

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

## Prueba de hipótesis

Pruebas a-posteriori

reg3 <- lm(diastema$duration ~ diastema$species)
anova(reg3)
## Analysis of Variance Table
## 
## Response: diastema$duration
##                   Df  Sum Sq Mean Sq F value    Pr(>F)    
## diastema$species   2 0.80213 0.40107  1162.8 < 2.2e-16 ***
## Residuals        898 0.30974 0.00034                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Prueba de hipótesis

Pruebas a-posteriori

pairwise.t.test(diastema$duration, diastema$species, p.adj="bonferroni")
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  diastema$duration and diastema$species 
## 
##             diastema hylaeformis
## hylaeformis <2e-16   -          
## vocator     <2e-16   <2e-16     
## 
## P value adjustment method: bonferroni

Prueba de hipótesis

Pruebas a-posteriori

reg4 <- aov(duration ~ species, data = diastema)

TukeyHSD(reg4)

Prueba de hipótesis

Pruebas a-posteriori

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = duration ~ species, data = diastema)
## 
## $species
##                             diff         lwr         upr p adj
## hylaeformis-diastema -0.01625793 -0.02043615 -0.01207970     0
## vocator-diastema     -0.07326968 -0.07683909 -0.06970026     0
## vocator-hylaeformis  -0.05701175 -0.06185379 -0.05216971     0

Prueba de hipótesis

Pruebas a-posteriori (paquete ‘multcomp’)

reg4T <- glht(reg4, linfct = mcp(species = "Tukey"))

confint(reg4T)

Prueba de hipótesis

Pruebas a-posteriori (paquete ‘multcomp’)

## 
##   Simultaneous Confidence Intervals
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: aov(formula = duration ~ species, data = diastema)
## 
## Quantile = 2.3382
## 95% family-wise confidence level
##  
## 
## Linear Hypotheses:
##                             Estimate lwr      upr     
## hylaeformis - diastema == 0 -0.01626 -0.02042 -0.01210
## vocator - diastema == 0     -0.07327 -0.07682 -0.06971
## vocator - hylaeformis == 0  -0.05701 -0.06183 -0.05219

Prueba de hipótesis

Pruebas a-posteriori (paquete ‘multcomp’)

##                           Estimate         lwr         upr
## hylaeformis - diastema -0.01625793 -0.02041944 -0.01209642
## vocator - diastema     -0.07326968 -0.07682481 -0.06971454
## vocator - hylaeformis  -0.05701175 -0.06183442 -0.05218908
## attr(,"conf.level")
## [1] 0.95
## attr(,"calpha")
## [1] 2.33821

Prueba de hipótesis

Pruebas a-posteriori (paquete ‘multcomp’)

plot(reg4T)

Prueba de hipótesis

Pruebas a-posteriori (paquete ‘multcomp’)

dif.letras <- cld(reg4T)

dif.letras
##    diastema hylaeformis     vocator 
##         "c"         "b"         "a"

Prueba de hipótesis

Pruebas a-posteriori

boxplot(diastema$duration ~ diastema$species, 
        col = sample(cols, 1), ylim = c(0, 0.2))

text(x = 1:4, y = 0.197, labels = dif.letras$mcletters$Letters)

Prueba de hipótesis

Pruebas a-posteriori

ggplot(diastema, aes(x = species, y = duration)) + 
  geom_boxplot(fill = sample(cols, 1)) + 
  annotate(geom = "text", x = 1:3, y = 0.2, 
           label = dif.letras$mcletters$Letters) 

Prueba de hipótesis

Pruebas a-posteriori

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = diastema$duration ~ diastema$species)
## 
## $`diastema$species`
##                             diff         lwr         upr p adj
## hylaeformis-diastema -0.01625793 -0.02043615 -0.01207970     0
## vocator-diastema     -0.07326968 -0.07683909 -0.06970026     0
## vocator-hylaeformis  -0.05701175 -0.06185379 -0.05216971     0

Prueba de hipótesis

  • Kruskal-Wallis (ANOVA no paramétrico)
    • Compara promedios entre grupos
  • Supuestos
    • Independencia de las observaciones

Prueba de hipótesis

Kruskal-Wallis (ANOVA no paramétrico)

kw <- kruskal.test(diastema$duration ~ diastema$species)

kw
## 
##  Kruskal-Wallis rank sum test
## 
## data:  diastema$duration by diastema$species
## Kruskal-Wallis chi-squared = 490.59, df = 2, p-value < 2.2e-16

Prueba de hipótesis

  • Chi-cuadrado (tabla de contingencia)
    • Evalúa si el número de muestras para la combinación de niveles en 2 variables categóricas sigue una distribución determinada
  • Supuestos
    • Independencia de las observaciones
    • N para cada celda debe ser >= 5

Prueba de hipótesis

Chi-cuadrado

diast2$dur.cat <- ifelse(diast2$duration > 0.1, 1, 0)
set.seed(10)
diast2$dur.cat[sample(1:nrow(diast2), 100)] <- 
  sample(c(0, 1), size = 50, replace = TRUE)

species.dur <- table(diast2$dur.cat, diast2$species)
species.dur
##    
##     diastema hylaeformis vocator
##   0      187           0     191
##   1      376           0      12

Prueba de hipótesis

Chi-cuadrado

species.dur <- species.dur[,-2]

chisq.test(species.dur)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  species.dur
## X-squared = 218.76, df = 1, p-value < 2.2e-16

Prueba de hipótesis

Chi-cuadrado con simulación del valor de P con Monte Carlo

chisq.test(species.dur, simulate.p.value = TRUE, B = 1000)
## 
##  Pearson's Chi-squared test with simulated p-value (based on 1000
##  replicates)
## 
## data:  species.dur
## X-squared = 221.19, df = NA, p-value = 0.000999

Ejercicio 1

Que prueba alternativa existe para tablas de contingencia y como se corre en R?

Prueba de hipótesis

  • Prueba de Mantel
    • Estima la correlación entre matrices de distancia (similitud)
    • Utiliza aleatorización para calcular valores de P
  • Supuestos
    • Matrices creadas con métodos de distancia equivalentes

Prueba de hipótesis

Prueba de Mantel

data(varespec)

data(varechem)

veg.dist <- vegdist(varespec) # Bray-Curtis
env.dist <- vegdist(scale(varechem), "euclid")

Prueba de hipótesis

Prueba de Mantel (paquete ‘vegan’)

mantel(veg.dist, env.dist)
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = veg.dist, ydis = env.dist) 
## 
## Mantel statistic r: 0.3047 
##       Significance: 0.001 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.124 0.150 0.185 0.220 
## Permutation: free
## Number of permutations: 999

Prueba de hipótesis

Métodos de aleatorización (Monte Carlo)

  • Calculan un estadístico luego de aleatorizar los datos
  • Repiten el procedimiento muchas veces hasta generar una distribución nula
  • Comparan el valor observado con la distribución nula para calcular el valor de P
  • No necesariamente se generan todas las combinaciones posibles en la aleatorización

Prueba de hipótesis

Métodos de aleatorización (Monte Carlo)

Prueba T simple:

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

Métodos de aleatorización (Monte Carlo) T pareada por aleatorización:

reps <- 10000
dif.prom <- NULL
sleep1 <- sleep

for(i in  1:reps) {
  sleep1$group <- sample(sleep1$group, nrow(sleep1))
m1 <-  mean(sleep1$extra[sleep1$group == "1"])
m2 <-  mean(sleep1$extra[sleep1$group == "2"])
  dif.prom[i] <- m1 - m2
}

m1 <-  mean(sleep$extra[sleep$group == "1"])
m2 <-  mean(sleep$extra[sleep$group == "2"])

Prueba de hipótesis

Métodos de aleatorización (Monte Carlo)

T pareada por aleatorización:

obs <- m1 - m2
obs
## [1] -1.58
pval <- length(dif.prom[dif.prom <= obs])/reps
pval
## [1] 0.0443

Prueba de hipótesis

Métodos de aleatorización (Monte Carlo)

T pareada por aleatorización:

hist(dif.prom, col = sample(cols, 1), main = NULL, 
     breaks = 20, xlim = c(range(c(obs, dif.prom))))
abline(v = obs, col = "red", lty = 3, lwd = 4)

Prueba de hipótesis

Métodos de aleatorización (Monte Carlo)

T pareada por aleatorización:

sleep2 <- sleep
  
  sleep2$extra[sleep2$group == 1] <- rnorm(n = 10, mean = 10)
  sleep2$extra[sleep2$group == 2] <- rnorm(n = 10, mean = 12)

Prueba de hipótesis

Métodos de aleatorización (Monte Carlo)

T pareada por aleatorización:

reps <- 1000
dif.prom <- NULL
sleep1 <- sleep2

for(i in  1:reps) {
  sleep1$group <- sample(sleep1$group, nrow(sleep1))
m1 <-  mean(sleep1$extra[sleep1$group == "1"])
m2 <-  mean(sleep1$extra[sleep1$group == "2"])
  dif.prom[i] <- m1 - m2
  }

m1 <-  mean(sleep2$extra[sleep2$group == "1"])
m2 <-  mean(sleep2$extra[sleep2$group == "2"])

Prueba de hipótesis

Métodos de aleatorización (Monte Carlo)

T pareada por aleatorización:

obs <- m1 - m2

pval <- length(dif.prom[dif.prom <= obs])/reps


obs 
## [1] -2.127808
pval
## [1] 0

Prueba de hipótesis

Métodos de aleatorización (Monte Carlo)

T pareada por aleatorización:

hist(dif.prom, col = sample(cols, 1), main = NULL, 
     breaks = 20, xlim = c(range(c(obs, dif.prom))))
abline(v = obs, col = "red", lty = 3, lwd = 4)

Prueba de hipótesis

Métodos de aleatorización (Monte Carlo)

T pareada por aleatorización:

sleep2 <- sleep
  
  sleep2$extra[sleep2$group == 1] <- rnorm(n = 10, mean = 10)
  sleep2$extra[sleep2$group == 2] <- rnorm(n = 10, mean = 12)

Prueba de hipótesis

Métodos de aleatorización (Monte Carlo)

T pareada por aleatorización:

reps <- 10000
dif.prom <- NULL
sleep1 <- sleep2

for(i in  1:reps) {
  sleep1$group <- sample(sleep1$group, nrow(sleep1))
m1 <-  mean(sleep1$extra[sleep1$group == "1"])
m2 <-  mean(sleep1$extra[sleep1$group == "2"])
  dif.prom[i] <- m1 - m2
  }

m1 <-  mean(sleep2$extra[sleep2$group == "1"])
m2 <-  mean(sleep2$extra[sleep2$group == "2"])

Prueba de hipótesis

Métodos de aleatorización (Monte Carlo)

T pareada por aleatorización:

obs <- m1 - m2

pval <- length(dif.prom[dif.prom <= obs])/reps

obs 
## [1] -1.701389
pval
## [1] 0.0015
rereps <- 1000
pvals <- NULL

for(i in 1:1000)
{
  samp.dif <- sample(dif.prom, rereps)
pvals[i] <- length(samp.dif[samp.dif <= obs])/rereps 
  
}

quantile(pvals, c(0.025, 0.975))
##  2.5% 97.5% 
## 0.000 0.004

Prueba de hipótesis

Métodos de aleatorización (Monte Carlo)

T pareada por aleatorización:

hist(dif.prom, col = sample(cols, 1), main = NULL, 
     breaks = 20, xlim = c(range(c(obs, dif.prom))))
abline(v = obs, col = "red", lty = 3, lwd = 4)

Prueba de hipótesis

Regresión logística

log.mod <- glm(diast2$dur.cat ~ diast2$duration, family="binomial")

summary(log.mod)

Prueba de hipótesis

Regresión logística

## 
## Call:
## glm(formula = diast2$dur.cat ~ diast2$meanfreq, family = "binomial")
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.0974  -0.8342   0.4682   0.8167   2.6699  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      10.4499     0.7985   13.09   <2e-16 ***
## diast2$meanfreq  -2.6050     0.2003  -13.00   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1061.77  on 765  degrees of freedom
## Residual deviance:  799.85  on 764  degrees of freedom
## AIC: 803.85
## 
## Number of Fisher Scoring iterations: 4

Prueba de hipótesis

Regresión logística

plot(diast2$meanfreq, diast2$dur.cat, xlab="Frecuencia",ylab="Duracion (prob)",
     col = sample(cols, 1), pch = 20, cex = 2)

curve(plogis(log.mod$coefficients[1] + 
               (log.mod$coefficients[2] *x)), 
      add=TRUE, col = "black", lwd = 3)

Prueba de hipótesis

Regresión logística

Prueba de hipótesis

Regresión logística

ggplot(diast2, aes(x=meanfreq, y= dur.cat)) + 
  geom_point(size = 3, col = sample(cols, 1)) + 
  stat_smooth(method="glm", 
    method.args=list(family="binomial"), col = "red2") + 
  theme_classic() + labs(x = "Frecuencia", 
                         y = "Duracion (probabilidad)")

Prueba de hipótesis

Regresión logística

Prueba de hipótesis

Regresión logística por aleatorización

reps <- 1000
dif.prom <- NULL
datos <- diast2

coef.sim <- sapply(1:reps, function(x)
{
  datos$dur.cat <- sample(datos$dur.cat, nrow(datos)) 
  log.mod <- glm(datos$dur.cat ~ datos$meanfreq, family="binomial")
  cf <- as.numeric(coef(log.mod)[2])
  return(cf)
  }
  )

Prueba de hipótesis

Regresión logística por aleatorización

log.mod <- glm(diast2$dur.cat ~ diast2$meanfreq, family="binomial")

obs <- coef(log.mod)[2]

pval <- length(coef.sim[coef.sim <= obs])/reps
pval

hist(coef.sim, col = sample(cols, 1), main = NULL, 
     breaks = 20, xlim = c(range(c(obs, coef.sim))))
abline(v = obs, col = "red", lty = 3, lwd = 4)

Prueba de hipótesis

Regresión logística por aleatorización

## [1] "p = 1"

Prueba de hipótesis

Regresión logística por aleatorización

diast2$meanfreq <- rnorm(nrow(diast2))

log.mod <- glm(diast2$dur.cat ~ diast2$meanfreq, family="binomial")

obs <- coef(log.mod)[2]

pval <- length(coef.sim[coef.sim <= obs])/reps
pval

hist(coef.sim, col = sample(cols, 1), main = NULL, 
     breaks = 20, xlim = c(range(c(obs, coef.sim))))
abline(v = obs, col = "red", lty = 3, lwd = 4)

Prueba de hipótesis

Regresión logística por aleatorización

## [1] "p = 0.847"

Ejercicio 2

Calcule un valor de P para la prueba de chi-cuadrado utilizando aleatorización

species.dur <- table(diast2$dur.cat, diast2$species)

chisq.test(species.dur)
diast4 <- diast2
diast4 <- droplevels(diast4)
reps <- 10000
chi.ale <- NULL 


for(i in 1:reps) {
  diast4$dur.cat <- sample(diast4$dur.cat)
 lc <- table(diast4$dur.cat, diast4$species)
  chi.ale[i] <- chisq.test(lc)$statistic
  }

obs <- chisq.test(species.dur)$statistic

pvals <- NULL
rereps <- 100

for(i in 1:1000)
{
  samp.chi <- sample(chi.ale, rereps)
pvals[i] <- length(samp.chi[samp.chi >= obs])/rereps 
  
}

quantile(pvals, c(0.025, 0.975))
##  2.5% 97.5% 
##     0     0
pval <- length(chi.ale[chi.ale >= obs])/reps 

pval
## [1] 0

Ejercicio 3

Calcule un valor de P para la prueba de correlación

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

cor(x = x2, y = y2)

cor.test(x = x2, y = y2)
reps <- 1000
c.reps <- NULL
for(i in 1:reps){
  x3 <- sample(x2)
  y3 <- sample(y2)
  creps <- cor.test(x3,y3)$estimate
  c.reps[i] <- creps
}

c.reps
##    [1] -1.015953e-02 -2.261267e-02 -1.111459e-02  2.880492e-02
##    [5]  3.436904e-02  3.768925e-02  4.171066e-02  2.803483e-02
##    [9]  4.392359e-02  2.800918e-02  2.501211e-02 -4.718466e-03
##   [13] -3.811913e-02 -1.127598e-02  7.171447e-03  2.011460e-02
##   [17] -5.479223e-02 -8.742714e-03 -2.157783e-02  2.367969e-02
##   [21] -2.061514e-02  3.806788e-02  8.069575e-02 -8.538483e-04
##   [25] -2.138330e-02  3.429724e-02 -8.445972e-03  1.327269e-03
##   [29] -8.150053e-04  1.228991e-02  7.850104e-02 -7.810280e-03
##   [33] -2.470054e-02  2.943678e-02  1.394471e-03 -1.805195e-02
##   [37] -1.724253e-02  2.632799e-03 -2.017357e-02 -5.174630e-04
##   [41] -1.476360e-02  9.081523e-03  1.458509e-02  2.143630e-02
##   [45]  7.150365e-02 -2.071055e-02 -8.198404e-03  2.471508e-02
##   [49]  2.838905e-02 -1.257428e-02  1.090595e-03 -3.221435e-02
##   [53] -8.542908e-03  2.971168e-02  2.744995e-02 -1.501305e-02
##   [57]  3.282222e-02 -4.018017e-02  3.065225e-03  2.743530e-02
##   [61]  3.497663e-02 -1.841111e-02 -6.323639e-03 -4.057748e-02
##   [65] -1.556743e-02 -1.308125e-02  6.349919e-03  2.542335e-02
##   [69]  2.777603e-02  3.902090e-02  1.079499e-02 -3.085947e-02
##   [73] -9.274109e-03  1.067688e-02 -8.528155e-02 -2.214772e-02
##   [77]  1.573147e-02  6.119947e-03 -1.631183e-02  2.593794e-02
##   [81] -5.478298e-02  2.624223e-02 -5.877859e-03  7.987307e-03
##   [85] -7.549016e-03 -1.359460e-02  6.405611e-03  3.824298e-02
##   [89] -1.487210e-02 -1.371926e-02  3.957992e-02  2.526918e-02
##   [93]  2.972728e-02  6.398930e-02 -8.663822e-03  3.932800e-02
##   [97] -3.101221e-02 -3.574420e-02  1.391167e-02 -5.006281e-03
##  [101]  3.450836e-02  5.339876e-03  5.108617e-02 -3.898307e-02
##  [105] -6.275977e-03  1.025326e-02  3.551235e-02  2.388270e-02
##  [109]  1.775279e-02 -5.022460e-02 -5.284387e-02  6.145855e-03
##  [113]  4.667264e-02  1.979461e-02  1.212743e-02 -3.291580e-02
##  [117] -6.160451e-02 -2.840035e-02 -3.808629e-03  5.070669e-03
##  [121]  1.276014e-02 -3.755358e-03  2.512825e-02 -4.704245e-02
##  [125] -1.201259e-02 -1.822888e-02 -4.180626e-03  6.530994e-02
##  [129] -2.569255e-02  3.599867e-02 -1.352518e-02  1.519200e-02
##  [133]  3.118046e-02  5.011823e-04 -2.896494e-02  1.535140e-02
##  [137]  1.020970e-02 -3.339649e-02 -5.805950e-02 -1.781430e-02
##  [141] -2.829633e-02 -2.053614e-03  1.904476e-02 -1.142590e-02
##  [145]  5.405626e-02 -6.085077e-03  3.524271e-02 -1.026068e-02
##  [149]  1.424594e-02  2.587803e-02 -6.771109e-03 -6.027765e-02
##  [153] -2.763801e-02 -5.934621e-02  8.174809e-03  5.658655e-02
##  [157]  5.480690e-02  4.707039e-02 -1.641950e-02 -4.105597e-02
##  [161]  8.639446e-02 -5.498639e-03  2.079214e-02 -2.089352e-02
##  [165] -2.341750e-02  5.647165e-02  2.895156e-02 -1.081575e-02
##  [169] -2.591591e-02 -3.541572e-02  1.681437e-02  4.380359e-02
##  [173]  3.491005e-02 -2.577211e-02 -7.603820e-02 -1.737532e-02
##  [177]  2.360779e-02  1.308455e-03 -9.932196e-03  4.040915e-03
##  [181]  8.415134e-02 -3.477658e-02 -9.677272e-03  2.733735e-03
##  [185] -5.323560e-02 -1.854715e-02 -1.155057e-02  1.152276e-01
##  [189]  1.896657e-03  4.787789e-03 -1.367061e-02  5.997051e-02
##  [193]  3.228635e-02  1.056015e-02  1.115334e-02  6.372680e-03
##  [197] -1.981266e-02  1.910001e-02  1.907692e-02 -4.406889e-02
##  [201] -1.540997e-02 -1.624485e-02  6.444029e-02  4.061388e-02
##  [205] -1.989690e-02 -8.501343e-02  1.394737e-02  2.280909e-02
##  [209] -3.859639e-02 -2.049566e-02  3.426250e-02 -4.501110e-02
##  [213]  8.140540e-02  8.605410e-03  3.149573e-02  3.356641e-02
##  [217]  3.528298e-02  2.916576e-03 -7.335771e-03 -3.647200e-02
##  [221]  3.744639e-02 -6.675822e-02  5.691454e-02  7.829772e-02
##  [225]  7.043087e-03 -1.391884e-02 -4.218040e-03  1.095678e-02
##  [229] -2.020839e-02 -4.950051e-02 -9.220357e-03  7.827905e-03
##  [233] -4.939828e-02  7.023112e-02  3.740012e-02 -1.266584e-02
##  [237]  2.143050e-02  4.853656e-02 -8.275441e-02  1.670382e-02
##  [241]  2.318216e-03 -9.593372e-03  2.016290e-03  3.122174e-03
##  [245]  4.537027e-02 -2.064375e-02  1.656531e-02  2.640164e-02
##  [249] -3.752374e-02  4.360290e-02  1.742491e-02  3.831721e-02
##  [253] -9.331136e-03 -2.330500e-02  8.082145e-03  8.753069e-03
##  [257]  8.605796e-02  1.018877e-02  2.733933e-02 -2.090856e-02
##  [261] -2.119001e-02  5.656339e-03 -1.314337e-02 -3.163789e-02
##  [265]  1.478108e-02  1.433034e-02 -4.193429e-03  2.444234e-02
##  [269]  4.337905e-02  1.144175e-02  2.430110e-02 -4.290647e-02
##  [273] -1.471899e-02  3.032169e-02 -2.214183e-03 -3.933583e-02
##  [277]  5.360219e-02  2.320916e-02  2.830184e-02 -8.652810e-03
##  [281]  2.943642e-02  3.889973e-02 -6.874890e-02 -1.622866e-02
##  [285]  5.187652e-02  3.635395e-02 -2.733985e-02 -6.138145e-03
##  [289] -1.958376e-02  2.619990e-02 -1.168696e-02 -5.558403e-03
##  [293]  3.516986e-02  4.161468e-02  6.723417e-02 -3.201782e-02
##  [297] -1.771501e-02 -3.992149e-02 -3.719268e-04  4.020031e-02
##  [301] -4.410596e-02 -3.055990e-02  2.118375e-03  1.610494e-02
##  [305]  1.565593e-02  2.640250e-02 -7.301965e-02 -2.650433e-03
##  [309]  5.094741e-02  3.364617e-02 -5.127683e-02  7.571211e-03
##  [313]  7.613924e-03 -3.863249e-03 -7.360515e-02  1.655258e-02
##  [317] -5.178176e-02 -4.182803e-02  1.258020e-02 -1.553909e-02
##  [321] -2.052912e-02  5.139980e-02 -1.616628e-02  4.790788e-03
##  [325] -1.806102e-02  9.916662e-03  2.699214e-02 -4.791170e-02
##  [329]  2.182731e-02  4.207977e-02  1.179925e-02  9.387650e-03
##  [333]  7.242128e-02 -1.512524e-02  7.019338e-03 -1.544852e-02
##  [337]  1.210781e-02  2.672451e-02 -3.205714e-02  3.850534e-02
##  [341] -1.607597e-02  8.534058e-02  1.315561e-02 -4.154977e-02
##  [345]  5.397665e-02 -1.045633e-02  2.219873e-02  5.620619e-03
##  [349]  6.394042e-02  4.767246e-04 -5.469729e-03 -4.163881e-02
##  [353] -1.994256e-02 -1.659646e-02  4.190347e-04  6.478057e-02
##  [357] -2.759843e-03  9.006659e-04 -1.320092e-02  4.765952e-02
##  [361] -2.535416e-03  3.720745e-03 -2.353053e-02  2.736643e-02
##  [365] -6.680803e-03 -8.523085e-03  4.586664e-02  4.078219e-03
##  [369] -2.963291e-02  5.245026e-02  5.274194e-02  5.378518e-02
##  [373] -6.318087e-02  4.908770e-02  1.316947e-02  1.269329e-03
##  [377] -1.450019e-02 -1.316618e-02  1.726789e-02 -9.471549e-03
##  [381] -3.569752e-02  3.021875e-02  6.776896e-03 -4.188548e-02
##  [385] -4.991912e-02  9.817693e-03  4.271031e-02 -5.690585e-02
##  [389]  5.507842e-03  2.416101e-02  1.984642e-02 -5.008449e-02
##  [393]  3.897383e-02  1.420635e-02 -1.852663e-02  6.033576e-02
##  [397]  3.678052e-02 -3.534581e-02  2.372654e-02 -5.009173e-03
##  [401]  5.704259e-03  3.439027e-03 -2.542819e-02  5.575728e-02
##  [405] -1.581904e-02 -2.499582e-02 -7.646558e-02 -2.938401e-02
##  [409] -6.254203e-03  2.658171e-02 -6.234265e-02 -3.829221e-02
##  [413]  6.513127e-03 -4.436000e-03 -4.730673e-02  4.596974e-02
##  [417]  3.444068e-02  4.813955e-02  1.658634e-02  8.786996e-03
##  [421]  2.173132e-02  4.405766e-02 -3.427448e-02 -3.621272e-02
##  [425] -2.699434e-02  1.650346e-02  2.478289e-02  3.985393e-02
##  [429] -1.802739e-02  3.605464e-02 -8.759139e-04 -1.137716e-02
##  [433] -4.560475e-03  1.681028e-03 -4.225100e-02 -1.948377e-02
##  [437]  3.173514e-02  1.028904e-02 -3.329038e-02  7.073988e-02
##  [441]  2.963628e-02 -6.137745e-02  1.187119e-02 -2.852749e-02
##  [445]  3.532933e-02 -3.764954e-02 -3.129952e-03  4.164864e-02
##  [449]  2.892018e-03 -2.969373e-02 -3.334556e-03 -1.438928e-03
##  [453] -5.987062e-03  7.387946e-02 -7.087791e-03  9.369771e-02
##  [457] -3.113293e-02 -2.633287e-02 -1.408613e-02 -4.766069e-02
##  [461] -5.573465e-03  1.120391e-03  2.234373e-02 -1.802187e-02
##  [465]  4.808691e-02  2.663974e-02 -2.416121e-02 -7.624378e-03
##  [469] -1.821649e-02 -5.249595e-02 -5.145299e-03  4.592822e-02
##  [473] -1.437253e-02  5.362163e-04 -1.246484e-03  1.267363e-02
##  [477]  2.029175e-02  1.396485e-03 -3.020636e-02 -2.087199e-02
##  [481]  2.943661e-02 -4.121701e-02  6.332692e-02 -4.185553e-02
##  [485] -9.847802e-03  6.344084e-03  1.739975e-02 -2.928550e-02
##  [489] -7.680745e-04  1.832996e-02  1.697156e-02  6.023965e-02
##  [493]  1.172690e-02  7.226164e-02  5.023776e-02  1.365324e-02
##  [497] -1.593915e-02 -3.943216e-02 -2.789986e-02 -5.335730e-03
##  [501] -2.902580e-02 -1.699539e-02 -1.526130e-02  2.414853e-02
##  [505] -4.920750e-02  2.935894e-02 -4.927568e-02  1.543185e-03
##  [509]  1.339865e-02 -1.882344e-02  2.654630e-02  2.260118e-02
##  [513] -1.573939e-03 -5.446232e-04 -2.928079e-02 -3.891943e-02
##  [517] -3.399635e-03  1.616649e-02 -2.457493e-03 -1.216016e-02
##  [521] -5.524330e-02 -1.961082e-02  4.722162e-02 -4.100664e-02
##  [525]  4.245821e-02 -3.589436e-02 -1.666057e-02  1.064037e-03
##  [529]  1.445741e-02 -1.179406e-02  4.911573e-02 -4.498152e-02
##  [533] -1.513266e-02  5.161893e-02  9.906717e-04 -3.080124e-02
##  [537] -3.078022e-02 -4.984510e-03 -2.793813e-02 -3.243323e-02
##  [541] -1.406635e-03  5.577392e-02  5.410491e-02  2.136039e-02
##  [545] -3.495766e-02 -4.192116e-04  2.556798e-03  1.425041e-02
##  [549]  2.334347e-02  3.912987e-02  6.162056e-02 -4.483896e-02
##  [553]  1.211629e-02  4.597238e-02  5.052588e-02 -4.138523e-03
##  [557] -3.243196e-02  1.758643e-02  4.234962e-02  1.298736e-02
##  [561] -2.486783e-03  3.098891e-04  5.640070e-02  5.770337e-02
##  [565] -4.852519e-03  5.691275e-02 -1.086546e-02  2.025411e-02
##  [569]  6.106300e-02  4.787227e-02 -5.675293e-03  8.471620e-02
##  [573]  1.988951e-02  1.830076e-02 -4.948987e-03 -7.610167e-03
##  [577]  7.679663e-02 -1.982148e-03  5.203113e-02  1.361164e-02
##  [581]  1.275663e-02 -2.020389e-02 -3.798591e-03  3.522183e-02
##  [585] -3.131389e-03  3.359982e-02 -4.210215e-02 -2.371798e-02
##  [589]  2.169612e-02 -1.112838e-02  5.442029e-02  1.005370e-02
##  [593]  1.610279e-02 -1.565699e-02  1.937763e-02 -3.925667e-03
##  [597] -1.083714e-04  1.249336e-02 -7.935656e-04 -4.845839e-02
##  [601] -3.233368e-03 -4.969029e-03 -4.782686e-03  7.348828e-02
##  [605] -4.769390e-02 -1.819939e-02  2.517381e-03  3.353162e-02
##  [609]  1.814400e-04 -4.120134e-02 -1.508279e-02 -6.329781e-02
##  [613] -2.416283e-02  5.077750e-03 -2.015331e-02  4.618135e-03
##  [617]  1.137071e-02 -3.346705e-02  8.118157e-03 -6.613164e-03
##  [621]  2.397905e-02  1.399499e-02 -5.589548e-03  2.362799e-02
##  [625] -2.387353e-02 -1.490649e-03 -1.952017e-02 -4.971954e-02
##  [629] -2.592255e-02 -2.054143e-02  1.013066e-04 -1.156568e-02
##  [633]  3.777362e-02 -8.972774e-03  2.810852e-02  2.995264e-02
##  [637]  3.263252e-02  7.184229e-02  4.174954e-02 -2.387313e-02
##  [641] -2.275688e-02 -2.897735e-04  5.951520e-02 -4.675142e-03
##  [645] -2.410984e-03 -1.949353e-02  1.527847e-02 -1.923451e-02
##  [649]  5.326102e-02  1.002748e-03  5.971534e-03 -4.678326e-02
##  [653] -3.846530e-02  1.967918e-02 -2.471415e-02 -1.825343e-02
##  [657]  1.032650e-02  3.945189e-03 -2.728292e-02 -3.892063e-02
##  [661]  1.217109e-02 -2.164449e-02 -2.335803e-02  2.090088e-02
##  [665]  1.482631e-02  3.997444e-02 -3.906199e-02 -3.003848e-02
##  [669] -3.368607e-02  4.595900e-02  4.155091e-02 -4.244513e-02
##  [673]  2.126371e-03  2.755743e-02  1.911537e-02  1.080055e-02
##  [677]  6.422529e-02 -1.662921e-02  1.012902e-02  3.848191e-02
##  [681] -8.064491e-05  3.693178e-02 -2.161839e-02 -2.860945e-02
##  [685] -2.254896e-02 -5.605723e-02 -2.071207e-02  1.967903e-02
##  [689]  4.523779e-03  1.837598e-02 -6.658514e-03 -2.096881e-02
##  [693] -9.614334e-03 -5.938547e-02 -2.372919e-02  3.760562e-02
##  [697] -4.917922e-02 -7.235948e-03 -1.083055e-02 -8.612905e-03
##  [701]  1.643139e-02 -6.434857e-02  3.476493e-02 -1.804783e-02
##  [705] -4.830361e-03 -1.315496e-02  6.910657e-03 -6.221000e-03
##  [709] -2.285467e-02 -7.595207e-04 -1.416169e-02  3.848935e-02
##  [713] -6.138591e-02 -5.435910e-02 -3.267669e-02 -2.086433e-02
##  [717]  5.328232e-02 -9.405048e-03 -3.225993e-02 -5.222192e-02
##  [721] -2.672973e-02  5.398832e-02 -3.840715e-02  7.023252e-03
##  [725] -3.697042e-02  5.795064e-02 -1.452002e-02 -1.447936e-02
##  [729] -7.230405e-03 -1.894983e-02  2.583143e-02 -3.650368e-02
##  [733] -4.166297e-03  3.516035e-02  5.133896e-02 -3.167801e-02
##  [737] -3.046989e-02  1.700286e-02 -3.790456e-02 -4.251962e-02
##  [741] -4.181170e-04  4.681447e-06 -3.179922e-02  2.559439e-02
##  [745] -2.559521e-02 -1.293442e-02  3.583909e-02  8.068202e-03
##  [749] -1.333969e-02  9.017700e-04 -5.603520e-02 -3.159744e-03
##  [753]  2.354644e-02  2.041108e-02 -2.258406e-02 -4.233968e-02
##  [757] -1.395410e-02  6.669981e-02 -1.797025e-02 -3.790076e-02
##  [761]  2.176142e-02  4.484144e-02  1.360978e-02  2.667411e-02
##  [765] -4.182889e-03  8.680848e-03  5.326164e-02 -7.012843e-03
##  [769] -5.381627e-03  2.520212e-02  8.729780e-02 -7.185118e-03
##  [773]  8.088558e-02 -4.630044e-02 -1.960607e-03 -2.195906e-02
##  [777] -1.779805e-02  3.371455e-02 -1.626543e-02  1.897931e-02
##  [781] -2.635172e-03  2.744734e-02  1.832008e-02  1.120207e-02
##  [785]  3.953930e-02  2.612655e-02 -9.940347e-03  4.344631e-02
##  [789] -1.404491e-01 -2.030419e-02  4.090368e-02 -9.206473e-02
##  [793]  3.892443e-02  5.911531e-02  1.072916e-02  2.081321e-02
##  [797] -3.869494e-02  1.743958e-02 -3.619372e-02  9.497982e-03
##  [801]  9.527769e-03 -3.418978e-02 -2.369592e-02  2.074769e-02
##  [805] -8.143880e-02  4.392439e-03 -6.445949e-03 -3.311062e-02
##  [809] -2.044082e-02  2.361157e-02 -2.327714e-02  1.406167e-02
##  [813] -8.519149e-04 -4.272997e-02  4.155070e-02  1.046826e-03
##  [817]  1.968644e-02 -9.962085e-03  3.761997e-02  2.793736e-02
##  [821]  1.335212e-02  2.574487e-02 -1.169535e-02 -2.685682e-02
##  [825]  2.211602e-02  1.614968e-02 -2.197058e-02 -1.110665e-02
##  [829]  2.889168e-02  2.576079e-03 -9.932759e-02  1.075344e-01
##  [833]  5.578109e-02  2.136735e-02  1.041943e-02  2.077717e-02
##  [837] -7.877688e-03  3.615430e-02 -2.614016e-02  2.786043e-02
##  [841] -7.406263e-03  3.242590e-02  7.781003e-03 -4.450991e-02
##  [845]  1.249749e-02  4.724716e-02 -5.142153e-02 -2.627645e-02
##  [849]  3.676705e-02 -1.384697e-02 -5.631522e-02  3.153971e-02
##  [853] -5.637565e-02  3.498632e-03  4.717273e-02  4.807012e-02
##  [857]  3.317858e-02  7.057534e-02 -7.338660e-04  4.935598e-02
##  [861]  7.989569e-02 -5.744489e-02  1.079751e-02  2.131994e-02
##  [865] -6.020267e-03  1.888126e-02  3.133678e-02 -2.824448e-02
##  [869]  2.088946e-02 -3.281464e-03 -5.176662e-02  3.171977e-03
##  [873]  1.666406e-02 -2.640448e-02 -1.515748e-02  2.801146e-02
##  [877] -7.482238e-03  4.168101e-03  2.005338e-02 -4.833638e-02
##  [881] -2.351412e-02  8.223680e-04 -2.831295e-02 -1.625549e-02
##  [885] -2.779990e-02 -4.284739e-02 -2.389612e-02  4.759841e-02
##  [889]  4.979815e-02  3.602854e-02 -3.812607e-02  8.829662e-03
##  [893] -1.544905e-02 -1.909993e-02  2.225670e-02  5.438077e-02
##  [897]  2.580379e-02  8.626487e-03  1.357861e-02  1.813987e-02
##  [901] -7.225976e-03 -5.941783e-02 -4.794708e-03  2.975585e-02
##  [905] -4.811812e-03 -4.198804e-02 -4.875349e-03  7.079470e-02
##  [909] -2.533706e-03  2.904029e-02  1.841588e-02 -2.823517e-02
##  [913]  3.202347e-02 -1.073192e-02 -2.322092e-03 -3.194054e-02
##  [917]  1.129558e-02 -3.951336e-02 -3.826654e-02  4.869067e-03
##  [921] -1.905794e-02  2.445711e-02 -1.395665e-03  3.073580e-03
##  [925] -3.625166e-03  3.495589e-02 -2.621830e-02  4.522680e-02
##  [929]  1.888062e-02  2.898510e-02  9.385867e-03 -1.255304e-02
##  [933] -4.101060e-03  7.270235e-03  2.261018e-02 -4.664092e-02
##  [937] -1.403038e-02  1.605414e-02 -4.234241e-02 -2.013423e-02
##  [941]  4.769423e-02  7.734680e-03  2.870226e-02  4.760386e-02
##  [945]  2.303481e-03 -6.315486e-02  2.245288e-02 -5.332358e-02
##  [949] -4.068768e-02 -3.491548e-02 -1.512635e-03 -2.198016e-02
##  [953]  1.475229e-02 -5.697055e-02  2.369785e-03 -2.188450e-02
##  [957] -7.214936e-02 -2.744246e-02 -1.398417e-02 -1.090718e-02
##  [961]  5.440762e-02 -3.094797e-03  4.784968e-03 -2.698914e-02
##  [965]  1.245283e-02  5.673139e-03  1.061822e-02  4.898375e-03
##  [969] -4.928850e-02 -1.231308e-02  1.230108e-02 -3.417847e-02
##  [973] -5.629633e-02  5.170395e-02 -3.532076e-02 -1.308372e-02
##  [977]  4.667539e-03  4.430163e-03  1.596862e-02  1.398784e-02
##  [981] -2.942175e-02 -2.727172e-02 -1.684678e-02 -6.744343e-02
##  [985]  6.282984e-03  2.250439e-02 -3.860157e-02 -5.874641e-02
##  [989]  1.134730e-01 -4.762697e-04  5.598025e-02 -9.298620e-03
##  [993]  5.338028e-02  1.302108e-03  1.444238e-03  4.408899e-02
##  [997] -3.744561e-02  3.796613e-02  3.435719e-03  2.420929e-02
cor.obs <- cor.test(x2, y2)$estimate

pval.cor <- length(c.reps[c.reps<=cor.obs])/reps
pval.cor
## [1] 0