Tigermed-data-about-depression-clinical-trial

intro

basic idea

Repeated measurements are taken on the same individuals over time.

The primary objective of most longitudinal studies is to characterize changes in a particular outcome variable of interest over time. However, we are not only interested in change itself but also in the factors that may influence this change. These factors could be treatment interventions, such as dummy variables abc, and we want to know whether the treatment has an effect on change over time. Does the intervention lead to greater improvements in the outcome over time? It could also be exposures like sex.

There are two aspects of the data that complicate the analysis, meaning that we cannot use many of the routine regression techniques we are already familiar with. The first is that if you take repeated measurements over time, those measurements will be correlated and are usually positively correlated. For example, if my systolic blood pressure is high today, it provides more information about future measurements. The second is that when you take repeated measurements over time, the variability of those measurements is often heterogeneous over time. This means that the variability of the outcome at the start of the longitudinal study can often be different from the variability of the outcome later on.

Both correlation and variability play important roles in this model, and therefore covariance matrix is particularly significant.

Introduce occasions or time, which are not the same.

This study period is short, just 28 days, which is quicker when computing the unstructure matrix

The data is balanced, meaning that each individual has been measured at the same occasions.

further improvement

The mixed model can be more flexible for smoothing and semi-parametric regression.

compareGroups

# 将SITEID分为5组
library(dplyr)

Attaching package: 'dplyr'
The following objects are masked from 'package:stats':

    filter, lag
The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union
library(dplyr)
library(compareGroups)
library(readxl)
data0 <- read_excel("3.xlsx")

# 将SITEID转换为数值型变量
data0 <- data0 %>%
  mutate(SITEID_num = as.numeric(factor(SITEID)))

# 使用ntile函数分组
data0 <- data0 %>%
  mutate(SITEID_group = paste0("Group", ntile(SITEID_num, 5)))


data0$SITEID_group <- as.factor(data0$SITEID_group)



data0$madrsc <- ifelse(data0$BASE > 40, "MADRS SCORE ≥ 40"," MADRS SCORE < 40")

data0$TRT01A <- as.factor(data0$TRT01A)
data0$TRT01A <- relevel(data0$TRT01A, ref = "Dummy A")
data0$SEX <- as.factor(data0$SEX)
data0$RACE <- as.factor(data0$RACE)
data0$RACE <- relevel(data0$RACE, ref = "WHITE")



# 选择变量并创建表格
data_table1 <- data0 %>%
  select(madrsc, AGE, SEX, RACE, TRT01A, PCHG, CHG, BASE, AVAL,SITEID_group )

# 检查变量类型
str(data_table1)
tibble [3,243 × 10] (S3: tbl_df/tbl/data.frame)
 $ madrsc      : chr [1:3243] "MADRS SCORE ≥ 40" "MADRS SCORE ≥ 40" "MADRS SCORE ≥ 40" "MADRS SCORE ≥ 40" ...
 $ AGE         : num [1:3243] 37 37 37 37 37 37 37 53 53 53 ...
 $ SEX         : Factor w/ 2 levels "F","M": 1 1 1 1 1 1 1 1 1 1 ...
 $ RACE        : Factor w/ 8 levels "WHITE","AMERICAN INDIAN OR ALASKA NATIVE",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ TRT01A      : Factor w/ 3 levels "Dummy A","Dummy B",..: 1 1 1 1 1 1 1 2 2 2 ...
 $ PCHG        : num [1:3243] NA -19.57 -8.7 -6.52 -8.7 ...
 $ CHG         : num [1:3243] NA -9 -4 -3 -4 -3 -3 NA -17 -10 ...
 $ BASE        : num [1:3243] 46 46 46 46 46 46 46 35 35 35 ...
 $ AVAL        : num [1:3243] 46 37 42 43 42 43 43 35 18 25 ...
 $ SITEID_group: Factor w/ 5 levels "Group1","Group2",..: 1 1 1 1 1 1 1 1 1 1 ...
# 检查样本量
table(data_table1$madrsc)

 MADRS SCORE < 40  MADRS SCORE ≥ 40 
             2525               716 
# 创建比较表格
table1 <- compareGroups(madrsc~., data = data_table1)
Warning: glm.fit: algorithm did not converge
Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
# 如果仍然出现警告，可以尝试使用其他方法来处理变量的尺度差异
# 或者检查数据是否存在数据分离问题等
table1 <- createTable(table1, show.all = TRUE, hide.no = "no",
                      show.p.overall=TRUE)
table1

--------Summary descriptives table by 'madrsc'---------

________________________________________________________________________________________________________ 
                                                  [ALL]      MADRS SCORE < 40 MADRS SCORE ≥ 40 p.overall 
                                                 N=3241          N=2525            N=716                 
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
AGE                                            45.2 (14.0)     45.5 (14.5)      44.4 (12.3)      0.045   
SEX:                                                                                             0.006   
    F                                         1928 (59.5%)    1470 (58.2%)      458 (64.0%)              
    M                                         1313 (40.5%)    1055 (41.8%)      258 (36.0%)              
RACE:                                                                                              .     
    WHITE                                     2788 (86.0%)    2177 (86.2%)      611 (85.3%)              
    AMERICAN INDIAN OR ALASKA NATIVE            6 (0.19%)       6 (0.24%)        0 (0.00%)               
    ASIAN                                      97 (2.99%)      83 (3.29%)        14 (1.96%)              
    BLACK OR AFRICAN AMERICAN                  241 (7.44%)     185 (7.33%)       56 (7.82%)              
    MULTIPLE                                   60 (1.85%)      46 (1.82%)        14 (1.96%)              
    NATIVE HAWAIIAN OR OTHER PACIFIC ISLANDER  14 (0.43%)       7 (0.28%)        7 (0.98%)               
    NOT REPORTED                               28 (0.86%)      21 (0.83%)        7 (0.98%)               
    UNKNOWN                                     7 (0.22%)       0 (0.00%)        7 (0.98%)               
TRT01A:                                                                                          0.429   
    Dummy A                                   1610 (49.7%)    1252 (49.6%)      358 (50.0%)              
    Dummy B                                    819 (25.3%)     650 (25.7%)      169 (23.6%)              
    Dummy C                                    812 (25.1%)     623 (24.7%)      189 (26.4%)              
PCHG                                          -24.64 (29.4)   -24.71 (29.8)    -24.39 (27.9)     0.802   
CHG                                           -8.82 (10.4)    -8.31 (9.82)     -10.61 (12.2)    <0.001   
BASE                                           35.4 (6.36)     33.1 (5.09)      43.6 (2.26)      0.000   
AVAL                                           27.9 (11.1)     26.0 (10.0)      34.5 (12.1)     <0.001   
SITEID_group:                                                                                   <0.001   
    Group1                                     649 (20.0%)     486 (19.2%)      163 (22.8%)              
    Group2                                     649 (20.0%)     516 (20.4%)      133 (18.6%)              
    Group3                                     647 (20.0%)     430 (17.0%)      217 (30.3%)              
    Group4                                     648 (20.0%)     552 (21.9%)       96 (13.4%)              
    Group5                                     648 (20.0%)     541 (21.4%)      107 (14.9%)              
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
summary(table1)



---Available data----

_________________________________________________________________________________________ 
             [ALL]  MADRS SCORE < 40 MADRS SCORE ≥ 40      method       select Fact OR/HR 
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
AGE          3241        2525              716        continuous-normal  ALL       1      
SEX          3241        2525              716           categorical     ALL       --     
RACE         3241        2525              716           categorical     ALL       --     
TRT01A       3241        2525              716           categorical     ALL       --     
PCHG         2765        2153              612        continuous-normal  ALL       1      
CHG          2765        2153              612        continuous-normal  ALL       1      
BASE         3241        2525              716        continuous-normal  ALL       1      
AVAL         3241        2525              716        continuous-normal  ALL       1      
SITEID_group 3241        2525              716           categorical     ALL       --     
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
?compareGroups
library(dplyr)
library(compareGroups)
data0 <- read_excel("3.xlsx")

data0$madrsc <- ifelse(data0$BASE > 40, "MADRS SCORE ≥ 40"," MADRS SCORE < 40")

data0$AGE <- scale(data0$AGE)
data0$AVAL <- scale(data0$AVAL)
data0$BASE <- scale(data0$BASE)
data0$TRT01A <- as.factor(data0$TRT01A)
data0$TRT01A <- relevel(data0$TRT01A, ref = "Dummy A")
data0$SEX <- as.factor(data0$SEX)
data0$RACE <- as.factor(data0$RACE)
data0$RACE <- relevel(data0$RACE, ref = "WHITE")


data0$PCHG <- scale(data0$PCHG)

data_table1 <- data0 %>%
  select(madrsc, AGE, SEX, RACE, TRT01A, PCHG, CHG, BASE, AVAL)

table1 <- compareGroups(madrsc~., data = data_table1)
Warning: glm.fit: algorithm did not converge
Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
table1 <- createTable(table1, show.all = TRUE, hide.no = "no",
                      show.p.overall=TRUE)
table1

--------Summary descriptives table by 'madrsc'---------

_______________________________________________________________________________________________________ 
                                                 [ALL]      MADRS SCORE < 40 MADRS SCORE ≥ 40 p.overall 
                                                 N=3241         N=2525            N=716                 
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
AGE                                           0.00 (1.00)     0.02 (1.03)      -0.06 (0.88)     0.045   
SEX:                                                                                            0.006   
    F                                         1928 (59.5%)   1470 (58.2%)      458 (64.0%)              
    M                                         1313 (40.5%)   1055 (41.8%)      258 (36.0%)              
RACE:                                                                                             .     
    WHITE                                     2788 (86.0%)   2177 (86.2%)      611 (85.3%)              
    AMERICAN INDIAN OR ALASKA NATIVE           6 (0.19%)       6 (0.24%)        0 (0.00%)               
    ASIAN                                      97 (2.99%)     83 (3.29%)        14 (1.96%)              
    BLACK OR AFRICAN AMERICAN                 241 (7.44%)     185 (7.33%)       56 (7.82%)              
    MULTIPLE                                   60 (1.85%)     46 (1.82%)        14 (1.96%)              
    NATIVE HAWAIIAN OR OTHER PACIFIC ISLANDER  14 (0.43%)      7 (0.28%)        7 (0.98%)               
    NOT REPORTED                               28 (0.86%)     21 (0.83%)        7 (0.98%)               
    UNKNOWN                                    7 (0.22%)       0 (0.00%)        7 (0.98%)               
TRT01A:                                                                                         0.429   
    Dummy A                                   1610 (49.7%)   1252 (49.6%)      358 (50.0%)              
    Dummy B                                   819 (25.3%)     650 (25.7%)      169 (23.6%)              
    Dummy C                                   812 (25.1%)     623 (24.7%)      189 (26.4%)              
PCHG                                          0.00 (1.00)     0.00 (1.01)      0.01 (0.95)      0.802   
CHG                                           -8.82 (10.4)   -8.31 (9.82)     -10.61 (12.2)    <0.001   
BASE                                          0.00 (1.00)    -0.37 (0.80)      1.29 (0.36)      0.000   
AVAL                                          0.00 (1.00)    -0.17 (0.90)      0.60 (1.09)     <0.001   
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
summary(table1)



---Available data----

___________________________________________________________________________________ 
       [ALL]  MADRS SCORE < 40 MADRS SCORE ≥ 40      method       select Fact OR/HR 
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
AGE    3241        2525              716        continuous-normal  ALL       1      
SEX    3241        2525              716           categorical     ALL       --     
RACE   3241        2525              716           categorical     ALL       --     
TRT01A 3241        2525              716           categorical     ALL       --     
PCHG   2765        2153              612        continuous-normal  ALL       1      
CHG    2765        2153              612        continuous-normal  ALL       1      
BASE   3241        2525              716        continuous-normal  ALL       1      
AVAL   3241        2525              716        continuous-normal  ALL       1      
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

+++++see mean

# Load required libraries
library(readxl)  # For reading Excel files
library(dplyr)   # For data manipulation
library(tidyr)   # For data manipulation

Attaching package: 'tidyr'
The following object is masked _by_ '.GlobalEnv':

    table1
# Read the data from the Excel file
data <- read_excel("3.xlsx")

