r-GLMM的最佳结构-具有1个纵向的3个因子

Question

我想对以下问题进行重复测量分析/纵向数据：

"在a区域和B区域分别分析了16棵树。在每个区域，冬季和夏季分别分析了8棵树和8棵树，但它们不是同一棵树。考虑到淀粉对每棵树直径在五个不同深度的感知。">

tree    Region  season  depth   starch
1       A       W       1       0.07
1       A       W       2       0.10
1       A       W       3       0.13
1       A       W       4       0.16
1       A       W       5       0.11
2       A       W       1       0.07
2       A       W       2       0.10
2       A       W       3       0.13
2       A       W       4       0.16
2       A       W       5       0.11
...     ...     ...     ...     ...
17      B       S       1       0.06
...     ...     ...     ...     ...

我认为，在R中具有伽玛分布的广义线性混合模型(GLMM)上进行拟合。作为GLMM的初学者，我想问某人我如何在R中进行拟合，以了解区域、季节和深度因子是否会对响应变量产生不同的影响。

如果我运行，这将是正确的

require(lme4)
Mod1=glmer(starch~Region*season*depth+(1|tree),data=data,family="gamma")
summary(Mod1)  
?

如果没有，最好的处理方法是什么？如果有人能给我一个方向或至少一个推荐信，我将不胜感激。

感谢你的帮助@Ben Bolker和@flies。张贴的稿件帮助很大。然后，我将确认是否有可能将深度视为定性和Region * stage * depth相互作用。这样做：

data = within (data, {
Region = factor (Region) 
season = factor (season) 
depth = factor (depth)  })
require (lme4) 
Mod1 = glmer (starch~Region*season*depth+(1|tree),data=data,family=Gamma(link="log"))
summary (Mod1)
library (car)
Anova(mod1)

获得以下结果：

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation)
glmerMod]
Family:Gamma(log)
Formula: starch ~Region*season*depth+(1|tree)
Date: data
AIC     BIC    logLik  deviance df.resid
-1358.4 -1290.7 701.2   -1402.4   138
Scaled residuals:
Min 1Q Median 3Q Max
-2.3398 -0.6699 -0.1065 0.6683 3.2020
Random effects:
Groups Name Variance Std.Dev.
tree (Intercept) 7.171e-05 0.008468
Residual 6.020e-04 0.024536
Number of obs: 160, groups: tree, 32
Fixed effects:
Estimate     Std. Error    t value    Pr (> | z |)
(Intercept) -2.593064      0.009621      -269.51     <2e-16 ***
RegionRP     0.260453      0.013607        19.14    <2e-16 ***
seasonV     -0.193693      0.013607       -14.23   <2e-16 ***
depth2       0.409813      0.011894        34.46    <2e-16 ***
depth3       0.594269      0.011893        49.97    <2e-16 ***
depth4       0.779051      0.011893        65.50    <2e-16 ***
depth5       0.432146      0.011893        36.34    <2e-16 ***
RegionRP:seasonV    0.088320    0.019243    4.59    4.44e-06 ***
localRP:depth2     -0.065211    0.016820    -3.88   0.000106 ***
localRP:depth3     -0.130185    0.016819    -7.74   9.92e-15 ***
localRP:depth4     -0.190743    0.016820   -11.34   <2e-16 ***
localRP:depth5     -0.067266    0.016820    -4.00   6.35e-05 ***
seasonV:depth2      0.031624    0.016821     1.88   0.060103.
seasonV:depth3      0.139424    0.016820     8.29   <2e-16 ***
seasonV:depth4      0.147717    0.016820     8.78   <2e-16 ***
seasonV:depth5      0.107589    0.016820     6.40   1.59e-10 ***
RegionRP:seasonV:depth2 -0.018490  0.023787 -0.78   0.436970
RegionRP:seasonV:depth3 -0.113810  0.023786 -4.78   1.71e-06 ***
RegionRP:seasonV:depth4 -0.112337  0.023787 -4.72   2.33e-06 ***
RegionRP:seasonV:depth5 -0.121932  0.023787 -5.13   2.96e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '' 1
Analysis of Deviance Table (Type II Wald chisquare tests)
Response: starch
Chisq      Df  Pr (> Chisq)
Region                 872.9486   1   <2.2e-16 ***
season                 282.9125   1   <2.2e-16 ***
depth                16726.2395   4   <2.2e-16 ***
Region:season            1.5641   1    0.2111
Region:depth           521.4171   4   <2.2e-16 ***
season:depth            85.5213   4   <2.2e-16 ***
Region:season:depth     49.1586   4    5.41e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '' 1

是否可以进行上述分析？给定估计参数的数量，是否应该考虑连续深度和加法模型？

Answer 1

