如何在我的一组数据集中执行ANCOVA？

Question

我无法找到使用 ANCOVA 进行亚组分析的适当方法。在我的分析中，我想比较两组之间的 6 个月抑郁评分 (six_PHQ)：根据婚姻状况 (Marital_Status) 进行干预与对照 (Study.Arm)。我还控制基线分数 (base_PHQ) 作为协变量。这些是数据框 PHQ1 中的列。

  dput(PHQ1[1:4, ])
structure(list(Patient = c("ASSC-011-JB", "ASSC-013-ME", "ASSC-019-JM", 
"ASSC-053-RH"), **Study.Arm = c("Control", "Control", "Control"**, 
**"Control")**, DOB = c("1968-01-11", "1939-06-12", "1955-03-24", 
"1944-05-04"), Enrolment.Date = c("2021-04-23", "2021-03-19", 
"2021-03-25", "2021-06-29"), Age_Enrollment = c(53L, 81L, 66L, 
77L), Age_Today = c(56L, 84L, 68L, 79L), Gender = c("m", "f", 
"f", "f"), **Marital_Status = structure(c(2L, 2L, 1L, 2L)**, levels = c("Married", 
"Not_Married"), class = "factor"), Highest_Education = c(4L, 
4L, 3L, 2L), Employment_Type = c(5L, 2L, 2L, 1L), Income_Annual = c(7L, 
2L, 1L, 1L), Canadian_Born = c(1L, 2L, 1L, 1L), Ethnicity = c("White", 
"White", "White", "White"), Number.of.Chronic.conditions = c(2L, 
5L, 5L, 6L), **base_PHQ = c(0L, 3L, 1L, 4L), six_PHQ = c(3L, 1L, 
0L, 2L)**, Change.Score = c(-3L, 2L, 1L, 2L), Clinical.Change = c("Minimal Increase", 
"Minimal Decrease", "Minimal Decrease", "Minimal Decrease"), 
    Change.of.5 = c("No Change", "No Change", "No Change", "No Change"
    )), row.names = c(NA, 4L), class = "data.frame").```

   
I tried using the `   group_by   ` function and making a lm model within that grouping, but this is not quite what I was hoping for:

`   Ancova_Marital <- PHQ1 %>%
   group_by(Marital_Status) %>%
do(model = lm(six_PHQ ~ Study.Arm + base_PHQ, data = PHQ1)) %>%
   summarise(
     Marital_Status = first(Marital_Status),
     beta = coef(model)[2],  # Slope coefficient (effect of covariate)
     p_value = summary(model)$coef[2, 4] # p-value for the slope coefficient 
   )    `

`   print(Ancova_Marital)
  Marital_Status   beta p_value
  <fct>           <dbl>   <dbl>
1 Married        -0.231   0.734
2 Not_Married    -0.231   0.734   `


With this method, I don't know how to take it a step further and use the ANCOVA code: `    Anova([aov model], type=3)    `. Apologies if this is unclear. I am new to using R and stack overflow. Thank you in advance. 

Answer 1

我不是 100% 确定我明白您到底在寻找什么，但这里有一个基于

Marital_Status

拆分数据并为每个子组执行 ANCOVA 的建议：

library(tidyverse)

# Creating a sample dataset
n <- 100 # number of observations
set.seed(0) # seed for reproducibility
data <- tibble(
  six_PHQ = sample(0:30, n, replace = T),
  Study.Arm = sample(c("Control", "Intervention"), n, replace = T),
  base_PHQ = sample(0:30, n, replace = T),
  Marital_Status = sample(c("Married", "Not_Married"), n, replace = T)
)

# Grouping the data by "Marital_Status" and creating nested data frames
data_subgroups <- data %>% 
  group_by(Marital_Status) %>% 
  nest()

# Custom function to perform ANCOVA analysis
perform_ANCOVA <- function(df) {
  model <- lm(six_PHQ ~ Study.Arm + base_PHQ, data = df)
  anova_result <- anova(model)
  return(anova_result)
}

# Applying ANCOVA analysis to each subgroup
result <- data_subgroups %>% 
  mutate(ANCOVA_result = map(data, perform_ANCOVA))

这是您要找的吗？

即使我不知道您正在分析的研究，您也可能会考虑将所有内容放入一个模型中，并且可能添加与

Marital_Status

的交互而不是分组。