作者：黄相源，复旦大学儿科医院硕士
Syntax
√输入1：
# 交互作用模型中分组效应及标准误的计算
rm(list=ls())
library(gmodels)
library(ggplot2)
### 生成模拟数据
set.seed(191114)
x 1000, mean=10, sd=5)
group 1:2
inter group
y 3*x + 4*group + 5*inter + rnorm(1000, mean=0, sd=5)
# group=1时, x变化一个单位，y值变化3+5个单位
# group=2时，x变化一个单位，y值变化3+5*2个单位
df group=group, inter=inter, y=y)
head(df)√结果1：
x group inter y
1 1.439778 1 1.439778 14.14025
2 1.603191 2 3.206383 24.15817
3 13.397201 1 13.397201 113.65427
4 5.260796 2 10.521591 81.34142
5 1.706308 1 1.706308 18.59737
6 7.986918 2 15.973836 107.56454√输入2：
plot(df$x, df$y, col=df$group)√结果2：
√输入3：
### 1.使用gmodels包下estimable计算各组效应
model1 group+inter)
summary(model1)√结果3：
Call:
lm(formula = y ~ x + group + inter, data = df)
Residuals:
Min 1Q Median 3Q Max
-14.5072 -3.2593 0.1235 3.2859 14.3020
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.11982 1.04732 -1.069 0.285
x 3.16932 0.09293 34.104 2e-16 ***
group 4.78550 0.65818 7.271 7.22e-13 ***
inter 4.88762 0.05847 83.586 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.73 on 996 degrees of freedom
Multiple R-squared: 0.9941, Adjusted R-squared: 0.9941
F-statistic: 5.605e+04 on 3 and 996 DF, p-value: 2.2e-16√输入4：
# group=1时，x对y的效应量及标准误
estimable(model1, c(0,1,0,1))√结果4：
Estimate Std. Error t value DF Pr(>|t|)
(0 1 0 1) 8.056947 0.04170063 193.2092 996 0√输入5：
# group=2时，x对y的效应量及标准误
estimable(model1, c(0,1,0,2))√结果5：
Estimate Std. Error t value DF Pr(>|t|)
(0 1 0 2) 12.94457 0.04099109 315.7899 996√输入6：
# 同时显示两组中的效应量
cm 'Group 1' = c(0, 1, 0, 1),
'Group 2' = c(0, 1, 0, 2))
estimable(model1, cm)√结果6：
Estimate Std. Error t value DF Pr(>|t|)
Group 1 8.056947 0.04170063 193.2092 996 0
Group 2 12.944571 0.04099109 315.7899 996 0√输入7：
### 2.操纵数据以计算各组效应
# group=1时的效应量及标准误
df_temp1 df_temp1$group $group-1
df_temp1$inter $group*df_temp1$x
summary(lm(data=df_temp1, y~x+group+inter))√结果7：
Call:
lm(formula = y ~ x + group + inter, data = df_temp1)
Residuals:
Min 1Q Median 3Q Max
-14.5072 -3.2593 0.1235 3.2859 14.3020
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.66568 0.47035 7.794 1.63e-14 ***
x 8.05695 0.04170 193.209 2e-16 ***
group 4.78550 0.65818 7.271 7.22e-13 ***
inter 4.88762 0.05847 83.586 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.73 on 996 degrees of freedom
Multiple R-squared: 0.9941, Adjusted R-squared: 0.9941
F-statistic: 5.605e+04 on 3 and 996 DF, p-value: 2.2e-16√输入8：
# group=2时的效应量及标准误
df_temp2 df_temp2$group $group-2
df_temp2$inter $group*df_temp2$x
summary(lm(data=df_temp2, y~x+group+inter))√结果8：
Call:
lm(formula = y ~ x + group + inter, data = df_temp2)
Residuals:
Min 1Q Median 3Q Max
-14.5072 -3.2593 0.1235 3.2859 14.3020
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.45118 0.46040 18.356 2e-16 ***
x 12.94457 0.04099 315.790 2e-16 ***
group 4.78550 0.65818 7.271 7.22e-13 ***
inter 4.88762 0.05847 83.586 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.73 on 996 degrees of freedom
Multiple R-squared: 0.9941, Adjusted R-squared: 0.9941
F-statistic: 5.605e+04 on 3 and 996 DF, p-value: 2.2e-16√输入9：
### 3.手动计算各组效应
model1 <- lm(data=df, y~x+group+inter)summary(model1)
# 模型中 b0=-1.12, b1=3.17, b2=4.79, b3=4.89
# 各参数标准误为seb0= 1.05, seb1=0.09, seb2=0.66, seb3=0.06coef <- summary(model1)$coefficientsb1 <- coef[2,1]b3 <- coef[4,1]seb1 <- coef[2,2]seb3 <- coef[4,2]
# 当group=n时，y随x变化的效应量effect=b1+n*b3
# effect的标准误等于 seb1^2 + n^2*seb3^2 + 2*n*cov(b1, b3) 的平方根
# cov(b1, b3)为系数b1和b3的协方差cov_matrix <- vcov(model1)covb1b3 <- cov_matrix[2,4]
# group=1时effect1 <- b1 + b3seeffect1 <- sqrt(seb1^2 + 1^2*seb3^2 + 2*1*covb1b3)cat("Effect: ", effect1, " , Standard Error: ", seeffect1)√结果9：
Effect: 8.056947 , Standard Error: 0.04170063√输入10：
# group=2时
effect2 2*b3
seeffect2 sqrt(seb1^2 + 2^2*seb3^2 + 2*2*covb1b3)
cat("Effect: ", effect2, " , Standard Error: ", seeffect2)√结果10：
Effect: 12.94457 , Standard Error: 0.04099109本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。
