Obiettivo dell'analisi: prevedere qualcosa
A seconda del contesto e degli obiettivi, ci interesserà avere un modello che
Valutazione delle variabili presenti nel dataset ed eventuale pulizia dello stesso. In particolare si dovrà
Selezionare le variabili “sensate”: andranno rimosse le variabili non informative, legate ad un identificativo, vuote o delle quali non si conosce il significato etc.
Ricodifica delle variabili:
Trattamento/rimozione di dati mancanti
Rimozione variabili concomitanti o leaker
perc = 2/3
set.seed(42)
indici <- sample(1:nrow(data), round(nrow(data)*perc))
train <- data[indici, ]
test <- data[-indici, ]
# Variabili risposta
# y.train <-
# y.test <-
#
# x.train <- train[, -which(names(train) == "")]
# x.test <- test[, -which(names(test) == " ")]
#
mat.train <- model.matrix(~. , data = x.train)[, -1]
mat.test <- model.matrix(~., data = x.test)[, -1]
# Trasformata
# Logaritmica: attezione agli zeri
y.train.log <- log(y.train)
y.test.log <- log(y.test)
# Matrice per il confronto tra modelli
confronto <- matrix(NA, nrow = 1 , ncol =4,
dimnames = list(c(), c("model", "est. scale", "err.original", "err.log")))
err <- function(prediction) errori(prediction, actual= y.test, log = T, actual.log= y.test.log)
#--------------------------------------#
# Divisione stima/verifica per selezione modello all'interno del training test
set.seed(42)
cb1 <- sample(1:NROW(train), floor(NROW(train)/2))
cb2 <- setdiff(1:NROW(train), cb1)
\[ \hat{y} = Sy \] \[ y_i - \hat{y}_{-i} = \frac{(y_i - \hat{y}_i)}{1 - S_{ii}} \]
Approssimazione tramite
\[ GCV = \frac{1}{n} \sum_{i = 1}^n [\frac{(y_i - \hat{f}(x_i)}{1 - tr(S)/n}]^2 \]
dove \( tr(S) \) approsima il numero effettivo di parametri.
Strategia generale: in linea di massima i dati a disposizione sono abbastanza da permettere la suddivisione in training/test set; nella fase di stima dei vari modelli scelgo il migliore per tipo con uno dei criteri
Nel caso ci sia un numero elevato di variabili esplicative, conviene selezionarne un sottoinsieme, soprattuto se devo stimare modelli il cui costo computazionale è strettamente connesso al numero di variabili. Ci sono diverse possibilità.
fit.lm <- lm(y.train~. ,data = x.train)
summary(fit.lm)
##
## Call:
## lm(formula = y.train ~ ., data = x.train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -35.30 -8.87 -0.73 7.65 303.78
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.47e+02 1.87e+00 78.58 < 2e-16 ***
## date -1.51e-03 1.68e-04 -8.99 < 2e-16 ***
## dowMonday 1.20e+00 9.20e-01 1.30 0.193
## dowSaturday 1.20e+00 9.16e-01 1.31 0.191
## dowSunday -1.55e+00 9.42e-01 -1.65 0.099 .
## dowThursday -4.63e-01 9.06e-01 -0.51 0.609
## dowTuesday 2.47e-01 9.19e-01 0.27 0.788
## dowWednesday -5.34e-01 9.08e-01 -0.59 0.557
## pm10tmean 1.04e-01 1.68e-02 6.18 7.4e-10 ***
## o3tmean 6.15e-02 3.18e-02 1.94 0.053 .
## so2tmean 2.21e-01 1.07e-01 2.06 0.040 *
## no2tmean -3.21e-04 5.33e-02 -0.01 0.995
## cotmean -8.20e-04 1.21e-03 -0.68 0.500
## tmpd -4.21e-01 4.49e-02 -9.38 < 2e-16 ***
## dptp 8.78e-02 4.36e-02 2.01 0.044 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.8 on 3212 degrees of freedom
## Multiple R-squared: 0.182, Adjusted R-squared: 0.178
## F-statistic: 50.9 on 14 and 3212 DF, p-value: <2e-16
pred.lm <- predict(fit.lm, newdata = x.test)
# Nel caso non siano ammesse previsioni negative, o altro determino un valore t
# t = 0.5
# pred.lm <- pmax(predict(fit.lm, newdata = x.test), t)
confronto[1, ] <- cbind("linear", "original", err(pred.lm))
# Considero le variabili significative nel primo modello di regressione
# lineare e andare avanti con quelle
p.val <- summary(fit.lm)$coefficients[,4]
names(which(p.val < 0.05))
# Tengo conto delle variabili qualitative a "mano"
# è plausibile che ci possano essere solo alcuni livelli di tali variabili significativi, ma i test per le singole variabili sono fatti in relazione ad un livello di base.
v.qual <-
v.quant <- names(train)[names(train) %in% names(which(p.val < 0.05))]
library(leaps)
best.sub <- leaps(x = mat.train, y=y.train)
id <- best.sub$which[which.min(best.sub$Cp),]
fit.best <- lm(y.train~., data=data.frame(mat.train[,id]))
summary(fit.best)
#Previsione
pred.best <- predict(fit.best, newdata=data.frame(mat.test[,id]))
plot(best.sub$size, best.sub$Cp,xlab="Numero di parametri", ylab="Cp")
abline(v=best.sub$size[which.min(best.sub$Cp)], lty="dashed")
confronto <- rbind(confronto, cbind("best.sub", "original", err(pred.best)))
library(lars)
# Scelta del modello
# Criterio di "penalizzazione": simile all'AIC, la procedura mi fornisce la statistica Cp che penalizza
# in base al numero di parametri
step.models <- lars(x = mat.train, y = y.train, type = "stepwise")
pred.step <-predict.lars(step.models, newx= mat.test)
plot(step.models$Cp)
J <- step.models$df[step.models$Cp==min(step.models$Cp)]
# Secondo migliore
K <- step.models$df[step.models$Cp==min(step.models$Cp[-J])]
# Da valutare anche ad occhio
pred.step <- pred.step$fit[,J]
confronto <- rbind(confronto, c("step","original" , err(pred.step)))
step.cv <- cv.lars(x = mat.train, y = y.train, type = "stepwise")
# Minimo
min.step.cv <- step.cv$index[which.min(step.cv$cv)]
# Regola un errore standard
one.se <- step.cv$index[which.min(abs(step.cv$cv[-min.step.cv] - step.cv$cv[min.step.cv] + step.cv$cv.error[min.step.cv]))]
fit.step <- lars(y.train, x=mat.train, type="step")
#Previsione
pred.step <- predict(fit.step, newx=mat.test, s= one.se, mode="step")$fit
confronto <- rbind(confronto, c("step","original" , err(pred.step)))
# Stime dei coefficienti entrati nel modello
coef(fit.step)[one.se, ]
# Stime dei coefficienti entrati nel modello
# coef(fit.step)[min.step.cv, ]
#Rappresentazione grafica per la scelta del parametro lambda
# plot(prova.step ,xvar="df",plottype="Cp");plot(prova.step)
# Tramite CV
# Non ne vale la pena
# library(foreach)
# library(doMC)
# registerDoMC(cores = 2)
library(glmnet)
# Cv sulla base del mean squared error
glmnet.cv <- cv.glmnet(mat.train, y.train , nfolds = 10, alpha = 1,
family = "gaussian")
plot(glmnet.cv)
# Lambda con minimo errore
glmnet.cv$lambda.min
# lambda con regole di un err. standard
glmnet.cv$lambda.1se
pred.glmnet1 <- predict(glmnet.cv, newx= mat.test, s="lambda.min")
pred.glmnet2 <- predict(glmnet.cv, newx=mat.test, s="lambda.1se")
confronto <- rbind(confronto, c("lasso.min","original" , err(pred.glmnet1)))
confronto <- rbind(confronto, c("lasso.1se","original" , err(pred.glmnet2)))
f.lasso <- lars(mat.train[cb1, ], y.train[cb1 ], type = "lasso")
p.lasso <- predict.lars(f.lasso, newx = mat.train[cb2, ] )
matrice.err<-tasso.lars(p.lasso$fit, y.train[cb2])
plot(matrice.err[, 2], type = "l")
# min
J <- matrice.err[matrice.err[,2]==min(matrice.err[,2], na.rm=T),1]
# oppure scegli ad occhio
fit.lasso <- lars(mat.train, y.train, type = "lasso", max.steps = J)
pred.lasso <- predict.lars(fit.lasso, newx = mat.test)$fit[, J]
# Coefficienti
beta.lasso <- coef(fit.lasso)[J,]
# coefficienti a zero
varname<-attr(mat.train, "dimnames")[[2]]
varname[beta.lasso==0]
# e quali no
varname[beta.lasso!=0]
confronto <- rbind(confronto, c("lasso(tasso)","original" , err(pred.lasso)))
Selezione variabili per modello GAM: 1. Se ho poche variabili posso usare una procedura di tipo stepwise
library(gam)
fit.gam <- step.gam(gam(y.train~1, data=x.train), scope=gam.scope(x.train))
plot(fit.gam)
pred.gam<-predict(fit.gam, x.test)
confronto <- rbind(confronto, c("gam","original" , err(pred.gam)))
fit.gam.sel <- gam(as.formula(paste("VAR_RISPOSTA ~", paste(paste("s(", paste(v.quant, collapse=")+s("),")+"),
paste(v.qual, collapse="+"), collapse = "+"))), data = train)
plot(fit.gam.sel, ask = T)
pred.gam.sel<-predict(fit.gam.sel, x.test)
confronto <- rbind(confronto, c("gam.sel","original" , err(pred.gam.sel)))
library(polspline)
# Indicare gli indici delle variabili categoriali!
# Se uso una selezione mi basta
# which(names(x.train) %in% v.qual)
fit.mars <- polymars(y.train, x.train, factors = )
tab <- fit.mars$model
tab$pred1 <- c("cost", names(x.train))[tab$pred1 + 1]
tab
## pred1 knot1 pred2 knot2 coefs SE
## 1 cost NA 0 NA 1.085e+02 5.716e+00
## 2 tmpd NA 0 NA 2.926e-01 1.223e-01
## 3 tmpd 71.5 0 NA -1.002e-01 2.670e-01
## 4 date NA 0 NA 2.529e-03 6.535e-04
## 5 pm10tmean NA 0 NA 1.149e-01 1.357e-02
## 6 date 8969.0 0 NA 1.615e-02 5.710e-03
## 7 tmpd 31.0 0 NA -3.397e-01 7.053e-02
## 8 date NA 8 NA -5.083e-05 1.246e-05
## 9 date 9677.0 0 NA 2.873e-02 4.380e-03
## 10 date 9212.0 0 NA -4.071e-02 8.418e-03
## 11 tmpd 80.5 0 NA -4.818e+00 2.589e+00
## 12 date NA 8 80.5 9.745e-04 3.261e-04
## 13 tmpd 62.5 0 NA -1.896e-02 4.356e-01
## 14 date 10613.0 0 NA -1.092e-02 3.373e-03
## 15 date NA 8 62.5 5.363e-05 4.617e-05
pred.mars <- predict(fit.mars, x.test)
confronto <- rbind(confronto, c("mars","original" , err(pred.mars)))
# Divisione in due insiemi
library(tree)
fit.trees <- tree(VAR_RISP ~. , data=train[cb1,],
control=tree.control(nobs=length(cb1), minsize=2, mindev=0.001))
prune <- prune.tree(fit.trees, newdata=train[cb2 ,])
plot(prune)
J <- prune$size[prune$dev == min(prune$dev)]
J
fit.tree <- prune.tree(fit.trees, best = J)
plot(fit.tree)
text(fit.tree, cex = 0.6)
pred.tree <- predict(fit.tree, newdata = x.test)
confronto <- rbind(confronto, c("tree","original" , err(pred.tree)))
library(rpart)
fit.rpart <- rpart(y.train ~., data = x.train, control = rpart.control(xval = 10))
plot(fit.rpart)
text(fit.rpart, cex = 0.7)
pred.rpart <- predict(fit.rpart, x.test)
confronto <- rbind(confronto, c("tree(rpart)","original" , err(pred.rpart)))
Per le reti neurali ho scritto una funzione in cui vengono provati diversi valori per il weigth-decay e numero di strati latenti. L'insieme di training viene suddiviso ulteriormente in due parti, la prima per la stima, la secondo per il calcolo dell'errore quadratico medio (in scala originale e versione logaritmica). Viene restituita una matrice contenente, per ogni modello stimato, gli errori nelle due scale, il valore del weight-decay, il numero di strati latenti e un indicatore di covergenza del processo di stima.
library(nnet)
decay<- 10^(seq(-4, -2, length=5))
layers <- seq(1, 3, by = 1)
err.mat <- neural.reg(formula = as.formula(VAR_RISP)~. , response = y.train, data= train,
index.train= cb1, index.test = cb2,
decay = decay, layers = layers , maxit = 500)
# Rete con errore minore
net.or <- err.mat[which.min(err.mat[, 1]), ]
fit.net1 <- nnet(y.train ~. , data =x.train, size = as.numeric(net.or[4]) ,
decay = as.numeric(net.or[3]), maxit = 1000, linout = T)
## # weights: 17
## initial value 44438557.373929
## final value 747830.372322
## converged
pred1 <- predict(fit.net1, x.test)
confronto <- rbind(confronto, c("nnet1", "original", err(pred1)))
Procedura bootstrap, in cui per ogni campione viene stimato il modello. Bisogna decidere il numero di replicazioni bootstrap. La media di predittori “indipendenti” porta ad una riduzione in termini di varianza per le previsioni ottenute (legge dei grandi numeri); un alto numero di replicazioni quindi non dovrebbe comportare problemi di sovrdattamento, ma comporta un maggiore costo computazionale.
Funziona bene per procedure instabili e può portare ad un notevole miglioramento in termini di errore quadratico medio; per procedure stabili è inutile (e.g. regressione lineare)
Di default vengono usati gli alberi di classificazione, ma in generale di possono considerare altri modelli, e si potrebbe sfruttare l'errore di predizione ottenuto dai compioni out-of-bag per la selezione. (analogamente per tarare il numero di replicazioni bootstrap)
library(ipred)
b <- 100
fit.bag<-bagging(y.train~., data= x.train ,nbagg= b,coob=T)
print(fit.bag)
##
## Bagging regression trees with 100 bootstrap replications
##
## Call: bagging.data.frame(formula = y.train ~ ., data = x.train, nbagg = b,
## coob = T)
##
## Out-of-bag estimate of root mean squared error: 13.84
pred.bag<-predict(fit.bag, x.test)
confronto <- rbind(confronto , c("bagging","original", err(pred.bag)))
###Foreste casuali
Usa gli alberi come modello base, e non solo vengono campionate le unità statistiche ma sono “campionate” anche le variabili. Nella stima dell'albero, ad ogni nodo viene scelto casualmente un certo numero di variabili esplicative.
Parametri di regolazione
Scelto un B abbastanza grande, sarebbe da scegliere il valore di q ottimale con i soliti metodi, ma sarebbe una cosa infinita. Il default di R per la regressione è dato da \( p/3 \).
library (randomForest)
set.seed(42)
fit.rf <- randomForest(y.train ~. , data =x.train, nodesize=1, ntree = 1000, mtry=3)
# Convergenza errore
plot(fit.rf)
pred.rf <-predict(fit.rf, x.test)
confronto <- rbind(confronto, c("rf", "original", err(pred.rf)))
library(gbm)
# Mi accontento del pacchetto gbm e lascio tutto a se stesso
fit.boost <- gbm(y.train ~.,data=x.train, n.trees=1000, cv.folds=4, distribution="gaussian")
summary(fit.boost)
## var rel.inf
## tmpd tmpd 97.85052
## dptp dptp 2.05236
## o3tmean o3tmean 0.09712
## dow dow 0.00000
## pm10tmean pm10tmean 0.00000
## so2tmean so2tmean 0.00000
## no2tmean no2tmean 0.00000
## cotmean cotmean 0.00000
pred.boost <- predict(fit.boost, newdata = x.test)
## Using 1000 trees...
confronto <- rbind(confronto , c("boosting","original", err(pred.boost)))
confronto
## model est. scale err.original err.log
## [1,] "linear" "original" "322391.856930402" "21.8303950125345"
## [2,] "best.sub" "original" "323395.347702883" "21.9033015865722"
## [3,] "step" "original" "322800.591228009" "21.859089172397"
## [4,] "step" "original" "322622.406124466" "21.8470148911844"
## [5,] "lasso.min" "original" "324563.907106534" "22.0018319024122"
## [6,] "lasso.1se" "original" "364683.180505321" "25.2828994066419"
## [7,] "lasso(tasso)" "original" "322391.856930401" "21.8303950125345"
## [8,] "mars" "original" "308452.945884283" "21.1747069169727"
## [9,] "tree(rpart)" "original" "318239.718091554" "21.8201618122354"
## [10,] "nnet1" "original" "386108.328780649" "26.899535713664"
## [11,] "bagging" "original" "311257.675668858" "21.3912721210362"
## [12,] "rf" "original" "291858.399891" "20.1315136105897"
## [13,] "boosting" "original" "339395.606214689" "23.4454159496182"