2012-05-03 29 views
1

私は生息地の分布モデルを生成するために生物の視覚トランスセクトデータを分析しようとしています。いったん生物が目撃されると、所与の時間間隔でポイントデータが収集されるので、それらは追跡される。これらの「自己相関」の自己相関のために、私はPirotta et al。のものと同様のGAM-GEEアプローチを利用したいと考えています。 2011年、パッケージ「yags」と「スプライン」を使用しています(http://www.int-res.com/abstracts/meps/v436/p257-272/)。彼らのRスクリプトはここに表示されています(http://www.int-res.com/articles/suppl/m436p257_supp/m436p257_supp1-code.r)。私はこのコードを限定的に成功させ、収束しないモデルの複数の問題を使用しました。gamm4 RパッケージでGAM-GEEを実行していますか?

以下

私のデータの構造である:ここ

> str(dat2) 

'data.frame': 10792 obs. of 4 variables: 

$ dist_slag  : num 26475 26340 25886 25400 24934 ... 
$ Depth   : num -10.1 -10.5 -16.6 -22.2 -29.7 ... 
$ dolphin_presence: int 0 0 0 0 0 0 0 0 0 0 ... 


$ block   : int 1 1 1 1 1 1 1 1 1 1 ... 


> head(dat2) 

    dist_slag Depth dolphin_presence block 
1 26475.47 -10.0934    0  1 
2 26340.47 -10.4870    0  1 
3 25886.33 -16.5752    0  1 
4 25399.88 -22.2474    0  1 



5 24934.29 -29.6797    0  1 
6 24519.90 -26.2370    0  1 

Iは、自己相関しかしながら、各ブロック内

> summary(dat2$block) 
    Min. 1st Qu. Median Mean 3rd Qu. Max. 
    1.00 39.00 76.00 73.52 111.00 148.00 

が存在するグループの数を示す(私のブロック変数の要約であります私はSimon Wood教授のパッケージと関数にもっと精通しているので、パッケージ 'gamm4'を使いたいと思っています。そして、gamm4が最も適切かもしれないと思われます。トランセクトに沿って)なぜ私はgamm4がgammよりも適切だと思うのですか? gammのヘルプでは、それが要因の中に、自己相関のために、次の例を提供します。

## more complicated autocorrelation example - AR errors 
## only within groups defined by `fac' 
e <- rnorm(n,0,sig) 
for (i in 2:n) e[i] <- 0.6*e[i-1]*(fac[i-1]==fac[i]) + e[i] 
y <- f + e 
b <- gamm(y~s(x,k=20),correlation=corAR1(form=~1|fac)) 

この例に続いて、次は私が(出力を調べることによって、しかし、私のデータセットの

b <- gamm4(dolphin_presence~s(dist_slag)+s(Depth),random=(form=~1|block), family=binomial(),data=dat) 

を使用するコードです私は結果をどのように解釈するのか分からない、あるいはグループ内の自己相関が考慮されているとは思わない。

> summary(b$gam) 

Family: binomial 
Link function: logit 

Formula: 
dolphin_presence ~ s(dist_slag) + s(Depth) 

Parametric coefficients: 


      Estimate Std. Error z value Pr(>|z|) 
    (Intercept) -13.968  5.145 -2.715 0.00663 ** 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Approximate significance of smooth terms: 


       edf Ref.df Chi.sq p-value  
s(dist_slag) 4.943 4.943 70.67 6.85e-14 *** 
s(Depth)  6.869 6.869 115.59 < 2e-16 *** 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 



R-sq.(adj) = 0.317glmer.ML score = 10504 Scale est. = 1   n = 10792 
> 

> summary(b$mer) 
Generalized linear mixed model fit by the Laplace approximation 


    AIC BIC logLik deviance 
10514 10551 -5252 10504 
Random effects: 
Groups Name   Variance Std.Dev. 
Xr  s(dist_slag) 1611344 1269.39 
Xr.0 s(Depth)  98622 314.04 
Number of obs: 10792, groups: Xr, 8; Xr.0, 8 



Fixed effects: 
       Estimate Std. Error z value Pr(>|z|) 
X(Intercept)  -13.968  5.145 -2.715 0.00663 ** 
Xs(dist_slag)Fx1 -35.871  33.944 -1.057 0.29063 
Xs(Depth)Fx1  3.971  3.740 1.062 0.28823 


--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Correlation of Fixed Effects: 
      X(Int) X(_)F1 
Xs(dst_s)F1 0.654  
Xs(Dpth)Fx1 -0.030 0.000 
> 

"ブロック"変数の一意の値内で自己相関が実際に考慮されるようにするにはどうすればよいですか? "summary(b $ mer)"の出力を解釈する最も簡単な方法は何ですか?

結果は、同じ変数とパラメータを使用して、通常のgam(パッケージmgcv)と "correlation = ..."という用語が異なる場合、異なるものが発生していることを示します。

しかし、私は相関項(シーズン)のために別の変数を使用する場合、私は同じ出力を得る:私はちょうどそれが正しくため、各値内の相関関係を考慮していることを確認したい

> dat2 <- data.frame(dist_slag = dat$dist_slag, Depth = dat$Depth, dolphin_presence = dat$dolphin_presence, 

+ block = dat$block, season=dat$season) 
> head(dat2) 
     dist_slag Depth dolphin_presence block season 
1 26475.47 -10.0934    0  1  F 
2 26340.47 -10.4870    0  1  F 

3 25886.33 -16.5752    0  1  F 
4 25399.88 -22.2474    0  1  F 
5 24934.29 -29.6797    0  1  F 
6 24519.90 -26.2370    0  1  F 

> summary(dat2$season) 

    F S 
3224 7568 


> b <- gamm4(dolphin_presence~s(dist_slag)+s(Depth),correlation=corAR1(1, form=~1 | season), family=binomial(),data=dat2) 
> summary(b$gam) 

Family: binomial 
Link function: logit 


Formula: 
dolphin_presence ~ s(dist_slag) + s(Depth) 

Parametric coefficients: 
      Estimate Std. Error z value Pr(>|z|) 
    (Intercept) -13.968  5.145 -2.715 0.00663 ** 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 


Approximate significance of smooth terms: 
       edf Ref.df Chi.sq p-value  
s(dist_slag) 4.943 4.943 70.67 6.85e-14 *** 
s(Depth)  6.869 6.869 115.59 < 2e-16 *** 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 


R-sq.(adj) = 0.317glmer.ML score = 10504 Scale est. = 1   n = 10792 
> summary(b$mer) 
Generalized linear mixed model fit by the Laplace approximation 
    AIC BIC logLik deviance 

10514 10551 -5252 10504 
Random effects: 
Groups Name   Variance Std.Dev. 
Xr  s(dist_slag) 1611344 1269.39 
Xr.0 s(Depth)  98622 314.04 
Number of obs: 10792, groups: Xr, 8; Xr.0, 8 


Fixed effects: 
       Estimate Std. Error z value Pr(>|z|) 
X(Intercept)  -13.968  5.145 -2.715 0.00663 ** 
Xs(dist_slag)Fx1 -35.871  33.944 -1.057 0.29063 
Xs(Depth)Fx1  3.971  3.740 1.062 0.28823 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Correlation of Fixed Effects: 
      X(Int) X(_)F1 
Xs(dst_s)F1 0.654  
Xs(Dpth)Fx1 -0.030 0.000 
> 

を"ブロック"変数。自己相関がブロックの単一値内に存在することができると言うモデルを定式化するにはどうすればよいですか?ブロック間の独立性を仮定します。別のノートで

、私はまた、(2よりも多くの変数を持つ)大きなモデルのモデルが完了した後、次の警告メッセージを受信して​​います:

Warning message: 
In mer_finalize(ans) : false convergence (8) 

答えて

3
  • gamm4はいた、lme4の上に構築されていますではないはとは対照的に、correlationのパラメータを許容します(mgcv::gammの下にあるパッケージ)。mgcv::gammは、バイナリデータを扱いますが、一般にgamm4/lme4のようにLaplace/GHQ近似値よりも精度が低いPQLを使用します。 correlation引数が無視されているという警告が表示されないということは残念ですが(残念ですが)引数を使用して簡単な例を試してみると、警告が表示されますが、余分な引数はgamm4のどこかで飲み込まれています)。
  • 自己相関構造(「自己相関はブロックの各単一値内に存在することができますが、ブロック間の独立性を仮定します)」は、nlme(したがってmgcv::gamm)にコーディングされた相関構造です。
  • 私はmcgv::gammを使用し、既知の構造を持ついくつかのシミュレートされたデータで可能な限り試してみることをお勧めします(または上記の補足資料で提供されているデータセットを使用して、あなたの定性的結論を再現できるかどうかを確認してください)。代替アプローチ)。
  • StackOverflowのはいいですが、より多くの混合モデルの専門知識が[email protected]
