所以当你看this教程的时候，R并没有给予你平均值，就像你看SPSS的时候一样。也许你想尝试一个不同的库来给你更多的信息？你也可以像this教程一样手动操作。他们只是简单的得到每组的平均值：

PlantGrowth %>%
  group_by(group) %>%
  get_summary_stats(weight, type = "mean_sd")

请注意，您的分类变量需要是二进制的，以单独识别每个类别的影响。（如果需要，创建虚拟变量）

我生成了一个随机数据集来增加样本量，作为执行ANCOVA来生成可估计的效应和均值的说明性示例。（checking model assumptions, choosing the best model, testing independence between covariates and factor variable）。例如，我们需要验证协变量（Control1,2）和处理（ind.v）是相互独立的。因为仅当协变量和治疗效应对应答相互独立时，将协变量项纳入模型才有意义变量。另外，另一点是决定哪个模型是最佳模型。我们需要使用完整模型，还是从完整模型中删除一些没有任何意义的项，将它们包含到模型中。另一点是，看看?aov，当数据平衡且处理具有相同的重复次数时，函数是合适的。否则，函数如lm()将更合适。我假设问题所有者提供的完整数据集是balanced，因此我将考虑aov()作为ANOVA。此外，我还想显示ANOVA基础模型（无协变量）和问题所有者指定的完整模型，以比较两个模型估计的平均值和治疗之间的比较。希望这些能有所帮助。

set.seed(123)
x = rnorm(160)

# to generate data in ranges between 0-1.
min.x = min(x)
max.x = max(x)

dep.variable = (x - min.x)/(max.x - min.x)
ind.v = rep(1:4, each=40)
Control1=sample(1:2,size = 160,replace = T)
Control2=sample(1:7,size = 160,replace = T)

# make data frame
df<-data.frame(dep.variable,ind.v,Control1,Control2)

# convert ind.v to factor
df$ind.v <- factor(df$ind.v)

# center the covariates and save centered values in new columns
df$Ctrl1s<-scale(df$Control1, center = T, scale = F)
df$Ctrl2s<-scale(df$Control2, center = T, scale = F)

# run anova with and without covariates
aov.noCov <- aov (dep.variable ~ ind.v , data = df)
aov.Cov <- aov (dep.variable ~ ind.v*Ctrl1s*Ctrl2s, data = df)

#partitioning the sums of squares from the linear model using ANOVA here 
#We use type II sums of squares to address the presence of covariates

library(car)
Anova(aov.noCov, type="II") # anova for model without covariates
                  
Response: dep.variable
          Sum Sq  Df F value Pr(>F)
ind.v     0.0247   3   0.188 0.9045
Residuals 6.8224 156 

Anova(aov.Cov, type="II") # anova for model with covariates

Anova Table (Type II tests)
Response: dep.variable
                    Sum Sq  Df F value Pr(>F)
ind.v               0.0254   3  0.1850 0.9064
Ctrl1s              0.0275   1  0.5995 0.4400
Ctrl2s              0.0602   1  1.3138 0.2536
ind.v:Ctrl1s        0.0592   3  0.4302 0.7317
ind.v:Ctrl2s        0.0511   3  0.3714 0.7738
Ctrl1s:Ctrl2s       0.0014   1  0.0315 0.8594
ind.v:Ctrl1s:Ctrl2s 0.0095   3  0.0691 0.9763
Residuals           6.6005 144

为了在控制和不控制协变量的情况下获得调整后的平均值，我们使用来自包lsmeans的lsmeans()。

library(lsmeans)
lsmc1<-lsmeans(aov.noCov, ~ ind.v) # get adjusted means with no covariate
lsmc2<-lsmeans(aov.Cov, ~ ind.v) # lsmeans controlling for covariates

lsmc1
ind.v lsmean     SE  df lower.CL upper.CL
 1      0.524 0.0331 156    0.458    0.589
 2      0.512 0.0331 156    0.447    0.577
 3      0.515 0.0331 156    0.450    0.581
 4      0.490 0.0331 156    0.425    0.555

lsmc2
 ind.v lsmean     SE  df lower.CL upper.CL
 1      0.512 0.0375 144    0.438    0.586
 2      0.513 0.0361 144    0.441    0.584
 3      0.503 0.0377 144    0.429    0.578
 4      0.492 0.0346 144    0.424    0.561

两个模型的估计均值不同，因为我们在第二个模型中包含协变量。最后，我们定义事后比较以进行成对比较。我们在这里使用Scheffé或Bonferroni程序，因为Tukey’s方法不适用于ANCOVA（Kutner et al. 2005）。

