Scarica Business Data Analytics e più Schemi e mappe concettuali in PDF di E-Business solo su Docsity!
PARTE 1
- Introduzione allo statistical learning
- Inferenza (test) per due popolazioni. One-way ANOVA per dati univariati.
- Analisi multivariata: esplorazione, quantificazione della dipendenza, covarianza e correlazione, matrice di varianza e covarianza.
- Statistical learning non supervisionato
- Riduzione dimensionale: Analisi delle Componenti Principali
- Classificazione non supervisionata: Clustering gerarchico e K-means
- Statistical learning supervisionato:
- Modelli parametrici di Regressione Lineare e Lineare Generalizzata
- Regressione lineare semplice e multipla. Stima dei coefficienti, valutazione dell’accuratezza del modello. Predittori qualitativi. PRESS.
- Selezione del modello e regolarizzazione: subset selection, metodi di penalizzazione, Ridge regression e LASSO.
- Regressione logistica binaria e multinomiale. Sensitività, Specificità, Curva ROC, AUC.
- Regressione e classificazione non parametriche o KNN, CART, Random Forest.
TEST DI IPOTESI
TEST DI IPOTESI DI 1 POPOLAZIONE
1) Verifica di normalità
Test di Shapiro (Per numerosità<5000) shapiro.test(col)
H0 : X~N vs H1: dati NON gaussiani P-value alto (>0.05): accetto H0, i dati provengono da una distribuzione normale. P-value basso: rifiutiamo l'ipotesi di normalità dei dati.
Graficamente t = seq(min(col), max(col), length.out = 1000) f = dnorm(t, mean = mu_hat, sd = sd_hat) # densità normale a cui dò in input t, media e sd quartz() hist(col, breaks = 15, probability = T) lines(t, f, type = "l", lwd = 3)
Confronto graficamente se i dati hanno una distribuzione simile alla rispettiva normale.
Graficamente quartz() qqnorm(col) qqline(col, lwd=2)
Confronto graficamente i quantili empirici con quelli teorici. In caso di normalità perfetta, i dati giacciono esattamente sulla retta.
2) Formulazione ipotesi e calcolo p-value
CASO 1 : popolazione distribuita normalmente: t-test Uso funzione t.test t.test(dati$Colesterolo, mu = 200, alternative = "less") UNILATERO A SINISTRA t.test(dati$Colesterolo, mu = 200, alternative = "greater") UNILATERO A DESTRA t.test(dati$Colesterolo, mu = 200, alternative = "two.sided") BILATERO H0 test bilatero: mu=200. Se p-value alto, accetto H0, quindi mu è circa uguale a 200.
Oppure a mano: Test unilatero a sinistra (H0: mu >=200 vs H1: mu<200) mu_0 = 200 n = length(dati$Colesterolo) media_campionaria = mean(dati$Colesterolo) dev_std_campionaria = sd(dati$Colesterolo) t = (media_campionaria - mu_0)*sqrt(n)/dev_std_campionaria pval = pt(t, df = n-1) pval
P-value alto: accetto H0. P-value basso: rifiuto H
TEST DI IPOTESI DI 2 POPOLAZIONI
1) Verifica di normalità
Test di Shapiro (Per numerosità<5000) shapiro.test(colesterolo_uomini) shapiro.test(colesterolo_donne)
H0 : X~N vs H1: dati NON gaussiani P-value alto: accetto H0, i dati provengono da una distribuzione distribuita normalmente. P-value basso: rifiutiamo l'ipotesi di normalità dei dati.
2) Formulazione ipotesi e calcolo p-value
Caso a varianze non note Le VARIANZE nelle due sottopopolazioni sono UGUALI? O DIVERSE? Esistono dei testi approssimati per l'uguaglianza tra varianze di due campioni, come ad esempio il test di Fisher ( basta saper leggere il p-value ) var.test(colesterolo_uomini, colesterolo_donne) t.test(colesterolo_uomini,colesterolo_donne, var.equal = T) Oppure a mano: X = colesterolo_donne Y = colesterolo_uomini nX = length(X) nY = length(Y) s2x = var(X) s2Y = var(Y) s2p = 1/(nX + nY - 2) * ((nX-1)s2x + (nY-1)s2Y) sp = sqrt(s2p) t = (mean(X) - mean(Y)) / sp / sqrt (1/nX + 1/nY) pval = 2*(1 - pt(t, df = nX + nY -2)) Risultati coerenti pval
P-value alto: accetto H0, le due medie sono uguali P-value basso: rifiuto H0, le due medie sono diverse
Caso a varianze note var_M = 150 var_F = 110 differenza_medie = mean(colesterolo_uomini) - mean(colesterolo_donne) n_uomini = length(colesterolo_uomini) n_donne = length(colesterolo_donne) z = differenza_medie / sqrt(var_M/n_uomini + var_F/n_donne) pvalue = 2*(1-pnorm(abs(z))) pvalue
P-value alto: accetto H0, le due medie sono uguali P-value basso: rifiuto H0, le due medie sono diverse
ANOVA TEST (Per g>2 popolazioni)
Verifica che tutti i g gruppi abbiano la stessa media (H0), altrimenti ne esiste almeno uno con media diversa (H1).
1) Verifica bilanciamento classi
Specie in presenza di pochi dati, è bene accertarsi che le classi siano ragionevolmente bilanciate. table(dati_polli$feed)
2) Verifica normalità
tapply(dati_polli$weight, dati_polli$feed, shapiro.test) La funzione tapply suddivide i valori di una variabile quantitativa per i livelli di una variabile qualitativa e, per ogni livello, applica una particolare funzione (ad esempio: mean, max, shapiro.test,...) sui dati. Per ogni sottogruppo verifico che la distribuzione sia normale (p-value>0.05)
3) Verifica di omoschedasticità
Ovvero che le varianze in ogni gruppo possano essere assunte uguali. Test di Bartlett bartlett.test(dati_polli$weight ~ dati_polli$feed) p-value alto: tutti i gruppi hanno varianza simile p-value basso: almeno uno dei gruppi ha varianza diversa dagli altri.
4) ANOVA
Ricordiamo la parametrizzazione mu_j = mu + tau_j , con tau_1 = 0 Quindi il modello è X_ji = mu + tau_j + epsilon_ij, epsilon_ij ~ N(0, sigma^2) H0: tau_j = 0 per ogni j VS H1: esiste j tale che tau_j!=
Anova_polli = aov(weight ~ feed, data = dati_polli) summary(Anova_polli) Numero di livelli del fattore feed: numero di gruppi: g Df del feed: gradi di libertà della variabile feed: g- Dimensione del dataset: numero di osservazioni: n Df dei Residuals: gradi di libertà dei residui: n-g Il p-value è indicato nel summary, nella cella (feed; Pr(>F)). p-value alto: accetto H0 (e la variabile categorica feed non dà distinzioni tra i gruppi). p-value basso: rifiuto H0.
Stima della comune varianza S , cioè la varianza stimata per tutte le popolazioni. S = sum(Anova_polli$residuals ^ 2) / Anova_polli$df.residual S
Anova_polli$df.residuals sono i gradi di libertà dei residui.
Intervalli di confidenza alpha = 0.05 (Intervallo di confidenza al 95%)
source("Confronti_multipli_da_anova.R") intervalli_anova = Confronti_multipli_da_anova(Anova_polli, alpha = alpha) names(intervalli_anova) intervalli = intervalli_anova$intervals
PCA
Riduce il numero di variabili che descrivono un insieme di dati, limitando il più possibile la perdita di informazioni. library(GGally) library(ggfortify)
1) Isolo componenti numeriche
turisti_obs = turisti[ ,-c(1,2)]
2) Verifica correlazione tra le variabili
Se c’è correlazione forte tra le variabili, allora ha senso fare PCA. Graficamente library(GGally) quartz() ggpairs(turisti_obs) Con matrice di correlazione cor(turisti_obs) Se tutti gli elementi della matrice sono vicini a 1, allora la struttura di correlazione è abbastanza forte, quindi c’è omogeneità, si può fare PCA.
3) Verifica covarianza tra le variabili
Matrice di covarianza cov(turisti_obs) Se le varianze assumono valori molto diversi, allora c’è forte disomogeneità.
4) Se c’è forte disomogeneità => Standardizzazione dei dati
Se c’è forte disomogeneità, data da matrice di correlazione e/o di covarianza, allora bisogna standardizzare il dataframe (ogni dato avrà media=0 e varianza=1). turisti_obs_scaled = data.frame(scale(turisti_obs)) Si nota che così matrice di correlazione e di covarianza sono uguali.
5) Calcolo componenti principali PC
I principal component sono i vettori lungo cui si sviluppano maggiormente i dati, estrapolati dal campione X, costituito da n osservazioni distribuite normalmente con media mu e matrice di covarianza sig. PC = princomp(X) Deviazioni standard di ogni PC var_props = (PC$sdev^2) / sum(PC$sdev^2) Percentuale di varianza (percentuale di variabilità) di ogni PC cumsum(var_props) Restituisce la percentuale di varianza cumulata summary(PC) PC, var_props, cumsum(var_props)
Loadings : quanto quella variabile influisce su quel principal component. PC$loadings Dove non c'è valore vuol dire che è 0 Scores : valore dei dati proiettati lungo le 3 direzioni principali PC$scores
6) Scree Plot
quartz() plot(cumsum(PC$sd^2)/sum(PC$sd^2), type='b', axes=F, xlab='number of components', ylab='contribution to the total variance', ylim=c(0,1), lwd=2) abline(h=1, col='red', lwd=1.5) abline(h= x , lty=2, col='blue') box() axis(2,at=0:10/10,labels=0:10/10) axis(1,at=1:ncol(turisti),labels=1:ncol(turisti),las=2)
Posso così visualizzare il contributo di ogni componente alla variabilità totale. Posso ottenere così il numero minimo di componenti principali da considerare per spiegare almeno x% della variabilità del dataset (prendo le componenti fino a quella che supera la linea blu).
Varianze in ordine decrescente variances = diag(cov(turisti_obs)) sort(variances, decreasing = TRUE)
7) Altre visualizzazioni grafiche
library(ggfortify) autoplot(PC_tour_scaled) Visualizza gli scores sul piano PC1-PC2, cioè tutti i dati
Biplot : autoplot(PC_tour_scaled, loadings = TRUE, loadings.label = TRUE) Mostra sul piano PC1-PC2 come ogni variabile influisce su PC1 e PC
Per specificare quale PC vogliamo sugli assi: autoplot(PC_tour_scaled, x = 2, y = 3, loadings = TRUE, loadings.label = TRUE)
Grafico loadings delle prime 3 PC par(mar = c(1,4,0,2), mfrow = c(3,1)) for(i in 1:3)barplot(load.food[,i], ylim = c(-1, 1))
8) Aggiunta di variabili categoriche
Costruisco un dataframe con gli scores, che ha come colonne le PC e, in più, le variabili categoriche. Scores = data.frame(PC_tour_scaled$scores) names(Scores) = paste0("Comp",1:8) Scores = cbind(Scores,turisti_lab)
Boxplot di PC1 rispetto alla variabile categorica dei mesi. quartz() boxplot(Scores$Comp1 ~ Scores$Month)
4) Selezione di K
Grafico della within sum of squares (tot.withinss) tot.withinss = NULL for(k in 1:10) tot.withinss = c(tot.withinss, kmeans(nidi, k, nstart=20)$tot.withinss)
quartz() plot(1:10, tot.withinss, type='b', xlab='K') Si cerca un gomito in corrispodenza del K dopo cui il guadagno è trascurabile.
Si noti che a volte, si preferisce guardare al grafico della tot.withinss normalizzandolo rispetto alla total sum of squares ovvero la variabilità totale. Si ottiene una curva del tutto equivalente ai fini interpretativi, che ha però il pregio di assumere valori tra 0 e 1: la si può interpretare come frazione della variabilità totale dovuta alla variabilità intra-cluster. ratio = NULL for(i in 1:10){ kmeans_out = kmeans(data,i,nstart=20) ratio = c(ratio, kmeans_out$tot.withinss/kmeans_out$totss) } quartz() plot(1:10, ratio, type='b', xlab='K', ylab = "tot.withinss/totss", lwd=2) L'andamento della curva è esattamente identico: ogni suo punto è semplicemente stato diviso per una costante, la variabilità totale. Ricordiamo infatti che, mentre sia betweenss che tot.withinss sono proprietà della clusterizzazione, la loro somma, cioè la variabilità totale, dipende solo dai dati. Questa curva può essere interpretata come frazione della variabilità totale dovuta alla variabilità intra-cluster.
Grafico della silhouette per k = 3 e k = 2 Si ricorda che:
- Silhouette di circa 1 significa che i cluster sono molto densi e ben separati
- Silhouette di circa 0 significa che i cluster si sovrappongono
- Silhouette minore di 0 significa che i cluster potrebbero essere sbagliati/non corretti
quartz() par(mfrow = c(2,2)) for ( i in 2:5){ km.output = kmeans(data, centers = i, nstart = 20) silhouettekm = silhouette(km.output$cluster , dist(data)) plot(silhouettekm, col = 1:i, main=paste("Silhouette per k =", i, sep=" ")) abline(v = mean(silhouettekm[,3]), lwd=3, col='red4') }
Devo prendere il K che ha average silhouette width più vicina a 1.
CLUSTERING GERARCHICO
library(cluster)
1) Definisco dist
Il comando dist calcola una matrice metrica (triangolare inferiore), in forma vettoriale per evitare ridondanze, seguendo una distanza specificata dall'utente. Il default è la distanza euclidea.
metric = dist(pioppi) Sarebbe equivalente a metric = dist(pioppi, method="euclidean") Altrimenti si può anche usare, ad esempio, il metodo manhattan.
dist(pioppi[1:7, ]) R non rappresenta la matrice metrica in forma matriciale, per evitare di ripetere elementi che sono uguali per definizione, ma questo non è un problema, poichè utilizzeremo tale oggetto passandolo direttamente agli algoritmi di clustering gerarchico.
2) hclust
Clustering gerarchico a partire da una matrice metrica: comando hclust. Occorre anche passare il metodo di linkage desiderato: Complete Linkage: distanza tra cluster è distanza tra i punti più lontani hc.complete = hclust(metric, method="complete") Average Linkage: distanza tra cluster è distanza tra i centroidi hc.average = hclust(metric, method="average") Single Linkage: distanza tra cluster è distanza tra i punti più vicini hc.single = hclust(metric, method="single")
3) Dendrogrammi
quartz() par(mfrow=c(1,3)) plot(hc.complete, main="Complete Linkage", xlab="", sub="", cex=.9, labels = F) plot(hc.average, main="Average Linkage", xlab="", sub="", cex=.9, labels = F) plot(hc.single, main="Single Linkage", xlab="", sub="", cex=.9, labels = F)
Mostra come funziona graficamente il clustering gerarchico: in basso ci sono i singoli dati che vengono raggruppati a due a due, a loro volta raggruppati e così via. Le zone dove sono più lunghi i rettangoli sono zone di stabilità. Il numero di gruppi che si vogliono identificare è tendenzialmente in una zona di stabilità molto alta.
Si deve scegliere il tipo di clustering con zona di stabilità più alta.
MODELLI LINEARI
library(GGally) library(MASS) library(car) library(rgl)
REGRESSIONE LINEARE SEMPLICE
1) lm
first_model = lm(Income ~ Education, data = data) summary(first_model) Il regressore è Education.
2) Graficamente
quartz() par(mfrow = c(2,2)) plot(first_model)
3) Info che si possono ottenere da lm
names(first_model) first_model$coefficients Valori stimati dei beta first_model$residuals Residui : epsilon_hat first_model$rank p + 1 first_model$fitted.values Valori fittati: y_hat first_model$df.residual Gradi di liberta' del modello: n - p - 1 first_model$model Il data frame dei dati sigma2_hat = sum(first_model$residuals^2)/first_model$df.residual Stima della varianza
vcov(first_model) Matrice di covarianza dei beta rstandard(first_model) Residui standardizzati
4) Intervalli di confidenza per i beta
alpha = 0. confint(first_model, level = 1 - alpha)
5) Correggere per la nostra dimensionalità (p+1)
confint(first_model, level = 1 - alpha/2) In generale: 1 - alpha/(p+1)
6) Intervalli di confidenza e predizione per una nuova osservazione x
x_0 = 16. Intervallo per la media: predict(first_model, data.frame(Education = x_0), interval = "confidence", level = 1 - alpha)
Intervallo di predizione: predict(first_model, data.frame(Education = x_0), interval = "prediction", level = 1- alpha)
7) Riportiamo graficamente gli intervalli puntuali di predizione e confidenza
(point-wise intervals!)
X0 = data.frame(Education = seq(10, 22, length=100)) Creiamo 100 "nuovi" punti e interpoliamo Conf = predict(first_model, X0, interval='confidence') Pred = predict(first_model, X0, interval='prediction')
quartz() plot(data$Education, data$Income, pch=19, xlab = 'Education', ylab = 'Income') lines(X0[,1], Conf[,'fit'], lwd=1.5) lines(X0[,1], Conf[,'lwr'], lty=2, col='red', lwd=3) lines(X0[,1], Conf[,'upr'], lty=2, col='red', lwd=3)
lines(X0[,1], Pred[,'lwr'], lty=3, col='deepskyblue', lwd=3) lines(X0[,1], Pred[,'upr'], lty=3, col='deepskyblue', lwd=3)
legend('topleft', legend = c('Fit', 'Confidence Interval', 'Prediction Interval'), lwd=2, lty = c(1,2,3) , col = c('black', 'red', 'deepskyblue') )
REGRESSIONE LINEARE MULTIPLA
library(GGally) library(MASS) library(car) library(rgl) Adesso facciamo regressione lineare multipla mettendo anche un altro regressore: l'anzianità. Quindi il modello e' y_i = beta_0 + beta_1Education + beta_2Seniority + eps_i
1) lm
second_model = lm(Income ~ Education + Seniority, data = data) summary(second_model)
I regressori sono Education e Seniority. Ha senso aggiungere un regressore se il p-value relativo alla sua significatività (Pr(>|t|)) è basso. I regressori con (Pr(>|t|)) alto non sono significativi: si possono togliere.
2) Controllo la normalità dei residui
shapiro.test(second_model$residuals)
3) Diagnostica del modello
quartz() par(mfrow = c(2,2)) plot(second_model)
Controllo che omoschedasticità sia ok e che gli effetti di leva non siano forti. Bisogna sempre fare diagnostica sui residui e sulla retta di regressione, non limitarsi a guardare SOLO R2.
- Grafico: quartz() plot(speed2, dist, main='Scatterplot Speed^2 vs. Dist', lwd=2, xlab='Speed2', ylab='dist') matplot(X.new, IC, add=T, type='l', col=c('black','red','red'), lwd=2, lty=2) matplot(X.new, IP, add=T, type='l', col=c('black','blue','blue'), lwd=2, lty=2) legend('topleft', legend=c('Regression line','Confidence intervals','Prediction intervals'), col=c('black','red','blue'), lwd=2, lty=c(1,2,2), cex=1.2)
MODELLO DI REGRESSIONE LINEARE MULTIPLA CON UN PREDITTORE
CATEGORICO
Modello: Average_Score = beta_0 + beta_1 * Years_Service + Eps, Eps ~ N (0, sigma^2) Si vuole cioè includere una variabile categorica nel modello lineare.
1) lm
result = lm(Average_Score ~ Years_Service, data = work) summary(result)
2) Esistono graficamente differenze tra gli elementi di una variabile categorica?
work$Sex = as.factor(work$Sex) Importante farlo: la categorica deve diventare di tipo factor str(work) quartz() plot(work$Years_Service, work$Average_Score, main='Scatterplot di Y vs X', lwd=2, pch=19, xlab='Years of Service', ylab='Average Score', col = ifelse(work$Sex == "Male", "blue","red"))
3) Definizione delle dummy variables
Costruiamo una variabile dummy , ovvero una variabile binaria (solo valori 0 o 1) che codifica l'appartenenza ad una categoria. La variabile avrà infatti valore: 1 se l'osservazione è della categoria maschi, 0 altrimenti (femmine) - categoria di riferimento E' chiaro quindi che per una categorica con K livelli ci serviranno sempre K-1 variabili dummies.
Definisco la dummy e la inserisco nel dataset: work$dummy = ifelse(work$Sex == "Male",1,0) work$dummy
Il nostro modello con l'inclusione della dummy avrà la seguente forma: Average_Score = beta0 + beta1dummy + beta2Years_Service + beta3dummyYears_Service + eps
Notate che possiamo raccogliere Years_Service: Average_Score = beta0 + beta1dummy + (beta2 + beta3dummy)*Years_Service + eps Ovvero:
Average_Score = B0[g] + B1[g]*Years_Service + eps with B0[g]=beta0 se l'unità è nel gruppo per cui dummy=0 (females) B0[g]=beta0+beta1 se l'unità è nel gruppo per cui dummy=1 (males) B1[g]=beta2 se l'unità è nel gruppo per cui dummy=0 (females) B1[g]=beta2+beta3 se l'unità è nel gruppo per cui dummy=1 (males)
Fittiamo il modello: model = lm(Average_Score ~ dummy + Years_Service + dummy:Years_Service, data = work) ) Basta fare questo: model = lm(Average_Score ~ Sex + Years_Service + Sex:Years_Service, data = work) summary(model)
La dummy è rilevante se ha p-value basso. L’effetto della dummy sulla pendenza della retta si determina dal p-value di beta (dummy:Years_Service). Se è basso allora ha effetto rilevante. Se è alto allora si può togliere.
Eliminazione di una variabile non rilevante model = lm(Average_Score ~ dummy + Years_Service, data = work) summary(model)
Interpretazione del modello: modello per le donne: Y = 7.035 + 0.097 X + eps modello per gli uomini: Y = 7.035 - 2.59 + 0.097 X + eps = 4.44 + 0.097 X + eps
4) Costruzione grafica
coefficienti = model$coefficients beta_0_donne = coefficienti[1] beta_1_donne = beta_1_uomini = coefficienti[3] beta_0_uomini = coefficienti[1] + coefficienti[2]
5) Diagnostica dei residui
quartz(width = 14, height = 8) par(mfrow=c(2,2)) plot(model)
shapiro.test(model$residuals)
MODELLO DI REGRESSIONE LINEARE MULTIPLA CON 2 PREDITTORI
CATEGORICI
Guarda 6.0, Esempio 2.
STEPWISE SELECTION
1) lm senza variabili categoriche + diagnosi dei residui
first = lm(depth ~ lat + long, data = data) summary(first)
quartz() par(mfrow = c(2,2)) plot(first)
shapiro.test(first$residuals)
2) Aggiungiamo variabile categorica
data$zone = as.factor(data$zone) second = lm(depth ~ lat + long + zone + zone:lat + zone:long, data = data) summary(second)
Diagnosi dei residui quartz() par(mfrow = c(2,2)) plot(second)
3) Modello polinomiale (molto ricco) + diagnosi residui
fourth = lm(depth ~ lat + long + I(lat^2) + I(long^2) + I(latlong) + zone + zone:lat + zone:long + zone:I(lat^2) + zone:I(long^2) + zone : I(latlong) , data = data) summary(fourth)
shapiro.test(fourth$residuals) quartz(width = 14, height = 8) par(mfrow = c(2,2)) plot(fourth)
Considerazioni:
- modello lineare si adatta bene anche per fittare due paraboloidi
- attenzione: non si faccia inferenza (quindi previsione) in punti dello spazio in cui non ci sono osservazioni
- summary complicati all'aumentare dei predittori: come trovare il giusto compromesso?
4) StepAIC
Il processo di selezione e riduzione del modello puo' essere reso automatico con alcuni algoritmi! Utilizziamo stepAIC che procede con la stepwise selection nelle sue tre forme, a seconda dell'opzione direction specificata, e restituisce il modello che, in termini di AIC, non poteva essere ulteriormente migliorato (la regola è quella di preferire i modelli con l'AIC più basso). Si ricorda agli studenti che tra gli algoritmi di scelta delle variabili abbiamo:
- Backward elimination (
"backward"): elimina predittori in sequenza dal modello iniziale. Produce una sequenza di modelli di complessita' decrescente fino ad ottenere il modello ottimo. - Forward selection (
"forward"): aggiunge predittori in sequenza, usando quelli disponibili nell'argomento data di lm. Produce una sequenza di modelli di complessità crescente fino a raggiungere il modello ottimo. - Stepwise regression (
"both"): una ricerca forward-backward che, ad ogni step, decide se includere o escludere un predittore. A differenza delle modalità precedenti, un predittore che era stato incluso/escluso in uno step precedente, può essere escluso/incluso in uno step successivo.
Stepwise forward selection Definisco il modello minimo con solo l'intercetta, passando però tutto il dataset alla funzione questo permetterà a stepAIC di sapere dove prendere le variabili da aggiungere fino a raggiungere l'upper bound. miminal_model = lm(depth ~ 1, data = data)
forward = stepAIC(miminal_model, direction="forward", scope=list(lower=miminal_model, upper=fourth)) forward$anova summary(forward) quartz(width = 14, height = 8) par(mfrow = c(2,2)) plot(forward)
Parte dal modello base aggiungendo tutti i termini presenti in fourth tranne:
- l'interazione tra zona sismica e quadrato della latitudine (zone:I(lat^2))
- prodotto tra latitudine e longitudine (lat:long)
Stepwise backward selection backward = stepAIC(fourth, details = T, direction = "backward") backward$anova summary(backward) quartz(width = 14, height = 8) par(mfrow = c(2,2)) plot(backward)
Rimuove dal modello fourth:
- l'interazione tra zona sismica e quadrato della latitudine (I(lat^2):zone)
- prodotto tra latitudine e longitudine (lat:long)
Stepwise mixed selection both = stepAIC(fourth, details = T, direction = "both") both$anova summary(both) quartz(width = 14, height = 8) par(mfrow = c(2,2)) plot(both)
