In [1]:
# rm(list=ls()) 
options(OutDec = ",")
#==============================================================================
# Descricao
#
# Dados sobre tratamento de esgoto. 
# MANOVA: comparacoes emparelhadas
#
#==============================================================================
x       <- read.table("../dados/d04_esgoto.txt",header=T)
n       <- dim(x)[1]
p       <- 2
D       <- x[,1:2]-x[,3:4]
colnames(D) <- c("Dif_OBD","Dif_SS")
print(D)

   Dif_OBD Dif_SS
1      -19     12
2      -22     10
3      -18     42
4      -27     15
5       -4     -1
6      -10     11
7      -14     -4
8       17     60
9        9     -2
10       4     10
11     -19     -7
In [2]:
#==============================================================================
# Estatísticas
#==============================================================================
Dbarra  <- apply(D,2,mean)   # vetor de medias amostrais
S       <- var(D)            # matriz de covariancias amostrais
R       <- cor(D)            # matriz de correlacoes amostrais

print("==== Médias ====")
print(Dbarra)
print("==== Covariâncias ====")
print(S)
print("==== Correlações ====")
print(R)

[1] "==== Médias ===="
  Dif_OBD    Dif_SS 
-9,363636 13,272727 
[1] "==== Covariâncias ===="
          Dif_OBD    Dif_SS
Dif_OBD 199,25455  88,30909
Dif_SS   88,30909 418,61818
[1] "==== Correlações ===="
          Dif_OBD    Dif_SS
Dif_OBD 1,0000000 0,3057682
Dif_SS  0,3057682 1,0000000
In [3]:
#==============================================================================
# Teste de hipotese
#==============================================================================
alpha   <- 0.05
T2      <- n*Dbarra%*%solve(S)%*%Dbarra
Fcrit   <- (n-1)*p*qf(1-alpha,p,n-p)/(n-p)
print(c(T2 > Fcrit)) # Logo, rejeita-se a hipótese nula.

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

[1] TRUE
[1] "==== Valor p ===="
           [,1]
[1,] 0,02082779
In [4]:
#==============================================================================
# Intervalos-T2 simultaneos
#==============================================================================
mu1.li  <- Dbarra[1] - sqrt(Fcrit*S[1,1]/n)
mu1.ls  <- Dbarra[1] + sqrt(Fcrit*S[1,1]/n)
mu2.li  <- Dbarra[2] - sqrt(Fcrit*S[2,2]/n)
mu2.ls  <- Dbarra[2] + sqrt(Fcrit*S[2,2]/n)
print(c(mu1.li,mu1.ls))  # IC para mu_D1
print(c(mu2.li,mu2.ls))  # IC para mu_D2

  Dif_OBD   Dif_OBD 
-22,45327   3,72600 
   Dif_SS    Dif_SS 
-5,700119 32,245574
In [5]:
#==============================================================================
# Intervalos de Bonferroni
#==============================================================================
tcrit    <- qt(1-alpha/(2*p),n-1)
mu1t.li  <- Dbarra[1] - tcrit*sqrt(S[1,1]/n)
mu1t.ls  <- Dbarra[1] + tcrit*sqrt(S[1,1]/n)
mu2t.li  <- Dbarra[2] - tcrit*sqrt(S[2,2]/n)
mu2t.ls  <- Dbarra[2] + tcrit*sqrt(S[2,2]/n)
print(c(mu1t.li,mu1t.ls))  # IC para mu_D1
print(c(mu2t.li,mu2t.ls))  # IC para mu_D2

   Dif_OBD    Dif_OBD 
-20,573107   1,845835 
   Dif_SS    Dif_SS 
-2,974903 29,520358