summary(contrast(lsmc2, method="pairwise", adjust="Scheffe"),
        infer=c(T,T))
        
 contrast         estimate     SE  df lower.CL upper.CL t.ratio p.value
 ind.v1 - ind.v2 -0.000195 0.0521 144   -0.147    0.147  -0.004  1.0000
 ind.v1 - ind.v3  0.008977 0.0531 144   -0.141    0.159   0.169  0.9987
 ind.v1 - ind.v4  0.020152 0.0510 144   -0.124    0.164   0.395  0.9843
 ind.v2 - ind.v3  0.009172 0.0522 144   -0.139    0.157   0.176  0.9986
 ind.v2 - ind.v4  0.020347 0.0500 144   -0.121    0.162   0.407  0.9829
 ind.v3 - ind.v4  0.011175 0.0512 144   -0.134    0.156   0.218  0.9973

在这里查看p值，我们观察到配对比较均不显著，因为它是随机数据集。这只是尝试显示存在协变量时ANCOVA的步骤。

aov()是lm()的 Package 器，专门用于提供方差分析所需的常见输出。您可以像使用lm()输出一样使用coef()获得系数，或者直接使用lm()：

x = data.frame(
  y = c(
    rnorm(10, 5, 2),
    rnorm(10, 8, 2)
  ),
  g = rep(letters[1:2], each = 10)
)

aov(y ~ g, data = x)

Call:
   aov(formula = y ~ g, data = x)

Terms:
                       g Residuals
Sum of Squares  33.06216  81.16610
Deg. of Freedom        1        18

Residual standard error: 2.123494
Estimated effects may be unbalanced

使用coef()获得组平均值：
一个二个一个一个
或者直接调用lm()：

lm(y ~ g, data = x)

Call:
lm(formula = y ~ g, data = x)

Coefficients:
(Intercept)           gb  
      5.785        2.571

请注意，第一组'a'的均值在(Intercept)项中，'b'的效应是相对于此截距的差值。实际组均值是'a'的截距与'b'的两个系数之和。
对于两个以上的因子，效应会受到混淆，并且无法通过模型输出获得组平均值（不是aov就是lm）.这样做的原因，我怀疑，这些边际平均值只有在因素之间没有相互作用的情况下才有意义。有一个专门的软件包来处理这个问题，叫做emmeans。我们必须使用模拟数据集，因为最小示例对于水平数的观测值仍然太少：

# sim parameters
n = 10
g_effects = c(
  'a' = 5,
  'b' = 8
)
ctrl_effects = c(
  'c' = 2,
  'd' = 4
)

effects = outer(g_effects, ctrl_effects, `+`)
effects

c  d
a  7  9
b 10 12

# generate data
X = data.frame(
  y = rep(effects, each = n) + rnorm(n * 4),
  g = rep(rownames(effects), each = n),
  ctrl = rep(colnames(effects), each = n * 2)
)

X

y g ctrl
1   6.519858 a    c
2   6.704983 a    c
3   7.542411 a    c
4   3.915727 a    c
5   7.564559 a    c
6   7.840707 a    c
7   5.137290 a    c
8   7.696412 a    c
9   8.300361 a    c
10  5.942079 a    c
11  9.682047 b    c
12 10.551709 b    c
13  9.208265 b    c
14  9.249308 b    c
15  9.994033 b    c
16 10.969805 b    c
17 10.547202 b    c
18  8.613042 b    c
19 10.893566 b    c
20 10.513579 b    c
21  9.447953 a    d
22 10.432833 a    d
23 10.443755 a    d
24  8.896230 a    d
25  8.385421 a    d
26  8.250635 a    d
27  8.673966 a    d
28  9.684629 a    d
29 10.183822 a    d
30  8.721494 a    d
31 11.734383 b    d
32 10.844294 b    d
33 11.906789 b    d
34 11.890356 b    d
35 11.941118 b    d
36 11.708217 b    d
37 14.436147 b    d
38 12.601567 b    d
39 10.799885 b    d
40 14.598940 b    d

现在，我们拟合方差分析并求出边际平均值：

library('emmeans')

m = aov(y ~ g + ctrl, data = X)
emmeans(m, spec = 'g')

g emmean    SE df lower.CL upper.CL
 a   8.01 0.247 37     7.51     8.51
 b  11.13 0.247 37    10.63    11.63

Results are averaged over the levels of: ctrl 
Confidence level used: 0.95