#r #function
#r #функция
Вопрос:
У меня есть такие данные, как:
data = data.frame(Y1 = runif(10000),
E1=runif(10000),
E2=runif(10000),
E3=runif(10000),
AQ=runif(10000),
WE=runif(10000),
SZ=runif(10000),
PO=runif(10000),
LL=runif(10000),
SCHOOL=sample(1:2, r = T),
CLASS = sample(1:4, r = T))
Мои цели — соответствовать этим простым регрессионным моделям:
Y1 = Ei AQ
Y1 = Ei AQ WE SZ
Y1 = Ei AQ WE SZ PO LL
для каждой комбинации ШКОЛЫ и КЛАССА это дает 3 * 2 * 4 = 24 МОДЕЛИ.
Из всех моделей я просто хочу сохранить коэффициент и значение p, а также доверительный интервал для всех независимых переменных.
Комментарии:
1. Что вы пробовали, что не сработало? Можете ли вы опубликовать свои попытки / проблемы?
2. Внимание:
SCHOOL=sample(1:2, r = T)вполне может привести к тому, что столбец будет полностью заполнен 2 или 1, поскольку выборкой может стать векторc(2,2)илиc(1,1). Если вам нужны случайные значения по всему набору данных, используйтеSCHOOL=sample(1:2, 10000, r = T). То же самое касается класса
Ответ №1:
Вы можете передавать регрессионные модели в lm() в виде символов. Итак, вам просто нужно перебирать имена ваших переменных и обновлять модель для каждого запуска:
results = list()
models = list()
model = "Y1 ~ "
for(i in 2:length(names(data))){
model = paste(model, ' ', names(data)[i])
models[i] = model
results[i] = list(lm(model,data=data))
}
Этот код расширяет уравнение регрессии на каждой итерации и сохраняет уравнение модели и результаты регрессии в списках.
РЕДАКТИРОВАТЬ: Извините, я неправильно прочитал первоначальный пост и предложенные модели регрессии, которые не были востребованы.
Приведенный ниже код вычисляет три различные модели регрессии для всех E-переменных и для всех классов и школ по отдельности:
counter = 1
models = c()
results = list()
for(s in unique(data$SCHOOL)){
for(c in unique(data$CLASS)){
subset = data[data$SCHOOL==s amp; data$CLASS==c,]
for(ei in 1:3){
for(dv in c("AQ", "AQ WE SZ", "AQ WE SZ PO LL")){
model = paste("Y1 ~ E",ei," ",dv,sep="")
models[counter] = paste(model,"for SCHOOL=",s," and CLASS=",c)
r = lm(model,data=subset)
results[counter]=list(r)
counter=counter 1
}
}
}
}
Третья модель в этом цикле будет:
> models[[3]]
[1] "Y1 ~ E1 AQ WE SZ PO LL for SCHOOL= 2 and CLASS= 3"
> summary(results[[3]])
Call:
lm(formula = model, data = subset)
Residuals:
Min 1Q Median 3Q Max
-0.53038 -0.24119 -0.00201 0.24940 0.54257
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.445346 0.024633 18.079 < 2e-16 ***
E1 0.035352 0.019676 1.797 0.072504 .
AQ 0.017344 0.019958 0.869 0.384919
WE -0.002737 0.019820 -0.138 0.890174
SZ 0.067423 0.020079 3.358 0.000797 ***
PO -0.029509 0.019897 -1.483 0.138188
LL 0.018326 0.019483 0.941 0.346988
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.2863 on 2493 degrees of freedom
Multiple R-squared: 0.007309, Adjusted R-squared: 0.00492
F-statistic: 3.059 on 6 and 2493 DF, p-value: 0.005518
Комментарии:
1. большое вам спасибо, это действительно очень полезно. Однако у меня есть вопрос. Ваш код включает E1, E2 и E3 в модели, но я надеюсь подогнать модели, как в вопросе, где E1, E2 и E3 представлены в отдельных моделях, а модели оцениваются группами.
2. Извините, я неправильно понял вашу проблему. Правка исправила это. Теперь она перебирает разные подмножества данных и разные модели.
Ответ №2:
Вот один из вариантов с tidyverse
library(dplyr) # >= 1.0.0
library(tidyr)
library(purrr)
library(broom)
lst1 <- list("AQ", c("AQ", "WE", "SZ"),
c("AQ", "WE", "SZ", "PO", "LL"))
nm1 <- grep("^E\d $", names(data), value = TRUE)
fmlst <- do.call(c, lapply(lst1, function(vec)
lapply(nm1, function(nm) reformulate(c(nm, vec), response = 'Y1'))))
data %>%
nest_by(SCHOOL, CLASS) %>%
summarise(lmmodels = map(fmlst, ~ lm(.x, data = data)),
tidyout = map(lmmodels, tidy))
# A tibble: 36 x 4
# Groups: SCHOOL, CLASS [4]
# SCHOOL CLASS lmmodels tidyout
# <int> <int> <list> <list>
# 1 1 2 <lm> <tibble [3 × 5]>
# 2 1 2 <lm> <tibble [3 × 5]>
# 3 1 2 <lm> <tibble [3 × 5]>
# 4 1 2 <lm> <tibble [5 × 5]>
# 5 1 2 <lm> <tibble [5 × 5]>
# 6 1 2 <lm> <tibble [5 × 5]>
# 7 1 2 <lm> <tibble [7 × 5]>
# 8 1 2 <lm> <tibble [7 × 5]>
# 9 1 2 <lm> <tibble [7 × 5]>
#10 1 4 <lm> <tibble [3 × 5]>
# … with 26 more rows
Комментарии:
1. большое вам спасибо, это действительно превосходно; как мне получить доступ к результатам или сохранить их? Я пытаюсь индексировать [], но не могу получить это .. в частности, я надеюсь на оценку, значение pvalue и доверительные интервалы в одном data.frame
2. @bvowe. Это
listстолбец. Таким образом, вы можете перебирать список с помощьюlapplyилиmapизpurrr, т.е. если ‘out’ является результатом, тоlapply(out$tidyout, function(x) x$yourcolumnname)3. я пробовал, TEST=data %>% nest_by(ШКОЛА, КЛАСС) %>% summarise(lmmodels = map(fmlst, ~ lm(.x, data = данные)), tidyout = map(lmmodels, tidy)) lapply (ПРОВЕРКА $ tidyout, функция (x) x $ коэффициенты)
4. @bvowe вы можете сделать
map(TEST$tidyout, ~ .x$estimate)илиmap(TEST$tidyout, ~ .x$p.value)
Ответ №3:
Я бы решил это так
library(tidyverse)
data <- data.frame(
Y1 = runif(10000),
E1 = runif(10000),
E2 = runif(10000),
E3 = runif(10000),
AQ = runif(10000),
WE = runif(10000),
SZ = runif(10000),
PO = runif(10000),
LL = runif(10000),
SCHOOL = sample(1:2, size = 10000, r = T),
CLASS = sample(1:4, size = 10000, r = T)
)
xx <- data %>%
nest_by(SCHOOL, CLASS) %>%
mutate(
model1 = list(lm(Y1 ~ E1 AQ, data = data)),
model2 = list(lm(Y1 ~ E1 AQ WE SZ, data = data)),
model3 = list(lm(Y1 ~ E1 AQ WE SZ PO LL, data = data))
) %>%
mutate(across(contains("model"), .fns = ~ list(broom::tidy(.x))))
xx$model1
#> [[1]]
#> # A tibble: 3 x 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.511 0.0215 23.8 5.04e-103
#> 2 E1 -0.0188 0.0277 -0.680 4.97e- 1
#> 3 AQ -0.0152 0.0279 -0.546 5.85e- 1
#>
#> [[2]]
#> # A tibble: 3 x 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.496 0.0209 23.7 2.05e-103
#> 2 E1 0.0239 0.0273 0.875 3.81e- 1
#> 3 AQ -0.0162 0.0274 -0.591 5.54e- 1
#>
#> [[3]]
#> # A tibble: 3 x 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.489 0.0217 22.6 9.05e-95
#> 2 E1 -0.00647 0.0288 -0.225 8.22e- 1
#> 3 AQ 0.0315 0.0290 1.09 2.76e- 1
#>
#> [[4]]
#> # A tibble: 3 x 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.521 0.0221 23.6 1.21e-102
#> 2 E1 -0.0264 0.0278 -0.949 3.43e- 1
#> 3 AQ -0.00902 0.0287 -0.315 7.53e- 1
#>
#> [[5]]
#> # A tibble: 3 x 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.521 0.0213 24.5 1.49e-108
#> 2 E1 0.00735 0.0282 0.260 7.95e- 1
#> 3 AQ -0.0257 0.0280 -0.917 3.59e- 1
#>
#> [[6]]
#> # A tibble: 3 x 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.490 0.0218 22.5 3.92e-94
#> 2 E1 0.0359 0.0285 1.26 2.09e- 1
#> 3 AQ -0.00878 0.0288 -0.305 7.61e- 1
#>
#> [[7]]
#> # A tibble: 3 x 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.529 0.0217 24.4 3.04e-107
#> 2 E1 -0.0498 0.0287 -1.74 8.28e- 2
#> 3 AQ -0.0341 0.0289 -1.18 2.38e- 1
#>
#> [[8]]
#> # A tibble: 3 x 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.532 0.0229 23.3 5.57e-99
#> 2 E1 -0.0760 0.0296 -2.56 1.04e- 2
#> 3 AQ -0.00109 0.0286 -0.0382 9.70e- 1
xx$model2
#> [[1]]
#> # A tibble: 5 x 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.507 0.0293 17.3 5.37e-60
#> 2 E1 -0.0183 0.0277 -0.660 5.09e- 1
#> 3 AQ -0.0143 0.0279 -0.514 6.07e- 1
#> 4 WE -0.0288 0.0283 -1.02 3.09e- 1
#> 5 SZ 0.0378 0.0287 1.32 1.89e- 1
#>
#> [[2]]
#> # A tibble: 5 x 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.504 0.0288 17.5 8.12e-62
#> 2 E1 0.0238 0.0273 0.871 3.84e- 1
#> 3 AQ -0.0164 0.0274 -0.599 5.49e- 1
#> 4 WE -0.0135 0.0275 -0.489 6.25e- 1
#> 5 SZ -0.00165 0.0280 -0.0589 9.53e- 1
#>
#> [[3]]
#> # A tibble: 5 x 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.486 0.0308 15.8 3.03e-51
#> 2 E1 -0.00598 0.0288 -0.208 8.36e- 1
#> 3 AQ 0.0316 0.0290 1.09 2.77e- 1
#> 4 WE 0.0135 0.0287 0.469 6.39e- 1
#> 5 SZ -0.00838 0.0286 -0.293 7.69e- 1
#>
#> [[4]]
#> # A tibble: 5 x 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.516 0.0299 17.3 4.36e-60
#> 2 E1 -0.0255 0.0279 -0.916 3.60e- 1
#> 3 AQ -0.00908 0.0287 -0.317 7.52e- 1
#> 4 WE -0.00593 0.0278 -0.213 8.31e- 1
#> 5 SZ 0.0157 0.0280 0.559 5.76e- 1
#>
#> [[5]]
#> # A tibble: 5 x 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.512 0.0293 17.5 1.93e-61
#> 2 E1 0.00714 0.0282 0.253 8.00e- 1
#> 3 AQ -0.0261 0.0280 -0.933 3.51e- 1
#> 4 WE 0.0464 0.0282 1.64 1.01e- 1
#> 5 SZ -0.0283 0.0283 -1.00 3.16e- 1
#>
#> [[6]]
#> # A tibble: 5 x 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.500 0.0299 16.7 1.73e-56
#> 2 E1 0.0356 0.0286 1.25 2.13e- 1
#> 3 AQ -0.00939 0.0289 -0.325 7.45e- 1
#> 4 WE -0.0184 0.0290 -0.633 5.27e- 1
#> 5 SZ -0.000915 0.0291 -0.0314 9.75e- 1
#>
#> [[7]]
#> # A tibble: 5 x 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.494 0.0307 16.1 6.12e-53
#> 2 E1 -0.0497 0.0287 -1.73 8.35e- 2
#> 3 AQ -0.0297 0.0290 -1.02 3.06e- 1
#> 4 WE 0.0210 0.0295 0.714 4.76e- 1
#> 5 SZ 0.0452 0.0287 1.58 1.15e- 1
#>
#> [[8]]
#> # A tibble: 5 x 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.516 0.0306 16.9 2.41e-57
#> 2 E1 -0.0760 0.0296 -2.57 1.04e- 2
#> 3 AQ -0.00122 0.0286 -0.0425 9.66e- 1
#> 4 WE 0.0396 0.0289 1.37 1.70e- 1
#> 5 SZ -0.00644 0.0292 -0.220 8.26e- 1
xx$model3
#> [[1]]
#> # A tibble: 7 x 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.515 0.0355 14.5 3.63e-44
#> 2 E1 -0.0184 0.0278 -0.662 5.08e- 1
#> 3 AQ -0.0143 0.0279 -0.513 6.08e- 1
#> 4 WE -0.0286 0.0283 -1.01 3.12e- 1
#> 5 SZ 0.0374 0.0288 1.30 1.94e- 1
#> 6 PO -0.00520 0.0282 -0.185 8.53e- 1
#> 7 LL -0.0117 0.0280 -0.419 6.75e- 1
#>
#> [[2]]
#> # A tibble: 7 x 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.493 0.0348 14.1 2.24e-42
#> 2 E1 0.0244 0.0273 0.893 3.72e- 1
#> 3 AQ -0.0174 0.0274 -0.635 5.26e- 1
#> 4 WE -0.0124 0.0275 -0.452 6.52e- 1
#> 5 SZ -0.00223 0.0280 -0.0795 9.37e- 1
#> 6 PO 0.0317 0.0275 1.15 2.50e- 1
#> 7 LL -0.0103 0.0272 -0.380 7.04e- 1
#>
#> [[3]]
#> # A tibble: 7 x 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.469 0.0366 12.8 2.04e-35
#> 2 E1 -0.00503 0.0288 -0.175 8.61e- 1
#> 3 AQ 0.0302 0.0290 1.04 2.98e- 1
#> 4 WE 0.0132 0.0287 0.460 6.46e- 1
#> 5 SZ -0.00696 0.0286 -0.244 8.08e- 1
#> 6 PO -0.0194 0.0281 -0.691 4.89e- 1
#> 7 LL 0.0547 0.0281 1.95 5.18e- 2
#>
#> [[4]]
#> # A tibble: 7 x 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.525 0.0372 14.1 4.15e-42
#> 2 E1 -0.0258 0.0279 -0.923 3.56e- 1
#> 3 AQ -0.0101 0.0288 -0.353 7.24e- 1
#> 4 WE -0.00707 0.0278 -0.254 7.99e- 1
#> 5 SZ 0.0155 0.0280 0.552 5.81e- 1
#> 6 PO 0.0171 0.0279 0.613 5.40e- 1
#> 7 LL -0.0327 0.0282 -1.16 2.47e- 1
#>
#> [[5]]
#> # A tibble: 7 x 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.501 0.0366 13.7 5.28e-40
#> 2 E1 0.00872 0.0282 0.309 7.57e- 1
#> 3 AQ -0.0249 0.0280 -0.888 3.75e- 1
#> 4 WE 0.0455 0.0283 1.61 1.08e- 1
#> 5 SZ -0.0258 0.0283 -0.913 3.61e- 1
#> 6 PO -0.0287 0.0281 -1.02 3.07e- 1
#> 7 LL 0.0456 0.0284 1.60 1.09e- 1
#>
#> [[6]]
#> # A tibble: 7 x 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.514 0.0363 14.2 2.57e-42
#> 2 E1 0.0361 0.0285 1.26 2.06e- 1
#> 3 AQ -0.00899 0.0288 -0.312 7.55e- 1
#> 4 WE -0.0167 0.0290 -0.576 5.65e- 1
#> 5 SZ 0.000892 0.0291 0.0307 9.76e- 1
#> 6 PO -0.0644 0.0288 -2.24 2.56e- 2
#> 7 LL 0.0325 0.0291 1.12 2.64e- 1
#>
#> [[7]]
#> # A tibble: 7 x 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.484 0.0366 13.2 1.74e-37
#> 2 E1 -0.0508 0.0288 -1.76 7.85e- 2
#> 3 AQ -0.0302 0.0291 -1.04 3.00e- 1
#> 4 WE 0.0214 0.0295 0.725 4.68e- 1
#> 5 SZ 0.0459 0.0287 1.60 1.10e- 1
#> 6 PO 0.0203 0.0287 0.708 4.79e- 1
#> 7 LL 0.000378 0.0288 0.0131 9.90e- 1
#>
#> [[8]]
#> # A tibble: 7 x 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.503 0.0376 13.4 3.27e-38
#> 2 E1 -0.0761 0.0297 -2.57 1.04e- 2
#> 3 AQ -0.00115 0.0286 -0.0400 9.68e- 1
#> 4 WE 0.0400 0.0289 1.38 1.67e- 1
#> 5 SZ -0.00594 0.0293 -0.203 8.39e- 1
#> 6 PO 0.0208 0.0288 0.721 4.71e- 1
#> 7 LL 0.00232 0.0295 0.0786 9.37e- 1
Создано 2020-08-13 пакетом reprex (версия 0.3.0)