Esercitazioni Analisi statistica multivariata

# se plottiamo un modello lineare, il parametro "which" fornisce dei grafici pre-stabiliti:

# which=1 => residui VS valori teorici

# which=2 => Q-Q plot dei residui

# which=3 => scale-location: valuta l'omoschedasticità dei residui

# which=4 => residui VS punti di leva: vedi sotto

summary(hatvalues(m1))

# hatvalues: resituisce il vettore n x 1 che corrisponde alla diagonale della matrice di

proiezione H.

# Serve per l'analisi dei punti di leva

length(m1$coe_icients)/nrow(data)

# la media degli hatvalues corrisponde al numero di coe_icienti (regressori+intercetta) fratto

la

# numerosità del campione

par(mfrow=c(2,2))

plot(hatvalues(m1), ylab="Leverage", xlab="Indice", type="h" )

abline(h=3*mean(hatvalues(m1)), col="red")

# traccia una linea orizzontale, che ci permette di vedere quali punti hanno valore H

# superiore a 3 volte la media.

plot(x=hatvalues(m1), y=rstudent(m1), ylab="R Student", xlab="Leverage" )

# rstandard: come rstandard ma restituisce i residui studentizzati.

# Il grafico di sopra plotta gli hatvalues con i residui studentizzati, per valutare la loro potenza

# di leva.

abline(v=2*mean(hatvalues(m1)))

# traccia una linea verticale, che ci permette di vedere quali punti hanno valore H superiore

# a 2 volte la media. Dato che i valori di leva sono sull'asse x, serve una abline verticale

abline(h=c(-2,2))

# traccia due linee orizzontali, che ci permettono di vedere quali residui studentizzati hanno

# valore superiore a 2 o inferiori a -2

plot(cooks.distance(m1), ylab="Cook's D", xlab="Indice", type="h" )

# cooks.distance: resituisce il vettore n x 1 contenente le distanze di Cook per ciascuna u.s.

abline(h=4/(nrow(data)-length(m1$coe_icients) ))

# un criterio per trovare i punti influenti è quello di considerare quelli con distanza di Cook

# superiore a 4/(n°unità - n°regressori - 1)

which( hatvalues(m1)> 2*mean(hatvalues(m1)))

which( abs(rstudent(m1))>2 )

which(cooks.distance(m1)> 4/(nrow(data)-length(m1$coe_icients)-1))

ES. 2

rm(list=ls())

setwd("H:/Il mio Drive/sgi modelli statistici 2023/es2_final")

data<-read.csv("airfare.csv", sep=",", header=T)

lm1<-lm( fare ~ concen+dist+passen ,data)

summary(lm1)

alm1<-anova(lm1)

# anova con un modello: valuta la significatività dei parametri.

alm1

# varianza del modello

alm1$`Sum Sq`[4]/alm1$Df[4]

# uguale a

sum(lm1$residuals^2)/lm1$df.residual

lm0<-lm( fare ~ 1 ,data)

# modello di grado 0

anova(lm1, lm0)

# anova con due modelli: valuta se c'è una di_erenza significativa nella capacità esplicativa

# dei due modelli.

se<-sqrt(vcov(lm1)["dist","dist"])

# calcola lo standard error del modello

# vcov: matrice di varianze e covarianze di un modello lineare. Da non confondere con var,

che è la

# matrice var/cov di una distribuzione di dati.

lm1$coe_icients["dist"] +c(-1,1)*qt(0.90, lm1$df.residual)*se

lm1$coe_icients["dist"] +c(-1,1)*qt(0.95, lm1$df.residual)*se

lm1$coe_icients["dist"] +c(-1,1)*qt(0.99, lm1$df.residual)*se

# creazione degli intervalli di confidenza

# qt: funzione che determina i quantili di una T di Student. Richiede due parametri in entrata,

# il valore di probabilità associato al quantile e i gradi di libertà della distribuzione.

confint(lm1, level=0.9)

# confint: funzione che determina gli intervalli di confidenza per i parametri di un modello

summary(lm1, cov=T)$cov

