1 Correlación: teoría

Este tema es parte del Curso Fundamental de Posgrado: Análisis Estadístico en Ciencias Biológicas Utilizando R, semestre 2017-1, de la Universidad Nacional Autónoma de México, impartido entre Agosto-Diciembre 2016 en el Centro de Ciencias Genómicas. Para más información consultar la página del curso en: http://www.ccg.unam.mx/~vinuesa/R4biosciences/.

Una presentación web de este documento está disponible aquí: http://www.ccg.unam.mx/~vinuesa/R4biosciences/docs/Tema8_correlacion_presentacionR.html#

La parte teórica está basada en (Crawley 2015), (Everitt and Hothorn 2014) y (A. P. Field, Miles, and Field 2012).

Este documento está aún en construcción y es generado con R (R Core Team 2016), rstudio (RStudio Team 2016), knitr (Xie 2016), rmarkdown (Allaire et al. 2016), pandoc (MacFerlane 2016) y \(LaTeX\).

1.1 Introducción: el concepto de correlación

La correlación es una medida de la relación (covariación) lineal entre dos variables cuantitativas contínuas (x, y). La manera más sencilla de saber si dos variables están correlacionadas es determinar si co-varían (varían conjuntamente). Es importante hacer notar que esta covariación no implica necesariamente causalidad, la correlación puede ser fortuita, como en el caso clásico de la correlación entre entre el número de venta de helados e incendios, debido al efecto de una tercera variable, la temperatura ambiental.

La correlación es en esencia una medida normalizada de asociación o covariación lineal entre dos variables. Esta medida o índice de correlación \(r\) puede variar entre -1 y +1, ambos extremos indicando correlaciones perfectas, negativa y positiva respectivamente. Un valor de \(r\) = 0 indica que no existe relación lineal entre las dos variables. Una correlación positiva indica que ambas variables varían en el mismo sentido. Una correlación negativa significa que ambas variables varían en sentidos opuestos. Lo interesante del índice de correlación es que \(r\) es en sí mismo una medida del tamaño del efecto, que suele interpretarse de la siguiente manera:

  • correlación despreciable: \(r\) < |0.1|
  • correlación baja: |0.1| < \(r\) <= |0.3|
  • correlación mediana : |0.3| < \(r\) <= |0.5|
  • correlación fuerte o alta: \(r\) > |0.5|

La siguiente figura muesta ejemplos de pares de variables con correlación positiva moderada, negativa fuerte, así como correlación despreciable. ¿Sabrías identificar cada caso? Estos son datos que provienen del set state.x77 del paquete \(datasets\) que viene con la instalación de base de R. Se trata de una matriz de 50 filas y 8 columnas con estadísticas para 50 estados de EU relativos al tamaño de la población, renta per cápita, % de analfabetismo, expectancia de vida, tasa de asesinato, % de graduados de preparatoria, número promedio de días con heladas y área del estado en millas cuadradas. Esta información la pueden ver con el comando \(help("state.x77")\).

Analizemos visualmente las relaciones entre un subconjunto de las variables de state.x77 para afianzar los conceptos clave. Nótese que en las Figs. C-D, además de los gráficos de dispersión (“scatterplots”), se muestra la recta de regresión correspondiente al ajuste de un modelo lineal a los datos, con el fin de visualizar mejor la desviación de los puntos con respecto al modelo lineal.

1.2 Definiciones formales

La correlación se define en términos de la varianza (\(s^2\)) de las variables x e y, así como de la covarianza \(cov\) de x,y. Es por tanto una medida de la variación conjunta de ambas variables (\(cov(x,y)\)).

1.2.1 varianza \((s^2)\)

La varianza de una muestra representa el promedio de la desviación de los datos con respecto a la media \[\displaystyle Varianza(s^2) = \frac{\sum(x_i-\bar{x})^2} {N-1} = \frac{\sum(x_i-\bar{x})(xi-\bar{x})} {N-1}\\\]

1.2.2 covarianza \(cov(x,y)\)

La covarianza entre dos variables \(x\) e \(y\) es una medida de la relación “promedio” éstas. Es la desviación promedio del producto cruzado entre ellas: \[\displaystyle cov(x,y) = \frac{\sum(x_i-\bar{x})(yi-\bar{y})} {N-1}\\\]

Veamos un ejemplo de dos variables que co-varían (respuesta ~ dosis en 5 pacientes), es decir, que están correlacionadas.
Los datos: dosis=(8,9,10,13,15); resp=(5,4,4,6,8);

