Acrescenta conteúdo ao script de Monte Carlo.

parent f433011b
#=======================================================================
# 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 do p-valor em um teste cuja estatística de teste tem
# distribuição exata. # distribuição exata.
simul <- function() { simul_diff_medias <- function(dist = "normal", n = 20) {
x <- rnorm(20) amos <- switch(dist,
y <- rnorm(20) "normal" = {
t.test(x, y, var.equal = TRUE)$p.value 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 # wilcox.test(x, y)$p.value
} }
simul()
simul_diff_medias()
# P-valores do teste sob H0. # P-valores do teste sob H0.
p <- replicate(1000, simul()) p <- replicate(1000, simul_diff_medias())
# Distribuição do teste é uniforme quando é exato. # Distribuição do teste é uniforme quando é exato.
plot(ecdf(p)) plot(ecdf(p))
...@@ -22,15 +37,14 @@ abline(a = 0, b = 1, col = 2) ...@@ -22,15 +37,14 @@ abline(a = 0, b = 1, col = 2)
g <- gl(5, 4) g <- gl(5, 4)
simul <- function() { simul_f_anova <- function() {
y <- rnorm(g) y <- rnorm(g)
# plot(y ~ as.integer(g))
m0 <- lm(y ~ g) m0 <- lm(y ~ g)
anova(m0)[1, 5] anova(m0)[1, 5]
} }
# P-valores do teste sob H0. # P-valores do teste sob H0.
p <- replicate(1000, simul()) p <- replicate(1000, simul_f_anova())
# Distribuição do teste é uniforme quando é exato. # Distribuição do teste é uniforme quando é exato.
plot(ecdf(p)) plot(ecdf(p))
...@@ -41,23 +55,25 @@ abline(a = 0, b = 1, col = 2) ...@@ -41,23 +55,25 @@ abline(a = 0, b = 1, col = 2)
g <- gl(2, 5) g <- gl(2, 5)
simul <- function() { simul_chi_deviance <- function() {
y <- rpois(g, lambda = 5) y <- rpois(g, lambda = 5)
m0 <- glm(y ~ g, family = poisson) m0 <- glm(y ~ g, family = poisson)
anova(m0, test = "Chi")[2, 5] anova(m0, test = "Chi")[2, 5]
} }
# P-valores do teste sob H0. # P-valores do teste sob H0.
p <- replicate(1000, simul()) p <- replicate(1000, simul_chi_deviance())
# Distribuição do teste é uniforme quando é exato. # Distribuição do teste é uniforme quando é exato.
plot(ecdf(p)) plot(ecdf(p), cex = 0.1)
abline(a = 0, b = 1, col = 2) abline(a = 0, b = 1, col = 2)
#----------------------------------------------------------------------- #=======================================================================
# Cobertura de intervalos de confiança. # Cobertura de intervalos de confiança.
#-----------------------------------------------------------------------
# Gerando dados de um modelo de regressão linear simples. # Gerando dados de um modelo de regressão linear simples.
x <- 0:10 x <- 0:10
X <- cbind(1, x) X <- cbind(1, x)
beta <- c(3, 0.5) beta <- c(3, 0.5)
...@@ -68,13 +84,13 @@ plot(y ~ x) ...@@ -68,13 +84,13 @@ plot(y ~ x)
m0 <- lm(y ~ x) m0 <- lm(y ~ x)
confint(m0) confint(m0)
simul <- function() { simul_confint <- function() {
y <- X %*% beta + rnorm(x, 0, 0.3) y <- X %*% beta + rnorm(x, 0, 0.3)
m0 <- lm(y ~ x) m0 <- lm(y ~ x)
confint(m0) confint(m0)
} }
ic <- replicate(100, simul()) ic <- replicate(100, simul_confint())
str(ic) str(ic)
# Com i = 1 é o intercepto e i = 2 é a taxa. # Com i = 1 é o intercepto e i = 2 é a taxa.
...@@ -87,12 +103,12 @@ segments(1:ncol(ic[i, , ]), ...@@ -87,12 +103,12 @@ segments(1:ncol(ic[i, , ]),
abline(h = beta[i]) abline(h = beta[i])
# Aumentando o número de simulações # Aumentando o número de simulações
ic <- replicate(1000, simul()) ic <- replicate(1000, simul_confint())
str(ic) str(ic)
# Conta quantos IC contem o valor paramétrico. # Conta quantos IC contem o valor paramétrico real.
sum(ic[1, 1, ] < 3 & ic[1, 2, ] > 3) sum(ic[1, 1, ] < 3 & ic[1, 2, ] > 3)/dim(ic)[3]
sum(ic[2, 1, ] < 0.5 & ic[2, 2, ] > 0.5) sum(ic[2, 1, ] < 0.5 & ic[2, 2, ] > 0.5)/dim(ic)[3]
#----------------------------------------------------------------------- #-----------------------------------------------------------------------
# Taxa de cobertura em modelo não linear. # Taxa de cobertura em modelo não linear.
...@@ -137,7 +153,7 @@ layout(1) ...@@ -137,7 +153,7 @@ layout(1)
# plot(profile(n0)) # plot(profile(n0))
# layout(1) # layout(1)
simul <- function() { simul_confint_nls <- function() {
y <- A * x/(V + x) + rnorm(x, 0, 0.3) y <- A * x/(V + x) + rnorm(x, 0, 0.3)
n0 <- nls(y ~ A * x/(V + x), n0 <- nls(y ~ A * x/(V + x),
start = list(A = 10, V = 2)) start = list(A = 10, V = 2))
...@@ -145,7 +161,7 @@ simul <- function() { ...@@ -145,7 +161,7 @@ simul <- function() {
# confint(n0) # confint(n0)
} }
ic <- replicate(100, simul()) ic <- replicate(100, simul_confint_nls())
str(ic) str(ic)
# Com i = 1 é o intercepto e i = 2 é a taxa. # Com i = 1 é o intercepto e i = 2 é a taxa.
...@@ -159,7 +175,7 @@ abline(h = c(A, V)[i]) ...@@ -159,7 +175,7 @@ abline(h = c(A, V)[i])
# Aumentando para ter uma estimativa mais precisa. # Aumentando para ter uma estimativa mais precisa.
n <- 1000 n <- 1000
ic <- replicate(n, simul()) ic <- replicate(n, simul_confint_nls())
str(ic) str(ic)
# Taxas de cobertura. # Taxas de cobertura.
...@@ -169,4 +185,103 @@ sum(ic[2, 1, ] < V & ic[2, 2, ] > V)/n ...@@ -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 # TODO: substitua o IC de Wald pelo IC de perfil de verossimilhança e
# repita tudo. Tire suas próprias conclusões. # 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()
#----------------------------------------------------------------------- #-----------------------------------------------------------------------
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment