Diferentes maneiras de modelar interações entre preditores contínuos e categóricos no GAM

A pergunta a seguir se baseia na discussão encontrada nesta página . Dada uma variável de resposta y, uma variável explicativa contínua xe um fator fac, é possível definir um Modelo Aditivo Geral (GAM) com uma interação entre xe facusando o argumento by=. De acordo com o arquivo de ajuda ?gam.models no pacote R mgcv, isso pode ser feito da seguinte maneira:

gam1 <- gam(y ~ fac +s(x, by = fac), ...)

O @GavinSimpson aqui sugere uma abordagem diferente:

gam2 <- gam(y ~ fac +s(x) +s(x, by = fac, m=1), ...)

Eu tenho brincado com um terceiro modelo:

gam3 <- gam(y ~ s(x, by = fac), ...)

Minhas principais perguntas são: alguns desses modelos estão errados ou são simplesmente diferentes? Neste último caso, quais são suas diferenças? Baseado no exemplo que discutirei abaixo, acho que consegui entender algumas das diferenças, mas ainda estou perdendo alguma coisa.

Como exemplo, vou usar um conjunto de dados com espectros de cores para flores de duas espécies diferentes de plantas medidas em locais diferentes.

rm(list=ls())
# install.packages("RCurl")
library(RCurl) # allows accessing data from URL
df <- read.delim(text=getURL("https://raw.githubusercontent.com/marcoplebani85/datasets/master/flower_color_spectra.txt"))
library(mgcv)

Para maior clareza, cada linha da figura acima representa o espectro médio de cores previsto para cada local com um GAM separado, com density~s(wl)base em amostras de ~ 10 flores. As áreas cinzas representam 95% de IC para cada GAM.

Meu objetivo final é modelar o efeito (potencialmente interativo) Taxone o comprimento wlde onda na refletância (referidos densityno código e no conjunto de dados), enquanto consideramos Localityum efeito aleatório em um GAM de efeito misto. No momento, não adicionarei a parte de efeito misto ao meu prato, que já está cheio o suficiente com a tentativa de entender como modelar interações.

Começarei com o mais simples dos três GAMs interativos:

gam.interaction0 <- gam(density ~ s(wl, by = Taxon), data = df) 
# common intercept, different slopes
plot(gam.interaction0, pages=1)

summary(gam.interaction0)

Produz:

Family: gaussian 
Link function: identity 

Formula:
density ~ s(wl, by = Taxon)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  28.3490     0.1693   167.4   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                      edf Ref.df     F p-value    
s(wl):TaxonSpeciesA 8.938  8.999 884.3  <2e-16 ***
s(wl):TaxonSpeciesB 8.838  8.992 325.5  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.523   Deviance explained = 52.4%
GCV = 284.96  Scale est. = 284.42    n = 9918

A parte paramétrica é a mesma para as duas espécies, mas splines diferentes são ajustados para cada espécie. É um pouco confuso ter uma parte paramétrica no resumo dos GAMs, que não são paramétricos. @IsabellaGhement explica:

Se você observar as plotagens dos efeitos suaves estimados correspondentes ao seu primeiro modelo, notará que eles estão centralizados em torno de zero. Portanto, você precisa 'mudar' essas suavizações (se a interceptação estimada for positiva) ou para baixo (se a interceptação estimada for negativa) para obter as funções suaves que você pensava estar estimando. Em outras palavras, você precisa adicionar a interceptação estimada aos suaves para obter o que realmente deseja. Para o seu primeiro modelo, o 'turno' é assumido como sendo o mesmo para os dois suaves.

Se movendo:

gam.interaction1 <- gam(density ~ Taxon +s(wl, by = Taxon, m=1), data = df)
plot(gam.interaction1,pages=1)

summary(gam.interaction1)

Dá:

Family: gaussian 
Link function: identity 

Formula:
density ~ Taxon + s(wl, by = Taxon, m = 1)

Parametric coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    40.3132     0.1482   272.0   <2e-16 ***
TaxonSpeciesB -26.0221     0.2186  -119.1   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                      edf Ref.df    F p-value    
s(wl):TaxonSpeciesA 7.978      8 2390  <2e-16 ***
s(wl):TaxonSpeciesB 7.965      8  879  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.803   Deviance explained = 80.3%
GCV = 117.89  Scale est. = 117.68    n = 9918

Agora, cada espécie também tem sua própria estimativa paramétrica.

O próximo modelo é o que tenho problemas para entender:

gam.interaction2 <- gam(density ~ Taxon + s(wl) + s(wl, by = Taxon,  m=1), data = df)
plot(gam.interaction2, pages=1)

Não tenho uma ideia clara do que esses gráficos representam.