+0

教授Bolker-で、おそらくそこにあるが、あなたの迅速かつ徹底的な対応をありがとうございました!私はRのエコロジカルモデルに関するあなたのテキストを楽しんだ。私は、私のデータとPirotta et alのデータの両方にmgcv:gammを使用しようと試みるだろう。私は式の相関構造の一部がgamm4で試みたのと同じであると仮定していますか?すなわち、 "相関= corAR1(1、フォーム=〜1 |ブロック)"もう一度ありがとう! – user1367925

+0

私はそれが正しい相関仕様(corar1(form =〜1 | block)で十分かもしれない)だと思います。 –

+0

残念ながら、次のコードを実行すると、後続のエラーが発生します。 > fit1 < - gamm(dolphin_presence_s(深さ、k = 4)+ s(dist_slag、k = 4)、相関= corAR1(形式=〜1 |ブロック)、family = binomial()、data = dat PQL反復の最大数:COEF '20 反復1 反復2 反復3 反復4 反復5 エラー< - corAR1'(' * TMP * '、値=値[parMap [I])。 : 外部関数呼び出し(arg 1)のNA/NaN/Infこれがなぜ発生するかについてのアイデアはありますか?私はまた、より多くのvarsで "convergence error code = 1"のエラーを取得します。 – user1367925

関連する問題