编写循环/函数以生成各种混合模型结果

Question

我正在用 R 编写循环或函数，但我还没有真正理解如何做到这一点。目前，我需要编写一个循环/函数（不确定哪个更好）来在同一数据框中创建多个混合模型的结果。

示例数据集如下所示：

dataset <- read.table(text = 
"ID A_2 B_2 C_2 A_1 B_1 C_1 chkgp
M1  10  20  60  30  54  33  Treatment
M1  20  50  40  33  31  44  Placebo
M2  40  80  40  23  15  66  Placebo
M2  30  90  40  67  67  66  Treatment
M3  30  10  20  22  89  77  Treatment
M3  40  50  30  44  50  88  Placebo
M4  40  30  40  42  34  99  Treatment
M4  30  40  50  33  60  80  Placebo",header = TRUE, stringsAsFactors = FALSE)

我建立了一个模型，以Variable _2为因变量，variable _1为自变量，请有“lme4”包来运行混合模型

dep_vars<-grep("_2$",colnames(dataset),value = T) #This selects all variables ending in `_2` which should all be dependent variables.


#This removes the `_2` from the dependent variables which should give you the common stem which can be used to select both dependent and independent variables from your data frame.

reg_vars<-gsub("_2$","",dep_vars)

## To check that we have exact the correct variable which _2
dep_vars


[1] "A_2" "B_2" "C_2"

创建一个循环获取所有结果

full_results <- lapply(reg_vars, function(i) summary(lmer(paste0("log(",i,"_2)~",i,"_1+chkgp+(1|ID)"),data=dataset)))

查看第一个模型结果的摘要

full_results[1]
[[1]]
Linear mixed model fit by REML ['lmerMod']
Formula: log(A_2) ~ A_1 + chkgp + (1 | ID)
   Data: dataset

REML criterion at convergence: 16.9

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.16981 -0.24161  0.04418  0.37744  0.95925 

Random effects:
 Groups   Name        Variance Std.Dev.
 ID       (Intercept) 0.1314   0.3625  
 Residual             0.1188   0.3446  
Number of obs: 8, groups:  ID, 4

Fixed effects:
                Estimate Std. Error t value
(Intercept)     3.293643   0.411924   7.996
A_1             0.004512   0.009844   0.458
chkgpTreatment -0.276792   0.253242  -1.093

Correlation of Fixed Effects:
            (Intr) A_1   
A_1         -0.795       
chkgpTrtmnt -0.068 -0.272

问题：我想得到每个模型的 chkgpTreatment std.error t 值、P 值、上 CI 和下 CI 的结果，并将其存储在这样的数据框中

Depend_var  Indep_var mean difference
                                       p.value Upper ci Lower_ci
A_2 A_1 chkgpTreatment          
B_2 B_1 chkgpTreatment          
C_2 C_1 chkgpTreatment

Answer 1

我认为这会产生您正在寻找的输出。为了获得 p 值，我必须安装 lmerTest 包。此解决方案还使用 purrr 和 dplyr 包，因此如果您想使用它，则需要安装它们。

library(lmerTest)
library(purrr)
library(dplyr)

dataset <- read.table(
   text =
      "ID A_2 B_2 C_2 A_1 B_1 C_1 chkgp
M1  10  20  60  30  54  33  Treatment
M1  20  50  40  33  31  44  Placebo
M2  40  80  40  23  15  66  Placebo
M2  30  90  40  67  67  66  Treatment
M3  30  10  20  22  89  77  Treatment
M3  40  50  30  44  50  88  Placebo
M4  40  30  40  42  34  99  Treatment
M4  30  40  50  33  60  80  Placebo",
   header = TRUE,
   stringsAsFactors = FALSE
)
#This selects all variables ending in `_2` which should all be dependent
#variables.
dep_vars <-
   grep("_2$", colnames(dataset), value = T) 


#This removes the `_2` from the dependent variables which should give you the
#common stem which can be used to select both dependent and independent
#variables from your data frame.

reg_vars <- gsub("_2$", "", dep_vars)

## To check that we have exact the correct variable which _2
dep_vars

full_results <-
   map(reg_vars, function(i) {
      summary(lmerTest::lmer(paste0("log(", i, "_2)~", i, "_1+chkgp+(1|ID)"),
                   data = dataset))
   })



results <- map_dfr(full_results, function(.x) {
   data.frame( 
      # extract indep_var name and replace "1" with "2"
      depend_var = paste0(substr(row.names(.x$coefficients)[2], 1, 2), "2"),
      # extract depend_var name
      indep_var = row.names(.x$coefficients)[2],
      # Get the coefficient associated with chkgpTrtmnt
      mean_difference = .x$coefficients[3, 1],
      # Get std. error
      se = .x$coefficients[3, 2],
      # Get p-value
      p.value = .x$coefficients[3, 5],
      # Calculate the CI by +/- 1.96 * the standard error
      lower_ci = (.x$coefficients[3, 1] - (1.96 * .x$coefficients[3, 4])),
      upper_ci = (.x$coefficients[3, 1] + (1.96 * .x$coefficients[3, 4]))
   )
})