预测从`LM 'MLM' 线性模型对象（）`的

Question

我有三个数据集：

响应 - 5（样品）矩阵X 10（因变量）

预测 - 5矩阵（样品）×2（自变量）

TEST_SET - 10（样品基质）×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") 

我正在使用训练集多元线性模型限定的S的响应和预测集的结合，我想用这个模型进行了测试集的预测：

training_dataframe <- data.frame(predictors, response) 
fit <- lm(response ~ predictors, data = training_dataframe) 
predictions <- predict(fit, data.frame(test_set)) 

然而，预测结果是非常奇怪：

predictions 

第一关矩阵的尺寸是5×10，这是响应变量中的采样数量除以DV的数量。

我对R中这种类型的分析技术不是很熟练，但不应该得到10 x 10的矩阵，这样我就可以在我的test_set中预测每行了吗？

这个问题的任何帮助，将不胜感激， 马丁

Answer 1

您踏入河中支持不佳一部分，你有模型类是“传销”，即“多元线性模型”，这不是标准的“lm”类。当有多个（独立的）响应变量用于一组常用协变量/预测变量时，您可以得到它。虽然lm()函数可以适合这样的模型，但对于“mlm”类，predict方法很差。如果你看看methods(predict)，你会看到一个predict.mlm*。通常对于具有“lm”级别的线性模型，当您拨打predict时会调用predict.lm;但是对于“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()为了保护矩阵输入，数据很混乱。如果您不使用I()，那么令人惊讶的是，这不会导致lm问题，但predict.mlm将很难获得正确的预测矩阵。

那么，我会建议在这种情况下使用“列表”而不是“数据框”。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。