ME.Junio.2018.Rmd

---
title: "Modelizacion Estad�stica. Segunda parte. Evaluaci�n extraordinaria. Junio 2018"
author: "F. Javier P�rez Ram�rez"
date: "2 de mayo de 2018"
output: 
  pdf_document:
    number_sections: true
    toc: true
    toc_depth: 3
header-includes:
  - \usepackage{xcolor}
  - \usepackage{framed}
  
  - \usepackage{fancyhdr}
  - \pagestyle{fancy}
  - \fancyhead[RO,RE]{Mod. Estad�stica. 2� parte.}
  - \fancyfoot{}
  - \fancyfoot[RE,LO]{F. Javier P�rez. 2017-2018}
  - \fancyfoot[LE,RO]{\thepage} 
bibliography: BibliografiaME.bib
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, warning=FALSE, message = FALSE)

#    include = FALSE prevents code and results from appearing in the finished file. R Markdown still runs the code in the chunk, and the results can be used by other chunks.
#    echo = FALSE prevents code, but not the results from appearing in the finished file. This is a useful way to embed figures.
#    message = FALSE prevents messages that are generated by code from appearing in the finished file.
#    warning = FALSE prevents warnings that are generated by code from appearing in the finished.
#    fig.cap = "..." adds a caption to graphical results.

```

***

\pagebreak

# Descripci�n del fichero de datos.
**Indicar los individuos analizados y describir las variables indicando el tipo (cualitativa o cuantitativa) y el resumen estad�stico m�s adecuado en cada caso.**

Se carga el fichero "coches.xls" y se comprueba que los datos se han le�do correctamente. A continuaci�n se muestra un resumen de los datos y su descripci�n.

```{r include=FALSE}
# Borrar toda las variables
### rm(list=ls())
```


```{r coches}
library(readxl)
coches <- read_excel("coches.xls")
# Se guarda una copia de los datos originales
coches.original <- coches
# head(coches)
```

Se tiene un fichero con datos relativos a veh�culos:

| Variable | Tipo | Descripci�n                                        |
|----------|------|----------------------------------------------------|
| `MARCA`  | Cualitativa | Marca del veh�culo |
| `MODELO` | Cualitativa | Modelo del veh�culo |
| `TIPO`   | Cualitativa | Tipo de veh�culo: `Automobile`" o "`Truck`" |
| `VENTAS`    | Cuantitativa | N�mero de unidades vendidas |
| `REVENTA`   | Cuantitativa | Precio de reventa del coche a los 4 a�os  |
| `PRECIO`    | Cuantitativa | Precio del veh�culo |
| `MOTOR`     | Cuantitativa | Tama�o del  motor (en c.c.) |
| `CV`        | Cuantitativa | Potencia del motor en CV |
| `PISADA`    | Cuantitativa | Base del neum�tico |
| `ANCHURA  ` | Cuantitativa | Ancho del veh�culo (en pulgadas) |
| `LONGITUD`  | Cuantitativa | Longitud del veh�culo (en pulgadas) |
| `PESO`      | Cuantitativa | Peso del veh�culo (en TM) |
| `CAPACIDAD` | Cuantitativa | Capacidad del tanque de combustible (ltr.) |
| `MPG`       | Cuantitativa | Consumo de combustible (Millas por gal�n) |


[comment_tipo_var]: <> (* Variables cualitativas: `MARCA`, `MODELO` y `TIPO`. La variable `TIPO` s�lo tiene dos posibles valores en este fichero: "`Automobile`" y "`Truck`"
* Variables cuantitativas: `VENTAS`, `REVENTA`, `PRECIO`, `MOTOR`, `CV`, `PISADA`, `ANCHURA`, `LONGITUD`, `PESO`, `CAPACIDAD` y `MPG`.)

[comment_color]: <> (\colorlet{shadecolor}{red!90} 
\begin{shaded}
Revisar lo de la variable booleana
\end{shaded})


## Inspecci�n de los valores de cada campo.

```{r}
par(mfrow=c(2,6))
boxplot(coches$VENTAS, xlab="VENTAS")
boxplot(coches$REVENTA, xlab="REVENTA")
boxplot(coches$PRECIO, xlab="PRECIO")
boxplot(coches$MOTOR, xlab="MOTOR")
boxplot(coches$CV, xlab="CV")
boxplot(coches$PISADA, xlab="PISADA")
boxplot(coches$ANCHURA, xlab="ANCHURA")
boxplot(coches$LONGITUD, xlab="LONGITUD")
boxplot(coches$PESO, xlab="PESO")
boxplot(coches$CAPACIDAD, xlab="CAPACIDAD")
boxplot(coches$MPG, xlab="MPG")
```

Observando el *boxplot* de cada una de las variables, parece que hubiera tres *outliers* en: 'VENTAS', 'MOTOR' Y 'MPG'. Se analizan a continuaci�n.

```{r include=FALSE}
# library(sandwich)
# library(carData)
# library(car)
library(RcmdrMisc)
numSummary(coches[,c("VENTAS", "REVENTA", "PRECIO", "MOTOR", "CV", "PISADA", "ANCHURA", "LONGITUD", "PESO", "CAPACIDAD", "MPG")], statistics=c("mean", "sd", "quantiles"), quantiles=c(0,.25,.5,.75, 1))
```


### Variable `VENTAS`.

```{r}
#  summary(coches$VENTAS)
sprintf("Desviaci�n estandar para la variable VENTAS: %.1f", sd(coches$VENTAS))
sprintf("Desviacion VENTAS m�ximas respecto a la media: %.1f", 
        (max(coches$VENTAS) - mean(coches$VENTAS)) / sd(coches$VENTAS))
sprintf("Desviacion ventas 'penultimo' elemento respecto a la media: %.1f", (276747 - mean(coches$VENTAS)) / sd(coches$VENTAS))
```
La media de veh�culos vendidos es de 58.619 y su error est�ndar es de 74.372. Sin embargo existe un veh�culo (el 'Ford' 'F-Series') cuyas ventas son 540.561 que va mucho m�s all� de lo esperable puesto que su desviaci�n est�ndar respecto a la media es de 6.5. M�s a�n la ventas de este modelo pr�cticamente duplican las del anterior que son de 276.747 ventas, cuya desviaci�n est�ndar respecto a la media es de 2.9; siendo este �ltimo un valor admisible.

Sin embargo se puede comprobar que s� se han vendido  esas unidades de ese modelo en Estados Unidos. Por tanto ese dato es correcto (Fuente: `http://carsalesbase.com/us-car-sales-data/ford/ford-f-series/`)

### Variable `MOTOR`.
En el `MOTOR` parece que hubiera un *outlier* en su valor m�ximo, que son 8.000 cc. Pero ese modelo existe tal y como se describe en la p�gina web  `https://es.wikipedia.org/wiki/Dodge_Viper`.


### Variable `MPG`.
As� mismo tambi�n es real que el Chevrolet Metro recorre la distancia de 45 'MPG'. Este dato se puede ver en la web  `https://www.fueleconomy.gov/feg/pdfs/guides/FEG2000.pdf`.


```{r}
summary(coches$MPG)
sprintf("Desviaci�n estandar para la variable MPG: %.1f", sd(coches$MPG))
sprintf("Desviacion MPG m�ximo respecto a la media: %.1f", (max(coches$MPG) - mean(coches$MPG)) / sd(coches$MPG))
sprintf("Desviacion MPG 'penultimo' respecto a la media: %.1f", (33 - mean(coches$MPG)) / sd(coches$MPG))
```

## Resumen de los datos.


```{r}
summary(coches)
```

Por tanto se puede decir que en principio todos los datos de la tabla son correctos.

\pagebreak

# Realizar An�lisis Factorial y An�lisis Cluster

**Objetivo: Resumir las variables del fichero en un n�mero m�s reducido de factores y realizar una tipolog�a de los coches, en base a esos factores.**

Se nos proporciona una *matriz de datos* $X$ con $n = 120$ individuos y 14 columnas que ser�n reducidas a $p = 11$ variables cuantitativas. El objetivo del *An�lisis Factorial* es resumir la informaci�n contenida en esa matriz de datos.

## Identificar las variables para el AF.

**Identificar las variables a incluir en el An�lisis Factorial.**

### Normalidad de las variables.


```{r include = FALSE}
library(nortest)
# Pruebas de normalidad Shapiro-Wilk
shp.coches <- rbind(shapiro.test(coches$VENTAS),
                    shapiro.test(coches$REVENTA),
                    shapiro.test(coches$PRECIO),
                    shapiro.test(coches$MOTOR),
                    shapiro.test(coches$CV),
                    shapiro.test(coches$PISADA),
                    shapiro.test(coches$ANCHURA),
                    shapiro.test(coches$LONGITUD),
                    shapiro.test(coches$PESO),
                    shapiro.test(coches$CAPACIDAD),
                    shapiro.test(coches$MPG))

print(shp.coches)

# Pruebas de normalidad Lilliefors
lillie.coches <- rbind(lillie.test(coches$VENTAS),
                    lillie.test(coches$REVENTA),
                    lillie.test(coches$PRECIO),
                    lillie.test(coches$MOTOR),
                    lillie.test(coches$CV),
                    lillie.test(coches$PISADA),
                    lillie.test(coches$ANCHURA),
                    lillie.test(coches$LONGITUD),
                    lillie.test(coches$PESO),
                    lillie.test(coches$CAPACIDAD),
                    lillie.test(coches$MPG))

print(lillie.coches)

# Pruebas de normalidad Kolmogorov-Smirnov
ks.coches <- rbind(ks.test(coches$VENTAS, pnorm, mean(coches$VENTAS), sd(coches$VENTAS)),
                   ks.test(coches$REVENTA, pnorm, mean(coches$REVENTA), sd(coches$REVENTA)),
                   ks.test(coches$PRECIO, pnorm, mean(coches$PRECIO), sd(coches$PRECIO)),
                   ks.test(coches$MOTOR, pnorm, mean(coches$MOTOR), sd(coches$MOTOR)),
                   ks.test(coches$CV, pnorm, mean(coches$CV), sd(coches$CV)),
                   ks.test(coches$PISADA, pnorm, mean(coches$PISADA), sd(coches$PISADA)),
                   ks.test(coches$ANCHURA, pnorm, mean(coches$ANCHURA), sd(coches$ANCHURA)),
                   ks.test(coches$LONGITUD, pnorm, mean(coches$LONGITUD), sd(coches$LONGITUD)),
                   ks.test(coches$PESO, pnorm, mean(coches$PESO), sd(coches$PESO)),
                   ks.test(coches$CAPACIDAD, pnorm, mean(coches$CAPACIDAD), sd(coches$CAPACIDAD)),
                   ks.test(coches$MPG, pnorm, mean(coches$MPG), sd(coches$MPG)))
                   
print(ks.coches)
```

Antes de realizar ning�n an�lisis, se estudia la normalidad de las variables. Para ello s�lo se incluir�n las once variables cuantitativas mencionadas anteriormente: `VENTAS`, `REVENTA`, `PRECIO`, `MOTOR`, `CV`, `PISADA`, `ANCHURA`, `LONGITUD`, `PESO`, `CAPACIDAD` y `MPG`.

Entonces no se usan las variables cualitativas.
```{r echo=FALSE}
# coches$MARCA <- NULL
# coches$MODELO <- NULL
# coches$TIPO <- NULL
```
Se estudia la normalidad de las variables con el test de *Shapiro-Wilk* que es el m�s estricto. En este caso se plantea como hip�tesis nula que la muestra proviene de una poblaci�n normalmente distribuida. Para ello se hace la prueba con cada una de las variables.

Este test nos muestra que solamente hay dos variables con distribuci�n normal: `LONGITUD` y `PESO`,  cuyo *p-valor* > 0.05. Las dem�s no est�n normalmente distribuidas. Por tanto se prueba la normalidad del resto de las variables con el test de *Lilliefors*. Pero con este test, tomando un nivel de significaci�n del 5%, tampoco se obtienen mejores resultados, porque en el mejor de los casos el p-valor vale 0.027. Por �ltimo se prueba con el test de *Kolmogorov-Smirnov*. Se obtiene que todas las variables tienen normalidad excepto `VENTAS`, `REVENTA` y `PRECIO`.



| Variable  | Shapiro-Wilk (p-valor) | Lilliefors (p-valor) | Kolmogorov-Smirnov  (p-valor) |
| ----------|------------------------|----------------------|---------------|
| VENTAS    | 6.8192e-15 | 4.9031e-15 | **2.8197e-05** |
| REVENTA   | 3.5684e-13 | 9.3104e-14 | **6.9501e-05** |
| PRECIO    | 4.1681e-11 | 5.9375e-09 | **0.002153** |
| MOTOR     | 9.6331e-06 | 0.003366   | 0.1586 |
| CV        | 2.0126e-05 | 0.009568   | 0.2274 |
| PISADA    | 6.5897e-06 | 0.002656   | 0.1463 |
| ANCHURA   | 0.009119   | 0.02722   | 0.3281 |
| LONGITUD  | 0.8386     |           |           |
| PESO      | 0.662      |           |           |
| CAPACIDAD | 1.3447e-05 | 0.0007571 | 0.09565 |
| MPG       | 6.1219e-05 | 0.0002374 | 0.06494 |

Se intenta encontrar un transformaci�n de cada una de esas tres variables para conseguir normalidad.

```{r}
par(mfrow=c(1,3))
hist(coches$VENTAS, xlab = "VENTAS", main = "VENTAS")
hist(coches$REVENTA, xlab = "REVENTA", main = "REVENTA")
hist(coches$PRECIO, xlab = "PRECIO", main = "PRECIO")
```

La forma del histograma de la variable `VENTAS` nos puede sugerir que para conseguir la normalidad podemos transformarla a logaritmo neperiano.  Las variables `REVENTA` y  `PRECIO` tienen una distribuci�n mas compleja y nos pueden recordar a una *distribuci�n de Poisson* o quiz�s a una *distribuci�n Gamma*. Pero si se hace una transformaci�n de ese estilo, despu�s la interpretaci�n de los datos es bastante m�s ardua. Por tanto se prueba a transformar cada una de estas tres variables a un sencillo *logaritmo neperiano*.  

```{r echo=TRUE}
coches$lnVENTAS <- log(coches$VENTAS)
coches$lnREVENTA <- log(coches$REVENTA)
coches$lnPRECIO <- log(coches$PRECIO)
```

```{r include = FALSE}
# Pruebas de normalidad Shapiro-Wilks
shp.ln.coches <- rbind(shapiro.test(coches$lnVENTAS),
                    shapiro.test(coches$lnREVENTA),
                    shapiro.test(coches$lnPRECIO))
print(shp.ln.coches)

# Pruebas de normalidad Lilliefors
lillie.ln.coches <- rbind(lillie.test(coches$lnVENTAS),
                    lillie.test(coches$lnREVENTA),
                    lillie.test(coches$lnPRECIO))
print(lillie.ln.coches)

# Pruebas de normalidad Kolmogorov-Smirnov
ks.ln.coches <- rbind(ks.test(coches$lnVENTAS, pnorm, mean(coches$lnVENTAS), sd(coches$lnVENTAS)),
                   ks.test(coches$lnREVENTA, pnorm, mean(coches$lnREVENTA), sd(coches$lnREVENTA)),
                   ks.test(coches$lnPRECIO, pnorm, mean(coches$lnPRECIO), sd(coches$lnPRECIO)))
print(ks.ln.coches)
```

Tabla con los test de normalidad: *p-valores* de las transformaciones. 

