La siguiente pregunta se basa en la discusión que se encuentra en esta página . Dada una variable de respuesta y, una variable explicativa continua xy un factor fac, es posible definir un Modelo Aditivo General (GAM) con una interacción entre xy facusando 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 preguntas principales son: ¿algunos de estos modelos son incorrectos o simplemente son 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 utilizar un conjunto de datos con espectros de color para 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 en la figura de arriba representa el espectro de color medio predicho para cada ubicación con un GAM de forma separado density~s(wl)basado en muestras de ~ 10 flores. Las áreas grises representan un IC del 95% para cada GAM.
Mi objetivo final es modelar el efecto (potencialmente interactivo) Taxony la longitud wlde onda en la reflectancia (referido como densityen el código y el conjunto de datos) mientras que considero Localitycomo un efecto aleatorio en un GAM de efectos mixtos. Por el momento no agregaré la parte de efectos mixtos a mi plato, que ya está lo suficientemente lleno para tratar de entender cómo modelar las interacciones.
Comenzaré con el más simple de los tres GAM interactivos:
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 = 9918La parte paramétrica es la misma para ambas especies, pero se ajustan diferentes splines para cada especie. Es un poco confuso tener una parte paramétrica en el resumen de GAM, que no son paramétricos. @IsabellaGhement explica:
Si observa las gráficas de los efectos suaves estimados (suavizados) correspondientes a su primer modelo, notará que están centrados alrededor de cero. Por lo tanto, debe 'desplazar' esos suavizados hacia arriba (si la intersección estimada es positiva) o hacia abajo (si la intersección estimada es negativa) para obtener las funciones uniformes que pensó que estaba estimando. En otras palabras, debe agregar la intercepción estimada a los suavizados para obtener lo que realmente desea. Para su primer modelo, se supone que el 'cambio' es el mismo para ambos suavizados.
Hacia 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 = 9918Ahora, cada especie también tiene su propia estimación paramétrica.
El siguiente modelo es el que tengo problemas para 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 = 9918La parte paramétrica de gam.interaction2es aproximadamente la misma que para gam.interaction1, pero ahora hay tres estimaciones para términos suaves, que no puedo interpretar.
Gracias de antemano a cualquiera que se tome el tiempo para ayudarme a comprender las diferencias en los tres modelos.
fuente
gam1más algo para elSampleIDefecto más que necesita hacer algo sobre el problema de la varianza no constante; Estos datos no parecen estar distribuidos condicionalmente gaussianos debido al límite inferior.Respuestas:
gam1ygam2están bien son modelos diferentes, aunque intentan hacer lo mismo, que son suavizaciones específicas del grupo de modelos.La
gam1formahace esto al estimar un suavizador separado para cada nivel de
f(suponiendo quefsea un factor estándar), y de hecho, también se estima un parámetro de suavidad separado para cada suavizado.La
gam2formalogra el mismo objetivo que
gam1(de modelar la relación fluida entrexyypara cada nivel def) pero lo hace estimando un efecto suave global o promedio dexony(els(x)término) más un término de diferencia suave (el segundos(x, by = f, m = 1)término). Como la penalización aquí está en 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 promedio.gam3formarestá mal independientemente de lo bien que pueda encajar en una situación particular. La razón por la que digo que está mal es que cada suavizado especificado por laY en cada uno de los grupos definidos por Y para el grupo actual y la media general (modelo de intercepción), más el buen efecto de Y .
s(x, by = f)parte se centra en cero debido a la restricción de suma a cero impuesta para la identificabilidad del modelo. Como tal, no hay nada en el modelo que explique la media def. Solo existe la media general dada por la intercepción del modelo. Esto significa que ahora el más suave, que está centrado en torno a cero y que ha tenido la función de base plana eliminada de la expansión de la basex(como se confunde con la intercepción del modelo) ahora es responsable de modelar tanto la diferencia en la media dexsobreSin embargo, ninguno de estos modelos es apropiado para sus datos; ignorando por ahora la distribución incorrecta de la respuesta (
densityno puede ser negativa y hay un problema de heterogeneidad que un no gaussianofamilysolucionaría o resolvería), no ha tenido en cuenta la agrupación por flor (SampleIDen su conjunto de datos).Si su objetivo es modelar
Taxoncurvas específicas, un modelo de la forma sería un punto de partida:donde agregué un efecto aleatorio
SampleIDy aumenté el tamaño de la expansión de base para losTaxonsuavizados específicos.Este modelo,
m1modela las observaciones como resultado de unwlefecto suave dependiendo de la especie (Taxon) de la que proviene la observación (elTaxontérmino paramétrico solo establece la mediadensitypara cada especie y se necesita como se discutió anteriormente), más una intercepción aleatoria. En conjunto, las curvas para flores individuales surgen de versiones desplazadas de lasTaxoncurvas específicas, con la cantidad de desplazamiento dada por la intercepción aleatoria. Este modelo supone que todos los individuos tienen la misma forma de liso que el liso para el particular delTaxonque proviene la flor individual.Otra versión de este modelo es la
gam2forma de arriba pero con un efecto aleatorio agregadoEste modelo se ajusta mejor, pero no creo que esté resolviendo el problema en absoluto, ver más abajo. Una cosa que creo que sugiere es que el valor predeterminado
kes potencialmente demasiado bajo para lasTaxoncurvas específicas en estos modelos . Todavía hay mucha variación suave residual que no estamos modelando si nos fijamos en los gráficos de diagnóstico.Es muy probable que este modelo sea demasiado restrictivo para sus datos; Algunas de las curvas en su diagrama de los suavizados individuales no parecen ser versiones simples desplazadas de las
Taxoncurvas promedio. Un modelo más complejo permitiría también suavizaciones específicas para cada individuo. Tal modelo puede estimarse utilizando la base de interacciónfso factor-liso . Todavía queremosTaxoncurvas específicas, pero también queremos tener un suavizado separado para cada unaSampleID, pero a diferencia de losbysuavizados, sugeriría que inicialmente desee que todas esasSampleIDcurvas específicas tengan la misma ondulación. En el mismo sentido que la intersección aleatoria que incluimos anteriormente, elfsbase agrega una intercepción aleatoria, pero también incluye una spline "aleatoria" (uso las comillas de miedo como en una interpretación bayesiana del GAM, todos estos modelos son solo variaciones de efectos aleatorios).Este modelo se ajusta a sus datos como
Tenga en cuenta que he aumentado
kaquí, en caso de que necesitemos más ondulación en los suavizadosTaxonespecíficos. Todavía necesitamos elTaxonefecto paramétrico por las razones explicadas anteriormente.Ese modelo tarda mucho tiempo en ajustarse a un solo núcleo
gam();bam()lo más probable es que sea mejor para ajustar este modelo, ya que aquí hay una cantidad relativamente grande de efectos aleatorios.Si comparamos estos modelos con una versión de AIC corregida por selección de parámetros de suavidad, vemos cuán dramáticamente mejor este último modelo
m3se compara con los otros dos, aunque utiliza un orden de magnitud más grados de libertadSi observamos los suavizados de este modelo, obtenemos una mejor idea de cómo se ajustan los datos:
(Tenga en cuenta que esto se produjo utilizando
draw(m3)ladraw()función de mi paquete gratia . Los colores en el diagrama inferior izquierdo son irrelevantes y no ayudan aquí).Cada
SampleIDcurva ajustada se construye a partir de la intersección o el término paramétricoTaxonSpeciesBmás uno de los dosTaxonsuavizados específicos, dependiendo de a qué perteneceTaxoncada unoSampleID, más su propioSampleIDsuavizado específico.Tenga en cuenta que todos estos modelos siguen siendo incorrectos, ya que no tienen en cuenta la heterogeneidad; Los modelos gamma o Tweedie con un enlace de registro serían mis opciones para llevar esto más lejos. Algo como:
Pero estoy teniendo problemas con este modelo en este momento, lo que podría indicar que es demasiado complejo con múltiples suavidades
wlincluidas.Una forma alternativa es utilizar el enfoque de factor ordenado, que realiza una descomposición similar a ANOVA en los suavizados:
Taxontérmino paramétrico se retienes(wl)es un suave que representará el nivel de referencias(wl, by = Taxon)tendrá una diferencia separada suave para cada otro nivel. En su caso, solo tendrá uno de estos.Este modelo está equipado como
m3,pero la interpretación es diferente; el primero
s(wl)se referiráTaxonAy el suave implicado pors(wl, by = fTaxon)será una diferencia suave entre el suave paraTaxonAy el deTaxonB.fuente
SampleIDes un espectrograma de una sola flor, cada uno de una planta diferente, por lo que no creoSampleIDque deba especificarse como aleatorio (pero corrígeme si me equivoco). De hecho, he usado un modelo similar al suyom3con unTaxonfactor ordenado, pero que especifica+ s(Locality, bs="re") + s(Locality, wl, bs="re")como aleatorio. Analizaré los problemas que plantea sobre la distribución de los residuos y la heterocedasticidad. ¡Salud!SampleIDcomo aleatoria los datos de una sola flor que probablemente estén relacionados y más si la función completa se relaciona con la flor, por lo que, en cierto sentido, las funciones (suaves) son aleatorias. También podría necesitar un efecto aleatorio simple para la planta si hubiera varias flores por planta y múltiples plantas por taxón en el estudio (use elbs = 're'"suave" que mencioné anteriormente en la respuesta.m3confamily = Gamma(link = 'log')ofamily = tw()me estaba poniendo problemas reales con mgcv no ser capaz de encontrar buenos valores de partida y otros errores que causan mgcv a mierda, que no he llegado a la parte inferior de todavía. Ciertamente, a partir de los datos que proporcionó, un modelo gaussiano no es correcto. Obtuve un gaussiano con un enlace de registro que encaja y me ayudó, pero tampoco captura toda la heterogeneidad.Esto es lo que Jacolien van Rij escribe en su página de tutorial:
Las variables categóricas deben especificarse como factores, factores ordenados o factores binarios con las funciones R apropiadas. Para comprender cómo interpretar los resultados y lo que cada modelo puede y no puede decirnos, consulte la página del tutorial de Jacolien van Rij directamente. Su tutorial también explica cómo ajustar los GAM de efectos mixtos. Para comprender el concepto de interacciones en el contexto de los GAM, esta página de tutorial de Peter Laurinec también es útil. Ambas páginas proporcionan mucha más información para ejecutar GAM correctamente en diferentes escenarios.
fuente