1.2.2.1 Cálculo “a mano” de la covarianza de dos variables

\(\displaystyle cov(dosis,resp) = \frac{\sum(x_i-\bar{x})(yi-\bar{y})} {N-1} =\frac{(-0.4)(-3)+(-1.4)(-2)+(-1.4)(-1)+(0.6)(2)+(2.6)(4)}{4} = \frac{17}{4} = 4.25\\\)

Un valor de covarianza positivo indica que ambas variables se desvían de la media en la misma dirección, mientras que uno negativo indica que las desviaciones acontecen en sentidos opuestos.

1.2.2.2 Cálculo en R de las medias y desviaciones para cada variable, así como el coeficiente de covariación

# genermos los vectores dosis y resp
dosis <- c(5,4,4,6,8)
resp <- c(8,9,10,13,15)
# calculemos la dosis y respuesta medias
dosis.mean <- mean(dosis); resp.mean <- mean(resp);  
cat("dosis media =", dosis.mean, "; resp.mean =", resp.mean)
## dosis media = 5.4 ; resp.mean = 11
# cálculo de las desviaciones de cada dosis y respuesa con respecto a sus  
# valores promedio
dosis.dev <- dosis - mean(dosis); resp.dev <- resp - mean(resp);  
cat("dosis.dev =", dosis.dev, "; resp.dev =", resp.dev)
## dosis.dev = -0.4 -1.4 -1.4 0.6 2.6 ; resp.dev = -3 -2 -1 2 4
# Cálculo del coef. de covariación
Covar <- sum((dosis.dev)*(resp.dev))/(length(dosis)-1); cat("cov =", Covar)
## cov = 4.25
# cálculo de la covariación entre dosis y respuesta con cov(x,y)
cov(dosis,resp)
## [1] 4.25

El problema de usar la covarianza como medida de relación entre variables estriba en que depende de la escala de las medidas usadas. Es decir, la covarianza no es una medida estandarizada. Por tanto la covarianza no puede ser usada para comparar las relaciones entre variables medidas en diferentes unidades.

1.2.3 El coeficiente de correlación de Pearson \(r\) = coef. de covariación estandarizado

Para resolver el problema de dependencia de la escala o unidades de las mediciones (valores), necesitamos una unidad a la cual pueda convertirse cualquier medida. Esta unidad de medida libre de escala es la desviación estándar (s ó \(\sigma\)). Al igual que la varianza, mide la desviación promedio de los datos con respecto a la media aritmética por no ser otra cosa que la \(\sqrt{varianza}\) ó \(\sqrt{s^2}\). Al dividir cualquier distancia de la media por la desviación estándar, obtendremos una distancia en unidades de desviación estándar.

Por tanto, para normalizar la covarianza la tenemos que dividir por la desviación estándar. Como la covarianza se calcula para dos variables cov(x,y), tenemos que calcular la desviación estándar para cada variable, multiplicándolas entre ellas, es decir:

\[\displaystyle \text{Coef. de correlación de Pearson} (r) = \frac{cov(x,y)} {s_xs_y} =\frac{\sum(x_i-\bar{x})(yi-\bar{y})} {(N-1)s_xs_y}\\\]

1.2.3.1 Cálculo del coeficiente de correlación de Pearson en R:

# uso básico de la función cor()
cor(dosis,resp)
## [1] 0.8711651
# comprobamos el resultado usando la fórmula de r indicada arriba
cov(dosis,resp)/(sd(resp)*sd(dosis))
## [1] 0.8711651

1.2.4 Correlaciones parciales

Permiten evaluar la correlación entre dos variables (Var.1 y Var.2) considerando el efecto (varianza) de una tercera (Var.3) o más variables.

Eliminando la varianza compartida por las variables de interés con la o las variables auxiliares, obtenemos una medida de \(r\) que refleja los efectos de las variables de interés primario.

En R podemos hacer análisis de correlación parcial usando la función \(pcor()\) del paquete \(ggm\). Veremos en la prática el uso de las funciones \(ggm::pcor()\) y $ggm::cpor.test().

1.2.5 Supuestos hechos por el estadístico de correlación de Pearson \(r\)

  • Ambas deben ser variables cuantitativas contínuas (medidas de intervalo)
  • Si queremos correr tests de significancia, las variables deben estar normalmente distribuidas

