0

私は、2 つの生息地 (侵襲されたものと侵略されていないもの) と 3 つの異なる柱頭タイプ (ウェット、ドライ、セミドライ) の間の花粉の沈着の違いをテストしています。これは、疑似複製と非独立性に対処するために、サイトごとのサンプル数と種の数が不均衡であり、データの非正規分布がガンマ誤差分布に適合するネストされたランダム構造を持つことになったコミュニティ アプローチです。

最適なモデルを見つけるために、尤度比検定を使用しました。これは、固定効果の相互作用を持つモデルがよりよく適合することを示しています。

 > m1b<-glmer(nb~habitat*stigmatype+(1|sitecode/stigmaspecies), family=Gamma(link=log))
 > m2b<-glmer(nb~habitat+stigmatype+(1|sitecode/stigmaspecies), family=Gamma(link=log))
 > anova(m1b,m2b)
Data: 
Models:
m2b: nb ~ habitat + stigmatype + (1 | sitecode/stigmaspecies)
m1b: nb ~ habitat * stigmatype + (1 | sitecode/stigmaspecies)
Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
m2b  7 3032.8 3061.3 -1509.4   3018.8                           
m1b  9 3030.1 3066.7 -1506.0   3012.1 6.6672      2    0.03566 *

そこから、固定項の p 値を解釈する方法について少し混乱しています。以下の出力を見て、生息地と柱頭のタイプの p 値を交互作用項の独立した結果として解釈できますか? 言い換えれば、変数の生息地は、侵略されていない生息地と侵略された生息地 (切片) が異なるように、それ自体に大きな影響を与えていると言えますか? スティグマ型も同じ思考？または、相互作用がわずかに重要であるため、固定値を個別に解釈できなくなりましたか? そして、事後テストだけが、実際の違いがどこにあるのかを教えてくれるでしょうか?

m1b<-glmer(nb~habitat*stigmatype+(1|sitecode/stigmaspecies), family=Gamma(link=log))
summary(m1b)

Generalized linear mixed model fit by maximum likelihood ['glmerMod']
Family: Gamma ( log )
Formula: nb ~ habitat * stigmatype + (1 | sitecode/stigmaspecies) 

  AIC       BIC    logLik  deviance 
 3030.101  3066.737 -1506.050  3012.101 

 Random effects:
 Groups                 Name        Variance  Std.Dev. 
 stigmaspecies:sitecode (Intercept) 5.209e+00 2.2822436
 sitecode               (Intercept) 2.498e-07 0.0004998
 Residual                           2.070e+00 1.4388273
 Number of obs: 433, groups: stigmaspecies:sitecode, 109; sitecode, 20

 Fixed effects:
                                 Estimate Std. Error t value Pr(>|z|)    
 (Intercept)                            2.3824     0.4080   5.839 5.26e-09 ***
 habitatnon-invaded                    -1.8270     0.6425  -2.843  0.00446 ** 
 stigmatypesemidry                     -1.7531     0.7573  -2.315  0.02061 *  
 stigmatypewet                         -1.7210     0.8944  -1.924  0.05434 .  
 habitatnon-invaded:stigmatypesemidry   2.0774     1.1440   1.816  0.06938 
 habitatnon-invaded:stigmatypewet       1.3120     1.4741   0.890  0.37346

ありがとうございます！

4

1 に答える 1

0

データをグラフ化してみてください。侵略されたサイトは、侵略されていないサイトより常に高い/低いように見えますか、それとも明確なパターンはありませんか? 全体像がはっきりしない場合は、事後テストをお勧めします。このコードを使用して、同様の事後分析を行ったところです。

levels(data$habitat)
levels(data$stigmatype)#I believe 1 is the first value that is returned when you run levels(), 0 is the next, and -1 the last

testInteractions(m1b, custom=list(habitat='invaded', stigmatype=c(1,0,-1), adjustment="none"))
testInteractions(m1b, custom=list(habitat='non-invaded', stigmatype=c(1,0,-1), adjustment="none"))

testInteractions(m1b, custom=list(stigmatype='wet', habitat=c(1,-1), adjustment="none"))
testInteractions(m1b, custom=list(stigmatype='dry', habitat=c(1,-1), adjustment="none"))
testInteractions(m1b, custom=list(stigmatype='semidry', habitat=c(1,-1), adjustment="none"))
于 2014-09-19T04:22:45.293 に答える