mgcv：修復GAMM中的平滑參數和模型嵌套的有效性

Question

我目前正在嘗試對來自mgcv包的協變量x（Age）和性別0-1變量與gamm之間的交互進行建模。在爲每個性別指定一個主要模型（讓我們稱之爲M0）之後，我想測試一個更簡單的假設，即性別之間的差異是線性的（而不是任意平滑的）。我有以下兩個問題：

1. 當試圖鳥巢模型正確，我想從M0走出了0性別流暢平滑參數，並在簡單的模型規範使用它。但我得到了以下錯誤消息：

錯誤gamm.setup（GP，pterms = pTerms，數據= MF，繩結=結，parametric.only = FALSE，：

GAMM不能處理連接的平滑參數（可能是從使用`id或自適應潤滑肌膚的）

第二個問題是比較笨之一。在本模型甚至嵌套時我從光滑去爲每個性別0性別平滑和1性別線性差異？

下面是一個例子。我模擬了一些隨機數據，所以數據不會顯示我實際數據的行爲，但問題依然存在。

library(mgcv) 

### Simulate random data 
x <- rnorm(100, mean = 10, sd = 1.5) 
y <- rnorm(100, mean = 1, sd = 0.025) 

id <- sample(1:10, size = 100, replace = T) 
id <- as.factor(id) 

gender <- sample(c(0,1), size = 100, replace = T) 


### Specify main model, M0 
ctrl <- list(niterEM=0,optimMethod="L-BFGS-B", msMaxIter = 100) 
M0 <- gamm(y~s(x, by = as.factor(gender)) + gender, 
      random=list(id=~1+x), control=ctrl) 

### Specify model with linear difference between gender0 and gender1 
M1 <- gamm(y~s(x) + gender:x + gender, 
      random=list(id=~1+x), control=ctrl) 

### Testing 
anova(M0$lme, M1$lme) 

### Problems: 
sp0 <- M0$gam$sp[1] 
M1 <- gamm(y~s(x, sp = sp0) + gender:x + gender, 
      random=list(id=~1+x), control=ctrl) 

錯誤gamm.setup（GP，pterms = pTerms，數據= MF，繩結=疙瘩，parametric.only = FALSE，：

GAMM不能處理從連接平滑參數（可能使用'身份證」或自適應潤滑肌膚）的

有什麼想法？在此先感謝。

Answer 1

關於gamm錯誤

這是一件非常有趣的事情......那麼，我應該首先解釋邏輯。

原則它是違法的gamm修復平滑參數，因爲gamm將把光滑如隨機效應的波浪組成部分，其中的差額通過lme估計（你有高斯數據）。如果你試圖修正平滑參數，那實際上是說你想修正隨機效應的方差。那麼，lme不允許你這樣做（我懷疑這種嘗試在混合建模中是否合法。）

因此gamm將在平滑參數禁用任何可能的限制，其中包括：

1. 下界平滑參數的，經由min.sp;
2. 在s()連接平滑共享相同的平滑參數，經由id;
3. 直接指定經由sp平滑參數，經由s()。

前兩個是完全檢查。 gamm沒有min.sp的說法，如gam;即使你通過...指定它，有沒有機會，它會習慣（如後面這是一個的過程中gamm.setup傳遞到gam.setupNULL，所以你指定的min.sp被忽略）。 id的規格也會被你看到的錯誤信息檢測到，但是你當然沒有指定id，所以上面的錯誤在這裏沒有報告正確的問題，因此是一個錯誤。

第三個，卻沒有被直接通過gamm檢查。理想情況下，只要gamm/gam公式被解釋（通過interpret.gam），sp應重置爲-1（如果不是很容易），並且應該發出關於此的警告消息。但是，這部分不見了。因此，目前你能做的最好的事情就是不指定sp。

---

您是否有有效的模型嵌套？

現在讓我們談談您對嵌套的關注。 嵌套與基礎設置有關，而不是選擇平滑參數。只要您具有相同的基礎組（相同的基礎類型，相同的數字和/或「節點」的位置），模型矩陣將是相同的。例如，在你的模型M0和M1，你有s()與mgcv默認bs = 'tp', k = 10相同的配置。所以s()的設計矩陣在你的兩個模型中是相同的。 by = factor(gender)只是將s()複製到您的主要型號M0中的各個級別gender。也許它不容易看到，但實際上你的M1嵌套在M0。

讓我們考慮一個小例子來驗證這一點。爲了簡單起見，我不會使用從mgcv，但使用poly(x, degree = 2)（想象它是s(x)）。讓我們先產生一些玩具數據：

x <- 1:10 
f <- gl(2, 5, labels = c("M", "F")) 

由於f不是一個有序的因素，s(x, by = factor(f))通過複製s(x)爲f各級產生設計矩陣：

## original design matrix for `s(x)` 
X0 <- poly(x, 2) 
## design matrix for `f`, without contrasting 
Xf <- model.matrix(~ f + 0) 
## design matrix for level `M` 
X1 <- X0 * Xf[, 1] 
## design matrix for level `F` 
X2 <- X0 * Xf[, 2] 
## design matrix for `s(x, by = f)` "please, imagine it as `poly`" 
X <- cbind(X1, X2) 

#    1   2   1   2 
# [1,] -0.49543369 0.52223297 0.00000000 0.00000000 
# [2,] -0.38533732 0.17407766 0.00000000 0.00000000 
# [3,] -0.27524094 -0.08703883 0.00000000 0.00000000 
# [4,] -0.16514456 -0.26111648 0.00000000 0.00000000 
# [5,] -0.05504819 -0.34815531 0.00000000 0.00000000 
# [6,] 0.00000000 0.00000000 0.05504819 -0.34815531 
# [7,] 0.00000000 0.00000000 0.16514456 -0.26111648 
# [8,] 0.00000000 0.00000000 0.27524094 -0.08703883 
# [9,] 0.00000000 0.00000000 0.38533732 0.17407766 
#[10,] 0.00000000 0.00000000 0.49543369 0.52223297 

你的第二個模型M1，只有一個平滑術語s(x)，其設計矩陣爲X0。

下面是我們如何可以看到您的M1嵌套在M0：

1. 作爲X1 + X2 = X0的s(x)設計矩陣和s(x, by = f)具有相同的列跨度，從而s(x)嵌套在s(x, by = f);
2. 因爲x嵌套在s(x),x:f嵌套在s(x, by = f)。

---

我會做

什麼雖然你的模型已經很好地嵌套，M0沒有與你M1可比解釋的主力機型。您的主要型號M0將以每個級別的獨立平滑度結束，而您的M1重點關注兩組之間的差異。

如果我們能夠控制M0接受「參考平滑+差異平滑」的形式，那將是一件好事。然後，如果差異平滑發現一條線，而實際上沒有擬合，我們已經知道沒有證據表明組之間存在非線性差異。

在mgcv中，如果您的因子是有序的，則會構建差異平滑。因此，我建議你用適合您的主要型號：

gender1 <- ordered(gender) ## create an ordered factor 
s(x) + s(x, by = gender1) + gender 

如果估計結果顯示爲一條線的區別順利s(x, by = gender1)，你知道你可以改爲適合簡單的模型

s(x) + gender:x + gender 

即使不使用F-test。

請注意，爲了構建「差異」平滑，具有排序因子by非常重要。如果你

s(x) + s(x, by = gender) + gender ## note, it is "gender" in "by" 

s(x)和s(x, by = gender)是完全線性相關。得到的模型矩陣將是等級缺陷。

---

一些後續

忘在我的例子，我在第一比較參數化作爲s(x, by = as.factor(gender))和s(x) + s(x, by = gender)由AIC相同的模型，包括（請回想gender是0-1的數值變量）。這些模型在統計上是等價的，但在這些情況下明顯估計平滑參數是不同的，因此AIC有所不同。

哦，是的。您的gender是二元的，因此數值爲by也是構建差異平滑的好主意。但要小心。數字by不產生居中平滑。因此，s(x) + s(x, by = gender)將隱含有一個攔截列，與模型攔截混淆。你應該去s(x) + s(x, by = gender) - 1。