Написание цикла / функции для генерации различных результатов смешанной модели

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 каждой модели и сохранить его во фрейме данных, например

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]))
   )
})