Curso Udemy
Jose Antonio Constantino Jaime
25/1/2022
#1 Generalidades de R
##1.1 Cargar librerias
library("XML")
library("tidyverse")#Consultar donde esta ubicado un archivo
getwd()## [1] "C:/Users/Antonio Constantino/Desktop/Rstudio/Cursos/curso youtube"#Como cambiar esa direccion
setwd("C:/Users/Antonio Constantino/Desktop/Rstudio/Cursos/curso youtube/")Como se lee un Csv
Dentro de la funcion read , sep= indica por que simbolo separa los datos Nota de oro : cuando se tiene los datos separados por un “;” entoces se usa la funcion read.ccv2 y en el caso que este separado por tabuladores , entonces se debe usar dento de la funcion como “sep=”t”
auto1 <-read.csv("../../r-course-master/data/tema1/auto-mpg.csv",header = TRUE,sep = ",")Cuando no tiene encabezado el archivo
sin_encabezado <- read.csv("../../r-course-master/data/tema1/auto-mpg-noheader.csv")Lo que veras es que pone una x de un primer dato ( la coma y el cuato lo que te dice es que solo te muestre las primeras 4 filas )
head(sin_encabezado,4)## X1 X28 X4 X140 X90 X2264 X15.5 X71 chevrolet.vega.2300
## 1 2 19 3 70 97 2330 13.5 72 mazda rx2 coupe
## 2 3 36 4 107 75 2205 14.5 82 honda accord
## 3 4 28 4 97 92 2288 17.0 72 datsun 510 (sw)
## 4 5 21 6 199 90 2648 15.0 70 amc gremlin#Como introduzco el encabezado
#Se va a reescribir
sin_encabezado <- read.csv("../../r-course-master/data/tema1/auto-mpg-noheader.csv",header = F,col.names = c("numero","MPG","Cilindros","desplazamiento","Caballos_de_potencia","peso","aceleracion","ano","modelo"))
# aqui lo que haces es :header( carnal no tienes encabezados ),colnames = ( te voy a poner nombre c( estos son los nombres )) Ahora veremos como se observan este dato
head(sin_encabezado,4)## numero MPG Cilindros desplazamiento Caballos_de_potencia peso aceleracion ano
## 1 1 28 4 140 90 2264 15.5 71
## 2 2 19 3 70 97 2330 13.5 72
## 3 3 36 4 107 75 2205 14.5 82
## 4 4 28 4 97 92 2288 17.0 72
## modelo
## 1 chevrolet vega 2300
## 2 mazda rx2 coupe
## 3 honda accord
## 4 datsun 510 (sw)#ya quedo Ahora : que chingao hago cuando falta un dato Pasan 2 cosas Cuando tengo datos de tipo numero :asi aparecen como N/A ( Not Available “No disponible”) Y categoricas ( stirngs o letras ):se tiene que poner (na.sting) Tambien puede pasar que el string lo este tomando como factor ( letras como numero (por ejemplo el numero de cilindros de un carro)) Hagamoslo
sin_encabezado <- read.csv("../../r-course-master/data/tema1/auto-mpg-noheader.csv",stringsAsFactors = TRUE,na.strings = "null",header = F,col.names = c("numero","MPG","Cilindros","desplazamiento","Caballos_de_potencia","peso","aceleracion","ano","modelo"))as.character: es para cambiar variables categoricas a caracteres
#Esto escribe el nombre de la cabecera
names(auto1)## [1] "No" "mpg" "cylinders" "displacement" "horsepower"
## [6] "weight" "acceleration" "model_year" "car_name"Como ver el numero de filas que tiene el DF
nrow(auto1)## [1] 398#Como extraer datos de una pagina web o “scraping”
library(XML)
#esto es un valor
#Nota una coleccion de valores estructurado es lo que se llama un data frame ( como solo un valor de una celda de excel)
url <- "../../r-course-master/data/tema1/cd_catalog.xml"
#esto es una varibale que casualmete se llama igual que la funcion
#esta variable es un "apuntador" que localiza el docmuento y debe pasar por todod los nodos
xmldoc<- xmlParse(url)
#Se debe obtener la raiz
rootnode<- xmlRoot(xmldoc)
#Aqui puedes ver una parte de la estuctura
rootnode[2]## $CD
## <CD>
## <TITLE>Hide your heart</TITLE>
## <ARTIST>Bonnie Tyler</ARTIST>
## <COUNTRY>UK</COUNTRY>
## <COMPANY>CBS Records</COMPANY>
## <PRICE>9.90</PRICE>
## <YEAR>1988</YEAR>
## </CD>
##
## attr(,"class")
## [1] "XMLInternalNodeList" "XMLNodeList"Creemos un data frame desde un xml
cds_data <- xmlSApply(rootnode,function(x)xmlSApply(x,xmlValue))
#xmlValue Estrae informacion del nodo raiz
# la funcion xmlApply hace una matriz por este motivo en las variables de la derecha se tiene un DF pero sin el boton de play
View(cds_data)#Ahora tenemos que trasponerlo
cds_catalog <- data.frame(t(cds_data),row.names = NULL)
cds_catalog[1:5,]## TITLE ARTIST COUNTRY COMPANY PRICE YEAR
## 1 Empire Burlesque Bob Dylan USA Columbia 10.90 1985
## 2 Hide your heart Bonnie Tyler UK CBS Records 9.90 1988
## 3 Greatest Hits Dolly Parton USA RCA 9.90 1982
## 4 Still got the blues Gary Moore UK Virgin records 10.20 1990
## 5 Eros Eros Ramazzotti EU BMG 9.90 1997Otras funciones que no se explican en este curso pero que
Nota : Los corchetes se leen asi [Filas,Columnas], fin del comunicado jejeje son utiles para XMLson: xpathSApply(): getNodeSet():
#Ahora agamos un ejercicio para obtener valores de wikipedia
#population_url<-readHTMLTable("../../r-course-master/data/tema1/WorldPopulation-wiki.htm")
#tables<- (population_url)
#View(tables)Esto te crea a una lista de tablas que estan dentro de la paguina ahora seleccionaremos una tabla dentro de ese listado
#most_populated<-tables[[6]]
#head(most_populated)Ahora : si no quieres la lista y despues extarer la tabla , osea , ya sabes dentro del xml cual es el numero de tabla que quieres
#custom_table <- readHTMLTable(population_url,which = 6)
#head(custom_table)Esto te puede ayudar a identificar si la paguina tiene demaciado peso y no se vuelva lenta tu computadora entonces puede usarla
Ficheros tipo json
library("jsonlite")##
## Attaching package: 'jsonlite'## The following object is masked from 'package:purrr':
##
## flattendat.1<- fromJSON("../../r-course-master/data/tema1/students.json")
dat.2<- fromJSON("../../r-course-master/data/tema1/student-courses.json")Para obtener valores de Internet en formato jason
urljason <- "http://www.floatrates.com/daily/usd.json"Aqui solo lo extrae : ahora lo que hace falta es meterlo a una lista
library("curl")## Using libcurl 7.64.1 with Schannel##
## Attaching package: 'curl'## The following object is masked from 'package:readr':
##
## parse_datedfdeurl <- fromJSON(urljason)Ahora debes tener cuidado con este tema ya que depende mucho de la estuctura del JSON por que dentro de cada llave puede tener catergorias
Lo leeremos y lo sobre escibiremos a la mismo DF ( Esta madre aun no me sale )
#Aqui tienes que hacer una conversion
dfdeurl <- dfdeurl$eur$code
#Esto me crea un valor , pero quiero un DFPara convertir un un dataframe a JSON
#conversion.df.JSON<- toJSON(mos_populated)
#head(mos_populated)Leer archivo con Formato con ancho fijo (.TXT) fwf
los widths=ancho ( los numeros dentro de la c son los numeros de carateres por columna )
student<- read.fwf("C:/Users/Antonio Constantino/Desktop/Rstudio/curso de youtube/curso youtube/../../r-course-master/data/tema1/student-fwf.txt",
width=c(4,15,20,15,4),
col.names = c("id","name","email","carrera","año"))ahora si tuviera encabezados Sep=separador
student1 <- read.fwf("C:/Users/Antonio Constantino/Desktop/Rstudio/curso de youtube/curso youtube/../../r-course-master/data/tema1/student-fwf.txt",width=c(4,15,20,15,4),header = F)
head(student1)## V1 V2 V3 V4 V5
## 1 1044 Fran Morton fran.m@quama.edu Statistics 2015
## 2 1035 Kato Mullins tempor@giat.com Marketing 2014
## 3 1046 Jordan Keller Cum@gmail.com Chemistry 2014
## 4 1084 Macey Simpson eu@famesacturpis.com Finance 2016
## 5 1010 Deanna Ortega Dea.Ortega@velit.ca Chemistry 2016
## 6 1096 Joe Eaton je@euismod.org MIS 2014Notaras que pone una V en el lufar de la fila que no tiene nombre
Te quieres saltar algunas lineas agregas skip Si quieres excluir una columna agregas los que quieras eliminar Nota Importante :Si usas agregar las columnas aparte de agregar el menos debes borrar el nombre de la columna
student2<- read.fwf("C:/Users/Antonio Constantino/Desktop/Rstudio/curso de youtube/curso youtube/../../r-course-master/data/tema1/student-fwf.txt",width=c(4,15,-20,15,4),header = F,sep = "\t",skip = 2)
head(student2)## V1 V2 V3 V4
## 1 1046 Jordan Keller Chemistry 2014
## 2 1084 Macey Simpson Finance 2016
## 3 1010 Deanna Ortega Chemistry 2016
## 4 1096 Joe Eaton MIS 2014
## 5 1051 Yoko Washington Marketing 2017
## 6 1074 Ebony Farrell Chemistry 2017Aqui se salto 2 y elimino la fila del correo
Como se crean las variables
Clientes <- c("juan gabriel","ricardo","pedro ")
fechas <- as.Date(c("2017-12-27","2017-11-1","2018-1-1"))
pago<- c(315,192.55,40.15)crear data frame
pedidos <- data.frame(Clientes,fechas,pago)Como se guarda el data frame
este es formato .Rdata
save(pedidos,file = "../../r-course-master/data/tema1/pedidos.Rdata ")#m Como se guarda el data frame este es formato .rds
saveRDS(pedidos,file = "../../r-course-master/data/tema1/pedidos.rds")Asi de borra una variable o un DF
# remove (pedidos )Como se Carga una varibale ya trabajada
No es lo mismo guardarlo que “mandarlo llamar” cargarlo para trabajarlo en R
Para formatos .Rdata
load("../../r-course-master/data/tema1/pedidos.Rdata ")Para formatos RDS, necesitas crear una variable para que tenga un nombre ( por decirlo asi ), te lo guarda para datos en bruto como si tubieras un excel sin nombre y necesitas cargarlo
orders <-readRDS("../../r-course-master/data/tema1/pedidos.rds")Como se guarada toda la sesión
#No es tan conveniente
save.image("../curso youtube/Trabajo.Rdata")Cargar datos que ya viene con R( precargados )
data("iris")
data("cars")
View(cars)Nota : el que guarda todo en archivo es el formato Rdata
Valiar que estoy sobre escribiendo variables
#Atach Cuando faltan datos dentro de un DF
library(readr)
datosfaltantes <-read.csv("../../r-course-master/data/tema1/missing-data.csv",na.strings = "")
#na,string : conviene los vacios datos los cambia a N/A
View(datosfaltantes)Para Limpiar los datos , donde no hay N/A: Te elimina estas filas donde hay N/A
datoslimpios<- na.omit(datosfaltantes)Eliminar N/A de filas de solo de una columna seleccionada
la columna Ingresos = Income
df.elim.igresos<- datosfaltantes[!is.na(datosfaltantes$Income),]
View(df.elim.igresos)Validar si tiene valores faltantes
complete.cases(datosfaltantes)## [1] TRUE TRUE TRUE FALSE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE
## [13] FALSE TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [25] TRUE TRUE TRUETodos los faltantes te los muestra como FALSE
Reemplazar datos que faltan a un N/A
Paso 1 cuando tiene ceros ,convertirlos en N/A
datosfaltantes1 <-read.csv("../../r-course-master/data/tema1/missing-data.csv",na.strings = "")
datosfaltantes1$Income[datosfaltantes1$Income==0]<- NA
datosfaltantes1$Income[datosfaltantes1$Phone_type==0]<- NAAjustar N/A o cero al promedio
Pero tiene N/A asi que nos va a amandar error , veamoslo
mean(datosfaltantes$Income)## [1] NAEntonces lo que le pediremos es : Borralos y te saldra
mean(datosfaltantes$Income,na.rm = T)## [1] 57872Desviacion estandar : eliminando los N/A
sd(datosfaltantes$Income,na.rm = T)## [1] 33170.09Para no perder esos datos , se reemplaza por el promedio (pero existe un algoritmo mejor )
#ingresaremos una columna donde los que faltan se vana ingresar los datos pendientes (N/A)
datosfaltantes$columna_Incertada <-
ifelse(is.na(datosfaltantes$Income),
mean(datosfaltantes$Income,na.rm= TRUE),datosfaltantes$Income)Algoritmo para agregar datos faltantes (Funcion parte 1)
Ahora si el N/A es una variable categorica ( en vez de numeros son letras ) Abusado que esta bien perro
pero es basada en distribuciones gussianas , es muy impotante que se limpie o se asigne valores
# X es un vector que contiene N/A
rand.impute <- function(x){ #x viene a ser la columna del data frame con la que trabajaremos
missing <- is.na(x) #missing va a ser un vector del mismo tamano que x pero que contiene TRUE or false dependiendo si es X es igual NA o no
n.missing <- sum(missing)# esta variable indica cuantos valores son NA dentro de x
# Nota : cuando se pide sumar valores booleanos solo suma los FALSOS
x.obs <- x[!missing]# x.obs son los valores x que NO tienen un valor N/A dentro de x
imputed <- x #imputed es el valor que vamos a devolver, le damos el valor x por defecto
imputed[missing] <- sample(x.obs, n.missing, replace = TRUE) # imputed en la posicion que falten los missing values, va a ser un muestreo aleatorio de las observaviones de x que no esten missing, de tama?o n.missing, para que sea representativa respecto los missing values, y con la possibilidad de remplazo para que cambie los NA por valores aleatorios
return(imputed) # devuelveme el valor imputed
}Funcion parte 2
rand.impute.data.frame <- function(df,cols){
names <- names(df) #damos al vector name, los nombres de las columnas
for (col in cols) {
name <- paste(names[col], ".imputed", sep = "") #Anadir el nombre imputed despues del column name para identificar las nuevas columnas
df[name] = rand.impute(df[,col])#llamaremos la funcion de sustitucion de valores con valores aeatorios
}
df #recordar devolver el data frame con el que estamos trabajando
}ahora aplicamos la funcion al Df para que agregue las columnas
Leyendo esta sentencia es: al objeto “datos faltantes” le vas a aplicar la funcion de agregar columnas en los dtaos faltantes de forma aleatorias al DF y entre parentesis le dices cuales columnas va a sumar (de la 1 a la 2)
Otra forma de explicarlo :ahora se manda llamar la funcion Ahora metemos en una variable llamada data ,las columnas del df “datos flatantes” + las variables categoricas y numericas faltantes distribucion gaussiana ( solo aplica para los que no tenia valores) los demas los deja igual
datosfaltantes1<-rand.impute.data.frame(datosfaltantes1,c(1,2))
head(datosfaltantes1)## Income Phone_type Car_type Income.imputed Phone_type.imputed
## 1 89800 Android Luxury 89800 Android
## 2 47500 Android Non-Luxury 47500 Android
## 3 45000 iPhone Luxury 45000 iPhone
## 4 44700 <NA> Luxury 44700 iPhone
## 5 59500 iPhone Luxury 59500 iPhone
## 6 NA Android Non-Luxury 83000 AndroidComo Borrar una Funcion, Variable , puntero y DF
# rm("Nombre de la funcion")Evitar duplicaciones de datos
Llenado de datos para crear un DF
family.salary = c(40000, 60000, 50000, 80000, 60000, 70000, 60000)
family.size = c(4, 3, 2, 2, 3, 4, 3)
family.car = c("Lujo", "Compacto", "Utilitario", "Lujo",
"Compacto", "Compacto", "Compacto")Creamos un DF
family <- data.frame(family.salary, family.size, family.car)La funcion unique borra los registros que estan duplicados y esto lo guardamos en otra variable
family.unique <- unique(family)Esta funcion nos manda una variable de tipo booleanos ( verdaero o falso ) solo manda los que estan duplicados el “original , no lo detecta”
duplicated(family)## [1] FALSE FALSE FALSE FALSE TRUE FALSE TRUESi quiero saber cuales son las filas duplicadas
family[duplicated(family),]## family.salary family.size family.car
## 5 60000 3 Compacto
## 7 60000 3 Compacto#Ojo: los corchetes le estan diciendo [dame esto,"todo"( es dejarlo en blanco)]Reescalar datos
Primero se instala un paquete ( Aqui esta comentado por que solo se tiene que instalar una ves recoordemos)
#install.packages("scales")Acivar el paquete y cargamos los datos
library("scales")##
## Attaching package: 'scales'## The following object is masked from 'package:purrr':
##
## discard## The following object is masked from 'package:readr':
##
## col_factorstudents <-read.csv( "../../r-course-master/data/tema1/data-conversion.csv")
head(students)## Age State Gender Height Income
## 1 23 NJ F 61 5000
## 2 13 NY M 55 1000
## 3 36 NJ M 66 3000
## 4 31 VA F 64 4000
## 5 58 NY F 70 30000
## 6 29 TX F 63 10000Ceramos una nueva variable reescalar es una transformacion lineal
students$Income.rescaled<- rescale(students$Income)
head(students)## Age State Gender Height Income Income.rescaled
## 1 23 NJ F 61 5000 0.07407407
## 2 13 NY M 55 1000 0.00000000
## 3 36 NJ M 66 3000 0.03703704
## 4 31 VA F 64 4000 0.05555556
## 5 58 NY F 70 30000 0.53703704
## 6 29 TX F 63 10000 0.16666667La funcion rescale es un resumen de la formula los resultados son exactamente lo mismo
(students$Income - min(students$Income))/(max(students$Income)-min(students$Income))## [1] 0.07407407 0.00000000 0.03703704 0.05555556 0.53703704 0.16666667
## [7] 0.90740741 1.00000000 0.01851852 0.35185185En el caso que esten en valores porcentuales aqui lo que se hace es agregar “to=c(min, max)”
rescale(students$Income, to= c(0,100))## [1] 7.407407 0.000000 3.703704 5.555556 53.703704 16.666667
## [7] 90.740741 100.000000 1.851852 35.185185# El resultado es el mismo pero multiplicado por 100 en este caso Re-esclado para varias columnas/Variables
rescaled.many<-function(dataframe ,cols){
names<- names(dataframe)
for (col in cols){
name<- paste(names[col],"rescaled",sep = ".")
dataframe[name]<- rescale(dataframe[,col])
}
cat(paste("Hemos reescalado", length(cols),"variable(s)"))
dataframe
}Se obtiene para cada variable un reescalamiento
students<- rescaled.many(students,c(1,4))## Hemos reescalado 2 variable(s)head(students)## Age State Gender Height Income Income.rescaled Age.rescaled Height.rescaled
## 1 23 NJ F 61 5000 0.07407407 0.2222222 0.4000000
## 2 13 NY M 55 1000 0.00000000 0.0000000 0.0000000
## 3 36 NJ M 66 3000 0.03703704 0.5111111 0.7333333
## 4 31 VA F 64 4000 0.05555556 0.4000000 0.6000000
## 5 58 NY F 70 30000 0.53703704 1.0000000 1.0000000
## 6 29 TX F 63 10000 0.16666667 0.3555556 0.5333333Normalizando y Estandarizando un DF
Normalizar (campana de Gauss):Hacer que algo se ajuste a una norma, una regla o un modelo común Nota Importante : esta tecnica solo sirve cuando el DF contiene solo Numeros
housing<-read.csv("../../r-course-master/data/tema1/BostonHousing.csv")
housing.z<- scale(housing)Ente mas cerca del cero :son los mas cencanos a la media Ente mas lejos del cero :son mas alejados de la media
Si quisiera ajustar la media ( que no se calculara con respecto a esta)
housing.mean<- scale(housing,center = T,scale = F)
#scale es la desviacion estandar Ahora hagamos una funcion que agrege esa normalizacion y no sobre escriba el DF ( para conservar los datos )
scale.many <- function(dataframe,cols){
names <- names(dataframe)
for (col in cols){
name<-paste(names[col],"z",sep = ".")
dataframe[name]<- scale(dataframe[,cols])
}
cat(paste("hemos normalizado",length(cols),"variables"))
dataframe
}Ahora ejecutaremos la funcion
housing <- scale.many(housing,c(1,3,5:8))## hemos normalizado 6 variablesCategorizando Informacion numerica
Ejemplo:Meter numeros en categorias ( 10 = pobre y 100 = rico )
students3 <-read.csv("../../r-course-master/data/tema1/data-conversion.csv")
bp<-c(-Inf,10000,31000,Inf)
names<-c("low","averange", "High")
#siempre es una categoria menos que los parametros y hace intervalos de tipo (abierto y cerrado]
# Inf a 10000, 10000 a 31000,31000 a Inf
students3$Income.cat<- cut(students3$Income,breaks=bp,labels=names)
head(students3)## Age State Gender Height Income Income.cat
## 1 23 NJ F 61 5000 low
## 2 13 NY M 55 1000 low
## 3 36 NJ M 66 3000 low
## 4 31 VA F 64 4000 low
## 5 58 NY F 70 30000 averange
## 6 29 TX F 63 10000 lowTambien puedo decirle que no ponga etiquetas (no las defino yo) observa como sumo una columna con los valores en intervalo y no en etiquetas
students3$Income.cat2<-cut(students3$Income,breaks=bp)
head(students3)## Age State Gender Height Income Income.cat Income.cat2
## 1 23 NJ F 61 5000 low (-Inf,1e+04]
## 2 13 NY M 55 1000 low (-Inf,1e+04]
## 3 36 NJ M 66 3000 low (-Inf,1e+04]
## 4 31 VA F 64 4000 low (-Inf,1e+04]
## 5 58 NY F 70 30000 averange (1e+04,3.1e+04]
## 6 29 TX F 63 10000 low (-Inf,1e+04]Tambien puedo pedirle a R corte como el quiera solo diciendo en cuantas categorias cortar
students3$Income.cat3<-cut(students3$Income,breaks=4,labels = c("1","2","3","4"))
head(students3)## Age State Gender Height Income Income.cat Income.cat2 Income.cat3
## 1 23 NJ F 61 5000 low (-Inf,1e+04] 1
## 2 13 NY M 55 1000 low (-Inf,1e+04] 1
## 3 36 NJ M 66 3000 low (-Inf,1e+04] 1
## 4 31 VA F 64 4000 low (-Inf,1e+04] 1
## 5 58 NY F 70 30000 averange (1e+04,3.1e+04] 3
## 6 29 TX F 63 10000 low (-Inf,1e+04] 1Variables ficticias
Instalamos la libreria
#install.packages("dummies")
library("dummies")## dummies-1.5.6 provided by Decision PatternsComo una conversion de variables a todo el Df
students4 <-read.csv("../../r-course-master/data/tema1/data-conversion.csv")
students.dummy<- dummy.data.frame(students4,sep=".", all = T)## Warning in model.matrix.default(~x - 1, model.frame(~x - 1), contrasts = FALSE):
## non-list contrasts argument ignored
## Warning in model.matrix.default(~x - 1, model.frame(~x - 1), contrasts = FALSE):
## non-list contrasts argument ignorednames(students.dummy)## [1] "Age" "State.NJ" "State.NY" "State.TX" "State.VA" "Gender.F" "Gender.M"
## [8] "Height" "Income"#all detecta los valores numericos y los respeta , no los hace como los estados Como una conversion de variables a una columna ( femenino y Masculino )
students4$nuevo <- dummy(students4$Gender)## Warning in model.matrix.default(~x - 1, model.frame(~x - 1), contrasts = FALSE):
## non-list contrasts argument ignoredSelecionamos los datos que queremos aplicarle la funcion dummy a algunas columas en especial
students4<- dummy.data.frame(students4,names = c("State"))## Warning in model.matrix.default(~x - 1, model.frame(~x - 1), contrasts = FALSE):
## non-list contrasts argument ignored#Este ejemplo lo que hace es convierte una dummy en varias columans
students4$Gender.dummy<-dummy(students4$Gender)## Warning in model.matrix.default(~x - 1, model.frame(~x - 1), contrasts = FALSE):
## non-list contrasts argument ignored#cuando solo tienes 2 ( hombre y mujer ) te crea una columna y no hace 2 columnas : asigna para F=1 y M=0Formas de Eliminar información que falta
bostonprincipal<- read.csv("../../r-course-master/data/tema1/housing-with-missing-value.csv",header = T,stringsAsFactors = F)Si queremos una vista rapida de los datos ( estadiasticos basicos ) Enterga :Minimo,Primer cuartil ,mediana,media,3er quartil, maximo y numero de N/A por columna
summary(bostonprincipal)## X crim zn indus
## Min. : 1.0 Min. : 0.00632 Min. : 0.00 Min. : 0.46
## 1st Qu.:127.2 1st Qu.: 0.08205 1st Qu.: 0.00 1st Qu.: 5.19
## Median :253.5 Median : 0.25651 Median : 0.00 Median : 9.69
## Mean :253.5 Mean : 3.61352 Mean : 11.36 Mean :11.14
## 3rd Qu.:379.8 3rd Qu.: 3.67708 3rd Qu.: 12.50 3rd Qu.:18.10
## Max. :506.0 Max. :88.97620 Max. :100.00 Max. :27.74
##
## chas nox rm age
## Min. :0.00000 Min. :0.3850 Min. :3.561 Min. : 2.90
## 1st Qu.:0.00000 1st Qu.:0.4490 1st Qu.:5.886 1st Qu.: 45.02
## Median :0.00000 Median :0.5380 Median :6.208 Median : 77.50
## Mean :0.06917 Mean :0.5547 Mean :6.285 Mean : 68.57
## 3rd Qu.:0.00000 3rd Qu.:0.6240 3rd Qu.:6.623 3rd Qu.: 94.08
## Max. :1.00000 Max. :0.8710 Max. :8.780 Max. :100.00
##
## dis rad tax ptratio
## Min. : 1.130 Min. : 1.000 Min. :187.0 Min. :12.60
## 1st Qu.: 2.100 1st Qu.: 4.000 1st Qu.:279.0 1st Qu.:17.40
## Median : 3.207 Median : 5.000 Median :330.0 Median :19.10
## Mean : 3.795 Mean : 9.515 Mean :408.2 Mean :18.47
## 3rd Qu.: 5.188 3rd Qu.:24.000 3rd Qu.:666.0 3rd Qu.:20.20
## Max. :12.127 Max. :24.000 Max. :711.0 Max. :22.00
## NA's :40 NA's :40
## b lstat medv
## Min. : 0.32 Min. : 1.73 Min. : 5.00
## 1st Qu.:375.38 1st Qu.: 6.95 1st Qu.:17.02
## Median :391.44 Median :11.36 Median :21.20
## Mean :356.67 Mean :12.65 Mean :22.53
## 3rd Qu.:396.23 3rd Qu.:16.95 3rd Qu.:25.00
## Max. :396.90 Max. :37.97 Max. :50.00
## Primera solucion n.omit ( eliminar las filas ) Pero esto borra datos que pueden ser muy importantes para el analsis y se mueven algunos datos
boston1<- bostonprincipal
#na.omit(boston1)Puedes compararlos
summary(boston1)## X crim zn indus
## Min. : 1.0 Min. : 0.00632 Min. : 0.00 Min. : 0.46
## 1st Qu.:127.2 1st Qu.: 0.08205 1st Qu.: 0.00 1st Qu.: 5.19
## Median :253.5 Median : 0.25651 Median : 0.00 Median : 9.69
## Mean :253.5 Mean : 3.61352 Mean : 11.36 Mean :11.14
## 3rd Qu.:379.8 3rd Qu.: 3.67708 3rd Qu.: 12.50 3rd Qu.:18.10
## Max. :506.0 Max. :88.97620 Max. :100.00 Max. :27.74
##
## chas nox rm age
## Min. :0.00000 Min. :0.3850 Min. :3.561 Min. : 2.90
## 1st Qu.:0.00000 1st Qu.:0.4490 1st Qu.:5.886 1st Qu.: 45.02
## Median :0.00000 Median :0.5380 Median :6.208 Median : 77.50
## Mean :0.06917 Mean :0.5547 Mean :6.285 Mean : 68.57
## 3rd Qu.:0.00000 3rd Qu.:0.6240 3rd Qu.:6.623 3rd Qu.: 94.08
## Max. :1.00000 Max. :0.8710 Max. :8.780 Max. :100.00
##
## dis rad tax ptratio
## Min. : 1.130 Min. : 1.000 Min. :187.0 Min. :12.60
## 1st Qu.: 2.100 1st Qu.: 4.000 1st Qu.:279.0 1st Qu.:17.40
## Median : 3.207 Median : 5.000 Median :330.0 Median :19.10
## Mean : 3.795 Mean : 9.515 Mean :408.2 Mean :18.47
## 3rd Qu.: 5.188 3rd Qu.:24.000 3rd Qu.:666.0 3rd Qu.:20.20
## Max. :12.127 Max. :24.000 Max. :711.0 Max. :22.00
## NA's :40 NA's :40
## b lstat medv
## Min. : 0.32 Min. : 1.73 Min. : 5.00
## 1st Qu.:375.38 1st Qu.: 6.95 1st Qu.:17.02
## Median :391.44 Median :11.36 Median :21.20
## Mean :356.67 Mean :12.65 Mean :22.53
## 3rd Qu.:396.23 3rd Qu.:16.95 3rd Qu.:25.00
## Max. :396.90 Max. :37.97 Max. :50.00
## Segunda solucion: para eliminar solo de alguna columna
drop_na<-c("rad")
boston2<- bostonprincipal[complete.cases(bostonprincipal[,!names(bostonprincipal)%in% drop_na]),]Elimianr columna de un DF
Esto modifica el DF, elimina la columna
bostonsin1col<- bostonprincipal
bostonsin1col$rad=NULL Eliminar diferentes columnas a la vez
eliminarcol<-c("rad","ptratio")
bostonsin2col<- bostonprincipal[,!names(bostonprincipal)%in% eliminarcol]Formas de completar información que falta
#install.packages("Hmisc")Activamos la libreria y cargamos los datos
library("Hmisc")## Loading required package: lattice## Loading required package: survival## Loading required package: Formula##
## Attaching package: 'Hmisc'## The following objects are masked from 'package:dplyr':
##
## src, summarize## The following objects are masked from 'package:base':
##
## format.pval, unitsbostoncopy<- bostonprincipalImpute :Es para rellenar valores NA con la media =mean
bostoncopy$rad<-impute(bostoncopy$rad,mean)
bostoncopy$ptratio<-impute(bostoncopy$ptratio,mean)
summary(bostoncopy)##
## 40 values imputed to 9.515021
##
##
## 40 values imputed to 18.4676## X crim zn indus
## Min. : 1.0 Min. : 0.00632 Min. : 0.00 Min. : 0.46
## 1st Qu.:127.2 1st Qu.: 0.08205 1st Qu.: 0.00 1st Qu.: 5.19
## Median :253.5 Median : 0.25651 Median : 0.00 Median : 9.69
## Mean :253.5 Mean : 3.61352 Mean : 11.36 Mean :11.14
## 3rd Qu.:379.8 3rd Qu.: 3.67708 3rd Qu.: 12.50 3rd Qu.:18.10
## Max. :506.0 Max. :88.97620 Max. :100.00 Max. :27.74
## chas nox rm age
## Min. :0.00000 Min. :0.3850 Min. :3.561 Min. : 2.90
## 1st Qu.:0.00000 1st Qu.:0.4490 1st Qu.:5.886 1st Qu.: 45.02
## Median :0.00000 Median :0.5380 Median :6.208 Median : 77.50
## Mean :0.06917 Mean :0.5547 Mean :6.285 Mean : 68.57
## 3rd Qu.:0.00000 3rd Qu.:0.6240 3rd Qu.:6.623 3rd Qu.: 94.08
## Max. :1.00000 Max. :0.8710 Max. :8.780 Max. :100.00
## dis rad tax ptratio
## Min. : 1.130 Min. : 1.000 Min. :187.0 Min. :12.60
## 1st Qu.: 2.100 1st Qu.: 4.000 1st Qu.:279.0 1st Qu.:17.40
## Median : 3.207 Median : 5.000 Median :330.0 Median :18.60
## Mean : 3.795 Mean : 9.515 Mean :408.2 Mean :18.47
## 3rd Qu.: 5.188 3rd Qu.: 9.515 3rd Qu.:666.0 3rd Qu.:20.20
## Max. :12.127 Max. :24.000 Max. :711.0 Max. :22.00
## b lstat medv
## Min. : 0.32 Min. : 1.73 Min. : 5.00
## 1st Qu.:375.38 1st Qu.: 6.95 1st Qu.:17.02
## Median :391.44 Median :11.36 Median :21.20
## Mean :356.67 Mean :12.65 Mean :22.53
## 3rd Qu.:396.23 3rd Qu.:16.95 3rd Qu.:25.00
## Max. :396.90 Max. :37.97 Max. :50.00Nota : ya no tiene N/A
Para este caso lo que se va a reemplazar con la mediana
bostonmedian<-bostonprincipal
bostonmedian$rad<-impute(bostonmedian$rad,median)
bostonmedian$ptratio<-impute(bostonmedian$ptratio,median)Nota Importante : se puede usar con la media , mediana , un numero etc,etc ,etc
Matriz donde vemos los campos vacios
#install.packages("mice")
#install.packages("VIM")
library(mice)##
## Attaching package: 'mice'## The following object is masked from 'package:stats':
##
## filter## The following objects are masked from 'package:base':
##
## cbind, rbindlibrary(VIM)## Loading required package: colorspace## Loading required package: grid## VIM is ready to use.## Suggestions and bug-reports can be submitted at: https://github.com/statistikat/VIM/issues##
## Attaching package: 'VIM'## The following object is masked from 'package:datasets':
##
## sleepEntrega una matriz donde te muestra los datos completos y los casos donde faltan valores
md.pattern(bostonprincipal)
## X crim zn indus chas nox rm age dis tax b lstat medv rad ptratio
## 431 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0
## 35 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1
## 35 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1
## 5 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 2
## 0 0 0 0 0 0 0 0 0 0 0 0 0 40 40 80Grafiquemoslo Col= color labels= pone el nombre de la columna cex.axis=Tamaño de letra Gap= distancia entre graficos ylab=Nombre de la grafica en el eje y
aggr(bostonprincipal,col=c('blue','red',numbers=T,shortVar=T),labels=names(bostonprincipal),cey.axis=0.5,cex.axis=0.6, gap=3,ylab=c("histograma de NA","Matriz"))
# Combinando y Separando datos Lo interesante de este tema es lo que tiene valor Cargamos libreria tidyr
#install.packages("tidyr")
library(tidyr)Carga de datos
crime.data <- read.csv("../../r-course-master/data/tema1/USArrests.csv", stringsAsFactors = F)
head(crime.data)## X Murder Assault UrbanPop Rape
## 1 Alabama 13.2 236 58 21.2
## 2 Alaska 10.0 263 48 44.5
## 3 Arizona 8.1 294 80 31.0
## 4 Arkansas 8.8 190 50 19.5
## 5 California 9.0 276 91 40.6
## 6 Colorado 7.9 204 78 38.7Como Ingresar una columna
crime.data<-cbind(state=rownames(crime.data),crime.data)Rownames ( enumera las filas) # Como se pueden resumir datos Array asociativo /diccionario gater = Función que unifica columnas Tablas de doble entrada
crime.data1<- gather(crime.data,key="crime_type", value = "arrest_estimate",Murder:UrbanPop)Resume en 2 variables
crime.data2<- gather(crime.data,key="crime_type", value = "arrest_estimate",-state)lo contrario a gether
crime.data3<-spread(crime.data2,key="crime_type", value="arrest_estimate")Es como Juntar valores
crime.data4<- unite(crime.data,col="Mordet_Assailt",Murder,Assault,sep="_")Ahora si te pasan datos que estan juntos con caracteres se separan con la funcion separate
crime.data5<-separate(crime.data4,col="Mordet_Assailt",
into = c("Murder","Assault"),
sep = "_")Modelos predictivos para imputar datos
Paso 1:pedrice valores que se parezcan con un metodo Paso 2:acompletar los datos paso 3:Sustituir en el Df original los datos completos
Borrarlos no es una opción Tecnicas de imputación multivariable
library(mice)
housingdata <- read.csv("../../r-course-master/data/tema1/housing-with-missing-value.csv",header = T,stringsAsFactors = F)columns<-c("ptratio", "rad")
imputed_data<-mice(housingdata[,names(housingdata)%in% columns],
m=5,# numero de imputaciones
maxit = 3,#Numero de iteraciones
method="pmm",#Predice el promedio de los datos
seed=2018#Semilla
)##
## iter imp variable
## 1 1 rad ptratio
## 1 2 rad ptratio
## 1 3 rad ptratio
## 1 4 rad ptratio
## 1 5 rad ptratio
## 2 1 rad ptratio
## 2 2 rad ptratio
## 2 3 rad ptratio
## 2 4 rad ptratio
## 2 5 rad ptratio
## 3 1 rad ptratio
## 3 2 rad ptratio
## 3 3 rad ptratio
## 3 4 rad ptratio
## 3 5 rad ptratioMetodos =method pmm= comparación predictiva de medios logreg =regresión logistica “solo se pude con 2” polyreg = regresión logistica politomica polr= modelo de probabilidades proporcionales
summary(imputed_data)## Class: mids
## Number of multiple imputations: 5
## Imputation methods:
## rad ptratio
## "pmm" "pmm"
## PredictorMatrix:
## rad ptratio
## rad 0 1
## ptratio 1 0Ahora con el modelo hecho acompletaremos los datos
complete.cases<-mice::complete(imputed_data)asi se indica cuando los metodos se confunden o R o entiende de que libreria (mice::complete)la libreria + el metodo Ahora se sustituye los valores na del Df original por los que estimaste con el metodo mice en otras palabras: mete en esta columna lo que esta en esta columna
housingdata$ptratio<-complete.cases$ptratio
housingdata$rad<-complete.cases$radValidamos cuantos NA existen
anyNA(housingdata)## [1] FALSEOtra forma de inferir los NA
library(Hmisc)
housingdata1 <- read.csv("../../r-course-master/data/tema1/housing-with-missing-value.csv",header = T,stringsAsFactors = F)
impute_arg<- aregImpute(~ptratio+rad,data = housingdata1,n.impute = 5)## Iteration 1
Iteration 2
Iteration 3
Iteration 4
Iteration 5
Iteration 6
Iteration 7
Iteration 8 Veamos que nos arroja
impute_arg##
## Multiple Imputation using Bootstrap and PMM
##
## aregImpute(formula = ~ptratio + rad, data = housingdata1, n.impute = 5)
##
## n: 506 p: 2 Imputations: 5 nk: 3
##
## Number of NAs:
## ptratio rad
## 40 40
##
## type d.f.
## ptratio s 2
## rad s 1
##
## Transformation of Target Variables Forced to be Linear
##
## R-squares for Predicting Non-Missing Values for Each Variable
## Using Last Imputations of Predictors
## ptratio rad
## 0.244 0.263con este metodo de donde obtengo los valores que quiero imputar
impute_arg$imputed$ptratio## [,1] [,2] [,3] [,4] [,5]
## 43 18.6 16.9 21.0 18.5 18.8
## 61 18.9 17.9 19.1 19.2 20.9
## 99 18.4 17.8 18.2 18.7 19.7
## 103 14.7 21.0 17.9 21.2 14.7
## 104 17.6 21.0 18.7 21.0 16.6
## 113 19.2 18.7 16.9 14.7 14.7
## 117 16.4 21.2 18.6 20.2 20.9
## 118 20.9 21.0 17.9 21.0 16.1
## 125 19.1 17.9 15.3 18.6 12.6
## 136 21.2 21.0 21.0 21.0 17.8
## 153 14.7 18.6 17.9 20.9 21.0
## 156 18.4 18.6 19.1 18.6 17.8
## 164 18.7 18.7 17.8 16.8 14.7
## 168 15.2 15.2 18.0 19.2 20.1
## 173 15.2 21.0 17.8 14.7 17.9
## 197 15.2 21.0 18.0 17.6 14.7
## 214 17.6 16.6 18.0 18.6 20.1
## 215 14.7 17.6 19.1 17.6 18.6
## 219 14.7 16.0 17.8 18.9 14.7
## 226 17.8 21.0 18.4 14.7 19.1
## 247 19.7 19.0 17.8 20.9 19.7
## 272 18.3 18.6 19.1 17.6 17.6
## 291 21.0 21.0 18.8 17.0 17.6
## 292 16.9 22.0 17.4 17.6 17.9
## 295 18.3 14.4 16.4 21.2 18.5
## 303 16.1 19.2 21.0 13.0 16.6
## 313 19.2 19.2 18.0 21.2 14.4
## 314 21.2 19.2 18.0 18.6 18.5
## 327 13.0 18.3 19.1 17.0 19.2
## 332 20.2 18.0 13.6 17.8 21.0
## 364 20.2 20.2 20.2 20.2 20.2
## 391 20.2 20.2 20.2 20.2 20.2
## 392 20.2 20.2 20.2 20.2 20.2
## 403 20.2 20.2 20.2 20.2 20.2
## 417 20.2 20.2 20.2 20.2 20.2
## 437 17.8 16.9 17.9 12.6 18.6
## 442 21.0 21.0 17.8 13.0 15.3
## 458 20.2 20.2 20.2 20.2 20.2
## 465 20.2 20.2 20.2 20.2 20.2
## 500 17.4 19.6 17.8 18.7 18.7impute_arg$imputed$rad## [,1] [,2] [,3] [,4] [,5]
## 29 4 5 1 4 1
## 61 2 4 6 5 5
## 83 7 2 7 2 5
## 85 4 4 4 4 2
## 101 5 4 3 4 5
## 131 4 8 4 4 4
## 133 4 5 4 4 4
## 139 4 4 3 5 4
## 156 3 4 4 3 6
## 166 5 5 5 5 3
## 175 1 5 4 6 6
## 176 4 1 8 8 5
## 185 6 3 5 2 8
## 197 5 4 4 6 5
## 203 5 4 5 8 5
## 233 2 4 8 3 7
## 236 8 2 3 3 5
## 243 5 5 6 8 5
## 261 5 5 5 5 4
## 263 5 5 1 5 7
## 273 3 4 4 5 4
## 279 6 4 6 3 6
## 296 6 5 4 1 6
## 310 3 4 3 4 3
## 329 4 8 8 4 1
## 337 24 24 24 24 24
## 345 3 3 4 3 4
## 362 24 24 24 24 24
## 370 24 24 24 24 24
## 376 24 24 24 24 24
## 407 24 24 24 24 24
## 419 24 24 24 24 24
## 423 5 24 24 24 24
## 437 3 5 2 5 4
## 442 4 4 6 2 5
## 452 24 24 24 24 24
## 464 24 24 24 24 24
## 486 24 24 24 24 5
## 491 24 24 24 24 24
## 497 4 5 7 5 8Solo falatria sustituirlo con las columnas que se tienen en el DF Original
Detección de outliners atravez de box pot
Los valores que se salen de lo comun es decir (1 o 2 y aparece un 73 de la nada ) Cambiar ruta del archivo
ozone.data<- read.csv("../../r-course-master/data/tema1/ozone.csv",header = T,stringsAsFactors = F)Nota importante : antes de esimar los NA seria importente detectar esto para que los metodos no muevan el resultado
outliner_values<-boxplot(ozone.data$pressure_height,main="presure Heigh",boxwex=0.5)
summary(ozone.data$pressure_height)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 5320 5700 5770 5753 5830 5950 12Comparando los datos
boxplot(pressure_height~Month,
data = ozone.data,
main="presure Height per Month"
)
boxplot(ozone_reading~Month,
data = ozone.data,
main="ozone reading per onth"
)$out # el out muestra un mensaje en consola ## [1] 11.06 9.93 22.89 24.29 29.79mtext("hola")
# Enmascarando Outlinerscon trasnformaciones y cappings Lo que hacemos es los que esten arriba del percentil 95 se sustitute con la media , asi los datos no se ven tan alejados
impute_outliners<- function(x,removeNA=T){
quantiles<- quantile(x,c(0.05,0.95),na.rm=removeNA)
x[x<quantiles[1]]<- mean(x,na.rm=removeNA)
x[x>quantiles[2]]<- mean(x,na.rm=removeNA)
x
}Validemos si acomodo los datos que estaban separados Aplicamos la función
imputed.data.ozone<-impute_outliners(ozone.data$pressure_height)Ordenar 2 graficos
par(mfrow=c(1,2))
boxplot(ozone.data$pressure_height,main="presión con Outliners")
boxplot(imputed.data.ozone,main="presión sin Outliners")
# Outliners , pero reemplazando datos Nota: es diferente imputar los datos a reemplazarlos
replace.outliners<- function(x,removeNA=T){
qrts<-quantile(x,probs=c(0.25,0.75),na.rm=removeNA)
caps<-quantile(x,probs=c(0.05,0.95),na.rm=removeNA)
iqr<-qrts[2]-qrts[1]#Rango Intercuartilico
h<- 1.5*iqr
x[x<qrts[2]-h]<-caps[1]
x[x<qrts[1]+h]<-caps[2]
x
}Ahora esa selección de valores fuera de esos rangos , se van a sutituir con la función de arriba ya cargada
capped_presure_heigh<-replace.outliners(ozone.data$pressure_height)Grafiquemoslo 5
par(mfrow=c(1,2))
boxplot(ozone.data$pressure_height,main="presión con Outliners")
boxplot(capped_presure_heigh,main="Presión sin Outliners")
# Resumen de datos y estructura del DF
estructura<-read.csv("../../r-course-master/data/tema2/auto-mpg.csv",header = T,stringsAsFactors = F)
anyNA(estructura)#validamos si tiene NA's el DF## [1] FALSELos cilindros de los carros no nos sirven como valores , por eso los vamos a comvertir en factores
estructura$cylinders<-factor(estructura$cylinders,
levels= c(3,4,6,8),
labels= c("3cil","4cil","6cil","8cil"))summary da valores estadisticos
summary(estructura)## No mpg cylinders displacement horsepower
## Min. : 1.0 Min. : 9.00 3cil: 4 Min. : 68.0 Min. : 46.0
## 1st Qu.:100.2 1st Qu.:17.50 4cil:204 1st Qu.:104.2 1st Qu.: 76.0
## Median :199.5 Median :23.00 6cil: 84 Median :148.5 Median : 92.0
## Mean :199.5 Mean :23.51 8cil:103 Mean :193.4 Mean :104.1
## 3rd Qu.:298.8 3rd Qu.:29.00 NA's: 3 3rd Qu.:262.0 3rd Qu.:125.0
## Max. :398.0 Max. :46.60 Max. :455.0 Max. :230.0
## weight acceleration model_year car_name
## Min. :1613 Min. : 8.00 Min. :70.00 Length:398
## 1st Qu.:2224 1st Qu.:13.82 1st Qu.:73.00 Class :character
## Median :2804 Median :15.50 Median :76.00 Mode :character
## Mean :2970 Mean :15.57 Mean :76.01
## 3rd Qu.:3608 3rd Qu.:17.18 3rd Qu.:79.00
## Max. :5140 Max. :24.80 Max. :82.00Aqui lo que vas a ver es como R esta leyendo el archivo previamnte cargado
str(estructura)## 'data.frame': 398 obs. of 9 variables:
## $ No : int 1 2 3 4 5 6 7 8 9 10 ...
## $ mpg : num 28 19 36 28 21 23 15.5 32.9 16 13 ...
## $ cylinders : Factor w/ 4 levels "3cil","4cil",..: 2 1 2 2 3 2 4 2 3 4 ...
## $ displacement: num 140 70 107 97 199 115 304 119 250 318 ...
## $ horsepower : int 90 97 75 92 90 95 120 100 105 150 ...
## $ weight : int 2264 2330 2205 2288 2648 2694 3962 2615 3897 3755 ...
## $ acceleration: num 15.5 13.5 14.5 17 15 15 13.9 14.8 18.5 14 ...
## $ model_year : int 71 72 82 72 70 75 76 81 75 76 ...
## $ car_name : chr "chevrolet vega 2300" "mazda rx2 coupe" "honda accord" "datsun 510 (sw)" ...Tambien puedes verlo en un zoom de cada columna lo que vas a ver es: cuantos niveles ( etiquetas y que posicion con respecto a esa etiqueta tiene ) 2=4cil ,1=4cil y asi sucesivamente
str(estructura$cylinders)## Factor w/ 4 levels "3cil","4cil",..: 2 1 2 2 3 2 4 2 3 4 ...Estadisticos y Medidas(solo valores numericos)
#install.packages(c("modeest","raster","moments"))
library(modeest)## Registered S3 method overwritten by 'rmutil':
## method from
## print.response httrlibrary(raster)## Loading required package: sp##
## Attaching package: 'raster'## The following object is masked from 'package:dplyr':
##
## selectlibrary(moments)##
## Attaching package: 'moments'## The following object is masked from 'package:modeest':
##
## skewnessCreamos un vector desde una DF
x<-estructura$mpgmean(x)#promedio (medida aritmetica)## [1] 23.51457median(x)#mediana ## [1] 23mfv(x)#moda es el valor que mas se repite ## [1] 13quantile(x)#percentiles## 0% 25% 50% 75% 100%
## 9.0 17.5 23.0 29.0 46.6Medidas de Dispercion
como se distribuyen los datos (que tan lejos estan de la media)
var(x)#Varianza## [1] 61.08961sd(x)#desviacion tipica(se usa intervalos de confianza )## [1] 7.815984cv(x)#Coheficiente de variacion ## [1] 33.2389Momento de # Asimetria de Fisher Asimetria de fisher(que tan chueca esta la campana de gauss), si sale negativo es asimetrica negativa , y si es positiva , pues es asimetrica positiva
skewness(x)## [1] 0.4553419Curtosos
Que tan chaparra esta la campana de gauss Leptocurtica = sin colas ,“una curva flaquita” Mesocurtica=Esta es perfecta platicurtica=La curva es casi toda curvas
moments::kurtosis(x)## [1] 2.480575#Histograma
hist(x)
# Obtener un subconnjunto de los datos Cargamos los datos
dat<- read.csv("../../r-course-master/data/tema2/auto-mpg.csv",header = T,stringsAsFactors = F)Indice por posicion
dat[1:5,8:9]## model_year car_name
## 1 71 chevrolet vega 2300
## 2 72 mazda rx2 coupe
## 3 82 honda accord
## 4 72 datsun 510 (sw)
## 5 70 amc gremlinind.posic<-dat[c(1:5,c(8,9))]
ind.posic[1:3,]## No mpg cylinders displacement horsepower model_year car_name
## 1 1 28 4 140 90 71 chevrolet vega 2300
## 2 2 19 3 70 97 72 mazda rx2 coupe
## 3 3 36 4 107 75 82 honda accordO tienes la opcion de seleccionar por nombre de la columna: Se lle mas o menos asi :De la fila 1 a la 5 de las columnas que se llaman “etc”
dat[1:5,c("model_year","car_name")]## model_year car_name
## 1 71 chevrolet vega 2300
## 2 72 mazda rx2 coupe
## 3 82 honda accord
## 4 72 datsun 510 (sw)
## 5 70 amc gremlinObtener el Minimo y maximo de una fila
dat[dat$mpg==max(dat$mpg)|dat$mpg==min(dat$mpg),]## No mpg cylinders displacement horsepower weight acceleration model_year
## 190 190 9.0 8 304 193 4732 18.5 70
## 269 269 46.6 4 86 65 2110 17.9 80
## car_name
## 190 hi 1200d
## 269 mazda glcObtener valores con rangos
dat[dat$mpg<30 & dat$cylinders==6,c("car_name","mpg")]## car_name mpg
## 5 amc gremlin 21.0
## 9 chevroelt chevelle malibu 16.0
## 21 plymouth volare 20.5
## 29 amc hornet 18.0
## 37 amc concord 19.4
## 41 dodge aspen 6 20.6
## 43 ford mustang 18.0
## 48 dodge aspen 19.1
## 50 chevrolet citation 23.5
## 52 ford torino 500 19.0
## 69 chevrolet chevelle malibu classic 16.0
## 72 ford maverick 15.0
## 73 pontiac firebird 19.0
## 75 chevrolet nova 18.0
## 77 chevrolet malibu 20.5
## 85 ford maverick 21.0
## 110 plymouth valiant 18.0
## 115 buick century 22.4
## 117 pontiac lemans v6 21.5
## 124 ford pinto 18.0
## 137 ford fairmont (auto) 20.2
## 147 mercury capri v6 21.0
## 149 amc gremlin 18.0
## 157 plymouth satellite custom 16.0
## 161 datsun 810 maxima 24.2
## 174 buick skyhawk 21.0
## 177 amc pacer d/l 17.5
## 178 ford granada l 22.0
## 182 plymouth duster 23.0
## 183 chevrolet nova 15.0
## 184 chevrolet nova 22.0
## 201 amc hornet 19.0
## 202 toyota cressida 25.4
## 203 ford maverick 24.0
## 207 plymouth satellite sebring 18.0
## 209 amc hornet sportabout (sw) 18.0
## 215 buick century special 20.6
## 219 amc matador 18.0
## 223 amc gremlin 20.0
## 225 oldsmobile omega brougham 26.8
## 238 dodge aspen 18.6
## 239 mercury zephyr 6 19.8
## 247 plymouth duster 20.0
## 258 peugeot 604sl 16.2
## 270 buick skylark 20.5
## 273 ford granada 18.5
## 274 buick regal sport coupe (turbo) 17.7
## 276 ford granada gl 20.2
## 277 buick century 17.0
## 282 chevrolet concours 17.5
## 284 amc matador 15.0
## 287 amc concord dl 6 20.2
## 292 mercury monarch 15.0
## 294 toyota mark ii 20.0
## 305 chrysler lebaron salon 17.6
## 308 plymouth valiant 22.0
## 310 pontiac ventura sj 18.5
## 312 toyota mark ii 19.0
## 315 ford granada ghia 18.0
## 324 ford maverick 18.0
## 325 amc hornet 18.0
## 333 ford maverick 21.0
## 334 plymouth volare custom 19.0
## 336 datsun 810 22.0
## 337 amc pacer 19.0
## 341 pontiac phoenix lj 19.2
## 343 chevrolet chevelle malibu 17.0
## 353 chevrolet citation 28.8
## 356 dodge aspen se 20.0
## 363 amc concord d/l 18.1
## 366 plymouth duster 22.0
## 370 buick century limited 25.0
## 373 chevrolet nova custom 16.0
## 374 plymouth valiant custom 19.0
## 377 volvo 264gl 17.0
## 380 amc hornet 22.5
## 381 mercury zephyr 20.8
## 382 amc matador 16.0
## 385 plymouth fury 18.0
## 394 mercedes-benz 280s 16.5
## 398 amc gremlin 19.0Subset ( subconjunto )
subset(dat,mpg>30 & cylinders==6,select = c("car_name","mpg"))#si no se incluye el select lo que se hace es : Te entrega las columnas que selecciones , de lo contrario , te entrega todas las columnas ## car_name mpg
## 12 volvo diesel 30.7
## 300 oldsmobile cutlass ciera (diesel) 38.0
## 364 datsun 280-zx 32.7Nota importante : Una recoemndacion es , cuando se trabaja con solo unas columnas es mejor trabajar con nombre de la columna
Ejemplo c1:2,2:3 # Nomenclatura & = y ( este y este ) | = o ( este o este ) ! = Lo contrario “==” = igual
#Excuir columnas
dat[1:5,c(-1,-9)]## mpg cylinders displacement horsepower weight acceleration model_year
## 1 28 4 140 90 2264 15.5 71
## 2 19 3 70 97 2330 13.5 72
## 3 36 4 107 75 2205 14.5 82
## 4 28 4 97 92 2288 17.0 72
## 5 21 6 199 90 2648 15.0 70#si se quiere usar con los nombres de las columnas no se puede
#dat[1:5,-c("car_name",mpg)]
#se tiene que hacer de esta forma
excluir.col<-dat[dat$mpg ,!names(dat)%in% c("No","car_name")]# evito el dato escrito en las columnas
excluir.col[1:3,]## mpg cylinders displacement horsepower weight acceleration model_year
## 28 15 8 383 170 3563 10.0 70
## 19 14 8 302 140 4638 16.0 74
## 36 14 8 351 153 4154 13.5 71Como estrarer valores exactos
dat[dat$mpg %in% c(15,20),c("car_name","mpg")]## car_name mpg
## 23 chevrolet vega 20
## 24 audi 100ls 20
## 28 dodge challenger se 15
## 55 chevrolet monza 2+2 20
## 72 ford maverick 15
## 91 mercury cougar brougham 15
## 98 ford galaxie 500 15
## 108 buick skylark 320 15
## 134 chevrolet monte carlo s 15
## 141 amc ambassador dpl 15
## 170 amc matador (sw) 15
## 172 dodge coronet custom 15
## 183 chevrolet nova 15
## 187 plymouth fury iii 15
## 194 toyota carina 20
## 223 amc gremlin 20
## 226 dodge dart custom 15
## 235 chevrolet bel air 15
## 247 plymouth duster 20
## 250 volvo 245 20
## 284 amc matador 15
## 285 chevrolet monte carlo 15
## 292 mercury monarch 15
## 294 toyota mark ii 20
## 356 dodge aspen se 20Diviciones con split /unsplit
Hace subconjuntos y crea listas de DF
car.list<-split(dat,dat$cylinders)como accedemos a este Df
car.list[[1]]## No mpg cylinders displacement horsepower weight acceleration model_year
## 2 2 19.0 3 70 97 2330 13.5 72
## 199 199 18.0 3 70 90 2124 13.5 73
## 251 251 23.7 3 70 100 2420 12.5 80
## 365 365 21.5 3 80 110 2720 13.5 77
## car_name
## 2 mazda rx2 coupe
## 199 maxda rx3
## 251 mazda rx-7 gs
## 365 mazda rx-4compara como se diferencia de usar un corchete doble de uno sencillo
str(dat[1])## 'data.frame': 398 obs. of 1 variable:
## $ No: int 1 2 3 4 5 6 7 8 9 10 ...str(dat[[1]])## int [1:398] 1 2 3 4 5 6 7 8 9 10 ...Particion de data frame con variable numerica
#install.packages("caret")
library(caret)##
## Attaching package: 'caret'## The following object is masked from 'package:survival':
##
## cluster## The following object is masked from 'package:purrr':
##
## liftdat.casa<- read.csv("../../r-course-master/data/tema2/BostonHousing.csv")metV: Valor mediano de casas ocupadas con mil dolares # Training practicion con el 80% de los datos
training.ids<- createDataPartition(dat.casa$MEDV,p=.80, list=F)
data.training<-dat.casa[training.ids,]#
data.validation <- dat.casa[-training.ids,]Se interpreta:crea una particion( array) con el 80% de los datos y mo ,me lo des como lista
Crearemos una particion con el 70%
training.ids2<- createDataPartition(dat.casa$MEDV,p=.70, list=F)
data.training2<-dat.casa[training.ids2,]
temp<- dat.casa[-training.ids2,]
data.validation2 <-createDataPartition(temp$MEDV,p=.5, list=F)
data.testing<- temp[-data.validation2,]Particion de data frame con variable categorica
data2<-read.csv("../../r-course-master/data/tema2/boston-housing-classification.csv", head =T)Paricion 70-30
training.ids.3<- createDataPartition(data2$MEDV_CAT,p=0.7,list = F)
data.training.3<-data2[training.ids.3,]
data.validation3<-data2[-training.ids.3,]#Funciones de 2 con pariciones definidas
rda.cb.partition2<- function(dataframe,target.index,prob){
library(caret)
training.ids<- createDataPartition(dataframe[,target.index],p=prob, list = F)
list(train=dataframe[training.ids,],val=dataframe[-training.ids,])
}Funcion de 3 particiones
rda.cb.partition3<- function(dataframe,target.index,
prob.training,prob.val){
library(caret)
training.ids<- createDataPartition(dataframe[,target.index],p=prob.training, list = F)
train.data<-dataframe[training.ids,]
temp<- dataframe[-training.ids,]
validation.ids<-createDataPartition(temp[,target.index],p=prob.val,list=F)
list(train=train.data,val=validation.ids,test=temp[-validation.ids,])
}veamos como hace las particiones de la funcion de 2 particiones
data.particion2<-rda.cb.partition2(dat.casa,14,0.8)Para partir con la funcion de 3
data.particion3<- rda.cb.partition3(dat.casa,14,0.7,0.5)Como se accede a las propiedades de una lista
data.particion2.1<-data.particion2[[1]]
data.particion2.1[1:3,]## CRIM ZN INDUS CHAS NOX RM AGE DIS RAD TAX PTRATIO B LSTAT
## 2 0.02731 0 7.07 0 0.469 6.421 78.9 4.9671 2 242 17.8 396.90 9.14
## 3 0.02729 0 7.07 0 0.469 7.185 61.1 4.9671 2 242 17.8 392.83 4.03
## 4 0.03237 0 2.18 0 0.458 6.998 45.8 6.0622 3 222 18.7 394.63 2.94
## MEDV
## 2 21.6
## 3 34.7
## 4 33.4data.particion3.1<-data.particion3[[1]]
data.particion3.1[1:3,]## CRIM ZN INDUS CHAS NOX RM AGE DIS RAD TAX PTRATIO B LSTAT
## 1 0.00632 18 2.31 0 0.538 6.575 65.2 4.0900 1 296 15.3 396.90 4.98
## 4 0.03237 0 2.18 0 0.458 6.998 45.8 6.0622 3 222 18.7 394.63 2.94
## 5 0.06905 0 2.18 0 0.458 7.147 54.2 6.0622 3 222 18.7 396.90 5.33
## MEDV
## 1 24.0
## 4 33.4
## 5 36.2obtener una muestra del DF Original ( muestreo aleatorio )
muestra<-sample(dat.casa$CRIM,40,replace = F)# replase = False significa que esos 40 que tome valores quer no estan repoetidos ( Nota :no puedes pedir mas de los que tiene tu DF)
# puede validar el numero de filas con nrowHistogramas ,boxplots y scatterplots
auto.data<- read.csv("../../r-course-master/data/tema2/auto-mpg.csv", stringsAsFactors = F)
#sobre escribiremos a cartor la columna de cilindros
auto.data$cylinders<-factor(auto.data$cylinders,
levels= c(3,4,5,6,8),labels=c("3cil","4cil","5cil","6cil","8cil"))
attach(auto.data)#carga el paquete para trabajar y con esto te evitas referenciar el DF , solo indicas la columna ## The following object is masked from package:ggplot2:
##
## mpgHistograma
hist(acceleration,
col =rainbow(16),#Color del Histograma o rainbow para arcohiris o solo col="Blue"
xlab="Aceleracion",# Nombre del eje X
ylab = "Frecuencia",# Nombre del eje y
main ="Histograma de la aceleraciones",# Titulo que quieres poner
breaks = 16)#numero de diviciones que deseas
Boxplot ( Caja de gatos con bigotes )
Se puede hacer que te mida todo el Df de un jalon
#boxplot(auto.data)boxplot(mpg,
xlab="Millas por Galeon")
Hacer un boxplot en funcion de otra variable
boxplot(mpg~model_year, #muestra las millas por galeon por ano
xlab="millas por galeon" )
Hacer otro boxplot en fincion de otra varioable distinta a la mencionada arriba
boxplot(mpg~cylinders,
xlab="Millas por Galeon en funcion de los cilindros" )
# Scatter Plot (solo variables numericas)
plot(mpg~horsepower)
linearmodel1<-lm(mpg~horsepower)#Crea una linea de regresion lineal
abline(linearmodel1)# sobrepone la linea sobre la grafica anterior 
Es mas funcional tener una matriz de scaterplot
pairs(~mpg+displacement+horsepower+weight)
Como pintar una grafica
plot(mpg~horsepower,#Grafica de nube de puntos
type="n") # No pintes nada ( no pongas los puntos )
linearmodel<-lm(mpg~horsepower)#Crea una linea de regresion lineal
abline(linearmodel)# sobrepone la linea sobre la grafica anterior
# Agregar colores para una variable ( en este caso los cilindros )
with(subset(auto.data,# toma es te subconjunto
cylinders=="8cil"),# que se llama asi
points(horsepower,mpg,# de estas variables que mandate llamar en el grafico
col="red"))# y los coloreas de este color
with(subset(auto.data,cylinders=="6cil"), points(horsepower,mpg,col="blue"))
with(subset(auto.data,cylinders=="5cil"), points(horsepower,mpg,col="green"))
with(subset(auto.data,cylinders=="4cil"), points(horsepower,mpg))# si no se manda llamar algun color , lo pone negro por defoult
with(subset(auto.data,cylinders=="3cil"),points(horsepower,mpg,,col="yellow"))
Histograma con linea de tendencia
Nota el Histograma tiene que estar en frecuencias acumuladas
hist(mpg,prob=T)
lines(density(mpg))
Combinar 2 graficas en una sola imagen con par()
old.par<- par()#guardas que esta era la configuracion inicial de los parametros :ojo esto es importante
par(mfrow=c(1,2))# Cmabias la configuracion ( por columnas, una al lado de otra )
# para mfcol ( los tienes un gafico uno sobre otro )
with(auto.data,{
plot(mpg~weight, main="Peso vs Consumo")
plot(mpg~acceleration,main="Aceleracion Vs el Consumo ")
}
)
par(old.par)# aqui regresas a decirle que solo es para este grafico ( esto que te mandas en rojo, solo es ten cuidado , pero no esta mal )## Warning in par(old.par): el parámetro del gráfico "cin" no puede ser
## especificado## Warning in par(old.par): el parámetro del gráfico "cra" no puede ser
## especificado## Warning in par(old.par): el parámetro del gráfico "csi" no puede ser
## especificado## Warning in par(old.par): el parámetro del gráfico "cxy" no puede ser
## especificado## Warning in par(old.par): el parámetro del gráfico "din" no puede ser
## especificado## Warning in par(old.par): el parámetro del gráfico "page" no puede ser
## especificadoLattise
bwplot xyplot densityplot splom
library(lattice)
bwplot(~auto.data$mpg|auto.data$cylinders,#asi se representan valores(x) con factores(y) de manera rapida
main="MPG segun los cilindros ",
xlab = "Millas por Galeon")
# Cajas de gatos con vigotes de menara rapida con valores numeros y vategorias
xyplot(mpg~weight|cylinders,data = auto.data,
main="Peso Vs Consumo Vs Cilindros",
xlab = "Peso en KG",
ylab="Consumo (MPG)")
#Se puede agregar el valor aspect=1 a la funcion para que guarde la proporcion Comparar a travez de representaciones
alquiler<- read.csv("../../r-course-master/data/tema2/daily-bike-rentals.csv", header = T)
alquiler$season<-factor(alquiler$season, levels = c(1,2,3,4),labels=c("Invierno","Primavera","Verano","Otono"))
alquiler$workingday<-factor(alquiler$workingday, levels = c(0,1),labels=c("festivo ","De trabajo"))
alquiler$weathersit<-factor(alquiler$weathersit, levels = c(1,2,3),labels=c("Despejado","Nublado", "LLuvia/Nieve "))
alquiler$dteday<- as.Date(alquiler$dteday, format = "%Y-%m-%d")# cambiar la columna a dia en el formatro especializado
attach(alquiler)## The following object is masked _by_ .GlobalEnv:
##
## temp## The following object is masked from package:mice:
##
## windspeed4 graficas en una sola imagen
winter <- subset(alquiler,season=="Invierno")$cnt
spring<-subset(alquiler,season=="Primavera")$cnt
summer<-subset(alquiler,season=="Verano")$cnt
fall<-subset(alquiler,season=="Otono")$cntpar(mfrow=c(2,2))#2 filas y 2 columnas hist(winter,prob=T, xlab = “alquiler de bicicletas en Invierno”, main = ““)
par(mfrow=c(2,2))#2 filas y 2 columnas
hist(winter,prob=T, xlab = "alquiler de bicicletas en Invierno", main = "")
lines(density(winter))
abline(v=mean(winter), col="blue")
abline(v=median(winter), col="red")
hist(spring,prob=T, xlab = "alquiler de bicicletas en Primavera", main = "")
lines(density(spring))
abline(v=mean(spring), col="blue")
abline(v=median(spring), col="red")
hist(summer,prob=T, xlab = "alquiler de bicicletas en Verano", main = "")
lines(density(summer))
abline(v=mean(summer), col="blue")
abline(v=median(summer), col="red")
hist(fall,prob=T, xlab = "alquiler de bicicletas en otono", main = "")
lines(density(fall))
abline(v=mean(fall), col="blue")
abline(v=median(fall), col="red")
# Visualizacion de datos
#install.packages("beanplot")
library(beanplot)
beanplot(alquiler$cnt~alquiler$season, col = c("blue","red","green"))
#ctm es la cantidad de bicicletas que se alqularon Posible Causalidad ( Esto se presenta a direccion )
Esto puede ayudar validar hipotesis Analisis de causalidad
bwplot(cnt~weathersit, data = alquiler,
layout=c(1,1),
xlab = "Pronostico del tiempo",
ylab = "Frecuencias"
)
Ahora si queriamos comparar cuantos datos se tienen por distribucion
bwplot(cnt~weathersit, data = alquiler,
layout=c(1,1),
xlab = "Pronostico del tiempo",
ylab = "Frecuencias",
panel = function(x,y,...){
panel.bwplot(x,y,...)
panel.stripplot(x,y,jitter.data = T,...,col="blue")#anade /sobrepone los datos de los puntos
}
)
# Tecnica de Validacion cruzada Se explicara mas a detalle mas adelante # Q-Q Plot Graficos Cuartil- Cuartil ( validar si es normal una muestra o datos )
s<-seq(0.01,0.99,0.01)# crea una secuencia desde el 0.001 hasta el 0.99 en contadores de 0.001
qn<-qnorm(s)#cuartiles normales
df<-data.frame(p=s,q=qn)
sample<-rnorm(200)# se hizo una distribucion normal con 200 datos
quantile(sample,probs = s)## 1% 2% 3% 4% 5% 6%
## -1.911985523 -1.907636194 -1.722155533 -1.660761013 -1.455566078 -1.421415812
## 7% 8% 9% 10% 11% 12%
## -1.395348406 -1.354863151 -1.284591665 -1.250437500 -1.224555833 -1.182464018
## 13% 14% 15% 16% 17% 18%
## -1.161734550 -1.098535783 -1.046278868 -1.002362907 -0.945845418 -0.877679720
## 19% 20% 21% 22% 23% 24%
## -0.827154004 -0.799978937 -0.765828909 -0.733186183 -0.704539796 -0.672544002
## 25% 26% 27% 28% 29% 30%
## -0.618693565 -0.595425049 -0.580098186 -0.557597207 -0.518925713 -0.480378379
## 31% 32% 33% 34% 35% 36%
## -0.442520622 -0.403188988 -0.399510601 -0.386718482 -0.336680094 -0.306370977
## 37% 38% 39% 40% 41% 42%
## -0.251621677 -0.221982820 -0.195486862 -0.139503108 -0.087338992 -0.078656962
## 43% 44% 45% 46% 47% 48%
## -0.070741933 -0.049643005 0.001783342 0.024594170 0.036941386 0.059434399
## 49% 50% 51% 52% 53% 54%
## 0.087524221 0.103231341 0.110023048 0.125529223 0.147124309 0.166610699
## 55% 56% 57% 58% 59% 60%
## 0.183045674 0.216728447 0.235484293 0.246049550 0.260361740 0.295438080
## 61% 62% 63% 64% 65% 66%
## 0.312726100 0.323398672 0.351307497 0.397524013 0.421949682 0.447452355
## 67% 68% 69% 70% 71% 72%
## 0.465805467 0.488676388 0.557884141 0.603094547 0.616598823 0.630499579
## 73% 74% 75% 76% 77% 78%
## 0.701442327 0.751856787 0.829118991 0.891602646 0.920840932 0.952813717
## 79% 80% 81% 82% 83% 84%
## 0.967495272 0.985640922 0.997784678 1.011041585 1.125196823 1.175981055
## 85% 86% 87% 88% 89% 90%
## 1.250664737 1.278877292 1.310043392 1.346990786 1.378243725 1.390780416
## 91% 92% 93% 94% 95% 96%
## 1.408825439 1.471364986 1.542513636 1.550994193 1.582807037 1.684863625
## 97% 98% 99%
## 1.808690558 1.896880097 2.328570115qqnorm (tiene que seguir aproximadamente una linea recta )
qqnorm(trees$Height)
# como diferenciar cuando es una distribucion normal de una uniforme Se ve como una ojiva y recordemos que tiene que ser una linea recta
qqnorm(randu$x)
qqplot otra forma de hacer grafica quantil-cuantil para distribuciones uniformes
#randu es una distribucion uniforme( todos tienen la misma probabilidad de salir )
n<-length(randu$x)# obtiene el largo de los datos ( son 400)
y<- qunif(ppoints(n))#qunif distribucion uniforme
#ppoints genera un nuemro de probabilidades , en este caso , ya sacamos el largo por eso es n
qqplot(y,randu$x)
representacion de ji-cuadrado
Asi se ve una distribucion de ji cuadrado
chi3<- qchisq(ppoints(30)#qchiq= cuantil chi square
,df=3)# df son los grados de libertad
n30<- qnorm(ppoints(30))
qqplot(n30,chi3)
distribucion de cuauchy
se ve como una onda y tiene que ser una linea recta
cauchy<-qcauchy(ppoints(30))
qqplot(n30,cauchy)
Ahora , si comparas la misma
x<- seq(-3,3,0.01)# aqui solo generas numeros
par(mfrow=c(3,2))# configuracion par que sean 2 graficas
plot(x,dnorm(x),
main="Distibucion Normal")#funcion de distribucion
plot(x,pnorm(x),
main="Probabilidad de una normal")# funcion de probabilidad acumulativa
plot(x,dchisq(x,df=3),
#Grados de libertad
main="Distibucion chi")#funcion de distribucion de chi
plot(x,pchisq(x,df=3),
main="probabilidad chi")#funcion de probabilidad de Chi
plot(x,dcauchy(x),
#Grados de libertad
main="Distibucion chi")#funcion de distribucion de cauchi
plot(x,pcauchy(x),
main="probabilidad chi")#funcion de probabilidad de cauchi
Clasificacion en data science
Matices de confucion
calificaciones<-read.csv("../../r-course-master/data/tema3/college-perf.csv",header = T)
calificaciones$Perf<-ordered(calificaciones$Perf,levels=c("Low","Medium","High"))# odenamos los valores ne este orden , en el DF no se ve pero ya lo hizo
calificaciones$Pred<-ordered(calificaciones$Pred,levels=c("Low","Medium","High"))# mismo comentario de arriba pero a otr columna
table<-table(calificaciones$Perf,calificaciones$Pred,dnn = c("Actual"," Prediccion"))
table# matriz de confucion ## Prediccion
## Actual Low Medium High
## Low 1150 84 98
## Medium 166 1801 170
## High 35 38 458si se quiere saber en porcentaje
prop.table(table)*100#en porcentaje ( pero esto es del total) ## Prediccion
## Actual Low Medium High
## Low 28.750 2.100 2.450
## Medium 4.150 45.025 4.250
## High 0.875 0.950 11.450round(prop.table(table,1)*100,2)#esto es por filas## Prediccion
## Actual Low Medium High
## Low 86.34 6.31 7.36
## Medium 7.77 84.28 7.96
## High 6.59 7.16 86.25round(prop.table(table,2)*100,2)#por columnas ## Prediccion
## Actual Low Medium High
## Low 85.12 4.37 13.50
## Medium 12.29 93.66 23.42
## High 2.59 1.98 63.09Matiz de confusión en una grafica
barplot(table,leyend=T,
Xlab="Nota predecida por el modelo")## Warning in plot.window(xlim, ylim, log = log, ...): "leyend" is not a graphical
## parameter## Warning in plot.window(xlim, ylim, log = log, ...): "Xlab" is not a graphical
## parameter## Warning in axis(if (horiz) 2 else 1, at = at.l, labels = names.arg, lty =
## axis.lty, : "leyend" is not a graphical parameter## Warning in axis(if (horiz) 2 else 1, at = at.l, labels = names.arg, lty =
## axis.lty, : "Xlab" is not a graphical parameter## Warning in title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...):
## "leyend" is not a graphical parameter## Warning in title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...): "Xlab"
## is not a graphical parameter## Warning in axis(if (horiz) 1 else 2, cex.axis = cex.axis, ...): "leyend" is not
## a graphical parameter## Warning in axis(if (horiz) 1 else 2, cex.axis = cex.axis, ...): "Xlab" is not a
## graphical parameter
mosaicplot(table,main="Eficiencia del modelo")
summary(table)# el p-valor indica dependencia entre el valor actual y el predicho , es decir el modelo es eficiente## Number of cases in table: 4000
## Number of factors: 2
## Test for independence of all factors:
## Chisq = 4449, df = 4, p-value = 0Analisis de Componentes Principales
Tranforma los datos de muchas dimenciones y busca los basicos
crimenes<-read.csv("../../r-course-master/data/tema3/USArrests.csv")
rownames(crimenes)<-crimenes$X #cambia como el identificador a una columna en particular
crimenes$X<- NULL# Eliminamos la fila que ya no necesitamos
apply(crimenes,
2,#por columnas
var)# la funcion varianza## Murder Assault UrbanPop Rape
## 18.97047 6945.16571 209.51878 87.72916Podemos observar que la que tiene mayor varianza es la de assault
Normaliza y estandarizar los datos
acp<-prcomp(crimenes,# analisis de componentes principales
center = T,#Resta la media de cada variable
scale. = T)#y divide por la desviación estandar
acp## Standard deviations (1, .., p=4):
## [1] 1.5748783 0.9948694 0.5971291 0.4164494
##
## Rotation (n x k) = (4 x 4):
## PC1 PC2 PC3 PC4
## Murder -0.5358995 0.4181809 -0.3412327 0.64922780
## Assault -0.5831836 0.1879856 -0.2681484 -0.74340748
## UrbanPop -0.2781909 -0.8728062 -0.3780158 0.13387773
## Rape -0.5434321 -0.1673186 0.8177779 0.08902432Graficar la lista Nota:usar la tecnica del codo la cual consiste en donde la tecnica del “codo del brazo”, por ejemplo, en este caso lo que se mostraria, serian las 2 variables leyendolas de izquierda a derecha
plot(acp,type="l")
summary(acp)## Importance of components:
## PC1 PC2 PC3 PC4
## Standard deviation 1.5749 0.9949 0.59713 0.41645
## Proportion of Variance 0.6201 0.2474 0.08914 0.04336
## Cumulative Proportion 0.6201 0.8675 0.95664 1.00000cuando se observa el summary de acp en la caracteristica cumulative proportion significa la cantidad de datos que se representan #representacion de acp
biplot(acp,scale=0)
si se observa hacia donde se mueve mas es rape en el espacio PC2 #Analizando la diferentes componentes principales
pc1<- apply(acp$rotation[,1]*crimenes,1,sum)
pc2<- apply(acp$rotation[,1]*crimenes,1,sum)#de la lista sumamos el vector de rotacion , recuerda que son matrices
crimenes$pc1<-pc1#sumamos una nueva fila a dataset
crimenes$pc2<-pc2curva ROC ( caracteristicas operativas del receptor )
Es una tecnica de clasificacion: Sensibilidad vs Especificidad Entre mas sensibilidad es mejor y entre menos especificidad es mejor Cargamos los paquetes
#install.packages("ROCR")
library(ROCR)
curvaroc1<- read.csv("../../r-course-master/data/tema3/roc-example-1.csv")# variables numericas
curvaroc2<- read.csv("../../r-course-master/data/tema3/roc-example-2.csv")# variables categoricasPara el DF de variables numericas tpr=true, positive rate, fpr=False Positive Rate
#nota : las clases del DF significa: 0 fallo y 1 exito
pred1<-prediction(curvaroc1$prob,curvaroc1$class)#Con respecto a este objeto
perf1<-performance(pred1,"tpr","fpr")
#Objeto de probabilidad con respecto ala predicción
plot(perf1)
lines(par()$usr[1:2],par()$usr[3:4])
Cuales son las variables que son los que se tienen que perfecionar
probcortes<-data.frame(cut=perf1@alpha.values[[1]],#valor de corte
pfr=perf1@x.values[[1]],#toma los valores de perfilamiento del eje x
tpr=perf1@y.values)#toma los valores de perfilamiento del eje x
colnames(probcortes)<-c("cut","pfr","tpr")# renombraremos las columnas
head(probcortes)#La cabeza de los datos ## cut pfr tpr
## 1 Inf 0 0.00000000
## 2 0.9917340 0 0.01851852
## 3 0.9768288 0 0.03703704
## 4 0.9763148 0 0.05555556
## 5 0.9601505 0 0.07407407
## 6 0.9351574 0 0.09259259tail(probcortes)#Cola de los datos ## cut pfr tpr
## 96 0.10426897 0.8913043 1
## 97 0.07292866 0.9130435 1
## 98 0.07154785 0.9347826 1
## 99 0.04703280 0.9565217 1
## 100 0.04652589 0.9782609 1
## 101 0.00112760 1.0000000 1Se debe encontar un equilibrio donde la empresa se sienta comoda con los verdederos positivos
probcortes[probcortes$tpr>.80,]#dame todos los valores verdaderos positivos por arriba del 80%## cut pfr tpr
## 55 0.49815058 0.2173913 0.8148148
## 56 0.49616956 0.2173913 0.8333333
## 57 0.47840739 0.2391304 0.8333333
## 58 0.47754679 0.2608696 0.8333333
## 59 0.46323419 0.2826087 0.8333333
## 60 0.45227354 0.2826087 0.8518519
## 61 0.44950615 0.3043478 0.8518519
## 62 0.43443516 0.3260870 0.8518519
## 63 0.43274841 0.3478261 0.8518519
## 64 0.42845777 0.3695652 0.8518519
## 65 0.40701592 0.3695652 0.8703704
## 66 0.40272660 0.3695652 0.8888889
## 67 0.40248242 0.3913043 0.8888889
## 68 0.40140505 0.3913043 0.9074074
## 69 0.37646732 0.4130435 0.9074074
## 70 0.36254324 0.4347826 0.9074074
## 71 0.35547851 0.4565217 0.9074074
## 72 0.34872470 0.4782609 0.9074074
## 73 0.33814262 0.5000000 0.9074074
## 74 0.33528819 0.5217391 0.9074074
## 75 0.33041954 0.5434783 0.9074074
## 76 0.31878806 0.5652174 0.9074074
## 77 0.31438300 0.5869565 0.9074074
## 78 0.31164060 0.6086957 0.9074074
## 79 0.30761685 0.6086957 0.9259259
## 80 0.23584000 0.6304348 0.9259259
## 81 0.23310990 0.6304348 0.9444444
## 82 0.22243908 0.6521739 0.9444444
## 83 0.20977357 0.6739130 0.9444444
## 84 0.19511299 0.6739130 0.9629630
## 85 0.18730064 0.6739130 0.9814815
## 86 0.18210618 0.6956522 0.9814815
## 87 0.17060700 0.7173913 0.9814815
## 88 0.15536249 0.7391304 0.9814815
## 89 0.14907332 0.7391304 1.0000000
## 90 0.14288827 0.7608696 1.0000000
## 91 0.13868194 0.7826087 1.0000000
## 92 0.13090790 0.8043478 1.0000000
## 93 0.12741177 0.8260870 1.0000000
## 94 0.12676979 0.8478261 1.0000000
## 95 0.11107864 0.8695652 1.0000000
## 96 0.10426897 0.8913043 1.0000000
## 97 0.07292866 0.9130435 1.0000000
## 98 0.07154785 0.9347826 1.0000000
## 99 0.04703280 0.9565217 1.0000000
## 100 0.04652589 0.9782609 1.0000000
## 101 0.00112760 1.0000000 1.0000000Para el DF de variables numericas
predict2<-prediction(curvaroc2$prob,curvaroc2$class,label.ordering = c("non-buyer","buyer"))
# nota importante :deben cohincidri exactamente los nombres con las etiquetas del DF
#Probabilidad mas baja de compra a la mas alta ( en ese orden )
perf2<-performance(predict2,"tpr","fpr")
plot(perf2)
lines(par()$usr[1:2],par()$usr[3:4])# pinta la linea recta # Arboles de clasificacion
#install.packages(c("rpart","rpart.plot","caret"))
library(rpart)#construye los arboles
library(rpart.plot)#grafico de esos arboles## Warning: package 'rpart.plot' was built under R version 4.1.3library(caret)
library(ggplot2)
library(lattice)Cargamos los datos
banknote<-read.csv("../../r-course-master/data/tema3/banknote-authentication.csv")#caret:: set.seed(2018)# semilla
training.ids<-caret::createDataPartition(banknote$class,p=0.7,list=F)#revisar Particion de data frame con variable numerica
mod<-rpart(class~.,# especifica las variables dependientes e independientes ( clase depende de todas las anteriores(.))
data = banknote[training.ids,],#De donde va a tomar los datos ( de solo los datos de entrenamiento)
method = "class",#Clasifique casos
control = rpart.control(minsplit = 20,#Solo debe tener mas elementos en cada clase
cp=0.01))#ajustar la complejidad
mod # Asi se ve el modelo de forma (analitica )## n= 961
##
## node), split, n, loss, yval, (yprob)
## * denotes terminal node
##
## 1) root 961 411 0 (0.572320499 0.427679501)
## 2) variance>=0.320165 511 52 0 (0.898238748 0.101761252)
## 4) curtosis>=-4.3856 487 34 0 (0.930184805 0.069815195)
## 8) variance>=1.5922 349 3 0 (0.991404011 0.008595989) *
## 9) variance< 1.5922 138 31 0 (0.775362319 0.224637681)
## 18) curtosis>=-2.2525 116 11 0 (0.905172414 0.094827586) *
## 19) curtosis< -2.2525 22 2 1 (0.090909091 0.909090909) *
## 5) curtosis< -4.3856 24 6 1 (0.250000000 0.750000000)
## 10) variance>=2.3098 7 1 0 (0.857142857 0.142857143) *
## 11) variance< 2.3098 17 0 1 (0.000000000 1.000000000) *
## 3) variance< 0.320165 450 91 1 (0.202222222 0.797777778)
## 6) skew>=7.5653 78 14 0 (0.820512821 0.179487179)
## 12) variance>=-4.8358 64 0 0 (1.000000000 0.000000000) *
## 13) variance< -4.8358 14 0 1 (0.000000000 1.000000000) *
## 7) skew< 7.5653 372 27 1 (0.072580645 0.927419355)
## 14) skew>=5.86535 20 8 1 (0.400000000 0.600000000)
## 28) variance>=-2.38228 8 0 0 (1.000000000 0.000000000) *
## 29) variance< -2.38228 12 0 1 (0.000000000 1.000000000) *
## 15) skew< 5.86535 352 19 1 (0.053977273 0.946022727)
## 30) curtosis>=6.21865 101 18 1 (0.178217822 0.821782178)
## 60) skew>=-4.72105 18 0 0 (1.000000000 0.000000000) *
## 61) skew< -4.72105 83 0 1 (0.000000000 1.000000000) *
## 31) curtosis< 6.21865 251 1 1 (0.003984064 0.996015936) *Como verlo de manera grafica
prp(mod,type=2,#cada nodo quede etiquetado y sea exactamente a la mitad( investigar/ variar)
extra=104,#mostrar la probabilidad de cada uno de los casos por en base acada rama
nn=T,# numero del nodo
fallen.leaves=T,#Nodos finales los muestra en la parte final
faclen=4,# abreviar los nombres de las clases
varlen=8,# abreviar la longitud de los nombres de las variables
shadow.col="blue")#sombras 
# Podar arboles de Clasificacion Seleccionar solo algunas ramas muestra las componentes principales
mod$cptable## CP nsplit rel error xerror xstd
## 1 0.65206813 0 1.00000000 1.00000000 0.03731631
## 2 0.12165450 1 0.34793187 0.35036496 0.02692085
## 3 0.03406326 2 0.22627737 0.24817518 0.02323235
## 4 0.02919708 3 0.19221411 0.20681265 0.02141697
## 5 0.02189781 4 0.16301703 0.18491484 0.02035518
## 6 0.01581509 6 0.11922141 0.12408759 0.01690839
## 7 0.01216545 10 0.05596107 0.10705596 0.01576549
## 8 0.01000000 11 0.04379562 0.09975669 0.01524342se debe seleccionar minimo xerror y sumada alguna la fila de xstd( desviacion estandar) debe dar un valor minimo de xerror para este ejemplo tomamos la numero 6( desviacion estandar wey) # Creamos el modelo de clasificacion recortado
mod.pruned<-prune(mod,mod$cptable[6,"CP"])Ver el modelo de clasificacion de arboles “Podado”
prp(mod.pruned,type=2,#cada nodo quede etiquetado y sea exactamente a la mitad( investigar/ variar)
extra=104,#mostrar la probabilidad de cada uno de los casos por en base acada rama
nn=T,# numero del nodo
fallen.leaves=T,#Nodos finales los muestra en la parte final
faclen=4,# abreviar los nombres de las clases
varlen=8,# abreviar la longitud de los nombres de las variables
shadow.col="blue")#sombras 
Esto es una Predicion hecha con el 70% de los datos
Con los datos totales coteja que tan bien se apega a los valores de el DF original
pred.pruned<-predict(mod,banknote[-training.ids,],type = "class")# generea la prediccion
table(banknote[-training.ids,]$class,pred.pruned,dnn=c("Actual","Predicho"))## Predicho
## Actual 0 1
## 0 208 4
## 1 8 191Se ve excelente
#comparando con el arbol “podado” Se han movido los datos y subieron los errores
pred.pruned.recort<-predict(mod.pruned,banknote[-training.ids,],type = "class")# generea la prediccion
table(banknote[-training.ids,]$class,pred.pruned.recort,dnn=c("Actual","Predicho"))## Predicho
## Actual 0 1
## 0 194 18
## 1 7 192Diagramas ROC para arboles
pred.pruned.predict<-predict(mod.pruned,banknote[-training.ids,],type = "prob")# genera modelo probabilistico
pred.prob<-ROCR::prediction(pred.pruned.predict[,2],banknote[-training.ids,"class"])
perf.pruned.pred<-ROCR::performance(pred.prob,"tpr","fpr")
plot(perf.pruned.pred)
Esto si es un pinche buen modelo
Random Forest
#install.packages("randomForest")
library(randomForest)## randomForest 4.7-1## Type rfNews() to see new features/changes/bug fixes.##
## Attaching package: 'randomForest'## The following object is masked from 'package:dplyr':
##
## combine## The following object is masked from 'package:ggplot2':
##
## marginlibrary(ROCR)
library(caret)
library(ggplot2)Cargamos los datos
random.forest<-read.csv("../../r-course-master/data/tema3/banknote-authentication.csv")Cambiemos a factores la columna de class
random.forest$class<-factor(random.forest$class)Hagamos una semilla y particion de datos
set.seed(2018)
trainining.random.fores<-createDataPartition(random.forest$class,p=.70,list=F)Apliquemos la tecnica de Random Forest
mod<-randomForest(x=random.forest[trainining.random.fores,1:4],# todas las variables independientes ( todo el Df)
y=random.forest[trainining.random.fores,5],# esta es la variable que queremos predecir
ntree = 500,#numero de arboles
keep.forest =TRUE)#se queda con todos los arboles intermedios ( para predicciones futuras)Predecimos los valores del Df con los valores predichos por el modelo
pred.random.forest<-predict(mod,random.forest[-trainining.random.fores,])
table(random.forest[-trainining.random.fores,"class"],pred.random.forest,dnn=c("Actual","Predicho"))## Predicho
## Actual 0 1
## 0 225 3
## 1 0 183Diagrama ROC para el Random Forest
prob.random.forest<-predict(mod,random.forest[-trainining.random.fores,],type = "prob")# genero probabilidades
head(prob.random.forest)## 0 1
## 4 1.000 0.000
## 6 1.000 0.000
## 7 1.000 0.000
## 9 1.000 0.000
## 10 0.958 0.042
## 12 1.000 0.000library(ROCR)
pred.random.forest.roc<-ROCR::prediction(prob.random.forest[,2],random.forest[-trainining.random.fores,"class"])#Predecir
perf.random.forest<-performance(pred.random.forest.roc,"tpr","fpr")# desempeno de las prediciiones
plot(perf.random.forest)
Este modelo es la leche # maquinas de soporte vectorial ( Suport Verctor Machine ) Hiper plano de separacion se trata de encontrar la linea que separa ambos conjuntos de datos
#install.packages("e1071")
library(e1071)##
## Attaching package: 'e1071'## The following objects are masked from 'package:moments':
##
## kurtosis, moment, skewness## The following object is masked from 'package:raster':
##
## interpolate## The following object is masked from 'package:modeest':
##
## skewness## The following object is masked from 'package:Hmisc':
##
## imputelibrary(caret)
banknote.svm<-read.csv("../../r-course-master/data/tema3/banknote-authentication.csv")
banknote.svm$class =factor(banknote.svm$class)Plantamos la semilla
set.seed(2018)
t.ids<-createDataPartition(banknote.svm$class,p=.70,list = F)
#La fuccion svm toma una media y desviacion estandar ( unificada )
# el parametro puede darle diferentes pesos a las clases en las cuales se quieres orrer la funcion svm ( class.weights =c("0"=0.03),"1"=0.7))
# cost =1 (por default)es el parametro donde se le dice al modelo qe margen de error podemos darle
mod.svm<-svm(class~.,#tipo de modelo ( calsificacion o de regresion )cuando es clasificacion solo separa
data=banknote.svm[t.ids,],
cost=1)Como calculo el margen de error mas eficiente Investigar un poco mas, ya que aqui nos sale un error
tuned<-tune.svm(class~.,data=banknote.svm[t.ids,],gama=10^(-6:-1),10^(1:2))
summary(tuned)##
## Error estimation of 'svm' using 10-fold cross validation: 0Grafiquemos los datos del svm
table(banknote.svm[t.ids,"class"],fitted(mod.svm)#Trae los valores dentro de si , vienen intrinsecos
,dnn=c("Actual","Predicho"))## Predicho
## Actual 0 1
## 0 534 0
## 1 0 427NO existe error alguno
Se fue catalogado de manera correcta
pred.smv<-predict(mod.svm,banknote.svm[-t.ids,])# predice
table(banknote.svm[-t.ids,"class"],pred.smv,dnn=c("Actual","Predicho"))## Predicho
## Actual 0 1
## 0 228 0
## 1 0 183Grafica de 2 variables ( en este caso con los valores de entrenamiento)
Grafiquemoslo para ver como se observa como tenemos 4 variables seleccionamos 2 , pero esto es para los valores de entrenamiento
plot(mod.svm,data=banknote.svm[t.ids,],skew~variance)
# Grafica del SVM sin los valores de entrenamiento
plot(mod.svm,data=banknote.svm[-t.ids,],skew~variance)
Generemos la interpolacion de bayes
pred.bayes<-predict(mod.bayes,ep[-t.id.bayes,])# prediccion
# tabla de confucion
tab.bayes<-table(ep[-t.id.bayes,]$Purchase,pred.bayes,dnn=c("Actual", "Predicho"))
confusionMatrix(tab.bayes)## Confusion Matrix and Statistics
##
## Predicho
## Actual No Yes
## No 1 1
## Yes 2 0
##
## Accuracy : 0.25
## 95% CI : (0.0063, 0.8059)
## No Information Rate : 0.75
## P-Value [Acc > NIR] : 0.9961
##
## Kappa : -0.5
##
## Mcnemar's Test P-Value : 1.0000
##
## Sensitivity : 0.3333
## Specificity : 0.0000
## Pos Pred Value : 0.5000
## Neg Pred Value : 0.0000
## Prevalence : 0.7500
## Detection Rate : 0.2500
## Detection Prevalence : 0.5000
## Balanced Accuracy : 0.1667
##
## 'Positive' Class : No
## K-Nearest Neighbors
Es una metodologia no parametrica de clasificacion y regresion
#install.packages("class")# regularmente instalado pero se puede actualizar
library(class)
library(caret)cargamos el archivo
vacaciones<-read.csv("../../r-course-master/data/tema3/vacation-trip-classification.csv")
anyNA(vacaciones)# existen NA## [1] FALSENormalizaremos el DF con 2 nuevas variables ya que estan super diferentes los valores de entrada
vacaciones$income.z<-scale(vacaciones$Income)#crea una columna normalizada
vacaciones$Family_size.z<-scale(vacaciones$Family_size)#crea una columna normalizadaSe hara de la siguiente forma 50% del Df para entrenar 25% para validar 25% para testear
set.seed(2018)
#50% de los datos
t.ids.kneigh<-createDataPartition(vacaciones$Result,p=0.5, list = F)#Data de entrenamiento
train.kneigh<-vacaciones[t.ids.kneigh,]# Data de entrenamiento
temp.kneigh<-vacaciones[-t.ids.kneigh,]# es el restante es decir todos los que no son el conjunto de entrenamiento
# Aqui se encuentra la partición del 25% que esta ( 50% del 50%)
validation.id.kneigh<-createDataPartition(temp.kneigh$Result, p=0.5, list=F)#Datos de validacion
validation.kneigh<-temp.kneigh[validation.id.kneigh,]
test.kneigh<- temp.kneigh[-t.ids.kneigh,]Empecemos a predecir , para este momento es iterativo , es decir manualmente se varia k=1,2,3,4..n
pred1.kneigh<-knn(train.kneigh[,4:5],validation.kneigh[,4:5],train.kneigh[,3],k=1)
errormatrix1<-table(validation.kneigh$Result,pred1.kneigh,dnn=c("Actual","Predict"))
errormatrix1## Predict
## Actual Buyer Non-buyer
## Buyer 3 2
## Non-buyer 5 0Generareremos una predidcion 2 y sale diferente y observamos que con K=5 nos manda un error de tan solo 3 datos que no predice
pred2.kneigh<-knn(train.kneigh[,4:5],test.kneigh[,4:5],train.kneigh[,3],k=5)
errormatrix2<- table(test.kneigh$Result,pred2.kneigh,dnn =c("Actual","Predict"))
errormatrix2## Predict
## Actual Buyer Non-buyer
## Buyer 4 2
## Non-buyer 1 4Funcion que elija el mejor numero de vecinos en la metodologia K-Nearest Neighbors
knn.automate<-function(tr.prediction, val.prediction,tr.target,val.target,start.k,end.k){
for(k in start.k:end.k){
pred.k<-knn(tr.prediction,val.prediction,tr.target,k)
tab<- table(val.target,pred.k,dnn = c("Actual","Predicho"))
cat(paste("Matriz de confusion para k=",k,"\n"))
cat("===============================\n")
print(tab)
cat("-------------------------------\n")
}
}Como chingaos se utiliza esta funcion recuerda que el nombre de las variables puede variar dependeindo de como tu llamaste tus variables
knn.automate(train.kneigh[,4:5],# el conjunto de entrenamiento en este caso normalizado
validation.kneigh[,4:5],#las columnas de validacion
train.kneigh[,3],# le decimos que de la columna 3 que es por la categoria que queremos clasificar
validation.kneigh[,3],#tambien le decimos que lo queremos por esta variable de este conjunto de datos
1,8#Aqui le decimos de que numero a que numero va a iterar
)## Matriz de confusion para k= 1
## ===============================
## Predicho
## Actual Buyer Non-buyer
## Buyer 3 2
## Non-buyer 4 1
## -------------------------------
## Matriz de confusion para k= 2
## ===============================
## Predicho
## Actual Buyer Non-buyer
## Buyer 3 2
## Non-buyer 4 1
## -------------------------------
## Matriz de confusion para k= 3
## ===============================
## Predicho
## Actual Buyer Non-buyer
## Buyer 3 2
## Non-buyer 5 0
## -------------------------------
## Matriz de confusion para k= 4
## ===============================
## Predicho
## Actual Buyer Non-buyer
## Buyer 3 2
## Non-buyer 4 1
## -------------------------------
## Matriz de confusion para k= 5
## ===============================
## Predicho
## Actual Buyer Non-buyer
## Buyer 3 2
## Non-buyer 1 4
## -------------------------------
## Matriz de confusion para k= 6
## ===============================
## Predicho
## Actual Buyer Non-buyer
## Buyer 3 2
## Non-buyer 1 4
## -------------------------------
## Matriz de confusion para k= 7
## ===============================
## Predicho
## Actual Buyer Non-buyer
## Buyer 3 2
## Non-buyer 1 4
## -------------------------------
## Matriz de confusion para k= 8
## ===============================
## Predicho
## Actual Buyer Non-buyer
## Buyer 3 2
## Non-buyer 0 5
## -------------------------------ya existe un algoritmo que no es a “mano” se puede usar la funcion de carte::trainControl
training.control<-trainControl(method = "repeatedcv",#Valodacion cruzada repetida
number = 10,#numero de remuestreos
repeats = 3)#repeticiones de la validacion cruzada
cater.knn.fit<-train(Result~Family_size.z+Income,data = train.kneigh,
method="knn",#metodo por el cual se va a hacer
trControl=training.control,#el de arriba
preProcess=c("center","scale"),#Centrar y restar la media
tuneLength=10)#longitud del ajuste## Warning in knn3Train(train = structure(c(1.31914046975573, -0.860309002014607, :
## k = 21 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.31914046975573, -0.860309002014607, :
## k = 23 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.28892483532644, -0.937399880237408, :
## k = 21 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.28892483532644, -0.937399880237408, :
## k = 23 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.44871854585396, -0.84508581841481, :
## k = 21 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.44871854585396, -0.84508581841481, :
## k = 23 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.35320178424367, -0.789367707475475, :
## k = 21 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.35320178424367, -0.789367707475475, :
## k = 23 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.79845760605179, -0.83005735663929, :
## k = 21 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.79845760605179, -0.83005735663929, :
## k = 23 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.48638756616766, -0.772921534407183, :
## k = 21 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.48638756616766, -0.772921534407183, :
## k = 23 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(-0.772921534407183,
## -0.772921534407183, : k = 21 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(-0.772921534407183,
## -0.772921534407183, : k = 23 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.31914046975573, -0.860309002014607, :
## k = 21 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.31914046975573, -0.860309002014607, :
## k = 23 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.23421938529559, -0.987375508236473, :
## k = 19 exceeds number 18 of patterns## Warning in knn3Train(train = structure(c(1.23421938529559, -0.987375508236473, :
## k = 21 exceeds number 18 of patterns## Warning in knn3Train(train = structure(c(1.23421938529559, -0.987375508236473, :
## k = 23 exceeds number 18 of patterns## Warning in knn3Train(train = structure(c(1.28892483532644, -0.937399880237408, :
## k = 21 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.28892483532644, -0.937399880237408, :
## k = 23 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.44871854585396, -0.84508581841481, :
## k = 21 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.44871854585396, -0.84508581841481, :
## k = 23 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.31914046975573, -0.860309002014607, :
## k = 21 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.31914046975573, -0.860309002014607, :
## k = 23 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.31914046975573, -0.860309002014607, :
## k = 21 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.31914046975573, -0.860309002014607, :
## k = 23 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.28892483532644, -0.937399880237408, :
## k = 21 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.28892483532644, -0.937399880237408, :
## k = 23 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.28892483532644, -0.937399880237408, :
## k = 21 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.28892483532644, -0.937399880237408, :
## k = 23 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.35320178424367, -0.789367707475475, :
## k = 21 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.35320178424367, -0.789367707475475, :
## k = 23 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.31914046975573, -0.860309002014607, :
## k = 21 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.31914046975573, -0.860309002014607, :
## k = 23 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(-0.84508581841481, -0.84508581841481, :
## k = 21 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(-0.84508581841481, -0.84508581841481, :
## k = 23 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.75609839671568, -0.913171166292152, :
## k = 21 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.75609839671568, -0.913171166292152, :
## k = 23 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.46926177338491, -0.734630886692453, :
## k = 19 exceeds number 18 of patterns## Warning in knn3Train(train = structure(c(1.46926177338491, -0.734630886692453, :
## k = 21 exceeds number 18 of patterns## Warning in knn3Train(train = structure(c(1.46926177338491, -0.734630886692453, :
## k = 23 exceeds number 18 of patterns## Warning in knn3Train(train = structure(c(1.31914046975573, -0.860309002014607, :
## k = 21 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.31914046975573, -0.860309002014607, :
## k = 23 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.26324401467532, -0.902317153339516, :
## k = 19 exceeds number 18 of patterns## Warning in knn3Train(train = structure(c(1.26324401467532, -0.902317153339516, :
## k = 21 exceeds number 18 of patterns## Warning in knn3Train(train = structure(c(1.26324401467532, -0.902317153339516, :
## k = 23 exceeds number 18 of patterns## Warning in knn3Train(train = structure(c(1.28892483532644, -0.937399880237408, :
## k = 21 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.28892483532644, -0.937399880237408, :
## k = 23 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(2.05211309734947, -0.836046076697931, :
## k = 21 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(2.05211309734947, -0.836046076697931, :
## k = 23 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.31914046975573, -0.860309002014607, :
## k = 21 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.31914046975573, -0.860309002014607, :
## k = 23 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.31914046975573, -0.860309002014607, :
## k = 21 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.31914046975573, -0.860309002014607, :
## k = 23 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.28892483532644, -0.937399880237408, :
## k = 21 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.28892483532644, -0.937399880237408, :
## k = 23 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(-0.772921534407183,
## -0.772921534407183, : k = 21 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(-0.772921534407183,
## -0.772921534407183, : k = 23 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.35320178424367, -0.789367707475475, :
## k = 21 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.35320178424367, -0.789367707475475, :
## k = 23 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.44871854585396, -0.84508581841481, :
## k = 21 exceeds number 19 of patterns## Warning in knn3Train(train = structure(c(1.44871854585396, -0.84508581841481, :
## k = 23 exceeds number 19 of patternscater.knn.fit## k-Nearest Neighbors
##
## 21 samples
## 2 predictor
## 2 classes: 'Buyer', 'Non-buyer'
##
## Pre-processing: centered (2), scaled (2)
## Resampling: Cross-Validated (10 fold, repeated 3 times)
## Summary of sample sizes: 19, 19, 19, 19, 19, 19, ...
## Resampling results across tuning parameters:
##
## k Accuracy Kappa
## 5 0.7944444 0.58000000
## 7 0.7722222 0.53000000
## 9 0.6833333 0.35000000
## 11 0.7333333 0.45000000
## 13 0.6500000 0.28333333
## 15 0.7444444 0.48333333
## 17 0.7111111 0.43000000
## 19 0.5277778 0.04666667
## 21 0.5055556 0.01000000
## 23 0.4722222 -0.04333333
##
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was k = 5.Como lo menciona la salida , este ejercicio el algoritmo encuentra que el mas optimo es de K=5
ahora que ya conocemos el valor mas optimo de K , correremos el valor mas eficiente
pred5<- knn(train.kneigh[,4:5],validation.kneigh[,4:5],train.kneigh[,3],k=5,prob = T)
pred5## [1] Buyer Buyer Buyer Non-buyer Non-buyer Non-buyer Non-buyer
## [8] Buyer Non-buyer Non-buyer
## attr(,"prob")
## [1] 1.0000000 0.6000000 0.8333333 0.6000000 0.6000000 0.6666667 0.6000000
## [8] 0.8333333 0.6000000 0.6000000
## Levels: Buyer Non-buyerRedes Neuronales
Las cajas negras
#install.packages("nnet")
library(nnet)
library(caret)cargamos los datos
bn<-read.csv("../../r-course-master/data/tema3/banknote-authentication.csv")
bn$class<-factor(bn$class)# solo para corroborar que lo toma como factor Para seleccionar el rango que se va a utilizar dentro de la red neuronal sigue a grosomodo la fomrmula rango*maximo (|variables|)~1
apply(bn,2,max)## variance skew curtosis entropy class
## " 6.8248000" " 12.95160000" "17.92740000" " 2.4495000" "1"el maximo de las varibles que vemos es 17, por lo tanto se debe buscar un numero que se acerque a uno
podemos hacer una ecuacion rango= 1/ maximo de los numeros que nos mostro la funcion apply
Creamos particion
t.id.redn<-createDataPartition(bn$class,p=0.70,list=F)
mod.redn<-nnet::nnet(class~.,data=bn[t.id.redn,],
size=3,# especifica el numero de capas , puede ser el numero de columnas
maxit=10000,#Numero de iteraciones maximas
decay=.001,#Se utiliza para el overfiting ( no se especializa en los datos que entrenamos )se debe dejar rango para las validacioens
rang=0.05,#ver nota superior
na.action=na.omit)#Para que omita los NA## # weights: 19
## initial value 665.779419
## iter 10 value 445.463295
## iter 20 value 61.023803
## iter 30 value 32.802513
## iter 40 value 17.213270
## iter 50 value 2.682251
## iter 60 value 2.006200
## iter 70 value 1.788214
## iter 80 value 1.581753
## iter 90 value 1.522893
## iter 100 value 1.477252
## iter 110 value 1.454335
## iter 120 value 1.354654
## iter 130 value 1.258758
## iter 140 value 1.191089
## iter 150 value 1.012912
## iter 160 value 0.922466
## iter 170 value 0.855131
## iter 180 value 0.825344
## iter 190 value 0.752956
## iter 200 value 0.662098
## iter 210 value 0.632855
## iter 220 value 0.626612
## iter 230 value 0.625244
## iter 240 value 0.624736
## iter 250 value 0.624530
## iter 260 value 0.624512
## iter 270 value 0.624510
## iter 280 value 0.624509
## iter 290 value 0.624509
## final value 0.624509
## convergedHagamos la prediccion
pred.redn<-predict(mod.redn,newdata = bn[-t.id.redn,],type = "class")
table(bn[-t.id.redn,]$class,pred.redn,dnn = c("Actual","Predicho"))## Predicho
## Actual 0 1
## 0 228 0
## 1 0 183#curva Roc de una red neuronal
library(ROCR)
pred2.redn<-predict(mod.redn,newdata = bn[-t.id.redn,],type = "raw")
perf.redn<-performance(prediction(pred2.redn,bn[-t.id.redn,"class"]),"tpr","fpr")
plot(perf.redn)
Esta madre es perfecta
Analisis del discriminante lineal(lda)
Se trata de ver si las medias o varianzas son parecidas entre ellas Se parece mucho a la prueba de ANoVA la varible dependiente :categorica la varible independiente: continua
#install.packages("MASS")
library(MASS)##
## Attaching package: 'MASS'## The following object is masked _by_ '.GlobalEnv':
##
## housing## The following objects are masked from 'package:raster':
##
## area, select## The following object is masked from 'package:dplyr':
##
## selectlibrary(caret)cargamos datos
#seet.seed(2018)
bn.disc.lin<- read.csv("../../r-course-master/data/tema3/banknote-authentication.csv")
bn.disc.lin$class<-factor(bn.disc.lin$class)
t.id.disc.lin<-createDataPartition(bn.disc.lin$class,p=0.7,list=F)modelo lineal de dicriminacion
mod.disc.lin<-lda(bn.disc.lin[t.id.disc.lin,1:4],#parametro de entrada
bn.disc.lin[t.id.disc.lin,5] )#parametro de salida , con la indicación de la columna
# otra forma de escribirlo
# mod.disc.lin.2<-lda(class~.,data = bn.disc.lin[t.id.disc.lin,])agregamos una columana al DF y creamos la tabla para los datos de entrenamiento
bn.disc.lin[t.id.disc.lin,"Pred"]<-predict(mod.disc.lin,bn.disc.lin[t.id.disc.lin,1:4])$class# valor predicho y se mete desde el una columna nueva
table(bn.disc.lin[t.id.disc.lin,"class"],
bn.disc.lin[t.id.disc.lin,"Pred"],dnn = c("Actual","Predicho"))## Predicho
## Actual 0 1
## 0 515 19
## 1 0 427Excluimos los valores de etrenamiento y predecimos la matriz de confusion
bn.disc.lin[-t.id.disc.lin,"Pred"]<-predict(mod.disc.lin,bn.disc.lin[-t.id.disc.lin,1:4])$class
table(bn.disc.lin[-t.id.disc.lin,"class"],
bn.disc.lin[-t.id.disc.lin,"Pred"],dnn = c("Actual","Predicho"))## Predicho
## Actual 0 1
## 0 215 13
## 1 0 183Regresión logistica (gml)
Resultado de una variable categorica en fucnión de variables predictoras calculo integral (supongo resuelve las integrales a travez de metodos numericos )
library(caret)
library(lattice)
#cargamos datos
boston_Reg_log<-read.csv("../../r-course-master/data/tema3/boston-housing-logistic.csv")
boston_Reg_log$CLASS<-factor(boston_Reg_log$CLASS,levels=c(0,1))
set.seed(2018)
# creamos la particion aleatoria
t_id_reg_log<-createDataPartition(boston_Reg_log$CLASS,p=0.7,list = F)
#creamos el modelo
mod_reg_log<- glm(CLASS~.,data=boston_Reg_log[t_id_reg_log,],family=binomial)#si es vaerdadero o falso
summary(mod_reg_log)##
## Call:
## glm(formula = CLASS ~ ., family = binomial, data = boston_Reg_log[t_id_reg_log,
## ])
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.17116 -0.47682 0.09069 0.44375 2.83196
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 23.873324 3.936114 6.065 1.32e-09 ***
## NOX -21.564510 4.658177 -4.629 3.67e-06 ***
## DIS -0.368276 0.143892 -2.559 0.0105 *
## RAD 0.174650 0.071923 2.428 0.0152 *
## TAX -0.005127 0.003755 -1.365 0.1722
## PTRATIO -0.717774 0.133575 -5.374 7.72e-08 ***
## B 0.009092 0.003918 2.320 0.0203 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 353.03 on 254 degrees of freedom
## Residual deviance: 167.07 on 248 degrees of freedom
## AIC: 181.07
##
## Number of Fisher Scoring iterations: 6Investigar fisher scoring
Creamos una columana de datos de las probabilidades
boston_Reg_log[-t_id_reg_log,"Prob_success"]<-predict(mod_reg_log,newdata = boston_Reg_log[-t_id_reg_log,],
type="response")# mide las probabilidades de los objetos clasificados
# creamos otra columna donde si esta por encima de 50% de prob, entonces es un exito , de lo contrario es un fracaso
boston_Reg_log[-t_id_reg_log,"pred_50"]<-ifelse(boston_Reg_log[-t_id_reg_log,"Prob_success"]>=0.5,1,0)Creamos la matriz de confusion para la regresion logistica
table(boston_Reg_log[-t_id_reg_log,"CLASS"],boston_Reg_log[-t_id_reg_log,"pred_50"],dnn=c("Actual","Predicho"))## Predicho
## Actual 0 1
## 0 45 6
## 1 6 51Analizar en tiempo real de TWiter
Se deja pendiente por validacion de perfil de desarrollados por parte de Twiter
#La raiz del error cuadratico medio Se trata de restar Para este ejemplo se toma el precio=price y lo predicho por alguno de los algotirmos antes vistos
dat_error_medio<-read.csv("../../r-course-master/data/tema4/rmse.csv")
# Paso 1 se resta la columna del precio menos la columana predicha y este resultado al cuadrado
sqrt(mean((dat_error_medio$price-dat_error_medio$pred)^2))## [1] 2.934995#graficamos los datos de prediccion vs precio real
plot(dat_error_medio$price,dat_error_medio$pred,xlab = "Actual",ylab = "Predicho")
abline(0,1)
Hagamos una funcion para error cuadratico medio
rmse.fuction<-function(actual,predicted){
return(sqrt(mean((actual-predicted)^2)))
}Como se usa esta vaina
rmse.fuction(dat_error_medio$price,dat_error_medio$pred)## [1] 2.934995regresion K nearest neiighbors
#install.packages("FNN")
library(FNN)##
## Attaching package: 'FNN'## The following objects are masked from 'package:class':
##
## knn, knn.cvlibrary(dummies)
library(scales)
library(caret)
edu<-read.csv("../../r-course-master/data/tema4/education.csv")Se va a hacer dummy la variable región
dms<-dummy(edu$region,sep="_")## Warning in model.matrix.default(~x - 1, model.frame(~x - 1), contrasts = FALSE):
## non-list contrasts argument ignoredLe agregamos esta varibel al DF edu
edu<-cbind(edu,dms)Cambiemos a factores los datos que son necesarios estandarizar el resto de
edu$urban.s<-rescale(edu$urban)
edu$income.z<-rescale(edu$income)
edu$under18.z<-rescale(edu$under18)Creamos la perticion de los datos
set.seed(2018)
library(caret)
t_id_kneigh.regresion<-caret::createDataPartition(edu$expense,p=0.6, list=F)conjunto de entrenamiento
tr_knn_regresion<-edu[t_id_kneigh.regresion,]# meten todos lo que fueron seleccionados
temp<-edu[-t_id_kneigh.regresion,] #Aqui entra el resto de los datosDatos de validacion
v_ids_kneigh.regresion<-createDataPartition(temp$expense,p=0.5,list=F)#Se Creo otra partición de los datos de temp
val.kneigh.regresion<-temp[v_ids_kneigh.regresion,]
test_kneight_regresion<-temp[-v_ids_kneigh.regresion,]Determinemos la K que es mas eficiente
training_control<-trainControl(method = "repeatedcv",#Valodacion cruzada repetida
number = 10,#numero de remuestreos
repeats = 3)#repeticiones de la validacion cruzada
cater.knn.fit<-train(expense~.,data =edu,
method="knn",#metodo por el cual se va a hacer
trControl=training_control,#el de arriba
tuneLength=10)#longitud del ajuste
cater.knn.fit## k-Nearest Neighbors
##
## 50 samples
## 12 predictors
##
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 3 times)
## Summary of sample sizes: 45, 46, 43, 45, 44, 46, ...
## Resampling results across tuning parameters:
##
## k RMSE Rsquared MAE
## 5 47.66062 0.4239637 37.55211
## 7 46.00489 0.4819194 36.71542
## 9 46.63363 0.4587353 36.98325
## 11 45.77538 0.4655686 36.77972
## 13 46.30492 0.4597074 36.87004
## 15 45.32533 0.4888495 36.30289
## 17 45.44537 0.4919484 36.72798
## 19 46.36769 0.4783885 37.50904
## 21 46.42500 0.4896844 37.61464
## 23 46.60530 0.4987181 37.80011
##
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was k = 15.Regresion 1
reg1<-knn.reg(tr_knn_regresion[,7:12],# conjunto de entrenamiento
val.kneigh.regresion[,7:12],#Conjunto de validacion
tr_knn_regresion$expense, # variable ( columna) que se desea predecir
k=25, # iterar para el valor mas cercano
algorithm = "brute")#este es por defaultregrecion predicha
reg.test<-knn.reg(tr_knn_regresion[,7:12],test_kneight_regresion[,7:12],tr_knn_regresion$expense,k=25,algorithm = "brute")Cargamos la funcion que en la sesión pasada generamos
rmse_reg_kneight<-rmse.fuction(reg1$pred,edu$expense)
rmse_reg_kneight## [1] 61.36413Regresion K nearest neiighbors sin particiones de test
creamos la regresion de donde se van a tomar los valores sin test
reg_knn_sin.test<-knn.reg(tr_knn_regresion[,7:12],#valores a de entrenamiento
test = NULL,# sin valores de prueba
y=tr_knn_regresion$expense,#respuesta de cada observacion de entrenamiento
k=3,# numero de vecinos
algorithm = "brute")Error medio
rmse_reg_knn<- sqrt(mean((reg_knn_sin.test$residuals)^2))
rmse_reg_knn## [1] 40.9394funcion que calcula la regresion lineal
regresion_knn<-function(tr_predictor,val_predictors,
tr_target,val_target,k) {
library(FNN)
result<-knn.reg(tr_predictor,val_predictors,tr_target,k,algorithm = "brute")
rmserror<-sqrt(mean((tr_target-result$pred)^2))
cat(paste("RMSE para k=",toString(k),":",rmserror,"\n",sep=" "))
rmserror
}Como se usa esta funcion de regresion knn con particiones de validacion (test)
regresion_knn(tr_knn_regresion[,7:12],val.kneigh.regresion[,7:12]
,tr_knn_regresion$expense,val.kneigh.regresion$expense,25)## Warning in tr_target - result$pred: longitud de objeto mayor no es múltiplo de
## la longitud de uno menor## RMSE para k= 25 : 48.4137955545731## [1] 48.4138Funcion que grafica los valores de error medio con base a lo calculado anteriormente RMSE
error_knn_cuadratico_regresion <-function(tr_predictors, val_predictors,
tr_target, val_target, start_k, end_k){
rms_errors <- vector()
for(k in start_k:end_k){
rms_error <- regresion_knn(tr_predictors, val_predictors,
tr_target, val_target, k)
rms_errors <- c(rms_errors, rms_error)
}
plot(rms_errors, type = 'o', xlab = "k", ylab = "RMSE")
}como se usa esta pinche funcion de RMSE
error_knn_cuadratico_regresion(tr_knn_regresion[,7:12],val.kneigh.regresion[,7:12],tr_knn_regresion$expense,val.kneigh.regresion$expense,1,25)## Warning in tr_target - result$pred: longitud de objeto mayor no es múltiplo de
## la longitud de uno menor## RMSE para k= 1 : 69.0455013016779## Warning in tr_target - result$pred: longitud de objeto mayor no es múltiplo de
## la longitud de uno menor## RMSE para k= 2 : 64.9357855885335## Warning in tr_target - result$pred: longitud de objeto mayor no es múltiplo de
## la longitud de uno menor## RMSE para k= 3 : 66.3299569911246## Warning in tr_target - result$pred: longitud de objeto mayor no es múltiplo de
## la longitud de uno menor## RMSE para k= 4 : 64.8152572942899## Warning in tr_target - result$pred: longitud de objeto mayor no es múltiplo de
## la longitud de uno menor## RMSE para k= 5 : 63.6186096987352## Warning in tr_target - result$pred: longitud de objeto mayor no es múltiplo de
## la longitud de uno menor## RMSE para k= 6 : 62.3610608450676## Warning in tr_target - result$pred: longitud de objeto mayor no es múltiplo de
## la longitud de uno menor## RMSE para k= 7 : 60.011111853709## Warning in tr_target - result$pred: longitud de objeto mayor no es múltiplo de
## la longitud de uno menor## RMSE para k= 8 : 59.4158371077842## Warning in tr_target - result$pred: longitud de objeto mayor no es múltiplo de
## la longitud de uno menor## RMSE para k= 9 : 57.8842635174964## Warning in tr_target - result$pred: longitud de objeto mayor no es múltiplo de
## la longitud de uno menor## RMSE para k= 10 : 57.3507029163549## Warning in tr_target - result$pred: longitud de objeto mayor no es múltiplo de
## la longitud de uno menor## RMSE para k= 11 : 57.1439027354162## Warning in tr_target - result$pred: longitud de objeto mayor no es múltiplo de
## la longitud de uno menor## RMSE para k= 12 : 56.6673138749805## Warning in tr_target - result$pred: longitud de objeto mayor no es múltiplo de
## la longitud de uno menor## RMSE para k= 13 : 55.613860163701## Warning in tr_target - result$pred: longitud de objeto mayor no es múltiplo de
## la longitud de uno menor## RMSE para k= 14 : 54.0501866064668## Warning in tr_target - result$pred: longitud de objeto mayor no es múltiplo de
## la longitud de uno menor## RMSE para k= 15 : 53.570560945355## Warning in tr_target - result$pred: longitud de objeto mayor no es múltiplo de
## la longitud de uno menor## RMSE para k= 16 : 52.7946669226045## Warning in tr_target - result$pred: longitud de objeto mayor no es múltiplo de
## la longitud de uno menor## RMSE para k= 17 : 52.3436748243534## Warning in tr_target - result$pred: longitud de objeto mayor no es múltiplo de
## la longitud de uno menor## RMSE para k= 18 : 51.0961651607614## Warning in tr_target - result$pred: longitud de objeto mayor no es múltiplo de
## la longitud de uno menor## RMSE para k= 19 : 51.0052869862921## Warning in tr_target - result$pred: longitud de objeto mayor no es múltiplo de
## la longitud de uno menor## RMSE para k= 20 : 51.0679594694658## Warning in tr_target - result$pred: longitud de objeto mayor no es múltiplo de
## la longitud de uno menor## RMSE para k= 21 : 49.8425056400732## Warning in tr_target - result$pred: longitud de objeto mayor no es múltiplo de
## la longitud de uno menor## RMSE para k= 22 : 49.1411825647657## Warning in tr_target - result$pred: longitud de objeto mayor no es múltiplo de
## la longitud de uno menor## RMSE para k= 23 : 48.5951736558941## Warning in tr_target - result$pred: longitud de objeto mayor no es múltiplo de
## la longitud de uno menor## RMSE para k= 24 : 48.0292263769835## Warning in tr_target - result$pred: longitud de objeto mayor no es múltiplo de
## la longitud de uno menor## RMSE para k= 25 : 48.4137955545731
#Regresion lineal Es
library(caret)
auto <- read.csv("../../r-course-master/data/tema4/auto-mpg.csv")
auto$cylinders <- factor(auto$cylinders,
levels = c(3,4,5,6,8),
labels = c("3c", "4c", "5c", "6c", "8c"))Realizamos la particion
set.seed(2018)#muestreo
t_id_reg_lin<- createDataPartition(auto$mpg, p = 0.7, list = F)# creamos particion
names(auto)#imprimimos el nombre de las columnas## [1] "No" "mpg" "cylinders" "displacement" "horsepower"
## [6] "weight" "acceleration" "model_year" "car_name"Creamos el modelo Ya hace las varibles dummy de por default la funcion lm
mod <- lm(mpg ~ .#mpg en funcion de todos
, data = auto[t_id_reg_lin,# con los datos de entrenamieto
-c(1,8,9)])#quitamos las columnas que no necesitamos(solo la 1,8 y 9)
mod##
## Call:
## lm(formula = mpg ~ ., data = auto[t_id_reg_lin, -c(1, 8, 9)])
##
## Coefficients:
## (Intercept) cylinders4c cylinders5c cylinders6c cylinders8c
## 37.284202 6.231475 8.248195 2.131026 4.568171
## displacement horsepower weight acceleration
## 0.002245 -0.057543 -0.004665 0.050745Este es el modelo que predice el Millas Por galeon #mpg = 37.284202 + # + 6.2314754c + 8.2481955c + 2.1310266c + 4.5681718c + # + 0.002245 * displacement -0.057543 * horsepower + # - 0.004665 * weight + 0.050745 * acceleration
Residual es el error que estoy cometiendo
Recuerda que los * que pinta summary (son mas importantes para el modelo lineal) residual standar error =Rmse Si el p-value es chiquito es mejor la estimacion
summary(mod)##
## Call:
## lm(formula = mpg ~ ., data = auto[t_id_reg_lin, -c(1, 8, 9)])
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.0606 -2.4686 -0.4435 1.9821 16.0907
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.2842024 3.6497412 10.216 < 2e-16 ***
## cylinders4c 6.2314753 2.4926855 2.500 0.01301 *
## cylinders5c 8.2481946 3.8091396 2.165 0.03123 *
## cylinders6c 2.1310256 2.7759570 0.768 0.44335
## cylinders8c 4.5681710 3.2054454 1.425 0.15527
## displacement 0.0022449 0.0108924 0.206 0.83687
## horsepower -0.0575428 0.0202773 -2.838 0.00489 **
## weight -0.0046652 0.0009999 -4.665 4.84e-06 ***
## acceleration 0.0507454 0.1443575 0.352 0.72547
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.092 on 271 degrees of freedom
## Multiple R-squared: 0.7304, Adjusted R-squared: 0.7224
## F-statistic: 91.75 on 8 and 271 DF, p-value: < 2.2e-16Haremos una caja de gato con bigotes de los residuos (de los errores)
boxplot(mod$residuals)
RMSE, aunque le summary lo arroja lo podemos comparar
sqrt(mean((mod$fitted.values - auto[t_id_reg_lin,]$mpg)^2))## [1] 4.026021Predecimos los valores de consumo de Millas por Galeon
pred <- predict(mod, auto[-t_id_reg_lin, -c(1,8,9)])
sqrt(mean((pred - auto[-t_id_reg_lin,]$mpg)^2))## [1] 3.894627Interpretemos los parametros
par(mfrow=c(2,2))#recuerda que esto es para graficar 2 por 2
plot(mod)
Analicemos los modelos
los numeros de los plots son el numero de dato de la izquierda del DF original
1._se contesta la pregunta ,Existe una relacion que es lineal? 2._El de la parte superior derecha ( sirve para ver si los errores siguen una distribucion normal) 3._ El de la parte inferior izquierda(Representa los residuos estandarizados sobre los valores ajustados)spred location ( si todos los valores tienen la misma varianza “homosedasticos”)si estan ordenados , es bueno
4._la de la derecha abajo Posibles sujetos influyentes (los outliners pueden influenciar los valores )(distancia de cook’s distance) arriba derecha y abajo derecha son los puntos criticos
Si se quiere indicar al modelo que utilice un dato especifico , se le pude que sea otro el de referencia
Nota importante : el que tengas mas valores categoricos se recomienda ponerlo como principal
auto <- within(auto,
cylinders <- relevel(cylinders, ref="4c"))
# se vuelve a correr el modelo
mod <- lm(mpg ~. , data = auto[t_id_reg_lin, -c(1,8,9)])
summary(mod)##
## Call:
## lm(formula = mpg ~ ., data = auto[t_id_reg_lin, -c(1, 8, 9)])
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.0606 -2.4686 -0.4435 1.9821 16.0907
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 43.5156778 2.8644210 15.192 < 2e-16 ***
## cylinders3c -6.2314753 2.4926855 -2.500 0.013014 *
## cylinders5c 2.0167192 2.9414057 0.686 0.493532
## cylinders6c -4.1004497 1.0580642 -3.875 0.000134 ***
## cylinders8c -1.6633044 1.8447503 -0.902 0.368048
## displacement 0.0022449 0.0108924 0.206 0.836872
## horsepower -0.0575428 0.0202773 -2.838 0.004886 **
## weight -0.0046652 0.0009999 -4.665 4.84e-06 ***
## acceleration 0.0507454 0.1443575 0.352 0.725467
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.092 on 271 degrees of freedom
## Multiple R-squared: 0.7304, Adjusted R-squared: 0.7224
## F-statistic: 91.75 on 8 and 271 DF, p-value: < 2.2e-16pred <- predict(mod, auto[-t_id_reg_lin, -c(1,8,9)])
sqrt(mean((pred-auto[-t_id_reg_lin,]$mpg)^2))## [1] 3.894627par(mfrow=c(2,2))
plot(mod)
Se puede mejorar el analisis , esto sirve para determinar que columnas (variables )son importates y cuales no (stepAIC)
library(MASS)
mod##
## Call:
## lm(formula = mpg ~ ., data = auto[t_id_reg_lin, -c(1, 8, 9)])
##
## Coefficients:
## (Intercept) cylinders3c cylinders5c cylinders6c cylinders8c
## 43.515678 -6.231475 2.016719 -4.100450 -1.663304
## displacement horsepower weight acceleration
## 0.002245 -0.057543 -0.004665 0.050745summary(mod)##
## Call:
## lm(formula = mpg ~ ., data = auto[t_id_reg_lin, -c(1, 8, 9)])
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.0606 -2.4686 -0.4435 1.9821 16.0907
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 43.5156778 2.8644210 15.192 < 2e-16 ***
## cylinders3c -6.2314753 2.4926855 -2.500 0.013014 *
## cylinders5c 2.0167192 2.9414057 0.686 0.493532
## cylinders6c -4.1004497 1.0580642 -3.875 0.000134 ***
## cylinders8c -1.6633044 1.8447503 -0.902 0.368048
## displacement 0.0022449 0.0108924 0.206 0.836872
## horsepower -0.0575428 0.0202773 -2.838 0.004886 **
## weight -0.0046652 0.0009999 -4.665 4.84e-06 ***
## acceleration 0.0507454 0.1443575 0.352 0.725467
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.092 on 271 degrees of freedom
## Multiple R-squared: 0.7304, Adjusted R-squared: 0.7224
## F-statistic: 91.75 on 8 and 271 DF, p-value: < 2.2e-16step.model <- MASS::stepAIC(mod, direction="backward")## Start: AIC=797.96
## mpg ~ cylinders + displacement + horsepower + weight + acceleration
##
## Df Sum of Sq RSS AIC
## - displacement 1 0.71 4539.2 796.00
## - acceleration 1 2.07 4540.5 796.08
## <none> 4538.5 797.96
## - horsepower 1 134.87 4673.3 804.16
## - weight 1 364.53 4903.0 817.59
## - cylinders 4 619.82 5158.3 825.80
##
## Step: AIC=796
## mpg ~ cylinders + horsepower + weight + acceleration
##
## Df Sum of Sq RSS AIC
## - acceleration 1 1.76 4541.0 794.11
## <none> 4539.2 796.00
## - horsepower 1 136.06 4675.2 802.27
## - weight 1 475.09 5014.3 821.87
## - cylinders 4 629.80 5169.0 824.38
##
## Step: AIC=794.11
## mpg ~ cylinders + horsepower + weight
##
## Df Sum of Sq RSS AIC
## <none> 4541.0 794.11
## - horsepower 1 289.06 4830.0 809.39
## - cylinders 4 642.38 5183.3 823.16
## - weight 1 605.26 5146.2 827.14summary(step.model)##
## Call:
## lm(formula = mpg ~ cylinders + horsepower + weight, data = auto[t_id_reg_lin,
## -c(1, 8, 9)])
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.9108 -2.4854 -0.4535 1.9444 16.2137
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 44.3060838 1.5539812 28.511 < 2e-16 ***
## cylinders3c -6.4358865 2.3970441 -2.685 0.0077 **
## cylinders5c 2.0004828 2.9289747 0.683 0.4952
## cylinders6c -4.0112530 0.8426183 -4.760 3.14e-06 ***
## cylinders8c -1.4973133 1.3818662 -1.084 0.2795
## horsepower -0.0611080 0.0146588 -4.169 4.12e-05 ***
## weight -0.0044161 0.0007321 -6.032 5.24e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.078 on 273 degrees of freedom
## Multiple R-squared: 0.7302, Adjusted R-squared: 0.7243
## F-statistic: 123.2 on 6 and 273 DF, p-value: < 2.2e-16Lo que te arroja los variables de la ultima iteracion por lo que se puede concluir que (cylinders, hoursepower y wheight) son las varibles para la obtener un modelo completo
Regrecion con Arboles de clasificacion
Instalamos y cargamos las librerias
#install.packages("rpart.plot")
library(rpart)
library(rpart.plot)
library(caret)Cargamos el documento ( cambiar dirección)
bh <- read.csv("../../r-course-master/data/tema4/BostonHousing.csv")Creamos una paticion
t.id <-caret::createDataPartition(bh$MEDV, p = .7, list = F)bfit <- rpart(MEDV ~., data = bh[t.id,])
bfit## n= 356
##
## node), split, n, deviance, yval
## * denotes terminal node
##
## 1) root 356 28822.1900 22.40618
## 2) RM< 6.805 293 11281.7400 19.51433
## 4) LSTAT>=14.91 112 1896.9530 14.34196
## 8) CRIM>=5.76921 58 591.3710 11.96552 *
## 9) CRIM< 5.76921 54 626.2083 16.89444 *
## 5) LSTAT< 14.91 181 4534.2900 22.71492
## 10) LSTAT>=5.245 172 2639.5900 22.13140
## 20) LSTAT>=9.715 89 566.7760 20.60674 *
## 21) LSTAT< 9.715 83 1644.0860 23.76627 *
## 11) LSTAT< 5.245 9 716.8800 33.86667 *
## 3) RM>=6.805 63 3694.3560 35.85556
## 6) RM< 7.437 44 1094.5590 32.39545 *
## 7) RM>=7.437 19 853.1011 43.86842 *Grafiquemoslo
par(mfrow=c(1,1))
prp(bfit, type = 2, nn=T,
fallen.leaves = T, faclen = 4,
varlen = 8, shadow.col = "gray")
Detalle de tabla
bfit$cptable## CP nsplit rel error xerror xstd
## 1 0.48039697 0 1.0000000 1.0033907 0.09825971
## 2 0.16829040 1 0.5196030 0.5827505 0.06549191
## 3 0.06060246 2 0.3513126 0.3935690 0.05315030
## 4 0.04086502 3 0.2907102 0.3487069 0.05377561
## 5 0.02357120 4 0.2498451 0.3007673 0.04694361
## 6 0.01487496 5 0.2262740 0.2785457 0.04780752
## 7 0.01000000 6 0.2113990 0.2642192 0.04714776Grafiquemos lo que dice la tabla
cuando se tiene muchas clasificaciones se puede elejir alguno de los factores que sea representativo
plotcp(bfit)
Subes a la tablay tomas el valor que sea mas factible (y recortas el arbol)
bfitpruned <- prune(bfit, cp = 0.02343779)Volvemos a representar el arbol
par(mfrow=c(1,1))
prp(bfitpruned, type = 2, nn=T,
fallen.leaves = T, faclen = 4,
varlen = 8, shadow.col = "gray")
Calculamos las predicciones
Pedeciremos los datos de entrenamiento para el arbol no recortado, todos los datos de entrenamiento
preds <- predict(bfit, bh[t.id,])
#Calculamos el error medio para los datos de validacion
sqrt(mean((preds - bh[t.id,]$MEDV)^2))## [1] 4.137042Pedeciremos los datos deentrenamiento para el arbol no recortado, todos los datos que no son de entrenamiento
preds <- predict(bfit, bh[-t.id, ])
sqrt(mean((preds - bh[-t.id,]$MEDV)^2))## [1] 5.689694Pedeciremos los datos de entrenamiento para el arbol recortado, todos los datos entrenamiento
preds <- predict(bfitpruned, bh[t.id,])
sqrt(mean((preds - bh[t.id,]$MEDV)^2))## [1] 4.280118Pedeciremos los datos de validacion para el arbol recortado, todos los datos que no son de entrenamiento
preds <- predict(bfitpruned, bh[-t.id, ])
sqrt(mean((preds - bh[-t.id,]$MEDV)^2))## [1] 5.898846#predictores categoricos
#Cargamos los datos
ed <- read.csv("../../r-course-master/data/tema4/education.csv")
#convertimos a factores
ed$region <- factor(ed$region)
#partimos el data frame
t.id <- caret::createDataPartition(ed$expense, p = 0.7, list = F)
#Predeciren función de todo la variable
fit <- rpart::rpart(expense ~ region+urban+income+under18, data = ed[t.id,])
#pintemos el arbol que diseñamos
rpart.plot::prp(fit, type = 2, nn=T, fallen.leaves = T, faclen=4, varlen=8, shadow.col = "gray")
#Bagging Combina de forma conjunta lo de diferentes arboles
#install.packages("ipred")
library(ipred)## Warning: package 'ipred' was built under R version 4.1.3##
## Attaching package: 'ipred'## The following object is masked from 'package:raster':
##
## cv# Construimos el bagguing a partir de los datos de entrenamiento
bagging.fit <- bagging(MEDV~.,#predice esta variable en función de las demas
data=bh[t.id, ])# vas a tomar del DF inicial y las filas de entrenamiento
# Predeciremos apartir de el conjunto de entrenamiento
prediction.t <- predict(bagging.fit, bh[t.id,])
#RMSE de bagging
sqrt(mean((prediction.t-bh[t.id,]$MEDV)^2))## [1] 2.561796#Predecirmos con los valores a partir de los de validación ( los que nos on de entrenamiento)
prediction.v <- predict(bagging.fit, bh[-t.id,])
sqrt(mean((prediction.v-bh[-t.id,]$MEDV)^2))## [1] 6.936346#Boosting
#install.packages("gbm")
library(gbm)#gradient busted machine ## Warning: package 'gbm' was built under R version 4.1.3## Loaded gbm 2.1.8gbmfit <- gbm(MEDV~., #predice en funciónde todas las variables
data = bh, #datos tomalos de aqui
distribution = "gaussian")# distribucion gaussina
#Predice en función de
prediction.t <- predict(gbmfit, bh)## Using 100 trees...Random Forest para regresiones
library(randomForest)
library(caret)
bh <- read.csv("../../r-course-master/data/tema4/BostonHousing.csv")
#elejimos semilla
set.seed(2018)
#Creamos datos de entrenamiento
t.id <- createDataPartition(bh$MEDV, p = 0.7, list = F)Creamos el modelo
mod <- randomForest(x = bh[t.id, 1:13],#predictores
y = bh[t.id, 14],#las variables a predecir
ntree = 1000, #numero de arboles
xtest = bh[-t.id, 1:13],#son los predictors que se usan en la particion
ytest = bh[-t.id, 14],#las variables de salida en la misma particio de validacion
importance = T,#Si tiene que computar o no las puntuaciones de varibales predictoras
keep.forest = T)#mantener los arboles del modelo
mod##
## Call:
## randomForest(x = bh[t.id, 1:13], y = bh[t.id, 14], xtest = bh[-t.id, 1:13], ytest = bh[-t.id, 14], ntree = 1000, importance = T, keep.forest = T)
## Type of random forest: regression
## Number of trees: 1000
## No. of variables tried at each split: 4
##
## Mean of squared residuals: 11.77997
## % Var explained: 86.72
## Test set MSE: 10.05
## % Var explained: 86.51Parametros que se pueden introducir en la tecnica de Random Forest
mtry = m/3, donde m = # de predictores nodesize = 5 minimos tamaño que debe tener un nodo termina maxnodes ( maximo de nodos terminales )
Se ven la importancia para cada variable %purity el dato mayor es la que tine mayor importancia
mod$importance## %IncMSE IncNodePurity
## CRIM 8.1353944 1795.0927
## ZN 0.9797628 240.7091
## INDUS 4.7049920 1536.7625
## CHAS 0.5831449 194.5591
## NOX 8.1279708 1623.9652
## RM 41.3092319 10179.2089
## AGE 3.8170358 802.6095
## DIS 6.4509516 1953.6404
## RAD 1.0707028 263.1413
## TAX 4.4276314 1188.9411
## PTRATIO 8.9288034 2591.7216
## B 1.8193482 667.7755
## LSTAT 55.9722090 8001.2091Grafiquemoslo soble los datos predichos ( con los datos de entrenamiento)
plot(bh[t.id,]$MEDV, predict(mod, newdata = bh[t.id,]),
xlab = "Actual", ylab = "Predichos")
abline(0,1)
Grafiquemos los datos en caso de los datos de predicccion
plot(bh[-t.id,]$MEDV, predict(mod, newdata = bh[-t.id,]),
xlab = "Actual", ylab = "Predichos")
abline(0,1)
Regresion con Redes neuronales
library(nnet)
library(caret)
library(devtools)## Warning: package 'devtools' was built under R version 4.1.3## Loading required package: usethis## Warning: package 'usethis' was built under R version 4.1.3bh <- read.csv("../../r-course-master/data/tema4/BostonHousing.csv")
set.seed(2018)
t.id <- createDataPartition(bh$MEDV, p= 0.7, list = F)
summary(bh$MEDV)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5.00 17.02 21.20 22.53 25.00 50.00creamos la red neuronal ( parara en 290 iteraciones)
fit <- nnet(MEDV/50 ~., data=bh[t.id, ],
size = 6, decay = 0.1,
maxit = 1000, linout=T)## # weights: 91
## initial value 1711.376818
## iter 10 value 14.279679
## iter 20 value 9.849751
## iter 30 value 7.511037
## iter 40 value 6.285160
## iter 50 value 3.989395
## iter 60 value 3.214674
## iter 70 value 2.455499
## iter 80 value 1.950545
## iter 90 value 1.782552
## iter 100 value 1.735083
## iter 110 value 1.686014
## iter 120 value 1.632630
## iter 130 value 1.516829
## iter 140 value 1.490113
## iter 150 value 1.438664
## iter 160 value 1.385513
## iter 170 value 1.362224
## iter 180 value 1.349518
## iter 190 value 1.322281
## iter 200 value 1.307838
## iter 210 value 1.304445
## iter 220 value 1.303188
## iter 230 value 1.302788
## iter 240 value 1.302763
## iter 250 value 1.302370
## iter 260 value 1.301663
## iter 270 value 1.300696
## iter 280 value 1.298021
## iter 290 value 1.293863
## iter 300 value 1.290511
## iter 310 value 1.279228
## iter 320 value 1.264173
## iter 330 value 1.241595
## iter 340 value 1.211158
## iter 350 value 1.193082
## iter 360 value 1.187331
## iter 370 value 1.184634
## iter 380 value 1.181169
## iter 390 value 1.179331
## iter 400 value 1.177550
## iter 410 value 1.170499
## iter 420 value 1.169408
## iter 430 value 1.159805
## iter 440 value 1.149498
## iter 450 value 1.147599
## iter 460 value 1.147367
## final value 1.147343
## convergedDescargueos el recurso
source_url("https://gist.githubusercontent.com/fawda123/7471137/raw/466c1474d0a505ff044412703516c34f1a4684a5/nnet_plot_update.r")## i SHA-1 hash of file is 74c80bd5ddbc17ab3ae5ece9c0ed9beb612e87efGrafiquemos la red neuronal
plot(fit, max.sp = T)
Calculemos el valor medio de la red neuronal de los datos de entrenamiento
sqrt(mean((fit$fitted.values*50# esto se hace por que estaban escalados (x/50)
-bh[t.id,"MEDV"])^2))## [1] 2.354379Predigamos los valores de validación con base al modelo
# calculamos la predicción
pred <- predict(fit, bh[-t.id,])
# calculemos el error medio
sqrt(mean((pred*50 - bh[-t.id,"MEDV"])^2))## [1] 4.082578#k-fold cross validation( para evitar el over-fiting ) La técnica que hace un muestreo de los datos dejando siempre un conjunto no vacío y de más de un elemento como testing de los mismos usando el resto como conjunto de entrenamiento, de modo que en algún momento todo elemento del data set habrá sido parte del conjunto de testing
bh <- read.csv("../../r-course-master/data/tema4/BostonHousing.csv")Función que nos ayuda de validacion cruzada (Kross Validaion )
kfold.crossval.reg <- function(df, nfolds){
fold <- sample(1:nfolds,#diviciones
nrow(df),#valor maximo es el numero de columnas
replace = T)# con remplazo
cat(fold)
mean.sqr.errs <- sapply(1:nfolds,
kfold.cval.reg.iter,
df, fold)
list("MSE "= mean.sqr.errs,
"Overall_Mean_Sqr_Error" = mean(mean.sqr.errs),
"Std_Mean_Sqr_Error" = sd(mean.sqr.errs))
}Esta fucnión es tal vez ala mas importante como lo hemos visto
kfold.cval.reg.iter <- function(k, df, fold){
tr.ids <- !fold %in% c(k)# los que no estan el el fold
test.ids <- fold %in% c(k)
mod <- lm(MEDV ~., data = df[tr.ids,])#este es Lineal , pero puede ser cualquiera de los que se han visto
pred <- predict(mod, df[test.ids,])#predecimos los valores con los de entrenamiento
sqr.errs <- (pred - df[test.ids,"MEDV"])^2# calculamos los valores con los de entreanamiento
mean(sqr.errs)
}como se ocupa la funcion
res <- kfold.crossval.reg(bh, 5)## 2 5 5 4 3 1 5 1 4 4 1 4 4 4 1 4 1 5 3 4 5 3 3 1 4 2 4 5 3 2 1 2 5 5 2 3 1 5 4 5 1 2 5 4 4 3 4 5 2 3 3 2 4 3 3 5 4 2 1 2 4 2 1 4 2 4 4 4 2 4 2 1 5 3 1 1 5 3 2 3 2 3 1 4 2 5 1 1 4 3 4 5 1 1 1 1 2 5 2 1 4 4 5 5 3 5 2 3 4 5 2 3 1 2 2 5 3 5 4 5 1 5 1 5 3 4 4 5 1 1 5 5 1 3 2 5 3 3 5 3 5 3 2 3 1 3 2 4 2 5 2 3 3 3 1 4 3 3 1 2 1 3 2 2 3 1 5 4 3 5 3 4 4 2 4 5 1 3 4 5 2 4 3 2 4 1 1 5 3 1 1 2 5 3 5 1 2 1 2 1 2 4 1 1 1 5 5 5 1 1 3 3 3 1 2 5 4 2 2 2 4 2 1 2 3 4 4 5 1 3 5 4 3 5 4 5 3 5 1 4 4 4 1 5 2 3 3 3 5 4 5 2 2 3 2 5 4 5 3 4 1 2 3 4 4 3 4 2 2 5 4 5 5 1 4 4 1 2 2 5 5 1 4 2 4 3 5 4 1 5 3 3 4 1 3 1 3 5 4 5 2 5 5 2 2 5 3 3 1 5 2 4 4 5 5 3 3 4 4 2 4 5 3 3 5 4 2 3 5 5 5 5 2 5 1 3 1 1 2 2 5 5 1 3 1 5 4 4 2 1 3 4 5 1 1 3 5 5 2 2 5 2 1 3 2 3 5 4 5 1 1 2 5 2 2 5 4 2 2 3 4 3 2 4 3 1 4 4 2 4 3 1 1 4 3 4 4 2 3 3 4 1 1 3 3 3 4 5 4 2 3 2 5 2 5 2 5 5 4 1 5 4 2 1 2 3 5 4 2 2 4 1 1 1 4 4 1 4 5 4 3 5 2 4 4 1 4 3 4 1 1 3 5 5 1 2 1 1 3 1 2 1 4 2 2 5 5 3 4 4 3 5 2 1 4 4 1 5 2 2 4 4 5 2 3 4 2 5 4 2 2 3 4 4 2 1 4 5 4 3 5 5 5 3 2 2res## $`MSE `
## [1] 26.53045 29.62104 19.72100 15.59297 29.42232
##
## $Overall_Mean_Sqr_Error
## [1] 24.17756
##
## $Std_Mean_Sqr_Error
## [1] 6.249624Función que
loocv.reg <- function(df){
mean.sqr.errors <- sapply(1:nrow(df),
loocv.reg.iter, df)
list("MSE"=mean.sqr.errors,
"overall_mean_sqr_errors" = mean(mean.sqr.errors),
"sd_mean_sqr_errors" = sd(mean.sqr.errors))
}#Leave one out cross validatión La técnica que hace un muestreo de los datos dejando siempre un elemento como testing de los mismos usando el resto como conjunto de entrenamiento, de modo que en algún momento todo elemento del data set habrá sido probado en el mismo es
loocv.reg.iter <- function(k, df){
mod <- lm(MEDV~., data = df[-k,])
pred <- predict(mod, df[k,])
sqr.error <- (pred - df[k,"MEDV"])^2
sqr.error
}Vemos el resultado de la funcion locv.reg(bh)
res <- loocv.reg(bh)
res## $MSE
## 1 2 3 4 5 6
## 3.729797e+01 1.199961e+01 1.748783e+01 2.371101e+01 7.041998e+01 1.221852e+01
## 7 8 9 10 11 12
## 1.066078e-02 6.083233e+01 2.760620e+01 4.354925e-04 1.718911e+01 7.572428e+00
## 13 14 15 16 17 18
## 6.522578e-01 7.377726e-01 1.208910e+00 3.733962e-01 6.923949e+00 3.547527e-01
## 19 20 21 22 23 24
## 1.697014e+01 4.355471e-02 1.201640e+00 3.821703e+00 4.135365e-01 4.993998e-01
## 25 26 27 28 29 30
## 6.331468e-03 2.721472e-01 1.331116e+00 8.652999e-03 1.362318e+00 1.574726e-02
## 31 32 33 34 35 36
## 1.611519e+00 1.313453e+01 2.067570e+01 1.443976e+00 4.475714e-02 2.463602e+01
## 37 38 39 40 41 42
## 5.592485e+00 4.552686e+00 3.286605e+00 3.312408e-01 5.014142e-01 2.130980e+00
## 43 44 45 46 47 48
## 9.679261e-03 8.527549e-03 3.100071e+00 8.032073e+00 1.844215e-01 2.157248e+00
## 49 50 51 52 53 54
## 3.108574e+01 4.973245e+00 2.566659e+00 1.243208e+01 7.237428e+00 4.330235e-01
## 55 56 57 58 59 60
## 1.389389e+01 1.940582e+01 2.618686e-02 2.472834e+00 2.411174e+00 2.266668e+00
## 61 62 63 64 65 66
## 7.126901e-01 6.707352e+00 3.332020e+00 6.237148e+00 1.001863e+02 4.957845e+01
## 67 68 69 70 71 72
## 3.986293e+01 8.085272e-01 4.762654e-04 1.357698e-02 1.056106e+00 1.892730e-03
## 73 74 75 76 77 78
## 3.244263e+00 4.352050e-01 2.105917e+00 6.779739e+00 8.939215e+00 6.729250e+00
## 79 80 81 82 83 84
## 3.953892e-03 4.682032e+00 1.670276e-01 9.770162e+00 1.549403e+00 4.721297e+00
## 85 86 87 88 89 90
## 7.947959e-01 1.439888e+00 1.121866e-01 1.391472e+01 5.172308e+01 4.657649e+00
## 91 92 93 94 95 96
## 2.092745e+01 3.007600e+01 3.806832e+01 1.761120e+01 4.348412e+01 5.188964e-02
## 97 98 99 100 101 102
## 1.140565e+01 8.993158e+00 7.964098e+01 9.323476e-01 8.796581e+00 8.471550e-01
## 103 104 105 106 107 108
## 1.561794e+00 1.054589e+00 1.843013e+00 9.606184e-01 5.567635e+00 1.272343e-01
## 109 110 111 112 113 114
## 8.422781e+00 1.431191e-01 1.141503e+00 1.421192e+01 4.012310e+00 4.185680e+00
## 115 116 117 118 119 120
## 4.570305e+01 4.663181e+00 4.841275e+00 2.068423e+01 4.232158e-03 2.282826e+00
## 121 122 123 124 125 126
## 7.941337e-03 5.336048e+00 3.731817e-03 1.002866e+00 3.516524e+00 1.325469e+00
## 127 128 129 130 131 132
## 1.355923e+00 1.089893e+00 9.276111e-01 6.229524e-02 7.355950e-01 3.796536e-02
## 133 134 135 136 137 138
## 9.110577e+00 7.297429e+00 5.771485e+00 7.366954e-01 2.418407e+00 5.385606e+00
## 139 140 141 142 143 144
## 2.782415e-01 1.916950e+00 1.958247e-01 1.187180e+02 1.648491e+00 1.312498e+01
## 145 146 147 148 149 150
## 1.051908e+01 3.492333e+00 5.403145e-02 4.135578e+01 7.243438e+01 3.865230e-01
## 151 152 153 154 155 156
## 4.774235e-01 1.861299e+00 2.725144e+01 4.887157e+00 3.267753e+01 2.424213e+01
## 157 158 159 160 161 162
## 3.113948e-01 6.897296e+01 2.403866e+01 5.712380e+00 3.625518e+01 1.883310e+02
## 163 164 165 166 167 168
## 1.003511e+02 7.590807e+01 4.601966e+00 1.531457e-01 1.776428e+02 5.385503e-01
## 169 170 171 172 173 174
## 7.141359e+00 1.992025e+01 2.800766e+01 2.847223e+01 1.597189e-02 3.089412e+01
## 175 176 177 178 179 180
## 1.593424e+01 1.814710e+00 6.004608e+00 2.117228e+01 2.417781e+00 1.889265e+01
## 181 182 183 184 185 186
## 2.700943e+01 7.380785e+01 1.683769e+01 2.388951e+00 1.426146e+01 2.407323e+01
## 187 188 189 190 191 192
## 2.082895e+02 2.136408e+00 7.143840e+00 1.544500e-01 4.038493e+01 4.588249e-02
## 193 194 195 196 197 198
## 1.256855e+01 1.061280e+00 6.230196e+00 8.896089e+01 8.516689e+00 6.006171e+00
## 199 200 201 202 203 204
## 1.170881e-02 2.496221e+01 5.508817e+00 2.875440e+01 2.919602e+01 4.615213e+01
## 205 206 207 208 209 210
## 5.129553e+01 8.475094e-03 5.215106e-01 2.201520e+01 8.825574e-01 9.969343e+00
## 211 212 213 214 215 216
## 5.206310e-01 5.555872e+00 1.244453e-01 8.517094e+00 1.843480e+02 2.440186e-01
## 217 218 219 220 221 222
## 1.180092e+01 1.224423e-01 1.264901e+01 4.821454e+01 4.536309e+01 4.762079e+00
## 223 224 225 226 227 228
## 2.311447e+01 1.282296e-01 4.341674e+01 1.114443e+02 2.034105e-04 6.545355e-01
## 229 230 231 232 233 234
## 1.348988e+02 7.382200e-02 3.471605e-02 2.588687e+00 1.403737e+01 1.299645e+02
## 235 236 237 238 239 240
## 7.888862e+00 1.634536e+00 2.682745e+01 1.525101e+00 2.281558e+01 2.668342e+01
## 241 242 243 244 245 246
## 2.865031e+01 1.363589e+01 3.776324e+00 1.412930e+01 1.715674e+00 2.743474e+01
## 247 248 249 250 251 252
## 1.891622e+01 4.323744e-01 1.062594e+01 4.655989e+00 3.867855e-02 6.070924e-02
## 253 254 255 256 257 258
## 2.333755e+01 1.857514e+02 4.552659e+00 6.702129e-01 4.504348e+01 4.929602e+01
## 259 260 261 262 263 264
## 2.498378e-01 2.562539e+01 1.085240e+00 3.742820e+01 6.618498e+01 1.265035e+01
## 265 266 267 268 269 270
## 4.714732e-01 3.236633e+01 2.956476e-01 9.015552e+01 1.837542e+01 2.719936e+01
## 271 272 273 274 275 276
## 1.479861e+00 4.164324e+00 1.721232e+01 8.352803e-02 1.497046e+01 3.314929e+00
## 277 278 279 280 281 282
## 6.268429e+00 3.281041e+00 1.615043e+00 4.566495e-02 4.540504e+01 1.176022e+00
## 283 284 285 286 287 288
## 3.475458e+01 3.231028e+01 3.856496e-01 2.956488e+01 3.221712e-06 1.517697e+01
## 289 290 291 292 293 294
## 2.489079e+01 4.593853e+00 2.655275e+01 9.144115e+00 1.689299e+01 3.928269e+00
## 295 296 297 298 299 300
## 7.862971e+00 2.158638e-02 7.318452e-02 6.199754e-01 4.571617e+01 8.848605e+00
## 301 302 303 304 305 306
## 3.740617e+01 4.908793e+01 6.331184e+00 9.313059e-02 8.707894e+00 5.865127e+00
## 307 308 309 310 311 312
## 4.925755e+00 2.134519e+01 3.490444e+01 1.097246e+01 6.295753e+00 2.339483e+01
## 313 314 315 316 317 318
## 1.532674e+01 1.582541e+01 2.870317e+00 1.912256e+01 3.450631e-02 2.052398e+00
## 319 320 321 322 323 324
## 1.430690e+00 1.069377e-01 1.198376e+00 3.176837e+00 6.198795e+00 9.227381e-01
## 325 326 327 328 329 330
## 1.416033e-02 4.804102e-03 4.765176e-01 8.336834e+00 3.781088e+00 2.954469e+00
## 331 332 333 334 335 336
## 3.441631e+00 8.528804e+00 1.587045e+01 3.698186e-03 7.654037e-01 2.415336e-01
## 337 338 339 340 341 342
## 4.512952e-01 6.365567e-01 2.542731e+00 5.219794e+00 7.668080e+00 5.806470e+00
## 343 344 345 346 347 348
## 3.408684e+01 1.489448e+01 7.243904e+00 9.687327e-01 6.213789e+00 4.995692e+00
## 349 350 351 352 353 354
## 9.689524e+00 2.109834e+01 6.324526e+00 1.385283e+01 3.218312e+00 2.582157e+01
## 355 356 357 358 359 360
## 1.685384e+01 1.801083e+01 3.683016e+00 1.141446e+00 2.744186e-01 1.208874e+01
## 361 362 363 364 365 366
## 5.772775e+00 9.752284e-01 7.091829e+00 1.291493e+01 2.852837e+02 2.149911e+02
## 367 368 369 370 371 372
## 4.275165e+01 1.779938e+02 7.873947e+02 3.351657e+02 2.640801e+02 6.543936e+02
## 373 374 375 376 377 378
## 6.339490e+02 6.339202e+01 1.869516e+02 1.116278e+02 1.532786e+01 5.002861e+01
## 379 380 381 382 383 384
## 7.834563e+00 4.540748e+01 3.271900e+01 5.936606e+01 4.666056e+00 5.998654e-01
## 385 386 387 388 389 390
## 3.213960e+01 7.774770e-01 2.024656e+01 3.378258e+00 1.470464e+01 7.556259e+00
## 391 392 393 394 395 396
## 4.588322e+00 3.579287e+01 3.790106e-02 4.228283e+01 2.818196e+01 5.328189e+01
## 397 398 399 400 401 402
## 4.761296e+01 6.304668e+01 2.668367e+00 2.232802e+01 4.152712e+01 1.158000e+02
## 403 404 405 406 407 408
## 3.896707e+01 2.289148e+01 1.396653e+00 1.449052e+01 1.583904e+01 6.525032e+01
## 409 410 411 412 413 414
## 1.275883e+01 6.108433e+01 6.495700e-02 6.170145e-02 2.925406e+02 2.156831e+01
## 415 416 417 418 419 420
## 1.483805e+02 6.146177e+00 3.742888e+01 1.289712e+01 1.072386e+01 4.155393e+01
## 421 422 423 424 425 426
## 8.599456e+00 1.575544e+01 5.345426e+00 5.494603e-02 9.396361e+00 2.815285e+00
## 427 428 429 430 431 432
## 4.101458e+01 1.140076e+01 1.113683e+01 1.329283e+01 1.418733e+01 2.253565e+01
## 433 434 435 436 437 438
## 3.126912e+01 7.832180e+00 1.898620e+01 1.566939e-03 2.581520e+01 1.540751e-02
## 439 440 441 442 443 444
## 1.363672e+01 7.302687e-02 5.062740e+00 3.949323e-02 1.197830e-01 7.300796e+00
## 445 446 447 448 449 450
## 5.229042e-01 3.231554e-02 7.953288e+00 3.147774e+01 1.200647e+01 1.829087e+01
## 451 452 453 454 455 456
## 1.057631e+01 1.831199e+01 6.343967e+00 2.318438e+01 1.557890e-01 3.159516e+00
## 457 458 459 460 461 462
## 1.753985e-04 4.202941e-01 5.355962e+00 2.282126e+00 7.235331e+00 6.317140e+00
## 463 464 465 466 467 468
## 7.758864e-02 5.171410e+00 1.204679e+00 4.230099e+00 2.308054e+01 4.803341e+00
## 469 470 471 472 473 474
## 4.610735e+00 2.372174e+00 7.392626e-02 1.200459e+01 5.737959e-01 1.851058e+01
## 475 476 477 478 479 480
## 6.970657e+00 8.239081e+00 1.522962e+01 2.160582e-01 2.176554e+01 2.202696e-01
## 481 482 483 484 485 486
## 2.296869e-01 1.214821e+01 1.347174e+01 5.440092e-01 1.376131e+00 1.091587e+00
## 487 488 489 490 491 492
## 3.195204e-01 5.479531e-01 1.304814e+01 1.753278e+00 2.335405e+01 2.967204e-02
## 493 494 495 496 497 498
## 2.040458e+01 1.419016e+00 1.572809e+01 4.111324e+01 3.327803e+01 6.650704e-01
## 499 500 501 502 503 504
## 9.817067e-03 9.329858e-01 1.369270e+01 1.336045e+00 3.277817e+00 1.463799e+01
## 505 506
## 1.785779e+01 1.136972e+02
##
## $overall_mean_sqr_errors
## [1] 23.72575
##
## $sd_mean_sqr_errors
## [1] 65.33334este conjunto de funciones hace la cantidad de daros que se tiene , pues ese numero , hace regresiones lineales para cada uno de los valores
Reduciendo datos con Clustering o analisis componantes principales
Merodos partitivos Aglomeramos o separamos ( divisitibo)( dendograma ) Metodos gerargicos #clustering gerargico
protein <- read.csv("../../r-course-master/data/tema5/protein.csv")Escalamos los valores numericos
data <- as.data.frame(scale(protein[,-1]))pero se elimino la columna , no pasa nada, lo incluimos
data$Country = protein$CountryLe ponemos numero a las filas
rownames(data) = data$Countrycluster gerargico ( este es “el mas eficiente”)
hc <- hclust(dist# medir distancias
(data,#Data frame
method = "euclidean"),#metodo para distancia
method = "ward.D2")#es de minima varianza ## Warning in dist(data, method = "euclidean"): NAs introducidos por coerciónhc##
## Call:
## hclust(d = dist(data, method = "euclidean"), method = "ward.D2")
##
## Cluster method : ward.D2
## Distance : euclidean
## Number of objects: 25veamoslo graficamente
plot(hc, hang = -0.01, cex = 0.7)
hc2 <- hclust(dist(data, method = "euclidean"),
method = "single")## Warning in dist(data, method = "euclidean"): NAs introducidos por coerciónplot(hc2, hang=-0.01, cex = 0.7)
Las diferencias de metodos de aglomeracion
hc3 <- hclust(dist(data, method = "euclidean"),
method = "complete")## Warning in dist(data, method = "euclidean"): NAs introducidos por coerciónplot(hc3, hang=-0.01, cex = 0.7)
hc3$merge## [,1] [,2]
## [1,] -18 -25
## [2,] -2 -14
## [3,] -6 -20
## [4,] -3 -24
## [5,] -12 4
## [6,] -5 -7
## [7,] -4 1
## [8,] -15 3
## [9,] -10 -13
## [10,] -21 2
## [11,] -9 -22
## [12,] -8 8
## [13,] -16 6
## [14,] 5 10
## [15,] -17 -19
## [16,] -1 7
## [17,] -23 13
## [18,] -11 17
## [19,] 11 14
## [20,] 9 16
## [21,] 12 19
## [22,] 18 21
## [23,] 15 20
## [24,] 22 23hc4 <- hclust(dist(data, method = "euclidean"),
method = "average")## Warning in dist(data, method = "euclidean"): NAs introducidos por coerciónplot(hc4, hang=-0.01, cex = 0.7)
Esta matriz lo que hace es te da las posiciones de como aglomera cada uno de los valores
hc4$merge## [,1] [,2]
## [1,] -18 -25
## [2,] -2 -14
## [3,] -6 -20
## [4,] -3 -24
## [5,] -12 4
## [6,] -15 3
## [7,] -5 -7
## [8,] -4 1
## [9,] -21 2
## [10,] -10 -13
## [11,] 5 9
## [12,] -8 6
## [13,] -16 7
## [14,] -9 -22
## [15,] 11 13
## [16,] -1 8
## [17,] -17 -19
## [18,] 14 15
## [19,] 12 18
## [20,] -11 16
## [21,] -23 20
## [22,] 10 21
## [23,] 17 22
## [24,] 19 23En esta tabla podemos ver como se ven las relaciones numericas entre los nuemros de aglomeracion
d <- dist(data, method = "euclidean")## Warning in dist(data, method = "euclidean"): NAs introducidos por coerciónd ## Albania Austria Belgium Bulgaria Czechoslovakia Denmark
## Austria 6.467966
## Belgium 6.270545 2.583921
## Bulgaria 2.914078 5.165346 5.520376
## Czechoslovakia 5.419246 2.246156 2.340354 4.164048
## Denmark 6.993021 3.182556 2.676462 6.363126 3.551469
## E Germany 6.738023 2.721647 2.228452 5.701727 1.981974 2.914117
## Finland 6.197500 4.291913 3.693850 6.132320 4.199483 2.775939
## France 6.648173 3.783072 2.313117 5.861471 3.551333 3.860951
## Greece 4.486144 5.443265 4.949126 3.966231 5.133578 5.899629
## Hungary 4.926246 3.463473 4.210760 3.526341 2.900155 5.311659
## Ireland 7.122631 2.889077 1.746375 6.548645 3.328471 2.982614
## Italy 4.243320 3.918642 3.919850 3.020464 3.527166 5.037933
## Netherlands 6.333064 1.184075 2.370619 5.448173 2.312485 2.673810
## Norway 5.760859 4.085139 3.120846 5.562267 3.717820 2.101509
## Poland 6.201022 2.947240 3.094715 4.674394 2.218277 4.052664
## Portugal 6.969670 6.882349 5.956965 6.329434 5.817515 6.187499
## Romania 2.835084 4.902056 5.017853 1.991640 3.754924 5.833256
## Spain 5.869554 5.152500 4.213951 5.103852 4.373624 5.388490
## Sweden 5.959496 3.095772 2.725222 5.688855 3.433509 1.456405
## Switzerland 5.400868 2.321822 2.471070 4.725284 2.765509 3.360708
## UK 6.261715 3.950516 2.051293 6.109584 4.038158 3.662983
## USSR 4.580364 4.387765 3.331584 4.038101 2.863559 4.387004
## W Germany 6.698441 1.733344 1.494659 5.914431 2.302033 2.521632
## Yugoslavia 3.101508 5.740029 5.906917 2.100768 4.575402 6.706354
## E Germany Finland France Greece Hungary Ireland Italy
## Austria
## Belgium
## Bulgaria
## Czechoslovakia
## Denmark
## E Germany
## Finland 4.292653
## France 3.998525 4.843663
## Greece 5.923580 5.801368 4.790900
## Hungary 3.875142 5.686633 5.243799 4.332346
## Ireland 3.197670 3.411083 3.322193 6.013141 5.078293
## Italy 4.554132 5.180912 4.007816 2.266429 3.323988 5.105872
## Netherlands 2.668781 3.567546 3.592490 5.434954 3.679952 2.470828 4.132090
## Norway 3.449373 2.174253 4.132575 4.877926 5.173575 3.804993 4.216275
## Poland 2.855449 4.352077 3.793482 4.652913 3.207114 3.939571 3.286878
## Portugal 5.536763 6.859720 5.961973 5.042438 6.006109 7.445740 4.914201
## Romania 5.042962 5.340870 5.825068 3.815704 2.604885 5.907937 3.277615
## Spain 4.308394 5.777546 4.690822 3.266241 4.090129 5.568623 3.029450
## Sweden 3.225851 2.177177 4.023670 5.255041 4.917562 3.009951 4.359957
## Switzerland 3.764376 3.728661 2.555926 4.329069 4.062262 2.968257 3.095835
## UK 4.140596 4.092941 2.710334 4.871947 5.396215 2.375613 4.411907
## USSR 3.601734 3.657170 4.466362 4.336835 3.615445 4.109018 3.752108
## W Germany 1.994169 3.851187 3.094365 5.653969 4.113548 1.905283 4.361052
## Yugoslavia 5.823788 6.108572 6.647119 4.143284 3.194559 6.811222 3.774715
## Netherlands Norway Poland Portugal Romania Spain
## Austria
## Belgium
## Bulgaria
## Czechoslovakia
## Denmark
## E Germany
## Finland
## France
## Greece
## Hungary
## Ireland
## Italy
## Netherlands
## Norway 3.545298
## Poland 2.922854 3.907456
## Portugal 6.714237 5.055744 5.075801
## Romania 4.893325 4.936560 4.168287 5.934467
## Spain 5.129482 4.386674 3.582102 3.091414 4.472014
## Sweden 2.532559 1.585203 4.049537 6.186614 5.122213 5.069184
## Switzerland 2.003929 3.518374 3.239766 6.453545 4.594953 4.827802
## UK 3.707378 3.741913 4.742940 6.891621 5.716957 4.969019
## USSR 4.091736 3.436253 3.074917 5.349661 2.905607 3.823926
## W Germany 1.341817 3.477531 3.159113 6.474562 5.365969 4.852100
## Yugoslavia 5.800653 5.700872 4.733978 6.141142 1.039579 4.814148
## Sweden Switzerland UK USSR W Germany
## Austria
## Belgium
## Bulgaria
## Czechoslovakia
## Denmark
## E Germany
## Finland
## France
## Greece
## Hungary
## Ireland
## Italy
## Netherlands
## Norway
## Poland
## Portugal
## Romania
## Spain
## Sweden
## Switzerland 2.834669
## UK 3.302204 2.995056
## USSR 3.974365 4.000245 4.222185
## W Germany 2.602745 2.408679 3.051419 4.105817
## Yugoslavia 6.010580 5.491393 6.605472 3.536161 6.286440Calculemos las distanciasentre los paises
alb<-data["Albania",-10]
aus<-data["Austria",-10]
sqrt(sum((alb-aus)^2))#Errror medio ## [1] 6.136051sum(abs(alb-aus))#este se calcula para la distancia de manhatan## [1] 15.97134Como hacemos el clustering divicitivo
#install.packages("cluster")
library(cluster)## Warning: package 'cluster' was built under R version 4.1.3dv <- diana(data, metric = "euclidean")
par(mfrow=c(1,2))
plot(dv)
Ver como se comport la separacionde cada uno de esos valores
par(mfrow=c(1,1))
#Se recorta el arbol en el numero indicado
fit <- cutree(hc, k=4)
#Asi se ve como se van categorizando los datos
fit## Albania Austria Belgium Bulgaria Czechoslovakia
## 1 2 2 1 3
## Denmark E Germany Finland France Greece
## 2 3 2 2 4
## Hungary Ireland Italy Netherlands Norway
## 3 2 4 2 2
## Poland Portugal Romania Spain Sweden
## 3 4 1 4 2
## Switzerland UK USSR W Germany Yugoslavia
## 2 2 3 2 1#Esta tabla se ejecutan los datos se agruparon por grupo (parte superior) y como se fueron acomodando los paises(el nuemro de abajo )
table(fit)## fit
## 1 2 3 4
## 4 12 5 4#graficamos los datos
plot(hc, hang = -0.01, cex = 0.7)
#Le pedimos que nos pinte las segmentaciones
rect.hclust(hc, k=4, border="red")
# Clustering con K-means
Cuidado con los outliners por que se mueve sumamente feo
protein <- read.csv("../../r-course-master/data/tema5/protein.csv")
rownames(protein) = protein$Country # asignamos los valores como identificador para cada fila Eliminamos la columna donde se mostraban estos numeros
protein$Country = NULLAhora que todo los valores son numericos , entoces escalamos los valores
protein.scaled = as.data.frame(scale(protein))Cargamos e instalamos desde paquete desde git_hub Minimiza el uso a traves de analisis de componentes principales
library(devtools)
devtools::install_github("kassambara/factoextra")#esta libreria nos ayuda cuando tenemos muchas variables ## WARNING: Rtools is required to build R packages, but is not currently installed.
##
## Please download and install Rtools 4.0 from https://cran.r-project.org/bin/windows/Rtools/.## Skipping install of 'factoextra' from a github remote, the SHA1 (1689fc74) has not changed since last install.
## Use `force = TRUE` to force installation#Se puede validar en el github kassambra (https://github.com/kassambara/factoextra)Creamos los valores que nos ayudan a hacer los K-means devuelve los grupos a los cuales pertenece cada uno de esos valores Que porcentaje esta informacion esta resumida dentro de la informacion (between_SS( suma de cada cluster ) / total_SS(suma de los cuadrados totales ) = 59.8 % )
km <- kmeans(protein.scaled,4)
km## K-means clustering with 4 clusters of sizes 10, 8, 5, 2
##
## Cluster means:
## RedMeat WhiteMeat Eggs Milk Fish Cereals Starch
## 1 -0.1218974 0.6101653 0.3972740 0.5711137 0.2486383 -0.6175974 0.3573866
## 2 -0.6096362 -0.6553728 -0.8934192 -0.7300065 -0.6859595 1.2382474 -0.8956083
## 3 1.5990065 0.2988565 0.9341308 0.6091128 -0.1422470 -0.5948180 0.3451473
## 4 -0.9494848 -1.1764767 -0.7480204 -1.4583242 1.8562639 -0.3779572 0.9326321
## Nuts Fr.Veg
## 1 -0.8823165 -0.38028649
## 2 1.0401967 -0.06153324
## 3 -0.3484949 0.10200104
## 4 1.1220326 1.89256278
##
## Clustering vector:
## Albania Austria Belgium Bulgaria Czechoslovakia
## 2 1 3 2 1
## Denmark E Germany Finland France Greece
## 1 1 1 3 2
## Hungary Ireland Italy Netherlands Norway
## 2 3 2 1 1
## Poland Portugal Romania Spain Sweden
## 1 4 2 4 1
## Switzerland UK USSR W Germany Yugoslavia
## 3 3 2 1 2
##
## Within cluster sum of squares by cluster:
## [1] 37.502503 39.167085 12.069794 4.300578
## (between_SS / total_SS = 56.9 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault"Hora lo que hacemos es
aggregate(protein.scaled,# de este DF
by = list(cluster = km$cluster),# de la lista de cluster que esta contenida en el metodo de calculo anterior
mean)#dame la media para cada cluster ( conjunto de datos )## cluster RedMeat WhiteMeat Eggs Milk Fish Cereals
## 1 1 -0.1218974 0.6101653 0.3972740 0.5711137 0.2486383 -0.6175974
## 2 2 -0.6096362 -0.6553728 -0.8934192 -0.7300065 -0.6859595 1.2382474
## 3 3 1.5990065 0.2988565 0.9341308 0.6091128 -0.1422470 -0.5948180
## 4 4 -0.9494848 -1.1764767 -0.7480204 -1.4583242 1.8562639 -0.3779572
## Starch Nuts Fr.Veg
## 1 0.3573866 -0.8823165 -0.38028649
## 2 -0.8956083 1.0401967 -0.06153324
## 3 0.3451473 -0.3484949 0.10200104
## 4 0.9326321 1.1220326 1.89256278Visualicemos como se ve de manera grafica este
library(factoextra)## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBafviz_cluster(km, data = protein.scaled)
#mini Bach K-means para segmentar imagenes intalamos paquetes y cargamos
#install.packages(c("OpenImageR", "ClusterR"))
library(OpenImageR)### Warning: package 'OpenImageR' was built under R version 4.1.3library(ClusterR)### Warning: package 'ClusterR' was built under R version 4.1.3## Loading required package: gtools##
## Attaching package: 'gtools'## The following object is masked from 'package:e1071':
##
## permutationsCargamos la imagen pequena
imagename <- "../../r-course-master/data/tema5/bodegon.jpg"Hacemos que leea la imagen
img <- readImage(imagename)Usamos un algoritmo que redimenciona el tamano y cantidad de pixeles dependiendo del algoritmo que se desee usar
img.resize <- resizeImage(img, 350, 350,
method = "bilinear")#lo que hacemos es redimencionar los pixeles que se desea
imageShow(img.resize)
Vectorizamos la imagen
img.vector <- apply(img.resize, 3,# se compacta de 3 en 3 los pixeles
as.vector)Veamos la dimencion del vector
dim(img.vector)#Elementos y pixeles ## [1] 122500 3Predice el centroide de cadauno de los trosos el los que se dividira la imagen
kmmb<-MiniBatchKmeans(img.vector,
clusters = 10,#aprox numero de colores
batch_size = 20, #
num_init = 5,#Numero inicial
max_iters = 100,#Maximo de iteraciones
init_fraction = 0.2,#
initializer = "kmeans++",#metodo
early_stop_iter = 10, # que no se pare nunca antes de la iteracion 10
verbose = F)
kmmb## KMeans Cluster
## Call: NULL
## Data cols: 3
## Centroids: 10
## BSS/SS:
## SS: = 15136.32 (WSS) + (BSS)Predice los centroides Mira cuales de los 10 colores tiene mas cerca
prmb <- predict_MBatchKMeans(img.vector,
kmmb$centroids)#de esta variable toma los centroides por alguna extrana rason esto no me lo compila
library(OpenImageR)
get.cent.mb <- kmmb$centroids#jala los centroides
new.img <- get.cent.mb[prmb,]#crea la nueva imagen
#dim(new.img) <- c(nrow(new.img), ncol(img),3)
#imageShow(new.img)#caracteristicas de clustering estandar 1. Tiende a crear clusters convexos 2.Tiende a crear clusters redondos 3. clasifica todos los puntos aunque sean outliners
#Tecnicas de Particionar alrededor de los k medoides medoide: es el que mas se parece a los datos de ese grupo de datos (geometricamente el de en medio )
protein<-read.csv("../../r-course-master/data/tema5/protein.csv")Trasformaciones al DF Nota importante: el escalado es fundamental esn estas tecnicas con la finalidad de evitar inconvenientes con los outliners
rownames(protein) = protein$Country#ponemos nombre a las filas
protein$Country = NULL#eliminamos esas filas
protein.scaled <- as.data.frame(scale(protein))#escalamos el DFCargamos las librerias
library(cluster)
library(factoextra)Medoides (algoritmo pam ) Caracteristica y ventaja de este algoritmo ( lo que hace es toma un valor existente dentro de Df y lo toma para cada uno)
km <- pam(protein.scaled,#DF
4)# numero de medoides
km## Medoids:
## ID RedMeat WhiteMeat Eggs Milk Fish Cereals
## Romania 18 -1.0839304 -0.43204252 -1.2848772 -0.8461152 -0.9651632 1.5810786
## W Germany 24 0.4696633 1.24631815 1.0415022 0.2375653 -0.2598064 -1.2435778
## Sweden 20 0.0215113 -0.02598752 0.5046454 1.0679178 0.9451781 -1.1615716
## Spain 19 -0.8150392 -1.21708219 0.1467409 -1.1979595 0.7982288 -0.2777275
## Starch Nuts Fr.Veg
## Romania -0.7196689 1.1220326 -0.7406162
## W Germany 0.5654541 -0.7916675 -0.1862628
## Sweden -0.3524909 -0.8420280 -1.1840990
## Spain 0.8714358 1.4241958 1.6985391
## Clustering vector:
## Albania Austria Belgium Bulgaria Czechoslovakia
## 1 2 2 1 2
## Denmark E Germany Finland France Greece
## 3 2 3 2 4
## Hungary Ireland Italy Netherlands Norway
## 1 2 4 2 3
## Poland Portugal Romania Spain Sweden
## 2 4 1 4 3
## Switzerland UK USSR W Germany Yugoslavia
## 2 2 1 2 1
## Objective function:
## build swap
## 1.847172 1.839215
##
## Available components:
## [1] "medoids" "id.med" "clustering" "objective" "isolation"
## [6] "clusinfo" "silinfo" "diss" "call" "data"Veamoslo graficamente
fviz_cluster(km)
Medoide ( algoritmo clara) Extrae muestras y utiliza el partition arround mediods(pAm) Tecnica de particion para datos enormes ( ayuda para problemas de memoria)
1._Hace muetreo del data set original y elabora un cluster a la muestra 2._Devuelve el mejor clustering de todos los que elabora
clarafit <- clara(protein.scaled, 4,
samples = 6)#tamano de la muestra
clarafit## Call: clara(x = protein.scaled, k = 4, samples = 6)
## Medoids:
## RedMeat WhiteMeat Eggs Milk Fish Cereals
## Romania -1.0839304 -0.43204252 -1.2848772 -0.8461152 -0.9651632 1.5810786
## W Germany 0.4696633 1.24631815 1.0415022 0.2375653 -0.2598064 -1.2435778
## Sweden 0.0215113 -0.02598752 0.5046454 1.0679178 0.9451781 -1.1615716
## Spain -0.8150392 -1.21708219 0.1467409 -1.1979595 0.7982288 -0.2777275
## Starch Nuts Fr.Veg
## Romania -0.7196689 1.1220326 -0.7406162
## W Germany 0.5654541 -0.7916675 -0.1862628
## Sweden -0.3524909 -0.8420280 -1.1840990
## Spain 0.8714358 1.4241958 1.6985391
## Objective function: 1.839215
## Clustering vector: Named int [1:25] 1 2 2 1 2 3 2 3 2 4 1 2 4 2 3 2 4 1 ...
## - attr(*, "names")= chr [1:25] "Albania" "Austria" "Belgium" "Bulgaria" "Czechoslovakia" "Denmark" "E Germany" ...
## Cluster sizes: 6 11 4 4
## Best sample:
## [1] Albania Austria Belgium Bulgaria Czechoslovakia
## [6] Denmark E Germany Finland France Greece
## [11] Hungary Ireland Italy Netherlands Norway
## [16] Poland Portugal Romania Spain Sweden
## [21] Switzerland UK USSR W Germany Yugoslavia
##
## Available components:
## [1] "sample" "medoids" "i.med" "clustering" "objective"
## [6] "clusinfo" "diss" "call" "silinfo" "data"Veamos este algoritmo
fviz_cluster(clarafit)
#Validando los resultados de un cluster
Como puedo saber que porcenteje estaq bien catalogado y cuando no esto sirve para elejir o variar el metodo de clustering Validacion de estabilidad (consistencia del cluster )= va quitando columnas (variables )y ve si se meueve mucho la casificacion Cargamos los paquetes
#install.packages(c("fpc", "NbClust"))
library(factoextra)
library(cluster)
library(fpc)## Warning: package 'fpc' was built under R version 4.1.3library(NbClust)volvemos a trabajar con el dataset de proteinas explicado en las lineas de codigo 2502 aprox
protein <- read.csv("../../r-course-master/data/tema5/protein.csv")
rownames(protein) = protein$Country
protein$Country = NULL
protein.scaled <- as.data.frame(scale(protein))Determinar numero de clusters que se desean nota ver “Conclusion”, muestra k=3 , pero el numero de diviciones es 7 que correponde a esa linea
nb <- NbClust(protein.scaled,#DF escalado
distance = "euclidean",#Tipo de distancia
min.nc = 2,# minimo de cluster
max.nc = 12, #Maximo de cluster
method = "ward.D2",# metodo aglomerativo
index = "all")#paramerto por defecto ## Warning in pf(beale, pp, df2): NaNs produced
## *** : The Hubert index is a graphical method of determining the number of clusters.
## In the plot of Hubert index, we seek a significant knee that corresponds to a
## significant increase of the value of the measure i.e the significant peak in Hubert
## index second differences plot.
## 
## *** : The D index is a graphical method of determining the number of clusters.
## In the plot of D index, we seek a significant knee (the significant peak in Dindex
## second differences plot) that corresponds to a significant increase of the value of
## the measure.
##
## *******************************************************************
## * Among all indices:
## * 6 proposed 2 as the best number of clusters
## * 7 proposed 3 as the best number of clusters
## * 3 proposed 6 as the best number of clusters
## * 1 proposed 7 as the best number of clusters
## * 1 proposed 10 as the best number of clusters
## * 5 proposed 12 as the best number of clusters
##
## ***** Conclusion *****
##
## * According to the majority rule, the best number of clusters is 3
##
##
## *******************************************************************generemos una visualizacion mas sencilla El analisis
fviz_nbclust(nb) + theme_minimal()## Warning in if (class(best_nc) == "numeric") print(best_nc) else if
## (class(best_nc) == : la condición tiene longitud > 1 y sólo el primer elemento
## será usado## Warning in if (class(best_nc) == "matrix") .viz_NbClust(x, print.summary, : la
## condición tiene longitud > 1 y sólo el primer elemento será usado## Warning in if (class(best_nc) == "numeric") print(best_nc) else if
## (class(best_nc) == : la condición tiene longitud > 1 y sólo el primer elemento
## será usado## Warning in if (class(best_nc) == "matrix") {: la condición tiene longitud > 1 y
## sólo el primer elemento será usado## Among all indices:
## ===================
## * 2 proposed 0 as the best number of clusters
## * 1 proposed 1 as the best number of clusters
## * 6 proposed 2 as the best number of clusters
## * 7 proposed 3 as the best number of clusters
## * 3 proposed 6 as the best number of clusters
## * 1 proposed 7 as the best number of clusters
## * 1 proposed 10 as the best number of clusters
## * 5 proposed 12 as the best number of clusters
##
## Conclusion
## =========================
## * According to the majority rule, the best number of clusters is 3 .
Ahora que tenemos el numero de clusters mas eficiente generamos kmeans
km.res <- kmeans(protein.scaled, 3)Creamos una sirueta : que es otro tipo de analisis
sil.km <- silhouette(km.res$cluster, #del objeto creado
dist(protein.scaled))#calcule la distancia del Df inicial Del objeto que creamos le la varibale de la sirueta le hacemos un summary
sil.sum <- summary(sil.km)
sil.sum## Silhouette of 25 units in 3 clusters from silhouette.default(x = km.res$cluster, dist = dist(protein.scaled)) :
## Cluster sizes and average silhouette widths:
## 4 10 11
## 0.4187964 0.1936819 0.2068237
## Individual silhouette widths:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.03801 0.14773 0.23808 0.23548 0.32609 0.46940visualicemos este objeto (la sirueta) Para medir la eficiencia del objeto anchura (analisis de DAnn),conectividad , separacion ,
fviz_silhouette(sil.km)## cluster size ave.sil.width
## 1 1 4 0.42
## 2 2 10 0.19
## 3 3 11 0.21
25 paises , pues 25 columnas mirar la grafica de arriba entre mas cerca de -1 esta peor clusterizado y entre mas cercano a 1 mejor esta clasificado
la linea intermedia es la linea roja punteada es el promedio de que tan bien se clasifico entre la k elejida
“ver teoria de clusters” nota : es importante notar que si todos salen positivos es muy bueno
fviz_cluster(km.res, data = protein.scaled)
dd <- dist(protein.scaled, method = "euclidean")#calculo de distancias
km_stats <- fpc::cluster.stats(dd, km.res$cluster)#calculo de los estados con base a las distancias recien hechas
km_stats## $n
## [1] 25
##
## $cluster.number
## [1] 3
##
## $cluster.size
## [1] 4 10 11
##
## $min.cluster.size
## [1] 4
##
## $noisen
## [1] 0
##
## $diameter
## [1] 2.633487 6.612010 4.499548
##
## $average.distance
## [1] 1.940134 3.792570 2.654427
##
## $median.distance
## [1] 2.028172 3.627695 2.581981
##
## $separation
## [1] 2.39223 2.71661 2.39223
##
## $average.toother
## [1] 4.196834 4.904011 4.398396
##
## $separation.matrix
## [,1] [,2] [,3]
## [1,] 0.000000 3.259916 2.39223
## [2,] 3.259916 0.000000 2.71661
## [3,] 2.392230 2.716610 0.00000
##
## $ave.between.matrix
## [,1] [,2] [,3]
## [1,] 0.000000 5.134783 3.344152
## [2,] 5.134783 0.000000 4.820093
## [3,] 3.344152 4.820093 0.000000
##
## $average.between
## [1] 4.550228
##
## $average.within
## [1] 2.995397
##
## $n.between
## [1] 194
##
## $n.within
## [1] 106
##
## $max.diameter
## [1] 6.61201
##
## $min.separation
## [1] 2.39223
##
## $within.cluster.ss
## [1] 114.7435
##
## $clus.avg.silwidths
## 1 2 3
## 0.4187964 0.1936819 0.2068237
##
## $avg.silwidth
## [1] 0.2354826
##
## $g2
## NULL
##
## $g3
## NULL
##
## $pearsongamma
## [1] 0.5320644
##
## $dunn
## [1] 0.3618008
##
## $dunn2
## [1] 0.8817642
##
## $entropy
## [1] 1.020961
##
## $wb.ratio
## [1] 0.658296
##
## $ch
## [1] 9.707054
##
## $cwidegap
## [1] 2.062677 2.932773 2.194416
##
## $widestgap
## [1] 2.932773
##
## $sindex
## [1] 2.39223
##
## $corrected.rand
## NULL
##
## $vi
## NULLSuma de los cuadrados internos del cluster
km_stats$within.cluster.ss# este es el total de todos los centros ## [1] 114.7435km_stats$clus.avg.silwidths#promedio de cada K## 1 2 3
## 0.4187964 0.1936819 0.2068237#indice de Dunn entre mas alejado este un cluster ( ente mas grande es este numero es mejor , significa que el equivocarse al clasificar se hace mas raro , minimiza el reisgo de equivocarse)
km_stats$dunn# indice de dun (distancia mas pequena de un dato entre dos clusters /distancia maxima entre elementos dentro de un cluster)## [1] 0.3618008Vamos a comparar la diferencia de los indices de Dunn entre metodos
kmed <- pam(protein.scaled, 3)
sil.kmed <- silhouette(kmed$clustering,
dist(protein.scaled))
fviz_cluster(kmed, data = protein.scaled)#graficamos las categorias del calculo de pam 
fviz_silhouette(sil.kmed)# calculamos la cirueta del paquete PAM ## cluster size ave.sil.width
## 1 1 6 0.31
## 2 2 15 0.38
## 3 3 4 0.21
kmed_stats <- cluster.stats(dd, kmed$clustering)
kmed_stats$within.cluster.ss## [1] 105.9863kmed_stats$clus.avg.silwidths## 1 2 3
## 0.3062056 0.3811136 0.2063242kmed_stats$dunn## [1] 0.5678917Comparamos los 2 clusters
res.com <- cluster.stats(dd, km.res$cluster,
kmed$clustering)
res.com$corrected.rand#Indice de correcion aleatoria ( simiitudes )entre -1 y 1( conhincidencia casi exacta )( este lo queremos mas hacia el 1)## [1] 0.5243425res.com$vi#Variacion entre custers (distancia entre custers )( lo mas pequeno posible )## [1] 0.6171538#tecnicas avanzadas de clustering Un [unto de la frontera en un cluster es quel que :tiene unvecindario inferior al minimo estipulado a la distancia menor o igual a la indicada #clustering basado en densidades ( cuando tenemos ruido y outliners )
library(fpc)
library(factoextra)
data("multishapes", package = "factoextra")# esta precargado como se ve el DF
par(mfrow=c(1,1))
head(multishapes)## x y shape
## 1 -0.8037393 -0.8530526 1
## 2 0.8528507 0.3676184 1
## 3 0.9271795 -0.2749024 1
## 4 -0.7526261 -0.5115652 1
## 5 0.7068462 0.8106792 1
## 6 1.0346985 0.3946550 1Obtenemos la columna 1(X) y 2(y)
dataPoints <- multishapes[,1:2]
plot(multishapes)
plot(dataPoints)# esta es la que vamos a trabajar 
intentemos meter estos valores dentro de un k-means Esta agrupando horrible (jajajaja)
km <- kmeans(dataPoints, 5)
fviz_cluster(km, data= dataPoints)
eliminamos los ruidos
dsFit <- dbscan(dataPoints,
eps = 0.15,#epsilon :Radio que debe de existir entre valores
MinPts = 5)#minimo de puntos:no se cree un cluster abajo del numero indicado
dsFit## dbscan Pts=1100 MinPts=5 eps=0.15
## 0 1 2 3 4 5
## border 31 24 1 5 7 1
## seed 0 386 404 99 92 50
## total 31 410 405 104 99 51veamos como lo categoriza
fviz_cluster(dsFit, dataPoints, geom = "point")
Clustering basado en modelos
BIC:Criterio de informacion bayesiano clasification : uncertanty :regiones de incertidumbre ( el algoritmo no sabe donde poner estos datos ) density:
#install.packages("mclust")
library(mclust)## Warning: package 'mclust' was built under R version 4.1.3## Package 'mclust' version 5.4.9
## Type 'citation("mclust")' for citing this R package in publications.##
## Attaching package: 'mclust'## The following object is masked _by_ '.GlobalEnv':
##
## banknote## The following object is masked from 'package:VIM':
##
## diabetes## The following object is masked from 'package:purrr':
##
## mapmclust <- Mclust(dataPoints)
plot(mclust)
que hizo Mclust Modelo Gaussiano mixto finito
summary(mclust)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust VEV (ellipsoidal, equal shape) model with 9 components:
##
## log-likelihood n df BIC ICL
## -2062.757 1100 45 -4440.653 -4577.743
##
## Clustering table:
## 1 2 3 4 5 6 7 8 9
## 134 187 112 132 128 142 108 104 53Reducir dimenciones con ACP(Analisis de componentes principales )
#install.packages("corrplot")
library(corrplot)#correlaciones ## Warning: package 'corrplot' was built under R version 4.1.3## corrplot 0.92 loadedbh <- read.csv("../../r-course-master/data/tema5/BostonHousing.csv")Hacemos el objeto para las correlaciones
corr <- cor(bh[,-14])
corr## CRIM ZN INDUS CHAS NOX
## CRIM 1.00000000 -0.20046922 0.40658341 -0.055891582 0.42097171
## ZN -0.20046922 1.00000000 -0.53382819 -0.042696719 -0.51660371
## INDUS 0.40658341 -0.53382819 1.00000000 0.062938027 0.76365145
## CHAS -0.05589158 -0.04269672 0.06293803 1.000000000 0.09120281
## NOX 0.42097171 -0.51660371 0.76365145 0.091202807 1.00000000
## RM -0.21924670 0.31199059 -0.39167585 0.091251225 -0.30218819
## AGE 0.35273425 -0.56953734 0.64477851 0.086517774 0.73147010
## DIS -0.37967009 0.66440822 -0.70802699 -0.099175780 -0.76923011
## RAD 0.62550515 -0.31194783 0.59512927 -0.007368241 0.61144056
## TAX 0.58276431 -0.31456332 0.72076018 -0.035586518 0.66802320
## PTRATIO 0.28994558 -0.39167855 0.38324756 -0.121515174 0.18893268
## B -0.38506394 0.17552032 -0.35697654 0.048788485 -0.38005064
## LSTAT 0.45562148 -0.41299457 0.60379972 -0.053929298 0.59087892
## RM AGE DIS RAD TAX PTRATIO
## CRIM -0.21924670 0.35273425 -0.37967009 0.625505145 0.58276431 0.2899456
## ZN 0.31199059 -0.56953734 0.66440822 -0.311947826 -0.31456332 -0.3916785
## INDUS -0.39167585 0.64477851 -0.70802699 0.595129275 0.72076018 0.3832476
## CHAS 0.09125123 0.08651777 -0.09917578 -0.007368241 -0.03558652 -0.1215152
## NOX -0.30218819 0.73147010 -0.76923011 0.611440563 0.66802320 0.1889327
## RM 1.00000000 -0.24026493 0.20524621 -0.209846668 -0.29204783 -0.3555015
## AGE -0.24026493 1.00000000 -0.74788054 0.456022452 0.50645559 0.2615150
## DIS 0.20524621 -0.74788054 1.00000000 -0.494587930 -0.53443158 -0.2324705
## RAD -0.20984667 0.45602245 -0.49458793 1.000000000 0.91022819 0.4647412
## TAX -0.29204783 0.50645559 -0.53443158 0.910228189 1.00000000 0.4608530
## PTRATIO -0.35550149 0.26151501 -0.23247054 0.464741179 0.46085304 1.0000000
## B 0.12806864 -0.27353398 0.29151167 -0.444412816 -0.44180801 -0.1773833
## LSTAT -0.61380827 0.60233853 -0.49699583 0.488676335 0.54399341 0.3740443
## B LSTAT
## CRIM -0.38506394 0.4556215
## ZN 0.17552032 -0.4129946
## INDUS -0.35697654 0.6037997
## CHAS 0.04878848 -0.0539293
## NOX -0.38005064 0.5908789
## RM 0.12806864 -0.6138083
## AGE -0.27353398 0.6023385
## DIS 0.29151167 -0.4969958
## RAD -0.44441282 0.4886763
## TAX -0.44180801 0.5439934
## PTRATIO -0.17738330 0.3740443
## B 1.00000000 -0.3660869
## LSTAT -0.36608690 1.0000000Grafiquemos las correlaciones Rojo =poco correlacionado azul =muy correlacionado
corrplot(corr, method = "circle")
#scale = T, matriz de correlaciones #scale = F, matriz de covarianzas
aqui vamos a correr el ejercicio donde podemos observar cuales son
bh.acp <- prcomp(bh[,-14], scale = T)#calculo de analisis de con cuantas variables que % del DF puedo explicar (cumulative Proportion )
summary(bh.acp)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 2.4752 1.1972 1.11473 0.92605 0.91368 0.81081 0.73168
## Proportion of Variance 0.4713 0.1103 0.09559 0.06597 0.06422 0.05057 0.04118
## Cumulative Proportion 0.4713 0.5816 0.67713 0.74310 0.80732 0.85789 0.89907
## PC8 PC9 PC10 PC11 PC12 PC13
## Standard deviation 0.62936 0.5263 0.46930 0.43129 0.41146 0.25201
## Proportion of Variance 0.03047 0.0213 0.01694 0.01431 0.01302 0.00489
## Cumulative Proportion 0.92954 0.9508 0.96778 0.98209 0.99511 1.00000Grafiquemoslo lo que vamos a observar son las componentes que explican el DF
plot(bh.acp)
plot(bh.acp, type = "lines")
Grafiquemos con color estas son las columnas que tenemos que utilizar
biplot(bh.acp, col = c("gray", "red"))
head(bh.acp$x,5)
La matriz de rotacion Cual es la que explica mejor los datos
#bh.acp$rotationDesviacion estandar conclusion ( las primeras 7 son suficientes para entender los datos ) alijerar ( disminuir el tamano del df)
bh.acp$sdev## [1] 2.4752472 1.1971947 1.1147272 0.9260535 0.9136826 0.8108065 0.7316803
## [8] 0.6293626 0.5262541 0.4692950 0.4312938 0.4114644 0.2520104#Series temporales y su utilidad Cargamos los datos
AAPL = read.csv("../../r-course-master/data/tema6/AAPL.csv", stringsAsFactors = F)
AMZN = read.csv("../../r-course-master/data/tema6/AMZN.csv", stringsAsFactors = F)
FB = read.csv("../../r-course-master/data/tema6/FB.csv", stringsAsFactors = F)
GOOG = read.csv("../../r-course-master/data/tema6/GOOG.csv", stringsAsFactors = F)
str(AAPL)## 'data.frame': 2519 obs. of 7 variables:
## $ Date : chr "2007-12-31" "2008-01-02" "2008-01-03" "2008-01-04" ...
## $ Open : num 28.5 28.5 27.9 27.4 25.9 ...
## $ High : num 28.6 28.6 28.2 27.6 26.2 ...
## $ Low : num 28.2 27.5 27.5 25.6 24.3 ...
## $ Close : num 28.3 27.8 27.8 25.7 25.4 ...
## $ Adj.Close: num 25.4 25 25 23.1 22.7 ...
## $ Volume : int 134833300 269794700 210516600 363958000 518048300 380954000 453470500 370743800 308071400 275112600 ...Filtramos el dia que no nos sirve
AAPL=AAPL[AAPL$Date>="2008-01-01",]
AMZN=AMZN[AMZN$Date>="2008-01-01",]
GOOG=GOOG[GOOG$Date>="2008-01-01",]Hacemos el casting a fecha de las mismas…
AAPL$Date = as.Date(AAPL$Date)
AMZN$Date = as.Date(AMZN$Date)
FB$Date = as.Date(FB$Date)
GOOG$Date = as.Date(GOOG$Date)Grafiquemos
library(ggplot2)
ggplot(AAPL, aes(Date#eje X
, Close#eje Y
)) +
geom_line(aes(color="Apple")) +
geom_line(data=AMZN, aes(color = "Amazon"))+
geom_line(data=FB, aes(color="Facebook"))+
geom_line(data=GOOG, aes(color="Google"))+
labs(color="Legend")+
scale_color_manual("",
breaks = c("Apple", "Amazon", "Facebook", "Google" ),
values = c("green", "yellow", "blue", "red"))+
ggtitle("Comparaciones de cierre de stocks")+
theme(plot.title = element_text(lineheight = 0.7, face = "bold"))
# cargar datos en tiempo real
#install.packages("quantmod")
library(quantmod)## Warning: package 'quantmod' was built under R version 4.1.3## Loading required package: xts## Warning: package 'xts' was built under R version 4.1.3## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric##
## Attaching package: 'xts'## The following objects are masked from 'package:dplyr':
##
## first, last## Loading required package: TTR## Warning: package 'TTR' was built under R version 4.1.3## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo##
## Attaching package: 'quantmod'## The following object is masked from 'package:Hmisc':
##
## Lagse conecta y es un objeto de tipo xts
getSymbols("AAPL")## 'getSymbols' currently uses auto.assign=TRUE by default, but will
## use auto.assign=FALSE in 0.5-0. You will still be able to use
## 'loadSymbols' to automatically load data. getOption("getSymbols.env")
## and getOption("getSymbols.auto.assign") will still be checked for
## alternate defaults.
##
## This message is shown once per session and may be disabled by setting
## options("getSymbols.warning4.0"=FALSE). See ?getSymbols for details.## [1] "AAPL"Grafica volumen en millones
barChart(AAPL)
chartSeries(AAPL, TA = "NULL")#tecnival indicatior 
candleChart(AAPL, up.col = "black",
dn.col="red", theme = "white")
Aqui muestra los porcentajes de cada uno de los valores ( variaciones en la media)
chartSeries(AAPL[,4],TA = "addMACD()")
Obtengamos la de Netflix
getSymbols("BTC")## [1] "BTC"barChart(BTC)
chartSeries(BTC, TA = "NULL")
chartSeries(BTC[,4],TA = "addMACD()")