Estudiando Agricultura en Taiwán: r = 0.816; graficar para pensar

Existen conocidos ejemplos de científicos que empleaban el dibujo para pensar —Ramón y Cajal, se nos viene a la mente.

# Fuente: https://gist.github.com/amoeba/7576126
library(ggplot2)
theme_set(theme_bw(base_size=18))
anscombe_m <- data.frame()

# Código para agrupar datos y gráficos 
for(i in 1:4)
  anscombe_m <- rbind(anscombe_m, data.frame
                      (set = i, 
                        x  = anscombe[,i], 
                        y  = anscombe[,i+4]))

# Gráfico
ggplot(anscombe_m, aes(x, y)) + geom_point(size=5, color="red", fill="purple", shape=21) + 
  geom_smooth(method="lm", fill=NA, fullrange=TRUE) + facet_wrap(~set, ncol=2)

Graficar datos obtenidos de experimentos, aunque menos artísticos que los creados por Cajal, puede encerrar abundante información que ayude a visualizar y pensar cómo proceder con el análisis posterior.

Gráficos de diagnósticos de una regresión lineal, generados en R.

par(mfrow = c(2, 2))
plot(modelo1, main = "Modelo 1")

En 1973, Frank Anscombe publicó un artículo científico para ilustrar la importancia de graficar datos, en donde las correlaciones son idénticas (r = 0.816) pero los gráficos son diferentes, hoy conocido como “Cuarteto de Anscombe”.

library(ggplot2)

gg <- ggplot(anscombe.data, aes(x = x, y = y))
gg <- gg + geom_point(color = "black")
gg <- gg + facet_wrap(~Set, ncol = 2)
gg <- gg + geom_smooth(formula = y ~ x, method = "lm", se = FALSE, data = anscombe.data)
gg

Cógidos y gráficos producidos en R, replicando “Graph in Statistical Analysis”.

# Cuarteto de Anscombe para Regresión Lineal Simple
# Fuente: Cuarteto de Anscombe de 'Idéntica' Regresión Lineal Simple 
# https://rstudio-pubs-static.s3.amazonaws.com/52381_36ec82827e4b476fb968d9143aec7c4f.html

#### Exposición de Datos
knitr::kable(anscombe)

x1	x2	x3	x4	y1	y2	y3	y4
10	10	10	8	8.04	9.14	7.46	6.58
8	8	8	8	6.95	8.14	6.77	5.76
13	13	13	8	7.58	8.74	12.74	7.71
9	9	9	8	8.81	8.77	7.11	8.84
11	11	11	8	8.33	9.26	7.81	8.47
14	14	14	8	9.96	8.10	8.84	7.04
6	6	6	8	7.24	6.13	6.08	5.25
4	4	4	19	4.26	3.10	5.39	12.50
12	12	12	8	10.84	9.13	8.15	5.56
7	7	7	8	4.82	7.26	6.42	7.91
5	5	5	8	5.68	4.74	5.73	6.89

anscombe.1 <- data.frame(x = anscombe[["x1"]], y = anscombe[["y1"]], Set = "Anscombe Set 1")
anscombe.2 <- data.frame(x = anscombe[["x2"]], y = anscombe[["y2"]], Set = "Anscombe Set 2")
anscombe.3 <- data.frame(x = anscombe[["x3"]], y = anscombe[["y3"]], Set = "Anscombe Set 3")
anscombe.4 <- data.frame(x = anscombe[["x4"]], y = anscombe[["y4"]], Set = "Anscombe Set 4")


#### Estadísticas descriptivas
# Media
anscombe.data <- rbind(anscombe.1, anscombe.2, anscombe.3, anscombe.4)
aggregate(cbind(x, y) ~ Set, anscombe.data, mean)

##              Set x        y
## 1 Anscombe Set 1 9 7.500909
## 2 Anscombe Set 2 9 7.500909
## 3 Anscombe Set 3 9 7.500000
## 4 Anscombe Set 4 9 7.500909

# Desviación típica
aggregate(cbind(x, y) ~ Set, anscombe.data, sd)

##              Set        x        y
## 1 Anscombe Set 1 3.316625 2.031568
## 2 Anscombe Set 2 3.316625 2.031657
## 3 Anscombe Set 3 3.316625 2.030424
## 4 Anscombe Set 4 3.316625 2.030579

#### Correlación entre variables x y
library(plyr)

# Función para agrupar datos
correlation <- function(data) {
  x <- data.frame(r = cor(data$x, data$y)) 
  return(x)
}

# Idéntica correlación para las cuatro series de datos: r = 0.816
ddply(.data = anscombe.data, .variables = "Set", .fun = correlation)

##              Set         r
## 1 Anscombe Set 1 0.8164205
## 2 Anscombe Set 2 0.8162365
## 3 Anscombe Set 3 0.8162867
## 4 Anscombe Set 4 0.8165214