# Convert AVISIT to a numeric factor to sort by the day of the visit
data$AVISIT <- factor(data$AVISIT, levels = unique(data$AVISIT)[order(as.numeric(gsub("Day \\d+ \\(", "", gsub(" \\(DB\\)", "", data$AVISIT)) ))])
Warning in order(as.numeric(gsub("Day \\d+ \\(", "", gsub(" \\(DB\\)", "", :
NAs introduced by coercion
# Group by AVISIT and TRT01A, then calculate mean of AVAL for each group
result <- data %>%
  group_by(AVISIT, TRT01A) %>%
  summarise(mean_AVAL = mean(AVAL, na.rm = TRUE)) %>%
  ungroup() %>%  # Remove grouping to allow sorting
  arrange(AVISIT) %>%  # Sort by AVISIT (now a numeric factor)
  pivot_wider(names_from = TRT01A, values_from = mean_AVAL)  # Pivot the table
`summarise()` has grouped output by 'AVISIT'. You can override using the
`.groups` argument.
# Display the result
print(result)
# A tibble: 7 × 4
  AVISIT        `Dummy A` `Dummy B` `Dummy C`
  <fct>             <dbl>     <dbl>     <dbl>
1 Baseline (DB)      35.0      35.6      35.5
2 Day 2 (DB)         24.0      25.5      24.9
3 Day 8 (DB)         27.9      29.2      27.6
4 Day 15 (DB)        27.5      29.1      27.9
5 Day 22 (DB)        26.2      27.2      25.9
6 Day 28 (DB)        25.4      27.5      26.0
7 Endpoint (DB)      25.7      27.8      26.1
#data$mean_AVAL <- mean(data$AVAL, na.rm = TRUE)


# ggplot2 for visualization --3 lines in a plot and y-axis is AVAL, x-axis is AVISIT
library(ggplot2)
#ggplot(data$mean_AVAL ,aes(x = AVISIT)) +
 # geom_line(aes(y = `Dummy A`, color = "Dummy A")) +
  #geom_line(aes(y = `Dummy B`, color = "Dummy B")) +
  #geom_line(aes(y = `Dummy C`, color = "Dummy C")) +
  #labs(title = "Mean AVAL by AVISIT and Treatment Group",
    #   x = "Visit",
     #  y = "Mean AVAL") +
#  theme_minimal() +
  #s#cale_color_manual(values = c("blue", "red", "green")) +
  #theme(legend.title = element_blank())
# Load required libraries


# For plotting, keep the data in long format
plot_data <- data %>%
  group_by(AVISIT, TRT01A) %>%
  summarise(mean_AVAL = mean(AVAL, na.rm = TRUE), .groups = "drop")

# Create the plot
ggplot(plot_data, aes(x = AVISIT, y = mean_AVAL, color = TRT01A, group = TRT01A)) +
  geom_line() +
  geom_point() +
  labs(title = "Mean AVAL by AVISIT and Treatment Group",
       x = "Visit",
       y = "Mean AVAL") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_color_manual(values = c("blue", "red", "green"))

# Optionally, you can write the result to a new Excel file
# library(writexl)
# write_xlsx(result, "result.xlsx")

main result

library(readxl)
library(mmrm)
Warning: package 'mmrm' was built under R version 4.4.3
library(ggplot2)
library(dplyr)


admadrs <- read_excel("3.xlsx", 
                      sheet = "Sheet1", col_names = TRUE, 
                      col_types = NULL, na = "", skip = 0)
# 创建分类变量 base_situation
admadrs$base_situation <- cut(
  admadrs$BASE,
  breaks = c( 9, 19, 39, Inf),
  labels = c("10-19", "20-39", "40+"),

)



# 将 base_situation 转换为因子类型
admadrs$base_situation <- as.factor(admadrs$base_situation)

admadrs$USUBJID <- factor(admadrs$USUBJID)
admadrs$AVISITN <- admadrs$AVISITN - 20000
admadrs$AVISITN <- factor(admadrs$AVISITN)
admadrs$CHG <- as.numeric(admadrs$CHG)

admadrs$base_situation <- relevel(admadrs$base_situation, ref = "40+")

# Delete a duplicate point in time
admadrs <- admadrs[!duplicated(admadrs[c("USUBJID", "AVISIT")]), ]

# scale data 
admadrs$AGE <- scale(admadrs$AGE)
admadrs$AVAL <- scale(admadrs$AVAL)




# 更新模型以包括 base_situation
fit <- mmrm(
  formula = AVAL ~ AVISITN:TRT01A +TRT01A + SITEID+ RACE + AGE + base_situation + us(AVISITN | USUBJID),
  data = admadrs)

# 获取模型摘要
model_summary <- summary(fit)

model_summary
mmrm fit

Formula:     
AVAL ~ AVISITN:TRT01A + TRT01A + SITEID + RACE + AGE + base_situation +  
    us(AVISITN | USUBJID)
Data:        admadrs (used 3214 observations from 476 subjects with maximum 7 
timepoints)
Covariance:  unstructured (28 variance parameters)
Method:      Satterthwaite
Vcov Method: Asymptotic
Inference:   REML

Model selection criteria:
     AIC      BIC   logLik deviance 
  3852.3   3969.0  -1898.2   3796.3 

Coefficients: 
                                                Estimate Std. Error         df
(Intercept)                                     1.806053   0.375454 415.500000
TRT01ADummy B                                   0.001719   0.041133 416.000000
TRT01ADummy C                                   0.048026   0.041752 416.900000
SITEIDUS10002                                   0.076913   0.124708 415.700000
SITEIDUS10004                                  -0.221056   0.172061 415.500000
SITEIDUS10005                                   0.238437   0.362054 415.300000
SITEIDUS10006                                  -0.142527   0.123867 415.700000
SITEIDUS10008                                  -0.555081   0.189226 416.700000
SITEIDUS10009                                  -0.108724   0.134404 415.500000
SITEIDUS10013                                   0.098598   0.120475 415.300000
SITEIDUS10014                                   0.120042   0.135360 417.100000
SITEIDUS10015                                   0.154076   0.121021 415.400000
SITEIDUS10016                                  -0.281098   0.261486 421.100000
SITEIDUS10017                                  -0.126686   0.158631 416.200000
SITEIDUS10019                                  -0.076086   0.134493 415.300000
SITEIDUS10020                                  -0.359461   0.139698 415.400000
SITEIDUS10021                                  -0.156841   0.116689 415.300000
SITEIDUS10023                                  -0.622644   0.258943 415.500000
SITEIDUS10025                                   0.030272   0.158685 415.300000
SITEIDUS10027                                  -0.178228   0.260325 415.300000
SITEIDUS10029                                   0.035942   0.108453 415.500000
SITEIDUS10030                                  -0.061651   0.101481 415.600000
SITEIDUS10031                                  -0.054005   0.151611 415.300000
SITEIDUS10033                                   0.003187   0.155735 433.100000
SITEIDUS10034                                   0.038837   0.171523 415.300000
SITEIDUS10036                                   0.301038   0.125752 415.300000
SITEIDUS10039                                   0.569455   0.361001 415.300000
SITEIDUS10040                                   0.063481   0.190898 415.500000
SITEIDUS10041                                   0.129966   0.157449 415.500000
SITEIDUS10043                                  -0.173603   0.116430 415.800000
SITEIDUS10044                                  -0.083537   0.215919 415.800000
SITEIDUS10045                                  -0.014193   0.188651 415.300000
SITEIDUS10046                                   0.183459   0.159190 415.700000
SITEIDUS10049                                   0.111136   0.158066 415.300000
SITEIDUS10051                                  -0.070450   0.114373 415.300000
SITEIDUS10053                                   0.392623   0.385587 415.400000
SITEIDUS10055                                   0.012793   0.122778 415.400000
SITEIDUS10056                                  -0.386136   0.129926 416.100000
SITEIDUS10057                                   0.068405   0.136772 415.600000
SITEIDUS10058                                  -0.128815   0.120894 415.400000
SITEIDUS10062                                  -0.006573   0.088860 415.600000
SITEIDUS10068                                   0.069483   0.106422 415.300000
SITEIDUS10069                                   0.260331   0.215821 415.300000
SITEIDUS10070                                  -0.093592   0.087621 415.600000
SITEIDUS10071                                  -0.280235   0.160520 415.500000
SITEIDUS10073                                  -0.087786   0.361084 415.200000
SITEIDUS10074                                  -0.125308   0.111749 415.900000
SITEIDUS10075                                  -0.003134   0.088734 415.500000
SITEIDUS10076                                  -0.361377   0.148825 415.800000
RACEASIAN                                      -0.283532   0.381605 415.400000
RACEBLACK OR AFRICAN AMERICAN                  -0.566955   0.373347 415.400000
RACEMULTIPLE                                   -0.648127   0.392952 415.500000
RACENATIVE HAWAIIAN OR OTHER PACIFIC ISLANDER  -0.437063   0.450218 415.500000
RACENOT REPORTED                               -0.271027   0.414808 415.500000
RACEUNKNOWN                                    -0.484391   0.521374 415.500000
RACEWHITE                                      -0.463115   0.367678 415.400000
AGE                                             0.021668   0.018133 417.800000
base_situation10-19                            -2.339846   0.187831 415.600000
base_situation20-39                            -0.857016   0.040751 415.600000
TRT01ADummy A:AVISITN2                         -1.006297   0.059276 468.500000
TRT01ADummy B:AVISITN2                         -0.934629   0.082854 467.300000
TRT01ADummy C:AVISITN2                         -0.981017   0.083632 472.400000
TRT01ADummy A:AVISITN8                         -0.651015   0.056054 463.500000
TRT01ADummy B:AVISITN8                         -0.569553   0.078328 462.100000
TRT01ADummy C:AVISITN8                         -0.716786   0.079153 469.200000
TRT01ADummy A:AVISITN15                        -0.694693   0.059041 464.500000
TRT01ADummy B:AVISITN15                        -0.583654   0.082304 459.400000
TRT01ADummy C:AVISITN15                        -0.697462   0.082969 464.100000
TRT01ADummy A:AVISITN22                        -0.822214   0.065033 470.600000
TRT01ADummy B:AVISITN22                        -0.742937   0.090897 468.700000
TRT01ADummy C:AVISITN22                        -0.838812   0.091607 474.000000
TRT01ADummy A:AVISITN28                        -0.842654   0.064190 470.300000
TRT01ADummy B:AVISITN28                        -0.694904   0.089841 470.400000
TRT01ADummy C:AVISITN28                        -0.843249   0.090264 471.200000
TRT01ADummy A:AVISITN999                       -0.847770   0.064133 470.400000
TRT01ADummy B:AVISITN999                       -0.702906   0.089752 470.300000
TRT01ADummy C:AVISITN999                       -0.850659   0.090181 471.200000
                                              t value Pr(>|t|)    
(Intercept)                                     4.810 2.11e-06 ***
TRT01ADummy B                                   0.042  0.96668    
TRT01ADummy C                                   1.150  0.25070    
SITEIDUS10002                                   0.617  0.53774    
SITEIDUS10004                                  -1.285  0.19959    
SITEIDUS10005                                   0.659  0.51054    
SITEIDUS10006                                  -1.151  0.25054    
SITEIDUS10008                                  -2.933  0.00354 ** 
SITEIDUS10009                                  -0.809  0.41902    
SITEIDUS10013                                   0.818  0.41359    
SITEIDUS10014                                   0.887  0.37568    
SITEIDUS10015                                   1.273  0.20368    
SITEIDUS10016                                  -1.075  0.28299    
SITEIDUS10017                                  -0.799  0.42497    
SITEIDUS10019                                  -0.566  0.57189    
SITEIDUS10020                                  -2.573  0.01042 *  
SITEIDUS10021                                  -1.344  0.17965    
SITEIDUS10023                                  -2.405  0.01663 *  
SITEIDUS10025                                   0.191  0.84880    
SITEIDUS10027                                  -0.685  0.49395    
SITEIDUS10029                                   0.331  0.74050    
SITEIDUS10030                                  -0.608  0.54384    
SITEIDUS10031                                  -0.356  0.72187    
SITEIDUS10033                                   0.020  0.98368    
SITEIDUS10034                                   0.226  0.82098    
SITEIDUS10036                                   2.394  0.01711 *  
SITEIDUS10039                                   1.577  0.11546    
SITEIDUS10040                                   0.333  0.73965    
SITEIDUS10041                                   0.825  0.40959    
SITEIDUS10043                                  -1.491  0.13671    
SITEIDUS10044                                  -0.387  0.69904    
SITEIDUS10045                                  -0.075  0.94007    
SITEIDUS10046                                   1.152  0.24980    
SITEIDUS10049                                   0.703  0.48239    
SITEIDUS10051                                  -0.616  0.53825    
SITEIDUS10053                                   1.018  0.30915    
SITEIDUS10055                                   0.104  0.91706    
SITEIDUS10056                                  -2.972  0.00313 ** 
SITEIDUS10057                                   0.500  0.61724    
SITEIDUS10058                                  -1.066  0.28726    
SITEIDUS10062                                  -0.074  0.94107    
SITEIDUS10068                                   0.653  0.51418    
SITEIDUS10069                                   1.206  0.22841    
SITEIDUS10070                                  -1.068  0.28608    
SITEIDUS10071                                  -1.746  0.08158 .  
SITEIDUS10073                                  -0.243  0.80803    
SITEIDUS10074                                  -1.121  0.26279    
SITEIDUS10075                                  -0.035  0.97184    
SITEIDUS10076                                  -2.428  0.01560 *  
RACEASIAN                                      -0.743  0.45790    
RACEBLACK OR AFRICAN AMERICAN                  -1.519  0.12963    
RACEMULTIPLE                                   -1.649  0.09983 .  
RACENATIVE HAWAIIAN OR OTHER PACIFIC ISLANDER  -0.971  0.33222    
RACENOT REPORTED                               -0.653  0.51387    
RACEUNKNOWN                                    -0.929  0.35339    
RACEWHITE                                      -1.260  0.20853    
AGE                                             1.195  0.23277    
base_situation10-19                           -12.457  < 2e-16 ***
base_situation20-39                           -21.030  < 2e-16 ***
TRT01ADummy A:AVISITN2                        -16.977  < 2e-16 ***
TRT01ADummy B:AVISITN2                        -11.280  < 2e-16 ***
TRT01ADummy C:AVISITN2                        -11.730  < 2e-16 ***
TRT01ADummy A:AVISITN8                        -11.614  < 2e-16 ***
TRT01ADummy B:AVISITN8                         -7.271 1.53e-12 ***
TRT01ADummy C:AVISITN8                         -9.056  < 2e-16 ***
TRT01ADummy A:AVISITN15                       -11.766  < 2e-16 ***
TRT01ADummy B:AVISITN15                        -7.091 5.04e-12 ***
TRT01ADummy C:AVISITN15                        -8.406 5.23e-16 ***
TRT01ADummy A:AVISITN22                       -12.643  < 2e-16 ***
TRT01ADummy B:AVISITN22                        -8.173 2.82e-15 ***
TRT01ADummy C:AVISITN22                        -9.157  < 2e-16 ***
TRT01ADummy A:AVISITN28                       -13.127  < 2e-16 ***
TRT01ADummy B:AVISITN28                        -7.735 6.36e-14 ***
TRT01ADummy C:AVISITN28                        -9.342  < 2e-16 ***
TRT01ADummy A:AVISITN999                      -13.219  < 2e-16 ***
TRT01ADummy B:AVISITN999                       -7.832 3.23e-14 ***
TRT01ADummy C:AVISITN999                       -9.433  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Covariance estimate:
         0      2      8     15     22     28    999
0   0.1287 0.0874 0.0944 0.0802 0.0866 0.0829 0.0830
2   0.0874 0.8628 0.5757 0.5341 0.5563 0.5582 0.5632
8   0.0944 0.5757 0.7903 0.6412 0.6584 0.6251 0.6278
15  0.0802 0.5341 0.6412 0.8388 0.7806 0.7324 0.7313
22  0.0866 0.5563 0.6584 0.7806 1.0243 0.8904 0.8913
28  0.0829 0.5582 0.6251 0.7324 0.8904 1.0050 1.0005
999 0.0830 0.5632 0.6278 0.7313 0.8913 1.0005 1.0039

draw ci

library(ggplot2)
library(dplyr)
# 提取系数估计值和标准误差，并筛选出包含 "SITE" 的系数
coef_table <- as.data.frame(model_summary$coefficients)


# 计算置信区间
coef_table$LowerCI <- coef_table$Estimate - 1.96 * coef_table$`Std. Error`
coef_table$UpperCI <- coef_table$Estimate + 1.96 * coef_table$`Std. Error`


# 美化P值显示
format_pvalue <- function(p) {
  ifelse(p < 0.001, "<0.001", 
         ifelse(p < 0.01, sprintf("%.3f", p),
                sprintf("%.2f", p)))
}

result_table <- data.frame(
  variable = rownames(model_summary$coefficients),
  Coefficient = round(coef_table$Estimate, 3),
  CI = paste0("(", round(coef_table$LowerCI, 3), ", ", round(coef_table$UpperCI, 3), ")"),
  `P-value` = sapply(coef_table$`Pr(>|t|)`, format_pvalue)
)

result_table
                                        variable Coefficient               CI
1                                    (Intercept)       1.806    (1.07, 2.542)
2                                  TRT01ADummy B       0.002  (-0.079, 0.082)
3                                  TRT01ADummy C       0.048   (-0.034, 0.13)
4                                  SITEIDUS10002       0.077  (-0.168, 0.321)
5                                  SITEIDUS10004      -0.221  (-0.558, 0.116)
6                                  SITEIDUS10005       0.238  (-0.471, 0.948)
7                                  SITEIDUS10006      -0.143    (-0.385, 0.1)
8                                  SITEIDUS10008      -0.555 (-0.926, -0.184)
9                                  SITEIDUS10009      -0.109  (-0.372, 0.155)
10                                 SITEIDUS10013       0.099  (-0.138, 0.335)
11                                 SITEIDUS10014       0.120  (-0.145, 0.385)
12                                 SITEIDUS10015       0.154  (-0.083, 0.391)
13                                 SITEIDUS10016      -0.281  (-0.794, 0.231)
14                                 SITEIDUS10017      -0.127  (-0.438, 0.184)
15                                 SITEIDUS10019      -0.076   (-0.34, 0.188)
16                                 SITEIDUS10020      -0.359 (-0.633, -0.086)
17                                 SITEIDUS10021      -0.157  (-0.386, 0.072)
18                                 SITEIDUS10023      -0.623  (-1.13, -0.115)
19                                 SITEIDUS10025       0.030  (-0.281, 0.341)
20                                 SITEIDUS10027      -0.178  (-0.688, 0.332)
21                                 SITEIDUS10029       0.036  (-0.177, 0.249)
22                                 SITEIDUS10030      -0.062  (-0.261, 0.137)
23                                 SITEIDUS10031      -0.054  (-0.351, 0.243)
24                                 SITEIDUS10033       0.003  (-0.302, 0.308)
25                                 SITEIDUS10034       0.039  (-0.297, 0.375)
26                                 SITEIDUS10036       0.301   (0.055, 0.548)
27                                 SITEIDUS10039       0.569  (-0.138, 1.277)
28                                 SITEIDUS10040       0.063  (-0.311, 0.438)
29                                 SITEIDUS10041       0.130  (-0.179, 0.439)
30                                 SITEIDUS10043      -0.174  (-0.402, 0.055)
31                                 SITEIDUS10044      -0.084   (-0.507, 0.34)
32                                 SITEIDUS10045      -0.014  (-0.384, 0.356)
33                                 SITEIDUS10046       0.183  (-0.129, 0.495)
34                                 SITEIDUS10049       0.111  (-0.199, 0.421)
35                                 SITEIDUS10051      -0.070  (-0.295, 0.154)
36                                 SITEIDUS10053       0.393  (-0.363, 1.148)
37                                 SITEIDUS10055       0.013  (-0.228, 0.253)
38                                 SITEIDUS10056      -0.386 (-0.641, -0.131)
39                                 SITEIDUS10057       0.068    (-0.2, 0.336)
40                                 SITEIDUS10058      -0.129  (-0.366, 0.108)
41                                 SITEIDUS10062      -0.007  (-0.181, 0.168)
42                                 SITEIDUS10068       0.069  (-0.139, 0.278)
43                                 SITEIDUS10069       0.260  (-0.163, 0.683)
44                                 SITEIDUS10070      -0.094  (-0.265, 0.078)
45                                 SITEIDUS10071      -0.280  (-0.595, 0.034)
46                                 SITEIDUS10073      -0.088   (-0.796, 0.62)
47                                 SITEIDUS10074      -0.125  (-0.344, 0.094)
48                                 SITEIDUS10075      -0.003  (-0.177, 0.171)
49                                 SITEIDUS10076      -0.361  (-0.653, -0.07)
50                                     RACEASIAN      -0.284  (-1.031, 0.464)
51                 RACEBLACK OR AFRICAN AMERICAN      -0.567  (-1.299, 0.165)
52                                  RACEMULTIPLE      -0.648  (-1.418, 0.122)
53 RACENATIVE HAWAIIAN OR OTHER PACIFIC ISLANDER      -0.437  (-1.319, 0.445)
54                              RACENOT REPORTED      -0.271  (-1.084, 0.542)
55                                   RACEUNKNOWN      -0.484  (-1.506, 0.538)
56                                     RACEWHITE      -0.463  (-1.184, 0.258)
57                                           AGE       0.022  (-0.014, 0.057)
58                           base_situation10-19      -2.340 (-2.708, -1.972)
59                           base_situation20-39      -0.857 (-0.937, -0.777)
60                        TRT01ADummy A:AVISITN2      -1.006  (-1.122, -0.89)
61                        TRT01ADummy B:AVISITN2      -0.935 (-1.097, -0.772)
62                        TRT01ADummy C:AVISITN2      -0.981 (-1.145, -0.817)
63                        TRT01ADummy A:AVISITN8      -0.651 (-0.761, -0.541)
64                        TRT01ADummy B:AVISITN8      -0.570 (-0.723, -0.416)
65                        TRT01ADummy C:AVISITN8      -0.717 (-0.872, -0.562)
66                       TRT01ADummy A:AVISITN15      -0.695  (-0.81, -0.579)
67                       TRT01ADummy B:AVISITN15      -0.584 (-0.745, -0.422)
68                       TRT01ADummy C:AVISITN15      -0.697  (-0.86, -0.535)
69                       TRT01ADummy A:AVISITN22      -0.822  (-0.95, -0.695)
70                       TRT01ADummy B:AVISITN22      -0.743 (-0.921, -0.565)
71                       TRT01ADummy C:AVISITN22      -0.839 (-1.018, -0.659)
72                       TRT01ADummy A:AVISITN28      -0.843 (-0.968, -0.717)
73                       TRT01ADummy B:AVISITN28      -0.695 (-0.871, -0.519)
74                       TRT01ADummy C:AVISITN28      -0.843  (-1.02, -0.666)
75                      TRT01ADummy A:AVISITN999      -0.848 (-0.973, -0.722)
76                      TRT01ADummy B:AVISITN999      -0.703 (-0.879, -0.527)
77                      TRT01ADummy C:AVISITN999      -0.851 (-1.027, -0.674)
   P.value
1   <0.001
2     0.97
3     0.25
4     0.54
5     0.20
6     0.51
7     0.25
8    0.004
9     0.42
10    0.41
11    0.38
12    0.20
13    0.28
14    0.42
15    0.57
16    0.01
17    0.18
18    0.02
19    0.85
20    0.49
21    0.74
22    0.54
23    0.72
24    0.98
25    0.82
26    0.02
27    0.12
28    0.74
29    0.41
30    0.14
31    0.70
32    0.94
33    0.25
34    0.48
35    0.54
36    0.31
37    0.92
38   0.003
39    0.62
40    0.29
41    0.94
42    0.51
43    0.23
44    0.29
45    0.08
46    0.81
47    0.26
48    0.97
49    0.02
50    0.46
51    0.13
52    0.10
53    0.33
54    0.51
55    0.35
56    0.21
57    0.23
58  <0.001
59  <0.001
60  <0.001
61  <0.001
62  <0.001
63  <0.001
64  <0.001
65  <0.001
66  <0.001
67  <0.001
68  <0.001
69  <0.001
70  <0.001
71  <0.001
72  <0.001
73  <0.001
74  <0.001
75  <0.001
76  <0.001
77  <0.001
# 导出到Excel
#library(writexl)
#write_xlsx(result_table, "mmrm_results.xlsx")

main result with controling sex

library(readxl)
library(mmrm)
library(ggplot2)
library(dplyr)


admadrs <- read_excel("3.xlsx", 
                      sheet = "Sheet1", col_names = TRUE, 
                      col_types = NULL, na = "", skip = 0)
# 创建分类变量 base_situation
admadrs$base_situation <- cut(
  admadrs$BASE,
  breaks = c( 9, 19, 39, Inf),
  labels = c("10-19", "20-39", "40+"),

)



# 将 base_situation 转换为因子类型
admadrs$base_situation <- as.factor(admadrs$base_situation)

admadrs$USUBJID <- factor(admadrs$USUBJID)
admadrs$AVISITN <- admadrs$AVISITN - 20000
admadrs$AVISITN <- factor(admadrs$AVISITN)
admadrs$CHG <- as.numeric(admadrs$CHG)

admadrs$base_situation <- relevel(admadrs$base_situation, ref = "40+")

# Delete a duplicate point in time
admadrs <- admadrs[!duplicated(admadrs[c("USUBJID", "AVISIT")]), ]

# scale data 
admadrs$AGE <- scale(admadrs$AGE)
admadrs$AVAL <- scale(admadrs$AVAL)




# 更新模型以包括 base_situation
fit <- mmrm(
  formula = AVAL ~ AVISITN:TRT01A +TRT01A +SEX+ SITEID+ RACE + AGE + base_situation + us(AVISITN | USUBJID),
  data = admadrs)

# 获取模型摘要
model_summary <- summary(fit)

model_summary
mmrm fit

Formula:     AVAL ~ AVISITN:TRT01A + TRT01A + SEX + SITEID + RACE + AGE +  
    base_situation + us(AVISITN | USUBJID)
Data:        admadrs (used 3214 observations from 476 subjects with maximum 7 
timepoints)
Covariance:  unstructured (28 variance parameters)
Method:      Satterthwaite
Vcov Method: Asymptotic
Inference:   REML

Model selection criteria:
     AIC      BIC   logLik deviance 
  3856.3   3972.9  -1900.1   3800.3 

Coefficients: 
                                                Estimate Std. Error         df
(Intercept)                                    1.801e+00  3.756e-01  4.145e+02
TRT01ADummy B                                  1.714e-03  4.114e-02  4.150e+02
TRT01ADummy C                                  4.943e-02  4.178e-02  4.159e+02
SEXM                                          -3.458e-02  3.568e-02  4.151e+02
SITEIDUS10002                                  8.039e-02  1.248e-01  4.147e+02
SITEIDUS10004                                 -2.118e-01  1.723e-01  4.145e+02
SITEIDUS10005                                  2.211e-01  3.625e-01  4.143e+02
SITEIDUS10006                                 -1.343e-01  1.241e-01  4.147e+02
SITEIDUS10008                                 -5.671e-01  1.897e-01  4.157e+02
SITEIDUS10009                                 -1.067e-01  1.344e-01  4.145e+02
SITEIDUS10013                                  1.027e-01  1.206e-01  4.143e+02
SITEIDUS10014                                  1.153e-01  1.354e-01  4.161e+02
SITEIDUS10015                                  1.455e-01  1.213e-01  4.145e+02
SITEIDUS10016                                 -2.962e-01  2.618e-01  4.200e+02
SITEIDUS10017                                 -1.249e-01  1.587e-01  4.152e+02
SITEIDUS10019                                 -6.638e-02  1.350e-01  4.143e+02
SITEIDUS10020                                 -3.592e-01  1.397e-01  4.144e+02
SITEIDUS10021                                 -1.470e-01  1.174e-01  4.144e+02
SITEIDUS10023                                 -6.167e-01  2.590e-01  4.146e+02
SITEIDUS10025                                  2.665e-02  1.587e-01  4.143e+02
SITEIDUS10027                                 -1.780e-01  2.604e-01  4.143e+02
SITEIDUS10029                                  4.001e-02  1.086e-01  4.146e+02
SITEIDUS10030                                 -5.887e-02  1.016e-01  4.146e+02
SITEIDUS10031                                 -5.395e-02  1.517e-01  4.143e+02
SITEIDUS10033                                  4.711e-04  1.558e-01  4.321e+02
SITEIDUS10034                                  3.307e-02  1.716e-01  4.143e+02
SITEIDUS10036                                  3.019e-01  1.258e-01  4.143e+02
SITEIDUS10039                                  5.580e-01  3.613e-01  4.143e+02
SITEIDUS10040                                  7.344e-02  1.912e-01  4.145e+02
SITEIDUS10041                                  1.468e-01  1.584e-01  4.145e+02
SITEIDUS10043                                 -1.695e-01  1.165e-01  4.149e+02
SITEIDUS10044                                 -9.636e-02  2.164e-01  4.147e+02
SITEIDUS10045                                 -1.581e-02  1.888e-01  4.143e+02
SITEIDUS10046                                  1.849e-01  1.592e-01  4.147e+02
SITEIDUS10049                                  1.036e-01  1.582e-01  4.143e+02
SITEIDUS10051                                 -8.296e-02  1.149e-01  4.143e+02
SITEIDUS10053                                  4.168e-01  3.866e-01  4.144e+02
SITEIDUS10055                                  1.125e-02  1.228e-01  4.145e+02
SITEIDUS10056                                 -3.738e-01  1.305e-01  4.151e+02
SITEIDUS10057                                  6.903e-02  1.368e-01  4.146e+02
SITEIDUS10058                                 -1.250e-01  1.210e-01  4.144e+02
SITEIDUS10062                                  1.243e-03  8.932e-02  4.147e+02
SITEIDUS10068                                  7.517e-02  1.066e-01  4.144e+02
SITEIDUS10069                                  2.575e-01  2.159e-01  4.143e+02
SITEIDUS10070                                 -9.618e-02  8.767e-02  4.146e+02
SITEIDUS10071                                 -2.820e-01  1.605e-01  4.145e+02
SITEIDUS10073                                 -1.095e-01  3.614e-01  4.142e+02
SITEIDUS10074                                 -1.269e-01  1.118e-01  4.149e+02
SITEIDUS10075                                 -6.080e-03  8.878e-02  4.145e+02
SITEIDUS10076                                 -3.652e-01  1.489e-01  4.149e+02
RACEASIAN                                     -2.722e-01  3.818e-01  4.145e+02
RACEBLACK OR AFRICAN AMERICAN                 -5.534e-01  3.736e-01  4.144e+02
RACEMULTIPLE                                  -6.407e-01  3.931e-01  4.146e+02
RACENATIVE HAWAIIAN OR OTHER PACIFIC ISLANDER -4.206e-01  4.507e-01  4.146e+02
RACENOT REPORTED                              -2.520e-01  4.152e-01  4.145e+02
RACEUNKNOWN                                   -4.518e-01  5.227e-01  4.146e+02
RACEWHITE                                     -4.490e-01  3.680e-01  4.144e+02
AGE                                            2.227e-02  1.815e-02  4.169e+02
base_situation10-19                           -2.334e+00  1.881e-01  4.146e+02
base_situation20-39                           -8.527e-01  4.104e-02  4.146e+02
TRT01ADummy A:AVISITN2                        -1.006e+00  5.927e-02  4.685e+02
TRT01ADummy B:AVISITN2                        -9.346e-01  8.285e-02  4.673e+02
TRT01ADummy C:AVISITN2                        -9.810e-01  8.363e-02  4.725e+02
TRT01ADummy A:AVISITN8                        -6.510e-01  5.605e-02  4.635e+02
TRT01ADummy B:AVISITN8                        -5.695e-01  7.832e-02  4.621e+02
TRT01ADummy C:AVISITN8                        -7.166e-01  7.915e-02  4.692e+02
TRT01ADummy A:AVISITN15                       -6.947e-01  5.904e-02  4.645e+02
TRT01ADummy B:AVISITN15                       -5.836e-01  8.230e-02  4.594e+02
TRT01ADummy C:AVISITN15                       -6.973e-01  8.296e-02  4.642e+02
TRT01ADummy A:AVISITN22                       -8.222e-01  6.503e-02  4.706e+02
TRT01ADummy B:AVISITN22                       -7.430e-01  9.089e-02  4.687e+02
TRT01ADummy C:AVISITN22                       -8.386e-01  9.160e-02  4.741e+02
TRT01ADummy A:AVISITN28                       -8.426e-01  6.419e-02  4.703e+02
TRT01ADummy B:AVISITN28                       -6.949e-01  8.984e-02  4.704e+02
TRT01ADummy C:AVISITN28                       -8.431e-01  9.026e-02  4.712e+02
TRT01ADummy A:AVISITN999                      -8.478e-01  6.413e-02  4.704e+02
TRT01ADummy B:AVISITN999                      -7.029e-01  8.975e-02  4.703e+02
TRT01ADummy C:AVISITN999                      -8.505e-01  9.018e-02  4.712e+02
                                              t value Pr(>|t|)    
(Intercept)                                     4.797 2.25e-06 ***
TRT01ADummy B                                   0.042  0.96678    
TRT01ADummy C                                   1.183  0.23741    
SEXM                                           -0.969  0.33304    
SITEIDUS10002                                   0.644  0.51993    
SITEIDUS10004                                  -1.229  0.21977    
SITEIDUS10005                                   0.610  0.54227    
SITEIDUS10006                                  -1.082  0.27979    
SITEIDUS10008                                  -2.989  0.00297 ** 
SITEIDUS10009                                  -0.793  0.42801    
SITEIDUS10013                                   0.851  0.39515    
SITEIDUS10014                                   0.851  0.39523    
SITEIDUS10015                                   1.200  0.23089    
SITEIDUS10016                                  -1.131  0.25857    
SITEIDUS10017                                  -0.787  0.43159    
SITEIDUS10019                                  -0.492  0.62330    
SITEIDUS10020                                  -2.570  0.01050 *  
SITEIDUS10021                                  -1.251  0.21150    
SITEIDUS10023                                  -2.381  0.01772 *  
SITEIDUS10025                                   0.168  0.86671    
SITEIDUS10027                                  -0.684  0.49467    
SITEIDUS10029                                   0.368  0.71276    
SITEIDUS10030                                  -0.580  0.56242    
SITEIDUS10031                                  -0.356  0.72223    
SITEIDUS10033                                   0.003  0.99759    
SITEIDUS10034                                   0.193  0.84725    
SITEIDUS10036                                   2.399  0.01686 *  
SITEIDUS10039                                   1.544  0.12327    
SITEIDUS10040                                   0.384  0.70103    
SITEIDUS10041                                   0.927  0.35462    
SITEIDUS10043                                  -1.455  0.14654    
SITEIDUS10044                                  -0.445  0.65633    
SITEIDUS10045                                  -0.084  0.93330    
SITEIDUS10046                                   1.161  0.24619    
SITEIDUS10049                                   0.655  0.51274    
SITEIDUS10051                                  -0.722  0.47075    
SITEIDUS10053                                   1.078  0.28162    
SITEIDUS10055                                   0.092  0.92704    
SITEIDUS10056                                  -2.865  0.00438 ** 
SITEIDUS10057                                   0.505  0.61416    
SITEIDUS10058                                  -1.033  0.30204    
SITEIDUS10062                                   0.014  0.98890    
SITEIDUS10068                                   0.705  0.48119    
SITEIDUS10069                                   1.193  0.23357    
SITEIDUS10070                                  -1.097  0.27325    
SITEIDUS10071                                  -1.756  0.07974 .  
SITEIDUS10073                                  -0.303  0.76212    
SITEIDUS10074                                  -1.135  0.25699    
SITEIDUS10075                                  -0.068  0.94543    
SITEIDUS10076                                  -2.453  0.01458 *  
RACEASIAN                                      -0.713  0.47634    
RACEBLACK OR AFRICAN AMERICAN                  -1.481  0.13936    
RACEMULTIPLE                                   -1.630  0.10385    
RACENATIVE HAWAIIAN OR OTHER PACIFIC ISLANDER  -0.933  0.35130    
RACENOT REPORTED                               -0.607  0.54414    
RACEUNKNOWN                                    -0.864  0.38792    
RACEWHITE                                      -1.220  0.22316    
AGE                                             1.227  0.22045    
base_situation10-19                           -12.409  < 2e-16 ***
base_situation20-39                           -20.780  < 2e-16 ***
TRT01ADummy A:AVISITN2                        -16.977  < 2e-16 ***
TRT01ADummy B:AVISITN2                        -11.280  < 2e-16 ***
TRT01ADummy C:AVISITN2                        -11.730  < 2e-16 ***
TRT01ADummy A:AVISITN8                        -11.615  < 2e-16 ***
TRT01ADummy B:AVISITN8                         -7.272 1.53e-12 ***
TRT01ADummy C:AVISITN8                         -9.054  < 2e-16 ***
TRT01ADummy A:AVISITN15                       -11.767  < 2e-16 ***
TRT01ADummy B:AVISITN15                        -7.091 5.04e-12 ***
TRT01ADummy C:AVISITN15                        -8.405 5.26e-16 ***
TRT01ADummy A:AVISITN22                       -12.644  < 2e-16 ***
TRT01ADummy B:AVISITN22                        -8.174 2.80e-15 ***
TRT01ADummy C:AVISITN22                        -9.155  < 2e-16 ***
TRT01ADummy A:AVISITN28                       -13.128  < 2e-16 ***
TRT01ADummy B:AVISITN28                        -7.735 6.34e-14 ***
TRT01ADummy C:AVISITN28                        -9.341  < 2e-16 ***
TRT01ADummy A:AVISITN999                      -13.220  < 2e-16 ***
TRT01ADummy B:AVISITN999                       -7.832 3.22e-14 ***
TRT01ADummy C:AVISITN999                       -9.432  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Covariance estimate:
         0      2      8     15     22     28    999
0   0.1287 0.0888 0.0953 0.0803 0.0868 0.0823 0.0825
2   0.0888 0.8656 0.5780 0.5356 0.5578 0.5590 0.5641
8   0.0953 0.5780 0.7920 0.6422 0.6594 0.6253 0.6281
15  0.0803 0.5356 0.6422 0.8390 0.7808 0.7318 0.7309
22  0.0868 0.5578 0.6594 0.7808 1.0246 0.8898 0.8910
28  0.0823 0.5590 0.6253 0.7318 0.8898 1.0037 0.9993
999 0.0825 0.5641 0.6281 0.7309 0.8910 0.9993 1.0028

draw ci

library(ggplot2)
library(dplyr)
# 提取系数估计值和标准误差，并筛选出包含 "SITE" 的系数
coef_table <- as.data.frame(model_summary$coefficients)


# 计算置信区间
coef_table$LowerCI <- coef_table$Estimate - 1.96 * coef_table$`Std. Error`
coef_table$UpperCI <- coef_table$Estimate + 1.96 * coef_table$`Std. Error`


# 美化P值显示
format_pvalue <- function(p) {
  ifelse(p < 0.001, "<0.001", 
         ifelse(p < 0.01, sprintf("%.3f", p),
                sprintf("%.2f", p)))
}

result_table <- data.frame(
  variable = rownames(model_summary$coefficients),
  Coefficient = round(coef_table$Estimate, 3),
  CI = paste0("(", round(coef_table$LowerCI, 3), ", ", round(coef_table$UpperCI, 3), ")"),
  `P-value` = sapply(coef_table$`Pr(>|t|)`, format_pvalue)
)

result_table
                                        variable Coefficient               CI
1                                    (Intercept)       1.801   (1.065, 2.538)
2                                  TRT01ADummy B       0.002  (-0.079, 0.082)
3                                  TRT01ADummy C       0.049  (-0.032, 0.131)
4                                           SEXM      -0.035  (-0.105, 0.035)
5                                  SITEIDUS10002       0.080  (-0.164, 0.325)
6                                  SITEIDUS10004      -0.212   (-0.55, 0.126)
7                                  SITEIDUS10005       0.221  (-0.489, 0.932)
8                                  SITEIDUS10006      -0.134  (-0.377, 0.109)
9                                  SITEIDUS10008      -0.567 (-0.939, -0.195)
10                                 SITEIDUS10009      -0.107   (-0.37, 0.157)
11                                 SITEIDUS10013       0.103  (-0.134, 0.339)
12                                 SITEIDUS10014       0.115   (-0.15, 0.381)
13                                 SITEIDUS10015       0.146  (-0.092, 0.383)
14                                 SITEIDUS10016      -0.296  (-0.809, 0.217)
15                                 SITEIDUS10017      -0.125  (-0.436, 0.186)
16                                 SITEIDUS10019      -0.066  (-0.331, 0.198)
17                                 SITEIDUS10020      -0.359 (-0.633, -0.085)
18                                 SITEIDUS10021      -0.147  (-0.377, 0.083)
19                                 SITEIDUS10023      -0.617 (-1.124, -0.109)
20                                 SITEIDUS10025       0.027  (-0.284, 0.338)
21                                 SITEIDUS10027      -0.178  (-0.688, 0.332)
22                                 SITEIDUS10029       0.040  (-0.173, 0.253)
23                                 SITEIDUS10030      -0.059   (-0.258, 0.14)
24                                 SITEIDUS10031      -0.054  (-0.351, 0.243)
25                                 SITEIDUS10033       0.000  (-0.305, 0.306)
26                                 SITEIDUS10034       0.033  (-0.303, 0.369)
27                                 SITEIDUS10036       0.302   (0.055, 0.549)
28                                 SITEIDUS10039       0.558   (-0.15, 1.266)
29                                 SITEIDUS10040       0.073  (-0.301, 0.448)
30                                 SITEIDUS10041       0.147  (-0.164, 0.457)
31                                 SITEIDUS10043      -0.170  (-0.398, 0.059)
32                                 SITEIDUS10044      -0.096  (-0.521, 0.328)
33                                 SITEIDUS10045      -0.016  (-0.386, 0.354)
34                                 SITEIDUS10046       0.185  (-0.127, 0.497)
35                                 SITEIDUS10049       0.104  (-0.206, 0.414)
36                                 SITEIDUS10051      -0.083  (-0.308, 0.142)
37                                 SITEIDUS10053       0.417  (-0.341, 1.175)
38                                 SITEIDUS10055       0.011  (-0.229, 0.252)
39                                 SITEIDUS10056      -0.374  (-0.63, -0.118)
40                                 SITEIDUS10057       0.069  (-0.199, 0.337)
41                                 SITEIDUS10058      -0.125  (-0.362, 0.112)
42                                 SITEIDUS10062       0.001  (-0.174, 0.176)
43                                 SITEIDUS10068       0.075  (-0.134, 0.284)
44                                 SITEIDUS10069       0.258  (-0.166, 0.681)
45                                 SITEIDUS10070      -0.096  (-0.268, 0.076)
46                                 SITEIDUS10071      -0.282  (-0.597, 0.033)
47                                 SITEIDUS10073      -0.109  (-0.818, 0.599)
48                                 SITEIDUS10074      -0.127  (-0.346, 0.092)
49                                 SITEIDUS10075      -0.006   (-0.18, 0.168)
50                                 SITEIDUS10076      -0.365 (-0.657, -0.073)
51                                     RACEASIAN      -0.272  (-1.021, 0.476)
52                 RACEBLACK OR AFRICAN AMERICAN      -0.553  (-1.286, 0.179)
53                                  RACEMULTIPLE      -0.641   (-1.411, 0.13)
54 RACENATIVE HAWAIIAN OR OTHER PACIFIC ISLANDER      -0.421  (-1.304, 0.463)
55                              RACENOT REPORTED      -0.252  (-1.066, 0.562)
56                                   RACEUNKNOWN      -0.452  (-1.476, 0.573)
57                                     RACEWHITE      -0.449   (-1.17, 0.272)
58                                           AGE       0.022  (-0.013, 0.058)
59                           base_situation10-19      -2.334 (-2.702, -1.965)
60                           base_situation20-39      -0.853 (-0.933, -0.772)
61                        TRT01ADummy A:AVISITN2      -1.006  (-1.122, -0.89)
62                        TRT01ADummy B:AVISITN2      -0.935 (-1.097, -0.772)
63                        TRT01ADummy C:AVISITN2      -0.981 (-1.145, -0.817)
64                        TRT01ADummy A:AVISITN8      -0.651 (-0.761, -0.541)
65                        TRT01ADummy B:AVISITN8      -0.570 (-0.723, -0.416)
66                        TRT01ADummy C:AVISITN8      -0.717 (-0.872, -0.561)
67                       TRT01ADummy A:AVISITN15      -0.695  (-0.81, -0.579)
68                       TRT01ADummy B:AVISITN15      -0.584 (-0.745, -0.422)
69                       TRT01ADummy C:AVISITN15      -0.697  (-0.86, -0.535)
70                       TRT01ADummy A:AVISITN22      -0.822  (-0.95, -0.695)
71                       TRT01ADummy B:AVISITN22      -0.743 (-0.921, -0.565)
72                       TRT01ADummy C:AVISITN22      -0.839 (-1.018, -0.659)
73                       TRT01ADummy A:AVISITN28      -0.843 (-0.968, -0.717)
74                       TRT01ADummy B:AVISITN28      -0.695 (-0.871, -0.519)
75                       TRT01ADummy C:AVISITN28      -0.843  (-1.02, -0.666)
76                      TRT01ADummy A:AVISITN999      -0.848 (-0.973, -0.722)
77                      TRT01ADummy B:AVISITN999      -0.703 (-0.879, -0.527)
78                      TRT01ADummy C:AVISITN999      -0.850 (-1.027, -0.674)
   P.value
1   <0.001
2     0.97
3     0.24
4     0.33
5     0.52
6     0.22
7     0.54
8     0.28
9    0.003
10    0.43
11    0.40
12    0.40
13    0.23
14    0.26
15    0.43
16    0.62
17    0.01
18    0.21
19    0.02
20    0.87
21    0.49
22    0.71
23    0.56
24    0.72
25    1.00
26    0.85
27    0.02
28    0.12
29    0.70
30    0.35
31    0.15
32    0.66
33    0.93
34    0.25
35    0.51
36    0.47
37    0.28
38    0.93
39   0.004
40    0.61
41    0.30
42    0.99
43    0.48
44    0.23
45    0.27
46    0.08
47    0.76
48    0.26
49    0.95
50    0.01
51    0.48
52    0.14
53    0.10
54    0.35
55    0.54
56    0.39
57    0.22
58    0.22
59  <0.001
60  <0.001
61  <0.001
62  <0.001
63  <0.001
64  <0.001
65  <0.001
66  <0.001
67  <0.001
68  <0.001
69  <0.001
70  <0.001
71  <0.001
72  <0.001
73  <0.001
74  <0.001
75  <0.001
76  <0.001
77  <0.001
78  <0.001
# 导出到Excel
#library(writexl)
#write_xlsx(result_table, "mmrm_results.xlsx")

lm

# 使用 lm 拟合线性模型
fit_lm <- lm(
  formula = AVAL ~ AVISITN:TRT01A +TRT01A + SITEID+ RACE + AGE + base_situation,
  data = admadrs
)

summary(fit_lm)

Call:
lm(formula = AVAL ~ AVISITN:TRT01A + TRT01A + SITEID + RACE + 
    AGE + base_situation, data = admadrs)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.7612 -0.4780  0.1187  0.5814  2.1366 

Coefficients:
                                              Estimate Std. Error t value
(Intercept)                                    1.71542    0.37155   4.617
TRT01ADummy B                                 -0.02879    0.09610  -0.300
TRT01ADummy C                                  0.04646    0.09678   0.480
SITEIDUS10002                                  0.02089    0.11608   0.180
SITEIDUS10004                                  0.12149    0.15902   0.764
SITEIDUS10005                                  0.24366    0.33050   0.737
SITEIDUS10006                                 -0.17280    0.11428  -1.512
SITEIDUS10008                                 -0.34969    0.19240  -1.817
SITEIDUS10009                                  0.24909    0.12509   1.991
SITEIDUS10013                                  0.56116    0.11016   5.094
SITEIDUS10014                                  0.17056    0.12858   1.326
SITEIDUS10015                                  0.03213    0.11271   0.285
SITEIDUS10016                                 -0.30379    0.30909  -0.983
SITEIDUS10017                                 -0.04176    0.15110  -0.276
SITEIDUS10019                                 -0.15188    0.12295  -1.235
SITEIDUS10020                                 -0.38872    0.12855  -3.024
SITEIDUS10021                                  0.27363    0.10667   2.565
SITEIDUS10023                                 -0.29919    0.25405  -1.178
SITEIDUS10025                                  0.31938    0.14496   2.203
SITEIDUS10027                                 -0.33069    0.23768  -1.391
SITEIDUS10029                                  0.03777    0.10003   0.378
SITEIDUS10030                                  0.01878    0.09352   0.201
SITEIDUS10031                                 -0.26735    0.13852  -1.930
SITEIDUS10033                                  0.19930    0.13631   1.462
SITEIDUS10034                                  0.30245    0.15667   1.931
SITEIDUS10036                                 -0.36938    0.11498  -3.213
SITEIDUS10039                                  1.00520    0.32949   3.051
SITEIDUS10040                                  0.12201    0.17742   0.688
SITEIDUS10041                                  0.15199    0.14525   1.046
SITEIDUS10043                                 -0.28370    0.11056  -2.566
SITEIDUS10044                                  0.26565    0.20138   1.319
SITEIDUS10045                                 -0.13825    0.17227  -0.803
SITEIDUS10046                                 -0.09002    0.14859  -0.606
SITEIDUS10049                                  0.28183    0.14436   1.952
SITEIDUS10051                                 -0.09597    0.10495  -0.914
SITEIDUS10053                                 -0.55402    0.35330  -1.568
SITEIDUS10055                                 -0.42579    0.11334  -3.757
SITEIDUS10056                                 -0.31387    0.12409  -2.529
SITEIDUS10057                                  0.47336    0.12656   3.740
SITEIDUS10058                                  0.15978    0.11139   1.434
SITEIDUS10062                                 -0.12712    0.08208  -1.549
SITEIDUS10068                                 -0.30722    0.09731  -3.157
SITEIDUS10069                                  0.79177    0.19707   4.018
SITEIDUS10070                                 -0.02214    0.08083  -0.274
SITEIDUS10071                                 -0.13676    0.14748  -0.927
SITEIDUS10073                                 -0.51585    0.35468  -1.454
SITEIDUS10074                                 -0.15776    0.10491  -1.504
SITEIDUS10075                                 -0.07545    0.08233  -0.916
SITEIDUS10076                                  0.20851    0.13704   1.522
RACEASIAN                                     -0.08511    0.37347  -0.228
RACEBLACK OR AFRICAN AMERICAN                 -0.38612    0.36620  -1.054
RACEMULTIPLE                                  -0.70294    0.38410  -1.830
RACENATIVE HAWAIIAN OR OTHER PACIFIC ISLANDER -0.24321    0.43215  -0.563
RACENOT REPORTED                               0.08491    0.40171   0.211
RACEUNKNOWN                                    0.60616    0.49442   1.226
RACEWHITE                                     -0.43075    0.36130  -1.192
AGE                                            0.07715    0.01666   4.631
base_situation10-19                           -1.69274    0.17404  -9.726
base_situation20-39                           -0.82355    0.03766 -21.871
AVISITN2:TRT01ADummy A                        -1.00458    0.07913 -12.696
AVISITN8:TRT01ADummy A                        -0.65815    0.07929  -8.301
AVISITN15:TRT01ADummy A                       -0.69466    0.07974  -8.711
AVISITN22:TRT01ADummy A                       -0.82964    0.08014 -10.352
AVISITN28:TRT01ADummy A                       -0.86057    0.07967 -10.802
AVISITN999:TRT01ADummy A                      -0.84786    0.07867 -10.777
AVISITN2:TRT01ADummy B                        -0.92254    0.11080  -8.326
AVISITN8:TRT01ADummy B                        -0.57217    0.11103  -5.153
AVISITN15:TRT01ADummy B                       -0.58461    0.11103  -5.265
AVISITN22:TRT01ADummy B                       -0.75172    0.11206  -6.708
AVISITN28:TRT01ADummy B                       -0.72849    0.11229  -6.488
AVISITN999:TRT01ADummy B                      -0.70291    0.11031  -6.372
AVISITN2:TRT01ADummy C                        -0.96668    0.11200  -8.631
AVISITN8:TRT01ADummy C                        -0.72971    0.11249  -6.487
AVISITN15:TRT01ADummy C                       -0.70641    0.11223  -6.294
AVISITN22:TRT01ADummy C                       -0.87271    0.11358  -7.684
AVISITN28:TRT01ADummy C                       -0.88039    0.11251  -7.825
AVISITN999:TRT01ADummy C                      -0.85249    0.11102  -7.679
                                              Pr(>|t|)    
(Intercept)                                   4.05e-06 ***
TRT01ADummy B                                 0.764480    
TRT01ADummy C                                 0.631244    
SITEIDUS10002                                 0.857208    
SITEIDUS10004                                 0.444936    
SITEIDUS10005                                 0.461023    
SITEIDUS10006                                 0.130631    
SITEIDUS10008                                 0.069240 .  
SITEIDUS10009                                 0.046537 *  
SITEIDUS10013                                 3.71e-07 ***
SITEIDUS10014                                 0.184772    
SITEIDUS10015                                 0.775630    
SITEIDUS10016                                 0.325749    
SITEIDUS10017                                 0.782303    
SITEIDUS10019                                 0.216802    
SITEIDUS10020                                 0.002515 ** 
SITEIDUS10021                                 0.010359 *  
SITEIDUS10023                                 0.239007    
SITEIDUS10025                                 0.027647 *  
SITEIDUS10027                                 0.164221    
SITEIDUS10029                                 0.705780    
SITEIDUS10030                                 0.840871    
SITEIDUS10031                                 0.053694 .  
SITEIDUS10033                                 0.143797    
SITEIDUS10034                                 0.053631 .  
SITEIDUS10036                                 0.001328 ** 
SITEIDUS10039                                 0.002302 ** 
SITEIDUS10040                                 0.491681    
SITEIDUS10041                                 0.295448    
SITEIDUS10043                                 0.010332 *  
SITEIDUS10044                                 0.187200    
SITEIDUS10045                                 0.422316    
SITEIDUS10046                                 0.544646    
SITEIDUS10049                                 0.050997 .  
SITEIDUS10051                                 0.360574    
SITEIDUS10053                                 0.116949    
SITEIDUS10055                                 0.000175 ***
SITEIDUS10056                                 0.011475 *  
SITEIDUS10057                                 0.000187 ***
SITEIDUS10058                                 0.151531    
SITEIDUS10062                                 0.121524    
SITEIDUS10068                                 0.001609 ** 
SITEIDUS10069                                 6.01e-05 ***
SITEIDUS10070                                 0.784170    
SITEIDUS10071                                 0.353804    
SITEIDUS10073                                 0.145940    
SITEIDUS10074                                 0.132736    
SITEIDUS10075                                 0.359513    
SITEIDUS10076                                 0.128221    
RACEASIAN                                     0.819740    
RACEBLACK OR AFRICAN AMERICAN                 0.291789    
RACEMULTIPLE                                  0.067334 .  
RACENATIVE HAWAIIAN OR OTHER PACIFIC ISLANDER 0.573612    
RACENOT REPORTED                              0.832604    
RACEUNKNOWN                                   0.220285    
RACEWHITE                                     0.233258    
AGE                                           3.78e-06 ***
base_situation10-19                            < 2e-16 ***
base_situation20-39                            < 2e-16 ***
AVISITN2:TRT01ADummy A                         < 2e-16 ***
AVISITN8:TRT01ADummy A                         < 2e-16 ***
AVISITN15:TRT01ADummy A                        < 2e-16 ***
AVISITN22:TRT01ADummy A                        < 2e-16 ***
AVISITN28:TRT01ADummy A                        < 2e-16 ***
AVISITN999:TRT01ADummy A                       < 2e-16 ***
AVISITN2:TRT01ADummy B                         < 2e-16 ***
AVISITN8:TRT01ADummy B                        2.72e-07 ***
AVISITN15:TRT01ADummy B                       1.49e-07 ***
AVISITN22:TRT01ADummy B                       2.33e-11 ***
AVISITN28:TRT01ADummy B                       1.01e-10 ***
AVISITN999:TRT01ADummy B                      2.14e-10 ***
AVISITN2:TRT01ADummy C                         < 2e-16 ***
AVISITN8:TRT01ADummy C                        1.01e-10 ***
AVISITN15:TRT01ADummy C                       3.52e-10 ***
AVISITN22:TRT01ADummy C                       2.05e-14 ***
AVISITN28:TRT01ADummy C                       6.87e-15 ***
AVISITN999:TRT01ADummy C                      2.13e-14 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.8545 on 3137 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.2873,    Adjusted R-squared:  0.2701 
F-statistic: 16.64 on 76 and 3137 DF,  p-value: < 2.2e-16
AIC(fit_lm)
[1] 8188.17
BIC(fit_lm)
[1] 8662.041
logLik(fit_lm)
'log Lik.' -4016.085 (df=78)

draw subgroups

≥18 to 24 young adults

# 加载必要的包
library(readxl)
library(mmrm)
library(ggplot2)

?mmrm
# 读取数据
admadrs <- read_excel("/Users/qiujiahui/Desktop/≥18 to 24 young adults.xlsx", 
                      sheet = "Sheet1", col_names = TRUE, 
                      col_types = NULL, na = "", skip = 0)
head(admadrs)

# 删除重复的时间点
admadrs <- admadrs[!duplicated(admadrs[c("USUBJID", "AVISIT")]), ]

# 将相关列转换为因子
admadrs$USUBJID <- factor(admadrs$USUBJID)
admadrs$AVISITN <- factor(admadrs$AVISITN)
admadrs$ADY <- as.factor(admadrs$ADY)
admadrs$SITEID <- as.factor(admadrs$SITEID)
admadrs$TRTSEQA <- as.factor(admadrs$TRTSEQA)
admadrs$BASE <- as.numeric(admadrs$BASE)

# 将 TRT01A 转换为因子，并重新编码参考水平
admadrs$TRT01A <- as.factor(admadrs$TRT01A)
admadrs$TRT01A <- relevel(admadrs$TRT01A, ref = "Dummy A")

#admadrs$AVISITN <- relevel(admadrs$AVISITN, ref = "")

# 对数值型变量进行归一化处理


str(admadrs$BASE)

# 拟合模型
fit <- mmrm(
  formula = AVAL ~ AVISITN:TRT01A + TRT01A+ SEX + RACE + SITEID  + us(AVISITN | USUBJID),
  data = admadrs
)

# 查看模型结果
model_summary <- summary(fit)
print(model_summary)

≥18 to 24 young adults 分析


# 加载必要的包
library(emmeans)
library(ggplot2)
library(dplyr)

############## 年龄组分析 ##############
cat("\n### 年龄组分析 ###\n")

# 重新计算边际均值
emm_young <- emmeans(fit, ~ AVISITN | TRT01A)
print(emm_young)

# 基线和时间点设置
baseline_level <- "0"  # 假设基线是0，请调整为实际值

# 获取AVISITN的水平
avisitn_levels <- levels(admadrs$AVISITN)
time_points <- avisitn_levels[avisitn_levels != baseline_level]

# 如果AVISITN是数值而不是因子
if(length(avisitn_levels) == 0) {
  cat("AVISITN似乎不是因子或没有水平，尝试获取唯一值...\n")
  avisitn_values <- sort(unique(admadrs$AVISITN))
  avisitn_values <- as.character(avisitn_values)
  baseline_level <- as.character(min(as.numeric(avisitn_values)))
  time_points <- avisitn_values[avisitn_values != baseline_level]
}

cat("基线:", baseline_level, "\n")
cat("比较时间点:", time_points, "\n")

# 执行对比分析
tryCatch({
  # 对每个治疗组，计算各时间点与基线的对比
  comparisons_young <- contrast(emm_young, method = "revpairwise", by = "TRT01A")
  cat("\n年轻成人组的时间点比较:\n")
  print(comparisons_young)
  
  # 保存结果
  results_young <- as.data.frame(comparisons_young)
  results_young$AgeGroup <- "Young Adults (18-24)"
  
}, error = function(e) {
  cat("年轻成人组的对比分析出错:", conditionMessage(e), "\n")
  results_young <- data.frame()
})

middle-aged adults (≥24 to 55 years)

# 加载必要的包
library(readxl)
library(mmrm)
library(ggplot2)

?mmrm
# 读取数据
admadrs <- read_excel("/Users/qiujiahui/Desktop/≥24 to 55 middle-aged adults.xlsx", 
                      sheet = "Sheet1", col_names = TRUE, 
                      col_types = NULL, na = "", skip = 0)
head(admadrs)

# 删除重复的时间点
admadrs <- admadrs[!duplicated(admadrs[c("USUBJID", "AVISIT")]), ]

# 将相关列转换为因子
admadrs$USUBJID <- factor(admadrs$USUBJID)
admadrs$AVISITN <- factor(admadrs$AVISITN)
admadrs$ADY <- as.factor(admadrs$ADY)
admadrs$SITEID <- as.factor(admadrs$SITEID)
admadrs$TRTSEQA <- as.factor(admadrs$TRTSEQA)
admadrs$BASE <- as.numeric(admadrs$BASE)

# 将 TRT01A 转换为因子，并重新编码参考水平
admadrs$TRT01A <- as.factor(admadrs$TRT01A)
admadrs$TRT01A <- relevel(admadrs$TRT01A, ref = "Dummy A")

#admadrs$AVISITN <- relevel(admadrs$AVISITN, ref = "")

# 对数值型变量进行归一化处理


str(admadrs$BASE)

# 拟合模型
fit <- mmrm(
  formula = AVAL ~ AVISITN:TRT01A + TRT01A+ SEX + RACE + SITEID  + us(AVISITN | USUBJID),
  data = admadrs
)

# 查看模型结果
model_summary <- summary(fit)
print(model_summary)

middle-aged adults (≥24 to 55 years)分析


############## 中年成人分析 ##############
# 读取中年成人数据并拟合模型
# (假设这里已经完成)

cat("\n### 中年成人组分析 ###\n")

# 重新计算边际均值
emm_middle <- emmeans(fit, ~ AVISITN | TRT01A)
print(emm_middle)

# 执行对比分析
tryCatch({
  # 对每个治疗组，计算各时间点与基线的对比
  comparisons_middle <- contrast(emm_middle, method = "revpairwise", by = "TRT01A")
  cat("\n中年成人组的时间点比较:\n")
  print(comparisons_middle)
  
  # 保存结果
  results_middle <- as.data.frame(comparisons_middle)
  results_middle$AgeGroup <- "Middle-Aged Adults (24-55)"
  
}, error = function(e) {
  cat("中年成人组的对比分析出错:", conditionMessage(e), "\n")
  results_middle <- data.frame()
})

≥55 older adults

# 加载必要的包
library(readxl)
library(mmrm)
library(ggplot2)

?mmrm
# 读取数据
admadrs <- read_excel("/Users/qiujiahui/Desktop/≥55 older adults.xlsx", 
                      sheet = "Sheet1", col_names = TRUE, 
                      col_types = NULL, na = "", skip = 0)
head(admadrs)

# 删除重复的时间点
admadrs <- admadrs[!duplicated(admadrs[c("USUBJID", "AVISIT")]), ]

# 将相关列转换为因子
admadrs$USUBJID <- factor(admadrs$USUBJID)
admadrs$AVISITN <- factor(admadrs$AVISITN)
admadrs$ADY <- as.factor(admadrs$ADY)
admadrs$SITEID <- as.factor(admadrs$SITEID)
admadrs$TRTSEQA <- as.factor(admadrs$TRTSEQA)
admadrs$BASE <- as.numeric(admadrs$BASE)

# 将 TRT01A 转换为因子，并重新编码参考水平
admadrs$TRT01A <- as.factor(admadrs$TRT01A)
admadrs$TRT01A <- relevel(admadrs$TRT01A, ref = "Dummy A")

#admadrs$AVISITN <- relevel(admadrs$AVISITN, ref = "")

# 对数值型变量进行归一化处理


str(admadrs$BASE)

# 拟合模型
fit <- mmrm(
  formula = AVAL ~ AVISITN:TRT01A + TRT01A+ SEX + RACE + SITEID  + us(AVISITN | USUBJID),
  data = admadrs
)

# 查看模型结果
model_summary <- summary(fit)
print(model_summary)

≥55 older adults 分析

############## 老年成人分析 ##############
# 读取老年成人数据并拟合模型
# (假设这里已经完成)

cat("\n### 老年成人组分析 ###\n")

# 重新计算边际均值
emm_older <- emmeans(fit, ~ AVISITN | TRT01A)
print(emm_older)

# 执行对比分析
tryCatch({
  # 对每个治疗组，计算各时间点与基线的对比
  comparisons_older <- contrast(emm_older, method = "revpairwise", by = "TRT01A")
  cat("\n老年成人组的时间点比较:\n")
  print(comparisons_older)
  
  # 保存结果
  results_older <- as.data.frame(comparisons_older)
  results_older$AgeGroup <- "Older Adults (55+)"
  
}, error = function(e) {
  cat("老年成人组的对比分析出错:", conditionMessage(e), "\n")
  results_older <- data.frame()
})

年龄组疗效

############## 比较不同年龄组疗效差异 ##############
cat("\n### 比较不同年龄组疗效差异 ###\n")

# 在合并数据之前确保AgeGroup标签正确
results_young$AgeGroup <- "Young Adults (18-24)"
results_middle$AgeGroup <- "Middle-Aged Adults (24-55)"
results_older$AgeGroup <- "Older Adults (55+)"

# 合并数据
combined_results <- rbind(results_young, results_middle, results_older)

# 提取对比信息，并简化时间点标签
combined_results$TimeComparison <- gsub(" - .*$", "", combined_results$contrast)
# 使用更简洁的时间点标签，例如T1, T2, T3...
combined_results$TimePoint <- factor(combined_results$TimeComparison, 
                                    levels = unique(combined_results$TimeComparison),
                                    labels = paste0("T", seq_along(unique(combined_results$TimeComparison))))

# 创建改进的交互图
p <- ggplot(combined_results, aes(x = TimePoint, y = estimate, color = TRT01A, 
                               group = interaction(TRT01A, AgeGroup), 
                               linetype = AgeGroup)) +
  geom_line() +
  geom_point() +
  geom_errorbar(aes(ymin = estimate - SE, ymax = estimate + SE), width = 0.2) +
  facet_wrap(~ TRT01A) +
  labs(title = "Comparison of Treatment Effects Across Age Groups",
       subtitle = "Changes compared to baseline",
       x = "Time Point",
       y = "Estimated Difference",
       color = "Treatment Group",
       linetype = "Age Group") +
  theme_minimal() +
  # 改进x轴标签显示
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, size = 8),  # 倾斜标签并减小字体
    legend.position = "bottom",  # 将图例放在底部
    legend.box = "vertical"  # 垂直排列图例
  )

# 打印图表
print(p)

# 添加时间点说明表格
time_point_table <- data.frame(
  TimePointLabel = paste0("T", seq_along(unique(combined_results$TimeComparison))),
  OriginalTimePoint = unique(combined_results$TimeComparison)
)

# 打印时间点对照表
cat("\nTime Point Reference Table:\n")
print(time_point_table)

# 创建一个函数提取特定对比和治疗组的结果
get_result <- function(results_df, contrast_pattern, trt) {
  subset(results_df, grepl(contrast_pattern, contrast) & TRT01A == trt)
}

# 定义要比较的对比模式（如0 vs 2, 0 vs 8等）
contrast_patterns <- paste0(baseline_level, " - ", time_points)

# 创建数据框用于年龄组间比较
# 年轻成人 vs 中年成人
young_vs_middle_table <- data.frame()
# 年轻成人 vs 老年成人
young_vs_older_table <- data.frame()
# 中年成人 vs 老年成人
middle_vs_older_table <- data.frame()

for(pattern in contrast_patterns) {
  for(trt in unique(combined_results$TRT01A)) {
    # 获取三个年龄组的结果
    young_result <- get_result(results_young, pattern, trt)
    middle_result <- get_result(results_middle, pattern, trt)
    older_result <- get_result(results_older, pattern, trt)
    
    # 年轻vs中年
    if(nrow(young_result) > 0 && nrow(middle_result) > 0) {
      diff_estimate <- young_result$estimate - middle_result$estimate
      diff_se <- sqrt(young_result$SE^2 + middle_result$SE^2)
      z_stat <- diff_estimate / diff_se
      p_value <- 2 * (1 - pnorm(abs(z_stat)))
      
      young_vs_middle_table <- rbind(young_vs_middle_table, data.frame(
        Contrast = pattern,
        TRT01A = trt,
        Young_Est = young_result$estimate,
        Young_SE = young_result$SE,
        Middle_Est = middle_result$estimate,
        Middle_SE = middle_result$SE,
        Diff_Est = diff_estimate,
        Diff_SE = diff_se,
        Z_stat = z_stat,
        P_value = p_value
      ))
    }
    
    # 年轻vs老年
    if(nrow(young_result) > 0 && nrow(older_result) > 0) {
      diff_estimate <- young_result$estimate - older_result$estimate
      diff_se <- sqrt(young_result$SE^2 + older_result$SE^2)
      z_stat <- diff_estimate / diff_se
      p_value <- 2 * (1 - pnorm(abs(z_stat)))
      
      young_vs_older_table <- rbind(young_vs_older_table, data.frame(
        Contrast = pattern,
        TRT01A = trt,
        Young_Est = young_result$estimate,
        Young_SE = young_result$SE,
        Older_Est = older_result$estimate,
        Older_SE = older_result$SE,
        Diff_Est = diff_estimate,
        Diff_SE = diff_se,
        Z_stat = z_stat,
        P_value = p_value
      ))
    }
    
    # 中年vs老年
    if(nrow(middle_result) > 0 && nrow(older_result) > 0) {
      diff_estimate <- middle_result$estimate - older_result$estimate
      diff_se <- sqrt(middle_result$SE^2 + older_result$SE^2)
      z_stat <- diff_estimate / diff_se
      p_value <- 2 * (1 - pnorm(abs(z_stat)))
      
      middle_vs_older_table <- rbind(middle_vs_older_table, data.frame(
        Contrast = pattern,
        TRT01A = trt,
        Middle_Est = middle_result$estimate,
        Middle_SE = middle_result$SE,
        Older_Est = older_result$estimate,
        Older_SE = older_result$SE,
        Diff_Est = diff_estimate,
        Diff_SE = diff_se,
        Z_stat = z_stat,
        P_value = p_value
      ))
    }
  }
}

# 输出比较结果
cat("\n年轻成人 vs. 中年成人组疗效比较:\n")
if(nrow(young_vs_middle_table) > 0) {
  print(young_vs_middle_table)
} else {
  cat("没有足够数据进行比较\n")
}

cat("\n年轻成人 vs. 老年成人组疗效比较:\n")
if(nrow(young_vs_older_table) > 0) {
  print(young_vs_older_table)
} else {
  cat("没有足够数据进行比较\n")
}

cat("\n中年成人 vs. 老年成人组疗效比较:\n")
if(nrow(middle_vs_older_table) > 0) {
  print(middle_vs_older_table)
} else {
  cat("没有足够数据进行比较\n")
}

# 添加显著性标记到图表
# 创建一个数据框，包含所有显著差异的p值<0.05的结果
significant_diffs <- rbind(
  subset(young_vs_middle_table, P_value < 0.05),
  subset(young_vs_older_table, P_value < 0.05),
  subset(middle_vs_older_table, P_value < 0.05)
)

if(nrow(significant_diffs) > 0) {
  cat("\n显著性差异总结 (p < 0.05):\n")
  print(significant_diffs)
}
# 创建年龄组显著性热图
if(nrow(young_vs_middle_table) > 0 && nrow(young_vs_older_table) > 0 && nrow(middle_vs_older_table) > 0) {
  # 准备热图数据
  heatmap_data <- NULL
  
  # 处理各组间比较的p值
  for(trt in unique(combined_results$TRT01A)) {
    for(pattern in contrast_patterns) {
      # 提取各组比较的p值
      y_m_row <- subset(young_vs_middle_table, Contrast == pattern & TRT01A == trt)
      y_o_row <- subset(young_vs_older_table, Contrast == pattern & TRT01A == trt)
      m_o_row <- subset(middle_vs_older_table, Contrast == pattern & TRT01A == trt)
      
      if(nrow(y_m_row) > 0 && nrow(y_o_row) > 0 && nrow(m_o_row) > 0) {
        new_row <- data.frame(
          TimePoint = pattern,
          Treatment = trt,
          Young_vs_Middle_p = y_m_row$P_value,
          Young_vs_Older_p = y_o_row$P_value,
          Middle_vs_Older_p = m_o_row$P_value
        )
        heatmap_data <- rbind(heatmap_data, new_row)
      }
    }
  }
  
  if(!is.null(heatmap_data)) {
    # 创建长格式数据用于热图
    heatmap_long <- tidyr::pivot_longer(
      heatmap_data, 
      cols = c(Young_vs_Middle_p, Young_vs_Older_p, Middle_vs_Older_p),
      names_to = "Comparison", 
      values_to = "P_Value"
    )
    
    # 生成热图
    p_heatmap <- ggplot(heatmap_long, aes(x = TimePoint, y = Comparison, fill = -log10(P_Value))) +
      geom_tile() +
      scale_fill_gradient(low = "white", high = "red") +
      facet_wrap(~ Treatment) +
      labs(
        title = "Statistical Significance of Age Group Differences",
        subtitle = "-log10(p-value): higher values = more significant",
        x = "Time Point", 
        y = "Group Comparison",
        fill = "-log10(p)"
      ) +
      theme_minimal() +
      theme(axis.text.x = element_text(angle = 45, hjust = 1))
    
    print(p_heatmap)
  }
}

# 创建各年龄组在各治疗组的治疗效应折线图
# 将结果按年龄组和治疗方案分组
age_group_effects <- ggplot(combined_results, 
                           aes(x = TimePoint, y = estimate, 
                               group = AgeGroup, color = AgeGroup)) +
  geom_line() +
  geom_point() +
  geom_errorbar(aes(ymin = estimate - SE, ymax = estimate + SE), width = 0.2) +
  facet_wrap(~ TRT01A) +
  labs(title = "Treatment Effects by Age Group",
       x = "Time Point", 
       y = "Estimated Effect",
       color = "Age Group") +
  theme_minimal() +
  theme(legend.position = "bottom")

print(age_group_effects)

# 总结分析
cat("\n## 年龄组治疗效应总结分析\n")
cat("1. 基于以上分析，我们可以观察到不同年龄组对治疗的反应存在以下差异：\n")

# 检查是否有显著差异
if(exists("significant_diffs") && nrow(significant_diffs) > 0) {
  cat("   - 存在显著的年龄组间差异(p < 0.05)，特别是在以下时间点和治疗组：\n")
  for(i in 1:nrow(significant_diffs)) {
    row <- significant_diffs[i,]
    cat(paste0("     * 时间点 ", row$Contrast, "，治疗组 ", row$TRT01A, 
               " (p = ", round(row$P_value, 4), ")\n"))
  }
} else {
  cat("   - 在分析的时间点和治疗组中未发现显著的年龄组间差异\n")
}

# 药效学差异总结
cat("\n2. 药效学差异总结：\n")
cat("   - young adults (18-24 years)：\n")
young_effect <- mean(subset(results_young, TRT01A != "Dummy A")$estimate)
cat(paste0("     Mean therapeutic effect: ", round(young_effect, 2), "\n"))

cat("   - middle-aged adults (24-55 years)：\n")
middle_effect <- mean(subset(results_middle, TRT01A != "Dummy A")$estimate)
cat(paste0("     Mean therapeutic effect: ", round(middle_effect, 2), "\n"))

cat("   - older adults (>= 55 years)：\n")
older_effect <- mean(subset(results_older, TRT01A != "Dummy A")$estimate)
cat(paste0("     Mean therapeutic effect: ", round(older_effect, 2), "\n"))

# 整体结论
cat("\n3. 临床意义和建议：\n")
# 比较三组平均效应
groups <- c("年轻成人", "中年成人", "老年成人")
effects <- c(young_effect, middle_effect, older_effect)
best_group_idx <- which.max(abs(effects))
lowest_group_idx <- which.min(abs(effects))

cat(paste0("   - ", groups[best_group_idx], "组显示出最强的治疗反应\n"))
cat(paste0("   - ", groups[lowest_group_idx], "组显示出最弱的治疗反应\n"))

# 最终建议
cat("\n4. 临床应用建议：\n")
cat("   - 根据不同年龄组的反应差异，可能需要考虑针对特定年龄组调整治疗方案\n")
cat("   - 特别关注对治疗反应较弱的年龄组，可能需要增加治疗强度或辅助治疗\n")
cat("   - 应进一步研究年龄相关的药代动力学和药效学差异，以优化治疗策略\n")

主要疗效分析（Day28）

library(emmeans)
mmrm() registered as emmeans extension
Welcome to emmeans.
Caution: You lose important information if you filter this package's results.
See '? untidy'
library(readxl)
library(dplyr)
library(ggplot2)
library(emmeans)
library(mmrm)
library(tidyr)
library(viridis)
Loading required package: viridisLite
library(grid)
library(gridExtra)

Attaching package: 'gridExtra'
The following object is masked from 'package:dplyr':

    combine
# 获取Day28的最小二乘均值
emm_day28 <- emmeans(
  fit, 
  specs = ~ TRT01A | AVISITN,
  at = list(AVISITN = "28") # "28"是AVISITN中Day28对应的水平（减20000后）
)

# 治疗组与安慰剂对比
contrast_specs <- list(
  "Esk56_vs_Placebo" = c(-1, 1, 0),
  "Esk84_vs_Placebo" = c(-1, 0, 1)
)

contrast_results <- contrast(
  emm_day28,
  method = contrast_specs,
  adjust = "none"
) %>% 
  summary(infer = TRUE) # 包含置信区间和p值

# 显示对比结果
print(contrast_results)
AVISITN = 28:
 contrast         estimate     SE  df lower.CL upper.CL t.ratio p.value
 Esk56_vs_Placebo   -0.841 0.0992 651   -1.036   -0.647  -8.482  <.0001
 Esk84_vs_Placebo   -0.798 0.0791 596   -0.954   -0.643 -10.098  <.0001

Results are averaged over the levels of: SEX, SITEID, RACE, base_situation 
Confidence level used: 0.95

疗效趋势可视化及森林图

library(ggplot2)
# 获取各时间点预测值
#admadrs$AVISITN <- droplevels(admadrs$AVISITN)
emm_trend <- emmeans(
  fit,
  specs = ~ TRT01A | AVISITN, 
  at = list(AVISITN = levels(admadrs$AVISITN)),# 使用所有时间点
  nuisance = "RACE") %>% # 将RACE设为无关因素保证行数不超过10000
  summary()

# 绘制疗效趋势图
(p_trend <- ggplot(emm_trend, aes(x = AVISITN, y = emmean, 
                                 color = TRT01A, group = TRT01A)) +
  geom_line(linewidth = 1.2) +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), 
                width = 0.2, 
                position = position_dodge(0.1)) +
  geom_point(size = 3) +
  scale_color_viridis(discrete = TRUE, 
                     labels = c("Placebo", "56mg", "84mg")) +
  labs(
    title = "Treatment Effect Over Time",
    x = "Study Visit (Day)", 
    y = "Adjusted Mean Change from Baseline (LS Means)", 
    color = "Treatment Group"
  ) +
  theme_minimal() +
  theme(
    legend.position = "top",
    plot.title = element_text(hjust = 0.5, face = "bold"),
    axis.text.x = element_text(angle = 45, hjust = 1)
  ))

