Intuição por trás de interações com produtos tensores em GAMs (pacote MGCV em R)

30

Modelos aditivos generalizados são aqueles em que por exemplo. as funções são suaves e devem ser estimadas. Geralmente por splines penalizados. MGCV é um pacote em R que faz isso, e o autor (Simon Wood) escreve um livro sobre seu pacote com exemplos de R. Ruppert et al. (2003) escrevem um livro muito mais acessível sobre versões mais simples da mesma coisa.

y=α+f1(x1)+f2(x2)+ei

Minha pergunta é sobre interações nesse tipo de modelo. E se eu quiser fazer algo parecido com o seguinte: se estivéssemos em terra OLS (onde é apenas um beta) , Não teria problema em interpretar \ hat {f} _3 . Se estimarmos por splines penalizados, também não tenho problemas com a interpretação no contexto aditivo.

y=α+f1(x1)+f2(x2)+f3(x1×x2)+ei
f 3ff^3

Mas o pacote MGCV no GAM tem essas coisas chamadas "suavização do produto tensorial". Eu pesquiso no Google "produto tensor" e meus olhos imediatamente olham para tentar ler as explicações que encontro. Ou eu não sou inteligente o suficiente ou a matemática não é explicada muito bem, ou ambas.

Em vez de codificar

normal = gam(y~s(x1)+s(x2)+s(x1*x2))

um produto tensorial faria a mesma coisa (?)

what = gam(y~te(x1,x2))

quando eu faço

plot(what)

ou

vis.gam(what)

Eu recebo uma saída muito legal. Mas não tenho idéia do que está acontecendo dentro da caixa preta te(), nem como interpretar a saída legal acima mencionada. Na outra noite, tive um pesadelo que estava dando em um seminário. Mostrei a todos um gráfico interessante, eles me perguntaram o que significava e eu não sabia. Então eu descobri que não estava com roupas.

Alguém poderia ajudar a mim e à posteridade, dando um pouco de mecânica e intuição sobre o que está acontecendo por baixo do capô aqui? Idealmente, falando um pouco sobre a diferença entre o caso normal de interação aditiva e o caso tensor? Pontos de bônus por dizer tudo em inglês simples antes de passar para a matemática.

generic_user
fonte
exemplo simples, extraído do livro do autor do pacote: dados da biblioteca (mgcv) (árvores) ct5 <- gam (Volume ~ te (altura, circunferência, k = 5), família = gama (link = log), dados = árvores) ct5 vis.gam (CT5) trama (CT5, too.far = 0,15)
generic_user

Respostas:

30

Vou tentar responder a isso em três etapas: primeiro, vamos identificar exatamente o que queremos dizer com suavidade univariada. A seguir, descreveremos uma suavização multivariada (especificamente, uma suavidade de duas variáveis). Finalmente, farei a minha melhor tentativa de descrever um produto tensorial suave.

1) Univariado suave

Digamos que temos alguns dados de resposta que conjecturamos é uma função desconhecida de uma variável preditora mais algum erro . O modelo seria:f x εyfxε

y=f(x)+ε

Agora, para se ajustar a esse modelo, precisamos identificar a forma funcional de . A maneira como fazemos isso é identificando funções básicas, que são superpostas para representar a função na sua totalidade. Um exemplo muito simples é uma regressão linear, na qual as funções são apenas e , a interceptação. Aplicando a expansão da base, temosf β 2 x β 1ffβ2xβ1

y=β1+β2x+ε

Em forma de matriz, teríamos:

Y=Xβ+ε

Onde é um vetor de coluna n por 1, é uma matriz de modelo n por 2, é um vetor de coluna 2 por 1 dos coeficientes do modelo e é um vetor de erros n por 1 da coluna . tem duas colunas porque existem dois termos em nossa expansão de base: o termo linear e o intercepto.X β ε XYXβεX

