En vue de se préparer à l’examen (qui aura lieu vendredi prochain), je mets en ligne quelques exercices de séries temporelles, ainsi que quelques questions à choix multiples. On regarde tout ça mercredi.
Daily Archives: 18/03/2016
Classification on the German Credit Database
In our data science course, this morning, we’ve use random forrest to improve prediction on the German Credit Dataset. The dataset is
> url="http://freakonometrics.free.fr/german_credit.csv" > credit=read.csv(url, header = TRUE, sep = ",")
Almost all variables are treated a numeric, but actually, most of them are factors,
> str(credit) 'data.frame': 1000 obs. of 21 variables: $ Creditability : int 1 1 1 1 1 1 1 1 1 1 ... $ Account.Balance : int 1 1 2 1 1 1 1 1 4 2 ... $ Duration : int 18 9 12 12 12 10 8 ... $ Purpose : int 2 0 9 0 0 0 0 0 3 3 ...
(etc). Let us convert categorical variables as factors,
> F=c(1,2,4,5,7,8,9,10,11,12,13,15,16,17,18,19,20) > for(i in F) credit[,i]=as.factor(credit[,i])
Let us now create our training/calibration and validation/testing datasets, with proportion 1/3-2/3
> i_test=sample(1:nrow(credit),size=333) > i_calibration=(1:nrow(credit))[-i_test]
The first model we can fit is a logistic regression, on selected covariates
> LogisticModel <- glm(Creditability ~ Account.Balance + Payment.Status.of.Previous.Credit + Purpose + Length.of.current.employment + Sex...Marital.Status, family=binomial, data = credit[i_calibration,])
Based on that model, it is possible to draw the ROC curve, and to compute the AUC (on ne validation dataset)
> fitLog <- predict(LogisticModel,type="response", + newdata=credit[i_test,]) > library(ROCR) > pred = prediction( fitLog, credit$Creditability[i_test]) > perf <- performance(pred, "tpr", "fpr") > plot(perf) > AUCLog1=performance(pred, measure = "auc")@y.values[[1]] > cat("AUC: ",AUCLog1,"\n") AUC: 0.7340997
An alternative is to consider a logistic regression on all explanatory variables
> LogisticModel <- glm(Creditability ~ ., + family=binomial, + data = credit[i_calibration,])
We might overfit, here, and we should observe that on the ROC curve
> fitLog <- predict(LogisticModel,type="response", + newdata=credit[i_test,]) > pred = prediction( fitLog, credit$Creditability[i_test]) > perf <- performance(pred, "tpr", "fpr") > plot(perf) > AUCLog2=performance(pred, measure = "auc")@y.values[[1]] > cat("AUC: ",AUCLog2,"\n") AUC: 0.7609792
There is a slight improvement here, compared with the previous model, where only five explanatory variables were considered.
Consider now some regression tree (on all covariates)
> library(rpart) > ArbreModel <- rpart(Creditability ~ ., + data = credit[i_calibration,])
We can visualize the tree using
> library(rpart.plot) > prp(ArbreModel,type=2,extra=1)
The ROC curve for that model is
> fitArbre <- predict(ArbreModel, + newdata=credit[i_test,], + type="prob")[,2] > pred = prediction( fitArbre, credit$Creditability[i_test]) > perf <- performance(pred, "tpr", "fpr") > plot(perf) > AUCArbre=performance(pred, measure = "auc")@y.values[[1]] > cat("AUC: ",AUCArbre,"\n") AUC: 0.7100323
As expected, a single has a lower performance, compared with a logistic regression. And a natural idea is to grow several trees using some boostrap procedure, and then to agregate those predictions.
> library(randomForest) > RF <- randomForest(Creditability ~ ., + data = credit[i_calibration,]) > fitForet <- predict(RF, + newdata=credit[i_test,], + type="prob")[,2] > pred = prediction( fitForet, credit$Creditability[i_test]) > perf <- performance(pred, "tpr", "fpr") > plot(perf) > AUCRF=performance(pred, measure = "auc")@y.values[[1]] > cat("AUC: ",AUCRF,"\n") AUC: 0.7682367
Here this model is (slightly) better than the logistic regression. Actually, if we create many training/validation samples, and compare the AUC, we can observe that – on average – random forests perform better than logistic regressions,
> AUC=function(i){ + set.seed(i) + i_test=sample(1:nrow(credit),size=333) + i_calibration=(1:nrow(credit))[-i_test] + LogisticModel <- glm(Creditability ~ ., + family=binomial, + data = credit[i_calibration,]) + summary(LogisticModel) + fitLog <- predict(LogisticModel,type="response", + newdata=credit[i_test,]) + library(ROCR) + pred = prediction( fitLog, credit$Creditability[i_test]) + AUCLog2=performance(pred, measure = "auc")@y.values[[1]] + RF <- randomForest(Creditability ~ ., + data = credit[i_calibration,]) + fitForet <- predict(RF, + newdata=credit[i_test,], + type="prob")[,2] + pred = prediction( fitForet, credit$Creditability[i_test]) + AUCRF=performance(pred, measure = "auc")@y.values[[1]] + return(c(AUCLog2,AUCRF)) + } > A=Vectorize(AUC)(1:200) > plot(t(A))
Prévision de l’intérêt pour le chocolat, avec des modèles ARIMA
Hier en cours, on s’était lancé sur la prévision de la série obtenue sur Google Trends pour le mot clé “chocolat“. La série est la suivante
> url = "http://freakonometrics.free.fr/report-chocolat.csv" > report = read.csv(url,skip=6, + header=FALSE,nrows=636) > X=report$V2 > plot(X,type="l")
Comme on l’avait évoqué, c’est plus simple de construire un modèle ARIMA sur des séries mensuelles qu’hebdomadaire. Donc on va mensualiser la série. Ah oui, on va aussi enlever les deux dernières années pour faire une comparaison de modèles (en confrontant les prévision avec les observations)
> date1=substr(as.character(report$V1),1,10) > date2=substr(as.character(report$V1),14,23) > d1=as.Date(date1,"%Y-%m-%d") > d2=as.Date(date2,"%Y-%m-%d") > Vect_J=Vect_M=Vect_X=Vect_A= NULL > for(i in 1:length(d1)){ + jours = seq(d1[i],d2[i],length=7) + mois = as.POSIXlt(jours)$mon+1 + annees = as.POSIXlt(jours)$year+1900 + donnees = rep(X[i],7) + Vect_J = c(Vect_J,jours) + Vect_M = c(Vect_M,mois) + Vect_A = c(Vect_A,annees) + Vect_X = c(Vect_X,donnees) + } > base = data.frame(J = Vect_J,M = Vect_M, + A = Vect_A,X = Vect_X) > base$AM = base$A*100+base$M > Y = tapply(base$X,as.factor(base$AM),mean) > Z = ts(as.numeric(Y[1:(146-24)]), + start=c(2004,1),frequency=12) > plot(Z)
Compte tenu de l’amplitude croissante des cycles, il peut être intéressant, ici, de passer au logarithme de la série. Et comme on observe une rupture dans la tendance (linéaire au début, et relativement constante ensuite), on peut ajuste une régression linéaire avec une rupture (mais continue)
> plot(as.numeric(Z),type="l") > X=log(as.numeric(Z)) > T=1:length(X) > db=data.frame(X,T) > trend=lm(X~T+I((T-80)*(T>80)),data=db) > plot(X,type="l") > lines(predict(trend),col="red")
C’est le résidu de cette série qu’on va chercher à modéliser,
> Y=residuals(trend) > plot(Y,type="l")
Pour modéliser cette série, on peut commencer par regarder sa série d’autocorrélation
> plot(acf(Y,lag=36),lwd=5)
et d’autocorrélation partielles,
> plot(pacf(Y,lag=36),lwd=5)
La série est stationnaire, mais il y a une composante cyclique forte (peut être une racine unité saisonnière ?). On peut tenter un premier modèle AR, ou ARMA (avec un AR les résidus ne sont pas un bruit blanc).
> model1 = arima(Y,order=c(12,0,12), + seasonal = list(order = c(0, 0, 0), period = 12))
La série des résidus est ici un bruit blanc
> E=residuals(model1) > plot(acf(E,lag=36),lwd=5)
On peut alors utiliser ce modèle pour faire des prévisions, sur la série initiale (avant de prendre le logarithme, et d’ajuster une tendance)
> futur=length(X)+1:24 > Yp=predict(model1,n.ahead=24)$pred+ + predict(trend,newdata=data.frame(T=futur)) > Xp=exp(as.numeric( + Yp+.5*predict(model1,n.ahead=24)$se^2)) > plot(as.numeric(Z),type="l",xlim=c(1,150)) > lines(futur,Xp,col="red") > Y2=tapply(base$X,as.factor(base$AM),mean) > obs_reel=as.numeric(Y2[(146-24)+1:24]) > lines(futur,obs_reel,col="blue")
avec en rouge notre modèle, et en bleu les vraies observations. On peut alors calculer la somme des carrés des erreurs, sur ces 24 observations,
> sum( (obs_reel-Xp)^2 ) [1] 190.9722
Mais on peut tenter d’autres modèles, par exemple en changeant la tendance (on ne va pas le faire) ou en changeant le modèle ARIMA (en intégrant une racine unité saisonnière)
> model3=arima(Y,order=c(0,0,2), + seasonal = list(order = c(0, 1, 1), period = 12)) > E=residuals(model3) > model3 Coefficients: ma1 ma2 sma1 0.2246 0.3034 -0.9999 s.e. 0.0902 0.0925 0.3503 sigma^2 estimated as 0.002842: log likelihood = 152.37, aic = -296.75
On vérifie que les résidus sont effectivement un bruit blanc
> plot(acf(E,lag=36),lwd=5) > Box.test(E,lag=12) Box-Pierce test data: E X-squared = 6.326, df = 12, p-value = 0.8988
On peut alors faire une prévision sur notre série brute,
> futur=length(X)+1:24 > Yp=predict(model3,n.ahead=24)$pred+ + predict(trend,newdata=data.frame(T=futur)) > Xp=exp(as.numeric( + Yp+.5*predict(model3,n.ahead=24)$se^2)) > plot(as.numeric(Z),type="l",xlim=c(1,150)) > lines(futur,as.numeric(Xp),col="red") > Y2=tapply(base$X,as.factor(base$AM),mean) > obs_reel=as.numeric(Y2[(146-24)+1:24]) > lines(futur,obs_reel,col="blue")
La somme des carrés des erreurs est ici un peu plus faible
> sum( (obs_reel-Xp)^2 ) [1] 173.8138
C’est à dire que ce second modèle est meilleur que le précédant, sur les deux dernières années. C’est probablement celui que l’on devrait utiliser pour faire des prévisions pour les années à venir (à moins qu’en continuant l’analyse on n’en trouve un meilleur).