Cursor Euskaltel
Cursor Euskaltel
Cursor Euskaltel
1. Información y pre-requisitos 5
5. Distribuciones de probabilidad en R 59
5.1. Distribución binomial Bin(n, p) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5.2. Distribución de Poisson P ois(λ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
5.3. Distribution Exponencial Exp(λ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3
4 ÍNDICE GENERAL
11.Extras 185
11.1. El paquete caret . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 185
Capítulo 1
Información y pre-requisitos
Requisitos:
• Última versión del software R (https://www.r-project.org)
• Rstudio (https://www.rstudio.com).
Horario:
• Día 11 de marzo 10:30 a 14:30 horas
• Días 12 y 13 de marzo de 09.00 a 13.00 horas
• Día 14 de marzo 09.00 a 12.00 horas
5
6 CAPÍTULO 1. INFORMACIÓN Y PRE-REQUISITOS
Capítulo 2
2.2. Rstudio
RStudio es un IDE (Integrated Development Environment, o Entorno de Desarrollo Integrado) de código
abierto para R, que permite interactuar con R de manera muy simple.
Entre otras ventajas, Rstudio utiliza diferentes colores para las distintas clases de objetos de R, permite
autocompletar código, incluye un sistema de menús de ayuda muy completo, cuenta con un potente sistema
para la gestión, descarga y construcción de librerías, dispone de un depurador de código que detecta posibles
errores de sintaxis, es multiplataforma (existen versiones para Windows, Linux y Mac).
A la izquierda, la consola donde se ejecutan los comandos de R.
A la derecha, en la parte superior, tenemos una ventana que muestra nuestro entorno (environment) de
trabajo, en el que iremos viendo las variables y funciones que vayamos cargando, creando, etc. Obsérvese
que esta ventana tiene algunos iconos que permiten guardar el contenido de la memoria, cargar el contenido
de la memoria de una sesión de trabajo anterior, importar archivos de datos que se hayan guardado como
texto, y limpiar el contenido de la memoria.
7
8 CAPÍTULO 2. INTRODUCCIÓN AL SOFTWARE ESTADÍSTICO R
A la derecha, en la parte inferior, se muestra el contenido de nuestro directorio home donde R arranca por
defecto. Observemos que esta ventana tiene varias pestañas:
Files: archivos en el directorio actual.
Plots: en esta ventana se irán mostrando los gráficos que generemos con el programa.
Packages: permite ver qué librerías (colecciones de funciones que extienen la funcionalidad de R)
tenemos instaladas; asimismo nos permite descargar e instalar nuevas librerías.
Help: permite acceder a ayuda sobre ‘R.
Viewer: Permite acceder a contenido web local.
Ejemplo:
install.packages("DAAG")
# Varios paquetes con el comando c("< Paquete 1 >","< Paquete 2 >")
install.packages(c("DAAG","psych","gdata","foreign","Hmisc","xlsx",
"MASS","calibrate","corrplot","fields","RColorBrewer",
"ggplot2","lattice","visreg",
"car","InformationValue","ROCR"))
Desde Rstudio: Files > Click en ... > More > Set as workind directory
Ver los últimos comandos utilizados en la consola
history() # mostrar los últimos 25 comandos
history(max.show=Inf) # mostrar todos los comandos anteriores
## [1] 2 3 5 7 11 13 17
scan("ex.txt", skip = 1, nlines = 1) # only 1 line after the skipped one
## [1] 2 3 5 7
unlink("ex.data") # tidy up
Existen diferentes formatos (.txt, .csv, .xls, .xlsx, SAS, Stata, etc…)
Algunas librerías de R para importar datos:
library(gdata)
library(foreign)
O también
library(Hmisc)
mydata = spss.get("mydata.por", use.value.labels=TRUE) # SPSS
2.9. Vectores
Crear dos vectores
weight<-c(60,72,57,90,95,72)
class(weight)
## [1] "numeric"
height<-c(1.75,1.80,1.65,1.90,1.74,1.91)
Resumen de un vector
summary(weight)
Cuantiles y percentiles
quantile(weight) # por defecto cuantil 25%, 50% y 75%
Si σxy < 0m hay dependencia inversa o negativa, es decir, a grandes valores de x corresponden pequeños
valores de y.
1∑
n
Cov(x, y) = (xi − x̄)(yi − ȳ)
n i=1
cov(weight,height)
## [1] 0.6773333
El coeficiente de correlación mide la relacion lineal (positiva o negativa) entre dos variables. Formalmente es
el cociente entre la covarianza y el producto de las desviaciones típicas de ambas variables. Siendo σx y σy
las desviaciones estandar y σx y la covarianza entre x e y.
σxy
ρxy =
σx σy
cor(weight,height)
## [1] 0.437934
## [1] "character"
table(sex)
## sex
## FEMALE MALE
## 2 4
## [1] "data.frame"
str(Dat) # Ver la estructura del data.frame
3. función by o colMeans
# 'by' divide los datos en factores y realiza
# los cálculos para cada grupo
by(Dat[,3:5],sex, colMeans)
4. función aggregate
# otra opción
aggregate(Dat[,3:5], by=list(sex),mean)
2.13. VECTORES LÓGICOS 15
Lista de 2 vectores
zz <- list(x,y) # crea una lista
unlist(zz) # deshace la lista convirtiéndola en un vector concatenado
## [1] 2 3 5 2 7 1 10 15 12
Subconjunto de vectores
x[c(1,3,4)]
## [1] 2 5 2
x[-c(2,6)] # simbolo - omite los elementos
## [1] 2 5 2 7
Secuencias
seq(1,9) # ó 1:9
## [1] 1 2 3 4 5 6 7 8 9
seq(1,9,by=1)
## [1] 1 2 3 4 5 6 7 8 9
seq(1,9,by=0.5)
## [1] 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0
seq(1,9,length=20)
rep(c(2,3,5), 4)
rep(1:2,c(10,15))
## [1] 1 2 3 4 5 6 7 8 9 10 11 12
dim(x)<-c(3,4) # 3 filas y 4 columnas
X <- matrix(1:12,nrow=3,byrow=TRUE)
X
## D E F G
## A 1 4 7 10
## B 2 5 8 11
## C 3 6 9 12
colnames(X) <- month.abb[4:7]
X
## C 3 6 9 12
Concatenar filas y columnas rbind(), cbind()
Y <- matrix(0.1*(1:12),3,4)
2.16. Factores
gender<-c(rep("female",691),rep("male",692))
class(gender)
## [1] "character"
# cambiar vector a factor (por ejemplo a una categoria)
gender<- factor(gender)
levels(gender)
## female male
## 691 692
table(gender)
## gender
## female male
## 691 692
status<- c(0,3,2,1,4,5) # Crear vector numerico,
# transformarlo a niveles.
fstatus <- factor(status, levels=0:5)
levels(fstatus) <- c("student","engineer",
"unemployed","lawyer","economist","dentist")
max(a[b])
## [1] 4
sum(a[b])
## [1] 5
## [1] NA
El argumento na.rm=TRUE excluye los valores NA en el cálculo de algunos valores
sum(a,na.rm=TRUE)
## [1] 10
a <- c(1,2,3,4,NA)
is.na(a) # YES or NO
## [1] 1 2 3 4
## attr(,"na.action")
## [1] 5
## attr(,"class")
## [1] "omit"
NA en data frames:
2.19. TRABAJANDO CON DATA FRAMES 19
require(graphics)
?airquality
pairs(airquality, panel = panel.smooth, main = "airquality data")
airquality data
0 100 250 60 80 0 10 20 30
0 100
Ozone
0 200
Solar.R
5 15
Wind
90
Temp
60
9
Month
7
5
20
Day
0
0 50 150 5 10 20 5 6 7 8 9
ok <- complete.cases(airquality)
airquality[ok,]
## mpg am
## Mazda RX4 21.0 1
## Mazda RX4 Wag 21.0 1
## Datsun 710 22.8 1
## Hornet 4 Drive 21.4 0
## Hornet Sportabout 18.7 0
## Valiant 18.1 0
## Duster 360 14.3 0
## Merc 240D 24.4 0
## Merc 230 22.8 0
## Merc 280 19.2 0
## Merc 280C 17.8 0
## Merc 450SE 16.4 0
## Merc 450SL 17.3 0
## Merc 450SLC 15.2 0
## Cadillac Fleetwood 10.4 0
## Lincoln Continental 10.4 0
## Chrysler Imperial 14.7 0
## Fiat 128 32.4 1
## Honda Civic 30.4 1
## Toyota Corolla 33.9 1
## Toyota Corona 21.5 0
## Dodge Challenger 15.5 0
## AMC Javelin 15.2 0
## Camaro Z28 13.3 0
## Pontiac Firebird 19.2 0
## Fiat X1-9 27.3 1
## Porsche 914-2 26.0 1
## Lotus Europa 30.4 1
2.19. TRABAJANDO CON DATA FRAMES 21
## 25% 50%
## 352.00000 39.60853 84.20792 0.00000 2.42875 4.00000
## 75%
## 18.97500 472.00000 -341.00000
22 CAPÍTULO 2. INTRODUCCIÓN AL SOFTWARE ESTADÍSTICO R
Capítulo 3
23
24 CAPÍTULO 3. ANÁLISIS DE DATOS BÁSICO EN R
Scatterplot Example
30
Miles Per Gallon
25
20
15
10
2 3 4 5
Car Weight
Matriz scatterplot
pairs(~mpg+disp+drat+wt,data=mtcars,
main="Simple Scatterplot Matrix")
mpg 30
20
10
300
disp
100
5.0
4.0
drat
3.0
2 3 4 5
wt
-1.bb 4 6 8
Piechart o diagrama de tarta
pie(tab)
8
Ejercicio:
1. El data.frame VADeaths contiene las tasas de mortalidad por cada 1000 habitantes en Virginia (EEUU)
en 1940
Las tasas de mortalidad se miden cada 1000 habitantes por año. Se encuentran clasificadas por grupo
de edad (filas) y grupo de población (columnas). Los grupos de edad son: 50-54, 55-59, 60-64, 65-69,
70-74 y los grupos de población: Rural/Male, Rural/Female, Urban/Male and Urban/Female.
data(VADeaths)
VADeaths
• Result:
• Resultado:
Crear una tabla de conteos para cada species y realiza un gráfico descriptivo.
• Resultado:
##
## Acacia mabellae C. fraseri Acmena smithii B. myrtifolia
## 16 12 26 11
25
20
15
10
5
0
Realiza un gráfico que relacione la biomasa de la madera (wood) y el di?metro a la altura del pecho
(dbh). Utiliza tambi?n la escala logarítmica.
Acmena <- subset(rainforest, species == "Acmena smithii")
28 CAPÍTULO 3. ANÁLISIS DE DATOS BÁSICO EN R
7
6
1000
log(wood)
wood
5
500
4
3
0
dbh log(dbh)
Histogram of Acmena$dbh
8
6
Frequency
4
2
0
10 20 30 40 50
-1.bb Acmena$dbh
4. Crea un vector de n?meros enteros positivos impares the longitud 100 y calcula los valores entre 60 y
80.
Result:
## [1] 61 63 65 67 69 71 73 75 77 79
3.2. Scatterplots
library(MASS)
data("mammals")
?mammals
head(mammals)
## body brain
## Arctic fox 3.385 44.5
## Owl monkey 0.480 15.5
## Mountain beaver 1.350 8.1
## Cow 465.000 423.0
## Grey wolf 36.330 119.5
## Goat 27.660 115.0
attach(mammals)
species <- row.names(mammals)
x <- body
y <- brain
library(calibrate)
# scatterplot
30 CAPÍTULO 3. ANÁLISIS DE DATOS BÁSICO EN R
Asian elephant
brain weight in gr
3000
Human
1000
Giraffe
Horse
Okapi
Chimpanzee
Cowseal
Donkey
Gorilla
Grey
Pig
Baboon
Rhesus
Sheep
Brazilian monkey
tapir
Jaguar
Grey
Goat
Patas
Roe
Verbet
Kangaroo
Red
Genet
Owl
Slow
Rabbit
Tree
Rockfox
Arctic
Raccoon
Cat
Echidna wolf
deer
Giant
Chinchilla
RatN.A.
Guinea
Galago
Ground
Water
MoleE.
Golden monkey
armadillo
fox
hyrax−b
Yellow−bellied
monkey
loris
hyrax
Phalanger
Mountainhyrax−a
Nine−banded
African
opossum
Arctic
European
Tenrec rat
Desert
shrew
Star−nosed
Mouse
Musk
Big
Little
Lesser
brown beaver
giant
ground
pig
squirrel
opossum
hedgehog
American
shrew
brown bat marmot
armadillo
pouched
squirrel
hedgehog
mole
hamster
mole
short−tailed
bat rat
shrew
0
En escala logarítmica
plot(log(x),log(y), xlab = "log body weight in kgr", ylab = "log brain weight in gr",
main="log Body vs log Brain weight \n for 62 Species of Land Mammals")
textxy(log(x),log(y),labs=species,col = "blue",cex=0.85)
3.3. MÁS OPCIONES GRÁFICAS 31
8 Asian African
elephantelephant
Human
Giraffe
Horse
Okapi
log brain weight in gr
Chimpanzee Cow
GreyDonkey
Gorilla
seal
6
Baboon
Rhesus monkey
Sheep Pig
JaguarBrazilian tapir
Patas
Goat
Roe Grey
monkey
deer wolf
Giant armadillo
Verbet
Red
Arctic fox
fox Kangaroo
Raccoon
4
Cat
Echidna hyrax−b
Owl monkey Genet Rock Yellow−bellied marmot
Rock hyrax−a Slow
Treeloris
Rabbithyrax
Phalanger Nine−banded armadillo
Chinchilla Mountain
African
N.A. beaver
giant
opossum pouched rat
Arctic ground squirrel Guinea pig
2
Ground Galago
squirrel Water opossum
European
Mole hedgehog
rat Tenrec
Tree shrew
Desert hedgehogRat
E. American
Star−nosed mole
Goldenmole
hamster
0
Big brownMouse
Musk shrew
bat
Little brown bat
−2
−4 −2 0 2 4 6 8
Una vez realizado un plot, el comando points permite añadir nuevas observaciones.
set.seed(1234)
x <- rnorm(10,sd=5,mean=20)
y <- 2.5*x - 1.0 + rnorm(10,sd=9,mean=0)
cor(x,y)
## [1] 0.7512194
plot(x,y,xlab="Independent",ylab="Dependent",main="Random plot")
x1 <- runif(8,15,25)
y1 <- 2.5*x1 - 1.0 + runif(8,-6,6)
points(x1,y1,col=2)
32 CAPÍTULO 3. ANÁLISIS DE DATOS BÁSICO EN R
Random plot
60
50
Dependent
40
30
20
10 15 20 25
Independent
con la leyenda
set.seed(1234)
x2 <- runif(8,15,25)
y2 <- 2.5*x2 - 1.0 + runif(8,-6,6)
plot(x,y,xlab="Independent",ylab="Dependent",main="Random plot")
points(x1,y1,col=2,pch=3)
points(x2,y2,col=4,pch=5)
legend("topleft",c("Original","one","two"),col=c(1,2,4),pch=c(1,3,5))
3.3. MÁS OPCIONES GRÁFICAS 33
Random plot
60 Original
one
two
50
Dependent
40
30
20
10 15 20 25
Independent
2
2
1
1
y−label
0
0
−1
−2
0 1 2 3 4 5 6 7 0 20 40 60 80 100
third plot
5
60
30
4
Frequency
Frequency
rexp(100)
40
20
3
2
20
10
1
0
0
−4 −2 0 1 2 3 0 1 2 3 4 5 6 7 −2 −1 0 1 2
4 5 6 7 −2 0 2 4
0
−2
6
v
4
12
w
8
x
2
−2
18
y
12
−2 0 1 8 10 12 14 12 16 20
Gráfico de correlaciones
La función corrplot de la librer?a corrplot permite visualizar una matriz de correlaciones calculada me-
diante la función cor
library(corrplot)
M <- cor(d)
corrplot(M, method="circle",type="upper")
36 CAPÍTULO 3. ANÁLISIS DE DATOS BÁSICO EN R
w
u
y
1
u 0.8
0.6
v 0.4
0.2
w 0
−0.2
x −0.4
−0.6
y −0.8
−1
par(mfrow=c(2,2))
image(x,y,f)
contour(x,y,f)
contour(x,y,f,nlevels=4)
image(x,y,f,col=grey.colors(100))
contour(x,y,f,nlevels=4,add=TRUE,col="red")
3.3. MÁS OPCIONES GRÁFICAS 37
0.
0
0.
−0 .4
6
6
−0
.6
−0
0.6
−0 .8
0.
−0 .4
2
0.
−0 .6
4
−0
.2 .8
0
5
0
0.
5
−0 0
0.
2 −0
6
.6 .8
0. −0 −0
4 .4 .4
0. −0
2 .6 −0
0.
0 .8 0
6
4
0. −0
4
0. 4 −0 .2
6 0. .4 −0
2 .8
0. 0.8 −0 −0
4 .2 .6
0.2 −0 −0
.8 − .6
y
0.2
3
3
0 0 −0
−0 .4 −0
0.6 .2 .8
0.4
0.2 −0.6 0.8 0.8
0.6 −0. −
2 −0 0.8 −0.4
.6
2
0.4
2
−0.4 −0
0.2 0.8 .8
−0.2 0
−0.4 0.8
−0.6 0
1
0.6 0.4 −0.8
1
0
0
0 1 2 3 4 5 6 0 1 2 3 4 5 6
0 0
6
0 0 −0
6
−0 0. .5
0. .5 5
5
−0
−0 .5
.5 0. 0
0. 5
5
5 0 0
5
0. 0 −0
5 −0 0.5 .5
.5
0 0
−0 −0
.5 .5
4
0.5 0.5
4
0 − 0 −0
0.5 0 .5 .5
0.5
0 0
y
−0. −0.
5
3
0.5 5
0 −0.5 0.5 0
−0.5
0.5 0.5
0 0
−0.5
2
−0.5
2
−0.5 0
1
−0.5
0.5 0 −0.5
0 −0.5
0.5
0.5
0
0.5
0
0 1 2 3 4 5 6 0 1 2 3 4 5 6
o representar imágenes
library(fields)
data(lennon)
image(lennon,col=grey(seq(0,1,l=256)))
3.4. TABLAS DE CLASIFICACIÓN CRUZADA O DE CONTIGENCIA 39
1.0
0.8
0.6
0.4
0.2
0.0
## Sex
## F M
## 80 66
table(Sex,Age)
## Age
## Sex F0 F1 F2 F3
## F 10 32 19 19
## M 17 14 21 14
# or xtabs
xtabs(~Sex+Age,data=quine)
## Age
3.5. DATOS CUALITATIVOS 41
## Sex F0 F1 F2 F3
## F 10 32 19 19
## M 17 14 21 14
xtabs(~Sex+Age+Eth,data=quine)
## , , Eth = A
##
## Age
## Sex F0 F1 F2 F3
## F 5 15 9 9
## M 8 5 11 7
##
## , , Eth = N
##
## Age
## Sex F0 F1 F2 F3
## F 5 17 10 10
## M 9 9 10 7
Cálculos sobre tablas de contigencia
tapply(Days,Age,mean)
## F0 F1 F2 F3
## 14.85185 11.15217 21.05000 19.60606
tapply(Days,list(Sex,Age),mean)
## F0 F1 F2 F3
## F 18.70000 12.96875 18.42105 14.00000
## M 12.58824 7.00000 23.42857 27.21429
tapply(Days,list(Sex,Age),function(x) sqrt(var(x)/length(x)))
## F0 F1 F2 F3
## F 4.208589 2.329892 5.299959 2.940939
## M 3.768151 1.418093 3.766122 4.569582
Tabla de contigencia
xtabs(~treatment+improved)
## improved
## treatment none some marked
## placebo 29 7 7
## treated 13 7 21
42 CAPÍTULO 3. ANÁLISIS DE DATOS BÁSICO EN R
De manera gráfica,
spineplot(improved ~ treatment)
1.0
0.8
some
0.6
improved
0.4
none
0.2
0.0
placebo treated
treatment
El conjunto de datos de R, UCBAdmissionscontiene los datos agregadps de los solicitantes a universidad de
Berkeley a los seis departamentos más grandes en 1973 clasificados por sexo y admisión.
data("UCBAdmissions")
?UCBAdmissions
apply(UCBAdmissions, c(2,1), sum)
## Admit
## Gender Admitted Rejected
## Male 1198 1493
## Female 557 1278
prop.table(apply(UCBAdmissions, c(2,1), sum))
## Admit
## Gender Admitted Rejected
## Male 0.2646929 0.3298719
## Female 0.1230667 0.2823685
ftable(UCBAdmissions)
## Dept A B C D E F
## Admit Gender
## Admitted Male 512 353 120 138 53 22
## Female 89 17 202 131 94 24
## Rejected Male 313 207 205 279 138 351
## Female 19 8 391 244 299 317
Con ftable podemos presentar la información con mayor claridad
ftable(round(prop.table(UCBAdmissions), 3),
row.vars="Dept", col.vars = c("Gender", "Admit"))
3.5. DATOS CUALITATIVOS 43
Resulta más intereseante mostrar la información por género Gender y Dept combinados (dimensiones 2 y
3 del array). Nótese que las tasas de admisión por male y female son más o menos similares en todos los
departamentos, excepto en “A”, donde las tasas de las mujeres es mayor.
# prop.table(UCBAdmissions, c(2,3))
ftable(round(prop.table(UCBAdmissions, c(2,3)), 2),
row.vars="Dept", col.vars = c("Gender", "Admit"))
## Gender
## Admit Male Female
## Admitted 1198 557
## Rejected 1493 1278
Gráficamente
spineplot(margin.table(UCBAdmissions, c(3, 2)),
main = "Applications at UCB")
44 CAPÍTULO 3. ANÁLISIS DE DATOS BÁSICO EN R
Applications at UCB
Female
1.0
0.8
0.6
Gender
Male
0.4
0.2
0.0
A B C D E F
Dept
spineplot(margin.table(UCBAdmissions, c(3, 1)),
main = "Admissions at UCB")
Admissions at UCB
1.0
Rejected
0.8
0.6
Admit
Admitted
0.4
0.2
0.0
A B C D E F
Dept
Estos datos ilustran la denominada paradoja de Simpson. Este hecho ha sido analizado como un posible
caso de discriminación por sexo en las tasas de admisión en Berkeley. De los 2691 hombres que solicitaron
3.6. DATOS CUANTITATIVOS 45
se admitidos, 1198 (44.5 %) fueron admitidos, comparado con las 1835 mujeres de las cuales tan sólo 557
(30.4 %) fueron admitidas. Se podr?a por tanto concluir que los hombres tienes tasas de admisión mayores
que las mujeres. Wikipedia: Gender Bias UC Berkeley.
1000
[Marker Y]
[Marker Y]
600
600
200
200
0
0 20 60 100 0 10 30 50
[Marker X] [Marker X]
## eruptions waiting
## 1 3.600 79
## 2 1.800 54
## 3 3.333 74
## 4 2.283 62
## 5 4.533 85
## 6 2.883 55
Consideremos los datos del geyser Old Faithful en el parque nacional de Yellowstone, EEUU.
plot(faithful)
46 CAPÍTULO 3. ANÁLISIS DE DATOS BÁSICO EN R
90
80
waiting
70
60
50
eruptions
## [1] 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5
La función cut nos permite divider el rango en los intervalos que especifiquemos, con el argumento
right=FALSE, consideramos el intervalo cerrado por la derecha.
duration.cut = cut(duration, breaks, right=FALSE)
## duration.cut
## [1.5,2) [2,2.5) [2.5,3) [3,3.5) [3.5,4) [4,4.5) [4.5,5) [5,5.5)
## 51 41 5 7 30 73 61 4
Con hist podemos realizarlo de manera autom?tica:
freq <- hist(duration)
freq
3.6. DATOS CUANTITATIVOS 47
## $breaks
## [1] 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5
##
## $counts
## [1] 55 37 5 9 34 75 54 3
##
## $density
## [1] 0.40441176 0.27205882 0.03676471 0.06617647 0.25000000 0.55147059
## [7] 0.39705882 0.02205882
##
## $mids
## [1] 1.75 2.25 2.75 3.25 3.75 4.25 4.75 5.25
##
## $xname
## [1] "duration"
##
## $equidist
## [1] TRUE
##
## attr(,"class")
## [1] "histogram"
freq <- hist(duration,breaks = breaks)
Histogram of duration
60
Frequency
40
20
0
2 3 4 5
duration
hist(duration,50)
48 CAPÍTULO 3. ANÁLISIS DE DATOS BÁSICO EN R
Histogram of duration
12
10
8
Frequency
6
4
2
0
duration
Estimación de densidad construye una estimación dada una distribucion de probabilidad para una muestra
dada.
require(graphics)
d <- density(faithful$eruptions)
d
##
## Call:
## density.default(x = faithful$eruptions)
##
## Data: faithful$eruptions (272 obs.); Bandwidth 'bw' = 0.3348
##
## x y
## Min. :0.5957 Min. :0.0002262
## 1st Qu.:1.9728 1st Qu.:0.0514171
## Median :3.3500 Median :0.1447010
## Mean :3.3500 Mean :0.1813462
## 3rd Qu.:4.7272 3rd Qu.:0.3086071
## Max. :6.1043 Max. :0.4842095
plot(d)
3.6. DATOS CUANTITATIVOS 49
density.default(x = faithful$eruptions)
0.5
0.4
0.3
Density
0.2
0.1
0.0
1 2 3 4 5 6
70
60
50
Duration in minutes
h2
##
## ----------------------------
50 CAPÍTULO 3. ANÁLISIS DE DATOS BÁSICO EN R
Frecuencias relativas
duration.relfreq <- duration.freq / nrow(faithful)
tab <- cbind(duration.freq, duration.relfreq)
apply(tab,2,sum)
## duration.freq duration.relfreq
## 272 1
0.005
90
0.015
80
0.02
waiting
70
0.01
0.005
60
0.01
0.015
50
eruptions
detach("faithful")
Gráficos persp
persp(Dens2d,phi=30,theta=20,d=5,xlab="eruptions",ylab="waiting",zlab="",shade=.2,col="lightblue",expand
0.025
0.020
0.015
0.010
90
0.005 80
ng
70
iti
60
wa
2
3 50
erupti 4
ons 5
52 CAPÍTULO 3. ANÁLISIS DE DATOS BÁSICO EN R
Capítulo 4
4.1. Condicionales
Comparaciones
equal: ==
"hola" == "hola"
## [1] TRUE
"hola" == "Hola"
## [1] FALSE
1 == 2-1
## [1] TRUE
not equal: !=
a <- c(1,2,4,5)
b <- c(1,2,3,5)
a == b
## [1] TRUE TRUE TRUE FALSE TRUE TRUE FALSE TRUE TRUE TRUE
mayor/menor que o igual: >= <=
set.seed(2)
a <- rnorm(10)
53
54 CAPÍTULO 4. INTRODUCCIÓN A LA PROGRAMACIÓN BÁSICA CON R
b <- rnorm(10)
a >= b
## [1] FALSE FALSE TRUE FALSE FALSE TRUE FALSE FALSE TRUE FALSE
which
set.seed(3)
which(a>b)
## [1] 3 6 9
LETTERS
## [1] "A" "B" "C" "D" "E" "F" "G" "H" "I" "J" "K" "L" "M" "N" "O" "P" "Q"
## [18] "R" "S" "T" "U" "V" "W" "X" "Y" "Z"
which(LETTERS=="R")
## [1] 18
which.min o which.max
set.seed(4)
a <- rnorm(10)
a
## [1] 7
which.max(a)
## [1] 9
is.na
a[2] <- NA
is.na(a)
## [1] FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
which(is.na(a))
## [1] 2
## [1] 4 5 6
or: |
z = 1:6
(z > 2) & (z < 5)
## [1] 3 4
not: !
x <- c(TRUE,FALSE,0,6)
y <- c(FALSE,TRUE,FALSE,TRUE)
!x
## [1] FALSE
|| vs |
x||y
## [1] TRUE
x|y
4.3. if
if(cond1=true) { cmd1 } else { cmd2 }
if(1==0) {
print(1)
} else {
print(2)
}
## [1] 2
4.4. ifelse
ifelse(test, true_value, false_value)
x <- 1:10 # Creates sample data
ifelse(x<5 | x>8, x, 0)
## [1] 1 2 3 4 0 0 0 0 9 10
4.5.1. for
Sintaxis
for(variable in sequence) {
statements
}
for (j in 1:5)
{
print(j^2)
}
## [1] 1
## [1] 4
## [1] 9
## [1] 16
## [1] 25
## [1] 1 4 9 16 25
## [1] 14.563
Ejemplo:
x <- 1:10
z <- NULL
for(i in seq(along=x)) {
if (x[i]<5) {
z <- c(z,x[i]-1)
} else {
stop("values need to be <5")
}
}
## Error: values need to be <5
z
## [1] 0 1 2 3
4.5. LOOPS O BUCLES 57
4.5.2. while
Similar al bucle for, pero las iteraciones están controladas por una condición.
z <- 0
while(z < 5) {
z <- z + 2
print(z)
}
## [1] 2
## [1] 4
## [1] 6
58 CAPÍTULO 4. INTRODUCCIÓN A LA PROGRAMACIÓN BÁSICA CON R
Capítulo 5
Distribuciones de probabilidad en R
El paquete stats de R (que se instala por defecto al instalar R, y se carga en memoria siempre que iniciamos
sesión) implementa numerosas funciones para la realización de cálculos asociados a distintas distribuciones
de probabilidad.
Entre las utilizadas más comúnmente podemos citar:
Distribuciones discretas:
Distribución Nombre en R
Binomial binom
Poisson pois
Geométrica geom
Hipergeométrica hyper
Distribuciones contínuas:
Distribución Nombre en R
Uniforme unif
Normal norm
t-Student t
F-Fisher F
Chi-cuadrado chisq
Examinamos algunas de las operaciones básicas asociadas con las distribuciones de probabilidad.
Hay un gran número de distribuciones de probabilidad disponibles, pero sólo observamos unas pocas.
Para obtener una lista completa de las distribuciones disponibles en R puede utilizar el siguiente co-
mando:
help("Distributions")
Para cada distribución hay cuatro comandos. Los comandos para cada distribución están precedidos
de una letra para indicar la funcionalidad:
59
60 CAPÍTULO 5. DISTRIBUCIONES DE PROBABILIDAD EN R
( )
n k
f (k, n, p) = Pr(X = k) = p (1 − p)n−k , k = 0, 1, 2, ..., n
k
k ( )
∑ n
F (k; n, p) = Pr(X ≤ k) = pi (1 − p)n−i
i=0
i
0.6
0.10
P(X=k)
0.2
0.00
0.0
0 10 20 30 40 0 10 20 30 40
x x
5.1. DISTRIBUCIÓN BINOMIAL BIN (N, P ) 61
## [1] 0.1328756
Para encontrar la probabilidad de tener cuatro o menos respuestas correctas mediante intentos aleatorios,
aplicamos la función dbinom' conk=0,1,2,4‘.
prob <- NULL
for(k in 0:4){
prob <- c(prob,dbinom(k,n,p))
prob
}
prob
## [1] 0.9274445
# or simply
sum(dbinom(0:4,n,p))
## [1] 0.9274445
Alternativamente, podemos usar la función de probabilidad acumulada para la distribución binomial ‘pbi-
nom’.
pbinom(4,size=n,prob=0.2)
## [1] 0.9274445
Solución: La probabilidad de que cuatro o menos preguntas sean contestadas correctamente al azar en un
cuestionario de opción múltiple de doce preguntas es del 92,7 %.
¿Cuál es la probabilidad de que 2 ó 3 preguntas sean respondidas correctamente?
sum(dbinom(2:3,n,p))
## [1] 0.519691
Pregunta:
Supongamos que la empresa A fabricó un producto B con una probabilidad de 0,005 de ser defectuoso.
Suponga que el producto B se envía en una caja que contiene 25 B artículos. ¿Cuál es la probabilidad de
que una caja elegida al azar contenga exactamente un producto defectuoso?
62 CAPÍTULO 5. DISTRIBUCIONES DE PROBABILIDAD EN R
Solución:
( )
25
P (X = 1) = 0,0051 (1 − 0,005)(25−1)
1
p=0.005
choose(25,1)*p^1*(1-p)^(25-1)
## [1] 0.1108317
# or
dbinom(1,25,0.005)
## [1] 0.1108317
¿Cuál es la probabilidad de que una caja elegida al azar no contenga más de un artículo defectuoso?
Solución:
P (X ≤ 1) = P (X = 0) + P (X = 1)
sum(dbinom(0:1,25,p))
## [1] 0.9930519
# or
pbinom(1,25,0.005)
## [1] 0.9930519
λk e−λ
Pr(k eventos en el intervalo) =
k!
e−λ λx
P (X ≤ x | λ) = for x = 0, 1, 2, ...
x!
5.2. DISTRIBUCIÓN DE POISSON P OIS(λ) 63
1.0
0.3
0.8
0.6
0.2
P(X ≤ k)
P(X=k)
λ=1 λ=1
λ=4 λ=4
λ = 10 λ = 10
0.4
0.1
0.2
0.0
0.0
0 10 20 30 40 0 10 20 30 40
x x
Pregunta:
Supongamos que el número de plantas individuales de una especie dada que esperamos encontrar en un
cuadrado de un metro cuadrado sigue la distribución de Poisson con una media de λ = 10. Encuentra la
probabilidad de encontrar exactamente 12 dólares por persona.
Pregunta:
Si hay doce coches cruzando un puente por minuto en promedio, encuentre la probabilidad de tener diecisiete
o más coches cruzando el puente en un minuto en particular.
## [1] 0.2578387
64 CAPÍTULO 5. DISTRIBUCIONES DE PROBABILIDAD EN R
Se puede demostrar que la distribución Binomial puede ser aproximada con la función de masa de probabi-
lidad de Poisson cuando n es grande. Usando la distribución de Poisson, la media λ = np
lambda <- n*p
sum(dpois(0:3,lambda))
## [1] 0.2650259
Es importante tener en cuenta que la aproximación de Poisson a la distribución binomial funciona bien sólo
cuando n es grande y p es pequeño. En general, la aproximación funciona bien si n ≥ 20 y p ≤ 0,05, o si
n ≥ 100 y p ≤ 0,10.
f (x; λ) = λ exp(−λx)
donde λ > 0 es la tasa del evento (también conocida como parámetro de tasa, tasa de llegada, tasa de
mortalidad, tasa de fracaso, tasa de transición). La variable exponencial x ∈ [0,
inf ty]
5.3. DISTRIBUTION EXPONENCIAL EXP (λ) 65
1.0
1.5
0.8
1.0
0.6
λ = 0.5 λ = 0.5
P(X ≤ x)
λ=1 λ=1
λ = 1.5 λ = 1.5
0.4
0.5
0.2
0.0
0.0
0 1 2 3 4 5 0 1 2 3 4 5
x x
{
1 − e−λx x≥0
F (x) = Pr(X ≤ x) =
0 x<0
## [1] 0.2231302
b.
pexp(5,rate=1/10,lower.tail = FALSE)
## [1] 0.6065307
1 (x−µ)2
f (x|µ, σ 2 ) = √ e− 2σ 2 ,
2σ 2 π
dónde
σ 2 es la variación.
El proceso para estandarizar la distribución Normal consiste en transformar la variable Normal N (µ, σ)
en N (0, 1), es decir
X −µ
Z= ∼ N (0, 1)
σ
5.4. DISTRIBUTION NORMAL N (µ, σ 2 ) 67
1.0
0.8
µ = −2
σ2 = 0.5
0.8
µ=0
σ2 = 5
µ=0
σ2 = 1
0.6
µ=0
σ2 = 0.2
0.6
µ = −2
P(X ≤ x)
σ2 = 0.5
µ=0
σ2 = 5
0.4
µ=0
0.4
σ2 = 1
µ=0
σ2 = 0.2
0.2
0.2
0.0
0.0
−4 −2 0 2 4 −4 −2 0 2 4
x x
Pregunta:
X es una variable normalmente distribuida con una media de µ = 30 y una desviación estándar de σ = 4.
Encontrar
a) P (x < 40)
b) P (x > 21)
c) P (30 < x < 35)
Solucion:
a) Para x = 40, la z estandarizada es (40 − 30)/4 = 2,5 y por tanto:
pnorm(2.5) # or
## [1] 0.9937903
pnorm(40,mean=30,sd=4,lower.tail=TRUE)
## [1] 0.9937903
b) P (x > 21)
68 CAPÍTULO 5. DISTRIBUCIONES DE PROBABILIDAD EN R
pnorm(21,mean=30,sd=4,lower.tail = FALSE)
## [1] 0.9877755
c) P (30 < x < 35)
pnorm(35,mean=30,sd=4,lower.tail = TRUE)-pnorm(30,mean=30,sd=4,lower.tail = TRUE)
## [1] 0.3943502
Pregunta:
El ingreso a una determinada universidad se determina mediante un examen nacional. Los resultados de
esta prueba se distribuyen normalmente con una media de 500 y una desviación estándar de 100. Tom
quiere ser admitido en esta universidad y sabe que debe obtener mejores resultados que al menos el 70 %
de los estudiantes que tomaron el examen. Tom toma el examen y saca 585 puntos. ¿Será admitido en esta
universidad?
Solución:
N = 1000
hist(rnorm(N,500,100),20,col="grey")
abline(v=585,col=2)
hist-rnorm-03-1.bb
Histogram of rnorm(N, 500, 100)
200
150
Frequency
100
50
0
## [1] 0.8023375
Tom obtuvo una puntuación mejor que el 80.23 % de los estudiantes que tomaron el examen y será admitido
en esta universidad.
5.5. DISTRIBUCIÓN UNIFORME U (A, B) 69
1 a≤x≤b
f (x) = a − b
0 elsewhere
La función runif() puede ser usada para simular variables aleatorias uniformes independientes de n.
Por ejemplo, podemos generar 5 números aleatorios uniformes en [0, 1] como sigue:
set.seed(1234)
runif(5)
Histogram of U
15
Frequency
10
5
0
1 2 3 4 5
U
70 CAPÍTULO 5. DISTRIBUCIONES DE PROBABILIDAD EN R
0.10
0.05
0.00
0 1 2 3 4 5 6
71
72 CAPÍTULO 6. MODELOS LINEALES Y ANÁLISIS DE LA VARIANZA
y = f (x1 , x2 , ..., xp ) + ϵ,
donde f es una función desconocida y ϵ es el error en esta representación. Puesto que normalmente no
tenemos suficientes datos para intentar estimar f directamente (problema inverso), normalmente tenemos
que asumir que tiene alguna forma restringida.
La forma más simple y común es el modelo lineal (LM).
y = β0 + β1 x1 + β2 xz + ϵ,
ŷ = β̂0 + β̂1 x
set.seed(1)
n <- 50
x <- seq(1,n)
beta0 <- 15
beta1 <- 0.5
20
10
0 10 20 30 40 50
Un procedimiento matemático para encontrar la curva que mejor se ajusta a un conjunto dado de puntos
minimizando la suma de los cuadrados de los residuos de los puntos de la línea ajustada. Ilustración de los
mínimos cuadrados de ajuste
74 CAPÍTULO 6. MODELOS LINEALES Y ANÁLISIS DE LA VARIANZA
40
y = 29.36
30
y^ − y
y
y^=27.8
20
10
x=25
0 10 20 30 40 50
x
Podemos calcular directamente las cantidades de interés, es decir, la solución ordinaria de mínimos cuadrados:
∑
n
mı́n = (yi − ŷi )2
β0 ,β1
i=1
Donde
∑n
xi yi
β̂1 = ∑i=1
n
x2
y β̂0 = ȳ − β̂1 x̄
i=1 i
β̂ = (X ′ X)−1 X ′ y
0
50
1
40
Boston$medv
30
20
10
0 20 40 60 80
Boston$crim
plot(Boston[Boston$chas==0,c("crim","medv")],xlim=range(Boston$crim),ylim=range(Boston$medv),col="blue",
points(Boston[Boston$chas==1,c("crim","medv")],col="red",pch=2)
legend("topright",c("CHAS = 0", "CHAS = 1"), col=c(4,2),pch=c(1,2))
76 CAPÍTULO 6. MODELOS LINEALES Y ANÁLISIS DE LA VARIANZA
50 CHAS = 0
CHAS = 1
40
30
medv
20
10
0 20 40 60 80
crim
attach(Boston)
80
60
40
20
0
40
20
0
0 1
CHAS
50
40
30
crim
20
10
0 1
CHAS
library(ggplot2)
qplot(crim,medv,data=Boston, colour=factor(chas))
80 CAPÍTULO 6. MODELOS LINEALES Y ANÁLISIS DE LA VARIANZA
50
40
30
factor(chas)
medv
0
1
20
10
0 25 50 75
crim
qplot(crim,medv,data=Boston, colour=tax)
6.3. DEFINICIÓN DE MODELOS EN R 81
50
40
tax
30 700
600
medv
500
400
300
200
20
10
0 25 50 75
crim
library(lattice)
xyplot(medv~crim,groups=factor(chas),auto.key = TRUE)
82 CAPÍTULO 6. MODELOS LINEALES Y ANÁLISIS DE LA VARIANZA
0
1
50
40
30
medv
20
10
0 20 40 60 80
crim
xyplot(medv~crim|factor(chas),auto.key = TRUE)
6.3. DEFINICIÓN DE MODELOS EN R 83
0 20 40 60 80
0 1
50
40
30
medv
20
10
0 20 40 60 80
crim
names(Boston)
Comenzaremos usando la función lm() para encajar un modelo de regresión lineal simple, con medv como
respuesta y lstat como predictor. La sintaxis básica de lm() es lm(y~x,data), donde y es la respuesta, x
es el predictor, y los datos son el conjunto de datos en el que se guardan estas dos variables.
84 CAPÍTULO 6. MODELOS LINEALES Y ANÁLISIS DE LA VARIANZA
Si escribimos lm.fit, se obtiene información básica sobre el modelo. Para información más detallada, usamos
summary(lm.fit). Esto nos da valores de p y errores estándar para los coeficientes, así como la estadística
R2 y el estadístico para el modelo.
lm.fit
##
## Call:
## lm(formula = medv ~ lstat, data = Boston)
##
## Coefficients:
## (Intercept) lstat
## 34.55 -0.95
summary(lm.fit)
##
## Call:
## lm(formula = medv ~ lstat, data = Boston)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.168 -3.990 -1.318 2.034 24.500
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 34.55384 0.56263 61.41 <2e-16 ***
## lstat -0.95005 0.03873 -24.53 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.216 on 504 degrees of freedom
## Multiple R-squared: 0.5441, Adjusted R-squared: 0.5432
## F-statistic: 601.6 on 1 and 504 DF, p-value: < 2.2e-16
Podemos usar la función names() para averiguar qué otra información se almacena en lm.fit. Aunque
podemos extraer estas cantidades por nombre - por ejemplo lm.fit$coefficients - es más seguro usar las
funciones del extractor como coef() para acceder a ellas.
names(lm.fit)
## (Intercept) lstat
## 34.5538409 -0.9500494
lm.fit[[1]]
## (Intercept) lstat
6.3. DEFINICIÓN DE MODELOS EN R 85
## 34.5538409 -0.9500494
coef(lm.fit)
## (Intercept) lstat
## 34.5538409 -0.9500494
Para obtener un intervalo de confianza para las estimaciones del coeficiente, podemos usar el comando
confint().
confint(lm.fit, level = 0.95)
## 2.5 % 97.5 %
## (Intercept) 33.448457 35.6592247
## lstat -1.026148 -0.8739505
La función predict() puede ser usada para producir intervalos de confianza e intervalos de predicción para
la predicción de medv para un valor dado de lstat.
CI <- predict(object = lm.fit, newdata = data.frame(lstat = c(5, 10, 15)),
interval = "confidence")
CI
Por ejemplo, el intervalo de confianza del 95 % asociado con un valor “stat” de 10 es (24,474132, 25,6325627) y
el intervalo de predicción del 95 % es (12,8276263, 37,2790683). Como se esperaba, los intervalos de confianza
86 CAPÍTULO 6. MODELOS LINEALES Y ANÁLISIS DE LA VARIANZA
y predicción se centran en el mismo punto (un valor predicho de 25,0533473 para medv cuando lstat es
igual a 10), pero estos últimos son sustancialmente más amplios.
Ahora trazaremos medv y lstat junto con la línea de regresión de los mínimos cuadrados usando las funciones
plot() y abline().
plot(medv ~ lstat, data = Boston)
abline(lm.fit)
50
40
30
medv
20
10
10 20 30
lstat
# Or using ggplot2
library(ggplot2)
ggplot(data = Boston, aes(x = lstat, y = medv)) +
geom_point() +
geom_smooth(method = "lm") +
theme_bw()
6.3. DEFINICIÓN DE MODELOS EN R 87
50
40
30
medv
20
10
0 10 20 30
lstat
Hay alguna evidencia de no linealidad en la relación entre lstat y medv. Esta cuestión se discutirá más
adelante.
plot(medv ~ lstat, data = Boston,pch=15,cex=.65,col="lightgrey")
abline(lm.fit, lwd = 3, col = "red")
88 CAPÍTULO 6. MODELOS LINEALES Y ANÁLISIS DE LA VARIANZA
50
40
30
medv
20
10
10 20 30
lstat
50
40
30
medv
20
10
10 20 30
lstat
6.3. DEFINICIÓN DE MODELOS EN R 89
Opciones pch
plot(1:20, 1:20, pch = 1:20)
20
15
1:20
10
5
5 10 15 20
1:20
50
40
30
medv
20
10
0 10 20 30
lstat
# thicker line
ggplot(data = Boston, aes(x = lstat, y = medv)) +
geom_point() +
geom_smooth(method = "lm", color = "red", size = 2) +
theme_bw()
6.3. DEFINICIÓN DE MODELOS EN R 91
50
40
30
medv
20
10
0 10 20 30
lstat
4
372
373 372
373
268 268
20
3
Standardized residuals
2
10
Residuals
1
0
0
−1
−10
−2
−20
0 5 10 15 20 25 30 −3 −2 −1 0 1 2 3
372
373
268
4
Standardized residuals
Standardized residuals
1.5
215
413 375
2
1.0
0
0.5
−2
Cook's distance
0.0
Alternatively, we can compute the residuals from a linear regression fit using the residuals() function. The
function rstudent() will return the studentized residuals, and we can use this function to plot the residuals
against the fitted values.
par(mfrow = c(1, 2))
plot(predict(lm.fit), residuals(lm.fit))
plot(predict(lm.fit), rstudent(lm.fit))
6.3. DEFINICIÓN DE MODELOS EN R 93
4
20
3
2
10
residuals(lm.fit)
rstudent(lm.fit)
1
0
0
−1
−10
0 5 10 15 20 25 30 −2 0 5 10 15 20 25 30
predict(lm.fit) predict(lm.fit)
La biblioteca car tiene una función residualPlots para evaluar los residuos (calcula una prueba de curva-
tura para cada una de las parcelas añadiendo un término cuadrático y probando que la cuadrática sea cero).
Ver ?residualPlots.
library(car)
residualPlots(lm.fit)
94 CAPÍTULO 6. MODELOS LINEALES Y ANÁLISIS DE LA VARIANZA
20
20
10
10
Pearson residuals
Pearson residuals
0
0
−10
−10
10 20 30 0 5 10 15 20 25 30
Sobre la base de las gráficos de los residuos, hay alguna evidencia de no linealidad. Las estadísticas de
apalancamiento pueden ser calculadas para cualquier número de predictores usando la función hatvalues.
La función influenceIndexPlot del paquete car crea cuatro gráficos de diagnóstico que incluyen un gráfico
de los valores de sombrero.
plot(hatvalues(lm.fit))
6.3. DEFINICIÓN DE MODELOS EN R 95
0.025
hatvalues(lm.fit)
0.015
0.005
Index
which.max(hatvalues(lm.fit))
## 375
## 375
influenceIndexPlot(lm.fit, id = 5)
Diagnostic Plots
Cook's distance
375
413
−2 1 3 0.00 0.06
Studentized residuals
372
373
hat−valuesBonferroni p−value
0.6
373
0.005 0.0250.0
372
375 415
Index
96 CAPÍTULO 6. MODELOS LINEALES Y ANÁLISIS DE LA VARIANZA
##
## Call:
## lm(formula = medv ~ lstat + age, data = Boston)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.981 -3.978 -1.283 1.968 23.158
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 33.22276 0.73085 45.458 < 2e-16 ***
## lstat -1.03207 0.04819 -21.416 < 2e-16 ***
## age 0.03454 0.01223 2.826 0.00491 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.173 on 503 degrees of freedom
## Multiple R-squared: 0.5513, Adjusted R-squared: 0.5495
## F-statistic: 309 on 2 and 503 DF, p-value: < 2.2e-16
El conjunto de datos Boston contiene 13 variables, por lo que sería engorroso tener que escribirlas todas
para poder realizar una regresión utilizando todos los predictores. En su lugar, podemos utilizar la siguiente
abreviatura:
ls.fit <- lm(medv ~ ., data = Boston)
summary(ls.fit)
##
## Call:
## lm(formula = medv ~ ., data = Boston)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.595 -2.730 -0.518 1.777 26.199
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.646e+01 5.103e+00 7.144 3.28e-12 ***
## crim -1.080e-01 3.286e-02 -3.287 0.001087 **
## zn 4.642e-02 1.373e-02 3.382 0.000778 ***
## indus 2.056e-02 6.150e-02 0.334 0.738288
## chas 2.687e+00 8.616e-01 3.118 0.001925 **
## nox -1.777e+01 3.820e+00 -4.651 4.25e-06 ***
6.3. DEFINICIÓN DE MODELOS EN R 97
Podemos acceder a los componentes individuales de un objeto de resumen por nombre (escriea ?summary.lm
para ver qué hay disponible). Por lo tanto summary(lm.fit)$r.sq nos da los R2 , y summary(lm.fit)$sigma
nos da hatsigma.
Si queremos realizar una regresión usando todas las variables pero excepto una, podemos eliminarla usando
-. Por ejemplo, en la salida de regresión anterior, age tiene un alto valor p. Así que tal vez queramos hacer
una regresión excluyendo este predictor. La siguiente sintaxis resulta en una regresión usando todos los
predictores excepto age.
ls.fit1 <- lm(medv ~ . - age, data = Boston)
summary(ls.fit1)
##
## Call:
## lm(formula = medv ~ . - age, data = Boston)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.6054 -2.7313 -0.5188 1.7601 26.2243
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 36.436927 5.080119 7.172 2.72e-12 ***
## crim -0.108006 0.032832 -3.290 0.001075 **
## zn 0.046334 0.013613 3.404 0.000719 ***
## indus 0.020562 0.061433 0.335 0.737989
## chas 2.689026 0.859598 3.128 0.001863 **
## nox -17.713540 3.679308 -4.814 1.97e-06 ***
## rm 3.814394 0.408480 9.338 < 2e-16 ***
## dis -1.478612 0.190611 -7.757 5.03e-14 ***
## rad 0.305786 0.066089 4.627 4.75e-06 ***
## tax -0.012329 0.003755 -3.283 0.001099 **
## ptratio -0.952211 0.130294 -7.308 1.10e-12 ***
## black 0.009321 0.002678 3.481 0.000544 ***
## lstat -0.523852 0.047625 -10.999 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
98 CAPÍTULO 6. MODELOS LINEALES Y ANÁLISIS DE LA VARIANZA
6.3.2.1. Interacciones
Es fácil incluir términos de interacción en un modelo lineal usando la función lm(). La sintaxis lstat:black
indica a R que incluya un término de interacción entre lstat y black. La sintaxis lstat*age incluye
simultáneamente lstat,age, y el término de interacción lstat ×age como predictores; es una abreviatura
de lstat + age + lstat:age.
summary(lm(medv ~ lstat*age, data = Boston))
##
## Call:
## lm(formula = medv ~ lstat * age, data = Boston)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.806 -4.045 -1.333 2.085 27.552
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 36.0885359 1.4698355 24.553 < 2e-16 ***
## lstat -1.3921168 0.1674555 -8.313 8.78e-16 ***
## age -0.0007209 0.0198792 -0.036 0.9711
## lstat:age 0.0041560 0.0018518 2.244 0.0252 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.149 on 502 degrees of freedom
## Multiple R-squared: 0.5557, Adjusted R-squared: 0.5531
## F-statistic: 209.3 on 3 and 502 DF, p-value: < 2.2e-16
##
## Call:
## lm(formula = medv ~ lstat + I(lstat^2), data = Boston)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.2834 -3.8313 -0.5295 2.3095 25.4148
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 42.862007 0.872084 49.15 <2e-16 ***
## lstat -2.332821 0.123803 -18.84 <2e-16 ***
6.3. DEFINICIÓN DE MODELOS EN R 99
50
40
medv
30
20
10
10 20 30
-1.bb lstat
El valor p cercano a cero asociado con el término cuadrático sugiere que conduce a un modelo mejorado.
Usamos la función anova() para cuantificar aún más hasta qué punto el ajuste cuadrático es superior al
ajuste lineal.
anova(lm.fit, lm.fit2)
Aquí el Modelo 1 (lm.fit) representa el submodelo lineal que contiene sólo un predictor, lstat, mientras
que el Modelo 2 (lm.fit2) corresponde al modelo cuadrático más grande que tiene dos predictores, lstat
y lstat2 . La función anova() realiza una prueba de hipótesis comparando los dos modelos. La hipótesis
nula es que los dos modelos se ajustan a los datos igualmente bien, y la hipótesis alternativa es que el
modelo completo es superior. Aquí la estadística F es 135.1998221 y el valor p asociado es virtualmente
100 CAPÍTULO 6. MODELOS LINEALES Y ANÁLISIS DE LA VARIANZA
cero. Esto proporciona una evidencia muy clara de que el modelo que contiene los predictores lstat y lstat2
es muy superior al modelo que sólo contiene el predictor lstat. Esto no es sorprendente, ya que antes vimos
evidencia de no linealidad en la relación entre medv y lstat. Si escribimos
par(mfrow = c(2,2))
plot(lm.fit2)
Residuals vs Fitted Normal Q−Q
30
4
20
Standardized residuals
10
Residuals
2
0
0
−10
−2
−20
15 20 25 30 35 40 −3 −2 −1 0 1 2 3
268 1
4
Standardized residuals
Standardized residuals
0.5
1.5
205
215
2
1.0
0
0.5
415
−2
Cook's distance
0.0
0.5
entonces vemos que cuando el término lstat2 se incluye en el modelo, hay poco patrón discernible en los
residuos.
Para crear un ajuste cúbico, podemos incluir un predictor de la forma I(X^3). Sin embargo, este enfoque
puede empezar a ser engorroso para los polinomios de orden superior. Un mejor enfoque implica usar la
función poly() para crear el polinomio dentro de lm(). Por ejemplo, el siguiente comando produce un ajuste
polinómico de quinto orden:
lm.fit5 <- lm(medv ~ poly(lstat, 5), data = Boston)
summary(lm.fit5)
##
## Call:
## lm(formula = medv ~ poly(lstat, 5), data = Boston)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.5433 -3.1039 -0.7052 2.0844 27.1153
##
6.3. DEFINICIÓN DE MODELOS EN R 101
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 22.5328 0.2318 97.197 < 2e-16 ***
## poly(lstat, 5)1 -152.4595 5.2148 -29.236 < 2e-16 ***
## poly(lstat, 5)2 64.2272 5.2148 12.316 < 2e-16 ***
## poly(lstat, 5)3 -27.0511 5.2148 -5.187 3.10e-07 ***
## poly(lstat, 5)4 25.4517 5.2148 4.881 1.42e-06 ***
## poly(lstat, 5)5 -19.2524 5.2148 -3.692 0.000247 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.215 on 500 degrees of freedom
## Multiple R-squared: 0.6817, Adjusted R-squared: 0.6785
## F-statistic: 214.2 on 5 and 500 DF, p-value: < 2.2e-16
Esto sugiere que incluir términos polinómicos adicionales, hasta el quinto orden, conduce a una mejora
en el ajuste del modelo. Sin embargo, una investigación adicional de los datos revela que ningún término
polinómico más allá del quinto orden tiene valores p significativos en un ajuste de regresión.
library(ggplot2)
ggplot(data = Boston, aes(x = lstat, y = medv)) +
geom_point() +
theme_bw() +
stat_smooth(method = "lm", formula = y ~ poly(x, 5))
50
40
30
medv
20
10
0
0 10 20 30
lstat
102 CAPÍTULO 6. MODELOS LINEALES Y ANÁLISIS DE LA VARIANZA
Por supuesto, no estamos restringidos a usar transformaciones polinómicas de los predictores. Aquí probamos
una transformación logarítmica.
summary(lm(medv ~ log(rm), data = Boston))
##
## Call:
## lm(formula = medv ~ log(rm), data = Boston)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.487 -2.875 -0.104 2.837 39.816
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -76.488 5.028 -15.21 <2e-16 ***
## log(rm) 54.055 2.739 19.73 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.915 on 504 degrees of freedom
## Multiple R-squared: 0.4358, Adjusted R-squared: 0.4347
## F-statistic: 389.3 on 1 and 504 DF, p-value: < 2.2e-16
ggplot(data = Boston, aes(x = log(rm), y = medv)) +
geom_point() +
theme_bw() +
stat_smooth(method = "lm")
40
medv
20
8.7)
ventas <- c(22.1, 10.4, 9.3, 18.5, 12.9, 7.2, 11.8, 13.2, 4.8, 10.6, 8.6, 17.4,
9.2, 9.7, 19, 22.4, 12.5, 24.4, 11.3, 14.6, 18, 12.5, 5.6, 15.5, 9.7, 12,
15, 15.9, 18.9, 10.5, 21.4, 11.9, 9.6, 17.4, 9.5, 12.8, 25.4, 14.7, 10.1,
21.5, 16.6, 17.1, 20.7, 12.9, 8.5, 14.9, 10.6, 23.2, 14.8, 9.7, 11.4, 10.7,
22.6, 21.2, 20.2, 23.7, 5.5, 13.2, 23.8, 18.4, 8.1, 24.2, 15.7, 14, 18,
9.3, 9.5, 13.4, 18.9, 22.3, 18.3, 12.4, 8.8, 11, 17, 8.7, 6.9, 14.2, 5.3,
11, 11.8, 12.3, 11.3, 13.6, 21.7, 15.2, 12, 16, 12.9, 16.7, 11.2, 7.3, 19.4,
22.2, 11.5, 16.9, 11.7, 15.5, 25.4, 17.2, 11.7, 23.8, 14.8, 14.7, 20.7,
19.2, 7.2, 8.7, 5.3, 19.8, 13.4, 21.8, 14.1, 15.9, 14.6, 12.6, 12.2, 9.4,
15.9, 6.6, 15.5, 7, 11.6, 15.2, 19.7, 10.6, 6.6, 8.8, 24.7, 9.7, 1.6, 12.7,
5.7, 19.6, 10.8, 11.6, 9.5, 20.8, 9.6, 20.7, 10.9, 19.2, 20.1, 10.4, 11.4,
10.3, 13.2, 25.4, 10.9, 10.1, 16.1, 11.6, 16.6, 19, 15.6, 3.2, 15.3, 10.1,
7.3, 12.9, 14.4, 13.3, 14.9, 18, 11.9, 11.9, 8, 12.2, 17.1, 15, 8.4, 14.5,
7.6, 11.7, 11.5, 27, 20.2, 11.7, 11.8, 12.6, 10.5, 12.2, 8.7, 26.2, 17.6,
22.6, 10.3, 17.3, 15.9, 6.7, 10.8, 9.9, 5.9, 19.6, 17.3, 7.6, 9.7, 12.8,
25.5, 13.4)
El modelo lineal múltiple que se obtiene empleando las variables tv, radio y periodico como predictores
de ventas es el siguiente:
modelo <- lm(ventas ~ tv + radio + periodico, data = datos)
summary(modelo)
##
## Call:
## lm(formula = ventas ~ tv + radio + periodico, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.8277 -0.8908 0.2418 1.1893 2.8292
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.938889 0.311908 9.422 <2e-16 ***
## tv 0.045765 0.001395 32.809 <2e-16 ***
## radio 0.188530 0.008611 21.893 <2e-16 ***
## periodico -0.001037 0.005871 -0.177 0.86
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.686 on 196 degrees of freedom
## Multiple R-squared: 0.8972, Adjusted R-squared: 0.8956
## F-statistic: 570.3 on 3 and 196 DF, p-value: < 2.2e-16
De acuerdo al p-value obtenido para el coeficiente parcial de regresión de periodico, esta variable no
contribuye de forma significativa al modelo. Como resultado de este análisis se concluye que las variables tv
y radio están asociadas con la cantidad de ventas.
##
## Call:
## lm(formula = ventas ~ tv + radio, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.7977 -0.8752 0.2422 1.1708 2.8328
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.92110 0.29449 9.919 <2e-16 ***
## tv 0.04575 0.00139 32.909 <2e-16 ***
## radio 0.18799 0.00804 23.382 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.681 on 197 degrees of freedom
## Multiple R-squared: 0.8972, Adjusted R-squared: 0.8962
## F-statistic: 859.6 on 2 and 197 DF, p-value: < 2.2e-16
25
20
ventas
15
10
40
30
5
radio
20
50 100 10
150
tv 200 250 0
El modelo lineal a partir del cual se han obtenido las conclusiones asume que el efecto sobre las ventas debido
a un incremento en el presupuesto de uno de los medios de comunicación es independiente del presupuesto
gastado en los otros. Por ejemplo, el modelo lineal considera que el efecto promedio sobre las ventas debido
a aumentar en una unidad el presupuesto de anuncios en TV es siempre de 0.04575, independientemente de
la cantidad invertida en anuncios por radio. Sin embargo, la representación gráfica muestra que el modelo
tiende a sobrevalorar las ventas cuando el presupuesto es muy alto en uno de los medios pero muy bajo en el
otro. Por contra, los valores de ventas predichos por el modelo están por debajo de las ventas reales cuando el
presupuesto está repartido de forma equitativa entre ambos medios. Este comportamiento sugiere que existe
interacción entre los predictores, por lo que el efecto de cada uno de ellos sobre la variable respuesta depende
en cierta medida del valor que tome el otro predictor.
Tal y como se ha definido previamente, un modelo lineal con dos predictores sigue la ecuación:
106 CAPÍTULO 6. MODELOS LINEALES Y ANÁLISIS DE LA VARIANZA
Y = β0 + β1 X1 + β2 X2 + ϵ
De acuerdo a esta definición, el incremento de una unidad en el predictor X1 produce un incremento promedio
de la variable Y de β1 . Modificaciones en el predictor X2 no alteran este hecho, y lo mismo ocurre con X2
respecto a X1 . Para que el modelo pueda contemplar la interacción se introduce un tercer predictor, que se
construye con el producto de los predictores X1 y X2 .
Y = β0 + β1 X1 + β2 X2 + β3 X1 X2 + ϵ
Y = β0 + (β1 + β3 X2 )X1 + β2 X2 + ϵ
El efecto de X1 sobre Y ya no es constante, sino que depende del valor que tome X2 .
En R se puede introducir interacción entre predictores de dos formas, indicando los predictores individuales
y entre cuales se quiere evaluar la interacción, o bien de forma directa.
modelo_interaccion <- lm(formula = ventas ~ tv + radio + tv:radio, data = datos)
summary(modelo_interaccion)
##
## Call:
## lm(formula = ventas ~ tv + radio + tv:radio, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.3366 -0.4028 0.1831 0.5948 1.5246
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.750e+00 2.479e-01 27.233 <2e-16 ***
## tv 1.910e-02 1.504e-03 12.699 <2e-16 ***
## radio 2.886e-02 8.905e-03 3.241 0.0014 **
## tv:radio 1.086e-03 5.242e-05 20.727 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9435 on 196 degrees of freedom
## Multiple R-squared: 0.9678, Adjusted R-squared: 0.9673
## F-statistic: 1963 on 3 and 196 DF, p-value: < 2.2e-16
# O tambien
# lm(formula = ventas ~ tv * radio, data = datos) es equivalente.
6.3. DEFINICIÓN DE MODELOS EN R 107
25
20
ventas
15
40
10 30
radio
20
50 100 10
150
tv 200 250 0
Los resultados muestran una evidencia clara de que la interacción tv:radio es significativa y de que el modelo
que incorpora la interacción (Adjusted R-squared = 0.9673) es superior al modelo que solo contemplaba el
efecto de los predictores por separado (Adjusted R-squared = 0.8956).
Se puede emplear un ANOVA para realizar un test de hipótesis y obtener un p-value que evalúe la hipótesis
nula de que ambos modelos se ajustan a los datos igual de bien.
anova(modelo, modelo_interaccion)
Regresión logística
Una regresión logística se utiliza típicamente cuando hay una variable de resultado dicotómica (como ganar
o perder), y una variable predictiva continua que está relacionada con la probabilidad o “odds” de la variable
de resultado. También se puede utilizar con predictores categóricos y con múltiples predictores.
Si usamos una regresión lineal para modelar una variable dicotómica (como Y ), es posible que el modelo
resultante no restrinja los Y pronosticados dentro de 0 y 1. Además, otros supuestos de regresión lineal como
la normalidad del error pueden ser violados. Así que en vez de eso, modelamos las probabilidades del evento
de $log (logit) o logit, donde, p es la probabilidad del evento.
pi
zi = ln( ) = β0 + β1 x1 + ... + βp xp
1 − pi
La ecuación anterior puede ser modelada usando glm() mediante el argumento family="binomial. Pero
estamos más interesados en la probabilidad del evento que en las probabilidades logarítmicas del evento. Por
lo tanto, los valores pronosticados del modelo anterior, es decir, las probabilidades logarítmicas del evento,
pueden convertirse a probabilidad de evento de la siguiente manera:
1
pi = 1 −
1 + exp(zi )
109
110 CAPÍTULO 7. REGRESIÓN LOGÍSTICA
1.0
0.8
0.6
p
0.4
0.2
0.0
−6 −4 −2 0 2 4 6
log(p/1−p)
El objetivo de este tipo de análisis es minimizar el riesgo y maximizar el beneficio del banco.
Para minimizar las pérdidas desde la perspectiva del banco, el banco necesita una regla de decisión sobre
a quién dar la aprobación del préstamo y a quién no. Los perfiles demográficos y socioeconómicos de un
solicitante son considerados por los administradores de préstamos antes de que se tome una decisión sobre
su solicitud de préstamo.
Los datos German Credit Data contiene datos sobre 20 variables y la clasificación de si un solicitante es
considerado un Bueno o un Malo riesgo de crédito para 1000 solicitantes de préstamos.
Se espera que un modelo predictivo elaborado a partir de estos datos sirva de guía al gerente del banco para
tomar una decisión sobre la aprobación de un préstamo a un posible solicitante en función de su perfil.
Ver descripción (en inglés)
# German Credit Data
gcreditdata <- read.table("https://archive.ics.uci.edu/ml/machine-learning-databases/statlog/german/germ
colnames(gcreditdata)<-c("account.status","months",
7.1. EJEMPLO: DATOS DE CRÉDITO EN ALEMANIA (CREDIT SCORING) 111
"credit.history","purpose","credit.amount",
"savings","employment","installment.rate","personal.status",
"guarantors","residence","property","age","other.installments",
"housing","credit.cards","job","dependents","phone","foreign.worker","credit.rating")
head(gcreditdata)
##
## 1 2
## 700 300
gcreditdata$credit.rating <- ifelse(gcreditdata$credit.rating==1,1,0)
table(gcreditdata$credit.rating) # Good = 1 / Bad = 0
##
## 0 1
## 300 700
levels(gcreditdata$account.status) <- c("<0DM","<200DM",">200DM","NoStatus")
levels(gcreditdata$credit.history) <- c("No","Allpaid","Allpaidtillnow","Delayinpaying","Critical")
levels(gcreditdata$purpose) <- c("car(new)","car(used)","furniture/equipment","radio/television","domest
"repairs","education","vacation-doesnotexist?","retraining","business","others")
112 CAPÍTULO 7. REGRESIÓN LOGÍSTICA
##
## Call:
## glm(formula = credit.rating ~ ., family = binomial, data = gcreditdata.train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.8353 -0.6178 0.3415 0.6729 2.0514
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.259e-01 1.371e+00 0.384 0.701257
## account.status<200DM 3.532e-01 2.719e-01 1.299 0.193876
## account.status>200DM 8.198e-01 4.626e-01 1.772 0.076353 .
## account.statusNoStatus 1.487e+00 2.782e-01 5.344 9.08e-08 ***
## months -3.687e-02 1.156e-02 -3.189 0.001428 **
## credit.historyAllpaid -1.347e-01 6.539e-01 -0.206 0.836779
## credit.historyAllpaidtillnow 5.296e-01 4.975e-01 1.065 0.287038
## credit.historyDelayinpaying 7.156e-01 5.575e-01 1.284 0.199266
## credit.historyCritical 1.439e+00 5.073e-01 2.837 0.004559 **
## purposecar(used) 1.838e+00 4.856e-01 3.786 0.000153 ***
## purposefurniture/equipment 1.095e+00 9.231e-01 1.186 0.235728
## purposeradio/television 8.879e-01 3.289e-01 2.700 0.006944 **
## purposedomesticappliances 9.928e-01 3.107e-01 3.196 0.001395 **
## purposerepairs 5.399e-01 9.006e-01 0.599 0.548881
## purposeeducation -4.889e-01 6.587e-01 -0.742 0.457933
## purposevacation-doesnotexist? -1.800e-01 4.648e-01 -0.387 0.698616
## purposeretraining 2.093e+00 1.266e+00 1.653 0.098417 .
## purposebusiness 5.772e-01 3.985e-01 1.448 0.147489
## credit.amount -1.107e-04 5.456e-05 -2.028 0.042526 *
## savingsA62 5.219e-01 3.548e-01 1.471 0.141326
## savingsA63 8.722e-01 5.710e-01 1.528 0.126623
## savingsA64 1.746e+00 6.533e-01 2.673 0.007514 **
## savingsA65 1.354e+00 3.350e-01 4.042 5.30e-05 ***
7.1. EJEMPLO: DATOS DE CRÉDITO EN ALEMANIA (CREDIT SCORING) 113
Variables significativas:
sig.var<- summary(m1)$coeff[-1,4] <0.01
names(sig.var)[sig.var == TRUE]
##
## ytest 0 1
## 0 48 51
## 1 25 176
error1<- sum(result1[1,2], result1[2,1])/sum(result1)
error1
## [1] 0.2533333
0.6
0.4
0.2
0.0
## AUC: 0.7699382 n
##
## 0 1
## 24720 7841
Claramente, existe un sesgo de clase, una condición observada cuando la proporción de eventos es mucho
menor que la proporción de no eventos. Por lo tanto, debemos muestrear las observaciones en proporciones
aproximadamente iguales para obtener mejores modelos.
Crear Muestras de Entrenamiento y Pruebas
Una manera de abordar el problema del sesgo de clase es muestrear los 0 y 1 para los trainingData (muestra
de entrenamiento) en proporciones iguales. Al hacerlo, pondremos el resto de los inputData no incluidos para
la formación en testData (muestra de validación). Como resultado, el tamaño de la muestra de entrenamiento
será menor que el de la validación, lo que está bien, porque hay un gran número de observaciones (>10K).
# Create Training Data
input_ones <- inputData[which(inputData$ABOVE50K == 1), ] # all 1's
input_zeros <- inputData[which(inputData$ABOVE50K == 0), ] # all 0's
testData <- rbind(test_ones, test_zeros) # row bind the 1's and 0's
Cuando usamos la función de predicción en este modelo, predecirá las probabilidades de la variable Y . Para
convertirlo en una probabilidad de predicción que esté entre 0 y 1, usamos el plogis().
El puntaje de probabilidad de la predicción de corte por defecto es de 0,5 o la proporción de 1’s y 0’s en
los datos de entrenamiento. Pero a veces, afinar el corte de probabilidad puede mejorar la precisión tanto
en las muestras de desarrollo como en las de validación. La función InformationValue::optimalCutoff
proporciona formas de encontrar el punto de corte óptimo para mejorar la predicción de 1, 0, 1 y 0 y reducir
el error de clasificación. Permite calcular la puntuación óptima que minimiza el error de clasificación para el
modelo anterior.
library(InformationValue)
optCutOff <- optimalCutoff(testData$ABOVE50K, predicted)[1]
optCutOff
## [1] 0.89
Error de clasificación
El error de clasificación errónea es el desajuste porcentual de los valores predefinidos frente a los reales,
independientemente de que sean 1 o 0. Cuanto menor sea el error de clasificación, mejor será su modelo.
misClassError(testData$ABOVE50K, predicted, threshold = optCutOff)
## [1] 0.0892
ROC
ROC Curve
1.00
0.75
Sensitivity (TPR)
0.50
AUROC: 0.8915
0.25
0.00
## [1] 0.3442414
specificity(testData$ABOVE50K, predicted, threshold = optCutOff)
## [1] 0.9800853
Los números anteriores se calculan a partir de la muestra de validación que no se utilizó para la formación
del modelo. Así que, una tasa de detección de la verdad de 34.42 % en los datos de prueba es bueno.
Matriz de Confusión Las columnas son reales, mientras que las filas son predicciones.
confusionMatrix(testData$ABOVE50K, predicted, threshold = optCutOff)
## 0 1
## 0 18849 1543
## 1 383 810
118 CAPÍTULO 7. REGRESIÓN LOGÍSTICA
Questions:
Ajustar un modelo logístico con pclass como variable explicativa. ¿Cuál es la interpretación del modelo
ajustado?
Encontrar el mejor modelo de regresión logística posible basado en todas las variables disponibles.
model <- glm(Survived ~.,family=binomial(link='logit'),data=train)
summary(model)
##
## Call:
## glm(formula = Survived ~ ., family = binomial(link = "logit"),
## data = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.6064 -0.5954 -0.4254 0.6220 2.4165
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.137627 0.594998 8.635 < 2e-16 ***
## Pclass -1.087156 0.151168 -7.192 6.40e-13 ***
## Sexmale -2.756819 0.212026 -13.002 < 2e-16 ***
## Age -0.037267 0.008195 -4.547 5.43e-06 ***
7.3. EJEMPLO: DATOS DE LOS SUPERVIVIENTES DEL TITANIC 119
##
## Call:
## glm(formula = Survived ~ as.factor(Pclass), family = binomial,
## data = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.3787 -0.7515 -0.7515 0.9887 1.6747
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.4616 0.1474 3.131 0.00174 **
## as.factor(Pclass)2 -0.5455 0.2138 -2.551 0.01074 *
120 CAPÍTULO 7. REGRESIÓN LOGÍSTICA
##
## Call:
## glm(formula = Survived ~ Pclass + Sex + Age + SibSp, family = binomial(logit),
## data = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.6595 -0.6125 -0.4247 0.6149 2.4302
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.05604 0.50130 10.086 < 2e-16 ***
## Pclass -1.14391 0.12585 -9.089 < 2e-16 ***
## Sexmale -2.75564 0.20471 -13.461 < 2e-16 ***
## Age -0.03725 0.00812 -4.588 4.48e-06 ***
## SibSp -0.33075 0.10892 -3.037 0.00239 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
7.3. EJEMPLO: DATOS DE LOS SUPERVIVIENTES DEL TITANIC 121
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1065.39 on 799 degrees of freedom
## Residual deviance: 713.43 on 795 degrees of freedom
## AIC: 723.43
##
## Number of Fisher Scoring iterations: 5
anova(mod2,test="Chisq")
##
## Call:
## glm(formula = Survived ~ Pclass + Sex + Pclass:Sex + Age + SibSp,
## family = binomial(logit), data = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.1993 -0.6265 -0.4770 0.4485 2.3093
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 7.606411 0.960804 7.917 2.44e-15 ***
## Pclass -2.108360 0.316024 -6.672 2.53e-11 ***
## Sexmale -5.887480 0.920417 -6.397 1.59e-10 ***
## Age -0.038063 0.008498 -4.479 7.50e-06 ***
## SibSp -0.310269 0.109370 -2.837 0.004556 **
## Pclass:Sexmale 1.254202 0.338241 3.708 0.000209 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1065.39 on 799 degrees of freedom
122 CAPÍTULO 7. REGRESIÓN LOGÍSTICA
En los pasos anteriores, evaluamos brevemente el ajuste del modelo, ahora nos gustaría ver cómo lo está ha-
ciendo el modelo al predecir y en un nuevo conjunto de datos. Estableciendo el parámetro type='response',
R producirá probabilidades en la forma de P (y = 1|X). Nuestro límite de decisión será de 0,5.
Si P (y = 1|X) > 0,5 entonces y = 1 en caso contrario y = 0. Tener en cuenta que para algunas aplicaciones
diferentes límites de decisión podría ser una mejor opción.
fitted.results <- predict(mod3,newdata=test,type='response')
fitted.results <- ifelse(fitted.results > 0.5,1,0)
La precisión de 0.8075 en el conjunto de test es un buen resultado. Sin embargo, hay que tener en cuenta que
este resultado depende en cierta medida de la división manual de los datos que hice anteriormente, por lo
tanto, si deseamos una puntuación más precisa, sería mejor realizar algún tipo de validación cruzada, como
la validación cruzada k-fold.
1.0
0.8
True positive rate
0.6
0.4
0.2
0.0
## [1] 0.8543155
124 CAPÍTULO 7. REGRESIÓN LOGÍSTICA
Capítulo 8
8.1. Suavizado
Suavizado Scatterplot
Consiste en resaltar la tendencia subyacente en los datos. La tendencia subyacente sería una función como:
f (x) = E(y|x)
yi = f (xi ) + ϵi , E(ϵi ) = 0,
en cuyo caso el problema se denomina a menudo como regresión no-paramétrica, donde f (·) es una función
suave no conocida que se estima a partir de los datos (xi , yi )
Existen varios métodos para suavizar una gráfica de dispersión, incluyendo splines, regresión kernel, medias
móviles, loess o su versión ponderada lowess.
El conjunto de datos lidar consta de 221 observaciones de un experimento de detección y alcance de luz
(LIDAR).
library(SemiPar)
#library(car)
data(lidar)
plot(logratio~range,data=lidar)
125
126 CAPÍTULO 8. MODELOS ADITIVOS GENERALIZADOS
0.0
−0.2
−0.4
logratio
−0.6
−0.8
range
y = β0 + β1 X1 + β2 X2 + βp X p + ϵ ϵ ∼ N (0, σ 2 )
library(lattice)
cols. <- rainbow(5)
range.seq <- seq(min(lidar$range),max(lidar$range),l=200)
plot(logratio~range,data=lidar,cex=.5,pch=19,col="grey",xlim=c(380,725),ylim=c(-0.90,0.05))
for(i in 1:5){
lines(range.seq,predict(lm(logratio~poly(range,i),data=lidar),data.frame(range=range.seq)),lwd=2,col=c
}
legend("bottomleft",c("p = 1","p = 2","p = 3","p = 4","p = 5"),col=cols.,lty=1:5,lwd=4)
8.1. SUAVIZADO 127
0.0
−0.2
logratio
−0.4
p=1
−0.6
p=2
p=3
−0.8
p=4
p=5
range
Podemos tratar estas relaciones no lineales de manera efectiva a través de técnicas más flexibles.
Regresión local:
LOESS y LOWESS son dos métodos de regresión no-paramétrica muy relacionados que combinan múltiples
modelos de regresión en un modo basado en k-neast-neighbor.
attach(lidar)
span <- 1/c(1,2,3,4,5)
plot(logratio~range,data=lidar,cex=.5,pch=19,col="grey",xlim=c(380,725),ylim=c(-0.90,0.05))
cols <- rainbow(length(span))
for(i in 1:length(span)){
lines(loess.smooth(range,logratio,span=span[i], degree=2),col=cols.[i],lwd=2.5,lty=i)
}
legend("bottomleft",paste("span = ",round(span,3)),col=cols.,lwd=2.5,lty=1:6)
128 CAPÍTULO 8. MODELOS ADITIVOS GENERALIZADOS
0.0
−0.2
logratio
−0.4
span = 1
−0.6
span = 0.5
span = 0.333
−0.8
span = 0.25
span = 0.2
range
∑n ( )
x −x
j=1yi K i b j
ŷi = ∑ ( ) ,
n xi −xj
j=1 K b
donde b es un parámetro de ancho de banda (bandwidth), y K una función del núcleo/kernel, como en la
estimación de densidad. Hay diferentes opciones para el kernel K (por ejemplo: Gaussian), el parámetro de
elección crítico es el ancho de banda b.
En R la función ksmooth en library(stats) permite el suavizado tipo kernel (otras opciones son locpoly
en library(KernSmooth))
attach(lidar)
b <- c(10,20,50,120,1000)
plot(logratio~range,data=lidar,cex=.5,pch=19,col="grey",xlim=c(380,725),ylim=c(-0.90,0.05))
cols. <- rainbow(length(b))
for(i in 1:length(b)){
lines(ksmooth(range,logratio,x.points=range,"normal", bandwidth=b[i]),col=cols.[i],lwd=2.5,lty=i)
}
legend("bottomleft",paste("bandwidth = ",b),col=cols.,lwd=2.5,lty=1:5)
8.2. GAMS 129
0.0
−0.2
logratio
−0.4
bandwidth = 10
−0.6
bandwidth = 20
bandwidth = 50
−0.8
bandwidth = 120
bandwidth = 1000
range
8.2. GAMs
Los GAMs (del inglés generalized additive models) son una generalización de los GLMs para incorporar
formas no lineales de los predictores (plines, Polinomios, o funciones Step, etc…).
y ∼ ExpoFam(µ, σ 2 , ...)
µ = E[y]
Ahora usamos las funciones suaves f (·) de nuestras variables predictoras, que pueden tomar formas más
flexibles. Se supone que los valores observados son de alguna distribución de la familia exponencial, y µ está
relacionado con los predictores a través de una función enlace (o link).
En este curso se estudiará la implementación de GAMs en R usando el paquete mgcv y funciones de base
tipo spline.
Cargaremos la librería mgcv, y la función gam:
La función gam
La sintaxis es:
library(mgcv)
fit.gam <- gam(y ~ s(x))
En general, la sintaxis es
gam(formula,method="...",select="...",family=gaussian())
Los argumentos principales para esta función se deben especificar como una formula:
La primera opción es la base utilizada para representar los términos suaves s(x) (ver ?s o
?smooth.terms). El tipo de función base se puede modificar utilizando bs dentro de s(x,bs="ps"):
130 CAPÍTULO 8. MODELOS ADITIVOS GENERALIZADOS
m El orden de la penalización.
k la dimensión de la base utilizada para representar el término de suavizado. k debe ser mayor que m).
by una variable numérica o factorial de la misma dimensión que cada covariable.
sp cualquier parámetro de suavizado suministrado para el término de suavizado.
etc …
summary(fit.gam)
library(mgcv)
gam1 <- gam(wage~s(age,k=6)+s(year,k=6)+education,data = Wage)
summary(gam1)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## wage ~ s(age, k = 6) + s(year, k = 6) + education
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 85.456 2.153 39.692 < 2e-16 ***
## education2. HS Grad 10.978 2.428 4.521 6.4e-06 ***
## education3. Some College 23.530 2.558 9.197 < 2e-16 ***
## education4. College Grad 38.159 2.543 15.005 < 2e-16 ***
## education5. Advanced Degree 62.559 2.760 22.668 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(age) 4.435 4.840 46.69 < 2e-16 ***
## s(year) 1.101 1.195 11.11 0.000379 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.29 Deviance explained = 29.2%
## GCV = 1240.7 Scale est. = 1236.3 n = 3000
par(mfrow=c(1,2)) #to partition the Plotting Window
plot(gam1,se = TRUE)
8.2. GAMS 131
10
10
0
0
−10
−10
s(age,4.44)
s(year,1.1)
−20
−20
−30
−30
−40
−40
20 40 60 80 2003 2005 2007 2009
age year
s(year,1.1)
−10
−10
−40
−40
age year
education
Partial for education
40
0
education
300 300
250 250
f(year)
f(age) 200 200
150 150
100 100
50 50
0 0
age year
f(education)
300
250
200
150
100
50
education
El procedimiento de ajuste es el mismo que se muestra en las secciones anteriores. Simplemente los efectos
paramétricos están incluidos en la parte de efectos fijos, X. Un caso interesado es cuando la parte paramétrica
incluye una variable factorial con dos o más niveles. Como vimos en el caso lineal, podemos considerar
diferentes situaciones: términos suaves paralelos (efecto aditivo) o términos suaves no paralelos (efecto de
interacción). Además, podemos considerar la misma cantidad de alisamiento o una cantidad diferente.
Datos de cebollas
Para ilustrar un caso simple de un modelo semi-paramétrico, consideramos los data(onions) en la
library(SemiPar). Los datos consisten en 84 observaciones de un experimento que involucra la producción
de cebollas en dos localidades del sur de Australia. Las variables son:
dens densidad de área de las plantas (plantas por metro cuadrado)
yield rendimiento de la cebolla (gramos por planta).
location localización: 0=Purnong Landing, 1=Virginia.
La siguiente figura muestra que para una mayor densidad de plantas el rendimiento de Purnong Landing de
onios es mayor.
library(SemiPar)
data(onions)
attach(onions)
Purnong Landing
Virginia
5.0
log(yield)
4.5
4.0
3.5
50 100 150
dens
El modelo lineal es:
donde
{
0 si la i-ésima observación es de Purnong Landing
locationij =
1 si la i-ésima observación es de Virginia
La figura sugiere un término no lineal para dens, por lo tanto, un modelo mejor sería:
summary(fit1)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## log(yield) ~ L + s(dens, k = 20, m = 2, bs = "ps")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.85011 0.01688 287.39 <2e-16 ***
## LVirginia -0.33284 0.02409 -13.82 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(dens) 4.568 19 72.76 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.946 Deviance explained = 94.9%
## -REML = -54.242 Scale est. = 0.011737 n = 84
La siguiente figura muestra las curvas ajustadas para cada ubicación con el modelo fit1, asumimos que
ambas curvas son paralelas.
L.P <- rep(levels(L)[1],100)
L.V <- rep(levels(L)[2],100)
dens.g <- seq(min(dens),max(dens),l=100)
fit1.P <- predict(fit1,newdata=data.frame(L=L.P,dens=dens.g),se.fit=TRUE)
fit1.V <- predict(fit1,newdata=data.frame(L=L.V,dens=dens.g),se.fit=TRUE)
plot(dens,log(yield),col=points.cols[location+1],pch=16,cex=.55)
lines(dens.g,fit1.P$fit,col=2,lwd=2)
lines(dens.g,fit1.P$fit+1.96*fit1.P$se.fit,col=2,lwd=2,lty=2)
lines(dens.g,fit1.P$fit-1.96*fit1.P$se.fit,col=2,lwd=2,lty=2)
lines(dens.g,fit1.V$fit,col=4,lwd=2)
lines(dens.g,fit1.V$fit+1.96*fit1.V$se.fit,col=4,lwd=2,lty=2)
lines(dens.g,fit1.V$fit-1.96*fit1.V$se.fit,col=4,lwd=2,lty=2)
8.3. MODELOS SEMI-PARAMÉTRICOS 135
5.5
5.0
log(yield)
4.5
4.0
3.5
50 100 150
dens
Asumir curvas paralelas para ambas localizaciones implica que la disminución del rendimiento a medida que
aumenta la densidad es la misma en ambas localizaciones. En lugar de ajustar un modelo de efectos aditivos,
podemos ajustar un modelo de efectos con interacción, como por ejemplo:
donde
{
0 si la i-ésima observación es de Purnong Landing
L(j) =
1 si la i-ésima observación es de Virginia
##
## Family: gaussian
## Link function: identity
##
## Formula:
## log(yield) ~ L + s(dens, k = 20, m = 2, bs = "ps", by = L)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.84407 0.01603 302.12 <2e-16 ***
## LVirginia -0.33003 0.02271 -14.54 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
136 CAPÍTULO 8. MODELOS ADITIVOS GENERALIZADOS
## [1] -125.2307
AIC(fit2)
## [1] -131.1844
fit1$sp
## s(dens)1 s(dens)2
## 164.66124239 0.00147493
fit2$sp
plot(dens,log(yield),col=points.cols[location+1],pch=16,cex=.55)
lines(dens.g,fit2.P$fit,col=2,lwd=2)
lines(dens.g,fit2.P$fit+1.96*fit1.P$se.fit,col=2,lwd=2,lty=2)
lines(dens.g,fit2.P$fit-1.96*fit1.P$se.fit,col=2,lwd=2,lty=2)
lines(dens.g,fit2.V$fit,col=4,lwd=2)
lines(dens.g,fit2.V$fit+1.96*fit2.V$se.fit,col=4,lwd=2,lty=2)
lines(dens.g,fit2.V$fit-1.96*fit2.V$se.fit,col=4,lwd=2,lty=2)
8.3.1. Ejemplo:
Continuémos con el ejemplo Wage
gam2 <- gam(wage ~ jobclass + s(age,k=6)+s(year,k=6)+education,data = Wage)
gam3 <- gam(wage ~ jobclass + s(age,by=jobclass,k=6)+s(year,k=6)+education,data = Wage)
summary(gam2)
##
## Family: gaussian
## Link function: identity
##
## Formula:
8.3. MODELOS SEMI-PARAMÉTRICOS 137
s(dens,3.1):LPurnong Landing
1.0
1.0
s(dens,4.73):LVirginia
0.0
0.0
−1.0
dens dens
5.5
log(yield)
4.5
3.5
50 100 150
dens
##
## Family: gaussian
## Link function: identity
##
## Formula:
## wage ~ jobclass + s(age, by = jobclass, k = 6) + s(year, k = 6) +
## education
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 84.039 2.186 38.453 < 2e-16 ***
## jobclass2. Information 4.284 1.352 3.169 0.00154 **
## education2. HS Grad 10.733 2.425 4.427 9.9e-06 ***
## education3. Some College 22.796 2.564 8.892 < 2e-16 ***
## education4. College Grad 36.940 2.571 14.367 < 2e-16 ***
## education5. Advanced Degree 60.621 2.822 21.482 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(age):jobclass1. Industrial 2.980 3.607 33.95 < 2e-16 ***
## s(age):jobclass2. Information 4.551 4.901 20.27 < 2e-16 ***
## s(year) 1.000 1.000 14.27 0.000162 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.293 Deviance explained = 29.6%
## GCV = 1237.2 Scale est. = 1231.3 n = 3000
8.3. MODELOS SEMI-PARAMÉTRICOS 139
par(mfrow=c(2,2))
plot(gam2,all.terms=TRUE,scheme=1)
s(age,4.41)
s(year,1)
−10
−10
−40
−40
20 30 40 50 60 70 80 2003 2005 2007 2009
age year
jobclass education
30 60
30 60
0
jobclass education
plot(gam3,all.terms=TRUE,scheme=1)
s(age,4.55):jobclass2. Information
s(age,2.98):jobclass1. Industrial
−10
−10
−50
−50
20 30 40 50 60 70 80 20 30 40 50 60 70 80
age age
jobclass
Partial for jobclass
30 60
s(year,1)
−10
−50
year jobclass
140 CAPÍTULO 8. MODELOS ADITIVOS GENERALIZADOS
education
30 60
0
1. < HS Grad 4. College Grad
education
Ejercicio: Selecciona un modelo GAM.
100
Ozone
0
200
Solar.R
0
15
Wind
5
90
Temp
60
Month 9
7
5
Day
15
0
0 50 150 5 10 20 5 6 7 8 9
Una simple gráfica de dispersión muestra que algunas de las variables tienen una relación no lineal.
Consideremos un modelo lineal para la variable de respuesta Ozone
airq.lm <- lm(Ozone ~ Temp + Wind + Solar.R, data=airquality)
summary(airq.lm)
##
## Call:
## lm(formula = Ozone ~ Temp + Wind + Solar.R, data = airquality)
##
## Residuals:
8.4. EJEMPLO: CALIDAD DEL AIRE 141
40
40
20
20
20
Partial for Solar.R
Partial for Temp
0
−20
−20
−20
−40
−40
−40
60 70 80 90 5 10 15 20 0 50 150 250
También podemos graficar los residuos del modelo para comprobar las hipótesis del modelo.
par(mfrow=c(1,2))
plot(airq.lm,which=1:2)
142 CAPÍTULO 8. MODELOS ADITIVOS GENERALIZADOS
100
5
117 117
4
Standardized residuals
3
62
30 62
50
30
Residuals
2
1
0
0
−2 −1
−50
−20 20 60 100 −2 −1 0 1 2
La falta de normalidad en los residuos se debe a la asimetría de la variable de respuesta Ozone, podemos
aplicar una transformación log para conseguir asimetría, es decir:
par(mfrow=c(1,2))
hist(airquality$Ozone, main="Ozone")
hist(log(airquality$Ozone), main ="log(Ozone)")
Ozone log(Ozone)
20
30
15
Frequency
Frequency
20
10
10
5
0
0 50 100 150 0 1 2 3 4 5
airquality$Ozone log(airquality$Ozone)
8.4. EJEMPLO: CALIDAD DEL AIRE 143
Ahora podemos ajustar un modelo lineal del log(Ozono), ¿el modelo se ve más adecuado?
lairq.lm <- lm(log(Ozone)~ Temp + Wind + Solar.R, data=airquality)
summary(lairq.lm)
##
## Call:
## lm(formula = log(Ozone) ~ Temp + Wind + Solar.R, data = airquality)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.06193 -0.29970 -0.00231 0.30756 1.23578
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.2621323 0.5535669 -0.474 0.636798
## Temp 0.0491711 0.0060875 8.077 1.07e-12 ***
## Wind -0.0615625 0.0157130 -3.918 0.000158 ***
## Solar.R 0.0025152 0.0005567 4.518 1.62e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5086 on 107 degrees of freedom
## (42 observations deleted due to missingness)
## Multiple R-squared: 0.6644, Adjusted R-squared: 0.655
## F-statistic: 70.62 on 3 and 107 DF, p-value: < 2.2e-16
plot(lairq.lm)
Residuals vs Fitted
24
117
1
Residuals
0
−1
−2
21
Fitted values
lm(log(Ozone) ~ Temp + Wind + Solar.R)
144 CAPÍTULO 8. MODELOS ADITIVOS GENERALIZADOS
Normal Q−Q
3 24
117
2
Standardized residuals
1
0
−4 −3 −2 −1
21
−2 −1 0 1 2
Theoretical Quantiles
lm(log(Ozone) ~ Temp + Wind + Solar.R)
Scale−Location
21
2.0
Standardized residuals
24
1.5
117
1.0
0.5
0.0
Fitted values
lm(log(Ozone) ~ Temp + Wind + Solar.R)
8.4. EJEMPLO: CALIDAD DEL AIRE 145
Residuals vs Leverage
2 24
48
Standardized residuals
0
−2
0.5
−4
Cook's distance 21
Leverage
lm(log(Ozone) ~ Temp + Wind + Solar.R)
Ajustando un modelo lineal podemos concluir que puede haber alguna heterokedasticidad que no se tenga
en cuenta en el modelo. Sin embargo, la causa más posible es que la relación entre la variable de respuesta
y las covariables está lejos de ser lineal.
Ajustemos un modelo GAM, en primer lugar con una sola variable, es decir:
library(mgcv)
airq.gam1 <- gam(log(Ozone) ~ s(Wind,bs="ps",m=2,k=10),
method="REML", select=TRUE,data=airquality)
summary(airq.gam1)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## log(Ozone) ~ s(Wind, bs = "ps", m = 2, k = 10)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.4185 0.0655 52.19 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Wind) 2.565 9 6.453 2.76e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.336 Deviance explained = 35%
146 CAPÍTULO 8. MODELOS ADITIVOS GENERALIZADOS
Los resultados muestran que el efecto suave s(Wind) es significativo (con un valor aproximado de p cercano
a cero y edf 2,56). La siguiente figura muestra el efecto de viento suave. El incremento de la velocidad del
viento disminuye los niveles de ozono (dramáticamente hasta 10mph). Tener en cuenta que incluyendo el
intercepto, los efectos suaves se centran en cero.
plot(airq.gam1,residuals=TRUE,scheme=1)
2
1
s(Wind,2.56)
0
−1
−2
−3
5 10 15 20
Wind
gam.check(airq.gam1)
8.4. EJEMPLO: CALIDAD DEL AIRE 147
1
0
deviance residuals
−1
−2
−3
−2 −1 0 1 2
theoretical quantiles
Resids vs. linear pred.
1
0
residuals
−1
−2
−3
linear predictor
148 CAPÍTULO 8. MODELOS ADITIVOS GENERALIZADOS
Histogram of residuals
35
30
25
Frequency
20
15
10
5
0
−3 −2 −1 0 1
Residuals
8.4. EJEMPLO: CALIDAD DEL AIRE 149
2
1
0
Fitted Values
##
## Method: REML Optimizer: outer newton
## full convergence after 8 iterations.
## Gradient range [-1.455961e-06,1.183099e-06]
## (score 128.6926 & scale 0.4976891).
## Hessian positive definite, eigenvalue range [0.4379927,57.51548].
## Model rank = 10 / 10
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(Wind) 9.00 2.56 0.69 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Nótese que, en el modelo airq.gam1 usamos k=10 nodos, la variable Wind tiene valores únicos de 31, y por
lo tanto 10 nodos deberían ser suficientes. Ajustamos un modelo para los residuos del modelo gam ajustado
resids <- residuals(airq.gam1)
resids.gam <- gam(resids~s(Wind,k=20,m=2),method="REML",
select=TRUE,data=airq.gam1$model)
summary(resids.gam)
##
## Family: gaussian
150 CAPÍTULO 8. MODELOS ADITIVOS GENERALIZADOS
−0.001
−0.003
5 10 15 20
Wind
Ahora añadimos la funcion Temp:
airq.gam2 <- gam(log(Ozone)~s(Wind,bs="ps",m=2,k=10)+s(Temp,bs="ps",m=2,k=10),
method="REML",select=TRUE,data=airquality)
summary(airq.gam2)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## log(Ozone) ~ s(Wind, bs = "ps", m = 2, k = 10) + s(Temp, bs = "ps",
## m = 2, k = 10)
8.4. EJEMPLO: CALIDAD DEL AIRE 151
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.41852 0.04963 68.89 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Wind) 2.068 9 1.353 0.000981 ***
## s(Temp) 3.990 9 10.496 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.619 Deviance explained = 63.9%
## -REML = 101.64 Scale est. = 0.28568 n = 116
1
0
0
s(Temp,3.99)
s(Wind,2.07)
−1
−1
−2
−2
−3
−3
5 10 15 20 60 70 80 90
Wind Temp
Podemos comparar ambos modelos utilizando la función anova:
anova(airq.gam1,airq.gam2)
## [1] 255.0315
AIC(airq.gam2)
## [1] 195.6561
Incluimos la variable Solar.R que contiene valores faltantes (NA’s).
sum(is.na(airquality$Solar.R))
## [1] 7
Para comparar el modelo anterior airq.gam2 y el nuevo modelo que incluye Solar.R utilizaremos los mismos
datos, por lo que omitiremos los valores que faltan y volveremos a ajustar los modelos.
new.airquality <- na.omit(airquality)
airq.gam22=gam(log(Ozone)~s(Wind,bs="ps",m=2,k=10)+s(Temp,bs="ps",m=2,k=10),
method="REML",select=TRUE,data=new.airquality)
airq.gam3=gam(log(Ozone)~s(Wind,bs="ps",m=2,k=10)+s(Temp,bs="ps",m=2,k=10)+s(Solar.R,bs="ps",m=2,k=20),
method="REML",select=TRUE,data=new.airquality)
summary(airq.gam22)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## log(Ozone) ~ s(Wind, bs = "ps", m = 2, k = 10) + s(Temp, bs = "ps",
## m = 2, k = 10)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.41593 0.05128 66.61 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Wind) 2.1879 9 1.919 1e-04 ***
## s(Temp) 0.9874 9 8.601 7.73e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.611 Deviance explained = 62.2%
## -REML = 95.215 Scale est. = 0.29194 n = 111
summary(airq.gam3)
##
## Family: gaussian
## Link function: identity
##
8.4. EJEMPLO: CALIDAD DEL AIRE 153
## Formula:
## log(Ozone) ~ s(Wind, bs = "ps", m = 2, k = 10) + s(Temp, bs = "ps",
## m = 2, k = 10) + s(Solar.R, bs = "ps", m = 2, k = 20)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.41593 0.04586 74.49 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Wind) 2.318 9 2.255 2.44e-05 ***
## s(Temp) 1.852 9 6.128 1.12e-12 ***
## s(Solar.R) 2.145 19 1.397 1.23e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.689 Deviance explained = 70.7%
## -REML = 86.106 Scale est. = 0.23342 n = 111
AIC(airq.gam22)
## [1] 185.7769
AIC(airq.gam3)
## [1] 166.0423
par(mfrow=c(2,2))
plot(airq.gam3,residuals=TRUE)
s(Temp,1.85)
s(Wind,2.32)
0
−2
−2
5 10 15 20 60 70 80 90
Wind Temp
s(Solar.R,2.14)
0
−2
Solar.R
154 CAPÍTULO 8. MODELOS ADITIVOS GENERALIZADOS
Capítulo 9
Análisis multivariante
155
156 CAPÍTULO 9. ANÁLISIS MULTIVARIANTE
Recodificamos las pruebas relativas a 3 carreras hurdles, run200m y run800m, restando el valor más alto en
cada carrera, cada uno de los 35 tiempos de los atletas.
heptathlon$hurdles <- max(heptathlon$hurdles) - heptathlon$hurdles
heptathlon$run200m <- max(heptathlon$run200m) - heptathlon$run200m
heptathlon$run800m <- max(heptathlon$run800m) - heptathlon$run800m
Diagrama de dispersion
score <- which(colnames(heptathlon) == "score")
plot(heptathlon[,-score])
1.50 1.70 0 2 4 36 42
hurdles
2
0
1.50 1.85
highjump
15
shot
10
0 2 4
run200m
7.0
longjump
5.0
36 44
javelin
30
run800m
0
0 1 2 3 10 13 16 5.0 6.0 7.0 0 20 40
Matriz de correlación
round(cor(heptathlon[,-score]),3)
(PNG) de Papúa Nueva Guinea - eliminaremos esta observación para ver si la matriz de correlación es
significativamente diferente:
heptathlon <- heptathlon[-which(rownames(heptathlon)=="Launa (PNG)"),]
plot(heptathlon[,-score])
1.70 1.80 0 2 4 36 42
3.5
hurdles
1.5
highjump
1.70
15
shot
10
0 2 4
run200m
7.0
longjump
5.5
36 44
javelin
20 35
run800m
Eliminando al atleta de Papúa Nueva Guinea, las correlaciones cambian sustancialmente y en el diagrama
de dispersión matricial no se observan valores extremos.
round(cor(heptathlon[,-score]),3)
## $sdev
## [1] 2.0793370 0.9481532 0.9109016 0.6831967 0.5461888 0.3374549 0.2620420
##
## $rotation
## PC1 PC2 PC3 PC4 PC5
158 CAPÍTULO 9. ANÁLISIS MULTIVARIANTE
summary
summary(heptathlon_pca)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Standard deviation 2.0793 0.9482 0.9109 0.68320 0.54619 0.33745
## Proportion of Variance 0.6177 0.1284 0.1185 0.06668 0.04262 0.01627
## Cumulative Proportion 0.6177 0.7461 0.8646 0.93131 0.97392 0.99019
## PC7
## Standard deviation 0.26204
## Proportion of Variance 0.00981
## Cumulative Proportion 1.00000
Comprobamos la variabilidad explicada:
head(heptathlon_pca$x)
heptathlon_pca
4
3
Variances
2
1
0
El biplot es una representación gráfica de datos multivariantes. Así como una gráfica de dispersión muestra
la distribución combinada de dos variables, una biplot representa tres o más variables.
Si ordenamos de mayor a menor la variable “puntuación” tenemos las tres medallas de oro, plata y bronce.
head(heptathlon[order(heptathlon$score,decreasing = TRUE),],3)
El biplot, nos muestra los atletas proyectados en sus dos primeros componentes principales, pero también
las flechas nos dan información sobre las varianzas y covarianzas de las variables (direcciones de máxima
variabilidad).
biplot(heptathlon_pca)
9.1. ANÁLISIS DE COMPONENTES PRINCIPALES 161
−4 −2 0 2 4
0.4
4
run800m
Hui−Ing (TAI)
Choubenkova (URS)
0.2
John (GDR)
2
Mulliner (GB)
Bouraga (URS)
Behmer (GDR) Hautenauve (BEL)
Dimitrova
Fleming(BUL)
(AUS)
run200m Greiner (USA) Kytola (FIN) Jeong−Mi (KOR)
Schulz (GDR)
Sablovskaite (URS)
hurdles
Geremias (BRA)
0.0
0
longjumpshot Hagger (GB)
Joyner−Kersee (USA)
PC2
Wijnsma (HOL)
javelin Ruotsalainen
Lajbnerova (CZE) (FIN)
−0.2
−2
highjump Scheider (SWI)
Braun (FRG)
−0.4
−4
Yuping (CHN)
Brown (USA)
PC1
Por ejemplo, el ganador Joyner-Kersee (USA) acumula puntuaciones más altas en el longjump,hurdlesp
y run200m.
Podemos analizar la correlación entre la puntuación de la variable y PC1. Esto indica que la correlación es
muy negativa y muy fuerte con el score.
cor(heptathlon$score, heptathlon_pca$x[,1])
## [1] -0.9931168
162 CAPÍTULO 9. ANÁLISIS MULTIVARIANTE
plot(heptathlon$score, heptathlon_pca$x[,1])
4
heptathlon_pca$x[, 1]
2
0
−2
−4
heptathlon$score
objeto prcomp
class(heptathlon_pca)
## [1] "prcomp"
−4 −2 0 2 4
0.4
4
run800m
Hui−Ing (TAI)
Choubenkova (URS)
0.2
John (GDR)
2
Mulliner (GB)
Bouraga (URS)
Behmer (GDR) Hautenauve (BEL)
Dimitrova
Fleming(BUL)
(AUS)
run200m Greiner (USA) Kytola (FIN) Jeong−Mi (KOR)
Schulz (GDR)
Sablovskaite (URS)
hurdles
Geremias (BRA)
0.0
0
longjumpshot Hagger (GB)
Joyner−Kersee (USA)
PC2
Wijnsma (HOL)
−2
highjump Scheider (SWI)
Braun (FRG)
−0.4
−4
Yuping (CHN)
Brown (USA)
PC1
Es un método alternativo más adecuado a la regresión logística cuando la variable cualitativa tiene más de
dos niveles (K ≥ 2). Supone también un modelo más estable cuando el tamaño muestral n es pequeño y la
distribución de los predictores es aproximadamente normal en cada una de sus clases. El propósito del LDA
es encontrar la combinación lineal de las variables originales que permita la mejor separación entre grupos
de un set de datos. El LDA está basado en el clasificador Bayesiano.
164 CAPÍTULO 9. ANÁLISIS MULTIVARIANTE
Ilustración
library(ggplot2)
ggplot(data = data.frame(x = -4:4), aes(x = x)) +
# distribución de densidad de k = 1
stat_function(fun = dnorm,
args = list(mean = -1.25, sd = 1),
color = "green3") +
# distribución de densidad de k = 2
stat_function(fun = dnorm,
args = list(mean = 1.25, sd = 1),
color = "firebrick") + # límite de decisión de Bayes
geom_vline(xintercept = 0,
linetype = "longdash") +
theme_bw()
0.4
0.3
0.2
y
0.1
0.0
−4 −2 0 2 4
x
En este caso, el clasificador de Bayes asignará la observación a la clase k = 1 (verde) si x < 0, y a la clase
k = 2 (rojo) si x > 0.
9.2.1. Ejemplo:
Vamos a predecir si el rendimiento del combustible (gas mileage) de un automóvil es alto o bajo en función
del resto de predictores del set de datos.
library(ISLR)
head(Auto, 3)
attach(Auto) # Vector de “0”s (rendimiento bajo) con la misma longitud que la variable mpg
200
7
300
6
5
50 100
4
100
3
0 1 0 1 0 1
82
4500
20
78
3000
15
74
10
1500
70
0 1 0 1 0 1
60
count mpg01
40
0
20 1
0
100 200 300 400
displacement
40
30 mpg01
count
20 0
10 1
0
50 100 150 200
horsepower
30
mpg01
count
20
0
10 1
0
2000 3000 4000 5000
weight
modelo.lda
## Call:
## lda(mpg01 ~ cylinders + displacement + horsepower + weight, data = Auto.train)
##
## Prior probabilities of groups:
## 0 1
## 0.514377 0.485623
##
## Group means:
## cylinders displacement horsepower weight
## 0 6.677019 269.8323 128.71429 3601.379
9.3. K-MEDIAS 171
9.3. K-medias
K-means Clustering es un algoritmo de aprendizaje no supervisado que intenta agrupar datos basados en su
similitud. El aprendizaje sin supervisión significa que no hay resultados que predecir, y el algoritmo sólo trata
de encontrar patrones en los datos. En k-means clustering, tenemos que especificar el número de clusters
en los que queremos que se agrupen los datos. El algoritmo asigna aleatoriamente cada observación a un
cúmulo, y encuentra el centroide de cada cúmulo. Luego, el algoritmo itera a través de dos pasos:
1. Reasigne los puntos de datos al cluster cuyo centroide está más cerca.
2. Calcular el nuevo centroide de cada cúmulo.
Estos dos pasos se repiten hasta que la variación dentro del conglomerado no pueda reducirse más. La
variación dentro del conglomerado se calcula como la suma de la distancia euclídea entre los puntos de datos
y sus respectivos centros de conglomerados.
Ejemplo:
library(rattle)
data(wine)
head(wine)
## $names
## [1] "cluster" "centers" "totss" "withinss"
## [5] "tot.withinss" "betweenss" "size" "iter"
## [9] "ifault"
##
172 CAPÍTULO 9. ANÁLISIS MULTIVARIANTE
## $class
## [1] "kmeans"
Una pregunta fundamental es cómo determinar el valor del parámetro k. Si consideramos el porcentaje de
varianza explicado en función del número de grupos:
Uno debe elegir un número de grupos para que la incorporación de otro grupo no dé un mejor modelado
de los datos. Más precisamente, si el porcentaje de varianza explicado por los clusters se traza de acuerdo
al número de grupos, los primeros clusters añadirán mucha información (explican mucha varianza), pero en
algún momento la ganancia marginal disminuirá, dando un ángulo en el gráfico. En este punto se elige el
número de racimos, de ahí el “criterio del codo”.
wssplot <- function(data, nc=15, seed=1234){
wss <- (nrow(data)-1)*sum(apply(data,2,var))
for (i in 2:nc){
set.seed(seed)
wss[i] <- sum(kmeans(data, centers=i)$withinss)}
plot(1:nc, wss, type="b", xlab="Number of Clusters",
ylab="Within groups sum of squares")}
wssplot(wine.stand, nc=6)
2200
Within groups sum of squares
1800
1400
1000
1 2 3 4 5 6
Number of Clusters
1 159
178 176 2 4
167
170
151 160 19
174150 177
173
154 169 17 6 15
157152
156 165 34854
1640
2
175 153
149 46 18 50 59 53
158168 162
148 49 56
3132
57 111
145 58 7 47 43 14
Component 2
166161 164
172 26
5 37 20
55 3
9
138
147 146
163 35 27
12 29
13 10
4152
4821
143
136
139 155
141
133 144
132134
142
140 44 122
22 74
137 131 42 38 36453096 51
171 84 6997 3 28 33 2 23
0
135 25
24
123
6271113 124 66 7939
73 103 80 72
119 130128
108
61
78 63
127 121
125
82
85 111
11075
7099
93
106 89
92 65
129
115 118
112 86 100 6467
−2
91 87 120
114
102 126
107109105
77 101
9495
7688
83 68 98
90
104
60 117
81
116
−4
−4 −2 0 2 4
Component 1
These two components explain 55.41 % of the point variability.
Sabiendo que hay tres tipos de wine$Type wines, podemos calcular la matriz de confusión.
table(wine$Type)
##
## 1 2 3
## 59 71 48
table(wine[,1],k.means.fit$cluster)
##
## 1 2 3
## 1 0 59 0
## 2 3 3 65
## 3 48 0 0
El criterio de desviación mínima de Ward minimiza la desviación total dentro del grupo de empresas.
H.fit <- hclust(d, method="ward")
## The "ward" method has been renamed to "ward.D"; note new "ward.D2"
Height
##
##
##
##
##
174
0 20 40 60 80 100 120
159
160
154
176
177
groups
class(H.fit)
168
172
165
173
157
149
175
3 0 0 48
2 7 58 6
1 58 1 0
1 2 3
## [1] "hclust"
170
178
167
169
153
151
150
152
162
141
143
140
163
84
144
table(wine[,1],groups)
139
166
119
147
148
161
156
174
158
137
138
155
135
136
62
61
69
145
142
146
plot(H.fit) # display dendogram
164
17197
131
133
132
134
49
563
31
34
16
54
rect.hclust(H.fit, k=3, border="red")
17
18
14
15
43
59
534
11
32
506
19
36
24
25
33
66
27
58
37
35
382
28
39
52
12
d
1347
2379
30
1
21
41
57
La opción plot devuelve un hclust que muestra el dendograma:
10
45
55
22
Cluster Dendrogram
42
44
20
408
46
96
70
79
122
74
72
26
295
106
93
108
89
91
92
120
13071
78
65
113
87
83
88
114
90
115
100
116
128
80
123 60
63
77
76
101
94
95
81
98
82
86
104
109
102
68
107
105
11775
64
99
51
67
111
124
125
121
85
110
73
118
103
112
126
127
CAPÍTULO 9. ANÁLISIS MULTIVARIANTE
129
Capítulo 10
En esta sección describimos una clase de métodos de aprendizaje que se desarrollaron por separado en
diferentes campos: estadística e inteligencia artificial basadas en modelos esencialmente idénticos.
Las redes neuronales artificiales se inspiran en el comportamiento conocido del cerebro humano (principal-
mente el referido a las neuronas y sus conexiones), trata de crear modelos artificiales que solucionen problemas
difíciles de resolver mediante técnicas algorítmicas convencionales.
La idea central es extraer las combinaciones lineales de las variables de entrada como características derivadas
y, a continuación, modelar el objetivo como una función no lineal de estas características. El resultado es un
poderoso método de aprendizaje, con amplias aplicaciones en muchos campos.
la neurona artificial pretende mimetizar las características más importantes de la neurona biólogica. En
general, recibe las señales de entrada de las neuronas v ecinas ponderadas por los pesos de las conexiones. La
suma de estas señales ponderadas proporciona la entrada total o neta de la neurona y, mediante la aplicación
de una función matemática - denominada función de salida - , sobre la entrada neta, se calcula un valor de
salida, el cual es enviado a otras neuronas.
Tanto los valores de entrada a la neurona como su salida pueden ser señales excitatorias (cuando el valor es
positivo) o inhibitorias (cuando el valor es negativo).
175
176 CAPÍTULO 10. INTRODUCCION A LAS REDES NEURONALES ARTIFICIALES
es predecir el valor medio de las viviendas ocupadas por sus propietarios (‘medv) utilizando todas las demás
variables continuas disponibles.
set.seed(500)
library(MASS)
data <- Boston
Primero tenemos que comprobar que no hay datos faltantes, de lo contrario tenemos que arreglar el conjunto
de datos.
apply(data,2,function(x) sum(is.na(x)))
##
## Call:
## glm(formula = medv ~ ., data = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -14.9143 -2.8607 -0.5244 1.5242 25.0004
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 43.469681 6.099347 7.127 5.50e-12 ***
## crim -0.105439 0.057095 -1.847 0.065596 .
## zn 0.044347 0.015974 2.776 0.005782 **
## indus 0.024034 0.071107 0.338 0.735556
## chas 2.596028 1.089369 2.383 0.017679 *
## nox -22.336623 4.572254 -4.885 1.55e-06 ***
## rm 3.538957 0.472374 7.492 5.15e-13 ***
## age 0.016976 0.015088 1.125 0.261291
## dis -1.570970 0.235280 -6.677 9.07e-11 ***
## rad 0.400502 0.085475 4.686 3.94e-06 ***
## tax -0.015165 0.004599 -3.297 0.001072 **
## ptratio -1.147046 0.155702 -7.367 1.17e-12 ***
## black 0.010338 0.003077 3.360 0.000862 ***
## lstat -0.524957 0.056899 -9.226 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 23.26491)
##
178 CAPÍTULO 10. INTRODUCCION A LAS REDES NEURONALES ARTIFICIALES
La función sample(x,size) simplemente produce un vector del tamaño especificado de muestras selecciona-
das aleatoriamente desde el vector x. Por defecto el muestreo es sin reemplazo: el índice es esencialmente un
vector aleatorio de muestras.
Como se trata de un problema de regresión, vamos a utilizar el error cuadrático medio (ECM) como medida
de lo lejos que están nuestras predicciones de los datos reales.
No hay una regla fija en cuanto a cuántas capas y neuronas usar, aunque hay varias reglas generales más
o menos aceptadas. Normalmente, si es necesario, una capa oculta es suficiente para un gran número de
aplicaciones.
En cuanto al número de neuronas, debería estar entre el tamaño de la capa de entrada y el tamaño de la
capa de salida, normalmente 2/3 del tamaño de entrada. No hay garantía de que ninguna de estas reglas se
ajuste mejor a su modelo de modo que es recomendable probar con diferentes combinaciones.
Para este ejemplo, vamos a usar 2 capas ocultas con esta configuración: 13:5:3:1.
La capa de entrada tiene 13 entradas,
las dos capas ocultas tienen 5 y 3 neuronas y
la capa de salida tiene, por supuesto, una sola salida ya que estamos haciendo regresión.
Vamos a entrar en la red:
library(neuralnet)
n <- names(train_)
f <- as.formula(paste("medv ~", paste(n[!n %in% "medv"], collapse = " + ")))
nn <- neuralnet(f,data=train_,hidden=c(5,3),linear.output=TRUE)
10.4. PREDICCIÓN CON UNA RED NEURONAL 179
Nota:
La expresión y~. no es aceptada por la función neuralnet().
El argumento oculto acepta un vector con el número de neuronas para cada capa oculta, mien-
tras que el argumento linear.output se usa para especificar si queremos hacer una regresión
linear.output=TRUE o clasificación linear.output=FALSE.
La librería neuralnet proporciona una buena herramienta para representar gráficamente el modelo con los
pesos en cada conexión:
plot(nn)
Las líneas negras muestran las conexiones entre cada capa y los pesos de cada conexión, mientras que las líneas
azules muestran el término de sesgo añadido en cada paso. El sesgo puede ser pensado como el intercepto de
un modelo lineal.
La red neuronal es esencialmente una caja negra, por lo que no podemos decir mucho sobre el ajuste, los
pesos y el modelo. Basta decir que el algoritmo de entrenamiento ha convergido y por lo tanto el modelo
está listo para ser utilizado.
## [1] 15.75184
Podemos ahora comparar los dos ECM:
print(paste(MSE.lm,MSE.nn))
40
50
30
40
pr.nn_
pr.lm
30
20
20
10
10
NN LM
10 20 30 40 50 10 20 30 40 50
Inspeccionando visualmente podemos ver que las predicciones hechas por la red neural están (en general)
más concentradas alrededor de la línea (una alineación perfecta con la línea indicaría un ECM de 0 y por lo
tanto una predicción perfecta ideal) que las hechas por el modelo lineal.
plot(test$medv,pr.nn_,col='red',main='Real vs predicted NN',pch=18,cex=0.7)
points(test$medv,pr.lm,col='blue',pch=18,cex=0.7)
abline(0,1,lwd=2)
legend('bottomright',legend=c('NN','LM'),pch=18,col=c('red','blue'))
10.5. EVALUANDO LA PREDICCIÓN MEDIANTE VALIDACIÓN CRUZADA 181
Real vs predicted NN
50
40
pr.nn_
30
20
NN
10
LM
10 20 30 40 50
test$medv
## [1] 23.8356
set.seed(450)
cv.error <- NULL
k <- 10
for(i in 1:k){
index <- sample(1:nrow(data),round(0.9*nrow(data)))
182 CAPÍTULO 10. INTRODUCCION A LAS REDES NEURONALES ARTIFICIALES
nn <- neuralnet(f,data=train.cv,hidden=c(5,2),linear.output=T)
Calculamos el ECM promedio y graficamos los resultados como una gráfica de caja.
mean(cv.error)
## [1] 10.32698
cv.error
6 8 10 12 14 16 18
MSE CV
Como se puede ver, la media del ECM de la red neural (10.33) es inferior al del modelo lineal, aunque parece
haber cierto grado de variación en laos ECMs de la validación cruzada. Esto puede depender de la división de
10.5. EVALUANDO LA PREDICCIÓN MEDIANTE VALIDACIÓN CRUZADA 183
los datos o de la inicialización aleatoria de los pesos en la red neuronal. Al ejecutar la simulación en diferentes
momentos con diferentes semillas, puede obtener una estimación puntual más precisa para el ECM medio.
1. 2 capas ocultas de 5 y 3.
2. 1 capa oculta de 8
Boston.nn.5.3 <- neuralnet(f
, data=train_
, hidden=c(5,3)
, linear.output=TRUE)
Como alternativa a la función plot.nn(), la librería NeuralNetTools incluye funciones gráficas más elegantes.
Este elegante gráfico resuelve el problema del desorden visual utilizando el grosor de la línea para representar
la magnitud del peso y el color de la línea para representar el signo de peso (negro = positivo, gris = negativo).
library(NeuralNetTools)
plotnet(Boston.nn.5.3)
crim B1 B2 B3
I1
zn I2
indus I3
chas I4
nox I5 H1
rm I6 H2 H1
age I7 H3 H2 O1 medv
dis I8 H4 H3
rad I9 H5
tax I10
ptratio I11
black I12
lstat I13
plotnet(Boston.nn.8)
184 CAPÍTULO 10. INTRODUCCION A LAS REDES NEURONALES ARTIFICIALES
crim B1 B2
I1
zn I2
indus I3
chas H1
I4
nox H2
I5
rm H3
I6
age H4
I7 O1 medv
dis H5
I8
rad H6
I9
tax H7
I10
ptratio H8
I11
black I12
lstat I13
Capítulo 11
Extras
185