# ROGER GUEVARA (roger.guevara@inecol.edu.mx)
# INSTITUTO DE ECOLOGIA A.C RED DE BIOLOGÍA EVOLUTIVA
# ESTADISTICA COMPUTACIONAL EN R
# Correlación-------
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"
# La correlación es una mediad de cambio entre dos variable numéricas. Y
# esta mediad se obtiene de la covarianza entre las dos variables
# ponderado por el producto de las desviaciones estándares.
# La covarianza se calcula en la misma forma que la varianza salvo que es
# el la sumatoria del producto de las diferencias de cada pareja de datos
# respecto a su media.
# Aquí calculamos las medias. Usamos el argumento << na.rm= TRUE>> para
# omitir los valores ausentes del calculo. De no usar este argumento R no
# procesará la función
M_F <- mean(FLORES, na.rm = TRUE)
M_NF <- mean(NUM_FLORES, na.rm = TRUE)
# De la misma forma calculamos las desviaciones estándares de cada
# variable.
S_F <- sd(FLORES, na.rm = TRUE)
S_NF <- sd(NUM_FLORES, na.rm = TRUE)
# Y con esto obtenemos los productos de las diferencias de las medias
# entre las dos variables.
SPXY <- (FLORES - M_F) * (NUM_FLORES - M_NF)
# La suma de de SPXY ponderada por los grados de libertad es la
# covarianza.
length(SPXY)
## [1] 218
sum(is.na(SPXY))
## [1] 4
COV <- sum(SPXY, na.rm = TRUE)/213
COV
## [1] 1.315
# La función << var() >> cuando toma como argumento una matriz, tabla o
# data.frame con dos o mas columnas numéricas el resultado es una matriz
# de varianzas y covarianzas. Las varianzas están en la diagonal y las
# covarianzas fuera de la diagonal. La matriz es en espejo asi que por
# encima y debajo de la diagonal existe la misma información
VARCOVAR <- var(cbind(FLORES, NUM_FLORES), na.rm = TRUE)
var(DATOS[, 9:13], na.rm = TRUE)
## TALLO HOJAS FLORES CAPSULAS SEMILLAS
## TALLO 77.153 36.228 2.1790 22.0563 46.476
## HOJAS 36.228 18.578 1.0238 10.9561 23.694
## FLORES 2.179 1.024 0.0793 0.6191 1.313
## CAPSULAS 22.056 10.956 0.6191 7.7681 16.679
## SEMILLAS 46.476 23.694 1.3127 16.6786 37.050
VARCOVAR[2, 1]/(S_F * S_NF)
## [1] 0.9066
# Usamos la función << cor() >>
cor(FLORES, NUM_FLORES, use = "na.or.complete")
## [1] 0.9073
# Entendiendo el ámbito de los coeficientes de correlación
# -------------------------
# Cuando calculamos la correlación sobre la misma variable la congruencia
# de cambio entre las dos variables es perfecta y en el cálculo de la
# covarianza obtenemos. 1) En el numerado (X-X_MEDIA)*(X-X_MEDIA) =
# (X-X_MEDIA)^2 que es la suma de los cuadrados de X y entre los grados de
# libertad tenemos la varianza de X. En el denominador tenemos la
# desviación estándar de X que multiplica a la desviación estándar de X,
# es decir la desviación estándar de X al cuadrado, lo que también nos da
# la varianza de X. Por lo tanto la correlación de dos variables
# idénticas da 1.
# Por el contrario cuando X y Y son exactamente inversas el numerador del
# cálculo de la covarianza es la varianza misma pero con signo negativo ya
# que en la suma de los producto de XY todos los valores serán negativos
# (o iguales a cero si coinciden con el valor de la media)
# Considere estos dos vectores
X <- c(1, 2, 3, 4, 5)
Y <- c(5, 4, 3, 2, 1)
# Ambos vectores tienen la misma media ya que contienen los mismo valores
M_X <- mean(X)
M_Y <- mean(Y)
# Y en las parejas de datos de las diferencias que se multiplican siempre
# uno de los elementos será negativo excepto cuando algún datos coinciden
# con el valor de la media
cbind((X - M_X), (Y - M_Y))
## [,1] [,2]
## [1,] -2 2
## [2,] -1 1
## [3,] 0 0
## [4,] 1 -1
## [5,] 2 -2
# Esto resulta en el estimado de la varianza pero con signo negativo, y
# por lo tanto al dividirlo por le producto de las desviaciones estándares
# de cada vector de datos, que al contener exactamente los mismos valores
# dan como resultado la varianza de cualquiera de los vectores, y por lo
# tanto el cociente resultante del cálculo de la correlación será -1
# Cuando no existe relación entre las variables el numerador tiende a cero
# y en consecuencia también el coeficiente de correlación.
# Para conocer si el coeficiente de correlación es significativo debemos
# dividir este valor entre su error estándar, y esto lo podemo obetner
# directamente con la función <<cor.test() >>
cor.test(FLORES, NUM_FLORES, na.rm = TRUE)
##
## Pearson's product-moment correlation
##
## data: FLORES and NUM_FLORES
## t = 31.42, df = 212, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.8803 0.9284
## sample estimates:
## cor
## 0.9073
# Regresión lineal simple ------
# La regresión lineal simple es un procedimiento mediante el cual se busca
# ajustar una recta a un conjunto de datos definidos por dos vectores
# numéricos en las que las observaciones se presentan como parejas de
# datos. La recta que se ajusta debe minimizar la distancia de los datos a
# la recta. Para poder definir una recta necesitamos dos elementos. Una
# inclinación, pendiente y un punto cualquiera, que en este caso de la
# regresión es el punto en el que la recta interceptaría el eje Y. Una
# propiedad de la recta que minimiza la distancia a los datos es que está
# debe pasar por el punto definido por las dos medias. Por lo tanto de la
# expresión siguiente donde X y Y representan las medias de las variables
# solo debemos estimar << a >> o << b >> y despejar para conocer el valor
# del otro parámetro.
# Y = a + bX
# El mas fácil de estimar es la pendiente que es una medida de cambio
# entre las dos variables, es decir por cada unidad en el eje X cuanto
# cambia en el eje Y. Entonces podríamos calcular todos los cocientes
# entre las parejas de datos y una mediad de tendencia central puede ser
# nuestro estimado de pendiente. Aquí calculamos dos, la media y la
# mediana
P_MEDIA <- mean(FLORES/NUM_FLORES, na.rm = TRUE)
P_MEDIANA <- median(FLORES/NUM_FLORES, na.rm = TRUE)
# Despejando de la formula Y = a + bX, donde << Y >> y << X >> son los
# valores promedios tenemos que a = Y - bX, y calculamos los interceptos a
# partir de los estimados de la pendiente
I_M <- M_F - P_MEDIA * M_NF
I_Q2 <- M_F - P_MEDIANA * M_NF
# Aquí creamos un gráfico de dispersión usando la sintaxis con el tilde
# que R interpreta como <<modelado en función de ... >> y graficará el
# peso de las flores << FLORES >> en función del número de flores <<
# NUM_FLORES >>
plot(FLORES ~ NUM_FLORES)
# Y usamos la función << abline() >> para pintar las lineas de regresión
# estimadas. Esta función usa dos argumentos, un valor de intercepto y un
# valor de pendiente.
abline(I_M, P_MEDIA, col = "yellow")
abline(I_Q2, P_MEDIANA, col = "pink")
# Con la función << points() >> colocamos el punto definido por los
# promedios de las dos variables, en color rojo y que despliegue un punto
# sólido << pch=19>>
points(M_NF, M_F, col = "red", pch = 19)
# Note como las dos rectas pasan por el punto definido por los valores
# promedio, pero son malos predictores de los datos ya que por una parte
# subestima los valores pequeños y por otra parte sobre estima los valores
# grandes. Además esta forma de calcular la pendiente tiene el defecto de
# que nunca puede producir un valor de pendiente negativa que describiría
# una relación inversa entre dos variables
# Usando la suma de los productos de XY, que como ya vimos es el producto
# de los desvíos respecto de las medias y dividido por la suma de
# cuadrados de X podemos estimar la pendiente, teniendo la posibilidad de
# que el valor resultante sea positivo o negativo según se la relación
# entre las variables
SPXY
## [1] 3.166598 0.964412 0.797900 2.901482 -1.265774 0.589761 0.839528
## [8] 0.704552 0.349668 0.328040 0.761296 0.814226 0.165482 3.270412
## [15] 1.157668 6.962505 0.764366 1.766226 0.262180 0.475621 -0.004751
## [22] -0.004751 0.180645 -0.006425 -0.085123 -0.043262 0.048459 0.327621
## [29] 1.797528 -0.039634 1.735900 1.040924 1.014133 1.157668 1.971342
## [36] 5.000505 2.338040 0.032831 1.391156 0.910970 0.023621 0.492412
## [43] NA 3.475528 0.031947 0.765854 0.377296 -0.021448 1.104459
## [50] 5.394552 0.032087 1.391156 2.194505 0.717482 0.704552 0.569156
## [57] 0.988273 0.137110 0.125528 2.149994 0.006970 0.607528 5.375900
## [64] 2.603994 4.845947 7.361435 4.016784 1.916505 2.729900 3.817389
## [71] 3.680645 0.839528 0.701110 1.350970 0.922784 1.797528 2.675342
## [78] 0.108738 0.732738 1.402598 2.176459 1.566831 1.878412 5.464273
## [85] 3.817389 -0.009774 1.346319 0.241528 0.191156 0.572366 -0.051634
## [92] 0.392924 -0.006425 0.187994 2.586691 0.645900 0.080366 -0.001402
## [99] 1.402598 0.530784 1.040924 1.454226 0.606226 0.684273 0.454040
## [106] 0.717482 0.001900 1.216040 0.959342 4.364366 1.089296 1.274412
## [113] 1.674273 1.402598 0.645900 1.797528 0.185435 0.964412 NA
## [120] 1.047668 1.402598 1.551017 5.316180 2.071900 2.168645 1.505854
## [127] 1.089296 -0.051634 0.807435 7.361435 0.865807 2.176459 0.147063
## [134] 0.765854 0.523994 NA 0.964412 5.498505 1.247714 0.669482
## [141] -0.003076 1.454226 0.964412 2.808273 0.176366 0.922784 0.088784
## [148] 3.954133 0.669482 0.637854 0.881156 1.006040 1.427761 1.797528
## [155] 1.797528 1.242877 0.922784 1.797528 1.797528 1.402598 0.589761
## [162] 1.402598 1.551017 1.797528 1.402598 1.089296 1.006040 1.350970
## [169] 1.489389 1.047668 1.181249 1.454226 1.454226 1.735900 1.797528
## [176] 0.989575 0.631389 1.735900 1.454226 1.247714 1.242877 1.797528
## [183] 1.489389 0.937947 1.299342 1.674273 1.454226 0.881156 1.489389
## [190] 0.922784 1.402598 1.350970 1.299342 0.797900 0.542970 1.144459
## [197] 0.384831 0.479435 1.299342 0.922784 0.783063 1.181249 0.964412
## [204] 1.350970 0.606226 0.673017 0.669482 0.783063 0.589761 -0.117913
## [211] 0.436180 NA 0.094970 0.922784 0.989575 0.036924 1.299342
## [218] 0.473296
SPX <- (NUM_FLORES - M_NF)^2
PEND <- sum(SPXY, na.rm = TRUE)/sum(SPX, na.rm = TRUE)
# Otra vez usando la fórmula donde a = Y - bX obtenemos el intercepto
# calculado para esta nueva pendiente
I <- M_F - PEND * M_NF
# Al pintar la linea que describen este intercepto y pendiente vemis que
# es un mucho mejor modelo de los datos ya que no tiende a sub o sobre
# estimar en ninguna parte de la linea de regresión.
plot(FLORES ~ NUM_FLORES)
abline(I, PEND, col = "pink")
# La función << lm() >> contiene toda la codificación para poder llevar
# acabo una regresión lineal y utiliza la notación con tilde, peso de las
# flores modelado en función del número de flores.
M <- lm(FLORES ~ NUM_FLORES)
# El objeto que se genera el del tipo lm,
class(M)
## [1] "lm"
# y la misma función << summary() >> que habíamos utilizado para describir
# un vector de datos nos sirve para extraer los detalles del modelo de
# regresión.
summary(M)
##
## Call:
## lm(formula = FLORES ~ NUM_FLORES)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5018 -0.0697 -0.0150 0.0385 0.6789
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.05439 0.01397 -3.89 0.00013 ***
## NUM_FLORES 0.04968 0.00158 31.42 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.119 on 212 degrees of freedom
## (4 observations deleted due to missingness)
## Multiple R-squared: 0.823, Adjusted R-squared: 0.822
## F-statistic: 987 on 1 and 212 DF, p-value: <2e-16
# Volvamos a graficar los datos
plot(FLORES ~ NUM_FLORES)
# Y agreguemos la línea de regresión esta vez utilizando otra notación.
# Como es una regresión lineal simple vasta con especificar el nombre del
# modelo que al contener información del intercepto y la pendiente la
# función << abline() >> sabrá que hace!
abline(M, col = "red")
# Una función muy útil cuando ajustamos modelos es la función << predict()
# >>
PRED <- predict(M)
# Esta función nos muestra cuanto es el estimado de Y para cada valor de X
# observado.
PRE
## Error: object 'PRE' not found
# La mayor utilidad de esta función radica en que nos puede dar la
# predicción para cualquier valor de X, no necesariamente un valor de los
# datos observados. Así podemos solicitarle la predicción para una
# secuencia de valores que cubran el ámbito de los datos para poder pintar
# la linea de regresión. Esto parece algo absurdo ya que como vimos la
# función << abline() >> justo hace esto. Aquí el uso de << predict() >>
# es redundante, pero en modelos mas complejos con mas de una pendiente y
# mas de un intercepto la función << predict() >> es la única vía para
# pintar la linea de regresión. Para usar la función debemos creear nuevos
# datos, y estos deben tener el mismo nombre que la variable usada en el
# modelo. Esto lo logramos creando un data.frame en el que la variable
# utilizada en el modelo la igualamos a una secuencia de datos. Usualmente
# unos 30 valores son suficientes para tener lineas sin brincos o
# escalones.
SEQ <- seq(0, 20, length = 30)
NUEVOS_DATOS <- data.frame(NUM_FLORES = SEQ)
# Ahora si, teniendo los nuevos datos usamos la función << predict() >>
# utilizando el nombre del modelo como argumento, y el data.frame con los
# nuevos datos, y tendremos el estimado para estos nuevos 30 valores
PRED <- predict(M, NUEVOS_DATOS)
PRED
## 1 2 3 4 5 6 7 8
## -0.05439 -0.02013 0.01414 0.04840 0.08266 0.11692 0.15119 0.18545
## 9 10 11 12 13 14 15 16
## 0.21971 0.25398 0.28824 0.32250 0.35677 0.39103 0.42529 0.45955
## 17 18 19 20 21 22 23 24
## 0.49382 0.52808 0.56234 0.59661 0.63087 0.66513 0.69940 0.73366
## 25 26 27 28 29 30
## 0.76792 0.80218 0.83645 0.87071 0.90497 0.93924
# Note como PRED es un objeto de clase numérica
class(PRED)
## [1] "numeric"
# Con esto podemos añadir la linea graficando cada valor en la secuencia
# contra el valor estimado.
lines(SEQ, PRED, col = "red")
# Además cada estimado de Y lo podemos solicitar acompañado de su error
# estándar lo cual nos permitirá trazar el intervalo de confianza. La
# solicitud de error estándar la hacemos con el argumento << se=TRUE >>
# dentro de a misma función << predict() >>
PRED <- predict(M, NUEVOS_DATOS, se = TRUE)
# Ahora note como el objeto << PRED >> ya no es numérico, ahora es del
# tipo << 'list' >>
class(PRED)
## [1] "list"
# con la función << ls() >> podemos conocer los componente de esta lista.
ls(PRED)
## [1] "df" "fit" "residual.scale" "se.fit"
# Podemos generar dos objetos a partir de esta predicción, uno
# correspondiente al componente << $fit >> que contiene las predicciones
# del moldeo, es decir la recta de regresión.
LIN_REG <- PRED$fit
# Y un segundo objeto a partir del componente << $se.fit >> que contiene
# los errores estándares de las predicciones
EE_REG <- PRED$se.fit
# Con esta información podemos incluir los límites del intervalo de
# confianza de la linea de regresión de la siguiente manera
lines(SEQ, LIN_REG + 1.96 * EE_REG, col = "red", lty = 2)
lines(SEQ, LIN_REG - 1.96 * EE_REG, col = "red", lty = 2)
# Note que sumamos y restamos de las predicciones 1.96 veces el error
# estándar con base en la distribución normal. Sin embargo, es preferible
# ser conservador y utilizar la distribución t en lugar de la z para
# calcular los intervalos de confianza. Lo que se requiere es saber cuanto
# vale los cuantiles de t ,positivo y negativo (es el mismo valor pero con
# signo opuesto), entre los cuales se encuentra el 95% del área bajo la
# curva de la distribución.
# Esto lo logramos con la función probabilística de t para cuantiles, <<
# qt() >> que usa como argumentos la probabilidad para la cual queremos
# saber el valor del cuantil y los grados de libertad. En este caso
# usaremos el valor de probabilidad de 0.025, y al ser a dos colas
# tendremos el 0.05. Los grados de libertad son 212.
qt(0.025, 212)
## [1] -1.971
# El cuantil de t es igual a -1.971217 y los límites del intervalo de
# confianza se calculan y grafican de la siguiente forma
T <- qt(0.025, 212)
lines(SEQ, LIN_REG + T * EE_REG, col = "green", lty = 2)
lines(SEQ, LIN_REG - T * EE_REG, col = "green", lty = 2)
# En este caso la diferencia entre el cálculo del intervalo de confianza
# con z o con t tiene una mínima diferencia ya que con 212 grados de
# libertad la distribución de t es prácticamente igual a la normal.