##############################
# Leitura dos dados
##############################
dados=read.table("cov.R",h=T)
attach(dados)
dados$maq<-as.factor(dados$maq) 
str(dados)

# Gráfico que relaciona a variável resposta e a covariável
with(dados, plot(resp ~ cov, col = c(1, 2, 4)[maq], pch = (1:3)[maq], xlab="Diâmetro (cm)", ylab="Comprimento (cm)"))
##############################


##############################
# Estimação dos parâmetros

a=3
n=5

# Matriz delineamento
col1=c(rep(1,n*a))
col2=c(rep(1,n), rep(0,2*n))
col3=c(rep(0,n), rep(1,n), rep(0,n))
col4=c(rep(0,2*n), rep(1,n))
Z=cbind(col1,col2,col3,col4)

# Vetor de covariáveis e variável dependente
X=dados$cov
Y=dados$resp

library(MASS)

# Matriz P
P=Z%*%ginv(crossprod(Z))%*%t(Z)
I=diag(a*n)

# Soma de quadrados e produtos dos resíduos
Rx2=t(X)%*%(I-P)%*%X
Rxy=t(X)%*%(I-P)%*%Y
Ry2=t(Y)%*%(I-P)%*%Y

# Estimativa de Beta
beta=solve(Rx2)%*%Rxy

# Estimativa de Alfa
alfa=ginv(crossprod(Z))%*%t(Z)%*%Y - ginv(crossprod(Z))%*%t(Z)%*%X%*%beta

# Coeficiente de correlação
r=Rxy/sqrt(Rx2*Ry2)
tc=(r/sqrt(1-r^2))*sqrt(a*(n-1)-1)
1-pt(tc,(a*(n-1)-1))
##############################


##############################
# Teste para efeito de tratamentos ajustado

# Soma de Quadrados de Resíduo (modelo reduzido)
SQR.aj=Ry2 - ((Rxy^2)/Rx2)

# Soma de Quadrados Total (modelo completo)
modelo.r=lm(resp~cov)
h=model.matrix(modelo.r)
H=h%*%solve(crossprod(h))%*%t(h)
SQT.aj=t(Y)%*%(I-H)%*%Y

# Soma de quadrados Tratamento ajustado
SQA=SQT.aj-SQR.aj
F1=(SQA/(a-1))/(SQR.aj/ ((a*(n-1))-1))
##############################


##############################
# Teste sobre a inclinação
SQRL=(((Rxy)^2)/Rx2)
F2=SQRL/(SQR.aj/(a*(n-1)-1))
##############################

##############################
# Teste para homogeneidade de inclinações
X2=matrix(c(dados[1:5,2], rep(0,a*n), dados[6:10,2], rep(0,a*n), dados[11:15,2]),a*n, byrow=F)
R1x2=t(X2)%*%(I-P)%*%X2
R1xy=t(X2)%*%(I-P)%*%Y
R1y2=t(Y)%*%(I-P)%*%Y
beta1=solve(R1x2)%*%R1xy
SQF=R1y2-(t(R1xy)%*%beta1)
F3=((SQR.aj-SQF)/2)/(SQF/(a*(n-2)))
qf(0.95,2,9)
1-pf(F3,2,9)


##############################
# UTILIZANDO FUNÇÕES PRONTAS
##############################

# Teste sobre a inclinação da reta de regressão
dados.cov<-aov(resp ~ maq + cov, data = dados) 
summary(dados.cov)
##############################


##############################
# Teste sobre o efeito do tratamento ajustado 
dados.trat<-aov(resp ~ cov + maq, data = dados) 
summary(dados.trat)
##############################


##############################
# Análise de Covariância
require(car) 
Anova(dados.cov, type = "III")
# ou
Anova(dados.trat, type = "III")
##############################