2

私は現在、パッケージから共変量x(年齢) と性別 0-1 変数の間の相互作用をモデル化しようとしています。性別ごとに滑らかな項を持つメイン モデル ( と呼びましょう) を指定した後、性別間の差が線形である (任意に滑らかではない) という単純な仮説をテストしたいと思います。次の 2 つの質問があります。gammmgcvM0

  1. モデルを適切にネストしようとするとき、 から 0-gender Smooth のスムージング パラメータを取得しM0、より単純なモデル仕様で使用したいと考えています。しかし、次のエラー メッセージが表示されます。

gamm.setup(gp、pterms = pTerms、データ = mf、ノット = ノット、parametric.only = FALSE、のエラー:

gamm はリンクされた平滑化パラメータを処理できません (おそらく「id」または適応平滑化の使用による)

  1. 2 番目の質問は、より愚かな質問です。性別ごとのスムーズから性別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 はリンクされた平滑化パラメータを処理できません (おそらく「id」または適応平滑化の使用による)

何かご意見は？前もって感謝します。

4

1 に答える 1

3

gammエラーについて

これは非常に興味深いことです...まず、ロジックを説明する必要があります。

原則として、でスムージング パラメータを修正することは違法ですgammgamm。これは、 がスムースの波状コンポーネントをランダム効果として扱い、その分散を推定するlmeためです (ガウス データがある場合)。平滑化パラメーターを修正しようとすると、それは基本的に、ランダム効果の分散を修正したいということです。まあ、lmeこれを行うことはできません (また、そのような試みが混合モデリングで合法であるかどうかも疑問です)。

したがってgamm、次のようなスムージング パラメーターの制約をすべて無効にします。

  1. 平滑化パラメータの下限min.sp。
  2. idinを介して、同じスムージング パラメータを共有するリンクされたスムースs()。
  3. sp、経由で直接指定された平滑化パラメータs()。

最初の 2 つは完全にチェックされます。;のような引数gammはありません。経由で指定しても、使用される可能性はありません（後でduringに渡されるため、指定は無視されます）。の指定は、表示されたエラー メッセージによっても検出されますが、もちろん指定しなかったため、上記のエラーはここで正しい問題を報告していないため、バグです。min.spgam...NULLgam.setupgamm.setupmin.spidid

3 番目のものは、実際には を介し​​て直接チェックされていませんgamm。理想的には、gamm/gam式が (によってinterpret.gam)解釈されたらすぐに、そうでない場合はspにリセットし-1、これに関する警告メッセージを発行する必要があります。ただし、この部分が欠けています。したがって、現時点でできる最善のことは、指定しないことspです。

有効なモデルのネストがありますか?

次に、ネスティングに関する懸念に目を向けましょう。ネスティングは、平滑化パラメーターの選択ではなく、基本設定に関連しています。基底のセットが同じである限り (基底の種類、「結び目」の数および/または位置が同じ)、モデル マトリックスは同じになります。たとえば、モデルM0とでは、デフォルトのとM1同じ構成になっています。したがって、 の計画行列は 2 つのモデルで同じです。はこれをメイン モデルのすべてのレベルに複製します。わかりにくいかもしれませんが、実際には にネストされています。s()mgcvbs = 'tp', k = 10s()by = factor(gender)s()genderM0M1M0

これを確認するために、小さな例を考えてみましょう。簡単にするために、s(x)fromを使用せずに( であると想像して)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

2 番目のモデルには、計画行列が である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. 2 つのグループの違いに焦点を当てながら、メイン モデルM0はレベルごとに独立したスムーズになります。M1

M0「参照平滑＋差分平滑」の形を許容するように制御できればよい。次に、スムーズな差が直線になった場合、実際にフィッティングすることなくM1、グループ間に非線形差の証拠がないことをすでに知っています。

では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))最初に、AICと同じモデルをパラメーター化して比較したことを例に含めるのを忘れていましたs(x) + s(x, by = gender)(リコールgenderは 0 ～ 1 の数値変数です)。これらのモデルは統計的に同等ですが、平滑化パラメーターはケースによって明らかに異なって推定されるため、AIC は少し異なります。

ああ、そうです。あなたgenderはバイナリなので、数値byもスムーズな違いを構築するための良い考えです。しかし、これは慎重に行ってください。Numericalbyは、中央のスムーズを生成しません。したがって、s(x) + s(x, by = gender)暗黙的に切片列があり、モデルの切片と混同されます。あなたは一緒に行くべきs(x) + s(x, by = gender) - 1です。

于 2016-11-09T12:41:40.150 に答える