从 lm() 预测 'mlm' 线性模型对象
Prediction of 'mlm' linear model object from `lm()`
我有三个数据集：
响应 - 5（样本）x 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 中的每一行都有预测吗？
非常感谢对此问题的任何帮助，
马丁
您正在进入 R 中支持不佳的部分。您拥有的模型 class 是 "mlm"，即 "multiple linear models"，这不是标准 "lm" class。当您对一组共同的协变量/预测变量有几个（独立的）响应变量时，您就会得到它。虽然 lm() 函数可以拟合这样的模型，但是 predict 方法对 "mlm" class 很差。如果您查看 methods(predict)，您会看到一个 predict.mlm*。通常对于具有 "lm" class 的线性模型，当您调用 predict 时会调用 predict.lm；但是对于 "mlm" class 会调用 predict.mlm*。
predict.mlm*太原始了。它不允许 se.fit，即它不能产生预测误差、置信度/预测区间等，尽管这在理论上是可能的。它只能计算预测均值。如果是这样，我们为什么要使用 predict.mlm* 呢？！预测均值可以通过一个简单的矩阵-矩阵乘法得到（在标准中"lm" class这是一个矩阵-向量乘法），所以我们可以自己做。
考虑这个小的重现示例。
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
注意我使用data.frame()时的I()。当我们想要获得一个矩阵的数据帧时，这是必须的。你可以比较一下两者的区别：
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().
嗯，我建议在这种情况下使用 "list" 而不是 "data frame"。 data lm 中的参数] 以及 predict 中的 newdata 参数允许列表输入。 "list" 是比数据框更通用的结构，它可以毫无困难地容纳任何数据结构。我们可以做到：
## 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 相比有很大的不同。
我有三个数据集：
响应 - 5（样本）x 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 中的每一行都有预测吗？
非常感谢对此问题的任何帮助， 马丁
您正在进入 R 中支持不佳的部分。您拥有的模型 class 是 "mlm"，即 "multiple linear models"，这不是标准 "lm" class。当您对一组共同的协变量/预测变量有几个（独立的）响应变量时，您就会得到它。虽然 lm() 函数可以拟合这样的模型，但是 predict 方法对 "mlm" class 很差。如果您查看 methods(predict)，您会看到一个 predict.mlm*。通常对于具有 "lm" class 的线性模型，当您调用 predict 时会调用 predict.lm；但是对于 "mlm" class 会调用 predict.mlm*。
predict.mlm*太原始了。它不允许 se.fit，即它不能产生预测误差、置信度/预测区间等，尽管这在理论上是可能的。它只能计算预测均值。如果是这样，我们为什么要使用 predict.mlm* 呢？！预测均值可以通过一个简单的矩阵-矩阵乘法得到（在标准中"lm" class这是一个矩阵-向量乘法），所以我们可以自己做。
考虑这个小的重现示例。
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
注意我使用data.frame()时的I()。当我们想要获得一个矩阵的数据帧时，这是必须的。你可以比较一下两者的区别：
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().
嗯，我建议在这种情况下使用 "list" 而不是 "data frame"。 data lm 中的参数] 以及 predict 中的 newdata 参数允许列表输入。 "list" 是比数据框更通用的结构，它可以毫无困难地容纳任何数据结构。我们可以做到：
## 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 相比有很大的不同。