summary(lm1, cor=T)$cor

library(ellipse)

plot(ellipse(lm1, c(3,4) ), type="l" )

# ellipse: traccia un ellisse di confidenza che contiene al 95% i valori congiunti dei parametri

# dei regressori scelti. In questo caso, dato che il secondo argomento è c(3,4), verranno scelti

il

# terzo e quarto regressore, quindi "dist" e "passen".

points(m1$coe_icients["dist"],m1$coe_icients["passen"],pch=19 )

abline(v=confint(lm1, level=0.95)["dist",], lty=2 )

abline(h=confint(lm1, level=0.95)["passen",], lty=2 )

# si possono usare gli estremi di un intervallo di confidenza per plottare una abline

lm2<-lm(fare ~ concen,data)

anova(lm1,lm2)

data<-read.csv("dati.txt", sep = "\t")

data<-data[,-1]

summary(data)

cor(data)

pairs.panels(data, cex.cor=2)

data2<-data[,c("peso","bicipite")]

pairs.panels(data2, lm=T)

hist(data$peso, main="peso")

hist(data$bicipite, main="bicipite")

lin<-lm(peso~bicipite, data2)

summary(lin)

plot(lin, which=1)

linlog<-lm(peso~I(log(bicipite)),data2)

# per trasformare una covariata, bisogna mettere la trasformazione dentro I(). Esempio: se

volessi

# considerare x al quadrato, metterei I(x^2)

summary(linlog)

##se la cicrf del bicipite varia del 1% ->2.31 lb

plot(linlog, which=1)

loglog<-lm(I(log(peso))~I(log(bicipite)),data2)

summary(loglog)

##variazione % di Y / variazione % di X

#se X aumenta del 1% ->Y aumenta del 1.3%

plot(loglog, which=1)

loglin<-lm(I(log(peso))~bicipite,data2)

summary(loglin)

##variazione % di Y / variazione di X-> moltip 100

#se X aumenta di 1 cm ->Y aumenta del 4.07%

plot(loglin, which=1)

ES 3

data<-read.csv("nazioni.csv", sep=";")

library(psych)

numeric<-unlist(lapply(data,is.numeric))

Corr<-cor(data[,numeric])

require(ggm)

Corr.Par<-parcor(Corr)

require(corrplot)

par(mfrow=c(2,1))

corrplot(Corr)

corrplot(Corr.Par)

par(mfrow=c(1,1))

table(data$relig)

data$relig<-factor(data$relig)

# factor: fattorizza una variabile qualitativa in livelli, ciascuno dei quali diventa una variabile

# dicotomica autonoma.

m1<-lm(pil~.-nazione,data)

# il . indica tutte le variabili (come il * su SQL). Si possono escludere variabili del modello

# mettendole dopo il segno -

summary(m1)

require(car)

vif(m1)

# vif: calcola il variance inflation factor (il reciproco della tolleranza) per ciascun

# regressore del modello

m2<-lm(pil~.-nazione-vitamas,data)

summary(m2)

vif(m2)

data$relig <- relevel(data$relig,"Prot")

# relevel: rifattorizza una variabile divisa per livelli. Il livello di riferimento diventerà

# il livello corrispondente alla stringa data come argomento.

m3<-lm(pil~.-nazione-vitamas-relig,data )

anova(m3,m2)

m2<- lm(pil~densita + urbana + vitafem + alfabet + relig,data)

n<-length(m2$fitted.values)

k<-length(m2$coe_icients)

SSE<- sum(m2$residuals^2)

n*log(SSE/n)+2*k

extractAIC(m2)[2]

# extractAIC: restituisce due valori, il primo è il numero di gradi di libertà del modello,

# il secondo è l'e_ettivo valore dell'AIC

#AIC "classico"

#n*log(SSE/n)+2*(k+1)+n+n*log(2*pi)

#AIC(m2)

step(m2)

# step: procedura stepwise

extractAIC(update(m2, . ~ . - alfabet))[2]

# update: accetta in entrata due parametri, un modello e una formula. Restituisce il modello