| Variable  | Shapiro-Wilk (p-valor) | Lilliefors (p-valor) | Kolmogorov-Smirnov (p-valor) |
|-----------|------------------------|----------------------|---------------|
| lnVENTAS  | 0.0004844 | **0.1210**   | 0.5450 |
| lnREVENTA | 0.0003796 | 0.004019     | **0.1685** |
| lnPRECIO  | 0.01470   | 0.02674      | **0.3260** |

Donde el test de *Lilliefors* ya nos indica que el `ln(VENTAS)` posee una distribuci�n normal. Y el test de *Kolmogorov-Smirnov* nos muestra que las variables `ln(REVENTA)` y `ln(PRECIO)` tambi�n tienen una distribuci�n normal.

Por tanto el resto de los test y an�lisis se realizar� con este **grupo de variables cuantitativas que est�n normalizadas:  `MOTOR`, `CV`, `PISADA`, `ANCHURA`, `LONGITUD`, `PESO`, `CAPACIDAD`, `MPG` y `lnVENTAS`, `lnREVENTA`, `lnPRECIO`**.

```{r echo=TRUE}
# Se crea un nuevo data frame con las variables normalizadas.
# Se incluyen: "MARCA", "MODELO", (NO: "VENTAS") (NO: "REVENTA"), "TIPO",
# (NO: "PRECIO"), "MOTOR", "CV", "PISADA", "ANCHURA",
# "LONGITUD", "PESO", "CAPACIDAD", "MPG", "lnVENTAS",
# "lnREVENTA", "lnPRECIO" 
names(coches)
coches_ln <- coches[,c(1:2,5,7:17)]
names(coches_ln)
```

```{r include=FALSE}
# Se elimina el data frame coches para que no haya problemas.
# coches <- NULL
```


## Comprobar condiciones para AF.

**Comprobar que se verifican las condiciones para realizar un An�lisis Factorial.**

### Determinante de la matriz de correlaciones
Para comprobar si se verifican las condiciones para realizar un AF, hay que hallar el determinante de la matriz de correlaciones. Si el determinante es muy peque�o, significa que las correlaciones entre variables son altas y que puede ser adecuado realizar un an�lisis factorial.[@urkaregi2017apuntes]

```{r echo=TRUE}
# Se halla la matriz de correlaciones de los datos cuantitativos, 
# desde MOTOR hasta lnPRECIO
matriz.cor.coches_ln <- cor(coches_ln[,4:14])
# Se calcula el determinante de la matriz de correlaciones
det.cor.coches_ln <- det(matriz.cor.coches_ln)
# Se muestra su valor
det.cor.coches_ln
```

El determinante de la matriz de correlaciones vale  1.09e-06, que es un valor muy peque�o. O sea que se deduce que las correlaciones entre las variables son altas. Por tanto seg�n este test **ser�a adecuado un an�lisis factorial**.

### Test de de esfericidad de Bartlett
No obstante se puede usar adem�s el Test de Bartlett para comprobar si existen correlaciones en las variables. Este test contrasta la hip�tesis de que la matriz de correlaciones es la matriz identidad. Si es as�, la consecuencia es que las correlaciones entre variables son nulas.
```{r}
library(psych)
```

```{r, echo=TRUE}
# Se toma n (numero de muestras) igual al numero de filas de nuestros datos
cortest.bartlett(matriz.cor.coches_ln, n=nrow(coches_ln))
```

Se obtiene un *p-valor* = 1.22e-292. Con lo cual s� se podr�a rechazar la hip�tesis nula, $H_0$. Por tanto se puede decir que hay correlaciones entre las variables y se confirma una vez m�s que **se puede hacer un an�lisis factorial**.  

### El �ndice Kaiser-Meyer-Olkin (KMO).
Este �ndice compara los coeficientes de correlaci�n y los coeficientes de correlaci�n parciales observados:

$$KMO = \frac{\sum_{i \neq j}\sum_{j \neq i} r_{}ij^2}{\sum_{i \neq j} r_{}ij^2 + \sum_{i \neq j}a_{}ij^2}$$

donde $r_{ij}$ es el *coeficiente de correlaci�n entre las variables $i$ y $j$* y $a_{ij}$ es el *coeficiente de correlaci�n parcial entre las variables $i$ y $j$*.

El coeficiente de correlaci�n parcial entre dos variables representa la relaci�n entre las dos variables, eliminada la influencia del resto de variables. Si las variables tienen factores comunes, los coeficientes de correlaci�n parciales ser�n peque�os, dado que se eliminan los efectos lineales del resto de variables. En consecuencia, para que sea adecuada la utilizaci�n del an�lisis factorial, los coeficientes de correlaci�n parciales tienen que ser peque�os.

Si la suma de los cuadrados de los coeficientes de correlaci�n parciales es peque�a en relaci�n con la suma de los cuadrados de los coeficientes de correlaci�n, el �ndice KMO estar� pr�ximo a 1. Si los valores del �ndice *KMO* son peque�os, no ser� adecuado realizar un an�lisis factorial, dado que las correlaciones entre pares de variables no se pueden explicar por medio de otras variables.

```{r}
kmo.coches_ln <- KMO(matriz.cor.coches_ln)
sprintf("KMO = %.2f", kmo.coches_ln$MSA)
```

Se obtiene que **$KMO$ = 0.82**, que es **alto**, mayor que 0.80. Por tanto s� ser�a **adecuado realizar un an�lisis factorial**.

###  Medida de adecuaci�n de la muestra (MSA)
La **medida de adecuaci�n de la muestra (MSA)** es similar al �ndice KMO, pero en este caso s�lo se incluyen los coeficientes relativos a la variable que queremos comprobar. Si los �ndices MSA son peque�os, no llevaremos a cabo un an�lisis factorial. Si aparecen �ndices *MSA* peque�os s�lo en algunas variables, podemos plantearnos el eliminar esas variables del an�lisis.
```{r}
print(kmo.coches_ln)
```


Los valores *MSA* son en general altos. El valor m�s bajo lo tiene la variable `lnPRECIO` = 0.74, seguido de `CV` = 0.78 y `PISADA` = 0.78. El resto son superiores a 0.8.

O sea que este �ltimo test corrobora de nuevo la idea de que **se puede realizar un an�lisis factorial**. En todo caso se podr�an eliminar las variables con los *MSA* m�s bajos.


## An�lisis Factorial.
**Realizar el An�lisis Factorial y valorar razonadamente las posibles mejoras de este an�lisis**

M�s adelante, observando las comunalidades se justificar� que se pueden plantear al menos tres opciones de an�lisis:

|  Variables | 1: Todas las variables | 2: Todas las variables sin `lnVENTAS` | 3: Todas las variables sin `lnVENTAS` ni `MPG` |
| --------|:--------:|:--------:|:--------:|
|  `MOTOR`     |  x |  x |  x  |
|  `CV`        |  x |  x |  x  |
|  `PISADA`    |  x |  x |  x  |
|  `ANCHURA`   |  x |  x |  x  |
|  `LONGITUD`  |  x |  x |  x  |
|  `PESO`      |  x |  x |  x  |
|  `CAPACIDAD` |  x |  x |  x  |
|  `MPG`       |  x |  x |     |
|  `lnVENTAS`  |  x |    |     |
|  `lnREVENTA` |  x |  x |  x  |
|  `lnPRECIO`  |  x |  x |  x  |

### Opci�n 1. Todas las variables cuantitativas.

Lista de variables: `MOTOR`, `CV`, `PISADA`, `ANCHURA`, `LONGITUD`, `PESO`, `CAPACIDAD`, `MPG`, `lnVENTAS`, `lnREVENTA` y `lnPRECIO`. 

Se hallan los valores propios de la matriz de correlaciones.

```{r}
# eigen(matriz.cor.coches)$values
sprintf("%.2f", eigen(matriz.cor.coches_ln)$values)
```

Hay dos valores propios mayores que 1. Por lo tanto nos quedaremos con 2 factores para realizar el an�lisis.

Se pueden usar diferentes m�todos para calcular los factores. Se elegir� el que explique m�s varianza que en este caso es el m�todo de *Componentes Principales*.

```{r include=FALSE}
# An�lisis Factorial inicial
lista_af = cbind(rbind("minres", "wls", "gls", "pa", "ml", "principal"),
                 rbind(fa(coches_ln[,4:14], nfactors=2, fm="minres")$Vaccounted[3,2],
                       fa(coches_ln[,4:14], nfactors=2, fm="wls")$Vaccounted[3,2],
                       fa(coches_ln[,4:14], nfactors=2, fm="gls")$Vaccounted[3,2],
                 fa(coches_ln[,4:14], nfactors=2, fm="pa")$Vaccounted[3,2],
                 fa(coches_ln[,4:14], nfactors=2, fm="ml")$Vaccounted[3,2],
                 principal(coches_ln[,4:14], nfactors=2, rotate=F, scores=F)$Vaccounted[3,2]))

print(lista_af)
```


| M�todo                                 | Varianza explicada(%) |
|----------------------------------------|------|
| "minres", residuos m�nimos (OLS)       | 75 % |
| "wls", m�nimos cuadrados ponderados    | 76 % |
| "gls", m�nimos cuadrados generalizados | 76 % |
| "pa", factor principal                 | 75 % |
| "ml", M�xima verosimilitud             | 75 % |
| Componentes Principales                | 79 % |


Inicialmente se realiza el An�lisis de Componentes Principales (ACP) sin rotar

```{r}
# M�todo de componentes principales
cp.coches_ln_sin_rotar <- principal(coches_ln[,4:14], nfactors=2, rotate=F, scores=F)
print(cp.coches_ln_sin_rotar)
# plot(cp.coches_ln_sin_rotar)
```

Los dos primeros valores propios explican el 79 % de la varianza. Las comunalidades m�s peque�as corresponden a `lnVENTAS` (0.61) seguido de `MPG` (0.71), `CAPACIDAD` (0.73), `ANCHURA` (0.75), `LONGITUD` (0.77) y `MOTOR` (0.78). Las otras cinco variables tienen comunalidades superiores a 0.80.

Observando las comunalidades se pueden plantear al menos tres opciones de an�lisis:

|  Opci�n | 1: Todas las variables | 2: Todas las variables sin `lnVENTAS` | 3: Todas las variables sin `lnVENTAS` ni `MPG` |
| --------|:--------:|:--------:|:--------:|
|  `MOTOR`     |  x |  x |  x  |
|  `CV`        |  x |  x |  x  |
|  `PISADA`    |  x |  x |  x  |
|  `ANCHURA`   |  x |  x |  x  |
|  `LONGITUD`  |  x |  x |  x  |
|  `PESO`      |  x |  x |  x  |
|  `CAPACIDAD` |  x |  x |  x  |
|  `MPG`       |  x |  x |     |
|  `lnVENTAS`  |  x |    |     |
|  `lnREVENTA` |  x |  x |  x  |
|  `lnPRECIO`  |  x |  x |  x  |


A continuaci�n se muestra muestra el An�lisis de Componentes principales rotados con varimax

```{r}
library(GPArotation)
# Componentes principales rotados
cp.coches_ln <- principal(coches_ln[,4:14], nfactors=2, rotate="varimax", scores=TRUE)
cp.coches_ln
```


|  Variable | 1: Todas las variables. Valor comunalidades ($h^2$) |
| --------|:--------:|
|  `lnPRECIO`  |  0.91 |
|  `PESO`      |  0.88 |
|  `lnREVENTA` |  0.87 |
|  `PISADA`    |  0.85 |
|  `CV`        |  0.84 |
|  `MOTOR`     |  0.78 |
|  `LONGITUD`  |  0.77 |
|  `ANCHURA`   |  0.75 |
|  `CAPACIDAD` |  0.73 |
|  `MPG`       |  0.71 |
|  `lnVENTAS`  |  0.61 |



### Opci�n 2. Todas las variables cuantitativas sin `lnVENTAS`.

Todas las variables cuantitativas sin `lnVENTAS`.

Lista de variables: `MOTOR`, `CV`, `PISADA`, `ANCHURA`, `LONGITUD`, `PESO`, `CAPACIDAD`, `MPG` , `lnREVENTA` y `lnPRECIO`. 



```{r include=FALSE}
# Se crea la matriz con los datos
coches_ln_sin_ventas <-  coches_ln[,c(1:11,13:14)]
names(coches_ln_sin_ventas)

# Se halla la matriz de correlaciones de los datos
matriz.cor.coches_ln_sin_ventas <- cor(coches_ln_sin_ventas[,c(4:13)])
# Se calcula el determinante de la matriz de correlaciones
det.cor.coches_ln_sin_ventas <- det(matriz.cor.coches_ln_sin_ventas)
# Se muestra su valor
det.cor.coches_ln_sin_ventas


library(psych)
# Se toma n (numero de muestras) igual al numero de filas de nuestros datos
cortest.bartlett(matriz.cor.coches_ln_sin_ventas, n=nrow(coches_ln_sin_ventas))


kmo.coches_ln_sin_ventas <- KMO(matriz.cor.coches_ln_sin_ventas)
```

```{r}
# Valor del determinante de la matriz de correlaciones
sprintf("Determinante matriz correlaciones = %.2e", det.cor.coches_ln_sin_ventas)

# Indice KMO
sprintf("KMO = %.2f", kmo.coches_ln_sin_ventas$MSA)
print(kmo.coches_ln_sin_ventas)
```

Se ha calculado la nueva matriz de datos, la matriz de correlaciones y su determinante. Tambi�n se ha obtenido que el �ndice *KMO* = 0.82, que es un valor aceptable.

Se hallan los valores propios de la matriz de correlaciones.

```{r}
# eigen(matriz.cor.coches)$values
sprintf("%.2f", eigen(matriz.cor.coches_ln_sin_ventas)$values)
```

Hay dos valores propios mayores que 1. Por lo tanto nos quedaremos con 2 factores para realizar el an�lisis.

Se pueden usar diferentes m�todos para calcular los factores. Se elegir� el que explique m�s varianza que en este caso es el m�todo de *Componentes Principales*.

```{r include=FALSE}
# An�lisis Factorial inicial
num_factores = 2

lista_af_sin_ventas = cbind(rbind("minres", "wls", "gls", "pa", "ml", "principal"),
                 rbind(fa(coches_ln_sin_ventas[,c(4:13)], nfactors=num_factores, fm="minres")$Vaccounted[3,num_factores],
                       fa(coches_ln_sin_ventas[,c(4:13)], nfactors=num_factores, fm="wls")$Vaccounted[3,num_factores],
                       fa(coches_ln_sin_ventas[,c(4:13)], nfactors=num_factores, fm="gls")$Vaccounted[3,num_factores],
                       fa(coches_ln_sin_ventas[,c(4:13)], nfactors=num_factores, fm="pa")$Vaccounted[3,num_factores],
                       fa(coches_ln_sin_ventas[,c(4:13)], nfactors=num_factores, fm="ml")$Vaccounted[3,num_factores],
                       principal(coches_ln_sin_ventas[,c(4:13)], nfactors=num_factores, rotate=F, scores=F)$Vaccounted[3,num_factores]))

print(lista_af_sin_ventas)
```


| M�todo                                 | Varianza explicada(%) |
|----------------------------------------|------|
| "minres", residuos m�nimos (OLS)       | 78 % |
| "wls", m�nimos cuadrados ponderados    | 80 % |
| "gls", m�nimos cuadrados generalizados | 80 % |
| "pa", factor principal                 | 78 % |
| "ml", M�xima verosimilitud             | 78 % |
| Componentes Principales                | 82 % |


