CE-092: Extensões de modelos de regressão

Tópicos Avançados de Splines: Estabilidade Numérica, Bases bs()/ns(), Intervalos e Extrapolação

Introdução

Os materiais anteriores da série regressao-polinomial-eolico a regressao-partes23-eolico, e o aprofundamento em regressao-partes-texto construíram uma sequencia completa: médias, retas, parábolas e cúbicas por partes, utilizando o conjunto de dados eólicos e com dado simulado de nó conhecido que foi utilizado para validar estimação e bootstrap. Outros tópicos ficaram visíveis com dados que sugeriam uma curvatura mais pronunciada e extrapolar para fora do domínio observado, ilustrado nos scripts medias2splines.R e xySpline.R exploram.

Este material explora o conteúdo desses dois scripts e trata:

1. Recapitulação rápida da escada de restrições, agora com uma curva que tem curvatura real. Painéis mostrando o comportamento de funções ajustadas e suas derivadas ilustram os efeitos das restrições;
2. Instabilidade numérica da base de potências truncadas, quantificada pelo número de condição de \(X'X\);
3. Uso prático de splines::bs(): nós explícitos vs. automáticos (df=), e o papel de Boundary.knots;
4. Splines cúbicos naturais (ns()): a restrição extra que impõem, o custo em ajuste, e por que compensa fora da amostra;
5. As funções base de splines::bs() e splines::ns() são comentadas;
6. Intervalos de confiança vs. de predição;
7. Extrapolação: onde todas as diferenças entre os modelos aparecem com clareza.

Dados

Dados simulados com curvatura real, bem diferente da curva eólica, quase linear, gerados com \(n=50\), \(x\sim U(0,10)\) e \(Y_i = 0{,}2x_i + \cos(x_i+1) + \varepsilon_i\), \(\varepsilon_i\sim N(0,0{,}5^2)\), a mesma função geradora utilizada em medias2splines.R.

library(splines)  ## bs(), ns() -- carregado aqui (não em um chunk no meio do
                   ## documento) porque o efeito colateral de attach do pacote
                   ## não é reproduzido de forma confiável por chunks em cache

set.seed(201802)
df <- data.frame(x = sort(runif(50, 0, 10)))
df <- transform(df, Ytrue = 0.2 * x + cos(x + 1))
df <- transform(df, Y = Ytrue + rnorm(50, mean = 0, sd = 0.5))

nd <- data.frame(x = seq(0, 10, length = 2001))
nd <- transform(nd, Ytrue = 0.2 * x + cos(x + 1))
## NA exatamente nos nós, para não conectar segmentos separados nos gráficos abaixo
nd$x[nd$x == 3.5] <- NA
nd$x[nd$x == 7]   <- NA

with(df, plot(Y ~ x, pch = 20, col = "gray30"))
with(nd, lines(Ytrue ~ x, lty = 3))

Dados simulados (curva verdadeira tracejada)

fitFun <- function(m) c(logLik = as.numeric(logLik(m)), df = attr(logLik(m), "df"),
                         AIC = AIC(m), BIC = BIC(m))

Recapitulando a escada de restrições, com uma curva de verdade

Com dois nós fixos em \(x=3{,}5\) e \(x=7\) (três regiões), a mesma escada de 0 a 3 restrições de regressao-partes23-eolico se aplica aqui. Não repetimos a dedução algébrica que já feita naquele material e nos concentramos no ajuste e a comparação, para servir de base à visualização de derivadas a seguir.

df <- transform(df, x1 = ifelse(x < 3.5, 0, x - 3.5), x2 = ifelse(x < 7, 0, x - 7))
nd  <- transform(nd,  x1 = ifelse(x < 3.5, 0, x - 3.5), x2 = ifelse(x < 7, 0, x - 7))
df <- transform(df, I1 = ifelse(x < 3.5, 0, 1), I2 = ifelse(x < 7, 0, 1))
nd  <- transform(nd,  I1 = ifelse(x < 3.5, 0, 1), I2 = ifelse(x < 7, 0, 1))

fit1 <- lm(Y ~ x, data = df)                                   # reta global
fit3 <- lm(Y ~ poly(x, 3), data = df)                           # cúbica global
fit4 <- lm(Y ~ poly(x, 4), data = df)                           # quártica global

## 0 restrições ("livre"): jump permitido, cúbica independente em cada região
m3_livre <- lm(Y ~ poly(x, 3) + I1 + poly(x1, 3) + I2 + poly(x2, 3), data = df)
## 1 restrição (continuidade da função)
m3_r1 <- lm(Y ~ poly(x, 3) + poly(x1, 3) + poly(x2, 3), data = df)
## 2 restrições (continuidade da função e da 1a derivada)
m3_r2 <- lm(Y ~ poly(x, 3, raw = TRUE) + poly(x1, 3, raw = TRUE)[, -1] +
                poly(x2, 3, raw = TRUE)[, -1], data = df)
## 3 restrições (spline cúbico "de fato")
m3_r3 <- lm(Y ~ poly(x, 3) + I(x1^3) + I(x2^3), data = df)

nd <- transform(nd, y1 = predict(fit1, newdata = nd), y3 = predict(fit3, newdata = nd),
                 y4 = predict(fit4, newdata = nd))
nd <- transform(nd, y3.livre = predict(m3_livre, newdata = nd),
                 y3.r1 = predict(m3_r1, newdata = nd),
                 y3.r2 = predict(m3_r2, newdata = nd),
                 y3.r3 = predict(m3_r3, newdata = nd))

sapply(list(livre = m3_livre, r1 = m3_r1, r2 = m3_r2, r3 = m3_r3, cubica = fit3, quartica = fit4), fitFun)

##            livre        r1        r2        r3    cubica  quartica
## logLik -35.41825 -37.38099 -38.69683 -40.06912 -66.56633 -39.59515
## df      13.00000  11.00000   9.00000   7.00000   5.00000   6.00000
## AIC     96.83650  96.76198  95.39365  94.13824 143.13266  91.19029
## BIC    121.69280 117.79424 112.60186 107.52240 152.69278 102.66243

Função e suas derivadas, para cada nível de restrição

O painel a seguir é o coração da comparação: cada coluna mostra a curva ajustada (topo), sua 1ª derivada numérica (meio) e sua 2ª derivada numérica (base). As quebras que desaparecem de uma coluna para a seguinte são exatamente as restrições impostas — e ficam muito mais visíveis aqui do que na curva eólica (suave demais para revelar quinas de derivada), porque a curva geradora (\(\cos\)) tem curvatura de verdade.

dxx <- with(nd, unique(diff(x)))[1]

deriv1 <- function(y) c(diff(y) / dxx, NA)
deriv2 <- function(y1) c(diff(y1) / dxx, NA)

nd <- transform(nd, y3.livre.d1 = deriv1(y3.livre)); nd <- transform(nd, y3.livre.d2 = deriv2(y3.livre.d1))
nd <- transform(nd, y3.r1.d1 = deriv1(y3.r1));       nd <- transform(nd, y3.r1.d2 = deriv2(y3.r1.d1))
nd <- transform(nd, y3.r2.d1 = deriv1(y3.r2));       nd <- transform(nd, y3.r2.d2 = deriv2(y3.r2.d1))
nd <- transform(nd, y3.r3.d1 = deriv1(y3.r3));       nd <- transform(nd, y3.r3.d2 = deriv2(y3.r3.d1))
nd <- transform(nd, y3.d1 = deriv1(y3));             nd <- transform(nd, y3.d2 = deriv2(y3.d1))
nd <- transform(nd, y4.d1 = deriv1(y4));             nd <- transform(nd, y4.d2 = deriv2(y4.d1))

par(mfcol = c(3, 6), mar = c(3, 3, 1.5, 0.5), mgp = c(1.7, 0.7, 0))

for (nome in c("livre", "r1", "r2", "r3")) {
    with(df, plot(Y ~ x, pch = 20, cex = 0.6, col = "gray40", main = nome))
    with(nd, lines(get(paste0("y3.", nome)) ~ x))
    plot(get(paste0("y3.", nome, ".d1")) ~ x, data = nd, type = "l", ylab = "1a deriv.")
    plot(get(paste0("y3.", nome, ".d2")) ~ x, data = nd, type = "l", ylab = "2a deriv.")
}
with(df, plot(Y ~ x, pch = 20, cex = 0.6, col = "gray40", main = "cúbica global"))
with(nd, lines(y3 ~ x)); plot(y3.d1 ~ x, data = nd, type = "l", ylab = "1a deriv.")
plot(y3.d2 ~ x, data = nd, type = "l", ylab = "2a deriv.")
with(df, plot(Y ~ x, pch = 20, cex = 0.6, col = "gray40", main = "quártica global"))
with(nd, lines(y4 ~ x)); plot(y4.d1 ~ x, data = nd, type = "l", ylab = "1a deriv.")
plot(y4.d2 ~ x, data = nd, type = "l", ylab = "2a deriv.")

Função, 1ª e 2ª derivadas (numéricas): livre, 1, 2, 3 restrições, cúbica e quártica globais

par(mfrow = c(1, 1))

Note como a 1ª derivada só fica sem saltos a partir de “r2” (2 restrições) e a 2ª derivada só fica sem saltos em “r3” — e como as duas últimas colunas (polinômios globais) não têm quinas em lugar nenhum, mas em compensação não conseguem acompanhar a curvatura local dos dados tão bem quanto a spline cúbica (compare os valores de logLik na tabela acima: a cúbica global tem log-verossimilhança pior que a spline “r3” com o mesmo número de parâmetros).

Instabilidade numérica da base de potências truncadas

Os quatro modelos acima foram construídos com potências de \(x\), \(x_1=(x-3{,}5)_+\) e \(x_2=(x-7)_+\) — a base de potências truncadas. Ela é a mais fácil de escrever à mão, mas tem um problema conhecido: colunas como \(x, x^2, x^3\) (e suas versões truncadas) ficam fortemente correlacionadas entre si, o que deixa \(X'X\) mal condicionada. Isso não afeta os valores ajustados (o modelo continua sendo exatamente o mesmo, matematicamente), mas afeta a precisão numérica do cálculo de \(\hat\beta=(X'X)^{-1}X'y\).

O número de condição de \(X'X\) mede esse mal-condicionamento — quanto maior, mais instável. Comparamos, para o mesmo modelo (mesmo nível de restrição, portanto mesmo espaço de colunas), duas parametrizações:

## Mesmo modelo (1 restrição), duas bases: potências cruas vs. poly() ortogonal
X_1r_raw  <- model.matrix(Y ~ x + I(x^2) + I(x^3) +
                               I(x1) + I(x1^2) + I(x1^3) +
                               I(x2) + I(x2^2) + I(x2^3), data = df)
X_1r_orth <- model.matrix(Y ~ poly(x, 3) + poly(x1, 3) + poly(x2, 3), data = df)

## Mesmo modelo (3 restrições = spline cúbico "de fato"), duas bases: crua vs. bs()
X_3r_raw <- model.matrix(Y ~ poly(x, 3, raw = TRUE) + I(x1^3) + I(x2^3), data = df)
X_3r_bs  <- model.matrix(Y ~ bs(x, knots = c(3.5, 7)), data = df)

data.frame(
    modelo = c("1 restrição, base crua", "1 restrição, poly() ortogonal",
               "3 restrições, base truncada crua", "3 restrições, bs()"),
    numero_condicao = c(kappa(crossprod(X_1r_raw)), kappa(crossprod(X_1r_orth)),
                         kappa(crossprod(X_3r_raw)), kappa(crossprod(X_3r_bs)))
)

##                             modelo numero_condicao
## 1           1 restrição, base crua    1.586891e+08
## 2    1 restrição, poly() ortogonal    1.537788e+07
## 3 3 restrições, base truncada crua    4.835958e+07
## 4               3 restrições, bs()    4.626399e+02

A diferença é enorme: para o modelo de 3 restrições (o spline cúbico final), o número de condição cai de ~48 milhões (base crua) para ~463 (bs()) — mais de 5 ordens de grandeza. E, apesar disso, os dois ajustes são idênticos:

m_raw <- lm(Y ~ poly(x, 3, raw = TRUE) + I(x1^3) + I(x2^3), data = df)
m_bs  <- lm(Y ~ bs(x, knots = c(3.5, 7)), data = df)

c(logLik_raw = as.numeric(logLik(m_raw)), logLik_bs = as.numeric(logLik(m_bs)))

## logLik_raw  logLik_bs 
##  -40.06912  -40.06912

max(abs(fitted(m_raw) - fitted(m_bs)))

## [1] 2.9865e-14

Ou seja: a instabilidade numérica é um problema de como calcular \(\hat\beta\), não de o que o modelo representa. poly() (base ortogonal) já ajuda bastante (comparado à base crua), mas é bs() — a base B-spline do pacote splines — que resolve o problema de forma robusta, e é por isso que ela é a escolha padrão na prática a partir de poucos nós.

bs(): nós explícitos, nós automáticos e Boundary.knots

fit_bs_k <- lm(Y ~ bs(x, knots = c(3.5, 7), Boundary.knots = c(0, 10)), data = df)

Em vez de escolher os nós manualmente, bs(x, df = ...) os posiciona automaticamente nos quantis de \(x\):

fit_bs_df5 <- lm(Y ~ bs(x, df = 5, Boundary.knots = c(0, 10)), data = df)
attr(with(df, bs(x, df = 5, Boundary.knots = c(0, 10))), "knots")

## [1] 3.544345 6.239567

Com df=5 (e grau 3), bs() escolheu automaticamente 2 nós internos, nos quantis de x — próximos, mas não iguais, aos \(3{,}5\) e \(7\) escolhidos “à mão”.

Boundary.knots merece atenção: se omitido, bs() usa o intervalo amostral de \(x\) (range(df$x) \(\approx\) 0.28, 9.95), não o domínio pretendido do problema (\([0,10]\), por construção da simulação). Definir Boundary.knots explicitamente como o domínio pretendido evita que a base mude de forma sutil dependendo de onde a amostra por acaso começou e terminou — e, como veremos na Seção “Extrapolação”, também importa para o que acontece fora da amostra.

ns(): splines cúbicos naturais

ns() (natural splines) adiciona uma restrição extra às já vistas: força a 2ª e a 3ª derivadas a serem zero além dos nós de fronteira — ou seja, a curva é exatamente linear fora do intervalo definido por Boundary.knots. Isso custa 2 graus de liberdade a mais que bs() para o mesmo conjunto de nós.

fit_ns_k <- lm(Y ~ ns(x, knots = c(3.5, 7), Boundary.knots = c(0, 10)), data = df)
fit_ns_df5 <- lm(Y ~ ns(x, df = 5, Boundary.knots = c(0, 10)), data = df)
attr(with(df, ns(x, df = 5, Boundary.knots = c(0, 10))), "knots")

## [1] 1.802739 3.857011 5.184644 8.104730

data.frame(
    modelo = c("bs 2 nós", "bs df=5", "ns 2 nós", "ns df=5", "cúbica global"),
    gl     = sapply(list(fit_bs_k, fit_bs_df5, fit_ns_k, fit_ns_df5, fit3), function(m) length(coef(m))),
    logLik = round(sapply(list(fit_bs_k, fit_bs_df5, fit_ns_k, fit_ns_df5, fit3), function(m) as.numeric(logLik(m))), 2),
    AIC    = round(sapply(list(fit_bs_k, fit_bs_df5, fit_ns_k, fit_ns_df5, fit3), AIC), 2)
)

##          modelo gl logLik    AIC
## 1      bs 2 nós  6 -40.07  94.14
## 2       bs df=5  6 -38.92  91.84
## 3      ns 2 nós  4 -65.70 141.39
## 4       ns df=5  6 -40.24  94.47
## 5 cúbica global  4 -66.57 143.13

Este resultado é instrutivo: ns() com os mesmos 2 nós de bs() tem log-verossimilhança bem pior (\(-65{,}7\) contra \(-40{,}1\)) — a restrição de linearidade nas bordas, com poucos nós, custa caro em ajuste. Mas ns(df=5) posiciona automaticamente 4 nós internos (mais que os 2 de bs(df=5)) e recupera um ajuste quase idêntico ao de bs() (\(-40{,}2\) contra \(-38{,}9\)). Ou seja: ns() geralmente precisa de mais nós que bs() para o mesmo ajuste, porque parte da flexibilidade vai para a restrição de linearidade nas bordas — a Seção “Extrapolação” mostra por que, mesmo assim, costuma valer a pena.

Funções base: o que bs() e ns() calculam

As colunas de bs()/ns() são funções-base fixas de \(x\) — o mesmo tipo de matplot() de colunas usado em xySpline.R para as bases de potências truncadas. Aqui estão as bases correspondentes aos ajustes acima:

par(mfrow = c(1, 2))
xs <- seq(0, 10, length = 501)
matplot(xs, bs(xs, knots = c(3.5, 7), Boundary.knots = c(0, 10)), type = "l",
        main = "bs(), 2 nós", xlab = "x", ylab = "")
abline(v = c(3.5, 7), lty = 3, col = "gray")
matplot(xs, ns(xs, knots = c(3.5, 7), Boundary.knots = c(0, 10)), type = "l",
        main = "ns(), 2 nós", xlab = "x", ylab = "")
abline(v = c(3.5, 7), lty = 3, col = "gray")

Funções base de bs() (esquerda) e ns() (direita), mesmos nós (3.5, 7)

par(mfrow = c(1, 1))

Repare como, na base de ns() (direita), cada função-base vira uma reta fora do intervalo \([0,10]\) — a restrição de linearidade “embutida” na própria base, visível diretamente, sem precisar olhar coeficiente algum.

Intervalos de confiança vs. intervalos de predição

predict(..., interval = "confidence") dá a incerteza sobre a função média \(E[Y\mid x]\); interval = "prediction" soma a variância residual e dá a incerteza sobre uma nova observação individual — por isso a banda de predição é sempre mais larga.

nd_ic <- data.frame(x = seq(0, 10, length = 201))
nd_ic <- transform(nd_ic, x1 = ifelse(x < 3.5, 0, x - 3.5), x2 = ifelse(x < 7, 0, x - 7))
ic_r3   <- predict(m3_r3, newdata = nd_ic, interval = "confidence")
ic_p4   <- predict(fit4,  newdata = nd_ic, interval = "confidence")
ip_r3   <- predict(m3_r3, newdata = nd_ic, interval = "prediction")
ip_p4   <- predict(fit4,  newdata = nd_ic, interval = "prediction")

par(mfrow = c(1, 2))
with(df, plot(Y ~ x, pch = 20, col = "gray40", main = "Confiança"))
matlines(nd_ic$x, ic_r3, col = "darkgreen", lty = c(1, 2, 2))
matlines(nd_ic$x, ic_p4, col = "steelblue", lty = c(1, 2, 2))
legend("topleft", c("spline cúbica (r3)", "quártica global"), col = c("darkgreen", "steelblue"), lty = 1, bty = "n", cex = 0.8)

with(df, plot(Y ~ x, pch = 20, col = "gray40", main = "Predição"))
matlines(nd_ic$x, ip_r3, col = "darkgreen", lty = c(1, 2, 2))
matlines(nd_ic$x, ip_p4, col = "steelblue", lty = c(1, 2, 2))

Intervalos de confiança (esquerda) e de predição (direita): spline cúbica (r3) vs. quártica global

par(mfrow = c(1, 1))

Dentro do domínio observado, as bandas de confiança da spline cúbica e da quártica global são parecidas (ambas capturam bem a curvatura); a diferença relevante entre elas aparece só quando se sai do domínio observado — assunto da próxima seção.

Extrapolação: onde tudo isso importa

Estendendo a predição de \([0,10]\) para \([-1,11]\) — só 1 unidade além de cada lado — compara-se bs() (2 nós), bs() (1 nó), ns(df=5) e o polinômio global de grau 4:

fit_bs1 <- lm(Y ~ bs(x, knots = 5, Boundary.knots = c(0, 10)), data = df)

nd_ex <- data.frame(x = seq(-1, 11, length = 401))
nd_ex <- transform(nd_ex, Ytrue = 0.2 * x + cos(x + 1))

ex_bs2 <- predict(fit_bs_k,  newdata = nd_ex, interval = "confidence")

## Warning in bs(x, degree = 3L, knots = c(3.5, 7), Boundary.knots = c(0, 10:
## alguns valores de 'x' depois dos nós de fronteira podem causar bases má
## condicionadas

ex_bs1 <- predict(fit_bs1,   newdata = nd_ex, interval = "confidence")

## Warning in bs(x, degree = 3L, knots = 5, Boundary.knots = c(0, 10), intercept =
## FALSE): alguns valores de 'x' depois dos nós de fronteira podem causar bases má
## condicionadas

ex_ns5 <- predict(fit_ns_df5, newdata = nd_ex, interval = "confidence")
ex_p4  <- predict(fit4,      newdata = nd_ex, interval = "confidence")

with(df, plot(Y ~ x, pch = 20, col = "gray40", xlim = c(-1, 11), ylim = c(-5, 15),
              main = "Extrapolação: bs() x ns() x polinômio"))
with(nd_ex, lines(Ytrue ~ x, lty = 3))
matlines(nd_ex$x, ex_bs2, col = "tomato", lty = c(1, 2, 2))
matlines(nd_ex$x, ex_bs1, col = "purple", lty = c(1, 2, 2))
matlines(nd_ex$x, ex_ns5, col = "darkgreen", lty = c(1, 2, 2))
matlines(nd_ex$x, ex_p4,  col = "steelblue", lty = c(1, 2, 2))
abline(v = c(0, 10), lty = 3, col = "gray60")
legend("top", c("bs 2 nós", "bs 1 nó", "ns df=5", "poly grau 4"),
       col = c("tomato", "purple", "darkgreen", "steelblue"), lty = 1, bty = "n", ncol = 4, cex = 0.75)

Extrapolação além do domínio [0,10]: bs() e polinômio disparam; ns() se mantém estável

bs() emite, no processo, um aviso: “alguns valores de ‘x’ depois dos nós de fronteira podem causar bases mal condicionadas” — um alerta explícito do próprio pacote de que ele não foi desenhado para extrapolar. ns() não emite esse aviso, exatamente porque foi desenhado para isso. Comparando as previsões (e bandas) 1 unidade além de cada extremo do domínio:

nd_pts <- data.frame(x = c(-1, 0, 10, 11))
tab_ex <- rbind(
    data.frame(modelo = "bs 2 nós",   round(predict(fit_bs_k,   newdata = nd_pts, interval = "confidence"), 2)),
    data.frame(modelo = "bs 1 nó",    round(predict(fit_bs1,    newdata = nd_pts, interval = "confidence"), 2)),
    data.frame(modelo = "ns df=5",    round(predict(fit_ns_df5, newdata = nd_pts, interval = "confidence"), 2)),
    data.frame(modelo = "poly grau4", round(predict(fit4,       newdata = nd_pts, interval = "confidence"), 2))
)

## Warning in bs(x, degree = 3L, knots = c(3.5, 7), Boundary.knots = c(0, 10:
## alguns valores de 'x' depois dos nós de fronteira podem causar bases má
## condicionadas

## Warning in bs(x, degree = 3L, knots = 5, Boundary.knots = c(0, 10), intercept =
## FALSE): alguns valores de 'x' depois dos nós de fronteira podem causar bases má
## condicionadas

tab_ex$x <- rep(nd_pts$x, 4)
tab_ex[, c("modelo", "x", "fit", "lwr", "upr")]

##        modelo  x  fit   lwr   upr
## 1    bs 2 nós -1 8.65  4.58 12.72
## 2    bs 2 nós  0 2.16  0.80  3.52
## 3    bs 2 nós 10 3.06  2.28  3.84
## 4    bs 2 nós 11 9.33  6.37 12.30
## 11    bs 1 nó -1 7.51  4.91 10.10
## 21    bs 1 nó  0 1.96  0.91  3.02
## 31    bs 1 nó 10 2.78  2.10  3.45
## 41    bs 1 nó 11 6.98  5.22  8.74
## 12    ns df=5 -1 2.50  0.49  4.50
## 22    ns df=5  0 0.95 -0.11  2.02
## 32    ns df=5 10 2.51  1.85  3.18
## 42    ns df=5 11 4.10  2.88  5.32
## 13 poly grau4 -1 9.95  6.89 13.01
## 23 poly grau4  0 2.37  1.24  3.49
## 33 poly grau4 10 2.96  2.26  3.66
## 43 poly grau4 11 8.99  6.89 11.09

Em \(x=-1\) (só 1 unidade fora da amostra), a spline bs() de 2 nós prevê \(\hat Y\approx 8{,}6\) e o polinômio de grau 4 prevê \(\hat Y\approx 10{,}0\) — ambos claramente absurdos, considerando que os valores observados nas proximidades de \(x=0\) giram em torno de \(1\)–\(2\). O ns(df=5) prevê \(\hat Y\approx 2{,}5\): ainda uma extrapolação (não é mágica — é só uma reta a partir do que foi observado perto da borda), mas muito mais razoável, com banda de confiança bem mais estreita. O mesmo padrão se repete em \(x=11\). Esta é a razão prática mais forte para preferir ns() a bs() (ou a um polinômio global) sempre que houver alguma chance de predição fora — ou perto das bordas — do domínio observado.

Síntese

Tópico	Ferramenta / resultado
Instabilidade numérica	base truncada crua tem \(\kappa(X'X)\) ~5 ordens de grandeza maior que bs(), para o mesmo modelo
Nós de bs()/ns()	knots= (manual) ou df= (automático, quantis de \(x\)); Boundary.knots deve refletir o domínio pretendido, não só a amostra
ns() vs. bs()	ns() força linearidade além das bordas — custa ajuste com poucos nós, mas estabiliza a extrapolação
Funções base	visíveis diretamente via matplot(bs(...))/matplot(ns(...)) — a restrição de ns() aparece literalmente como retas fora do domínio
IC vs. IP	confiança = incerteza da média; predição = confiança + variância residual (sempre mais larga)
Extrapolação	bs() e polinômios globais podem divergir fortemente fora do domínio; ns() é a opção estável por construção

Com isso, fecha-se o conjunto de tópicos sinalizados como pendentes ao longo da série (regressao-partes-texto.qmd, regressao-partes23-eolico.qmd) — os únicos tópicos ainda não cobertos são a extensão para múltiplos nós livres em graus maiores e a conexão com P-splines/penalização, já introduzida conceitualmente (mas não implementada) em regressao-partes-texto.qmd.