diff --git a/_site.yml b/_site.yml
index 3473fc751e39fffeb2261068fe1d54dca1d73b9a..c4f899c126193faa17f055d8200ba7cf9a403f7c 100644
--- a/_site.yml
+++ b/_site.yml
@@ -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)"
diff --git a/tutoriais/08-jackknife.Rmd b/tutoriais/08-jackknife.Rmd
new file mode 100644
index 0000000000000000000000000000000000000000..cc4530bd7cdbf7635663901a752f2757b47a1f01
--- /dev/null
+++ b/tutoriais/08-jackknife.Rmd
@@ -0,0 +1,298 @@
+---
+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")
+}
+```