# 绘制森林图
contrast_df <- as.data.frame(contrast_results)
contrast_df$Comparison <- rownames(contrast_df)

(p_forest <- ggplot(contrast_df, 
                   aes(x = estimate, y = Comparison, 
                       xmin = lower.CL, xmax = upper.CL)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey50") +
  geom_errorbarh(height = 0.2) +
  geom_point(size = 3, color = "#440154") +
  scale_y_discrete(labels = c("56mg vs Placebo", "84mg vs Placebo")) +
  labs(
    title = "Treatment Comparisons at Day 28",
    x = "Treatment Difference (95% CI)", 
    y = ""
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    panel.grid.major.y = element_blank()
  ))

# 组合图形
grid.arrange(p_trend, p_forest, ncol = 1, heights = c(2, 1))

# 保存图形（PDF格式）
ggsave("treatment_analysis.pdf", 
       arrangeGrob(p_trend, p_forest, ncol = 1), 
       width = 10, height = 8)

# 保存统计结果
write.csv(contrast_results, "treatment_contrasts.csv", row.names = FALSE)

治疗组vs安慰剂的估计边际均值和差异计算及可视化

# 将Day 28对应的AVISITN水平赋值
day28_level <- "28" 

# 计算Day 28时的估计边际均值
em_day28 <- emmeans(fit, ~ TRT01A | AVISITN, at = list(AVISITN = day28_level))
print(em_day28)
AVISITN = 28:
 TRT01A  emmean    SE   df lower.CL upper.CL
 Dummy A -0.266 0.129  798   -0.519  -0.0126
 Dummy B -1.107 0.157 1000   -1.416  -0.7989
 Dummy C -1.064 0.147  979   -1.352  -0.7764

Results are averaged over the levels of: SEX, SITEID, RACE, base_situation 
Confidence level used: 0.95 
# 计算治疗组vs安慰剂在Day28时的差异
contrasts_day28 <- contrast(em_day28, method = "trt.vs.ctrl", ref = 1)
print(contrasts_day28)
AVISITN = 28:
 contrast          estimate     SE  df t.ratio p.value
 Dummy B - Dummy A   -0.841 0.0992 651  -8.482  <.0001
 Dummy C - Dummy A   -0.798 0.0791 596 -10.098  <.0001

Results are averaged over the levels of: SEX, SITEID, RACE, base_situation 
P value adjustment: dunnettx method for 2 tests 
# 计算所有时间点的估计边际均值
em_all <- emmeans(fit, ~ TRT01A | AVISITN, nuisance = "RACE")
print(em_all)
AVISITN = 0:
 TRT01A  emmean    SE   df lower.CL upper.CL
 Dummy A  0.304 0.103  416   0.1020   0.5053
 Dummy B  0.305 0.105  417   0.0994   0.5113
 Dummy C  0.353 0.106  417   0.1445   0.5616

AVISITN = 2:
 TRT01A  emmean    SE   df lower.CL upper.CL
 Dummy A -0.703 0.117  634  -0.9325  -0.4728
 Dummy B -1.396 0.149  883  -1.6881  -1.1032
 Dummy C -1.237 0.147  992  -1.5249  -0.9488

AVISITN = 8:
 TRT01A  emmean    SE   df lower.CL upper.CL
 Dummy A -0.631 0.132  823  -0.8899  -0.3720
 Dummy B -1.327 0.155 1028  -1.6310  -1.0222
 Dummy C -1.404 0.150 1035  -1.6972  -1.1103

AVISITN = 15:
 TRT01A  emmean    SE   df lower.CL upper.CL
 Dummy A -0.677 0.132  829  -0.9371  -0.4176
 Dummy B -1.419 0.160 1049  -1.7321  -1.1051
 Dummy C -1.467 0.184  763  -1.8276  -1.1055

AVISITN = 22:
 TRT01A  emmean    SE   df lower.CL upper.CL
 Dummy A -0.347 0.116  614  -0.5746  -0.1202
 Dummy B -1.188 0.152  884  -1.4870  -0.8897
 Dummy C -0.993 0.150  982  -1.2870  -0.6987

AVISITN = 28:
 TRT01A  emmean    SE   df lower.CL upper.CL
 Dummy A -0.266 0.129  798  -0.5193  -0.0126
 Dummy B -1.107 0.157 1000  -1.4157  -0.7989
 Dummy C -1.064 0.147  979  -1.3520  -0.7764

AVISITN = 999:
 TRT01A  emmean    SE   df lower.CL upper.CL
 Dummy A -0.413 0.130  805  -0.6673  -0.1587
 Dummy B -1.114 0.157  980  -1.4215  -0.8068
 Dummy C -1.214 0.184  758  -1.5761  -0.8519

Results are averaged over the levels of: 1 nuisance factors, SEX, SITEID, base_situation 
Confidence level used: 0.95 
# 各组在各时间点的均值比较
contrasts_all <- contrast(em_all, method = "trt.vs.ctrl", ref = 1)
print(contrasts_all)
AVISITN = 0:
 contrast          estimate     SE  df t.ratio p.value
 Dummy B - Dummy A  0.00171 0.0411 415   0.042  0.9972
 Dummy C - Dummy A  0.04943 0.0418 416   1.183  0.3944

AVISITN = 2:
 contrast          estimate     SE  df t.ratio p.value
 Dummy B - Dummy A -0.69297 0.0748 588  -9.270  <.0001
 Dummy C - Dummy A -0.53418 0.0923 670  -5.788  <.0001

AVISITN = 8:
 contrast          estimate     SE  df t.ratio p.value
 Dummy B - Dummy A -0.69563 0.0926 668  -7.511  <.0001
 Dummy C - Dummy A -0.77280 0.0796 596  -9.714  <.0001

AVISITN = 15:
 contrast          estimate     SE  df t.ratio p.value
 Dummy B - Dummy A -0.74127 0.0962 509  -7.705  <.0001
 Dummy C - Dummy A -0.78919 0.0970 516  -8.134  <.0001

AVISITN = 22:
 contrast          estimate     SE  df t.ratio p.value
 Dummy B - Dummy A -0.84093 0.0788 590 -10.675  <.0001
 Dummy C - Dummy A -0.64546 0.0991 657  -6.515  <.0001

AVISITN = 28:
 contrast          estimate     SE  df t.ratio p.value
 Dummy B - Dummy A -0.84138 0.0992 651  -8.482  <.0001
 Dummy C - Dummy A -0.79832 0.0791 596 -10.098  <.0001

AVISITN = 999:
 contrast          estimate     SE  df t.ratio p.value
 Dummy B - Dummy A -0.70119 0.0947 512  -7.401  <.0001
 Dummy C - Dummy A -0.80107 0.0953 515  -8.406  <.0001

Results are averaged over the levels of: 1 nuisance factors, SEX, SITEID, base_situation 
P value adjustment: dunnettx method for 2 tests 
# 将结果转换为数据框以便绘图
em_all_df <- as.data.frame(em_all)

# 创建一个时间轴变量，将AVISITN转换为数值以便绘图
em_all_df$Time <- as.numeric(as.character(em_all_df$AVISITN))

# 绘制各治疗组随时间的MADRS变化趋势图
ggplot(em_all_df, aes(x = AVISITN, y = emmean, color = TRT01A, group = TRT01A)) +
  geom_line(linewidth = 1) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = emmean - SE, ymax = emmean + SE), width = 0.2) +
  labs(
    title = "Change from Baseline in MADRS Total Score Over Time",
    x = "Day",
    y = "Estimated Mean Change from Baseline",
    color = "Treatment Group"
  ) +
  scale_color_viridis(discrete = TRUE, option = "D", 
                      labels = c("Placebo", "Esketamine 56 mg", "Esketamine 84 mg")) +
  theme_minimal() +
  theme(
    legend.position = "bottom",
    panel.grid.minor = element_blank(),
    axis.title = element_text(face = "bold"),
    plot.title = element_text(hjust = 0.5, face = "bold")
  )

