setwd("G:\\Meu Drive\\Dropbox\\2022\\Segundo Semestre\\LCE5801_Regrressão\\Aula_08")

# ajuste do modelo
m1=lm(Peso~PT+AC+AG, data=dados);anova(m1)
m2=lm(Peso~AC+AG+PT, data=dados);anova(m2)

pred=predict(m1);pred
# res?duos ordin?rios
res_ord=residuals(m1);res_ord
# res?duos estudentizados internamente
res_pad=rstandard(m1);res_pad
# res?duos estudentizados externamente
res_stud=rstudent(m1);res_stud

plot(dados)

# matriz X e vetor Y
X=model.matrix(m1);X
Y=as.matrix(dados$Peso);Y

# matriz H e diagonal principal
H=X%*%solve(t(X)%*%X)%*%t(X)
hii=diag(H);hii
n=length(Y)
p=ncol(X)
hbarra=p/n
pad_hii=(3*p)/n;pad_hii

# Pontos inconsistentes - residuo grande
res_stud
n=length(Y)
p=ncol(X)
alpha=0.05/(2*n)
t_alpha=qt(alpha,n-p-1,lower.tail=F);t_alpha

# Pontos de alavanca - leverage grande
# leverage - hii
pad_hii=(3*p)/n;pad_hii
hii

# DFBeta - estimativas dos par?metros

# DFFits - ajuste geral
pad_Fit=3*sqrt(p/(n-p));pad_Fit

# COVRATIO -raz?o de vari?ncias generalizadas
pad_Cov=cbind(1-(3*p/(n-p)),1+(3*p/(n-p)));pad_Cov

# Di - dist?ncia de Cooks - afastamento do vetor
# de estimativas com a retirada da obs i
pad_Cook=qf(0.5,p,n-p,lower.tail=F);pad_Cook


# Pontos influentes
Meddiag=influence.measures(m1);Meddiag
summary(Meddiag)


par(mar=c(2.5,2.5,2.5,2.5))
plot(m1)

# Gr?ficos de medidas de diagn?stico

  # gr?ficos de indice para hii
    plot(hii,ylim=c(0,2))
    abline(h=pad_hii,col="red")

  # gr?ficos de indice para DFFits
    dfit = dffits(m1)
    plot(dfit,ylim=c(-4,4))
    abline(h=pad_Fit,col="red")
    abline(h=-pad_Fit,col="red")

  # gr?ficos de indice para DFBetas
    dfb.1  = dfbetas(m1)[,1];
    dfb.PT = dfbetas(m1)[,2];
    dfb.AC = dfbetas(m1)[,3];
    dfb.AG = dfbetas(m1)[,4];
    par(mar=c(2.5,2.5,2.5,2.5))
    
    par(mfrow=c(2,2))
    plot(dfb.1,ylim=c(-3,3))
    abline(h=1,col="red")
    abline(h=-1,col="red")
    plot(dfb.PT,ylim=c(-3,3))
    abline(h=1,col="red")
    abline(h=-1,col="red")
    plot(dfb.AC,ylim=c(-3,3))
    abline(h=1,col="red")
    abline(h=-1,col="red")
    plot(dfb.AG,ylim=c(-3,3))
    abline(h=1,col="red")
    abline(h=-1,col="red")

  # gr?ficos de indice para Dcook
    par(mfrow=c(1,1))
    dcook =cooks.distance(m1)
    plot(dcook,ylim=c(0,3))
    abline(h=pad_Cook,col="red")


# Gr?ficos de res?duos

    # gr?fico res?duo versus Xfora

       # X1=PT fora
       m2=lm(Peso~AC+AG,data=dados)
       r2=residuals(m2)
       plot(dados$PT,r2)
       abline(h=0,col="red")
       m2.1=lm(r2~dados$PT)
       summary(m2.1)

       # X1=PT ajustada  e m2=lm(Peso~AC+AG)
       m5=lm(PT~AC+AG,data=dados)
       r5=residuals(m5)
       plot(r2,r5)
       abline(h=0,col="red")
       m5.1=lm(r2~r5)
       summary(m5.1)
       
       # X2=AC fora
       m3=lm(Peso~PT+AG,data=dados)
       r3=residuals(m3)
       m3.AC=lm(AC~PT+AG, data=dados)
       r3.ac=residuals(m3.AC)
       plot(r3,r3.ac)
       summary(lm(r3~r3.ac))
       
       # X3=AG fora
       m4=lm(Peso~PT+AC,data=dados)
       r4=residuals(m4)
       plot(dados$AG,r4)
       abline(h=0,col="red")
       m4.1=lm(r4~dados$AG)
       summary(m4.1)

    # gr?fico da vari?vel adicionada

       # X1=PT ajustada  e m2=lm(Peso~AC+AG)
       mod1=lm(Peso~PT, data=dados)
       mod2=lm(AC~PT, data=dados)
       r1=residuals(mod1)
       r2=residuals(mod2)
       plot(r2,r1)
       abline(h=0, col="red")
       m=lm(r1~r2)
       summary(m)
       plot(mod1)

       
       
       # X2 ajustada  e m3=lm(Peso~PT+AG)
       m6=lm(AC~PT+AG,data=dados)
       r6=residuals(m6)
       plot(r3,r6)
       abline(h=0,col="red")
       m6.1=lm(r3~r6)
       summary(m6.1)

       # X3 ajustada  e m4=lm(Peso~PT+AC)
       m7=lm(AG~PT+AC,data=dados)
       r7=residuals(m7)
       plot(r4,r7)
       abline(h=0,col="red")
       m7.1=lm(r4~r7)
       summary(m7.1)

   # gr?fico da vari?vel j? inclu?da

       # Fazer

   #gr?fico dos res?duos parciais (res?duo + componente)
   
       # Fazer    

    # teste de normalidade
    shapiro.test(res_stud)

    # gr?fico Normal de probabilidade
    qqnorm(res_stud,ylab="Residuos Estudentizados Externamente", main="")
    qqline(res_stud,col="red")

    # boxplot
    boxplot(res_stud)

# necessidade de transforma??o
library(MASS)

boxcox(m1,lambda = seq(-10, 30))
locator()