Inicialmente se realiza el An�lisis de Componentes Principales (ACP) sin rotar

```{r}
# M�todo de componentes principales
principal(coches_ln_sin_ventas[,c(4:13)], nfactors=num_factores, rotate=F, scores=F)
```

Los dos primeros valores propios explican el 82 % de la varianza. Las comunalidades m�s peque�as corresponden a `MPG` (0.71), `CAPACIDAD` (0.74), `ANCHURA` (0.76), `MOTOR` (0.79) y `LONGITUD` (0.80). Las otras cinco variables tienen comunalidades superiores a 0.80.

En la siguiente tabla se puede ver que con este an�lisis mejoran las comunalidades:


|  Variable | 1: Todas las variables. Valor comunalidades ($h^2$) | 2: Todas las variables sin `lnVENTAS`. Valor comunalidades ($h^2$) |
| --------|:--------:|:--------:|
|  `lnPRECIO`  |  0.91 |  0.93 |
|  `PESO`      |  0.88 |  0.88 |
|  `lnREVENTA` |  0.87 |  0.89 |
|  `PISADA`    |  0.85 |  0.86 |
|  `CV`        |  0.84 |  0.87 |
|  `MOTOR`     |  0.78 |  0.79 |
|  `LONGITUD`  |  0.77 |  0.80 |
|  `ANCHURA`   |  0.75 |  0.76 |
|  `CAPACIDAD` |  0.73 |  0.74 |
|  `MPG`       |  0.71 |  0.71 |
|  `lnVENTAS`  |  0.61 |  --- |



* Porcentaje de varianza explicada. Opci�n 1 (Todas las variables) = **79 %**.
* Porcentaje de varianza explicada. Opci�n 2 (Todas las variables sin `lnVENTAS`) = **82 %**.



### Opci�n 3. Sin `lnVENTAS` y sin `MPG`.

Todas las variables cuantitativas sin `lnVENTAS` y sin `MPG`.

Lista de variables: `MOTOR`, `CV`, `PISADA`, `ANCHURA`, `LONGITUD`, `PESO`, `CAPACIDAD` , `lnREVENTA` y `lnPRECIO`. 


```{r include=FALSE}
names(coches_ln)
# Se crea la matriz con los datos
coches_ln_sin_vtas_mpg <-  coches_ln[,c(4:10,13:14)]
names(coches_ln_sin_vtas_mpg)


# Se halla la matriz de correlaciones de los datos
matriz.cor.coches_ln_sin_vtas_mpg <- cor(coches_ln_sin_vtas_mpg)
# Se calcula el determinante de la matriz de correlaciones
det.cor.coches_ln_sin_vtas_mpg <- det(matriz.cor.coches_ln_sin_vtas_mpg)
# Se muestra su valor
det.cor.coches_ln_sin_vtas_mpg


library(psych)
# Se toma n (numero de muestras) igual al numero de filas de nuestros datos
cortest.bartlett(matriz.cor.coches_ln_sin_vtas_mpg, n=nrow(coches_ln_sin_vtas_mpg))


kmo.coches_ln_sin_vtas_mpg <- KMO(matriz.cor.coches_ln_sin_vtas_mpg)
sprintf("KMO = %.2f", kmo.coches_ln_sin_vtas_mpg$MSA)


print(kmo.coches_ln_sin_vtas_mpg)
```



```{r}
# Valor del determinante de la matriz de correlaciones
sprintf("Determinante matriz correlaciones = %.2e", det.cor.coches_ln_sin_vtas_mpg)

# Indice KMO
sprintf("KMO = %.2f", kmo.coches_ln_sin_vtas_mpg$MSA)
print(kmo.coches_ln_sin_vtas_mpg)
```

Se ha calculado la nueva matriz de datos, la matriz de correlaciones y su determinante. Tambi�n se ha obtenido que el �ndice *KMO* = 0.78, que no es un valor muy aceptable dado que es inferior a 0.80.

Se hallan los valores propios de la matriz de correlaciones.




Despu�s de crear la nueva matriz de datos, la matriz de correlaciones y su determinante, se hallan los valores propios de la matriz de correlaciones.

```{r}
# eigen(matriz.cor.coches)$values
sprintf("%.2f", eigen(matriz.cor.coches_ln_sin_vtas_mpg)$values)
```

Hay dos valores propios mayores que 1. Por lo tanto nos quedaremos con 2 factores para realizar el an�lisis.

Se pueden usar diferentes m�todos para calcular los factores. Se elegir� el que explique m�s varianza que en este caso es el m�todo de *Componentes Principales*.

```{r include=FALSE}
# An�lisis Factorial inicial
# An�lisis Factorial inicial
num_factores = 2

lista_af_sin_vtas_mpg = cbind(rbind("minres", "wls", "gls", "pa", "ml", "principal"),
                            rbind(fa(coches_ln_sin_vtas_mpg, nfactors=num_factores, fm="minres")$Vaccounted[3,num_factores],
                                  fa(coches_ln_sin_vtas_mpg, nfactors=num_factores, fm="wls")$Vaccounted[3,num_factores],
                                  fa(coches_ln_sin_vtas_mpg, nfactors=num_factores, fm="gls")$Vaccounted[3,num_factores],
                                  fa(coches_ln_sin_vtas_mpg, nfactors=num_factores, fm="pa")$Vaccounted[3,num_factores],
                                  fa(coches_ln_sin_vtas_mpg, nfactors=num_factores, fm="ml")$Vaccounted[3,num_factores],
                                  principal(coches_ln_sin_vtas_mpg, nfactors=num_factores, rotate=F, scores=F)$Vaccounted[3,num_factores]))

print(lista_af_sin_vtas_mpg)
```


| M�todo                                 | Varianza explicada(%) |
|----------------------------------------|------|
| "minres", residuos m�nimos (OLS)       | 80 % |
| "wls", m�nimos cuadrados ponderados    | 81 % |
| "gls", m�nimos cuadrados generalizados | 81 % |
| "pa", factor principal                 | 80 % |
| "ml", M�xima verosimilitud             | 79 % |
| Componentes Principales                | 84 % |


Inicialmente se realiza el An�lisis de Componentes Principales (ACP) sin rotar

```{r}
# M�todo de componentes principales
principal(coches_ln_sin_vtas_mpg, nfactors=num_factores, rotate=F, scores=F)
```

Los dos primeros valores propios explican el 84 % de la varianza. Las comunalidades m�s peque�as corresponden a `CAPACIDAD` (70), `ANCHURA` (0.78) y `MOTOR` (0.79). Las otras variables tienen comunalidades superiores a 0.80.


En el siguiente apartado se discute cual puede ser la mejor opci�n.


## An�lisis Factorial definitivo.

**Realizar el An�lisis Factorial definitivo (desde tu punto de vista)**

A continuaci�n se muestra un resumen de los datos obtenidos hasta ahora y se eval�a cual puede ser el mejor an�lisis.

Adem�s se a�ade un resumen de una cuarta opci�n: Todas las variables sin `lnVENTAS`, ni `MPG`, ni tampoco la variable `CAPACIDAD` cuyo $h^2$ tambi�n es bajo.
```{r include=FALSE}
# Se crea la matriz con los datos de la opci�n 4
coches_ln_sin_vtas_mpg_cap <-  coches_ln[,c(4:9,13:14)]
```

### �ndice KMO.
* �ndice KMO. Opci�n 1 (Todas las variables) = **0.82**.
* �ndice KMO. Opci�n 2 (Todas las variables sin `lnVENTAS`) = **0.82**.
* �ndice KMO. Opci�n 3 (Todas las variables sin `lnVENTAS` y sin `MPG`) = **0.78**.
* �ndice KMO. Opci�n 4 (Todas las variables sin `lnVENTAS`, sin `MPG` y sin `MPG`) = **0.76**.

El �ndice KMO muestra las correlaciones entre pares de variables. Se aprecia que el �ndice KMO es inferior a 0.80 en la tercera y cuarta opci�n. O sea que **para las opciones 3 y 4 no ser�a adecuado realizar un an�lisis factorial**.


### �ndices MSA.

```{r include=FALSE}
# Se crea la matriz con los datos
coches_ln_sin_vtas_mpg_cap <-  coches_ln[,c(4:9,13:14)]
# Se halla la matriz de correlaciones de los datos
matriz.cor.coches_ln_sin_vtas_mpg_cap <- cor(coches_ln_sin_vtas_mpg_cap)
kmo.coches_ln_sin_vtas_mpg_cap <- KMO(matriz.cor.coches_ln_sin_vtas_mpg_cap)


kmo.coches_ln$MSAi
kmo.coches_ln_sin_ventas$MSAi
kmo.coches_ln_sin_vtas_mpg$MSAi
kmo.coches_ln_sin_vtas_mpg_cap$MSAi
```


|  Variable | 1: Todas las variables. (Valor MSA) | 2: Todas las variables sin `lnVENTAS`. (Valor MSA) | 3: Todas las variables sin `lnVENTAS` ni `MPG`. (Valor MSA) | 4: Todas las variables sin `lnVENTAS` ni `MPG` ni `CAPACIDAD`. (Valor MSA) |
|-----------|:----:|:----:|:----:|:----:|
| ANCHURA   | 0.95 | 0.95 | 0.95 | 0.94 |
|    MPG    | 0.92 | 0.92 |  ---  | ---  |
| CAPACIDAD | 0.85 | 0.86 | 0.82 | ---  |
| lnVENTAS  | 0.85 |  --- |  --- |  --- |
| MOTOR     | 0.83 | 0.83 | 0.80 | 0.76 |
|   PESO    | 0.82 | 0.82 | 0.77 | 0.76 |
| LONGITUD  | 0.79 | 0.79 | 0.76 | 0.76 |
| lnREVENTA | 0.79 | 0.77 | 0.75 | 0.74 |
| CV        | 0.78 | 0.77 | 0.74 | 0.73 |
| PISADA    | 0.78 | 0.79 | 0.78 | 0.72 |
| lnPRECIO  | 0.74 | 0.72 | 0.70 | 0.69 |

Vemos que los �ndices MSA van disminuyendo seg�n se eliminan variables de manera que para las opciones 3 y 4 tiene demasiados valores por debajo de 0.80. Con lo cual desde el punto de vista de los �ndices MSA **para las opciones 3 y 4 no ser�a adecuado realizar un an�lisis factorial**.


### Comunalidades $h^2$.

La siguiente tabla muestra las comunalidades para las cuatro opciones.

|  Variable | 1: Todas las variables. Valor comunalidades ($h^2$) | 2: Todas las variables sin `lnVENTAS`. Valor comunalidades ($h^2$) | 3: Todas las variables sin `lnVENTAS` ni `MPG`. Valor comunalidades ($h^2$) | 4: Todas las variables sin `lnVENTAS` ni `MPG`  ni `CAPACIDAD`. Valor comunalidades ($h^2$) |
| --------|:--------:|:--------:|:--------:|:--------:|
|  `lnPRECIO`  |  0.91 |  0.93 |  0.94  |  0.94  |
|  `PESO`      |  0.88 |  0.88 |  0.86  |  0.82  |
|  `lnREVENTA` |  0.87 |  0.89 |  0.90  |  0.89  |
|  `PISADA`    |  0.85 |  0.86 |  0.88  |  0.90  |
|  `CV`        |  0.84 |  0.87 |  0.89  |  0.80  |
|  `MOTOR`     |  0.78 |  0.79 |  0.79  |  0.87  |
|  `LONGITUD`  |  0.77 |  0.80 |  0.82  |  0.87  |
|  `ANCHURA`   |  0.75 |  0.76 |  0.78  |  0.79  |
|  `CAPACIDAD` |  0.73 |  0.74 |  0.70  |  ---  |
|  `MPG`       |  0.71 |  0.71 |  ---   |  ---  |
|  `lnVENTAS`  |  0.61 |  ---  |  ---   |  ---  |


Para las tres primeras opciones se encuentran comunalidades menores de 0.80, pero adem�s en la opci�n 1 hay un $h^2$ = 0.61, que es muy baja respecto al resto. Las comunalidades aumentan seg�n se van eliminando variables, *pero se va perdiendo informaci�n*. El mayor cambio sucede en la opci�n 2, dado que all� se ha eliminado el valor de $h^2$ = 0.61 para `lnVENTAS`. O sea que en este caso la **opci�n 2** podr�a ser **la mejor**.


### Varianza total explicada

* Porcentaje de varianza explicada. Opci�n 1 (Todas las variables) = **79 %**.
* Porcentaje de varianza explicada. Opci�n 2 (Todas las variables sin `lnVENTAS`) = **82 %**.
* Porcentaje de varianza explicada. Opci�n 3 (Todas las variables sin `lnVENTAS` y sin `MPG`) = **84 %**.
* Porcentaje de varianza explicada. Opci�n 4 (Todas las variables sin `lnVENTAS`, sin `MPG` y sin `MPG`) = **86 %**.

**Cualquiera de las opciones explica gran cantidad de varianza** y la varianza explicada aumenta a medida que se eliminan variables del an�lisis; lo cual es de esperar.


### Diagramas de componentes y variables

Al ir eliminando las variables el contenido de los dos factores no var�a. Por tanto estos diagramas no nos dan informaci�n adicional para tomar la decisi�n final.


```{r echo=FALSE}
par(mfrow=c(1,2))

fa.diagram(cp.coches_ln, digits = 3, main = "Opci�n 1", rsize=.37)

cp.coches_ln_sin_ventas <- principal(coches_ln_sin_ventas[,c(4:13)], nfactors=num_factores, rotate="varimax", scores=TRUE)
fa.diagram(cp.coches_ln_sin_ventas, digits = 3, main = "Opci�n 2. Sin ventas", rsize=.37)

par(mfrow=c(1,2))
cp.coches_ln_sin_vtas_mpg <- principal(coches_ln_sin_vtas_mpg, nfactors=num_factores, rotate="varimax", scores=TRUE)
fa.diagram(cp.coches_ln_sin_vtas_mpg, digits = 3, main = "Opci�n 3. Sin ventas, mpg", rsize=.37)

cp.coches_ln_sin_vtas_mpg_cap <- principal(coches_ln_sin_vtas_mpg_cap, nfactors=num_factores, rotate="varimax", scores=TRUE)
fa.diagram(cp.coches_ln_sin_vtas_mpg_cap, digits = 3, main = "Opci�n 4. Sin ventas, mpg, capac.", rsize=.37)
```


**Conclusi�n:**
Decido tomar la **opci�n 2**, eliminando �nicamente la variable `log(VENTAS)`. 


* Lista de variables: `MOTOR`, `CV`, `PISADA`, `ANCHURA`, `LONGITUD`, `PESO`, `CAPACIDAD`, `MPG` , `lnREVENTA` y `lnPRECIO`. 
* El �ndice KMO = 0.82, que es un valor adecuado.
* Tiene los mejores valores para los �ndices MSA de entre las cuatro opciones.
* En esta opci�n se elimina la comunalidad m�s baja (0.61) sin perder un gran n�mero de variables y de informaci�n.
* Porcentaje de varianza explicada de todas las variables es del 79 %.

## Interpretar los factores obtenidos.

Se muestran los datos de las componentes principales rotadas por *varimax* para interpretar los factores.