O mesmo princípio se aplica à expansão da base no MGCV, embora as funções da base sejam muito mais sofisticadas. Especificamente, as funções de base individuais não precisam ser definidas sobre o domínio completo da variável independente . Esse é geralmente o caso ao usar bases baseadas em nós (consulte "exemplo baseado em nós"x) O modelo é então representado como a soma das funções básicas, cada uma das quais é avaliada em todos os valores da variável independente. No entanto, como mencionei, algumas dessas funções básicas assumem um valor zero fora de um determinado intervalo e, portanto, não contribuem para a expansão da base fora desse intervalo. Como exemplo, considere uma base de spline cúbica na qual cada função de base é simétrica em relação a um valor diferente (nó) da variável independente - em outras palavras, cada função de base tem a mesma aparência, mas é apenas deslocada ao longo do eixo da variável independente (isso é uma simplificação excessiva, pois qualquer base prática também incluirá um termo de interceptação e linear, mas espero que você entenda).

Para ser explícito, uma expansão básica da dimensão pode parecer com:i2

y=β1+β2x+β3f1(x)+β4f2(x)+...+βifi2(x)+ε

onde cada função é, talvez, uma função cúbica da variável independente .xfx

A equação da matriz ainda pode ser usada para representar nosso modelo. A única diferença é que agora é uma matriz n-por-i; isto é, possui uma coluna para cada termo na expansão da base (incluindo o termo de interceptação e linear). Como o processo de expansão da base nos permitiu representar o modelo na forma de uma equação matricial, podemos usar mínimos quadrados lineares para ajustar o modelo e encontrar os coeficientes . X βY=Xβ+εXβ

Este é um exemplo de regressão não compensada, e uma das principais forças do MGCV é sua estimativa de suavidade por meio de uma matriz de penalidade e parâmetro de suavização. Em outras palavras, em vez de:

β=(XTX)1XTY

temos:

β=(XTX+λS)1XTY

onde é uma matriz de penalidade quadrática por- e é um parâmetro de suavização escalar. Não vou entrar na especificação da matriz de penalidade aqui, mas seria suficiente dizer que, para qualquer base de expansão de alguma variável independente e definição de uma penalidade quadrática de "wiggliness" (por exemplo, uma penalidade de segunda derivada), uma pode calcular a matriz pena .i i λ SSiiλS

O MGCV pode usar vários meios de estimar o parâmetro de suavização ideal . Não vou entrar nesse assunto, já que meu objetivo aqui era fornecer uma ampla visão geral de como uma suavidade univariada é construída, o que acredito ter feito.λ

2) Multivariada suave

A explicação acima pode ser generalizada para várias dimensões. Vamos voltar ao nosso modelo que fornece a resposta como uma função dos preditores e . A restrição a duas variáveis ​​independentes evitará desordenar a explicação com notação arcana. O modelo é então:f x zyfxz

y=f(x,z)+ε

Agora, deve ser intuitivamente óbvio que vamos representar com uma expansão de base (isto é, uma superposição de funções de base), exatamente como fizemos no caso univariado de acima. Também deve ser óbvio que pelo menos uma, e quase certamente muitas mais, dessas funções básicas devem ser funções de e (se esse não fosse o caso, então implicitamente seria separável de modo que ). Uma ilustração visual de uma base de spline multidimensional pode ser encontrada aqui . Uma expansão de base bidimensional completa da dimensão pode ser algo como:f ( x ) x z f f ( x , z ) = f x ( x ) + f z ( z ) i - 3f(x,z)f(x)xzff(x,z)=fx(x)+fz(z)i3

y=β1+β2x+β3z+β4f1(x,z)+...+βifi3(x,z)+ε

Eu acho que é bastante claro que ainda podemos representar isso em forma de matriz com:

Y=Xβ+ε

simplesmente avaliando cada função básica em cada combinação única de e . A solução ainda é:zxz

β=(XTX)1XTY

