diff --git a/slides/08-jackknife.Rmd b/slides/08-jackknife.Rmd
index a4f4501b86e9687e059885b9fa7228813fae7300..980eee7b0986d096ba2dc6f9286916703ffcb28f 100644
--- a/slides/08-jackknife.Rmd
+++ b/slides/08-jackknife.Rmd
@@ -15,12 +15,12 @@ csl: ../config/ABNT-UFPR-2011-Mendeley.csl
## Ideia
-A ideia fundamentada no estimador da média
+A ideia é fundamentada no estimador da média
$$
\bar{X} = \frac{1}{n} \sum_{i=1}^{n} X_i.
$$
-A média com a $j$-ésima observação removida, $X_j$ é
+A média com a $j$-ésima observação removida, $X_j$, é
$$
\bar{X}_{-j} = \frac{1}{n - 1}
\left( \sum_{i=1}^{n} X_i \right) - X_j.
@@ -123,82 +123,21 @@ De acordo com @efron2016computerage
* O erro-padrão Jackknife é viciado para estimar o verdadeiro erro
padrão.
-## Exemplo 2 - regressão local
-
-```{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)
-
-n <- nrow(kidney)
-i <- 1:n
-# kidney$i <- i
-
-# Ajuste do modelo de regressão local.
-m0 <- with(kidney, lowess(x = age, y = tot, f = 1/3))
-fit0 <- unique(as.data.frame(m0))
-head(fit0)
-
-plot(tot ~ age, data = kidney)
-lines(m0)
-
-# Obtendo os erros padrões de predição Jackknife.
-res <- lapply(i,
- FUN = function(ii) {
- m1 <- with(kidney[-ii, ],
- lowess(x = age,
- y = tot,
- f = 1/3))
- fit1 <- unique(as.data.frame(m1))
- names(fit1)[2] <- "y1"
- fit <- merge(fit0, fit1, by = "x")
- fit$pv <- with(fit, n * y - (n - 1) * y1)
- return(fit[, c("x", "pv")])
- })
-res <- do.call(rbind, res)
-str(res)
-head(res)
-
-res <- aggregate(pv ~ x,
- data = res,
- FUN = function(y) {
- c(m = mean(y),
- sd = sd(y)/sqrt(n))
- })
-
-str(res)
-head(res)
-
-# ATTENTION: requer revisão.
-
-plot(tot ~ age, data = kidney)
-lines(m0)
-abline(v = 25, lty = 2, col = "gray")
-for (i in 1:nrow(res)) {
- segments(res$x[i],
- res$pv[i, 1] - 2 * res$pv[i, 2],
- res$x[i],
- res$pv[i, 1] + 2 * res$pv[i, 2])
-}
-
-res[res$x == 25, ]
-
-# TODO: continuar.
-```
-
-## Exemplo 3 - regressão não linear
+## Exemplo 2 - regressão não linear
```{r}
+# Carregando dados.
data(segreg, package = "alr3")
library(latticeExtra)
-xyplot(C ~ Temp, data = segreg, type = c("p", "smooth"))
+# Visualizando os dados.
+xyplot(C ~ Temp,
+ data = segreg,
+ type = c("p", "smooth"),
+ span = 0.4)
+# Encontrando valores iniciais.
start <- list(th0 = 75,
th1 = 0.5,
th2 = 50)
@@ -207,6 +146,7 @@ xyplot(C ~ Temp, data = segreg) +
0 * (x < th2)),
data = start)
+# Ajuste do modelo aos dados.
n0 <- nls(C ~ th0 + th1 * (Temp - th2) * (Temp >= th2) +
0 * (Temp < th2),
data = segreg,
@@ -214,8 +154,131 @@ n0 <- nls(C ~ th0 + th1 * (Temp - th2) * (Temp >= th2) +
summary(n0)
confint.default(n0)
-# TODO:
-# Usar Jackknife para obter erro padrão para os parâmetros.
+# Curva ajustada sobreposta aos dados.
+xyplot(C ~ Temp, data = segreg) +
+ layer(panel.curve(th0 + th1 * (x - th2) * (x >= th2) +
+ 0 * (x < th2),
+ col = "orange"),
+ data = as.list(coef(n0)))
+
+#-----------------------------------------------------------------------
+# 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)
+```
+
+## Exemplo 3 - 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")
+}
```
## Próxima aula