```{r}
# Se muestran las componentes principales rotadas
# Uso un alias para acortar la notacion
cargas <- cp.coches_ln_sin_ventas$loadings

plot(cargas[,2]~cargas[,1],  xlim=c(-1, 1.2),ylim=c(-1, 1), pch = 16,
     xlab = list('CPR 1', col="red"), ylab = list('CPR 2', col="blue"), 
     main = 'Componentes principales rotadas',
     col=ifelse(colnames(cp.coches_ln_sin_ventas$Vaccounted)=="RC1","red","blue"))
text(cargas[,2]~cargas[,1], labels = row.names(cargas), pos = 4)
abline(v=0, h=0)
abline(v=1, h=1, lty=3)
abline(v=-1, h=-1, lty=3)

# Se muestran los factores y sus variables
fa.diagram(cp.coches_ln_sin_ventas$loadings, digits = 3)

# Resumen de las componentes principales
cp.coches_ln_sin_ventas$weights
cp.coches_ln_sin_ventas$Vaccounted
```

* El factor `RC1` explica el 51 % de la varianza. Est� relacionado con `PISADA` (Base del neum�tico), `LONGITUD` del veh�culo, `ANCHURA`  del veh�culo, `PESO` del veh�culo, `CAPACIDAD` del tanque de combustible y tiene correlaci�n negativa con `MPG` (millas recorridas por gal�n de combustible). Por tanto este primer factor se podr�a definir como *"Tama�o y aspecto externo del veh�culo"*.
* El factor `RC2` explica el 49 % de la varianza. Est� relacionado con `log(PRECIO)`  del veh�culo, `log(REVENTA)` precio de reventa del veh�culo a los 4 a�os, `CV` (potencia del motor en CV) y con `MOTOR` (tama�o del  motor). Se podr�a definir como *"Potencia y precio del veh�culo"*. Aunque hay que se�alar que debido a la presencia de los logaritmos en este factor, los veh�culos m�s baratos tienen mayor preponderancia que los caros.

## Coches en el plano factorial.
**Representar los coches en el plano factorial**

Se representan el plano factorial de los coches de tres formas diferentes: 

* Plano1: destacando los precios de los coches en tres grupos "baratos", "precio medio" y "caros" (sin logaritmos en este caso).
* Plano2: destacando el tama�o del motor.
* Plano3: destacando el consumo (MPG).

Y se aprecia que en general, los coches m�s baratos tienen un motor m�s peque�o y consumen menos combustible.

```{r include=FALSE}
# Se guardan las componentes rotadas en la matriz de datos
coches_ln_sin_ventas$CPR1 <- cp.coches_ln_sin_ventas$scores[,1]
coches_ln_sin_ventas$CPR2 <- cp.coches_ln_sin_ventas$scores[,2]

library(ggplot2)
library(qgraph)
```
```{r include=FALSE}
# Se muestra la marca del coche.
# Para poner un n�mero de columna sustituir
# coches_ln_sin_ventas$MARCA por rownames(coches_ln_sin_ventas)
ggplot(coches_ln_sin_ventas, aes(x=CPR1, y=CPR2)) +
  geom_point() + geom_rug() +
  xlim(-2.5, 3.5) + ylim(-2.1,4.2)+
  geom_vline(xintercept = 0) + geom_hline(yintercept = 0) +
  geom_text(aes(label=ifelse((CPR2 > 1 | CPR2 < -1 | CPR1 > 0.8 | CPR1 < -1.1) ,as.character(coches_ln_sin_ventas$MARCA),'')),hjust=0,vjust=0, size=3.3)  +
  scale_colour_brewer(palette = "Set1")
```


```{r}
# Precio medio
precio_med <- function() 
{
  ifelse(exp(coches_ln_sin_ventas$lnPRECIO) < mean(exp(coches_ln_sin_ventas$lnPRECIO)),"baratos","caros")
}

precio_med_3 <- function() 
{
  margen = 0.1
  
  txt1 = paste("baratos", 100*(1 - 2*margen) / 2, "%")
  txt2 = paste("medio", 2*100*margen, "%")
  txt3 = paste("caros", 100*(1 - 2*margen) / 2, "%")
  
  ifelse(exp(coches_ln_sin_ventas$lnPRECIO) < ((1 - margen) * mean(exp(coches_ln_sin_ventas$lnPRECIO))),txt1,
         ifelse(exp(coches_ln_sin_ventas$lnPRECIO) > ((1 + margen) * mean(exp(coches_ln_sin_ventas$lnPRECIO))),txt3,txt2))
}

# ggplot(coches_ln_sin_ventas, aes(x=CPR1, y=CPR2, color=precio_med_3())) +
#   geom_point() + geom_rug() +
#   xlim(-2.5, 3.5) + ylim(-2.1,4.2)+
#   geom_vline(xintercept = 0) + geom_hline(yintercept = 0) +
#   geom_text(aes(label=ifelse((CPR2 > 1 | CPR2 < -1 | CPR1 > 0.8 | CPR1 < -1.1),as.character(rownames(coches_ln_sin_ventas)),'')),hjust=0,vjust=0, size=3.3)  +
#   scale_colour_brewer(palette = "Set1")


# library(ggrepel)
# Cambiar la funcion geom_text() por geom_text_repel()
ggplot(coches_ln_sin_ventas, aes(x=CPR1, y=CPR2, color=precio_med_3())) +
  geom_point() + geom_rug() +
  xlim(-2.5, 3.5) + ylim(-2.1,4.2)+
  geom_vline(xintercept = 0) + geom_hline(yintercept = 0) +
  geom_text(aes(label=ifelse((CPR2 > 1 | CPR2 < -1 | CPR1 > 0.8 | CPR1 < -1.1),as.character(paste(rownames(coches_ln_sin_ventas), coches_ln_sin_ventas$MARCA)),'')),hjust=0,vjust=0, size=3.3)  +
  scale_colour_brewer(palette = "Set1")
```

```{r include=FALSE}
#Tipo de coche
ggplot(coches_ln_sin_ventas, aes(x=CPR1, y=CPR2, color=TIPO)) +
  geom_point() + geom_rug() +
  xlim(-2.5, 3.5) + ylim(-2.1,4.2)+
  geom_vline(xintercept = 0) + geom_hline(yintercept = 0) +
  geom_text(aes(label=ifelse((CPR2 > 1 | CPR2 < -1 | CPR1 > 0.8 | CPR1 < -1.1),as.character(paste(rownames(coches_ln_sin_ventas), coches_ln_sin_ventas$MARCA)),'')),hjust=0,vjust=0, size=3.3)  +
  scale_colour_brewer(palette = "Set1")
```



```{r include=FALSE}
pisada_med <- function() 
{
  ifelse(coches_ln_sin_ventas$PISADA < mean(coches_ln_sin_ventas$PISADA),"estrecho","ancho")
}

pisada_med_3 <- function() 
{
  margen = 0.05
  
  txt1 = paste("estrecho", 100*(1 - 2*margen) / 2, "%")
  txt2 = paste("medio", 2*100*margen, "%")
  txt3 = paste("ancho", 100*(1 - 2*margen) / 2, "%")
  
  ifelse(coches_ln_sin_ventas$PISADA < ((1 - margen) * mean(coches_ln_sin_ventas$PISADA)),txt1,
         ifelse(coches_ln_sin_ventas$PISADA > ((1 + margen) * mean(coches_ln_sin_ventas$PISADA)),txt3,txt2))
}


ggplot(coches_ln_sin_ventas, aes(x=CPR1, y=CPR2, color=pisada_med_3())) +
  geom_point() + geom_rug() +
  xlim(-2.5, 3.5) + ylim(-2.1,4.2)+
  geom_vline(xintercept = 0) + geom_hline(yintercept = 0) +
  geom_text(aes(label=ifelse((CPR2 > 1 | CPR2 < -1 | CPR1 > 0.8 | CPR1 < -1.1),as.character(paste(rownames(coches_ln_sin_ventas), coches_ln_sin_ventas$MARCA)),'')),hjust=0,vjust=0, size=3.3)  +
  scale_colour_brewer(palette = "Set1")
```



```{r}
motor_med <- function() 
{
  ifelse(coches_ln_sin_ventas$MOTOR < mean(coches_ln_sin_ventas$MOTOR),"grande","peque�o")
}

motor_med_3 <- function() 
{
  margen = 0.1
  
  txt1 = paste("peque�o", 100*(1 - 2*margen) / 2, "%")
  txt2 = paste("medio", 2*100*margen, "%")
  txt3 = paste("grande", 100*(1 - 2*margen) / 2, "%")
  
  ifelse(coches_ln_sin_ventas$MOTOR < ((1 - margen) * mean(coches_ln_sin_ventas$MOTOR)),txt1,
         ifelse(coches_ln_sin_ventas$MOTOR > ((1 + margen) * mean(coches_ln_sin_ventas$MOTOR)),txt3,txt2))
}


ggplot(coches_ln_sin_ventas, aes(x=CPR1, y=CPR2, color=motor_med_3())) +
  geom_point() + geom_rug() +
  xlim(-2.5, 3.5) + ylim(-2.1,4.2)+
  geom_vline(xintercept = 0) + geom_hline(yintercept = 0) +
  geom_text(aes(label=ifelse((CPR2 > 1 | CPR2 < -1 | CPR1 > 0.8 | CPR1 < -1.1),as.character(paste(rownames(coches_ln_sin_ventas), coches_ln_sin_ventas$MARCA)),'')),hjust=0,vjust=0, size=3.3)  +
  scale_colour_brewer(palette = "Set1")
```


```{r include=FALSE}
# Precio reventa al cabo de 4 a�os
reventa_med <- function() 
{
  ifelse(exp(coches_ln_sin_ventas$lnREVENTA)<mean(exp(coches_ln_sin_ventas$lnREVENTA)),"bajo","alto")
}

reventa_med_3 <- function() 
{
  margen = 0.1
  
  txt1 = paste("bajo", 100*(1 - 2*margen) / 2, "%")
  txt2 = paste("medio", 2*100*margen, "%")
  txt3 = paste("alto", 100*(1 - 2*margen) / 2, "%")
  
  ifelse(exp(coches_ln_sin_ventas$lnREVENTA) < ((1 - margen) * mean(exp(coches_ln_sin_ventas$lnREVENTA))),txt1,
         ifelse(exp(coches_ln_sin_ventas$lnREVENTA) > ((1 + margen) * mean(exp(coches_ln_sin_ventas$lnREVENTA))),txt3,txt2))
}


ggplot(coches_ln_sin_ventas, aes(x=CPR1, y=CPR2, color=reventa_med_3())) +
  geom_point() + geom_rug() +
  xlim(-2.5, 3.5) + ylim(-2.1,4.2)+
  geom_vline(xintercept = 0) + geom_hline(yintercept = 0) +
  geom_text(aes(label=ifelse((CPR2 > 1 | CPR2 < -1 | CPR1 > 0.8 | CPR1 < -1.1),as.character(rownames(coches_ln_sin_ventas)),'')),hjust=0,vjust=0, size=3.3)  +
  scale_colour_brewer(palette = "Set1")
```


```{r}
# Consumo MPG
mpg_med <- function() 
{
  ifelse(coches_ln_sin_ventas$MPG < mean(coches_ln_sin_ventas$MPG),"bajo","alto")
}

mpg_med_3 <- function() 
{
  margen = 0.1
  
  txt1 = paste("alto", 100*(1 - 2*margen) / 2, "%")
  txt2 = paste("medio", 2*100*margen, "%")
  txt3 = paste("bajo", 100*(1 - 2*margen) / 2, "%")
  
  ifelse(coches_ln_sin_ventas$MPG < ((1 - margen) * mean(coches_ln_sin_ventas$MPG)),txt1,
         ifelse(coches_ln_sin_ventas$MPG > ((1 + margen) * mean(coches_ln_sin_ventas$MPG)),txt3,txt2))
}


ggplot(coches_ln_sin_ventas, aes(x=CPR1, y=CPR2, color=mpg_med_3())) +
  geom_point() + geom_rug() +
  xlim(-2.5, 3.5) + ylim(-2.1,4.2)+
  geom_vline(xintercept = 0) + geom_hline(yintercept = 0) +
  geom_text(aes(label=ifelse((CPR2 > 1 | CPR2 < -1 | CPR1 > 0.8 | CPR1 < -1.1),as.character(paste(rownames(coches_ln_sin_ventas), coches_ln_sin_ventas$MARCA)),'')),hjust=0,vjust=0, size=3.3)  +
  scale_colour_brewer(palette = "Set1")
```


## Caracter�sticas m�s destacables de los coches.

**Se�alar las caracter�sticas m�s destacables de los coches en base a los factores obtenidos.**

En general, las caracter�sticas m�s destacables de los coches en base a los factores son:

* Seg�n el factor `RC1`: los veh�culos de m�s tama�o y peso consumen m�s combustible.
* Seg�n el factor `RC2`: los veh�culos con motor m�s grande y mayor potencia tienen un precio m�s alto y su reventa tambi�n es a precio m�s alto.


## Cluster jer�rquico.

**Realizar un cluster jer�rquico, tomando como variables las coordenadas factoriales obtenidas en el an�lisis anterior.**
 
 
Antes de realizar un *an�lisis cluster* deben tomarse tres decisiones:
1. Selecci�n de las variables relevantes para identificar a los grupos. Esto ya se ha hecho en los apartados previos.
2. Elecci�n de la medida de proximidad entre los individuos. En este caso se elegir� la distancia eucl�dea al cuadrado, que es la requieren los m�todos *centroide*, de la *mediana* y  de *Ward*.
3. Elecci�n del criterio para agrupar individuos. En este caso a�n est� por decidir


 Existen diferentes m�todos de aglomeraci�n, entre ellos: *"distancias m�nimas"*, *"distancias m�ximas"*, *"promedio entre grupos"*, *"media ponderada"*, *"centroide"*, *"mediana"* o *"ward"*. Se puede probar con cada uno de ellos y ver c�mo es el dendogr�ma para decidir por el m�s adecuado a nuestros prop�sitos. 
 
 Adem�s se puede buscar el n�mero �ptimo de clusters (`k`) usando la funci�n `daisy` (*Dissimilarity Matrix Calculation*) que computa las desigualdades en las distancias entre parejas. En este caso se ha usado la m�trica *euclidean* para ese c�lculo.

A continuaci�n se muestra un resumen de todo ello:
 
 
 
 
```{r}
# Funcion para clacular el k optimo
# Creamos una funcion para calcular la matriz de distancia dos a dos con los grupos
# daisy: Compute all the pairwise dissimilarities (distances) between observations in the data set.
grpdist <- function(X, metrica="gower"){
  # metrica = "euclidean", "gower"
  # Se carga la libreria "cluster"
  require(cluster)
  gr <- as.data.frame(as.factor(X))
  distgr <- daisy(gr, metrica) # euclidean, gower
  distgr
  }
```