1.2.6 El coeficiente de correlación no paramétrico de Kendall \(\tau\)

El coeficiente de correlación \(\tau\) de Kendall es no paramétrico, es decir, se puede usar cuando se viola el supesto de distribución normal de las variables a comparar. La correlación \(\tau\) de Kendall es particularmente adecuada cuando tenemos un set de datos pequeño con muchos valores en el mismo rango o clase. Se puede usar por ejemplo con datos categóricos codificados binariamene (0,1). Estudios estadísticos han demostrado que el coeficiente de correlación \(\tau\) de Kendall es un mejor estimador de la correlación en la población que el coeficiente de correlación no paramétrico de Spearman \(\rho\), por lo que se recomienda usar \(\tau\) para análisis de datos no paramétricos. En R se puede estimar la correlación \(\tau\) o \(\rho\) cambiando el valor del argumento method=“kendall” o method=“spearman” de la función \(cor\), en la que por defecto method=“pearson”. Por ejemplo: \(cor(x,y, method="kendall")\)

1.2.7 El coeficiente de determinación \(R^2\)

El coeficiente de correlación elevado al cuadrado es el coeficiente de determinación, \(R^2\), que mide la cantidad de variación en una variable que es compartida por otra. Vimos en el ejemplo anterior que la \(r\) para cor(dosis,resp) era de 0.8711651, y por tanto \(R^2 = 0.7589286\). Por tanto podemos decir que la respuesta comparte un ~76% de la variación mostrada por la dosis. Tengan en cuenta de nuevo que compartir variabilidad no implica necesariamente causalidad.

# Cálculo del coeficiente de determinación, expresado como porcentaje, de la correlación  
# entre dosis y respuesta.
cat("R^2(dosis,resp) =", round( cor(dosis,resp)^2 * 100, 1), "%.")
## R^2(dosis,resp) = 75.9 %.

1.3 Significancia del coeficiente de correlación \((r)\)

Se trata de probar la hipótesis de que \(r\neq0\), es decir, buscamos rechazar la \(H_0: r = 0\).
Tenemos dos maneras de retar la \(H_0\): i) usando \(z-scores\); ii) mediante estadístico \(t\).

1.3.1 Cálculo de la significancia de \(r\) usando \(z-scores\)

Los \(z-scores\) son útiles ya que conocemos la probabilidad de ocurrencia de un valor de \(z\) determinado, siempre y cuando la distribución de la que proviene sea normal. Dado que es sabido que \(r\) tiene una distribución muestral no normal, tenemos que hacer una transformación de Fisher para normalizarla:

\[\displaystyle z_r = \frac{1}{2}log_e\left(\frac{1+r}{1-r}\right)\\\]

La \(z_r\) resultante tiene un error estándar de:
\[\displaystyle SE_{z_r} = \frac{1}{\sqrt{N-3}}\]

  • Para nuestro ejemplo tenemos que \(r\) = .87, \(z_r = 1.33\) y \(SE_{z_r} = .71\), como se demustra abajo aplicando las dos fórmulas arriba enunciadas:
# el valor de r previamente calculado
r <- round(cor(dosis,resp), 2) # .87

# Cálculo del z-score para r
Zr <- 0.5*(log((1+r)/(1-r)))
cat("z-score para r (Zr): ", Zr)
## z-score para r (Zr):  1.33308
# Cálculo del error estándar de zr
SEzr <- 1/sqrt(5-3)
cat("error estándar de zr (SEzr):", SEzr)
## error estándar de zr (SEzr): 0.7071068

Ahora podemos transformar esta \(r\) ajustada a un \(z-score\), tal y como vimos anteriormente para los scores crudos. Recuerden que para calcular un \(z-score\) que represente el tamaño de la correlación con respecto a un valor particular, simplemente calculamos el \(z-score\) contra dicho valor usando el \(SE_{z_r}\) asociado. Este valor generalmente va a ser 0, ya que queremos probar que \(H_1 \neq H_0\), es decir, trataremos de rechazar la \(H_0: r = 0\). Por tanto \[\displaystyle z = \frac{z_r}{SE_{z_r}}\]

  • Cálculo de la significancia de \(r\) usando \(z\)
