Ce devoir maison est à envoyer par mail aude.sportisse@upmc.fr pour le lundi 28 septembre (soir) en format .Rmd et .html.

Il comporte deux exercices, on pourra utiliser les fiches 5 et 6.

Exercice 1

A partir du jeu de données icecream, nous allons étudier la consommation de glace aux Etats-Unis sur une période de 30 semaines du 18 Mars 1950 to 11 Juillet 1953. Les variables sont la consommation (Consumption en pintes par habitant), le salaire hebdomadaire (Income en dollars), le prix des glaces (Price en dollars), la température (Temp en degré fahrenheit) et la catégorie socio-professionnelle (sc).

  1. Charger le jeu de données icecream avec le nom des colonnes en utilisant la fonction read.table et en specifiant correctement les arguments sep, row.names et header. Précisez la nature des variables (qualitative ou quantitative) et faites une analyse rapide des données.
data <- read.table("icecream.txt", sep=",", header = TRUE)
summary(data)
##       cons            income          price             temp          sc    
##  Min.   :0.2560   Min.   :76.00   Min.   :0.2600   Min.   :24.00   bc  :10  
##  1st Qu.:0.3113   1st Qu.:79.25   1st Qu.:0.2672   1st Qu.:32.00   prof: 8  
##  Median :0.3515   Median :83.50   Median :0.2770   Median :43.00   wc  :12  
##  Mean   :0.3594   Mean   :84.60   Mean   :0.2749   Mean   :47.83            
##  3rd Qu.:0.3912   3rd Qu.:89.25   3rd Qu.:0.2815   3rd Qu.:63.00            
##  Max.   :0.5480   Max.   :96.00   Max.   :0.2920   Max.   :72.00
  1. Créer un jeu de données comprenant seulement les variables quantitatives. Découper aléatoirement ce jeu de données en deux échantillons: un échantillon d’apprentissage (avec 70% des données) et un échantillon de test.
data_quanti <- data[,1:4]
obs_train <- sample(1:dim(data_quanti)[1],ceiling(0.7*dim(data_quanti)[1]),replace=FALSE)
obs_test <- setdiff(1:dim(data_quanti)[1],obs_train)
data_train <- data_quanti[obs_train,]
data_test <- data_quanti[obs_test,]
  1. Sur le jeu de données d’apprentissage, représenter les nuages de points de cons en fonction de(s) variable(s) quantitative(s). Effectuer ensuite la régression linéaire de la variable cons en fonction de toute(s) le(s) variable(s) quantitative(s) du données. Ce modèle sera appelé modele1. Afficher un résumé et interpréter le résultat (donner notamment les variables significatives et en quel sens elles le sont).
modele1 <- lm(cons~., data = data_quanti)
summary(modele1)
## 
## Call:
## lm(formula = cons ~ ., data = data_quanti)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.068710 -0.014802 -0.002724  0.015850  0.070768 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.2955881  0.2560958   1.154  0.25891    
## income       0.0034569  0.0011946   2.894  0.00761 ** 
## price       -1.4175283  0.8041789  -1.763  0.08970 .  
## temp         0.0033673  0.0004531   7.431 6.85e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03768 on 26 degrees of freedom
## Multiple R-squared:  0.7059, Adjusted R-squared:  0.6719 
## F-statistic:  20.8 on 3 and 26 DF,  p-value: 4.399e-07

Testées une par une, les seules variables significatives au niveau 5% sont temp et income.

  1. Effectuer une procédure de sélection de variable par le critère d’information bayésien (BIC). Quelles sont les variables sélectionnées ?
library(leaps)
choix <- regsubsets(cons~., data = data_train, nbest = 1)
plot(choix, scale = "bic")

Le critère est optimisé pour la ligne en haut du graphique. Nous conservons donc, pour le critère BIC, le modèle à 2 variables (en plus de la constante) : temp et income

  1. Effectuer une nouvelle régression linéaire avec uniquement les variables retenues en question précédente, ce modèle sera appelé modele2. Refaire ensuite la procédure de sélection de variables.
modele2 <- lm(cons ~ temp + income, data = data_train)
summary(modele2)
## 
## Call:
## lm(formula = cons ~ temp + income, data = data_train)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.066099 -0.013556  0.004955  0.015474  0.080946 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.180840   0.136107  -1.329  0.20056    
## temp         0.003845   0.000532   7.226 1.01e-06 ***
## income       0.004166   0.001411   2.952  0.00853 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03558 on 18 degrees of freedom
## Multiple R-squared:  0.7454, Adjusted R-squared:  0.7171 
## F-statistic: 26.35 on 2 and 18 DF,  p-value: 4.494e-06
choixbis <- regsubsets(cons ~ temp + income, data = data_train, nbest = 1)
plot(choixbis, scale = "bic")

  1. Dans un vecteur Y1, stocker les valeurs de \(Y\) prédites par le modèle modele1 pour chaque individu de l’échantillon de test. Construire de même un vecteur Y2 à partir de modele2. Utiliser pour cela la fonction predict.