# 保存图片
ggsave("MADRS_change_over_time.png", width = 10, height = 6, dpi = 300)

# 创建Day 28时的柱状图比较
day28_df <- subset(em_all_df, AVISITN == day28_level)

# 将各组在各时间点的均值比较结果转换为数据框以便绘图
contrasts_day28_df <- as.data.frame(contrasts_day28)

# 绘制Day 28时的柱状图
p1 <- ggplot(day28_df, aes(x = TRT01A, y = emmean, fill = TRT01A)) +
  geom_bar(stat = "identity", width = 0.7) +
  geom_errorbar(aes(ymin = emmean - SE, ymax = emmean + SE), width = 0.2) +
  labs(
    title = "Change from Baseline in MADRS Total Score at Day 28",
    x = "Treatment Group",
    y = "Estimated Mean Change from Baseline",
    fill = "Treatment Group"
  ) +
  scale_fill_viridis(discrete = TRUE, option = "D",
                    labels = c("Placebo", "Esketamine 56 mg", "Esketamine 84 mg")) +
  theme_minimal() +
  theme(
    legend.position = "none",
    axis.title = element_text(face = "bold"),
    plot.title = element_text(hjust = 0.5, face = "bold")
  )

# 绘制治疗组vs安慰剂的差异图
p2 <- ggplot(contrasts_day28_df, aes(x = contrast, y = estimate, fill = contrast)) +
  geom_bar(stat = "identity", width = 0.7) +
  geom_errorbar(aes(ymin = estimate - SE, ymax = estimate + SE), width = 0.2) +
  labs(
    title = "Treatment Effect vs Placebo at Day 28",
    x = "Contrast",
    y = "Estimated Difference",
    fill = "Contrast"
  ) +
  scale_fill_viridis(discrete = TRUE, option = "E") +
  theme_minimal() +
  theme(
    legend.position = "none",
    axis.title = element_text(face = "bold"),
    plot.title = element_text(hjust = 0.5, face = "bold")
  ) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red")