```{r}

# Distancia euclidea al cuadrado entre los elementos de la matriz
coches_ln_sin_ventas.dist<-dist(model.matrix(~-1 + CPR1+CPR2, coches_ln_sin_ventas))^2
# Lista de metodos a testear
lista_mtdos = c("ward.D", "complete", "single", "average", "mcquitty", "median", "centroid")
# data frame auxialiar
kt <- data.frame(k=1:nrow(coches_ln_sin_ventas), r=0)

# Bucle para testear todos los metodos con la distancia elegida
for (mtdo_usado in lista_mtdos)
{
  HClust_aux <- hclust(coches_ln_sin_ventas.dist , method= mtdo_usado)
  # plot(HClust_aux, main= "Cluster Dendrograma para coches", xlab= "Veh�culos estudiados",sub=paste("Metodo=", mtdo_usado,"; Distancia=euclidea al cuadrado", sep = ""))

  for(i in 2:(nrow(coches_ln_sin_ventas)-1)){
    gr <- cutree(HClust_aux, i)
    distgr <- grpdist(gr, "euclidean")
    mt <- cor(coches_ln_sin_ventas.dist, distgr, method="pearson")
    kt[i,2] <- mt
    }

  cat('\nM�todo: "', mtdo_usado,'". k �ptimo = ',k.best <- which.max(kt$r),'\nComposici�n de los clusters:\n', sep = "")

  print(summary(as.factor(cutree(HClust_aux, k = k.best))))

}
```
 
 
 * Est� claro que los **peores m�todos** para nuestro an�lisis son: "single" (single linkage o distancias m�nimas), "average" (promedio entre grupos), "mcquitty" (Similarity Analysis by Reciprocal Pairs for Discrete and Continuous Data), "median" (WPGMC) y "centroid" (UPGMC). En estos casos se crea un cluster muy grande con casi todos los elementos (> 90%) y el resto de clusters tiene menos del 10 % de los elementos. Esto impide una buena clasificaci�n.
 * Los **mejores m�todos** para nuestro an�lisis pueden ser "ward" y "complete" ( complete linkage method o distancias m�ximas). Ambos m�todos generan tres primeros grupos similares en tama�o. Pero el m�todo complete genera un cluster de un solo elemento. Por tanto, decido usar el **m�todo ward**  con las distancias euclideas al cuadrado.


## Decidir el n�mero de cluster y caracterizar estos cluster.

A�n est� por decicir el criterio para agrupar individuos y  el n�mero de clusters.

Se muestra como ejemplo una divisi�n cluster por el m�todo *Ward* o tambi�n conocido como *"p�rdida de inercia m�nima"*

```{r}
## Realizamos un analisis Cluster
# the agglomeration method to be used. 
# This should be (an unambiguous abbreviation of) one of 
# "ward.D", "ward.D2", "single", "complete", "average" (= UPGMA), 
# "mcquitty" (= WPGMA), "median" (= WPGMC) or "centroid" (= UPGMC).
mtdo_usado = "ward.D" # "ward", "single", "complete"

coches_ln_sin_ventas.dist<-dist(model.matrix(~-1 + CPR1+CPR2, coches_ln_sin_ventas))^2
HClust_coches_ln_sin_ventas <- hclust(coches_ln_sin_ventas.dist , method= mtdo_usado)
plot(HClust_coches_ln_sin_ventas, main= "Cluster Dendrograma para coches", xlab= "Veh�culos estudiados", sub=paste("Metodo=", mtdo_usado,"; Distancia=euclidea al cuadrado", sep = ""))
```


Se pueden apreciar 4 grupos en el dendograma. Por otro lado se busca el �ptimo usando la funci�n `daisy` (*Dissimilarity Matrix Calculation*) tal y como se mencion� anteriormente.

```{r message=FALSE, warning=FALSE}
kt <- data.frame(k=1:nrow(coches_ln_sin_ventas), r=0)
for(i in 2:(nrow(coches_ln_sin_ventas)-1)){
  gr <- cutree(HClust_coches_ln_sin_ventas, i)
  # Se halla las distancias entre parejas
  distgr <- grpdist(gr, "euclidean")
  mt <- cor(coches_ln_sin_ventas.dist, distgr, method="pearson")
  kt[i,2] <- mt}
# Se muestra la distancia optima
cat('Busqueda iterativa mediante distancias euclideas =',k.best <- which.max(kt$r))
```


Nos da como resultado 5 clusters

Se muestra la composici�n de los de los datos con 4 y con 5 clusters:
```{r}
cat("Composicion del dendogr�ma con 4 clusters\n")
summary(as.factor(cutree(HClust_coches_ln_sin_ventas, k = 4))) # Cluster Sizes
cat("Composicion del dendogr�ma con 5 clusters\n")
summary(as.factor(cutree(HClust_coches_ln_sin_ventas, k = 5))) # Cluster Sizes
```

Hay poca variaci�n entre la opci�n con 4 o con 5 clusters. Sencillamente desaparecen 3 elementos del tercer grupo y se crea un quinto grupo. Se realizar�n los dos an�lisis y se comparar�n resultados.


### Opci�n 1: Cuatro clusters.

```{r include=FALSE}
# k = 4
num_clusters = 4
```

```{r}
(centroides4<-by(model.matrix(~-1 + CPR1+CPR2, coches_ln_sin_ventas), as.factor(cutree(HClust_coches_ln_sin_ventas, k = num_clusters)), colMeans))
```

* En el primer grupo (CPR1 $\approx$ -0.70 , CPR2 $\approx$ -0.74): 
    + Tama�o: bajo  e intermedio.
    + Precio: bajo  e intermedio.
    + Motor: bajo  e intermedio.
* En el segundo grupo (CPR1 $\approx$ 0.10, CPR2 $\approx$ 0.69):
    + Tama�o: medio-bajo y medio-alto.
    + Precio: medio y alto.
    + Motor: medio.
* En el tercer grupo (CPR1 $\approx$ 1.24, CPR2 $\approx$ -0.31):
    + Tama�o: medio y gran tama�o.
    + Precio: desde intermedio-bajo hasta alto.
    + Motor: tama�o peque�o y medio-alto
* En el cuarto grupo (CPR1 $\approx$ -1.31, CPR2 $\approx$ 2.7). Se situar�an los veh�culos de tama�o muy peque�o, con precio y motor muy grande. Posiblemente veh�culos *"deportivos"*:
    + Tama�o: desde muy peque�o a tama�o medio.
    + Precio: desde alto hasta muy alto.
    + Motor: tama�o desde grande hasta muy grande.



```{r}
library("RcmdrMisc")
# Se asigna el numero de cluster a cada individuo
coches_ln_sin_ventas$num_cluster <- assignCluster(model.matrix(~-1 + CPR1+CPR2, coches_ln_sin_ventas), coches_ln_sin_ventas,cutree(HClust_coches_ln_sin_ventas, k = num_clusters))
```


```{r include=FALSE}
library(ggplot2)
library(factoextra)
fviz_dend(HClust_coches_ln_sin_ventas, k = num_clusters,rect = T, show_labels = T)
```

Gr�ficos de los cuatro clusters

```{r}
ggplot(coches_ln_sin_ventas, aes(x=CPR1, y=CPR2, color=num_cluster)) +
  geom_point() + geom_rug() +
  geom_vline(xintercept = 0) + geom_hline(yintercept = 0) +
  geom_text(aes(label=ifelse((CPR2 > 1 | CPR2 < -1 | CPR1 > 0.8 | CPR1 < -1.1) ,as.character(paste(rownames(coches_ln_sin_ventas), coches_ln_sin_ventas$MARCA)),'')),hjust=0,vjust=0, size=3.3)  +
  scale_colour_brewer(palette = "Set1")
```


### Opci�n 2: Cinco clusters.

```{r include=FALSE}
# k = 5
num_clusters = 5
```

```{r}
(centroides5<-by(model.matrix(~-1 + CPR1+CPR2, coches_ln_sin_ventas), as.factor(cutree(HClust_coches_ln_sin_ventas, k = num_clusters)), colMeans))
```


* En el primer grupo (CPR1 $\approx$ -0.70 , CPR2 $\approx$ -0.74): 
    + Tama�o: bajo  e intermedio.
    + Precio: bajo  e intermedio.
    + Motor: bajo  e intermedio.
* En el segundo grupo (CPR1 $\approx$ 0.10, CPR2 $\approx$ 0.69):
    + Tama�o: medio-bajo y medio-alto.
    + Precio: medio y alto.
    + Motor: medio.
* En el tercer grupo (CPR1 $\approx$ 1.03, CPR2 $\approx$ -0.29):
    + Tama�o: medio y gran tama�o.
    + Precio: desde intermedio-bajo hasta alto.
    + Motor: tama�o peque�o y medio-alto
* En el cuarto grupo (CPR1 $\approx$ -1.32, CPR2 $\approx$ 2.70). Se situar�an los veh�culos de tama�o muy peque�o, con precio y motor muy grande. Posiblemente veh�culos *"deportivos"*:
    + Tama�o: desde muy peque�o a tama�o medio.
    + Precio: desde alto hasta muy alto.
    + Motor: tama�o desde grande hasta muy grande.
* En el quinto grupo (CPR1 $\approx$ 3.14, CPR2 $\approx$ -0.50). Pueden corresponder a algunos veh�culos `TIPO` "truck" con motores peque�os y precio econ�mico en su sector:
    + Tama�o: desde muy grande.
    + Precio: medio.
    + Motor: tama�o medio-bajo.




En este caso los cuatro primeros cluster son pr�cticamente id�nticos, sin embargo aparece un nuevo cluster con s�lo 3 elementos que caracterizan a los *trucks* "econ�micos" que se comercializan en Estados Unidos, lo cual parece interesante. Por tanto pienso que la **segunda opci�n con cinco clusters** es la que mejor puede describir nuestro conjunto.

```{r include=FALSE}
library("RcmdrMisc")
# Se asigna el numero de cluster a cada individuo
coches_ln_sin_ventas$num_cluster <- assignCluster(model.matrix(~-1 + CPR1+CPR2, coches_ln_sin_ventas), coches_ln_sin_ventas,cutree(HClust_coches_ln_sin_ventas, k = num_clusters))
```



Gr�ficos de los 5 clusters

```{r}
ggplot(coches_ln_sin_ventas, aes(x=CPR1, y=CPR2, color=num_cluster)) +
  geom_point() + geom_rug() +
  geom_vline(xintercept = 0) + geom_hline(yintercept = 0) +
  geom_text(aes(label=ifelse((CPR2 > 1 | CPR2 < -1 | CPR1 > 0.8 | CPR1 < -1.1) ,as.character(paste(rownames(coches_ln_sin_ventas), coches_ln_sin_ventas$MARCA)),'')),hjust=0,vjust=0, size=3.3)  +
  scale_colour_brewer(palette = "Set1")
```


Dendogr�ma de los cinco clusters.


```{r}
# num_clusters
library(ggplot2)
library(factoextra)
fviz_dend(HClust_coches_ln_sin_ventas, k = num_clusters,rect = T, show_labels = T)
```


Gr�fico de silueta para los cinco clusters.


```{r}
library(cluster)
k <- num_clusters
cutg <- cutree(HClust_coches_ln_sin_ventas, k=num_clusters)
sil <- silhouette(cutg, coches_ln_sin_ventas.dist)
plot(sil, main="Silhouette plot - Chord - average", cex.names=0.8, col=2:(k+1), nmax=nrow(coches_ln_sin_ventas))
```


Seg�n el gr�fico de siluetas, la coherencia global no es muy alta (0.61) y el grupo 3 es el menos coherente. Los 5 clusters tienen todos los casos bien clasificados (salvo el grupo 3).





```{r include=FALSE}
# Relacionamos cada veh�culo a cada cluster
Tu_C1_5<-subset(coches_ln_sin_ventas,coches_ln_sin_ventas$num_cluster=="1")
Tu_C2_5<-subset(coches_ln_sin_ventas,coches_ln_sin_ventas$num_cluster=="2")
Tu_C3_5<-subset(coches_ln_sin_ventas,coches_ln_sin_ventas$num_cluster=="3")
Tu_C4_5<-subset(coches_ln_sin_ventas,coches_ln_sin_ventas$num_cluster=="4")
Tu_C5_5<-subset(coches_ln_sin_ventas,coches_ln_sin_ventas$num_cluster=="5")
```


```{r include=FALSE}
# Cluster 1
# row.names(Tu_C1_5)
paste(Tu_C1_5$MARCA, Tu_C1_5$MODELO)

## library(RcmdrMisc)
#numSummary(Tu_C1_5[,c("lnREVENTA", "lnPRECIO", "MOTOR", "CV", "PISADA", "ANCHURA", "LONGITUD", "PESO", "CAPACIDAD", "MPG")], #statistics=c("mean", "sd", "quantiles"), quantiles=c(0,.25,.5,.75, 1))

```


```{r include=FALSE}
#Cluster 2
# row.names(Tu_C2_5)
paste(Tu_C2_5$MARCA, Tu_C2_5$MODELO)

```


```{r include=FALSE}
# Cluster 3
# row.names(Tu_C3_5)
paste(Tu_C3_5$MARCA, Tu_C3_5$MODELO)
```


```{r include=FALSE}
#Cluster 4
# row.names(Tu_C4_5)
paste(Tu_C4_5$MARCA, Tu_C4_5$MODELO)

```


```{r include=FALSE}
# Cluster 5
# row.names(Tu_C5_5)
paste(Tu_C5_5$MARCA, Tu_C5_5$MODELO)

```




\pagebreak



# Realizar un An�lisis de Correspondencias.

**Objetivo: Dibujar las marcas de coches en un plano factorial (biplot).**

El an�lisis de correspondencias permite cuantificar de forma objetiva atributos cualitativos. Es una t�cnica descriptiva para representar tablas de contingencia, es decir, tablas donde recogemos las frecuencias de aparici�n de dos o m�s variables en un conjunto de elementos. Constituye el equivalente de componentes principales y coordenadas principales para variables cualitativas.[@pena2013analisis]


## Seleccionar las variables Cualitativas o construirlas.

**Seleccionar las variables (o construirlas) para construir la tabla de contingencia (se puede tambi�n construir una nueva variable cruzando dos variables cualitativas).**

En general, una tabla de contingencia es un conjunto de n�meros positivos dispuestos en una matriz, donde el n�mero en cada casilla representa la frecuencia absoluta observada para esa combinaci�n de las dos variables.[@pena2013analisis]

En nuestro caso a partir de la variable cuantitativa `PRECIO` se va a construir la variable cualitativa `nivelPrecio` que contiene el precio de los veh�culos por escalones de 20.000 USD. Por ejemplo:

| `MARCA`   | `PRECIO` | `nivelPrecio` | Interpretaci�n |
|-----------|---------:|-------------|----------------------|
| Chevrolet | 9235  | **USD.20.K** | Precio menor de 20.000 USD |
| Hyundai   | 9699  | **USD.20.K** | Precio menor de 20.000 USD |
| Ford      | 26935 | **USD.40.K** | Precio entre de 20.000 USD y 40.000 USD |
| Cadillac  | 44475 | **USD.60.K** | Precio entre de 40.000 USD y 60.000 USD |
| Porsche   | 71020 | **USD.80.K** | Precio entre de 60.000 USD y 80.000 USD |
| Mercedes-Benz | 82600 | **USD.100.K** | Precio entre de 80.000 USD y 100.000 USD |

**Nota**:
�Por qu� elegir un escal�n de 20.000 USD en vez de elegir uno de 10.000 USD?
Sencillamente porque despu�s de haber probado ambas opciones, este escal�n da pie a hacer m�s interpretaciones y m�s diversas durante el an�lisis de correspondencias. Porque la calidad que nos da el escal�n de 20.000 USD ser� mejor que el de 10.000 USD para todos los tramos y despu�s el porcentaje de varianza explicada tambi�n ser� mayor y todo esto dar� m�s "*juego*" en el desarrollo de este trabajo.



## Tabla de contingencia y An�lisis de Correspondencias.
**Construir la tabla de contingencia y realizar el An�lisis de Correspondencias.**

### Construir la tabla de contingencia


