含3个水平和r个相互作用的线性模型中拦截率和斜率的计算

Question

我试图在R中计算线性模型的斜率，我已经将我的数据集聚合到

agg_df <- aggregate(cbind(rate.output, crit.intercept) ~ lvl + treatment, data = d, FUN = mean)

这是Bolker在这个thread中推荐的

这使得我的数据有了一个可复制的例子：

lvl <- as.factor(rep(c(1, 2, 3), 3)) 
treatment <- as.factor(c(rep(c("green"), 3), rep(c("purple"), 3)), 
rep(c("red"), 3)) 
o2 <- c(0.035941608, 0.042206981, 0.023556132, 
     0.016169792, 0.041431159, 0.054221145, 0.007571207, 0.008033468, 0.012353746)  
df <- as.data.frame(cbind(lvl, treatment, o2))

然后，我使用交互作用项运行线性模型：

o2 <- lm(rate.output ~ treatment*lvl, data = agg_df) |>
   summary()

这将返回(我认为在“估计”之后添加了正确的处理方法和级别)：

Coefficients:
                      Estimate <br>
(Intercept)           0.035942  Green Level 1 
treatmentpurple      -0.019772  Purple Level 1
treatmentRed         -0.028370  Red Level 1
lvl2                  0.006265  Green Level 2
lvl3                 -0.012385  Green Level 3
treatmentpurple:lvl2  0.018996  Purple Level 2
treatmentRed:lvl2    -0.005803  Red Level 2
treatmentpurple:lvl3  0.050437  Purple Level 3
treatmentRed:lvl3     0.017168  Red Level 3

然后，我想计算我的不同处理方法和它们的斜率的截距。我认为这是如何做到的一个例子，但给了我错误的结果：

用紫色处理lvl2的斜率：

0.035942+(-0.019772)+0.018996 == 0.035166

有人告诉我，我可以通过这样做来控制我的计算：

o2_purp <- agg_df[agg_df$treatment=="purple",]  
fit_p <- lm(rate.output ~ lvl, data = o2_purp) |> summary() 

Output from fit_p model: 
Coefficients: 
            Estimate Std. Error t value Pr(>|t|) 
(Intercept)  0.035942        NaN     NaN      NaN 
lvl2         0.006265        NaN     NaN      NaN 
lvl3        -0.012385        NaN     NaN      NaN 

这告诉我紫色治疗2级的斜率是0.02526。和上面的方程不一样。我哪里出问题了？如何计算治疗紫和治疗红2级和3级的斜率？

在我计算了斜率后，我想测试不同的治疗水平均值之间的差异，例如，图基后特别分析。

谢谢你抽出时间回答我的问题。

编辑:我这样做的原因是为了比较我在不同的测试介质中接触到的无脊椎动物的正常摄氧量。在测试中，处理被添加为一个梯度，因此我想测试级别之间是否存在差异，这就是为什么我试图计算不同级别的斜率。希望这能澄清一点。

Answer 1

lvl <- as.factor(rep(c(1, 2, 3), 3))
treatment <- as.factor(rep(c("green", "purple","red"), each=3))
o2 <- c(0.035941608, 0.042206981, 0.023556132, 0.016169792, 0.041431159, 0.054221145, 0.007571207, 0.008033468, 0.012353746)
df <- data.frame(lvl, treatment, o2)

mod <- lm(o2 ~ treatment*lvl, data = df)

summary(mod)
Call:
lm(formula = o2 ~ treatment * lvl, data = df)

Residuals:
ALL 9 residuals are 0: no residual degrees of freedom!

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)
(Intercept)           0.035942        NaN     NaN      NaN
treatmentpurple      -0.019772        NaN     NaN      NaN
treatmentred         -0.028370        NaN     NaN      NaN
lvl2                  0.006265        NaN     NaN      NaN
lvl3                 -0.012385        NaN     NaN      NaN
treatmentpurple:lvl2  0.018996        NaN     NaN      NaN
treatmentred:lvl2    -0.005803        NaN     NaN      NaN
treatmentpurple:lvl3  0.050437        NaN     NaN      NaN
treatmentred:lvl3     0.017168        NaN     NaN      NaN

Residual standard error: NaN on 0 degrees of freedom
Multiple R-squared:      1, Adjusted R-squared:    NaN 
F-statistic:   NaN on 8 and 0 DF,  p-value: NA

# To caluclate the slope of Purple treatment lvl2
mod$coefficients["(Intercept)"] + mod$coefficients["treatmentpurple"] + mod$coefficients["treatmentpurple:lvl2"]
0.03516579 
mod2 = lm(o2 ~ treatment*lvl, data = df, subset = df$treatment == 'purple')

Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) : 
contrasts can be applied only to factors with 2 or more level

正如您所看到的，您甚至不能运行第二次回归，因为您需要两个或两个以上的级别才能在回归中使用一个因子变量。在这两种情况下，您都不应该期望得到相同的结果，因为样本大小和规范是不同的。