z <- Zr/SEzr
cat("valor de z: ", z, "## ¿Será significativo?")
## valor de z:  1.885259 ## ¿Será significativo?
pvalue1sided = pnorm(-abs(z))
cat("valor p de una cola: ", pvalue1sided)
## valor p de una cola:  0.02969742
pvalue2sided <- 2*pnorm(-abs(z))
cat("valor p de dos colas: ", pvalue2sided, "## <<< Generalmente usamos este caso")
## valor p de dos colas:  0.05939484 ## <<< Generalmente usamos este caso

Pueden usar la función \(convert.z.score()\) para convertir \(z-scores\) a \(p-values\), mostrada abajo.

# Vean cómo se importa código de en un archivo almacenado en el disco local
source("./code/func_convert.z.score.R") 

# con este comando podemos imprimir a STDOUT (consola) el contenido del archivo guardado  
# en disco
cat(readLines('./code/func_convert.z.score.R'), sep = '\n')
## # función para convertir z-scores a p-values, tomada de https://www.biostars.org/p/17227/
## convert.z.score<-function(z, one.sided=NULL) {
##   if(is.null(one.sided)) {
##     pval = pnorm(-abs(z));
##     pval = 2 * pval
##   } else if(one.sided=="-") {
##     pval = pnorm(z);
##   } else {
##     pval = pnorm(-z);                                                                                 
##   }
##   return(pval);
## }
# Con estas líneas llamamos a la función que hemos importado con source
pval.2s <- convert.z.score(z)
pval.1s <- convert.z.score(z, 1)

# imprimimos el contenido de las variables que almacenan el resultado devuelto  
# por la función
cat("2-sided pval: ", pval.2s, "| 1-sided pval: ", pval.1s)
## 2-sided pval:  0.05939484 | 1-sided pval:  0.02969742

1.3.2 Cálculo de la significancia de \(r\) mediante el estadístico-\(t\) y la funcion \(cor.test()\)

En R el test de significancia de \(r\) está implementado en la función \(cor.test()\), basado en el estadístico-\(t\). con \(N-2\) grados de libertad, que pueden obtenerse directamente de \(r\):

\[\displaystyle t_r = \frac{r \sqrt{(N-2)}}{\sqrt{1-r^2}}\\\]

#>>> Cálculo a mano de la significancia de r
# Cálculo de r
r <- cor(dosis,resp)

# Cálculo del estadístico-t
tr <- r*(sqrt((5-2)))/sqrt((1-r**2))
cat("valor del estadístico-t: ", tr)
## valor del estadístico-t:  3.073181
# cálculo de la p del estadístico
ptr <- 2*pt(tr, 3, lower.tail = FALSE)
cat("valor de la p para el estadístico-t: ", ptr)
## valor de la p para el estadístico-t:  0.05442624
#>>> Cálculo automático con la función de R cor.test()
(cor.test(dosis, resp))
## 
##  Pearson's product-moment correlation
## 
## data:  dosis and resp
## t = 3.0732, df = 3, p-value = 0.05443
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.0479747  0.9914236
## sample estimates:
##       cor 
## 0.8711651

1.3.3 Análisis de potencia y significancia estadística de \(r\)

El análisis de potencia es un aspecto importante del diseño experimental. Nos permite determinar el tamaño muestral (la famosa \(n\)) requerido para detectar un efecto de un determinado tamaño con un grado determinado de confianza. De modo complementario, tambíen nos permite determinar la probabilidad de detectar un efecto de un tamaño determinado, dados un nivel de confiaza y tamaño de muestra predeterminados. Así por ejemplo, si el nivel de confianza o estima del tamaño del efecto no son satisfactorios para un \(n\) dada, lo único aconsejable es incrementar la \(n\) o abandonar el experimento.

Las siguientes cuatro magnitudes están íntimamente relacionadas:
tamaño muestal (número de observaciones)
tamaño del efecto (medida de la fuerza de un fenómeno; la magnitud del estadístico, \(r\), \(t\), \(chi^2\) …) wikipedia
nivel de significancia = P(error Tipo I) = probabilidad de observar un effecto inexistente (usualmente \(\alpha = .05\))
potencia = 1 - P(error Tipo II) = probabilidad de observar un effecto real (un valor frecuentemente usado es 0.8)

Dadas cualesquiera tres, podemos determinar la cuarta magnitud.

Recordemos los tipos de error: - Error de tipo I: ocurre cuando estimamos que existe un efecto en la población, cuando en realidad no hay tal. - Error de tipo II: ocurre cuando estimamos que no existe un efecto en la población, cuando en realidad sí lo hay.

