Прогнозирование объекта линейной модели mlm из lm()
У меня есть три набора данных:
ответ - матрица 5(выборки) х 10(зависимые переменные)
предикторы - матрица 5(выборки) x 2(независимые переменные)
test_set - матрица из 10(выборок) x 10(зависимые переменные, определенные в ответе)
response <- matrix(sample.int(15, size = 5*10, replace = TRUE), nrow = 5, ncol = 10)
colnames(response) <- c("1_DV","2_DV","3_DV","4_DV","5_DV","6_DV","7_DV","8_DV","9_DV","10_DV")
predictors <- matrix(sample.int(15, size = 7*2, replace = TRUE), nrow = 5, ncol = 2)
colnames(predictors) <- c("1_IV","2_IV")
test_set <- matrix(sample.int(15, size = 10*2, replace = TRUE), nrow = 10, ncol = 2)
colnames(test_set) <- c("1_IV","2_IV")
Я делаю многомерную линейную модель с использованием обучающего набора, определенного как комбинация наборов ответов и предикторов, и я хотел бы использовать эту модель для прогнозирования тестового набора:
training_dataframe <- data.frame(predictors, response)
fit <- lm(response ~ predictors, data = training_dataframe)
predictions <- predict(fit, data.frame(test_set))
Тем не менее, результаты для прогнозов действительно странные:
predictions
Сначала размеры матрицы составляют 5 x 10, то есть количество выборок в переменной отклика по числу DV.
Я не очень разбираюсь в анализе такого типа в R, но разве мне не нужно получать матрицу 10 x 10, чтобы у меня были прогнозы для каждой строки в моем test_set?
Любая помощь с этим вопросом будет принята с благодарностью, Мартин
1 ответ
Вы входите в плохо поддерживаемую часть в R. У вас есть класс модели "mlm", то есть "несколько линейных моделей", который не является стандартным классом "lm". Вы получаете это, когда у вас есть несколько (независимых) переменных ответа для общего набора ковариат / предикторов. Хотя lm() функция может соответствовать такой модели, predict Метод плох для класса "MLM". Если вы посмотрите на methods(predict), вы бы увидели predict.mlm*, Обычно для линейной модели с классом "lm", predict.lm называется, когда вы звоните predict; но для класса "MLM" predict.mlm* называется.
predict.mlm* слишком примитивно Не позволяет se.fitон не может давать ошибки прогнозирования, доверительные интервалы / интервалы прогнозирования и т. д., хотя это теоретически возможно. Это может только вычислить среднее значение предсказания. Если так, почему мы хотим использовать predict.mlm* совсем?! Среднее значение предсказания может быть получено с помощью тривиального умножения матрицы на матрицу (в стандартном классе "lm" это умножение матрицы на вектор), поэтому мы можем сделать это самостоятельно.
Посмотрите на этот маленький пример.
set.seed(0)
## 2 response of 10 observations each
response <- matrix(rnorm(20), 10, 2)
## 3 covariates with 10 observations each
predictors <- matrix(rnorm(30), 10, 3)
fit <- lm(response ~ predictors)
class(fit)
# [1] "mlm" "lm"
beta <- coef(fit)
# [,1] [,2]
#(Intercept) 0.5773235 -0.4752326
#predictors1 -0.9942677 0.6759778
#predictors2 -1.3306272 0.8322564
#predictors3 -0.5533336 0.6218942
Когда у вас есть набор данных прогноза:
# 2 new observations for 3 covariats
test_set <- matrix(rnorm(6), 2, 3)
сначала нужно заполнить столбец перехвата
Xp <- cbind(1, test_set)
Затем сделайте эту матрицу умножения
pred <- Xp %*% beta
# [,1] [,2]
#[1,] -2.905469 1.702384
#[2,] 1.871755 -1.236240
Возможно, вы заметили, что я даже не использовал фрейм данных здесь. Да, это не нужно, поскольку у вас есть все в матричной форме. Для тех волшебников R, возможно, используя lm.fit или даже qr.solve это более просто.
Но как полный ответ, это необходимо, чтобы продемонстрировать, как использовать predict.mlm чтобы получить желаемый результат.
## still using previous matrices
training_dataframe <- data.frame(response = I(response), predictors = I(predictors))
fit <- lm(response ~ predictors, data = training_dataframe)
newdat <- data.frame(predictors = I(test_set))
pred <- predict(fit, newdat)
# [,1] [,2]
#[1,] -2.905469 1.702384
#[2,] 1.871755 -1.236240
Обратите внимание I() когда я использую data.frame(), Это необходимо, когда мы хотим получить матрицу данных. Вы можете сравнить разницу между:
str(data.frame(response = I(response), predictors = I(predictors)))
#'data.frame': 10 obs. of 2 variables:
# $ response : AsIs [1:10, 1:2] 1.262954.... -0.32623.... 1.329799.... 1.272429.... 0.414641.... ...
# $ predictors: AsIs [1:10, 1:3] -0.22426.... 0.377395.... 0.133336.... 0.804189.... -0.05710.... ...
str(data.frame(response = response, predictors = predictors))
#'data.frame': 10 obs. of 5 variables:
# $ response.1 : num 1.263 -0.326 1.33 1.272 0.415 ...
# $ response.2 : num 0.764 -0.799 -1.148 -0.289 -0.299 ...
# $ predictors.1: num -0.2243 0.3774 0.1333 0.8042 -0.0571 ...
# $ predictors.2: num -0.236 -0.543 -0.433 -0.649 0.727 ...
# $ predictors.3: num 1.758 0.561 -0.453 -0.832 -1.167 ...
Без I() чтобы защитить матричный ввод, данные являются беспорядочными. Удивительно, что это не вызовет проблем с lm, но predict.mlm будет трудно получить правильную матрицу для прогноза, если вы не используете I(),
Ну, я бы рекомендовал использовать "список" вместо "фрейма данных" в этом случае. data аргумент в lm также newdata аргумент в predict позволяет вводить список. "Список" - это более общая структура, чем фрейм данных, который может без труда содержать любую структуру данных. Мы можем:
## still using previous matrices
training_list <- list(response = response, predictors = predictors)
fit <- lm(response ~ predictors, data = training_list)
newdat <- list(predictors = test_set)
pred <- predict(fit, newdat)
# [,1] [,2]
#[1,] -2.905469 1.702384
#[2,] 1.871755 -1.236240
Возможно, в самом конце, я должен подчеркнуть, что всегда безопасно использовать интерфейс формулы, а не матричный интерфейс. Я буду использовать R встроенный набор данных trees в качестве воспроизводимого примера.
fit <- lm(cbind(Girth, Height) ~ Volume, data = trees)
## use the first two rows as prediction dataset
predict(fit, newdata = trees[1:2, ])
# Girth Height
#1 9.579568 71.39192
#2 9.579568 71.39192
Возможно, вы до сих пор помните мое высказывание, что predict.mlm* слишком примитивен, чтобы поддерживать se.fit, Это шанс проверить это.
predict(fit, newdata = trees[1:2, ], se.fit = TRUE)
#Error in predict.mlm(fit, newdata = trees[1:2, ], se.fit = TRUE) :
# the 'se.fit' argument is not yet implemented for "mlm" objects
Упс... Как насчет доверительных интервалов / интервалов прогнозирования (фактически без возможности вычислить стандартную ошибку невозможно произвести эти интервалы)? Что ж, predict.mlm* просто проигнорирую это.
predict(fit, newdata = trees[1:2, ], interval = "confidence")
# Girth Height
#1 9.579568 71.39192
#2 9.579568 71.39192
Так что это так отличается по сравнению с predict.lm,