analisis.discriminante.fnc
Copia, Pega y Adapta
modelo.ad= analisis.discriminante.fnc(datos, paso.a.paso=F,
grupo='diagnostico', variables=2:6, grafica=T)
datos = puntuacion.discriminante.fnc(datos, modelo=modelo.ad)
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:
Estimación paso a paso para la selección de las variables discriminantes.
Tabla de medias de las p variables seleccionadas en los J niveles del factor incluido en el argumento grupo.
Manova omnibus.
Funciones discriminantes estimadas.
Valor de Wilks para cada variable si fuese eliminada de la o las funciones.
Prueba de homogeneidad de las varianzas de MBox.
Centroides, coeficientes no típicos, típicos y estructura para las p variables seleccionadas en las funciones estimadas.
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.
Crearemos una nueva variable que indique el número del caso o registro que llamaremos registro.
Extraeremos una muestra aleatoria por cada nivel de la variable de agrupación Species.
Sobre esa muestra estimaremos el modelo discriminante.
Extraeremos de la base de datos iris los registros que no han participado en la muestra aleatoria anterior. La llamaremos arbitrariamente resto.
Estimaremos el grupo de pertenencia de cada registro de resto utilizando el modelo estimado a partir de la muestra.
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