# 15. prednáška 23. listopadu 2009 # rekurzivní a parciální rezidua, # grafická regresní diagnostika # rm(list=ls()) data(primates,package="DAAG") # podá potrebná data primates attach(primates) plot(Brainwt~Bodywt,pch=16) # závislost váhy mozku na váze tela u primátu identify(Bodywt,Brainwt,rownames(primates)) plot(log(Brainwt)~log(Bodywt),pch=16) # budeme potrebovat data usporádaná podle hodnoty x (ord = order(Bodywt)) primates[ord,] brain=Brainwt[ord] body=Bodywt[ord] n = rep(NA,3) plot(log(brain)~log(body),pch=16) # jiné měřítko: plot(log(brain)~log(body),ylim=c(1,8),pch=16) abline(a <- lm(log(brain)~log(body))) # # rekurzivní rezidua # abline(a2 <- lm(log(brain)~log(body),subset=1:2),col=2) summary(a2) abline(v=log(body[3]),lty=3,col=4) print(pred <- predict(a2,newdata=list(body=body[3]),se.fit=TRUE,scale=1)) # POZOR, volbou scale=1 dosazujeme jako odhad sigma jedničku, # takže opravdu počítáme \sqrt((x_{3.})'(X_2'X_2)^{-1}x_{3.}) points(log(body[3]),pred$fit,col=2,pch=16) lines(rep(log(body[3]),2),c(log(brain[3]),pred$fit),col=4) print(n[1] <- (log(brain[3])-pred$fit)/sqrt(1+pred$se.fit^2)) # abline(a3 <- lm(log(brain)~log(body),subset=1:3),col=2) summary(a3) abline(v=log(body[4]),lty=3) print(pred <- predict(a3,newdata=list(body=body[4]),se.fit=TRUE,scale=1)) points(log(body[4]),pred$fit,col=2,pch=16) lines(rep(log(body[4]),2),c(log(brain[4]),pred$fit),col=4) print(n[2] <- (log(brain[4])-pred$fit)/sqrt(1+pred$se.fit^2)) # abline(a4 <- lm(log(brain)~log(body),subset=1:4),col=2) summary(a4) abline(v=log(body[5]),lty=3) print(pred <- predict(a4,newdata=list(body=body[5]),se.fit=TRUE,scale=1)) points(log(body[5]),pred$fit,col=2,pch=16) lines(rep(log(body[5]),2),c(log(brain[5]),pred$fit),col=4) print(n[3] <- (log(brain[5])-pred$fit)/sqrt(1+pred$se.fit^2)) abline(a5 <- lm(log(brain)~log(body),subset=1:5),col=2) summary(a5) summary(a) n # přehled rekurzivních reziduí # library(strucchange) recresid(log(brain)~log(body)) # predpokládá usporádání recresid(log(Brainwt)~log(Bodywt)) # NA USPOŘÁDÁNÍ ZÁVISÍ!!! # Harvey-Collier psi-test: t.test(n) library(lmtest) harvtest(log(brain)~log(body)) harvtest(log(Brainwt)~log(Bodywt)) # pozor na poradí!!! harvtest(log(Brainwt)~log(Bodywt),order.by=Bodywt) detach(primates) # # parciální rezidua # Policie=read.csv2("data/Policie.csv") summary(Policie) attach(Policie) # # parciální rezidua # plot(fat~height) abline(aH<-lm(fat~height)) summary(aHW<-lm(fat~height+weight)) source("parcres.R") #parcres parcres(aHW,1) summary(aH) summary(aHW) # pseudomodel parciální rezidua pro height ~ height : summary(aHWheight <- lm(resid(aHW)-coef(aHW)["height"]*height~height)) parcres(aHW,2) summary(update(aHW,.~.-height)) summary(aHW) # v knihovne car: library(car) cr.plot(aHW,height,smooth=FALSE,main="Parciální rezidua") cr.plot(aHW,height,main="Parciální rezidua") # výpocet ve std. knihovne stats print(partRes <- resid(aHW,type="partial")) dim(partRes) apply(partRes,2,sum) # mají nulové prumery colnames(partRes) str(partRes) # strucne jak objekt vypadá attr(partRes,"constant") # co to je? mean(fat) # průměr Y partRes+attr(partRes,"constant") # aby měly stejný prumer jako y: apply(partRes+attr(partRes,"constant"),2,mean) # opravdu jej mají # ještě nakreslit plot(height,partRes[,"height"]+attr(partRes,"constant")) abline(lm(partRes[,"height"]+attr(partRes,"constant")~height)) # # používané grafy reziduí # plot(fitted(aHW),resid(aHW)) abline(h=0,lty=3) # občas: plot(fitted(aHW),fat) abline(0,1,lty=3) # proti nějaké doprovodné proměnné plot(react,resid(aHW)) abline(h=0,lty=3) # konstantnost rozptylu: plot(fitted(aHW),resid(aHW)); abline(h=0,lty=3) plot(fitted(aHW),resid(aHW)^2) abline(a<-lm(resid(aHW)^2~fitted(aHW))) summary(a) # naivní pseudotest bptest(aHW,~fitted(aHW)) # knihovna lmtest bptest(aHW,~fitted(aHW),studentize=FALSE) ncv.test(aHW) # knihovna car # normalita (normovaná, obyčejná, nikoliv studentizovaná!) # NASTAV HISTORII!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! qqnorm(rstandard(aHW)) qqline(rstandard(aHW)) # jednoduchá aproximace plot(qnorm(((1:50)-3/8)/(50+1/4)),sort(rstandard(aHW))) shapiro.test(rstandard(aHW)) cor(qnorm(((1:50)-3/8)/(50+1/4)),sort(rstandard(aHW)))^2 # # autokorelace # musí mít význam pořadí pozorování, zejm. čas data(Hartnagel) # kriminalita žen (data z knihovny car) Hartnagel[sort(sample(nrow(Hartnagel),3)),] # fertilita, zaměstnanost žen, podíl vysokoškolaček, # počet tresných činů žen, .., mužů, .. summary(a <- lm(fconvict ~ tfr + partic + degrees + mconvict, data=Hartnagel)) plot(u<-resid(a)) # závislost reziduí na pořadí plot(u[-1]~u[-length(u)]) # grafy, které nabízí erko # # plot(aHW) plot(aHW,which=1:6) # standardne which = c(1:3,5) #