O cálculo da matriz de penalidade da segunda derivada é praticamente o mesmo que no caso univariado, exceto que, em vez de integrar a segunda derivada de cada função de base em relação a uma única variável, integramos a soma de todas as segundas derivadas (incluindo parciais) com relação a para todas as variáveis ​​independentes. Os detalhes do exposto acima não são especialmente importantes: o ponto é que ainda podemos construir a matriz de penalidade e usar o mesmo método para obter o valor ideal do parâmetro de suavização e, dado esse parâmetro de suavização, o vetor de coeficientes ainda é:λSλ

β=(XTX+λS)1XTY

Agora, esse liso bidimensional tem uma penalidade isotrópica : isso significa que um único valor de se aplica nas duas direções. Isso funciona bem quando e estão aproximadamente na mesma escala, como um aplicativo espacial. Mas e se substituirmos a variável espacial pela variável temporal ? As unidades de podem ser muito maiores ou menores que as unidades de , e isso pode prejudicar a integração de nossas segundas derivadas, porque algumas dessas derivadas contribuirão desproporcionalmente para a integração geral (por exemplo, se medirmos em nanossegundos e x z z t t x t x t x xλxzzttxtxem anos-luz, a integral da segunda derivada em relação a pode ser muito maior que a integral da segunda derivada em relação a e, portanto, a "ondulação" ao longo da direção pode ser amplamente não-penalizada). O slide 15 da "caixa de ferramentas suave" que vinculei tem mais detalhes sobre esse tópico.txx

Vale a pena notar que não decompusemos as funções de base em bases marginais de e . A implicação aqui é que suavizações multivariadas devem ser construídas a partir de bases que suportam múltiplas variáveis. Os suportes de produto tensor suportam a construção de bases multivariadas a partir de bases marginais univariadas, como explico abaixo.zxz

3) O produto tensor suaviza

O produto tensor suaviza o problema das respostas de modelagem a interações de várias entradas com unidades diferentes. Vamos supor que temos uma resposta que é uma função da variável espacial e variável temporal . Nosso modelo é então:f x tyfxt

y=f(x,t)+ε

O que gostaríamos de fazer é construir uma base bidimensional para as variáveis e . Isso será muito mais fácil se pudermos representar como:t fxtf

f(x,t)=fx(x)ft(t)

No sentido algébrico / analítico, isso não é necessariamente possível. Mas lembre-se, estamos discretizando os domínios de e (imagine uma "tridimensional" bidimensional definida pela localização dos nós nos eixos e ), de modo que a função "verdadeira" seja representada pela superposição de funções de base . Assim como assumimos que uma função univariada muito complexa pode ser aproximada por uma função cúbica simples em um intervalo específico de seu domínio, podemos assumir que a função não separável possa ser aproximada pelo produto de funções mais simples extxtff(x,t)fx(x)ft(t) em um intervalo - desde que nossa escolha de dimensões básicas torne esses intervalos suficientemente pequenos!

A expansão de nossa base, dada uma base dimensional em e dimensional em , ficaria assim:ixjt

y=β1+β2x+β3fx1(x)+β4fx2(x)+...+βifx(i3)(x)+βi+1t+βi+2tx+βi+3tfx1(x)+βi+4tfx2(x)+...+β2itfx(i3)(x)+β2i+1ft1(t)+β2i+2ft1(t)x+β2i+3ft1(t)fx1(x)+βi+4ft1(t)fx2(x)+...+β2ift1(t)fx(i3)(x)++βijft(j3)(t)fx(i3)(x)+ε

O que pode ser interpretado como um produto tensorial. Imagine que avaliamos cada função base em e , construindo assim as matrizes do modelo n-por-i e n-por-j e , respectivamente. Poderíamos então calcular o -by- tensor produto destas duas matrizes modelo e reorganizar em colunas, de tal modo que cada coluna representada uma combinação única . Lembre-se de que as matrizes do modelo marginal tinham colunas e , respectivamente. Esses valores correspondem às suas respectivas dimensões base. Nossa nova base de duas variáveis ​​deve ter a dimensãoxtXTn2ij XTijijije, portanto, o mesmo número de colunas em sua matriz de modelo.