```{r include=FALSE}

# se crea una copia de la matriz de datos original
coches.corresp <- coches.original

# Se crean una nueva columna segun el nivel de precio del vehiculo con xxxxx niveles
escalon_precio = 20000

# VENTAS vs PRECIO
# Se suman todas las ventas hechas en cada modelo una marca para cada escalon de precios

miles = trunc(escalon_precio/1000)



# Se halla el numero total de nuevas columnas
numNuevasColumnas =ceiling(1 + (max(coches.original$PRECIO) - min(coches.original$PRECIO) )/escalon_precio)

# Se crea una lista con el nombre de las nuevas columnas
lista_columnas <- vector(mode="character", length=0)

# Se crean las columnas
for (nFila in 1:numNuevasColumnas) {

coches.corresp$nombreFalsoDecolumna <- 0
nombreCol = paste("USD.", miles * nFila , ".K", sep = "")
colnames(coches.corresp)[colnames(coches.corresp)=="nombreFalsoDecolumna"] <- nombreCol

lista_columnas <- c(lista_columnas, nombreCol) 
}


# Se meten las ventas en cada columna
for (nFila in 1:nrow(coches.corresp)) {

numero_col_destino = as.integer(trunc(1 + coches.corresp[nFila,"PRECIO"]/escalon_precio))
coches.corresp[nFila, lista_columnas[numero_col_destino]] <- coches.corresp[nFila,"VENTAS"]
}

## # Se suma el contenido de todas las columnas que coinciden en MARCA y nivelPrecio
## tab.suma <- aggregate(coches.corresp[,sapply(coches.corresp,is.numeric)], c(coches.corresp["MARCA"], coches.corresp["nivelPrecio"]),sum)
# Se suma el contenido de todas las columnas que coinciden en MARCA
coches.corresp <- aggregate(coches.corresp[,sapply(coches.corresp,is.numeric)],coches.corresp["MARCA"],sum)

### rm(tab.cont.coches)

# Se crea la tabla de contingencia
# Se crea la tabla con la primera columna: MARCA
tab.cont.coches <- coches.corresp["MARCA"]
# Se a�aden el resto de columnas
for (estaCol in 1:numNuevasColumnas){
  tab.cont.coches = append(tab.cont.coches, coches.corresp[lista_columnas[estaCol]])
}


tab.cont.coches <- as.data.frame(tab.cont.coches)
# Se usa la MARCA como nombre de fila
rownames(tab.cont.coches) <- tab.cont.coches$MARCA
tab.cont.coches$MARCA <- NULL

## write.table(tab.cont.coches, "mydata.txt", sep="\t")

tab.cont.coches
```




Se crea una tabla de contingencia que contiene el n�mero de `VENTAS` por `MARCA` en cada `nivelPrecio`. 

As� por ejemplo: la `MARCA` `Acura` no vende ning�n veh�culo de menos de 20.000 USD, vende 70.417 veh�culos cuyo precio est� entre 20.000 USD y 40.000 USD, vende 8.588  veh�culos cuyo precio est� entre 40.000 USD y 60.000 USD y no vende ning�n veh�culo por encima de 60.000 USD.


| `MARCA`  | USD.20.K | USD.40.K | USD.60.K | USD.80.K  | USD.100.K |
|----------|---------:|---------:|---------:|----------:|----------:|
| Acura | 0 | 70417 | 8588 | 0 | 0 |
| Audi | 0 | 39177 | 0 | 1380 | 0 |
| BMW | 0 | 26758 | 0 | 0 | 0 |
| Buick | 0 | 242019 | 0 | 0 | 0 |
| Cadillac | 0 | 81450 | 15943 | 0 | 0 |
| Chevrolet | 402021 | 26402 | 17947 | 0 | 0 |
| Chrysler | 40160 | 77385 | 0 | 0 | 0 |
| Dodge | 631929 | 175981 | 0 | 916 | 0 |
| Ford | 571760 | 1275205 | 0 | 0 | 0 |
| Honda | 430587 | 162087 | 0 | 0 | 0 |
| Hyundai | 137326 | 0 | 0 | 0 | 0 |
| Infiniti | 0 | 23713 | 0 | 0 | 0 |
| Jeep | 55557 | 237596 | 0 | 0 | 0 |
| Lexus | 0 | 36770 | 6375 | 0 | 0 |
| Lincoln | 0 | 13798 | 48911 | 0 | 0 |
| Mercedes-Benz | 0 | 18392 | 27602 | 16774 | 3311 |
| Mercury       | 108836 | 129163 | 0 | 0 | 0 |
| Mitsubishi    | 124389 | 56506 | 0 | 0 | 0 |
| Nissan        | 42643 | 237829 | 0 | 0 | 0 |
| Oldsmobile    | 1112 | 59068 | 0 | 0 | 0 |
| Plymouth      | 62129 | 0 | 0 | 0 | 0 |
| Pontiac       | 131097 | 199865 | 0 | 0 | 0 |
| Porsche       | 0 | 0 | 8982 | 3146 | 0 |
| Saturn | 110389 | 0 | 0 | 0 | 0 |
| Toyota | 532991 | 132260 | 9835 | 0 | 0 |
| Volkswagen | 108647 | 51102 | 0 | 0 | 0 |



### Realizar el An�lisis de Correspondencias.

```{r}
# # row percentages
# prop.table(tab.cont.coches, 1)
# # column percentages
# prop.table(tab.cont.coches, 2) 

library(ca)
# Se realiza el analisis de correspondencias.
tab.cont.coches.ac <- ca(tab.cont.coches)
# tab.cont.coches.ac
summary(tab.cont.coches.ac)
```


## Interpretar los factores obtenidos.

Como era de esperar para una matriz de $26 \times 5$ se han obtenido 4 valores propios con los que se puede explicar el 100% de la varianza. Con los dos primeros factores se explica el 84.3 % de la varianza.

Se utilizan las contribuciones absolutas (`ctr`) para interpretar las dimensiones (o factores). Respecto a las columnas, que contiene los tramos de precios de los veh�culos:

* En la primera dimensi�n las contribuciones m�s altas corresponden a:
    + "< USD 60 K", con contribuci�n negativa. Su contribuci�n relativa `cor` es muy alta (= 869); o sea que est� bien representada.
    + "< USD 80 K", tambi�n con contribuci�n negativa.  Su contribuci�n relativa `cor` es alta (= 594); o sea que tambi�n  est� bien representada.
    
  Por tanto la primera dimensi�n representa a los veh�culos vendidos m�s caros (entre 40.000 y 80.000 USD) en el lado negativo y los m�s baratos a penas est�n representados.
  
* En la segunda dimensi�n las contribuciones absolutas m�s altas corresponden a:
    + "< 40 K USD", con contribuci�n positiva. Su contribuci�n relativa `cor` es muy alta (= 999); o sea que  est� bien representada.
    + "< 20 K USD", con contribuci�n negativa.  Su contribuci�n relativa `cor` es alta; o sea que tambi�n  est� bien representada.
    
  O sea que en la segunda dimensi�n est�n representados los veh�culos vendidos m�s baratos situ�ndose los de menor precio (< 20.000 USD) en la parte negativa y los que oscilan entre 20.000 USD y 40.000 USD est�n ubicados en la parte positiva de este factor.

* La modalidad `USD.100.K` (veh�culos entre 80.000 y 100.000 USD) es la que peor est� representada, dado que su par�metro de calidad (`qlt`) es el menor de todas

* Respecto a las marcas, todas est�n bien representadas en el plano factorial dado que su calidad (`qlt`) es alta en todos los casos.


En el �ltimo apartado de este ejercicio se relacionan el nivel de ventas de cada marca con los precios de los veh�culos.


## Representar las marcas de coches en el primer plano factorial.

Se muestra el gr�fico sim�trico en el que las filas y columnas se representan en coordenadas principales, es decir, se realiza el calibrado de manera que la inercia es igual a la inercia principal (valor propio o cuadrado del valor singular). En este gr�fico se aprecia el nivel de `VENTAS` de cada `MARCA` asociada a un rango de `PRECIO`.

```{r}
# symmetric map
## plot(tab.cont.coches.ac) # Por defecto map = "symmetric"
## plot(tab.cont.coches.ac, map = "rowprincipal", main="rowprincipal")
## plot(tab.cont.coches.ac, map = "colprincipal", main="colprincipal")

## # asymmetric map 
## plot(tab.cont.coches.ac, mass = TRUE, contrib = "absolute", map = "rowgreen", arrows = c(FALSE, TRUE))

fviz_ca_biplot(tab.cont.coches.ac, repel=TRUE, map = "symmetric")

```


```{r}

```


## Explicar las conclusiones.

Se pueden realizar las siguientes asociaciones a cada `MARCA`:

* `Porche` se caracteriza por vender veh�culos caros entre 60.000 y 80.000 USD. Dado que se sit�a muy cerca del punto `USD.80.K` y dicho punto presenta una alta calidad (`qlt`= 60,4%).
* `Lincoln` se caracteriza por vender veh�culos entre 40.000 y 60.000 USD. Dado que se sit�a muy cerca del punto `USD.60.K` y dicho punto tambi�n presenta una alta calidad.
* `Mercedes-Benz` se sit�a entre el punto 'USD.60.K y 'USD.80.K`. Ambos puntos presentan alta calidad. O sea que esta marca se caracteriza por vender entre los rangos de 40.000 y 80.000 USD.

En el otro extremo del plano factorial, 

* Los veh�culos de precio medio-alto (`USD.40.K`) se encuentran en la parte superior. Por ejemplo: `BMW`, `Buick`, `Audi`, `Oldsmobile`, etc.
* Los veh�culos de precio bajo (`USD.20.K`)se encuentran en la parte inferior. Por ejemplo: `Saturn`, `Hyundai`, `Plymouth`, `Chevrolet`, `Toyota`, etc. Seg�n se asciende por este eje, el nivel de `PRECIOS` asociados a las `VENTAS` de esa `MARCA` tambi�n aumenta.

Con este an�lisis de correspondencias se llegan a conclusiones adicionales a las que se han obtenido en el an�lisis factorial. Aqu� se distribuyen las marcas seg�n los precios con m�s detalle.


# An�lisis de Correspondencias M�ltiples y An�lisis Cluster.

**Realizar un An�lisis de Correspondencias y, a continuaci�n, un An�lisis Cluster**

**Objetivo: Similar al An�lisis Factorial, pero incluyendo todas las variables (o la mayor�a).**


## Dividir las variables cuantitativas en cuartiles.

**Dividir las variables cuantitativas en cuartiles (pod�is elegir otra divisi�n si os parece m�s adecuada).**

```{r include=FALSE}

library("gtools")
# Funcion para dividir en cuartiles
cutN <- function(X , n = 4){
     cut(X ,
         include.lowest = TRUE ,
         breaks = quantile(X , probs = (0:n)/n ,na.rm = TRUE),
         labels = FALSE)}
```


```{r}
df_coches <- coches.original
# names(df_coches)

num_divisiones = 4 # 7, 6, 5, 4

df_coches$VENTAS.q<-factor(cutN(df_coches$VENTAS, n = num_divisiones))
df_coches$REVENTA.q<-factor(cutN(df_coches$REVENTA, n = num_divisiones))
df_coches$PRECIO.q<-factor(cutN(df_coches$PRECIO, n = num_divisiones))
df_coches$MOTOR.q<-factor(cutN(df_coches$MOTOR, n = num_divisiones))
df_coches$CV.q<-factor(cutN(df_coches$CV), n = num_divisiones)
df_coches$PISADA.q<-factor(cutN(df_coches$PISADA), n = num_divisiones)
df_coches$ANCHURA.q<-factor(cutN(df_coches$ANCHURA), n = num_divisiones)
df_coches$LONGITUD.q<-factor(cutN(df_coches$LONGITUD), n = num_divisiones)
df_coches$PESO.q<-factor(cutN(df_coches$PESO), n = num_divisiones)
df_coches$CAPACIDAD.q<-factor(cutN(df_coches$CAPACIDAD), n = num_divisiones)
df_coches$MPG.q<-factor(cutN(df_coches$MPG, n = num_divisiones))
```


Las variables cualitativas divididas en cuartiles tienen la siguiente distribuci�n:

```{r}
# Se crea una nueva matriz con las siguientes 14 columnas:
# "MARCA", "MODELO", "TIPO", "VENTAS.q", "REVENTA.q", "PRECIO.q", "MOTOR.q", "CV.q"       
# "PISADA.q", "ANCHURA.q", "LONGITUD.q", "PESO.q", "CAPACIDAD.q", "MPG.q"
## df_coches.ac<-df_coches[,c(15:25)]
###df_coches.ac<-df_coches[,c(5, 15:25)]

# Se incluyen: MARCA, MODELO, TIPO y  los cuartiles
df_coches.ac<-df_coches[,c(1:2, 5, 15:25)]
# Se declaran MARCA, MODELO y TIPO como un factor
df_coches.ac$MARCA <- factor(df_coches.ac$MARCA)
df_coches.ac$MODELO <- factor(df_coches.ac$MODELO)
df_coches.ac$TIPO <- factor(df_coches.ac$TIPO)

# Se reordena la matriz. Primero las variables activas, despues las ilustrativas y al final MARCA y MODELO
#  - Primero las variables activas: VENTAS, REVENTA, PRECIO. Que estan en la posicion (4:6)
#  - Despues las variables ilustrativas: en la posicion 3 esta el TIPO y en (7:14) los cuartiles
#  - Al final MARCA y MODELO: en la posicion (1:2)
# De tal forma que ahora el orden queda: 
#  - Variables activas (1:3). VENTAS, REVENTA, PRECIO
#  - Variables ilustrativas: (4:12). TIPO y cuartiles
df_coches.ac <- df_coches.ac[c(4:6,3,7:14,1:2)]


## names(df_coches.ac)

# Se muestran los datos de las 11 variables cuanitativas dividudas en cuartiles.
summary(df_coches.ac[c(1:3, 5:12)])
```
```{r include=FALSE}
0.02 * 120
```

El criterio para definir "*efectivo peque�o*" es: $2 \% n =   0,02*n = 0,02*120 = 2,4$. O sea, un m�nimo de 3 respuestas. Por tanto nuestros datos cumplen adecuadamente el criterio.


## Variables activas y variables ilustrativas.

**Definir las variables activas y las ilustrativas que os parezcan m�s adecuadas para realizar un An�lisis de Correspondencias M�ltiples.**

Se tomar�n como **variables activas**, las m�s subjetivas, las que pueden variar por formar parte la decisi�n de las personas compradoras: el n�mero de `VENTAS.q` y el precio de `REVENTA.q`. As� mismo se puede considerar `PRECIO.q` como una variable subjetiva, porque el valor de salida del veh�culo se puede variar a lo largo de la campa�a de ventas. Se tomar�n como **variables ilustrativas** los datos objetivos de los veh�culos. Por otro lado, la `MARCA`y el `MODELO` no se incluyen en el an�lisis porque se usar�n como el nombre de cada fila o �ndice.

* Variables activas: `VENTAS.q`, `REVENTA.q`, `PRECIO.q`.
* Variables ilustrativas: `TIPO`, `MOTOR.q`, `CV.q`, `PISADA.q`, `ANCHURA.q`, `LONGITUD.q`, `PESO.q`, `CAPACIDAD.q`, `MPG.q`.


## An�lisis de Correspondencias M�ltiples.

**Realizar el An�lisis de Correspondencias M�ltiples y valorar razonadamente las posibles mejoras de este an�lisis.**

Se realiza un an�lisis de correspondencias m�ltiples considerando inicialmente solo las variables activas.

```{r}
library(ca)
#  VENTAS.q, REVENTA.q, PRECIO.q
ACM_coches <- mjca(df_coches.ac[,1:3])
#  VENTAS.q, REVENTA.q
## ACM_coches <- mjca(df_coches.ac[,c(1:2)])
# Todas las variables menos MARCA y MODELOS
## ACM_coches <- mjca(df_coches.ac[,1:12])