summary(gam.interaction2)

Dá:

Family: gaussian 
Link function: identity 

Formula:
density ~ Taxon + s(wl) + s(wl, by = Taxon, m = 1)

Parametric coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    40.3132     0.1463   275.6   <2e-16 ***
TaxonSpeciesB -26.0221     0.2157  -120.6   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                      edf Ref.df     F p-value    
s(wl)               8.940  8.994 30.06  <2e-16 ***
s(wl):TaxonSpeciesA 8.001  8.000 11.61  <2e-16 ***
s(wl):TaxonSpeciesB 8.001  8.000 19.59  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.808   Deviance explained = 80.8%
GCV = 114.96  Scale est. = 114.65    n = 9918

A parte paramétrica de gam.interaction2é quase a mesma que para gam.interaction1, mas agora existem três estimativas para termos suaves, que não posso interpretar.

Agradeço antecipadamente a qualquer um que dedique algum tempo para me ajudar a entender as diferenças nos três modelos.

gam1e gam2estão bem; eles são modelos diferentes, embora estejam tentando fazer a mesma coisa, que é suavizações específicas do grupo de modelos.

O gam1formulário

y ~ f + s(x, by = f)

faz isso estimando um nível mais suave separado para cada nível de f(supondo que esse fseja um fator padrão) e, de fato, um parâmetro de suavidade separado também é estimado para cada nível suave.

O gam2formulário

y ~ f + s(x) + s(x, by = f, m = 1)

