最佳结构 - 3个因素用1纵向

Question

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

“有在的区域被分析的16个树和在B地区有16个。在每个地区，冬季分析8棵树，夏季分析8棵树，但它们不是同一棵树。考虑到在对每棵树直径五个不同深度的starch's知觉。”

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 
...  ...  ...  ...  ... 

我认为在获得适合于与伽玛分布广义线性混合模型（GLMM）在R上我GLMM是一种初学者，我想问问某人我如何在R中表现出这种适合性，以便知道区域，季节和深度因素是否会在响应变量中产生不同的效果。

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

如果不是，将它程序的最佳方式：

，如果我跑这将是正确的吗？如果有人能给我一个方向或至少有一个参考，我非常感激。

谢谢你的帮助@Ben Bolker和@flies。发布的贡献帮助了很多。 然后我会确认是否可以将深度视为定性和交互。这样做：

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;见哈勒尔的回归建模策略 ），除非您的测量结果非常精确，并且存在少量生物变异。完整的固定效应模型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 

难道上面的分析可以执行？鉴于估计参数的数量，您是否应该考虑连续深度和相加模型？