Como incluir um termo de interação no GAM?

24

O código a seguir avalia a semelhança entre duas séries temporais:

set.seed(10)
RandData <- rnorm(8760*2)
America <- rep(c('NewYork','Miami'),each=8760)

Date = seq(from=as.POSIXct("1991-01-01 00:00"), 
           to=as.POSIXct("1991-12-31 23:00"), length=8760)

DatNew <- data.frame(Loc = America,
                     Doy = as.numeric(format(Date,format = "%j")),
                     Tod = as.numeric(format(Date,format = "%H")),
                     Temp = RandData,
                     DecTime = rep(seq(1, length(RandData)/2) / (length(RandData)/2),
                                   2))
require(mgcv)
mod1 <- gam(Temp ~ Loc + s(Doy) + s(Doy,by = Loc) +
  s(Tod) + s(Tod,by = Loc),data = DatNew, method = "ML")

Aqui, gamé usado para avaliar como a temperatura em Nova York e Miami varia da temperatura média (de ambos os locais) em diferentes momentos do dia. O problema que tenho agora é que preciso incluir um termo de interação que mostre como a temperatura de cada local varia ao longo do dia em diferentes dias do ano. Eventualmente, espero exibir todas essas informações em um gráfico (para cada local). Então, para Miami, espero ter um gráfico que mostre como a temperatura varia em relação à média em diferentes horários do dia e em diferentes épocas do ano (gráfico 3D?)

KatyB
fonte
2
Você pode encontrar a resposta para esta pergunta stats.stackexchange.com/questions/18937/… relevante.
21412 jbowman

Respostas:

18

O "a" em "gam" significa "aditivo", o que significa que não há interações. Portanto, se você ajustar as interações, não estará mais adaptando um modelo de gam.

Dito isto, existem maneiras de obter alguma interação, como termos dentro dos termos aditivos em um gam, você já está usando um deles usando o byargumento para s. Você pode tentar estender isso para que o argumento byseja uma matriz com uma função (sin, cos) de doy ou tod. Você também pode ajustar splines de suavização em um modelo linear regular que permita interações (isso não fornece o ajuste que o gam faz, mas ainda pode ser útil).

Você também pode considerar a regressão de busca de projeção como outra ferramenta de ajuste. Modelos Loess ou mais paramétricos (com pecado e / ou cos) também podem ser úteis.

Parte da decisão sobre quais ferramentas usar é qual pergunta você está tentando responder. Você está apenas tentando encontrar um modelo para prever datas e horários futuros? você está tentando testar se determinados preditores são significativos no modelo? você está tentando entender a forma do relacionamento entre um preditor e o resultado? Algo mais?

Greg Snow
fonte
3
x1 1,x2
y=f1 1(x1 1)+f2(x2)+f3(x1 1x2)+ε
gam
y=f1 1(x1 1)+f2(x2)+f3(x1 1)x2+f4(x2)x1 1+ε
by
11
@ Macro, isso depende se você deseja dividir os cabelos ou não (bem, tecnicamente, o que você escreveu talvez fosse "am's", já que você não está realmente usando a parte g).
Greg Neve
2
@ Macro, tenho certeza que mais foi dito sobre o assunto, mas veja a página 4 dos GAMs do artigo deste Venable, Exegeses on Linear Models . Não está claro exatamente como os efeitos principais não aditivos e os efeitos de interação são identificados simultaneamente.
Andy W
muito obrigado por seus comentários. Meu principal objetivo aqui é não prever valores futuros, mas simplesmente ver como cada série temporal se altera da média de ambas as séries. Por exemplo, o resultado do mod1 mostra como a série temporal em cada local varia da média primeiro na escala de tempo do dia do ano (Doy) e depois na hora do dia (Tod). A partir disso, gostaria de ver como cada série varia em função de Doy e Tod, e espero que as séries temporais diferam bastante durante o período de verão.
24512 KatyB
2
@AndyW É importante observar que os GAM percorreram um longo caminho desde que a Venable os comentou - desenvolvimentos na última década com os splines penalizados sensu Simon Wood (conforme implementado no mgcv ) e os tratamentos totalmente bayesianos abordam questões como seleção de suavidade, interações e como ajustá-los (produtos tensores de bases marginais sendo uma abordagem) na estrutura do modelo aditivo. Estou razoavelmente confiante com as objeções de Venable & Cox aos GAMs, conforme descrito na exegeses foram em grande parte abordadas por esses desenvolvimentos recentes na teoria do GAM.
Reintegrar Monica - G. Simpson
25