atinge o mesmo objetivo que gam1(de modelar a relação suave entre xe ypara cada nível de f), mas o faz estimando um efeito suave global ou médio de xon y(o s(x)termo) mais um termo de diferença suave (o segundo s(x, by = f, m = 1)termo). Como a penalidade aqui está na primeira derivada ( m = 1) for this difference smoother, it is penalising departure from a flat line, which when added to the global or average smooth term (s (x) `) reflete um desvio do efeito global ou médio.

gam3 Formato

y ~ s(x, by = f)

está errado, independentemente de quão bem possa caber em uma situação específica. A razão pela qual digo que está errado é que cada suavização especificada pela s(x, by = f)peça é centrada em torno de zero devido à restrição de soma a zero imposta para a identificabilidade do modelo. Como tal, não há nada no modelo que responda pela média de em cada um dos grupos definidos por . Existe apenas a média geral dada pela interceptação do modelo. Isso significa que agora, mais suave, centrada em torno de zero e cuja função de base plana foi removida da expansão de base (como é confundida com a interceptação do modelo) agora é responsável por modelar a diferença na média de da corrente grupo e a média geral (interceptação de modelo), além do efeito suave de$Y$Y Yfx$Y$xem .$Y$

Nenhum desses modelos é apropriado para seus dados; ignorando por enquanto a distribuição incorreta da resposta ( densitynão pode ser negativo e há um problema de heterogeneidade que um não-gaussiano familyiria corrigir ou resolver), você não levou em consideração o agrupamento por flor ( SampleIDno seu conjunto de dados).

Se seu objetivo é modelar Taxoncurvas específicas, um modelo do formulário seria um ponto de partida:

m1 <- gam(density ~ Taxon + s(wl, by = Taxon, k = 20) + s(SampleID, bs = 're'),
          data = df, method = 'REML')

onde adicionei um efeito aleatório SampleIDe aumentei o tamanho da expansão da base para os Taxonsuavizadores específicos.

Esse modelo m1modela as observações como provenientes de um wlefeito suave , dependendo de quais espécies ( Taxon) a observação vem (o Taxontermo paramétrico apenas define a média densitypara cada espécie e é necessário conforme discutido acima), além de uma interceptação aleatória. Tomadas em conjunto, as curvas para flores individuais surgem de versões deslocadas das Taxoncurvas específicas, com a quantidade de mudança dada pela interceptação aleatória. Este modelo pressupõe que todos os indivíduos tenham a mesma forma de suavidade que a suavizada para o particular de Taxononde vem a flor.

Outra versão deste modelo é a gam2forma acima, mas com um efeito aleatório adicional

m2 <- gam(density ~ Taxon + s(wl) + s(wl, by = Taxon, m = 1) + s(SampleID, bs = 're'),
          data = df, method = 'REML')

Esse modelo se encaixa melhor, mas acho que não está resolvendo o problema, veja abaixo. Uma coisa que acho que sugere é que o padrão ké potencialmente muito baixo para as Taxoncurvas específicas desses modelos . Ainda há muita variação suave residual que não estamos modelando, se você observar os gráficos de diagnóstico.

É provável que esse modelo seja muito restritivo para seus dados; algumas das curvas em sua plotagem de suavizações individuais não parecem ser versões simples deslocadas das Taxoncurvas médias. Um modelo mais complexo também permitiria suavizações específicas de cada indivíduo. Esse modelo pode ser estimado usando a base de interação suavefs ou fator . Ainda queremos Taxoncurvas específicas, mas também queremos ter uma suavização separada para cada uma SampleID, mas, diferentemente das bysuavizações, sugiro que inicialmente você deseje que todas essas SampleIDcurvas específicas tenham a mesma ondulação. No mesmo sentido que a interceptação aleatória que incluímos anteriormente, ofs A base adiciona uma interceptação aleatória, mas também inclui um spline "aleatório" (eu uso as aspas do susto, como em uma interpretação bayesiana do GAM, todos esses modelos são apenas variações nos efeitos aleatórios).

Este modelo é ajustado para seus dados como

m3 <- gam(density ~ Taxon + s(wl, by = Taxon, k = 20) + s(wl, SampleID, bs = 'fs'), 
          data = df, method = 'REML')

Observe que eu tenho um aumento kaqui, caso precisemos de mais oscilação nos suavizadores Taxonespecíficos. Ainda precisamos do Taxonefeito paramétrico pelas razões explicadas acima.

Esse modelo leva muito tempo para encaixar em um único núcleo com gam()- bam()provavelmente será melhor para encaixar esse modelo, pois há um número relativamente grande de efeitos aleatórios aqui.

Se compararmos esses modelos com uma versão do AIC corrigida pela seleção de parâmetros de suavidade, veremos quão melhor é este modelo, m3comparado com os outros dois, mesmo que use uma ordem de magnitude em graus de liberdade

> AIC(m1, m2, m3)
          df      AIC
m1  190.7045 67264.24
m2  192.2335 67099.28
m3 1672.7410 31474.80

Se observarmos as suavidades deste modelo, teremos uma idéia melhor de como ele está ajustando os dados:

(Observe que isso foi produzido usando draw(m3)a draw()função do meu pacote gratia . As cores no gráfico inferior esquerdo são irrelevantes e não ajudam aqui.)

SampleIDA curva ajustada de cada uma é construída a partir do termo de interceptação ou do parâmetro paramétrico, TaxonSpeciesBmais um dos dois Taxonsuavizadores específicos, dependendo de qual Taxoncada um SampleIDpertence, além da sua SampleIDsuavidade específica.

Observe que todos esses modelos ainda estão errados, pois não explicam a heterogeneidade; modelos gama ou Tweedie com um link de log seriam minhas escolhas para levar isso adiante. Algo como:

m4 <- gam(density ~ Taxon + s(wl, by = Taxon) + s(wl, SampleID, bs = 'fs'), 
          data = df, method = 'REML', family = tw())

Mas estou tendo problemas com o ajuste deste modelo no momento, o que pode indicar que é muito complexo com vários suaves wlincluídos.

Uma forma alternativa é usar a abordagem de fatores ordenados, que faz uma decomposição do tipo ANOVA nos suaves:

• Taxon termo paramétrico é retido
• s(wl)é suave que representará o nível de referência
• s(wl, by = Taxon)terá uma diferença separada suave para o outro nível. No seu caso, você terá apenas um deles.

Este modelo é montado como m3,

df <- transform(df, fTaxon = ordered(Taxon))
m3 <- gam(density ~ fTaxon + s(wl) + s(wl, by = fTaxon) +
            s(wl, SampleID, bs = 'fs'), 
          data = df, method = 'REML')

mas a interpretação é diferente; o primeiro s(wl)se referirá TaxonAe o suave implícito s(wl, by = fTaxon)será uma diferença suave entre o suave para TaxonAe o de TaxonB.

Isto é o que Jacolien van Rij escreve em sua página de tutorial:

Como configurar a interação depende do tipo de preditor de agrupamento:

• com fator incluem diferença de interceptação: Group + s(Time, by=Group)
• com fator ordenado incluem diferença de interceptação e referência suave: Group + s(Time) + s(Time, by=Group)
• com preditor binário incluem referência suave: s(Time) + s(Time, by=IsGroupChildren)

As variáveis ​​categóricas devem ser especificadas como fatores, fatores ordenados ou fatores binários com as funções R apropriadas. Para entender como interpretar as saídas e o que cada modelo pode e não pode nos dizer, consulte a página de tutorial de Jacolien van Rij diretamente. Seu tutorial também explica como ajustar GAMs de efeito misto. Para entender o conceito de interações no contexto de GAMs, esta página de tutorial de Peter Laurinec também é útil. Ambas as páginas fornecem muitas informações adicionais para executar GAMs corretamente em diferentes cenários.