df <- read.csv("http://www.leg.ufpr.br/~paulojus/data/dadosExt01.csv")
n <- nrow(df)
x <- df$x
y <- df$y
k <- 5 # nó (knot) usado na simulação
plot(x, y, pch = 19, col = "gray40", main = "Dados simulados")
Regression Splines para Suavização
Parte 1 — Spline Linear com Um Nó: Modelo, Restrição de Continuidade e Estimação do Nó
1 Introdução
Este material apresenta como construir e ajustar regression splines no caso mais simples possível: uma spline linear com um único nó (knot). O objetivo é mostrar como restrições de suavização (continuidade da função e, em extensões, de suas derivadas) são incorporadas na representação matricial \(Y = X\beta + \varepsilon\) do modelo de regressão linear eindicar como isso se generaliza para graus polinomiais maiores.
Os tópicos abordados aqui são:
- Modelo segmentado sem restrição de continuidade;
- Modelo com restrição de continuidade, obtido por reparametrização (redução de colunas de \(X\));
- Formas alternativas de impor a mesma restrição — via multiplicadores de Lagrange (equações normais aumentadas) e via pseudo-observação ponderada;
- Estimação do nó \(k\) quando ele não é conhecido a priori;
- Incerteza de \(\hat k\) via bootstrap — implementado manualmente e via pacote
boot; - Exercícios propostos.
A extensão para splines quadráticas e cúbicas com 1 nó, a generalização para múltiplos nós, a instabilidade numérica da base de potências truncadas e a base de B-splines serão abordados em materiais posteriores.
1.1 Dados utilizados
2 Modelo sem restrição de continuidade
Considere um nó \(k\) que divide o domínio de \(x\) em duas regiões, e a variável indicadora \(IND_i = I(x_i > k)\). O modelo por partes pode ser escrito em uma única equação:
\[y_i = \beta_1 + \beta_2 x_i + \beta_3 IND_i + \beta_4 (x_i IND_i) + \varepsilon_i\]
Para \(x>k\), o intercepto passa a ser \(\beta_1+\beta_3\) e a inclinação \(\beta_2+\beta_4\). A matriz \(X\) tem 4 colunas: \([1,\,x,\,IND,\,xIND]\), e a solução de mínimos quadrados é a usual, \(\hat\beta = (X'X)^{-1}X'y\). O tamanho do salto (gap) no nó é \(\beta_3+\beta_4 k\).
IND <- as.numeric(x > k)
X <- cbind(Intercepto = 1, x = x, IND = IND, xIND = x * IND)
(beta <- drop(solve(crossprod(X), crossprod(X, y))))Intercepto x IND xIND
0.9312854 1.8989314 11.9870281 -2.6128125 ## Conferindo com lm()
(fit <- lm(y ~ x + IND + I(x * IND)))
Call:
lm(formula = y ~ x + IND + I(x * IND))
Coefficients:
(Intercept) x IND I(x * IND)
0.9313 1.8989 11.9870 -2.6128 ## predição
x_pred <- seq(0, 10, length.out = 101)
IND_pred <- as.numeric(x_pred > k)
X_pred <- model.matrix(~ x_pred + IND_pred + I(x_pred * IND_pred))
y_pred <- X_pred %*% beta
plot(x, y, pch = 19, col = "gray40", main = "Dados simulados")
lines(x_pred, y_pred, col = "steelblue", lwd = 2)
## Tamanho do salto no nó com esses betas
(gap <- beta[3] + beta[4] * k) IND
-1.077034 O gap calculado confirma numericamente que, sem a restrição, o ajuste é descontínuo no nó.
3 Modelo com restrição de continuidade
Continuidade em \(x=k\) exige \(\beta_3+\beta_4 k=0 \Rightarrow \beta_3=-\beta_4 k\). Substituindo de volta:
\[y_i = \beta_1 + \beta_2 x_i + \beta_3 (x_i-k)_+ + \varepsilon_i, \qquad (x-k)_+ = \max(x-k,0)\]
A matriz \(X\) passa a ter apenas 3 colunas: \([1,\,x,\,(x-k)_+]\). Perdeu-se exatamente 1 grau de liberdade pela restrição imposta.
xk_pp <- pmax(x - k, 0) # (x-k)_+
X_r <- cbind(Intercepto = 1, x = x, `(x-k)+` = xk_pp)
(beta_r <- drop(solve(crossprod(X_r), crossprod(X_r, y))))Intercepto x (x-k)+
1.242356 1.712845 -2.556061 ## Conferindo com lm()
(fit_r <- lm(y ~ x + xk_pp))
Call:
lm(formula = y ~ x + xk_pp)
Coefficients:
(Intercept) x xk_pp
1.242 1.713 -2.556 ## predição
X_pred_r <- model.matrix(~ x_pred + I((x_pred - k) * IND_pred))
y_pred_r <- X_pred_r %*% beta_r
plot(x, y, pch = 19, col = "gray40", main = "Dados simulados")
lines(x_pred, y_pred, col = "steelblue", lwd = 2)
lines(x_pred, y_pred_r, col = "tomato", lwd = 2)
legend("topright", legend = c("sem restrição", "com restrição"),
col = c("steelblue", "tomato"), lwd = 2, bty = "n")
3.1 Visualização comparativa
par(mfrow = c(1, 2))
plot(x, y, pch = 19, col = "gray40", main = "Sem restrição")
xs1 <- seq(0, k, length.out = 50)
xs2 <- seq(k, 10, length.out = 50)
lines(xs1, cbind(1, xs1, 0, 0) %*% beta, col = "blue", lwd = 2)
lines(xs2, cbind(1, xs2, 1, xs2) %*% beta, col = "red", lwd = 2)
abline(v = k, lty = 2, col = "gray")
plot(x, y, pch = 19, col = "gray40", main = "Com restrição\n(contínuo)")
xs <- seq(0, 10, length.out = 100)
Xs <- cbind(1, xs, pmax(xs - k, 0))
lines(xs, Xs %*% beta_r, col = "darkgreen", lwd = 2)
abline(v = k, lty = 2, col = "gray")
par(mfrow = c(1, 1))
4 Formas alternativas de impor a restrição
A reparametrização acima (redução de colunas de \(X\)) não é a única forma de impor \(R\beta=r\). Duas alternativas clássicas são apresentadas a seguir. Ambas devem produzir o mesmo \(\hat\beta\) restrito.
4.1 Mínimos quadrados restritos via Lagrange (linhas nas equações normais)
Com o modelo “livre” de 4 colunas, a restrição de continuidade é \(R\beta=r\), com \(R=[0\;\;0\;\;1\;\;k]\) e \(r=0\). Minimizando \((y-X\beta)'(y-X\beta)\) sujeito a \(R\beta=r\) via multiplicador de Lagrange \(\lambda\), chega-se ao sistema aumentado
\[\begin{bmatrix} X'X & R' \\ R & 0\end{bmatrix}\begin{bmatrix}\beta\\\lambda\end{bmatrix} = \begin{bmatrix}X'y\\ r\end{bmatrix}\]
com solução fechada equivalente
\[\hat\beta_R = \hat\beta - (X'X)^{-1}R'\big[R(X'X)^{-1}R'\big]^{-1}(R\hat\beta - r)\]
O código a seguir mostra as duas formas equivalentes: resolver o sistema de equações normais aumentado ou aplicar a correção em \(\hat\beta\).
R <- t(c(0,0,1,k))
r <- 0
XtX <- crossprod(X)
Xty <- crossprod(X, y)
XtXR <- rbind(XtX, R)
XtXr <- cbind(XtXR, c(R, 0))
drop(solve(XtXr, c(Xty, r)))
## Intercepto x IND xIND
## 1.242356 1.712845 12.780304 -2.556061 -2.609633
beta_r
## Intercepto x (x-k)+
## 1.242356 1.712845 -2.556061
beta - drop(solve(XtX, crossprod(R, solve(R %*% solve(XtX, t(R)), (R%*%beta - r)))))
## Intercepto x IND xIND
## 1.242356 1.712845 12.780304 -2.5560614.2 Pseudo-observação ponderada (linha extra em \(X\))
Adiciona-se uma linha extra em \(X\) (e em \(y\)) representando a restrição como uma “observação” adicional, com peso \(w\) muito grande:
\[X_{aug} = \begin{bmatrix} X \\ \sqrt{w}\,R\end{bmatrix}, \qquad y_{aug} = \begin{bmatrix} y \\ \sqrt{w}\,r\end{bmatrix}\]
Quando \(w\to\infty\), essa observação domina o ajuste e força \(R\hat\beta \to r\). Este mecanismo é a base conceitual das P-splines (Eilers & Marx), em que uma penalidade de rugosidade é adicionada da mesma forma, ponderada por \(\sqrt\lambda\).
Aqui também mostramos duar formas de implementar a ponderação
## (pseudo)Observação ponderada
w <- 1e8
Xaug <- rbind(X, sqrt(w)*c(0,0,1,k))
yaug <- c(y, sqrt(w)*0)
drop(solve(crossprod(Xaug), crossprod(Xaug, yaug)))
## Intercepto x IND xIND
## 1.242356 1.712845 12.780304 -2.556061
beta_r
## Intercepto x (x-k)+
## 1.242356 1.712845 -2.556061
## ou usando lm() com argumento weights
Xaug1 <- rbind(X, c(0,0,1,k))
yaug1 <- c(y, 0)
coef(lm(yaug1 ~ Xaug1 - 1, weights = c(rep(1, n), w)))
## Xaug1Intercepto Xaug1x Xaug1IND Xaug1xIND
## 1.242356 1.712845 12.780304 -2.5560615 Estimação do nó \(k\)
Na base restrita \(X=[1,\,x,\,(x-k)_+]\), o parâmetro \(k\) aparece dentro da própria coluna, ou seja, para cada candidato de \(k\), \(X\) muda inteiramente. O modelo é linear em \(\beta\) apenas condicional a \(k\) fixo, o que torna a estimação conjunta um problema de mínimos quadrados não lineares em \(k\).
5.1 Perfil de RSS (busca em grade)
Os passos são:
- Define-se um grid de valores candidatos para \(k\);
- Para cada \(k\) candidato:
- Constroi-se as variáveis do modelo;
- Ajusta-se o modelo;
- Calcula-se a soma de quadrados dos resíduos (RSS(k));
- Toma-se o valor de \(k\) da menor RSS(k).
\[RSS(k) = \big(y - X(k)\hat\beta(k)\big)'\big(y - X(k)\hat\beta(k)\big), \qquad \hat k = \arg\min_k RSS(k)\]
rss_perfil <- function(k_cand, x, y) {
xk_cand <- pmax(x - k_cand, 0)
sum(residuals(lm(y ~ x + xk_cand))^2)
}
## grade de candidatos a k evitando bordas
k_grid <- seq(quantile(x, 0.10), quantile(x, 0.90), length.out = 200)
RSS_k <- sapply(k_grid, rss_perfil, x = x, y = y)
(k_hat <- k_grid[which.min(RSS_k)])[1] 4.597513plot(k_grid, RSS_k, type = "l", lwd = 2, col = "steelblue",
xlab = "k candidato", ylab = "RSS(k)",
main = "Perfil de RSS para k")
abline(v = k_hat, lty = 2, col = "tomato")
points(k_hat, min(RSS_k), pch = 19, col = "tomato")
legend("topleft", legend = sprintf("k_hat = %.3f", k_hat),
bty = "n", text.col = "tomato")
xk_pp_hat <- pmax(x - k_hat, 0)
X_r_hat <- cbind(Intercepto = 1, x = x, `(x-k)+` = xk_pp_hat)
(beta_r_hat <- drop(solve(crossprod(X_r_hat), crossprod(X_r_hat, y))))Intercepto x (x-k)+
0.9121986 1.9109675 -2.6407503 x_pred_hat <- seq(0, 10, length.out = 101)
X_pred_r_hat <- cbind(1, x_pred_hat, pmax(x_pred_hat - k_hat, 0))
y_pred_r_hat <- X_pred_r_hat %*% beta_r_hat
plot(x, y, pch = 19, col = "gray40", main = "Nó estimado vs nó fixo (k=5)")
lines(x_pred, y_pred_r, col = "tomato", lwd = 2, lty = 2) # k fixo = 5
lines(x_pred_hat, y_pred_r_hat, col = "darkgreen", lwd = 2) # k estimado
abline(v = k, lty = 3, col = "tomato")
abline(v = k_hat, lty = 3, col = "darkgreen")
legend("topleft", legend = c("k fixo = 5", sprintf("k_hat = %.3f", k_hat)),
col = c("tomato", "darkgreen"), lwd = 2, lty = c(2,1), bty = "n")
5.2 Conferência com o pacote segmented
O pacote segmented (Muggeo, 2003) expande \((x-k)_+\) em série de Taylor em torno de um valor inicial \(k_0\):
\[(x-k)_+ \approx (x-k_0)_+ - (k-k_0)\,I(x>k_0)\]
o que permite estimar a correção \(\Delta k\) via um modelo linear auxiliar, iterando até convergência o que corresponde a um Gauss-Newton “disfarçado” de regressão linear.
fit_lin <- lm(y ~ x)
fit_seg <- segmented::segmented(fit_lin, seg.Z = ~x, psi = list(x = 5))
summary(fit_seg)
***Regression Model with Segmented Relationship(s)***
Call:
segmented.lm(obj = fit_lin, seg.Z = ~x, psi = list(x = 5))
Estimated Break-Point(s):
Est. St.Err
psi1.x 4.597 0.202
Coefficients of the linear terms:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.9120 0.4270 2.136 0.0396 *
x 1.9111 0.1539 12.414 1.43e-14 ***
U1.x -2.6408 0.1863 -14.176 NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.83 on 36 degrees of freedom
Multiple R-Squared: 0.8776, Adjusted R-squared: 0.8674
Boot restarting based on 6 samples. Last fit:
Convergence attained in 2 iterations (rel. change 1.0356e-11)6 Incerteza de \(\hat k\): bootstrap
O intervalo de confiança de \(\hat k\) não sai da fórmula usual de mínimos quadrados (que assume \(X\) fixo), já que \(X\) depende do próprio \(\hat k\) estimado.
O bootstrap segue os seguintes passos:
6.1 Bootstrap manual (não paramétrico, pares)
set.seed(202602)
B <- 1000
k_boot <- numeric(B)
## mesma grade usada no perfil original para as réplicas
k_grid <- seq(quantile(x, 0.10), quantile(x, 0.90), length.out = 200)
for (b in 1:B) {
idx_b <- sample(seq_along(x), size = n, replace = TRUE)
x_b <- x[idx_b]
y_b <- y[idx_b]
RSS_b <- sapply(k_grid, rss_perfil, x = x_b, y = y_b)
k_boot[b] <- k_grid[which.min(RSS_b)]
}
summary(k_boot) Min. 1st Qu. Median Mean 3rd Qu. Max.
3.875 4.477 4.598 4.590 4.718 5.280 (erro_padrao_k <- sd(k_boot))[1] 0.1695747(IC_normal <- k_hat + c(-1, 1) * qnorm(0.975) * erro_padrao_k)[1] 4.265152 4.929873(IC_percentil <- quantile(k_boot, probs = c(0.025, 0.975))) 2.5% 97.5%
4.275343 4.918678 par(mfrow = c(1, 2))
hist(k_boot, breaks = 30, col = "lightblue", border = "white",
main = "Distribuição bootstrap de k_hat", xlab = "k_hat*")
abline(v = k_hat, col = "darkgreen", lwd = 2)
abline(v = IC_percentil, col = "tomato", lwd = 2, lty = 2)
legend("topright", legend = c("k_hat (dados originais)", "IC percentil 95%"),
col = c("darkgreen", "tomato"), lwd = 2, lty = c(1,2), bty = "n", cex = 0.8)
boxplot(k_boot, horizontal = TRUE, col = "lightblue",
main = "Boxplot de k_hat*", xlab = "k_hat*")
abline(v = k, col = "gray40", lty = 3)
points(k_hat, 1, pch = 19, col = "darkgreen")
legend("topright", legend = "k verdadeiro", lty = 3, col = "gray40", bty = "n")
par(mfrow = c(1, 1))
Nota
Melhoria futura: a grade k_grid usada dentro do loop bootstrap é fixa, calculada sobre os dados originais. Tecnicamente seria mais correto recalculá-la dinamicamente em cada réplica (por exemplo, quantile(x_b, c(.10, .90))), já que o range de x_b varia entre reamostras. Fica registrado como possível aprimoramento.
6.2 Bootstrap via pacote boot (com IC BCa)
k_stat <- function(dados, indices) {
x_b <- dados$x[indices]
y_b <- dados$y[indices]
RSS_b <- sapply(k_grid, rss_perfil, x = x_b, y = y_b)
k_grid[which.min(RSS_b)]
}
set.seed(202602) # mesma seed do bootstrap manual
boot_out <- boot::boot(data = df, statistic = k_stat, R = 1000)
boot_out
ORDINARY NONPARAMETRIC BOOTSTRAP
Call:
boot::boot(data = df, statistic = k_stat, R = 1000)
Bootstrap Statistics :
original bias std. error
t1* 4.597513 -0.007105794 0.1910292ci_out <- boot::boot.ci(boot_out, type = c("norm", "basic", "perc", "bca"))
ci_outBOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 1000 bootstrap replicates
CALL :
boot::boot.ci(boot.out = boot_out, type = c("norm", "basic",
"perc", "bca"))
Intervals :
Level Normal Basic
95% ( 4.230, 4.979 ) ( 4.236, 4.959 )
Level Percentile BCa
95% ( 4.236, 4.959 ) ( 4.196, 4.919 )
Calculations and Intervals on Original Scaleplot(boot_out)
6.3 Comparando visualmente os diferentes tipos de intervalo de confiança
O objeto retornado por boot.ci() contém os limites de cada tipo de intervalo em posições diferentes da matriz correspondente. Extraímos os quatro tipos calculados acima (Normal, Básico, Percentil e BCa) e os plotamos lado a lado, junto com \(\hat k\) e o valor verdadeiro \(k=5\) usado na simulação dos dados:
## extraindo os limites inferior/superior de cada tipo de IC
lim_norm <- ci_out$normal[2:3]
lim_basic <- ci_out$basic[4:5]
lim_perc <- ci_out$percent[4:5]
lim_bca <- ci_out$bca[4:5]
ci_tab <- rbind(
Normal = lim_norm,
Basico = lim_basic,
Percentil = lim_perc,
BCa = lim_bca
)
colnames(ci_tab) <- c("liminf", "limsup")
knitr::kable(round(ci_tab, 4), caption = "Limites dos IC 95% para k, por tipo")| liminf | limsup | |
|---|---|---|
| Normal | 4.2302 | 4.9790 |
| Basico | 4.2362 | 4.9588 |
| Percentil | 4.2362 | 4.9588 |
| BCa | 4.1961 | 4.9187 |
## grafico comparativo
plot(NA, xlim = range(c(ci_tab, k, k_hat)), ylim = c(0.5, nrow(ci_tab) + 0.5),
yaxt = "n", xlab = "k", ylab = "",
main = "IC bootstrap (95%) para k, por tipo")
axis(2, at = 1:nrow(ci_tab), labels = rownames(ci_tab), las = 1)
for (i in 1:nrow(ci_tab)) {
segments(ci_tab[i, 1], i, ci_tab[i, 2], i, lwd = 4, col = "steelblue")
points(k_hat, i, pch = 19, col = "darkgreen")
}
abline(v = k, col = "tomato", lty = 2, lwd = 2)
legend("bottomright",
legend = c("k_hat (estimativa pontual)", "k verdadeiro = 5"),
col = c("darkgreen", "tomato"),
pch = c(19, NA), lty = c(NA, 2), lwd = c(NA, 2), bty = "n")
Leitura do gráfico: todos os quatro intervalos devem cobrir o valor verdadeiro \(k=5\) (linha tracejada vermelha) neste exemplo. Vale observar se o IC BCa fica deslocado ou assimétrico em relação aos demais — isso é esperado quando a distribuição bootstrap de \(\hat k\) (vista no histograma anterior) não é simétrica, e é exatamente a situação em que o BCa tende a ser mais confiável que o IC normal ou percentil simples.
7 Síntese
| Etapa | Ferramenta “na unha” | Pacote/extensão equivalente |
|---|---|---|
| Ajuste com restrição de continuidade | crossprod + solve na base \((x-k)_+\) |
lm(y ~ x + xk_pp) |
| Restrição alternativa (Lagrange / pseudo-obs.) | equações normais aumentadas / linha extra em \(X\) | base conceitual das P-splines |
| Estimação de \(k\) | perfil de RSS por busca em grade | segmented::segmented() |
| Incerteza de \(\hat k\) | bootstrap de pares manual | boot::boot() + boot.ci(type="bca") |
Esta primeira parte cobriu o caso de um único nó, com spline linear. Materiais posteriores tratam a extensão para splines quadráticas e cúbicas, a generalização para múltiplos nós, a instabilidade numérica da base de potências truncadas, e a transição para B-splines.
8 Exercícios propostos
Exercício 1 — Formas alternativas de restrição, outro nó. Repita a comparação entre reparametrização, Lagrange e pseudo-observação ponderada (Seção “Formas alternativas de impor a restrição”), agora fixando \(k=3\) em vez de \(k=5\). Construa a matriz \(R\) correspondente, obtenha \(\hat\beta_R\) pelas três formas e monte uma tabela comparativa como a apresentada no texto. Os três resultados continuam coincidindo? O que muda na matriz \(R\) quando o nó muda?
Exercício 2 — Sensibilidade da estimação de \(k\) ao ruído. Gere um novo conjunto de dados simulado, com a mesma função geradora usada em dadosExt01.csv (spline linear com quebra em \(k=5\)), mas aumentando o desvio-padrão do ruído (por exemplo, sd = 3 em vez de sd = 1). Refaça o perfil de \(RSS(k)\) e a estimativa \(\hat k\) para esses novos dados. O vale do perfil de RSS fica mais raso ou mais estreito? Refaça também o bootstrap manual (Seção “Bootstrap manual”) e compare o erro-padrão de \(\hat k\) com o obtido para os dados originais — a incerteza aumentou como esperado?
Exercício 3 — Regressão polinomial global, para comparação. Usando os mesmos dados (df), ajuste modelos de regressão polinomial global (sem nó, sem restrição de continuidade por partes) de graus 1, 2, 3 e 4 — por exemplo, lm(y ~ poly(x, grau, raw = TRUE)). Sobreponha as quatro curvas ajustadas e a curva da spline linear restrita (com \(k=5\)) em um único gráfico. Compare visualmente o comportamento das curvas, especialmente nas extremidades do domínio de \(x\). Calcule e compare o \(R^2\) e o RSS de cada ajuste — o polinômio de grau mais alto necessariamente produz um ajuste “melhor”? Que problema conhecido de polinômios globais de grau alto (ex.: oscilações espúrias, fenômeno de Runge) fica evidente nessa comparação, e por que a spline (local, por partes) tende a evitá-lo?
Exercício 4 — Modelo linear com platô (linear-plateau), nó fixo e estimado. Um modelo alternativo, comum em aplicações (ex.: resposta a dose, crescimento), é o linear com platô: a resposta cresce linearmente até um ponto de corte \(c\) e depois permanece constante (inclinação zero) — ou seja, \(y_i=\beta_1+\beta_2\min(x_i,c)+\varepsilon_i\), com matriz \(X=[1,\,\min(x,c)]\) (apenas 2 colunas). Note a diferença em relação à spline linear com restrição de continuidade vista neste material: lá a inclinação após o nó era livre (\(\beta_2+\beta_3\)); aqui ela é forçada a ser exatamente zero após \(c\).
Ajuste esse modelo aos dados (
df) com o ponto de corte fixo em \(c=5\), construindo \(X\) explicitamente e obtendo \(\hat\beta\) por mínimos quadrados (confira comlm()). Compare visualmente o ajuste com o da spline linear restrita já vista (mesmo \(k=5\)).Agora trate \(c\) como desconhecido e estime-o pela mesma estratégia de perfil de RSS usada para \(\hat k\) (Seção “Estimação do nó \(k\)”), testando uma grade de candidatos \(c\) e escolhendo o que minimiza \(RSS(c)=\sum_i(y_i-\hat\beta_1-\hat\beta_2\min(x_i,c))^2\). Compare \(\hat c\) com o \(\hat k\) obtido para a spline linear restrita — fazem sentido próximos, iguais, ou diferentes, dado que os modelos representam formas distintas depois do ponto de corte?
(Opcional) Repita o bootstrap manual de pares para obter o erro-padrão e um IC para \(\hat c\), e compare com o obtido para \(\hat k\).
Exercício 5 (opcional) — Estabilidade do IC BCa. Repita o bootstrap via pacote boot (Seção “Bootstrap via pacote boot”) com pelo menos 3 sementes (set.seed) diferentes e \(R=1000\) réplicas cada vez. Monte um gráfico como o da Seção “Comparando visualmente os diferentes tipos de intervalo de confiança”, mas agora comparando o IC BCa obtido em cada semente. O intervalo permanece razoavelmente estável entre execuções? O que isso sugere sobre o número de réplicas \(B\) necessário para um resultado confiável em relatórios ou publicações?