analisis.discriminante.fnc

Análisis discriminante

Partiremos de la base de datos iris (disponible en la librería MASS) sobre la que llevaremos a cabo un análisis discriminante para intentar encontrar combinaciones lineales de 4 variables cuantitativas que maximicen la distancia entre 3 grupos de especies vegetales.

head(iris)

Sepal.Length Sepal.Width Petal.Length Petal.Width Species

1 5.1 3.5 1.4 0.2 setosa

2 4.9 3.0 1.4 0.2 setosa

3 4.7 3.2 1.3 0.2 setosa

4 4.6 3.1 1.5 0.2 setosa

5 5.0 3.6 1.4 0.2 setosa

6 5.4 3.9 1.7 0.4 setosa

Solicitamos la distribución de frecuencias de la variable de agrupación Species.

frecuencias.fnc(iris, variable='Species')

#--------------------------------------------------------------------------------

# TABLA DE FRECUENCIAS

#--------------------------------------------------------------------------------

$Species

$Species$n.total

[1] 150

$Species$tabla

setosa versicolor virginica

50 50 50

modelo1= analisis.discriminante.fnc(iris, grupo='Species')

Observa que hemos asignado la salida de la función al objeto modelo1. Esto tendría especial utilidad si quisiéramos utilizar el modelo estimado tanto para extraer las puntuaciones discriminantes para cada sujeto como para llevar a cabo predicciones de grupo de pertenencias para sujetos no incluidos en el proceso de estimación.

Por defecto si solo explicitas el argumento obligatorio grupo la función analizará todas las variables de la base de datos (en este ejemplo iris) y obtendrás:

1. Estimación paso a paso para la selección de las variables discriminantes.

2. Tabla de medias de las p variables seleccionadas en los J niveles del factor incluido en el argumento grupo.

3. Manova omnibus.

4. Funciones discriminantes estimadas.

5. Valor de Wilks para cada variable si fuese eliminada de la o las funciones.

6. Prueba de homogeneidad de las varianzas de MBox.

7. Centroides, coeficientes no típicos, típicos y estructura para las p variables seleccionadas en las funciones estimadas.

8. Tabla de contingencia de frecuencias absolutas y relativas de la clasificación de los J grupos a partir de las funciones discriminantes.

Veamos la llamada anterior incluyendo el nuevo argumento grafica=T.

modelo1=analisis.discriminante.fnc(iris, variables=1:4,

paso.a.paso=T, grupo='Species', grafica=T)

En la salida posterior podrás ver lo que "cuentan" estas gráficas. En la de la izquierda puedes ver que la primera función separa de forma muy clara a la especie setosa del resto (virgínica y versicolor), mientras que la segunda apenas aporta nada en dicha separación.

Con el argumento grafica=T obtendremos una gráfica por cada función discriminante con los histogramas de las puntuaciones discriminantes en cada uno de los J grupos.

Estimación Bootstrap

Si incluimos el argumento bootstrap=T, la función llevará a cabo la extracción de 1000 muestras bootstrap para calcular el sesgo en la estimación de los coeficientes no típicos, así como el error típico bootstrap el cual se utilizará para estimar el estadístico "t" así como el intervalo de confianza con corrección de sesgo (adjusted bootstrap percentile BCa method) Efron (1987))

analisis.discriminante.fnc(iris, variables=1:4,

paso.a.paso=T, grupo='Species', boot=T)

#------------------------------------------------------------------

# Extrayendo 1000 muestras bootstrap

#------------------------------------------------------------------

$Bootstrap

Original Sesgo ET.boot t.boot IC.BCa.inf IC.BCa.sup

Petal.LengthLD1 -2.201212 -0.059461 0.33141 -6.642026 -2.81476 -1.4582

Sepal.WidthLD1 1.534473 0.066349 0.38397 3.996358 0.80602 2.2607

Petal.WidthLD1 -2.810460 -0.097774 0.49104 -5.723475 -3.72158 -1.8067

Sepal.LengthLD1 0.829378 0.016132 0.29701 2.792456 0.20618 1.3879

Petal.LengthLD2 0.931921 -1.004522 1.09814 0.848636 -0.93315 2.5729

Sepal.WidthLD2 -2.164521 2.283580 2.20273 -0.982652 -3.46220 2.0637

Petal.WidthLD2 -2.839188 3.051409 2.96426 -0.957805 -5.49121 2.6505