Y1 <- predict(modele1, newdata = data_test)
Y2 <- predict(modele2, newdata = data_test)
  1. Pour \(i=1,\dots,N_t\) où \(N_t\) est le nombre d’observation de l’échantillon de test, on note \(\hat Y_i^j\) la prévision par le modèle \(j\) du \(i\)-ème individu de l’échantillon test, et \(Y_i\) la valeur de \(Y\) observée sur le \(i\)-ème individu de l’échantillon test. Calculer \[\textrm{EQM}1=\frac{1}{N_t}\sum_{i=1}^{N_t}(\hat Y_i^1-Y_i)^2\hspace{1cm}\textrm{et}\hspace{1cm}\textrm{EQM}2=\frac{1}{N_t}\sum_{i=1}^{N_t}(\hat Y_i^2-Y_i)^2.\] Interpréter.
Y <- data_test$cons
Nt <- dim(data_test)[1]
EQM1 <- sum((Y1 - Y)^2/Nt)
EQM1
## [1] 0.001625836
EQM2 <- sum((Y2 - Y)^2/Nt)
EQM2
## [1] 0.002232533
  1. Représenter cons en fonction de(s) variable(s) qualitative(s). Effectuer une analyse de la variance. Conclure sur l’effect de(s) variable(s) qualitative(s) sur la consommation de glace.
library(ggplot2)
ggplot(data)+aes(x=sc,y=cons)+geom_boxplot() 

modele<-lm(cons~sc,data=data)
anova(modele)
## Analysis of Variance Table
## 
## Response: cons
##           Df   Sum Sq   Mean Sq F value Pr(>F)
## sc         2 0.009525 0.0047623  1.1085 0.3446
## Residuals 27 0.115999 0.0042962

La p-value est supérieur à 0.1, on en conclue que la catégorie socio-professionnelle n’a pas d’impact sur la consommation de glace.

Exercice 2

  1. Simuler 1000 échantillons \((x_i,Y_i)_{1\leq i\leq 100}\) suivant le modèle \[Y_i=1+4x_i+\varepsilon_i,\qquad \varepsilon_i\sim\mathcal N(0,1),\] avec les \(x_i \sim \mathcal{N}(3,1)\).
x <- rnorm(100,mean=3,sd=1)
n <- length(x)
B <- 100
matY <- matrix(0, ncol = n, nrow = B)
beta <- matrix(nrow = B, ncol = 2)
for (i in 1:B){
    matY[i, ] <- 1+4*x + rnorm(n)
}
  1. Pour les 1000 échantillons simulés, stocker les estimateurs des moindres carrés \(\hat \beta_0\) et \(\hat \beta_1\) renvoyés par la fonction lm dans une matrice de taille \(1000\times 2\). Stocker également les moindres carrés \(\hat \beta_0\) et \(\hat \beta_1\) renvoyés par la fonction lm dans une autre matrice de taille \(1000\times 2\) en les calculant à la main avec la formule \(\hat{\beta}=(X'X)^{-1}X'Y\). Vérifier que les valeurs sont bien les mêmes.
X <- matrix(c(rep(1,n),x),ncol=2)
beta_cal <- matrix(nrow = B, ncol = 2)
for (i in 1:B){
  beta[i, ] <- coef(lm(matY[i, ]~x))
  beta_cal[i, ] <- solve(t(X)%*%X)%*%t(X)%*%matY[i, ]
}
  1. Rappeler la moyenne et l’écart type théorique de \(\hat\beta_0\). Calculer l’écart-type. Comparer les valeurs théoriques avec la valeurs observées.
X <- cbind(rep(1, n), x)
inverse <- solve(t(X)%*%X)
sigma0 <- sqrt(inverse[1, 1])
sigma0
## [1] 0.299118

Ce qui coïncide approximativement avec la valeur observée :

sd(beta[, 1])
## [1] 0.2974473
  1. Superposer la densité théorique de \(\hat \beta_0\) et un estimateur de la densité obtenu à partir de l’échantillon obtenu précédemment.

On peut voir aussi que les distributions empirique et théorique de \(\hat\beta_0\) sont cohérentes.

plot3 <- ggplot(data.frame(beta[, 1]), aes(x = beta[, 1])) + geom_line(stat = "density")
plot3 <- plot3 + stat_function(fun = dnorm, args = list(mean = 1, sd = sigma0), col = "red") + xlab("")
plot3

  1. Reprendre la question 1 en simulant augmentant le bruit \(\epsilon\sim\mathcal{N}(0,10)\). Comme dans la question 2, pour les 1000 échantillons simulés, stocker les estimateurs des moindres carrés \(\hat \beta_0\) et \(\hat \beta_1\) renvoyés par la fonction lm dans une matrice de taille \(1000\times 2\). Que remarque-t-on sur les valeurs estimées \(\hat{\beta}\) ? Expliquer pourquoi (avec un graphique idéalement).
x <- rnorm(100,mean=3,sd=1)
n <- length(x)
B <- 100
matY <- matrix(0, ncol = n, nrow = B)
beta <- matrix(nrow = B, ncol = 2)
for (i in 1:B){
    matY[i, ] <- 1+4*x + rnorm(n,sd=10)
    beta[i, ] <- coef(lm(matY[i, ]~x))
}

Le niveau de bruit est très élevé.

plot(1+4*x+ rnorm(n,sd=10))
points(1:100,1+4*x,col="red")