#1 - Chargement des donnees jura=read.table("http://informatique-mia.inra.fr/biosp/sites/ciam.inra.fr.biosp/files/jura.txt",header=T) #Analyse descriptives dim(jura) boxplot(jura) pairs(jura) #3 test de coorélation---------- cor.test(jura$Ni,jura$Co) n=nrow(jura) r=cor(jura$Ni,jura$Co) r tobs=r/sqrt(1-r^2)*sqrt(n-2) df=n-2 pvalue=2*pt(-abs(tobs),df) pvalue # 4 la regression #4.1 res1=lm(Ni~Co,jura) summary(res1) n=nrow(jura) #Question1 y=jura$Ni x=jura$Co plot(x,y) abline(coefficients(res1)) points(x,fitted(res1),pch=16,col=4) #Question2 b=cov(x,y)/var(x) a=mean(y)-b*mean(x) sigma=sqrt(1/(n-2)*sum((residuals(res1))^2)) #residual standard error sigmab=sigma^2/sum((x-mean(x))^2) #std error de Co tobs=b/sqrt(sigmab) #tvalue pour Co pvalue=2*pt(-abs(tobs),n-2) R2=cor(x,y)^2 R2 #mais on trouve aussi R2 en calculant : SCR=sum((residuals(res1))^2) SCT=sum((y-mean(y))^2) 1-SCR/SCT #ou encore 1-var(residuals(res1))/var(y) #ou encore var(fitted(res1))/var(y) # La statistique F R2/(1-R2)*(n-2) 1-pf(R2/(1-R2)*(n-2),1,n-2) #4.2 plot(x,residuals(res1)) abline(h=0) #plot(y,residuals(res1)) qqnorm(residuals(res1)) qqline(residuals(res1)) #Le graphique 2 semble indiquer un une structure ds les resioducs.. #5. La regression multiple #5.1 corp<-function(y,x,z) { #Calcule de la corrélation partielle entre # y et x en retirant l effet de z n<-cor(y,x)-cor(y,z)*cor(x,z) d<-sqrt((1-cor(y,z)^2)*(1-cor(x,z)^2)) r<-n/d print(r) } #5.2 corp(jura$Ni,jura$Cd,jura$Co) corp(jura$Ni,jura$Cr,jura$Co) corp(jura$Ni,jura$Cu,jura$Co) corp(jura$Ni,jura$Pb,jura$Co) corp(jura$Ni,jura$Zn,jura$Co) # Cr semble la meilleure variable au sens de la correlation partielle cor(jura) #meme reponse avec la correlation classique #5.3 res2=lm(Ni~Co+Cr,jura) summary(res1) summary(res2) par(mfrow=c(1,2)) plot(res1) plot(jura$Co,residuals(res1)) plot(jura$Cr,residuals(res2),col=2) # 5.4 tmp=lm(Cr~Co,jura) cor(res1$res,tmp$res) corp(jura$Ni,jura$Cr,jura$Co) # les 2 correlations sont identiques: #res1$res represente ce qui n est pas expliqué par Co dans Ni #tmp$res represente ce qui n est pas explique par Co dans Cr #on cherche la variable, qui une fois que l on retire l information apportée par Co, a des résidus # fortement correles a la partie non expliquée par CO dans Ni. #et cette information est la plus fortment correlee aux residus de res1 #5.5 correlation partielle d ordre >1 e1=residuals(res2) e2Cd=residuals(lm(Cd~Co+Cr,jura)) e2Zn=residuals(lm(Zn~Co+Cr,jura)) cor(e1,e2Cd) cor(e1,e2Zn) res3=lm(Ni~Co+Cr+Zn,jura) summary(res3) plot(y,residuals(res3)) res4=lm(Ni~Co+Cr+Zn+Cd,jura) summary(res4) res5=lm(Ni~Co+Cr+Cd+Zn+Pb,jura) summary(res5) #6 Comparaison de modele #6.1 anova(res3,res5) n=nrow(jura) q=5 #M3 p=7 #M5 sumres3=sum((residuals(res3)^2)) sumres5=sum((residuals(res5)^2)) fobs=(sumres3-sumres5)/sumres5 *((n-6)/(2)) 1-pf(fobs,2,n-7) #R2M3=var(fitted(res3))/var(y) #R2M5=var(fitted(res5))/var(y) #fobs= (n-p)/q*(R2M3-R2M5)/(1-R2M5) #6.2 AIC(res1) AIC(res2) AIC(res3) AIC(res4) AIC(res5) BIC(res1) BIC(res2) BIC(res3) BIC(res4) BIC(res5) AIC(res3) -2*logLik(res3)+2*5 BIC(res3) -2*logLik(res3)+log(n)*5 #6.3 install.packages("MuMIn") options(na.action = "na.fail") toto <- lm(Ni ~ ., data = jura) dd <- dredge(toto,rank=AIC) dd subset(dd, delta < 4) dd <- dredge(toto,rank=BIC) dd subset(dd, delta < 4)