Melhora os exemplos.

parent 7d07fb83
......@@ -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
......
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