Регрессия с многомерными векторами и извлечением коэффициентов
Я хочу создать 1000 выборок из 200 двумерных нормально распределенных векторов.
set.seed(42) # for sake of reproducibility
mu <- c(1, 1)
S <- matrix(c(0.56, 0.4,
0.4, 1), nrow=2, ncol=2, byrow=TRUE)
bivn <- mvrnorm(200, mu=mu, Sigma=S)
так что я могу запустить регрессии OLS для каждой выборки и, следовательно, получить 1000 оценок. я пробовал это
library(MASS)
bivn_1000 <- replicate(1000, mvrnorm(200, mu=mu, Sigma=S), simplify=FALSE)
но я застрял там, потому что теперь я не знаю, как продолжить регрессию для каждой выборки.
Я был бы признателен за помощь, чтобы узнать, как запустить эти 1000 регрессий, а затем извлечь коэффициенты.
Попробуйте sapply(bivn_1000, function(x) summary(lm(x[, 1] ~ x[, 2]))$coef). Ваши данные двумерные, или я ошибаюсь?
да. Каждый образец содержит 200 двумерных векторов, и я пытаюсь создать 1000 образцов.
Пожалуйста, смотрите мой ответ ниже.
Ответы 1
Мы могли бы написать пользовательскую функцию регрессии.
regFun1 <- function(x) summary(lm(x[, 1] ~ x[, 2]))
который мы можем перебирать данные с помощью lapply:
l1 <- lapply(bivn_1000, regFun1)
Коэффициенты сохраняются внутри списка и могут быть извлечены следующим образом:
l1[[1]]$coefficients # for the first regression
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 0.5554601 0.06082924 9.131466 7.969277e-17
# x[, 2] 0.4797568 0.04255711 11.273246 4.322184e-23
Редактировать:
Если нам нужны только оценки без статистики, мы соответствующим образом корректируем вывод функции.
regFun2 <- function(x) summary(lm(x[, 1] ~ x[, 2]))$coef[, 1]
Поскольку нам может понадобиться вывод в матричной форме, мы используем sapply далее.
m2 <- t(sapply(bivn_1000, regFun2))
head(m2)
# (Intercept) x[, 2]
# [1,] 0.6315558 0.4389721
# [2,] 0.5514555 0.4840933
# [3,] 0.6782464 0.3250800
# [4,] 0.6350999 0.3848747
# [5,] 0.5899311 0.3645237
# [6,] 0.6263678 0.3825725
куда
dim(m2)
# [1] 1000 2
уверяет нас, что у нас есть 1000 оценок.
Да, это сработало. Но есть ли возможность получить 1000 оценщиков, по 1 на выборку?
Я не уверен, правильно ли я вас понял. В списке 1000 оценщиков. Попробуйте length(l1), что дает 1000.
Да, вы правы, я этого не заметил. Мне нужны только коэффициенты для первого столбца. Я попробовал этот цикл for(i in 1:1000){ a<-l1[[i]]$coefficients[,1] b_hat<-rbind() } Но он вставляет только последний коэффициент. Я тоже пробовал этоb_hat<-l1[[c(1:1000)]]$coefficients[,1] но произошла ошибка. Не могли бы вы помочь мне и с этим? Спасибо
Я только что попробовал, но получил описательную статистику и ряд NA :(