# 将两个图组合在一起
combined_plot <- grid.arrange(p1, p2, ncol = 2)

# 保存组合图
ggsave("Day28_treatment_effect.png", combined_plot, width = 12, height = 6, dpi = 300)

# 创建森林图展示各时间点的治疗效果
contrasts_all_df <- as.data.frame(contrasts_all)
contrasts_all_df$Time <- as.numeric(as.character(contrasts_all_df$AVISITN))

# 森林图
ggplot(contrasts_all_df, aes(x = AVISITN , y = estimate, color = contrast)) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = estimate - SE*1.96, ymax = estimate + SE*1.96), width = 0.2) +
  geom_line(aes(group = contrast), linewidth = 1) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  facet_wrap(~ contrast, ncol = 1) +
  labs(
    title = "Treatment Effect vs Placebo Over Time",
    x = "Day",
    y = "Estimated Difference vs Placebo",
    color = "Treatment Group"
  ) +
  scale_color_viridis(discrete = TRUE, option = "D") +
  theme_minimal() +
  theme(
    legend.position = "none",
    strip.background = element_rect(fill = "lightgrey"),
    strip.text = element_text(face = "bold"),
    axis.title = element_text(face = "bold"),
    plot.title = element_text(hjust = 0.5, face = "bold")
  )

