Capítulo 8

Regressão não paramétrica e modelos aditivos


A técnica estatística mais comumente usada, paramétrica ou não paramétrica, é a regressão. A premissa básica da regressão não paramétrica é a existência de uma função suave \(r(\cdot)\) que relaciona a resposta \(y\) e o preditor \(x\): \[ \tag{8.1} y_i=r(x_i)+\epsilon_i, \quad \mbox{para} \quad 1\leq i\leq n, \qquad \mbox{onde} \quad \epsilon_i\sim g\big(\mu,\sigma_\epsilon^2(x)\big)\cdot \] Frequentemente, assume-se que \(g\) segue uma distribuição normal e \(\sigma_\epsilon^2(x) = \sigma^2\), uma constante.

Neste capítulo, exploraremos a relação entre o estimador de densidade subjacente e o estimador de regressão. Muitos estimadores de regressão discutidos neste capítulo são suavizadores lineares, ou seja, combinações lineares das respostas observadas. Questões de robustez são relevantes, particularmente em dimensões mais elevadas.

Suavizadores não lineares não paramétricos serão discutidos. Para problemas de regressão em alta dimensão, modelos aditivos da forma \[ r(y_i) = \sum_{j=1}^d r_j(x_{ij}) \] têm recebido muita atenção, ver C. J. Stone (1985), Wahba (1990) e Hastie and Tibshirani (1990).

Uma implementação baseada em estimativas de histograma deslocado médio (ASH) será apresentada. Outras técnicas promissoras, incluindo o pouco conhecido método de regressão modal e o procedimento de regressão \(L_1\) não paramétrico, serão introduzidas. Vários outros livros fornecem tratamentos abrangentes de vários aspectos deste amplo campo, ver Eubank (1988), Müller (1988), Härdle (1990) e Simonoff (1996).



8.1 Regressão kernel não paramétrica


8.1.1 O estimador Nadaraya-Watson


Há dois casos distintos de (8.1) que devem ser considerados, dependendo da estrutura probabilística dos dados \(\{(x_i,y_i) \, : \, 1\leq i\leq n\}\). O primeiro caso ocorre quando os dados \(\{x_i\}\) provêm de um delineamento fixo; isto é, os dados \(x_i\) não são aleatórios, mas escolhidos pelo experimentador.

O segundo caso é o delineamento aleatório, que ocorre quando os dados provêm de uma função de densidade de probabilidade conjunta \(f(x,y)\). A ênfase neste capítulo será no último caso, em consonância com o foco deste livro.

A função de regressão teórica é definida como sendo \[ \tag{8.2} r(x)=\mbox{E}\big( Y \, | \, X=x \big)=\int y\, f(x \, | \,y) \mbox{d}y=\dfrac{\displaystyle \int y \, f(x,y)\mbox{d}y}{\displaystyle \int f(x,y)\mbox{d}y}\cdot \]

Considere a construção de um estimador de regressão não paramétrico obtido pelo cálculo de (8.2) usando o estimador kernel de produto bivariado \[ \widehat{f}(x,y)=\dfrac{1}{n\, h_x h_y}\sum_{i=1}^n K\Bigg(\dfrac{x-x_i}{h_x} \Bigg)K\Bigg(\dfrac{y-y_i}{h_y} \Bigg)=\dfrac{1}{n}\sum_{i=1}^n K_{h_x}(x-x_i)K_{h_y}(y-y_i) \] em vez da densidade bivariada desconhecida \(f(x,y)\).

O denominador em (8.2) contém a função de densidade marginal, que aqui se torna \[ \int \widehat{f}(x,y)\mbox{d}y=\dfrac{1}{n}\sum_{i=1}^n K_{h_x}(x-x_i)\int K_{h_y}(y-y_i)\mbox{d}y=\dfrac{1}{n}K_{h_x}(x-x_i) \] pois \(\displaystyle \int K_h(t)\mbox{d}t = 1\).

A restrição para um núcleo de ordem 2, de que \(\displaystyle \int t K_h (t)\mbox{d}t = 0\), implica que \(\displaystyle \int y K_h (y - y_i )\mbox{d}y = y_i\); portanto, o numerador se torna \[ \int y \, \widehat{f}(x,y)\mbox{d}y=\dfrac{1}{n}\sum_{i=1}^n y_i \, K_{h_i}(x-x_i)\cdot \]

Portanto, o estimador de regressão kernel não paramétrico correspondente a (8.2) é \[ \tag{8.3} \widehat{r}(x)=\dfrac{\displaystyle \dfrac{1}{n} \sum_{i=1}^n y_i \, K_{h_x}(x-x_i)}{\displaystyle \dfrac{1}{n} \sum_{i=1}^n K_{h_x}(x-x_i)}=\sum_{i=1}^n \omega_{h_x}(x,x_i)\, y_i, \] onde \[ \omega_{h_x}(x,x_i)=\dfrac{K_{h_x}(x-x_i)}{\displaystyle \dfrac{1}{n} \sum_{i=1}^n K_{h_x}(x-x_i)}\cdot \] Este estimador foi proposto por Nadaraya (1964) e Watson (1964).

Duas observações devem ser feitas. Primeiro, o estimador de regressão kernel (8.3) é linear nas observações \(\{y_i\}\) e, portanto, é um suavizador linear. Essa característica é compartilhada por muitos outros estimadores de regressão não paramétricos. Segundo, a estimativa de regressão kernel é independente da escolha específica do parâmetro de suavização \(h_y\).

Essa característica não é totalmente inesperada, dado o resultado anterior no Exercício 16, Capítulo 6, de que \[ \int x \, \widehat{f}_h (x)\mbox{d}x = \overline{x} \] para qualquer escolha de \(h\). No entanto, certas deficiências no modelo, como a falta de robustez da regressão kernel, decorrem dessa característica.

Finalmente, suponha que os dados não sejam aleatórios, mas provenham de um delineamento fixo \(\{x_i = i/n, \, i = 1,\cdots, n\}\), e que a função de regressão seja periódica em [0,1]. Se a estimativa de regressão for calculada nos pontos de delineamento, então o denominador na equação (8.3) é constante e o vetor de pesos \(\omega_h (x_j, x_i )\) precisa ser calculado apenas uma vez.

Como a avaliação repetida do kernel pode acarretar grande parte da carga computacional, o delineamento fixo com espaçamento uniforme é particularmente vantajoso. De fato, é possível simular essa carga computacional eficiente mesmo com um delineamento aleatório, “arredondando” os pontos de delineamento \(\{x_i\}\) para uma malha fixa com espaçamento uniforme. Essa proposta lembra o algoritmo ASH e foi investigada por Härdle and Scott (1988), ver Seção 8.4.1.


8.1.2 Estimadores polinomiais locais de mínimos quadrados


O uso de polinômios de ordem superior em regressão paramétrica para aproximar uma classe maior de curvas de regressão possíveis merece alguma consideração. Aqui, o grau do polinômio desempenha um papel correspondente ao parâmetro de suavização na regressão não paramétrica.

No entanto, os ajustes polinomiais paramétricos podem apresentar oscilações rápidas à medida que a ordem do polinômio aumenta. Assim, o uso de polinômios de ordem cada vez mais alta para regressão “não paramétrica” não parece prático. Se a curva de regressão verdadeira for suave, então um ajuste polinomial de baixa ordem deve ser adequado localmente (ver C. J. Stone 1977). Esta observação nada mais é do que uma reformulação do teorema de Taylor. Dois estimadores não paramétricos bem conhecidos emergem dessa linha de raciocínio.


8.1.2.1 Ajuste local de constantes


Globalmente, o melhor ajuste de regressão constante, polinômio de grau 0, a um diagrama de dispersão é \(\widehat{r}̂(x) =\overline{y}\), que é a estimativa que minimiza o critério dos mínimos quadrados: \[ \tag{8.5} \overline{y}=\arg \min_a \sum_{i=1}^n (a-y_i)^2; \] aqui, a notação \(\arg \min_a\) indica que a constante \(a = \overline{y}\) é a escolha (argumento) que minimiza o critério.

Considere um ajuste constante local em \(x\) aos dados. Aqui, “local” pode significar incluir apenas os dados \((x_i, y_i)\) para os quais \(x_i\in (x-h, x+h)\) na soma em (8.5), ou pode significar incluir apenas os \(k\) pontos de projeto mais próximos de \(x\).

É conveniente introduzir a função kernel \(K_h(x-x_i)\) para indicar precisamente quais termos estão incluídos e os pesos, se houver, desses termos. Assim, o melhor ajuste constante local é \[ \tag{8.6} \widehat{r}(x)=\arg \min_a \sum_{i=1}^n \Big( K_h(x-x_i)\times (a-y_i)^2\Big)\cdot \]

Minimizar o lado direito de (8.6) em relação a \(a\) leva a \[ \tag{8.7} \widehat{r}(x)=\dfrac{\displaystyle \sum_{i=1}^n y_i\, K_h(x-x_i)}{\displaystyle \sum_{j=1}^n K_h(x-x_j)}, \] que é precisamente o estimador de kernel de Nadaraya-Watson (8.3).

Este resultado pontual pode ser estendido a toda a função de regressão, observando que \[ \widehat{r}(\cdot) =\arg \min_{a(\cdot)} \int \sum_{i=1}^n K_h(x-x_i)\times \big(a(x)-y_i \big)^2 \mbox{d}x, \] visto que o integrando é minimizado por \(\widehat{a} (x) = \widehat{r}(x)\) em (8.7) para cada \(x\).


8.1.2.2 Ajuste polinomial local


Os ajustes de constantes locais na seção anterior são polinômios de ordem zero. O uso de ajustes polinomiais locais lineares ou quadráticos é intuitivamente atraente. O uso de polinômios de ordem superior localmente resulta em uma ordem de viés diferente, como ocorre com kernels de ordem superior.

