0

我正在测试两个栖息地（入侵和非入侵）和三种不同的柱头类型（湿、干和半干）之间花粉沉积的差异。这是一种社区方法，每个站点的样本和物种数量不平衡，数据的非正态分布，最终具有嵌套的随机结构，符合伽马误差分布，以处理伪复制和非独立性。

为了找出最佳模型，我使用了似然比检验，这表明具有固定效应相互作用的模型更适合：

 > 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

尝试绘制数据。看起来被入侵的站点总是高于/低于未入侵的站点，还是没有明确的模式？如果图片不是很清晰，我推荐post-hoc测试。我刚刚使用此代码进行了类似的事后分析：

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 回答