diff --git a/scripts/ce089-07.R b/scripts/ce089-07.R
index d91dbf43b9ab191c721dd797b90ad83f24e4ca5b..a5534011a1fdd8d5d0fd65ba5ca488144d86d977 100644
--- a/scripts/ce089-07.R
+++ b/scripts/ce089-07.R
@@ -1,17 +1,32 @@
+#=======================================================================
+# Avaliação da distribuição sob H_0 de testes de hipótese.
+
#-----------------------------------------------------------------------
# Distribuição do p-valor em um teste cuja estatÃstica de teste tem
# distribuição exata.
-simul <- function() {
- x <- rnorm(20)
- y <- rnorm(20)
- t.test(x, y, var.equal = TRUE)$p.value
+simul_diff_medias <- function(dist = "normal", n = 20) {
+ amos <- switch(dist,
+ "normal" = {
+ list(x = rnorm(n),
+ y = rnorm(n))
+ },
+ "exponencial" = {
+ list(x = rexp(n, rate = 2),
+ y = rexp(n, rate = 2))
+ },
+ "poisson" = {
+ list(x = rpois(n, lambda = 4),
+ y = rpois(n, lambda = 4))
+ })
+ t.test(amos$x, amos$y, var.equal = TRUE)$p.value
# wilcox.test(x, y)$p.value
}
-simul()
+
+simul_diff_medias()
# P-valores do teste sob H0.
-p <- replicate(1000, simul())
+p <- replicate(1000, simul_diff_medias())
# Distribuição do teste é uniforme quando é exato.
plot(ecdf(p))
@@ -22,15 +37,14 @@ abline(a = 0, b = 1, col = 2)
g <- gl(5, 4)
-simul <- function() {
+simul_f_anova <- function() {
y <- rnorm(g)
- # plot(y ~ as.integer(g))
m0 <- lm(y ~ g)
anova(m0)[1, 5]
}
# P-valores do teste sob H0.
-p <- replicate(1000, simul())
+p <- replicate(1000, simul_f_anova())
# Distribuição do teste é uniforme quando é exato.
plot(ecdf(p))
@@ -41,23 +55,25 @@ abline(a = 0, b = 1, col = 2)
g <- gl(2, 5)
-simul <- function() {
+simul_chi_deviance <- function() {
y <- rpois(g, lambda = 5)
m0 <- glm(y ~ g, family = poisson)
anova(m0, test = "Chi")[2, 5]
}
# P-valores do teste sob H0.
-p <- replicate(1000, simul())
+p <- replicate(1000, simul_chi_deviance())
# Distribuição do teste é uniforme quando é exato.
-plot(ecdf(p))
+plot(ecdf(p), cex = 0.1)
abline(a = 0, b = 1, col = 2)
-#-----------------------------------------------------------------------
+#=======================================================================
# Cobertura de intervalos de confiança.
+#-----------------------------------------------------------------------
# Gerando dados de um modelo de regressão linear simples.
+
x <- 0:10
X <- cbind(1, x)
beta <- c(3, 0.5)
@@ -68,13 +84,13 @@ plot(y ~ x)
m0 <- lm(y ~ x)
confint(m0)
-simul <- function() {
+simul_confint <- function() {
y <- X %*% beta + rnorm(x, 0, 0.3)
m0 <- lm(y ~ x)
confint(m0)
}
-ic <- replicate(100, simul())
+ic <- replicate(100, simul_confint())
str(ic)
# Com i = 1 é o intercepto e i = 2 é a taxa.
@@ -87,12 +103,12 @@ segments(1:ncol(ic[i, , ]),
abline(h = beta[i])
# Aumentando o número de simulações
-ic <- replicate(1000, simul())
+ic <- replicate(1000, simul_confint())
str(ic)
-# Conta quantos IC contem o valor paramétrico.
-sum(ic[1, 1, ] < 3 & ic[1, 2, ] > 3)
-sum(ic[2, 1, ] < 0.5 & ic[2, 2, ] > 0.5)
+# Conta quantos IC contem o valor paramétrico real.
+sum(ic[1, 1, ] < 3 & ic[1, 2, ] > 3)/dim(ic)[3]
+sum(ic[2, 1, ] < 0.5 & ic[2, 2, ] > 0.5)/dim(ic)[3]
#-----------------------------------------------------------------------
# Taxa de cobertura em modelo não linear.
@@ -137,7 +153,7 @@ layout(1)
# plot(profile(n0))
# layout(1)
-simul <- function() {
+simul_confint_nls <- function() {
y <- A * x/(V + x) + rnorm(x, 0, 0.3)
n0 <- nls(y ~ A * x/(V + x),
start = list(A = 10, V = 2))
@@ -145,7 +161,7 @@ simul <- function() {
# confint(n0)
}
-ic <- replicate(100, simul())
+ic <- replicate(100, simul_confint_nls())
str(ic)
# Com i = 1 é o intercepto e i = 2 é a taxa.
@@ -159,7 +175,7 @@ abline(h = c(A, V)[i])
# Aumentando para ter uma estimativa mais precisa.
n <- 1000
-ic <- replicate(n, simul())
+ic <- replicate(n, simul_confint_nls())
str(ic)
# Taxas de cobertura.
@@ -169,4 +185,103 @@ sum(ic[2, 1, ] < V & ic[2, 2, ] > V)/n
# TODO: substitua o IC de Wald pelo IC de perfil de verossimilhança e
# repita tudo. Tire suas próprias conclusões.
+#=======================================================================
+# Curva de poder de teste.
+
+#-----------------------------------------------------------------------
+# Teste t versus teste de Wilcox.
+
+# Várias opções de distribuição para comparar com fuga das suposições.
+simul_dados <- function(mu1 = 5, mu2 = mu1, sigma = 1, n = 10, dist = "normal") {
+ amos <- switch(dist,
+ "normal" = {
+ list(y1 = rnorm(n, mu1, sigma),
+ y2 = rnorm(n, mu2, sigma))
+ },
+ "exponencial" = {
+ list(y1 = rexp(n, rate = mu1),
+ y2 = rexp(n, rate = mu2))
+ },
+ "uniforme" = {
+ list(y1 = runif(n, min = mu1, max = mu1 + 1),
+ y2 = runif(n, min = mu2, max = mu2 + 1))
+ },
+ "poisson" = {
+ list(y1 = rpois(n, lambda = mu1),
+ y2 = rpois(n, lambda = mu2))
+ })
+ return(amos)
+}
+
+simul_t <- function(y1, y2) {
+ t.test(y1, y2, var.equal = TRUE)$p.value
+}
+
+simul_wilcox <- function(y1, y2) {
+ wilcox.test(y1, y2)$p.value
+}
+
+d_grid <- seq(0, 3, length.out = 30)
+p_grid <- sapply(d_grid,
+ FUN = function(d) {
+ alpha <- 0.05
+ N <- 1000
+ r <- replicate(N, {
+ # Avaliação pareada do poder.
+ dts <- simul_dados(mu1 = 5, mu2 = 5 + d)
+ c(simul_t(dts$y1, dts$y2) < alpha,
+ simul_wilcox(dts$y1, dts$y2) < alpha)
+ }, simplify = TRUE)
+ rowSums(r)/N
+ })
+
+p_grid <- t(p_grid)
+
+matplot(d_grid,
+ p_grid,
+ xlab = expression("Diferença entre médias:" ~ mu[1] - mu[2]),
+ ylab = expression(Pr("Rejeitar" ~ H[0] ~ "|" ~ mu[1] - mu[2])),
+ type = "l",
+ ylim = c(0, 1))
+abline(h = 0.05, col = 2)
+legend("right",
+ legend = c("Teste t", "Teste de Wilcox"),
+ lty = 1:2,
+ col = 1:2,
+ bty = "n")
+grid()
+
+#--------------------------------------------
+# Usando distribuição uniforme (fuga de pressuposto pro teste t).
+
+d_grid <- seq(0, 1.2, length.out = 30)
+p_grid <- sapply(d_grid,
+ FUN = function(d) {
+ alpha <- 0.05
+ N <- 1000
+ r <- replicate(N, {
+ # Avaliação pareada do poder.
+ dts <- simul_dados(n = 6, mu1 = 5, mu2 = 5 + d, dist = "uniforme")
+ c(simul_t(dts$y1, dts$y2) < alpha,
+ simul_wilcox(dts$y1, dts$y2) < alpha)
+ }, simplify = TRUE)
+ rowSums(r)/N
+ })
+
+p_grid <- t(p_grid)
+
+matplot(d_grid,
+ p_grid,
+ xlab = expression("Diferença entre médias:" ~ mu[1] - mu[2]),
+ ylab = expression(Pr("Rejeitar" ~ H[0] ~ "|" ~ mu[1] - mu[2])),
+ type = "l",
+ ylim = c(0, 1))
+abline(h = 0.05, col = 2)
+legend("right",
+ legend = c("Teste t", "Teste de Wilcox"),
+ lty = 1:2,
+ col = 1:2,
+ bty = "n")
+grid()
+
#-----------------------------------------------------------------------