El paquete gapminder contiene un fichero de datos de
población, esperanza de vida y renta per cápita de los países del mundo
entre 1952 y 2007. Instala el paquete y lleva a cabo los siguientes
gráficos:
Un histograma de la esperanza de vida en 2007 de los países de Europa.
Diagramas de caja con las esperanzas de vida de cada continente en el año 1952.
Un diagrama de dispersión de la renta per cápita y la esperanza de vida de cada país en el año 2007.
Mejora el gráfico anterior representando cada punto de un color diferente en función del continente al que pertenece cada país y representando la renta per cápita en una escala logarítmica.
Solución
(a) Histograma de la esperanza de vida en 2007 de los países de Europa.
library(gapminder)
View(gapminder)
head(gapminder)## # A tibble: 6 × 6
## country continent year lifeExp pop gdpPercap
## <fct> <fct> <int> <dbl> <int> <dbl>
## 1 Afghanistan Asia 1952 28.8 8425333 779.
## 2 Afghanistan Asia 1957 30.3 9240934 821.
## 3 Afghanistan Asia 1962 32.0 10267083 853.
## 4 Afghanistan Asia 1967 34.0 11537966 836.
## 5 Afghanistan Asia 1972 36.1 13079460 740.
## 6 Afghanistan Asia 1977 38.4 14880372 786.dim(gapminder)## [1] 1704 6# Para seleccionar los países de Europa
# defino un vector de Trues y Falses indicando
# qué país de la muestra es europeo:
EuropeTrue = gapminder$continent == "Europe"
# Selecciono las filas de la matriz de datos
# correspondientes a 2007:
True2007 = gapminder$year == "2007" Observación: No podemos dar a ningún objeto de R un nombre que comience con un número, de ahí que haya usado el nombre True2007 y no 2007True.
# Para evaluar la esperanza de vida sólo en los países europeos:
EspVidaEuropa = gapminder$lifeExp[EuropeTrue]
# Evalúo la esperanza de vida en los países europeos y en 2007:
EspVidaEuropa2007 = gapminder$lifeExp[EuropeTrue&True2007]
# Histograma de la esperanza de vida en 2007 de
# los países de Europa:
hist(EspVidaEuropa2007,freq=F,xlab="Esperanza de vida")
Vemos que hay dos grupos: los países cuya esperanza de vida está entre 72 y 78 y aquéllos cuya esperanza de vida está entre 78 y 82.
(b) Diagramas de caja con las esperanzas de vida de cada continente en el año 1952.
True1952 = gapminder$year == "1952"
boxplot(gapminder$lifeExp[True1952] ~ gapminder$continent[True1952], xlab="Continente",ylab="Esperanza de vida")
La tilde ~ en R siempre significa “en función de”.
(c) Un diagrama de dispersión de la renta per cápita y la esperanza de vida de cada país en el año 2007.
GDP2007 = gapminder$gdpPercap[True2007]
EspVida2007 = gapminder$lifeExp[True2007]
plot(GDP2007,EspVida2007,
xlab="Renta per capita",
ylab="Esperanza de vida",
main="Año 2007")
Si escribimos en la consola de R
?pointso equivalentemente
help("points")o en la pestaña de Help, en el recuadro de la lupa,
tecleamos points y pulsamos Intro, podemos ver
parámetros gráficos con los que mejorar el aspecto de los gráficos de R.
Por ejemplo, pch permite cambiar el tipo de punto
dibujado:
plot(GDP2007,EspVida2007,pch=19,
xlab="Renta per capita",
ylab="Esperanza de vida",
main="Año 2007")
Vemos que la esperanza de vida aumenta con el PIB per capita, pero no de manera lineal. De hecho, a mí me ha recordado a la curva de Michaelis-Menten, que describe la velocidad en una reacción enzimática en función de la concentración de sustrato y por eso he tratado de ajustar el mismo tipo de ecuación a estos datos: \[ y = 40 + 40*\frac{x}{1000+x}, \] donde \(x\) = “Renta per capita” e \(y\) = “Esperanza de vida”. ¡Da un ajuste muy bueno para haber determinado los parámetros de la curva a ojo! Es muy interesante observar que la “esperanza de vida” alcanza la saturación a los 80 años aproximadamente (es decir, superado un umbral de producto interior bruto per capita, da prácticamente lo mismo si éste es mayor o menor).
plot(GDP2007,EspVida2007,pch=19,
xlab="Renta per capita",
ylab="Esperanza de vida",
main="Año 2007")
k = 1000
x = 0:50000
y = 40+40*x/(k+x)
lines(x,y,col="red")
(d) Mejora el gráfico anterior representando cada punto de un color diferente en función del continente al que pertenece cada país y representando la renta per cápita en una escala logarítmica.
Veamos qué continentes aparecen en 2007:
continent2007 = gapminder$continent[True2007]
summary(continent2007)## Africa Americas Asia Europe Oceania
## 52 25 33 30 2El continente es una variable cualitativa, o lo que R llama un factor (esta denominzación también se usa en Estadística, en los Modelos Lineales y, más concretamente, en Análisis de la Varianza). Los posibles valores de un factor se denominan niveles:
levels(continent2007)## [1] "Africa" "Americas" "Asia" "Europe" "Oceania"Podemos usar el comando unclass para colorear según los
niveles de un factor:
plot(GDP2007,EspVida2007,pch=19,
xlab="Producto interior bruto per capita",
ylab="Esperanza de vida", main="Año 2007",
col=c("red","green","black","blue","pink")[unclass(continent2007)])
legend("bottomright",pch=16,
col=c("red","green","black","blue","pink"),
legend=levels(continent2007))
Si buscáis R colors en Google veréis la paleta de
colores de R y si tecleáis en la consola de R
colors()veréis la lista de posibles colores.
Representar la renta per cápita en una escala logarítmica quiere decir sustituir el producto interior bruto per capita en el gráfico por su logaritmo. En el problema de regresión de \(y\) sobre \(x\) transformar una o las dos variables (mediante funciones \(g\) y \(f\) repectivamente) puede “linealizar” la relación entre las variables transformadas. Esto quiere decir que, aunque las variables originales no tengan relación lineal, es posible que sí la tengan transformaciones de dichas variables, es decir, \[ g(y) \simeq a+bf(x). \]
plot(log(GDP2007),EspVida2007,pch=16,col=c("red","green","black","blue","pink")[unclass(continent2007)])
legend("topleft",pch=16,col=c("red","green","black","blue","pink"),legend=levels(continent2007))
Personalmente, creo que no ha linealizado bien la relación, pero podemos probar a ajustar una recta de regresión a los datos:
recta = lm(EspVida2007~log(GDP2007))
plot(log(GDP2007),EspVida2007,pch=16,col=c("red","green","black","blue","pink")[unclass(continent2007)])
legend("topleft",pch=16,col=c("red","green","black","blue","pink"),legend=levels(continent2007))
abline(recta)
summary(recta)##
## Call:
## lm(formula = EspVida2007 ~ log(GDP2007))
##
## Residuals:
## Min 1Q Median 3Q Max
## -25.947 -2.661 1.215 4.469 13.115
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.9496 3.8577 1.283 0.202
## log(GDP2007) 7.2028 0.4423 16.283 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.122 on 140 degrees of freedom
## Multiple R-squared: 0.6544, Adjusted R-squared: 0.652
## F-statistic: 265.2 on 1 and 140 DF, p-value: < 2.2e-16Es decir, denotando \(x\) = “Renta per capita” e \(y\) = “Esperanza de vida”, el modelo planteado sería: \[ y \simeq 4.9496 + 7.2028 \log(x), \] que es como el que aparece en Wikipedia en la entrada “Renta per capita”. “Deshaciendo la transformación” gráficamente,
plot(GDP2007,EspVida2007,pch=19,
xlab="Renta per capita", ylab="Esperanza de vida",
main="Año 2007")
k = 1000
x = 0:50000
y = 40+40*x/(k+x)
lines(x,y,col="red",lwd=2)
y2 = recta$coefficients[1] + recta$coefficients[2] * log(x)
lines(x,y2,col="blue",lwd=2)
Os dejo pensar qué transformación de las variables originales lleva a cabo la curva de Michaelis-Menten para linealizar la relación.
J.R. Berrendero me ha pasado amablemente su código para resolver el
ejercicio con ggplot2:
# Esperanzas de vida en 2007
df <- gapminder %>%
filter(year == 2007 & continent == "Europe")
g_europa <- ggplot(df) +
geom_histogram(aes(x = lifeExp),
bins = 8,
col = 'black',
fill = "olivedrab4") +
labs(title = "Esperanza de vida en países europeos (2007)",
x = " ", y = " ")
df <- gapminder %>%
filter(year == 2007 & continent == "Africa")
g_africa <- ggplot(df) +
geom_histogram(aes(x = lifeExp),
bins = 8,
col = 'black',
fill = "olivedrab4") +
labs(title = "Esperanza de vida en países africanos (2007)",
x = " ", y = " ")
g_europa / g_africa
# Diagramas de cajas
ggplot(gapminder %>% filter(year == 1952)) +
geom_boxplot(aes(x = continent, y = lifeExp), fill='olivedrab4') +
labs(title = 'Diagramas de cajas para cada continente', y = 'Esperanza de vida', x = 'Continente')
# Diagramas de dispersión
ggplot(gapminder %>% filter(year == 2007)) +
geom_point(aes(x = gdpPercap, y = lifeExp)) +
labs(title = 'Renta per cápita y esperanza de vida en 2007', y = 'Esperanza de vida', x = 'Renta per cápita')
ggplot(gapminder %>% filter(year == 2007)) +
geom_point(aes(x = gdpPercap, y = lifeExp, col = continent)) +
labs(title = 'Renta per cápita y esperanza de vida en 2007', y = 'Esperanza de vida', x = 'Renta per cápita') +
scale_x_log10()