demo3=read.table("C:/Documents and Settings/CHEIK/Bureau/ETE 2014/Bases de données individuelles/Mortality.csv",sep=";",dec=",",header=TRUE) head(demo3) ##################### # Question(2) # ##################### ## Analyse statistique du taux de mortalité avec chacune des variables explicatives**** # moyenne, variance, etc summary(demo3$txm99) summary(demo3$txpop65) summary(demo3$nbm000) summary(demo3$espvi) summary(demo3$revbh) # Boite à moustache par(mfrow=c(3,2)) boxplot(demo3$txm99) boxplot(demo3$txpop65) boxplot(demo3$nbm000) boxplot(demo3$espvi) boxplot(demo3$revbh) # répérage des observations aberrantes resul1=boxplot(demo3$txm99) vatxm99=resul1$out which(demo3$txm99%in%vatxm99) resul2=boxplot(demo3$revbh) varevbh=resul2$out which(demo3$revbh%in%varevbh) # Les nuages de points plot(demo3[,2:6]) par(mfrow=c(2,2)) plot(demo3$txpop65,demo3$txm99,xlab="%population âgé de 65 ans et plus",ylab="Taux de mortalité") plot(demo3$nbm000,demo3$txm99,xlab="Le nombre de medecins pour 1000 hbts",ylab="Taux de mortalité") plot(demo3$espvi,demo3$txm99,xlab="Espérance de vie",ylab="Taux de mortalité") plot(demo3$revbh,demo3$txm99,xlab="Le revenu brut par habitant",ylab="Taux de mortalité") # Test corrélation linéaire cor.test(demo3$txm99,demo3$txpop65,data=demo3) cor.test(demo3$txm99,demo3$nbm000,data=demo3) cor.test(demo3$txm99,demo3$espvi,data=demo3) cor.test(demo3$txm99,demo3$revbh,data=demo3) ##################################### # Question(3) # ##################################### ##Regression linéaire multiple ## test de boxcox ## ltxm99=log(demo3$txm99) demo3$ltxm99=ltxm99 reg0=lm(demo3$ltxm99~demo3$txpop65+demo3$nbm000+demo3$espvi+demo3$revbh,data=demo3) library(MASS) boxcox(reg0) ## Modèle complet reg1=lm(demo3$txm99~demo3$txpop65+demo3$nbm000+demo3$espvi+demo3$revbh,data=demo3) summary(reg1) # H0: scr1=0 (ou R2=0) # H1: scr1!=0 (ou R2!=0) # la stat. associée est celle de Fisher ################################ # Question (4) # ################################ ## Analyse des résidus par(mfrow=c(4,3)) plot(reg1,ask=FALSE) qqnorm(residuals(reg1)) qqline(residuals(reg1),col="red") plot(residuals(reg1),ylab="Residus",main="Représentation des résidus",col="blue") plot(density(residuals(reg1)),main="Distribution des résidus",ylab="Résidus") ## Test de normalité des résidus ## necessite des packages tseries ## Jarque Bera library(tseries) jarque.bera.test(reg1$residuals) qchisq(0.95,2) ## Shapiro shapiro.test(reg1$residuals) ######################### # Question(5) # ######################## plot(reg1,which=1) # les observations correspondants aux résidus >2: obs 87, obs 35 et obs 5 ######################### # Question(6) # ######################### demo31=demo3[reg1$residuals<2,] reg2=lm(demo31$txm99~demo31$txpop65+demo31$nbm000+demo31$espvi+demo31$revbh,data=demo31) ################################### # Question(7) # ################################### ## test d 'homoscedascité ## Breush et Pagan ## avec base de données demo3 ##A la main reg1=lm(demo3$txm99~demo3$txpop65+demo3$nbm000+demo3$espvi+demo3$revbh,data=demo3) regu2=lm(I(residuals(reg1)^2)~demo3$txpop65+demo3$nbm000+demo3$espvi+demo3$revbh,data=demo3) # R-squared r2aux=summary(regu2)$r.square # stat. de Breush et Pagan bp=99*r2aux pvalue=1-pchisq(bp,4) pvalue # sinon en cosidérant les carrés et produits des variables txpop652=demo3$txpop65^2 nbm0002=demo3$nbm00^2 espvi2=demo3$espvi^2 revbh2=demo3$revbh^2 prod1=demo3$txpop65*demo3$nbm00 prod2=demo3$txpop65*demo3$espvi prod3=demo3$txpop65*demo3$revbh prod4=demo3$nbm00*demo3$espvi prod5=demo3$nbm00*demo3$revbh prod6=demo3$espvi*demo3$revbh prod7=demo3$txpop65*demo3$nbm00*demo3$espvi prod8=demo3$txpop65*demo3$nbm00*demo3$revbh prod9=demo3$nbm00*demo3$espvi*demo3$revbh prod10=demo3$txpop65*demo3$nbm00*demo3$espvi*demo3$revbh demo3$txpop652=txpop652 demo3$nbm0002=nbm0002 demo3$espvi2=espvi2 demo3$revbh2=revbh2 demo3$prod1=prod1 demo3$prod2=prod2 demo3$prod3=prod3 demo3$prod4=prod4 demo3$prod5=prod5 demo3$prod6=prod6 demo3$prod7=prod7 demo3$prod8=prod8 demo3$prod9=prod9 demo3$prod10=prod10 head(demo3) # modèle du carré des résidus regu2=lm(I(residuals(reg1)^2)~demo3$txpop65+demo3$nbm000+demo3$espvi+demo3$revbh+txpop652+nbm0002+espvi2+revbh2+prod1+prod2+prod3+prod4+prod5+prod6+prod7+prod8+prod9+prod10,data=demo3) # R-squared r2aux=summary(regu2)$r.square # stat. de Breush et Pagan bp=99*r2aux pvalue=1-pchisq(bp,18) pvalue # automatique, nécessite les packeges AER library(AER) bptest(reg1,data=demo3) ####################### # Question(8) # ###################### sce1=deviance(reg1) sigmar=sqrt(sce1/94) txm99sigmar=(demo3$txm99/sigmar) txpop65sigmar=(demo3$txpop65/sigmar) nbm000sigmar=(demo3$nbm000/sigmar) espvisigmar=(demo3$espvi/sigmar) revbhsigmar=(demo3$revbh/sigmar) reg1c=lm(txm99sigmar~txpop65sigmar+nbm000sigmar+espvisigmar+revbhsigmar,data=demo3) summary(reg1c) library(tseries) jarque.bera.test(reg1c$residuals) ######################### # Question(9) # ######################## demo32=demo31[demo31$txm99<22.70,] newdemo3=demo32[demo32$revbh<60870,] reg2c=lm(newdemo3$txm99~newdemo3$txpop65+newdemo3$nbm000+newdemo3$espvi+newdemo3$revbh,data=newdemo3) jarque.bera.test(reg2c$residuals) plot(density(reg2c$residuals),main="Distribution des résidus",ylab="Résidus") curve(dnorm(x,mean=mean(residuals(reg2c)),sd=sd(residuals(reg2c))),add=TRUE,lty=2,col="red") ########################## # Question(10) # ########################## par(mfrow=c(2,2)) plot(newdemo3$txpop65,reg2c$residuals,xlab="txpop65",ylab="Résidus",main="Residu vs txpop65"); abline(h=0,col="red") plot(newdemo3$nbm000,reg2c$residuals,xlab="nbm000",ylab="Résidus",main="Residu vs nbm000"); abline(h=0,col="red") plot(newdemo3$espvi,reg2c$residuals,xlab="espvi",ylab="Résidus",main="Residu vs espvi"); abline(h=0,col="red") plot(newdemo3$revbh,reg2c$residuals,xlab="revbh",ylab="Résidus",main="Residu vs revbh"); abline(h=0,col="red") ########################## # Question (11) # ########################## Y=newdemo3$txm99 X1=newdemo3$txpop65 X2=newdemo3$nbm000 X3=newdemo3$espvi X4=newdemo3$revbh AIC(lm(Y~I(X1^2)+X2+X3+X4,data=newdemo3)) AIC(lm(Y~I(X2^2)+X1+X3+X4,data=newdemo3)) AIC(lm(Y~I(X3^2)+X1+X2+X4,data=newdemo3)) AIC(lm(Y~I(X4^2)+X1+X2+X3,data=newdemo3)) AIC(lm(Y~I(X1^2)+I(X2^2)+X3+X4,data=newdemo3)) AIC(lm(Y~I(X1^2)+I(X3^2)+X2+X4,data=newdemo3)) AIC(lm(Y~I(X1^2)+I(X4^2)+X2+X3,data=newdemo3)) AIC(lm(Y~I(X2^2)+I(X3^2)+X1+X4,data=newdemo3)) AIC(lm(Y~I(X2^2)+I(X4^2)+X1+X3,data=newdemo3)) AIC(lm(Y~I(X3^2)+I(X4^2)+X1+X2,data=newdemo3)) AIC(lm(Y~I(X1^2)+I(X3^2)+I(X2^2)+X4,data=newdemo3)) AIC(lm(Y~I(X1^2)+I(X2^2)+I(X4^2)+X3,data=newdemo3)) AIC(lm(Y~I(X1^2)+I(X4^2)+I(X3^2)+X2,data=newdemo3)) AIC(lm(Y~I(X2^2)+I(X4^2)+I(X3^2)+X1,data=newdemo3)) AIC(lm(Y~I((X1+X2+X3+X4)^2),data=newdemo3)) AIC(lm(Y~I((1/(X1+X2+X3+X4)^2)),data=newdemo3)) AIC(lm(Y~I(1/X1)+I(1/X2)+X3+X4,data=newdemo3)) AIC(lm(Y~I(1/X1)+I(1/X3)+X2+X4,data=newdemo3)) AIC(lm(Y~I(1/X1)+I(1/X4)+X2+X3,data=newdemo3)) AIC(lm(Y~I(1/X3)+I(1/X2)+X1+X4,data=newdemo3)) AIC(lm(Y~I(1/X2)+I(1/X4)+X1+X3,data=newdemo3)) AIC(lm(Y~I(1/X3)+I(1/X4)+X1+X2,data=newdemo3)) AIC(lm(Y~I(1/X1)+I(1/X2)+I(1/X3)+X4,data=newdemo3)) AIC(lm(Y~I(1/X1)+I(1/X2)+I(1/X4)+X3,data=newdemo3)) AIC(lm(Y~I(1/X3)+I(1/X2)+I(1/X4)+X1,data=newdemo3)) AIC(lm(Y~I(1/X3)+I(1/X2)+I(1/X4)+I(1/X1),data=newdemo3)) AIC(lm(Y~I((1/(X1+X2+X3+X4))),data=newdemo3)) AIC(lm(Y~I((1/(X1+X2)))+X3+X4,data=newdemo3)) AIC(lm(Y~I((1/(X2+X3+X4)))+X1,data=newdemo3)) AIC(lm(Y~I((1/(X1+X2)))+I((1/(X3+X4))),data=newdemo3)) AIC(lm(Y~I((1/(X1+X3)))+I((1/(X2+X4))),data=newdemo3)) AIC(lm(Y~I((1/(X1+X4)))+I((1/(X2+X3))),data=newdemo3)) AIC(lm(Y~I((1/(X1+X3)))+X2+X4,data=newdemo3)) AIC(lm(Y~I((1/(X1+X4)))+X2+X3,data=newdemo3)) AIC(lm(Y~I((1/(X2+X3)))+X1+X4,data=newdemo3)) AIC(lm(Y~I((1/(X2+X4)))+X1+X3,data=newdemo3)) AIC(lm(Y~I((1/(X1+X2+X3)))+I(1/X4),data=newdemo3)) AIC(lm(Y~I((1/(X1+X2+X4)))+I(1/X3),data=newdemo3)) AIC(lm(Y~I((1/(X3+X2+X4)))+I(1/X1),data=newdemo3)) AIC(lm(Y~I((1/(X3+X4)))+X1+X2,data=newdemo3)) AIC(lm(Y~I((1/(X1+X2+X3)))+X4,data=newdemo3)) AIC(lm(Y~I((1/(X1+X2+X4)))+X3,data=newdemo3)) AIC(lm(Y~I((1/(X1+X3+X4)))+X2,data=newdemo3)) AIC(lm(Y~I(1/X1)+X2+X3+X4,data=newdemo3)) AIC(lm(Y~I(1/X2)+X1+X3+X4,data=newdemo3)) AIC(lm(Y~I(1/X3)+X1+X2+X4,data=newdemo3)) AIC(lm(Y~I(1/X4)+X1+X2+X3,data=newdemo3)) AIC(lm(Y~X1+X2+log(X3)+X4,data=newdemo3)) ################################# # Question (12) # ################################# reg3=lm(Y~X1+X2+log(X3)+X4,data=newdemo3) aovreg3=aov(Y~X1+X2+log(X3)+X4,data=newdemo3) ## analyse de la variance aovreg3 confint(reg3,level=0.95) ## intervalles de confiance des coefficients estimés ## Méthode de stepwise ## Modèles réduits reg4=lm(Y~X1+X2+log(X3),data=newdemo3) summary(reg4) ## analyse de la variance du modèle # H0: scr4=0 # H1: scr4!=0 # la stat. associée est celle de Fisher anova(reg4) aov(reg4) reg5=lm(Y~X1+log(X3),data=demo3) summary(reg5) ## analyse de la variance du modèle # H0: scr5=0 # H1: scr5!=0 # la stat. associée est celle de Fisher anova(reg5) aov(reg5) ## Test de Fisher de reg3 vs reg1 # H0: beta2=beta4=0 (reg5) # H1: beta2!=0 ou beta4!=0 (reg3) anova(reg3,reg5) n=length(newdemo3$txm99) sce1=deviance(reg1) # somme des carrés des résidus du modèle 1 sce3=deviance(reg3) # somme des carrés des résidus du modèle 3 ####### à la main #### # sce3 X1=cbind(1,newdemo3$txpop65,newdemo3$nbm000,log(newdemo3$espvi),newdemo3$revbh) B1=solve(t(X1)%*%X1)%*%t(X1)%*%newdemo3$txm99 epsilon1=newdemo3$txm99-X1%*%B1 sce1=sum(epsilon1^2) sce1 # sce5 X3=cbind(1,newdemo3$txpop65,log(newdemo3$espvi)) B3=solve(t(X3)%*%X3)%*%t(X3)%*%newdemo3$txm99 epsilon3=newdemo3$txm99-X3%*%B3 sce3=sum(epsilon3^2) sce3 # scr3 : somme des carrés de la régression 3 (modèle 3) scr1=sum((X1%*%B1-mean(demo3$txm99))^2) # scr5 : somme des carrés de la régression 5 (modèle 5) scr3=sum((X3%*%B3-mean(demo3$txm99))^2) dl1=n-5 dl3=n-3 Fcal=((sce3-sce1)/(dl3-dl1))/(sce1/dl1) Fcal qf(0.95,dl3-dl1,dl1) ## Test de normalité des résidus du modèle 3 jarque.bera.test(reg3$residuals) shapiro.test(residuals(reg3)) plot(density(reg3$residuals),main="Distribution des résidus",ylab="Résidus") curve(dnorm(x,mean=mean(residuals(reg3)),sd=sd(residuals(reg3))),add=TRUE,lty=2,col="red") ## Test de normalité des résidus du modèle 4 jarque.bera.test(reg4$residuals) shapiro.test(residuals(reg4)) plot(density(reg4$residuals),main="Distribution des résidus",ylab="Résidus") curve(dnorm(x,mean=mean(residuals(reg4)),sd=sd(residuals(reg4))),add=TRUE,lty=2,col="red") ## Test de normalité des résidus du modèle 5 jarque.bera.test(reg5$residuals) shapiro.test(residuals(reg5)) plot(density(reg5$residuals),main="Distribution des résidus",ylab="Résidus") curve(dnorm(x,mean=mean(residuals(reg5)),sd=sd(residuals(reg5))),add=TRUE,lty=2,col="red") par(mfrow=c(2,2)) plot(newdemo3$txpop65,reg3$residuals,xlab="txpop65",ylab="Résidus",main="Residu vs txpop65"); abline(h=0,col="red") plot(newdemo3$nbm000,reg3$residuals,xlab="nbm000",ylab="Résidus",main="Residu vs nbm000"); abline(h=0,col="red") plot(log(newdemo3$espvi),reg3$residuals,xlab="log(espvi)",ylab="Résidus",main="Residu vs espvi"); abline(h=0,col="red") plot(newdemo3$revbh,reg3$residuals,xlab="revbh",ylab="Résidus",main="Residu vs revbh"); abline(h=0,col="red") ## Test de Fisher de reg5 vs reg4 # H0: beta2=beta4=0 (reg5) # H1: beta2!=0 ou beta4!=0 (reg4) anova(reg5,reg4) # camparaison de variance dans les deux modèles avec et sans log ## critère BIC BIC(reg3,reg4,reg5) ## R2 ajustés summary(reg3) # R2 ajusté = 0.9778 summary(reg4) # R2 ajusté = 0.9779 summary(reg5) # R2 ajusté = 0.9753 ############################## # Question (13) # ############################# ## Test de Chow: homogéneité de l'échantillon # Effet continent plot(demo3$txpop65[demo3$dumaf==1],demo3$txm99[demo3$dumaf==1],xlim=c(range(demo3$txpop65)[1],range(demo3$txpop65)[2]),ylim=c(range(demo3$txm99)[1],range(demo3$txm99)[2]),col="Blue", main="visualisation de l'effet continent",xlab="txpop65",ylab="txm99") points(demo3$txpop65[demo3$dumaf==0],demo3$txm99[demo3$dumaf==0],col="Red") legend(15,22,legend=c("Pays d'Afrique","Autres pays"),bty="o",lwd=2,col=c("blue","red")) plot(demo3$nbm000[demo3$dumaf==1],demo3$txm99[demo3$dumaf==1],xlim=c(range(demo3$nbm000)[1],range(demo3$nbm000)[2]),ylim=c(range(demo3$txm99)[1],range(demo3$txm99)[2]),col="Blue", main="visualisation de l'effet continent",xlab="nbm00",ylab="txm99") points(demo3$nbm000[demo3$dumaf==0],demo3$txm99[demo3$dumaf==0],col="Red") plot(demo3$espvi[demo3$dumaf==1],demo3$txm99[demo3$dumaf==1],xlim=c(range(demo3$espvi)[1],range(demo3$espvi)[2]),ylim=c(range(demo3$txm99)[1],range(demo3$txm99)[2]),col="Blue", main="visualisation de l'effet continent",xlab="espvi",ylab="txm99") points(demo3$espvi[demo3$dumaf==0],demo3$txm99[demo3$dumaf==0],col="Red") plot(demo3$revbh[demo3$dumaf==1],demo3$txm99[demo3$dumaf==1],xlim=c(range(demo3$revbh)[1],range(demo3$revbh)[2]),ylim=c(range(demo3$txm99)[1],range(demo3$txm99)[2]),col="Blue", main="visualisation de l'effet continent",xlab="revbh",ylab="txm99") points(demo3$revbh[demo3$dumaf==0],demo3$txm99[demo3$dumaf==0],col="Red") ## 1) subdivision de l'échantillon en deux sous échantillons newdemo31=newdemo3[newdemo3$dumaf==1,] ## pays Afrique uniquement newdemo32=newdemo3[newdemo3$dumaf==0,] ## autres pays # regression avec le modèle 3, échatillon pays d'Afrique reg31=lm(newdemo31$txm99~newdemo31$txpop65+newdemo31$nbm000+log(newdemo31$espvi)+newdemo31$revbh,data=newdemo31) summary(reg31) AIC(reg31) BIC(reg31) jarque.bera.test(reg31$residuals) shapiro.test(residuals(reg31)) # regression avec le modèle 3, échantillon autres pays reg32=lm(newdemo32$txm99~newdemo32$txpop65+newdemo32$nbm000+log(newdemo32$espvi)+newdemo32$revbh,data=newdemo32) summary(reg32) AIC(reg32) BIC(reg32) jarque.bera.test(reg32$residuals) shapiro.test(residuals(reg32)) # 2) calcul de la statistique de Fisher associée SCE=deviance(reg3) # somme des carrés des résidus du modèle 3, l'ensemble de l'échatillon SCE1=deviance(reg31) # somme des carrés des résidus du modèle 3, échantillon Afrique SCE2=deviance(reg32) # somme des carrés des résidus du modèle 3, échantillon autres pays SCE12=SCE1+SCE2 Fcal=(((SCE-SCE12)/(89-84))/(SCE12/84)) Fcal qf(0.95,5,84) # sous excel INVERSE.LOI.F(0,05;5;84) # on compare Fcal à qf # niveau de significativité associé # sous excel =loi.F(Fcal;5;84) ############################# # Question (14) # ############################# ## Prédictions avec intervalles de confiance # Avec modèle 3 newdemo33=data.frame(newdemo3[1:84,]) # ensemble de pays prevtxm_3=predict(reg3, newdemo33, level = 0.9, interval= "confidence") prevtxm_3 newdemo311=data.frame(newdemo31[1:20,]) ## avec le modèle pays Afrique prev1txm=predict(reg31, newdemo311, level = 0.9, interval= "confidence") newdemo312=data.frame(newdemo32[1:50,]) ## avec le modèle autres pays prev2txm=predict(reg32, newdemo312, level = 0.9, interval= "confidence") ### Autres I=c(1,trunc(nrow(newdemo3)/2),nrow(newdemo3)) baseestim=newdemo3[-I,] basetest=newdemo3[I,] ### estimations avec la base diminuée de 3 obs. reg3_91=lm(baseestim$txm99~baseestim$txpop65+baseestim$nbm000+log(baseestim$espvi)+baseestim$revbh,data=baseestim) reg4_91=lm(baseestim$txm99~baseestim$txpop65+baseestim$nbm000+log(baseestim$espvi),data=baseestim) reg5_91=lm(baseestim$txm99~baseestim$nbm000+log(baseestim$espvi),data=baseestim) # prediction avec le modèle 3 prevtxm_3=predict(reg3_91,newdata=basetest,interval="pred",level=0.95) # prediction avec le modèle 4 prevtxm_4=predict(reg4_91, newdata=basetest,interval="pred",level=0.95) # prediction avec le modèle 5 prevtxm_5=predict(reg5_91, newdata=basetest,interval="pred",level=0.95)