# 保存森林图
ggsave("Treatment_effect_forest_plot.png", width = 10, height = 8, dpi = 300)

# 创建热图展示p值
# 提取所有时间点的p值
p_values <- contrasts_all_df
p_values$sig <- ifelse(p_values$p.value < 0.001, "***",
                      ifelse(p_values$p.value < 0.01, "**",
                            ifelse(p_values$p.value < 0.05, "*", "ns")))

# 处理热图所需数据
p_values_wide <- p_values %>%
  select(AVISITN, contrast, p.value) %>%
  pivot_wider(names_from = contrast, values_from = p.value)

# 绘制p值热图
p_values_long <- p_values %>%
  select(Time, contrast, p.value, sig)

ggplot(p_values_long, aes(x = factor(Time), y = contrast, fill = -log10(p.value))) +
  geom_tile(color = "white") +
  geom_text(aes(label = sig), color = "white", size = 5) +
  scale_fill_viridis(option = "C", name = "-log10(p-value)") +
  labs(
    title = "Statistical Significance of Treatment Effect vs Placebo",
    x = "Day",
    y = "Treatment Group",
    caption = "* p<0.05, ** p<0.01, *** p<0.001, ns: not significant"
  ) +
  theme_minimal() +
  theme(
    axis.title = element_text(face = "bold"),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.position = "right"
  )

