我在75%的数据集上拟合了一个线性回归模型,包括~11000个观测值和143个变量:
gl.fit <- lm(y[1:ceiling(length(y)*(3/4))] ~ ., data= x[1:ceiling(length(y)*(3/4)),]) #3/4 for training
,我的R ^ 2为0.43。然后,我尝试使用其余数据预测我的测试数据:
ytest=y[(ceiling(length(y)*(3/4))+1):length(y)]
x.test <- cbind(1,x[(ceiling(length(y)*(3/4))+1):length(y),]) #The rest for test
yhat <- as.matrix(x.test)%*%gl.fit$coefficients #Calculate the predicted values
我现在想计算测试数据的R ^ 2值。有没有简单的方法来计算?
谢谢
这里有几个问题。首先,这不是使用lm(...)
的好方法。 lm(...)
用于数据框,公式表达式引用df中的列。所以,假设你的数据是两个向量x
和y
,
set.seed(1) # for reproducible example
x <- 1:11000
y <- 3+0.1*x + rnorm(11000,sd=1000)
df <- data.frame(x,y)
# training set
train <- sample(1:nrow(df),0.75*nrow(df)) # random sample of 75% of data
fit <- lm(y~x,data=df[train,])
现在,fit
拥有基于训练集的模型。以这种方式使用lm(...)
可以让您在没有所有矩阵乘法的情况下生成预测。
第二个问题是R平方的定义。 conventional definition是:
1 - SS.residuals / SS.total
对于训练集和训练集,
SS.total = SS.regression + SS.residual
所以
SS.regression = SS.total - SS.residual,
因此
R.sq = SS.regression / SS.total
所以R.sq是由模型解释的数据集中的可变性部分,并且总是在0和1之间。
你可以在下面看到这个。
SS.total <- with(df[train,],sum((y-mean(y))^2))
SS.residual <- sum(residuals(fit)^2)
SS.regression <- sum((fitted(fit)-mean(df[train,]$y))^2)
SS.total - (SS.regression+SS.residual)
# [1] 1.907349e-06
SS.regression/SS.total # fraction of variation explained by the model
# [1] 0.08965502
1-SS.residual/SS.total # same thing, for model frame ONLY!!!
# [1] 0.08965502
summary(fit)$r.squared # both are = R.squared
# [1] 0.08965502
但这不适用于测试集(例如,当您从模型进行预测时)。
test <- -train
test.pred <- predict(fit,newdata=df[test,])
test.y <- df[test,]$y
SS.total <- sum((test.y - mean(test.y))^2)
SS.residual <- sum((test.y - test.pred)^2)
SS.regression <- sum((test.pred - mean(test.y))^2)
SS.total - (SS.regression+SS.residual)
# [1] 8958890
# NOT the fraction of variability explained by the model
test.rsq <- 1 - SS.residual/SS.total
test.rsq
# [1] 0.0924713
# fraction of variability explained by the model
SS.regression/SS.total
# [1] 0.08956405
在这个人为的例子中没有太大区别,但很可能有一个R平方。值小于0(以这种方式定义)。
例如,如果模型是具有测试集的非常差的预测器,那么残差实际上可能大于测试集中的总变化。这相当于说使用它的平均值比使用从训练集派生的模型更好地建模测试集。
我注意到你使用数据的前四分之三作为训练集,而不是随机抽样(如本例所示)。如果y
对x
的依赖性是非线性的,并且x
是有序的,那么你可以得到一个带有测试集的负R-sq。
关于下面的OP评论,使用测试集评估模型的一种方法是通过比较模型内和模型外均方误差(MSE)。
mse.train <- summary(fit)$sigma^2
mse.test <- sum((test.pred - test.y)^2)/(nrow(df)-length(train)-2)
如果我们假设训练和测试集都是正态分布,具有相同的方差并且具有遵循相同模型公式的均值,那么该比率应该具有带有(n.train-2)和(n.test-)的F分布。 2)自由度。如果MSE基于F检验显着不同,则该模型不能很好地拟合测试数据。
你有没有绘制你的test.y和pred.y vs x?仅这一点就会告诉你很多。
计算测试数据的R平方有点棘手,因为你必须记住你的基线是什么。您的基线预测是您的训练数据的平均值。
因此,扩展上面的@jlhoward提供的示例:
SS.test.total <- sum((test.y - mean(df[train,]$y))^2)
SS.test.residual <- sum((test.y - test.pred)^2)
SS.test.regression <- sum((test.pred - mean(df[train,]$y))^2)
SS.test.total - (SS.test.regression+SS.test.residual)
# [1] 11617720 not 8958890
test.rsq <- 1 - SS.test.residual/SS.test.total
test.rsq
# [1] 0.09284556 not 0.0924713
# fraction of variability explained by the model
SS.test.regression/SS.test.total
# [1] 0.08907705 not 0.08956405
更新:miscTools::rSquared()
函数假设R平方是在训练模型的同一数据集上计算的,因为它计算
yy <- y - mean(y)
在这里的第184行幕后:https://github.com/cran/miscTools/blob/master/R/utils.R
如果你想要一个函数,miscTools
包有一个rSquared
函数。
require(miscTools)
r2 <- rSquared(ytest, resid = ytest-yhat)
当您对(非)样本使用R2度量时,您会忽略对R2的解释的某些方面:
如果你想使用R,我会推荐功能modelr::rsquare
。请注意,这使用了测试样本中的SSR总数,而不是训练样本(正如一些人似乎所倡导的那样)。
在这里,我举一个例子,我们的列车数据只有3分,因此我们有很高的风险,我们有一个糟糕的模型,因此样本表现不佳,实际上,你可以看到R2是负的!
library(modelr)
train <- mtcars[c(1,3,4),]
test <- mtcars[-c(1,3,4),]
mod <- lm(carb ~ drat, data = train)
计算列车数据:
## train
y_train <- train$carb
SSR_y_train <- sum((y_train-mean(y_train))^2)
cor(fitted(mod), y_train)^2
#> [1] 0.2985092
rsquare(mod, train)
#> [1] 0.2985092
1-sum(residuals(mod)^2)/SSR_y_train
#> [1] 0.2985092
计算测试数据:
## test
pred_test <- predict(mod, newdata = test)
y_test <- test$carb
SSR_y_test <- sum((y_test-mean(y_test))^2)
cor(pred_test, y_test)^2
#> [1] 0.01737236
rsquare(mod, test)
#> [1] -0.6769549
1- 28* var(pred_test-y_test)/SSR_y_train
#> [1] -19.31621
1- 28* var(pred_test-y_test)/SSR_y_test
#> [1] -0.6769549