summary(ACM_coches)
```
El an�lisis de correspondencias m�ltiples (ACM) ha generado una matriz de Burt, que es la matriz compuesta por todas las tablas resultantes de cruzar todas las variables dos a dos [@greenacre2008practica].

En todas las tablas, las inercias de todos los ejes principales ser�n iguales a 1. Y, por tanto, la inercia total de la tabla $Z_q$ ser� igual a su dimensionalidad, es decir, igual a $J_q -1$. La inercia de $Z$ ser� la media de la inercias de las tablas que la componen:

$$inercia (Z) = \frac{1}{Q} \sum _q (Z_q) = \frac{1}{Q} \sum _q (J_q - 1) = \frac {J-Q}{Q}$$

Dado que $J-Q$ es la dimensionalidad de $Z$, la inercia media por dimension ser� $1/Q$. Se utiliza este valor $1/Q$ como umbral para decidir para qu� ejes es interesante el ACM [@greenacre2008practica]. Similar al umbral de 1 de los valores propios en el an�lisis de componentes principales.

En nuestro caso $Q = 3$, por tanto $\frac{1}{Q} = \frac{1}{3} = 0.333$ y nuestra inercia principal vale 0.591 que est� por encima de $1/3$.

Con dos factores se explica el 85.3 % de la varianza. Las contribuciones m�s altas para cada factor son:

* Factor 1: `REVENTA.q:4` (positivo), `PRECIO.q:4` (positivo),  `VENTAS.q:1` (positivo)
* Factor 2: `PRECIO.q:3` (positivo), `REVENTA.q:3` (positivo),  `REVENTA.q:1` (negativo), `PRECIO.q:1`  (negativo)

Este an�lisis puede mostrar el comportamiento para los precios, el n�mero de ventas y el precio de reventas m�s altos y m�s bajos. No as� las para las intermedias.

M�s adelante se explicar�n con detalle los factores.


```{r}
# plot(ACM_coches)
# plot(ACM_coches, mass = c(FALSE, TRUE))
# plot(ACM_coches, mass = c(FALSE, TRUE), contrib="absolute")

library(ggplot2)
library(ggrepel)

# Se calculan las distancias al cuadrado
dt_act<-cbind.data.frame(ACM_coches$colcoord[,1],
                     ACM_coches$colcoord[,2],
                     ACM_coches$levelnames)  ### paste(rownames(df_coches.ac), df_coches.ac$MARCA)
# Se crea el objeto con el grafico de las variables activas
var_act_grob <- ggplot(dt_act, aes(x=dt_act[,1], y=dt_act[,2])) +
  geom_point() + geom_rug() +
  # xlim(-5.5, 2) + ylim(-2.5,10)+
  geom_text_repel(aes(label=dt_act[,3]),hjust=0, vjust=0,size=3.3) +
  scale_colour_brewer(palette = "Set1") +
  geom_vline(xintercept = 0) + geom_hline(yintercept = 0) + 
  ggtitle("Variables activas en el plano factorial") +
  xlab(paste("Dimension 1 (", round(100 * ACM_coches$inertia.e[1], digits = 1) , "% )", sep = "")) +
  ylab(paste("Dimension 2 (", round(100 * ACM_coches$inertia.e[2], digits = 1) , "% )", sep = ""))

var_act_grob
```

### Posibles mejoras para el ACM.

#### Opci�n 1.
Si se aumenta el n�mero de dimensiones a 3 se mejora ligeramente el resultado final tal y como se puede ver en la siguiente tabla resumen. Al aparecer una nueva dimensi�n, las dos primeras no var�an, como es de esperar. Pero la tercera, en este caso, muestra altas contribuciones para `PRECIO.q:3` (positivo), `REVENTA.q:3`(positivo), `REVENTA.q:1` (negativo) y `PRECIO.q:1` (negativo). Estas variables no estaban bien representadas en ninguna de las otras dos dimensiones.

Pero tambi�n hay que hacer notar que al a�adir esta nueva dimensi�n el porcentaje de varianza explicada mejora muy poco, solo un 3.8 %. por tanto no merece la pena a�adir esta dimensi�n complicando las interpretaciones.

```{r include=FALSE}
89.1 - 85.3
```



```{r}
ACM_coches_nd3 <- mjca(df_coches.ac[,1:3])
summary(ACM_coches_nd3)
```

##### Opci�n 2.

Otra posible mejora ser�a tratar todas las variables como ilustrativas (ver tabla con el resumen m�s adelante). Esto nos dar�a mayor precisi�n en la interpretaci�n de los factores. Puesto que para este caso las contribuciones m�s altas para cada factor son :

* Factor 1: `MOTOR.q:1` (positivo), `PESO.q:1` (positivo), `CAPACIDAD.q:1` (positivo), `CV.q:1` (positivo), `MPG.q:4` (positivo), `PRECIO.q:1` (positivo), `ANCHURA.q:1` (positivo), `REVENTA.q:1` (positivo), `LONGITUD.q:1` (positivo), `PESO.q:4` (negativo), `MPG.q:1` (negativo), `PISADA.q:1` (positivo), `PISADA.q:4` (negativo) y `MOTOR.q:4` (negativo). Que claramente situa en el eje positivo a los veh�culos m�s peque�os, con menos motor y con mayor autonom�a de combustible. Y pone en el negativo a los m�s grandes, con m�s peso y menor autonom�a (menos `MPG`).
* Factor 2: `MOTOR.q:4` (negativo), `PESO.q:4` (negativo), `CV.q:2` (positivo), `CAPACIDAD.q:2` (positivo), `MOTOR.q:2` (positivo), `MPG.q:1` (negativo), `PRECIO.q:2` (positivo), `PESO.q:2` (positivo), `MPG.q:3` (positivo), `MPG.q:2` (positivo), `PRECIO.q:4` (negativo),  `REVENTA.q:4` (negativo). Que situa en el lado negativo a los veh�culos con mayor motor, menor auton�m�a (`MPG.q:1`), m�s precio y m�s reventa. Y posiciona en el lado positivo a los que tienen potencia, capacidad de tanque, motor y consumo medio-bajo.

Este an�lisis es m�s certero a la hora clasificar las caracteristicas principales de los veh�culos pero no nos muestra ninguna posible relaci�n entre las decisiones de las personas compradoras y las caracter�sticas de los veh�culos.


```{r}
# library(ca)

# Todas las variables menos MARCA y MODELOS
ACM_coches.todas <- mjca(df_coches.ac[,1:12])

summary(ACM_coches.todas)
```

## An�lisis de Correspondencias M�ltiples definitivo.

**Realizar el An�lisis de Correspondencias M�ltiples definitivo (desde tu punto de vista)**

Por lo expuesto anteriormente, el grupo de variables elegido es el inicial dado que el an�lisis puede resultar m�s interesante. 

Se crea un modelo con la siguiente distribuci�n de variables:

* Variables activas: `VENTAS.q`, `REVENTA.q`, `PRECIO.q`.
* Variables ilustrativas: `TIPO`, `MOTOR.q`, `CV.q`, `PISADA.q`, `ANCHURA.q`, `LONGITUD.q`, `PESO.q`, `CAPACIDAD.q`, `MPG.q`.


```{r echo = TRUE}
# Hay que meter la matriz bien ordenada:
#  - Primero, en las posiciones (1:3), las variables activas: VENTAS.q, REVENTA.q, PRECIO.q.
#  - Despues, en las posiciones (4:12) las variables ilustrativas: 
#      en la posicion 4 esta el TIPO y en (5:12) los cuartiles
# No se incluyen las variables MARCA y MODELO (13:14). Se usar�n como etiquetas.
ACM_coches.ac.ilu <- mjca(df_coches.ac[1:12], lambda="Burt", supcol=c(4:12), subsetcol="1:3", nd=2)
summary(ACM_coches.ac.ilu)
```

```{r}
# Variables ilustrativas en el plano factorial

#### strAux = colnames(df_coches.ac)[5:12]
#### lista_var_ilustrativas <- vector("character")
#### for (estaCol in strAux)
####     for (n in 1:4)
####       lista_var_ilustrativas <- append(lista_var_ilustrativas, paste(estaCol,":", n, sep=""))
#### 
### lista_var_ilustrativas

### lista_var_ilustrativas <- ACM_coches.ac.ilu$levelnames[c(15:46)]

### ACM_coches.ac.ilu$levelnames[c(15:46)]
### ACM_coches.ac.ilu$levelnames[ACM_coches.ac.ilu$colsup]

library(ggplot2)
library(ggrepel)
##  geom_text(aes(label=dt[,3]),hjust=0, vjust=0,size=3.3) +

dt<-cbind.data.frame(ACM_coches.ac.ilu$colcoord[ACM_coches.ac.ilu$colsup,1],
                     ACM_coches.ac.ilu$colcoord[ACM_coches.ac.ilu$colsup,2],
                     ACM_coches.ac.ilu$levelnames[ACM_coches.ac.ilu$colsup])  ### paste(rownames(df_coches.ac), df_coches.ac$MARCA)
ggplot(dt, aes(x=dt[,1], y=dt[,2])) +
  geom_point() + geom_rug() +
  # xlim(-5.5, 2) + ylim(-2.5,10)+
  geom_text_repel(aes(label=dt[,3]),hjust=0, vjust=0,size=3.3) +
  scale_colour_brewer(palette = "Set1") +
  geom_vline(xintercept = 0) + geom_hline(yintercept = 0) + 
  ggtitle("Variables ilustrativas en el plano factorial") +
  xlab(paste("Dimension 1 (", round(100 * ACM_coches$inertia.e[1], digits = 1) , "% )", sep = "")) +
  ylab(paste("Dimension 2 (", round(100 * ACM_coches$inertia.e[2], digits = 1) , "% )", sep = ""))

```



## Interpretar los factores obtenidos.

Tal y como se dijo anteriormente con dos factores se explica el 85.3 % de la varianza. 

Las contribuciones absolutas m�s altas (`ctr`) para cada factor son:

* Factor 1: `REVENTA.q:4` (positivo), `PRECIO.q:4` (positivo), `VENTAS.q:1` (positivo). Todas estas variables muestran una alta calidad (`qlt`). Adem�s est�n bien representadas porque sus contribuciones relativas (`cor`) son altas para cada una de ellas. Por tanto este factor representa a los veh�culos m�s caros que son los que menos ventas tienen, los de precios m�s altos y que en reventa tambi�n alcanzan un alto precio. Se situan en el lado positivo del eje.
* Factor 2: `PRECIO.q:3` (positivo), `REVENTA.q:3` (positivo), `REVENTA.q:1` (negativo), `PRECIO.q:1`  (negativo). Tambi�n presentan una alta calidad. Adem�s est�n bien representadas porque sus contribuciones relativas (`cor`) son altas para cada una de ellas. Esta dimensi�n representa a los veh�culos de precio medio-alto (lado positivo del eje) y a los de bajo precio (lado negativo del eje).

Este an�lisis puede mostrar el comportamiento para los precios, ventas y reventas m�s altas y m�s bajas. No as� las para las intermedias.



Para interpretar las variables ilustrativas en el plano factorial tenemos que tener en cuenta la calidad (`qlt`) de la representaci�n de cada variable ilustrativa del veh�culo. No podremos sacar conclusiones sobre las caracter�sticas de los veh�culos con mala calidad.

* Datos con `qlt` baja (< 50%) de los que no se pueden sacar conclusiones: `CV.q:2`, `MPG.q:3`, `LONGITUD.q:4`, `ANCHURA.q:3`, `PESO.q:2`, `MOTOR.q:2`, `MPG.q:2`, `CAPACIDAD.q:2`, `ANCHURA.q:2`

A continuaci�n se muestra el plano factorial sin las etiquetas de esas variables y se continua con el an�lisis.


```{r}
# Se elimina la etiqueta de esa variables ilustrativas

# Se hace una copia de las lista de variables ilustrativas
lista_var_ilustrativas_qlt_alta <- ACM_coches.ac.ilu$levelnames[ACM_coches.ac.ilu$colsup]


### lista_var_ilustrativas_qlt_alta
### which(ACM_coches.ac.ilu$colcor[ACM_coches.ac.ilu$colsup,1] + ACM_coches.ac.ilu$colcor[ACM_coches.ac.ilu$colsup,2] < 0.5)

qlt_baja = 0.5

# qlt = k1(cor) + k2(cor)
# qlt = ACM_coches.ac.ilu$colcor[ACM_coches.ac.ilu$colsup,1] + ACM_coches.ac.ilu$colcor[ACM_coches.ac.ilu$colsup,2]

# Se eliminan solo las etiquetas de las variables ilustrativas con qlt baja
lista_var_ilustrativas_qlt_alta[which(ACM_coches.ac.ilu$colcor[ACM_coches.ac.ilu$colsup,1] + ACM_coches.ac.ilu$colcor[ACM_coches.ac.ilu$colsup,2] < qlt_baja)] = ""

## lista_var_ilustrativas_qlt_alta

dt<-cbind.data.frame(ACM_coches.ac.ilu$colcoord[ACM_coches.ac.ilu$colsup,1],
                     ACM_coches.ac.ilu$colcoord[ACM_coches.ac.ilu$colsup,2],
                     lista_var_ilustrativas_qlt_alta)  ### paste(rownames(df_coches.ac), df_coches.ac$MARCA)
ggplot(dt, aes(x=dt[,1], y=dt[,2])) +
  geom_point() + geom_rug() +
  # xlim(-5.5, 2) + ylim(-2.5,10)+
  geom_text_repel(aes(label=dt[,3]),hjust=0, vjust=0,size=3.3) +
  scale_colour_brewer(palette = "Set1") +
  geom_vline(xintercept = 0) + geom_hline(yintercept = 0) + 
  ggtitle(paste("Variables ilustrativas con qlt >", qlt_baja * 100, "% en el plano factorial")) +
  xlab(paste("Dimension 1 (", round(100 * ACM_coches$inertia.e[1], digits = 1) , "% )", sep = "")) +
  ylab(paste("Dimension 2 (", round(100 * ACM_coches$inertia.e[2], digits = 1) , "% )", sep = ""))