# 保存热图
ggsave("Treatment_significance_heatmap.png", width = 10, height = 6, dpi = 300)

# 模型诊断图
# 提取残差
residuals <- residuals(fit)
fitted_values <- fitted(fit)

# 创建诊断数据框
diag_df <- data.frame(
  Fitted = fitted_values,
  Residuals = residuals
)

# 绘制残差图
p3 <- ggplot(diag_df, aes(x = Fitted, y = Residuals)) +
  geom_point(alpha = 0.5) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  geom_smooth(method = "loess", se = FALSE, color = "blue") +
  labs(
    title = "Residuals vs Fitted Values",
    x = "Fitted Values",
    y = "Residuals"
  ) +
  theme_minimal() +
  theme(
    axis.title = element_text(face = "bold"),
    plot.title = element_text(hjust = 0.5, face = "bold")
  )

# 绘制残差的QQ图
p4 <- ggplot(diag_df, aes(sample = Residuals)) +
  stat_qq() +
  stat_qq_line() +
  labs(
    title = "Normal Q-Q Plot of Residuals",
    x = "Theoretical Quantiles",
    y = "Sample Quantiles"
  ) +
  theme_minimal() +
  theme(
    axis.title = element_text(face = "bold"),
    plot.title = element_text(hjust = 0.5, face = "bold")
  )

# 组合诊断图
diag_plot <- grid.arrange(p3, p4, ncol = 2)
`geom_smooth()` using formula = 'y ~ x'

# 保存诊断图
ggsave("Model_diagnostics.png", diag_plot, width = 12, height = 6, dpi = 300)

数据可视化（有图例和解释）

# 各治疗组随时间的MADRS变化趋势图

# 输出emmeans结果数据
em_all_df
AVISITN = 0:
 TRT01A      emmean        SE      df   lower.CL   upper.CL Time
 Dummy A  0.3036347 0.1025859  415.71  0.1019829  0.5052865    0
 Dummy B  0.3053490 0.1047895  416.57  0.0993668  0.5113311    0
 Dummy C  0.3530675 0.1061115  416.77  0.1444871  0.5616480    0

AVISITN = 2:
 TRT01A      emmean        SE      df   lower.CL   upper.CL Time
 Dummy A -0.7026543 0.1170686  634.38 -0.9325430 -0.4727655    2
 Dummy B -1.3956254 0.1489953  883.30 -1.6880515 -1.1031993    2
 Dummy C -1.2368376 0.1467976  992.47 -1.5249070 -0.9487683    2

AVISITN = 8:
 TRT01A      emmean        SE      df   lower.CL   upper.CL Time
 Dummy A -0.6309696 0.1319135  823.38 -0.8898959 -0.3720433    8
 Dummy B -1.3266031 0.1551385 1028.50 -1.6310271 -1.0221791    8
 Dummy C -1.4037658 0.1495618 1035.19 -1.6972447 -1.1102869    8

AVISITN = 15:
 TRT01A      emmean        SE      df   lower.CL   upper.CL Time
 Dummy A -0.6773596 0.1323395  828.55 -0.9371198 -0.4175994   15
 Dummy B -1.4186303 0.1597666 1049.01 -1.7321288 -1.1051319   15
 Dummy C -1.4665458 0.1839236  762.97 -1.8276022 -1.1054894   15

AVISITN = 22:
 TRT01A      emmean        SE      df   lower.CL   upper.CL Time
 Dummy A -0.3473856 0.1156954  613.79 -0.5745923 -0.1201788   22
 Dummy B -1.1883160 0.1521635  884.23 -1.4869597 -0.8896723   22
 Dummy C -0.9928504 0.1499149  981.72 -1.2870410 -0.6986599   22

AVISITN = 28:
 TRT01A      emmean        SE      df   lower.CL   upper.CL Time
 Dummy A -0.2659051 0.1290646  798.18 -0.5192512 -0.0125591   28
 Dummy B -1.1072886 0.1571637  999.66 -1.4156972 -0.7988800   28
 Dummy C -1.0642296 0.1466599  979.15 -1.3520335 -0.7764257   28