Esse resultado foi demonstrado para o delineamento fixo por Härdle (1990) e para o delineamento aleatório por Fan (1992). O ajuste local de polinômios lineares foi especialmente defendido por Cleveland (1979) com seu procedimento LOWESS distribuído na linguagem S (Becker, Chambers, and Wilks 1988). Cleveland mostra que, por meio de uma organização cuidadosa do algoritmo, o trabalho de ajuste dos muitos polinômios locais pode ser minimizado adicionando e excluindo pontos de delineamento um de cada vez à soma dos quadrados. Wang (1990) considerou o ajuste polinomial local, mas com um critério de erro absoluto aplicado aos resíduos em vez de um critério de erro quadrático.


8.1.3 Erro quadrático médio pontual


As propriedades do erro quadrático médio (MSE) do estimador de Nadaraya-Watson são bastante complexas, pois o estimador é a razão entre duas variáveis aleatórias correlacionadas. O viés e a variância do estimador de Nadaraya-Watson na equação (8.3) podem ser obtidos usando aproximações para os erros padrão de funções de variáveis aleatórias (Stuart and Ord 1987).

Se o numerador e o denominador em (8.3) convergirem para uma constante (positiva), então a esperança assintótica da razão é a razão entre as esperanças assintóticas do numerador e do denominador, até a primeira ordem.

As propriedades do estimador kernel no denominador foram bem estudadas (ver equações (6.16) e (6.17)), com o resultado de que \[ \tag{8.8} \mbox{E}\big(\widehat{f}(x) \big)\approx f(x)+h^2 \sigma_K^2 \, f''(x)/2 \qquad \mbox{e} \qquad \mbox{Var}\big(\widehat{f}(x) \big)\approx \dfrac{R(K) f(x)}{nh}\cdot \]

Em seguida, considere a esperança do numerador: \[ \tag{8.9} \begin{array}{rcl} \displaystyle \mbox{E}\left(\sum_{i=1}^n y_i\, K_h(x-x_i) \right) & = & \displaystyle \int \int \nu \, \dfrac{1}{h} K\Bigg( \dfrac{x-u}{h}\Bigg) f(u,\nu)\, \mbox{d}u \,\mbox{d}\nu \\[0.8em] & = & \displaystyle \int \int \nu \, K(s) f(x-hs,\nu) \, \mbox{d}s \,\mbox{d}\nu, \end{array} \] após a mudança de variável \(s = (x-\nu)/h\).

Agora a densidade condicional satisfaz \(f(\nu \, | \, x-hs)f(x-hs) = f (x-hs, \nu)\). Portanto, a integral sobre \(\nu\) em (8.9) é igual a (ignorando \(K(s)\)) \[ f(x-hs) \int \nu \, f(\nu \, | \, x-hs) \, \mbox{d}\nu = f(x-hs) \, r(x-hs), \] como a integral é a média condicional e \(r(\cdot)\) é a verdadeira função de regressão definida em (8.1).

Continuando com (8.9), temos \[ \tag{8.10} \begin{array}{rcl} & = & \displaystyle \int K(s)f(x-hs) \, r(x-hs) \mbox{d}s \\[0.8em] & = & \displaystyle f(x) \, r(x) + h^2 \sigma_K^2 \Big(f'(x) r'(x) + f''(x) r(x)/2 + f(x)r''(x)/2 \Big) \end{array} \] após expandir \(f(x-hs)\) e \(r(x-hs)\) em séries de Taylor até a ordem \(h^2\), assumindo que \(K\) é um kernel de segunda ordem.

Assim, a esperança do estimador de Nadaraya-Watson com um delineamento aleatório é a razão entre as esperanças em (8.10) e (8.8), com a densidade \(f\) fatorada de cada uma: \[ \tag{8.11} \begin{array}{rcl} \mbox{E}\big(\widehat{f}(x) \big) & = & \displaystyle \dfrac{f(x)\times \Big(r(x)+h^2 \sigma_K^2 \big(f' r'/f + f'' r /(2f)+r''/2 \big) \Big)}{f(x)\times \Big(1+h^2\sigma_K^2 \, f''(x)/(2f) \Big)} \\[0.8em] & \approx & \displaystyle r(x)+\dfrac{1}{2}h^2 \sigma_K^2 \Bigg(r''(x)+ 2r' (x)\dfrac{f'(x)}{f(x)} \Bigg) \end{array} \] usando qualquer programa de manipulação simbólica, ou a aproximação \[ (1+h^2 c)^{-1} \approx (1-h^2 c) \] para \(h\approx 0\) no fator do denominador e multiplicando.

O viés no caso de projeto fixo é \(h^2 \sigma_K^2 r''(x)/2\), que aparece na equação (8.11), (veja o Exercício 1). O termo \(2r'(x)f'(x)/f(x)\) em (8.11) é pequeno se houver muitos pontos de dados na janela, ou seja, \(f(x)\) grande. Se houver mais pontos de projeto no intervalo \((x,x+h)\) do que em \((x-h, x)\), ou seja, \(f'(x) > 0\), então a média local será enviesada, pois a média incluirá mais respostas em \((x, x+h)\) do que em \((x-h, x)\). O viés será positivo se \(r'(x) >\) 0 e negativo caso contrário.

Uma interpretação semelhante existe quando \(f'(x) < 0\). A interação das várias possibilidades para os sinais de \(f'(x)\) e \(r'(x)\) é resumida nesse termo. Embora a expressão para o viés com um desenho fixo cometerá com o desenho aleatório quando \(f'(x) = 0\), observe que essas duas configurações não são idênticas. O desenho aleatório tem probabilidade zero de ser igualmente espaçado mesmo quando \(f = U(0, 1)\).

A variância de \(\widehat{r}(x)\) pode ser calculada usando a seguinte aproximação para a variância da razão de duas variáveis aleatórias (Stuart and Ord 1987): \[ \tag{8.12} \mbox{Var}\Bigg(\dfrac{U}{V} \Bigg)=\Bigg(\dfrac{\mbox{E}(U)}{\mbox{E}(V)} \Bigg)^2 \Bigg(\dfrac{\mbox{Var}(U)}{\big(\mbox{E}(U)\big)^2}+\dfrac{\mbox{Var}(V)}{\big(\mbox{E}(V)\big)^2}-\dfrac{2\, \mbox{Cov}(U,V)}{\mbox{E}(U)\mbox{E}(V)} \Bigg)\cdot \]

Cálculos semelhantes aos que seguem a equação (8.9) mostram que \[ \tag{8.13} \begin{array}{rcl} \displaystyle \mbox{Var}\Bigg(\dfrac{1}{n} \sum_{i=1}^n y_i\, K_h(x-x_i) \Bigg) & = & \displaystyle \dfrac{1}{n}\mbox{E}\big( K_h(x-x_i) \, y_i\big)^2 - O(n^{-1}) \\[0.8em] & \approx & \displaystyle \dfrac{R(K)f(x)}{nh}\big( \sigma_\epsilon^2 +r(x)^2\big), \end{array} \] utilizando o fato de \[ \int \nu^2 f(\nu \, | \, x-hs)\mbox{d}\nu = \sigma_\epsilon^2 (x-hs) + r(x-hs)^2 \] e \(\sigma_\epsilon^2(x)=\sigma_\epsilon^2\) para todo \(x\). A variância de \(\widehat{f}(x)\) é dada pela equação (8.8). Finalmente, \[ \tag{8.14} \mbox{Cov}\Bigg(\dfrac{1}{n}\sum_{i=1}^n K-h(x-x_i) \, Y_i-\dfrac{1}{n}\sum_{i=1}^n K_h(x-x_i) \Bigg) = \\[0.8em] \qquad \qquad \qquad \qquad =\dfrac{1}{n}\mbox{E}\big(K_h(x-x_i) \, Y_i\big)-O(n^{-1}) \approx \dfrac{R(K)f(x)r(x)}{nh}\cdot \]

Substituindo cuidadosamente as equações (8.8), (8.10), (8.13) e (8.14) em (8.12), obtém-se: \[ \mbox{Var}\big(\widehat{r}(x) \big)\approx \dfrac{R(K)\, \sigma_\epsilon^2}{n\, hf(x)}\cdot \]

Além do fator conhecido \(R(K)/(nh)\), a variância de \(\widehat{r}(x)\) inclui fatores relacionados à variância do ruído \(\sigma_\epsilon^2\) e à quantidade (relativa) de dados por meio de \(f(x)\).

Esses resultados podem ser reunidos no seguinte teorema (Rosenblatt 1969). Uma versão multivariada é apresentada por Mack and Müller (1989).


Teorema 8.1:

O \(\mbox{MSE}(x)\) assintótico do estimador de Nadaraya-Watson é \[ \tag{8.15} \mbox{AMSE}\big(\widehat{r}(x) \big)=\dfrac{R(K)\, \sigma_\epsilon^2}{n \, h f(x)}+\dfrac{1}{4}h^4 \sigma_K^4 \Big( r''(x)+2r'(x)\dfrac{f'(x)}{f(x)} \Bigg)^2\cdot \]


Demonstração. Scott (2015).


O fato de o estimador de Nadaraya-Watson decorrer diretamente do estimador kernel de produto em sua forma funcional merece uma análise mais detalhada. Para o caso da densidade bivariada, \(h^*_x = O(n^{-1/6})\); entretanto, ocorre que \(h^*_x = O(n^{-1/5})\) para este estimador específico, que é a taxa para a estimação da densidade univariada. Assim, o suavizamento adicional fornecido pela integração da função de densidade bivariada dispensa a necessidade de um parâmetro de suavização mais amplo. Diversos autores consideraram pesos diferentes dos de Nadaraya-Watson e obtiveram expressões de viés significativamente mais simples.

O esquema de pesos integrais ou de convolução desenvolvido por Gasser and Müller (1979) apresenta certas vantagens (Gasser and Engel 1990). Além de uma expressão de viés semelhante à do caso de projeto fixo, a forma do viés é favorável quando \(r(x)\) é linear, mas a variância é aumentada. Fan (1992) mostrou que a técnica de ajuste polinomial local alcança simultaneamente boas propriedades de viés e variância, mesmo com o projeto aleatório.


8.1.4 Seleção da largura de banda


