Adiciona tutorial de Jackknife.

parent d2cd0ed9
......@@ -41,8 +41,10 @@ navbar:
href: tutoriais/04-rel-va.html
- text: "GNA para Normal - Box Muller"
href: tutoriais/05-box-muller.html
- text: "Jackknife"
- text: "Jackknife (slides)"
href: slides/07-jackknife.pdf
- text: "Jackknife (tutorial)"
href: tutoriais/08-jackknife.html
- text: "----------"
- text: "Arquivos complementares"
- text: "GNA Uniformes (2015)"
......
---
title: "Aplicações de Jackknife"
author: Prof. Walmes M. Zeviani
date: '`r Sys.Date()`'
bibliography: ../config/Refs.bib
csl: ../config/ABNT-UFPR-2011-Mendeley.csl
---
# Exemplo com correlação
Neste exemplo, será usado Jackknife para fazer inferência sobre o
coeficiente de correlação de Pearson.
```{r}
# Fertilidade (y) vs percentual de homens ocupados na agricultura (x).
plot(Fertility ~ Agriculture, data = swiss)
# O coeficiente de correlação de Pearson possui método de
# inferência. Porém, depende do atendimento dos pressupostos, etc.
with(swiss, cor.test(x = Fertility, y = Agriculture))
# Estimativa com todas as observações.
rho <- with(swiss, cor(x = Fertility, y = Agriculture))
rho
n <- nrow(swiss) # Número de observações.
i <- 1:n # Índices.
# Estimativas parciais (partial estimates).
pe <- sapply(i,
FUN = function(ii) {
with(swiss[-ii, ],
cor(x = Fertility, y = Agriculture))
})
sort(pe)
# Pseudo valores.
pv <- n * rho - (n - 1) * pe
sort(pv)
# ATENÇÃO: alguns pseudo valores tem valores fora do espaço paramétrico
# da correlação que é [-1, 1].
# Estimativa pontual Jackknife.
mean(pv)
# Erro-padrão Jackknife.
sqrt(var(pv)/n)
# Intervalo de confiança Jackknife (supõe independência e normalidade).
mean(pv) + c(-1, 1) * qt(p = 0.975, df = n - 1) * sqrt(var(pv)/n)
# Intervalo de confiança para a correlação de Pearson.
with(swiss,
cor.test(x = Fertility, y = Agriculture)$conf.int)
```
Em qual dos dois tipos de intervalo de confiança é mais seguro (para não
usar o termo confiável)? O que é ser mais seguro?
```{r, include = FALSE}
library(mvtnorm)
pearson_jackknife <- function(x, y) {
n <- length(x)
i <- 1:n
rho <- cor(x, y)
pe <- sapply(i,
FUN = function(ii) {
cor(x[-ii], y[-ii])
})
pv <- n * rho - (n - 1) * pe
mean(pv) + c(0, -1, 1) * qt(p = 0.975, df = n - 1) * sqrt(var(pv)/n)
}
pearson <- function(x, y) {
ctest <- cor.test(x, y)
c(ctest$estimate, ctest$conf.int)
}
rho <- 0.5
i <- 0
res <- replicate(3000,
simplify = FALSE, {
xy <- rmvnorm(n = 50, mean = c(0, 0), sigma = matrix(c(1, rho, rho, 1), 2, 2))
# out <- data.frame(param = c("est", "lwr", "upr"),
# pearson = pearson(xy[, 1], xy[, 2]),
# jackknife = pearson_jackknife(xy[, 1], xy[, 2]))
out <- rbind(pearson = pearson(xy[, 1], xy[, 2]),
jackknife = pearson_jackknife(xy[, 1], xy[, 2]))
colnames(out) <- c("est", "lwr", "upr")
i <<- i + 1
cbind(id = i, as.data.frame(out))
})
res[[1]]
res <- do.call(rbind, res)
res$met <- c("pearson", "jackknife")
str(res)
library(tidyverse)
res <- res %>%
gather(key = "param", value = "valor", est, lwr, upr)
head(res)
ggplot(filter(res, id <= 1000), aes(x = id, y = valor, color = param)) +
geom_point() +
# stat_summary(aes(group = param), fun.y = mean, geom = "line") +
geom_smooth(se = FALSE, method = lm, formula = y ~ 1) +
facet_wrap(~met)
ggplot(filter(res, id <= 50), aes(x = id, y = valor, color = met)) +
geom_point() +
facet_wrap(~param)
res %>%
group_by(param, met) %>%
summarise(m = mean(valor))
res %>%
spread(key = param, value = valor) %>%
mutate(range = upr - lwr,
vicio = est - rho,
eqm = (est - rho)^2,
cober = (upr >= rho) & (lwr <= rho)) %>%
group_by(met) %>%
summarise(range = mean(range),
cober = mean(cober),
vicio = mean(vicio),
eqm = mean(eqm))
```
# Exemplo com regressão não linear
```{r}
library(tidyverse)
# Carregando dados.
data(segreg, package = "alr3")
# Visualizando os dados.
ggplot(segreg, aes(x = Temp, y = C)) +
geom_point() +
geom_smooth(se = FALSE, span = 0.4)
# Encontrando valores iniciais.
start <- list(th0 = 75,
th1 = 0.5,
th2 = 50)
ggplot(segreg, aes(x = Temp, y = C)) +
geom_point() +
geom_smooth(se = FALSE, span = 0.4)
plot(C ~ Temp, data = segreg)
with(start, {
curve(th0 + th1 * (x - th2) * (x >= th2) + 0 * (x < th2), add = TRUE)
})
# Ajuste do modelo aos dados.
n0 <- nls(C ~ th0 + th1 * (Temp - th2) * (Temp >= th2) +
0 * (Temp < th2),
data = segreg,
start = start)
summary(n0)
confint.default(n0)
# Curva ajustada sobreposta aos dados.
plot(C ~ Temp, data = segreg)
with(as.list(coef(n0)), {
curve(th0 + th1 * (x - th2) * (x >= th2) + 0 * (x < th2), add = TRUE)
})
#-----------------------------------------------------------------------
# Aplicando o método Jackknife.
# Pseudo valores.
n <- nrow(segreg)
i <- 1:n
pv <- sapply(i,
FUN = function(ii) {
# Reajusta o modelo.
n1 <- update(n0,
data = segreg[-ii, ],
start = coef(n0))
# Obtém os pseudo valores.
n * coef(n0) - (n - 1) * coef(n1)
})
pv
# Estimativas pontuais de Jackknife e erros padrões.
jkest <- cbind("JK_Estimate" = apply(pv, MARGIN = 1, FUN = mean),
"JK_StdError" = apply(pv, MARGIN = 1, FUN = sd)/sqrt(n))
jkest
# Obtendo o IC.
me <- outer(X = c(lwr = -1, upr = 1) *
qt(0.975, df = n - length(coef(n0))),
Y = jkest[, 2],
FUN = "*")
t(sweep(me, MARGIN = 2, STATS = jkest[, 1], FUN = "+"))
#-----------------------------------------------------------------------
# Usando a função nlstools::nlsJack().
library(nlstools)
jk <- nlsJack(n0)
summary(jk)
ggplot(segreg, aes(x = Temp, y = C)) +
# geom_point() +
geom_text(label = 1:nrow(segreg))
```
# Exemplo com regressão polinomial
```{r}
# Carrega os dados da web.
url <- "https://web.stanford.edu/~hastie/CASI_files/DATA/kidney.txt"
kidney <- read.table(url, header = TRUE, sep = " ")
str(kidney)
# Visualiza os dados.
plot(tot ~ age, data = kidney)
# Ajuste do modelo polinomial.
m0 <- lm(tot ~ poly(age, degree = 5), data = kidney)
summary(m0)
# Checando os pressupostos.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# Predição com banda de confiança.
pred <- with(kidney,
data.frame(age = seq(min(age), max(age), by = 1)))
fit0 <- predict(m0, newdata = pred, interval = "confidence")
pred <- cbind(pred, fit0)
# Gráfico da curva ajustada com a banda de confiança.
plot(tot ~ age, data = kidney)
with(pred, matlines(age,
cbind(fit, lwr, upr),
lty = c(1, 2, 2),
col = 1))
# Fazendo Jackknife.
n <- nrow(kidney)
i <- 1:n
# Grid de valores para obter os preditos.
grid <- data.frame(age = seq(20, 85, by = 5))
grid$fit0 <- predict(m0, newdata = grid)
# Obtém os pseudo valores para cada ponto predito.
res <- lapply(i,
FUN = function(ii) {
m1 <- update(m0, data = kidney[-ii, ])
fit1 <- predict(m1, newdata = grid)
pv <- n * grid$fit0 - (n - 1) * fit1
return(pv)
})
# Desmonta lista para matriz.
res <- do.call(cbind, res)
# Obtém estimativa pontual e erro-padrão Jackknife.
res <- apply(res,
MARGIN = 1,
FUN = function(x) {
x <- na.omit(x)
m <- mean(x)
s <- sd(x)/sqrt(n)
c(m = m, s = s)
})
str(res)
# Colocando todos os resultados juntos.
plot(tot ~ age, data = kidney)
with(pred, matlines(age,
cbind(fit, lwr, upr),
lty = c(1, 2, 2),
col = 1))
qtl <- qt(0.975, df = df.residual(m0))
for (i in 1:nrow(grid)) {
points(grid$age[i], res["m", i], pch = 4, col = "red")
arrows(grid$age[i],
res["m", i] - qtl * res["s", i],
grid$age[i],
res["m", i] + qtl * res["s", i],
code = 3,
angle = 90,
length = 0.05,
col = "red")
}
```
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