如何使用R比较三个(或更多)组的回归系数？

Question

有时，您的研究可能会预测，回归系数的大小可能会因组而异。例如，您可能认为身高预测体重的回归系数在三个年龄组(年轻、中年、老年人)中会有所不同。下面，我们有一个数据文件，其中包括3个虚构的年轻人，3个虚构的中年人和3个虚构的老年人，以及他们的身高和体重。变量age表示年龄组，编码为1(年轻人)、2(中年)和3(老年人)。

那么，我如何使用R比较三个(或更多)组的回归系数(主要是斜率)？

示例数据：

age height weight 
 1  56 140   
 1  60 155   
 1  64 143     
 2  56 117   
 2  60 125   
 2  64 133      
 3  74 245   
 3  75 241   
 3  82 269   

Answer 1

为了确定回归系数是否“在三个年龄组之间不同”，我们可以在R中使用anova函数。例如，使用问题中的数据，并在末尾的注释中重复显示：

fm1 <- lm(weight ~ height, DF)
fm3 <- lm(weight ~ age/(height - 1), DF)

给出以下在2.7%水平上显着的结果，因此我们会得出结论，如果我们使用5%的截止值，则各组的回归系数存在差异，但如果我们使用1%的截止值，则不存在差异。

> anova(fm1, fm3)
 Analysis of Variance Table

Model 1: weight ~ height
Model 2: weight ~ age/(height - 1)
  Res.Df     RSS Df Sum of Sq      F  Pr(>F)  
1      7 2991.57                              
2      3  149.01  4    2842.6 14.307 0.02696 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

注1: fm3上面的有6个系数，每组有一个截距和斜率。如果你想要4个系数，一个共同的截距和独立的斜率，那么使用

lm(weight ~ age:height, DF)

注2:我们也可以比较一个模型，在这个模型中水平的子集是相同的。例如，我们可以将1岁和2岁相同的模型与它们都相同(fm1)和完全不同(Fm3)的模型进行比较：

fm2 <- lm(weight ~ age/(height - 1), transform(DF, age = factor(c(1, 1, 3)[age])))
anova(fm1, fm2, fm3)

如果你做了大量的测试，你可能会偶然得到一些有意义的测试，所以你会想要降低p值的截止值。

注释3:这里有一些关于lm公式的注释：http://science.nature.nps.gov/im/datamgmt/statistics/r/formulas/

注意4:我们使用这个作为输入：

Lines <- "age height weight
1 56 140
1 60 155
1 64 143
2 56 117
2 60 125
2 64 133
3 74 245
3 75 241
3 82 269"
DF <- read.table(text = Lines, header = TRUE)
DF$age <- factor(DF$age)

Answer 2

在CrossValidated中有一个很好的答案。但简而言之，

require(emmeans)
data <- data.frame(age = factor(c(1,1,1,2,2,2,3,3,3)), 
               height = c(56,60,64,56,60,64,74,75,82), 
               weight = c(140,155,142,117,125,133,245,241,269))

model <- lm(weight ~ height*age, data)
anova(model) #check the results

Analysis of Variance Table

Response: weight
           Df  Sum Sq Mean Sq  F value    Pr(>F)    
height      1 25392.3 25392.3 481.5984 0.0002071 ***
age         2  2707.4  1353.7  25.6743 0.0129688 *  
height:age  2   169.0    84.5   1.6027 0.3361518    
Residuals   3   158.2    52.7                       
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


slopes <- emtrends(model, 'age', var = 'height') #gets each slope
slopes
age height.trend   SE df lower.CL upper.CL
1           0.25 1.28  3    -3.84     4.34
2           2.00 1.28  3    -2.09     6.09
3           3.37 1.18  3    -0.38     7.12

Confidence level used: 0.95 


pairs(slopes) #gets their comparisons two by two
contrast estimate   SE df t.ratio p.value
1 - 2       -1.75 1.82  3 -0.964  0.6441 
1 - 3       -3.12 1.74  3 -1.790  0.3125 
2 - 3       -1.37 1.74  3 -0.785  0.7363 

P value adjustment: tukey method for comparing a family of 3 estimates