# ROGER GUEVARA (roger.guevara@inecol.edu.mx)
# INSTITUTO DE ECOLOGIA A.C RED DE BIOLOGÍA EVOLUTIVA
# ESTADISTICA COMPUTACIONAL EN R
# Continuamos con el mismo ejemplo de las siete calificaciones que vimos
# la sesión anterior
CAL <- c(10, 8, 7, 8, 10, 8, 7)
PROM <- mean(CAL)
DE <- sd(CAL)
# Lo que haremos ahora e repetir el mismo procedimiento de la
# verosimilitud que usamos para estimar la media cuando se maximiza la
# función de verosimilitud. La única diferencia es que cambiaremos
# arbitrariamente el valor de la desviación estándar. El código siguiente
# es el mismo que y habíamos elaborado, salvo que definidos el objeto <<
# DE2 >> que contiene diez valores entre 0.1 y 5 que se usarán un turno
# para sustituir le valor de la desviación estándar (1.253566) al momento
# del cálculo de la verosimilitud. También cambiaremos << DE >> dentro de
# la función << dnomr() >> a << DE2 >> para que el cálculo de la
# verosimilitud se haga con el valor de desviación estándar
# correspondiente, y anexamos un segundo bucle para hacer esto todo
# automáticamente.
for (DE2 in seq(0.1, 5, length = 10)) {
SEC <- seq(6, 10, length = 10000)
i <- 0
VERO <- numeric()
for (PROM2 in SEC) {
i <- i + 1
VERO[i] <- sum(dnorm(CAL, PROM2, DE2, log = T))
}
MAXIMO <- which(VERO == max(VERO))
print(c(DE2, SEC[MAXIMO]))
}
## [1] 0.100 8.286
## [1] 0.6444 8.2858
## [1] 1.189 8.286
## [1] 1.733 8.286
## [1] 2.278 8.286
## [1] 2.822 8.286
## [1] 3.367 8.286
## [1] 3.911 8.286
## [1] 4.456 8.286
## [1] 5.000 8.286
# Aquí se puede observar que el valor estimado de la media (columna 2) por
# el método de máxima verosimilitud en los diez casos es siempre el mismo
# independientemente del valor de la desviación estándar que se haya
# utilizado (columna 1).
# La Desviación estándar -----------------------------------
# La definición de la desviación estándar es un promedio de la distancia
# de cada observación a la media.
# Gráficamente representamos las siete calificaciones como sigue con la
# función << plot() >>, en la que el primer argumento, la secuencia de 1 a
# 7, son las coordenadas en el eje X. El segundo argumento << CAL >> son
# las coordenadas en el eje Y. Utilizamos el argumento << pch= >> para
# modificar el tipo de símbolo que se despliega en el gráfico, y los
# argumentos << xlab= >> y << ylab= >> para dar nombres a los ejes. El
# texto para nombrar los ejes debe estar entre comillas para que R no lo
# considere como objeto.
plot(1:7, CAL, pch = 19, xlab = "Número de dato", ylab = "Calificaciones")
abline(h = PROM, col = "red")
# Podemos calcular la diferencia entre las calificaciones y el promedio
DIFERENCIA <- CAL - PROM
# y con la función << lines() >> trazamos unas lineas verticales de cada
# dato a la media, que representan las diferencias que existen entre los
# datos y la media
for (i in 1:7) {
lines(c(i, i), c(CAL[i], PROM), lty = 2)
}
# y si estas diferencias las sumamos inevitablemente nos va a dar como
# resultado cero. Aquí no da cero exactamente pero es un valor muy pequeño
# a la menos decimoquinta potencia.
sum(DIFERENCIA)
## [1] -5.329e-15
# En consecuencia el promedio de las diferencias de los datos respecto a
# la media será también cero.
# Podríamos usar los valores absolutos de las diferencias para tener un
# estimado de la desviación promedio de los datos respecto a la media.
mean(abs(DIFERENCIA))
## [1] 0.9796
# Esta solución es poco elegante, por lo que en la historia se favoreció
# el método de elevar las diferencias al cuadrado lo cual elimina los
# signos negativos y nos permite tener un estimado de las diferencias de
# los datos respecto a la media pero elevadas al cuadrado. Sumando estas
# diferencias al cuadrado tenemos lo que se conoce como la SUMA DE
# CUADRADOS o LA SUMA DE LAS DESVIACIONES AL CUADRADO
sum(DIFERENCIA^2)
## [1] 9.429
# Y obteniendo el promedio, al dividir por siete (el número de datos)
# obtenemos lo que se conoce como la varianza
sum(DIFERENCIA^2)/7
## [1] 1.347
# Y finalmente sacando la raíz cuadrada de este promedio tenemos un
# estimado de la distancia promedio de los datos a la media en la misma
# escala en que están los datos, y se conoce como la desviación estándar
sqrt(mean(DIFERENCIA^2))
## [1] 1.161
# El estimado de la varianza y la desviación estándar de una muestra está
# subestimado mediante este procedimiento por lo que el cálculo del
# promedio de las desviaciones se hace sobre los grados de libertad, y no
# sobre el tamaño de la muestra. En este caso los grados de libertad son n
# - 1
# ERROR ESTANDAR Para esto usaremos ahora la base de datos << TOLOACHE.csv
# >>
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"
# y para la variable << DUREZA >> calcularemos su media, desviación
# estándar y error estándar. Al no existir una función que calcule el
# error estándar en R nosotros podemos generarla para ya contar con ella
# en cualquier momento. Para esto usamos la función << function() >> que
# tiene como argumentos una o mas variables que se deberá introducir para
# el cálculo. La sintaxis de esta función ademas incluye que entre las
# llaves << { } >> se colocan todas las instrucciones que sean necesarias
# para calcular lo que deseamos. En este caso sacar la raíz cuadrada de la
# varianza de una variable dividida por el número de datos
ee <- function(X) {
sqrt(var(X)/length(X))
}
# Ahora si calculamos los valores
MEDIA <- mean(DUREZA)
DE <- sd(DUREZA)
EE <- ee(DUREZA)
# Reflexionemos un poco sobre que es el error estándar. Ya sabemos como
# estimar el error estándar pero, ¿qué representa?. El error estándar de
# una muestra es un estimado de la desviación estándar que tendría una
# colección infinita de medias calculada cada una de una muestra aleatoria
# de la misma población y con el mismo número de observaciones que la
# muestra original. Es decir es como si hiciéramos muestreos repetidos la
# misma población. En la práctica no se puede hacer, pero aquí en el
# espacio virtual lo haremos fácilmente. Definimos primero un objeto
# llamado << MEDIAS >> en el que almacenaremos los estimados de la media
# de cada una de las muestras que generemos, y utilizamos un bucle para
# tomar diez mil muestras. Dentro del bucle tenemos la función << rnorm()
# >> que genera números aleatorios para una distribución normal con media,
# en este caso igual a la media de << DUREZA >>, y desviación estándar
# igual a la desviación estándar de << DUREZA >>. El primer argumento <<
# length(DUREZA) >> indica cuantos valores aleatorios queremos que se
# generen, en este caso 218 valores. Con esto dentro del bucle se obtiene
# una muestra aleatoria de la población definida y a esta muestra se le
# calcula la media y se almacena en el objeto << MEDIAS >>
MEDIAS <- numeric()
for (i in 1:10000) {
MUESTRA_A <- rnorm(length(DUREZA), MEDIA, DE)
MEDIAS[i] <- mean(MUESTRA_A)
}
# Ya concluido el bucle obtenemos la desviación estándar de las diez mil
# medias y este resultado lo comparamos con el cálculo del error estándar
# de los 218 datos de dureza y podemos observar que son muy similares, no
# idénticos por que el concepto de error estándar se base en una colección
# de infinitas medias. Si aumentamos el número de repeticiones hechas, por
# ejemplo a 50 mil observaremos que hay mayor coincidencia entre el error
# estándar estimado de los datos y la desviación estándar de la colección
# de medias.
sd(MEDIAS)
## [1] 0.8094
ee(DUREZA)
## [1] 0.8069
# Intervalo de confianza -----------------------------------
# Fundamentado en el teorema del límite central y la distribución
# probabilística normal sabemos que en el ámbito comprendido entre -1.96
# desviaciones estándares y +1.96 desviaciones estándares se incluye el
# 95% del área bajo la curva.
# Con el siguiente código generamos la información necesaria para pintar
# una distribución normal señalando la porción central que incluye 95% del
# área que bajo la curva. No se detalla como se hace esto por que no es
# relevante para los propósitos del curso
CUANTILES <- seq(-4, 4, length = 1000)
NORM <- dnorm(seq(-4, 4, length = 1000), 0, 1)
plot(CUANTILES, NORM, type = "l")
lines(CUANTILES[CUANTILES >= -1.96 & CUANTILES <= 1.96], NORM[CUANTILES >= -1.96 &
CUANTILES <= 1.96], col = "red", lwd = 3)
lines(c(-1.96, -1.96), c(0, dnorm(-1.96)), col = "red", lwd = 3, lty = 2)
lines(c(1.96, 1.96), c(0, dnorm(1.96)), col = "red", lwd = 3, lty = 2)
text(0, 0.15, "95%", col = "red")
# En consecuencia cuando nosotros tenemos un estimado de la desviación
# estándar de un conjunto infinito de medias, es decir el error estándar
# de una muestra que proviene de una población con media y desviación
# estándar distinto a cero y uno respectivamente debemos multiplicar el
# valor de error estándar por 1.96 y este valor restado y sumado de la
# media es el ámbito bajo la distribución normal de la población que
# incluye el 95% del área bajo la curva.
# Uno puede ingenuamente pensar que esto se traduce en que con 0.95 de
# probabilidad el valor real de la media de la población se encuentra
# dentro de este intervalo. Sin embargo, lo que realmente representa este
# intervalo calculado a partir de la muestra de datos es que 95% de los
# intervalos que se calcularan de muestras aleatorias de la misma
# población y con el mismo número de datos que nuestra muestra original
# contendrán al valor de la media real, y en consecuencia un 5% de dichos
# intervalos NO contendrá a la media de la población.
# Gráficamente podemos ilustrarlo de la manera siguiente. Con la función
# << rnorm() vamos a generar muestras aleatorias que contengan el mismo
# número de datos que tenemos en la variable << DUREZA >>, de una
# población que con la media y desviación estándar calculada de los datos.
# Para cada muestra extraeremos el error estándar, con la función que
# escribimos arriba, y también calculamos la media de cada muestra. En
# total vamos a generar 100 muestras y así tendremos una colección de 100
# errores estándares y 100 medias.
EES <- numeric()
MEDIAS <- numeric()
i <- 0
for (i in 1:100) {
MUESTRA_A <- rnorm(length(DUREZA), MEDIA, DE)
EES[i] <- ee(MUESTRA_A)
MEDIAS[i] <- mean(MUESTRA_A)
}
# Ya que tenemos estas colecciones de errores estándar y medias pintaremos
# en un gráfico los intervalos de confianza al 95% usando la función <<
# lines() >>. La lógica de esta función es que para pintar una linea
# cualquiera se requiere de definir dos punto en un plano. Así la función
# << lines() >> tiene dos argumentos fundamentales: 1) las coordenadas de
# los puntos en el eje X, y las coordenadas de los puntos en el eje Y. Las
# coordenadas en X para un intervalo de confianza serán los puntos
# definidos por la media ± 1.96*Error estándar. Y las coordenadas en Y
# para cada intervalo serán iguales, en este ejemplo es el objeto << i >>
# que tomará en secuencia los valores de uno a 100. Finalmente trazaremos
# una linea vertical en el valor de la media de la población.
plot(c(50, 58), c(0, 100), type = "n")
for (i in 1:100) {
lines(c(MEDIAS[i] + 1.96 * EES[i], MEDIAS[i] - 1.96 * EES[i]), c(i, i))
}
abline(v = MEDIA, col = "red")
# Como es evidente de los 100 intervalos de confianza calculados hay
# algunos de ellos que no incluyen el valor de la media de la población.
# En la realidad nunca tenemos certeza de si un intervalo de confianza que
# calculamos de nuestros datos incluya o no a la media de la población.
# Solo sabemos que con 0.95 de probabilidad ese intervalo si incluye a la
# media de la población.
# ELABORAR FUNCIÓN DESCRIPTIVA---------------------- Ahora elaboraremos
# una función que nos permita extraer la estadística descriptiva básica de
# una variable numérica lo cual nos facilitará hacer la exploración de los
# datos antes de realizar pruebas estadísticas. La función calcula la
# media, varianza, desviación estándar, error estándar, 1ro, 2do (mediana)
# y 3er cuartiles, mínimo y máximo ya sea de una única variable numérica o
# de una matriz que contenga dos o mas columnas de variables numéricas. El
# nombre de la función es << resumen >>
resumen <- function(X) {
if (is.data.frame(X)) {
MAT <- matrix(NA, 9, ncol((X)))
for (i in 1:ncol(X)) {
X2 <- X[, i]
MAT[, i] <- Y <- c(mean(X2, na.rm = T), var(X2, na.rm = T), sd(X2,
na.rm = T), sqrt(var(X2, na.rm = T)/length(X2)), quantile(X2,
c(0.25, 0.5, 0.75), na.rm = T), min(X2, na.rm = T), max(X2,
na.rm = T))
}
colnames(MAT) <- colnames(X)
rownames(MAT) <- c("Media", "Varianza", "Desv. Est.", "EE", "1er Q",
"Mediana", "3er Q", "Mínimo", "Máximo")
print(MAT)
}
if (is.vector(X)) {
Y <- c(mean(X), var(X, na.rm = T), sd(X, na.rm = T), sqrt(var(X, na.rm = T)/length(X)),
quantile(X, c(0.25, 0.5, 0.75), na.rm = T), min(X, na.rm = T), max(X,
na.rm = T))
names(Y) <- c("Media", "Varianza", "Desv. Est.", "Error Est", "1er Q",
"Mediana", "3er Q", "Mínimo", "Máximo")
print(Y)
}
}
# Como ejemplo extraemos la estadística descriptiva de la variable <<
# DUREZA >> y de un conjunto de variables numéricas que están en las
# columnas 9 a 13 de la base de datos << DATOS >>
resumen(DUREZA)
## Media Varianza Desv. Est. Error Est 1er Q Mediana
## 53.9617 141.9409 11.9139 0.8069 42.7550 48.4000
## 3er Q Mínimo Máximo
## 62.7700 21.4800 90.9100
resumen(DATOS[, 9:13])
## TALLO HOJAS FLORES CAPSULAS SEMILLAS
## Media 9.1171 4.7296 0.30167 2.5997 5.5383
## Varianza 77.0681 18.5819 0.07930 7.7554 36.9673
## Desv. Est. 8.7788 4.3107 0.28161 2.7849 6.0801
## EE 0.5946 0.2920 0.01907 0.1886 0.4118
## 1er Q 0.8100 0.5175 0.06000 0.0000 0.0000
## Mediana 3.9450 3.4850 0.20000 0.8100 0.9100
## 3er Q 16.9600 8.7275 0.48000 5.1200 11.6600
## Mínimo 0.0400 0.0200 0.01000 0.0000 0.0000
## Máximo 27.8200 12.9900 1.32000 9.0700 20.4400