Para duas variáveis ​​contínuas, você pode fazer o que quiser (seja uma interação ou não, deixarei outros discutirem conforme os comentários da resposta do @ Greg) usando:

mod1 <- gam(Temp ~ Loc + s(Doy, bs = "cc", k = 5) + 
                         s(Doy, bs = "cc", by = Loc, k = 5, m = 1) + 
                         s(Tod, bs = "cc", k = 5) + 
                         s(Tod, bs = "cc", by = Loc, k = 5, m = 1) +
                         te(Tod, Doy, by = Loc, bs = rep("cc",2)),
            data = DatNew, method = "ML")

O modelo mais simples deve ser aninhado no modelo mais complexo acima. Esse modelo mais simples é:

mod0 <- gam(Temp ~ Loc + s(Doy, bs = "cc", k = 5) + 
                         s(Doy, bs = "cc", by = Loc, k = 5, m = 1) + 
                         s(Tod, bs = "cc", k = 5) + 
                         s(Tod, bs = "cc", by = Loc, k = 5, m = 1),
            data = DatNew, method = "ML")

Observe duas coisas aqui:

  1. O tipo de base para cada suavidade é indicado. Nesse caso, esperamos que não haja descontinuidades no Temp entre 23:59 e 00:00 para Todnem entre Doy == 1e Doy == 365.25. Portanto, splines cúbicos cíclicos são apropriados, indicados aqui via bs = "cc".
  2. A dimensão base é declarada explicitamente ( k = 5). Isso corresponde à dimensão base padrão para cada suavização em um te()termo.

Juntos, esses recursos garantem que o modelo mais simples seja realmente aninhado no modelo mais complexo.

Para mais informações veja ?gam.modelsem mgcv .

Restabelecer Monica - G. Simpson
fonte
Relacionado ao seu segundo ponto - além da especificação de k, deve-se também fixar o número de nós (por exemplo fx=TRUE). Caso contrário, o modelo resultante mostra uma variação edfpara cada termo.
Marc na caixa
Eu devo atualizar esta resposta um pouco, considerando algumas novas funcionalidades no pacote mgcv em termos de splines para bases marginais. Dito isto, não concordo que você precise fixar os graus de liberdade para o spline. A chave é garantir que as bases dos modelos sejam aninhadas adequadamente. Então, diferenças entre modelos são possíveis, definindo alguns dos coeficientes para as funções básicas como zero, como aconteceria em um modelo linear com termos não spline.
Reinstate Monica - G. Simpson
3
Esperando que alguém ainda esteja assistindo este tópico e possa responder. Nesses modelos, como é que você deve especificar ambos s(Doy...)e s(Doy, by =Loc...)? Eu pensei que o primeiro seria aninhado no último e, portanto, não é necessário especificar?
ego_
3
Não, a primeira suavização é a função global e a by smooth representa diferenças específicas do site entre ela e a suavização global. O by smooths realmente precisa ser m = 1adicionado a eles para aplicar a penalidade na primeira derivada pela diferença de smooths.
Reponha Monica - G. Simpson
2
@ JoshuaRosenberg Se você quer dizer te(), depende do que você está incluindo no produto tensorial? As interações descritas aqui são interações fator-suave, mas te()implicariam duas ou mais variáveis ​​contínuas. Se você deseja termos globais e desvios específicos de um assunto, sim, te(DoY, Year, by = Loc, m = 1)pode ser usado em conjunto te(DoY, Year), embora existam outras maneiras de conseguir coisas semelhantes usando a interação fator-suave de te()efeito aleatório e termos contendo um spline de efeito aleatório.
Reintegrar Monica - G. Simpson