# rm(list=ls()) 
options(OutDec = ",")
#==============================================================================
# Descricao
#
# Dados sobre drogas aplicadas em caes
# Teste para a igualdade de tratamentos em experimentos de medidas repetidas
#
# HA_AP -> Tratamento 1: H ausente,  CO_2 alta pressao;
# HA_BP -> Tratamento 2: H ausente,  CO_2 baixa pressao;
# HP_AP -> Tratamento 3: H presente, CO_2 alta pressao; e
# HP_BP -> Tratamento 4: H presente, CO_2 baixa pressao.
#
#==============================================================================
silence <- suppressPackageStartupMessages # Para omitir mensagens de alertas
silence(library(Matrix))

#==============================================================================
# Lendo os dados
#==============================================================================
x       <- read.table("../dados/d05_caes.txt",header=T)
n       <- dim(x)[1]
q       <- 4
xbarra  <- apply(x,2,mean)  # vetor de medias amostrais
S       <- var(x)           # matriz de covariancias amostrais
R       <- cor(x)           # matriz de correlacoes amostrais

#==============================================================================
# Estatística
#==============================================================================
print("==== Médias ====")
print(xbarra) 
print("==== Covariâncias ====")
print(S)
print("==== Correlações ====")
print(R)

[1] "==== Médias ===="
   HA_AP    HA_BP    HP_AP    HP_BP 
368,2105 404,6316 479,2632 502,8947 
[1] "==== Covariâncias ===="
         HA_AP    HA_BP    HP_AP    HP_BP
HA_AP 2819,287 3568,415 2943,497 2295,357
HA_BP 3568,415 7963,135 5303,991 4065,459
HP_AP 2943,497 5303,991 6851,316 4499,640
HP_BP 2295,357 4065,459 4499,640 4878,988
[1] "==== Correlações ===="
          HA_AP     HA_BP     HP_AP     HP_BP
HA_AP 1,0000000 0,7531192 0,6697412 0,6188932
HA_BP 0,7531192 1,0000000 0,7180816 0,6522329
HP_AP 0,6697412 0,7180816 1,0000000 0,7782622
HP_BP 0,6188932 0,6522329 0,7782622 1,0000000

#==============================================================================
# Matriz de constrastes
#
# (mu_3 + mu_4) - (mu_1 + mu_2)
# (mu_1 + mu_3) - (mu_2 + mu_4)
# (mu_1 + mu_4) - (mu_2 + mu_3)
#
#==============================================================================
C       <- matrix(1,q-1,q)
C[1,1]  <- -1                
C[1,2]  <- -1
C[2,2]  <- -1
C[2,4]  <- -1
C[3,2]  <- -1
C[3,3]  <- -1
print(C)

print("==== Verificando o posto da matriz C. ====")
print("==== Deve ter posto q-1 ====")
print(rankMatrix(C)==(q-1))

     [,1] [,2] [,3] [,4]
[1,]   -1   -1    1    1
[2,]    1   -1    1   -1
[3,]    1   -1   -1    1
[1] "==== Verificando o posto da matriz C. ===="
[1] "==== Deve ter posto q-1 ===="
[1] TRUE

#==============================================================================
# Teste de hipotese
#==============================================================================
alpha   <- 0.05
Cxbarra <- C%*%(xbarra)
CS      <- C%*%S%*%t(C)
ICS     <- solve(CS)
Fcrit   <- (n-1)*(q-1)*qf(1-alpha,q-1,n-q+1)/(n-q+1)
T2      <- n*t(Cxbarra)%*%ICS%*%(Cxbarra)
print(T2 > Fcrit) # Logo, rejeita-se a hipótese nula.

# Valor p
print("==== Valor p ====")
print(1-pf(T2*(n-q+1)/((n-1)*(q-1)),q-1,n-q+1))

     [,1]
[1,] TRUE
[1] "==== Valor p ===="
             [,1]
[1,] 3,317767e-07

#==============================================================================
# Intervalos-T2 simultaneos
#==============================================================================
mu1.li  <- Cxbarra[1] - sqrt(Fcrit*CS[1,1]/n)
mu1.ls  <- Cxbarra[1] + sqrt(Fcrit*CS[1,1]/n)
mu2.li  <- Cxbarra[2] - sqrt(Fcrit*CS[2,2]/n)
mu2.ls  <- Cxbarra[2] + sqrt(Fcrit*CS[2,2]/n)
mu3.li  <- Cxbarra[3] - sqrt(Fcrit*CS[3,3]/n)
mu3.ls  <- Cxbarra[3] + sqrt(Fcrit*CS[3,3]/n)
print(c(mu1.li,mu1.ls))  # IC para mu_C1
print(c(mu2.li,mu2.ls))  # IC para mu_C2
print(c(mu3.li,mu3.ls))  # IC para mu_C3

[1] 135,6503 282,9813
[1] -114,72708   -5,37818
[1] -78,72858  53,14964

#==============================================================================
# Intervalos de Bonferroni
#==============================================================================
tcrit    <- qt(1-alpha/(2*(q-1)),n-1)
mu1t.li  <- Cxbarra[1] - tcrit*sqrt(CS[1,1]/n)
mu1t.ls  <- Cxbarra[1] + tcrit*sqrt(CS[1,1]/n)
mu2t.li  <- Cxbarra[2] - tcrit*sqrt(CS[2,2]/n)
mu2t.ls  <- Cxbarra[2] + tcrit*sqrt(CS[2,2]/n)
mu3t.li  <- Cxbarra[3] - tcrit*sqrt(CS[3,3]/n)
mu3t.ls  <- Cxbarra[3] + tcrit*sqrt(CS[3,3]/n)
print(c(mu1t.li,mu1t.ls))  # IC para mu_C1
print(c(mu2t.li,mu2t.ls))  # IC para mu_C2
print(c(mu3t.li,mu3t.ls))  # IC para mu_C2

[1] 150,5136 268,1180
[1] -103,6956  -16,4097
[1] -65,42422  39,84528