Как определить матрицу модели lm

Question

У меня есть измерения, полученные из 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))) 

Я заинтересован в количественной оценке, как measurement в каждом level зависит от group.

Я предполагаю линейную модель (lm) является подходящим подходом для этого, где group:level условия взаимодействия захватить эффекты Я заинтересован в.

Есть ли способ, чтобы указать, что lm будет только вычислить их условия взаимодействия: groupb:levelx, groupb:levely, и groupb:levelz? Я считаю, что это говорит мне, как каждый level влияет group «Ъ» (относительно group «а»), который я думаю, что я заинтересован в

Ближайший я получил это:.

lm(measurement ~ 0 + group * level - group, data = df) 

Но до сих пор вычисляет эффекты levelx, levely и levelz, что я не заинтересован.

Answer 1

Как @Lyzander упоминалось выше, вы должны прояснить немного больше, что вы хотите. Основываясь на том, что вы сказали «как измерение влияет на группу« 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 подхода согласуются друг с другом.

Я бы пошел со вторым подходом, так как это единственный, который дает вам информацию о том, является ли разница group (a vs b) в пределах level статистически значимой или нет. 1-й подход означает только сообщения. 3-й подход включает значения p, но они соответствуют сравнениям с базовым значением взаимодействия, а не сравнениям между группами a и b.

Вы упомянули «ТОЛЬКО вычислить эти условия взаимодействия: groupb: levelx, groupb: levely и groupb: levelz», что означает, что вы не получите других 3 взаимодействий a и x, y, z. Другими словами, вы вынуждаете свою модель включать эти 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 взаимодействий (б: х, б: у, б: г) с общей группой а. Вы не сравниваете v и b внутри x, y и z, но сравниваете x vs. y vs. z внутри группы b.

Answer 2

Основываясь на этом предложении: «Есть ли способ, чтобы указать лм, что будет только вычислять эти условия взаимодействия: GroupB: levelx, GroupB: levely и GroupB: levelz», я думаю, вы просто хотите:

lm(measurement ~ level +0, subset = group=="b", data = df)