Diferentes formas de modelar las interacciones entre los predictores continuos y categóricos en GAM

Question

La siguiente pregunta se basa en la discusión encontrada en esta página . Dada una variable de respuesta y una variable explicativa continua x y un factor fac es posible definir un Modelo Aditivo General (MAG) con una interacción entre x y fac usando el argumento by= . De acuerdo con el archivo de ayuda ?gam.models en el paquete R mgcv esto se puede lograr de la siguiente manera:

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

@GavinSimpson aquí sugiere un enfoque diferente:

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

He estado jugando con un tercer modelo:

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

Mis principales preguntas son: ¿algunos de estos modelos están mal, o son simplemente diferentes? En el último caso, ¿cuáles son sus diferencias? Basado en el ejemplo que voy a discutir a continuación, creo que podría entender algunas de sus diferencias, pero todavía me falta algo.

Como ejemplo voy a usar un conjunto de datos con espectros de color para las flores de dos especies de plantas diferentes medidas en diferentes lugares.

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 mayor claridad, cada línea de la figura de arriba representa el espectro de color medio predicho para cada lugar con un GAM de forma separada density~s(wl) basado en muestras de ~10 flores. Las áreas grises representan el 95% de CI para cada GAM.

Mi objetivo final es modelar el efecto (potencialmente interactivo) de Taxon y la longitud de onda wl sobre la reflectancia (denominada density en el código y el conjunto de datos) mientras se contabilizan Locality como un efecto aleatorio en un GAM de efecto mixto. Por el momento no voy a añadir la parte del efecto mixto a mi plato, que ya está bastante lleno con tratar de entender cómo modelar las interacciones.

Empezaré con la más simple de las tres GAM interactivas:

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

summary(gam.interaction0)

Produce:

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

La parte paramétrica es la misma para ambas especies, pero se colocan diferentes estrías para cada especie. Es un poco confuso tener una parte paramétrica en el resumen de las GAM, que no son paramétricas. @IsabellaGhement explica:

Si miras las gráficas de los efectos suaves estimados (smooths) correspondiente a su primer modelo, notará que son centrado en el cero. Así que tienes que "cambiar" esos alisados (si el intercepción estimada es positiva) o hacia abajo (si la intercepción estimada es negativo) para obtener las funciones suaves que pensabas que eras estimando. En otras palabras, necesitas añadir la intercepción estimada a las facilidades para llegar a lo que realmente quieres. Para su primer modelo, el Se supone que el "cambio" es el mismo para ambos tipos de lisos.

Sigamos adelante:

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

summary(gam.interaction1)

Da:

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

Ahora, cada especie también tiene su propia estimación paramétrica.

El siguiente modelo es el que me cuesta entender:

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

No tengo una idea clara de lo que representan estos gráficos.

summary(gam.interaction2)

Da:

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

La parte paramétrica de gam.interaction2 es más o menos lo mismo que para gam.interaction1 pero ahora hay tres estimaciones para términos suaves, que no puedo interpretar.

Gracias de antemano a todos los que se tomen el tiempo de ayudarme a entender las diferencias entre los tres modelos.

Answer 1

gam1 y gam2 están bien; son modelos diferentes, aunque intentan hacer lo mismo, que es modelar suavidades específicas del grupo.

El gam1 formulario

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

lo hace estimando un suavizador distinto para cada nivel de f (suponiendo que f es un factor estándar), y de hecho, se estima un parámetro de suavidad independiente para cada suavidad también.

El gam2 formulario

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

