ME731: Análise Discriminnate

Case 1:

O data set MEDIFIS contém informações sobre diversas medidas físicas de homens (1) e mulheres (0). Estamos interessados em classificar pessoas por sexo, segundo as medidas físicas que elas apresentam.

Os dados estão disponíveis aqui

library(MASS)
library(dplyr)
medifis <- read.table("https://raw.githubusercontent.com/jvega68/EA3/master/datos/medifis.dat")
glimpse(medifis)

Rows: 27
Columns: 8
$ V1 <int> 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, …
$ V2 <dbl> 159, 164, 172, 167, 164, 161, 168, 181, 18…
$ V3 <dbl> 49, 62, 65, 52, 51, 67, 48, 74, 74, 50, 65…
$ V4 <dbl> 36, 39, 38, 37, 36, 38, 39, 43, 41, 36, 36…
$ V5 <dbl> 68.0, 73.0, 75.0, 73.0, 71.0, 71.0, 72.5, …
$ V6 <dbl> 42.0, 44.0, 48.0, 41.5, 44.5, 44.0, 41.0, …
$ V7 <dbl> 57.0, 55.0, 58.0, 58.0, 54.0, 56.0, 54.5, …
$ V8 <dbl> 40, 44, 44, 44, 40, 42, 43, 47, 47, 41, 41…

colnames(medifis) <- c("sexo", "altura", "peso", "pe", "braco", "costas", "diam_craneo", "joelho")
glimpse(medifis)

Rows: 27
Columns: 8
$ sexo        <int> 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, …
$ altura      <dbl> 159, 164, 172, 167, 164, 161, 168, 181,…
$ peso        <dbl> 49, 62, 65, 52, 51, 67, 48, 74, 74, 50,…
$ pe          <dbl> 36, 39, 38, 37, 36, 38, 39, 43, 41, 36,…
$ braco       <dbl> 68.0, 73.0, 75.0, 73.0, 71.0, 71.0, 72.…
$ costas      <dbl> 42.0, 44.0, 48.0, 41.5, 44.5, 44.0, 41.…
$ diam_craneo <dbl> 57.0, 55.0, 58.0, 58.0, 54.0, 56.0, 54.…
$ joelho      <dbl> 40, 44, 44, 44, 40, 42, 43, 47, 47, 41,…

EAD

library(tidyr)
library(ggplot2)
medifis %>% pivot_longer(!sexo, names_to = "variaveis") %>% ggplot() + 
  geom_boxplot(aes(y = value, group = factor(sexo), fill = factor(sexo))) + 
  facet_wrap(.~variaveis, scales = "free_y")

Suposições:

biotools::boxM(medifis[,-1], medifis$sexo)

    Box's M-test for Homogeneity of Covariance Matrices

data:  medifis[, -1]
Chi-Sq (approx.) = 39.871, df = 28, p-value = 0.06789

mvnormalTest::mardia(medifis[medifis$sexo == 0, -1])

$mv.test

$uv.shapiro
            W      p-value UV.Normality
altura      0.9509 0.5383  Yes         
peso        0.8966 0.0845  Yes         
pe          0.8904 0.068   Yes         
braco       0.94   0.3821  Yes         
costas      0.9515 0.5477  Yes         
diam_craneo 0.9602 0.6951  Yes         
joelho      0.9475 0.4865  Yes         

mvnormalTest::mardia(medifis[medifis$sexo == 1, -1])

$mv.test

$uv.shapiro
            W      p-value UV.Normality
altura      0.9613 0.8028  Yes         
peso        0.9453 0.5702  Yes         
pe          0.9602 0.7868  Yes         
braco       0.9215 0.2984  Yes         
costas      0.9567 0.7361  Yes         
diam_craneo 0.9822 0.9909  Yes         
joelho      0.8495 0.0362  No          

Análise Discriminante

discrinante_linear <- lda(sexo ~ ., data = medifis)
discrinante_linear

Call:
lda(sexo ~ ., data = medifis)

Prior probabilities of groups:
        0         1 
0.5555556 0.4444444 

Group means:
    altura  peso       pe    braco   costas diam_craneo
0 161.7333 55.60 36.83333 70.03333 43.33333    56.63333
1 177.5833 74.25 41.66667 77.75000 49.00000    58.00000
    joelho
0 41.06667
1 45.62500

Coefficients of linear discriminants:
                     LD1
altura      -0.088142097
peso        -0.005872569
pe           0.650276334
braco        0.134222958
costas       0.123645380
diam_craneo -0.202387598
joelho       0.101000944

