Laboratorio 1 Statistica multivariata m

# LABORATORIO 1 STATISTICA MULTIVARIATA

# ESERCIZIO 1

#a) Caricare il dataset wine dal pacchetto "HDclassif". Esso contiene i

# risultati di un'analisi chimica effettuata su 178 vini coltivati in una

# regione italiana, dove i vini provengono da 3 diversi vitigni.

# L'analisi ha determinato le quantità di 13 componenti presenti in ciascuno

# dei 3 tipi di vino.

rm(list=ls())

install.packages("HDclassif")

library(HDclassif)

data("wine")

#b) Rappresentare il dataset graficamente, non considerando la variabile

# categoriale class.

str(wine)

# 'data.frame': 178 obs. of 14 variables:

# $ class: int 1 1 1 1 1 1 1 1 1 1 ...

# $ V1 : num 14.2 13.2 13.2 14.4 13.2 ...

# $ V2 : num 1.71 1.78 2.36 1.95 2.59 1.76 1.87 2.15 1.64 1.35 ...

# $ V3 : num 2.43 2.14 2.67 2.5 2.87 2.45 2.45 2.61 2.17 2.27 ...

pairs(wine[, -1]) # rimuoviamo la prima variabile

# Ottengo una È una matrice di dispersione: mostra, per ogni coppia di

# variabili, la relazione tra i loro valori. Serve a individuare rapidamente

# correlazioni, relazioni non lineari, gruppi e valori anomali.

# Una nuvola stretta e crescente, dal basso a sinistra verso l’alto a destra,

# indica una forte associazione positiva: quando una variabile aumenta, tende

# ad aumentare anche l’altra.

# Una nuvola decrescente indicherebbe un’associazione negativa. Una forma quasi

# circolare o molto dispersa indica invece un legame lineare debole.

# Forme curve, triangolari o “a ventaglio” possono indicare relazioni non

# lineari, limiti naturali delle variabili o varianza non costante.

#b) Calcolare la distanza di Mahalanobis, settare un valore di soglia e

# individuare eventuali outliers, rappresentandoli graficamente con la

# funzione pairs.

data <- wine[, -1]

n.obs <- nrow(data)

p <- ncol(data)

DM2 <- mahalanobis(data, center = apply(data,2,mean), cov = ((n.obs - 1)/n.obs)

* cov(data))

DM <- sqrt(DM2)

DM

soglia <- sqrt(qchisq(0.05, p, lower.tail = FALSE))

which(DM > soglia)

# [1] 14 60 69 70 72 74 75 79 96 97 111 116 122 159 160

out <- DM > soglia

pairs(data, col = ifelse(out, "red", "black"), pch = ifelse(out, 4, 19))

#c) Cosa succede se diminuiamo l'alpha del valore soglia? Perchè questo è un

# risultato atteso?

soglia2 <- sqrt(qchisq(0.01, p, lower.tail = FALSE))

which(DM > soglia2)

# [1] 60 70 74 96 111 122 159

out2 <- DM > soglia2

pairs(data, col = ifelse(out2, "red", "black"), pch = ifelse(out2, 4, 19))

#d) Dal grafico risultante dalla funzione "pairs", possiamo dedurre una

# non-normalità dei dati?

pairs(data)

#e) Analizzare la normalità delle singole variabili. Cosa possiamo dedurre?

par(mfrow = c(4,4))

for(j in 1:p){

qqnorm(data[ ,j],

xlab = "Quantili teorici",

ylab = "Quantili empirici",

main = colnames(data)[j],

cex.lab = 1.3)

qqline(data[, j])

}

#f) Costruire un QQPlot per confrontare i quantili empirici di u1,...,un con

# i quantili teorici della suddetta distribuzione beta.

# Distribuzione empirica

DM2 <- mahalanobis(data, center