dados=read.csv2("novilho.csv",h=T)

# ajuste do modelo
Peso=as.numeric(dados[,4])
PT=as.numeric(dados[,1])
AC=as.numeric(dados[,2])
AG=as.numeric(dados[,3])
m1=lm(Peso~PT+AC+AG)
anova(m1)
pred=predict(m1);pred
res_ord=residuals(m1);res_ord
res_pad=rstandard(m1);res_pad
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

# 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)

# 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(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
    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)
       r2=residuals(m2)
       plot(PT,r2)
       abline(h=0,col="red")
       m2.1=lm(r2~PT)
       summary(m2.1)

       # X2=AC fora
       m3=lm(Peso~PT+AG)
       r3=residuals(m3)
       plot(AC,r3)
       abline(h=0,col="red")
       m3.1=lm(r3~AC)
       summary(m3.1)

       # X3=AG fora
       m4=lm(Peso~PT+AC)
       r4=residuals(m4)
       plot(AG,r4)
       abline(h=0,col="red")
       m4.1=lm(r4~AG)
       summary(m4.1)

    # gráfico da variável adicionada

       # X1=PT ajustada  e m2=lm(Peso~AC+AG)
       m5=lm(PT~AC+AG)
       r5=residuals(m5)
       plot(r2,r5)
       abline(h=0,col="red")
       m5.1=lm(r2~r5)
       summary(m5.1)

       # X2 ajustada  e m3=lm(Peso~PT+AG)
       m6=lm(AC~PT+AG)
       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)
       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(-5, 5))
locator()