GAM 中連續和分類預測變量之間交互作用的不同建模方式
以下問題建立在此頁面上的討論之上。給定一個響應變量
y、一個連續解釋變量x和一個因子fac，可以定義一個通用加法模型 (GAM)，其中包含x和fac使用參數之間的交互by=。根據R包中的幫助文件 ，可以這樣完成：?gam.models``mgcvgam1 <- gam(y ~ fac +s(x, by = fac), ...)@GavinSimpson 在這裡提出了一種不同的方法：
gam2 <- gam(y ~ fac +s(x) +s(x, by = fac, m=1), ...)我一直在玩第三個模型：
gam3 <- gam(y ~ s(x, by = fac), ...)我的主要問題是：其中一些模型是錯誤的，還是只是不同？在後一種情況下，它們有什麼區別？ 根據我將在下面討論的示例，我想我可以理解它們的一些差異，但我仍然缺少一些東西。
作為一個例子，我將使用一個數據集，其中包含在不同位置測量的兩種不同植物物種的花朵的色譜。
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)這些是這兩個物種在局部級別的平均色譜（使用了滾動方式）：
每種顏色代表不同的物種。每一行代表不同的地方。
我的最終目標是模擬（可能交互的）
Taxon和波長wl對反射率百分比的影響（density在代碼和數據集中稱為），同時考慮Locality混合效應 GAM 中的隨機效應。目前我不會將混合效果部分添加到我的盤子中，這已經足夠了解如何對交互進行建模了。我將從三個交互式 GAM 中最簡單的一個開始：
gam.interaction0 <- gam(density ~ s(wl, by = Taxon), data = df) # common intercept, different slopes plot(gam.interaction0, pages=1)
summary(gam.interaction0)產生：
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兩個物種的參數部分相同，但每個物種都擬合了不同的樣條。在 GAM 的摘要中包含參數部分有點令人困惑，它是非參數的。@IsabellaGhement解釋說：
如果您查看與您的第一個模型相對應的估計平滑效果（平滑）圖，您會注意到它們以零為中心。因此，您需要向上（如果估計的截距為正）或向下（如果估計的截距為負）“移動”這些平滑，以獲得您認為正在估計的平滑函數。換句話說，您需要將估計的截距添加到平滑中以獲得您真正想要的。對於您的第一個模型，假定兩個平滑的“移位”相同。
繼續：
gam.interaction1 <- gam(density ~ Taxon +s(wl, by = Taxon, m=1), data = df) plot(gam.interaction1,pages=1)
summary(gam.interaction1)給出：
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現在，每個物種也有自己的參數估計。
下一個模型是我難以理解的模型：
gam.interaction2 <- gam(density ~ Taxon + s(wl) + s(wl, by = Taxon, m=1), data = df) plot(gam.interaction2, pages=1)
我不清楚這些圖表代表什麼。
summary(gam.interaction2)給出：
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的參數部分與
gam.interaction2的 大致相同gam.interaction1，但現在平滑項有三個估計值，我無法解釋。提前感謝任何願意花時間幫助我了解三種模型差異的人。
gam1並且gam2很好；它們是不同的模型，儘管它們試圖做同樣的事情，即模型組特定的平滑。
gam1表格_y ~ f + s(x, by = f)
f通過為每個平滑度（假設這是一個標準因子）估計一個單獨的平滑器來做到這一點f，實際上，每個平滑度也估計一個單獨的平滑度參數。
gam2表格_y ~ f + s(x) + s(x, by = f, m = 1)實現了與（對 的每個級別
gam1之間的平滑關係進行建模）相同的目標，但它是通過估計（項）加上平滑差分項（第二項）的全局或平均平滑效果來實現的。由於這裡的懲罰是在一階導數（s(x)）上反映了與全局或平均效應的偏差。xyfxys(x)s(x, by = f, m = 1)``m = 1) for this difference smoother, it is penalising departure from a flat line, which when added to the global or average smooth term (`
gam3形式y ~ s(x, by = f)無論它在特定情況下的適應程度如何，都是錯誤的。我說它是錯誤的原因是
s(x, by = f)由於對模型可識別性施加的總和為零約束，零件指定的每個平滑都以零為中心。因此，模型中沒有任何內容可以解釋 $ Y $ 在由 定義的每個組中f。模型截距僅給出整體平均值。這意味著以零為中心並且已經從 的基擴展中移除了平坦基函數x（因為它與模型截距混淆）現在負責對 $ Y $ 對於當前組和整體平均值（模型截距），加上xon的平滑效果 $ Y $ .但是，這些模型都不適合您的數據；現在忽略響應的錯誤分佈（
density不能是負數，並且存在非高斯family可以解決或解決的異質性問題），您沒有考慮按花分組（SampleID在您的數據集中）。如果您的目標是對
Taxon特定曲線進行建模，那麼形式的模型將是一個起點：m1 <- gam(density ~ Taxon + s(wl, by = Taxon, k = 20) + s(SampleID, bs = 're'), data = df, method = 'REML')我為特定平滑添加了隨機效果
SampleID並提高了基礎擴展的大小。Taxon該模型將觀察建模為
m1來自平滑wl效應，具體取決於Taxon觀察來自哪個物種 ( 總之，單個花朵的曲線來自特定曲線的偏移版本，偏移量由隨機截距給出。該模型假設所有個體都具有與單個花朵所來自的特定物體的光滑度相同的光滑形狀。Taxon``density``Taxon``Taxon該模型的另一個版本是
gam2上面的形式，但增加了隨機效應m2 <- gam(density ~ Taxon + s(wl) + s(wl, by = Taxon, m = 1) + s(SampleID, bs = 're'), data = df, method = 'REML')這個模型更適合，但我認為它根本不能解決問題，見下文。我認為它確實表明的一件事是，對於這些模型中的特定曲線，默認
k值可能太低。如果您查看診斷圖，仍然有很多剩餘平滑變化我們沒有建模。Taxon該模型很可能對您的數據過於嚴格；您的各個平滑圖中的某些曲線似乎不是
Taxon平均曲線的簡單移位版本。更複雜的模型也將允許個人特定的平滑。這樣的模型可以使用fs或因子平滑交互基礎進行估計。我們仍然需要Taxon特定的曲線，但我們也希望為每個曲線設置單獨的平滑SampleID，但與by平滑不同，我建議您最初希望所有這些SampleID特定曲線具有相同的擺動。在與我們之前包含的隨機截距相同的意義上，fs基礎添加了一個隨機截距，但還包括一個“隨機”樣條曲線（我在 GAM 的貝葉斯解釋中使用了嚇人的引號，所有這些模型都只是隨機效應的變體）。該模型適合您的數據
m3 <- gam(density ~ Taxon + s(wl, by = Taxon, k = 20) + s(wl, SampleID, bs = 'fs'), data = df, method = 'REML')請注意，我在
k這裡增加了，以防我們需要在特定的Taxon平滑中增加更多的擺動。Taxon由於上述原因，我們仍然需要參數效應。該模型需要很長時間才能適應單個內核
gam()-bam()因為這裡有相對大量的隨機效應，所以很可能會更好地適應這個模型。如果我們將這些模型與 AIC 的平滑參數選擇校正版本進行比較，我們會看到後一種模型 與其他兩個模型相比有多麼顯著，
m3即使它使用了一個數量級的自由度> AIC(m1, m2, m3) df AIC m1 190.7045 67264.24 m2 192.2335 67099.28 m3 1672.7410 31474.80如果我們看一下這個模型的平滑，我們會更好地了解它是如何擬合數據的：
（請注意，這是使用我的gratia
draw(m3)包中的函數生成的。左下圖中的顏色無關緊要，在這裡沒有幫助。）draw()每個
SampleID的擬合曲線是由截距或參數項TaxonSpeciesB加上兩個Taxon特定平滑之一（取決於Taxon每個平滑）SampleID，加上它自己的SampleID特定平滑。請注意，所有這些模型仍然是錯誤的，因為它們沒有考慮異質性；帶有日誌鏈接的 gamma 或 Tweedie 模型將是我更進一步的選擇。就像是：
m4 <- gam(density ~ Taxon + s(wl, by = Taxon) + s(wl, SampleID, bs = 'fs'), data = df, method = 'REML', family = tw())但我目前在擬合這個模型時遇到了問題，這可能表明它過於復雜，
wl包含多個平滑。另一種形式是使用有序因子方法，它對平滑進行類似 ANOVA 的分解：
Taxon保留參數項s(wl)是一個平滑，將代表參考水平s(wl, by = Taxon)將有一個單獨的差異平滑對於其他級別。在您的情況下，您將只有其中之一。這個模型適合像
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')但解釋不同；第一個
s(wl)是指TaxonA和 所暗示的平滑將是平滑 的和 的s(wl, by = fTaxon)之間的平滑區別。TaxonA``TaxonB