• 您需要family="Gamma"(引号是可选的，但必须是大写G)。我经常建议(1)使用对数链路(family=Gamma(link="log"))而不是默认的反向链路，和/或(2)在这种情况下使用对数线性混合模型(lmer(log(starch)~...))。两者在数值上都比默认的Gamma模型更稳定，参数也更容易解释
• 与上面的一些评论相反，默认情况下，glmer将把深度等数字变量视为连续预测器，这意味着它将假设(对数)-线性关系，并且只适合单个参数。你将无法检测到"哪个深度不同"，但你将能够(希望)检测到淀粉随着深度的变化而不断变化
• 如果你总共有32棵树，那么尝试拟合一个参数超过3个或最多4个的模型是有风险的(低功耗)(最大参数~n/10；请参阅Harrell的回归建模策略)，除非你的测量非常精确，并且有少量的生物变化。完全固定效果模型Region*season*depth为您提供了8个参数(2个区域X 2个季节X(截距、斜率)：不是上面讨论的20个，因为深度是连续的)
• 独立分析每个地区和季节的深度几乎与拟合具有所有交互作用的模型相同；由于每个地区/季节的组合中只有8棵树，因此很难建立可靠的模型
• 如果你愿意放弃你的互动，并拟合只有四个参数的加性模型Region + season + depth(截距+区域效应(假设所有季节和深度的常数)+季节效应(假设常数…)+深度效应(假设常量…)，再加上估计树木变异性的随机效应参数(在考虑固定效应后)-仍然有点太复杂，但也许你可以逃脱惩罚
• 不要忘记绘制数据，并检查模型的图形诊断

Answer 2

感谢您对@BenB和@flies的帮助。张贴的稿件帮助很大。

然后，我将确认是否有可能将深度视为定性和Region * stage * depth相互作用。这样做：

data = within (data, {
Region = factor (Region) 
season = factor (season) 
depth = factor (depth)  })
require (lme4) 
Mod1 = glmer (starch~Region*season*depth+(1|tree),data=data,family=Gamma(link="log"))
summary (Mod1)
library (car)
Anova(mod1)

获得以下结果：

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation)
glmerMod]
Family:Gamma(log)
Formula: starch ~Region*season*depth+(1|tree)
Date: data
AIC     BIC    logLik  deviance df.resid
-1358.4 -1290.7 701.2   -1402.4   138
Scaled residuals:
Min 1Q Median 3Q Max
-2.3398 -0.6699 -0.1065 0.6683 3.2020
Random effects:
Groups Name Variance Std.Dev.
tree (Intercept) 7.171e-05 0.008468
Residual 6.020e-04 0.024536
Number of obs: 160, groups: tree, 32
Fixed effects:
Estimate     Std. Error    t value    Pr (> | z |)
(Intercept) -2.593064      0.009621      -269.51     <2e-16 ***
RegionRP     0.260453      0.013607        19.14    <2e-16 ***
seasonV     -0.193693      0.013607       -14.23   <2e-16 ***
depth2       0.409813      0.011894        34.46    <2e-16 ***
depth3       0.594269      0.011893        49.97    <2e-16 ***
depth4       0.779051      0.011893        65.50    <2e-16 ***
depth5       0.432146      0.011893        36.34    <2e-16 ***
RegionRP:seasonV    0.088320    0.019243    4.59    4.44e-06 ***
localRP:depth2     -0.065211    0.016820    -3.88   0.000106 ***
localRP:depth3     -0.130185    0.016819    -7.74   9.92e-15 ***
localRP:depth4     -0.190743    0.016820   -11.34   <2e-16 ***
localRP:depth5     -0.067266    0.016820    -4.00   6.35e-05 ***
seasonV:depth2      0.031624    0.016821     1.88   0.060103.
seasonV:depth3      0.139424    0.016820     8.29   <2e-16 ***
seasonV:depth4      0.147717    0.016820     8.78   <2e-16 ***
seasonV:depth5      0.107589    0.016820     6.40   1.59e-10 ***
RegionRP:seasonV:depth2 -0.018490  0.023787 -0.78   0.436970
RegionRP:seasonV:depth3 -0.113810  0.023786 -4.78   1.71e-06 ***
RegionRP:seasonV:depth4 -0.112337  0.023787 -4.72   2.33e-06 ***
RegionRP:seasonV:depth5 -0.121932  0.023787 -5.13   2.96e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '' 1
Analysis of Deviance Table (Type II Wald chisquare tests)
Response: starch
Chisq      Df  Pr (> Chisq)
Region                 872.9486   1   <2.2e-16 ***
season                 282.9125   1   <2.2e-16 ***
depth                16726.2395   4   <2.2e-16 ***
Region:season            1.5641   1    0.2111
Region:depth           521.4171   4   <2.2e-16 ***
season:depth            85.5213   4   <2.2e-16 ***
Region:season:depth     49.1586   4    5.41e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '' 1

是否可以进行上述分析？给定估计参数的数量，是否应该考虑连续深度和加法模型？