# Lê os dados e disponibiliza-o
####IMPORTANTE : Alterar o diretório####
dados<-read.csv("coracao.csv",sep=";",dec=",") 
attach(dados)

# Define Grupo como fator (se necessário)
grupo<- factor(dados$Grupo)


# Modelo linear

fit.model<- lm(Coracao~grupo,  weights = 1/Peso)

summary(fit.model)

anovatab<-anova(fit.model)

# Diagnóstico

X <- model.matrix(fit.model)
n <- nrow(X)
p <- ncol(X)
H <- X%*%solve(t(X)%*%X)%*%t(X)
h <- diag(H)
lms <- summary(fit.model)
s <- lms$sigma
r <- resid(lms)
ts <- r/(s*sqrt(1-h))
di <- (1/p)*(h/(1-h))*(ts^2)
si <- lm.influence(fit.model)$sigma
tsi <- r/(si*sqrt(1-h))
a <- max(tsi)
b <- min(tsi)
ident <- diag(n)
epsilon <- matrix(0,n,100)
e <- matrix(0,n,100)
e1 <- numeric(n)
e2 <- numeric(n)
#
for(i in 1:100){
     epsilon[,i] <- rnorm(n,0,1)
     e[,i] <- (ident - H)%*%epsilon[,i]
     u <- diag(ident - H)
     e[,i] <- e[,i]/sqrt(u)
     e[,i] <- sort(e[,i]) }
#
for(i in 1:n){
     eo <- sort(e[i,])
     e1[i] <- (eo[2]+eo[3])/2
     e2[i] <- (eo[97]+eo[98])/2 }
#
med <- apply(e,1,mean)
faixa <- range(tsi,e1,e2)
#
par(mfrow=c(2,2))
qqnorm(tsi,xlab="Percentil da N(0,1)",
ylab="Residuo Studentizado", ylim=faixa, pch=16, main="")
par(new=T)
qqnorm(e1,axes=F,xlab="",ylab="",type="l",ylim=faixa,lty=1, main="")
par(new=T)
qqnorm(e2,axes=F,xlab="",ylab="", type="l",ylim=faixa,lty=1, main="")
par(new=T)
qqnorm(med,axes=F,xlab="",ylab="",type="l",ylim=faixa,lty=2, main="")
#------------------------------------------------------------#
plot(di,xlab="Índice", ylab="Distância de Cook", pch=16)
#identify(di, n=2)
#
plot(tsi,xlab="Índice", ylab="Resíduo Studentizado",
ylim=c(b-1,a+1), pch=16)
abline(2,0,lty=2)
abline(-2,0,lty=2)
#identify(tsi, n=1)
#
plot(fitted(fit.model),tsi,xlab="Valor Ajustado", 
ylab="Resíduo Studentizado", ylim=c(b-1,a+1), pch=16)
abline(2,0,lty=2)
abline(-2,0,lty=2)
#identify(fitted(fit.model),tsi, n=1)
par(mfrow=c(1,1))

# Teste para igualdade das variâncias

library(car)

leveneTest(tsi,grupo)

# Testes de normalidade

shapiro.test(tsi)

library(nortest)

ad.test(tsi)


# Comparacoes multiplas

testeF.CB <- function(fit.model,m.C)

{

v.beta <-  cbind(fit.model$coef)

n <- nrow(model.matrix(fit.model))

e.p <- nrow(v.beta)

e.q <- nrow(m.C)

m.cov.beta <- (vcov(fit.model))

e.F <- t(m.C%*%v.beta)%*%solve(m.C%*%m.cov.beta%*%t(m.C))%*%(m.C%*%v.beta)/e.q

e.pvalor <- 1-pf(e.F,e.q,n-e.p)

cat("Estatistica F = ",round(e.F,2),"\n")

cat("pvalor = ",round(e.pvalor,4),"\n")

cat("Matriz C :","\n")

print(m.C)

}


# Media do Grupo 2 (SHAM HR) igual a do Grupo 3 (I/R NR)

m.C23 <- cbind(0,1,-1,0,0,0)

# Media do Grupo 2 (SHAM HR) igual a do Grupo 4 (I/R HR)

m.C24 <- cbind(0,1,0,-1,0,0)

# Media do Grupo 2 (SHAM HR) igual a do Grupo 5 (I/R NR NAC)

m.C25 <- cbind(0,1,0,0,-1,0)

# Media do Grupo 2 (SHAM HR) igual a do Grupo 6 (I/R HR NAC)

m.C26 <- cbind(0,1,0,0,0,-1)

# Media do Grupo 3 (I/R NR) igual a do Grupo 4 (I/R HR)

m.C34 <- cbind(0,0,1,-1,0,0)

# Media do Grupo 3 (I/R NR) igual a do Grupo 5 (I/R NR NAC)

m.C35 <- cbind(0,0,1,0,-1,0)

# Media do Grupo 3 (I/R NR) igual a do Grupo 6 (I/R HR NAC)

m.C36 <- cbind(0,0,1,0,0,-1)

# Media do Grupo 4 (I/R HR) igual a do Grupo 5 (I/R NR NAC)

m.C45 <- cbind(0,0,0,1,-1,0)

# Media do Grupo 4 (I/R HR) igual a do Grupo 6 (I/R HR NAC)

m.C46 <- cbind(0,0,0,1,0,-1)

# Media do Grupo 5 (I/R NR NAC) igual a do Grupo 6 (I/R HR NAC)

m.C56 <- cbind(0,0,0,0,1,-1)

#testeF.CB(fit.model,m.C)