Los coeficientes de correlación son un ejemplo clásico de tamaños del efecto, por tanto podemos interpretar \(r\) sin la necesidad de \(valores p\), particularmente debido a que los \(valores p\) están relacionados con el tamaño de la muestra. Hay muchos estadísticos que recomiendan no obsesionarse con los famosos \(valores p\) para coeficientes de correlación, o coeficientes de regresión, ya que ambos son en sí estimadores de tamaños del efecto.

1.4 La importancia de visualizar gráficamente los datos antes de someterlos a análisis de correlación: lecciones del cuarteto de Anscombe

De lo explicado arriba debe quedar claro que el análisis de correlación sólo debe aplicarse entre pares de variables con relación lineal entre ellas. Es por tanto importante siempre graficar los datos para visualizar el tipo de relaciones entre las variables y detectar la presencia de valores atípicos o aberrantes (“outliers”).

El cuarteto de Anscombe comprende cuatro conjuntos de datos que tienen las mismas propiedades estadísticas (medias y varianzas), presentando los mismos valores de correlación entre pares de variables, con el mismo ajuste lineal (recta de regresión), pero que son marcadamente distintas al inspeccionarlas gráficamente.

Cada conjunto consiste de once puntos (x, y) y fueron construidos por el estadístico F. J. Anscombe. El cuarteto es una demostración de la importancia de visualizar gráficamente un conjunto de datos antes de someterlos a análisis estadísticos.

\(R\) trae el set de datos de Anscombe en el paquete \(datasets\). Exploremos sus características estadísticas, grafiquemos las relaciones entre los pares de variables x1-y1, x2-y2, x3-y3, x4-y4 y discutamos los resultados.

1.4.1 Código para generar las gráficas y estadísticas del cuarteto de Anscombe

Nota: este código puede ser un poco difícil de entender ya que contiene un bucle for y otros detalles sintácticos que no hemos visto en el curso. Muestro el código para los interesados en ejemplos para aprender estos elementos esenciales del lenguaje. Pero lo importante no es el código, sino las gráficas y análisis estadísticos resultantes, que discutiremos en la siguiente sección.

# visualizemos el conjunto de datos de Anscombe 
# recuerda: str(anscombe) nos muestra la estructura del objeto
anscombe
##    x1 x2 x3 x4    y1   y2    y3    y4
## 1  10 10 10  8  8.04 9.14  7.46  6.58
## 2   8  8  8  8  6.95 8.14  6.77  5.76
## 3  13 13 13  8  7.58 8.74 12.74  7.71
## 4   9  9  9  8  8.81 8.77  7.11  8.84
## 5  11 11 11  8  8.33 9.26  7.81  8.47
## 6  14 14 14  8  9.96 8.10  8.84  7.04
## 7   6  6  6  8  7.24 6.13  6.08  5.25
## 8   4  4  4 19  4.26 3.10  5.39 12.50
## 9  12 12 12  8 10.84 9.13  8.15  5.56
## 10  7  7  7  8  4.82 7.26  6.42  7.91
## 11  5  5  5  8  5.68 4.74  5.73  6.89
summary(anscombe)
##        x1             x2             x3             x4    
##  Min.   : 4.0   Min.   : 4.0   Min.   : 4.0   Min.   : 8  
##  1st Qu.: 6.5   1st Qu.: 6.5   1st Qu.: 6.5   1st Qu.: 8  
##  Median : 9.0   Median : 9.0   Median : 9.0   Median : 8  
##  Mean   : 9.0   Mean   : 9.0   Mean   : 9.0   Mean   : 9  
##  3rd Qu.:11.5   3rd Qu.:11.5   3rd Qu.:11.5   3rd Qu.: 8  
##  Max.   :14.0   Max.   :14.0   Max.   :14.0   Max.   :19  
##        y1               y2              y3              y4        
##  Min.   : 4.260   Min.   :3.100   Min.   : 5.39   Min.   : 5.250  
##  1st Qu.: 6.315   1st Qu.:6.695   1st Qu.: 6.25   1st Qu.: 6.170  
##  Median : 7.580   Median :8.140   Median : 7.11   Median : 7.040  
##  Mean   : 7.501   Mean   :7.501   Mean   : 7.50   Mean   : 7.501  
##  3rd Qu.: 8.570   3rd Qu.:8.950   3rd Qu.: 7.98   3rd Qu.: 8.190  
##  Max.   :10.840   Max.   :9.260   Max.   :12.74   Max.   :12.500
# calculemos las medias y varianza de cada columna  
# Noten el uso de la función sapply
sapply(anscombe, mean)
##       x1       x2       x3       x4       y1       y2       y3       y4 
## 9.000000 9.000000 9.000000 9.000000 7.500909 7.500909 7.500000 7.500909
sapply(anscombe, var)
##        x1        x2        x3        x4        y1        y2        y3 
## 11.000000 11.000000 11.000000 11.000000  4.127269  4.127629  4.122620 
##        y4 
##  4.123249
# veamos ahora las correlaciones entre las viables x1-y1, x2-y2 ...
cor(anscombe[,1:4], anscombe[,5:8])
##            y1         y2         y3         y4
## x1  0.8164205  0.8162365  0.8162867 -0.3140467
## x2  0.8164205  0.8162365  0.8162867 -0.3140467
## x3  0.8164205  0.8162365  0.8162867 -0.3140467
## x4 -0.5290927 -0.7184365 -0.3446610  0.8165214
# de la matriz anterior realmente necesitamos sólo su diagonal
diag(cor(anscombe[,1:4], anscombe[,5:8]))
## [1] 0.8164205 0.8162365 0.8162867 0.8165214
# Grafiquemos las relaciones entre las virables x1-y1, x2-y2 ...  
# El objetivo del código es:  
# 1) mediante par(mfrow=c(2,2), mar=c(6,4,1,1)) controlamos que las  
#    4 gráficas queden en un solo display, como una matriz de 2x2.  
# 2) Mediante dos bucles for anidados, obtenemos los índices para  
#    poder indexar el dataframe asncombe, atendiendo a los números  
#    de las columnas de las variables x,y que nos interesan: x=(1,2,3,4), y=(5,6,7,8).  
# 3) Lo que sigue es llamar a plot para cada par x,y  
# 4) para finalmente graficar la regresión lineal (explicado en el próximo tema).