consigue el mismo objetivo que gam1 (de modelar la relación fluida entre x y y para cada nivel de f ) pero lo hace estimando un efecto global o medio suave de x en y (el s(x) plazo) más un término de diferencia suave (el segundo s(x, by = f, m = 1) término). Como la penalización aquí es sobre la primera 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)`) refleja una desviación del efecto global o medio.

gam3 formulario

y ~ s(x, by = f)

es errónea, independientemente de lo bien que pueda encajar en una situación concreta. La razón por la que digo que es errónea es que cada liso especificado por el s(x, by = f) está centrada en cero debido a la restricción de suma a cero impuesta para la identificabilidad del modelo. Por lo tanto, no hay nada en el modelo que explique la media de $Y$ en cada uno de los grupos definidos por f . Sólo existe la media global dada por el intercepto del modelo. Esto significa que el suavizador, que está centrado en cero y al que se le ha eliminado la función de base plana de la expansión de base de x (ya que se confunde con el intercepto del modelo) es ahora responsable de modelar tanto la diferencia en la media de $Y$ para el grupo actual y la media global (intercepción del modelo), más el efecto suave de x en $Y$ .

Sin embargo, ninguno de estos modelos es apropiado para sus datos; ignorando por ahora la distribución errónea de la respuesta ( density no puede ser negativo y hay un problema de heterogeneidad que un no gaussiano family arreglaría o abordaría), no has tenido en cuenta la agrupación por flores ( SampleID en su conjunto de datos).

Si su objetivo es modelar Taxon curvas específicas, entonces un modelo de la forma sería un punto de partida:

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

donde he añadido un efecto aleatorio para SampleID y ha aumentado el tamaño de la expansión de la base para el Taxon lisos específicos.

Este modelo, m1 modeliza las observaciones como si procedieran de un wl efecto en función de las especies ( Taxon ) la observación proviene de (el Taxon El término paramétrico sólo establece la media density para cada especie y es necesario como se ha comentado anteriormente), más un intercepto aleatorio. En conjunto, las curvas de las flores individuales surgen de versiones desplazadas del Taxon curvas específicas, con la cantidad de desplazamiento dada por el intercepto aleatorio. Este modelo supone que todos los individuos tienen la misma forma de liso dada por el liso para el Taxon de la que proviene esa flor individual.

Otra versión de este modelo es el gam2 forma de arriba pero con un efecto aleatorio añadido

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

Este modelo se ajusta mejor, pero no creo que resuelva el problema en absoluto, véase más abajo. Una cosa que creo que sí sugiere es que el k es potencialmente demasiado bajo para el Taxon curvas específicas en estos modelos . Todavía hay una gran cantidad de variación suave residual que no estamos modelando si nos fijamos en los gráficos de diagnóstico.

Es más que probable que este modelo sea demasiado restrictivo para sus datos; algunas de las curvas en su gráfico de los alisados individuales no parecen ser simples versiones desplazadas del Taxon curvas medias. Un modelo más complejo permitiría también la existencia de suavizaciones específicas para cada individuo. Un modelo de este tipo puede estimarse utilizando el fs o base de interacción entre factores y lisos . Todavía queremos Taxon curvas específicas, pero también queremos tener un liso separado para cada SampleID pero a diferencia del by suaviza, yo sugeriría que inicialmente usted quiere que todos los SampleID -Curvas específicas para tener el mismo meneo. En el mismo sentido que el intercepto aleatorio que incluimos antes, el fs añade un intercepto aleatorio, pero también incluye un spline "aleatorio" (utilizo las comillas porque en una interpretación bayesiana del GAM, todos estos modelos son sólo variaciones de los efectos aleatorios).

Este modelo se ajusta a sus datos como

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

Tenga en cuenta que he aumentado k aquí, en caso de que necesitemos más meneo en el Taxon -suavidades específicas. Todavía necesitamos el Taxon efecto paramétrico por las razones explicadas anteriormente.

Ese modelo toma un largo tiempo para que quepa en un solo núcleo con gam() - bam() probablemente se ajustará mejor a este modelo, ya que hay un número relativamente grande de efectos aleatorios.

Si comparamos estos modelos con una versión del AIC corregida por la selección de parámetros de suavidad, veremos que este último modelo es mucho mejor, m3 es comparado con los otros dos a pesar de que utiliza un orden de magnitud más de grados de libertad

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

Si observamos los alisados de este modelo, podemos hacernos una mejor idea de cómo se ajusta a los datos:

(Nota: esto fue producido usando draw(m3) utilizando el draw() de mi gratia paquete. Los colores de la gráfica inferior izquierda son irrelevantes y no ayudan aquí).

Cada SampleID La curva ajustada se construye a partir del intercepto o del término paramétrico TaxonSpeciesB más uno de los dos Taxon -específicos, dependiendo de a qué Taxon cada SampleID pertenece, más es propio SampleID -especialmente suave.

Tenga en cuenta que todos estos modelos siguen siendo erróneos, ya que no tienen en cuenta la heterogeneidad; los modelos gamma o Tweedie con un enlace logarítmico serían mis opciones para llevar esto más lejos. Algo así como

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

Pero estoy teniendo problemas con el ajuste de este modelo en este momento, lo que podría indicar que es demasiado complejo con múltiples suavizados de wl incluido.

Una forma alternativa es utilizar el enfoque del factor ordenado, que realiza una descomposición similar a la de ANOVA en los suavizados:

• Taxon se mantiene el término paramétrico
• s(wl) es un liso que representará el referencia nivel
• s(wl, by = Taxon) tendrá un diferencia suave para cada uno de los otros niveles. En tu caso sólo tendrás uno de ellos.

Este modelo se ajusta 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')

pero la interpretación es diferente; la primera s(wl) se referirá a TaxonA y la suavidad implícita en s(wl, by = fTaxon) será una diferencia suave entre el suave para TaxonA y la de TaxonB .

Answer 2

Esto es lo que Jacolien van Rij escribe en su página de tutoriales:

La forma de establecer la interacción depende del tipo de agrupación predictor:

• con el factor incluye la diferencia de intercepción: Group + s(Time, by=Group)
• con factor ordenado incluyen diferencia de intercepción y suave de referencia: Group + s(Time) + s(Time, by=Group)
• con predictor binario incluyen referencia suave: s(Time) + s(Time, by=IsGroupChildren)

Las variables categóricas deben especificarse como factores, factores ordenados o factores binarios con las funciones R adecuadas. Para entender cómo interpretar los resultados y lo que cada modelo puede y no puede decirnos, véase Página del tutorial de Jacolien van Rij directamente. Su tutorial también explica cómo ajustar los GAM de efectos mixtos. Para entender el concepto de interacciones en el contexto de los GAM, este La página del tutorial de Peter Laurinec también es útil. Ambas páginas proporcionan mucha más información para ejecutar los GAM correctamente en diferentes escenarios.