classificacao <- predict(discrinante_linear, medifis[, -1])$class
table(classificacao, medifis$sexo)

             
classificacao  0  1
            0 15  0
            1  0 12

Significa que o modelo nunca erra?

Como podemos estimar a probabilidade de erro?

• Fórmula fechada (apenas valido neste caso em particular)
• Validação cruzada.

Coefficients of linear discriminants?

discrinante_linear$scaling

                     LD1
altura      -0.088142097
peso        -0.005872569
pe           0.650276334
braco        0.134222958
costas       0.123645380
diam_craneo -0.202387598
joelho       0.101000944

• Fisher (1938), sob outra abordagem que não assume normalidade, derivou uma forma de classificar (método conhecido como Análise Discriminante Linear de Fisher).
• A ideia é transformar \(\textbf{x} \in \mathbb{R}^p\) em \(z = \alpha'\textbf{x} \in \mathbb{R}\) de forma que \(z\) separe \(P_1\) de \(P_2\) o máximo que puder.
• Sejam \(z_{11}, \cdots, z_{1n_1}\) as transformações \(z = \alpha'\textbf{x}\) para \(\textbf{x} \in P_1\) e sejam \(z_{21}, \cdots, z_{2n_2}\) as transformações \(z = \alpha'\textbf{x}\) para \(\textbf{x} \in P_2\)
• Uma forma de medir a separação de ambos os grupos é atraves das suas médias \(\bar{z}_1\) e \(\bar{z}_2\). Assim queremos maximizar \[\phi = \dfrac{(\bar{z}_1 - \bar{z}_2)^2}{s^2_z} = \dfrac{(\alpha' [\bar{x}_1 - \bar{x}_2])^2}{\alpha' \textbf{S}_x \alpha} = \dfrac{\alpha' \overbrace{(\bar{x}_1 - \bar{x}_2)(\bar{x}_1 - \bar{x}_2)'}^{\textbf{B}}\alpha}{\alpha' \textbf{S}_x \alpha}\]

Derivando w.r.t \(\alpha\) e igualando a zero, temos: \[\dfrac{2\textbf{B}\alpha(\alpha'\textbf{S}_x \alpha) - \alpha' \textbf{B}\alpha 2 \textbf{S}_x \alpha}{(\alpha'\textbf{S}_x \alpha)^2} = 0\]

\[\textbf{B}\alpha \alpha'\textbf{S}_x\alpha = \alpha'\textbf{B}\alpha \textbf{S}_x \alpha \rightarrow \textbf{B}\alpha = \underbrace{\dfrac{\alpha'\textbf{B}\alpha}{\alpha'\textbf{S}_x \alpha}}_{\phi} \textbf{S}_x \alpha \rightarrow \textbf{S}_x^{-1}\textbf{B}\alpha = \phi \alpha\]

x <- medifis %>% filter(sexo == 0) %>% dplyr::select(!sexo)
y <- medifis %>% filter(sexo == 1) %>% dplyr::select(!sexo)
mu_1 <- matrix(apply(x, 2, mean), ncol = 1)
mu_2 <- matrix(apply(y, 2, mean), ncol = 1)
n_1 <- nrow(x)
n_2 <- nrow(y)
S <- (n_1 - 1)*cov(x) + (n_2 - 1)*cov(y)
S <- S/(n_1 + n_2 - 2)
B <- (mu_1 - mu_2) %*% t(mu_1 - mu_2)
eigen(solve(S) %*% B)$vectors[,1]

[1]  0.122807669+0i  0.008182203+0i -0.906024743+0i
[4] -0.187011759+0i -0.172274105+0i  0.281985000+0i
[7] -0.140723796+0i

discrinante_linear$scaling/sqrt(sum(discrinante_linear$scaling^2))

                     LD1
altura      -0.122807669
peso        -0.008182203
pe           0.906024743
braco        0.187011759
costas       0.172274105
diam_craneo -0.281985000
joelho       0.140723796

Observação:

• Validação cruzada é uma técnica para estimar a probabilidade de erro (pois não sempre podemos estimar esta probabilidade de forma analítica).
• Na prática, quando o tamanho amostral é pequeno, mesmo que \(\Sigma_1 \neq \Sigma_2\). Análise discriminante linear pode apresentar um melhor desempenho. Isto é devido a que análise discrinante quadrática implica estimar muitos parâmetros adicionais.
• Existem diversas extensões de análise discriminante (algumas delas serão apresentadas nos trabalhos da próxima semana)