Sepal.LengthLD2 -0.024102 0.020111 0.60601 -0.039772 -1.10872 1.2535

En la tabla de resultados Bootstrap podemos ver en primer lugar los valores originales de la estimación de los coeficientes no típicos, seguido del sesgo (diferencia entre el valor estimado y el promedio de las 1ooo replicaciones bootstrap), el error típico boostrap (desviación típica a partir de las 1000 muestras), el estadístico "t" (valor del parámetro dividido por el error típico bootstrap) y el intervalo de confianza corregido por el sesgo de asimetría (Efron 1987).

Validación cruzada (jacknife)

Con el argumento VC=T, la función lleva a cabo la asignación a los grupos excluyendo un sujeto en cada asignación (leave one out cross validation).

modelo1 = analisis.discriminante.fnc(iris, grupo='Species', VC=T)

A la salida del discriminante se añade la tabla de validación cruzada por jacknife.

$validacion.cruzada.jacknife

Predicho

Observado setosa versicolor virginica

setosa 50 0 0

versicolor 0 48 2

virginica 0 1 49

$prop.validacion.cruzada

Predicho

Observado setosa versicolor virginica

setosa 1.00 0.00 0.00

versicolor 0.00 0.96 0.04

virginica 0.00 0.02 0.98

Análisis discriminante Cuadrático

El análisis discriminante clásico (lineal) presenta un importante requisito relativo a la necesidad de cumplimiento del supuesto de igualdad de la matriz de varianzas y covarianzas entre los J niveles del factor de agrupación. Si se viola este supuesto el usuario puede añadir el argumento ADC=T indicando así que desea además la estimación del modelo de análisis discriminante cuadrático para el cual no se requiere dicho supuesto. Si se hace verdadero (T) este argumento, la función añadirá a la estimación lineal, la cuadrática. Esta podrá utilizarse para realizar las predicciones de las observaciones a los grupos en lugar del modelo lineal.

analisis.discriminante.fnc(iris, grupo='Species', ADC=T)

$modelo.adc

Call:

qda(eval(parse(text = modelo)), data = dat, prior = rep(1/n.group,

n.group))

Prior probabilities of groups:

setosa versicolor virginica

0.3333333 0.3333333 0.3333333

Group means:

Petal.Length Sepal.Width Petal.Width Sepal.Length

setosa 1.462 3.428 0.246 5.006

versicolor 4.260 2.770 1.326 5.936

virginica 5.552 2.974 2.026 6.588

$adc.pred

$adc.pred$variables

[1] "Species:qda.pred"

$adc.pred$tabla

qda.pred setosa versicolor virginica

Species

setosa 50 0 0

versicolor 0 48 2

virginica 0 1 49

$adc.pred$prop.fila

qda.pred setosa versicolor virginica

Species

setosa 1.00 0.00 0.00

versicolor 0.00 0.96 0.04

virginica 0.00 0.02 0.98

Puntuaciones discriminantes

Cuando estimamos un modelo discriminante lineal es posible extraer la puntuaciones discriminantes para cada sujeto en las funciones estimadas así como el grupo predicho para cada observación. Para ello utilizaremos la función: puntuacion.discriminante.fnc, cuyos únicos argumentos serán el nombre de la base de datos (iris) y el modelo lineal estimado (modelo1). Es muy importante reseñar aquí que el usuario puede utilizar esta función para extender la predicción del modelo a otra base de datos diferente a la utilizada para estimarlo. Obviamente estos nuevos datos deben contener las variables discriminantes utilizadas en la estimación.

iris.2 = puntuacion.discriminante.fnc(iris, modelo=modelo1)

*** Esta es la cabecera de tus nuevos datos. Se han incluido las puntuaciones discriminantes ***

Sepal.Length Sepal.Width Petal.Length Petal.Width Species LD1 LD2 grupo.predicho1

1 5.1 3.5 1.4 0.2 setosa 8.061800 -0.3004206 setosa

2 4.9 3.0 1.4 0.2 setosa 7.128688 0.7866604 setosa

3 4.7 3.2 1.3 0.2 setosa 7.489828 0.2653845 setosa

4 4.6 3.1 1.5 0.2 setosa 6.813201 0.6706311 setosa

5 5.0 3.6 1.4 0.2 setosa 8.132309 -0.5144625 setosa

6 5.4 3.9 1.7 0.4 setosa 7.701947 -1.4617210 setosa