aggiornato

# per la nuova formula. In questo caso, il . si riferisce non a tutte le variabili della distribuzione

# di dati, ma tutte le variabili presenti nella vecchia formula.

extractAIC(update(m2, . ~ . - urbana))[2]

extractAIC(update(m2, . ~ . - densita))[2]

extractAIC(update(m2, . ~ . - vitafem))[2]

extractAIC(update(m2, . ~ . - relig))[2]

m3<-update(m2, . ~ . - alfabet)

step(m2)

par(mfrow=c(1,1))

m4<-lm(formula = pil ~ vitafem + relig, data = data)

summary(m4)

plot(m4,which=1)

plot(m4,which=2)

plot(m4,which=3)

plot(m4, which=4)

plot(x=hatvalues(m1), y=rstudent(m1))

which(cooks.distance(m4) > 4/(nrow(data)-length(m4$coe_icients)-1))

which(abs(rstudent(m4))>2)

which(hatvalues(m4)>3*mean(hatvalues(m4)))

m4.2<-lm(formula = pil ~ vitafem + relig, data = data,

subset=(hatvalues(m4)<3*mean(hatvalues(m4)) & abs(rstudent(m4))<2) )

## l'argomento "subset=" accetta in entrata una a_ermazione logica (TRUE o FALSE) e

# la applica a tutte le righe, restituendo solo le righe che rispettano tale condizione.

new.data<-c(82, "Ort" )

names(new.data)<-c("vitafem","relig")

# names: assegna dei nomi alle caselle del vettore, invece di avere indice 1 e 2

new.data<-as.data.frame(t(new.data))

# as.data.frame: converte una struttura iterabile in un dataframe

new.data$vitafem<-as.numeric(new.data$vitafem)

predict(m4, newdata = new.data)

# predict: accetta in entrata un modello e una struttura di dati. Predice il valore della regressa

# (in questo caso "pil") date le osservazioni contenute nella struttura di dati. La struttura di

dati

# deve ovviamente contenere variabili chiamate come i regressori (nel nostro caso

"vitafem","relig")

m5<-lm(formula = I(log(pil)) ~ vitafem+ I(vitafem^2) + relig, data = data)

plot(m5,which=1)

library(leaps)

lmodc<-regsubsets(pil~densita + urbana + vitafem + alfabet + relig,data)

# regsubsets: sarò onesto capo, non ho 100% capito cosa sia, neanche 50% in realtà, ma so

che

# serve per valutare il modello più e_iciente (?)

lc<-summary(lmodc)

# la summary di una variabile risultata da regsubsets restituisce una variabile con vari

attributi.

# Quelli importanti sono "which", che ti dice nei modelli più e_icienti selezionati quali variabili

# vengono incluse e quali no, e "cp", che restituisce il Cp di Mallows del modello (state attenti

# quando cercate "Cp" su Google)

lc$which

lc$cp

# se da un modello con P variabili voglio tenerne K, un modello buono ha Cp pari a K

mallows.cp<-function(model, full.model ) {

var.full<- sum(full.model$residuals^2)/full.model$df.residual

dev.model<-sum(model$residuals^2)

dev.model/var.full-length(model$residuals)+2*length(model$coe_icients)

}

# dichiaro una funzione che determina il Cp di un modello. Assolutamente inutile dato che lo

fa

# già il summary di regsubsets, ma lo Spinellone voleva fare il ganassa

mallows.cp(lm(pil~ vitafem + relig,data) ,m2)

ES 4

setwd("H:\\Il mio Drive\\esercitazioni 2223")

setwd("C:\\Users\\DS064536\\Downloads")

data<-read.csv("studenti.csv", header=T)

summary(data)

table(data$anno)

library(psych)

pairs.panels(data[,c("altezza", "peso", "eta", "componenti")])

prop.table(table(data$fumo, data$genere))

# table: crea una tabella con le modalità di una variabile sulle righe e dell'alltra sulle colonne.

prop.table(table(data$fumo, data$genere), margin=1)

prop.table(table(data$fumo, d