8

私は、gam を使用して時系列データの季節性をモデル化することに大きな成功を収めました。私の最新のモデルは、季節的な変化に加えて、週ごとのパターンを明確に示しています。1 週間のパターン自体は 1 年を通して非常に安定していますが、その振幅は季節によっても異なります。したがって、理想的には、データを次のようにモデル化したいと考えています。

y ~ f(day in year) + g(day in year) * h(day in week)

ここでf、gとhは の循環平滑関数です。mgcv

gam(
  y ~ s(day_in_year, k=52, bs='cc') 
  + s(day_in_year, k=52, bs='cc'):s(day_in_week, k=5, bs='cc')
  , knots=list(
    day_in_year=c(0, 356)
    , day_in_week=c(0,7)
  )
  , data = data
)

残念ながら、これは機能せず、エラーがスローされますNA/NaN argument。whichを使用してみte(day_in_year, day_in_week, k=c(52, 5), bs='cc')ましたが、モデルが利用可能な年数が少ない特定の平日に当たる休日に過適合するため、自由度が多すぎます。

私がやろうとしている方法でモデルを指定することは可能ですか?

4

3 に答える 3

6

うわー、かなり古い質問です！

インタラクションについて

1 週間のパターン自体は 1 年を通して非常に安定していますが、その振幅は季節によっても異なります。

テンソル積スプライン基底の使用はte相互作用の正しい方法ですが、より適切なコンストラクターはti. あなたはそれteがあなたに多くのパラメータを返すと言った。もちろん。k = 52最初のマージンとk = 52 番目のマージンがあり、このテンソル項の係数が得られます52 * 5 - 1。しかし、これは相互作用が作成される方法にすぎません。

mgcvGAM 式では、:or*はパラメトリック項間の相互作用に対してのみ有効であることに注意してください。スムーズ間の相互作用はteまたはによって処理されtiます。

これがあなたが望んでいるものではない場合、あなたは「製品」に何を期待しますか? 2 つの周辺計画行列のアダマール積ですか? では、その意味は何ですか？ところで、アダマール積には同じ次元の行列が 2 つ必要です。ただし、2 つの余白の列数は同じではありません。

なぜ私が行列について話し続けるのか理解できない場合は、Simon の 2006 年の本を読む必要があります。GAM の推定については、現在はかなり時代遅れであると説明されていますが、GAM の構築/設定 (設計行列やペナルティ行列など) については、章で説明されています。 4 10年経っても全然変わらない。

わかりました、もう 1 つヒントをあげましょう。/計画行列の構築に使用される行単位のクロネッカー積は新しい発明ではありません。teti

平滑項s(x)は因子変数gによく似ています。単一の変数のように見えますが、多くの列を持つ計画行列として構築されます。これはgダミー行列でf(x)あり、基底行列です。したがって、2 つの滑らかな関数間の相互作用は、2 つの因子間の相互作用が構築されるのと同じ方法で構築されます。

g15 水準の因子と 10 水準の別の因子がある場合g2、それらの限界計画行列 (対比後) は 4 列と 9 列であり、交互作用g1:g2は 36 列になります。このような計画行列は、 と の計画行列の行単位のクロネッカー積g1ですg2。

オーバーフィッティングについて

おっしゃる通り、データは数年分しかありません。おそらく 2 年か 3 年ですか？k = 52その場合、モデルはforを使用して過度にパラメータ化されていday_in_yearます。たとえば に減らしてみてくださいk = 30。

オーバーフィッティングがまだ明らかな場合は、対処する方法がいくつかあります。

方法 1:滑らかさの選択に GCV を使用しています。試してみてくださいmethod = "REML"。GCV は常にデータをオーバーフィットする傾向があります。

方法 2: GCV にとどまり、手動で平滑化パラメーターを膨らませて、より重いペナルティを課します。gammaパラメータ ofgamはここで役立ちます。gamma = 1.5たとえば、試してみてください。

コードのタイプミス？

ノットの位置はday_in_year = c(0, 365)?

于 2017-05-17T07:41:25.240 に答える
0

私の質問は数年前のものであり、最終的に私にとって最適なソリューションを追加したいと考えています。

まず、私の質問で説明されているように、mgcv を使用してモデルを適合させることは不可能です。

しばらくの間、私は 2 段階のプロセスを使用しました。

model1 = gam(
  y ~ s(day_in_year, k=52, bs='cc') 
  + s(day_in_year, k=52, bs='cc')
  + as.factor(day_in_week)
  , knots=list(
    day_in_year=c(0, 366)
    , day_in_week=c(0,7)
  )
  , data = data
)
# get_weekday_offset gets the coefficients for each weekday and normalizes them to have mean 0
data$weekday_offset = get_weekday_offset(model1)[data$day_in_week]

model2 = gam(
  y ~ s(day_in_year, k=52, bs='cc') 
  + s(day_in_year, k=52, bs='cc')
  + s(day_in_year, k=52, bs='cc', by=weekday_offset)
  + as.factor(day_in_week)
  , knots=list(
    day_in_year=c(0, 366)
    , day_in_week=c(0,7)
  )
  , data = data
)

2 番目のモデルの形式は「y ~ f(年内の日) + g(年内の日) * h(曜日)」ですが、h は固定関数であり、2 番目のモデルのデータには当てはまりません。に最適model1です。これの明らかな欠点は、mgcv を 2 回実行する必要があり、最初の実行からのほぼすべての情報 (平日の係数を除く) が破棄されることです。

最終的にこのモデルを断念したのは、振幅の変化だけでなく、形状にも年々変化があったため、最終的にこのような外観に落ち着いたからです。

gam(
  y ~ s(day_in_year, k=52, bs='cc') 
  + s(day_in_year, k=10, bs='cc', by=as.factor(day_in_week)
  + as.factor(day_in_week)
  , knots=list(
    day_in_year=c(0, 366)
    , day_in_week=c(0,7)
  )
  , data = data
)

週次パターンの自由度を最大にしてデータを正当化したので、平日の を削除しs、平日と年間通算日の交互作用効果を追加しました ( s(day_in_year, k=12, bs='cc', by=as.factor(day_in_week))。これにより、別の季節性曲線が推定されます。ただし、季節項よりも自由度が大幅に低いため、項s(day_in_year, k=52, bs='cc')を試すときに発生した問題には遭遇しませんte。

于 2020-06-09T12:59:15.907 に答える