- R básico tiene funciones para las pruebas mas comunes de estadística:
- correlación
- regresión
- chi-cuadrado …
2019-11-10
Juego de datos:
diastema <- read.csv("diastema.csv")
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 ...
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
Mínimo
x <- diastema$meanfreq min(x)
## [1] 2.64822
Máximo
max(x)
## [1] 5.393333
Rango
range(x)
## [1] 2.648220 5.393333
Varianza
var(x)
## [1] 0.4473477
Desviación estándar
sd(x)
## [1] 0.6688406
Error estándar
sd(x)/sqrt(length(x))
## [1] 0.02228231
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
Resumen
summary(x)
## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 2.648 3.405 3.785 3.860 4.298 5.393
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 ##
Frecuencias (cuentas)
table(diastema$species)
## ## diastema hylaeformis vocator ## 563 135 203
Normalidad
shapiro.test(x)
## ## Shapiro-Wilk normality test ## ## data: x ## W = 0.97673, p-value = 8.306e-11
Normalidad
mi.paleta <- colorRampPalette(c("#E37222", "#66B9BF")) cols <- adjustcolor(mi.paleta(10), alpha.f = 0.4) hist(x, breaks = 8, col = sample(cols, 1))
Normalidad
Normalidad (ggplot)
library(ggplot2) ggplot(diastema, aes(x = meanfreq)) + geom_histogram(color = sample(cols, 1), fill = sample(cols, 1)) + theme_classic()
Normalidad
qqnorm(x, col = sample(cols, 1), pch = 20, cex= 2) qqline(x, col=2)
Normalidad (ggplot)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
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
Normalidad y transformación de variables
hist(log(x), breaks = 15, col = sample(cols, 1))
Normalidad y transformación de variables
qqnorm(log(x), col = sample(cols, 1), pch = 20, cex= 2) qqline(log(x), col=2)
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 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 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
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
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 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
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
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
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
Evaluar colinealidad de variables (paquete “corrplot”)
corrplot(corr = mat.cor)
Evaluar colinealidad de variables (paquete “corrplot”)
corrplot(corr = mat.cor, method = "ellipse", col = mi.paleta(200), tl.col = "black")
Evaluar colinealidad de variables (paquete “corrplot”)
corrplot(corr = mat.cor[1:6, 1:6], method = "number")
Evaluar colinealidad de variables (paquete “corrplot”)
corrplot.mixed(corr = mat.cor[1:6, 1:6], upper = "number", lower = "ellipse")
Evaluar colinealidad de variables
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
Regresión lineal
coef(reg)
## (Intercept) ## 0.09297865
confint(reg)
## 2.5 % 97.5 % ## (Intercept) 0.09068052 0.09527679
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)
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")
reg3 <- lm(diast2$duration ~ diast2$species)
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
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
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
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
Pruebas a-posteriori
reg4 <- aov(duration ~ species, data = diastema) TukeyHSD(reg4)
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
Pruebas a-posteriori (paquete ‘multcomp’)
reg4T <- glht(reg4, linfct = mcp(species = "Tukey")) confint(reg4T)
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
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
Pruebas a-posteriori (paquete ‘multcomp’)
plot(reg4T)
Pruebas a-posteriori (paquete ‘multcomp’)
dif.letras <- cld(reg4T) dif.letras
## diastema hylaeformis vocator ## "c" "b" "a"
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)
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)
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
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
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
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
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
Que prueba alternativa existe para tablas de contingencia y como se corre en R?
Prueba de Mantel
data(varespec) data(varechem) veg.dist <- vegdist(varespec) # Bray-Curtis env.dist <- vegdist(scale(varechem), "euclid")
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
Métodos de aleatorización (Monte Carlo)
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
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"])
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
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)
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)
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"])
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
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)
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)
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"])
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
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)
Regresión logística
log.mod <- glm(diast2$dur.cat ~ diast2$duration, family="binomial") summary(log.mod)
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
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)
Regresión logística
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)")
Regresión logística
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) } )
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)
Regresión logística por aleatorización
## [1] "p = 1"
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)
Regresión logística por aleatorización
## [1] "p = 0.847"
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
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