0

我有 9 条退化曲线，我想比较它们，并希望得到关于如何最好地做到这一点的建议。我最初的想法是比较非线性回归。我将首先解释这个问题，然后再详细介绍一下实验设计：

我的问题是：

  1. 如何比较我的 9 组之间的退化率？
  2. 我如何确定我的两个主要自变量（有机物类型和田间图）在多大程度上推动了降解率。

我将 3 种类型的有机物（X、Y、Z）放置在 3 个田间地块（A、B、C）中。每个小区放置12个有机质样品（每个小区36个样本，共108个样本）。我知道每个样品的原始有机物 (OM)（作为总值和干物质的百分比）含量。在每周间隔 3 个时间点（T1、T2、T3），我从每个地块中取出每种类型的 4 个样本，并再次测量有机物含量。

因此，对于 9 种组合（AX、AY、AZ、BX、BY、BZ、CX、CY、CZ）中的每一种，我有：在 T0 时对原始有机物进行 12 次测量，在后 3 次中对有机物进行 4 次测量点（T1、T2、T3）。

我希望我提供了足够的信息——如果我还没有，请询​​问。我非常感谢围绕此查询的任何帮助和建议。

谢谢你，圣诞快乐。

安德鲁。

链接到示例日期： https ://docs.google.com/spreadsheets/d/1a5w9BeeogprKAOwHi3WYSW7JF8EtjQOaaG9z-qzDRgw/pub?output=xlsx

4

1 回答 1

1

这里有一些快速让你开始的东西。由于您的时间点很少，我会使用线性模型。我假设 OM 的绝对差异在这里是明智的，即样本以某种有意义的方式归一化。您可能需要使用相对值（在这种情况下甚至可能需要 GLMM？）。

library(data.table)
DT <- fread("Untitled spreadsheet.csv")
setnames(DT, make.names(names(DT)))
DT[, DiffOM := OM.at.collection..g. - Original.OM..T0...g.]
DT <- DT[-1]

library(ggplot2)
p <- ggplot(DT, aes(x = Day.of.collection, y = DiffOM, color = Plot)) +
  geom_point() +
  facet_wrap(~ Sample.type, ncol = 1)
print(p)

情节1

有些人建议只在有大量组可用的情况下才拟合随机效应，但如果结果拟合看起来合理，我通常也相信只有少数组的模型。当然，在这种情况下，您不应该过于相信随机效应的方差估计。或者，您可以将Plot其视为固定效应，但您的模型将需要另外两个参数。然而，通常我们对情节差异不太感兴趣，更喜欢关注治疗效果。YMMV。

library(lmerTest)

fit1 <- lmer(DiffOM ~ Day.of.collection * Sample.type + (1 | Plot), data = DT)
fit2 <- lmer(DiffOM ~ Day.of.collection * Sample.type + (Day.of.collection | Plot), data = DT)
lme4:::anova.merMod(fit1, fit2)
#random slope doesn't really improve the model

fit3 <- lmer(DiffOM ~ Day.of.collection * Sample.type + (Sample.type | Plot), data = DT)
lme4:::anova.merMod(fit1, fit3)
#including the Sample type doesn't either

summary(fit1)
#apparently the interactions are far from significant

fit1a <- lmer(DiffOM ~ Day.of.collection + Sample.type + (1 | Plot), data = DT)
lme4:::anova.merMod(fit1, fit1a)
plot(fit1a)
#seems more or less okay with possibly exception of small degradation
#you could try a variance structure as implemented in package nlme

anova(fit1a)
#Analysis of Variance Table of type III  with  Satterthwaite 
#approximation for degrees of freedom
#                  Sum Sq Mean Sq NumDF DenDF F.value    Pr(>F)    
#Day.of.collection 3909.4  3909.4     1   102 222.145 < 2.2e-16 ***
#Sample.type        452.4   226.2     2   102  12.853 1.051e-05 ***
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

显然，采样前的降解率在样本类型之间是不同的（根据样本类型查看不同的截距），这意味着非线性速率（正如我们所期望的那样）。差异的线性模型意味着恒定的绝对降解率。

summary(fit1a)

newdat <- expand.grid(Day.of.collection = seq(28, 84, by = 1), 
                      Plot = c("A", "B", "C"), 
                      Sample.type = c("X", "Y", "Z"))
newdat$pred <- predict(fit1a, newdata = newdat)
newdat$pred0 <- predict(fit1a, newdata = newdat, re.form = NA)

p +
  geom_line(data = newdat, aes(y = pred, size = "subjects")) +
  geom_line(data = newdat, aes(y = pred0, size = "population", color = NULL)) +
  scale_size_manual(name = "Level",
                    values = c("subjects" = 0.5, "population" = 1.5))

情节2

于 2015-12-18T13:43:13.070 回答