A seleção de parâmetros de suavização por validação cruzada em regressão não paramétrica é mais fácil do que na estimação de densidade, até mesmo porque a curva de regressão passa pela nuvem de pontos de dados. Duas classes gerais de algoritmos surgiram com base em modificações simples do erro quadrático médio preditivo ingênuo, que é definido por \[ G(h)=\dfrac{1}{n}\sum_{i=1}^n \big( y_i-\widehat{r}(x_i)\big)^2\cdot \]

Denotando \(\widehat{r}_{-i}(\cdot)\) como o estimador de Nadaraya-Watson dos \(n-1\) pontos com o ponto \((x_i , y_i )\) omitido, Allen (1974), M. Stone (1974), Clark (1975) e Wahba and Wold (1975) propuseram escolher \[ \widehat{h}_{CV} =\arg \min_h \mbox{CV}(h), \qquad \mbox{onde} \qquad \mbox{CV}(h)=\dfrac{1}{n}\sum_{i=1}^n \big( y_i-\widehat{r}_{-i}(x_i)\big)^2\cdot \]

Ao omitir o \(i\)-ésimo ponto de dados, elimina-se a solução degenerada \(h = 0\), que corresponde à interpolação dos pontos de dados e a estimativa \(\mbox{CV}(h) = 0\), uma subestimação significativa do verdadeiro erro preditivo.

Alternativamente, a subestimação do erro preditivo pode ser evitada pela multiplicação por um fator simples: \[ \tag{8.16} \widehat{h}=\arg \min_h \Bigg( 1+\dfrac{2K(0)}{nh}\Bigg)\times G(h)\cdot \]

Pode-se demonstrar que essa escolha fornece um estimador consistente do erro preditivo. Essa fórmula específica em (8.16) é devida a Shibata (1981). Rice (1984) usa o fator \(\big(1-2K(0)/(nh)\big)^{-1}\). Härdle, Hall, and Marron (1988) mostraram como cinco algoritmos de validação cruzada têm a mesma série de Taylor como critério (8.16).


8.1.5 Suavização adaptativa


Uma breve revisão do Teorema 8.1 leva à conclusão de que a construção de curvas de regressão adaptativa assintoticamente ótimas é difícil, certamente mais difícil do que no caso da estimação de densidade.