El usuario podría ahora si lo desease llevar a cabo un Anova con contrastes poshoc par a par u ortogonales utilizando como variable dependiente la puntuación discriminante de cada sujeto.

Anova.fnc(iris.2, fac.inter='Species', vd='LD1',

poshoc='Species', grafica=F)

#--------------------------------------------------------------------------------

# Contrastes Poshoc

#--------------------------------------------------------------------------------

---------------------------------------------------------------------------

*** WARNING ***

Por defecto se lleva a cabo la correccion de la probabilidad asociada (p.val)

para contrastes poshoc de hochberg ajuste.p='hochberg'.

El numero de contrastes que se utiliza para la correccion son los realizados

dentro de cada subfamilia (menos conservador). Si deseas que se incluyan

los contrastes de la famila incluye el argumento: ajustep.tod.pares=T

---------------------------------------------------------------------------

$poshoc

$poshoc$Species

contrast estimate SE df t.ratio p.value

setosa - versicolor 9.43 0.2 147 47.163 <.0001

setosa - virginica 13.39 0.2 147 66.951 <.0001

versicolor - virginica 3.96 0.2 147 19.788 <.0001

P value adjustment: hochberg method for 3 tests

Creamos una familia ortogonal con un primer contraste que compara a la especie setosa versus las dos restantes y el segundo ortogonal al primero donde comparamos la especie virginica contra la versicolor.

mis.contrastes=list(set.vs.vv= c(-2,1,1), ver.vs.vir=c(0,1,-1))

Anova.fnc(iris.2, fac.inter='Species', vd='LD1',

poshoc='Species', contrastes=mis.contrastes)

#--------------------------------------------------------------------------------

# Contrastes Poshoc

#--------------------------------------------------------------------------------

$poshoc

$poshoc$`contrastes:Species`

set.vs.vv ver.vs.vir

setosa -2 0

versicolor 1 1

virginica 1 -1

$poshoc$Species

contrast estimate SE df t.ratio p.value

set.vs.vv -22.82 0.346 147 -65.884 <.0001

ver.vs.vir 3.96 0.200 147 19.788 <.0001

P value adjustment: hochberg method for 2 tests

Eficacia de un modelo predictivo

Partiendo de la misma base de datos iris vamos a realizar un procedimiento de estimación del modelo discriminante en la mitad de los datos y utilizaremos la otra mitad para probar la bondad del mismo.

1. Crearemos una nueva variable que indique el número del caso o registro que llamaremos registro.

2. Extraeremos una muestra aleatoria por cada nivel de la variable de agrupación Species.

3. Sobre esa muestra estimaremos el modelo discriminante.

4. Extraeremos de la base de datos iris los registros que no han participado en la muestra aleatoria anterior. La llamaremos arbitrariamente resto.

5. Estimaremos el grupo de pertenencia de cada registro de resto utilizando el modelo estimado a partir de la muestra.

6. Comprobaremos la bondad de ajuste del modelo mediante la tabla de contingencia de la variable Species y la variable predicha de Species.

frecuencias.fnc(iris, 'Species')

De esta manera vimos en la introducción anterior que había 50 registros por nivel de Species.

1.- Creamos la variable registro.

iris$registro= 1:150

2.- Extraemos una muestra aleatoria de tamaño 25 por nivel de Species.

muestra=extrae.muestra.fnc(iris, que.factor='Species', n=25, ID='registro')

3.- Estimamos el modelo discriminante sobre esta muestra y lo guardamos en ad.2.

ad.2 = analisis.discriminante.fnc(muestra, variables=1:4, grupo='Species')

4.- Extraemos de iris los registros que no forman parte de la muestra (Operadores Lógicos)

indice = iris$registro %in% muestra$registro

resto= iris[!indice, ]

5.- Aplicamos el modelo predictivo al resto de los datos

prediccion=puntuacion.discriminante.fnc(resto, modelo=ad.2)

6.- Comprobamos la bondad predictiva del modelo solicitando la tabla de contingencia del cruce de la variable observada Species y su versión predicha por el modelo grupo.predicho1.

frecuencias.fnc(prediccion, variables='Species:grupo.predicho1')

#--------------------------------------------------------------------------------

# TABLA DE FRECUENCIAS

#--------------------------------------------------------------------------------

$variables

[1] "Species:grupo.predicho1"

$tabla

grupo.predicho1 setosa versicolor virginica

Species

setosa 25 0 0

versicolor 0 25 0

virginica 0 2 23

Subir->