# ROGER GUEVARA (roger.guevara@inecol.edu.mx)
# INSTITUTO DE ECOLOGIA A.C RED DE BIOLOGÍA EVOLUTIVA
# ESTADISTICA COMPUTACIONAL EN R
source("~/desktop/curso R 2012/funciones.r")
# PRUEBAS DE HIPÓTESIS
# Datos de la velocidad de la velocidad luz estimados restando 299000 a
# cada valor. La idea es comparar si en promedio estos valores difieren
# estadísticamente de 1000.
VEL_LUZ <- c(850, 740, 900, 1070, 930, 850, 950, 980, 980, 880, 1000, 980, 930,
650, 760, 810, 1000, 1000, 960, 960)
# Para esto usaremos una prueba de t de una sola muestra en donde el
# parámetro de referencia es 1000. Esto se reduce a la diferencia del
# promedio de los datos y el parámetro de referencia dividido entre el
# error estándar de la diferencia de la media y el parámetro de
# referencia. Dado que el parámetro de referencia es fijo, no hay
# varianza, y por lo tanto el error estándar de la diferencia de la media
# de los datos y el parámetro de referencia es el mismo error estándar de
# los datos. El cociente de la diferencia entre su error estándar, se
# distribuye como t
MEDIAS <- mean(VEL_LUZ)
ERROR <- ee(VEL_LUZ)
T <- (MEDIAS - 1000)/ERROR
T
## [1] -3.879
# y con la con la función probabilística de t << pt() >> calculamos la
# probabilidad de observar este valor o uno mayor, y esto lo haremos a dos
# colas de la distribución. Este ejercicio tiene n-1 grados de libertad ya
# que solo estimamos un parámetro, la media de los datos.
pt(T, 19) * 2
## [1] 0.001011
# Con la función << t.test() >> que toma dos argumentos en el caso de la t
# de una muestra obtenemos el mismo resultado. Los argumentos son 1) el
# vector de datos, en este caos << VEL_LUZ >> y << mu= >> que es el
# parámetro de referencia.
t.test(VEL_LUZ, mu = 1000)
##
## One Sample t-test
##
## data: VEL_LUZ
## t = -3.879, df = 19, p-value = 0.001011
## alternative hypothesis: true mean is not equal to 1000
## 95 percent confidence interval:
## 859.9 958.1
## sample estimates:
## mean of x
## 909
# Prueba de t para dos muestras
DATOS <- read.table("~/desktop/curso R 2012/toloache.csv", header = TRUE, sep = ",")
attach(DATOS)
names(DATOS)
## [1] "NOTAS" "NUM" "INVERNADERO" "LUZ"
## [5] "MICORRIZA" "SUELO" "HERBIVORIA" "RAIZ"
## [9] "TALLO" "HOJAS" "FLORES" "CAPSULAS"
## [13] "SEMILLAS" "INMADUROS" "PESO_TOTAL" "NUM_FLORES"
## [17] "DIAS_A_LA_FLOR" "FRUT_INMADUROS" "FRUT_MADUROS" "TOT_FRUTOS"
## [21] "ABORTOS" "ALT_TOT" "HOJAS_TOT" "DUREZA"
## [25] "AREA" "PESO_HOJAS" "AFE" "LONG_HIFAS"
## [29] "HIFAS" "ESPORAS" "VESICULAS" "ARBUSCULOS"
## [33] "TOTAL"
# Haremos dos subconjuntos de los valores de << DUREZA >> como una función
# de la calidad de luz << LUZ >>
DUREZA_ML <- DUREZA[DATOS$LUZ == "ML"]
DUREZA_PL <- DUREZA[DATOS$LUZ == "PL"]
# Calculamos las medias, las varianzas y el tamaño de las muestras
M1 <- mean(DUREZA_ML)
M2 <- mean(DUREZA_PL)
VAR1 <- var(DUREZA_ML)
VAR2 <- var(DUREZA_PL)
N1 <- length(DUREZA_ML)
N2 <- length(DUREZA_PL)
# Con las varianzas y los tamaños de muestra calculamos el error estándar
# de la diferencia de las medias.
EEDIF <- sqrt(VAR1/N1 + VAR2/N2)
# Y la diferencia entre las medias dividido pro el error estándar de la
# diferencia de las medias nos da el valor de t, que tendrá n - 2 grados
# de libertad ya que se estimaron dos parámetros apartir de los datos.
T <- (M1 - M2)/EEDIF
pt(T, N1 + N2 - 2) * 2
## [1] 0.007088
# La función de prueba de de t << t.test() >> para comparar dos muestras
# toma tres argumentos, que son cada uno de las dos muestras y <<
# var.equal= >> definido como verdadero << TRUE >> fuerza a la función de
# t a asumir que hay homocedasticidad de varianzas, lo mismo que asumimos
# en el los cálculos que hicimos arriba
t.test(DUREZA_ML, DUREZA_PL, var.equal = TRUE)
##
## Two Sample t-test
##
## data: DUREZA_ML and DUREZA_PL
## t = -2.732, df = 216, p-value = 0.006822
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -7.484 -1.211
## sample estimates:
## mean of x mean of y
## 51.87 56.22
# De no especificar que las varianzas son iguales, la función << t.test()
# >> automáticamente lleva a cabo la llamada corrección de Welch que
# pondera (disminuye) los grados de libertad lo que causa que para
# alcanzar un resultado significativo con P = 0.05 el valor de t debe ser
# mayor que sin no usamos la corrección, en otras palabras, la diferencia
# entre las medias debe ser mayo para que rechacemos la hipótesis nula
t.test(DUREZA_ML, DUREZA_PL, var.equal = FALSE)
##
## Welch Two Sample t-test
##
## data: DUREZA_ML and DUREZA_PL
## t = -2.719, df = 207.4, p-value = 0.00711
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -7.500 -1.195
## sample estimates:
## mean of x mean of y
## 51.87 56.22
# Prueba de t-pareada
DATOS <- read.table("~/desktop/curso R 2012/paired.txt", sep = "\t", header = TRUE)
attach(DATOS)
# La prueba de t pareada, es sencillamente una prueba de t de una muestra
# en donde esa muestra son las diferencias entre las parejas de datos. La
# función << t.test() >> en este caso usa un tercer argumento << paired =
# TRUE >>
t.test(Downstream, Upstream, paired = TRUE)
##
## Paired t-test
##
## data: Downstream and Upstream
## t = 3.077, df = 8, p-value = 0.01519
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.02589 0.18077
## sample estimates:
## mean of the differences
## 0.1033
# y ahora probemos haciendo una prueba de una sola muestra con las
# diferencias entre las parejas de datos
DIFERENCIAS <- Downstream - Upstream
t.test(DIFERENCIAS, mu = 0)
##
## One Sample t-test
##
## data: DIFERENCIAS
## t = 3.077, df = 8, p-value = 0.01519
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.02589 0.18077
## sample estimates:
## mean of x
## 0.1033
# Noten como estos dos últimos resultados son idénticos
# Prueba de Wilcoxon-------------------- Para comparar dos medias se
# asignas rangos a los datos de las dos muestras de manera conjunta, y se
# suman los valores en rangos de cada muestra, y a la suma mas pequeña se
# les resta la suma de rangos esperada si las dos muestras son
# completamente distinta, en este caso si una muestra tuviese los primeros
# nueve valores jerárquicos. Este valor esperado de la suma de rangos para
# una muestra se obtiene con (n*(n+1))/2 siendo n el tamaño de muestra de
# con la suma mas pequeña.
# A continuación generamos una distribución probabilística para comparar
# dos muestras cada una con nueve datos como en el ejemplo que estamos
# viendo
SUMA <- numeric()
n <- 9
for (i in 1:10000) {
AZAR <- rnorm(18, 10, 3)
RANGOS <- rank(AZAR)
M1 <- RANGOS[1:9]
M2 <- RANGOS[10:18]
S1 <- sum(M1) - (n * (n + 1))/2
S2 <- sum(M2) - (n * (n + 1))/2
if (S1 < S2)
SUMA[i] <- S1 else SUMA[i] <- S2
}
hist(SUMA)
# Ahora obtenemos la suma de valores de rangos para cada una de las
# muestras
sum(rank(c(Downstream, Upstream))[1:9]) - (9 * 10)/2
## [1] 53.5
sum(rank(c(Downstream, Upstream))[10:18]) - (9 * 10)/2
## [1] 27.5
# y comparamos que probabilidad tendríamos de observar este valor o uno
# mas pequeño por azar contando todos los casos de la distribución
# empírica que obtuvimos aleatoriamente. El valor que comparamos es el mas
# pequeño de las dos muestras, 27.5
sum(SUMA <= 27.5)/i
## [1] 0.2501
# Ahora comparamos nuestro resultado con lo que obtenemos con la función
# << wilcox.test() >> en donde los argumentos son los objetos que definen
# cada una de las dos muestras.
wilcox.test(Downstream, Upstream)
## Warning: cannot compute exact p-value with ties
##
## Wilcoxon rank sum test with continuity correction
##
## data: Downstream and Upstream
## W = 53.5, p-value = 0.2694
## alternative hypothesis: true location shift is not equal to 0
# Ahora probemos utilizar la prueba paramétrica de t, pero usando los
# valores en rangos de as muestras
R_DS <- rank(c(Downstream, Upstream))[1:9]
R_US <- rank(c(Downstream, Upstream))[10:18]
t.test(R_DS, R_US)
##
## Welch Two Sample t-test
##
## data: R_DS and R_US
## t = 1.16, df = 15.99, p-value = 0.263
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -2.390 8.168
## sample estimates:
## mean of x mean of y
## 10.944 8.056
# Note la similitud de resultados entre estos procedimientos. Igualmente
# los procedimientos arrojan resultados muy similares si el análisis es
# pareado. Recuerden que esto significa que hay una sola muestra y la suma
# de valores en rangos con signo se espera que sea cero.
wilcox.test(Downstream, Upstream, paired = TRUE)
##
## Wilcoxon signed rank test
##
## data: Downstream and Upstream
## V = 44, p-value = 0.007812
## alternative hypothesis: true location shift is not equal to 0
t.test(R_DS, R_US, paired = TRUE)
##
## Paired t-test
##
## data: R_DS and R_US
## t = 3.956, df = 8, p-value = 0.004198
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 1.205 4.573
## sample estimates:
## mean of the differences
## 2.889
# Prueba de varianza---------------------- El cociente de dos varianza
# define la distribución F. Así, dadas dos muestras tomadas al azar de la
# misma población se esperaría que su cociente fuese 1, y mediante la
# distribución F podemos calcular la probabilidad de observar un valor
# igual o mayor al que resulte de dividir la varianza mayor por la
# varianza menor
V1 <- var(Downstream)
V2 <- var(Upstream)
CO <- V2/V1
# Con la distribución probabilística de F, usando la función << pf() >>
# colocamos tres argumentos. 1) el valor del cociente, 2) los grados de
# libertad del numerador, y 3) los grados de libertad del denominador. Y
# calculamos la probabilidad hacía la cola derecha de la distribución ya
# que en esta distribución la cola derecha se extiende hasta infinito
# mientras que la otra cola tiende a cero.
pf(CO, 8, 8, lower = FALSE)
## [1] 0.3265
# La probabilidad obtenida la debemos multiplicar por dos ya que en esta
# prueba no hay un criterio a priori para esperar que la varianza de una
# muestra sea mayor que la de la otra muestra.
# Y con la función << var.test() >> obtenemos este resultado directamente
var.test(Upstream, Downstream)
##
## F test to compare two variances
##
## data: Upstream and Downstream
## F = 1.389, num df = 8, denom df = 8, p-value = 0.6529
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.3134 6.1588
## sample estimates:
## ratio of variances
## 1.389
# Tablas de contingencia----------------------- Para poner a prueba si la
# frecuencia con que ocurren dos eventos son independientes se utilizan
# las tablas de contingencia. Estas tablas resumen los conteos realizados
# en donde cada conteo se puede asignar a la combinación de dos factores y
# cada factor tiene al menos dos niveles. En este ejemplo usaremos tres
# tipos de epífitas creciendo en tres tipos de forofitos.
# Primero generamos la matriz de observaciones
M <- matrix(c(50, 8, 10, 25, 20, 5, 2, 1, 12), 3, 3)
# Y para facilitar la identidad de las filas y columnas le agregamos
# nombres a las filas y columnas
rownames(M) <- c("Helechos", "Orqid.", "Bromelias")
colnames(M) <- c("Palmas", "Árboles", "Postes")
M
## Palmas Árboles Postes
## Helechos 50 25 2
## Orqid. 8 20 1
## Bromelias 10 5 12
# El principio de las tablas de contingencia es multiplicar las
# probabilidades de observar cada uno de los eventos por separado ya que
# esta es la hipótesis nula, de que los eventos ocurren
# independientemente. Entonces primero debemos obtener estas
# probabilidades que son en esencia la suma de frecuencias de cada fila y
# columna dividido por el total de la suma de frecuencias. Para esto
# usamos las funciones << colSums() >> y << rowSums() >>
colSums(M)
## Palmas Árboles Postes
## 68 50 15
colSums(M)
## Palmas Árboles Postes
## 68 50 15
# y las probabilidades son
PROB_FOROFITOS <- colSums(M)/sum(M)
PROB_EPIFITAS <- rowSums(M)/sum(M)
# Así podemos observar que la probabilidad de que una epífita cualquiera
# este en una palma es de 0.5112782, y la probabilidad de que una epífita
# sea un helecho es de 0.5789474
# De los productos de estas probabilidades marginales obtenemos la
# probabilidad de que dos eventos ocurren independientemente. Por ejemplo
# multiplicamos la probabilidad de que una epífita este en una palma por
# la probabilidad de que una epífita sea un helecho, 0.5112782 x 0.5789474
# de lo que obtenemos 0.2960032 que es la probabilidad de con que
# esperamos observar, solo por azar, el que un helecho este en una palma.
# Así al obtener las probabilidades de las diferentes combinaciones solo
# nos restaría multiplicarlas por la frecuencia total para obtener la
# frecuencia esperada.
PRO_ESP <- rep(PROB_FOROFITOS, each = 3) * rep(PROB_EPIFITAS, 3)
# y Podemos colocar las probabilidades en una matriz
PRO_ESP <- matrix(PRO_ESP, 3, 3)
rownames(PRO_ESP) <- c("Helechos", "Orqid.", "Bromelias")
colnames(PRO_ESP) <- c("Palmas", "Árboles", "Postes")
# Y la frecuencia esperada es multiplicar esta tabla por la frecuencia
# total de epífitas
FREC_ESPE <- PRO_ESP * sum(M)
# Y con la prueba de bondad de ajuste podemos calcular un valor de ji
# cuadrada usando la fórmula (OBSERVADO - ESPERADO)^2/ESPERADO
(M - FREC_ESPE)^2/FREC_ESPE
## Palmas Árboles Postes
## Helechos 2.871 0.5383 5.145
## Orqid. 3.143 7.5919 1.576
## Bromelias 1.049 2.6133 26.334
# Y el valor de ji-cuadrada es la suma de todas las celdas
sum((M - FREC_ESPE)^2/FREC_ESPE)
## [1] 50.86
# Finalmente con la distribución probabilistca de ji-cuadrada obtenemos la
# probabilidad esperada. Aqui debemos usar 1 menos la probabilidad
# calculada. Al igual que la dirtibución F, en la distribución de
# ji-cuadrada solo la cola derecha tiende a infinito por lo que se calcula
# la probabilidad sobre este lado de la distribución. Note que después de
# calcular cuatro frecuencias esperadas del interior de la tabla el resto
# ya no se estiman a partir de los datos ya que su valores quedan
# predefinidos para que las filas y columnas sumen lo mismo que en la
# tabla original.
pchisq(50.86188, 4, lower = FALSE)
## [1] 2.386e-10
# La función << chis.test() >> calcula esta prueba rápidamente y usa como
# argumento el objeto que contiene la matriz de frecuencias esperadas.
chisq.test(M)
## Warning: Chi-squared approximation may be incorrect
##
## Pearson's Chi-squared test
##
## data: M
## X-squared = 50.86, df = 4, p-value = 2.386e-10
# Si asignamos este análisis a un objeto podemos explorar las frecuencias
# esperadas, y los residuales estandarizados. Estos últimos nos servirán
# para explorar cuales celdas de la tabla contribuyen significativamente
# de lo esperado por azar.
CONTINGENCIA <- chisq.test(M)
## Warning: Chi-squared approximation may be incorrect
CONTINGENCIA$expected
## Palmas Árboles Postes
## Helechos 39.37 28.95 8.684
## Orqid. 14.83 10.90 3.271
## Bromelias 13.80 10.15 3.045
CONTINGENCIA$residuals
## Palmas Árboles Postes
## Helechos 1.694 -0.7337 -2.268
## Orqid. -1.773 2.7553 -1.256
## Bromelias -1.024 -1.6166 5.132
# De la tabla anterior las celdas con valores absolutos mayores a 1.96 son
# aquellas en las que el desvío de lo esperado es realmente significativo.
# Si son positivos los residuales indican que las frecuencias observadas
# son mayores de lo que se espera por azar. Por otra parte si los
# residuales son negativos indican que las frecuencias observadas son
# menores de lo que se espera por azar.
# La función << mosaicplot() >> que usa como argumento el objeto con las
# frecuencias esperadas nos da una representación gráfica de los
# residuales esperados, y usando el argumento << shade=TRUE >> el gráfico
# presenta colores rojos y azul en las celdas sub y sobre estimadas
# respectivamente. De cada color hay dos tonos. El mas claro de ellos
# indican que el desvío de las frecuencias es mayor 2 pero menor a 4
# desviaciones estándares,lo que significa que tienen una probabilidad
# entre 0.05 y 0.00006 y el tono mas oscuro indica aquellas celdas con
# probabilidad menor a 0.00006. Estas probabilidades se desprende de la
# distribución normal.
mosaicplot(t(M), shade = TRUE)