```


Las variables con `qlt` mayor del 60 % pueden ser aptas para sacar conclusiones: 

* En el cuadrante superior derecho:
    + Lado derecho. N�mero de `VENTAS` muy bajas (`VENTAS.q:1`), `PRECIO` de venta y `REVENTA` muy altos (`PRECIO.q:4` y `REVENTA.q:4`): Los veh�culos asociados a esta descrici�n son los que tienen `CV.q:4` (muy alto) y/o `MOTOR.q:4` (muy grande)  y/o `PESO.q:4` (mucho peso)  y/o `MPG.q:1` (muy poca autonom�a) y/o `CAPACIDAD.q:3`, `CAPACIDAD.q:4` (gran capacidad del tanque de combustible).
    + Lado superior del cuadrante. `PRECIO.q:3` (bastante alto), `REVENTA.q:3` (precio de reventa bastante alto). Este semieje representa a los veh�culos de precio medio-alto tanto en venta como en reventa:  Los veh�culos asociados a esta descrici�n son los que tienen `CV.q:3`, `MOTOR.q:3`, `LONGITUD.q:3`, `PESO.q:4`, `CAPACIDAD.q:3`, `ANCHURA.q:4`, `CAPACIDAD.q:4`, `MPG.q:1` y `PISADA.q:3`. O sea, veh�culos con un valor medio-alto en su `CV` o `MOTOR` o `LONGITUD` o `CAPACIDAD` de dep�sito o `ANCHURA` de neum�tico y  tambi�n los de gran `CAPACIDAD` de dep�sito y  muy poca autonom�a de combustible (`MPG` muy baja).
    
* En el cuadrante inferior izquierdo: 
    +  Lado izquierdo. Precio de `REVENTA` muy bajo (`REVENTA.q.1`) y `PRECIO` de venta muy bajo  (`PRECIO.q.1`). Este semieje representa a los veh�culos m�s econ�micos.  Los veh�culos asociados a esta descrici�n son los que tienen `CV.q:1`, `CAPACIDAD.q:1`, `MPG.q:4`, `PESO.q:1`, `MOTOR.q:1`, `PISADA.q:2` y `LONGITUD.q:1`. Los veh�culos con estas caracter�sticas est�n asociados a `CV.q:1`, `CAPACIDAD.q:1`, `MPG.q:4`, `PESO.q:1`, `MOTOR.q:1` y `ANCHURA.q:1`. O sea muy baja `CV`, `CAPACIDAD` de combustible, `PESO`, `MOTOR` y `ANCHURA` y adem�s muy alta autonom�a (`MPG` alta).

* Adem�s tambi�n se puede decir que los veh�culos `TIPO` `Truck` tienen un `PRECIO` m�s alto que los veh�culos `TIPO` `Automobile`. Dado que el segundo factor representa mayor precio cuanto mayor es su valor.



## Coches en el planos factorial.

**Representar los coches en los diferentes planos factoriales (o en el primero, si es suficiente)**

No hay m�s dos factores. O sea, un solo plano factorial. Por tanto se representan los coches en ese plano. Se evita poner la `MARCA` de cada veh�culo para no saturar el gr�fico.

Dado que las variables se han agrupado en cuartiles, ahora los veh�culos con caracter�sticas similares aparecen aglutinadas en un solo punto.

```{r}
# Se hace un ACP con todas las variables con el lambda ="indicator"
ACM_coches.ac.ilu.ind <- mjca(df_coches.ac[,1:12], lambda ="indicator", supcol=c(4:12), subsetcol="1:3",nd=2)
#summary(ACM_coches.ac.ilu.ind) #no son ajustado no se pueden interpretar
```

```{r}
library(ggplot2)
library(ggrepel)
##  geom_text(aes(label=dt[,3]),hjust=0, vjust=0,size=3.3) +

dt<-cbind.data.frame(ACM_coches.ac.ilu.ind$rowcoord[,1],
                     ACM_coches.ac.ilu.ind$rowcoord[,2],
                     rownames(df_coches.ac))  ### paste(rownames(df_coches.ac), df_coches.ac$MARCA)
ggplot(dt, aes(x=dt[,1], y=dt[,2])) +
  geom_point() + geom_rug() +
  # xlim(-5.5, 2) + ylim(-2.5,10)+
  geom_text_repel(aes(label=dt[,3]),hjust=0, vjust=0,size=3.3) +
  scale_colour_brewer(palette = "Set1") +
  geom_vline(xintercept = 0) + geom_hline(yintercept = 0)  +
  xlab("CPR1") +
  ylab("CPR2")

```



## Explicar las conclusiones y compararlas con los resultados del An�lisis Factorial.



FALTA    FALTA    FALTA    FALTA    FALTA    FALTA    FALTA    FALTA    
FALTA    FALTA    FALTA    FALTA    FALTA    FALTA    FALTA    FALTA    
FALTA    FALTA    FALTA    FALTA    FALTA    FALTA    FALTA    FALTA    
FALTA    FALTA    FALTA    FALTA    FALTA    FALTA    FALTA    FALTA    
FALTA    FALTA    FALTA    FALTA    FALTA    FALTA    FALTA    FALTA    
FALTA    FALTA    FALTA    FALTA    FALTA    FALTA    FALTA    FALTA    


## Cluster jer�rquico.

**Realizar un cluster jer�rquico, tomando como variables las coordenadas factoriales obtenidas en el ACM anterior.**

Antes de realizar un *an�lisis cluster* deben tomarse tres decisiones:
1. Selecci�n de las variables relevantes para identificar a los grupos. Esto ya se ha hecho en los apartados previos.
2. Elecci�n de la medida de proximidad entre los individuos. En este caso se elegir� la distancia eucl�dea al cuadrado, que es la requieren los m�todos *centroide*, de la *mediana* y  de *Ward*.
3. Elecci�n del criterio para agrupar individuos. En este caso a�n est� por decidir



## Decidir el n�mero de cluster y caracterizar estos cluster.

A�n est� por decicir el criterio para agrupar individuos y  el n�mero de clusters.

Se muestra como ejemplo una divisi�n cluster por el m�todo *Ward* o tambi�n conocido como *"p�rdida de inercia m�nima"*

```{r include=FALSE}

ACM_coches.def.ilu <- ACM_coches.ac.ilu

df_coches.def.ilu.ac <- df_coches.ac[,1:12]

### library(calibrate)
### 
### plot(ACM_coches.def.ilu$colcoord[c(13:46),1],ACM_coches.def.ilu$colcoord[c(13:46),2],xlim=c(-0.9, 1),ylim=c(-1.1, 1), xlab = "Variables ilustrativas en el plano factorial", ylab = "")
### textxy(ACM_coches.def.ilu$colcoord[c(13:46),1],ACM_coches.def.ilu$colcoord[c(13:46),2],labs= c(ACM_coches.def.ilu$levelnames[13:46]),cex=0.75)
### abline(h=0,v=0, lty=3)
### 
```


```{r}
df_coches.def.ilu.ac$C1 <- ACM_coches.ac.ilu.ind$rowcoord[,1]
df_coches.def.ilu.ac$C2 <- ACM_coches.ac.ilu.ind$rowcoord[,2]
# Distancia euclidea al cuadrado entre los elementos de la matriz
dist_coches<-dist(model.matrix(~-1 + C1+C2, df_coches.def.ilu.ac))^2
HClust.2 <- hclust(dist_coches , method= "ward.D")

plot(HClust.2, main= "Cluster Dendrograma para coches", xlab= 
  "Veh�culos estudiados", sub="Method=ward.d; Distancia=euclidea al cuadrado")
```

Parece que el n�mero de cluster adecuado es 3. Pero al igual que en el an�lisis de correspondencias simples, aqu�  tambi�n se puede buscar el n�mero �ptimo de clusters (`k`) usando la funci�n `daisy` (Dissimilarity Matrix Calculation) que computa las desigualdades en las distancias entre parejas. En este caso se ha usado la m�trica euclidean para ese c�lculo. Adem�s se van usando difentes m�todos de agrupaci�n.

A continuaci�n se muestra un resumen de todo ello:


```{r}

## # Lista de metodos a testear
## lista_mtdos = c("ward", "complete", "single", "average", "mcquitty", "median", "centroid")
# data frame auxialiar
kt <- data.frame(k=1:nrow(df_coches.def.ilu.ac), r=0)

# Bucle para testear todos los metodos con la distancia elegida
for (mtdo_usado in lista_mtdos)
{
  HClust_aux <- hclust(dist_coches , method= mtdo_usado)

  # plot(HClust_aux, main= "Cluster Dendrograma para coches", xlab= "Veh�culos estudiados",sub=paste("Metodo=", mtdo_usado,"; Distancia=euclidea al cuadrado", sep = ""))

  for(i in 2:(nrow(df_coches.def.ilu.ac)-1)){
    gr <- cutree(HClust_aux, i)
    distgr <- grpdist(gr)
    mt <- cor(dist_coches, distgr, method="pearson")
    kt[i,2] <- mt
    }
  # Calculamos mediante el estad�stico de Mantel el n�mero de clusters �ptimo:
  cat('\nM�todo: "', mtdo_usado,'". k �ptimo = ',k.best <- which.max(kt$r),'\nComposici�n de los clusters:\n', sep = "")

  print(summary(as.factor(cutree(HClust_aux, k = k.best))))

}

```


Todos los m�todos testeados nos muestran particiones muy similares. Con lo  cual, para tomar la decisi�n final, se representan los veh�culos en el plano factorial dentro de sus clusters.


```{r}

  library("RcmdrMisc")

# Funcion que genera el objeto con la figura
getGgplotACM <- function(hCluter_figura, k.optimo=3, titulo = "") {
  # Se crea una columna con el numero de cluster asignado a ese coche
  df_coches.def.ilu.ac$num_cluster <- assignCluster(model.matrix(~-1 +  C1+C2, df_coches.def.ilu.ac), 
                                                    df_coches.def.ilu.ac, cutree(hCluter_figura, k = k.optimo))
  
  ##### coches$num_cluster <- assignCluster(model.matrix(~-1 +  C1+C2, df_coches.def.ilu.ac), df_coches.def.ilu.ac, cutree(hCluter_figura, k = k.optimo))
  
  figura <- ggplot(df_coches.def.ilu.ac, aes(x=C1, y=C2, color=num_cluster)) +
    geom_point() + geom_rug() +
    geom_vline(xintercept = 0) + geom_hline(yintercept = 0) +
    geom_text_repel(aes(label=rownames(df_coches.def.ilu.ac)),hjust=0, vjust=0, size=3.3) +
    #geom_text(aes(label=rownames(df_coches.def.ilu.ac)),hjust=0, vjust=0, size=3.3) +
    theme(legend.position="top") + ggtitle(titulo)
  
  return(figura)
}

library(gridExtra)
library(grid)
library(ggplot2)
library(lattice)

# M�todo: "ward.D". k �ptimo = 3
HClust_aux <- hclust(dist_coches , method= "ward.D")
g1grob <- ggplotGrob(getGgplotACM(HClust_aux, k.optimo = 3, titulo = "M�todo ward.D"))

# M�todo: "complete". k �ptimo = 3
HClust_aux <- hclust(dist_coches , method= "complete")
g2grob <- ggplotGrob(getGgplotACM(HClust_aux, k.optimo = 3, titulo = "M�todo complete"))

# M�todo: "single". k �ptimo = 4
HClust_aux <- hclust(dist_coches , method= "single")
g3grob <- getGgplotACM(HClust_aux, k.optimo = 4, titulo = "M�todo single")

# M�todo: "average". k �ptimo = 3
HClust_aux <- hclust(dist_coches , method= "average")
g4grob <- getGgplotACM(HClust_aux, k.optimo = 3, titulo = "M�todo average")

#M�todo: "median". k �ptimo = 3
HClust_aux <- hclust(dist_coches , method= "median")
g5grob <- getGgplotACM(HClust_aux, k.optimo = 3, titulo = "M�todo median y M�todo mcquitty")

## # M�todo: "mcquitty". k �ptimo = 3
## HClust_aux <- hclust(dist_coches , method= "mcquitty")
## g6grob <- getGgplotACM(HClust_aux, k.optimo = 3, titulo = "M�todo mcquitty")

# M�todo: "centroid". k �ptimo = 3
HClust_aux <- hclust(dist_coches , method= "centroid")
g7grob <- getGgplotACM(HClust_aux, k.optimo = 3, titulo = "M�todo centroid")


texto1=textGrob("Composici�n de los clusters: 35 32 53", gp=gpar(fontsize=10))
texto2=textGrob("Composici�n de los clusters: 37 32 51", gp=gpar(fontsize=10))
texto3=textGrob("Composici�n de los clusters: 68 20 27  5", gp=gpar(fontsize=10))
texto4=textGrob("Composici�n de los clusters: 34 33 53", gp=gpar(fontsize=10))
texto5=textGrob("Composici�n de los clusters: 42 27 51", gp=gpar(fontsize=10))
texto7=textGrob("Composici�n de los clusters: 34 33 53", gp=gpar(fontsize=10))


# Se muestran la figuras
grid.arrange(
  grobs = list(g1grob, g2grob, texto1, texto2),
  heights = c(9, 1),
  layout_matrix = rbind(c(1, 2),
                        c(3, 4))
)

grid.arrange(
  grobs = list(g3grob, g4grob, texto3, texto4),
  heights = c(9, 1),
  layout_matrix = rbind(c(1, 2),
                        c(3, 4))
)

grid.arrange(
  grobs = list(g5grob, g7grob, texto5, texto7),
  heights = c(9, 1),
  layout_matrix = rbind(c(1, 2),
                        c(3, 4))
)

```

Llegados a este punto, ser�a conveniente consultar con personas especialistas en los datos que aqu� se tratan para decidir el nu�ro final de cluster con su composici�n. De la revisi�n de la bibliograf�a se deduce que, en caso de duda, el m�todo m�s recomendable, probablemente, sea el del **promedio entre grupos**, con la **distancia clidea al cuadrado** como medida de proximidad [@urkaregi2017apuntes].



```{r}
# La matriz dist_coches contiene las distancias eclideas al cuadrado.
#
# Se toma el metodo mas adecuado
mtdo_usado_acm = "median"
HClust.2 <- hclust(dist_coches , method = mtdo_usado_acm)

# Se elige el numero de clusters.
for(i in 2:(nrow(df_coches.def.ilu.ac)-1)){
  gr <- cutree(HClust.2, i)
  distgr <- grpdist(gr)
  mt <- cor(dist_coches, distgr, method="pearson")
  kt[i,2] <- mt
  }

num_clusters_acm <- which.max(kt$r)
```

Distribuci�n de los veh�culos en cada cluster:

```{r}
summary(as.factor(cutree(HClust.2, k = num_clusters_acm))) 
```



```{r}
library(factoextra)
# Dendograma
fviz_dend(HClust.2, k=num_clusters_acm,rect = T, show_labels = T)

# Mapa con vehiculos en el plano factorial.
# M�todo: "ward.D". k �ptimo = num_clusters_acm
coches_grob <- getGgplotACM(HClust.2, k.optimo = num_clusters_acm, titulo = paste("M�todo", mtdo_usado_acm)) +
  theme(legend.position="right")

# variables activas en el plano factorial
# var_act_grob se calculo mas arriba

# Se muestran la figuras
grid.arrange(
  grobs = list(var_act_grob, coches_grob),
  ## heights = c(9, 1),
  layout_matrix = rbind(c(1, 2))
)

```


```{r}
# Se crea una columna con el numero de cluster asignado a ese coche
coches$num_cluster_acm <- assignCluster(model.matrix(~-1 +  C1+C2, df_coches.def.ilu.ac), 
                                                  df_coches.def.ilu.ac, cutree(HClust.2, k = num_clusters_acm))


for (i in 1:num_clusters_acm) {
  cat(paste("\n Cluster n�mero:", i, "\t\n"))

  print(paste(rownames(coches)[which(coches$num_cluster_acm == i)], 
              coches$MARCA[which(coches$num_cluster_acm == i)], 
              coches$MODELO[which(coches$num_cluster_acm == i)],
              coches$TIPO[which(coches$num_cluster_acm == i)]))
}
```




## Compara los an�lisis cluster con An�lisis Factorial.

**Compara los cluster obtenidos con los del apartado de An�lisis Factorial, indicando, de manera razonada, cu�l te parece el an�lisis m�s adecuado**

FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   
FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   
FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   
FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   
FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   
FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   
FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   FALTA   



***

\pagebreak


---