library(car)

# Nastaveni generatoru nahodnych cisel
set.seed(100)

#	Nastaveni parametru a generovani nahodnych vyberu
k<-15;n<-2*k-1;sigma1<-3.5
xx<-seq(1,k,by=0.5)+sigma1*rnorm(n)
a<-5;b<-3;sigma2<-3.5;
yy<-a+b*xx+sigma2*rnorm(n)
x1<-c(xx,39)
y1<-c(yy,125)
x2<-x1
y2<-c(yy,11)


# Grafy body a regresnich primek
xlim<-range(c(x1,x2))
ylim<-range(c(y1,y2))
par(mfrow=c(1,2),mar=c(2,2,2,0.5)+0.05)
plot(x1[1:n],y1[1:n],pch=3,xlim=xlim,ylim=ylim)
points(x1[n+1],y1[n+1],pch=21,col="red",bg="orange")
abline(lm(y1~x1),col="dodgerblue",lwd=2)
abline(lm(y1[-(n+1)]~x1[-(n+1)]),col="violet",lwd=2,lty=2)
text(x1[n+1],y1[n+1],"leverage point and outlier   ",adj=c(1,0.5))
mtext("Data set 1")
plot(x2[1:n],y2[1:n],pch=3,xlim=xlim,ylim=ylim)
points(x2[n+1],y2[n+1],pch=21,col="red",bg="orange")
abline(lm(y2~x2),col="dodgerblue",lwd=2)
abline(lm(y2[-(n+1)]~x2[-(n+1)]),col="violet",lwd=2,lty=2)
text(x2[n+1],y2[n+1],"leverage point   ",adj=c(1,0.5))
mtext("Data set 2")


# Vypocet obou linearnich regresnich modelu
m1<-lm(y1~x1)
summary(m1)

m2<-lm(y2~x2)
summary(m2)


# Diagnosticke grafy pro model 1
par(mfrow=c(3,2),mar=c(5,5,1.5,0)+0.05)
plot(m1,which=1:6)


# Diagnosticke grafy pro model 1
par(mfrow=c(3,2),mar=c(5,5,1.5,0)+0.05)
plot(m2,which=1:6)


# Williamsuv graf, obsah bublinek je umerny Cookove vzdalenosti
k<-length(coef(m1))
n<-length(x1)
sr<-rstudent(m1)
Lr2<-abs(sr)>2
hii<-hatvalues(m1)
Lh2<-hii>2*k/n
L<-Lr2 | Lh2
par(mfrow=c(1,1),mar=c(5,5,1.5,0)+0.05)
influencePlot(m1)
text(hii[L],sr[L],as.character((1:n)[L]),adj=c(0.5,0.5),col="blue",cex=0.75)


# Cookova vzdalenost
Di<-cooks.distance(m1)
# Vypiseme podezrela pozorovani pro Dataset 1
cat(paste("meze pro hii:  2*k/n=",round(2*k/n,4),
          "   3*k/n=",round(3*k/n,4),"\n",sep=""))
cbind((1:n)[L],x1[L],y1[L],sr[L],hii[L],Di[L])


# Totez pro Dataset 2
k<-length(coef(m2))
n<-length(x2)
sr<-rstudent(m2)
Lr2<-abs(sr)>2
hii<-hatvalues(m2)
Lh2<-hii>2*k/n
L<-Lr2 | Lh2
par(mfrow=c(1,1),mar=c(5,5,1.5,0)+0.05)
influencePlot(m2)
text(hii[L],sr[L],as.character((1:n)[L]),adj=c(0.5,0.5),col="blue",cex=0.75)


# Cookova vzdalenost
Di<-cooks.distance(m2)
# Vypiseme podezrela pozorovani pro Dataset 2
cat(paste("meze pro hii:  2*k/n=",round(2*k/n,4),
          "   3*k/n=",round(3*k/n,4),"\n",sep=""))
cbind((1:n)[L],x2[L],y2[L],sr[L],hii[L],Di[L])


# Indexove grafy pro Dataset 1 a 2
infIndexPlot(m1,id.method="mahal",id.n=5)

infIndexPlot(m2,id.method="mahal",id.n=5)


# Simultanni test nulovosti stredni hodnoty rezidui pro Dataset 1 a 2
outlierTest(m1)

outlierTest(m2)


# Indexove grafu pro DFBETAS, DFFIT rezidua a COVRATIO pro Dataset 1
X<-x1
M<-m1
Y<-y1
IM<-influence.measures(M)
IDs<-2:4
txtY<-c("DFBETAS_2","DFFIT","COVRATIO")
shift<-c(0,0,1)
CutOff<-c(2/sqrt(length(X)),2*sqrt(length(coef(M)/length(X))),
          3*length(coef(M)/length(X)))
par(mfrow=c(1,3),mar=c(5,5,0.5,0.5)+0.05)
for(id in 1:3)
{
ID<-IDs[id]
plot(IM$infmat[,ID],type="p",xlab="index",ylab=txtY[id])
L<-abs(IM$infmat[,ID]-shift[id])>CutOff[id]
abline(h=shift[id])
abline(h=c(shift[id]-CutOff[id],shift[id]+CutOff[id]),lty=2)
if(sum(L)>0)
{
 text(as.numeric(rownames(IM$infmat)[L]),
     IM$infmat[L,ID],labels=rownames(IM$infmat)[L],
     cex=0.75,pos=2)
 pom<-data.frame(cbind(rownames(IM$infmat)[L],X[L],Y[L],IM$infmat[L,ID]))
 names(pom)<-c("index","x","y",colnames(IM$infmat)[ID])
  o<-order(-abs(IM$infmat[L,ID]))
 print(pom[o,],digits=4)
}
}


# Indexove grafu pro DFBETAS, DFFIT rezidua a COVRATIO pro Dataset 2
X<-x2
M<-m2
Y<-y2
IM<-influence.measures(M)
IDs<-2:4
txtY<-c("DFBETAS_2","DFFIT","COVRATIO")
shift<-c(0,0,1)
CutOff<-c(2/sqrt(length(X)),2*sqrt(length(coef(M)/length(X))),
          3*length(coef(M)/length(X)))
par(mfrow=c(1,3),mar=c(5,5,0.5,0.5)+0.05)
for(id in 1:3)
{
ID<-IDs[id]
plot(IM$infmat[,ID],type="p",xlab="index",ylab=txtY[id])
L<-abs(IM$infmat[,ID]-shift[id])>CutOff[id]
abline(h=shift[id])
abline(h=c(shift[id]-CutOff[id],shift[id]+CutOff[id]),lty=2)
if(sum(L)>0)
{
 text(as.numeric(rownames(IM$infmat)[L]),
     IM$infmat[L,ID],labels=rownames(IM$infmat)[L],
     cex=0.75,pos=2)
 pom<-data.frame(cbind(rownames(IM$infmat)[L],X[L],Y[L],IM$infmat[L,ID]))
 names(pom)<-c("index","x","y",colnames(IM$infmat)[ID])
  o<-order(-abs(IM$infmat[L,ID]))
 print(pom[o,],digits=4)
}
}

# Vykreslete Pregibonuv graf