如何在尊重可变对比度编码的同时使用 anova() 对 lm 和 lmer 对象进行显着性测试?
How can I use anova() for significance testing of lm and lmer objects while respecting variable contrast codings?
我对在 lm 或 lmer 对象上调用的 summary() 的输出中显示的显着性测试结果与在同一对象上调用的 anova() 的输出中显示的结果之间的关系感到困惑。具体来说,我不明白 (a) 为什么对于 df=1 的因素(应该可以比较结果),结果并不总是一致;和 (b) 为什么 summary() 尊重分配给每个因素的对比权重,但 anova() 不尊重。
这里有一个 lm 的例子:
data(iris)
## Apply treatment coding to Species, and fit model
contrasts(iris$Species) <- contr.treatment(length(levels(iris$Species)))
iris.lm.treatment <- lm(Sepal.Length ~ Petal.Length * Species, data=iris)
# check Petal.Length p-value in lm() output
coef(summary(iris.lm.treatment))["Petal.Length","Pr(>|t|)"]
[1] 0.05199902
# check Petal.Length p-value in anova() output
as.data.frame(anova(iris.lm.treatment))["Petal.Length","Pr(>F)"]
[1] 1.244558e-56
## Apply sum coding to Species, and fit model
contrasts(iris$Species) <- contr.sum(length(levels(iris$Species)))/2
iris.lm.sum <- lm(Sepal.Length ~ Petal.Length * Species, data=iris)
# check Petal.Length p-value in lm() output
coef(summary(iris.lm.sum))["Petal.Length","Pr(>|t|)"]
[1] 2.091453e-12
# check Petal.Length p-value in anova() output
as.data.frame(anova(iris.lm.sum))["Petal.Length","Pr(>F)"]
[1] 1.244558e-56
当 Species 的对比编码发生变化时,拟合 lm 中 Petal.Length 的显着性检验会发生变化——这是有道理的,因为模型评估每个因素时正交因素保持恒定为零。但是,方差分析结果中 Petal.Length 的显着性检验无论哪种方式都是相同的,并且与任一 lm 的结果都不匹配。
lmer 的行为(通过 lmerTest 包完成重要性测试)在相关方面令人困惑:
library(lmerTest)
data(ham)
## Apply treatment coding to variables, and fit model
contrasts(ham$Information) <- contr.treatment(length(levels(ham$Information)))
contrasts(ham$Product ) <- contr.treatment(length(levels(ham$Product )))
ham.lmer.treatment <- lmer(Informed.liking ~ Information * Product + (1 | Consumer) + (1 | Consumer:Product), data=ham)
# check Information p-value in lmer() output
coef(summary(ham.lmer.treatment))["Information2","Pr(>|t|)"]
[1] 0.4295516
# check Information p-value in anova() output
as.data.frame(anova(ham.lmer.treatment))["Information","Pr(>F)"]
[1] 0.04885354
## Apply sum coding to variables, and fit model
contrasts(ham$Information) <- contr.sum(length(levels(ham$Information)))/2
contrasts(ham$Product ) <- contr.sum(length(levels(ham$Product )))/2
ham.lmer.sum <- lmer(Informed.liking ~ Information * Product + (1 | Consumer) + (1 | Consumer:Product), data=ham)
# check Information p-value in lmer() output
coef(summary(ham.lmer.sum))["Information1","Pr(>|t|)"]
[1] 0.04885354
# check Information p-value in anova() output
as.data.frame(anova(ham.lmer.sum))["Information","Pr(>F)"]
[1] 0.04885354
在这里,变量编码似乎仍然会影响 summary() 输出中显示的结果,但不会影响 anova() 输出中显示的结果。但是,两个 anova() 结果与使用求和编码时获得的 lmer() 结果匹配。
在我看来,在这两种情况下,anova() 都忽略了使用的变量编码并使用了一些其他变量编码——在 lmer 的情况下,它似乎是求和编码——来评估显着性。我想知道如何执行使用分配的变量编码的统计测试。至少对于 lmer,我可以使用 contestMD() 来完成此操作;例如,
# test Information significance while respecting contrast weights
contestMD(ham.lmer.treatment, as.numeric(names(fixef(ham.lmer.treatment))=="Information2"))[,"Pr(>F)"]
[1] 0.4295516 # matches output from summary(ham.lmer.treatment)
但是,我想不出如何对 lm 进行等效测试(大概是使用 glht,但我想不出正确的函数调用)。所以,我的问题是:
从概念上讲,为什么 anova() 不遵守分配的变量编码? (大概这是所有预期的行为,但原因对我来说是不透明的。)
实际上,在 lm 对象上调用 anova() 时使用什么变量编码?
如何使用 lm 对象执行我想要的重要性测试? (我在上面使用了 df=1 的示例,因为它们可以在模型输出和 anova() 输出之间进行比较,但当然我真正想做的是测试 df>1 的效果。)
我仍然没有回答我的前两个问题,但在回答第三个问题时,我似乎可以通过创建子集模型来获得我想要的结果——每个模型都删除了一个因素——并将每个模型与使用 anova() 的完整模型。对于上面给出的示例 (iris.lm.treatment),我可以执行以下操作。 (在我的示例中,我遇到了麻烦,首先使用明确的数字预测变量重新拟合模型,否则我在使用 anova() 比较模型时会遇到困难。)
# create numeric columns with the same contrast codings as the nominal factor
Species.numeric <- as.data.frame(model.matrix(~ Species, data=iris))
# drop Intercept column
Species.numeric <- Species.numeric[,2:ncol(Species.numeric)]
# rename columns as Species.num1 & Species.num2 and append to iris
names(Species.numeric) <- paste0("Species.num", 1:ncol(Species.numeric))
iris <- cbind(iris, Species.numeric)
# re-fit lm with all numeric predictors
iris.lm.treatment.num <- lm(Sepal.Length ~ Petal.Length * (Species.num1 + Species.num2), data=iris)
# for each factor, create a subset model that has that factor removed
iris.lm.treatment.num.noPetalLength <- update(iris.lm.treatment.num, . ~ . - Petal.Length )
iris.lm.treatment.num.noSpecies <- update(iris.lm.treatment.num, . ~ . - (Species.num1 + Species.num2) )
iris.lm.treatment.num.noInteraction <- update(iris.lm.treatment.num, . ~ . - Petal.Length:(Species.num1 + Species.num2))
# use anova() to compare each subset model to the full model
anova(iris.lm.treatment.num.noPetalLength, iris.lm.treatment.num) # p = .052
anova(iris.lm.treatment.num.noSpecies, iris.lm.treatment.num) # p = 7.611e-06
anova(iris.lm.treatment.num.noInteraction, iris.lm.treatment.num) # p = .1895
花瓣长度的主效应产生的 p 值为 .052,与 iris.lm.treatment 中的结果相匹配。
我对在 lm 或 lmer 对象上调用的 summary() 的输出中显示的显着性测试结果与在同一对象上调用的 anova() 的输出中显示的结果之间的关系感到困惑。具体来说,我不明白 (a) 为什么对于 df=1 的因素(应该可以比较结果),结果并不总是一致;和 (b) 为什么 summary() 尊重分配给每个因素的对比权重,但 anova() 不尊重。
这里有一个 lm 的例子:
data(iris)
## Apply treatment coding to Species, and fit model
contrasts(iris$Species) <- contr.treatment(length(levels(iris$Species)))
iris.lm.treatment <- lm(Sepal.Length ~ Petal.Length * Species, data=iris)
# check Petal.Length p-value in lm() output
coef(summary(iris.lm.treatment))["Petal.Length","Pr(>|t|)"]
[1] 0.05199902
# check Petal.Length p-value in anova() output
as.data.frame(anova(iris.lm.treatment))["Petal.Length","Pr(>F)"]
[1] 1.244558e-56
## Apply sum coding to Species, and fit model
contrasts(iris$Species) <- contr.sum(length(levels(iris$Species)))/2
iris.lm.sum <- lm(Sepal.Length ~ Petal.Length * Species, data=iris)
# check Petal.Length p-value in lm() output
coef(summary(iris.lm.sum))["Petal.Length","Pr(>|t|)"]
[1] 2.091453e-12
# check Petal.Length p-value in anova() output
as.data.frame(anova(iris.lm.sum))["Petal.Length","Pr(>F)"]
[1] 1.244558e-56
当 Species 的对比编码发生变化时,拟合 lm 中 Petal.Length 的显着性检验会发生变化——这是有道理的,因为模型评估每个因素时正交因素保持恒定为零。但是,方差分析结果中 Petal.Length 的显着性检验无论哪种方式都是相同的,并且与任一 lm 的结果都不匹配。
lmer 的行为(通过 lmerTest 包完成重要性测试)在相关方面令人困惑:
library(lmerTest)
data(ham)
## Apply treatment coding to variables, and fit model
contrasts(ham$Information) <- contr.treatment(length(levels(ham$Information)))
contrasts(ham$Product ) <- contr.treatment(length(levels(ham$Product )))
ham.lmer.treatment <- lmer(Informed.liking ~ Information * Product + (1 | Consumer) + (1 | Consumer:Product), data=ham)
# check Information p-value in lmer() output
coef(summary(ham.lmer.treatment))["Information2","Pr(>|t|)"]
[1] 0.4295516
# check Information p-value in anova() output
as.data.frame(anova(ham.lmer.treatment))["Information","Pr(>F)"]
[1] 0.04885354
## Apply sum coding to variables, and fit model
contrasts(ham$Information) <- contr.sum(length(levels(ham$Information)))/2
contrasts(ham$Product ) <- contr.sum(length(levels(ham$Product )))/2
ham.lmer.sum <- lmer(Informed.liking ~ Information * Product + (1 | Consumer) + (1 | Consumer:Product), data=ham)
# check Information p-value in lmer() output
coef(summary(ham.lmer.sum))["Information1","Pr(>|t|)"]
[1] 0.04885354
# check Information p-value in anova() output
as.data.frame(anova(ham.lmer.sum))["Information","Pr(>F)"]
[1] 0.04885354
在这里,变量编码似乎仍然会影响 summary() 输出中显示的结果,但不会影响 anova() 输出中显示的结果。但是,两个 anova() 结果与使用求和编码时获得的 lmer() 结果匹配。
在我看来,在这两种情况下,anova() 都忽略了使用的变量编码并使用了一些其他变量编码——在 lmer 的情况下,它似乎是求和编码——来评估显着性。我想知道如何执行使用分配的变量编码的统计测试。至少对于 lmer,我可以使用 contestMD() 来完成此操作;例如,
# test Information significance while respecting contrast weights
contestMD(ham.lmer.treatment, as.numeric(names(fixef(ham.lmer.treatment))=="Information2"))[,"Pr(>F)"]
[1] 0.4295516 # matches output from summary(ham.lmer.treatment)
但是,我想不出如何对 lm 进行等效测试(大概是使用 glht,但我想不出正确的函数调用)。所以,我的问题是:
从概念上讲,为什么 anova() 不遵守分配的变量编码? (大概这是所有预期的行为,但原因对我来说是不透明的。)
实际上,在 lm 对象上调用 anova() 时使用什么变量编码?
如何使用 lm 对象执行我想要的重要性测试? (我在上面使用了 df=1 的示例,因为它们可以在模型输出和 anova() 输出之间进行比较,但当然我真正想做的是测试 df>1 的效果。)
我仍然没有回答我的前两个问题,但在回答第三个问题时,我似乎可以通过创建子集模型来获得我想要的结果——每个模型都删除了一个因素——并将每个模型与使用 anova() 的完整模型。对于上面给出的示例 (iris.lm.treatment),我可以执行以下操作。 (在我的示例中,我遇到了麻烦,首先使用明确的数字预测变量重新拟合模型,否则我在使用 anova() 比较模型时会遇到困难。)
# create numeric columns with the same contrast codings as the nominal factor
Species.numeric <- as.data.frame(model.matrix(~ Species, data=iris))
# drop Intercept column
Species.numeric <- Species.numeric[,2:ncol(Species.numeric)]
# rename columns as Species.num1 & Species.num2 and append to iris
names(Species.numeric) <- paste0("Species.num", 1:ncol(Species.numeric))
iris <- cbind(iris, Species.numeric)
# re-fit lm with all numeric predictors
iris.lm.treatment.num <- lm(Sepal.Length ~ Petal.Length * (Species.num1 + Species.num2), data=iris)
# for each factor, create a subset model that has that factor removed
iris.lm.treatment.num.noPetalLength <- update(iris.lm.treatment.num, . ~ . - Petal.Length )
iris.lm.treatment.num.noSpecies <- update(iris.lm.treatment.num, . ~ . - (Species.num1 + Species.num2) )
iris.lm.treatment.num.noInteraction <- update(iris.lm.treatment.num, . ~ . - Petal.Length:(Species.num1 + Species.num2))
# use anova() to compare each subset model to the full model
anova(iris.lm.treatment.num.noPetalLength, iris.lm.treatment.num) # p = .052
anova(iris.lm.treatment.num.noSpecies, iris.lm.treatment.num) # p = 7.611e-06
anova(iris.lm.treatment.num.noInteraction, iris.lm.treatment.num) # p = .1895
花瓣长度的主效应产生的 p 值为 .052,与 iris.lm.treatment 中的结果相匹配。