O uso dos pesos de Gasser-Müller simplifica um pouco as coisas, uma vez que o termo envolvendo \(r'\) não aparece na expansão do erro quadrático médio assintótico (\(\mbox{AMSE}\)). Müller (1988) descreve vários algoritmos para estimação adaptativa, assim como Staniswalis (1989) e Schucany (1995).

Um dos suavizadores adaptativos mais inovadores é o supersuavizador de Friedman (1984). O algoritmo constrói estimativas supersuavizadas, médias e subsuavizadas, denominadas “woofer”, “midrange” e “tweeter”, respectivamente. Com nove passagens pelos dados e o uso de validação cruzada local, uma estimativa adaptativa é construída.


8.2 Estimação linear não paramétrica geral


A estimativa de regressão kernel é uma combinação linear das respostas observadas na forma \(\pmb{\omega}^\top \pmb{y}\), que é um suavizador linear por definição. Na forma vetorial, os dados de regressão bivariada \(\{(x_i , y_i )\}\) serão denotados por \[ \pmb{x}=(x_1,\cdots,x_n)^\top \qquad \mbox{e} \qquad \pmb{y}=(y_1,\cdots,y_n)^\top\cdot \] Os pesos para outros estimadores lineares são derivados nas próximas três seções.


8.2.1 Regressão polinomial local


O estudo da regressão polinomial local é apenas ligeiramente mais complexo do que o caso global. Portanto, o caso global será considerado primeiro.

Qualquer regressão polinomial por mínimos quadrados pode ser facilmente demonstrada como “linear” no sentido descrito acima, mesmo quando os termos polinomiais são quadráticos ou cúbicos. Considere o caso da regressão linear com um termo de intercepto, que pode ser incorporado ao modelo adicionando uma coluna de uns aos pontos de projeto: \[ \pmb{X}=\begin{pmatrix} 1 & 1 & \cdots & 1 \\ x_1 & x_2 & \cdots & x_n \end{pmatrix}^\top \qquad \mbox{e} \qquad \beta=(\beta_0,\beta_1)^\top\cdot \] Então o problema de regressão pode ser escrito como \[ \tag{8.17} \pmb{y}=\pmb{X}\beta+\epsilon, \] do qual obtemos como estimativas \[ \widehat{\beta}=\big(\pmb{X}^\top \pmb{X} \big)^{-1} \pmb{X}^\top\pmb{y}, \] que é a conhecida estimativa de mínimos quadrados de \(\beta\).

Portanto, o melhor preditor de \(\pmb{y}\) em \(\pmb{x}\) é \[ \tag{8.18} \widehat{y}(x)=\widehat{r}(x)=\pmb{\omega}_x^\top\pmb{y} \] onde \[ \pmb{\omega}_x^\top \pmb{y}= \begin{pmatrix} 1 & x \end{pmatrix} \begin{pmatrix} \pmb{X}^\top \pmb{X}\end{pmatrix}^{-1}\pmb{X}^\top \pmb{y}\cdot \]

O preditor não só é uma combinação linear das respostas \(\{y_i\}\), como os próprios pesos se alinham em uma linha reta quando mostrados graficamente no plano \((x, y)\). Isso pode ser observado empiricamente na Figura 8.1, mas decorre da equação (8.18) ao examinar a forma como \(x\) entra na definição do vetor de pesos \(\pmb{\omega}_x\) (ver Exercício 3).

Figura 8.1: Suavizadores lineares aplicados ao conjunto de dados de fluxo de gás a 74.6 psia. No primeiro quadro, os 33 vetores de peso são plotados para um ajuste paramétrico de linha reta, com o círculo mostrando o peso no próprio ponto de dados. No segundo quadro, os dados brutos (“+”) são mostrados conectados por pequenos traços; as 33 estimativas em \(x = x_i\) do quadro 1 são mostradas (os símbolos “o” na linha reta conectados por traços grandes); as 33 estimativas “o” conectadas por uma linha contínua usando um modelo quadrático local com \(k = 7\), três pontos de cada lado de \(x_i\), exceto nos limites. No terceiro quadro, os 33 vetores de peso para os ajustes quadráticos locais são mostrados com o círculo mostrando o peso em \(x_i\). O quarto quadro mostra os 33 ajustes quadráticos locais sobrepostos aos dados.

Essa demonstração é bastante geral, pois a matriz de projeto \(\pmb{X}\) pode representar uma regressão polinomial de qualquer ordem, com múltiplas variáveis preditoras e com potências arbitrárias dessas variáveis. Assim, o modelo linear geral no contexto paramétrico também é um estimador linear de acordo com a definição não paramétrica.

Para algoritmos de regressão polinomial não paramétrica local, como o LOWESS, a demonstração acima é facilmente estendida. Em vez de um único modelo linear da forma (8.17) que vale para todo \(x\), um modelo linear local é formado para cada ponto \(x\). Apenas os pontos de projeto \(x_i\) “próximos” a \(x\) são incluídos no modelo linear local. O número de pontos incluídos no ajuste linear local desempenha o papel de parâmetro de suavização, embora o grau do polinômio também seja importante. O número de pontos incluídos no ajuste local pode ser designado diretamente como \(k\) ou indiretamente como uma fração \(f\) do tamanho da amostra. Por exemplo, o LOWESS escolhe os \(f\times n\) pontos mais próximos de \(x\); ou seja, \(k = fn\).


8.2.2 Suavização por spline


A construção de curvas suaves que passam pelos pontos de projeto \(\{(x_i , y_i ) \, , \, 1\leq i\leq n\}\) era um problema prático enfrentado pelos projetistas de cascos de navios. A solução encontrada foi utilizar um pedaço de madeira longo, fino e flexível, fixado aos pontos de projeto de forma a minimizar a deformação na madeira. Esse dispositivo mecânico foi denominado spline.

Assim, os matemáticos que consideravam esse problema de interpolação partiram de um modelo matemático da deformação por flexão dada pela fórmula \[ \int r'' (x)^2 \mbox{d}x\cdot \] A solução para o problema de interpolação \[ \min_r \int_{x_1}^{x_n} r''(x)^2 \mbox{d}x \qquad \mbox{restrito à} \qquad r(x_i)=y_i, \quad i=1,\dots,n \] é chamado de interpolação por spline.

Schoenberg (1964) mostrou que a interpolação por spline \(r(\cdot)\) é um polinômio cúbico definido por partes entre nós adjacentes \((x_i , x_{i+1} )\) de modo que \(r\) é contínuo em \((x_1 , x_n )\). Esta última afirmação é denotada por \(r\in C^2(x_1 , x_n )\).

A versão estatística desse problema consiste em encontrar a função de regressão mais suave com uma soma de quadrados dos resíduos fixa.

O spline de suavização é a solução para o problema variacional \[ \tag{8.19} \widehat{r}_\lambda (\cdot) = \arg \min_r S_\lambda (r) = \sum_{i=1}^n \big(y_i-r(x_i) \big)^2 +\lambda \int_{x_1}^{x_n} r''(x)^2 \mbox{d}x\cdot \]

Reinsch (1967) demonstrou que a solução de (8.19) também é um spline cúbico. O limite \[ \lim_{\lambda\to 0} \{\widehat{r}_\lambda\} \] é a interpolação spline, enquanto o limite de \[ \lim_{\lambda\to \infty} \{\widehat{r}_λ\} \] é o ajuste linear de mínimos quadrados aos dados, observando que a deformação por flexão de uma linha é 0.

O parâmetro \(\lambda\) desempenha o papel de parâmetro de suavização. Dado \(\lambda\), a quantidade de computação necessária para encontrar \(\widehat{r}_\lambda\) é \(O(n)\). Pacotes como o \(\mbox{IMSL}\) (1991) fornecem software para esse propósito.

Diversas propriedades do spline de suavização podem ser derivadas sem uma solução explícita. Para qualquer função \(q\in C^2 (x_1 , x_n )\) e constante \(\alpha\), \[ S_\lambda (\widehat{r}+\alpha q) \geq S_\lambda (\widehat{r}), \qquad \forall q\in C^2(x_1,x_n)\cdot \] Portanto, \(T(\alpha) = S_\lambda (\widehat{r}+\alpha q)\) possui um mínimo local em \(\alpha= 0\) para todo \(q\). Pela equação de Euler-Lagrange, a derivada de \(T\) deve se anular nesse ponto. Agora \[ \begin{array}{rcl} T'(\alpha) & = & \displaystyle \dfrac{\partial}{\partial \alpha} \Bigg( \sum_{i=1}^n \big(y_i-(\widehat{r}+\alpha q)(x_i) \big)^2 +\lambda \int_{x_1}^{x_n} \big((\widehat{r}+\alpha q)''(x) \big)^2 \mbox{d}x \Bigg) \\[0.8em] & = & \displaystyle \dfrac{\partial}{\partial \alpha} \Bigg( \sum_{i=1}^n \big(y_i-\widehat{r}(x_i)-\alpha q(x_i) \big)^2 +\lambda \int_{x_1}^{x_n} (\widehat{r}''^2+2\alpha \widehat{r}'' q'' +\alpha^2 q''^2\big) \mbox{d}x \Bigg)\cdot \end{array} \]

Realizando a diferenciação, a equação \(T'(\alpha) = 0\) quando \(\alpha = 0\) torna-se \[ \tag{8.20} -2 \sum_{i=1}^n \big(y_i-\widehat{r}(x_i) \big) q(x_i) +2\lambda \int_{x_1}^{x_n} \widehat{r}''(x) q''(x)\mbox{d}x=0\cdot \]

Agora, apresentaremos um breve argumento para demonstrar que \(\widehat{r}_\lambda\) é um suavizador linear. Suponha que existam dois vetores de resposta, \(\pmb{y}^{(1)}\) e \(\pmb{y}^{(2)}\), cada um com o mesmo vetor de projeto \(\pmb{x}\). Sejam \(\widehat{r}_\lambda^{(1)}\) e \(\widehat{r}_\lambda^{(2)}\) os respectivos splines de suavização. Então, se o spline de suavização for um suavizador linear, o spline de suavização para a soma dos dois vetores de resposta dos dados deverá ser a soma dos splines de suavização individuais; ou seja, \[ \tag{8.21} \big( \pmb{x},\pmb{y}^{(1)}+\pmb{y}^{(2)}\big) \qquad \mbox{implica que} \qquad \widehat{r}_\lambda^{(1+2)}=\widehat{r}_\lambda^{(1)}+\widehat{r}_\lambda^{(2)}\cdot \]

Como \(\big\{y_i^{(1)},y_i^{(2)},r_\lambda^{(1)},r_\lambda^{(2)}\big\}\) ao entrar linearmente na equação (8.20), a equação de Euler-Lagrange é válida para a solução proposta na equação (8.21) por aditividade (ver Exercício 8.2).

Como consequência da propriedade de suavização linear da spline de suavização, a solução para o conjunto de dados original \((\pmb{x},\pmb{y})\) pode ser considerada como a soma das soluções de \(n\) problemas individuais, como \[ \pmb{y}=\sum_{i=1}^n y_i \pmb{e}_i \qquad \mbox{onde} \qquad \pmb{e}_i=(0,\cdots,0,1,0,\cdots,0)^\top\cdot \]

Se \(\widehat{r}_\lambda^{(i)}\) é o spline de suavização para os dados \((\pmb{x},\pmb{e}_i)\), então a solução do spline para os dados originais \((\pmb{x},\pmb{y})\) é dado por \[ \widehat{r}_\lambda =\sum_{i=1}^n y_i \times \widehat{r}_\lambda ^{(i)}\cdot \]

Se o desenho for periódico, também chamado circular, e igualmente espaçado, então todos os \(n\) splines de suavização de componentes \(\widehat{r}_\lambda^{(i)}\) terão a mesma forma. Silverman (1984) demonstrou que a forma convergia para o kernel \[ \tag{8.22} K_s(t)=\dfrac{1}{2}e^{-|t|/\sqrt{2}}\sin\Bigg( \dfrac{\pi}{4}+\dfrac{|t|}{\sqrt{2}}\Bigg), \qquad \mbox{para} \qquad |t|<\infty\cdot \]

Pode-se demonstrar que \(K_s(t)\) é um kernel de ordem 4. Para dados provenientes de um delineamento aleatório, Silverman provou que o spline de suavização é um estimador de largura de banda adaptativo, e não fixo, com pesos \[ \omega_\lambda(x,x_i)\approx \dfrac{1}{f(x)h(x_i)}K_s\Bigg( \dfrac{x-x_i}{h(x_i)}\Bigg), \qquad \mbox{onde} \qquad h(x_i)\approx \Bigg(\dfrac{\lambda}{nf(x)} \Bigg)^{1/4}\cdot \]

Um kernel de quarta ordem deve ter largura de banda \(h(x_i) = O(n^{-1/9})\), o que ocorre se \(\lambda = O(n^{5/9})\) for escolhido. Infelizmente, o suavizador spline se adapta apenas à densidade de projeto \(f(x)\) e não à curvatura da curva desconhecida \(r(x)\), conforme exigido pelo Teorema 8.1.

Embora o estimador possa não fornecer uma estimativa adaptativa ótima, ele fornece estimativas visualmente atraentes na prática. A proposta original de Reinsch (1967) incluía uma função de ponderação em (8.19), que poderia ser projetada para fornecer propriedades adaptativas aprimoradas.


8.2.3 Kernels equivalentes


Uma ampla gama de algoritmos de regressão não paramétricos, incluindo LOWESS e splines de suavização, demonstraram ser suavizadores lineares. Para comparar essas estimativas com o estimador kernel original em (8.3), é natural comparar os vetores de peso \(\pmb{\omega}_x\) que multiplicam o vetor de resposta \(\pmb{y}\). Seguindo a definição dos pesos do kernel em (8.4), os vetores de peso para outros suavizadores lineares serão chamados de pesos de kernel equivalentes desse estimador. Para um conjunto de dados específico, esses pesos podem ser calculados exatamente para fins teóricos. Na prática, não é necessário obter os pesos explicitamente para calcular a estimativa.

Considere o conjunto de dados de precisão do fluxo de gás com \(n = 33\) leituras obtidas quando a pressão na linha de gás foi ajustada para 74.6 psia. A precisão do medidor é representada graficamente em função do logaritmo da taxa de fluxo de gás na tubulação. É evidente, a partir do painel superior direito da Figura 8.1, que a relação não é linear.

O primeiro painel mostra o conjunto de pesos equivalentes do kernel para o ajuste paramétrico de reta, e a estimativa é mostrada no quadro superior direito. Por exemplo, considere o ponto de dados mais à esquerda (1.76; 97.25). Os 33 pesos das 33 respostas usados para calcular a estimativa de regressão foram calculados de acordo com (8.18). Esses pesos se alinham em uma reta representada no quadro superior esquerdo da Figura 8.1, conectando os dois pontos extremos \((x,\omega) = (1.76; 0.145)\) e \((3.59; -0.081)\). Para identificar as 33 linhas, o peso \((x_i,\omega_i)\) é mostrado como um pequeno círculo na linha que representa o \(i\)-ésimo vetor de peso. Observe que os pesos não são locais, ou seja, os pesos não desaparecem fora de uma vizinhança de \(x_i\).

Um ajuste polinomial quadrático local também foi aplicado a esses dados. O ajuste foi baseado em \(k = 7\) pontos de projeto, \(\{x_{i-3},\cdots,x_i,\cdots,x_{i+3}\}\). Nos limites, os sete pontos mais próximos do limite são usados para o ajuste. Os 33 ajustes quadráticos locais são mostrados no quadro inferior direito da Figura 8.1, e os 33 vetores de peso equivalentes são mostrados no quadro inferior esquerdo. A pequena largura de banda permite que a estimativa acompanhe a pequena queda no meio da curva. Os vetores de peso equivalentes são 0 onde não são mostrados explicitamente. Por construção, os pesos são locais. Embora os kernels equivalentes sejam bastante irregulares e assumam valores negativos, parte dessa irregularidade resulta de sua capacidade de se adaptar corretamente nos limites (Fan 1992). Compare os limites dos kernels aqui com aqueles mostrados na Figura 6.10.

Os pesos do kernel de Nadaraya-Watson são mostrados na Figura 8.2 para o kernel biweight com uma largura de banda de \(h = 0.25\), que fornece um estimador mais suave do que o ajuste polinomial local. No entanto, um kernel de fronteira não foi empregado, como é evidente na estimativa na fronteira esquerda e na figura dos pesos do kernel.

Figura 8.2: Pesos e estimativa kernel Biweight Nadaraya-Watson para o conjunto de dados de fluxo de gás.

O ajuste do spline de suavização é mostrado na Figura 8.3. Observe que os pesos do kernel não são tão locais para esta pequena amostra porque a estimativa é muito mais suave (menos local). O comportamento na fronteira é razoável. Para uma discussão sobre o dual das larguras de banda equivalentes para suavizadores lineares, veja Hastie and Tibshirani (1990).

Figura 8.3: Pesos equivalentes nas estimativas do kernel e do spline de suavização para o conjunto de dados de fluxo de gás.


8.3 Robustez


Nos últimos anos, o esforço necessário para a entrada e verificação de dados tem sido maior do que o exigido para uma regressão completa e análise de variância. Embora as mudanças na computação tenham sido drásticas, os problemas de lidar com pontos de dados ruins ou influentes cresceram ainda mais rápido do que o tamanho dos conjuntos de dados.

Nesses casos, torna-se cada vez mais difícil identificar centenas de pontos influentes, particularmente em dimensões mais altas. Assim, a importância prática de procedimentos robustos que lidam “automaticamente” com outliers está em constante crescimento (Rousseeuw and Leroy 1987).


8.3.1 Estimadores resistentes


Como os suavizadores lineares são equivalentes a uma média ponderada local das respostas, esses estimadores não são resistentes a outliers no vetor de resposta. Seguindo a literatura bem desenvolvida em estimação robusta no contexto paramétrico (Huber 1964; Hampel 1974), surgiram duas propostas para melhorar a resistência do estimador kernel, com foco na caracterização dada na equação (8.6).

A função quadrática dos resíduos é substituída por uma função de influência, \(\rho(\cdot)\), que reduz a magnitude da contribuição de resíduos grandes ao critério de otimização. Uma escolha interessante é substituir a função de influência quadrática ilimitada, \(\rho(\epsilon) = \epsilon^2\), por uma função que tenha influência limitada; isto é, \(|\rho(\epsilon)| < c < \infty\).

Dada uma escolha para a função de influência robusta, o problema original é substituído por \[ \widehat{r}_\rho (\cdot)=\arg \min_{a(\cdot)} \int \sum_{i=1}^n K_h(x-x_i)\rho\big(a(x)-y_i\big) \mbox{d}x\cdot \] Este estimador foi discutido por Härdle (1984) e Härdle e Gasser (1984).

Um segundo algoritmo que reduz a influência de resíduos grandes foi proposto por Cleveland (1979). O procedimento LOWESS descrito na Seção 8.1.2 pode ser estendido para esse propósito. Em vez de substituir a função de influência quadrática dos resíduos \(\{\epsilon_i\}\), um fator multiplicativo adicional é introduzido termo a termo com base na magnitude relativa dos resíduos amostrais \(\{\widehat{\epsilon}_i\}\). O fator é pequeno para resíduos grandes. Cleveland sugere calcular uma estimativa de escala robusta \(\widehat{\sigma}\) e ponderar \(\{\delta_i\}\) por \[ \widehat{\sigma}=\mbox{mediana}\{|\widehat{\epsilon}_i|\} \qquad \mbox{o qual implica} \qquad \delta_i =\Big(1-\big(\widehat{\epsilon}_i/(6\widehat{\sigma}) \big)^2 \Big)_+^2\cdot \]

Por exemplo, ao ajustar um polinômio linear local, o polinômio ajustado resolve \[ (\widehat{a},\widehat{b})=\arg \min_{(a,b)} \sum_{i=1}^n \delta_i \times \omega_h(x,x_i)\big(y_i-(a+b\, x) \big)^2\cdot \]

A estimativa LOWESS é \(\widehat{a}+\widehat{b}x\). O procedimento é aplicado recursivamente até que as estimativas permaneçam inalteradas. O esforço computacional é uma ordem de magnitude maior do que a versão não robusta.


8.3.2 Regressão modal


É irônico que, embora o estimador de densidade kernel original seja bastante resistente a outliers, a função de regressão derivada da estimativa kernel não o seja. Existe algum estimador de regressão alternativo derivado do estimador kernel que seja resistente?

A resposta surge de uma breve reflexão sobre as três medidas de tendência central discutidas em livros didáticos básicos: média, mediana e moda. A proposta óbvia é ir além da média condicional do estimador kernel e considerar sua mediana condicional ou sua moda condicional. É claro que, se os erros forem normais ou simétricos, todas as três escolhas levam ao mesmo resultado, assintoticamente.

Um argumento forte pode ser feito para resumir os dados bivariados com um traço da moda condicional. Esse traço será chamado de curva de regressão modal ou traço modal. A definição formal é dada por \[ \mbox{Curva de regressão modal:} \quad \widehat{r}(x)=\arg(s) \max_y \widehat{f}(y\, | \, x), \] onde o plural em “\(\arg\)” indica todos os máximos locais.

Uma definição equivalente é: \[ \widehat{r}(x)=\arg(s) \max_y \widehat{f}(y,x)\cdot \]

Alguns exemplos revelam as vantagens práticas da moda condicional em comparação com os estimadores de média condicional, tanto robustos quanto não robustos. Considere o conjunto de dados de previsão de erupção do Old Faithful, apresentado na Figura 8.4, juntamente com o suavizador LOWESS. Devido ao comportamento multimodal (condicional) para \(x > 3\), a estimativa kernel resistente LOWESS é, na verdade, mais grosseira do que uma estimativa LOWESS não resistente.

Mas, como suavizador de diagramas de dispersão e resumidor de dados, qualquer suavizador linear é inadequado para resumir a estrutura. O quadro à direita da figura mostra o traço de regressão modal, bem como o ASH bivariado do qual foi derivado. Como medida de centro, a moda resume os valores condicionais “mais prováveis” em vez da média condicional. Por outro lado, quando a densidade condicional é simétrica, esses dois critérios coincidem. Neste exemplo, eles não coincidem.

Figura 8.4: Suavizadores de média condicional (LOWESS) e moda condicional do conjunto de dados de duração defasada do Old Faithful. A moda condicional é exibida com o símbolo “o” quando acima do contorno de 25% e com um “\(x\)” entre 5 e 25%. A linha de 45° também é mostrada.

Observe que a trajetória modal é muito mais suave do que a regressão usual. Essa suavidade pode parecer contraintuitiva à primeira vista. Compare as duas estimativas de regressão no intervalo \((1.; 2.5)\), no qual as densidades condicionais são unimodais, mas não simétricas.

A estimativa LOWESS segue fielmente a média condicional, que apresenta uma queda. No intervalo \((3.5; 5)\), as densidades condicionais são bimodais. A bifurcação da trajetória de regressão modal indica que uma única previsão, como a média condicional, não captura a estrutura dos dados. A moda condicional é quase meio minuto mais longa se a erupção anterior foi curta. Observe que a estimativa de regressão ordinária fornece previsões nesse intervalo onde há pouca massa de probabilidade, em torno de \(y = 3\) min.

Um segundo exemplo fornece mais evidências da utilidade dessa abordagem. A Figura 8.5 exibe o conjunto de dados de espessura da moeda de um centavo de dólar Lincoln dos EUA. A regressão kernel, novamente representada por LOWESS, rastreia as mudanças na espessura da moeda desde a Segunda Guerra Mundial. O Tesouro restaurou a espessura da moeda de um centavo à sua espessura pré-guerra na década de 1960, mas reduziu-a novamente na década de 1970.

Figura 8.5: Suavizadores da média condicional e da moda condicional do conjunto de dados de espessura de moedas de um centavo dos EUA.

Observe que, fora dos períodos de transição, as estimativas da média condicional e da moda condicional são bastante semelhantes. No entanto, a estimativa da moda condicional detecta e retrata essas mudanças abruptas. É interessante refletir sobre a percepção inicial comum de que as estimativas da média condicional parecem bastante razoáveis até serem comparadas com as estimativas da moda condicional.

Os traços da moda condicional se sobrepõem por vários anos em torno de 1957 e 1974. Essa sobreposição resulta do fato de existir uma condição de “limite” desconhecida nos dados e, portanto, a estimativa de densidade se comporta de maneira previsível. Seria interessante inserir kernels de limite interno neste problema de estimação (veja o Exercício 4).

Observe que a estimativa da moda condicional também detecta uma mudança rápida menor na espessura da moeda de um centavo por volta de 1985, mas não por um salto descontínuo. As estimativas kernel e da moda seguem direções diferentes em 1990.

Assim, a estimativa modal condicional, embora introduzida como um algoritmo de regressão alternativo resistente a outliers, tem potencial além dessa finalidade. O problema encontrado no conjunto de dados dos centavos é um exemplo do problema do “ponto de mudança”, onde ocorre um salto na curva de regressão (ver Siegmund 1988).

A exibição do traço modal condicional introduz algumas escolhas de design interessantes. As curvas acima foram calculadas fatia por fatia. As estimativas são mostradas com um símbolo “\(x\)” se a densidade for menor que 25% do valor máximo da densidade bivariada. Outras escolhas poderiam ser investigadas. Em particular, quando os dados contêm muitos outliers, uma consequência é o aparecimento de um traço modal adicional em cada um desses outliers. Os traços de outliers serão muito curtos, compostos por símbolos “\(x\)”, e servirão para identificar os outliers.

Os traços modais “centrais” não serão afetados por esses outliers; portanto, o traço modal é resistente a outliers. Finalmente, o cálculo propriamente dito das modas condicionais foi realizado através da combinação linear da estimativa ASH bivariada utilizando o kernel biweight. Tarter, Lock, and Mellin (1990) também desenvolveram um software para calcular modas condicionais. Uma ideia relacionada é discutida em Sager and Thisted (1982).


8.3.3 Regressão \(L_1\)


No modelo linear paramétrico geral, o critério \(L_1\) aplicado aos resíduos leva a uma estimativa particularmente resistente. Considere dados dispostos em uma linha reta, exceto por um outlier. Como o ajuste por mínimos quadrados passa pelo ponto \((\overline{x},\overline{y})\), fica claro que o ajuste será afetado pelo outlier.

Por outro lado, o ajuste \(L_1\) passará pelos \(n-1\) pontos da reta, ignorando o outlier. A razão é simples: qualquer deslocamento da reta acarreta um custo para cada um dos \(n-1\) pontos e um ganho apenas para o ponto discrepante. Esse argumento deixa claro que o ajuste \(L_1\) também será resistente a múltiplos outliers.

Wang (1990) investigou os aspectos teóricos e práticos do critério \(L_1\) em regressão não paramétrica. A ideia é essencialmente equivalente ao critério em (8.6), com a perda quadrática substituída pela perda absoluta. É bem conhecido que a solução para este problema é a solução de um problema de programação linear comum (Wagner 1959).

Se o \(i\)-ésimo resíduo for \(\epsilon_i = y_i-x_{(i)}\beta\) e \(\epsilon_i^+\) e \(\epsilon_i^-\) forem definidos como os valores absolutos das partes positiva e negativa de \(\epsilon_i\), respectivamente, então \[ \epsilon_i=\epsilon_i^+ -\epsilon_i^- \qquad \mbox{e} \qquad |\epsilon_i|=\epsilon_i^+ +\epsilon_i^-\cdot \] Assim, o problema de minimizar a soma dos resíduos absolutos torna-se \[ \min_{\epsilon_i^+,\epsilon_i^-} \sum_{i=1}^n \big( \epsilon_i^+ +\epsilon_i^-\big) \] restrito à \[ \pmb{y}=\pmb{X}\beta+\epsilon_i^+ -\epsilon_i^- \qquad \mbox{e} \qquad \epsilon_i^+,\epsilon_i^- \geq 0\cdot \]

Pacotes modernos de programação linear (PL) podem lidar com dezenas de milhares de pontos em tempo razoável. Para tornar o procedimento não paramétrico, o modelo \(L_1\) é aplicado localmente, com pesos de kernel nos valores absolutos dos critérios. O problema modificado permanece um problema de PL. Wang discute métodos eficientes para atualizar a série de programas lineares necessários para definir a curva de regressão não paramétrica \(L_1\).

Como exemplo, ajustes quadráticos locais foram calculados com os dados de acidentes de motocicleta usando os critérios \(L_1\) e \(L_2\) (ver Figura 8.6). Um kernel uniforme foi aplicado em um intervalo contendo os \(k\) vizinhos mais próximos.

Figura 8.6: (Esquerda) Ajustes quadráticos locais \(L_1\) e \(L_2\) aplicados aos dados da motocicleta, juntamente com um ajuste LOWESS com \(f = 0.25\). (Direita) Critérios de validação cruzada leave-one-out (LOO) normalizados, ou seja, o erro absoluto médio e o desvio padrão para os ajustes \(L_1\) e \(L_2\), respectivamente. Os valores brutos variam de \((16.2, 18.6)\) e \((18.5, 20.8)\), respectivamente. O melhor ajuste \(L_1\) ocorreu com 27 pontos em cada vizinhança local.

As curvas de validação cruzada leave-one-out são mostradas no quadro à direita da Figura 8.6, com uma modificação importante. Como um ajuste robusto produzirá grandes resíduos, usar a perda quadrática na função de validação cruzada LOO não será robusto por si só. Em vez disso, Wang and Scott (1994) sugerem usar a média dos resíduos absolutos para a função de validação cruzada leave-one-out (LOO-CV). Isso resultou em uma escolha baseada em dados de \(\widehat{k}=27\). A estimativa de mínimos quadrados locais foi calibrada em \(\widehat{k}=31\) usando a função LOO-CV usual; no entanto, a curva de risco não aumenta rapidamente para \(k > 31\). Uma curva LOWESS é incluída para referência, embora um parâmetro de suavização menor possa ser apropriado.

A característica especial que torna esta proposta interessante é a sua extensão imediata a múltiplas dimensões. Identificar pontos influentes torna-se problemático além de três dimensões (dois preditores) e esta abordagem parece bastante promissora, pois não requer interação (ou iteração) para eliminar os efeitos de outliers.

Um exemplo ilustrado na Figura 8.7, com uma superfície bivariada suave em uma malha \(61\times 61\) contaminada por 2.5% de ruído de grande variância \(N(0, 2.5^2)\) ou 97.5% de ruído \(N(0, 0.1^2)\), serve para ilustrar o poder da abordagem. A superfície verdadeira era uma porção de uma densidade normal bivariada padrão sobre \([-3, 3]^2\) multiplicada por \(2\pi\); portanto, a superfície varia de 0 a 1.

Figura 8.7: (Linha superior) Superfície de regressão que é uma normal bivariada 2π× em [−3, 3] × [−3, 3] em uma malha 61 × 61; estimativa não paramétrica quadrática completa local L1 com m x = 10 e m = 8; e estimativa não paramétrica quadrática completa local L2 com os mesmos parâmetros de suavização. (Linha inferior) Superfície contaminada com ruído da densidade da mistura 0,975 N(0, 0,12 ) + 0,025 N(0, 2,52 ), e superfícies residuais para as estimativas acima. Os gráficos de resíduos estão centrados em z = 0 e a escala vertical foi expandida por um fator de 3.

O modelo polinomial local era quadrático e um kernel uniforme foi introduzido no critério. O ruído é tão grande que a escala vertical da superfície dos dados brutos é truncada em \((-0.25, 1.25)\), sendo o intervalo vertical completo \((-5.17, 6.58)\). O critério robusto LOO-CV selecionou uma vizinhança retangular de \(21\times 17\). O ajuste \(L_1\) não é afetado pelos agrupamentos de outliers. O método dos mínimos quadrados locais apresenta mais ruído, mesmo neste exemplo com grande quantidade de dados. Se um modelo linear local for empregado (não mostrado), haverá um viés perceptível nos cantos, consulte Wang and Scott (1994) para mais detalhes e exemplos.


8.4 Regressão em várias dimensões


8.4.1 Suavização por kernel e WARPing


A extensão do estimador de Nadaraya-Watson para várias dimensões é direta. O estimador multivariado é novamente uma média local das respostas, com o kernel produto definindo o tamanho da vizinhança local e os pesos específicos das respostas.

Assim como no estimador de densidade kernel, a carga computacional pode ser bastante reduzida considerando-se um algoritmo de discretização semelhante ao ASH. Uma breve reflexão sugere que a discretização precisa ser feita apenas no espaço dos preditores, registrando-se somente a soma de todas as respostas em cada intervalo.

Para fins de notação, a soma das \(\nu_k\) respostas do \(k\)-ésimo intervalo com resposta média \(\overline{y}_k\) pode ser calculada por \(\nu_k\, \overline{y}_k\). No algoritmo de discretização de regressão (bivariada) (REG-BIN), a quantidade \(s\,y_k = \nu_k\, \overline{y}_k\) é tabelada, representando a soma das respostas do \(k\)-ésimo intervalo.

\[ \begin{array}{ll} & \mbox{REG-BIN}(x,y,n,a,b,nbin) \quad \mbox{Algoritmo:}\\[0.8em] & \delta=(b-a)/nbin \\[0.8em] & \mbox{for } k=1,nbin \; \{ \nu_k=0;sy_k=0 \,\} \\[0.8em] & \mbox{for } i=1,n \; \{ \\[0.8em] & \qquad k=(x_i-a)/\delta+1 \\[0.8em] & \qquad \mbox{if } (k\in [1,nbin]) \; \nu_k=\nu_k+1; st_k=sy_k+y_i \; \} \\[0.8em] & \mbox{return } (\{\nu_k\},\{sy_k\}) \end{array} \]

A versão de regressão univariada do ASH (REG-ASH) sobre a malha igualmente espaçada \(\{t_k\}\) correspondente ao estimador de Nadaraya-Watson (8.3) é calculada da seguinte forma: \[ \widehat{r}(x)=\dfrac{\displaystyle \sum_{i=1}^n K_h(x-x_i) \, y_i}{\displaystyle \sum_{i=1}^n K_h(x-x_i)}\approx \dfrac{\displaystyle \sum_{k=1}^{nbin} K_h(x-t_k)\, \nu_k \overline{y}_k}{\displaystyle \sum_{k=1}^{nbin} K_h(x-t_k) \,\nu_k}\cdot \]

Se a estimativa de regressão for calculada na mesma grade \(\{t_k\}\) que os dados agrupados, os pesos do kernel precisam ser calculados apenas uma vez. Quando um kernel de suporte finito é usado para calcular os pesos \(\omega_m(k)\), a estimativa de regressão em alguns intervalos pode ser igual a 0/0.

O algoritmo detecta essa condição e retorna “not-a-number” (NaN), “não é um número” traduzindo ao português. O ganho de velocidade é comparável ao observado com o ASH.

\[ \begin{array}{ll} & \mbox{REG-ASH}(m,\nu,sy,nbin,a,b,n,\omega_m) \quad \mbox{Algoritmo:}\\[0.8em] & \delta=(b-a)/nbin; h= m\delta \\[0.8em] & \mbox{for } k=1,nbin \; \{ f_k=0;r_k=0;t_k=a+(k-0.5)\delta \,\} \\[0.8em] & \mbox{for } k=1,nbin \; \{ \\[0.8em] & \qquad \mbox{if } (\nu_k=0) \mbox{ next } k \\[0.8em] & \qquad \mbox{for } i=\max(1,k-m+1),\min(nbin,k+m-1) \; \{ \\[0.8em] & \qquad \qquad f_i = f_i+\nu_k \omega_m(i-k)\\[0.8em] & \qquad \qquad r_i=r_i+sy_k\, \omega_m(i-k) \, \}\, \} \\[0.8em] & \mbox{for } k=1,nbin \, \{\mbox{ if } (f_k>0) r_k=r_k/f_k \, \mbox{ else } \, r_k=NaN \} \\[0.8em] & \mbox{return } (\pmb{x}=\{t_k\},\pmb{y}=\{r_k\}) \end{array} \]

Os métodos ASH e REG-ASH são casos especiais de uma abordagem geral para acelerar cálculos do tipo kernel. A ideia básica é arredondar os dados para uma malha fixa e realizar os cálculos com esses dados. A ideia de usar tais dados arredondados foi explorada por Scott (1981).

A abordagem WARPing é descrita por Härdle and Scott (1988). A abordagem da Transformada Rápida de Fourier (FFT) para regressão é descrita por Wand (1994).


8.4.2 Modelagem aditiva


Para mais do que algumas variáveis, além da crescente preocupação com a maldição da dimensionalidade e os efeitos de fronteira, o método kernel começa a perder sua interpretabilidade. Há uma necessidade de “restringir” os kernels multivariados de uma maneira que permita uma modelagem flexível, mas que mantenha a facilidade de interpretação da modelagem paramétrica multivariada. Os modelos de regressão aditiva atingem esse objetivo.

A superfície de regressão é modelada pela equação \[ r(\pmb{x})=r_0+\sum_{i=1}^d r_i(x_i)+\epsilon_i\cdot \]

A flexibilidade é alcançada ajustando suavizadores lineares unidimensionais para \(r_i\) em cada dimensão. O conjunto de superfícies multivariadas que podem ser geradas por este modelo não é muito complexo, mas é mais rico do que escolhas paramétricas como \(x_i\) ou \(x_i^2\) para \(r_i(x_i)\).

Uma flexibilidade adicional pode ser alcançada adicionando termos da forma \(r_{ij} (x_i, x_j)\). No entanto, a explosão combinatória ocorre rapidamente. O nível de cada função aditiva é arbitrário, embora uma constante \(r_0\) possa ser introduzida no modelo e cada função aditiva restringida a ter média 0.

O ajuste de tais modelos assume dois extremos. Por um lado, Wahba (1990) defendeu uma abordagem de função de penalidade que estima as \(d\) funções simultaneamente como splines. Um esquema iterativo mais simples foi proposto por Friedman and Stuetzle (1981), chamado de retroajuste. A ideia é começar com um conjunto inicial de estimativas de \(r_i\), possivelmente obtidas parametricamente, e iterar sobre cada uma das \(d\) variáveis, suavizando os resíduos não explicados pelos \(d-1\) preditores restantes.

Qualquer suavizador de diagrama de dispersão bivariado pode ser usado. Suavizadores de kernel simples podem ser usados, assim como splines de suavização. Friedman and Stuetzle (1981) usaram o super-suavizador no algoritmo de regressão por projeção (PPR), mais geral. Hastie and Tibshirani (1990) demonstram como o algoritmo de retroajuste imita o algoritmo de Gauss-Siedel para resolver equações lineares para certos suavizadores lineares. Em qualquer caso, o algoritmo geralmente converge rapidamente para uma vizinhança de uma solução razoável.

O algoritmo de retroajuste é facilmente implementado usando o estimador de regressão WARPing após a remoção da média de \(\{y_i\}\):

\[ \begin{array}{l} \mbox{Inicializar } r_j(x_j)=0 \\[0.8em] \mbox{Repita o processo em } j \mbox{ até a convergência } \, \{ \\[0.8em] \displaystyle \qquad \epsilon_j = y_j-\sum_{j\neq i} r_j(x_i) \\[0.8em] \displaystyle \qquad r_j = \mbox{alisador}(\epsilon_i) \, \} \\[0.8em] \end{array} \]

O modelo aditivo REG-ASH (RAM) foi aplicado a um exemplo de simulação de Wahba (1990) com \(\pmb{x}\in\mathbb{R}^4\), \(\epsilon_i\sim N(0,1)\), \(n= 100\) pontos de projeto amostrados uniformemente sobre o hipercubo \((0,1)^4\), e superfície aditiva verdadeira \[ \tag{8.23} r(\pmb{x})=10\times \sin(\pi x_2)+\exp(3\, x_3)+10^6\times x_4^{11}\times (1-x_4)^6+10^4 \times x_4^3 \times (1-x_4)^{10}, \] onde o valor verdadeiro da função aditiva \(r_1(x_1 ) = 0\).

Na Figura 8.8, a primeira linha exibe os gráficos de \((x_i ,\epsilon_i )\) para a primeira iteração juntamente com a suavização calculada; na segunda linha, a iteração final é exibida após cinco iterações do algoritmo. A variância residual inicial era de 57.2. Mesmo ao final da primeira iteração, a variação residual havia sido significativamente reduzida para 4.7 e era de 1.45 na solução final.

Figura 8.8: Iteração do modelo aditivo para dados simulados de (8.23). A função aditiva verdadeira é mostrada como uma linha pontilhada e a função aditiva estimada como uma linha contínua. A linha superior mostra o laço inicial e a linha inferior a iteração final (seis laços).

Um segundo exemplo mostrado na Figura 8.9 vem do conjunto completo de dados de fluxo de gás. O medidor foi testado em sete pressões diferentes com entre 24 e 51 fluxos em cada pressão. A primeira variável é o logaritmo da taxa de fluxo reescalado para (0,1). A segunda variável é o logaritmo da pressão (psia) também reescalado para (0,1). A variância residual inicial de 0.27 foi reduzida para 0.12. A relação com o logaritmo da pressão (psia) é quase perfeitamente linear, enquanto a relação com o logaritmo da taxa de fluxo é mais complexa. A precisão do medidor é afetada principalmente pela vazão, mas a pressão também exerce uma pequena influência.

Figura 8.9: Iteração do modelo aditivo para o conjunto de dados de fluxo de gás. As 214 medições foram realizadas em sete pressões diferentes na tubulação. Os ajustes aditivos para as iterações inicial e final são apresentados.


8.4.3 A maldição da dimensionalidade


Um modelo aditivo ainda é um estimador kernel, embora com kernels altamente restritos e sujeitos aos mesmos problemas assintóticos. O problema mais comum continua sendo o dos dados preditores com posto deficiente. Alguns resultados teóricos exigem que os dados se tornem densos sobre o hiper-retângulo que define o modelo aditivo. Essa suposição é claramente muito forte e, presumivelmente, necessária não na prática, mas para a demonstração de teoremas.

No entanto, o exemplo de Wahba na seção anterior pode ser modificado para ilustrar o efeito da maldição da dimensionalidade. Escolha \(\pmb{x}\sim U([0,1]^4)\), mas com correlações aos pares de 0.99. Na prática, isso é obtido aplicando-se a inversa da transformação de componentes principais a pontos amostrados uniformemente no hipercubo e, em seguida, reescalonando-os de volta para o hipercubo.

Os ajustes aditivos inicial e final são mostrados na Figura 8.10. As variâncias residuais inicial e final foram 6.6 e 2.1, respectivamente. Contudo, como os pontos de projeto estão concentrados ao longo da hiperdiagonal, as funções aditivas individuais não são identificáveis, uma vez que o espaço de projeto é singular. De fato, apenas uma função é necessária, defina \(x_2 = x_3 = x_4\) em (8.23). Em dimensões muito altas, tais singularidades são comuns e a maldição da dimensionalidade se aplica.

Figura 8.10: Iteração do modelo aditivo para dados simulados de (8.23). As correlações aos pares dos pontos de projeto são todas de 0.99. A verdadeira função aditiva é mostrada como uma linha pontilhada e a função aditiva estimada como uma linha contínua.

Um modelo aditivo fornecerá uma estimativa em toda a base do hiper-retângulo, mesmo onde não houver dados. Pode parecer simples verificar se existem dados próximos, mas isso nos leva de volta ao problema original de estimação de densidade e ao estimador de regressão kernel, para o qual o denominador pode ser usado para determinar a “precisão” da estimativa.

Finalmente, a interpretação das formas das funções \(r_i\) deve ser feita com cautela e somente após cuidadosa reflexão sobre as dificuldades de estimar a estrutura em modelos paramétricos. Ajustar o modelo aditivo permutando as variáveis e observar se a ordem afeta a forma dos ajustes é uma boa ferramenta de diagnóstico, assim como os métodos de reamostragem discutidos no Capítulo 9.

Como antes, se for possível projetar os dados originais em um subespaço conveniente, os algoritmos não paramétricos terão uma oportunidade aprimorada de realizar seu trabalho. As limitações dos modelos aditivos são discutidas por diversos autores na análise de Friedman (1991). Scott (1991) também levanta algumas ressalvas.

Componentes principais podem ser aplicados à matriz \(X\) antes da modelagem aditiva. Especificamente, os dados \((\pmb{X}, y)\) podem ser comprimidos para a forma \((\pmb{XP}, y)\), onde \(\pmb{P}\) é uma matriz de projeção. No entanto, isso ainda pode deixar um subespaço muito grande onde a resposta é constante, pois as componentes principais não incorporam nenhuma informação da resposta.

Uma nova ideia, a regressão inversa fatiada (SIR), proposta por Li (1991), incorpora a informação da resposta em uma abordagem projetada para dados normais, mas, assim como as componentes principais, possui aplicabilidade geral muito maior. O algoritmo SIR é uma extensão simples das ideias das componentes principais e componentes informativos. As respostas \(y_i\) são particionadas em \(k\) intervalos e a média do vetor \((x_i, y_i)\) de dimensão \((d+1)\) é calculada para cada intervalo, resultando no conjunto de \(k\) pontos. \[ (\overline{\pmb{x}}_j,\overline{y}_j), \qquad j=1,\cdots,k\cdot \]

A análise de componentes principais é aplicada à matriz de covariância dos \(k\) vetores \(\{\overline{\pmb{x}}_j\}\), e as direções com os maiores autovalores são retidas como preditores de \(y\).

As ideias básicas do SIR são ilustradas por um exemplo. Considere a superfície \[ r(x_1,x_2)=0.5625-0.1375 \, x_1 -0.2875 \, x_2 + 0.0125 \, x_1 \, x_2, \] sobre \(n = 100\) pontos de projeto amostrados de \(N(0, 0.4^2\pmb{I}_2 )\) com ruído \(\epsilon\sim N(0, 0.1^2 )\).

Na Figura 8.11, o contorno do plano que define a superfície real (e sua projeção) é exibido juntamente com os 100 pontos de dados. A escolha de \(k = 5\) intervalos foi feita de forma que cada intervalo tivesse 20 pontos. As seis marcas de escala no eixo \(y\) mostram as localizações dos intervalos.

Figura 8.11: Exemplo da técnica de redução de dimensionalidade SIR.

As médias condicionais para cada um dos cinco intervalos foram calculadas e são exibidas como losangos juntamente com uma linha que aponta para a projeção no espaço de projeto. A matriz de covariância dos cinco pontos no plano de projeto foi calculada e constatou-se que possui dois autovalores, dados por 0.365 e 0.034.

Assim, a direção SIR é o autovetor \((0.667, 0.745)^\top\) correspondente ao maior autovalor. A elipse correspondente à matriz de covariância é mostrada em torno das cinco médias. Observe que as respostas são essencialmente constantes na direção do outro autovetor. Para dados de dimensões mais elevadas que não seguem uma distribuição normal, essa ideia simples provou ser uma ferramenta poderosa para diagnóstico e exploração.

Modelos simples, como o algoritmo SIR, são atraentes devido à sua solução parcimoniosa. Um modelo simples semelhante é o modelo de árvore de regressão para a superfície de regressão que é constante por partes \[ \widehat{r}(\pmb{x})=\sum_{i=1}^p c_i \, \pmb{I}_{\{\pmb{x}\in \mathbb{N}\}}, \qquad \mbox{onde} \qquad \bigcup_{i=1}^p N_i=\mathbb{R}^d \quad \mbox{e} \quad N_i\cap N_j =\emptyset\cdot \] O chamado modelo de regressão CART tem tido muito sucesso na prática para problemas de classificação (ver Breiman et al. 1984).

A escolha do kernel em regressão é mais flexível do que em estimação de densidade. Em particular, o uso de kernels de ordem superior é mais atraente porque o problema da não negatividade desaparece. Deve-se ter cuidado nos limites para evitar anomalias na estimativa. Kernels de limite podem induzir artefatos com a mesma frequência com que os corrigem.

O uso direto de estimativas de kernel para regressão pode ou não ser a melhor escolha. Certamente, questões de resistência a outliers são de importância prática. Algoritmos de regressão de kernel ingênuos não oferecem proteção contra pontos influentes. O uso de regressão modal deve ser incentivado. Modelos aditivos fornecem a representação mais facilmente compreensível de superfícies de regressão multidimensionais, quando essas superfícies não são muito complexas. No entanto, a maldição da dimensionalidade sugere que as funções aditivas podem ser bastante imprecisas se os dados de projeto forem quase singulares.

Espera-se que o uso da norma \(L_1\) como uma medida robusta para ajuste polinomial local cresça drasticamente à medida que a disponibilidade de bons códigos de programação linear aumentar. No mínimo, o estimador \(L_1\) fornece uma excelente ferramenta de diagnóstico para dados com muitos valores discrepantes potenciais e pontos influentes.


8.5 Exercícios


  • 1- Considere o delineamento fixo com \(x_i = i/n\). Mostre que o erro quadrático médio (MSE) do estimador de Nadaraya-Watson neste caso é \[ \mbox{MSE}\big(\widehat{r}(x) \big)=\dfrac{\sigma_\epsilon^2 R(K)}{nh}+\dfrac{1}{4}h^4 \sigma_K^2 R(r'')\cdot \] Dica: Use a aproximação \(\displaystyle \sum_i r''(x_i)^2 \times (1/n)\approx \int_0^1 r''(x)^2 \mbox{d}x\).


  • 2- Verifique se o spline de suavização proposta na equação (8.21) satisfaz a equação de Euler-Lagrange, assumindo que \(\widehat{r}_\lambda^{(1)}\) e \(\widehat{r}_\lambda^{(2)}\) a satisfazem.


  • 3- Prove que o vetor de peso \(\pmb{\omega}_x\) na equação (8.18) é linear em \(x\).


  • 4- Construa uma modificação de limite interno para o conjunto de dados de moedas de um centavo exibido na Figura 8.5.


  • 5- Considere um modelo aditivo bivariado para uma função \(f (x, y)\), onde cada função aditiva é uma função constante por partes. Especificamente, \[ f_1(x)=\sum_i a_i \pmb{I}_{A_i}(x) \qquad \mbox{e} \qquad f_2(x)=\sum_j b_j \pmb{I}_{B_j}(x), \] onde \(\{A_i\}\) e \(\{B_j\}\) são intervalos ao longo dos eixos \(x\) e \(y\), respectivamente.

    1. Dados os vetores de dados \(\{(x_i, y_i, z_i), \, i = 1, n\}\), encontre as equações que correspondem às estimativas de mínimos quadrados de \(\{a_i\}\) e \(\{b_j\}\).

    2. Calcule e mostre graficamente esse estimador para o conjunto de dados de precisão do fluxo de gás, para diversas escolhas da malha.


8.6 Bibliografia


Allen, D. M. 1974. “The Relationship Between Variable Selection and Data Augmentation and a Method for Prediction.” Technometrics, no. 16: 125–27.
Becker, R. A., J. M. Chambers, and A. R. Wilks. 1988. The New s Language. Wadsworth & Brooks/Cole, Pacific Grove, CA.
Breiman, L., J. H. Friedman, A. Olshen, and C. J. Stone. 1984. CART: Classification and Regression Trees. Wadsworth, Belmont, CA.
Clark, R. M. 1975. “A Calibration Curve for Radiocarbon Dates.” Antiquity, no. 49: 251–66.
Cleveland, W. S. 1979. “Robust Locally Weighted Regression and Smoothing Scatterplots.” Journal of the American Statistical Association, no. 74: 829–36.
Eubank, R. L. 1988. Spline Smoothing and Nonparametric Regression. Marcel Dekker, New York.
Fan, J. Q. 1992. “Design-Adaptive Nonparametric Regression.” Journal of the American Statistical Association, no. 87: 998–1004.
Friedman, J. H. 1984. “A Variable Span Smoother.” Technical Report 5, Department of Statistics, Stanford University.
———. 1991. “Multivariate Adaptive Regression Splines (with Discussion).” Annals of Statistics, no. 19: 1–141.
Friedman, J. H., and W. Stuetzle. 1981. “Projection Pursuit Regression.” Journal of the American Statistical Association, no. 76: 817–23.
Gasser, T., and J. Engel. 1990. “The Choice of Weights in Kernel Regression Estimation.” Biometrika, no. 77: 377–81.
Gasser, T., and H. G. Müller. 1979. “Kernel Estimation of Regression Functions.” Springer-Verlag, Berlin.
Hampel, F. R. 1974. “The Influence Curve and Its Role in Robust Estimation.” Journal of the America Statistical Association, no. 69: 383–93.
Härdle, W. 1990. Applied Nonparametric Regression. Cambridge University Press, Cambridge, England.
Härdle, W., P. Hall, and J. S. Marron. 1988. “How Far Are Automatically Chosen Regression Smoothing Parameters from Their Optimum? (With Discussion).” Journal of the American Statistical Association, no. 83: 86–101.
Härdle, W., and D. W. Scott. 1988. “Smoothing in Low and High Dimensions by Weighted Averaging Using Rounded Points.” Computational Statistics, no. 7: 97–128.
Hastie, T. J., and R. J. Tibshirani. 1990. Generalized Additive Models. Chapman; Hall, London.
Huber, P. J. 1964. “Robust Estimation of a Location Parameter.” Annals of the Mathematical Statistics, no. 33: 73–101.
Li, K.-C. 1991. “Sliced Inverse Regression for Dimension Reduction (with Discussion).” Journal of the American Statistical Association, no. 86: 316–42.
Mack, Y. P., and H. G. Müller. 1989. “Convolution Type Estimators for Nonparametric Regression.” Statistical and Probility Letters, no. 7: 229–39.
Müller, H. G. 1988. Nonparametric Regression Analysis of Longitudinal Data. Springer-Verlag, Berlin.
Nadaraya, E. A. 1964. “On Estimating Regression.” Theory of Probability and Its Applications, no. 15: 134–37.
Reinsch, C. H. 1967. “Smoothing by Spline Functions.” Numerische Mathematik, no. 10: 177–83.
Rice, J. A. 1984. “Bandwidth Choice for Nonparametric Kernel Regression.” Annals of Statistics, no. 12: 1215–30.
Rosenblatt, M. 1969. “Conditional Probability Density and Regression Estimates.” Edited by P. R. Krishnaiah. Academic Press, New York.
Rousseeuw, P. J., and A. M. Leroy. 1987. Robust Regression and Outlier Detection. John Wiley, New York.
Sager, T. W., and R. A. Thisted. 1982. “Maximum Likelihood Estimation of Isotonic Modal Regression.” Annals of Statistics, no. 10: 690–707.
Schoenberg, I. 1964. “On Interpolation by Spline Functions and Its Minimum Properties.” International Series of Numerical Mathematics, no. 5: 109–29.
Schucany, W. R. 1995. “Adaptive Bandwidth Choice for Kernel Regression.” Journal of the American Statistical Association, no. 90: 535–40.
Scott, D. W. 1981. “Using Computer-Binned Data for Density Estimation.” Springer-Verlag, New York.
———. 1991. “On Estimation and Visualization of Higher Dimensional Surfaces.” Springer-Verlag, New York.
———. 2015. Multivariate Density Estimation. John Wiley & Sons. Inc.
Shibata, R. 1981. “An Optimal Selection of Regression Variables.” Biometrika, no. 68: 45–54.
Siegmund, D. 1988. “Confidence Sets in Change-Point Problems.” International Statistical Review, no. 56: 31–48.
Silverman, B. W. 1984. “Spline Smoothing: The Equivalent Variable Kernel Method.” Annals of Statistics, no. 12: 898–916.
Simonoff, J. S. 1996. Smoothing Methods in Statistics. Springer, New York.
Staniswalis, J. G. 1989. “Local Bandwidth Selection for Kernel Estimates.” Journal of the American Statistical Association, no. 84: 276–83.
Stone, C. J. 1977. “Consistent Nonparametric Regression (with Discussion).” Annals of Statistics, no. 5: 595–645.
———. 1985. “Additive Regression and Other Nonparametric Models.” Annals of Statistics, no. 13: 689–705.
Stone, M. 1974. “Cross-Validatory Choice and Assessment of Statistical Predictions (with Discussion).” Journal of the Royal Statistical Society, B, no. 36: 111–47.
Stuart, A., and J. K. Ord. 1987. Kendall’s Advanced Theory of Statistics, Volume 1. Oxford University Press, New York.
Tarter, M. E., M. D. Lock, and C. C. Mellin. 1990. Curves: Background and Program Description. Precision Data Group, Berkeley, CA.
Wagner, H. M. 1959. “Linear Programming Techniques for Regression Analysis.” Journal of the American Statistical Association, no. 54: 206–12.
Wahba, G. 1990. Spline Models for Observational Data. SIAM, Philadephia.
Wahba, G., and S. Wold. 1975. “A Completely Automatic French Curve: Fitting Spline Functions by Cross-Validation?” Communications in Statistics, no. 4: 1–17.
Wand, M. P. 1994. “Fast Computation of Multivariate Kernel Estimators.” Journal of Computational and Graphical Statistics, no. 3: 433–45.
Wang, F. T. 1990. “A New Method for Robust Nonparametric Regression.” PhD thesis, Department of Statistics, Rice University.
Wang, F. T., and D. W. Scott. 1994. “The \(L_1\) Method for Robust Nonparametric Regression.” Journal of the American Statistical Association, no. 89: 65–76.
Watson, G. S. 1964. “Smooth Regression Analysis.” Sankhyā Ser. A, no. 26: 359–72.