Adiciona Rmd com recortes sobre métodos intensivos.

parent 2d8830f1
---
title: "Métodos computacionalmente intensivos para inferência"
author: "[Walmes Marques Zeviani](http://lattes.cnpq.br/4410617539281650)"
date: '`r Sys.Date()`'
output:
html_document:
theme: yeti
---
```{r, include = FALSE}
opts_chunk$set(cache = TRUE,
tidy = FALSE,
fig.width = 7,
fig.height = 6,
fig.align = "center",
eval.after= "fig.cap",
warning = FALSE,
error = FALSE,
message = FALSE)
```
# Testes de aleatorização
## Uma senhora toma chá
```{r}
#-----------------------------------------------------------------------
# Uma senhora toma chá.
# Quantidade de maneiras de gerar dois grupos de 4 xícaras usando 8.
n_poss <- choose(8, 4)
n_poss
# De 4 xícaras selecionadas para um grupo, X representa o número de
# acertos.
n_acertar <- sapply(4:0,
FUN = function(i) {
# (formas de acertar x em 4) * (formas de errar x em 4).
choose(4, i) * choose(4, 4 - i)
})
n_acertar
# Pr(Acerto total) = 1/70 < 5%.
cumsum(n_acertar)/n_poss
# Comprovando por simulação.
v <- rep(0:1, each = 4)
mean(replicate(n = 100000,
expr = all(sample(v) == v)))
# Essa forma de fazer a conta é mais realista mas dá na mesma.
mean(replicate(n = 100000,
expr = all(sample(v) == sample(v))))
```
## Teste para a diferença de médias
```{r}
#-----------------------------------------------------------------------
# Exemplo com teste para a diferença de médias.
# Comprimentos de crânios de cães pré-históricos.
m <- c(120, 107, 110, 116, 114, 111, 113, 117, 114, 112)
f <- c(110, 111, 107, 108, 110, 105, 107, 106, 111, 111)
plot(ecdf(m), xlim = range(c(m, f)), col = "cyan")
lines(ecdf(f), col = "magenta")
rug(m, col = "cyan")
rug(f, col = "magenta")
# Diferença de média observada.
d <- mean(m) - mean(f)
d
# Todos as combinações possíveis.
choose(n = 20, k = 10)
#--------------------------------------------
# Com todas as combinações possíveis (exaustão).
# Para construir todas as combinações possíveis.
k <- combn(x = 1:20, m = 10)
dim(k)
# Vetor com os valores dos dois grupos concatenados.
mf <- c(m, f)
# Vetor cheio de zeros.
g <- integer(20)
# Calcula a diferença para todo arranjo possível.
D <- apply(k,
MARGIN = 2,
FUN = function(i) {
# Alguns elementos do vetor convertidos para 1.
g[i] <- 1L
# Cálculo da estatística de teste.
-diff(tapply(mf, g, FUN = mean))
})
# Histograma da distribuição da estatística sib H_0.
hist(D, col = "gray50")
rug(D)
abline(v = d, col = 2)
plot(ecdf(D), cex = 0)
rug(D)
abline(v = d, col = 2)
# P-valor do teste.
2 * sum(D >= d)/length(D)
#--------------------------------------------
# Com simulação (não necessáriamente exaustivo).
# Variáveis que indentifica os grupos.
g <- rep(1:2, each = 10)
# Apenas para conferir.
cbind(g, mf)
# Replicando a diferença para grupos formados por aleatorização.
D <- replicate(9999, {
gg <- sample(g)
-diff(tapply(mf, gg, FUN = mean))
})
# Tem que juntar a estatística observada com as simuladas.
D <- c(D, d)
hist(D, col = "gray50")
rug(D)
abline(v = d, col = 2)
plot(ecdf(D), cex = 0)
rug(D)
abline(v = d, col = 2)
# P-valor do teste.
2 * sum(D >= d)/length(D)
```
## Teste para a correlação
```{r}
#-----------------------------------------------------------------------
# Teste de aleatorização para a correlação.
# N = 5 para um par de medidas.
x <- c(4.1, 8.3, 2.9, 10.8, 9.5)
y <- c(3.7, 5.1, 1.0, 7.7, 8.9)
cbind(x, y)
plot(x, y)
# Estatística de teste na amostra original.
r0 <- cor(x, y, method = "pearson")
r0
# Todas as permutações possiveis: 5! = 120 do vetor x.
X <- gtools::permutations(n = length(x), r = length(x), v = x)
str(X)
head(X)
# As 120 correlações obtidas para cada arranjo.
r <- apply(X, MARGIN = 1, FUN = cor, y = y, method = "spearman")
# P-valor do teste.
2 * sum(r >= r0)/length(r)
```
## Índice de Moran
* <http://gis.stackexchange.com/questions/161887/significance-test-for-morans-i-using-monte-carlo-simulation>.
* <https://en.wikipedia.org/wiki/Moran's_I>.
* <http://rstudio-pubs-static.s3.amazonaws.com/268058_dd37b98f25a4496b9f0a7eb2fcf892cd.html>.
* <http://rspatial.org/analysis/rst/3-spauto.html>.
```{r}
#-----------------------------------------------------------------------
# Índice de Moran para medir dependência espacial.
# Coordenadas dos eventos em uma malha regular.
x <- 1:8
y <- 1:8
# Para criar a matriz de pesos.
ind <- expand.grid(i = 1:length(x), j = 1:length(y))
f <- function(i, j) {
u <- min(3, sum(abs(ind[i, ] - ind[j, ])))
c(0, 1, sqrt(1/2), 0)[u + 1]
}
w <- matrix(0, nrow(ind), nrow(ind))
for (i in 1:nrow(ind)) {
for (j in 1:nrow(ind)) {
w[i, j] <- f(i, j)
}
}
w <- w/sum(w)
image(w)
# Índice de Moran.
moran <- function(x, weights) {
# Tamanho da amostra.
n <- length(x)
# Valores normalizados.
z <- as.vector((x - mean(x))/sd(x))
# Índice de Moran.
as.vector(z %*% weights %*% (z * sqrt(n/(n - 1))))
}
# Teste de permutação com gráfico.
ppt <- function(z, w, N = 10000, stat, ...) {
# Índice de Moran por reamostragem.
sim <- replicate(N,
moran(sample(z), w))
# Determina o p-valor.
p.value <- mean((all <- c(stat, sim)) >= stat)
# Histograma da distribuição empírica sob H_0.
hist(sim,
sub = paste("p =", round(p.value, 4)),
xlim = range(all),
...)
abline(v = stat, col = "#903030", lty = 3, lwd = 2)
return(p.value)
}
# Observações simuladas.
set.seed(17)
par(mfrow = c(2, 3))
# Dados com dependência espacial ---------------------------------------
# Indução de autocorrelação por meio de uma tendência.
z <- matrix(rexp(length(x) * length(y),
outer(x, y^2)),
length(x))
image(log(z), main = "Com dependência")
# Índice de Moran com dados originais.
stat <- moran(z, w)
stat
hist(z)
ppt(z, w, stat = stat, main = "I de Moran", xlab = "I")
# Dados sem dependência espacial ---------------------------------------
# Geração de de um conjunto de dados sob hipótese nula.
z <- matrix(rnorm(length(x) * length(y), 0, 1/2), length(x))
image(z, main = "Sem dependência")
# Índice de Moran com dados originais.
stat <- moran(z, w)
stat
hist(z)
ppt(z, w, stat = stat, main = "I de Moran", xlab = "I")
```
# Jackknife
## Regressão não linear
```{r}
#-----------------------------------------------------------------------
# Aplicação de Jackknife em modelos de regressão não linear.
library(latticeExtra)
library(car)
library(alr3)
# Curva "palpite".
start <- list(th0 = 75, th1 = 0.5, th2 = 50)
xyplot(C ~ Temp, data = segreg) +
layer(panel.curve(th0 + th1 * (x - th2) * (x >= th2) +
0 * (x < th2), lty = 2),
data = start)
# Ajuste.
n0 <- nls(C ~ th0 + th1 * (Temp - th2) * (Temp >= th2) +
0 * (Temp < th2),
data = segreg,
start = start)
# Estimativas e medidas de ajuste.
summary(n0)
# Observados e preditos.
xyplot(C ~ Temp, data = segreg) +
layer(panel.curve(th0 + th1 * (x - th2) * (x >= th2) +
0 * (x < th2)),
data = as.list(coef(n0)))
#-----------------------------------------------------------------------
# Estimativas Jackknife para os parâmetros.
theta <- coef(n0) # Estimativas com todos os dados.
n <- nrow(segreg) # Número de observações.
# Matriz vazia para armazenar os pseudo-valores.
thetaj <- matrix(0, nrow = n, ncol = length(theta))
colnames(thetaj) <- names(theta)
# Ajustes deixando uma observação de fora (leave-one-out).
for (i in 1:n) {
# Reajusta o modelo, i.e. estima com leave-one-out.
n1 <- nls(C ~ th0 + th1 * (Temp - th2) * (Temp >= th2) +
0 * (Temp < th2),
data = segreg[-i, ],
start = coef(n0))
# Calcula os pseudo valores.
thetaj[i, ] <- n * theta - (n - 1) * coef(n1)
}
# Os primeiros pseudo valores.
head(thetaj)
# Estimativas pontuais.
jk <- colMeans(thetaj)
cbind(MLE = theta, Jack = jk, Vício = theta - jk)
# Erros padrões.
cbind(MLE = summary(n0)$coefficients[, "Std. Error"],
Jack = apply(thetaj, MARGIN = 2, sd)/sqrt(n))
# Pacote com função pronta para Jackknife para modelos não lineares.
library(nlstools)
# ls("package:nlstools")
# Aplica Jackknife sobre o modelo.
j0 <- nlsJack(n0)
summary(j0)
```
## Dose letal para 50% da população
```{r}
#-----------------------------------------------------------------------
# Inferência para a DL 50.
# library(labestData)
# data(PaulaTb3.12)
# str(PaulaTb3.12)
# help(PaulaTb3.12, h = "html")
# # Converte a resposta para binário.
# PaulaTb3.12$resp <- as.integer(PaulaTb3.12$resp) - 1
# dput(PaulaTb3.12)
PaulaTb3.12 <-
structure(list(vol = c(3.7, 3.5, 1.25, 0.75, 0.8, 0.7, 0.6, 1.1,
0.9, 0.9, 0.8, 0.55, 0.6, 1.4, 0.75, 2.3, 3.2, 0.85, 1.7, 1.8,
0.4, 0.95, 1.35, 1.5, 1.6, 0.6, 1.8, 0.95, 1.9, 1.6, 2.7, 2.35,
1.1, 1.1, 1.2, 0.8, 0.95, 0.75, 1.3), razao = c(0.825, 1.09,
2.5, 1.5, 3.2, 3.5, 0.75, 1.7, 0.75, 0.45, 0.57, 2.75, 3, 2.33,
3.75, 1.64, 1.6, 1.415, 1.06, 1.8, 2, 1.36, 1.35, 1.36, 1.78,
1.5, 1.5, 1.9, 0.95, 0.4, 0.75, 0.03, 1.83, 2.2, 2, 3.33, 1.9,
1.9, 1.625), resp = c(1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0,
1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1,
1, 1, 0, 0, 1)), .Names = c("vol", "razao", "resp"), row.names = c(NA,
39L), class = "data.frame")
layout(1)
plot(resp ~ vol, data = PaulaTb3.12)
m0 <- glm(resp ~ vol,
data = PaulaTb3.12,
family = binomial)
summary(m0)
# DL_50.
dl <- -coef(m0)[1]/coef(m0)[2]
# str(m0$family)$link
plot(resp ~ vol, data = PaulaTb3.12)
curve(m0$family$linkinv(coef(m0)[1] + coef(m0)[2] * x),
add = TRUE)
abline(v = dl, h = 0.5, lty = 2)
n <- nrow(PaulaTb3.12)
i <- 1:n
# Estimativas por Jackknife.
DL <- sapply(i,
FUN = function(i) {
m0 <- glm(resp ~ vol,
data = PaulaTb3.12[-i, ],
family = binomial)
# Estimativa Parcial.
DL <- -coef(m0)[1]/coef(m0)[2]
# Pseudo-valor.
n * dl - (n - 1) * DL
})
m <- mean(DL)
e <- sd(DL)/sqrt(n)
cbind(Est = rbind(MLE = dl,
Jack = m),
StdErr = rbind(MLE = car::deltaMethod(m0, g = "-(Intercept)/vol")["SE"],
Jack = e))
plot(ecdf(DL))
rug(DL)
plot(density(DL))
rug(DL)
# NOTE: alguma pista sobre porque a distribuição fica bimodal?
```
## Uma situação problemática
```{r}
#-----------------------------------------------------------------------
# Aplicação na Mediana (derivada não suave) -> problema.
# Ordena o vetor de valores.
x <- sort(precip)
# Calcula a mediana com os dados originais.
M <- median(x)
M
n <- length(x)
i <- 1:length(x)
# Estimativas da mediana por Jackknife.
y <- sapply(i,
FUN = function(i) {
ep <- median(x[-i])
pv <- n * M - (n - 1) * ep
return(pv)
})
stem(y)
# NOTE: Explique porque o Jackknife para a mediana só retornou 2
# valores? Considere que o tamanho da amostra é dado abaixo.
length(x)
```
# Bootstrap
## Regressão não linear
```{r}
#-----------------------------------------------------------------------
# Aplicação de bootstrap. Estudo de caso com modelos não lineares.
library(alr3)
str(turk0)
# Gráfico dos valores observados.
plot(Gain ~ A, data = turk0)
# Lista com os valores iniciais.
start <- list(Int = 625, Ass = 180, Mei = 0.1)
# Adiciona as curvas.
with(start, {
curve(Int + Ass * A/(Mei + A),
xname = "A", add = TRUE, col = 2)
curve(Int + Ass * (1 - 2^(-A/Mei)),
xname = "A", add = TRUE, col = 4)})
# Ajuste do primeiro modelo.
n0 <- nls(Gain ~ Int + Ass * A/(Mei + A),
data = turk0,
start = start)
# Ajuste do segundo modelo.
n1 <- nls(Gain ~ Int + Ass * (1 - 2^(-A/Mei)),
data = turk0,
start = start)
cbind(coef(n0), coef(n1))
#-----------------------------------------------------------------------
# Criando funções que fazem reamostragem dos dados e ajustam o modelo.
n <- nrow(turk0)
i <- 1:n
myboot0 <- function(formula, start) {
n0 <- nls(formula = formula,
data = turk0[sample(i, n, replace = TRUE), ],
start = start)
coef(n0)
}
#-----------------------------------------------------------------------
# Replicando.
# B = 999 reamostragens.
b0 <- replicate(999,
myboot0(Gain ~ Int + Ass * A/(Mei + A), start))
b1 <- replicate(999,
myboot0(Gain ~ Int + Ass * (1 - 2^(-A/Mei)), start))
# Conjunto de 1000 estimativas bootstrap.
b0 <- cbind(b0, coef(n0))
b1 <- cbind(b1, coef(n1))
str(b0)
# Calcula o vício absoluto (bootstrap - original).
rbind(n0 = rowMeans(b0) - coef(n0),
n1 = rowMeans(b1) - coef(n1))
# Variância do estimador.
rbind(n0 = apply(b0, MARGIN = 1, var),
n1 = apply(b1, MARGIN = 1, var))
# Função para calcular o erro quadrático médio.
eqm <- function(x) {
sum((x - tail(x, 1))^2)/length(x)
}
# Erro quadrático médio.
rbind(n0 = apply(b0, 1, FUN = eqm),
n1 = apply(b1, 1, FUN = eqm))
```
## Inferência para valores preditos
```{r}
#-----------------------------------------------------------------------
# Curvas ajustadas de onde é possível determinar uma banda de confiança.
bts <- cbind(rep(1:2, each = ncol(b0)),
as.data.frame(rbind(t(b0), t(b1))))
splom(bts[, -1], groups = bts[, 1])
par(mfrow = c(1, 2))
plot(Gain ~ A, data = turk0, type = "n")
apply(b0,
MARGIN = 2,
FUN = function(b) {
curve(b[1] + b[2] * x/(b[3] + x),
add = TRUE, col = rgb(0, 1, 1, 0.25), n = 51)
invisible()
})
points(Gain ~ A, data = turk0)
abline(v = 0.2, lty = 2)
plot(Gain ~ A, data = turk0, type = "n")
apply(b1,
MARGIN = 2,
FUN = function(b) {
curve(b[1] + b[2] * (1 - 2^(-x/b[3])),
add = TRUE, col = rgb(1, 1, 0, 0.25), n = 51)
invisible()
})
points(Gain ~ A, data = turk0)
abline(v = 0.2, lty = 2)
#-----------------------------------------------------------------------
# Estimativa do valor da função em x = 0.2.
p0 <- apply(b0, MARGIN = 2,
FUN = function(b){
b[1] + b[2] * 0.2/(b[3] + 0.2)
})
p1 <- apply(b1, MARGIN = 2,
FUN = function(b){
b[1] + b[2] * (1 - 2^(-0.2/b[3]))
})
cbind(Est = rbind(mean(p0), mean(p1)),
StdErr = rbind(sd(p0), sd(p1)))
```
## Usando o pacote `boot`
```{r}
#-----------------------------------------------------------------------
# Usando o pacote boot.
library(boot)
library(latticeExtra)
# dput(as.list(coef(n1)))
start <- list(Int = 622.958054146589,
Ass = 178.251908347242,
Mei = 0.097321927254495)
# Criar uma função com dois argumentos: o data.frame original e um vetor
# que vai representar o índice das linhas usado para reamostrar dos
# dados.
fitmodel <- function(dataset, index) {
n0 <- nls(Gain ~ Int + Ass * (1 - 2^(-A/Mei)),
data = dataset[index, ],
start = start)
c(coef(n0), s2 = deviance(n0)/df.residual(n0))
}
set.seed(321)
b0 <- boot(data = turk0,
statistic = fitmodel,
R = 1999)
class(b0)
methods(class = class(b0))
summary(b0)
# Como são obtidos.
theta <- b0$t
# Estimativas como conjunto de dados original.
b0$t0
# Vício.
apply(theta, MARGIN = 2, FUN = mean) - b0$t0
# Desvio-padrão boostrap.
apply(theta, MARGIN = 2, FUN = sd)
# Mediana bootstrap.
apply(theta, MARGIN = 2, FUN = median)
st <- summary(b0)
str(st)
# Gráfico de pares.
pairs(b0)
vcov(b0)
cov2cor(vcov(b0))
# Extrair os vetores de estativativas bootstrap.
str(b0)
colnames(b0$t) <- names(b0$t0)
best <- stack(as.data.frame(b0$t))
names(best)
st
densityplot(~values | ind,
data = best,
scales = "free",
as.table = TRUE) +
layer(panel.abline(v = st$original[which.packet()] +
c(0, st$bootBias[which.packet()]),
col = c(1, 2),
lty = 2))
# Intervalos de confiança.
confint(b0, type = "norm")
confint(b0, type = "basic")
confint(b0, type = "perc")
confint(b0, type = "bca")
```
# Monte Carlo
## Processo pontual
```{r}
#-----------------------------------------------------------------------
# Teste para independência de processo pontual.
plot(NULL, NULL, xlim = c(0, 1), ylim = c(0, 1), asp = 1)
lines(x = c(0, 1, 1, 0, 0), y = c(0, 0, 1, 1, 0))
# xy <- locator(n = 20, type = "p", pch = 19)
# dput(lapply(xy, round, digits = 3))
xy <- structure(list(x = c(0.204, 0.186, 0.529, 0.529, 0.385, 0.579,
0.918, 0.798, 0.634, 0.761, 0.704, 0.485,
0.291, 0.341, 0.402, 0.833, 0.972, 0.625,
0.732, 0.315),
y = c(0.829, 0.545, 0.526, 0.752, 0.674, 0.648,
0.792, 0.33, 0.121, 0.127, 0.332, 0.352,
0.188, 0.452, 0.221, 0.524, 0.957, 0.964,
0.755, 0.332)),
.Names = c("x", "y"))
#-----------------------------------------------------------------------
# Transforma a lista em matriz.
xy <- do.call(cbind, xy)
plot(xy,
NULL,
xlim = c(0, 1),
ylim = c(0, 1),
asp = 1,
axes = FALSE,
ann = FALSE)
lines(x = c(0, 1, 1, 0, 0),
y = c(0, 0, 1, 1, 0),
lty = 2)
# Calcula todas as distâncias (euclidianas).
d <- dist(xy)
d
# Obtém a menor distância.
m <- min(d)
# Número de pontos.
n <- nrow(xy)
# Geração de estatísticas por simulação do modelo assumido sob
# hipótese nula: localização uniforme no quadrado.
M <- replicate(9999, {
# As localizações de n pontos.
loc <- cbind(x = runif(n),
y = runif(n))
# A menor distância entre eles sob H_0.
min(dist(loc))
})
# Concatena as estatísticas simuladas com a observada.
M <- c(m, M)
# Gráfico da distribuição acumulada empírica sob H_0.
plot(ecdf(M))
abline(h = c(0.025, 0.975), lty = 2, col = 2)
abline(v = m, col = 2)
# Gráfico da distribuição acumulada empírica sob H_0.
plot(density(M))
abline(v = m, col = 2)
abline(v = quantile(M, c(0.025, 0.975)), lty = 2, col = 2)
rug(M)
# P-valor.
2 * sum(M > m)/length(M)
#-----------------------------------------------------------------------
# Moficando a estatística de teste.
xy <- structure(list(x = c(0.088, 0.326, 0.577, 0.846, 0.857, 0.568,
0.306, 0.077, 0.077, 0.328, 0.597, 0.883,
0.863, 0.64, 0.337, 0.088, 0.077, 0.346,
0.654, 0.619),
y = c(0.92, 0.922, 0.916, 0.935, 0.737, 0.674,
0.67, 0.665, 0.452, 0.502, 0.461, 0.454,
0.256, 0.26, 0.26, 0.219, 0.045, 0.06, 0.058,
0.439)),
.Names = c("x", "y"))
# Transforma a lista em matriz.
xy <- do.call(cbind, xy)
plot(xy,
NULL,
xlim = c(0, 1),
ylim = c(0, 1),
asp = 1,
axes = FALSE,
ann = FALSE)
lines(x = c(0, 1, 1, 0, 0),
y = c(0, 0, 1, 1, 0),
lty = 2)
# Todas as distância entre os pontos (a estatística de teste é um vetor
# e não um escalar).
d <- c(dist(xy))
d
# Gerando estatísticas de teste por simulação.
D <- replicate(499, {
dist(cbind(x = runif(n), y = runif(n)))
})
# Juntanto observado com simulado.
D <- cbind(d, D)
str(D)
plot(ecdf(D[, 1]), cex = 0)
for (i in 2:ncol(D)) {
lines(ecdf(D[, i]), cex = 0, col = "gray50")
}
lines(ecdf(D[, 1]), cex = 0, col = 2, lwd = 2)
```
## Amostragem por conjuntos ordenados
```{r}
#-----------------------------------------------------------------------
# Amostragem por Conjuntos Ordenados.