# Dados de uma amostra de n = 10 arvores de araucaria
# mensuradas em relação ao volume Y , área basal X1,
# area basal relativa X2 e altura em pés X3.
# Exemplo extra (NãO ESTÁ NA APOSTILA)

rm(list=ls(all=TRUE))
arv=read.table("arvores.txt",h=T);arv

# GRÁFICOS DE DISPERSÃO 2 A 2
pairs(arv)

# Montando um vetor de 1's para o intercepto beta0
b0=matrix(1,10,1);b0

# Montando a matriz de delineamneto (X), formada pelo vetor de 1's e pelas 
# variáveis explanatórias X1, X2 e X3
X=as.matrix(cbind(b0,arv[,-1]));X

# Montando o vetor da variável resposta
Y=as.matrix(arv$y);Y

# Matrizes necessárias para o processo de estimação dos parâmetros

# X'X
XtX=t(X)%*%X;
XtX

# X'X inversa
XtX_inv=solve(XtX);
XtX_inv

# vetor de parâmetros
theta=XtX_inv%*%(t(X)%*%Y);
theta

# Somas de quadrados
n=length(Y)
# Soma de Quadrados de Parâmetros
SQPAR=t(theta)%*%t(X)%*%Y;
SQPAR
# Correção
C=(sum(Y))^2/n;
C
# Soma de Quadrados de Regressão
SQREG=SQPAR-C;
SQREG
# Soma de Quadrados Total Corrigida
SQTOT=(t(Y)%*%Y)-C;
SQTOT
# Soma de Quadrados do Resíduo
SQRES=SQTOT-SQREG;
SQRES

# graus de liberdade associados às SQs
# número de parâmetros
p=ncol(X)
# gl regressão
glreg = p-1
glreg
# gl resíduo
glres = n-p
glres

# Quadrados médios
QMREG=SQREG/glreg
QMREG
QMRES=SQRES/glres
QMRES

# Estatística F
Fc=QMREG/QMRES
Fc

# Pr(>Fc)
pf(Fc,glreg,glres,lower.tail=F)

# Análise de Regressão Linear Multipla utilizando a função lm

# Soma de quadrado Sequencial 

# verificar a alteração nas somas de quadrados
# os números na frente do "mod" indica a ordem de entrada das variáveis no modelo

mod123<-lm(y~x1+x2+x3, data=arv) # X1, X2 e X3
anova(mod123)
summary(mod123)

mod213<-lm(y~x2+x1+x3, data=arv) # X2, X3 e X1
anova(mod213)
summary(mod213)

mod312<-lm(y~x3+x1+x2, data=arv) # X3, X1 e X2
anova(mod312)
summary(mod312)

# Soma de Quadrados Parcial
# Precisa instalar o pacote "car"
library(car)
  Anova(mod123,type="II")
  Anova(mod213,type="II")
  Anova(mod312,type="II")

# Outra forma de obter a Soma de Quadrados Parcial  

# precisa da diagonal principal da X'X inversa
diag(XtX_inv)
# precisa das estimativas dos par?metros
theta

SQParcialX1=(theta[2])^2/XtX_inv[2,2];SQParcialX1
SQParcialX2=(theta[3])^2/XtX_inv[3,3];SQParcialX2
SQParcialX3=(theta[4])^2/XtX_inv[4,4];SQParcialX3


# Princípio do Residuo Condicional para os testes das hipóteses

  # Exemplo 1, pg 86 apostila
  #  H0: Beta_j=0 (Usando Beta_2=0)

  ###############################################################
  # Modelo Completo y=b0+b1X1+b2X2+b3X3+e
  ###############################################################

  Xc = X
  XtXc = XtX
  XtX_inv_c = XtX_inv
  theta_c = theta
  # Somas de quadrado
  SQPARc = SQPAR;
  SQRESc = SQRES
  QMRESc = QMRES
  ###############################################################
  # Modelo Reduzido y=b0+b1X1+b3X3+e
  ###############################################################

  Xr = X[,-3];Xr
  XtXr = t(Xr)%*%Xr; XtXr
  XtX_inv_r = solve(XtXr);XtX_inv_r
  theta_r = XtX_inv_r%*%t(Xr)%*%Y;theta_r
  # Somas de quadrado
  SQPARr<-t(theta_r)%*%t(Xr)%*%Y;SQPARr
  SQRESr<-(t(Y)%*%Y)-SQPARr;SQRESr

  # ESTAT?STICA F
  Fcal=((SQRESr-SQRESc)/1)/ QMRESc;Fcal
  pf(Fcal,1,glres,lower.tail=F)


  # Exemplo 2, pg 87 apostila
  #  H0: b1=b2=b

  ##################################################################
  # Modelo reduzido   y=b0+bX1+bX2+b3X3=b0+b(X1+X2)+b3X3+e
  ##################################################################
  Xr = cbind(X[,1],X[,2]+X[,3],X[,4]);Xr
  XtXr = t(Xr)%*%Xr; XtXr
  XtX_inv_r = solve(XtXr);XtX_inv_r
  theta_r = XtX_inv_r%*%t(Xr)%*%Y;theta_r
  # Somas de quadrado
  SQPARr<-t(theta_r)%*%t(Xr)%*%Y;SQPARr
  SQRESr<-(t(Y)%*%Y)-SQPARr;SQRESr

  # ESTAT?STICA F
  Fcal=((SQRESr-SQRESc)/1)/ QMRESc;Fcal
  pf(Fcal,1,glres,lower.tail=F)