In [1]:
# rm(list=ls())
options(OutDec = ",")
#==============================================================================
# Descricao
#
# Teste de hipoteses - dados sobre transpiracao em mulheres
#
#==============================================================================
silence <- suppressPackageStartupMessages # Para omitir mensagens de alertas
silence(library(tseries)) # Para o teste de normalidade de Jarque-Bera
In [2]:
#==============================================================================
# Lendo os dados
#==============================================================================
x <- matrix(scan("../dados/d02_suor.txt"),,3,T)
colnames(x) <- c("suor","sodio","potassio")
print(head(x))
suor sodio potassio [1,] 3,7 48,5 9,3 [2,] 5,7 65,1 8,0 [3,] 3,8 47,2 10,9 [4,] 3,2 53,2 12,0 [5,] 3,1 55,5 9,7 [6,] 4,6 36,1 7,9
In [3]:
#==============================================================================
# Verificando normalidade
#==============================================================================
par(mfrow=c(2,2),bty="n")
hist(x[,1],prob=T,col="red",main="",xlab="Suor",ylab="Densidade")
hist(x[,2],prob=T,col="red",main="",xlab="Sodio",ylab="Densidade")
hist(x[,3],prob=T,col="red",main="",xlab="Potassio",ylab="Densidade")
In [4]:
#==============================================================================
# Verificando normalidade
#==============================================================================
par(mfrow=c(2,2),bty="n")
boxplot(x[,1],col="red",main="",xlab="Suor",ylim=c(0,10))
boxplot(x[,2],col="red",main="",xlab="Sodio",ylim=c(10,80))
boxplot(x[,3],col="red",main="",xlab="Potassio",ylim=c(6,16))
In [5]:
#==============================================================================
# Verificando normalidade
#==============================================================================
par(mfrow=c(2,2),bty="n")
qqnorm(x[,1],col="red",main="",pch=15,xlab="Percentil teórico",
ylab="Percentil amostral",xlim=c(-3,3),ylim=c(0,10))
qqline(x[,1])
qqnorm(x[,2],col="red",main="",pch=15,xlab="Percentil teórico",
ylab="Percentil amostral",xlim=c(-3,3),ylim=c(10,80))
qqline(x[,2])
qqnorm(x[,3],col="red",main="",pch=15,xlab="Percentil teórico",
ylab="Percentil amostral",,xlim=c(-3,3),ylim=c(6,16))
qqline(x[,3])
In [6]:
#==============================================================================
# Testando normalidade: teste de Jarque-Bera (parametrico)
#==============================================================================
jarque.bera.test(x[,1])
jarque.bera.test(x[,2])
jarque.bera.test(x[,3])
Jarque Bera Test data: x[, 1] X-squared = 0,64246, df = 2, p-value = 0,7253
Jarque Bera Test data: x[, 2] X-squared = 0,42693, df = 2, p-value = 0,8078
Jarque Bera Test data: x[, 3] X-squared = 0,95882, df = 2, p-value = 0,6191
In [7]:
#==============================================================================
# Testando normalidade: teste de Shapiro-Wilk (nao parametrico)
#==============================================================================
shapiro.test(x[,1])
shapiro.test(x[,2])
shapiro.test(x[,3])
Shapiro-Wilk normality test data: x[, 1] W = 0,97578, p-value = 0,8689
Shapiro-Wilk normality test data: x[, 2] W = 0,98584, p-value = 0,9862
Shapiro-Wilk normality test data: x[, 3] W = 0,96385, p-value = 0,6233
In [8]:
#==============================================================================
# Calculo das estatisticas
#==============================================================================
alpha <- 0.1 # nivel de significancia
p <- 3
n <- dim(x)[1] # 20
mu0 <- c(4,50,10)
xbarra <- apply(x,2,"mean")
S <- var(x)
IS <- solve(S) # inversa de S
T2 <- n*(xbarra-mu0) %*% IS %*% (xbarra-mu0)
RC <- (n-1)*p*qf(1-alpha,p,n-p)/(n-p)
print("==== Resultado do teste ====")
print(T2 > RC) # Rejeita-se H0.
qpv <- T2*(n-p)/((n-1)*p)
valorp <- 1-pf(qpv,p,n-p)
print("==== valor p do teste ====")
print(valorp)
[1] "==== Resultado do teste ===="
[,1]
[1,] TRUE
[1] "==== valor p do teste ===="
[,1]
[1,] 0,06492834
In [9]:
#==============================================================================
# Teste da razao de verossimilhanca
#==============================================================================
SIG0HAT <- matrix(0,3,3) # estimador sob H_0
for(i in 1:n)
SIG0HAT <- SIG0HAT + (x[i,]-mu0)%*%t(x[i,]-mu0)
SIG0HAT <- SIG0HAT/n
SIG1HAT <- (n-1)*S/n # estimador sob H_1
LAMBDA <- (det(SIG1HAT)/det(SIG0HAT))^(n/2)
T2B <- (n-1)*det(SIG0HAT)/det(SIG1HAT) - (n-1)
print(T2B)
[1] 9,738773