如何指定 lm 模型矩阵
How to specify lm model matrix
我有从 2 个组获得的测量值，每个组具有相同的 3 个水平。
这是我的例子 data.frame:
df <- data.frame(measurement = c(rnorm(10,1,1),rnorm(10,0.75,1),rnorm(10,1.25,1),
rnorm(10,0.5,1),rnorm(10,1.75,1),rnorm(10,0.25,1)),
group = as.factor(c(rep("a",30),rep("b",30))),
level = as.factor(rep(c(rep("x",10),rep("y",10),rep("z",10)),2)))
我有兴趣量化每个 level 中的 measurement 如何受到 group 的影响。
我想线性模型 (lm) 是合适的方法，其中 group:level 交互项捕捉了我感兴趣的效果。
有没有办法指定一个 lm 只计算这些交互项：groupb:levelx、groupb:levely 和 groupb:levelz？我相信这告诉我每个 level 如何受到 group "b" 的影响（相对于 group "a"），我认为这是我感兴趣的.
我得到的最接近的是：
lm(measurement ~ 0 + group * level - group, data = df)
但这仍然会计算 levelx、levely 和 levelz 的效果，我对此不感兴趣。
正如@Lyzander 上面提到的，您应该更清楚地说明您想要什么。根据您所说的 "how measurement is affected by group "b“（相对于组 "a"）对于每个级别”，我想有 3 种简单的方法可以做到这一点。
df <- data.frame(measurement = c(rnorm(10,1,1),rnorm(10,0.75,1),rnorm(10,1.25,1),
rnorm(10,0.5,1),rnorm(10,1.75,1),rnorm(10,0.25,1)),
group = as.factor(c(rep("a",30),rep("b",30))),
level = as.factor(rep(c(rep("x",10),rep("y",10),rep("z",10)),2)))
library(dplyr)
#### calculate stats (mean values) ---------------------------------------------
df %>% group_by(level, group) %>% summarise(MeanMeasurement = mean(measurement))
# level group MeanMeasurement
# (fctr) (fctr) (dbl)
# 1 x a 1.6708659
# 2 x b 0.8487751
# 3 y a 0.7977769
# 4 y b 1.4209206
# 5 z a 1.5484668
# 6 z b -0.3244225
#### build a model for each level ---------------------------------------------
summary(lm(measurement ~ group , data = df[df$level=="x",]))
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 1.6709 0.3174 5.264 5.27e-05 ***
# groupb -0.8221 0.4489 -1.831 0.0837 .
summary(lm(measurement ~ group , data = df[df$level=="y",]))
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 0.7978 0.2565 3.111 0.00604 **
# groupb 0.6231 0.3627 1.718 0.10295
summary(lm(measurement ~ group , data = df[df$level=="z",]))
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 1.5485 0.3549 4.363 0.000375 ***
# groupb -1.8729 0.5019 -3.731 0.001528 **
## build a model only with interactions ------------------------------------------
summary(lm(measurement ~ group : level , data = df))
# Coefficients: (1 not defined because of singularities)
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) -0.3244 0.3123 -1.039 0.303452
# groupa:levelx 1.9953 0.4416 4.518 3.43e-05 ***
# groupb:levelx 1.1732 0.4416 2.657 0.010354 *
# groupa:levely 1.1222 0.4416 2.541 0.013951 *
# groupb:levely 1.7453 0.4416 3.952 0.000227 ***
# groupa:levelz 1.8729 0.4416 4.241 8.76e-05 ***
# groupb:levelz NA NA NA NA
如果您检查统计数据（第一种方法）和模型的系数，您会发现所有这 3 种方法都相互一致。
我会选择第二种方法，因为它是唯一一种可以为您提供有关 level 中 group（a 与 b）的差异是否具有统计显着性的信息的方法.第一种方法只是报告手段。第三种方法包括 p 值，但它们对应于与基线交互值的比较，而不是组 a 和 b 之间的比较。
您提到了 "ONLY compute these interaction terms: groupb:levelx, groupb:levely, and groupb:levelz"，这意味着您不会获得 a 和 x,y,z 的其他 3 个交互项。换句话说，您强制您的模型包含这 3 种交互。
您可以像这样手动完成
df <- data.frame(measurement = c(rnorm(10,1,1),rnorm(10,0.75,1),rnorm(10,1.25,1),
rnorm(10,0.5,1),rnorm(10,1.75,1),rnorm(10,0.25,1)),
group = as.factor(c(rep("a",30),rep("b",30))),
level = as.factor(rep(c(rep("x",10),rep("y",10),rep("z",10)),2)))
library(dplyr)
df %>%
mutate(interactions = paste0(group,":",level),
interactions = ifelse(group=="a","a",interactions)) -> df2
summary(lm(measurement ~ interactions, data = df2))
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 0.9318 0.1831 5.089 4.36e-06 ***
# interactionsb:x -0.7803 0.3662 -2.131 0.03752 *
# interactionsb:y 0.2747 0.3662 0.750 0.45638
# interactionsb:z -1.1367 0.3662 -3.104 0.00299 **
但现在其他 3 个交互被组合在一起，每次您将 3 个交互（b:x、b:y、b:z）中的每一个与一般组 a 进行比较。您不比较 x、y 和 z 中的 a 与 b，而是比较 b 组中的 x、y 和 z。
根据这句话："Is there a way to specify an lm that will only compute these interaction terms: groupb:levelx, groupb:levely, and groupb:levelz?"，我想你只是想要：
lm( measurement ~ level +0, subset = group=="b", data = df)
我有从 2 个组获得的测量值，每个组具有相同的 3 个水平。
这是我的例子 data.frame:
df <- data.frame(measurement = c(rnorm(10,1,1),rnorm(10,0.75,1),rnorm(10,1.25,1),
rnorm(10,0.5,1),rnorm(10,1.75,1),rnorm(10,0.25,1)),
group = as.factor(c(rep("a",30),rep("b",30))),
level = as.factor(rep(c(rep("x",10),rep("y",10),rep("z",10)),2)))
我有兴趣量化每个 level 中的 measurement 如何受到 group 的影响。
我想线性模型 (lm) 是合适的方法，其中 group:level 交互项捕捉了我感兴趣的效果。
有没有办法指定一个 lm 只计算这些交互项：groupb:levelx、groupb:levely 和 groupb:levelz？我相信这告诉我每个 level 如何受到 group "b" 的影响（相对于 group "a"），我认为这是我感兴趣的.
我得到的最接近的是：
lm(measurement ~ 0 + group * level - group, data = df)
但这仍然会计算 levelx、levely 和 levelz 的效果，我对此不感兴趣。
正如@Lyzander 上面提到的，您应该更清楚地说明您想要什么。根据您所说的 "how measurement is affected by group "b“（相对于组 "a"）对于每个级别”，我想有 3 种简单的方法可以做到这一点。
df <- data.frame(measurement = c(rnorm(10,1,1),rnorm(10,0.75,1),rnorm(10,1.25,1),
rnorm(10,0.5,1),rnorm(10,1.75,1),rnorm(10,0.25,1)),
group = as.factor(c(rep("a",30),rep("b",30))),
level = as.factor(rep(c(rep("x",10),rep("y",10),rep("z",10)),2)))
library(dplyr)
#### calculate stats (mean values) ---------------------------------------------
df %>% group_by(level, group) %>% summarise(MeanMeasurement = mean(measurement))
# level group MeanMeasurement
# (fctr) (fctr) (dbl)
# 1 x a 1.6708659
# 2 x b 0.8487751
# 3 y a 0.7977769
# 4 y b 1.4209206
# 5 z a 1.5484668
# 6 z b -0.3244225
#### build a model for each level ---------------------------------------------
summary(lm(measurement ~ group , data = df[df$level=="x",]))
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 1.6709 0.3174 5.264 5.27e-05 ***
# groupb -0.8221 0.4489 -1.831 0.0837 .
summary(lm(measurement ~ group , data = df[df$level=="y",]))
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 0.7978 0.2565 3.111 0.00604 **
# groupb 0.6231 0.3627 1.718 0.10295
summary(lm(measurement ~ group , data = df[df$level=="z",]))
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 1.5485 0.3549 4.363 0.000375 ***
# groupb -1.8729 0.5019 -3.731 0.001528 **
## build a model only with interactions ------------------------------------------
summary(lm(measurement ~ group : level , data = df))
# Coefficients: (1 not defined because of singularities)
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) -0.3244 0.3123 -1.039 0.303452
# groupa:levelx 1.9953 0.4416 4.518 3.43e-05 ***
# groupb:levelx 1.1732 0.4416 2.657 0.010354 *
# groupa:levely 1.1222 0.4416 2.541 0.013951 *
# groupb:levely 1.7453 0.4416 3.952 0.000227 ***
# groupa:levelz 1.8729 0.4416 4.241 8.76e-05 ***
# groupb:levelz NA NA NA NA
如果您检查统计数据（第一种方法）和模型的系数，您会发现所有这 3 种方法都相互一致。
我会选择第二种方法，因为它是唯一一种可以为您提供有关 level 中 group（a 与 b）的差异是否具有统计显着性的信息的方法.第一种方法只是报告手段。第三种方法包括 p 值，但它们对应于与基线交互值的比较，而不是组 a 和 b 之间的比较。
您提到了 "ONLY compute these interaction terms: groupb:levelx, groupb:levely, and groupb:levelz"，这意味着您不会获得 a 和 x,y,z 的其他 3 个交互项。换句话说，您强制您的模型包含这 3 种交互。
您可以像这样手动完成
df <- data.frame(measurement = c(rnorm(10,1,1),rnorm(10,0.75,1),rnorm(10,1.25,1),
rnorm(10,0.5,1),rnorm(10,1.75,1),rnorm(10,0.25,1)),
group = as.factor(c(rep("a",30),rep("b",30))),
level = as.factor(rep(c(rep("x",10),rep("y",10),rep("z",10)),2)))
library(dplyr)
df %>%
mutate(interactions = paste0(group,":",level),
interactions = ifelse(group=="a","a",interactions)) -> df2
summary(lm(measurement ~ interactions, data = df2))
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 0.9318 0.1831 5.089 4.36e-06 ***
# interactionsb:x -0.7803 0.3662 -2.131 0.03752 *
# interactionsb:y 0.2747 0.3662 0.750 0.45638
# interactionsb:z -1.1367 0.3662 -3.104 0.00299 **
但现在其他 3 个交互被组合在一起，每次您将 3 个交互（b:x、b:y、b:z）中的每一个与一般组 a 进行比较。您不比较 x、y 和 z 中的 a 与 b，而是比较 b 组中的 x、y 和 z。
根据这句话："Is there a way to specify an lm that will only compute these interaction terms: groupb:levelx, groupb:levely, and groupb:levelz?"，我想你只是想要：
lm( measurement ~ level +0, subset = group=="b", data = df)