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.
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).
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
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,]
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.
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
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")
Y1 <- predict(modele1, newdata = data_test)
Y2 <- predict(modele2, newdata = data_test)
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
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.
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)
}
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, ]
}
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
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
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")