### Modelos de regresión lineal simple 
modelo1 <- lm(y ~ x, subset(anscombe.data, Set == "Anscombe Set 1"))
modelo2 <- lm(y ~ x, subset(anscombe.data, Set == "Anscombe Set 2"))
modelo3 <- lm(y ~ x, subset(anscombe.data, Set == "Anscombe Set 3"))
modelo4 <- lm(y ~ x, subset(anscombe.data, Set == "Anscombe Set 4"))


#### Resultados de modelos lineales
summary(modelo1)

## 
## Call:
## lm(formula = y ~ x, data = subset(anscombe.data, Set == "Anscombe Set 1"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.92127 -0.45577 -0.04136  0.70941  1.83882 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   3.0001     1.1247   2.667  0.02573 * 
## x             0.5001     0.1179   4.241  0.00217 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.237 on 9 degrees of freedom
## Multiple R-squared:  0.6665, Adjusted R-squared:  0.6295 
## F-statistic: 17.99 on 1 and 9 DF,  p-value: 0.00217

summary(modelo2)

## 
## Call:
## lm(formula = y ~ x, data = subset(anscombe.data, Set == "Anscombe Set 2"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.9009 -0.7609  0.1291  0.9491  1.2691 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)    3.001      1.125   2.667  0.02576 * 
## x              0.500      0.118   4.239  0.00218 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.237 on 9 degrees of freedom
## Multiple R-squared:  0.6662, Adjusted R-squared:  0.6292 
## F-statistic: 17.97 on 1 and 9 DF,  p-value: 0.002179

summary(modelo3)

## 
## Call:
## lm(formula = y ~ x, data = subset(anscombe.data, Set == "Anscombe Set 3"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.1586 -0.6146 -0.2303  0.1540  3.2411 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   3.0025     1.1245   2.670  0.02562 * 
## x             0.4997     0.1179   4.239  0.00218 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.236 on 9 degrees of freedom
## Multiple R-squared:  0.6663, Adjusted R-squared:  0.6292 
## F-statistic: 17.97 on 1 and 9 DF,  p-value: 0.002176

summary(modelo4)

## 
## Call:
## lm(formula = y ~ x, data = subset(anscombe.data, Set == "Anscombe Set 4"))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.751 -0.831  0.000  0.809  1.839 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   3.0017     1.1239   2.671  0.02559 * 
## x             0.4999     0.1178   4.243  0.00216 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.236 on 9 degrees of freedom
## Multiple R-squared:  0.6667, Adjusted R-squared:  0.6297 
## F-statistic:    18 on 1 and 9 DF,  p-value: 0.002165

#### Cuatro gráficos juntos
library(ggplot2)

gg <- ggplot(anscombe.data, aes(x = x, y = y))
gg <- gg + geom_point(color = "black")
gg <- gg + facet_wrap(~Set, ncol = 2)
gg <- gg + geom_smooth(formula = y ~ x, method = "lm", se = FALSE, data = anscombe.data)
gg

# Nota: 
# Observar los gráficos de dispersión con sus respectivas líneas de regresión en azul;
# luego, los correspondientes gráficos de diagnósticos. Para iniciar la interpretación, 
# el gráfico de "Residuals vs Fitted" o "Residuos vs Ajustados" debe asemejarse a un cielo 
# estrellado - con puntos totalmente aleatorios.


#### Gráficos de diagnósticos para los cuatro modelos
par(mfrow = c(2, 2))
plot(modelo1, main = "Modelo 1")

plot(modelo2, main = "Modelo 2")

plot(modelo3, main = "Modelo 3")

plot(modelo4, main = "Modelo 4")

#### Otra opción de gráfico
# Fuente: https://gist.github.com/amoeba/7576126
library(ggplot2)
theme_set(theme_bw(base_size=18))
anscombe_m <- data.frame()

# Código para agrupar datos y gráficos 
for(i in 1:4)
  anscombe_m <- rbind(anscombe_m, data.frame
                      (set = i, 
                        x  = anscombe[,i], 
                        y  = anscombe[,i+4]))

# Gráfico
ggplot(anscombe_m, aes(x, y)) + geom_point(size=5, color="red", fill="purple", shape=21) + 
  geom_smooth(method="lm", fill=NA, fullrange=TRUE) + facet_wrap(~set, ncol=2)

Hoy en día, hasta los paquetes más sofisticados de visualización de datos, siguen empleando el Cuarteto de Anscombe para enfatizar la importancia de graficar—ver sitio.

Archivo con cógidos completos: Cuarteto de Anscombe.r

Estudiando Agricultura en Taiwán

Páginas

mayo 09, 2017

r = 0.816; graficar para pensar

No hay comentarios:

Publicar un comentario