假设我有x
值、y
值和预期 y 值f
(来自一些非线性最佳拟合曲线)。
如何计算 R 中的 R^2?请注意,此函数不是线性模型,而是非线性最小二乘 ( nls
) 拟合,因此不是lm
拟合。
您只需使用该lm
函数来拟合线性模型:
x = runif(100)
y = runif(100)
spam = summary(lm(x~y))
> spam$r.squared
[1] 0.0008532386
请注意,r 平方不是为非线性模型定义的,或者至少非常棘手,引用自 R-help:
有一个很好的理由是,适合 R 的 nls 模型不提供 r-squared - r-squared 对于一般的 nls 模型没有意义。
r-squared 的一种思考方式是将拟合模型的残差平方和与仅由常数组成的平凡模型的残差平方和进行比较。在处理 nls 模型时,您不能保证这是嵌套模型的比较。如果模型没有嵌套,那么这种比较就没有多大意义。
所以答案是你可能一开始就不想这样做。
如果您想要同行评审的证据,请参阅这篇文章;并不是说您无法计算 R^2 值,只是它可能与线性模型案例中的含义不同/具有相同的理想属性。
听起来 f 是您的预测值。所以它们到实际值的距离除以 n * y 的方差
所以像
1-sum((y-f)^2)/(length(y)*var(y))
应该给你一个准 rsquared 值,只要你的模型相当接近线性模型并且 n 相当大。
作为对所提出问题的直接回答(而不是争论 R2/伪 R2 没有用),包nagelkerke
中的函数rcompanion
将报告 McFadden、Cox 和 Snell 提出的非线性最小二乘 (nls) 模型的各种伪 R2 值,和 Nagelkerke,例如
require(nls)
data(BrendonSmall)
quadplat = function(x, a, b, clx) {
ifelse(x < clx, a + b * x + (-0.5*b/clx) * x * x,
a + b * clx + (-0.5*b/clx) * clx * clx)}
model = nls(Sodium ~ quadplat(Calories, a, b, clx),
data = BrendonSmall,
start = list(a = 519,
b = 0.359,
clx = 2304))
nullfunct = function(x, m){m}
null.model = nls(Sodium ~ nullfunct(Calories, m),
data = BrendonSmall,
start = list(m = 1346))
nagelkerke(model, null=null.model)
该soilphysics
包还报告 Efron 的伪 R2 和调整后的模型的伪 R2 值为nls
1 - RSS/TSS:
pred <- predict(model)
n <- length(pred)
res <- resid(model)
w <- weights(model)
if (is.null(w)) w <- rep(1, n)
rss <- sum(w * res ^ 2)
resp <- pred + res
center <- weighted.mean(resp, w)
r.df <- summary(model)$df[2]
int.df <- 1
tss <- sum(w * (resp - center)^2)
r.sq <- 1 - rss/tss
adj.r.sq <- 1 - (1 - r.sq) * (n - int.df) / r.df
out <- list(pseudo.R.squared = r.sq,
adj.R.squared = adj.r.sq)
这也是由包中的函数pseudo R2
计算得出的。基本上,这个 R2 衡量你的合身程度与你只画一条穿过它们的水平线相比有多好。如果您的空模型是一个只允许拦截模型的模型,这对模型是有意义的。同样对于特定的其他非线性模型,它也很有意义。例如,对于使用严格递增样条曲线的骗局模型(样条曲线术语中的 bs="mpi"),最坏可能情况(例如,您的数据严格递减)的拟合模型将是一条平线,因此会导致一个accuracy
rcompanion
nls
R2
为零。调整后的 R2 也会惩罚拟合参数 nrs 较高的模型。使用调整后的 R2 值已经解决了上面链接的论文的许多批评,http ://www.ncbi.nlm.nih.gov/pmc/articles/PMC2892436/(如果有人发誓使用信息标准来做模型选择问题变成了使用哪一个 - AIC、BIC、EBIC、AICc、QIC 等)。
只是使用
r.sq <- max(cor(y,yfitted),0)^2
adj.r.sq <- 1 - (1 - r.sq) * (n - int.df) / r.df
我认为如果你有正常的高斯误差也是有意义的 - 即观察到的和拟合的 y 之间的相关性(剪裁为零,因此负关系意味着零预测能力)平方,然后针对拟合参数的 nr 进行调整调整后的版本。如果y
和yfitted
朝着相同的方向前进,这将是常规线性模型报告的R2
和值。adjusted R2
对我来说,这至少是完全合理的,所以我不同意完全拒绝模型pseudo R2
值的有用性,nls
正如上面的答案似乎暗示的那样。
对于非正常错误结构(例如,如果您使用带有非正常错误的 GAM),McFadden pseudo R2
类似地定义为
1-residual deviance/null deviance
非线性模型的另一个准 R 平方是对实际 y 值和预测 y 值之间的相关性进行平方。对于线性模型,这是常规的 R 平方。
modelr
包modelr::rsquare(nls_model, data)
nls_model <- nls(mpg ~ a / wt + b, data = mtcars, start = list(a = 40, b = 4))
modelr::rsquare(nls_model, mtcars)
# 0.794
rcompanion
这给出了与 Tom 从资源中描述的更长的方式基本相同的结果。
nagelkerke
功能更长的路nullfunct <- function(x, m){m}
null_model <- nls(mpg ~ nullfunct(wt, m),
data = mtcars,
start = list(m = mean(mtcars$mpg)))
nagelkerke(nls_model, null_model)[2]
# 0.794 or 0.796
lm(mpg ~ predict(nls_model), data = mtcars) %>% broom::glance()
# 0.795
就像他们说的,这只是一个近似值。
作为这个问题的替代方案,我多次使用以下过程:
向所有人致以最良好的祝愿。帕特里克。