# La información relevante (r y recta de regresión se imprime como texto sobre la gráfica)  
# Las líneas correspondientes están comentadas en el código
# Si les parece complicada mi solución, vean la que propone la documentación de R  
# tecleando help(anscombe); tiene detalles sintácticos muy interesantes

# par() controla los parámetros gráficos. guardamos en oldpar <- los ajustes originales  
# de par()
oldpar <- par()
# plotea 4 gráficas en 2 filas y 2 columnas, dando  márgenes entre ellas
par(mfrow=c(2,2), mar=c(4,4,1,1)) 
for(i in 1:4){ # 1,2,3,4
    j <- i+4  # 5,6,7,8
    # las siguientes variables capturan información que añadiremos al título de cada  
    # plot (..., main=h)
    
    # función de redondeo, a 2 dígitos
    r <- round(cor(anscombe[,i], anscombe[,j]), 2) 
    # pegamos diferentes elementos en una cadena
    h <- paste("Anscombe plot,,", i, " r =", r, sep = ' ')  
    
    # ajustemos un modelo lineal lm() y capturemos sus coeficientes:  
    # corte eje ordenadas (y) y pendiente. Capturamos los valores y los  
    # pegamos en la fórmula general de la recta: y = a + bx. 
    # lm() ajusta un modelo lineal a los datos
    fit <- lm(anscombe[,j] ~ anscombe[,i])
    
    # con fit$coefficients[[1]] extraemos del objeto fit (una lista),  
    # el primer [[1]] y segundo [[2]]  elemento almacenados en el vector  
    # llamado coefficients. Nótese el uso de [[]] para sacar elementos de una lista!!!
    regr <- paste("y =", round(fit$coefficients[[1]], 2), "+", 
                  round(fit$coefficients[[2]], 1), "x") 
    
    # estas líneas son para los "scatterplots" de los pares de variables
    plot(anscombe[,i], anscombe[,j], main = h, xlab = "", ylab = "")
    
    # y para graficar la recta de regresión
    abline(lm(anscombe[,j] ~ anscombe[,i]), lty=2)
    
    # con mtext() podemos escribir los coeficientes de la fórmula dentro del área  
    # de cada gráfica especificando el lado derecho (side=1), pegado a la derecha  
    # (adj = 1), y dos líneas abajo de la línea central (line = -2)
    mtext(regr, side = 1, adj = 1, line = -2)
}