AVISITN = 999:
 TRT01A      emmean        SE      df   lower.CL   upper.CL Time
 Dummy A -0.4129694 0.1295510  804.86 -0.6672671 -0.1586717  999
 Dummy B -1.1141616 0.1566069  980.15 -1.4214850 -0.8068383  999
 Dummy C -1.2140361 0.1844488  757.54 -1.5761277 -0.8519445  999

Results are averaged over the levels of: 1 nuisance factors, SEX, SITEID, base_situation 
Confidence level used: 0.95 
# 绘图
ggplot(em_all_df, aes(x = AVISITN, y = emmean, color = TRT01A, group = TRT01A)) +
  geom_line(linewidth = 1) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = emmean - SE, ymax = emmean + SE), width = 0.2) +
  labs(
    title = "Change from Baseline in MADRS Total Score Over Time",
    subtitle = "Mixed Model for Repeated Measures (MMRM) Analysis",
    x = "Study Day",
    y = "Estimated Mean Change from Baseline",
    color = "Treatment Group",
    caption = "Note: Error bars represent standard errors. Negative values indicate improvement in depressive symptoms."
  ) +
  scale_color_viridis(discrete = TRUE, option = "D", 
                      labels = c("Placebo", "Esketamine 56 mg", "Esketamine 84 mg")) +
  theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.background = element_rect(fill = "white", color = "grey80"),
    legend.title = element_text(face = "bold"),
    panel.grid.minor = element_blank(),
    axis.title = element_text(face = "bold"),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5),
    plot.caption = element_text(hjust = 0, face = "italic")
  ) +
  # 添加解释文本
  annotate("text", x = min(as.numeric(as.character(em_all_df$AVISITN))) + 2, 
           y = min(em_all_df$emmean) * 0.8, 
           label = "Lower values indicate greater\nimprovement in depression", 
           hjust = 0, size = 3.5, fontface = "italic")

# 保存图片
ggsave("MADRS_change_over_time_with_legend.png", width = 10, height = 6, dpi = 300)


# 输出Day 28时的估计边际均值
day28_df
    TRT01A AVISITN     emmean        SE       df   lower.CL    upper.CL Time
16 Dummy A      28 -0.2659051 0.1290646 798.1815 -0.5192512 -0.01255906   28
17 Dummy B      28 -1.1072886 0.1571637 999.6623 -1.4156972 -0.79888004   28
18 Dummy C      28 -1.0642296 0.1466599 979.1475 -1.3520335 -0.77642572   28
# Day 28时的柱状图
p1 <- ggplot(day28_df, aes(x = TRT01A, y = emmean, fill = TRT01A)) +
  geom_bar(stat = "identity", width = 0.7, color = "black") +
  geom_errorbar(aes(ymin = emmean - SE, ymax = emmean + SE), width = 0.2) +
  labs(
    title = "Change from Baseline in MADRS Total Score at Day 28",
    subtitle = "Estimated Marginal Means with Standard Errors",
    x = "Treatment Group",
    y = "Estimated Mean Change from Baseline",
    fill = "Treatment Group",
    caption = "Note: Error bars represent standard errors. Negative values indicate improvement in depressive symptoms."
  ) +
  scale_fill_viridis(discrete = TRUE, option = "D",
                    labels = c("Placebo", "Esketamine 56 mg", "Esketamine 84 mg")) +
  theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.background = element_rect(fill = "white", color = "grey80"),
    legend.title = element_text(face = "bold"),
    axis.title = element_text(face = "bold"),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5),
    plot.caption = element_text(hjust = 0, face = "italic"),
    axis.text.x = element_text(angle = 0, hjust = 0.5)
  ) +
  # 添加标签
  geom_text(aes(label = sprintf("%.1f", emmean)), position = position_dodge(width = 0.7), 
            vjust = -0.5, color = "black", fontface = "bold")


# 输出绘图数据
contrasts_day28_df
AVISITN = 28:
 contrast            estimate         SE     df t.ratio p.value
 Dummy B - Dummy A -0.8413835 0.09919806 650.97  -8.482  <.0001
 Dummy C - Dummy A -0.7983245 0.07905644 595.89 -10.098  <.0001

Results are averaged over the levels of: SEX, SITEID, RACE, base_situation 
P value adjustment: dunnettx method for 2 tests 
# 绘制治疗组vs安慰剂的差异图
p2 <- ggplot(contrasts_day28_df, aes(x = contrast, y = estimate, fill = contrast)) +
  geom_bar(stat = "identity", width = 0.7, color = "black") +
  geom_errorbar(aes(ymin = estimate - SE, ymax = estimate + SE), width = 0.2) +
  labs(
    title = "Treatment Effect vs Placebo at Day 28",
    subtitle = "Difference in Estimated Means with Standard Errors",
    x = "Treatment Comparison",
    y = "Estimated Difference vs Placebo",
    fill = "Treatment Comparison",
    caption = "Note: Values below zero favor active treatment over placebo.\nRed dashed line represents no difference from placebo."
  ) +
  scale_fill_viridis(discrete = TRUE, option = "E",
                    labels = c("Esketamine 56 mg vs Placebo", "Esketamine 84 mg vs Placebo")) +
  theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.background = element_rect(fill = "white", color = "grey80"),
    legend.title = element_text(face = "bold"),
    axis.title = element_text(face = "bold"),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5),
    plot.caption = element_text(hjust = 0, face = "italic"),
    axis.text.x = element_text(angle = 45, hjust = 1)
  ) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  # 添加p值和效应量标签
  geom_text(aes(label = sprintf("Diff: %.1f\np = %.3f", estimate, p.value)), 
            position = position_dodge(width = 0.7), vjust = -0.8, size = 3.5)

# 使用patchwork组合两个图
library(patchwork)
combined_plot <- p1 + p2 + plot_layout(ncol = 2)
combined_plot <- combined_plot + plot_annotation(
  title = "MADRS Score Change and Treatment Effect at Day 28",
  theme = theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 16))
)

# 保存组合图
ggsave("Day28_treatment_effect_with_legend.png", combined_plot, width = 12, height = 6, dpi = 300)

# 输出绘图数据
contrasts_all_df
AVISITN = 0:
 contrast            estimate         SE     df t.ratio p.value Time
 Dummy B - Dummy A  0.0017143 0.04113726 415.00   0.042  0.9972    0
 Dummy C - Dummy A  0.0494328 0.04177930 415.89   1.183  0.3944    0

AVISITN = 2:
 contrast            estimate         SE     df t.ratio p.value Time
 Dummy B - Dummy A -0.6929711 0.07475349 588.39  -9.270  <.0001    2
 Dummy C - Dummy A -0.5341834 0.09229578 670.43  -5.788  <.0001    2

AVISITN = 8:
 contrast            estimate         SE     df t.ratio p.value Time
 Dummy B - Dummy A -0.6956335 0.09261148 668.33  -7.511  <.0001    8
 Dummy C - Dummy A -0.7727962 0.07955698 595.88  -9.714  <.0001    8

AVISITN = 15:
 contrast            estimate         SE     df t.ratio p.value Time
 Dummy B - Dummy A -0.7412707 0.09620677 508.54  -7.705  <.0001   15
 Dummy C - Dummy A -0.7891862 0.09702519 515.94  -8.134  <.0001   15

AVISITN = 22:
 contrast            estimate         SE     df t.ratio p.value Time
 Dummy B - Dummy A -0.8409304 0.07877781 589.92 -10.675  <.0001   22
 Dummy C - Dummy A -0.6454649 0.09907534 656.61  -6.515  <.0001   22

AVISITN = 28:
 contrast            estimate         SE     df t.ratio p.value Time
 Dummy B - Dummy A -0.8413835 0.09919806 650.97  -8.482  <.0001   28
 Dummy C - Dummy A -0.7983245 0.07905644 595.89 -10.098  <.0001   28

AVISITN = 999:
 contrast            estimate         SE     df t.ratio p.value Time
 Dummy B - Dummy A -0.7011922 0.09474272 512.36  -7.401  <.0001  999
 Dummy C - Dummy A -0.8010667 0.09530024 514.97  -8.406  <.0001  999

Results are averaged over the levels of: 1 nuisance factors, SEX, SITEID, base_situation 
P value adjustment: dunnettx method for 2 tests 
# 创建森林图展示各时间点的治疗效果
ggplot(contrasts_all_df, aes(x = AVISITN, y = estimate, color = contrast)) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = estimate - SE*1.96, ymax = estimate + SE*1.96), width = 0.2) +
  geom_line(aes(group = contrast), linewidth = 1) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  facet_wrap(~ contrast, ncol = 1) +
  labs(
    title = "Treatment Effect vs Placebo Over Time",
    subtitle = "Estimated Differences with 95% Confidence Intervals",
    x = "Study Day",
    y = "Estimated Difference vs Placebo",
    color = "Treatment Comparison",
    caption = paste("Note: Points below zero favor active treatment over placebo.",
                   "Error bars represent 95% confidence intervals.",
                   "Non-overlapping confidence intervals with red dashed line (zero) indicate statistical significance.")
  ) +
  scale_color_viridis(discrete = TRUE, option = "D",
                     labels = c("Esketamine 56 mg vs Placebo", "Esketamine 84 mg vs Placebo")) +
  theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.background = element_rect(fill = "white", color = "grey80"),
    legend.title = element_text(face = "bold"),
    strip.background = element_rect(fill = "lightgrey"),
    strip.text = element_text(face = "bold"),
    axis.title = element_text(face = "bold"),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5),
    plot.caption = element_text(hjust = 0, face = "italic")
  ) +
  # 添加p值标签
  geom_text(aes(label = sprintf("p = %.3f", p.value)), 
            position = position_dodge(width = 0.5), vjust = -0.8, 
            color = "black", size = 3)

# 保存森林图
ggsave("Treatment_effect_forest_plot_with_legend.png", width = 10, height = 8, dpi = 300)




# 创建热图展示p值
# 提取所有时间点的p值
p_values <- contrasts_all_df
p_values$sig <- ifelse(p_values$p.value < 0.001, "***",
                     ifelse(p_values$p.value < 0.01, "**",
                          ifelse(p_values$p.value < 0.05, "*", "ns")))

# 数据准备
p_values_long <- p_values %>%
  select(Time, contrast, p.value, sig)

# 添加-log10(p值)列
p_values_long$neg_log10_p <- -log10(p_values_long$p.value)

p_values_long
   Time          contrast      p.value sig  neg_log10_p
1     0 Dummy B - Dummy A 9.971599e-01  ns  0.001235183
2     0 Dummy C - Dummy A 3.943676e-01  ns  0.404098792
3     2 Dummy B - Dummy A 2.484033e-10 ***  9.604842643
4     2 Dummy C - Dummy A 2.131253e-08 ***  7.671365020
5     8 Dummy B - Dummy A 1.882938e-13 *** 12.725163922
6     8 Dummy C - Dummy A 2.354221e-10 ***  9.628152722
7    15 Dummy B - Dummy A 8.864287e-11 *** 10.052356187
8    15 Dummy C - Dummy A 1.091197e-10 ***  9.962096784
9    22 Dummy B - Dummy A 2.467375e-10 ***  9.607764806
10   22 Dummy C - Dummy A 1.448419e-10 ***  9.839105738
11   28 Dummy B - Dummy A 2.220446e-16 *** 15.653559775
12   28 Dummy C - Dummy A 2.352738e-10 ***  9.628426432
13  999 Dummy B - Dummy A 1.000091e-10 ***  9.999960428
14  999 Dummy C - Dummy A 1.063928e-10 ***  9.973087829
# 绘制p值热图
ggplot(p_values_long, aes(x = factor(Time), y = contrast, fill =  -log10(p.value))) +
  geom_tile(color = "white") +
  # 根据显著性调整文本颜色
  geom_text(aes(label = sig, color = ifelse(p.value < 0.05, "white", "black")), size = 5) +
  # 手动设置文本颜色
  scale_color_identity() +
  # 设置小p值(高-log10(p值))为深色
  scale_fill_viridis(option = "C", name = "-log10(p-value)", direction = 1) +
  labs(
    title = "Statistical Significance of Treatment Effect vs Placebo",
    subtitle = "p-values Transformed to -log10 Scale for Visual Clarity",
    x = "Study Day",
    y = "Treatment Comparison",
    caption = paste("* p<0.05, ** p<0.01, *** p<0.001, ns: not significant",
                   "Darker colors represent smaller p-values (higher statistical significance)",
                   "The -log10 transformation converts p=0.05 to 1.3, p=0.01 to 2, p=0.001 to 3, etc.", sep = "\n")
  ) +
  theme_minimal() +
  theme(
    axis.title = element_text(face = "bold"),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5),
    plot.caption = element_text(hjust = 0, face = "italic"),
    legend.position = "right",
    legend.title = element_text(face = "bold"),
    legend.background = element_rect(fill = "white", color = "grey80")
  )

# 保存热图
ggsave("Treatment_significance_heatmap_with_legend.png", width = 10, height = 6, dpi = 300)



# 创建模型诊断图
# 提取残差
residuals <- residuals(fit)
fitted_values <- fitted(fit)

# 创建诊断数据框
diag_df <- data.frame(
  Fitted = fitted_values,
  Residuals = residuals
)


# 为QQ图准备数据
qq_data <- qqnorm(residuals, plot.it = FALSE)
qq_df <- data.frame(Theoretical = qq_data$x, Sample = qq_data$y)


# 绘制残差图
p3 <- ggplot(diag_df, aes(x = Fitted, y = Residuals)) +
  geom_point(alpha = 0.5, color = "blue") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  geom_smooth(method = "loess", se = TRUE, color = "red", fill = "pink", alpha = 0.3) +
  labs(
    title = "Residuals vs Fitted Values",
    subtitle = "Assessment of Homoscedasticity and Linearity",
    x = "Fitted Values",
    y = "Residuals",
    caption = "Note: Points should be randomly scattered around zero line (red dashed line).\nThe smooth curve (red) should be approximately horizontal."
  ) +
  theme_minimal() +
  theme(
    axis.title = element_text(face = "bold"),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5),
    plot.caption = element_text(hjust = 0, face = "italic")
  ) +
  # 添加解释性文本
  annotate("text", x = min(fitted_values) + 0.1*(max(fitted_values)-min(fitted_values)), 
           y = max(residuals)*0.8, 
           label = "Ideal pattern: Random scatter\naround the zero line with no pattern", 
           hjust = 0, size = 3.5, fontface = "italic")



# 绘制残差的QQ图
p4 <- ggplot(qq_df, aes(x = Theoretical, y = Sample)) +
  geom_point(alpha = 0.5, color = "blue") + 
  geom_abline(intercept = 0, slope = 10, color = "red") +
  labs(
    title = "Normal Q-Q Plot of Residuals",
    subtitle = "Assessment of Normality Assumption",
    x = "Theoretical Quantiles",
    y = "Sample Quantiles",
  ) +
  theme_minimal() +
  theme(
    axis.title = element_text(face = "bold"),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5),
    plot.caption = element_text(hjust = 0, face = "italic")
  ) 

# 使用patchwork组合诊断图
library(patchwork)
diag_plot <- p3 + p4 + plot_layout(ncol = 2)
diag_plot <- diag_plot + plot_annotation(
  title = "MMRM Model Diagnostic Plots",
  subtitle = "Assessing Model Assumptions for Valid Statistical Inference",
  caption = "These diagnostic plots help validate the assumptions of the Mixed Model for Repeated Measures (MMRM) analysis.",
  theme = theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 16),
               plot.subtitle = element_text(hjust = 0.5),
               plot.caption = element_text(hjust = 0, face = "italic"))
)

diag_plot
`geom_smooth()` using formula = 'y ~ x'