NOTA: Gostaria de salientar que, uma vez que construímos explicitamente as funções da base do produto tensorial, utilizando produtos de funções marginais, as bases do produto tensorial podem ser construídas a partir de bases marginais de qualquer tipo. Eles não precisam suportar mais de uma variável, ao contrário do multivariado suave discutido acima.

Na realidade, esse processo resulta em uma expansão geral da base da dimensão porque a multiplicação completa inclui a multiplicação de cada função de base pela interceptação x ( ), assim como a multiplicação de todas função base pela interceptação t ( ), mas devemos adicionar a interceptação novamente por si própria (então adicionamos 1). Isso é conhecido como aplicar uma restrição de identificabilidade.ijij+1tβx1jxβt1i

Então, podemos representar isso como:

y=β1+β2x+β3t+β4f1(x,t)+β5f2(x,t)+...+βijij+1fijij2(x,t)+ε

Onde cada uma das funções multivariadas é o produto de um par de funções marginais e . Novamente, é bastante claro que, tendo construído essa base, ainda podemos representar isso com a equação da matriz:fxt

Y=Xβ+ε

Qual (ainda) tem a solução:

β=(XTX)1XTY

Onde a matriz do modelo possui colunas . Quanto às matrizes de penalidade e , elas são construídas separadamente para cada variável independente da seguinte maneira:Xijij+1JxJt

Jx=βTIjSxβ

e,

Jt=βTStIiβ

Isso permite uma penalidade anisotrópica geral (diferente em cada direção) (Nota: as penalidades na segunda derivada de são somadas a cada nó no eixo vice-versa). Os parâmetros de suavização e agora podem ser estimados da mesma maneira que o único parâmetro de suavização foi para os suavizados univariados e multivariados. O resultado é que a forma geral de um produto tensorial suave é invariável ao redimensionamento de suas variáveis ​​independentes.xtλxλt

Eu recomendo a leitura de todas as vinhetas no site da MGCV, bem como " Modelos aditivos generalizados: e introdução ao R. ". Viva Simon Wood.

Josh
fonte
Boa resposta. Desde então, aprendi muito mais do que sabia há três anos. Mas não tenho certeza se teria entendido há 3 anos o que você escreveu hoje. Ou talvez eu tivesse. Penso que o ponto de partida é pensar em uma expansão de base em muitas dimensões como uma "rede" no espaço variável. Suponho que os tensores possam ser descritos como uma rede com padrões retangulares ... E talvez diferentes forças de "cisalhamento" puxando de cada direção.
generic_user 18/09/15
Em outra nota, eu o alertaria contra pensar no produto tensorial como representando algo espacial. Isso ocorre porque o produto tensorial real das funções marginais e incluirá toneladas de zeros que representam a avaliação das funções básicas fora do intervalo definido. O produto tensor real normalmente será muito escasso. xt
Josh
11
Obrigado por este ótimo resumo! Apenas uma observação: a equação após "Nossa expansão da base" não está completamente correta. Ele fornece as funções de base corretas, mas fornece uma parametrização onde os parâmetros correspondentes são da forma do produto ( ). βxiβtj
Javauh
11
@ Josh Ok, eu tentei. Não é fácil corrigi-lo e compreendê-lo ao mesmo tempo (e seguir a notação de outra pessoa). A propósito, o link para smooth-toolbox.pdf parece estar quebrado.
Javauh
11
Parece bom. Aparentemente, sua edição foi rejeitada, mas eu a substituí e a aprovei. Quando comecei a escrever essa resposta, não percebi o quão confusas as expansões seriam. Eu provavelmente deveria voltar e reescrevê-lo com a notação pi (produto) um dia desses.
Josh